- Timestamp:
- Feb 8, 2018, 9:54:05 PM (6 years ago)
- Branches:
- master, ctbto, dev
- Children:
- 02c8c50
- Parents:
- 6180177
- Location:
- python
- Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
python/FlexpartTools.py
rd69b677 r64cf353 1 # 2 # This software is licensed under the terms of the Apache Licence Version 2.0 3 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 4 # 5 # Functionality provided: Prepare input 3D-wind fields in hybrid coordinates + surface fields for FLEXPART runs 6 # 7 # Creation: October 2014 - Anne Fouilloux - University of Oslo 8 # Extension November 2015 - Leopold Haimberger - University of Vienna for: 9 # - using the WebAPI also for general MARS retrievals 10 # - job submission on ecgate and cca 11 # - job templates suitable for twice daily operational dissemination 12 # - dividing retrievals of longer periods into digestable chunks 13 # - retrieve also longer term forecasts, not only analyses and short term forecast data 14 # - conversion into GRIB2 15 # - conversion into .fp format for faster execution of FLEXPART 16 # 17 # Requirements: 18 # in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed 19 # ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/ 20 # dateutils 21 # matplotlib (optional, for debugging) 22 # 23 # 1 #!/usr/bin/env python 2 # -*- coding: utf-8 -*- 3 #************************************************************************ 4 # TODO AP 5 #AP 6 # - Functionality Provided is not correct for this file 7 # - localpythonpath should not be set in module load section! 8 # - Change History ist nicht angepasst ans File! Original geben lassen 9 # - def myerror muss angepasst werden da derzeit manuelle modifikation notwendig 10 # - def --init-- marsretrieval class, remove dataset and expect?! 11 # - the EIFlexpart class is jsut for EI ? So we should maybe create classes 12 # for each possible data set so the boundarys of the attributes can be 13 # validated! 14 #************************************************************************ 15 """ 16 @Author: Anne Fouilloux (University of Oslo) 17 18 @Date: October 2014 19 20 @ChangeHistory: 21 November 2015 - Leopold Haimberger (University of Vienna): 22 - using the WebAPI also for general MARS retrievals 23 - job submission on ecgate and cca 24 - job templates suitable for twice daily operational dissemination 25 - dividing retrievals of longer periods into digestable chunks 26 - retrieve also longer term forecasts, not only analyses and 27 short term forecast data 28 - conversion into GRIB2 29 - conversion into .fp format for faster execution of FLEXPART 30 31 February 2018 - Anne Philipp (University of Vienna): 32 - applied PEP8 style guide 33 - added documentation 34 35 @License: 36 (C) Copyright 2014 UIO. 37 38 This software is licensed under the terms of the Apache Licence Version 2.0 39 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 40 41 @Requirements: 42 - A standard python 2.6 or 2.7 installation 43 - dateutils 44 - matplotlib (optional, for debugging) 45 - ECMWF specific packages, all available from https://software.ecmwf.int/ 46 ECMWF WebMARS, gribAPI with python enabled, emoslib and 47 ecaccess web toolkit 48 49 @Description: 50 Further documentation may be obtained from www.flexpart.eu. 51 52 Functionality provided: 53 Prepare input 3D-wind fields in hybrid coordinates + 54 surface fields for FLEXPART runs 55 """ 56 # ------------------------------------------------------------------------------ 57 # MODULES 58 # ------------------------------------------------------------------------------ 24 59 import subprocess 25 from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter60 from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter 26 61 import traceback 27 62 import shutil 28 import os, errno,sys,inspect,glob 63 import os 64 import errno 65 import sys 66 import inspect 67 import glob 29 68 import datetime 30 #import re31 #import copy32 33 69 from string import atoi 34 from numpy import * 35 36 #sys.stdout=open('ECMWFDATA_'+str(os.getpid()),'w') 37 #from ecmwfapi import ECMWFDataServer 38 ecapi=True 70 from numpy import * 71 ecapi = True 39 72 try: 40 73 import ecmwfapi 41 74 except ImportError: 42 ecapi =False75 ecapi = False 43 76 from gribapi import * 44 77 from GribTools import GribTools 45 78 46 def interpret_args_and_control(*args,**kwargs): 47 48 parser = ArgumentParser(description='Retrieve FLEXPART input from ECMWF MARS archive', 79 def interpret_args_and_control(): 80 ''' 81 @Description: 82 Assigns the command line arguments and reads control file 83 content. Apply default values for non mentioned arguments. 84 85 @Input: 86 <nothing> 87 88 @Return: 89 args: instance of ArgumentParser 90 Contains the commandline arguments from script/program call. 91 92 c: instance of class Control 93 Contains all necessary information of a control file. The parameters 94 are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM, 95 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST, 96 RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, 97 SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, 98 ECFSDIR, MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR 99 For more information about format and content of the parameter see 100 documentation. 101 102 ''' 103 parser = ArgumentParser(description='Retrieve FLEXPART input from \ 104 ECMWF MARS archive', 49 105 formatter_class=ArgumentDefaultsHelpFormatter) 50 106 51 # the most important arguments 107 # the most important arguments 52 108 parser.add_argument("--start_date", dest="start_date", 53 109 help="start date YYYYMMDD") 54 parser.add_argument( 55 56 parser.add_argument( "--date_chunk", dest="date_chunk",default=None,57 58 110 parser.add_argument("--end_date", dest="end_date", 111 help="end_date YYYYMMDD") 112 parser.add_argument("--date_chunk", dest="date_chunk", default=None, 113 help="# of days to be retrieved at once") 114 59 115 # some arguments that override the default in the control file 60 parser.add_argument("--basetime", dest="basetime", help="base such as 00/12 (for half day retrievals)") 61 parser.add_argument("--step", dest="step", help="steps such as 00/to/48") 62 parser.add_argument("--levelist", dest="levelist",help="Vertical levels to be retrieved, e.g. 30/to/60") 63 parser.add_argument("--area", dest="area", help="area defined as north/west/south/east") 64 116 parser.add_argument("--basetime", dest="basetime", 117 help="base such as 00/12 (for half day retrievals)") 118 parser.add_argument("--step", dest="step", 119 help="steps such as 00/to/48") 120 parser.add_argument("--levelist", dest="levelist", 121 help="Vertical levels to be retrieved, e.g. 30/to/60") 122 parser.add_argument("--area", dest="area", 123 help="area defined as north/west/south/east") 124 65 125 # set the working directories 66 parser.add_argument("--inputdir", dest="inputdir", default=None,126 parser.add_argument("--inputdir", dest="inputdir", default=None, 67 127 help="root directory for storing intermediate files") 68 parser.add_argument("--outputdir", dest="outputdir", default=None,128 parser.add_argument("--outputdir", dest="outputdir", default=None, 69 129 help="root directory for storing output files") 70 130 parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts", 71 help="FLEXPART root directory (to find grib2flexpart and COMMAND file)\n\ 72 Normally ECMWFDATA resides in the scripts directory of the FLEXPART distribution") 73 131 help="FLEXPART root directory (to find grib2flexpart \ 132 and COMMAND file)\n\ Normally ECMWFDATA resides in \ 133 the scripts directory of the FLEXPART distribution") 134 74 135 # this is only used by prepareFLEXPART.py to rerun a postprocessing step 75 136 parser.add_argument("--ppid", dest="ppid", 76 help="Specify parent process id for rerun of prepareFLEXPART") 137 help="Specify parent process id for \ 138 rerun of prepareFLEXPART") 77 139 78 140 # arguments for job submission to ECMWF, only needed by submit.py 79 parser.add_argument("--job_template", dest='job_template',default="job.temp", 141 parser.add_argument("--job_template", dest='job_template', 142 default="job.temp", 80 143 help="job template file for submission to ECMWF") 81 144 #parser.add_argument("--remote", dest="remote", 82 #help="target for submission to ECMWF (ecgate or cca etc.)") 145 #help="target for submission to ECMWF \ 146 #(ecgate or cca etc.)") 83 147 parser.add_argument("--queue", dest="queue", 84 help="queue for submission to ECMWF (e.g. ecgate or cca )")85 86 87 parser.add_argument("--controlfile", dest="controlfile",default='CONTROL.temp',148 help="queue for submission to ECMWF \ 149 (e.g. ecgate or cca )") 150 parser.add_argument("--controlfile", dest="controlfile", 151 default='CONTROL.temp', 88 152 help="file with control parameters") 89 parser.add_argument("--debug", dest="debug", default=0,153 parser.add_argument("--debug", dest="debug", default=0, 90 154 help="Debug mode - leave temporary files intact") 155 91 156 args = parser.parse_args() 92 157 93 158 try: 94 c=Control(args.controlfile)159 c = Control(args.controlfile) 95 160 except IOError: 96 try: 97 c=Control(localpythonpath+args.controlfile) 98 99 except: 100 print 'Could not read control file "'+args.controlfile+'"' 101 print 'Either it does not exist or its syntax is wrong.' 102 print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information' 103 exit(1) 104 105 if args.start_date==None and getattr(c,'start_date')==None: 106 print 'start_date specified neither in command line nor in control file '+args.controlfile 107 print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information' 108 exit(1) 109 110 if args.start_date!=None: 111 c.start_date=args.start_date 112 if args.end_date!=None: 113 c.end_date=args.end_date 114 if c.end_date==None: 115 c.end_date=c.start_date 116 if args.date_chunk!=None: 117 c.date_chunk=args.date_chunk 118 119 if not hasattr(c,'debug'): 120 c.debug=args.debug 121 122 123 if args.inputdir==None and args.outputdir==None: 124 c.inputdir='../work' 125 c.outputdir='../work' 161 try: 162 c = Control(localpythonpath + args.controlfile) 163 except: 164 print('Could not read control file "' + args.controlfile + '"') 165 print('Either it does not exist or its syntax is wrong.') 166 print('Try "' + sys.argv[0].split('/')[-1] + 167 ' -h" to print usage information') 168 exit(1) 169 170 if args.start_date is None and getattr(c, 'start_date') is None: 171 print('start_date specified neither in command line nor \ 172 in control file ' + args.controlfile) 173 print('Try "' + sys.argv[0].split('/')[-1] + 174 ' -h" to print usage information') 175 exit(1) 176 177 if args.start_date is not None: 178 c.start_date = args.start_date 179 if args.end_date is not None: 180 c.end_date = args.end_date 181 if c.end_date is None: 182 c.end_date = c.start_date 183 if args.date_chunk is not None: 184 c.date_chunk = args.date_chunk 185 186 if not hasattr(c, 'debug'): 187 c.debug = args.debug 188 189 if args.inputdir is None and args.outputdir is None: 190 c.inputdir = '../work' 191 c.outputdir = '../work' 126 192 else: 127 if args.inputdir!=None:128 c.inputdir=args.inputdir129 if args.outputdir==None:130 c.outputdir=args.inputdir131 if args.outputdir!=None:132 c.outputdir=args.outputdir133 if args.inputdir==None:134 c.inputdir=args.outputdir135 136 if hasattr(c, 'outputdir')==False and args.outputdir==None:137 c.outputdir=c.inputdir193 if args.inputdir is not None: 194 c.inputdir = args.inputdir 195 if args.outputdir is None: 196 c.outputdir = args.inputdir 197 if args.outputdir is not None: 198 c.outputdir = args.outputdir 199 if args.inputdir is None: 200 c.inputdir = args.outputdir 201 202 if hasattr(c, 'outputdir') is False and args.outputdir is None: 203 c.outputdir = c.inputdir 138 204 else: 139 if args.outputdir!=None: 140 c.outputdir=args.outputdir 141 142 if args.area!=None: 143 afloat='.' in args.area 144 l=args.area.split('/') 145 if afloat: 146 for i in range(len(l)): 147 l[i]=str(int(float(l[i])*1000)) 148 c.upper,c.left,c.lower,c.right=l 149 150 # basetime aktiviert den ''operational mode'' 151 if args.basetime!=None: 152 c.basetime=args.basetime 153 #if int(c.basetime)==0: 154 #c.start_date=datetime.datetime.strftime(datetime.datetime.strptime(c.start_date,'%Y%m%d')-datetime.timedelta(days=1),'%Y%m%d') 155 #c.end_date=c.start_date 156 157 if args.step!=None: 158 l=args.step.split('/') 159 if 'to' in args.step.lower(): 160 if 'by' in args.step.lower(): 161 ilist=arange(int(l[0]),int(l[2])+1,int(l[4])) 162 c.step=['{:0>3}'.format(i) for i in ilist] 163 else: 164 myerror(None,args.step+':\n'+'please use "by" as well if "to" is used') 165 else: 166 c.step=l 167 168 if args.levelist!=None: 169 c.levelist=args.levelist 170 if 'to' in c.levelist: 171 c.level=c.levelist.split('/')[2] 172 else: 173 c.level=c.levelist.split('/')[-1] 174 175 176 if args.flexpart_root_scripts!=None: 177 c.flexpart_root_scripts=args.flexpart_root_scripts 178 179 return args,c 205 if args.outputdir is not None: 206 c.outputdir = args.outputdir 207 208 if args.area is not None: 209 afloat = '.' in args.area 210 l = args.area.split('/') 211 if afloat: 212 for i in range(len(l)): 213 l[i] = str(int(float(l[i]) * 1000)) 214 c.upper, c.left, c.lower, c.right = l 215 216 # basetime aktiviert den ''operational mode'' 217 if args.basetime is not None: 218 c.basetime = args.basetime 219 220 if args.step is not None: 221 l = args.step.split('/') 222 if 'to' in args.step.lower(): 223 if 'by' in args.step.lower(): 224 ilist = arange(int(l[0]), int(l[2]) + 1, int(l[4])) 225 c.step = ['{:0>3}'.format(i) for i in ilist] 226 else: 227 myerror(None, args.step + ':\n' + 228 'please use "by" as well if "to" is used') 229 else: 230 c.step = l 231 232 if args.levelist is not None: 233 c.levelist = args.levelist 234 if 'to' in c.levelist: 235 c.level = c.levelist.split('/')[2] 236 else: 237 c.level = c.levelist.split('/')[-1] 238 239 if args.flexpart_root_scripts is not None: 240 c.flexpart_root_scripts = args.flexpart_root_scripts 241 242 return args, c 243 180 244 181 245 def install_args_and_control(): 182 183 parser = ArgumentParser(description='Install ECMWFDATA software locally or on ECMWF machines', 246 ''' 247 @Description: 248 Assigns the command line arguments for installation and reads 249 control file content. Apply default values for non mentioned arguments. 250 251 @Input: 252 <nothing> 253 254 @Return: 255 args: instance of ArgumentParser 256 Contains the commandline arguments from script/program call. 257 258 c: instance of class Control 259 Contains all necessary information of a control file. The parameters 260 are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM, 261 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST, 262 RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, 263 SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, 264 ECFSDIR, MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR 265 For more information about format and content of the parameter see 266 documentation. 267 ''' 268 parser = ArgumentParser(description='Install ECMWFDATA software locally or \ 269 on ECMWF machines', 184 270 formatter_class=ArgumentDefaultsHelpFormatter) 185 parser.add_argument('--target',dest='install_target', 186 help="Valid targets: local | ecgate | cca , the latter two are at ECMWF") 187 parser.add_argument("--makefile", dest="makefile",help='Name of Makefile to use for compiling CONVERT2') 188 parser.add_argument("--ecuid", dest="ecuid",help='user id at ECMWF') 189 parser.add_argument("--ecgid", dest="ecgid",help='group id at ECMWF') 190 parser.add_argument("--gateway", dest="gateway",help='name of local gateway server') 191 parser.add_argument("--destination", dest="destination",help='ecaccess destination, e.g. leo@genericSftp') 271 272 parser.add_argument('--target', dest='install_target', 273 help="Valid targets: local | ecgate | cca , \ 274 the latter two are at ECMWF") 275 parser.add_argument("--makefile", dest="makefile", 276 help='Name of Makefile to use for compiling CONVERT2') 277 parser.add_argument("--ecuid", dest="ecuid", 278 help='user id at ECMWF') 279 parser.add_argument("--ecgid", dest="ecgid", 280 help='group id at ECMWF') 281 parser.add_argument("--gateway", dest="gateway", 282 help='name of local gateway server') 283 parser.add_argument("--destination", dest="destination", 284 help='ecaccess destination, e.g. leo@genericSftp') 192 285 193 286 parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts", 194 help="FLEXPART root directory on ECMWF servers (to find grib2flexpart and COMMAND file)\n\ 195 Normally ECMWFDATA resides in the scripts directory of the FLEXPART distribution, thus the:") 196 287 help="FLEXPART root directory on ECMWF servers \ 288 (to find grib2flexpart and COMMAND file)\n\ 289 Normally ECMWFDATA resides in the scripts directory \ 290 of the FLEXPART distribution, thus the:") 291 197 292 # arguments for job submission to ECMWF, only needed by submit.py 198 parser.add_argument("--job_template", dest='job_template',default="job.temp.o", 293 parser.add_argument("--job_template", dest='job_template', 294 default="job.temp.o", 199 295 help="job template file for submission to ECMWF") 200 296 201 parser.add_argument("--controlfile", dest="controlfile",default='CONTROL.temp', 297 parser.add_argument("--controlfile", dest="controlfile", 298 default='CONTROL.temp', 202 299 help="file with control parameters") 300 203 301 args = parser.parse_args() 204 302 205 303 try: 206 # if True: 207 c=Control(args.controlfile) 208 # else: 304 c = Control(args.controlfile) 209 305 except: 210 print 'Could not read control file "'+args.controlfile+'"' 211 print 'Either it does not exist or its syntax is wrong.' 212 print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information' 213 exit(1) 214 215 if args.install_target!='local': 216 if args.ecgid==None or args.ecuid==None or args.gateway==None or args.destination==None : 217 print 'Please enter your ECMWF user id and group id as well as the \nname of the local gateway and the ectrans destination ' 218 print 'with command line options --ecuid --ecgid --gateway --destination' 219 print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information' 220 print 'Please consult ecaccess documentation or ECMWF user support for further details' 221 sys.exit(1) 222 else: 223 c.ecuid=args.ecuid 224 c.ecgid=args.ecgid 225 c.gateway=args.gateway 226 c.destination=args.destination 227 306 print('Could not read control file "' + args.controlfile + '"') 307 print('Either it does not exist or its syntax is wrong.') 308 print('Try "' + sys.argv[0].split('/')[-1] + 309 ' -h" to print usage information') 310 exit(1) 311 312 if args.install_target != 'local': 313 if (args.ecgid is None or args.ecuid is None or args.gateway is None 314 or args.destination is None): 315 print('Please enter your ECMWF user id and group id as well as \ 316 the \nname of the local gateway and the ectrans \ 317 destination ') 318 print('with command line options --ecuid --ecgid \ 319 --gateway --destination') 320 print('Try "' + sys.argv[0].split('/')[-1] + 321 ' -h" to print usage information') 322 print('Please consult ecaccess documentation or ECMWF user support \ 323 for further details') 324 sys.exit(1) 325 else: 326 c.ecuid = args.ecuid 327 c.ecgid = args.ecgid 328 c.gateway = args.gateway 329 c.destination = args.destination 330 228 331 try: 229 c.makefile=args.makefile332 c.makefile = args.makefile 230 333 except: 231 pass 232 233 if args.install_target=='local': 234 if args.flexpart_root_scripts==None: 235 c.flexpart_root_scripts='../' 236 else: 237 c.flexpart_root_scripts=args.flexpart_root_scripts 238 239 if args.install_target!='local': 240 if args.flexpart_root_scripts==None: 241 c.ec_flexpart_root_scripts='${HOME}' 242 else: 243 c.ec_flexpart_root_scripts=args.flexpart_root_scripts 244 245 return args,c 334 pass 335 336 if args.install_target == 'local': 337 if args.flexpart_root_scripts is None: 338 c.flexpart_root_scripts = '../' 339 else: 340 c.flexpart_root_scripts = args.flexpart_root_scripts 341 342 if args.install_target != 'local': 343 if args.flexpart_root_scripts is None: 344 c.ec_flexpart_root_scripts = '${HOME}' 345 else: 346 c.ec_flexpart_root_scripts = args.flexpart_root_scripts 347 348 return args, c 349 246 350 247 351 def cleanup(c): 248 print "cleanup" 249 cleanlist=glob.glob(c.inputdir+"/*") 250 # cleanlist[-1:-1]=glob.glob(c.inputdir+"/flux*") 251 # cleanlist[-1:-1]=glob.glob(c.inputdir+"/*.grb") 252 # if c.ecapi==False and (c.ectrans=='1' or c.ecstorage=='1'): 253 # cleanlist[-1:-1]=glob.glob(c.inputdir+"/"+c.prefix+"*") 254 255 256 # silentremove(c.inputdir+"/VERTICAL.EC") 257 # silentremove(c.inputdir+"/date_time_stepRange.idx") 352 ''' 353 @Description: 354 Remove all files from intermediate directory 355 (inputdir from control file). 356 357 @Input: 358 c: instance of class Control 359 Contains all the parameters of control files, which are: 360 DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM, 361 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, 362 LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, 363 ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, 364 ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL, 365 GRIB2FLEXPART, FLEXPARTDIR 366 For more information about format and content of the parameter 367 see documentation. 368 369 @Return: 370 <nothing> 371 ''' 372 373 print("cleanup") 374 375 cleanlist = glob.glob(c.inputdir + "/*") 258 376 for cl in cleanlist: 259 if c.prefix not in cl: 260 silentremove(cl) 261 if c.ecapi==False and (c.ectrans=='1' or c.ecstorage=='1'): 262 silentremove(cl) 263 264 265 266 print "Done" 267 268 269 def myerror(c,message='ERROR'): 270 271 try: 272 target=c.mailfail 273 except AttributeError: 274 target=os.getenv('USER') 275 276 if(type(target) is not list): 277 target=[target] 278 print message 377 if c.prefix not in cl: 378 silentremove(cl) 379 if c.ecapi is False and (c.ectrans == '1' or c.ecstorage == '1'): 380 silentremove(cl) 381 382 print("Done") 383 384 return 385 386 387 def myerror(c, message='ERROR'): 388 ''' 389 @Description: 390 Print error message. 391 392 @Input: 393 c: instance of class Control 394 Contains all the parameters of control files, which are: 395 DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM, 396 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, 397 LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, 398 ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, 399 ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL, 400 GRIB2FLEXPART, FLEXPARTDIR 401 For more information about format and content of the parameter 402 see documentation. 403 404 message: string 405 Error message. Default value is "ERROR". 406 407 @Return: 408 <nothing> 409 ''' 410 # uncomment if user wants email notification directly from python 411 #try: 412 #target = c.mailfail 413 #except AttributeError: 414 #target = os.getenv('USER') 415 416 #if(type(target) is not list): 417 #target = [target] 418 419 print(message) 420 279 421 # uncomment if user wants email notification directly from python 280 422 #for t in target: 281 #p=subprocess.Popen(['mail','-s ECMWFDATA v7.0 ERROR', os.path.expandvars(t)], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE,bufsize=1) 282 #tr='\n'.join(traceback.format_stack()) 283 #pout=p.communicate(input=message+'\n\n'+tr)[0] 284 #print 'Email sent to '+os.path.expandvars(t) # +' '+pout.decode() 285 423 #p = subprocess.Popen(['mail','-s ECMWFDATA v7.0 ERROR', os.path.expandvars(t)], 424 # stdin = subprocess.PIPE, stdout = subprocess.PIPE, 425 # stderr = subprocess.PIPE, bufsize = 1) 426 #tr = '\n'.join(traceback.format_stack()) 427 #pout = p.communicate(input = message+'\n\n'+tr)[0] 428 #print 'Email sent to '+os.path.expandvars(t) # +' '+pout.decode() 429 286 430 exit(1) 287 288 def normalexit(c,message='Done!'):289 290 try:291 target=c.mailops292 293 294 if(type(target) is not list):295 target=[target]296 # Uncomment if user wants notification directly from python297 #for t in target:298 #p=subprocess.Popen(['mail','-s ECMWFDATA v7.0 normal exit', os.path.expandvars(t)], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE,bufsize=1)299 #pout=p.communicate(input=message+'\n\n')[0]300 #print pout.decode()301 302 except:303 pass304 305 print message306 431 307 432 return 308 433 434 435 def normalexit(c, message='Done!'): 436 ''' 437 @Description: 438 Print exit message. 439 440 @Input: 441 c: instance of class Control 442 Contains all the parameters of control files, which are: 443 DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM, 444 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, 445 LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, 446 ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, 447 ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL, 448 GRIB2FLEXPART, FLEXPARTDIR 449 For more information about format and content of the parameter 450 see documentation. 451 452 message: string 453 Message for exiting program. Default value is "Done!". 454 455 @Return: 456 <nothing> 457 458 ''' 459 # Uncomment if user wants notification directly from python 460 #try: 461 #target = c.mailops 462 #if(type(target) is not list): 463 #target = [target] 464 #for t in target: 465 #p = subprocess.Popen(['mail','-s ECMWFDATA v7.0 normal exit', 466 # os.path.expandvars(t)], 467 # stdin = subprocess.PIPE, 468 # stdout = subprocess.PIPE, 469 # stderr = subprocess.PIPE, bufsize = 1) 470 #pout = p.communicate(input = message+'\n\n')[0] 471 #print pout.decode() 472 #except: 473 #pass 474 475 print(message) 476 477 return 478 479 309 480 def product(*args, **kwds): 481 ''' 482 @Description: 483 484 @Input: 485 486 @Return: 487 488 ''' 310 489 # product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy 311 # product(range(2), repeat =3) --> 000 001 010 011 100 101 110 111490 # product(range(2), repeat = 3) --> 000 001 010 011 100 101 110 111 312 491 pools = map(tuple, args) * kwds.get('repeat', 1) 313 492 result = [[]] 314 493 for pool in pools: 315 result = [x +[y] for x in result for y in pool]494 result = [x + [y] for x in result for y in pool] 316 495 for prod in result: 317 496 yield tuple(prod) 318 497 319 ############################################################### 320 # utility to remove a file if it exists 321 # it does not fail if the file does not exist 322 ############################################################### 498 return 499 500 323 501 def silentremove(filename): 502 ''' 503 @Description: 504 Removes the file which name is passed to the function if 505 it exists. The function does not fail if the file does not 506 exist. 507 508 @Input: 509 filename: string 510 The name of the file to be removed without notification. 511 512 @Return: 513 <nothing> 514 ''' 324 515 try: 325 516 os.remove(filename) 326 except OSError as e: # this would be "except OSError, e:" before Python 2.6 327 if e.errno != errno.ENOENT: # errno.ENOENT = no such file or directory 328 raise # re-raise exception if a different error occured 517 except OSError as e: 518 # this would be "except OSError, e:" before Python 2.6 519 if e.errno is not errno.ENOENT: 520 # errno.ENOENT = no such file or directory 521 raise # re-raise exception if a different error occured 522 523 return 329 524 330 525 331 526 def init128(fn): 332 333 table128=dict() 527 ''' 528 @Description: 529 Opens and reads the grib table 128 file. 530 531 @Input: 532 fn: string 533 Path to file of ECMWF grib table number 128. 534 535 @Return: 536 table128: dictionary 537 Contains the ECMWF grib table 128 information. 538 ''' 539 table128 = dict() 334 540 with open(fn) as f: 335 336 337 if data[0]!='!':338 table128[data[0:3]]=data[59:64].strip()339 541 fdata = f.read().split('\n') 542 for data in fdata: 543 if data[0] != '!': 544 table128[data[0:3]] = data[59:64].strip() 545 340 546 return table128 341 547 342 def toparamId(pars,table): 343 344 cpar=pars.upper().split('/') 345 ipar=[] 548 549 def toparamId(pars, table): 550 ''' 551 @Description: 552 Transform parameter names to parameter ids 553 with ECMWF grib table 128. 554 555 @Input: 556 pars: string 557 Addpar argument from control file in the format of 558 parameter names instead of ids. 559 560 table: dictionary 561 Contains the ECMWF grib table 128 information. 562 563 @Return: 564 ipar: list of integer 565 List of addpar parameters from control file transformed to 566 parameter ids in the format of integer. 567 ''' 568 cpar = pars.upper().split('/') 569 ipar = [] 346 570 for par in cpar: 347 found=False 348 for k,v in table.iteritems(): 349 if par==k or par==v: 350 ipar.append(int(k)) 351 found=True 352 break 353 if found==False: 354 print 'Warning: par '+par+' not found in table 128' 571 found = False 572 for k, v in table.iteritems(): 573 if par == k or par == v: 574 ipar.append(int(k)) 575 found = True 576 break 577 if found is False: 578 print('Warning: par ' + par + ' not found in table 128') 579 355 580 return ipar 356 581 582 357 583 def dapoly(alist): 358 359 360 pya=(alist[3]-alist[0]+3.*(alist[1]-alist[2]))/6. 361 pyb=(alist[2]+alist[0])/2.-alist[1]-9.*pya/2. 362 pyc=alist[1]-alist[0]-7.*pya/2.-2.*pyb 363 pyd=alist[0]-pya/4.-pyb/3.-pyc/2. 364 sfeld=8.*pya+4.*pyb+2.*pyc+pyd 365 584 ''' 585 @Description: 586 587 @Input: 588 589 @Return: 590 591 ''' 592 pya = (alist[3] - alist[0] + 3. * (alist[1] - alist[2])) / 6. 593 pyb = (alist[2] + alist[0]) / 2. - alist[1] - 9. * pya / 2. 594 pyc = alist[1] - alist[0] - 7. * pya / 2. - 2. * pyb 595 pyd = alist[0] - pya / 4. - pyb / 3. - pyc / 2. 596 sfeld = 8. * pya + 4. * pyb + 2. * pyc + pyd 597 366 598 return sfeld 367 599 600 368 601 def darain(alist): 369 xa=alist[0] 370 xb=alist[1] 371 xc=alist[2] 372 xd=alist[3] 373 xa[xa<0.]=0. 374 xb[xb<0.]=0. 375 xc[xc<0.]=0. 376 xd[xd<0.]=0. 377 378 xac=0.5*xb 379 mask=xa+xc>0. 380 xac[mask]=xb[mask]*xc[mask]/(xa[mask]+xc[mask]) 381 xbd=0.5*xc 382 mask=xb+xd>0. 383 xbd[mask]=xb[mask]*xc[mask]/(xb[mask]+xd[mask]) 384 sfeld=xac+xbd 385 602 ''' 603 @Description: 604 Disaggregate a list of 4 precipitation values to generate a new value 605 for the second position of the 4 value list. This is used for 606 precipitation fields. 607 608 @Input: 609 alist: list of size 4, float 610 List of 4 precipitation values. 611 612 @Return: 613 sfeld: float 614 New value for the second position of the disaggregated 615 precipitation field. 616 ''' 617 xa = alist[0] 618 xb = alist[1] 619 xc = alist[2] 620 xd = alist[3] 621 xa[xa < 0.] = 0. 622 xb[xb < 0.] = 0. 623 xc[xc < 0.] = 0. 624 xd[xd < 0.] = 0. 625 626 xac = 0.5 * xb 627 mask = xa + xc > 0. 628 xac[mask] = xb[mask] * xc[mask] / (xa[mask] + xc[mask]) 629 xbd = 0.5 * xc 630 mask = xb + xd > 0. 631 xbd[mask] = xb[mask] * xc[mask] / (xb[mask] + xd[mask]) 632 sfeld = xac + xbd 633 386 634 return sfeld 387 635 636 388 637 class Control: 389 'class containing the information of the ECMWFDATA control file' 390 391 def __init__(self,filename): 392 393 #try: 394 with open(filename) as f: 395 fdata = f.read().split('\n') 396 for ldata in fdata: 397 data=ldata.split() 398 if len(data)>1: 399 if 'm_' in data[0].lower(): 400 data[0]=data[0][2:] 401 if data[0].lower()=='class': 402 data[0]='marsclass' 403 if data[0].lower()=='day1': 404 data[0]='start_date' 405 if data[0].lower()=='day2': 406 data[0]='end_date' 407 if data[0].lower()=='addpar': 408 if '/' in data[1]: 409 if data[1][0]=='/': 410 data[1]=data[1][1:] 411 dd=data[1].split('/') 412 data=[data[0]] 413 for d in dd: 414 data.append(d) 415 pass 416 if len(data)==2: 417 if '$' in data[1]: 418 # print data[1] 419 setattr(self,data[0].lower(),data[1]) 420 while '$' in data[1]: 421 i=data[1].index('$') 422 j=data[1].find('{') 423 k=data[1].find('}') 424 var=os.getenv(data[1][j+1:k]) 425 if var is not None: 426 data[1]=data[1][:i]+var+data[1][k+1:] 427 else: 428 myerror(None,'Could not find variable '+data[1][j+1:k]+' while reading '+filename) 429 430 # print data[0],data[1] 431 setattr(self,data[0].lower()+'_expanded',data[1]) 432 else: 433 if data[1].lower()!='none': 434 setattr(self,data[0].lower(),data[1]) 435 else: 436 setattr(self,data[0].lower(),None) 437 438 elif len(data)>2: 439 setattr(self,data[0].lower(),(data[1:])) 440 else: 441 pass 442 443 if not hasattr(self,'start_date'): 444 self.start_date=None 445 if not hasattr(self,'end_date'): 446 self.end_date=self.start_date 447 if not hasattr(self,'accuracy'): 448 self.accuracy=24 449 if not hasattr(self,'omega'): 450 self.omega='0' 451 if not hasattr(self,'cwc'): 452 self.cwc='0' 453 if not hasattr(self,'omegadiff'): 454 self.omegadiff='0' 455 if not hasattr(self,'etadiff'): 456 self.etadiff='0' 457 if not hasattr(self,'levelist'): 458 if not hasattr(self,'level'): 459 print 'Warning: neither levelist nor level specified in CONTROL file' 460 else: 461 self.levelist='1/to/'+self.level 462 else: 463 if 'to' in self.levelist: 464 self.level=self.levelist.split('/')[2] 465 else: 466 self.level=self.levelist.split('/')[-1] 467 468 if not hasattr(self,'maxstep'): 469 self.maxstep=0 470 for s in self.step: 471 if int(s)>self.maxstep: 472 self.maxstep=int(s) 473 else: 474 self.maxstep=int(self.maxstep) 475 if not hasattr(self,'prefix'): 476 self.prefix='EN' 477 if not hasattr(self,'makefile'): 478 self.makefile=None 479 if not hasattr(self,'basetime'): 480 self.basetime=None 481 if not hasattr(self,'date_chunk'): 482 self.date_chunk='3' 483 484 self.ecmwfdatadir=os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))+'/../' # script directory 485 # if not hasattr(self,'exedir'): 486 self.exedir=self.ecmwfdatadir+'src/' 487 if not hasattr(self,'grib2flexpart'): 488 self.grib2flexpart='0' 489 if not hasattr(self,'flexpart_root_scripts'): 490 self.flexpart_root_scripts=self.ecmwfdatadir 491 492 493 494 return 638 ''' 639 Class containing the information of the ECMWFDATA control file 640 641 Contains all the parameters of control files, which are: 642 DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM, 643 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, 644 LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, 645 ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, 646 ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL, 647 GRIB2FLEXPART, FLEXPARTDIR 648 For more information about format and content of the parameter 649 see documentation. 650 ''' 651 652 def __init__(self, filename): 653 ''' 654 @Description: 655 Initialises the instance of Control class and defines and 656 assign all controlfile variables. 657 658 @Input: 659 self: instance of Control class 660 Description see class documentation. 661 662 filename: string 663 Name of control file. 664 665 @Return: 666 <nothing> 667 ''' 668 with open(filename) as f: 669 fdata = f.read().split('\n') 670 for ldata in fdata: 671 data = ldata.split() 672 if len(data) > 1: 673 if 'm_' in data[0].lower(): 674 data[0] = data[0][2:] 675 if data[0].lower() == 'class': 676 data[0] = 'marsclass' 677 if data[0].lower() == 'day1': 678 data[0] = 'start_date' 679 if data[0].lower() == 'day2': 680 data[0] = 'end_date' 681 if data[0].lower() == 'addpar': 682 if '/' in data[1]: 683 # remove leading '/' sign from addpar content 684 if data[1][0] == '/': 685 data[1] = data[1][1:] 686 dd = data[1].split('/') 687 data = [data[0]] 688 for d in dd: 689 data.append(d) 690 pass 691 if len(data) == 2: 692 if '$' in data[1]: 693 setattr(self, data[0].lower(), data[1]) 694 while '$' in data[1]: 695 i = data[1].index('$') 696 j = data[1].find('{') 697 k = data[1].find('}') 698 var = os.getenv(data[1][j+1:k]) 699 if var is not None: 700 data[1] = data[1][:i] + var + data[1][k+1:] 701 else: 702 myerror(None, 'Could not find variable ' + 703 data[1][j+1:k] + ' while reading ' + 704 filename) 705 setattr(self, data[0].lower()+'_expanded', data[1]) 706 else: 707 if data[1].lower() != 'none': 708 setattr(self, data[0].lower(), data[1]) 709 else: 710 setattr(self, data[0].lower(), None) 711 elif len(data) > 2: 712 setattr(self, data[0].lower(), (data[1:])) 713 else: 714 pass 715 716 if not hasattr(self, 'start_date'): 717 self.start_date = None 718 if not hasattr(self, 'end_date'): 719 self.end_date = self.start_date 720 if not hasattr(self, 'accuracy'): 721 self.accuracy = 24 722 if not hasattr(self, 'omega'): 723 self.omega = '0' 724 if not hasattr(self, 'cwc'): 725 self.cwc = '0' 726 if not hasattr(self, 'omegadiff'): 727 self.omegadiff = '0' 728 if not hasattr(self, 'etadiff'): 729 self.etadiff = '0' 730 if not hasattr(self, 'levelist'): 731 if not hasattr(self, 'level'): 732 print(('Warning: neither levelist nor level \ 733 specified in CONTROL file')) 734 else: 735 self.levelist = '1/to/' + self.level 736 else: 737 if 'to' in self.levelist: 738 self.level = self.levelist.split('/')[2] 739 else: 740 self.level = self.levelist.split('/')[-1] 741 742 if not hasattr(self, 'maxstep'): 743 self.maxstep = 0 744 for s in self.step: 745 if int(s) > self.maxstep: 746 self.maxstep = int(s) 747 else: 748 self.maxstep = int(self.maxstep) 749 750 if not hasattr(self, 'prefix'): 751 self.prefix = 'EN' 752 if not hasattr(self, 'makefile'): 753 self.makefile = None 754 if not hasattr(self, 'basetime'): 755 self.basetime = None 756 if not hasattr(self, 'date_chunk'): 757 self.date_chunk = '3' 758 759 self.ecmwfdatadir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))+'/../' # script directory 760 #AP wieso nicht abfragen? Wieso immer setzen? 761 # if not hasattr(self,'exedir'): 762 self.exedir = self.ecmwfdatadir + 'src/' 763 if not hasattr(self, 'grib2flexpart'): 764 self.grib2flexpart = '0' 765 if not hasattr(self, 'flexpart_root_scripts'): 766 self.flexpart_root_scripts = self.ecmwfdatadir 767 768 return 769 495 770 def __str__(self): 496 497 attrs = vars(self) 498 # {'kids': 0, 'name': 'Dog', 'color': 'Spotted', 'age': 10, 'legs': 2, 'smell': 'Alot'} 499 # now dump this in some way or another 500 return ', '.join("%s: %s" % item for item in attrs.items()) 501 771 ''' 772 @Description: 773 Prepares a single string with all the comma seperated Control 774 class attributes including their values. The string has the format 775 #AP ????????????????? anschaun 776 777 @Input: 778 self: instance of Control class 779 Description see class documentation. 780 781 @Return: 782 string of Control class attributes with their values 783 ''' 784 attrs = vars(self) 785 # {'kids': 0, 'name': 'Dog', 'color': 'Spotted', 'age': 10, 'legs': 2, 'smell': 'Alot'} 786 # now dump this in some way or another 787 return ', '.join("%s: %s" % item for item in attrs.items()) 788 502 789 def tolist(self): 503 attrs=vars(self) 504 l=list() 505 for item in attrs.items(): 506 if '_expanded' in item[0]: 507 pass 508 elif 'exedir' in item[0]: 509 pass 510 elif 'flexpart_root_scripts' in item[0]: 511 pass 512 elif 'ecmwfdatadir' in item[0]: 513 pass 514 515 else: 516 if type(item[1]) is list: 517 stot='' 518 for s in item[1]: 519 stot+=s+' ' 520 521 l.append("%s %s" % (item[0],stot)) 522 else: 523 l.append("%s %s" % item ) 524 return sorted(l) 525 526 527 528 ############################################################## 529 # MARSretrieval class 530 ############################################################## 790 ''' 791 @Description: 792 Just generates a list of strings containing the attributes and 793 assigned values except the attributes "_expanded", "exedir", 794 "ecmwfdatadir" and "flexpart_root_scripts". 795 796 @Input: 797 self: instance of Control class 798 Description see class documentation. 799 800 @Return: 801 l: list 802 A sorted list of the all Control class attributes with 803 their values except the attributes "_expanded", "exedir", 804 "ecmwfdatadir" and "flexpart_root_scripts". 805 ''' 806 attrs = vars(self) 807 l = list() 808 for item in attrs.items(): 809 if '_expanded' in item[0]: 810 pass 811 elif 'exedir' in item[0]: 812 pass 813 elif 'flexpart_root_scripts' in item[0]: 814 pass 815 elif 'ecmwfdatadir' in item[0]: 816 pass 817 else: 818 if type(item[1]) is list: 819 stot = '' 820 for s in item[1]: 821 stot += s + ' ' 822 823 l.append("%s %s" % (item[0], stot)) 824 else: 825 #AP syntax error with doubled %s ??? 826 l.append("%s %s" % item ) 827 return sorted(l) 828 829 531 830 class MARSretrieval: 532 831 'class for MARS retrievals' 533 832 534 def __init__(self,server, marsclass="ei",type="",levtype="",levelist="", 535 repres="", date="",resol="",stream="",area="",time="",step="",expver="1", 536 number="",accuracy="", grid="", gaussian="", target="",param="",dataset="",expect=""): 537 # self.dataset=dataset 538 self.marsclass=marsclass 539 self.type=type 540 self.levtype=levtype 541 self.levelist=levelist 542 self.repres=repres 543 self.date=date 544 self.resol=resol 545 self.stream=stream 546 self.area=area 547 self.time=time 548 self.step=step 549 self.expver=expver 550 self.target=target 551 self.param=param 552 self.number=number 553 self.accuracy=accuracy 554 self.grid=grid 555 self.gaussian=gaussian 556 # self.expect=expect 557 self.server=server 833 def __init__(self, server, marsclass = "ei", type = "", levtype = "", 834 levelist = "", repres = "", date = "", resol = "", stream = "", 835 area = "", time = "", step = "", expver = "1", number = "", 836 accuracy = "", grid = "", gaussian = "", target = "", 837 param = "", dataset = "", expect = ""): 838 ''' 839 @Description: 840 841 @Input: 842 self 843 server 844 marsclass = "ei" 845 type = "" 846 levtype = "" 847 levelist = "" 848 repres = "" 849 date = "" 850 resol = "" 851 stream = "" 852 area = "" 853 time = "" 854 step = "" 855 expver = "1" 856 number = "" 857 accuracy = "" 858 grid = "" 859 gaussian = "" 860 target = "" 861 param = "" 862 dataset = "" 863 expect = "" 864 865 @Return: 866 <nothing> 867 ''' 868 869 # self.dataset = dataset 870 self.marsclass = marsclass 871 self.type = type 872 self.levtype = levtype 873 self.levelist = levelist 874 self.repres = repres 875 self.date = date 876 self.resol = resol 877 self.stream = stream 878 self.area = area 879 self.time = time 880 self.step = step 881 self.expver = expver 882 self.target = target 883 self.param = param 884 self.number = number 885 self.accuracy = accuracy 886 self.grid = grid 887 self.gaussian = gaussian 888 # self.expect = expect 889 self.server = server 890 891 return 558 892 559 893 def displayInfo(self): 560 attrs=vars(self) 561 for item in attrs.items(): 562 if item[0] in ('server'): 563 pass 564 else: 565 print item[0]+': '+str(item[1]) 566 567 return 894 ''' 895 @Description: 896 . 897 898 @Input: 899 self: 900 901 @Return: 902 <nothing> 903 ''' 904 attrs = vars(self) 905 for item in attrs.items(): 906 if item[0] in ('server'): 907 pass 908 else: 909 print(item[0] + ': ' + str(item[1])) 910 911 return 568 912 569 913 def dataRetrieve(self): 570 attrs=vars(self) 571 572 # self.server.retrieve(dicolist) 573 s='ret' 574 for k,v in attrs.iteritems(): 575 if k in ('server'): 576 continue 577 if k=='marsclass': 578 k='class' 579 if v=='': 580 continue 581 if k.lower()=='target': 582 target=v 583 else: 584 s=s+','+k+'='+str(v) 585 586 if self.server != False: 587 try: 588 self.server.execute(s,target) 589 except: 590 print 'MARS Request failed, have you already registered at apps.ecmwf.int?' 591 raise IOError 592 if os.stat(target).st_size==0: 593 print 'MARS Request returned no data - please check request' 594 raise IOError 914 ''' 915 @Description: 916 . 917 918 @Input: 919 self: 920 921 @Return: 922 <nothing> 923 ''' 924 attrs = vars(self) 925 # self.server.retrieve(dicolist) 926 s = 'ret' 927 for k, v in attrs.iteritems(): 928 if k in ('server'): 929 continue 930 if k == 'marsclass': 931 k = 'class' 932 if v == '': 933 continue 934 if k.lower() == 'target': 935 target = v 936 else: 937 s = s + ',' + k + '=' + str(v) 938 939 if self.server != False: 940 try: 941 self.server.execute(s,target) 942 except: 943 print('MARS Request failed, have you already registered at apps.ecmwf.int?') 944 raise IOError 945 if os.stat(target).st_size == 0: 946 print('MARS Request returned no data - please check request') 947 raise IOError 595 948 else: 596 597 s+=',target="'+target+'"' 598 p=subprocess.Popen(['mars'], stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE,bufsize=1) 599 pout=p.communicate(input=s)[0] 600 print pout.decode() 601 if 'Some errors reported' in pout.decode(): 602 print 'MARS Request failed - please check request' 603 raise IOError 604 605 if os.stat(target).st_size==0: 606 print 'MARS Request returned no data - please check request' 607 raise IOError 608 609 return 610 949 s += ',target = "' + target + '"' 950 p = subprocess.Popen(['mars'], stdin=subprocess.PIPE, 951 stdout=subprocess.PIPE, 952 stderr=subprocess.PIPE, bufsize = 1) 953 pout = p.communicate(input=s)[0] 954 print(pout.decode()) 955 if 'Some errors reported' in pout.decode(): 956 print('MARS Request failed - please check request') 957 raise IOError 958 959 if os.stat(target).st_size == 0: 960 print('MARS Request returned no data - please check request') 961 raise IOError 962 963 return 964 965 966 ############################################################## 611 967 ############################################################## 612 968 class EIFlexpart: 613 'class to retrieve Era Interim data for running FLEXPART' 614 ############################################################## 615 616 def __init__(self,c,fluxes=False): 617 # different mars types for retrieving reanalysis data for flexpart 618 619 self.types=dict() 620 try: 621 if c.maxstep>len(c.type): # Pure forecast mode 622 c.type=[c.type[1]] 623 c.step=['{:0>3}'.format(int(c.step[0]))] 624 c.time=[c.time[0]] 625 for i in range(1,c.maxstep+1): 626 c.type.append(c.type[0]) 627 c.step.append('{:0>3}'.format(i)) 628 c.time.append(c.time[0]) 629 except: 630 pass 631 632 self.inputdir=c.inputdir 633 self.basetime=c.basetime 634 self.dtime=c.dtime 635 self.mars={} 636 i=0 637 j=0 638 if fluxes==True and c.maxstep<24: 639 self.types[c.type[1]]={'times':'00/12','steps':'{}/to/12/by/{}'.format(c.dtime,c.dtime)} 640 i=1 641 for k in [0,12]: 642 for j in range(int(c.dtime),13,int(c.dtime)): 643 self.mars['{:0>3}'.format(i*int(c.dtime))]=[c.type[1],'{:0>2}'.format(k),'{:0>3}'.format(j)] 644 i+=1 645 else: 646 for ty,st,ti in zip(c.type,c.step,c.time): 647 btlist=range(24) 648 if c.basetime=='12': 649 btlist=[1,2,3,4,5,6,7,8,9,10,11,12] 650 # btday =[0,0,0,0,0,0,0,0,0, 0, 0, 0] 651 if c.basetime=='00': 652 btlist=[13,14,15,16,17,18,19,20,21,22,23,0] 653 # btday =[-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1] 654 # if c.step[0]=='00': 655 # btday[-1]=0 656 657 if mod(i,int(c.dtime))==0 and (i in btlist or c.maxstep>24): 658 if ty not in self.types.keys(): 659 self.types[ty]={'times':'','steps':''} 660 if ti not in self.types[ty]['times']: 661 if len(self.types[ty]['times'])>0: 662 self.types[ty]['times']+='/' 663 self.types[ty]['times']+=ti 664 if st not in self.types[ty]['steps']: 665 if len(self.types[ty]['steps'])>0: 666 self.types[ty]['steps']+='/' 667 self.types[ty]['steps']+=st 668 669 self.mars['{:0>3}'.format(j)]=[ty,'{:0>2}'.format(int(ti)),'{:0>3}'.format(int(st))] 670 j+=int(c.dtime) 671 672 i+=1 673 674 # Different grids need different retrievals 675 # SH=Spherical Harmonics, GG=Gaussian Grid, OG=Output Grid, ML=MultiLevel, SL=SingleLevel 676 self.params={'SH__ML':'','SH__SL':'', 677 'GG__ML':'','GG__SL':'', 678 'OG__ML':'','OG__SL':'', 679 'OG_OROLSM_SL':'','OG_acc_SL':''} 680 681 self.marsclass=c.marsclass 682 self.stream=c.stream 683 self.number=c.number 684 self.resol=c.resol 685 if 'N' in c.grid: # Gaussian output grid 686 self.grid=c.grid 687 self.area='G' 688 else: 689 self.grid='{}/{}'.format(int(c.grid)/1000.,int(c.grid)/1000.) 690 self.area='{}/{}/{}/{}'.format(int(c.upper)/1000.,int(c.left)/1000.,int(c.lower)/1000.,int(c.right)/1000.) 691 692 self.accuracy=c.accuracy 693 self.level=c.level 694 try: 695 self.levelist=c.levelist 696 except: 697 self.levelist='1/to/'+c.level 698 self.glevelist='1/to/'+c.level 699 try: 700 self.gaussian=c.gaussian 701 except: 702 self.gaussian='' 703 try: 704 self.dataset=c.dataset 705 except: 706 self.dataset='' 707 try: 708 self.expver=c.expver 709 except: 710 self.expver='1' 711 try: 712 self.number=c.number 713 except: 714 self.number='0' 715 716 self.outputfilelist=[] 717 718 719 # Now comes the nasty part that deals with the different scenarios we have: 720 # 1) Calculation of etadot on 721 # a) Gaussian grid 722 # b) Output grid 723 # c) Output grid using parameter 77 retrieved from MARS 724 # 3) Calculation/Retrieval of omega 725 # 4) Download also data for WRF 726 727 if fluxes==False: 728 self.params['SH__SL']=['LNSP','ML','1','OFF'] 729 # self.params['OG__SL']=["SD/MSL/TCC/10U/10V/2T/2D/129/172",'SFC','1',self.grid] 730 self.params['OG__SL']=["141/151/164/165/166/167/168/129/172",'SFC','1',self.grid] 731 self.params['OG_OROLSM__SL']=["160/27/28/173",'SFC','1',self.grid] 732 733 if len(c.addpar)>0: 734 if c.addpar[0]=='/': 735 c.addpar=c.addpar[1:] 736 self.params['OG__SL'][0]+='/'+'/'.join(c.addpar) 737 self.params['OG__ML']=['T/Q','ML',self.levelist,self.grid] 738 739 if c.gauss=='0' and c.eta=='1': # the simplest case 740 self.params['OG__ML'][0]+='/U/V/77' 741 elif c.gauss=='0' and c.eta=='0': # this is not recommended (inaccurate) 742 self.params['OG__ML'][0]+='/U/V' 743 elif c.gauss=='1' and c.eta=='0': # this is needed for data before 2008, or for reanalysis data 744 self.params['GG__SL']=['Q','ML','1','{}'.format((int(self.resol)+1)/2)] 745 self.params['SH__ML']=['U/V/D','ML',self.glevelist,'OFF'] 746 else: 747 print 'Warning: This is a very costly parameter combination, use only for debugging!' 748 self.params['GG__SL']=['Q','ML','1','{}'.format((int(self.resol)+1)/2)] 749 self.params['GG__ML']=['U/V/D/77','ML',self.glevelist,'{}'.format((int(self.resol)+1)/2)] 750 751 if c.omega=='1': 752 self.params['OG__ML'][0]+='/W' 753 754 try: # add cloud water content if necessary 755 if c.cwc=='1': 756 self.params['OG__ML'][0]+='/CLWC/CIWC' 757 except: 758 pass 759 760 try: # add vorticity and geopotential height for WRF if necessary 761 if c.wrf=='1': 762 self.params['OG__ML'][0]+='/Z/VO' 763 if '/D' not in self.params['OG__ML'][0]: 764 self.params['OG__ML'][0]+='/D' 765 # wrf_sfc='sp/msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4'.upper() 766 # 134/151/235/167/165/166/168/129/172/034/031/141/139/170/183/236/039/040/041/042 767 wrf_sfc='134/235/167/165/166/168/129/172/34/31/141/139/170/183/236/39/40/41/42'.upper() 768 lwrt_sfc=wrf_sfc.split('/') 769 for par in lwrt_sfc: 770 if par not in self.params['OG__SL'][0]: 771 self.params['OG__SL'][0]+='/'+par 772 773 except: 774 pass 775 else: 776 self.params['OG_acc_SL']=["LSP/CP/SSHF/EWSS/NSSS/SSR",'SFC','1',self.grid] 777 778 779 return 780 # add additional WRF specific parameters here 781 782 def create_namelist(self,c,filename): 783 784 # if inputdir=="": 785 # self.inputdir='.' 786 # else: 787 # self.inputdir=inputdir 788 789 self.inputdir=c.inputdir 790 area=asarray(self.area.split('/')).astype(float) 791 grid=asarray(self.grid.split('/')).astype(float) 792 793 # if area[3]<0: 794 # area[3]+=360 795 if area[1]>area[3]: 796 area[1]-=360 797 zyk=abs((area[3]-area[1]-360.)+grid[1])<1.e-6 798 maxl=int((area[3]-area[1])/grid[1])+1 799 maxb=int((area[0]-area[2])/grid[0])+1 800 801 f=open(self.inputdir+'/'+filename,'w') 802 f.write('&NAMGEN\n') 803 f.write(',\n '.join(['maxl='+str(maxl),'maxb='+str(maxb), 804 'mlevel='+self.level,'mlevelist='+'"'+self.levelist+'"', 805 'mnauf='+self.resol,'metapar='+'77', 806 'rlo0='+str(area[1]),'rlo1='+str(area[3]),'rla0='+str(area[2]),'rla1='+str(area[0]), 807 'momega='+c.omega,'momegadiff='+c.omegadiff,'mgauss='+c.gauss, 808 'msmooth='+c.smooth,'meta='+c.eta,'metadiff='+c.etadiff,'mdpdeta='+c.dpdeta])) 809 810 f.write('\n/\n') 811 f.close() 812 return 813 814 815 def retrieve(self, server, dates,times, inputdir=''): 816 self.dates=dates 817 self.server=server 818 819 if inputdir=="": 820 self.inputdir='.' 821 else: 822 self.inputdir=inputdir 823 824 # Retrieve Q not for using Q but as a template for a reduced gaussian grid one date and time is enough 825 # Take analysis at 00 826 qdate=self.dates 827 idx=qdate.find("/") 828 if (idx >0): 829 qdate=self.dates[:idx] 830 831 #QG= MARSretrieval(self.server, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream, type="an", levtype="ML", levelist="1", 832 #gaussian="reduced",grid='{}'.format((int(self.resol)+1)/2), resol=self.resol,accuracy=self.accuracy,target=self.inputdir+"/"+"QG.grb", 833 #date=qdate, time="00",expver=self.expver, param="133.128") 969 ''' 970 Class to retrieve Era Interim data for running FLEXPART 971 ''' 972 #AP change class name to ECFlexpart 973 def __init__(self, c, fluxes=False): 974 ''' 975 @Description: 976 Creates an object/instance of EIFlexpart with the 977 associated settings of its attributes for the retrieval. 978 979 @Input: 980 self: instance of EIFlexpart 981 The current object of the class. 982 983 c: instance of class Control 984 Contains all the parameters of control files, which are: 985 DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM, 986 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, 987 LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, 988 ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, 989 ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL, 990 GRIB2FLEXPART, FLEXPARTDIR 991 For more information about format and content of the parameter 992 see documentation. 993 994 fluxes: boolean, optional 995 Decides if a the flux parameter settings are stored or 996 the rest of the parameter list. 997 Default value is False. 998 999 @Return: 1000 <nothing> 1001 ''' 1002 # different mars types for retrieving reanalysis data for flexpart 1003 1004 self.types = dict() 1005 try: 1006 if c.maxstep > len(c.type): # Pure forecast mode 1007 c.type = [c.type[1]] 1008 c.step = ['{:0>3}'.format(int(c.step[0]))] 1009 c.time = [c.time[0]] 1010 for i in range(1, c.maxstep + 1): 1011 c.type.append(c.type[0]) 1012 c.step.append('{:0>3}'.format(i)) 1013 c.time.append(c.time[0]) 1014 except: 1015 pass 1016 1017 self.inputdir = c.inputdir 1018 self.basetime = c.basetime 1019 self.dtime = c.dtime 1020 self.mars = {} 1021 i = 0 1022 j = 0 1023 if fluxes is True and c.maxstep < 24: 1024 # only FC with start times at 00/12 (without 00/12) 1025 self.types[c.type[1]] = {'times': '00/12', 1026 'steps': '{}/to/12/by/{}'.format(c.dtime, 1027 c.dtime)} 1028 i = 1 1029 for k in [0, 12]: 1030 for j in range(int(c.dtime), 13, int(c.dtime)): 1031 self.mars['{:0>3}'.format(i * int(c.dtime))] = \ 1032 [c.type[1], '{:0>2}'.format(k), '{:0>3}'.format(j)] 1033 i += 1 1034 else: 1035 for ty, st, ti in zip(c.type, c.step, c.time): 1036 btlist = range(24) 1037 if c.basetime == '12': 1038 btlist = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] 1039 if c.basetime == '00': 1040 btlist = [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0] 1041 1042 if mod(i, int(c.dtime)) == 0 and \ 1043 (i in btlist or c.maxstep > 24): 1044 1045 if ty not in self.types.keys(): 1046 self.types[ty] = {'times': '', 'steps': ''} 1047 1048 if ti not in self.types[ty]['times']: 1049 if len(self.types[ty]['times']) > 0: 1050 self.types[ty]['times'] += '/' 1051 self.types[ty]['times'] += ti 1052 1053 if st not in self.types[ty]['steps']: 1054 if len(self.types[ty]['steps']) > 0: 1055 self.types[ty]['steps'] += '/' 1056 self.types[ty]['steps'] += st 1057 1058 self.mars['{:0>3}'.format(j)] = [ty, 1059 '{:0>2}'.format(int(ti)), 1060 '{:0>3}'.format(int(st))] 1061 j += int(c.dtime) 1062 1063 i += 1 1064 1065 # Different grids need different retrievals 1066 # SH = Spherical Harmonics, GG = Gaussian Grid, 1067 # OG = Output Grid, ML = MultiLevel, SL = SingleLevel 1068 self.params = {'SH__ML': '', 'SH__SL': '', 1069 'GG__ML': '', 'GG__SL': '', 1070 'OG__ML': '', 'OG__SL': '', 1071 'OG_OROLSM_SL': '', 'OG_acc_SL': ''} 1072 self.marsclass = c.marsclass 1073 self.stream = c.stream 1074 self.number = c.number 1075 self.resol = c.resol 1076 if 'N' in c.grid: # Gaussian output grid 1077 self.grid = c.grid 1078 self.area = 'G' 1079 else: 1080 self.grid = '{}/{}'.format(int(c.grid) / 1000., int(c.grid) / 1000.) 1081 self.area = '{}/{}/{}/{}'.format(int(c.upper) / 1000., 1082 int(c.left) / 1000., 1083 int(c.lower) / 1000., 1084 int(c.right) / 1000.) 1085 1086 self.accuracy = c.accuracy 1087 self.level = c.level 1088 try: 1089 self.levelist = c.levelist 1090 except: 1091 self.levelist = '1/to/' + c.level 1092 1093 self.glevelist = '1/to/' + c.level 1094 1095 try: 1096 self.gaussian = c.gaussian 1097 except: 1098 self.gaussian = '' 1099 1100 try: 1101 self.dataset = c.dataset 1102 except: 1103 self.dataset = '' 1104 1105 try: 1106 self.expver = c.expver 1107 except: 1108 self.expver = '1' 1109 1110 try: 1111 self.number = c.number 1112 except: 1113 self.number = '0' 1114 1115 self.outputfilelist = [] 1116 1117 1118 # Now comes the nasty part that deals with the different scenarios we have: 1119 # 1) Calculation of etadot on 1120 # a) Gaussian grid 1121 # b) Output grid 1122 # c) Output grid using parameter 77 retrieved from MARS 1123 # 3) Calculation/Retrieval of omega 1124 # 4) Download also data for WRF 1125 1126 if fluxes is False: 1127 self.params['SH__SL'] = ['LNSP', 'ML', '1', 'OFF'] 1128 # self.params['OG__SL'] = ["SD/MSL/TCC/10U/10V/2T/2D/129/172", 1129 # 'SFC','1',self.grid] 1130 self.params['OG__SL'] = ["141/151/164/165/166/167/168/129/172", \ 1131 'SFC', '1', self.grid] 1132 self.params['OG_OROLSM__SL'] = ["160/27/28/173", \ 1133 'SFC', '1', self.grid] 1134 1135 if len(c.addpar) > 0: 1136 if c.addpar[0] == '/': 1137 c.addpar = c.addpar[1:] 1138 self.params['OG__SL'][0] += '/' + '/'.join(c.addpar) 1139 self.params['OG__ML'] = ['T/Q', 'ML', self.levelist, self.grid] 1140 1141 if c.gauss == '0' and c.eta == '1': 1142 # the simplest case 1143 self.params['OG__ML'][0] += '/U/V/77' 1144 elif c.gauss == '0' and c.eta == '0': 1145 # this is not recommended (inaccurate) 1146 self.params['OG__ML'][0] += '/U/V' 1147 elif c.gauss == '1' and c.eta == '0': 1148 # this is needed for data before 2008, or for reanalysis data 1149 self.params['GG__SL'] = ['Q', 'ML', '1', \ 1150 '{}'.format((int(self.resol) + 1) / 2)] 1151 self.params['SH__ML'] = ['U/V/D', 'ML', self.glevelist, 'OFF'] 1152 else: 1153 print('Warning: This is a very costly parameter combination, \ 1154 use only for debugging!') 1155 self.params['GG__SL'] = ['Q', 'ML', '1', \ 1156 '{}'.format((int(self.resol) + 1) / 2)] 1157 self.params['GG__ML'] = ['U/V/D/77', 'ML', self.glevelist, \ 1158 '{}'.format((int(self.resol) + 1) / 2)] 1159 1160 if c.omega == '1': 1161 self.params['OG__ML'][0] += '/W' 1162 1163 try: 1164 # add cloud water content if necessary 1165 if c.cwc == '1': 1166 self.params['OG__ML'][0] += '/CLWC/CIWC' 1167 except: 1168 pass 1169 1170 try: 1171 # add vorticity and geopotential height for WRF if necessary 1172 if c.wrf == '1': 1173 self.params['OG__ML'][0] += '/Z/VO' 1174 if '/D' not in self.params['OG__ML'][0]: 1175 self.params['OG__ML'][0] += '/D' 1176 # wrf_sfc = 'sp/msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/ / 1177 # stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4'.upper() 1178 wrf_sfc = '134/235/167/165/166/168/129/172/34/31/141/ \ 1179 139/170/183/236/39/40/41/42'.upper() 1180 lwrt_sfc = wrf_sfc.split('/') 1181 for par in lwrt_sfc: 1182 if par not in self.params['OG__SL'][0]: 1183 self.params['OG__SL'][0] += '/' + par 1184 except: 1185 pass 1186 else: 1187 self.params['OG_acc_SL'] = ["LSP/CP/SSHF/EWSS/NSSS/SSR", \ 1188 'SFC', '1', self.grid] 1189 1190 # if needed, add additional WRF specific parameters here 1191 1192 return 1193 1194 1195 def create_namelist(self, c, filename): 1196 ''' 1197 @Description: 1198 Creates a namelist file in the temporary directory. 1199 It will contain values for maxl, maxb, mlevel, 1200 mlevelist, mnauf, metapar, rlo0, rlo1, rla0, rla1, 1201 momega, momegadiff, mgauss, msmooth, meta, metadiff, mdpdeta 1202 1203 @Input: 1204 self: instance of EIFlexpart 1205 The current object of the class. 1206 1207 c: instance of class Control 1208 Contains all the parameters of control files, which are: 1209 DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM, 1210 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, 1211 LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, 1212 ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, 1213 ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL, 1214 GRIB2FLEXPART, FLEXPARTDIR 1215 For more information about format and content of the parameter 1216 see documentation. 1217 1218 filename: string 1219 Name of the namelist file. 1220 1221 @Return: 1222 <nothing> 1223 ''' 1224 self.inputdir = c.inputdir 1225 area = asarray(self.area.split('/')).astype(float) 1226 grid = asarray(self.grid.split('/')).astype(float) 1227 1228 if area[1] > area[3]: 1229 area[1] -= 360 1230 zyk = abs((area[3] - area[1] - 360.) + grid[1]) < 1.e-6 1231 maxl = int((area[3] - area[1]) / grid[1]) + 1 1232 maxb = int((area[0] - area[2]) / grid[0]) + 1 1233 1234 f = open(self.inputdir + '/' + filename, 'w') 1235 f.write('&NAMGEN\n') 1236 f.write(',\n '.join(['maxl = ' + str(maxl), 'maxb = ' + str(maxb), 1237 'mlevel = ' + self.level, 1238 'mlevelist = ' + '"' + self.levelist + '"', 1239 'mnauf = ' + self.resol, 'metapar = ' + '77', 1240 'rlo0 = ' + str(area[1]), 'rlo1 = ' + str(area[3]), 1241 'rla0 = ' + str(area[2]), 'rla1 = ' + str(area[0]), 1242 'momega = ' + c.omega, 'momegadiff = ' + c.omegadiff, 1243 'mgauss = ' + c.gauss, 'msmooth = ' + c.smooth, 1244 'meta = ' + c.eta, 'metadiff = ' + c.etadiff, 1245 'mdpdeta = ' + c.dpdeta])) 1246 1247 f.write('\n/\n') 1248 f.close() 1249 return 1250 1251 def retrieve(self, server, dates, times, inputdir=''): 1252 ''' 1253 @Description: 1254 1255 1256 @Input: 1257 self: instance of EIFlexpart 1258 1259 server: instance of ECMWFService 1260 1261 dates: 1262 1263 times: 1264 1265 inputdir: string 1266 Default string is empty (''). 1267 1268 @Return: 1269 <nothing> 1270 ''' 1271 self.dates = dates 1272 self.server = server 1273 1274 if inputdir == "": 1275 self.inputdir = '.' 1276 else: 1277 self.inputdir = inputdir 1278 1279 # Retrieve Q not for using Q but as a template for a reduced gaussian 1280 # grid one date and time is enough 1281 # Take analysis at 00 1282 qdate = self.dates 1283 idx = qdate.find("/") 1284 if (idx > 0): 1285 qdate = self.dates[:idx] 1286 1287 #QG = MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = self.stream, type = "an", levtype = "ML", levelist = "1", 1288 #gaussian = "reduced",grid = '{}'.format((int(self.resol)+1)/2), resol = self.resol,accuracy = self.accuracy,target = self.inputdir+"/"+"QG.grb", 1289 #date = qdate, time = "00",expver = self.expver, param = "133.128") 834 1290 #QG.displayInfo() 835 1291 #QG.dataRetrieve() 836 1292 837 oro=False1293 oro = False 838 1294 for ftype in self.types: 839 for pk,pv in self.params.iteritems(): 840 if isinstance(pv,str): # or pk!='GG__SL' : 841 continue 842 mftype=''+ftype 843 mftime=self.types[ftype]['times'] 844 mfstep=self.types[ftype]['steps'] 845 mfdate=self.dates 846 mfstream=self.stream 847 mftarget=self.inputdir+"/"+ftype+pk+'.'+self.dates.split('/')[0]+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb" 848 if pk=='OG__SL': 849 pass 850 if pk=='OG_OROLSM__SL': 851 if oro==False: 852 mfstream='OPER' 853 mftype='AN' 854 mftime='00' 855 mfstep='000' 856 mfdate=self.dates.split('/')[0] 857 mftarget=self.inputdir+"/"+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb" 858 oro=True 859 else: 860 continue 861 862 if pk=='GG__SL' and pv[0]=='Q': 863 area="" 864 gaussian='reduced' 865 else: 866 area=self.area 867 gaussian=self.gaussian 868 869 if self.basetime==None: 870 MR= MARSretrieval(self.server, dataset=self.dataset, marsclass=self.marsclass, stream=mfstream, 871 type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian, 872 accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area, 873 date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0]) 874 875 MR.displayInfo() 876 MR.dataRetrieve() 877 # The whole else section is only necessary for operational scripts. It could be removed 878 else: # check if mars job requests fields beyond basetime. If yes eliminate those fields since they may not 879 # be accessible with user's credentials 880 sm1=-1 881 if 'by' in mfstep: 882 sm1=2 883 tm1=-1 884 if 'by' in mftime: 885 tm1=2 886 maxtime=datetime.datetime.strptime(mfdate.split('/')[-1]+mftime.split('/')[tm1],'%Y%m%d%H')+ \ 887 datetime.timedelta(hours=int(mfstep.split('/')[sm1])) 888 889 elimit=datetime.datetime.strptime(mfdate.split('/')[-1]+self.basetime,'%Y%m%d%H') 890 891 if self.basetime=='12': 892 if 'acc' in pk: 893 894 # Strategy: if maxtime-elimit>=24h reduce date by 1, 895 # if 12h<=maxtime-elimit<12h reduce time for last date 896 # if maxtime-elimit<12h reduce step for last time 897 # A split of the MARS job into 2 is likely necessary. 898 maxtime=elimit-datetime.timedelta(hours=24) 899 mfdate='/'.join(('/'.join(mfdate.split('/')[:-1]),datetime.datetime.strftime(maxtime,'%Y%m%d'))) 900 901 MR= MARSretrieval(self.server, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream, 902 type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian, 903 accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area, 904 date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0]) 905 906 MR.displayInfo() 907 MR.dataRetrieve() 908 909 maxtime=elimit-datetime.timedelta(hours=12) 910 mfdate=datetime.datetime.strftime(maxtime,'%Y%m%d') 911 mftime='00' 912 mftarget=self.inputdir+"/"+ftype+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb" 913 914 MR= MARSretrieval(self.server, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream, 915 type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian, 916 accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area, 917 date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0]) 918 919 MR.displayInfo() 920 MR.dataRetrieve() 921 else: 922 MR= MARSretrieval(self.server, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream, 923 type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian, 924 accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area, 925 date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0]) 926 927 MR.displayInfo() 928 MR.dataRetrieve() 929 else: 930 maxtime=elimit-datetime.timedelta(hours=24) 931 mfdate=datetime.datetime.strftime(maxtime,'%Y%m%d') 932 933 mftimesave=''.join(mftime) 934 935 if '/' in mftime: 936 times=mftime.split('/') 937 while int(times[0])+int(mfstep.split('/')[0])<=12 and pk!='OG_OROLSM__SL' and 'acc' not in pk: 938 times=times[1:] 939 if len(times)>1: 940 mftime='/'.join(times) 941 else: 942 mftime=times[0] 943 944 MR= MARSretrieval(self.server, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream, 945 type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian, 946 accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area, 947 date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0]) 948 949 MR.displayInfo() 950 MR.dataRetrieve() 951 952 if int(mftimesave.split('/')[0])==0 and int(mfstep.split('/')[0])==0 and pk!='OG_OROLSM__SL': 953 mfdate=datetime.datetime.strftime(elimit,'%Y%m%d') 954 mftime='00' 955 mfstep='000' 956 mftarget=self.inputdir+"/"+ftype+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb" 957 958 MR= MARSretrieval(self.server, dataset=self.dataset, marsclass=self.marsclass, stream=self.stream, 959 type=mftype, levtype=pv[1], levelist=pv[2],resol=self.resol, gaussian=gaussian, 960 accuracy=self.accuracy,grid=pv[3],target=mftarget,area=area, 961 date=mfdate, time=mftime,number=self.number,step=mfstep, expver=self.expver, param=pv[0]) 962 963 MR.displayInfo() 964 MR.dataRetrieve() 965 966 967 968 print "MARS retrieve done... " 969 970 def getFlexpartTime(self, type,step, time): 971 cstep='{:0>3}'.format(step) 972 ctime='{:0>2}'.format(int(time/100)) 973 1295 for pk, pv in self.params.iteritems(): 1296 if isinstance(pv, str): # or pk != 'GG__SL' : 1297 continue 1298 mftype = ''+ftype 1299 mftime = self.types[ftype]['times'] 1300 mfstep = self.types[ftype]['steps'] 1301 mfdate = self.dates 1302 mfstream = self.stream 1303 mftarget = self.inputdir+"/"+ftype+pk+'.'+self.dates.split('/')[0]+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb" 1304 if pk == 'OG__SL': 1305 pass 1306 if pk == 'OG_OROLSM__SL': 1307 if oro is False: 1308 mfstream = 'OPER' 1309 mftype = 'AN' 1310 mftime = '00' 1311 mfstep = '000' 1312 mfdate = self.dates.split('/')[0] 1313 mftarget = self.inputdir + "/" + pk + '.' + mfdate + \ 1314 '.' + str(os.getppid()) + '.' + \ 1315 str(os.getpid()) + ".grb" 1316 oro = True 1317 else: 1318 continue 1319 1320 if pk == 'GG__SL' and pv[0] == 'Q': 1321 area = "" 1322 gaussian = 'reduced' 1323 else: 1324 area = self.area 1325 gaussian = self.gaussian 1326 1327 if self.basetime is None: 1328 MR = MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = mfstream, 1329 type = mftype, levtype = pv[1], levelist = pv[2],resol = self.resol, gaussian = gaussian, 1330 accuracy = self.accuracy,grid = pv[3],target = mftarget,area = area, 1331 date = mfdate, time = mftime,number = self.number,step = mfstep, expver = self.expver, param = pv[0]) 1332 1333 MR.displayInfo() 1334 MR.dataRetrieve() 1335 # The whole else section is only necessary for operational scripts. 1336 # It could be removed 1337 else: 1338 # check if mars job requests fields beyond basetime. 1339 # If yes eliminate those fields since they may not 1340 # be accessible with user's credentials 1341 sm1 = -1 1342 if 'by' in mfstep: 1343 sm1 = 2 1344 tm1 = -1 1345 if 'by' in mftime: 1346 tm1 = 2 1347 maxtime = datetime.datetime.strptime(mfdate.split('/')[-1]+mftime.split('/')[tm1],'%Y%m%d%H')+ \ 1348 datetime.timedelta(hours = int(mfstep.split('/')[sm1])) 1349 1350 elimit = datetime.datetime.strptime(mfdate.split('/')[-1]+self.basetime,'%Y%m%d%H') 1351 1352 if self.basetime == '12': 1353 if 'acc' in pk: 1354 1355 # Strategy: if maxtime-elimit> = 24h reduce date by 1, 1356 # if 12h< = maxtime-elimit<12h reduce time for last date 1357 # if maxtime-elimit<12h reduce step for last time 1358 # A split of the MARS job into 2 is likely necessary. 1359 maxtime = elimit-datetime.timedelta(hours = 24) 1360 mfdate = '/'.join(('/'.join(mfdate.split('/')[:-1]),datetime.datetime.strftime(maxtime,'%Y%m%d'))) 1361 1362 MR = MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = self.stream, 1363 type = mftype, levtype = pv[1], levelist = pv[2],resol = self.resol, gaussian = gaussian, 1364 accuracy = self.accuracy,grid = pv[3],target = mftarget,area = area, 1365 date = mfdate, time = mftime,number = self.number,step = mfstep, expver = self.expver, param = pv[0]) 1366 1367 MR.displayInfo() 1368 MR.dataRetrieve() 1369 1370 maxtime = elimit-datetime.timedelta(hours = 12) 1371 mfdate = datetime.datetime.strftime(maxtime,'%Y%m%d') 1372 mftime = '00' 1373 mftarget = self.inputdir+"/"+ftype+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb" 1374 1375 MR = MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = self.stream, 1376 type = mftype, levtype = pv[1], levelist = pv[2],resol = self.resol, gaussian = gaussian, 1377 accuracy = self.accuracy,grid = pv[3],target = mftarget,area = area, 1378 date = mfdate, time = mftime,number = self.number,step = mfstep, expver = self.expver, param = pv[0]) 1379 1380 MR.displayInfo() 1381 MR.dataRetrieve() 1382 else: 1383 MR = MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = self.stream, 1384 type = mftype, levtype = pv[1], levelist = pv[2],resol = self.resol, gaussian = gaussian, 1385 accuracy = self.accuracy,grid = pv[3],target = mftarget,area = area, 1386 date = mfdate, time = mftime,number = self.number,step = mfstep, expver = self.expver, param = pv[0]) 1387 1388 MR.displayInfo() 1389 MR.dataRetrieve() 1390 else: 1391 maxtime = elimit-datetime.timedelta(hours = 24) 1392 mfdate = datetime.datetime.strftime(maxtime,'%Y%m%d') 1393 1394 mftimesave = ''.join(mftime) 1395 1396 if '/' in mftime: 1397 times = mftime.split('/') 1398 while int(times[0])+int(mfstep.split('/')[0])< = 12 and pk! = 'OG_OROLSM__SL' and 'acc' not in pk: 1399 times = times[1:] 1400 if len(times)>1: 1401 mftime = '/'.join(times) 1402 else: 1403 mftime = times[0] 1404 1405 MR = MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = self.stream, 1406 type = mftype, levtype = pv[1], levelist = pv[2],resol = self.resol, gaussian = gaussian, 1407 accuracy = self.accuracy,grid = pv[3],target = mftarget,area = area, 1408 date = mfdate, time = mftime,number = self.number,step = mfstep, expver = self.expver, param = pv[0]) 1409 1410 MR.displayInfo() 1411 MR.dataRetrieve() 1412 1413 if int(mftimesave.split('/')[0]) == 0 and int(mfstep.split('/')[0]) == 0 and pk! = 'OG_OROLSM__SL': 1414 mfdate = datetime.datetime.strftime(elimit,'%Y%m%d') 1415 mftime = '00' 1416 mfstep = '000' 1417 mftarget = self.inputdir+"/"+ftype+pk+'.'+mfdate+'.'+str(os.getppid())+'.'+str(os.getpid())+".grb" 1418 1419 MR = MARSretrieval(self.server, dataset = self.dataset, marsclass = self.marsclass, stream = self.stream, 1420 type = mftype, levtype = pv[1], levelist = pv[2],resol = self.resol, gaussian = gaussian, 1421 accuracy = self.accuracy,grid = pv[3],target = mftarget,area = area, 1422 date = mfdate, time = mftime,number = self.number,step = mfstep, expver = self.expver, param = pv[0]) 1423 1424 MR.displayInfo() 1425 MR.dataRetrieve() 1426 1427 print("MARS retrieve done... ") 1428 1429 return 1430 1431 1432 def getFlexpartTime(self, type, step, time): 1433 #AP remove this function, no longer needed 1434 ''' 1435 @Description: 1436 ???? 1437 #AP wozu? wird das überhaupt noch gebraucht? in create auskommentiert 1438 1439 @Input: 1440 self: instance of EIFlexpart 1441 The current object of the class. 1442 1443 type: string 1444 Type of the data field. E.g. FC - forecast or 1445 AN - analysis 1446 1447 step: integer 1448 Forecast step in hours. 1449 1450 time: integer 1451 Time in hours. 1452 1453 @Return: 1454 cflextime: 1455 1456 ''' 1457 cstep = '{:0>3}'.format(step) 1458 ctime = '{:0>2}'.format(int(time / 100)) 974 1459 ctype = str(type).upper() 975 myinfo = [ctype,ctime, cstep] 1460 1461 myinfo = [ctype, ctime, cstep] 976 1462 cflextime = None 977 1463 for t, marsinfo in self.mars.items(): 978 1464 if myinfo == marsinfo: 979 cflextime=t 1465 cflextime = t 1466 980 1467 return cflextime 981 1468 982 1469 def process_output(self, c): 983 984 print 'Postprocessing:\n Format: {}\n'.format(c.format) 985 if c.ecapi==False: 986 print 'ecstorage: {}\n ecfsdir: {}\n'.format(c.ecstorage,c.ecfsdir) 987 if not hasattr(c,'gateway'): 988 c.gateway=os.getenv('GATEWAY') 989 if not hasattr(c,'destination'): 990 c.destination=os.getenv('DESTINATION') 991 print 'ectrans: {}\n gateway: {}\n destination: {}\n '.format(c.ectrans, c.gateway,c.destination) 992 print 'Output filelist: \n',self.outputfilelist 993 994 if c.format.lower()=='grib2': 995 for ofile in self.outputfilelist: 996 p=subprocess.check_call(['grib_set','-s','edition=2,productDefinitionTemplateNumber=8',ofile,ofile+'_2']) 997 p=subprocess.check_call(['mv',ofile+'_2',ofile]) 998 999 if int(c.ectrans)==1 and c.ecapi==False: 1000 for ofile in self.outputfilelist: 1001 p=subprocess.check_call(['ectrans','-overwrite','-gateway',c.gateway, 1002 '-remote',c.destination,'-source',ofile]) 1003 print 'ectrans:',p 1004 if int(c.ecstorage)==1 and c.ecapi==False: 1005 for ofile in self.outputfilelist: 1006 p=subprocess.check_call(['ecp','-o',ofile,os.path.expandvars(c.ecfsdir)]) 1007 1008 #20131107 000000 EN13110700 ON DISC 1009 if c.outputdir!=c.inputdir: 1010 for ofile in self.outputfilelist: 1011 p=subprocess.check_call(['mv',ofile,c.outputdir]) 1012 1013 if c.grib2flexpart=='1': 1014 f=open(c.outputdir+'/'+'AVAILABLE','w') 1015 clist=[] 1016 for ofile in self.outputfilelist: # generate AVAILABLE file 1017 fname=ofile.split('/') 1018 if '.' in fname[-1]: 1019 l=fname[-1].split('.') 1020 timestamp=datetime.datetime.strptime(l[0][-6:]+l[1],'%y%m%d%H') 1021 timestamp+=datetime.timedelta(hours=int(l[2])) 1022 cdate=datetime.datetime.strftime(timestamp,'%Y%m%d') 1023 chms=datetime.datetime.strftime(timestamp,'%H%M%S') 1024 1025 else: 1026 cdate='20'+fname[-1][-8:-2] 1027 chms=fname[-1][-2:]+'0000' 1028 clist.append(cdate+' '+chms+' '*6+fname[-1]+' '*14+'ON DISC') 1029 clist.sort() 1030 f.write('\n'.join(clist)+'\n') 1031 f.close() 1032 1033 pwd=os.path.abspath(c.outputdir) 1034 f=open(pwd+'/pathnames','w') 1035 f.write(pwd+'/Options/\n') 1036 f.write(pwd+'/\n') 1037 f.write(pwd+'/\n') 1038 f.write(pwd+'/AVAILABLE\n') 1039 f.write('==================\n') 1040 f.close() 1041 1042 if not os.path.exists(pwd+'/Options'): 1043 os.makedirs(pwd+'/Options') 1044 f=open(os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))+'/../Options/COMMAND','r') 1045 lflist=f.read().split('\n') 1046 i=0 1047 for l in lflist: 1048 if 'LDIRECT' in l.upper(): 1049 break 1050 i+=1 1051 1052 clist.sort() 1053 lflist=lflist[:i+1]+[clist[0][:16],clist[-1][:16]]+lflist[i+3:] 1054 g=open(pwd+'/Options/COMMAND','w') 1055 g.write('\n'.join(lflist)+'\n') 1056 g.close() 1057 1058 os.chdir(c.outputdir) 1059 p=subprocess.check_call([os.path.expandvars(os.path.expanduser(c.flexpart_root_scripts))+'/../FLEXPART_PROGRAM/grib2flexpart', 1060 'useAvailable','.']) 1061 os.chdir(pwd) 1062 1063 def create(self, inputfiles, c,): 1064 1065 table128=init128(c.ecmwfdatadir+'/grib_templates/ecmwf_grib1_table_128') 1066 wrfpars=toparamId('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',table128) 1067 # index_keys=["date","time","stepRange"] 1068 index_keys=["date","time","step"] 1069 indexfile=c.inputdir+"/date_time_stepRange.idx" 1470 ''' 1471 @Description: 1472 1473 1474 @Input: 1475 self: instance of EIFlexpart 1476 The current object of the class. 1477 1478 c: instance of class Control 1479 Contains all the parameters of control files, which are: 1480 DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM, 1481 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, 1482 LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, 1483 ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, 1484 ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL, 1485 GRIB2FLEXPART, FLEXPARTDIR 1486 For more information about format and content of the parameter 1487 see documentation. 1488 1489 @Return: 1490 <nothing> 1491 1492 ''' 1493 1494 print 'Postprocessing:\n Format: {}\n'.format(c.format) 1495 if c.ecapi is False: 1496 print 'ecstorage: {}\n ecfsdir: {}\n'.format(c.ecstorage, c.ecfsdir) 1497 if not hasattr(c, 'gateway'): 1498 c.gateway = os.getenv('GATEWAY') 1499 if not hasattr(c, 'destination'): 1500 c.destination = os.getenv('DESTINATION') 1501 print 'ectrans: {}\n gateway: {}\n destination: {}\n '} 1502 .format(c.ectrans, c.gateway, c.destination) 1503 print 'Output filelist: \n', self.outputfilelist 1504 1505 if c.format.lower() == 'grib2': 1506 for ofile in self.outputfilelist: 1507 p = subprocess.check_call(['grib_set', '-s', 'edition = 2, \ 1508 productDefinitionTemplateNumber = 8', 1509 ofile, ofile + '_2']) 1510 p = subprocess.check_call(['mv', ofile + '_2', ofile]) 1511 1512 if int(c.ectrans) == 1 and c.ecapi is False: 1513 for ofile in self.outputfilelist: 1514 p = subprocess.check_call(['ectrans', '-overwrite', '-gateway', 1515 c.gateway, '-remote', c.destination, 1516 '-source', ofile]) 1517 print 'ectrans:',p 1518 if int(c.ecstorage) == 1 and c.ecapi is False: 1519 for ofile in self.outputfilelist: 1520 p = subprocess.check_call(['ecp', '-o', ofile, 1521 os.path.expandvars(c.ecfsdir)]) 1522 1523 # 20131107 000000 EN13110700 ON DISC 1524 if c.outputdir != c.inputdir: 1525 for ofile in self.outputfilelist: 1526 p = subprocess.check_call(['mv', ofile, c.outputdir]) 1527 1528 if c.grib2flexpart == '1': 1529 f = open(c.outputdir + '/' + 'AVAILABLE', 'w') 1530 clist = [] 1531 for ofile in self.outputfilelist: # generate AVAILABLE file 1532 fname = ofile.split('/') 1533 if '.' in fname[-1]: 1534 l = fname[-1].split('.') 1535 timestamp = datetime.datetime.strptime(l[0][-6:] + l[1], 1536 '%y%m%d%H') 1537 timestamp += datetime.timedelta(hours=int(l[2])) 1538 cdate = datetime.datetime.strftime(timestamp, '%Y%m%d') 1539 chms = datetime.datetime.strftime(timestamp, '%H%M%S') 1540 else: 1541 cdate = '20' + fname[-1][-8:-2] 1542 chms = fname[-1][-2:] + '0000' 1543 clist.append(cdate + ' ' + chms + ' '*6 + 1544 fname[-1] + ' '*14 + 'ON DISC') 1545 clist.sort() 1546 f.write('\n'.join(clist) + '\n') 1547 f.close() 1548 1549 pwd = os.path.abspath(c.outputdir) 1550 f = open(pwd + '/pathnames','w') 1551 f.write(pwd + '/Options/\n') 1552 f.write(pwd + '/\n') 1553 f.write(pwd + '/\n') 1554 f.write(pwd + '/AVAILABLE\n') 1555 f.write(' = == = == = == = == = == == = \n') 1556 f.close() 1557 1558 if not os.path.exists(pwd + '/Options'): 1559 os.makedirs(pwd+'/Options') 1560 f = open(os.path.expandvars( 1561 os.path.expanduser(c.flexpart_root_scripts)) + 1562 '/../Options/COMMAND', 'r') 1563 lflist = f.read().split('\n') 1564 i = 0 1565 for l in lflist: 1566 if 'LDIRECT' in l.upper(): 1567 break 1568 i += 1 1569 1570 clist.sort() 1571 lflist = lflist[:i+1] + \ 1572 [clist[0][:16], clist[-1][:16]] + \ 1573 lflist[i+3:] 1574 1575 with open(pwd + '/Options/COMMAND', 'w') as g: 1576 g.write('\n'.join(lflist) + '\n') 1577 1578 os.chdir(c.outputdir) 1579 p = subprocess.check_call([os.path.expandvars( 1580 os.path.expanduser(c.flexpart_root_scripts)) + 1581 '/../FLEXPART_PROGRAM/grib2flexpart', 1582 'useAvailable', '.']) 1583 os.chdir(pwd) 1584 1585 return 1586 1587 def create(self, inputfiles, c): 1588 ''' 1589 @Description: 1590 #AP WOZU und ist einrueckung richtig 1591 1592 @Input: 1593 self: instance of EIFlexpart 1594 The current object of the class. 1595 1596 inputfiles: list of strings 1597 A list of filenames. 1598 1599 c: instance of class Control 1600 Contains all the parameters of control files, which are: 1601 DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM, 1602 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, 1603 LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, 1604 ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, 1605 ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL, 1606 GRIB2FLEXPART, FLEXPARTDIR 1607 For more information about format and content of the parameter 1608 see documentation. 1609 1610 @Return: 1611 <nothing> 1612 ''' 1613 1614 table128 = init128(c.ecmwfdatadir + 1615 '/grib_templates/ecmwf_grib1_table_128') 1616 #AP wieso wrf 1617 wrfpars = toparamId('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/ \ 1618 stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4', 1619 table128) 1620 index_keys = ["date", "time", "step"] 1621 indexfile = c.inputdir + "/date_time_stepRange.idx" 1070 1622 silentremove(indexfile) 1071 grib=GribTools(inputfiles.files) 1072 iid=grib.index(index_keys=index_keys, index_file = indexfile) 1073 1074 print 'index done...' 1075 fdict={'10':None,'11':None,'12':None,'13':None,'16':None,'17':None,'19':None,'21':None,'22':None,'20':None} 1076 for f in fdict.keys(): 1077 silentremove(c.inputdir+"/fort."+f) 1078 1623 grib = GribTools(inputfiles.files) 1624 iid = grib.index(index_keys=index_keys, index_file=indexfile) 1625 print 'index done...' 1626 1627 fdict = {'10':None, '11':None, '12':None, '13':None, '16':None, 1628 '17':None, '19':None, '21':None, '22':None, '20':None} 1629 for f in fdict.keys(): 1630 silentremove(c.inputdir + "/fort." + f) 1631 1632 index_vals = [] 1633 for key in index_keys: 1634 key_vals = grib_index_get(iid, key) 1635 print key_vals 1636 index_vals.append(key_vals) 1637 1638 # creates index file 1639 for prod in product(*index_vals): 1640 for i in range(len(index_keys)): 1641 grib_index_select(iid, index_keys[i], prod[i]) 1642 1643 gid = grib_new_from_index(iid) 1644 # do convert2 program if gid at this time is not None, 1645 # therefore save in hid 1646 hid = gid 1647 cflextime = None 1648 for k, f in fdict.iteritems(): 1649 fdict[k] = open(c.inputdir + '/fort.' + k, 'w') 1650 if gid is not None: 1651 cdate = str(grib_get(gid, 'date')) 1652 time = grib_get(gid, 'time') 1653 type = grib_get(gid, 'type') 1654 step = grib_get(gid, 'step') 1655 # cflextime = self.getFlexpartTime(type,step, time) 1656 timestamp = datetime.datetime.strptime( 1657 cdate + '{:0>2}'.format(time/100), '%Y%m%d%H') 1658 timestamp += datetime.timedelta(hours=int(step)) 1659 # print gid,index_keys[i],prod[i],cdate,time,step,timestamp 1660 1661 cdateH = datetime.datetime.strftime(timestamp, '%Y%m%d%H') 1662 chms = datetime.datetime.strftime(timestamp, '%H%M%S') 1663 1664 if c.basetime is not None: 1665 slimit = datetime.datetime.strptime( 1666 c.start_date + '00', '%Y%m%d%H') 1667 bt = '23' 1668 if c.basetime == '00': 1669 bt = '00' 1670 slimit = datetime.datetime.strptime( 1671 c.end_date + bt, '%Y%m%d%H') - \ 1672 datetime.timedelta(hours=12-int(c.dtime)) 1673 if c.basetime == '12': 1674 bt = '12' 1675 slimit = datetime.datetime.strptime( 1676 c.end_date + bt, '%Y%m%d%H') - \ 1677 datetime.timedelta(hours=12-int(c.dtime)) 1678 1679 elimit = datetime.datetime.strptime(c.end_date + bt, '%Y%m%d%H') 1680 1681 if timestamp < slimit or timestamp > elimit: 1682 continue 1683 1684 try: 1685 if c.wrf == '1': 1686 if 'olddate' not in locals(): 1687 fwrf = open(c.outputdir + '/WRF' + cdate + 1688 '.{:0>2}'.format(time) + '.000.grb2', 'w') 1689 olddate = cdate[:] 1690 else: 1691 if cdate != olddate: 1692 fwrf = open(c.outputdir + '/WRF' + cdate + 1693 '.{:0>2}'.format(time) + '.000.grb2', 'w') 1694 olddate = cdate[:] 1695 except AttributeError: 1696 pass 1697 1698 # print 'cyear '+cyear+'/'+cmonth+'/'+'/EI'+cyear[2:4]+cmonth+cday+cflextime 1699 1700 savedfields = [] 1701 while 1: 1702 if gid is None: 1703 break 1704 paramId = grib_get(gid, 'paramId') 1705 gridtype = grib_get(gid, 'gridType') 1706 datatype = grib_get(gid, 'dataType') 1707 levtype = grib_get(gid, 'typeOfLevel') 1708 if paramId == 133 and gridtype == 'reduced_gg': 1709 # Relative humidity (Q.grb) is used as a template only 1710 # so we need the first we "meet" 1711 fout = open(c.inputdir + '/fort.18', 'w') 1712 grib_write(gid, fout) 1713 fout.close() 1714 elif paramId == 131 or paramId == 132: 1715 grib_write(gid, fdict['10']) 1716 elif paramId == 130: 1717 grib_write(gid, fdict['11']) 1718 elif paramId == 133 and gridtype != 'reduced_gg': 1719 grib_write(gid, fdict['17']) 1720 elif paramId == 152: 1721 grib_write(gid, fdict['12']) 1722 elif paramId == 155 and gridtype == 'sh': 1723 grib_write(gid, fdict['13']) 1724 elif paramId in [129, 138, 155] and levtype == 'hybrid' \ 1725 and c.wrf == '1': 1726 # print paramId,'not written' 1727 pass 1728 elif paramId == 246 or paramId == 247: 1729 # cloud liquid water and ice 1730 if paramId == 246: 1731 clwc = grib_get_values(gid) 1732 else: 1733 clwc += grib_get_values(gid) 1734 grib_set_values(gid, clwc) 1735 # grib_set(gid,'shortName','qc') 1736 grib_set(gid, 'paramId', 201031) 1737 grib_write(gid, fdict['22']) 1738 elif paramId == 135: 1739 grib_write(gid, fdict['19']) 1740 elif paramId == 77: 1741 grib_write(gid, fdict['21']) 1742 else: 1743 if paramId not in savedfields: 1744 grib_write(gid, fdict['16']) 1745 savedfields.append(paramId) 1746 else: 1747 print 'duplicate ' + str(paramId) + ' not written' 1748 1749 try: 1750 if c.wrf == '1': 1751 if levtype == 'hybrid': 1752 if paramId in [129, 130, 131, 132, 133, 138, 155]: 1753 grib_write(gid, fwrf) 1754 else: 1755 if paramId in wrfpars: 1756 grib_write(gid, fwrf) 1757 except AttributeError: 1758 pass 1759 1760 grib_release(gid) 1761 gid = grib_new_from_index(iid) 1762 1763 for f in fdict.values(): 1764 f.close() 1765 1766 # call for CONVERT2 1767 1768 if hid is not None: 1769 pwd = os.getcwd() 1770 os.chdir(c.inputdir) 1771 if os.stat('fort.21').st_size == 0 and int(c.eta) == 1: 1772 print 'Parameter 77 (etadot) is missing, most likely it is \ 1773 not available for this type or date/time\n' 1774 print 'Check parameters CLASS, TYPE, STREAM, START_DATE\n' 1775 myerror(c, 'fort.21 is empty while parameter eta is set \ 1776 to 1 in CONTROL file') 1777 1778 p = subprocess.check_call([os.path.expandvars( 1779 os.path.expanduser(c.exedir)) + '/CONVERT2'], shell=True) 1780 os.chdir(pwd) 1781 # create the corresponding output file fort.15 1782 # (generated by CONVERT2) 1783 # + fort.16 (paramId 167 and paramId 168) 1784 fnout = c.inputdir + '/' + c.prefix 1785 if c.maxstep > 12: 1786 suffix = cdate[2:8] + '.{:0>2}'.format(time/100) + \ 1787 '.{:0>3}'.format(step) 1788 else: 1789 suffix = cdateH[2:10] 1790 1791 fnout += suffix 1792 print "outputfile = " + fnout 1793 self.outputfilelist.append(fnout) # needed for final processing 1794 fout = open(fnout, 'wb') 1795 shutil.copyfileobj(open(c.inputdir + '/fort.15', 'rb'), fout) 1796 if c.cwc == '1': 1797 shutil.copyfileobj(open(c.inputdir + '/fort.22', 'rb'), fout) 1798 shutil.copyfileobj(open(c.inputdir + '/flux' + cdate[0:2] + 1799 suffix, 'rb'), fout) 1800 shutil.copyfileobj(open(c.inputdir + '/fort.16', 'rb'), fout) 1801 orolsm = glob.glob(c.inputdir + 1802 '/OG_OROLSM__SL.*.' + c.ppid + '*')[0] 1803 shutil.copyfileobj(open(orolsm, 'rb'), fout) 1804 fout.close() 1805 if c.omega == '1': 1806 fnout = c.outputdir + '/OMEGA' 1807 fout = open(fnout, 'wb') 1808 shutil.copyfileobj(open(c.inputdir + '/fort.25', 'rb'), fout) 1809 1810 try: 1811 if c.wrf == '1': 1812 fwrf.close() 1813 except: 1814 pass 1815 1816 grib_index_release(iid) 1817 1818 return 1819 1820 1821 def deacc_fluxes(self, inputfiles, c): 1822 ''' 1823 @Description: 1824 1825 1826 @Input: 1827 self: instance of EIFlexpart 1828 The current object of the class. 1829 1830 inputfiles: list of strings 1831 A list of filenames. 1832 1833 c: instance of class Control 1834 Contains all the parameters of control files, which are: 1835 DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM, 1836 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, 1837 LEVELIST, RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, 1838 ETADIFF, DPDETA, SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, 1839 ECSTORAGE, ECTRANS, ECFSDIR, MAILOPS, MAILFAIL, 1840 GRIB2FLEXPART, FLEXPARTDIR 1841 For more information about format and content of the parameter 1842 see documentation. 1843 1844 @Return: 1845 <nothing> 1846 ''' 1847 1848 table128 = init128(c.ecmwfdatadir + 1849 '/grib_templates/ecmwf_grib1_table_128') 1850 pars = toparamId(self.params['OG_acc_SL'][0], table128) 1851 index_keys = ["date", "time", "step"] 1852 indexfile = c.inputdir + "/date_time_stepRange.idx" 1853 silentremove(indexfile) # delete, if it exists already 1854 grib = GribTools(inputfiles.files) 1855 iid = grib.index(index_keys=index_keys, index_file=indexfile) 1856 print 'index done...' 1857 1858 # get the date, time and step values from indexfile 1079 1859 index_vals = [] 1080 1860 for key in index_keys: 1081 1861 key_vals = grib_index_get(iid,key) 1082 print key_vals 1083 1862 print key_vals 1863 # have to sort the steps, therefore convert to int first 1864 if key == 'step': 1865 l = [] 1866 for k in key_vals: 1867 l.append(int(k)) 1868 l.sort() 1869 # now convert back to str and overwrite key_vals 1870 key_vals = [] 1871 for k in l: 1872 key_vals.append(str(k)) 1084 1873 index_vals.append(key_vals) 1085 1874 1086 for prod in product(*index_vals): 1087 for i in range(len(index_keys)): 1088 grib_index_select(iid,index_keys[i],prod[i]) 1089 1090 1091 gid = grib_new_from_index(iid) 1092 hid = gid 1093 cflextime = None 1094 for k,f in fdict.iteritems(): 1095 fdict[k] = open(c.inputdir+'/fort.'+k,'w') 1096 if gid is not None: 1097 cdate = str(grib_get(gid, 'date')) 1098 time = grib_get(gid, 'time') 1099 type = grib_get(gid, 'type') 1100 step = grib_get(gid, 'step') 1101 # step = grib_get(gid, 'stepRange') 1102 # cflextime = self.getFlexpartTime(type,step, time) 1103 timestamp=datetime.datetime.strptime(cdate+'{:0>2}'.format(time/100),'%Y%m%d%H') 1104 timestamp+=datetime.timedelta(hours=int(step)) 1105 # print gid,index_keys[i],prod[i],cdate,time,step,timestamp 1106 1107 cdateH=datetime.datetime.strftime(timestamp,'%Y%m%d%H') 1108 chms=datetime.datetime.strftime(timestamp,'%H%M%S') 1109 1110 if c.basetime !=None: 1111 slimit=datetime.datetime.strptime(c.start_date+'00','%Y%m%d%H') 1112 bt='23' 1113 if c.basetime=='00': 1114 bt='00' 1115 slimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H')-datetime.timedelta(hours=12-int(c.dtime)) 1116 1117 if c.basetime=='12': 1118 bt='12' 1119 slimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H')-datetime.timedelta(hours=12-int(c.dtime)) 1120 1121 elimit=datetime.datetime.strptime(c.end_date+bt,'%Y%m%d%H') 1122 1123 if timestamp<slimit or timestamp>elimit: 1124 continue 1125 1126 1127 1128 try: 1129 if c.wrf=='1': 1130 if 'olddate' not in locals(): 1131 fwrf = open(c.outputdir+'/WRF'+cdate+'.{:0>2}'.format(time)+'.000.grb2','w') 1132 olddate=cdate[:] 1133 else: 1134 if cdate!=olddate: 1135 fwrf = open(c.outputdir+'/WRF'+cdate+'.{:0>2}'.format(time)+'.000.grb2','w') 1136 olddate=cdate[:] 1137 except AttributeError: 1138 pass 1139 1140 # print 'cyear '+cyear+'/'+cmonth+'/'+'/EI'+cyear[2:4]+cmonth+cday+cflextime 1141 1142 savedfields=[] 1143 while 1: 1144 if gid is None: break 1145 paramId = grib_get(gid, 'paramId') 1146 gridtype = grib_get(gid, 'gridType') 1147 datatype = grib_get(gid, 'dataType') 1148 levtype = grib_get(gid, 'typeOfLevel') 1149 if paramId == 133 and gridtype=='reduced_gg': 1150 # Relative humidity (Q.grb) is used as a template only so we need the first we "meet" 1151 fout=open(c.inputdir+'/fort.18','w') 1152 grib_write(gid,fout) 1153 fout.close() 1154 elif paramId == 131 or paramId == 132: 1155 grib_write(gid, fdict['10']) 1156 elif paramId == 130: 1157 grib_write(gid, fdict['11']) 1158 elif paramId == 133 and gridtype!='reduced_gg': 1159 grib_write(gid, fdict['17']) 1160 elif paramId == 152: 1161 grib_write(gid, fdict['12']) 1162 elif paramId == 155 and gridtype=='sh': 1163 grib_write(gid, fdict['13']) 1164 elif paramId in [129,138,155] and levtype=='hybrid' and c.wrf=='1': 1165 # print paramId,'not written' 1166 pass 1167 elif paramId == 246 or paramId == 247: # cloud liquid water and ice 1168 if paramId==246: 1169 clwc=grib_get_values(gid) 1170 else: 1171 clwc+=grib_get_values(gid) 1172 grib_set_values(gid,clwc) 1173 # grib_set(gid,'shortName','qc') 1174 grib_set(gid,'paramId',201031) 1175 grib_write(gid, fdict['22']) 1176 1177 elif paramId == 135: 1178 grib_write(gid, fdict['19']) 1179 elif paramId == 77: 1180 grib_write(gid, fdict['21']) 1181 else: 1182 if paramId not in savedfields: 1183 grib_write(gid, fdict['16']) 1184 savedfields.append(paramId) 1185 else: 1186 print 'duplicate '+str(paramId)+' not written' 1187 1188 try: 1189 if c.wrf=='1': 1190 if levtype=='hybrid': 1191 if paramId in [129,130,131,132,133,138,155]: 1192 grib_write(gid,fwrf) 1193 else: 1194 if paramId in wrfpars: 1195 grib_write(gid,fwrf) 1196 except AttributeError: 1197 pass 1198 1199 1200 1201 grib_release(gid) 1202 gid = grib_new_from_index(iid) 1203 # call for CONVERT2 1204 1205 for f in fdict.values(): 1206 f.close() 1207 1208 if hid is not None: 1209 pwd=os.getcwd() 1210 os.chdir(c.inputdir) 1211 if os.stat('fort.21').st_size==0 and int(c.eta)==1: 1212 print 'Parameter 77 (etadot) is missing, most likely it is not available for this type or date/time\n' 1213 print 'Check parameters CLASS, TYPE, STREAM, START_DATE\n' 1214 myerror(c,'fort.21 is empty while parameter eta is set to 1 in CONTROL file') 1215 1216 p=subprocess.check_call([os.path.expandvars(os.path.expanduser(c.exedir))+'/CONVERT2'],shell=True) 1217 os.chdir(pwd) 1218 # create the corresponding output file fort.15 (generated by CONVERT2) + fort.16 (paramId 167 and paramId 168) 1219 fnout=c.inputdir+'/'+c.prefix 1220 if c.maxstep>12: 1221 suffix=cdate[2:8]+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step) 1222 else: 1223 suffix=cdateH[2:10] 1224 1225 fnout+=suffix 1226 print "outputfile = " + fnout 1227 self.outputfilelist.append(fnout) # needed for final processing 1228 fout = open(fnout,'wb') 1229 shutil.copyfileobj(open(c.inputdir+'/fort.15','rb'), fout) 1230 if c.cwc=='1': 1231 shutil.copyfileobj(open(c.inputdir+'/fort.22','rb'), fout) 1232 shutil.copyfileobj(open(c.inputdir+'/flux'+cdate[0:2]+suffix,'rb'), fout) 1233 shutil.copyfileobj(open(c.inputdir+'/fort.16','rb'), fout) 1234 orolsm=glob.glob(c.inputdir+'/OG_OROLSM__SL.*.'+c.ppid+'*')[0] 1235 shutil.copyfileobj(open(orolsm,'rb'), fout) 1236 fout.close() 1237 if c.omega=='1': 1238 fnout=c.outputdir+'/OMEGA' 1239 fout = open(fnout,'wb') 1240 shutil.copyfileobj(open(c.inputdir+'/fort.25','rb'), fout) 1241 1242 1243 1244 try: 1245 if c.wrf=='1': 1246 fwrf.close() 1247 except: 1248 pass 1249 1250 grib_index_release(iid) 1251 1252 def deacc_fluxes(self, inputfiles, c): 1253 1254 table128=init128(c.ecmwfdatadir+'/grib_templates/ecmwf_grib1_table_128') 1255 pars=toparamId(self.params['OG_acc_SL'][0],table128) 1256 index_keys=["date","time","step"] 1257 indexfile=c.inputdir+"/date_time_stepRange.idx" 1258 silentremove(indexfile) 1259 grib=GribTools(inputfiles.files) 1260 iid=grib.index(index_keys=index_keys, index_file = indexfile) 1261 1262 print 'index done...' 1263 1264 index_vals = [] 1265 for key in index_keys: 1266 key_vals = grib_index_get(iid,key) 1267 print key_vals 1268 if key=='step': 1269 l=[] 1270 for k in key_vals: 1271 l.append(int(k)) 1272 l.sort() 1273 key_vals=[] 1274 for k in l: 1275 key_vals.append(str(k)) 1276 1277 index_vals.append(key_vals) 1278 1279 1280 valsdict={} 1281 svalsdict={} 1282 stepsdict={} 1283 for p in pars: 1284 valsdict[str(p)]=[] 1285 svalsdict[str(p)]=[] 1286 stepsdict[str(p)]=[] 1287 1288 for prod in product(*index_vals): 1289 for i in range(len(index_keys)): 1290 grib_index_select(iid,index_keys[i],prod[i]) 1291 1292 #for k,f in fdict.iteritems(): 1293 #fdict[k] = open('fort.'+k,'w') 1294 gid = grib_new_from_index(iid) 1295 hid = gid 1296 cflextime = None 1297 if gid is not None: 1298 cdate = grib_get(gid, 'date') 1299 #cyear = cdate[:4] 1300 #cmonth = cdate[4:6] 1301 #cday = cdate[6:8] 1302 time = grib_get(gid, 'time') 1303 type = grib_get(gid, 'type') 1304 step = grib_get(gid, 'step') 1305 # date+time+step-2*dtime (since interpolated value valid for step-2*dtime) 1306 sdate=datetime.datetime(year=cdate/10000,month=mod(cdate,10000)/100, 1307 day=mod(cdate,100),hour=time/100) 1308 fdate=sdate+datetime.timedelta(hours=step-2*int(c.dtime)) 1309 sdates=sdate+datetime.timedelta(hours=step) 1310 else: 1311 break 1312 1313 fnout=c.inputdir+'/' 1314 if c.maxstep>12: 1315 fnout+='flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step-2*int(c.dtime)) 1316 gnout=c.inputdir+'/flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step-int(c.dtime)) 1317 hnout=c.inputdir+'/flux'+sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100)+'.{:0>3}'.format(step) 1318 g=open(gnout,'w') 1319 h=open(hnout,'w') 1320 else: 1321 fnout+='flux'+fdate.strftime('%Y%m%d%H') 1322 gnout=c.inputdir+'/flux'+(fdate+datetime.timedelta(hours=int(c.dtime))).strftime('%Y%m%d%H') 1323 hnout=c.inputdir+'/flux'+sdates.strftime('%Y%m%d%H') 1324 g=open(gnout,'w') 1325 h=open(hnout,'w') 1326 print "outputfile = " + fnout 1327 f=open(fnout,'w') 1328 1329 # print 'cyear '+cyear+'/'+cmonth+'/'+'/EI'+cyear[2:4]+cmonth+cday+cflextime 1330 while 1: 1331 if gid is None: break 1875 valsdict = {} 1876 svalsdict = {} 1877 stepsdict = {} 1878 for p in pars: 1879 valsdict[str(p)] = [] 1880 svalsdict[str(p)] = [] 1881 stepsdict[str(p)] = [] 1882 #AP muss das hier wirklich eingerückt sein? 1883 for prod in product(*index_vals): 1884 for i in range(len(index_keys)): 1885 grib_index_select(iid, index_keys[i], prod[i]) 1886 1887 gid = grib_new_from_index(iid) 1888 hid = gid 1889 cflextime = None 1890 if gid is not None: 1891 cdate = grib_get(gid, 'date') 1892 #cyear = cdate[:4] 1893 #cmonth = cdate[4:6] 1894 #cday = cdate[6:8] 1895 time = grib_get(gid, 'time') 1896 type = grib_get(gid, 'type') 1897 step = grib_get(gid, 'step') 1898 # date+time+step-2*dtime (since interpolated value valid for step-2*dtime) 1899 sdate = datetime.datetime(year=cdate / 10000, 1900 month=mod(cdate, 10000) / 100, 1901 day=mod(cdate, 100), 1902 hour=time / 100) 1903 fdate = sdate + datetime.timedelta( 1904 hours=step - 2 * int(c.dtime)) 1905 sdates = sdate + datetime.timedelta(hours=step) 1906 else: 1907 break 1908 1909 if c.maxstep > 12: 1910 fnout = c.inputdir + '/flux' + \ 1911 sdate.strftime('%Y%m%d') + '.{:0>2}'.format(time/100) + \ 1912 '.{:0>3}'.format(step-2*int(c.dtime)) 1913 gnout = c.inputdir + '/flux' + \ 1914 sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100) + \ 1915 '.{:0>3}'.format(step-int(c.dtime)) 1916 hnout = c.inputdir + '/flux' + \ 1917 sdate.strftime('%Y%m%d')+'.{:0>2}'.format(time/100) + \ 1918 '.{:0>3}'.format(step) 1919 g = open(gnout, 'w') 1920 h = open(hnout, 'w') 1921 else: 1922 fnout = c.inputdir + '/flux' + fdate.strftime('%Y%m%d%H') 1923 gnout = c.inputdir + '/flux' + (fdate+datetime.timedelta( 1924 hours = int(c.dtime))).strftime('%Y%m%d%H') 1925 hnout = c.inputdir + '/flux' + sdates.strftime('%Y%m%d%H') 1926 g = open(gnout, 'w') 1927 h = open(hnout, 'w') 1928 print "outputfile = " + fnout 1929 f = open(fnout, 'w') 1930 1931 while 1: 1932 if gid is None: 1933 break 1332 1934 cparamId = str(grib_get(gid, 'paramId')) 1333 1935 step = grib_get(gid, 'step') 1334 1936 atime = grib_get(gid, 'time') 1335 ni=grib_get(gid, 'Ni') 1336 nj=grib_get(gid, 'Nj') 1337 if cparamId in valsdict.keys(): 1338 values = grib_get_values(gid) 1339 vdp=valsdict[cparamId] 1340 svdp=svalsdict[cparamId] 1341 sd=stepsdict[cparamId] 1342 1343 if cparamId=='142' or cparamId=='143': 1344 fak=1./1000. 1345 else: 1346 fak=3600. 1347 1348 values=(reshape(values,(nj,ni))).flatten()/fak 1349 vdp.append(values[:]) # save the accumulated values 1350 if step<=int(c.dtime): 1351 svdp.append(values[:]/int(c.dtime)) 1352 else: 1353 svdp.append((vdp[-1]-vdp[-2])/int(c.dtime)) 1354 1355 print cparamId,atime,step,len(values),values[0],std(values) 1356 #svdp.append(values[:]) # save the 1-hourly or 3-hourly specific values 1357 sd.append(step) 1358 if len(svdp)>=3: 1359 if len(svdp)>3: 1360 if cparamId=='142' or cparamId=='143': 1361 values=darain(svdp) 1362 else: 1363 values=dapoly(svdp) 1364 1365 if not (step==c.maxstep and c.maxstep>12 or sdates==elimit): 1366 vdp.pop(0) 1367 svdp.pop(0) 1368 else: 1369 if c.maxstep>12: 1370 values=svdp[1] 1371 else: 1372 values=svdp[0] 1373 1374 grib_set_values(gid, values) 1375 if c.maxstep>12: 1376 grib_set(gid,'step',max(0,step-2*int(c.dtime))) 1377 else: 1378 grib_set(gid,'step',0) 1379 grib_set(gid,'time',fdate.hour*100) 1380 grib_set(gid,'date',fdate.year*10000+fdate.month*100+fdate.day) 1381 grib_write(gid, f) 1382 1383 if c.basetime is not None: 1384 elimit=datetime.datetime.strptime(c.end_date+c.basetime,'%Y%m%d%H') 1385 else: 1386 elimit=sdate+datetime.timedelta(2*int(c.dtime)) 1387 1388 # squeeze out information of last two steps contained in svdp 1389 # if step+int(c.dtime)==c.maxstep and c.maxstep>12 or sdates+datetime.timedelta(hours=int(c.dtime))>=elimit: 1390 # Note that svdp[0] has not been popped in this case 1391 1392 1393 if step==c.maxstep and c.maxstep>12 or sdates==elimit: 1394 values=svdp[3] 1395 grib_set_values(gid, values) 1396 grib_set(gid,'step',0) 1397 truedatetime=fdate+datetime.timedelta(hours=2*int(c.dtime)) 1398 grib_set(gid,'time',truedatetime.hour*100) 1399 grib_set(gid,'date',truedatetime.year*10000+truedatetime.month*100+truedatetime.day) 1400 grib_write(gid, h) 1401 1402 #values=(svdp[1]+svdp[2])/2. 1403 if cparamId=='142' or cparamId=='143': 1404 values=darain(list(reversed(svdp))) 1405 else: 1406 values=dapoly(list(reversed(svdp))) 1407 1408 grib_set(gid,'step',0) 1409 truedatetime=fdate+datetime.timedelta(hours=int(c.dtime)) 1410 grib_set(gid,'time',truedatetime.hour*100) 1411 grib_set(gid,'date',truedatetime.year*10000+truedatetime.month*100+truedatetime.day) 1412 grib_set_values(gid, values) 1413 grib_write(gid, g) 1414 1415 grib_release(gid) 1416 1417 gid = grib_new_from_index(iid) 1418 1419 f.close() 1420 g.close() 1421 h.close() 1422 1423 1424 grib_index_release(iid) 1425 1937 ni = grib_get(gid, 'Ni') 1938 nj = grib_get(gid, 'Nj') 1939 if cparamId in valsdict.keys(): 1940 values = grib_get_values(gid) 1941 vdp = valsdict[cparamId] 1942 svdp = svalsdict[cparamId] 1943 sd = stepsdict[cparamId] 1944 1945 if cparamId == '142' or cparamId == '143': 1946 fak = 1./1000. 1947 else: 1948 fak = 3600. 1949 1950 values = (reshape(values, (nj, ni))).flatten() / fak 1951 vdp.append(values[:]) # save the accumulated values 1952 if step <= int(c.dtime): 1953 svdp.append(values[:] / int(c.dtime)) 1954 else: 1955 svdp.append((vdp[-1] - vdp[-2]) / int(c.dtime)) 1956 1957 print cparamId, atime, step, len(values), \ 1958 values[0], std(values) 1959 # save the 1/3-hourly or specific values 1960 # svdp.append(values[:]) 1961 sd.append(step) 1962 if len(svdp) >= 3: 1963 if len(svdp) > 3: 1964 if cparamId == '142' or cparamId == '143': 1965 values = darain(svdp) 1966 else: 1967 values = dapoly(svdp) 1968 1969 if not (step == c.maxstep and c.maxstep > 12 \ 1970 or sdates == elimit): 1971 vdp.pop(0) 1972 svdp.pop(0) 1973 else: 1974 if c.maxstep > 12: 1975 values = svdp[1] 1976 else: 1977 values = svdp[0] 1978 1979 grib_set_values(gid, values) 1980 if c.maxstep > 12: 1981 grib_set(gid, 'step', max(0, step-2*int(c.dtime))) 1982 else: 1983 grib_set(gid, 'step', 0) 1984 grib_set(gid, 'time', fdate.hour*100) 1985 grib_set(gid, 'date', fdate.year*10000 + 1986 fdate.month*100+fdate.day) 1987 grib_write(gid, f) 1988 1989 if c.basetime is not None: 1990 elimit = datetime.datetime.strptime(c.end_date + 1991 c.basetime, 1992 '%Y%m%d%H') 1993 else: 1994 elimit = sdate + datetime.timedelta(2*int(c.dtime)) 1995 1996 # squeeze out information of last two steps contained 1997 # in svdp 1998 # if step+int(c.dtime) == c.maxstep and c.maxstep>12 1999 # or sdates+datetime.timedelta(hours = int(c.dtime)) 2000 # >= elimit: 2001 # Note that svdp[0] has not been popped in this case 2002 2003 if (step == c.maxstep and c.maxstep > 12 2004 or sdates == elimit): 2005 values = svdp[3] 2006 grib_set_values(gid, values) 2007 grib_set(gid, 'step', 0) 2008 truedatetime = fdate + datetime.timedelta( 2009 hours=2*int(c.dtime)) 2010 grib_set(gid, 'time', truedatetime.hour * 100) 2011 grib_set(gid, 'date', truedatetime.year * 10000 + 2012 truedatetime.month * 100 + 2013 truedatetime.day) 2014 grib_write(gid, h) 2015 2016 #values = (svdp[1]+svdp[2])/2. 2017 if cparamId == '142' or cparamId == '143': 2018 values = darain(list(reversed(svdp))) 2019 else: 2020 values = dapoly(list(reversed(svdp))) 2021 2022 grib_set(gid, 'step',0) 2023 truedatetime = fdate + datetime.timedelta( 2024 hours=int(c.dtime)) 2025 grib_set(gid, 'time', truedatetime.hour * 100) 2026 grib_set(gid, 'date', truedatetime.year * 10000 + 2027 truedatetime.month * 100 + 2028 truedatetime.day) 2029 grib_set_values(gid, values) 2030 grib_write(gid, g) 2031 2032 grib_release(gid) 2033 2034 gid = grib_new_from_index(iid) 2035 2036 f.close() 2037 g.close() 2038 h.close() 2039 2040 grib_index_release(iid) 2041 2042 return 2043 -
python/GribTools.py
rd69b677 r64cf353 1 # 2 # (C) Copyright 2014 UIO. 3 # 4 # This software is licensed under the terms of the Apache Licence Version 2.0 5 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 6 # 7 # 8 # Creation: July 2014 - Anne Fouilloux - University of Oslo 9 # 10 # 11 1 #!/usr/bin/env python 2 # -*- coding: utf-8 -*- 3 #************************************************************************ 4 # TODO AP 5 #AP 6 # - GribTools name möglicherweise etwas verwirrend. 7 # - change self.filename in self.filenames!!! 8 # - 9 #************************************************************************ 10 """ 11 @Author: Anne Fouilloux (University of Oslo) 12 13 @Date: July 2014 14 15 @ChangeHistory: 16 February 2018 - Anne Philipp (University of Vienna): 17 - applied PEP8 style guide 18 - added documentation 19 - changed some naming 20 21 @License: 22 (C) Copyright 2014 UIO. 23 24 This software is licensed under the terms of the Apache Licence Version 2.0 25 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 26 27 @Requirements: 28 - A standard python 2.6 or 2.7 installation 29 - dateutils 30 - matplotlib (optional, for debugging) 31 - ECMWF specific packages, all available from https://software.ecmwf.int/ 32 ECMWF WebMARS, gribAPI with python enabled, emoslib and 33 ecaccess web toolkit 34 35 @Description: 36 Further documentation may be obtained from www.flexpart.eu. 37 38 Functionality provided: 39 Prepare input 3D-wind fields in hybrid coordinates + 40 surface fields for FLEXPART runs 41 """ 42 # ------------------------------------------------------------------------------ 43 # MODULES 44 # ------------------------------------------------------------------------------ 12 45 from gribapi import * 13 46 import traceback 14 import sys,os 15 16 17 ############################################################## 18 # GRIB utilities 19 ############################################################## 47 import sys, os 48 49 # ------------------------------------------------------------------------------ 50 # CLASS 51 # ------------------------------------------------------------------------------ 20 52 class GribTools: 21 'class for GRIB API with new methods' 22 def __init__(self,filename): 23 self.filename=filename 24 25 # get keyvalues for a given list of keynames 26 # a where statment can be given (list of key and list of values) 27 28 def getkeys(self,keynames,wherekeynames=[],wherekeyvalues=[]): 29 fileid=open(self.filename,'r') 30 31 return_list=[] 32 53 ''' 54 Class for GRIB utilities (new methods) based on GRIB API 55 ''' 56 # -------------------------------------------------------------------------- 57 # FUNCTIONS 58 # -------------------------------------------------------------------------- 59 def __init__(self, filenames): 60 ''' 61 @Description: 62 Initialise an object of GribTools and assign a list 63 of filenames. 64 65 @Input: 66 filenames: list of strings 67 A list of filenames. 68 69 @Return: 70 <nothing> 71 ''' 72 73 self.filename = filenames 74 75 return 76 77 78 def getkeys(self, keynames, wherekeynames=[], wherekeyvalues=[]): 79 ''' 80 @Description: 81 get keyvalues for a given list of keynames 82 a where statement can be given (list of key and list of values) 83 84 @Input: 85 keynames: list of strings 86 List of keynames. 87 88 wherekeynames: list of ??? 89 Default value is an empty list. 90 91 wherekeyvalues: list of ??? 92 Default value is an empty list. 93 94 @Return: 95 return_list: list of strings 96 List of keyvalues for given keynames. 97 ''' 98 99 fileid = open(self.filename, 'r') 100 101 return_list = [] 102 33 103 while 1: 34 gid_in = grib_new_from_file(fileid) 35 if gid_in is None: break 36 select=True 37 i=0 38 if len(wherekeynames) != len(wherekeyvalues): raise Exception("Give a value for each keyname!") 39 104 gid_in = grib_new_from_file(fileid) 105 106 if gid_in is None: 107 break 108 109 if len(wherekeynames) != len(wherekeyvalues): 110 raise Exception("Number of key values and key names must be \ 111 the same. Give a value for each keyname!") 112 113 select = True 114 i = 0 40 115 for wherekey in wherekeynames: 41 if not grib_is_defined(gid_in, wherekey): raise Exception("where Key was not defined") 42 select=select and (str(wherekeyvalues[i])==str(grib_get(gid_in, wherekey))) 43 i=i+1 116 if not grib_is_defined(gid_in, wherekey): 117 raise Exception("where key was not defined") 118 119 select = (select and (str(wherekeyvalues[i]) == 120 str(grib_get(gid_in, wherekey)))) 121 i += 1 122 44 123 if select: 45 124 llist = [] … … 47 126 llist.extend([str(grib_get(gid_in, key))]) 48 127 return_list.append(llist) 128 49 129 grib_release(gid_in) 130 50 131 fileid.close() 132 51 133 return return_list 52 53 # set keyvalues for a given list of keynames 54 # a where statment can be given (list of key and list of values) 55 # an input file must be given as an input for reading grib messages 56 # note that by default all messages are written out 57 # if you want to get only those meeting the where statement, use 58 # strict=true 59 def setkeys(self,fromfile,keynames,keyvalues, wherekeynames=[],wherekeyvalues=[], strict=False, filemode='w'): 60 fout=open(self.filename,filemode) 61 fin=open(fromfile) 62 134 135 136 def setkeys(self, fromfile, keynames, keyvalues, wherekeynames=[], 137 wherekeyvalues=[], strict=False, filemode='w'): 138 ''' 139 @Description: 140 Opens the file to read the grib messages and then write 141 them to a new output file. By default all messages are 142 written out. Also, the keyvalues of the passed list of 143 keynames are set or only those meeting the where statement. 144 (list of key and list of values). 145 146 @Input: 147 fromfile: string 148 Filename of the input file to read the grib messages from. 149 150 keynames: list of ??? 151 List of keynames. Default is an empty list. 152 153 keyvalues: list of ??? 154 List of keynames. Default is an empty list. 155 156 wherekeynames: list of ??? 157 Default value is an empty list. 158 159 wherekeyvalues: list of ??? 160 Default value is an empty list. 161 162 strict: boolean 163 Decides if everything from keynames and keyvalues 164 is written out the grib file (False) or only those 165 meeting the where statement (True). Default is False. 166 167 filemode: 168 Sets the mode for the output file. Default is "w". 169 170 @Return: 171 <nothing> 172 173 ''' 174 fout = open(self.filename, filemode) 175 fin = open(fromfile) 176 63 177 while 1: 64 gid_in = grib_new_from_file(fin) 65 if gid_in is None: break 66 67 select=True 68 i=0 69 if len(wherekeynames) != len(wherekeyvalues): raise Exception("Give a value for each keyname!") 70 178 gid_in = grib_new_from_file(fin) 179 180 if gid_in is None: 181 break 182 183 if len(wherekeynames) != len(wherekeyvalues): 184 raise Exception("Give a value for each keyname!") 185 186 select = True 187 i = 0 71 188 for wherekey in wherekeynames: 72 if not grib_is_defined(gid_in, wherekey): raise Exception("where Key was not defined") 73 select=select and (str(wherekeyvalues[i])==str(grib_get(gid_in, wherekey))) 74 i=i+1 189 if not grib_is_defined(gid_in, wherekey): 190 raise Exception("where Key was not defined") 191 192 select = (select and (str(wherekeyvalues[i]) == 193 str(grib_get(gid_in, wherekey)))) 194 i += 1 195 196 #AP is it secured that the order of keynames is equal to keyvalues? 75 197 if select: 76 i =0198 i = 0 77 199 for key in keynames: 78 200 grib_set(gid_in, key, keyvalues[i]) 79 i=i+1 201 i += 1 202 203 #AP this is a redundant code section 204 # delete the if/else : 205 # 206 # grib_write(gid_in, fout) 207 # 80 208 if strict: 81 209 if select: 82 grib_write(gid_in, fout)210 grib_write(gid_in, fout) 83 211 else: 84 grib_write(gid_in,fout) 212 grib_write(gid_in, fout) 213 #AP end 85 214 grib_release(gid_in) 215 86 216 fin.close() 87 217 fout.close() 88 218 89 # Add the content of a grib file but only messages 90 # corresponding to keys/values 91 # if selectWhere is False select fields that are different from keynames/keyvalues 92 def copy(self,filename_in, selectWhere=True, keynames=[], keyvalues=[],filemode='w'): 93 fin=open(filename_in) 94 fout=open(self.filename,filemode) 95 219 return 220 221 def copy(self, filename_in, selectWhere=True, 222 keynames=[], keyvalues=[], filemode='w'): 223 ''' 224 Add the content of another input grib file to the objects file but 225 only messages corresponding to keys/values passed to the function. 226 The selectWhere switch decides if to copy the keys equal to (True) or 227 different to (False) the keynames/keyvalues list passed to the function. 228 229 @Input: 230 filename_in: string 231 Filename of the input file to read the grib messages from. 232 233 selectWhere: boolean 234 Decides if to copy the keynames and values equal to (True) or 235 different to (False) the keynames/keyvalues list passed to the 236 function. Default is True. 237 238 keynames: list of ??? 239 List of keynames. Default is an empty list. 240 241 keyvalues: list of ??? 242 List of keynames. Default is an empty list. 243 244 filemode: 245 Sets the mode for the output file. Default is "w". 246 247 @Return: 248 <nothing> 249 ''' 250 251 fin = open(filename_in) 252 fout = open(self.filename, filemode) 253 96 254 while 1: 97 gid_in = grib_new_from_file(fin) 98 if gid_in is None: break 99 100 select=True 101 i=0 102 if len(keynames) != len(keyvalues): raise Exception("Give a value for each keyname!") 103 255 gid_in = grib_new_from_file(fin) 256 257 if gid_in is None: 258 break 259 260 if len(keynames) != len(keyvalues): 261 raise Exception("Give a value for each keyname!") 262 263 select = True 264 i = 0 104 265 for key in keynames: 105 if not grib_is_defined(gid_in, key): raise Exception("Key was not defined") 106 266 if not grib_is_defined(gid_in, key): 267 raise Exception("Key was not defined") 268 107 269 if selectWhere: 108 select=select and (str(keyvalues[i])==str(grib_get(gid_in, key))) 270 select = (select and (str(keyvalues[i]) == 271 str(grib_get(gid_in, key)))) 109 272 else: 110 select=select and (str(keyvalues[i])!=str(grib_get(gid_in, key))) 111 i=i+1 273 select = (select and (str(keyvalues[i]) != 274 str(grib_get(gid_in, key)))) 275 i += 1 276 112 277 if select: 113 grib_write(gid_in,fout) 278 grib_write(gid_in, fout) 279 114 280 grib_release(gid_in) 281 115 282 fin.close() 116 283 fout.close() 117 284 118 # Create index from a list of files if it does not exist or read it 119 def index(self,index_keys=["mars"], index_file = "my.idx"): 120 print "index to be done" 285 return 286 287 def index(self, index_keys=["mars"], index_file="my.idx"): 288 ''' 289 @Description: 290 Create index from a list of files if it does not exist or 291 read an index file. 292 293 @Input: 294 index_keys: list of ??? 295 Default is a list with a single entry string "mars". 296 297 index_file: string 298 Filename where the indices are stored. 299 Default is "my.idx". 300 301 @Return: 302 iid: integer 303 Grib index id. 304 ''' 305 print("... index will be done") 121 306 self.iid = None 122 307 123 308 if (os.path.exists(index_file)): 124 309 self.iid = grib_index_read(index_file) 125 print "Use existing index file: %s "%(index_file)310 print("Use existing index file: %s " % (index_file)) 126 311 else: 312 #AP does the for loop overwrite the iid all the time? 127 313 for file in self.filename: 128 print "Inputfile: %s "%(file)129 if self.iid ==None:130 self.iid = grib_index_new_from_file(file, index_keys)314 print("Inputfile: %s " % (file)) 315 if self.iid is None: 316 self.iid = grib_index_new_from_file(file, index_keys) 131 317 else: 132 grib_index_add_file(self.iid,file) 133 134 if self.iid != None: 135 grib_index_write(self.iid,index_file) 136 return self.iid 137 138 139 140 141 318 grib_index_add_file(self.iid, file) 319 #AP or does the if has to be in the for loop? 320 #AP would make more sense? 321 if self.iid is not None: 322 grib_index_write(self.iid, index_file) 323 324 return self.iid 325 326 327 328 329 -
python/UIOTools.py
r6180177 r64cf353 1 #!/usr/bin/env python 2 # -*- coding: utf-8 -*- 1 3 #************************************************************************ 2 4 # TODO AP 3 # 5 #AP 4 6 # - File name und Klassenname gleichsetzen? 5 7 # - checken welche regelmässigen methoden auf diese Files noch angewendet werden … … 13 15 14 16 @ChangeHistory: 15 Anne Philipp (University of Vienna) - February 2018: 16 Added documentation and applied pep8 style guides 17 February 2018 - Anne Philipp (University of Vienna): 18 - applied PEP8 style guide 19 - added documentation 17 20 18 21 @License: … … 21 24 This software is licensed under the terms of the Apache Licence Version 2.0 22 25 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 26 27 @Requirements: 28 - A standard python 2.6 or 2.7 installation 29 - dateutils 30 - matplotlib (optional, for debugging) 31 - ECMWF specific packages, all available from https://software.ecmwf.int/ 32 ECMWF WebMARS, gribAPI with python enabled, emoslib and 33 ecaccess web toolkit 34 35 @Description: 36 Further documentation may be obtained from www.flexpart.eu. 37 38 Functionality provided: 39 Prepare input 3D-wind fields in hybrid coordinates + 40 surface fields for FLEXPART runs 23 41 """ 24 25 26 42 # ------------------------------------------------------------------------------ 27 43 # MODULES … … 54 70 suffix: list of strings 55 71 Types of files to manipulate such as 56 [' .grib', 'grb', 'grib1', 'grib2', 'grb1','grb2']72 ['grib', 'grb', 'grib1', 'grib2', 'grb1', 'grb2'] 57 73 58 74 @Return: … … 83 99 <nothing> 84 100 ''' 101 #AP pathname zu path ändern 102 #AP is it possible for each possible file extension ? mabye regexx? 85 103 # Get the absolute path of the pathname parameter 86 104 pathname = os.path.abspath(pathname) -
python/getMARSdata.py
rd69b677 r64cf353 2 2 # 3 3 # This software is licensed under the terms of the Apache Licence Version 2.0 4 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 5 # 4 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 5 # 6 6 # Functionality provided: Prepare input 3D-wind fields in hybrid coordinates + surface fields for FLEXPART runs 7 7 # 8 8 # Creation: October 2014 - Anne Fouilloux - University of Oslo 9 # Extension November 2015 - Leopold Haimberger - University of Vienna for: 9 # Extension November 2015 - Leopold Haimberger - University of Vienna for: 10 10 # - using the WebAPI also for general MARS retrievals 11 11 # - job submission on ecgate and cca … … 15 15 # - conversion into GRIB2 16 16 # - conversion into .fp format for faster execution of FLEXPART 17 # 17 # 18 18 # 19 19 # Further documentation may be obtained from www.flexpart.eu 20 # 21 # Requirements: 20 # 21 # Requirements: 22 22 # in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed 23 23 # ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/ … … 26 26 # 27 27 # Get MARS GRIB fields from ECMWF for FLEXPART 28 # 28 # 29 29 30 30 #import socket … … 51 51 52 52 from FlexpartTools import MARSretrieval, EIFlexpart, silentremove, \ 53 Control,myerror,normalexit, interpret_args_and_control 53 Control, myerror, normalexit, \ 54 interpret_args_and_control 54 55 55 56 56 def getMARSdata(args, c):57 58 57 def getMARSdata(args, c): 58 59 59 60 if not os.path.exists(c.inputdir): 60 61 os.makedirs(c.inputdir) 61 print "start date %s "%(c.start_date)62 print "end date %s "%(c.end_date)62 print("start date %s " % (c.start_date)) 63 print("end date %s " % (c.end_date)) 63 64 64 65 # server = ECMWFDataServer()66 65 if ecapi: 67 66 server = ecmwfapi.ECMWFService("mars") … … 69 68 server = False 70 69 71 c.ecapi=ecapi 72 print 'ecapi:',c.ecapi 70 c.ecapi = ecapi 71 print 'ecapi:', c.ecapi 72 73 73 # Retrieve ERA interim data for running flexpart 74 #AP change this variant to correct format conversion with datetime 75 #AP import datetime and timedelta explicitly 76 syear = int(c.start_date[:4]) 77 smonth = int(c.start_date[4:6]) 78 sday = int(c.start_date[6:]) 79 start = datetime.date(year=syear, month=smonth, day=sday) 80 startm1 = start - datetime.timedelta(days=1) 81 if c.basetime == '00': 82 start = startm1 83 eyear = int(c.end_date[:4]) 84 emonth = int(c.end_date[4:6]) 85 eday = int(c.end_date[6:]) 86 end = datetime.date(year=eyear, month=emonth, day=eday) 87 if c.basetime == '00' or c.basetime == '12': 88 endp1 = end + datetime.timedelta(days=1) 89 else: 90 endp1 = end + datetime.timedelta(days=2) 74 91 75 syear=int(c.start_date[:4]) 76 smonth=int(c.start_date[4:6]) 77 sday=int(c.start_date[6:]) 78 start = datetime.date( year = syear, month = smonth, day = sday ) 79 startm1=start- datetime.timedelta(days=1) 80 if c.basetime=='00': 81 start=startm1 82 eyear=int(c.end_date[:4]) 83 emonth=int(c.end_date[4:6]) 84 eday=int(c.end_date[6:]) 85 86 end = datetime.date( year = eyear, month = emonth, day = eday ) 87 if c.basetime=='00' or c.basetime=='12': 88 endp1=end+ datetime.timedelta(days=1) 89 else: 90 endp1=end+ datetime.timedelta(days=2) 91 92 datechunk=datetime.timedelta(days=int(c.date_chunk)) 93 print 'removing content of '+c.inputdir 94 tobecleaned=glob.glob(c.inputdir+'/*_acc_*.'+str(os.getppid())+'.*.grb') 92 datechunk = datetime.timedelta(days=int(c.date_chunk)) 93 print 'removing content of ' + c.inputdir 94 tobecleaned = glob.glob(c.inputdir + '/*_acc_*.' + str(os.getppid()) + '.*.grb') 95 95 for f in tobecleaned: 96 96 os.remove(f) 97 97 98 98 times=None 99 99 if c.maxstep<24: … … 103 103 flexpart = EIFlexpart(c,fluxes=True) 104 104 if day+datechunk-datetime.timedelta(days=1)<endp1: 105 dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d") 105 dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d") 106 106 else: 107 dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d") 108 107 dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d") 108 109 109 print "retrieve " + dates + " in dir " + c.inputdir 110 110 try: … … 112 112 except IOError: 113 113 myerror(c,'MARS request failed') 114 114 115 115 day+=datechunk 116 116 else: … … 120 120 flexpart = EIFlexpart(c,fluxes=True) 121 121 if day+datechunk-datetime.timedelta(days=1)<end: 122 dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d") 122 dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d") 123 123 else: 124 dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d") 125 124 dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d") 125 126 126 print "retrieve " + dates + " in dir " + c.inputdir 127 127 flexpart.retrieve(server, dates, times, c.inputdir) 128 128 day+=datechunk 129 129 130 130 131 131 … … 136 136 times=None 137 137 while day<=end: 138 138 139 139 # we need to retrieve MARS data for this period (maximum one month) 140 140 flexpart = EIFlexpart(c) 141 141 if day+datechunk-datetime.timedelta(days=1)<end: 142 dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d") 142 dates= day.strftime("%Y%m%d") + "/to/" + (day+datechunk-datetime.timedelta(days=1)).strftime("%Y%m%d") 143 143 else: 144 dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d") 144 dates= day.strftime("%Y%m%d") + "/to/" + end.strftime("%Y%m%d") 145 145 print "retrieve " + dates + " in dir " + c.inputdir 146 146 147 147 flexpart.retrieve(server, dates, times, c.inputdir) 148 148 day+=datechunk 149 149 150 150 151 151 if __name__ == "__main__": 152 152 153 153 args,c=interpret_args_and_control() 154 154 getMARSdata(args,c) -
python/install.py
ra4b6cef r64cf353 1 1 #!/usr/bin/env python 2 # 3 # This software is licensed under the terms of the Apache Licence Version 2.0 4 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 5 # 6 # Functionality provided: Prepare input 3D-wind fields in hybrid coordinates + surface fields for FLEXPART runs 7 # 8 # Creation: October 2014 - Anne Fouilloux - University of Oslo 9 # Extension November 2015 - Leopold Haimberger - University of Vienna for: 10 # - using the WebAPI also for general MARS retrievals 11 # - job submission on ecgate and cca 12 # - job templates suitable for twice daily operational dissemination 13 # - dividing retrievals of longer periods into digestable chunks 14 # - retrieve also longer term forecasts, not only analyses and short term forecast data 15 # - conversion into GRIB2 16 # - conversion into .fp format for faster execution of FLEXPART 17 # 18 # 19 # Further documentation may be obtained from www.flexpart.eu 20 # 21 # Requirements: 22 # in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed 23 # ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/ 24 # dateutils 25 # matplotlib (optional, for debugging) 26 # 27 # 2 # -*- coding: utf-8 -*- 3 #************************************************************************ 4 # TODO AP 5 #AP 6 # - Functionality Provided is not correct for this file 7 # - localpythonpath should not be set in module load section! 8 # - create a class Installation and divide installation in 3 subdefs for 9 # ecgate, local and cca seperatly 10 # - Change History ist nicht angepasst ans File! Original geben lassen 11 #************************************************************************ 12 """ 13 @Author: Anne Fouilloux (University of Oslo) 14 15 @Date: October 2014 16 17 @ChangeHistory: 18 November 2015 - Leopold Haimberger (University of Vienna): 19 - using the WebAPI also for general MARS retrievals 20 - job submission on ecgate and cca 21 - job templates suitable for twice daily operational dissemination 22 - dividing retrievals of longer periods into digestable chunks 23 - retrieve also longer term forecasts, not only analyses and 24 short term forecast data 25 - conversion into GRIB2 26 - conversion into .fp format for faster execution of FLEXPART 27 28 February 2018 - Anne Philipp (University of Vienna): 29 - applied PEP8 style guide 30 - added documentation 31 32 @License: 33 (C) Copyright 2014 UIO. 34 35 This software is licensed under the terms of the Apache Licence Version 2.0 36 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 37 38 @Requirements: 39 - A standard python 2.6 or 2.7 installation 40 - dateutils 41 - matplotlib (optional, for debugging) 42 - ECMWF specific packages, all available from https://software.ecmwf.int/ 43 ECMWF WebMARS, gribAPI with python enabled, emoslib and 44 ecaccess web toolkit 45 46 @Description: 47 Further documentation may be obtained from www.flexpart.eu. 48 49 Functionality provided: 50 Prepare input 3D-wind fields in hybrid coordinates + 51 surface fields for FLEXPART runs 52 """ 53 # ------------------------------------------------------------------------------ 54 # MODULES 55 # ------------------------------------------------------------------------------ 28 56 import calendar 29 57 import shutil … … 44 72 from prepareFLEXPART import prepareFLEXPART 45 73 46 74 # ------------------------------------------------------------------------------ 75 # FUNCTIONS 76 # ------------------------------------------------------------------------------ 47 77 def main(): 48 78 ''' 49 79 ''' 50 # calledfromdir = os.getcwd()51 80 os.chdir(localpythonpath) 52 81 args, c = install_args_and_control() 53 # if c.outputdir[0]!='/':54 # c.outputdir=os.path.join(calledfromdir,c.outputdir)55 # c.inputdir=c.outputdir56 82 if args.install_target is not None: 57 83 install_via_gateway(c, args.install_target) … … 230 256 print('You should get an email with subject flexcompile \ 231 257 within the next few minutes') 258 232 259 else: 233 260 print('ERROR: unknown installation target ', target) -
python/plot_retrieved.py
rd69b677 r64cf353 1 1 #!/usr/bin/env python 2 2 # This software is licensed under the terms of the Apache Licence Version 2.0 3 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 4 # 3 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 4 # 5 5 # Functionality provided: Simple tool for creating maps and time series of retrieved fields. 6 # 7 # Requirements: 6 # 7 # Requirements: 8 8 # in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed 9 9 # ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/ … … 19 19 if localpythonpath not in sys.path: 20 20 sys.path.append(localpythonpath) 21 21 22 22 from matplotlib.pylab import * 23 23 import matplotlib.patches as mpatches … … 43 43 c.levels=asarray(c.levels,dtype='int') 44 44 c.area=asarray(c.area) 45 45 46 46 index_keys=["date","time","step"] 47 47 indexfile=c.inputdir+"/date_time_stepRange.idx" … … 67 67 68 68 index_vals.append(key_vals) 69 69 70 70 fdict=dict() 71 71 fmeta=dict() … … 101 101 fdatetimestep=fdatetime+datetime.timedelta(hours=step) 102 102 if len(fstamp)==0: 103 fstamp[key].append(fdatetimestamp) 103 fstamp[key].append(fdatetimestamp) 104 104 fmeta[key].append((paramId,parameterName,gtype,fdatetime,gtime,step,level)) 105 105 fdict[key].append(flipud(reshape(grib_get_values(gid),[gdict['Nj'],gdict['Ni']]))) … … 118 118 fmeta[key].append((paramId,parameterName,gtype,fdatetime,gtime,step,level)) 119 119 fdict[key].append(flipud(reshape(grib_get_values(gid),[gdict['Nj'],gdict['Ni']]))) 120 120 121 121 122 122 grib_release(gid) 123 123 gid = grib_new_from_index(iid) 124 124 125 125 for k in fdict.keys(): 126 126 fml=fmeta[k] 127 127 fdl=fdict[k] 128 128 129 129 for fd,fm in zip(fdl,fml): 130 130 ftitle=fm[1]+' {} '.format(fm[-1])+datetime.datetime.strftime(fm[3],'%Y%m%d%H')+' '+stats(fd) 131 131 pname='_'.join(fm[1].split())+'_{}_'.format(fm[-1])+datetime.datetime.strftime(fm[3],'%Y%m%d%H')+'.{:0>3}'.format(fm[5]) 132 132 plotmap(fd, fm,gdict,ftitle,pname+'.eps') 133 133 134 134 for k in fdict.keys(): 135 135 fml=fmeta[k] … … 143 143 lat=-20 144 144 lon=20 145 plottimeseries(fdl,fml,fsl,lat,lon, gdict, ftitle, pname+'.eps') 146 145 plottimeseries(fdl,fml,fsl,lat,lon, gdict, ftitle, pname+'.eps') 146 147 147 def plottimeseries(flist,fmetalist,ftimestamps,lat,lon,gdict,ftitle,filename): 148 148 t1=time.time() … … 158 158 plt.close(f) 159 159 print time.time()-t1,'s' 160 160 161 161 def plotmap(flist,fmetalist,gdict,ftitle,filename): 162 162 t1=time.time() 163 163 f=plt.figure(figsize=(12,6.7)) 164 mbaxes = f.add_axes([0.05, 0.15, 0.8, 0.7]) 164 mbaxes = f.add_axes([0.05, 0.15, 0.8, 0.7]) 165 165 m =Basemap(llcrnrlon=-180.,llcrnrlat=-90.,urcrnrlon=180,urcrnrlat=90.) 166 166 #if bw==0 : … … 168 168 #else: 169 169 #fill_color=rgb(0.85,0.85,0.85) 170 170 171 171 lw=0.3 172 172 m.drawmapboundary() … … 179 179 xleft=gdict['longitudeOfFirstGridPointInDegrees'] 180 180 if xleft>180.0: 181 xleft-=360. 181 xleft-=360. 182 182 x=linspace(xleft,gdict['longitudeOfLastGridPointInDegrees'],gdict['Ni']) 183 183 y=linspace(gdict['latitudeOfLastGridPointInDegrees'],gdict['latitudeOfFirstGridPointInDegrees'],gdict['Nj']) … … 186 186 s=m.contourf(xx,yy, flist) 187 187 title(ftitle,y=1.1) 188 cbaxes = f.add_axes([0.9, 0.2, 0.04, 0.6]) 188 cbaxes = f.add_axes([0.9, 0.2, 0.04, 0.6]) 189 189 cb=colorbar(cax=cbaxes) 190 190 191 191 savefig(c.outputdir+'/'+filename) 192 192 print 'created ',c.outputdir+'/'+filename … … 200 200 formatter_class=ArgumentDefaultsHelpFormatter) 201 201 202 # the most important arguments 202 # the most important arguments 203 203 parser.add_argument("--start_date", dest="start_date", 204 204 help="start date YYYYMMDD") … … 241 241 print 'Try "'+sys.argv[0].split('/')[-1]+' -h" to print usage information' 242 242 exit(1) 243 244 if args.levelist: 243 244 if args.levelist: 245 245 c.levels=args.levelist.split('/') 246 246 else: … … 250 250 else: 251 251 c.area='[0,0]' 252 252 253 253 c.paramIds=args.paramIds.split('/') 254 254 if args.start_step: … … 269 269 else: 270 270 c.outputdir=c.inputdir 271 271 272 272 return args,c 273 273 274 274 if __name__ == "__main__": 275 275 276 276 args,c=interpret_plotargs() 277 277 plot_retrieved(args,c) -
python/prepareFLEXPART.py
rd69b677 r64cf353 1 1 #!/usr/bin/env python 2 # 2 # 3 3 # This software is licensed under the terms of the Apache Licence Version 2.0 4 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 5 # 4 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 5 # 6 6 # Functionality provided: Prepare input 3D-wind fields in hybrid coordinates + surface fields for FLEXPART runs 7 7 # 8 8 # Creation: October 2014 - Anne Fouilloux - University of Oslo 9 # Extension November 2015 - Leopold Haimberger - University of Vienna for: 9 # Extension November 2015 - Leopold Haimberger - University of Vienna for: 10 10 # - using the WebAPI also for general MARS retrievals 11 11 # - job submission on ecgate and cca … … 15 15 # - conversion into GRIB2 16 16 # - conversion into .fp format for faster execution of FLEXPART 17 # 18 # Requirements: 17 # 18 # Requirements: 19 19 # in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed 20 # ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/ 20 # ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, 21 # all available from https://software.ecmwf.int/ 21 22 # dateutils 22 23 # matplotlib (optional, for debugging) 23 # 24 # 24 25 import calendar 25 26 import shutil … … 36 37 #from string import strip 37 38 from GribTools import GribTools 38 from FlexpartTools import EIFlexpart, Control, interpret_args_and_control, cleanup39 from FlexpartTools import EIFlexpart, Control, interpret_args_and_control, cleanup 39 40 40 41 hostname=socket.gethostname() … … 45 46 import ecmwfapi 46 47 except ImportError: 47 ecapi =False48 48 ecapi = False 49 49 50 50 51 def prepareFLEXPART(args,c): … … 53 54 54 55 namelist='fort.4' 55 56 56 57 if not args.ppid: 57 58 c.ppid=str(os.getppid()) … … 78 79 inputfiles.listFiles(c.inputdir, '*OG_acc_SL*.'+c.ppid+'.*') 79 80 if not os.path.exists(c.outputdir): 80 81 81 os.makedirs(c.outputdir) 82 82 83 flexpart = EIFlexpart(c,fluxes=True) 83 84 flexpart.create_namelist(c,'fort.4') … … 97 98 98 99 if int(c.debug)!=0: 99 100 print('Temporary files left intact') 100 101 else: 101 102 cleanup(c) 102 103 103 104 if __name__ == "__main__": 104 args, c=interpret_args_and_control()105 prepareFLEXPART(args, c)105 args, c = interpret_args_and_control() 106 prepareFLEXPART(args, c) 106 107 cleanup(c) 107 108 -
python/submit.py
rd69b677 r64cf353 1 1 #!/usr/bin/env python 2 # 3 # This software is licensed under the terms of the Apache Licence Version 2.0 4 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 5 # 6 # Functionality provided: Prepare input 3D-wind fields in hybrid coordinates + surface fields for FLEXPART runs 7 # 8 # Creation: October 2014 - Anne Fouilloux - University of Oslo 9 # Extension November 2015 - Leopold Haimberger - University of Vienna for: 10 # - using the WebAPI also for general MARS retrievals 11 # - job submission on ecgate and cca 12 # - job templates suitable for twice daily operational dissemination 13 # - dividing retrievals of longer periods into digestable chunks 14 # - retrieve also longer term forecasts, not only analyses and short term forecast data 15 # - conversion into GRIB2 16 # - conversion into .fp format for faster execution of FLEXPART 17 # 18 # Further documentation may be obtained from www.flexpart.eu 19 # 20 # Requirements: 21 # in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed 22 # ECMWF WebMARS, gribAPI with python enabled, emoslib, ecaccess web toolkit, all available from https://software.ecmwf.int/ 23 # dateutils 24 # matplotlib (optional, for debugging) 25 # 26 # 2 # -*- coding: utf-8 -*- 3 #************************************************************************ 4 # TODO AP 5 #AP 6 # - Change History ist nicht angepasst ans File! Original geben lassen 7 # - dead code ? what to do? 8 # - seperate operational and reanlysis for clarification 9 #************************************************************************ 10 """ 11 @Author: Anne Fouilloux (University of Oslo) 12 13 @Date: October 2014 14 15 @ChangeHistory: 16 November 2015 - Leopold Haimberger (University of Vienna): 17 - using the WebAPI also for general MARS retrievals 18 - job submission on ecgate and cca 19 - job templates suitable for twice daily operational dissemination 20 - dividing retrievals of longer periods into digestable chunks 21 - retrieve also longer term forecasts, not only analyses and 22 short term forecast data 23 - conversion into GRIB2 24 - conversion into .fp format for faster execution of FLEXPART 25 26 February 2018 - Anne Philipp (University of Vienna): 27 - applied PEP8 style guide 28 - added documentation 29 - minor changes in programming style for consistence 30 31 @License: 32 (C) Copyright 2014 UIO. 33 34 This software is licensed under the terms of the Apache Licence Version 2.0 35 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 36 37 @Requirements: 38 - A standard python 2.6 or 2.7 installation 39 - dateutils 40 - matplotlib (optional, for debugging) 41 - ECMWF specific packages, all available from https://software.ecmwf.int/ 42 ECMWF WebMARS, gribAPI with python enabled, emoslib and 43 ecaccess web toolkit 44 45 @Description: 46 Further documentation may be obtained from www.flexpart.eu. 47 48 Functionality provided: 49 Prepare input 3D-wind fields in hybrid coordinates + 50 surface fields for FLEXPART runs 51 """ 52 # ------------------------------------------------------------------------------ 53 # MODULES 54 # ------------------------------------------------------------------------------ 27 55 import calendar 28 56 import shutil 29 57 import datetime 30 58 import time 31 import os, sys,glob59 import os, sys, glob 32 60 import subprocess 33 61 import inspect 34 62 # add path to submit.py to pythonpath so that python finds its buddies 35 localpythonpath =os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))63 localpythonpath = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe()))) 36 64 sys.path.append(localpythonpath) 37 65 from UIOTools import UIOFiles 38 66 from string import strip 39 from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter67 from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter 40 68 from GribTools import GribTools 41 from FlexpartTools import EIFlexpart, Control, interpret_args_and_control, normalexit, myerror69 from FlexpartTools import EIFlexpart, Control, interpret_args_and_control, normalexit, myerror 42 70 from getMARSdata import getMARSdata 43 71 from prepareFLEXPART import prepareFLEXPART 72 # ------------------------------------------------------------------------------ 73 # FUNCTIONS 74 # ------------------------------------------------------------------------------ 75 def main(): 76 ''' 77 @Description: 78 Get the arguments from script call and initialize an object from 79 Control class. Decides from the argument "queue" if the local version 80 is done "queue=None" or the gateway version "queue=ecgate". 44 81 82 @Input: 83 <nothing> 45 84 46 47 def main(): 48 49 calledfromdir=os.getcwd() 50 # os.chdir(localpythonpath) 51 args,c=interpret_args_and_control() 52 if args.queue==None: 53 if c.inputdir[0]!='/': 54 c.inputdir=os.path.join(calledfromdir,c.inputdir) 55 if c.outputdir[0]!='/': 56 c.outputdir=os.path.join(calledfromdir,c.outputdir) 57 getMARSdata(args,c) 58 prepareFLEXPART(args,c) 85 @Return: 86 <nothing> 87 ''' 88 calledfromdir = os.getcwd() 89 args, c = interpret_args_and_control() 90 if args.queue is None: 91 if c.inputdir[0] != '/': 92 c.inputdir = os.path.join(calledfromdir, c.inputdir) 93 if c.outputdir[0] != '/': 94 c.outputdir = os.path.join(calledfromdir, c.outputdir) 95 getMARSdata(args, c) 96 prepareFLEXPART(args, c) 59 97 normalexit(c) 60 98 else: 61 submit(args.job_template,c,args.queue) 62 63 64 def submit(jtemplate,c,queue): 65 66 f=open(jtemplate) 67 lftext=f.read().split('\n') 68 insert_point=lftext.index('EOF') 69 f.close() 70 71 clist=c.tolist() 72 colist=[] 73 mt=0 99 submit(args.job_template, c, args.queue) 100 101 return 102 103 def submit(jtemplate, c, queue): 104 #AP divide in two submits , ondemand und operational 105 ''' 106 @Description: 107 Prepares the job script and submit it to the specified queue. 108 109 @Input: 110 jtemplate: string 111 Job template file for submission to ECMWF. It contains all necessary 112 module and variable settings for the ECMWF environment as well as 113 the job call and mail report instructions. 114 Default is "job.temp". 115 116 c: instance of class Control 117 Contains all necessary information of a control file. The parameters 118 are: DAY1, DAY2, DTIME, MAXSTEP, TYPE, TIME, STEP, CLASS, STREAM, 119 NUMBER, EXPVER, GRID, LEFT, LOWER, UPPER, RIGHT, LEVEL, LEVELIST, 120 RESOL, GAUSS, ACCURACY, OMEGA, OMEGADIFF, ETA, ETADIFF, DPDETA, 121 SMOOTH, FORMAT, ADDPAR, WRF, CWC, PREFIX, ECSTORAGE, ECTRANS, ECFSDIR, 122 MAILOPS, MAILFAIL, GRIB2FLEXPART, FLEXPARTDIR 123 For more information about format and content of the parameter see 124 documentation. 125 126 queue: string 127 Name of queue for submission to ECMWF (e.g. ecgate or cca ) 128 129 @Return: 130 <nothing> 131 ''' 132 133 # read template file and split from newline signs 134 with open(jtemplate) as f: 135 lftext = f.read().split('\n') 136 insert_point = lftext.index('EOF') 137 #AP es gibt mehrere EOFs überprüfen! 138 139 # put all parameters of control instance into a list 140 clist = c.tolist() # reanalysis (EI) 141 colist = [] # operational 142 mt = 0 143 144 #AP wieso 2 for loops? 145 #AP dieser part ist für die CTBTO Operational retrieves bis zum aktuellen Tag. 74 146 for elem in clist: 75 147 if 'maxstep' in elem: 76 mt =int(elem.split(' ')[1])77 148 mt = int(elem.split(' ')[1]) 149 78 150 for elem in clist: 79 151 if 'start_date' in elem: 80 elem ='start_date '+'${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}'152 elem = 'start_date ' + '${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}' 81 153 if 'end_date' in elem: 82 elem='start_date '+'${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}' 154 #AP Fehler?! Muss end_date heissen 155 elem = 'start_date ' + '${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}' 83 156 if 'base_time' in elem: 84 elem ='base_time '+'${MSJ_BASETIME}'85 if 'time' in elem and mt >24:86 elem ='time '+'${MSJ_BASETIME} {MSJ_BASETIME}'157 elem = 'base_time ' + '${MSJ_BASETIME}' 158 if 'time' in elem and mt > 24: 159 elem = 'time ' + '${MSJ_BASETIME} {MSJ_BASETIME}' 87 160 colist.append(elem) 88 89 lftextondemand=lftext[:insert_point]+clist+lftext[insert_point+2:] 90 lftextoper=lftext[:insert_point]+colist+lftext[insert_point+2:] 91 92 h=open('job.ksh','w') 93 h.write('\n'.join(lftextondemand)) 161 #AP end 162 163 #AP whats the difference between clist and colist ?! What is MSJ? 164 165 lftextondemand = lftext[:insert_point] + clist + lftext[insert_point + 2:] 166 lftextoper = lftext[:insert_point] + colist + lftext[insert_point + 2:] 167 168 with open('job.ksh', 'w') as h: 169 h.write('\n'.join(lftextondemand)) 170 171 #AP this is not used ?! what is it for? 172 #maybe a differentiation is needed 173 h = open('joboper.ksh', 'w') 174 h.write('\n'.join(lftextoper)) 94 175 h.close() 95 96 h=open('joboper.ksh','w') 97 h.write('\n'.join(lftextoper)) 98 h.close() 99 100 try: 101 p=subprocess.check_call(['ecaccess-job-submit','-queueName',queue,'job.ksh']) 176 #AP end 177 178 # submit job script to queue 179 try: 180 p = subprocess.check_call(['ecaccess-job-submit', '-queueName', 181 queue,' job.ksh']) 102 182 except: 103 print 'ecaccess-job-submit failed, probably eccert has expired'183 print('ecaccess-job-submit failed, probably eccert has expired') 104 184 exit(1) 105 #pout=p.communicate(input=s)[0]106 print 'You should get an email with subject flex.hostname.pid'107 185 186 print('You should get an email with subject flex.hostname.pid') 108 187 188 return 109 189 110 111 190 if __name__ == "__main__": 112 191 main() -
python/testsuite.py
rd69b677 r64cf353 2 2 3 3 # This software is licensed under the terms of the Apache Licence Version 2.0 4 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 5 # 4 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. 5 # 6 6 # Leopold Haimberger, Dec 2015 7 7 # 8 8 # Functionality provided: This script triggers the ECMWFDATA test suite. Call with 9 9 # testsuite.py [test group] 10 # 10 # 11 11 # 12 12 # Further documentation may be obtained from www.flexpart.eu 13 # 13 # 14 14 # Test groups are specified in testsuite.json 15 15 # in addition to a standard python 2.6 or 2.7 installation the following packages need to be installed … … 21 21 import json 22 22 import subprocess 23 23 24 try: 24 25 taskfile=open('testsuite.json') … … 41 42 os.makedirs('../test') 42 43 if len(sys.argv)>1: 43 groups=sys.argv[1:] 44 groups=sys.argv[1:] 44 45 else: 45 46 groups=['xinstall','default','ops','work','cv','fc']#,'hires'] … … 53 54 garglist=[] 54 55 for ttk,ttv in tv.iteritems(): 55 if isinstance(ttv,basestring): 56 if isinstance(ttv,basestring): 56 57 if ttk!='script': 57 58 garglist.append('--'+ttk) … … 64 65 arglist=[] 65 66 for tttk,tttv in ttv.iteritems(): 66 if isinstance(tttv,basestring): 67 if isinstance(tttv,basestring): 67 68 arglist.append('--'+tttk) 68 69 if '$' in tttv[0]: 69 arglist.append(os.path.expandvars(tttv)) 70 arglist.append(os.path.expandvars(tttv)) 70 71 else: 71 72 arglist.append(tttv) … … 86 87 print str(jobcounter-jobfailed)+' successful, '+str(jobfailed)+' failed' 87 88 print 'If tasks have been submitted via ECACCESS please check emails' 88 89 89 90 #print tasks
Note: See TracChangeset
for help on using the changeset viewer.