Skip to end of metadata
Go to start of metadata
compute_geopotential_on_ml.py
#!/usr/bin/env python
#
# Copyright 2015 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
#
# 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 params z and lnsp
#                 -l levelist (optional) - slash '/' 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')
#                 A fieldset of geopotential 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 import *
 
def main(t_q,z_lnsp,out_name,levelist):
 
    #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
 
    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 = codes_index_new_from_file(t_q,index_keys)
 
    #we need to get z and lnsp from the first level to do the calculations
    counter=0
    for date in codes_index_get(zlnsp,'date'):
        codes_index_select(zlnsp,'date',date)
        codes_index_select(iidtq,'date',date)
        for time in codes_index_get(zlnsp,'time'):
            codes_index_select(zlnsp,'time',time)
            codes_index_select(iidtq,'time',time)
         
            codes_index_select(zlnsp,'level',1)
    
            codes_index_select(zlnsp,'step',0)
            codes_index_select(zlnsp,'shortName','z')
            gid = codes_new_from_index(zlnsp)
             
            # surface geopotential
            values["z"] = codes_get_values(gid)
            z_h = values["z"]
            pv = codes_get_array(gid,'pv')
            levelSizeNV = codes_get(gid,'NV',int)/2 -1
            codes_release(gid)
 
            for step in codes_index_get(iidtq,'step'):
                z_h = values["z"]
                codes_index_select(iidtq,'step',step)
                codes_index_select(zlnsp,'step',step)
                for shortName in ["lnsp"]:
                   codes_index_select(zlnsp,'shortName',shortName)
                   gid = codes_new_from_index(zlnsp)
                   if codes_get(gid,"gridType",str) == "sh":
                       print(sys.argv[0]+' [ERROR] fields must be gridded, not spectral')
                       sys.exit(1)
                   values[shortName] = codes_get_values(gid)
                   pv = codes_get_array(gid,'pv')
                   levelSizeNV = codes_get(gid,'NV',int)/2 -1
                   codes_release(gid)
    
                # surface pressure
                sp = exp(values["lnsp"])
 
                # get the coefficients for computing the pressures
                # how many levels are we computing?
                codes_index_select(iidtq,'shortName',"t")
                levelSize=max(codes_index_get(iidtq,"level",int))
                if levelSize != levelSizeNV:
                    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
                # accumulating the height as we go.
                # See the IFS documentation:
                # https://software.ecmwf.int/wiki/display/IFS/CY41R1+Official+IFS+Documentation
                # part III
                # For speed and file I/O, we perform the computations with numpy vectors instead
                # of fieldsets.
 
                for lev in list(reversed(range(1,levelSize+1))):
 
                    # select the levelist and retrieve the vaules of t and q
                    # t_level: values for t
                    # q_level: values for q
                    codes_index_select(iidtq,'level',lev)
                    codes_index_select(iidtq,'shortName',"t")
                    gid = codes_new_from_index(iidtq)
                    t_level = codes_get_values(gid)
                    #gid_out will be used as output, cloning to get the attributes
                    gid_out = codes_clone(gid)
                    codes_release(gid)
                    codes_index_select(iidtq,'shortName',"q")
                    gid = codes_new_from_index(iidtq)
                    q_level = codes_get_values(gid)
                    codes_release(gid)
 
                    # compute moist temperature
                    t_level = t_level * (1.+0.609133*q_level)
 
                    # compute the pressures (on half-levels)
                    Ph_lev = A[lev-1] + (B[lev-1] * sp)
 
                    if lev == 1:
                        dlogP = log(Ph_levplusone/0.1)
                        alpha = log(2)
                    else:
                        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):
            print "[ERROR] file %s does not exist" %(fname)
            sys.exit(1)
    if args.levelist:
        if args.levelist != "all":
            levelist=args.levelist
    if args.output:
        out_name=args.output
 
    #calling main function
    main(args.t_q,args.z_lnsp,out_name,levelist)



  • No labels