- Created by Unknown User (uscs), last modified on Mar 21, 2017
You are viewing an old version of this page. View the current version.
Compare with Current View Page History
« Previous Version 21 Next »
Introduction
Some users are interested on geopotential (z) of the different model levels (ml). This Python script will allow you to compute it. The starting point is the geopotential and the pressure (z and lnsp) on the surface. The output will be written in GRIB format, providing the geopotential in m^2/s^2 for each level. (Note that you can get the height in meters if you divide the geopotential by the gravity of Earth (9.80665 m/s^2))
- It requires Python and ecCodes
- Input format has to be gridded, not spectral
- The script is based on the Metview macro mvl_geopotential_on_ml, code from Nils Wedi, IFS documentation:
For this use case, you need to retrieve in GRIB format both temperature (t) and Specific humidity (q) for each model level. Besides, you need both surface geopotential (z) and logarithm of surface pressure (lnsp) for model level = 1. Both files have to be synchronized using the same class, stream, date, time, step and gridType. The script will iterate through the date/time/step parameters, allowing you to get the values for multiple times.
If you are using MARS to retrieve the data, we have prepared a script that will allow you to get the data and then call the Python script automatically. You can set the date, time, step, class, grid resolution, type, area... The script will create and execute the MARS requests to get (t,q,z,lnsp). You can find more information about the script here: compute_geopotential_on_ml.ksh
Download
| Python program |
#!/usr/bin/env python3
'''
Copyright 2023 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
Xavi Abellan (03/12/2018) - compatibilty with Python 3
Xavi Abellan (27/03/2020) - Corrected steps in output file
Xavi Abellan (05/08/2020) - Better handling of levels
Xavi Abellan (31/01/2023) - Better handling of steps
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 :
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, codes_index_get, codes_get,
codes_index_select, codes_new_from_index, codes_set,
codes_index_add_file, codes_get_array, codes_get_values,
codes_index_release, codes_release, codes_set_values,
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',
default='all')
parser.add_argument('-o', '--output', help='name of the output file',
default='z_out.grib')
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()
# Handle levelist posibilities
if args.levelist == 'all':
args.levelist = range(1, 138)
elif "to" in args.levelist.lower():
if "by" in args.levelist.lower():
args.levelist = args.levelist.split('/')
args.levelist = list(range(int(args.levelist[0]),
int(args.levelist[2]) + 1,
int(args.levelist[4])))
else:
args.levelist = args.levelist.split('/')
args.levelist = list(range(int(args.levelist[0]),
int(args.levelist[2]) + 1))
else:
args.levelist = [int(l) for l in args.levelist.split('/')]
return args
def main():
'''Main function'''
args = parse_args()
print('Arguments: %s' % ", ".join(
['%s: %s' % (k, v) for k, v in vars(args).items()]))
fout = open(args.output, 'wb')
index_keys = ['date', 'time', 'shortName', 'level', 'step']
idx = codes_index_new_from_file(args.z_lnsp, index_keys)
codes_index_add_file(idx, args.t_q)
if 'u_v' in args:
codes_index_add_file(idx, args.u_v)
values = None
# iterate date
for date in codes_index_get(idx, 'date'):
codes_index_select(idx, 'date', date)
# iterate time
for time in codes_index_get(idx, 'time'):
codes_index_select(idx, 'time', time)
for step in codes_index_get(idx, 'step'):
codes_index_select(idx, 'step', step)
if not values:
values = get_initial_values(idx, keep_sample=True)
if 'height' in args:
values['height'] = args.height
values['gh'] = args.height * R_G + values['z']
if 'levelist' in args:
values['levelist'] = args.levelist
# surface pressure
try:
values['sp'] = get_surface_pressure(idx)
production_step(idx, step, values, fout)
except WrongStepError:
if step != '0':
raise
try:
codes_release(values['sample'])
except KeyError:
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, '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
else:
codes_release(gid)
return values
def check_max_level(idx, values):
'''Make sure we have all the levels required'''
# how many levels are we computing?
max_level = max(codes_index_get(idx, 'level', int))
if max_level != values['nlevels']:
print('%s [WARN] total levels should be: %d but it is %d' %
(sys.argv[0], values['nlevels'], max_level),
file=sys.stderr)
values['nlevels'] = max_level
def get_surface_pressure(idx):
'''Get the 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 gid is None:
raise WrongStepError()
if codes_get(gid, 'gridType', str) == 'sh':
print('%s [ERROR] fields must be gridded, not spectral' % sys.argv[0],
file=sys.stderr)
sys.exit(1)
# surface pressure
sfc_p = np.exp(codes_get_values(gid))
codes_release(gid)
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)
if gid is None:
raise MissingLevelError('T at level {} missing from input'.format(lev))
t_level = codes_get_values(gid)
codes_release(gid)
codes_index_select(idx, 'shortName', 'q')
gid = codes_new_from_index(idx)
if gid is None:
raise MissingLevelError('Q at level {} missing from input'.format(lev))
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 = np.log(ph_levplusone / 0.1)
alpha = np.log(2)
else:
dlog_p = np.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, step, values, fout):
'''Compute z at half & full level for the given level, based 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']
codes_set(values['sample'], 'step', int(step))
for lev in sorted(values['levelist'], reverse=True):
try:
z_h, z_f = compute_z_level(idx, lev, values, z_h)
# store the result (z_f) in a field and add to the output
codes_set(values['sample'], 'level', lev)
codes_set_values(values['sample'], z_f)
codes_write(values['sample'], fout)
except MissingLevelError as e:
print('%s [WARN] %s' % (sys.argv[0], e),
file=sys.stderr)
class WrongStepError(Exception):
''' Exception capturing wrong step'''
pass
class MissingLevelError(Exception):
''' Exception capturing missing levels in input'''
pass
if __name__ == '__main__':
main()
|
| Wrapper script |
#!/bin/bash
#
# Copyright 2023 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 (01/07/2015)
#
# Category : COMPUTATION
#
# OneLineDesc : Retrieves the files from mars and computes geopotential on model levels
#
# Description : This script will retrieve the necessary fields (t,q,z,lnsp) from MARS. Then it will
# call the Python script compute_geopotential_on_ml.py to calculate the geopotential
# on models levels
#
# Parameters : -d|--date <date> MARS formatted keyword. Default Yesterday (20151012)
# Examples: 20100501, 20100501/20100502/20100503, 20100501/to/20100503
# -t|--time <time> MARS formatted keyword time. Default 0
# Examples: 00, 00/12
# -s|--step <step> MARS formatted keyword step (only valid if type=fc). Default 0
# Examples: 00, 00/03/06/09
# -c|--class <class> MARS formatted keyword class. Default od
# Examples: od, ei
# -g|--grid <grid> MARS formatted keyword grid. Default 1.5/1.5
# Examples: 0.75/0.75, 128, 640
# -y|--type <type> MARS formatted keyword type. Default an
# Examples: fc, an
# -r|--stream <stream> MARS formatted keyword stream. Default oper
# Examples: oper, enfo
# -a|--area <area> MARS formatted keyword area. Default global
# Examples: europe, 75/60/10/-20
#
# Return Value : tq_ml_[date]_[time].grib - with all the levelist of t and q
# zlnsp_ml__[date]_[time].grib - levelist 1 for params z and lnsp
# compute_geopotential_on_ml.ksh.log - the logs will be stored here unless (-v)
#
# Dependencies : compute_geopotential_on_ml.py to calculate the geopotential
#
# Example Usage :
# ./compute_geopotential_on_ml.ksh 20150601
# ./compute_geopotential_on_ml.ksh -t 00/12 -d 20150101/to/20150110 -g 128 -a europe
#
usage(){
print " ./compute_geopotential_on_ml.ksh - Retrieves the files from mars and computes geopotential on model levels"
print
print " If you want to get more information about the format of the MARS keywords see:"
print " https://software.ecmwf.int/wiki/display/UDOC/MARS"
print
print "Usage: ./compute_geopotential_on_ml.ksh [ options ]"
print
print "options"
print " -v Enable verbose mode"
print " -k Do not clean after retrieve"
print
print " -d|--date <date> MARS formatted keyword. Default Yesterday ($DATE)"
print " Examples: 20100501, 20100501/20100502/20100503, 20100501/to/20100503"
print " -t|--time <time> MARS formatted keyword time. Default $TIME"
print " Examples: 00, 00/12"
print " -s|--step <step> MARS formatted keyword step (only valid if type=fc). Default $STEP"
print " Examples: 00, 00/03/06/09"
print " -c|--class <class> MARS formatted keyword class. Default $CLASS"
print " Examples: od, ei"
print " -g|--grid <grid> MARS formatted keyword grid. Default $GRID"
print " Examples: 0.75/0.75, 128, 640"
print " -y|--type <type> MARS formatted keyword type. Default $TYPE"
print " Examples: fc, an"
print " -r|--stream <stream> MARS formatted keyword stream. Default $STREAM"
print " Examples: oper, enfo"
print " -a|--area <area> MARS formatted keyword area. Default $AREA"
print " Examples: europe, 75/60/10/-20"
print " -l|--levelist <levelist> Slash separated list of levels to write in the output file. Default $LEVELIST"
print " Examples: 1, 91, 1/2/3/4"
print " -e|--extra <extra> (special) Extra keyword that you want to append to the MARS request."
print " Examples: \"expect=3\""
}
VERBOSE=0
DATE=`date +%Y%m%d -d "-1 day"`
TIME=0
CLASS="od"
GRID="1.5/1.5"
STEP=0
CLEAN=1
STREAM="oper"
TYPE="an"
LEVELIST="all"
AREA="global"
OPTS="h(help)v(verbose)k(keep)d:(date)t:(time)c:(class)g:(grid)y:(type)s:(step)r:(stream)a:(area)l:(levelist)e:(extra)"
while getopts "$OPTS" optchar
do
case "$optchar" in
d) DATE=$OPTARG;echo "date set to $OPTARG";;
t) TIME=$OPTARG;echo "time set to $OPTARG";;
c) CLASS=$OPTARG;echo "class set to $OPTARG";;
g) GRID=$OPTARG;echo "grid set to $OPTARG";;
y) TYPE=$OPTARG;echo "type set to $OPTARG";;
s) STEP=$OPTARG;echo "step set to $OPTARG";;
r) STREAM=$OPTARG;echo "stream set to $OPTARG";;
a) AREA=$OPTARG;echo "area set to $OPTARG";;
l) LEVELIST="$OPTARG";echo "levelist set to $OPTARG";;
e) EXTRA="$OPTARG,";echo "extra set to $OPTARG";;
v) VERBOSE=1;;
k) CLEAN=0;;
h) usage;
exit 0;;
?) usage;
exit 0;;
esac
done
EXECUTABLE=`basename $0`
IFS="/"
set -A DSTRING $DATE
set -A TSTRING $TIME
IFS=""
FNAME="${DSTRING[0]}_${TSTRING[0]}"
if [[ ! -f compute_geopotential_on_ml.py ]]; then
echo "[ERROR] you need compute_geopotential_on_ml.py in the same directory"
exit 1
fi
if [[ $TYPE == "fc" ]]; then
STEP="step=$STEP,"
else
STEP=""
fi
cat << EOF1 > tq.req
retrieve,
date= $DATE,
class = $CLASS,
time= $TIME,
stream = $STREAM,
levtype= ml,
grid= $GRID,
$STEP
$EXTRA
type= $TYPE,
area= $AREA,
param= t/q,
levelist= all,
target="tq_ml_${FNAME}.grib"
EOF1
cat << EOF2 > lnsp.req
retrieve,
date= $DATE,
class = $CLASS,
time= $TIME,
stream = $STREAM,
levtype= ml,
grid= $GRID,
$STEP
$EXTRA
type= $TYPE,
area= $AREA,
param=lnsp,
levelist=1,
target="lnsp_ml_${FNAME}.grib"
EOF2
cat << EOF3 > z.req
retrieve,
date= $DATE,
class = $CLASS,
time= $TIME,
stream = $STREAM,
levtype= ml,
grid= $GRID,
$EXTRA
area= $AREA,
param=z,
levelist=1,
type=an,
target="z_ml_${FNAME}.grib"
EOF3
if [ $VERBOSE -eq 0 ]; then
rm -f ${EXECUTABLE}.log
fi
for i in tq.req lnsp.req z.req; do
echo "executing mars request for $i ..."
if [ $VERBOSE -eq 1 ]; then
cmnd="`which mars` $i"
else
cmnd="`which mars` $i >>${EXECUTABLE}.log"
fi
echo $cmnd
eval $cmnd
ret_code=$?
if [ $ret_code != 0 ]; then
printf "Error : [%d] when executing command: '$cmnd'\n" $ret_code
exit $ret_code
fi
done
if [[ -f "lnsp_ml_${FNAME}.grib" && -f "z_ml_${FNAME}.grib" ]]; then
cat lnsp_ml_${FNAME}.grib z_ml_${FNAME}.grib > zlnsp_ml_${FNAME}.grib
rm -f lnsp_ml_${FNAME}.grib z_ml_${FNAME}.grib
if [ $CLEAN -eq 1 ]; then
rm -f z.req tq.req lnsp.req
fi
fi
cmnd="python3 compute_geopotential_on_ml.py tq_ml_${FNAME}.grib zlnsp_ml_${FNAME}.grib -l \"$LEVELIST\""
echo $cmnd
eval $cmnd
ret_code=$?
if [ $ret_code != 0 ]; then
printf "Error : [%d] when executing command: '$cmnd'\n" $ret_code
exit $ret_code
fi
if [ $CLEAN -eq 1 ]; then
rm -f z.req tq.req lnsp.req zlnsp_ml_${FNAME}.grib tq_ml_${FNAME}.grib
fi
|
Examples
This example will compute the geopotential on the 2015-10-08 time 00 operational analysis model levels (137). Below you can see the MARS user documentation request used to retrieve both files. You can set a different class/stream/type for the input data. The gribType and resolution can also be changed.
python compute_geopotential_on_ml.py tq_ml_20151008_00.grib zlnsp_ml_20151008_00.grib python compute_geopotential_on_ml.py tq_ml_20151008_00.grib zlnsp_ml_20151008_00.grib -o my_grib.grib
tq_ml_20151008_00.grib
retrieve, date= 20150801, class = od, time= 0, stream = oper, levtype= ml, grid= 1.5/1.5, type= an, param= t/q, levelist= all, target="tq_ml_20150801_00.grib"
zlnsp_ml_20151008_00.grib
retrieve, date= 20150801, class = od, time= 0, stream = oper, levtype= ml, grid= 1.5/1.5, type= an, param=lnsp, levelist=1, target="zlnsp_ml_20151008_00.grib" retrieve, type= an, param=z, levelist=1, target="zlnsp_ml_20151008_00.grib"
- No labels