- Timestamp:
- Oct 5, 2018, 3:35:18 PM (6 years ago)
- Branches:
- master, ctbto, dev
- Children:
- 5bad6ec
- Parents:
- 27fe969
- Location:
- source/python
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
source/python/classes/ControlFile.py
r4971f63 rca867de 57 57 import inspect 58 58 59 # software specific classes and modules from flex_extract 60 sys.path.append('../') 59 61 import _config 62 from mods.tools import my_error 60 63 61 64 # ------------------------------------------------------------------------------ … … 144 147 self.ectrans = 0 145 148 self.inputdir = _config.PATH_INPUT_DIR 146 self.outputdir = self.inputdir149 self.outputdir = None 147 150 self.ecmwfdatadir = _config.PATH_FLEXEXTRACT_DIR 148 151 self.exedir = _config.PATH_FORTRAN_SRC … … 177 180 <nothing> 178 181 ''' 179 from mods.tools import my_error180 181 # read whole CONTROL file182 182 183 183 try: … … 355 355 self.end_date = self.start_date 356 356 357 # basetime has only two possible values 358 if self.basetime: 359 if int(self.basetime) != 0 and int(self.basetime) != 12: 360 print('Basetime has an invalid value!') 361 print('Basetime = ' + str(self.basetime)) 362 sys.exit(1) 363 357 364 # assure consistency of levelist and level 358 365 if self.levelist is None and self.level is None: … … 407 414 self.flexpart_root_scripts = self.ecmwfdatadir 408 415 416 if not self.outputdir: 417 self.outputdir = self.inputdir 418 409 419 if not isinstance(self.mailfail, list): 410 420 if ',' in self.mailfail: -
source/python/classes/EcFlexpart.py
r27fe969 rca867de 74 74 # MODULES 75 75 # ------------------------------------------------------------------------------ 76 import os 77 import sys 78 import glob 79 import shutil 76 80 import subprocess 77 import shutil78 import os79 import glob80 81 from datetime import datetime, timedelta 81 82 import numpy as np … … 85 86 86 87 # software specific classes and modules from flex_extract 88 sys.path.append('../') 87 89 import _config 88 90 from GribTools import GribTools 89 91 from mods.tools import init128, to_param_id, silent_remove, product, my_error 90 92 from MarsRetrieval import MarsRetrieval 91 import mods.disaggregation 93 import mods.disaggregation as disaggregation 92 94 93 95 # ------------------------------------------------------------------------------ … … 287 289 'SFC', '1', self.grid] 288 290 289 # if needed, add additional WRF specific parameters here290 291 291 return 292 292 … … 382 382 383 383 384 def _mk_index_values(self, inputdir, inputfiles, keys): 385 ''' 386 @Description: 387 Creates an index file for a set of grib parameter keys. 388 The values from the index keys are returned in a list. 389 390 @Input: 391 keys: dictionary 392 List of parameter names which serves as index. 393 394 inputfiles: instance of UioFiles 395 Contains a list of files. 396 397 @Return: 398 iid: grib_index 399 This is a grib specific index structure to access 400 messages in a file. 401 402 index_vals: list 403 Contains the values from the keys used for a distinct selection 404 of grib messages in processing the grib files. 405 Content looks like e.g.: 406 index_vals[0]: ('20171106', '20171107', '20171108') ; date 407 index_vals[1]: ('0', '1200', '1800', '600') ; time 408 index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange 409 ''' 410 iid = None 411 index_keys = keys 412 413 indexfile = os.path.join(inputdir, _config.FILE_GRIB_INDEX) 414 silent_remove(indexfile) 415 grib = GribTools(inputfiles.files) 416 # creates new index file 417 iid = grib.index(index_keys=index_keys, index_file=indexfile) 418 419 # read the values of index keys 420 index_vals = [] 421 for key in index_keys: 422 #index_vals.append(grib_index_get(iid, key)) 423 #print(index_vals[-1]) 424 key_vals = grib_index_get(iid, key) 425 print(key_vals) 426 # have to sort the steps for disaggregation, 427 # therefore convert to int first 428 if key == 'step': 429 key_vals = [int(k) for k in key_vals] 430 key_vals.sort() 431 key_vals = [str(k) for k in key_vals] 432 index_vals.append(key_vals) 433 # index_vals looks for example like: 434 # index_vals[0]: ('20171106', '20171107', '20171108') ; date 435 # index_vals[1]: ('0', '1200') ; time 436 # index_vals[2]: (3', '6', '9', '12') ; stepRange 437 438 return iid, index_vals 439 440 384 441 def retrieve(self, server, dates, request, inputdir='.'): 385 442 ''' … … 427 484 t24h = timedelta(hours=24) 428 485 429 # dictionary which contains all parameter for the mars request 486 # dictionary which contains all parameter for the mars request, 430 487 # entries with a "None" will change in different requests and will 431 488 # therefore be set in each request seperately … … 674 731 table128 = init128(_config.PATH_GRIBTABLE) 675 732 pars = to_param_id(self.params['OG_acc_SL'][0], table128) 733 734 iid = None 735 index_vals = None 736 737 # get the values of the keys which are used for distinct access 738 # of grib messages via product 676 739 index_keys = ["date", "time", "step"] 677 indexfile = os.path.join(c.inputdir, _config.FILE_GRIB_INDEX) 678 silent_remove(indexfile) 679 grib = GribTools(inputfiles.files) 680 # creates new index file 681 iid = grib.index(index_keys=index_keys, index_file=indexfile) 682 683 # read values of index keys 684 index_vals = [] 685 for key in index_keys: 686 key_vals = grib_index_get(iid, key) 687 print(key_vals) 688 # have to sort the steps for disaggregation, 689 # therefore convert to int first 690 if key == 'step': 691 key_vals = [int(k) for k in key_vals] 692 key_vals.sort() 693 key_vals = [str(k) for k in key_vals] 694 index_vals.append(key_vals) 695 # index_vals looks for example like: 696 # index_vals[0]: ('20171106', '20171107', '20171108') ; date 697 # index_vals[1]: ('0', '1200') ; time 698 # index_vals[2]: (3', '6', '9', '12') ; stepRange 740 iid, index_vals = self._mk_index_values(c.inputdir, 741 inputfiles, 742 index_keys) 743 # index_vals looks like e.g.: 744 # index_vals[0]: ('20171106', '20171107', '20171108') ; date 745 # index_vals[1]: ('0', '1200', '1800', '600') ; time 746 # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange 699 747 700 748 valsdict = {} 701 749 svalsdict = {} 702 stepsdict = {}750 # stepsdict = {} 703 751 for p in pars: 704 752 valsdict[str(p)] = [] 705 753 svalsdict[str(p)] = [] 706 stepsdict[str(p)] = []754 # stepsdict[str(p)] = [] 707 755 708 756 print('maxstep: ', c.maxstep) 709 757 758 # "product" genereates each possible combination between the 759 # values of the index keys 710 760 for prod in product(*index_vals): 711 761 # e.g. prod = ('20170505', '0', '12') 712 762 # ( date ,time, step) 713 # per date e.g. time = 0, 1200 714 # per time e.g. step = 3, 6, 9, 12715 print('current prod: ', prod) 763 764 print('current product: ', prod) 765 716 766 for i in range(len(index_keys)): 717 767 grib_index_select(iid, index_keys[i], prod[i]) … … 720 770 gid = grib_new_from_index(iid) 721 771 722 # if there is data for this product combination 723 # prepare some date and time parameter before reading the data 724 if gid is not None: 725 cdate = grib_get(gid, 'date') 726 time = grib_get(gid, 'time') 727 step = grib_get(gid, 'step') 728 # date+time+step-2*dtime 729 # (since interpolated value valid for step-2*dtime) 730 sdate = datetime(year=cdate/10000, 731 month=(cdate % 10000)/100, 732 day=(cdate % 100), 733 hour=time/100) 734 fdate = sdate + timedelta(hours=step-2*int(c.dtime)) 735 sdates = sdate + timedelta(hours=step) 736 elimit = None 737 else: 738 break 772 # if there is no data for this specific time combination / product 773 # skip the rest of the for loop and start with next timestep/product 774 if not gid: 775 continue 776 777 # create correct timestamp from the three time informations 778 cdate = str(grib_get(gid, 'date')) 779 ctime = '{:0>2}'.format(grib_get(gid, 'time')/100) 780 cstep = '{:0>3}'.format(grib_get(gid, 'step')) 781 t_date = datetime.strptime(cdate + ctime, '%Y%m%d%H') 782 t_dt = t_date + timedelta(hours=int(cstep)) 783 t_m1dt = t_date + timedelta(hours=int(cstep)-int(c.dtime)) 784 t_m2dt = t_date + timedelta(hours=int(cstep)-2*int(c.dtime)) 785 t_enddate = None 739 786 740 787 if c.maxstep > 12: 741 788 fnout = os.path.join(c.inputdir, 'flux' + 742 sdate.strftime('%Y%m%d') + 743 '.{:0>2}'.format(time/100) + 789 t_date.strftime('%Y%m%d.%H') + 744 790 '.{:0>3}'.format(step-2*int(c.dtime))) 745 791 gnout = os.path.join(c.inputdir, 'flux' + 746 sdate.strftime('%Y%m%d') + 747 '.{:0>2}'.format(time/100) + 792 t_date.strftime('%Y%m%d.%H') + 748 793 '.{:0>3}'.format(step-int(c.dtime))) 749 794 hnout = os.path.join(c.inputdir, 'flux' + 750 sdate.strftime('%Y%m%d') + 751 '.{:0>2}'.format(time/100) + 795 t_date.strftime('%Y%m%d.%H') + 752 796 '.{:0>3}'.format(step)) 753 797 else: 754 798 fnout = os.path.join(c.inputdir, 'flux' + 755 fdate.strftime('%Y%m%d%H'))799 t_m2dt.strftime('%Y%m%d%H')) 756 800 gnout = os.path.join(c.inputdir, 'flux' + 757 (fdate + timedelta(hours=int(c.dtime))). 758 strftime('%Y%m%d%H')) 801 t_m1dt.strftime('%Y%m%d%H')) 759 802 hnout = os.path.join(c.inputdir, 'flux' + 760 sdates.strftime('%Y%m%d%H'))803 t_dt.strftime('%Y%m%d%H')) 761 804 762 805 print("outputfile = " + fnout) 763 f_handle = open(fnout, 'w')764 g_handle = open(gnout, 'w')765 h_handle = open(hnout, 'w')766 806 767 807 # read message for message and store relevant data fields 768 808 # data keywords are stored in pars 769 809 while 1: 770 if gid is None:810 if not gid: 771 811 break 772 812 cparamId = str(grib_get(gid, 'paramId')) 773 813 step = grib_get(gid, 'step') 774 atime = grib_get(gid, 'time')814 time = grib_get(gid, 'time') 775 815 ni = grib_get(gid, 'Ni') 776 816 nj = grib_get(gid, 'Nj') … … 779 819 vdp = valsdict[cparamId] 780 820 svdp = svalsdict[cparamId] 781 sd = stepsdict[cparamId]821 # sd = stepsdict[cparamId] 782 822 783 823 if cparamId == '142' or cparamId == '143': … … 793 833 svdp.append((vdp[-1] - vdp[-2]) / int(c.dtime)) 794 834 795 print(cparamId, atime, step, len(values),835 print(cparamId, time, step, len(values), 796 836 values[0], np.std(values)) 797 # save the 1/3-hourly or specific values 798 # svdp.append(values[:]) 799 sd.append(step) 837 800 838 # len(svdp) correspond to the time 801 839 if len(svdp) >= 3: … … 807 845 808 846 if not (step == c.maxstep and c.maxstep > 12 \ 809 or sdates == elimit):847 or t_dt == t_enddate): 810 848 vdp.pop(0) 811 849 svdp.pop(0) … … 817 855 818 856 grib_set_values(gid, values) 857 819 858 if c.maxstep > 12: 820 859 grib_set(gid, 'step', max(0, step-2*int(c.dtime))) 821 860 else: 822 861 grib_set(gid, 'step', 0) 823 grib_set(gid, 'time', fdate.hour*100) 824 grib_set(gid, 'date', fdate.year*10000 + 825 fdate.month*100+fdate.day) 826 grib_write(gid, f_handle) 862 grib_set(gid, 'time', t_m2dt.hour*100) 863 grib_set(gid, 'date', int(t_m2dt.strftime('%Y%m%d'))) 864 865 with open(fnout, 'w') as f_handle: 866 grib_write(gid, f_handle) 827 867 828 868 if c.basetime: 829 elimit = datetime.strptime(c.end_date + 830 c.basetime, '%Y%m%d%H') 869 t_enddate = datetime.strptime(c.end_date + 870 c.basetime, 871 '%Y%m%d%H') 831 872 else: 832 elimit = sdate + timedelta(2*int(c.dtime))873 t_enddate = t_date + timedelta(2*int(c.dtime)) 833 874 834 875 # squeeze out information of last two steps contained 835 876 # in svdp 836 877 # if step+int(c.dtime) == c.maxstep and c.maxstep>12 837 # or sdates+timedelta(hours = int(c.dtime))838 # >= elimit:878 # or t_dt+timedelta(hours = int(c.dtime)) 879 # >= t_enddate: 839 880 # Note that svdp[0] has not been popped in this case 840 881 841 882 if step == c.maxstep and c.maxstep > 12 or \ 842 sdates == elimit:883 t_dt == t_enddate: 843 884 844 885 values = svdp[3] 845 886 grib_set_values(gid, values) 846 887 grib_set(gid, 'step', 0) 847 truedatetime = fdate+ timedelta(hours=888 truedatetime = t_m2dt + timedelta(hours= 848 889 2*int(c.dtime)) 849 890 grib_set(gid, 'time', truedatetime.hour * 100) … … 851 892 truedatetime.month * 100 + 852 893 truedatetime.day) 853 grib_write(gid, h_handle) 894 with open(hnout, 'w') as h_handle: 895 grib_write(gid, h_handle) 854 896 855 897 #values = (svdp[1]+svdp[2])/2. … … 860 902 861 903 grib_set(gid, 'step', 0) 862 truedatetime = fdate+ timedelta(hours=int(c.dtime))904 truedatetime = t_m2dt + timedelta(hours=int(c.dtime)) 863 905 grib_set(gid, 'time', truedatetime.hour * 100) 864 906 grib_set(gid, 'date', truedatetime.year * 10000 + … … 866 908 truedatetime.day) 867 909 grib_set_values(gid, values) 868 grib_write(gid, g_handle) 869 870 grib_release(gid) 871 872 gid = grib_new_from_index(iid) 873 874 f_handle.close() 875 g_handle.close() 876 h_handle.close() 910 with open(gnout, 'w') as g_handle: 911 grib_write(gid, g_handle) 912 913 grib_release(gid) 914 915 gid = grib_new_from_index(iid) 877 916 878 917 grib_index_release(iid) … … 913 952 ''' 914 953 915 table128 = init128(_config.PATH_GRIBTABLE) 916 wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\ 917 stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4', 918 table128) 919 954 if c.wrf: 955 table128 = init128(_config.PATH_GRIBTABLE) 956 wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\ 957 stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4', 958 table128) 959 960 # these numbers are indices for the temporary files "fort.xx" 961 # which are used to seperate the grib fields to, 962 # for the Fortran program input 963 # 10: U,V | 11: T | 12: lnsp | 13: D | 16: sfc fields 964 # 17: Q | 18: Q , gaussian| 19: w | 21: etadot | 22: clwc+ciwc 965 fdict = {'10':None, '11':None, '12':None, '13':None, '16':None, 966 '17':None, '18':None, '19':None, '21':None, '22':None} 967 968 iid = None 969 index_vals = None 970 971 # get the values of the keys which are used for distinct access 972 # of grib messages via product 920 973 index_keys = ["date", "time", "step"] 921 indexfile = os.path.join(c.inputdir, _config.FILE_GRIB_INDEX) 922 silent_remove(indexfile) 923 grib = GribTools(inputfiles.files) 924 # creates new index file 925 iid = grib.index(index_keys=index_keys, index_file=indexfile) 926 927 # read values of index keys 928 index_vals = [] 929 for key in index_keys: 930 index_vals.append(grib_index_get(iid, key)) 931 print(index_vals[-1]) 932 # index_vals looks for example like: 933 # index_vals[0]: ('20171106', '20171107', '20171108') ; date 934 # index_vals[1]: ('0', '1200', '1800', '600') ; time 935 # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange 936 937 fdict = {'10':None, '11':None, '12':None, '13':None, '16':None, 938 '17':None, '19':None, '21':None, '22':None, '20':None} 939 974 iid, index_vals = self._mk_index_values(c.inputdir, 975 inputfiles, 976 index_keys) 977 # index_vals looks like e.g.: 978 # index_vals[0]: ('20171106', '20171107', '20171108') ; date 979 # index_vals[1]: ('0', '1200', '1800', '600') ; time 980 # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange 981 982 # "product" genereates each possible combination between the 983 # values of the index keys 940 984 for prod in product(*index_vals): 941 985 # e.g. prod = ('20170505', '0', '12') 942 986 # ( date ,time, step) 943 # per date e.g. time = 0, 1200 944 # per time e.g. step = 3, 6, 9, 12 945 # flag for Fortran program and file merging 946 convertFlag = False 947 print('current prod: ', prod) 948 # e.g. prod = ('20170505', '0', '12') 949 # ( date ,time, step) 950 # per date e.g. time = 0, 600, 1200, 1800 951 # per time e.g. step = 0, 3, 6, 9, 12 987 988 print('current product: ', prod) 989 952 990 for i in range(len(index_keys)): 953 991 grib_index_select(iid, index_keys[i], prod[i]) … … 956 994 gid = grib_new_from_index(iid) 957 995 958 # if there is data for this product combination 959 # prepare some date and time parameter before reading the data 960 if gid is not None: 961 # Combine all temporary data files into final grib file if 962 # gid is at least one time not None. Therefore set convertFlag 963 # to save information. The Fortran program is also 964 # only executed if convertFlag is True 965 convertFlag = True 966 # remove old fort.* files and open new ones 967 # they are just valid for a single product 968 for k, f in fdict.iteritems(): 969 fortfile = os.path.join(c.inputdir, 'fort.' + k) 970 silent_remove(fortfile) 971 fdict[k] = open(fortfile, 'w') 972 973 cdate = str(grib_get(gid, 'date')) 974 time = grib_get(gid, 'time') 975 step = grib_get(gid, 'step') 976 # create correct timestamp from the three time informations 977 # date, time, step 978 timestamp = datetime.strptime(cdate + '{:0>2}'.format(time/100), 979 '%Y%m%d%H') 980 timestamp += timedelta(hours=int(step)) 981 cdateH = datetime.strftime(timestamp, '%Y%m%d%H') 982 983 if c.basetime: 984 slimit = datetime.strptime(c.start_date + '00', '%Y%m%d%H') 985 bt = '23' 986 if c.basetime == '00': 987 bt = '00' 988 slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\ 989 - timedelta(hours=12-int(c.dtime)) 990 if c.basetime == '12': 991 bt = '12' 992 slimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H')\ 993 - timedelta(hours=12-int(c.dtime)) 994 995 elimit = datetime.strptime(c.end_date + bt, '%Y%m%d%H') 996 997 if timestamp < slimit or timestamp > elimit: 998 continue 999 else: 1000 pass 1001 1002 try: 1003 if c.wrf: 1004 if 'olddate' not in locals() or cdate != olddate: 1005 fwrf = open(os.path.join(c.outputdir, 1006 'WRF' + cdate + '.{:0>2}'.format(time) + 1007 '.000.grb2'), 'w') 1008 olddate = cdate[:] 1009 except AttributeError: 1010 pass 1011 1012 # helper variable to remember which fields were already used. 996 # if there is no data for this specific time combination / product 997 # skip the rest of the for loop and start with next timestep/product 998 if not gid: 999 continue 1000 1001 # remove old fort.* files and open new ones 1002 # they are just valid for a single product 1003 for k, f in fdict.iteritems(): 1004 fortfile = os.path.join(c.inputdir, 'fort.' + k) 1005 silent_remove(fortfile) 1006 fdict[k] = open(fortfile, 'w') 1007 1008 # create correct timestamp from the three time informations 1009 cdate = str(grib_get(gid, 'date')) 1010 ctime = '{:0>2}'.format(grib_get(gid, 'time')/100) 1011 cstep = '{:0>3}'.format(grib_get(gid, 'step')) 1012 timestamp = datetime.strptime(cdate + ctime, '%Y%m%d%H') 1013 timestamp += timedelta(hours=int(cstep)) 1014 cdate_hour = datetime.strftime(timestamp, '%Y%m%d%H') 1015 1016 # if the timestamp is out of basetime start/end date period, 1017 # skip this specific product 1018 if c.basetime: 1019 start_time = datetime.strptime(c.end_date + c.basetime, 1020 '%Y%m%d%H') - time_delta 1021 end_time = datetime.strptime(c.end_date + c.basetime, 1022 '%Y%m%d%H') 1023 if timestamp < start_time or timestamp > end_time: 1024 continue 1025 1026 if c.wrf: 1027 if 'olddate' not in locals() or cdate != olddate: 1028 fwrf = open(os.path.join(c.outputdir, 1029 'WRF' + cdate + '.' + ctime + '.000.grb2'), 'w') 1030 olddate = cdate[:] 1031 1032 # savedfields remembers which fields were already used. 1013 1033 savedfields = [] 1034 # sum of cloud liquid and ice water content 1035 scwc = None 1014 1036 while 1: 1015 if gid is None:1037 if not gid: 1016 1038 break 1017 1039 paramId = grib_get(gid, 'paramId') 1018 1040 gridtype = grib_get(gid, 'gridType') 1019 1041 levtype = grib_get(gid, 'typeOfLevel') 1020 if paramId == 133 and gridtype == 'reduced_gg': 1021 # Specific humidity (Q.grb) is used as a template only 1022 # so we need the first we "meet" 1023 with open(os.path.join(c.inputdir, 'fort.18'), 'w') as fout: 1024 grib_write(gid, fout) 1025 elif paramId == 131 or paramId == 132: 1042 if paramId == 77: # ETADOT 1043 grib_write(gid, fdict['21']) 1044 elif paramId == 130: # T 1045 grib_write(gid, fdict['11']) 1046 elif paramId == 131 or paramId == 132: # U, V wind component 1026 1047 grib_write(gid, fdict['10']) 1027 elif paramId == 130: 1028 grib_write(gid, fdict['11']) 1029 elif paramId == 133 and gridtype != 'reduced_gg': 1048 elif paramId == 133 and gridtype != 'reduced_gg': # Q 1030 1049 grib_write(gid, fdict['17']) 1031 elif paramId == 152: 1050 elif paramId == 133 and gridtype == 'reduced_gg': # Q, gaussian 1051 grib_write(gid, fdict['18']) 1052 elif paramId == 135: # W 1053 grib_write(gid, fdict['19']) 1054 elif paramId == 152: # LNSP 1032 1055 grib_write(gid, fdict['12']) 1033 elif paramId == 155 and gridtype == 'sh': 1056 elif paramId == 155 and gridtype == 'sh': # D 1034 1057 grib_write(gid, fdict['13']) 1035 elif paramId in [129, 138, 155] and levtype == 'hybrid' \ 1036 and c.wrf: 1037 pass 1038 elif paramId == 246 or paramId == 247: 1039 # cloud liquid water and ice 1040 if paramId == 246: 1041 clwc = grib_get_values(gid) 1058 elif paramId == 246 or paramId == 247: # CLWC, CIWC 1059 # sum cloud liquid water and ice 1060 if not scwc: 1061 scwc = grib_get_values(gid) 1042 1062 else: 1043 clwc += grib_get_values(gid)1044 grib_set_values(gid, clwc)1063 scwc += grib_get_values(gid) 1064 grib_set_values(gid, scwc) 1045 1065 grib_set(gid, 'paramId', 201031) 1046 1066 grib_write(gid, fdict['22']) 1047 elif paramId == 135: 1048 grib_write(gid, fdict['19']) 1049 elif paramId == 77: 1050 grib_write(gid, fdict['21']) 1067 elif c.wrf and paramId in [129, 138, 155] and \ 1068 levtype == 'hybrid': # Z, VO, D 1069 # do not do anything right now 1070 # these are specific parameter for WRF 1071 pass 1051 1072 else: 1052 1073 if paramId not in savedfields: 1074 # SD/MSL/TCC/10U/10V/2T/2D/Z/LSM/SDOR/CVL/CVH/SR 1075 # and all ADDPAR parameter 1053 1076 grib_write(gid, fdict['16']) 1054 1077 savedfields.append(paramId) … … 1074 1097 f.close() 1075 1098 1076 # call for Fortran program if flag is True 1077 if convertFlag: 1078 pwd = os.getcwd() 1079 os.chdir(c.inputdir) 1080 if os.stat('fort.21').st_size == 0 and c.eta: 1081 print('Parameter 77 (etadot) is missing, most likely it is \ 1082 not available for this type or date/time\n') 1083 print('Check parameters CLASS, TYPE, STREAM, START_DATE\n') 1084 my_error(c.mailfail, 'fort.21 is empty while parameter eta \ 1085 is set to 1 in CONTROL file') 1086 1087 # create the corresponding output file fort.15 1088 # (generated by Fortran program) + fort.16 (paramId 167 and 168) 1089 p = subprocess.check_call([os.path.join( 1090 c.exedir, _config.FORTRAN_EXECUTABLE)], shell=True) 1091 os.chdir(pwd) 1092 1093 # create final output filename, e.g. EN13040500 (ENYYMMDDHH) 1094 fnout = os.path.join(c.inputdir, c.prefix) 1095 if c.maxstep > 12: 1096 suffix = cdate[2:8] + '.{:0>2}'.format(time/100) + \ 1097 '.{:0>3}'.format(step) 1098 else: 1099 suffix = cdateH[2:10] 1100 fnout += suffix 1101 print("outputfile = " + fnout) 1102 self.outputfilelist.append(fnout) # needed for final processing 1103 1104 # create outputfile and copy all data from intermediate files 1105 # to the outputfile (final GRIB files) 1106 orolsm = os.path.basename(glob.glob( 1107 c.inputdir + '/OG_OROLSM__SL.*.' + c.ppid + '*')[0]) 1108 fluxfile = 'flux' + cdate[0:2] + suffix 1109 if not c.cwc: 1110 flist = ['fort.15', fluxfile, 'fort.16', orolsm] 1111 else: 1112 flist = ['fort.15', 'fort.22', fluxfile, 'fort.16', orolsm] 1113 1114 with open(fnout, 'wb') as fout: 1115 for f in flist: 1116 shutil.copyfileobj(open(os.path.join(c.inputdir, f), 1117 'rb'), fout) 1118 1119 if c.omega: 1120 with open(os.path.join(c.outputdir, 'OMEGA'), 'wb') as fout: 1121 shutil.copyfileobj( 1122 open(os.path.join(c.inputdir, 'fort.25'), 1123 'rb'), fout) 1099 # call for Fortran program to convert e.g. reduced_gg grids to 1100 # regular_ll and calculate detadot/dp 1101 pwd = os.getcwd() 1102 os.chdir(c.inputdir) 1103 if os.stat('fort.21').st_size == 0 and c.eta: 1104 print('Parameter 77 (etadot) is missing, most likely it is \ 1105 not available for this type or date/time\n') 1106 print('Check parameters CLASS, TYPE, STREAM, START_DATE\n') 1107 my_error(c.mailfail, 'fort.21 is empty while parameter eta \ 1108 is set to 1 in CONTROL file') 1109 1110 # Fortran program creates file fort.15 (with u,v,etadot,t,sp,q) 1111 p = subprocess.check_call([os.path.join( 1112 c.exedir, _config.FORTRAN_EXECUTABLE)], shell=True) 1113 os.chdir(pwd) 1114 1115 # create name of final output file, e.g. EN13040500 (ENYYMMDDHH) 1116 if c.maxstep > 12: 1117 suffix = cdate[2:8] + '.' + ctime + '.' + cstep 1124 1118 else: 1125 pass 1119 suffix = cdate_hour[2:10] 1120 fnout = os.path.join(c.inputdir, c.prefix + suffix) 1121 print("outputfile = " + fnout) 1122 # collect for final processing 1123 self.outputfilelist.append(os.path.basename(fnout)) 1124 1125 # create outputfile and copy all data from intermediate files 1126 # to the outputfile (final GRIB input files for FLEXPART) 1127 orolsm = os.path.basename(glob.glob(c.inputdir + 1128 '/OG_OROLSM__SL.*.' + c.ppid + '*')[0]) 1129 fluxfile = 'flux' + cdate[0:2] + suffix 1130 if not c.cwc: 1131 flist = ['fort.15', fluxfile, 'fort.16', orolsm] 1132 else: 1133 flist = ['fort.15', 'fort.22', fluxfile, 'fort.16', orolsm] 1134 1135 with open(fnout, 'wb') as fout: 1136 for f in flist: 1137 shutil.copyfileobj(open(os.path.join(c.inputdir, f), 'rb'), 1138 fout) 1139 1140 if c.omega: 1141 with open(os.path.join(c.outputdir, 'OMEGA'), 'wb') as fout: 1142 shutil.copyfileobj(open(os.path.join(c.inputdir, 'fort.25'), 1143 'rb'), fout) 1126 1144 1127 1145 if c.wrf: … … 1168 1186 .format(c.ectrans, c.gateway, c.destination)) 1169 1187 1170 print('Output filelist: \n')1188 print('Output filelist: ') 1171 1189 print(self.outputfilelist) 1172 1190 … … 1191 1209 if c.outputdir != c.inputdir: 1192 1210 for ofile in self.outputfilelist: 1193 p = subprocess.check_call(['mv', ofile, c.outputdir]) 1211 p = subprocess.check_call(['mv', 1212 os.path.join(c.inputdir, ofile), 1213 c.outputdir]) 1194 1214 1195 1215 return -
source/python/classes/MarsRetrieval.py
r295ff45 rca867de 62 62 # MODULES 63 63 # ------------------------------------------------------------------------------ 64 import os 65 import sys 64 66 import subprocess 65 import os 66 67 68 # software specific classes and modules from flex_extract 69 sys.path.append('../') 67 70 import _config 68 71 # ------------------------------------------------------------------------------ -
source/python/classes/UioFiles.py
r25b14be rca867de 48 48 # ------------------------------------------------------------------------------ 49 49 import os 50 import sys 50 51 import fnmatch 51 52 52 53 # software specific module from flex_extract 54 sys.path.append('../') 53 55 #import profiling 54 56 from mods.tools import silent_remove, get_list_as_string -
source/python/mods/get_mars_data.py
r27fe969 rca867de 48 48 import os 49 49 import sys 50 import datetime51 50 import inspect 51 from datetime import datetime, timedelta 52 53 # software specific classes and modules from flex_extract 54 sys.path.append('../') 55 import _config 56 from tools import my_error, normal_exit, get_cmdline_arguments, read_ecenv 57 from classes.EcFlexpart import EcFlexpart 58 from classes.UioFiles import UioFiles 59 52 60 try: 53 61 ecapi = True … … 55 63 except ImportError: 56 64 ecapi = False 57 58 # software specific classes and modules from flex_extract59 import _config60 from tools import my_error, normal_exit, get_cmdline_arguments, read_ecenv61 from classes.EcFlexpart import EcFlexpart62 from classes.UioFiles import UioFiles63 65 # ------------------------------------------------------------------------------ 64 66 # FUNCTION … … 147 149 # allerdings ist das relevant und ersichtlich an den NICHT FLUSS DATEN 148 150 149 150 # set start date of retrieval period 151 start = datetime.date(year=int(c.start_date[:4]), 152 month=int(c.start_date[4:6]), 153 day=int(c.start_date[6:])) 154 startm1 = start - datetime.timedelta(days=1) 155 156 # set end date of retrieval period 157 end = datetime.date(year=int(c.end_date[:4]), 158 month=int(c.end_date[4:6]), 159 day=int(c.end_date[6:])) 160 161 # set time period for one single retrieval 162 datechunk = datetime.timedelta(days=int(c.date_chunk)) 151 start = datetime.strptime(c.start_date, '%Y%m%d') 152 end = datetime.strptime(c.end_date, '%Y%m%d') 153 # time period for one single retrieval 154 datechunk = timedelta(days=int(c.date_chunk)) 163 155 164 156 if c.basetime == '00': 165 start = startm1 157 start = start - timedelta(days=1) 158 159 if c.maxstep <= 24: 160 startm1 = start - timedelta(days=1) 166 161 167 162 if c.basetime == '00' or c.basetime == '12': 168 # endp1 = end + datetime.timedelta(days=1)163 # endp1 = end + timedelta(days=1) 169 164 endp1 = end 170 165 else: 171 # endp1 = end + datetime.timedelta(days=2)172 endp1 = end + datetime.timedelta(days=1)166 # endp1 = end + timedelta(days=2) 167 endp1 = end + timedelta(days=1) 173 168 174 169 # -------------- flux data ------------------------------------------------ … … 243 238 # we only need to add datechunk - 1 days to retrieval 244 239 # for a period 245 delta_t_m1 = delta_t - datetime.timedelta(days=1)240 delta_t_m1 = delta_t - timedelta(days=1) 246 241 247 242 day = start -
source/python/mods/prepare_flexpart.py
r27fe969 rca867de 59 59 60 60 # software specific classes and modules from flex_extract 61 sys.path.append('../') 61 62 import _config 62 63 from classes.UioFiles import UioFiles 64 from classes.ControlFile import ControlFile 63 65 from tools import clean_up, get_cmdline_arguments, read_ecenv 64 66 from classes.EcFlexpart import EcFlexpart -
source/python/mods/tools.py
r27fe969 rca867de 60 60 # ------------------------------------------------------------------------------ 61 61 62 def none_or_str(value): 63 ''' 64 @Description: 65 Converts the input string into pythons None-type if the string 66 contains "None". 67 68 @Input: 69 value: string 70 String to be checked for the "None" word. 71 72 @Return: 73 None or value: 74 Return depends on the content of the input value. If it was "None", 75 then the python type None is returned. Otherwise the string itself. 76 ''' 77 if value == 'None': 78 return None 79 return value 80 81 def none_or_int(value): 82 ''' 83 @Description: 84 Converts the input string into pythons None-type if the string 85 contains "None". Otherwise it is converted to an integer value. 86 87 @Input: 88 value: string 89 String to be checked for the "None" word. 90 91 @Return: 92 None or int(value): 93 Return depends on the content of the input value. If it was "None", 94 then the python type None is returned. Otherwise the string is 95 converted into an integer value. 96 ''' 97 if value == 'None': 98 return None 99 return int(value) 100 62 101 def get_cmdline_arguments(): 63 102 ''' … … 79 118 80 119 # the most important arguments 81 parser.add_argument("--start_date", dest="start_date", default=None, 120 parser.add_argument("--start_date", dest="start_date", 121 type=none_or_str, default=None, 82 122 help="start date YYYYMMDD") 83 parser.add_argument("--end_date", dest="end_date", default=None, 123 parser.add_argument("--end_date", dest="end_date", 124 type=none_or_str, default=None, 84 125 help="end_date YYYYMMDD") 85 parser.add_argument("--date_chunk", dest="date_chunk", default=None, 126 parser.add_argument("--date_chunk", dest="date_chunk", 127 type=none_or_int, default=None, 86 128 help="# of days to be retrieved at once") 129 parser.add_argument("--controlfile", dest="controlfile", 130 type=none_or_str, default='CONTROL.temp', 131 help="file with CONTROL parameters") 132 133 # parameter for extra output information 134 parser.add_argument("--debug", dest="debug", 135 type=none_or_int, default=None, 136 help="debug mode - leave temporary files intact") 137 parser.add_argument("--request", dest="request", 138 type=none_or_int, default=None, 139 help="list all mars request in file mars_requests.dat \ 140 and skip submission to mars") 87 141 88 142 # some arguments that override the default in the CONTROL file 89 parser.add_argument("--basetime", dest="basetime", default=None, 90 help="base such as 00/12 (for half day retrievals)") 91 parser.add_argument("--step", dest="step", default=None, 143 parser.add_argument("--basetime", dest="basetime", 144 type=none_or_int, default=None, 145 help="base such as 00 or 12 (for half day retrievals)") 146 parser.add_argument("--step", dest="step", 147 type=none_or_str, default=None, 92 148 help="steps such as 00/to/48") 93 parser.add_argument("--levelist", dest="levelist", default=None, 149 parser.add_argument("--levelist", dest="levelist", 150 type=none_or_str, default=None, 94 151 help="Vertical levels to be retrieved, e.g. 30/to/60") 95 parser.add_argument("--area", dest="area", default=None, 152 parser.add_argument("--area", dest="area", 153 type=none_or_str, default=None, 96 154 help="area defined as north/west/south/east") 97 155 98 156 # set the working directories 99 parser.add_argument("--inputdir", dest="inputdir", default=None, 157 parser.add_argument("--inputdir", dest="inputdir", 158 type=none_or_str, default=None, 100 159 help="root directory for storing intermediate files") 101 parser.add_argument("--outputdir", dest="outputdir", default=None, 160 parser.add_argument("--outputdir", dest="outputdir", 161 type=none_or_str, default=None, 102 162 help="root directory for storing output files") 103 163 parser.add_argument("--flexpart_root_scripts", dest="flexpart_root_scripts", 104 default=None,164 type=none_or_str, default=None, 105 165 help="FLEXPART root directory (to find grib2flexpart \ 106 166 and COMMAND file)\n Normally flex_extract resides in \ … … 108 168 109 169 # this is only used by prepare_flexpart.py to rerun a postprocessing step 110 parser.add_argument("--ppid", dest="ppid", default=None, 170 parser.add_argument("--ppid", dest="ppid", 171 type=none_or_int, default=None, 111 172 help="specify parent process id for \ 112 173 rerun of prepare_flexpart") … … 114 175 # arguments for job submission to ECMWF, only needed by submit.py 115 176 parser.add_argument("--job_template", dest='job_template', 116 default="job.temp",177 type=none_or_str, default="job.temp", 117 178 help="job template file for submission to ECMWF") 118 parser.add_argument("--queue", dest="queue", default=None, 179 parser.add_argument("--queue", dest="queue", 180 type=none_or_str, default=None, 119 181 help="queue for submission to ECMWF \ 120 182 (e.g. ecgate or cca )") 121 parser.add_argument("--controlfile", dest="controlfile",122 default='CONTROL.temp',123 help="file with CONTROL parameters")124 parser.add_argument("--debug", dest="debug", default=None,125 help="debug mode - leave temporary files intact")126 parser.add_argument("--request", dest="request", default=None,127 help="list all mars request in file mars_requests.dat \128 and skip submission to mars")129 183 130 184 args = parser.parse_args()
Note: See TracChangeset
for help on using the changeset viewer.