Changeset ca867de in flex_extract.git


Ignore:
Timestamp:
Oct 5, 2018, 3:35:18 PM (13 months ago)
Author:
Anne Philipp <anne.philipp@…>
Branches:
dev
Children:
5bad6ec
Parents:
27fe969
Message:

refactored functions in EcFlexpart? and did some minor changes

Files:
1 added
10 edited

Legend:

Unmodified
Added
Removed
  • run/control/CONTROL.temp

    r2fb99de rca867de  
    1515RIGHT 45000
    1616LEVEL 60
    17 LEVELIST 55/to/60
     17LEVELIST 40/to/60
    1818RESOL 63
    1919GAUSS 1
  • run/jobscripts/job.ksh

    r25b14be rca867de  
    3333  module load grib_api/1.14.5
    3434  module load emos/437-r64
    35   export PATH=${PATH}:${HOME}/flex_extract_v7.1/python
     35  export PATH=${PATH}:${HOME}/flex_extract_v7.1/source/python
    3636  ;;
    3737  *cca*)
     
    4141  module load python
    4242  export SCRATCH=$TMPDIR
    43   export PATH=${PATH}:${HOME}/flex_extract_v7.1/python
     43  export PATH=${PATH}:${HOME}/flex_extract_v7.1/source/python
    4444  ;;
    4545esac
     
    4949cd python$$
    5050
    51 export CONTROL=$PWD/CONTROL
     51export CONTROL=CONTROL
    5252
    5353cat >$CONTROL<<EOF
    54 accuracy 24
     54accuracy 16
    5555addpar 186 187 188 235 139 39
    5656area
    5757basetime None
    58 controlfile CONTROL.test
     58controlfile CONTROL.temp
    5959cwc 0
    6060date_chunk 3
     
    6868ectrans 1
    6969ecuid km4a
    70 end_date 20160606
     70end_date 20120908
    7171eta 0
    7272etadiff 0
     
    7878grib2flexpart 0
    7979grid 5000
    80 inputdir ../work
     80inputdir ../../run/workspace/test
    8181install_target None
    8282job_template job.temp
    83 left -10000
     83left -15000
    8484level 60
    85 levelist 59/to/60
     85levelist 40/to/60
    8686logicals gauss omega omegadiff eta etadiff dpdeta cwc wrf grib2flexpart ecstorage ectrans debug request
    8787lower 30000
     
    9494omega 0
    9595omegadiff 0
    96 outputdir ../work
    97 prefix EItest_
    98 queue ecgate
     96outputdir ../../run/workspace/test
     97ppid 41511
     98prefix EI
     99queue local
    99100request 1
    100101resol 63
    101 right 10000
     102right 45000
    102103smooth 0
    103 start_date 20160606
     104start_date 20120908
    104105step 00 01 02 03 04 05 00 07 08 09 10 11 00 01 02 03 04 05 00 07 08 09 10 11
    105106stream OPER
    106107time 00 00 00 00 00 00 06 00 00 00 00 00 12 12 12 12 12 12 18 12 12 12 12 12
    107108type AN FC FC FC FC FC AN FC FC FC FC FC AN FC FC FC FC FC AN FC FC FC FC FC
    108 upper 40000
     109upper 75000
    109110wrf 0
    110111EOF
  • run/run.sh

    • Property mode changed from 100644 to 100755
    r2fb99de rca867de  
    11#!/bin/bash
    2  
    3 pyscript = ../python/submit.py
     2#
     3# @Author: Anne Philipp
     4#
     5# @Date: October, 4 2018
     6#
     7# @Description:
     8#
    49
     10
     11# -----------------------------------------------------------------
     12# AVAILABLE COMMANDLINE ARGUMENTS TO SET
     13#
     14# THE USER HAS TO SPECIFY THESE PARAMETER
     15#
     16
     17QUEUE=None
     18START_DATE='20120908'
     19END_DATE=None
     20DATE_CHUNK=None
     21BASETIME=None
     22STEP=None
     23LEVELIST=None
     24AREA=None
     25INPUTDIR='../../run/workspace/test'
     26OUTPUTDIR=None
     27FLEXPART_ROOT_SCRIPTS=None
     28PP_ID=None
     29JOB_TEMPLATE='job.temp'
     30CONTROLFILE='CONTROL.temp'
     31DEBUG=1
     32REQUEST=1
     33
     34# -----------------------------------------------------------------
     35#
     36# AFTER THIS LINE THE USER DOES NOT HAVE TO CHANGE ANYTHING !!!
     37#
     38# -----------------------------------------------------------------
     39
     40# PATH TO SUBMISSION SCRIPT
     41pyscript=../source/python/submit.py
     42
     43# INITIALIZE EMPTY PARAMETERLIST
     44parameterlist=""
     45
     46# CHECK FOR MORE PARAMETER
     47if [ -n "$START_DATE" ]; then
     48  parameterlist+=" --start_date=$START_DATE"
     49fi
     50if [ -n "$END_DATE" ]; then
     51  parameterlist+=" --end_date=$END_DATE"
     52fi
     53if [ -n "$DATE_CHUNK" ]; then
     54  parameterlist+=" --date_chunk=$DATE_CHUNK"
     55fi
     56if [ -n "$BASETIME" ]; then
     57  parameterlist+=" --basetime=$BASETIME"
     58fi
     59if [ -n "$STEP" ]; then
     60  parameterlist+=" --step=$STEP"
     61fi
     62if [ -n "$LEVELIST" ]; then
     63  parameterlist+=" --levelist=$LEVELIST"
     64fi
     65if [ -n "$AREA" ]; then
     66  parameterlist+=" --area=$AREA"
     67fi
     68if [ -n "$INPUTDIR" ]; then
     69  parameterlist+=" --inputdir=$INPUTDIR"
     70fi
     71if [ -n "$OUTPUTDIR" ]; then
     72  parameterlist+=" --outputdir=$OUTPUTDIR"
     73fi
     74if [ -n "$FLEXPART_ROOT_SCRIPTS" ]; then
     75  parameterlist+=" --flexpart_root_scripts=$FLEXPART_ROOT_SCRIPTS"
     76fi
     77if [ -n "$PP_ID" ]; then
     78  parameterlist+=" --ppid=$PP_ID"
     79fi
     80if [ -n "$JOB_TEMPLATE" ]; then
     81  parameterlist+=" --job_template=$JOB_TEMPLATE"
     82fi
     83if [ -n "$QUEUE" ]; then
     84  parameterlist+=" --queue=$QUEUE"
     85fi
     86if [ -n "$CONTROLFILE" ]; then
     87  parameterlist+=" --controlfile=$CONTROLFILE"
     88fi
     89if [ -n "$DEBUG" ]; then
     90  parameterlist+=" --debug=$DEBUG"
     91fi
     92if [ -n "$REQUEST" ]; then
     93  parameterlist+=" --request=$REQUEST"
     94fi
     95
     96# -----------------------------------------------------------------
     97# CALL INSTALLATION SCRIPT WITH DETERMINED COMMANDLINE ARGUMENTS
     98
     99$pyscript $parameterlist
     100
  • source/python/classes/ControlFile.py

    r4971f63 rca867de  
    5757import inspect
    5858
     59# software specific classes and modules from flex_extract
     60sys.path.append('../')
    5961import _config
     62from mods.tools import my_error
    6063
    6164# ------------------------------------------------------------------------------
     
    144147        self.ectrans = 0
    145148        self.inputdir = _config.PATH_INPUT_DIR
    146         self.outputdir = self.inputdir
     149        self.outputdir = None
    147150        self.ecmwfdatadir = _config.PATH_FLEXEXTRACT_DIR
    148151        self.exedir = _config.PATH_FORTRAN_SRC
     
    177180            <nothing>
    178181        '''
    179         from mods.tools import my_error
    180 
    181         # read whole CONTROL file
    182182
    183183        try:
     
    355355            self.end_date = self.start_date
    356356
     357        # basetime has only two possible values
     358        if self.basetime:
     359            if int(self.basetime) != 0 and int(self.basetime) != 12:
     360                print('Basetime has an invalid value!')
     361                print('Basetime = ' + str(self.basetime))
     362                sys.exit(1)
     363
    357364        # assure consistency of levelist and level
    358365        if self.levelist is None and self.level is None:
     
    407414            self.flexpart_root_scripts = self.ecmwfdatadir
    408415
     416        if not self.outputdir:
     417            self.outputdir = self.inputdir
     418
    409419        if not isinstance(self.mailfail, list):
    410420            if ',' in self.mailfail:
  • source/python/classes/EcFlexpart.py

    r27fe969 rca867de  
    7474# MODULES
    7575# ------------------------------------------------------------------------------
     76import os
     77import sys
     78import glob
     79import shutil
    7680import subprocess
    77 import shutil
    78 import os
    79 import glob
    8081from datetime import datetime, timedelta
    8182import numpy as np
     
    8586
    8687# software specific classes and modules from flex_extract
     88sys.path.append('../')
    8789import _config
    8890from GribTools import GribTools
    8991from mods.tools import init128, to_param_id, silent_remove, product, my_error
    9092from MarsRetrieval import MarsRetrieval
    91 import mods.disaggregation
     93import mods.disaggregation as disaggregation
    9294
    9395# ------------------------------------------------------------------------------
     
    287289                                        'SFC', '1', self.grid]
    288290
    289         # if needed, add additional WRF specific parameters here
    290 
    291291        return
    292292
     
    382382
    383383
     384    def _mk_index_values(self, inputdir, inputfiles, keys):
     385        '''
     386        @Description:
     387            Creates an index file for a set of grib parameter keys.
     388            The values from the index keys are returned in a list.
     389
     390        @Input:
     391            keys: dictionary
     392                List of parameter names which serves as index.
     393
     394            inputfiles: instance of UioFiles
     395                Contains a list of files.
     396
     397        @Return:
     398            iid: grib_index
     399                This is a grib specific index structure to access
     400                messages in a file.
     401
     402            index_vals: list
     403                Contains the values from the keys used for a distinct selection
     404                of grib messages in processing  the grib files.
     405                Content looks like e.g.:
     406                index_vals[0]: ('20171106', '20171107', '20171108') ; date
     407                index_vals[1]: ('0', '1200', '1800', '600') ; time
     408                index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
     409        '''
     410        iid = None
     411        index_keys = keys
     412
     413        indexfile = os.path.join(inputdir, _config.FILE_GRIB_INDEX)
     414        silent_remove(indexfile)
     415        grib = GribTools(inputfiles.files)
     416        # creates new index file
     417        iid = grib.index(index_keys=index_keys, index_file=indexfile)
     418
     419        # read the values of index keys
     420        index_vals = []
     421        for key in index_keys:
     422            #index_vals.append(grib_index_get(iid, key))
     423            #print(index_vals[-1])
     424            key_vals = grib_index_get(iid, key)
     425            print(key_vals)
     426            # have to sort the steps for disaggregation,
     427            # therefore convert to int first
     428            if key == 'step':
     429                key_vals = [int(k) for k in key_vals]
     430                key_vals.sort()
     431                key_vals = [str(k) for k in key_vals]
     432            index_vals.append(key_vals)
     433            # index_vals looks for example like:
     434            # index_vals[0]: ('20171106', '20171107', '20171108') ; date
     435            # index_vals[1]: ('0', '1200') ; time
     436            # index_vals[2]: (3', '6', '9', '12') ; stepRange
     437
     438        return iid, index_vals
     439
     440
    384441    def retrieve(self, server, dates, request, inputdir='.'):
    385442        '''
     
    427484        t24h = timedelta(hours=24)
    428485
    429         # dictionary which contains all parameter for the mars request
     486        # dictionary which contains all parameter for the mars request,
    430487        # entries with a "None" will change in different requests and will
    431488        # therefore be set in each request seperately
     
    674731        table128 = init128(_config.PATH_GRIBTABLE)
    675732        pars = to_param_id(self.params['OG_acc_SL'][0], table128)
     733
     734        iid = None
     735        index_vals = None
     736
     737        # get the values of the keys which are used for distinct access
     738        # of grib messages via product
    676739        index_keys = ["date", "time", "step"]
    677         indexfile = os.path.join(c.inputdir, _config.FILE_GRIB_INDEX)
    678         silent_remove(indexfile)
    679         grib = GribTools(inputfiles.files)
    680         # creates new index file
    681         iid = grib.index(index_keys=index_keys, index_file=indexfile)
    682 
    683         # read values of index keys
    684         index_vals = []
    685         for key in index_keys:
    686             key_vals = grib_index_get(iid, key)
    687             print(key_vals)
    688             # have to sort the steps for disaggregation,
    689             # therefore convert to int first
    690             if key == 'step':
    691                 key_vals = [int(k) for k in key_vals]
    692                 key_vals.sort()
    693                 key_vals = [str(k) for k in key_vals]
    694             index_vals.append(key_vals)
    695             # index_vals looks for example like:
    696             # index_vals[0]: ('20171106', '20171107', '20171108') ; date
    697             # index_vals[1]: ('0', '1200') ; time
    698             # index_vals[2]: (3', '6', '9', '12') ; stepRange
     740        iid, index_vals = self._mk_index_values(c.inputdir,
     741                                                inputfiles,
     742                                                index_keys)
     743        # index_vals looks like e.g.:
     744        # index_vals[0]: ('20171106', '20171107', '20171108') ; date
     745        # index_vals[1]: ('0', '1200', '1800', '600') ; time
     746        # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
    699747
    700748        valsdict = {}
    701749        svalsdict = {}
    702         stepsdict = {}
     750#        stepsdict = {}
    703751        for p in pars:
    704752            valsdict[str(p)] = []
    705753            svalsdict[str(p)] = []
    706             stepsdict[str(p)] = []
     754#            stepsdict[str(p)] = []
    707755
    708756        print('maxstep: ', c.maxstep)
    709757
     758        # "product" genereates each possible combination between the
     759        # values of the index keys
    710760        for prod in product(*index_vals):
    711761            # e.g. prod = ('20170505', '0', '12')
    712762            #             (  date    ,time, step)
    713             # per date e.g. time = 0, 1200
    714             # per time e.g. step = 3, 6, 9, 12
    715             print('current prod: ', prod)
     763
     764            print('current product: ', prod)
     765
    716766            for i in range(len(index_keys)):
    717767                grib_index_select(iid, index_keys[i], prod[i])
     
    720770            gid = grib_new_from_index(iid)
    721771
    722             # if there is data for this product combination
    723             # prepare some date and time parameter before reading the data
    724             if gid is not None:
    725                 cdate = grib_get(gid, 'date')
    726                 time = grib_get(gid, 'time')
    727                 step = grib_get(gid, 'step')
    728                 # date+time+step-2*dtime
    729                 # (since interpolated value valid for step-2*dtime)
    730                 sdate = datetime(year=cdate/10000,
    731                                  month=(cdate % 10000)/100,
    732                                  day=(cdate % 100),
    733                                  hour=time/100)
    734                 fdate = sdate + timedelta(hours=step-2*int(c.dtime))
    735                 sdates = sdate + timedelta(hours=step)
    736                 elimit = None
    737             else:
    738                 break
     772            # if there is no data for this specific time combination / product
     773            # skip the rest of the for loop and start with next timestep/product
     774            if not gid:
     775                continue
     776
     777            # create correct timestamp from the three time informations
     778            cdate = str(grib_get(gid, 'date'))
     779            ctime = '{:0>2}'.format(grib_get(gid, 'time')/100)
     780            cstep = '{:0>3}'.format(grib_get(gid, 'step'))
     781            t_date = datetime.strptime(cdate + ctime, '%Y%m%d%H')
     782            t_dt = t_date + timedelta(hours=int(cstep))
     783            t_m1dt = t_date + timedelta(hours=int(cstep)-int(c.dtime))
     784            t_m2dt = t_date + timedelta(hours=int(cstep)-2*int(c.dtime))
     785            t_enddate = None
    739786
    740787            if c.maxstep > 12:
    741788                fnout = os.path.join(c.inputdir, 'flux' +
    742                                      sdate.strftime('%Y%m%d') +
    743                                      '.{:0>2}'.format(time/100) +
     789                                     t_date.strftime('%Y%m%d.%H') +
    744790                                     '.{:0>3}'.format(step-2*int(c.dtime)))
    745791                gnout = os.path.join(c.inputdir, 'flux' +
    746                                      sdate.strftime('%Y%m%d') +
    747                                      '.{:0>2}'.format(time/100) +
     792                                     t_date.strftime('%Y%m%d.%H') +
    748793                                     '.{:0>3}'.format(step-int(c.dtime)))
    749794                hnout = os.path.join(c.inputdir, 'flux' +
    750                                      sdate.strftime('%Y%m%d') +
    751                                      '.{:0>2}'.format(time/100) +
     795                                     t_date.strftime('%Y%m%d.%H') +
    752796                                     '.{:0>3}'.format(step))
    753797            else:
    754798                fnout = os.path.join(c.inputdir, 'flux' +
    755                                      fdate.strftime('%Y%m%d%H'))
     799                                     t_m2dt.strftime('%Y%m%d%H'))
    756800                gnout = os.path.join(c.inputdir, 'flux' +
    757                                      (fdate + timedelta(hours=int(c.dtime))).
    758                                      strftime('%Y%m%d%H'))
     801                                     t_m1dt.strftime('%Y%m%d%H'))
    759802                hnout = os.path.join(c.inputdir, 'flux' +
    760                                      sdates.strftime('%Y%m%d%H'))
     803                                     t_dt.strftime('%Y%m%d%H'))
    761804
    762805            print("outputfile = " + fnout)
    763             f_handle = open(fnout, 'w')
    764             g_handle = open(gnout, 'w')
    765             h_handle = open(hnout, 'w')
    766806
    767807            # read message for message and store relevant data fields
    768808            # data keywords are stored in pars
    769809            while 1:
    770                 if gid is None:
     810                if not gid:
    771811                    break
    772812                cparamId = str(grib_get(gid, 'paramId'))
    773813                step = grib_get(gid, 'step')
    774                 atime = grib_get(gid, 'time')
     814                time = grib_get(gid, 'time')
    775815                ni = grib_get(gid, 'Ni')
    776816                nj = grib_get(gid, 'Nj')
     
    779819                    vdp = valsdict[cparamId]
    780820                    svdp = svalsdict[cparamId]
    781                     sd = stepsdict[cparamId]
     821 #                   sd = stepsdict[cparamId]
    782822
    783823                    if cparamId == '142' or cparamId == '143':
     
    793833                        svdp.append((vdp[-1] - vdp[-2]) / int(c.dtime))
    794834
    795                     print(cparamId, atime, step, len(values),
     835                    print(cparamId, time, step, len(values),
    796836                          values[0], np.std(values))
    797                     # save the 1/3-hourly or specific values
    798                     # svdp.append(values[:])
    799                     sd.append(step)
     837
    800838                    # len(svdp) correspond to the time
    801839                    if len(svdp) >= 3:
     
    807845
    808846                            if not (step == c.maxstep and c.maxstep > 12 \
    809                                     or sdates == elimit):
     847                                    or t_dt == t_enddate):
    810848                                vdp.pop(0)
    811849                                svdp.pop(0)
     
    817855
    818856                        grib_set_values(gid, values)
     857
    819858                        if c.maxstep > 12:
    820859                            grib_set(gid, 'step', max(0, step-2*int(c.dtime)))
    821860                        else:
    822861                            grib_set(gid, 'step', 0)
    823                             grib_set(gid, 'time', fdate.hour*100)
    824                             grib_set(gid, 'date', fdate.year*10000 +
    825                                      fdate.month*100+fdate.day)
    826                         grib_write(gid, f_handle)
     862                            grib_set(gid, 'time', t_m2dt.hour*100)
     863                            grib_set(gid, 'date', int(t_m2dt.strftime('%Y%m%d')))
     864
     865                        with open(fnout, 'w') as f_handle:
     866                            grib_write(gid, f_handle)
    827867
    828868                        if c.basetime:
    829                             elimit = datetime.strptime(c.end_date +
    830                                                        c.basetime, '%Y%m%d%H')
     869                            t_enddate = datetime.strptime(c.end_date +
     870                                                          c.basetime,
     871                                                          '%Y%m%d%H')
    831872                        else:
    832                             elimit = sdate + timedelta(2*int(c.dtime))
     873                            t_enddate = t_date + timedelta(2*int(c.dtime))
    833874
    834875                        # squeeze out information of last two steps contained
    835876                        # in svdp
    836877                        # if step+int(c.dtime) == c.maxstep and c.maxstep>12
    837                         # or sdates+timedelta(hours = int(c.dtime))
    838                         # >= elimit:
     878                        # or t_dt+timedelta(hours = int(c.dtime))
     879                        # >= t_enddate:
    839880                        # Note that svdp[0] has not been popped in this case
    840881
    841882                        if step == c.maxstep and c.maxstep > 12 or \
    842                            sdates == elimit:
     883                           t_dt == t_enddate:
    843884
    844885                            values = svdp[3]
    845886                            grib_set_values(gid, values)
    846887                            grib_set(gid, 'step', 0)
    847                             truedatetime = fdate + timedelta(hours=
     888                            truedatetime = t_m2dt + timedelta(hours=
    848889                                                             2*int(c.dtime))
    849890                            grib_set(gid, 'time', truedatetime.hour * 100)
     
    851892                                     truedatetime.month * 100 +
    852893                                     truedatetime.day)
    853                             grib_write(gid, h_handle)
     894                            with open(hnout, 'w') as h_handle:
     895                                grib_write(gid, h_handle)
    854896
    855897                            #values = (svdp[1]+svdp[2])/2.
     
    860902
    861903                            grib_set(gid, 'step', 0)
    862                             truedatetime = fdate + timedelta(hours=int(c.dtime))
     904                            truedatetime = t_m2dt + timedelta(hours=int(c.dtime))
    863905                            grib_set(gid, 'time', truedatetime.hour * 100)
    864906                            grib_set(gid, 'date', truedatetime.year * 10000 +
     
    866908                                     truedatetime.day)
    867909                            grib_set_values(gid, values)
    868                             grib_write(gid, g_handle)
    869 
    870                     grib_release(gid)
    871 
    872                     gid = grib_new_from_index(iid)
    873 
    874             f_handle.close()
    875             g_handle.close()
    876             h_handle.close()
     910                            with open(gnout, 'w') as g_handle:
     911                                grib_write(gid, g_handle)
     912
     913                grib_release(gid)
     914
     915                gid = grib_new_from_index(iid)
    877916
    878917        grib_index_release(iid)
     
    913952        '''
    914953
    915         table128 = init128(_config.PATH_GRIBTABLE)
    916         wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\
    917                             stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
    918                             table128)
    919 
     954        if c.wrf:
     955            table128 = init128(_config.PATH_GRIBTABLE)
     956            wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\
     957                                   stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
     958                                  table128)
     959
     960        # these numbers are indices for the temporary files "fort.xx"
     961        # which are used to seperate the grib fields to,
     962        # for the Fortran program input
     963        # 10: U,V | 11: T | 12: lnsp | 13: D | 16: sfc fields
     964        # 17: Q | 18: Q , gaussian| 19: w | 21: etadot | 22: clwc+ciwc
     965        fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
     966                 '17':None, '18':None, '19':None, '21':None, '22':None}
     967
     968        iid = None
     969        index_vals = None
     970
     971        # get the values of the keys which are used for distinct access
     972        # of grib messages via product
    920973        index_keys = ["date", "time", "step"]
    921         indexfile = os.path.join(c.inputdir, _config.FILE_GRIB_INDEX)
    922         silent_remove(indexfile)
    923         grib = GribTools(inputfiles.files)
    924         # creates new index file
    925         iid = grib.index(index_keys=index_keys, index_file=indexfile)
    926 
    927         # read values of index keys
    928         index_vals = []
    929         for key in index_keys:
    930             index_vals.append(grib_index_get(iid, key))
    931             print(index_vals[-1])
    932             # index_vals looks for example like:
    933             # index_vals[0]: ('20171106', '20171107', '20171108') ; date
    934             # index_vals[1]: ('0', '1200', '1800', '600') ; time
    935             # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
    936 
    937         fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
    938                  '17':None, '19':None, '21':None, '22':None, '20':None}
    939 
     974        iid, index_vals = self._mk_index_values(c.inputdir,
     975                                                inputfiles,
     976                                                index_keys)
     977        # index_vals looks like e.g.:
     978        # index_vals[0]: ('20171106', '20171107', '20171108') ; date
     979        # index_vals[1]: ('0', '1200', '1800', '600') ; time
     980        # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
     981
     982        # "product" genereates each possible combination between the
     983        # values of the index keys
    940984        for prod in product(*index_vals):
    941985            # e.g. prod = ('20170505', '0', '12')
    942986            #             (  date    ,time, step)
    943             # per date e.g. time = 0, 1200
    944             # per time e.g. step = 3, 6, 9, 12
    945             # flag for Fortran program and file merging
    946             convertFlag = False
    947             print('current prod: ', prod)
    948             # e.g. prod = ('20170505', '0', '12')
    949             #             (  date    ,time, step)
    950             # per date e.g. time = 0, 600, 1200, 1800
    951             # per time e.g. step = 0, 3, 6, 9, 12
     987
     988            print('current product: ', prod)
     989
    952990            for i in range(len(index_keys)):
    953991                grib_index_select(iid, index_keys[i], prod[i])
     
    956994            gid = grib_new_from_index(iid)
    957995
    958             # if there is data for this product combination
    959             # prepare some date and time parameter before reading the data
    960             if gid is not None:
    961                 # Combine all temporary data files into final grib file if
    962                 # gid is at least one time not None. Therefore set convertFlag
    963                 # to save information. The Fortran program is also
    964                 # only executed if convertFlag is True
    965                 convertFlag = True
    966                 # remove old fort.* files and open new ones
    967                 # they are just valid for a single product
    968                 for k, f in fdict.iteritems():
    969                     fortfile = os.path.join(c.inputdir, 'fort.' + k)
    970                     silent_remove(fortfile)
    971                     fdict[k] = open(fortfile, 'w')
    972 
    973                 cdate = str(grib_get(gid, 'date'))
    974                 time = grib_get(gid, 'time')
    975                 step = grib_get(gid, 'step')
    976                 # create correct timestamp from the three time informations
    977                 # date, time, step
    978                 timestamp = datetime.strptime(cdate + '{:0>2}'.format(time/100),
    979                                               '%Y%m%d%H')
    980                 timestamp += timedelta(hours=int(step))
    981                 cdateH = datetime.strftime(timestamp, '%Y%m%d%H')
    982 
    983                 if c.basetime:
    984                     slimit = datetime.strptime(c.start_date + '00', '%Y%m%d%H')
    985                     bt = '23'
    986                     if c.basetime == '00':
    987                         bt = '00'
    988                         slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\
    989                             - timedelta(hours=12-int(c.dtime))
    990                     if c.basetime == '12':
    991                         bt = '12'
    992                         slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\
    993                             - timedelta(hours=12-int(c.dtime))
    994 
    995                     elimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')
    996 
    997                     if timestamp < slimit or timestamp > elimit:
    998                         continue
    999             else:
    1000                 pass
    1001 
    1002             try:
    1003                 if c.wrf:
    1004                     if 'olddate' not in locals() or cdate != olddate:
    1005                         fwrf = open(os.path.join(c.outputdir,
    1006                                     'WRF' + cdate + '.{:0>2}'.format(time) +
    1007                                     '.000.grb2'), 'w')
    1008                         olddate = cdate[:]
    1009             except AttributeError:
    1010                 pass
    1011 
    1012             # helper variable to remember which fields were already used.
     996            # if there is no data for this specific time combination / product
     997            # skip the rest of the for loop and start with next timestep/product
     998            if not gid:
     999                continue
     1000
     1001            # remove old fort.* files and open new ones
     1002            # they are just valid for a single product
     1003            for k, f in fdict.iteritems():
     1004                fortfile = os.path.join(c.inputdir, 'fort.' + k)
     1005                silent_remove(fortfile)
     1006                fdict[k] = open(fortfile, 'w')
     1007
     1008            # create correct timestamp from the three time informations
     1009            cdate = str(grib_get(gid, 'date'))
     1010            ctime = '{:0>2}'.format(grib_get(gid, 'time')/100)
     1011            cstep = '{:0>3}'.format(grib_get(gid, 'step'))
     1012            timestamp = datetime.strptime(cdate + ctime, '%Y%m%d%H')
     1013            timestamp += timedelta(hours=int(cstep))
     1014            cdate_hour = datetime.strftime(timestamp, '%Y%m%d%H')
     1015
     1016            # if the timestamp is out of basetime start/end date period,
     1017            # skip this specific product
     1018            if c.basetime:
     1019                start_time = datetime.strptime(c.end_date + c.basetime,
     1020                                                '%Y%m%d%H') - time_delta
     1021                end_time = datetime.strptime(c.end_date + c.basetime,
     1022                                             '%Y%m%d%H')
     1023                if timestamp < start_time or timestamp > end_time:
     1024                    continue
     1025
     1026            if c.wrf:
     1027                if 'olddate' not in locals() or cdate != olddate:
     1028                    fwrf = open(os.path.join(c.outputdir,
     1029                                'WRF' + cdate + '.' + ctime + '.000.grb2'), 'w')
     1030                    olddate = cdate[:]
     1031
     1032            # savedfields remembers which fields were already used.
    10131033            savedfields = []
     1034            # sum of cloud liquid and ice water content
     1035            scwc = None
    10141036            while 1:
    1015                 if gid is None:
     1037                if not gid:
    10161038                    break
    10171039                paramId = grib_get(gid, 'paramId')
    10181040                gridtype = grib_get(gid, 'gridType')
    10191041                levtype = grib_get(gid, 'typeOfLevel')
    1020                 if paramId == 133 and gridtype == 'reduced_gg':
    1021                 # Specific humidity (Q.grb) is used as a template only
    1022                 # so we need the first we "meet"
    1023                     with open(os.path.join(c.inputdir, 'fort.18'), 'w') as fout:
    1024                         grib_write(gid, fout)
    1025                 elif paramId == 131 or paramId == 132:
     1042                if paramId == 77: # ETADOT
     1043                    grib_write(gid, fdict['21'])
     1044                elif paramId == 130: # T
     1045                    grib_write(gid, fdict['11'])
     1046                elif paramId == 131 or paramId == 132: # U, V wind component
    10261047                    grib_write(gid, fdict['10'])
    1027                 elif paramId == 130:
    1028                     grib_write(gid, fdict['11'])
    1029                 elif paramId == 133 and gridtype != 'reduced_gg':
     1048                elif paramId == 133 and gridtype != 'reduced_gg': # Q
    10301049                    grib_write(gid, fdict['17'])
    1031                 elif paramId == 152:
     1050                elif paramId == 133 and gridtype == 'reduced_gg': # Q, gaussian
     1051                    grib_write(gid, fdict['18'])
     1052                elif paramId == 135: # W
     1053                    grib_write(gid, fdict['19'])
     1054                elif paramId == 152: # LNSP
    10321055                    grib_write(gid, fdict['12'])
    1033                 elif paramId == 155 and gridtype == 'sh':
     1056                elif paramId == 155 and gridtype == 'sh': # D
    10341057                    grib_write(gid, fdict['13'])
    1035                 elif  paramId in [129, 138, 155] and levtype == 'hybrid' \
    1036                         and c.wrf:
    1037                     pass
    1038                 elif paramId == 246 or paramId == 247:
    1039                     # cloud liquid water and ice
    1040                     if paramId == 246:
    1041                         clwc = grib_get_values(gid)
     1058                elif paramId == 246 or paramId == 247: # CLWC, CIWC
     1059                    # sum cloud liquid water and ice
     1060                    if not scwc:
     1061                        scwc = grib_get_values(gid)
    10421062                    else:
    1043                         clwc += grib_get_values(gid)
    1044                         grib_set_values(gid, clwc)
     1063                        scwc += grib_get_values(gid)
     1064                        grib_set_values(gid, scwc)
    10451065                        grib_set(gid, 'paramId', 201031)
    10461066                        grib_write(gid, fdict['22'])
    1047                 elif paramId == 135:
    1048                     grib_write(gid, fdict['19'])
    1049                 elif paramId == 77:
    1050                     grib_write(gid, fdict['21'])
     1067                elif c.wrf and paramId in [129, 138, 155] and \
     1068                      levtype == 'hybrid': # Z, VO, D
     1069                    # do not do anything right now
     1070                    # these are specific parameter for WRF
     1071                    pass
    10511072                else:
    10521073                    if paramId not in savedfields:
     1074                        # SD/MSL/TCC/10U/10V/2T/2D/Z/LSM/SDOR/CVL/CVH/SR
     1075                        # and all ADDPAR parameter
    10531076                        grib_write(gid, fdict['16'])
    10541077                        savedfields.append(paramId)
     
    10741097                f.close()
    10751098
    1076             # call for Fortran program if flag is True
    1077             if convertFlag:
    1078                 pwd = os.getcwd()
    1079                 os.chdir(c.inputdir)
    1080                 if os.stat('fort.21').st_size == 0 and c.eta:
    1081                     print('Parameter 77 (etadot) is missing, most likely it is \
    1082                            not available for this type or date/time\n')
    1083                     print('Check parameters CLASS, TYPE, STREAM, START_DATE\n')
    1084                     my_error(c.mailfail, 'fort.21 is empty while parameter eta \
    1085                              is set to 1 in CONTROL file')
    1086 
    1087                 # create the corresponding output file fort.15
    1088                 # (generated by Fortran program) + fort.16 (paramId 167 and 168)
    1089                 p = subprocess.check_call([os.path.join(
    1090                     c.exedir, _config.FORTRAN_EXECUTABLE)], shell=True)
    1091                 os.chdir(pwd)
    1092 
    1093                 # create final output filename, e.g. EN13040500 (ENYYMMDDHH)
    1094                 fnout = os.path.join(c.inputdir, c.prefix)
    1095                 if c.maxstep > 12:
    1096                     suffix = cdate[2:8] + '.{:0>2}'.format(time/100) + \
    1097                              '.{:0>3}'.format(step)
    1098                 else:
    1099                     suffix = cdateH[2:10]
    1100                 fnout += suffix
    1101                 print("outputfile = " + fnout)
    1102                 self.outputfilelist.append(fnout) # needed for final processing
    1103 
    1104                 # create outputfile and copy all data from intermediate files
    1105                 # to the outputfile (final GRIB files)
    1106                 orolsm = os.path.basename(glob.glob(
    1107                     c.inputdir + '/OG_OROLSM__SL.*.' + c.ppid + '*')[0])
    1108                 fluxfile = 'flux' + cdate[0:2] + suffix
    1109                 if not c.cwc:
    1110                     flist = ['fort.15', fluxfile, 'fort.16', orolsm]
    1111                 else:
    1112                     flist = ['fort.15', 'fort.22', fluxfile, 'fort.16', orolsm]
    1113 
    1114                 with open(fnout, 'wb') as fout:
    1115                     for f in flist:
    1116                         shutil.copyfileobj(open(os.path.join(c.inputdir, f),
    1117                                                 'rb'), fout)
    1118 
    1119                 if c.omega:
    1120                     with open(os.path.join(c.outputdir, 'OMEGA'), 'wb') as fout:
    1121                         shutil.copyfileobj(
    1122                             open(os.path.join(c.inputdir, 'fort.25'),
    1123                                  'rb'), fout)
     1099            # call for Fortran program to convert e.g. reduced_gg grids to
     1100            # regular_ll and calculate detadot/dp
     1101            pwd = os.getcwd()
     1102            os.chdir(c.inputdir)
     1103            if os.stat('fort.21').st_size == 0 and c.eta:
     1104                print('Parameter 77 (etadot) is missing, most likely it is \
     1105                       not available for this type or date/time\n')
     1106                print('Check parameters CLASS, TYPE, STREAM, START_DATE\n')
     1107                my_error(c.mailfail, 'fort.21 is empty while parameter eta \
     1108                         is set to 1 in CONTROL file')
     1109
     1110            # Fortran program creates file fort.15 (with u,v,etadot,t,sp,q)
     1111            p = subprocess.check_call([os.path.join(
     1112                c.exedir, _config.FORTRAN_EXECUTABLE)], shell=True)
     1113            os.chdir(pwd)
     1114
     1115            # create name of final output file, e.g. EN13040500 (ENYYMMDDHH)
     1116            if c.maxstep > 12:
     1117                suffix = cdate[2:8] + '.' + ctime + '.' + cstep
    11241118            else:
    1125                 pass
     1119                suffix = cdate_hour[2:10]
     1120            fnout = os.path.join(c.inputdir, c.prefix + suffix)
     1121            print("outputfile = " + fnout)
     1122            # collect for final processing
     1123            self.outputfilelist.append(os.path.basename(fnout))
     1124
     1125            # create outputfile and copy all data from intermediate files
     1126            # to the outputfile (final GRIB input files for FLEXPART)
     1127            orolsm = os.path.basename(glob.glob(c.inputdir +
     1128                                        '/OG_OROLSM__SL.*.' + c.ppid + '*')[0])
     1129            fluxfile = 'flux' + cdate[0:2] + suffix
     1130            if not c.cwc:
     1131                flist = ['fort.15', fluxfile, 'fort.16', orolsm]
     1132            else:
     1133                flist = ['fort.15', 'fort.22', fluxfile, 'fort.16', orolsm]
     1134
     1135            with open(fnout, 'wb') as fout:
     1136                for f in flist:
     1137                    shutil.copyfileobj(open(os.path.join(c.inputdir, f), 'rb'),
     1138                                       fout)
     1139
     1140            if c.omega:
     1141                with open(os.path.join(c.outputdir, 'OMEGA'), 'wb') as fout:
     1142                    shutil.copyfileobj(open(os.path.join(c.inputdir, 'fort.25'),
     1143                                            'rb'), fout)
    11261144
    11271145        if c.wrf:
     
    11681186                  .format(c.ectrans, c.gateway, c.destination))
    11691187
    1170         print('Output filelist: \n')
     1188        print('Output filelist: ')
    11711189        print(self.outputfilelist)
    11721190
     
    11911209        if c.outputdir != c.inputdir:
    11921210            for ofile in self.outputfilelist:
    1193                 p = subprocess.check_call(['mv', ofile, c.outputdir])
     1211                p = subprocess.check_call(['mv',
     1212                                           os.path.join(c.inputdir, ofile),
     1213                                           c.outputdir])
    11941214
    11951215        return
  • source/python/classes/MarsRetrieval.py

    r295ff45 rca867de  
    6262# MODULES
    6363# ------------------------------------------------------------------------------
     64import os
     65import sys
    6466import subprocess
    65 import os
    66 
     67
     68# software specific classes and modules from flex_extract
     69sys.path.append('../')
    6770import _config
    6871# ------------------------------------------------------------------------------
  • source/python/classes/UioFiles.py

    r25b14be rca867de  
    4848# ------------------------------------------------------------------------------
    4949import os
     50import sys
    5051import fnmatch
    5152
    5253# software specific module from flex_extract
     54sys.path.append('../')
    5355#import profiling
    5456from mods.tools import silent_remove, get_list_as_string
  • source/python/mods/get_mars_data.py

    r27fe969 rca867de  
    4848import os
    4949import sys
    50 import datetime
    5150import inspect
     51from datetime import datetime, timedelta
     52
     53# software specific classes and modules from flex_extract
     54sys.path.append('../')
     55import _config
     56from tools import my_error, normal_exit, get_cmdline_arguments, read_ecenv
     57from classes.EcFlexpart import EcFlexpart
     58from classes.UioFiles import UioFiles
     59
    5260try:
    5361    ecapi = True
     
    5563except ImportError:
    5664    ecapi = False
    57 
    58 # software specific classes and modules from flex_extract
    59 import _config
    60 from tools import my_error, normal_exit, get_cmdline_arguments, read_ecenv
    61 from classes.EcFlexpart import EcFlexpart
    62 from classes.UioFiles import UioFiles
    6365# ------------------------------------------------------------------------------
    6466# FUNCTION
     
    147149    # allerdings ist das relevant und ersichtlich an den NICHT FLUSS DATEN
    148150
    149 
    150     # set start date of retrieval period
    151     start = datetime.date(year=int(c.start_date[:4]),
    152                           month=int(c.start_date[4:6]),
    153                           day=int(c.start_date[6:]))
    154     startm1 = start - datetime.timedelta(days=1)
    155 
    156     # set end date of retrieval period
    157     end = datetime.date(year=int(c.end_date[:4]),
    158                         month=int(c.end_date[4:6]),
    159                         day=int(c.end_date[6:]))
    160 
    161     # set time period for one single retrieval
    162     datechunk = datetime.timedelta(days=int(c.date_chunk))
     151    start = datetime.strptime(c.start_date, '%Y%m%d')
     152    end = datetime.strptime(c.end_date, '%Y%m%d')
     153    # time period for one single retrieval
     154    datechunk = timedelta(days=int(c.date_chunk))
    163155
    164156    if c.basetime == '00':
    165         start = startm1
     157        start = start - timedelta(days=1)
     158
     159    if c.maxstep <= 24:
     160        startm1 = start - timedelta(days=1)
    166161
    167162    if c.basetime == '00' or c.basetime == '12':
    168         # endp1 = end + datetime.timedelta(days=1)
     163        # endp1 = end + timedelta(days=1)
    169164        endp1 = end
    170165    else:
    171         # endp1 = end + datetime.timedelta(days=2)
    172         endp1 = end + datetime.timedelta(days=1)
     166        # endp1 = end + timedelta(days=2)
     167        endp1 = end + timedelta(days=1)
    173168
    174169    # --------------  flux data ------------------------------------------------
     
    243238    # we only need to add datechunk - 1 days to retrieval
    244239    # for a period
    245     delta_t_m1 = delta_t - datetime.timedelta(days=1)
     240    delta_t_m1 = delta_t - timedelta(days=1)
    246241
    247242    day = start
  • source/python/mods/prepare_flexpart.py

    r27fe969 rca867de  
    5959
    6060# software specific classes and modules from flex_extract
     61sys.path.append('../')
    6162import _config
    6263from classes.UioFiles import UioFiles
     64from classes.ControlFile import ControlFile
    6365from tools import clean_up, get_cmdline_arguments, read_ecenv
    6466from classes.EcFlexpart import EcFlexpart
  • source/python/mods/tools.py

    r27fe969 rca867de  
    6060# ------------------------------------------------------------------------------
    6161
     62def none_or_str(value):
     63    '''
     64    @Description:
     65        Converts the input string into pythons None-type if the string
     66        contains "None".
     67
     68    @Input:
     69        value: string
     70            String to be checked for the "None" word.
     71
     72    @Return:
     73        None or value:
     74            Return depends on the content of the input value. If it was "None",
     75            then the python type None is returned. Otherwise the string itself.
     76    '''
     77    if value == 'None':
     78        return None
     79    return value
     80
     81def none_or_int(value):
     82    '''
     83    @Description:
     84        Converts the input string into pythons None-type if the string
     85        contains "None". Otherwise it is converted to an integer value.
     86
     87    @Input:
     88        value: string
     89            String to be checked for the "None" word.
     90
     91    @Return:
     92        None or int(value):
     93            Return depends on the content of the input value. If it was "None",
     94            then the python type None is returned. Otherwise the string is
     95            converted into an integer value.
     96    '''
     97    if value == 'None':
     98        return None
     99    return int(value)
     100
    62101def get_cmdline_arguments():
    63102    '''
     
    79118
    80119    # the most important arguments
    81     parser.add_argument("--start_date", dest="start_date", default=None,
     120    parser.add_argument("--start_date", dest="start_date",
     121                        type=none_or_str, default=None,
    82122                        help="start date YYYYMMDD")
    83     parser.add_argument("--end_date", dest="end_date", default=None,
     123    parser.add_argument("--end_date", dest="end_date",
     124                        type=none_or_str, default=None,
    84125                        help="end_date YYYYMMDD")
    85     parser.add_argument("--date_chunk", dest="date_chunk", default=None,
     126    parser.add_argument("--date_chunk", dest="date_chunk",
     127                        type=none_or_int, default=None,
    86128                        help="# of days to be retrieved at once")
     129    parser.add_argument("--controlfile", dest="controlfile",
     130                        type=none_or_str, default='CONTROL.temp',
     131                        help="file with CONTROL parameters")
     132
     133    # parameter for extra output information
     134    parser.add_argument("--debug", dest="debug",
     135                        type=none_or_int, default=None,
     136                        help="debug mode - leave temporary files intact")
     137    parser.add_argument("--request", dest="request",
     138                        type=none_or_int, default=None,
     139                        help="list all mars request in file mars_requests.dat \
     140                        and skip submission to mars")
    87141
    88142    # some arguments that override the default in the CONTROL file
    89     parser.add_argument("--basetime", dest="basetime", default=None,
    90                         help="base such as 00/12 (for half day retrievals)")
    91     parser.add_argument("--step", dest="step", default=None,
     143    parser.add_argument("--basetime", dest="basetime",
     144                        type=none_or_int, default=None,
     145                        help="base such as 00 or 12 (for half day retrievals)")
     146    parser.add_argument("--step", dest="step",
     147                        type=none_or_str, default=None,
    92148                        help="steps such as 00/to/48")
    93     parser.add_argument("--levelist", dest="levelist", default=None,
     149    parser.add_argument("--levelist", dest="levelist",
     150                        type=none_or_str, default=None,
    94151                        help="Vertical levels to be retrieved, e.g. 30/to/60")
    95     parser.add_argument("--area", dest="area", default=None,
     152    parser.add_argument("--area", dest="area",
     153                        type=none_or_str, default=None,
    96154                        help="area defined as north/west/south/east")
    97155
    98156    # set the working directories
    99     parser.add_argument("--inputdir", dest="inputdir", default=None,
     157    parser.add_argument("--inputdir", dest="inputdir",
     158                        type=none_or_str, default=None,
    100159                        help="root directory for storing intermediate files")
    101     parser.add_argument("--outputdir", dest="outputdir", default=None,
     160    parser.add_argument("--outputdir", dest="outputdir",
     161                        type=none_or_str, default=None,
    102162                        help="root directory for storing output files")
    103163    parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts",
    104                         default=None,
     164                        type=none_or_str, default=None,
    105165                        help="FLEXPART root directory (to find grib2flexpart \
    106166                        and COMMAND file)\n Normally flex_extract resides in \
     
    108168
    109169    # this is only used by prepare_flexpart.py to rerun a postprocessing step
    110     parser.add_argument("--ppid", dest="ppid", default=None,
     170    parser.add_argument("--ppid", dest="ppid",
     171                        type=none_or_int, default=None,
    111172                        help="specify parent process id for \
    112173                        rerun of prepare_flexpart")
     
    114175    # arguments for job submission to ECMWF, only needed by submit.py
    115176    parser.add_argument("--job_template", dest='job_template',
    116                         default="job.temp",
     177                        type=none_or_str, default="job.temp",
    117178                        help="job template file for submission to ECMWF")
    118     parser.add_argument("--queue", dest="queue", default=None,
     179    parser.add_argument("--queue", dest="queue",
     180                        type=none_or_str, default=None,
    119181                        help="queue for submission to ECMWF \
    120182                        (e.g. ecgate or cca )")
    121     parser.add_argument("--controlfile", dest="controlfile",
    122                         default='CONTROL.temp',
    123                         help="file with CONTROL parameters")
    124     parser.add_argument("--debug", dest="debug", default=None,
    125                         help="debug mode - leave temporary files intact")
    126     parser.add_argument("--request", dest="request", default=None,
    127                         help="list all mars request in file mars_requests.dat \
    128                         and skip submission to mars")
    129183
    130184    args = parser.parse_args()
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG