- 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