#!/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_geopotentialwind_onspeed_mlheight
#
# Author (date) : Cristian Simarro (0920/1011/2015)
#
# Category : COMPUTATION
#
# OneLineDesc : Computes geopotential on model levels the u/v components of the wind at certain height
#
# Description : Computes geopotential on model levelscomputes the u/v components of the wind at certain height.
# Based on code from Nils Wedi, the IFS documentation:First, it calculates the geopotential of each model level. Once the requested height
# https://software.ecmwf.int/wiki/display/IFS/CY41R1+Official+IFS+Documentation
# part III. Dynamics and numerical procedures
#is found between two model levels, the program will vertically interpolate the u/v component values.
# Based on code from optimisedNils implementationWedi, bythe DominiqueIFS Lucas.documentation:
# ported to Python by Cristian Simarro
#https://software.ecmwf.int/wiki/display/IFS/CY41R1+Official+IFS+Documentation
# Parameters : tq.grib part III. Dynamics and numerical -procedures
# grib file with all the levelist of t and q
# optimised implementation by Dominique Lucas.
# zlnsp.grib ported to -Python gribby file with levelist 1 for params z and lnsp
#Cristian Simarro
#
# Parameters : -w wind -o outputheight in meters (optional)you -want nameto ofknow the output file (default='z_out.grib')
#
# Return Value : output (default='z_out.grib')
#u/v components
# tq.grib A fieldset of geopotential on model levels
#
# Dependencies : None
#
# Example Usage :
# - grib file with all the levelist of t and q
# uv.grib 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 checks and information printing
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):
- grib file with all the levelists of u/v
# zlnsp.grib - grib file with levelist 1 for params z and lnsp
# os.remove(out_name)
-o output fout = open(out_name,'w')
ftmp = open(t_q)
ftmp.close()
Rd = 287.06
RG = 9.80665
index_keys = ["date","time","shortName","level","step"]
values= {}
pv = {}
out = {}
gid_out = {}
values_plev = {}
values_lev = {}
h=int(h)
zlnsp = grib_index_new_from_file(z_lnsp,index_keys)
iidtq = grib_index_new_from_file(t_q,index_keys)
iiduv = grib_index_new_from_file(u_v,index_keys)
#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)
grib_index_select(iidtq,'date',date)
grib_index_select(iiduv,'date',date(optional) - name of the output file (default='uv_out_<wind>.grib')
#
# Return Value : output (default='uv_out_<wind>.grib')
# A fieldset the u/v components values for at the specified height
#
# Dependencies : None
#
# Example Usage :
# python compute_wind_speed_height.py tq_ml_20151116_0.grib zlnsp_ml_20151116_0.grib uv_ml_20151116_0.grib -w 100
from numpy import *
import sys,math,os
import argparse
from gribapi import *
PARAMS = ["u","v"]
def main(t_q,z_lnsp,u_v,out_name,h):
#some checks and information printing
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):
os.remove(out_name)
fout = open(out_name,'w')
ftmp = open(t_q)
ftmp.close()
Rd = 287.06
RG = 9.80665
index_keys = ["date","time","shortName","level","step"]
values= {}
pv = {}
out = {}
gid_out = {}
values_plev = {}
values_lev = {}
h=int(h)
zlnsp = grib_index_new_from_file(z_lnsp,index_keys)
iidtq = grib_index_new_from_file(t_q,index_keys)
iiduv = grib_index_new_from_file(u_v,index_keys)
counter=0
# iterate date
for date in grib_index_get(zlnsp,'date'):
grib_index_select(zlnsp,'date',date)
grib_index_select(iidtq,'date',date)
grib_index_select(iiduv,'date',date)
# iterate step
for time in grib_index_get(zlnsp,'time'):
grib_index_select(zlnsp,'time',time)
grib_index_select(iidtq,'time',time)
grib_index_select(iiduv,'time',time)
grib_index_select(zlnsp,'level',1)
grib_index_select(zlnsp,'step',0)
grib_index_select(zlnsp,'shortName','z')
gid = grib_new_from_index(zlnsp)
#gridType must be gridded, not spectral
if grib_get(gid,"gridType",str) == "sh":
print(sys.argv[0]+' [ERROR] fields must be gridded, not spectral')
sys.exit(1)
# surface geopotential
values["z"] = grib_get_values(gid)
for time in grib_index_get(zlnsp,'time'):
z_h = values["z"]
pv = grib_indexget_selectarray(zlnspgid,'timepv',time)
levelSizeNV = grib_index_selectget(iidtqgid,'timeNV',timeint)/2 -1
grib_index_select(iiduv,'time',time)
grib_index_select(zlnsp,'level',1)
release(gid)
grib_index_select(zlnsp,'step',0)
# iterate step all but grib_index_select(zlnsp,'shortName','z')
geopotential z wich is always step 0 (an)
for gidstep =in grib_new_from_index(zlnsp)
index_get(iidtq,'step'):
#gridType must be gridded,# notwe spectral
need to get z and lnsp from the first level to do if grib_get(gid,"gridType",str) == "sh":
the calculations
z_h = print(sys.argv[0]+' [ERROR] fields must be gridded, not spectral')
values["z"]
# heigh in geopotential
sys.exit(1)
my_z = #h*RG surface+ geopotentialz_h
values["z"] z_f_prev = grib_get_values(gid)z_h
z_h = values["z"]
pv = grib_getindex_arrayselect(gidiidtq,'pvstep',step)
levelSizeNV = grib_index_getselect(gidiiduv,'NVstep',int)/2 -1step)
for shortName in ["lnsp"]:
grib_release(gid)
for step in grib_index_getselect(iidtqzlnsp,'stepshortName',shortName):
z_h = values["z"]
grib_index_select(zlnsp,'step',step)
# heigh in geopotential
gid my_z = h*RG + z_h
= grib_new_from_index(zlnsp)
if z_f_prevgrib_get(gid,"gridType",str) = z_h= "sh":
print(sys.argv[0]+' [ERROR] fields must be gridded, not spectral')
grib_index_select(iidtq,'step',step)
grib_index_select(iiduv,'step',stepsys.exit(1)
for shortName in values["lnsp"]:shortName] = grib_get_values(gid)
pv = grib_indexget_selectarray(zlnspgid,'shortNamepv',shortName)
levelSizeNV = grib_index_selectget(zlnspgid,'stepNV',stepint)/2 -1
gid = grib_new_from_index(zlnsprelease(gid)
if grib_get(gid,"gridType",str) == "sh":# surface pressure
sp print(sys.argv[0]+' [ERROR] fields must be gridded, not spectral')
= exp(values["lnsp"])
# get sys.exit(1)
the coefficients for computing the pressures
values[shortName] = grib_get_values(gid)
# how many levels are we computing?
pv = grib_getindex_arrayselect(gidiidtq,'pvshortName',"t")
levelSizeNV = levelSize=max(grib_index_get(gid,'NV'iidtq,"level",int)/2 -1
)
grib_release(gid)
if levelSize != levelSizeNV:
# surface pressure
print(sys.argv[0]+' [WARN] total levels should be: '+str(levelSizeNV)+' but it sp = exp(values["lnsp"])
is '+str(levelSize))
A = # get the coefficients for computing the pressures
pv[0:levelSize+1]
B = pv[levelSize+1:]
# how many levels are we computing?
Ph_levplusone = A[levelSize] + (B[levelSize]*sp)
grib_index_select(iidtq,'shortName',"t")
# We want to integrate up into levelSize=max(grib_index_get(iidtq,"level",int))the atmosphere, starting at the ground
if# levelSizeso != levelSizeNV:
we start at the lowest level (highest number) and keep
print(sys.argv[0]+' [WARN] total levels should# be: '+str(levelSizeNV)+' but it is '+str(levelSize))accumulating the height as we go.
# See Athe =IFS pv[0:levelSize+1]documentation:
B = pv[levelSize+1:]# https://software.ecmwf.int/wiki/display/IFS/CY41R1+Official+IFS+Documentation
Ph_levplusone# = A[levelSize] + (B[levelSize]*sp)
part III
# For Wespeed wantand to integrate up intofile I/O, we perform the atmosphere,computations startingwith atnumpy thevectors groundinstead
# so we start at the lowest level (highest number) and keep
of fieldsets.
# accumulating the#initialize heightvalues asfor wethe go.output
#for Seeparam thein IFS documentationPARAMS:
# https://software.ecmwf.int/wiki/display/IFS/CY41R1+Official+IFS+Documentation grib_index_select(iiduv,'level',1)
# part III
grib_index_select(iiduv,'shortName',param)
# For speed and file I/O, we perform the computations with numpy vectors instead
gid = grib_new_from_index(iiduv)
gid_out[param]=grib_clone(gid)
# of fieldsets.
grib_release(gid)
out[param]=zeros(sp.size)
#initialize values for the output
found = [False for parami in ["u","v"]:range(sp.size)]
for lev grib_index_select(iiduv,'level',1)in list(reversed(range(1,levelSize+1))):
grib_index_select(iiduv,'shortName',param)# select the levelist and retrieve the vaules of t and q
gid = grib_new_from_index(iiduv)# t_level: values for t
gid_out[param]=grib_clone(gid)# q_level: values for q
grib_index_release(gidselect(iidtq,'level',lev)
out[param]=zeros(sp.sizegrib_index_select(iidtq,'shortName',"t")
found = [False for igid in range(sp.size)]
= grib_new_from_index(iidtq)
for lev in list(reversed(range(1,levelSize+1))): t_level = grib_get_values(gid)
# select the levelist and retrieve the vaules of t and qgrib_release(gid)
grib_index_select(iidtq,'shortName',"q")
# t_level: values for t
gid = grib_new_from_index(iidtq)
# q_level: values for q = grib_get_values(gid)
grib_index_select(iidtq,'level',lev)
release(gid)
# grib_index_select(iidtq,'shortName',"t")compute moist temperature
gidt_level = grib_new_from_index(iidtq)
t_level * (1.+0.609133*q_level)
# compute the t_level = grib_get_values(gidpressures (on half-levels)
grib_release(gid)
Ph_lev = A[lev-1] + (B[lev-1] * sp)
grib_index_select(iidtq,'shortName',"q")
if lev == 1:
gid dlogP = grib_new_from_index(iidtq)
log(Ph_levplusone/0.1)
q_levelalpha = grib_get_values(gidlog(2)
grib_release(gid)
else:
# compute moist temperature
dlogP = log(Ph_levplusone/Ph_lev)
t_level = t_level * (1.+0.609133*q_level)
dP # compute the pressures (on half-levels)
= Ph_levplusone-Ph_lev
Ph_levalpha = A[lev-1]. +- ((B[lev-1] * spPh_lev/dP)*dlogP)
ifTRd lev == 1: t_level*Rd
# z_f is the dlogPgeopotential = log(Ph_levplusone/0.1)
of this full level
# integrate alphafrom =previous log(2lower)
half-level z_h to the full level
else:
z_f = z_h + (TRd*alpha)
dlogP = log(Ph_levplusone/Ph_lev)
# z_h is the geopotential of 'half-levels'
dP = Ph_levplusone-Ph_lev
# integrate z_h to next half level
alpha = 1. - ((Ph_lev/dP) z_h=z_h+(TRd*dlogP)
TRdPh_levplusone = tPh_level*Rdlev
# z_f is the geopotential of this full retrieve u/v params for the current level
#for integrateparam from previous (lower) half-level z_h to the full levelin PARAMS:
z_f = z_h + (TRd*alpha)
grib_index_select(iiduv,'level',lev) #136
# z_h is the geopotential of 'half-levels'
grib_index_select(iiduv,'shortName',param)
# integrate z_h to next half level
gidp=grib_new_from_index(iiduv)
z_h=z_h+(TRd*dlogP)
values_lev[param] = grib_get_values(gidp)
Ph_levplusone = Ph_levgrib_release(gidp)
# store the resultif (z_f) in a field and add to the output fieldset
lev < levelSize):
# (add it to the front, not the end, because we are going 'backwards'
grib_index_select(iiduv,'level',lev+1) #137
# through the fields)
gidp=grib_new_from_index(iiduv)
for param in ["u","v"]:
values_plev[param] = grib_indexget_select(iiduv,'level',levvalues(gidp)
grib_index_select(iiduv,'shortName',paramrelease(gidp)
gidp=grib_new_from_index(iiduv)
else:
values_levplev[param] = grib_get_values(gidpzeros(sp.size)
# search if the provided grib_release(gidp)wind height converted to geopotential (my_z) is between
# the current if (lev < levelSize):level (z_f) and the previous one (z_f_prev)
for i grib_index_select(iiduv,'level',lev+1) #137in arange(z_f_prev.size):
if found[i]: continue
gidp=grib_new_from_index(iiduv)
if my_z[i] >= z_f_prev[i] and valuesmy_plevz[parami] =< grib_get_values(gidp)z_f[i]:
grib_release(gidp)found[i]=True
else:
# when found, interpolate vertically to get the value
values_plev[param] = zeros(sp.size)
# store it in out[param] to be written at the end
for i in arange(z_f_prev.size):
for param if found[i]: continue
in PARAMS:
if my_z out[param][i] >= z_f_prev( (values_lev[param][i] and* ( my_z[i] < -z_f_prev[i]:
)) + \
found[i]=True
#print "wind %d for point %d found between ml %d(%lf) and %d(%lf)\n" %(my_zvalues_plev[param][i],i,lev+1, * (z_f_prev[i],lev,z_f - my_z[i])
) ) / \
for param in ["u","v"]:
out[param][i] = ((values_lev[param][i] * ( myz_zf[i] - z_f_prev[i])) + (values_plev[param][i] * (z_f[i] - my_z[i]) )) / (z_f[i] -
# update z_f_prev[i])prev
z_f_prev=z_f
# simple error check
for i in arange(sp.size):
if not found[i]:
print "point ",i,"not found..."
# write the values in the fout file
for param in ["u","v"]PARAMS:
grib_set_values(gid_out[param],out[param])
grib_write(gid_out[param],fout)
grib_release(gid_out[param])
grib_index_release(iidtq)
grib_index_release(iiduv)
grib_index_release(zlnsp)
fout.close()
print("Done")
if __name__ == "__main__":
request_date=0
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,
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')
parser.add_argument('u_v', metavar='uv.grib', type=str,
help='grib file with u and v component of wind for the model levels')
args = parser.parse_args()
for fname in (args.t_q,args.z_lnsp,args.u_v):
if not os.path.isfile(fname):
print "[ERROR] file %s does not exist" %(fname)
sys.exit(1)
if args.wind:
wind = args.wind
out_name = 'uv_out_'+wind+'m.grib'
if args.output:
out_name=args.output
#calling main function
main(args.t_q,args.z_lnsp,args.u_v,out_name,wind)
|