Page History
...
Code Block | ||||
---|---|---|---|---|
| ||||
#!/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)) |