...
| Code Block | ||||
|---|---|---|---|---|
| ||||
#!/usr/bin/env python #''' # Copyright 20152019 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_geopotentialwind_onspeed_mlheight # # Author (date) : Cristian Simarro (0920/1011/2015) # # Categorymodified: Xavi :Abellan COMPUTATION # # OneLineDesc : Computes geopotential on model levels # # Description(03/12/2018) - compatibilty with Python 3 Category : Computes geopotential on model levels. # : COMPUTATION OneLineDesc : Computes the u/v components of the wind at certain height Description : Computes computes the u/v components of the wind at certain Based on code from Nilsheight. WediFirst, it calculates the geopotential IFS documentation: #of each model https://software.ecmwf.int/wiki/display/IFS/CY41R1+Official+IFS+Documentation # level. Once the requested height is found between two model part III. Dynamics and numerical procedures # levels, the program will vertically interpolate the u/v optimised implementation by Dominique Lucas. # component values. ported to Python by Cristian Simarro # # Parameters : tq.grib Based on code from Nils Wedi, the IFS documentation: https://www.ecmwf.int/en/forecasts/documentation-and-support/changes-ecmwf-model/ifs-documentation - grib filepart withIII. allDynamics theand levelistnumerical ofprocedures t and q # optimised implementation by Dominique zlnspLucas.grib - grib file withported levelistto 1Python forby params z and lnsp #Cristian Simarro Parameters : -w wind -o output (optional) - nameheight in ofmeters theyou outputwant file (default='z_out.grib') # # Return Value : output (default='z_out.grib') #to know A fieldset of geopotential on model levels # # Dependencies : None # # Example Usage : # u/v components python compute_geopotential_on_ml.py tq.grib zlnsp.grib from numpy import * import sys,math,os import argparse from gribapi import * def main(t_q,z_lnsp,u_v,out_name,h): #some- checksgrib andfile information printingwith all the levelist of print "Using as input files:\n ",t_q,z_lnsp,u_v print "The result will be stored in:\n ",out_name if os.path.exists(out_name): t os.remove(out_name)and q fout = open(out_name,'w') ftmp = open(t_q) ftmp.close() uv.grib Rd = 287.06 RG = 9.80665 index_keys- = ["date","time","shortName","level","step"] values= {}grib file with all the levelists of pv = {} out = {} gid_out = {} values_plev = {} values_lev = {}u/v h=int(h) zlnsp = grib_index_new_from_file(z_lnsp,index_keys) zlnsp.grib iidtq = grib_index_new_from_file(t_q,index_keys) iiduv =- grib_index_new_from_file(u_v,index_keys) file with levelist 1 for #we need to get z and lnsp from the first level to do the calculations counter=0 for date in grib_index_get(zlnsp,'date'): grib_index_select(zlnsp,'date',date) params z and lnsp grib_index_select(iidtq,'date',date) -o output grib_index_select(iiduv,'date',date) for time in grib_index_get(zlnsp,'time'): (optional) - name of the output file grib_index_select(zlnsp,'time',time) grib_index_select(iidtq,'time',time) grib_index_select(iiduv,'time',time) grib_index_select(zlnsp,'level',1(default='uv_out_<wind>.grib') Return Value : output (default='uv_out_<wind>.grib') grib_index_select(zlnsp,'step',0) A fieldset the u/v components values for at grib_index_select(zlnsp,'shortName','z') the specified gid = grib_new_from_index(zlnsp) height Dependencies : None Example Usage : #gridType must be gridded, not spectral compute_wind_speed_height.py tq.grib zlnsp.grib uv.grib if grib_get(gid,"gridType",str) == "sh": print(sys.argv[0]+' [ERROR] fields must be gridded, not spectral') -w 100 ''' from __future__ import print_function import sys import argparse import numpy as np from eccodes import (codes_index_new_from_file, codes_index_get, codes_get, codes_index_select, codes_new_from_index, codes_set, sys.exit(1) # surface geopotentialcodes_index_add_file, codes_get_array, codes_get_values, values["z"] = grib_get_values(gid) codes_index_release, z_h = values["z"]codes_release, codes_set_values, pv = grib_get_array(gid,'pv') codes_write) R_D levelSizeNV = grib_get(gid,'NV',int)/2 -1= 287.06 R_G = 9.80665 def parse_args(): ''' Parse program arguments using ArgumentParser''' parser grib_release(gid) = argparse.ArgumentParser( description='Python tool to calculate forthe wind stepat in grib_index_get(iidtq,'step'):certain height') parser.add_argument('-w', '--height', required=True, type=int, z_h = values["z"] help='height to calculate #the heigh in geopotentialwind components') parser.add_argument('-o', '--output', help='name of the output file') my_z = h*RG + z_hparser.add_argument('t_q', metavar='tq.grib', type=str, z_f_prev = z_h help=('grib file with temperature(t) and humidity(q)', grib_index_select(iidtq,'step',step) 'for the model levels')) grib_index_select(iiduv,'step',step) parser.add_argument('z_lnsp', metavar='zlnsp.grib', type=str, for shortName in ["lnsp"]: help=('grib file with geopotential(z) and Logarithm', grib_index_select(zlnsp,'shortName',shortName) grib_index_select(zlnsp,'step',step)'of surface pressure(lnsp) for the ml=1')) parser.add_argument('u_v', metavar='uv.grib', type=str, gid help= grib_new_from_index(zlnsp) ('grib file with u and v component of wind for', if grib_get(gid,"gridType",str) == "sh": 'the model levels')) print(sys.argv[0]+' [ERROR] fields must be gridded, not spectral')args = parser.parse_args() if not args.output: args.output = 'uv_out_%sm.grib' % args.height return sys.exit(1)args def main(): '''Main function''' args = parse_args() print('Arguments: %s' % values[shortName] = grib_get_values(gid)", ".join( ['%s: %s' % (k, v) for k, v pv = grib_get_array(gid,'pv') in vars(args).items()])) fout = open(args.output, 'wb') index_keys = ['date', 'time', levelSizeNV = grib_get(gid,'NV',int)/2 -1'shortName', 'level', 'step'] idx = grib_release(gidcodes_index_new_from_file(args.z_lnsp, index_keys) # surface pressure sp = exp(values["lnsp"]) # get the coefficients for computing the pressures # how many levels are we computing? grib_index_select(iidtq,'shortName',"t") levelSize=max(grib_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 codes_index_add_file(idx, args.t_q) if 'u_v' in args: codes_index_add_file(idx, # For speed and file I/O, we perform the computations with numpy vectors instead args.u_v) # iterate date for date in codes_index_get(idx, 'date'): codes_index_select(idx, 'date', date) # ofiterate fieldsets.step for time in codes_index_get(idx, 'time'): codes_index_select(idx, #initialize values for the output'time', time) values = for param in ["u","v"]:get_initial_values(idx, keep_sample=True) if 'height' in args: grib_index_select(iiduv,'level',1) values['height'] = args.height grib_index_select(iiduv,'shortName',param) values['gh'] = args.height * R_G + values['z'] gid = grib_new_from_index(iiduv) if 'levelist' in args: gid_out[param]=grib_clone(gid) values['levelist'] = args.levelist # iterate step all but geopotential z grib_release(gid) which is always step 0 (an) for step in out[param]=zeros(sp.size)codes_index_get(idx, 'step'): found = [False for i in range(sp.size)] codes_index_select(idx, 'step', step) # surface pressure for lev in list(reversed(range(1,levelSize+1))): try: # select the levelist and retrieve the vaules of t and q values['sp'] = get_surface_pressure(idx) # t_level:production_step(idx, values, for tfout) # q_level: values for q except WrongStepError: if step grib_index_select(iidtq,'level',lev) != '0': grib_index_select(iidtq,'shortName',"t") raise gid = grib_new_from_index(iidtq) try: codes_release(values['sample']) t_level = grib_get_values(gid) except KeyError: pass codes_index_release(idx) fout.close() def grib_release(gid)get_initial_values(idx, keep_sample=False): '''Get the values of surface z, pv and number of levels ''' gribcodes_index_select(iidtqidx, 'shortNamelevel',"q" 1) codes_index_select(idx, 'step', 0) codes_index_select(idx, 'shortName', 'z') gid = gribcodes_new_from_index(iidtqidx) values = {} # surface geopotential q_levelvalues['z'] = gribcodes_get_values(gid) values['pv'] = codes_get_array(gid, 'pv') values['nlevels'] = codes_get(gid, 'NV', int) // 2 grib_release(gid) - 1 check_max_level(idx, values) if keep_sample: # compute moist temperaturevalues['sample'] = gid else: codes_release(gid) return t_level = t_level * (1.+0.609133*q_level) values def check_max_level(idx, values): '''Make sure we have all the levels required''' # computehow the pressures (on half-levels) many levels are we computing? max_level = max(codes_index_get(idx, 'level', int)) if Ph_lev max_level != A[lev-1] + (B[lev-1] * sp) values['nlevels']: print('%s [WARN] total levels should be: %d but it is %d' % if lev == 1: (sys.argv[0], values['nlevels'], max_level), dlogP = log(Ph_levplusone/0.1file=sys.stderr) values['nlevels'] = max_level def get_surface_pressure(idx): '''Get the surface pressure for date-time-step''' codes_index_select(idx, alpha = log(2'level', 1) codes_index_select(idx, 'shortName', 'lnsp') gid else:= codes_new_from_index(idx) if gid is None: raise WrongStepError() if codes_get(gid, 'gridType', str) dlogP = log(Ph_levplusone/Ph_lev)== 'sh': print('%s [ERROR] fields must be gridded, not spectral' % sys.argv[0], dP = Ph_levplusone-Ph_levfile=sys.stderr) sys.exit(1) # surface pressure sfc_p = np.exp(codes_get_values(gid)) codes_release(gid) 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 return sfc_p def get_ph_levs(values, level): '''Return the presure at a given level and the next''' a_coef = values['pv'][0:values['nlevels'] + 1] b_coef = values['pv'][values['nlevels'] + 1:] ph_lev = a_coef[level - 1] + (b_coef[level - 1] * values['sp']) ph_levplusone = a_coef[level] + (b_coef[level] * values['sp']) return 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''' # select the levelist and retrieve the vaules of t and q # t_level: values for t # q_level: values for q codes_index_select(idx, 'level', lev) codes_index_select(idx, 'shortName', 't') gid = codes_new_from_index(idx) t_level = codes_get_values(gid) codes_release(gid) codes_index_select(idx, 'shortName', 'q') gid = codes_new_from_index(idx) 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, ph_levplusone = get_ph_levs(values, lev) if lev == 1: dlog_p # store the result (z_f) in a field and add to the output fieldset = np.log(ph_levplusone / 0.1) alpha = np.log(2) else: dlog_p #= (add it to the front, not the end, because we are going 'backwards' # through the fields) for param in ["u","v"]: grib_index_select(iiduv,'level',lev) grib_index_select(iiduv,'shortName',param) gidp=grib_new_from_index(iiduv) values_lev[param] = grib_get_values(gidp) grib_release(gidp) if (lev < levelSize): grib_index_select(iiduv,'level',lev+1) #137 gidp=grib_new_from_index(iiduv) values_plev[param] = grib_get_values(gidp) grib_release(gidpnp.log(ph_levplusone / ph_lev) alpha = 1. - ((ph_lev / (ph_levplusone - ph_lev)) * dlog_p) t_level = t_level * R_D # 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 + (t_level * alpha) # z_h is the geopotential of 'half-levels' # integrate z_h to next half level z_h = z_h + (t_level * dlog_p) return z_h, z_f def production_step(idx, values, fout): '''Produce u/v interpolated from ML for a given height''' # 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. out = {} params = ('u', 'v') # we need to get z and lnsp from the first level to do the # calculations z_h = values['z'] # height in geopotential z_f_prev = z_h # initialize values for the output for param in params: out[param] = np.zeros(values['sp'].size) found = [False for i in range(values['sp'].size)] for lev in list(reversed(list(range(1, values['nlevels'] + 1)))): z_h, z_f = compute_z_level(idx, lev, values, z_h) # retrieve u/v params for the else:current level for param in params: codes_index_select(idx, 'level', lev) values_plev[param] = zeros(sp.size) # 136 codes_index_select(idx, 'shortName', param) for igid in arange(z_f_prev.size):= codes_new_from_index(idx) values[param] = codes_get_values(gid) if found[i]: continue codes_release(gid) if lev < values['nlevels']: if my_z[i] >= zcodes_f_prev[i] and my_z[i] < z_f[i]: index_select(idx, 'level', lev + 1) # 137 gid = codes_new_from_index(idx) found[i]=True values['prev' + param] = codes_get_values(gid) #print "wind %d for point %d found between ml %d(%lf) and %d(%lf)\n" %(my_z[i],i,lev+1,z_f_prev[i],lev,z_f[i])codes_release(gid) else: values['prev' + param] = np.zeros(values['sp'].size) for param# in ["u","v"]: search if the provided wind height converted to # geopotential (my_z) is between the current level (z_f) # out[param][i] = ((values_lev[param][i] * ( my_z[i]-zand the previous one (z_f_prev[i])) + (values_plev[param][i] * (z_f[i] - my_z[i]) )) /for (z_f[i] -in range(z_f_prev[i].size): if found[i]: z_f_prev=z_f continue for i in arange(sp.size): if values['gh'][i] >= z_f_prev[i] and values['gh'][i] < z_f[i]: if not found[i]: = True # when found, interpolate vertically print "point ",i,"not found..."to get the for param# value and store it in ["u","v"]:out[param] to be written # at grib_set_values(gid_out[param],out[param])the end for param in params: grib_write(gid_out[param],fout) res grib_release(gid_out= (((float(values[param][i]) * grib_index_release(iidtq) grib_index_release(iiduv) grib_index_release(zlnsp) fout.close() print("Done") if __name__ == "__main__": request_date=0 (values['gh'][i] - z_f_prev[i])) + request_time=0 wind = 100 parser = argparse.ArgumentParser( description='Python tool to calculate the Z of the model levels') parser.add_argument("-w","--wind", help="height to calculate the wind components",required=True) parser.add_argument("-o","--output", help="name of the output file") parser.add_argument('t_q', metavar='tq.grib', type=str, (float(values['prev' + param][i]) * (z_f[i] - values['gh'][i]))) / help='grib file with temperature(t) and humidity(q) for the model levels') (z_f[i] - z_f_prev[i])) parser.add_argument('z_lnsp', metavar='zlnsp.grib', type=str, out[param][i] = res help='grib file# withupdate geopotential(z) and Logarithm of surface pressure(lnsp) for the ml=1')z_f_prev z_f_prev = z_f parser.add_argument('u_v', metavar='uv.grib', type=str, # simple error check for i in range(values['sp'].size): if not found[i]: help='grib file with u and v component of wind for the model levels') args = parser.parse_args() print('point ', i, 'not found...') for# write fnamethe in (args.t_q,args.z_lnsp,args.u_v): values in the fout file for ifparam not os.path.isfile(fname)in params: print "[ERROR] file %s does not exist" %(fname) sys.exit(1) if args.wind:codes_set(values['sample'], 'shortName', param) codes_set(values['sample'], 'typeOfLevel', 'heightAboveGround') codes_set(values['sample'], 'level', values['height']) wind = args.wind codes_set_values(values['sample'], out[param]) out_name = 'uv_out_'+wind+'m.grib' if args.output: out_name=args.output #calling main functioncodes_write(values['sample'], fout) class WrongStepError(Exception): ''' Exception capturing wrong step''' pass if __name__ == '__main__': main(args.t_q,args.z_lnsp,args.u_v,out_name,wind) ) |