Changeset 2fb99de in flex_extract.git for python/EcFlexpart.py


Ignore:
Timestamp:
Sep 20, 2018, 11:56:37 AM (6 years ago)
Author:
Anne Philipp <anne.philipp@…>
Branches:
master, ctbto, dev
Children:
5d42acd
Parents:
3232589
Message:

introduced config with path definitions and changed py files accordingly; Installation works; some tests were added for tarball making; Problems in submission to ecgate

File:
1 edited

Legend:

Unmodified
Added
Removed

  • python/EcFlexpart.py

    rdda0e9d r2fb99de  
    102102
    103103# software specific classes and modules from flex_extract
     104import _config
    104105from GribTools import GribTools
    105106from tools import init128, to_param_id, silent_remove, product, my_error
    … …  
    204205        self.accuracy = c.accuracy
    205206        self.level = c.level
    206 
    207         if c.levelist:
    208             self.levelist = c.levelist
    209         else:
    210             self.levelist = '1/to/' + c.level
    211 
     207        self.expver = c.expver
    212208        # for gaussian grid retrieval
    213209        self.glevelist = '1/to/' + c.level
    … …  
    217213        else:
    218214            self.gaussian = ''
    219 
    220         if hasattr(c, 'expver') and c.expver:
    221             self.expver = c.expver
    222         else:
    223             self.expver = '1'
    224 
    225         if hasattr(c, 'number') and c.number:
    226             self.number = c.number
    227         else:
    228             self.number = '0'
    229215
    230216        if 'N' in c.grid:  # Gaussian output grid
    … …  
    250236        # 4) Download also data for WRF
    251237
    252 
    253238        # Different grids need different retrievals
    254239        # SH = Spherical Harmonics, GG = Gaussian Grid,
    … …  
    274259            self.params['OG__ML'] = ['T/Q', 'ML', self.levelist, self.grid]
    275260
    276             if c.gauss == '0' and c.eta == '1':
     261            #if c.gauss == '0' and c.eta == '1':
     262            if not c.gauss and c.eta:
    277263                # the simplest case
    278264                self.params['OG__ML'][0] += '/U/V/77'
    279             elif c.gauss == '0' and c.eta == '0':
     265            #elif c.gauss == '0' and c.eta == '0':
     266            elif not c.gauss and not c.eta:
    280267            # this is not recommended (inaccurate)
    281268                self.params['OG__ML'][0] += '/U/V'
    282             elif c.gauss == '1' and c.eta == '0':
     269            #elif c.gauss == '1' and c.eta == '0':
     270            elif c.gauss and not c.eta:
    283271                # this is needed for data before 2008, or for reanalysis data
    284272                self.params['GG__SL'] = ['Q', 'ML', '1', \
    … …  
    286274                self.params['SH__ML'] = ['U/V/D', 'ML', self.glevelist, 'OFF']
    287275            else:
    288                 print 'Warning: This is a very costly parameter combination, \
    289                        use only for debugging!'
     276                print('Warning: This is a very costly parameter combination, \
     277                      use only for debugging!')
    290278                self.params['GG__SL'] = ['Q', 'ML', '1', \
    291279                                         '{}'.format((int(self.resol) + 1) / 2)]
    … …  
    293281                                         '{}'.format((int(self.resol) + 1) / 2)]
    294282
    295             if hasattr(c, 'omega') and c.omega == '1':
     283            if c.omega:
    296284                self.params['OG__ML'][0] += '/W'
    297285
    298286            # add cloud water content if necessary
    299             if hasattr(c, 'cwc') and c.cwc == '1':
     287            if c.cwc:
    300288                self.params['OG__ML'][0] += '/CLWC/CIWC'
    301289
    302290            # add vorticity and geopotential height for WRF if necessary
    303             if hasattr(c, 'wrf') and c.wrf == '1':
     291            if c.wrf:
    304292                self.params['OG__ML'][0] += '/Z/VO'
    305293                if '/D' not in self.params['OG__ML'][0]:
    … …  
    367355            f.write('&NAMGEN\n')
    368356            f.write(',\n  '.join(['maxl = ' + str(maxl), 'maxb = ' + str(maxb),
    369                                   'mlevel = ' + self.level,
    370                                   'mlevelist = ' + '"' + self.levelist + '"',
    371                                   'mnauf = ' + self.resol,
     357                                  'mlevel = ' + str(self.level),
     358                                  'mlevelist = ' + '"' + str(self.levelist)
     359                                                 + '"',
     360                                  'mnauf = ' + str(self.resol),
    372361                                  'metapar = ' + '77',
    373362                                  'rlo0 = ' + str(area[1]),
    … …  
    375364                                  'rla0 = ' + str(area[2]),
    376365                                  'rla1 = ' + str(area[0]),
    377                                   'momega = ' + c.omega,
    378                                   'momegadiff = ' + c.omegadiff,
    379                                   'mgauss = ' + c.gauss,
    380                                   'msmooth = ' + c.smooth,
    381                                   'meta = ' + c.eta,
    382                                   'metadiff = ' + c.etadiff,
    383                                   'mdpdeta = ' + c.dpdeta]))
     366                                  'momega = ' + str(c.omega),
     367                                  'momegadiff = ' + str(c.omegadiff),
     368                                  'mgauss = ' + str(c.gauss),
     369                                  'msmooth = ' + str(c.smooth),
     370                                  'meta = ' + str(c.eta),
     371                                  'metadiff = ' + str(c.etadiff),
     372                                  'mdpdeta = ' + str(c.dpdeta)]))
    384373
    385374            f.write('\n/\n')
    … …  
    387376        return
    388377
    389     def retrieve(self, server, dates, inputdir='.'):
     378    def retrieve(self, server, dates, request, inputdir='.'):
    390379        '''
    391380        @Description:
    … …  
    437426                    pass
    438427                if pk == 'OG_OROLSM__SL':
    439                     if oro is False:
     428                    if not oro:
    440429                        mfstream = 'OPER'
    441430                        mftype = 'AN'
    … …  
    457446
    458447    # ------  on demand path  --------------------------------------------------
    459                 if self.basetime is None:
     448                if not self.basetime:
    460449                    MR = MarsRetrieval(self.server,
    461450                                       marsclass=self.marsclass,
    … …  
    477466                                       param=pv[0])
    478467
    479                     MR.display_info()
    480                     MR.data_retrieve()
     468                    if request == 0:
     469                        MR.display_info()
     470                        MR.data_retrieve()
     471                    elif request == 1:
     472                        MR.print_info()
     473                    elif request == 2:
     474                        MR.print_info()
     475                        MR.display_info()
     476                        MR.data_retrieve()
     477                    else:
     478                        print('Failure')
    481479    # ------  operational path  ------------------------------------------------
    482480                else:
    … …  
    566564
    567565                            MR.display_info()
     566
    568567                            MR.data_retrieve()
    569568                        # --------------  non flux data ------------------------
    … …  
    662661                            MR.data_retrieve()
    663662
    664         print "MARS retrieve done... "
     663        if request == 0 or request == 2:
     664            print('MARS retrieve done ... ')
     665        elif request == 1:
     666            print('MARS request printed ...')
    665667
    666668        return
    … …  
    702704        '''
    703705
    704         print '\n\nPostprocessing:\n Format: {}\n'.format(c.format)
    705 
    706         if c.ecapi is False:
     706        print('\n\nPostprocessing:\n Format: {}\n'.format(c.format))
     707
     708        if not c.ecapi:
    707709            print('ecstorage: {}\n ecfsdir: {}\n'.
    708710                  format(c.ecstorage, c.ecfsdir))
    709             if not hasattr(c, 'gateway'):
    710                 c.gateway = os.getenv('GATEWAY')
    711             if not hasattr(c, 'destination'):
    712                 c.destination = os.getenv('DESTINATION')
     711            #if not hasattr(c, 'gateway'):
     712            #    c.gateway = os.getenv('GATEWAY')
     713            #if not hasattr(c, 'destination'):
     714            #    c.destination = os.getenv('DESTINATION')
    713715            print('ectrans: {}\n gateway: {}\n destination: {}\n '
    714716                  .format(c.ectrans, c.gateway, c.destination))
    715717
    716         print 'Output filelist: \n'
    717         print self.outputfilelist
     718        print('Output filelist: \n')
     719        print(self.outputfilelist)
    718720
    719721        if c.format.lower() == 'grib2':
    … …  
    724726                p = subprocess.check_call(['mv', ofile + '_2', ofile])
    725727
    726         if int(c.ectrans) == 1 and c.ecapi is False:
     728        if c.ectrans and not c.ecapi:
    727729            for ofile in self.outputfilelist:
    728730                p = subprocess.check_call(['ectrans', '-overwrite', '-gateway',
    … …  
    731733                #print('ectrans:', p)
    732734
    733         if int(c.ecstorage) == 1 and c.ecapi is False:
     735        if c.ecstorage and not c.ecapi:
    734736            for ofile in self.outputfilelist:
    735737                p = subprocess.check_call(['ecp', '-o', ofile,
    … …  
    742744        # prepare environment for the grib2flexpart run
    743745        # to convert grib to flexpart binary
    744         if c.grib2flexpart == '1':
     746        if c.grib2flexpart:
    745747
    746748            # generate AVAILABLE file
    … …  
    858860
    859861        index_keys = ["date", "time", "step"]
    860         indexfile = c.inputdir + "/date_time_stepRange.idx"
     862        indexfile = os.path.join(c.inputdir, _config.FILE_GRIB_INDEX)
    861863        silent_remove(indexfile)
    862864        grib = GribTools(inputfiles.files)
    … …  
    868870        for key in index_keys:
    869871            index_vals.append(grib_index_get(iid, key))
    870             print index_vals[-1]
     872            print(index_vals[-1])
    871873            # index_vals looks for example like:
    872874            # index_vals[0]: ('20171106', '20171107', '20171108') ; date
    … …  
    878880
    879881        for prod in product(*index_vals):
     882            # e.g. prod = ('20170505', '0', '12')
     883            #             (  date    ,time, step)
     884            # per date e.g. time = 0, 1200
     885            # per time e.g. step = 3, 6, 9, 12
    880886            # flag for Fortran program CONVERT2 and file merging
    881887            convertFlag = False
    882             print 'current prod: ', prod
     888            print('current prod: ', prod)
    883889            # e.g. prod = ('20170505', '0', '12')
    884890            #             (  date    ,time, step)
    … …  
    902908                # they are just valid for a single product
    903909                for k, f in fdict.iteritems():
    904                     silent_remove(c.inputdir + "/fort." + k)
    905                     fdict[k] = open(c.inputdir + '/fort.' + k, 'w')
     910                    fortfile = os.path.join(c.inputdir, 'fort.' + k)
     911                    silent_remove(fortfile)
     912                    fdict[k] = open(fortfile, 'w')
    906913
    907914                cdate = str(grib_get(gid, 'date'))
    … …  
    931938                    if timestamp < slimit or timestamp > elimit:
    932939                        continue
     940            else:
     941                pass
    933942
    934943            try:
    935                 if c.wrf == '1':
    936                     if 'olddate' not in locals():
    937                         fwrf = open(c.outputdir + '/WRF' + cdate +
    938                                     '.{:0>2}'.format(time) + '.000.grb2', 'w')
     944                if c.wrf:
     945                    if 'olddate' not in locals() or cdate != olddate:
     946                        fwrf = open(os.path.join(c.outputdir,
     947                                    'WRF' + cdate + '.{:0>2}'.format(time) +
     948                                    '.000.grb2'), 'w')
    939949                        olddate = cdate[:]
    940                     else:
    941                         if cdate != olddate:
    942                             fwrf = open(c.outputdir + '/WRF' + cdate +
    943                                         '.{:0>2}'.format(time) + '.000.grb2',
    944                                         'w')
    945                             olddate = cdate[:]
    946950            except AttributeError:
    947951                pass
    948952
    949             # helper variable to remember which fields are already used.
     953            # helper variable to remember which fields were already used.
    950954            savedfields = []
    951955            while 1:
    … …  
    958962                # Specific humidity (Q.grb) is used as a template only
    959963                # so we need the first we "meet"
    960                     with open(c.inputdir + '/fort.18', 'w') as fout:
     964                    with open(os.path.join(c.inputdir, 'fort.18'), 'w') as fout:
    961965                        grib_write(gid, fout)
    962966                elif paramId == 131 or paramId == 132:
    … …  
    971975                    grib_write(gid, fdict['13'])
    972976                elif  paramId in [129, 138, 155] and levtype == 'hybrid' \
    973                         and c.wrf == '1':
     977                        and c.wrf:
    974978                    pass
    975979                elif paramId == 246 or paramId == 247:
    … …  
    991995                        savedfields.append(paramId)
    992996                    else:
    993                         print 'duplicate ' + str(paramId) + ' not written'
     997                        print('duplicate ' + str(paramId) + ' not written')
    994998
    995999                try:
    996                     if c.wrf == '1':
    997                         if levtype == 'hybrid': # model layer
    998                             if paramId in [129, 130, 131, 132, 133, 138, 155]:
    999                                 grib_write(gid, fwrf)
    1000                         else: # sfc layer
    1001                             if paramId in wrfpars:
    1002                                 grib_write(gid, fwrf)
     1000                    if c.wrf:
     1001                        # model layer
     1002                        if levtype == 'hybrid' and \
     1003                           paramId in [129, 130, 131, 132, 133, 138, 155]:
     1004                            grib_write(gid, fwrf)
     1005                        # sfc layer
     1006                        elif paramId in wrfpars:
     1007                            grib_write(gid, fwrf)
    10031008                except AttributeError:
    10041009                    pass
    … …  
    10141019                pwd = os.getcwd()
    10151020                os.chdir(c.inputdir)
    1016                 if os.stat('fort.21').st_size == 0 and int(c.eta) == 1:
    1017                     print 'Parameter 77 (etadot) is missing, most likely it is \
    1018                            not available for this type or date/time\n'
    1019                     print 'Check parameters CLASS, TYPE, STREAM, START_DATE\n'
     1021                if os.stat('fort.21').st_size == 0 and c.eta:
     1022                    print('Parameter 77 (etadot) is missing, most likely it is \
     1023                           not available for this type or date/time\n')
     1024                    print('Check parameters CLASS, TYPE, STREAM, START_DATE\n')
    10201025                    my_error(c.mailfail, 'fort.21 is empty while parameter eta \
    10211026                             is set to 1 in CONTROL file')
    … …  
    10231028                # create the corresponding output file fort.15
    10241029                # (generated by CONVERT2) + fort.16 (paramId 167 and 168)
    1025                 p = subprocess.check_call(
    1026                     [os.path.expandvars(os.path.expanduser(c.exedir)) +
    1027                      '/CONVERT2'], shell=True)
     1030                p = subprocess.check_call([os.path.join(c.exedir, 'CONVERT2')],
     1031                                          shell=True)
    10281032                os.chdir(pwd)
    10291033
    10301034                # create final output filename, e.g. EN13040500 (ENYYMMDDHH)
    1031                 fnout = c.inputdir + '/' + c.prefix
     1035                fnout = os.path.join(c.inputdir, c.prefix)
    10321036                if c.maxstep > 12:
    10331037                    suffix = cdate[2:8] + '.{:0>2}'.format(time/100) + \
    … …  
    10361040                    suffix = cdateH[2:10]
    10371041                fnout += suffix
    1038                 print "outputfile = " + fnout
     1042                print("outputfile = " + fnout)
    10391043                self.outputfilelist.append(fnout) # needed for final processing
    10401044
    … …  
    10441048                    c.inputdir + '/OG_OROLSM__SL.*.' + c.ppid + '*')[0])
    10451049                fluxfile = 'flux' + cdate[0:2] + suffix
    1046                 if c.cwc != '1':
     1050                if not c.cwc:
    10471051                    flist = ['fort.15', fluxfile, 'fort.16', orolsm]
    10481052                else:
    … …  
    10511055                with open(fnout, 'wb') as fout:
    10521056                    for f in flist:
     1057                        shutil.copyfileobj(open(os.path.join(c.inputdir, f),
     1058                                                'rb'), fout)
     1059
     1060                if c.omega:
     1061                    with open(os.path.join(c.outputdir, 'OMEGA'), 'wb') as fout:
    10531062                        shutil.copyfileobj(
    1054                             open(c.inputdir + '/' + f, 'rb'), fout)
    1055 
    1056                 if c.omega == '1':
    1057                     with open(c.outputdir + '/OMEGA', 'wb') as fout:
    1058                         shutil.copyfileobj(
    1059                             open(c.inputdir + '/fort.25', 'rb'), fout)
    1060 
    1061         if hasattr(c, 'wrf') and c.wrf == '1':
     1063                            open(os.path.join(c.inputdir, 'fort.25'),
     1064                                 'rb'), fout)
     1065            else:
     1066                pass
     1067
     1068        if c.wrf:
    10621069            fwrf.close()
    10631070
    … …  
    11021109        pars = to_param_id(self.params['OG_acc_SL'][0], table128)
    11031110        index_keys = ["date", "time", "step"]
    1104         indexfile = c.inputdir + "/date_time_stepRange.idx"
     1111        indexfile = os.path.join(c.inputdir, _config.FILE_GRIB_INDEX)
    11051112        silent_remove(indexfile)
    11061113        grib = GribTools(inputfiles.files)
    … …  
    11121119        for key in index_keys:
    11131120            key_vals = grib_index_get(iid, key)
    1114             print key_vals
     1121            print(key_vals)
    11151122            # have to sort the steps for disaggregation,
    11161123            # therefore convert to int first
    … …  
    11331140            stepsdict[str(p)] = []
    11341141
    1135         print 'maxstep: ', c.maxstep
     1142        print('maxstep: ', c.maxstep)
    11361143
    11371144        for prod in product(*index_vals):
    … …  
    11401147            # per date e.g. time = 0, 1200
    11411148            # per time e.g. step = 3, 6, 9, 12
     1149            print('current prod: ', prod)
    11421150            for i in range(len(index_keys)):
    11431151                grib_index_select(iid, index_keys[i], prod[i])
    11441152
     1153            # get first id from current product
    11451154            gid = grib_new_from_index(iid)
     1155
     1156            # if there is data for this product combination
     1157            # prepare some date and time parameter before reading the data
    11461158            if gid is not None:
    11471159                cdate = grib_get(gid, 'date')
    … …  
    11611173
    11621174            if c.maxstep > 12:
    1163                 fnout = c.inputdir + '/flux' + \
    1164                     sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \
    1165                     '.{:0>3}'.format(step-2*int(c.dtime))
    1166                 gnout = c.inputdir + '/flux' + \
    1167                     sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \
    1168                     '.{:0>3}'.format(step-int(c.dtime))
    1169                 hnout = c.inputdir + '/flux' + \
    1170                     sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \
    1171                     '.{:0>3}'.format(step)
    1172                 g = open(gnout, 'w')
    1173                 h = open(hnout, 'w')
     1175                fnout = os.path.join(c.inputdir, 'flux' +
     1176                                     sdate.strftime('%Y%m%d') +
     1177                                     '.{:0>2}'.format(time/100) +
     1178                                     '.{:0>3}'.format(step-2*int(c.dtime)))
     1179                gnout = os.path.join(c.inputdir, 'flux' +
     1180                                     sdate.strftime('%Y%m%d') +
     1181                                     '.{:0>2}'.format(time/100) +
     1182                                     '.{:0>3}'.format(step-int(c.dtime)))
     1183                hnout = os.path.join(c.inputdir, 'flux' +
     1184                                     sdate.strftime('%Y%m%d') +
     1185                                     '.{:0>2}'.format(time/100) +
     1186                                     '.{:0>3}'.format(step))
    11741187            else:
    1175                 fnout = c.inputdir + '/flux' + fdate.strftime('%Y%m%d%H')
    1176                 gnout = c.inputdir + '/flux' + (fdate +
    1177                                                 timedelta(hours=int(c.dtime))
    1178                                                ).strftime('%Y%m%d%H')
    1179                 hnout = c.inputdir + '/flux' + sdates.strftime('%Y%m%d%H')
    1180                 g = open(gnout, 'w')
    1181                 h = open(hnout, 'w')
    1182 
    1183             print "outputfile = " + fnout
    1184             f = open(fnout, 'w')
     1188                fnout = os.path.join(c.inputdir, 'flux' +
     1189                                     fdate.strftime('%Y%m%d%H'))
     1190                gnout = os.path.join(c.inputdir, 'flux' +
     1191                                     (fdate + timedelta(hours=int(c.dtime))).
     1192                                     strftime('%Y%m%d%H'))
     1193                hnout = os.path.join(c.inputdir, 'flux' +
     1194                                     sdates.strftime('%Y%m%d%H'))
     1195
     1196            print("outputfile = " + fnout)
     1197            f_handle = open(fnout, 'w')
     1198            g_handle = open(gnout, 'w')
     1199            h_handle = open(hnout, 'w')
    11851200
    11861201            # read message for message and store relevant data fields
    … …  
    12431258                            grib_set(gid, 'date', fdate.year*10000 +
    12441259                                     fdate.month*100+fdate.day)
    1245                         grib_write(gid, f)
     1260                        grib_write(gid, f_handle)
    12461261
    12471262                        if c.basetime is not None:
    … …  
    12701285                                     truedatetime.month * 100 +
    12711286                                     truedatetime.day)
    1272                             grib_write(gid, h)
     1287                            grib_write(gid, h_handle)
    12731288
    12741289                            #values = (svdp[1]+svdp[2])/2.
    … …  
    12851300                                     truedatetime.day)
    12861301                            grib_set_values(gid, values)
    1287                             grib_write(gid, g)
     1302                            grib_write(gid, g_handle)
    12881303
    12891304                    grib_release(gid)
    … …  
    12911306                    gid = grib_new_from_index(iid)
    12921307
    1293             f.close()
    1294             g.close()
    1295             h.close()
     1308            f_handle.close()
     1309            g_handle.close()
     1310            h_handle.close()
    12961311
    12971312        grib_index_release(iid)
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG