Changeset ca867de in flex_extract.git for source/python/classes/EcFlexpart.py


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

refactored functions in EcFlexpart? and did some minor changes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG