Page History

#!/usr/bin/env python
#'''
# Copyright 20152018 ECMWF.
#
# This software is licensed under the terms of the Apache Licence Version 2.0
# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0
#
# In applying this licence, ECMWF does not waive the privileges and immunities 
granted to it by
# virtue of its status as an intergovernmental organisation 
nor does it submit to any jurisdiction.
#
# **************************************************************************
# Function      : compute_geopotential_on_ml
#
# Author (date) : Cristian Simarro (09/10/2015)
# modified:       Cristian Simarro (20/03/2017) - migrated to eccodes
#
# Category      : COMPUTATION
#
# OneLineDesc   : Computes geopotential on model levels
#
# DescriptionXavi Abellan    : Computes geopotential on model levels.
# (03/12/2018) - compatibilty with Python 3

Category      : COMPUTATION

OneLineDesc   : Computes geopotential on model levels

Description   : Computes geopotential on model levels.
                 Based on code from Nils Wedi, the IFS documentation:
#                 https://software.ecmwf.int/wiki/display/IFS/CY41R1+Official+IFS+Documentation
#                 part III. Dynamics and numerical procedures
#                 optimised implementation by Dominique Lucas.
#                 ported to Python by Cristian Simarro
#
# Parameters    : tq.grib                - grib file with all the levelist
 of t and q
#                 zlnsp.grib             - grib file with levelist 1 for paramsof zt and lnsp
#q
                zlnsp.grib             -l grib file with levelist (optional)1 -for slashparams
 '/' separated list of levelist to store in the output
#                 -o output   (optional) - name of the output file (default='z_out.grib')
#
# Return Value  : output (default='z_out.grib')
#    z and lnsp
                -l levelist (optional) - slash '/' separated A fieldsetlist of geopotentiallevelist
 on model levels
#
# Dependencies  : None
#
# Example Usage :
#                     python compute_geopotential_on_ml.py tq.grib zlnsp.grib
 
from numpy import *
import sys,math,os
import argparse
from eccodes importto *
store 
def main(t_q,z_lnsp,out_name,levelist):
 in the output
    #some checks and information printing
    print "Using as input files:\n   ",t_q,z_lnsp
    print "The result will be stored in:\n   ",out_name
 
 -o output   (optional) - name of the output file
           if levelist != "":
        print "Will only store these levels: "+levelist
        levelist_selected=levelist.split("/")
    fout = open(out_name,'w')
    ftmp = open(t_q)
    total_messages=codes_count_in_file(ftmp)/2
    ftmp.close()
    Rd = 287.06
    index_keys = ["date","time","shortName","level","step"]
 
    values= {}
    pv = {}
     
    zlnsp = codes_index_new_from_file(z_lnsp,index_keys)
    iidtq = default='z_out.grib')

Return Value  : output (default='z_out.grib')
                A fieldset of geopotential on model levels

Dependencies  : None

Example Usage :
                compute_geopotential_on_ml.py tq.grib zlnsp.grib
'''
from __future__ import print_function
import sys
import argparse
import numpy as np
from eccodes import (codes_index_new_from_file(t_q,index_keys)
 , codes_index_get, codes_get,
    #we need to  get z and lnsp from the first level to do the calculations
  codes_index_select,  counter=0codes_new_from_index, codes_set,
    for date in codes_index_get(zlnsp,'date'):
        codes_index_select(zlnsp,'date',date)
        codes_index_select(iidtq,'date',date)add_file, codes_get_array, codes_get_values,
        for time in codes_index_get(zlnsp,'time'):
          codes_index_release, codes_release, codes_indexset_select(zlnsp,'time',time)values,
            codes_index_select(iidtq,'time',time)
         
            codes_index_select(zlnsp,'level',1)
    
            codes_index_select(zlnsp,'step',0)
            codes_write)

R_D = 287.06
R_G = 9.80665


def parse_args():
    ''' Parse program arguments using ArgumentParser'''
    parser = argparse.ArgumentParser(
        description='Python tool to calculate the Z of the model levels')
    parser.add_argument('-l', '--levelist', help='levelist to store',
         codes_index_select(zlnsp,'shortName','z')
            gid = codes_new_from_index(zlnsp  default='')
    parser.add_argument('-o', '--output', help='name of the     
 output file',
           # surface geopotential
            values["z"] = codes_get_values(giddefault='z_out.grib')
    parser.add_argument('t_q', metavar='tq.grib', type=str,
         z_h = values["z"]
            pv help= codes_get_array(gid,'pv')('grib file with temperature(t) and humidity(q)',
            levelSizeNV = codes_get(gid,'NV',int)/2 -1
                 'for the model levels'))
    codesparser.add_release(gid)
 
argument('z_lnsp', metavar='zlnsp.grib', type=str,
                for step in codes_index_get(iidtq,'step'):      help=('grib file with geopotential(z) and Logarithm',
                z_h  = values["z"]
           'of surface pressure(lnsp) for  codes_index_select(iidtq,'step',stepthe ml=1'))
    args = parser.parse_args()
    if args.levelist     codes_index_select(zlnsp,'step',step)== 'all':
        args.levelist = ''
    return  for shortName in ["lnsp"]args


def main():
    '''Main function'''
    args   = parse_args()

    print('Arguments: %s' %  codes_index_select(zlnsp,'shortName',shortName)", ".join(
        ['%s: %s' % (k, v) for k, v in vars(args).items()]))

  gid  fout = codes_new_from_index(zlnspopen(args.output, 'wb')
    index_keys = ['date', 'time', 'shortName', 'level', 'step']

    idx = codes_index_new_from_file(args.z_lnsp, index_keys)
   if codes_get(gid,"gridType",str) == "sh"_index_add_file(idx, args.t_q)
    if 'u_v' in args:
        codes_index_add_file(idx, args.u_v)

    # iterate date
    for date   print(sys.argv[0]+' [ERROR] fields must be gridded, not spectral'in codes_index_get(idx, 'date'):
        codes_index_select(idx, 'date', date)
        # iterate step
        for time in   sys.exit(1)codes_index_get(idx, 'time'):
            codes_index_select(idx, 'time', time)
            values[shortName] = codesget_getinitial_values(gididx, keep_sample=True)
            if 'height' in args:
    pv    = codes_get_array(gid,'pv')
       values['height'] = args.height
          levelSizeNV = codes_get(gid,'NV',int)/2 -1
   values['gh'] = args.height * R_G + values['z']
          codes_release(gid)
      if 'levelist' in args:
                # surface pressurevalues['levelist'] = args.levelist
            # iterate step all sp = exp(values["lnsp"]but geopotential z which is always step 0 (an)
 
           for step in codes_index_get(idx, 'step'):
 # get the coefficients for computing the pressures
        codes_index_select(idx, 'step', step)
      #  how many levels are we computing?
    # surface pressure
                 codes_index_select(iidtq,'shortName',"t"values['sp'] = get_surface_pressure(idx)
                levelSize=max(codes_index_get(iidtq,"level",int))production_step(idx, values, fout)

        try:
        if levelSize != levelSizeNV:
    codes_release(values['sample'])
        except KeyError:
        print(sys.argv[0]+' [WARN] total levels should be: '+str(levelSizeNV)+' but it is '+str(levelSize))
                A = pv[0:levelSize+1]
                B = pv[levelSize+1:]
                Ph_levplusone = A[levelSize] + (B[levelSize]*sp)
 
                # We want to integrate up into the atmosphere, starting at the ground
                # so we start at the lowest level (highest number) and keep
         pass

    codes_index_release(idx)
    fout.close()


def get_initial_values(idx, keep_sample=False):
    '''Get the values of surface z, pv and number of levels '''
    codes_index_select(idx, 'level', 1)
    codes_index_select(idx, 'step', 0)
    codes_index_select(idx, 'shortName', 'z')
    gid = codes_new_from_index(idx)

    values = {}
    # surface geopotential
    values['z'] = codes_get_values(gid)
    values['pv'] = codes_get_array(gid, 'pv')
    values['nlevels'] = codes_get(gid, 'NV', int) // 2 - 1
    check_max_level(idx, values)
    if keep_sample:
        values['sample'] = gid
 # accumulating the height as we go. else:
        codes_release(gid)
    return    # See the IFS documentationvalues


def check_max_level(idx, values):
    '''Make sure we have all the levels required'''
     # https://software.ecmwf.int/wiki/display/IFS/CY41R1+Official+IFS+Documentation
           how many levels are we computing?
    max_level = max(codes_index_get(idx, 'level', int))
    if # part IIImax_level != values['nlevels']:
        print('%s [WARN] total levels should be: %d but #it Foris speed%d' and%
 file I/O, we perform the   computations with numpy vectors instead
  (sys.argv[0], values['nlevels'], max_level),
            # of fieldsetsfile=sys.stderr)
 
       values['nlevels'] = max_level


def get_surface_pressure(idx):
    '''Get the for lev in list(reversed(range(1,levelSize+1))):
 
     surface pressure for date-time-step'''
    codes_index_select(idx, 'level', 1)
    codes_index_select(idx, 'shortName', 'lnsp')
    gid = codes_new_from_index(idx)
    if codes_get(gid, 'gridType', str) == 'sh':
 # select the levelist and retrieve the vaules of t and q
     print('%s [ERROR] fields must be gridded, not spectral' % sys.argv[0],
              file=sys.stderr)
 # t_level: values for t    sys.exit(1)
    # surface pressure
    sfc_p = np.exp(codes_get_values(gid))
    codes_release(gid)
    #return q_level: values for qsfc_p


def get_ph_levs(values, level):
    '''Return the presure at a given level and the next'''
    a_coef   codes_index_select(iidtq,'level',lev)
  = values['pv'][0:values['nlevels'] + 1]
    b_coef = values['pv'][values['nlevels'] + 1:]
    ph_lev = a_coef[level - 1] + codes_index_select(iidtq,'shortName',"t"(b_coef[level - 1] * values['sp'])
    ph_levplusone = a_coef[level] + (b_coef[level] * values['sp'])
    return      gid = codes_new_from_index(iidtq)ph_lev, ph_levplusone


def compute_z_level(idx, lev, values, z_h):
    '''Compute z at half & full level for the given level, based on t/q/sp'''
   t_level = codes_get_values(gid)
        # select the levelist and retrieve the vaules of t and q
    # t_level: values for t
    # #gidq_outlevel: willvalues befor usedq
 as output, cloning to get the attributes codes_index_select(idx, 'level', lev)
       codes_index_select(idx, 'shortName', 't')
    gid = codes_new_from_index(idx)
       gid_outt_level = codes_cloneget_values(gid)
                    codes_release(gid)
                    codes_index_select(iidtqidx, 'shortName'," 'q"')
                    gid = codes_new_from_index(iidtqidx)
                    q_level = codes_get_values(gid)
                    codes_release(gid)
 
    # compute moist temperature
    t_level = t_level * (1. +  0.609133  # compute moist temperature* q_level)

    # compute the pressures (on half-levels)
           t_levelph_lev, ph_levplusone = t_level * (1.+0.609133*q_level)
 get_ph_levs(values, lev)

    if lev == 1:
        dlog_p = np.log(ph_levplusone / 0.1)
     #  compute thealpha pressures= (on half-levels)np.log(2)
    else:
        dlog_p = np.log(ph_levplusone / ph_lev)
        Ph_levalpha = A[lev-1]. +- (B[lev-1](ph_lev / (ph_levplusone - ph_lev)) * spdlog_p)
 
    t_level = t_level * R_D

    # z_f is the geopotential of this  if lev == 1:full level
    # integrate from previous (lower) half-level z_h to the
    # full level
    z_f = dlogPz_h =+ log(Ph_levplusone/0.1)t_level * alpha)

    # z_h is the geopotential of 'half-levels'
    # integrate z_h to next half level
    alphaz_h = log(2)
    z_h + (t_level * dlog_p)

    return z_h, z_f


def production_step(idx,         elsevalues, fout):
    '''Compute z at half & full level for the given level, based         dlogP = log(Ph_levplusone/Ph_lev)
                        dP    = Ph_levplusone-Ph_lev
                        alpha = 1. - ((Ph_lev/dP)*dlogP)
 
                    TRd = t_level*Rd
 
                    # z_f is the geopotential of this full level
                    # integrate from previous (lower) half-level z_h to the full level
                    z_f = z_h + (TRd*alpha)
 
                    # z_h is the geopotential of 'half-levels'
                    # integrate z_h to next half level
                    z_h=z_h+(TRd*dlogP)
 
                    Ph_levplusone = Ph_lev
                    # store the result (z_f) in a field and add to the output fieldset
                    # (add it to the front, not the end, because we are going 'backwards'
                    # through the fields)
                    if levelist == "" or str(lev) in levelist_selected:
                        codes_set(gid_out,"paramId",129)
                        codes_set(gid_out,'generatingProcessIdentifier',128)
                        codes_set(gid_out,'level', lev)
                        #codes_set(gid_out,"bitsPerValue",24)
                        codes_set_values(gid_out,z_f)
                        codes_write(gid_out,fout)
                    counter += 1
                    if counter >= int((total_messages+1)/20):
                        sys.stdout.write('.')
                        sys.stdout.flush()
                        counter=0
                    codes_release(gid_out)
            
    fout.close()
    print("Done")
 
if __name__ == "__main__":
    request_date=0
    request_time=0
    levelist = ""
    out_name = 'z_out.grib'
 
    parser = argparse.ArgumentParser(
        description='Python tool to calculate the Z of the model levels')
    parser.add_argument("-l","--levelist", help="levelist to store")
    parser.add_argument("-o","--output", help="name of the output file")
    parser.add_argument('t_q', metavar='tq.grib', type=str,
                   help='grib file with temperature(t) and humidity(q) for the model levels')
    parser.add_argument('z_lnsp', metavar='zlnsp.grib', type=str,
                   help='grib file with geopotential(z) and Logarithm of surface pressure(lnsp) for the ml=1')
    args = parser.parse_args()
    for fname in (args.t_q,args.z_lnsp):
        if not os.path.isfile(fname):on t/q/sp'''
    # We want to integrate up into the atmosphere, starting at the
    # ground so we start at the lowest level (highest number) and
    # keep accumulating the height as we go.
    # See the IFS documentation, part III
    # For speed and file I/O, we perform the computations with
    # numpy vectors instead of fieldsets.

    z_h = values['z']
    for lev in list(reversed(list(range(1, values['nlevels'] + 1)))):
        z_h, z_f = compute_z_level(idx, lev, values, z_h)
        # store the  print "[ERROR] file %s does not exist" %(fname)result (z_f) in a field and add to the output
        if values['levelist'] == '' or sys.exitstr(1lev)
 in values['levelist']:
  if args.levelist:
        if args.levelist != "all": codes_set(values['sample'], 'level', lev)
            levelist=args.levelistcodes_set_values(values['sample'], z_f)
    if args.output:
        out_name=args.output
 
    #calling main functioncodes_write(values['sample'], fout)


if __name__ == '__main__':
    main(args.t_q,args.z_lnsp,out_name,levelist))