From 0e576fcaf0dab8e7e586efa1a4b9b194ee90f8bb Mon Sep 17 00:00:00 2001
From: Anne Philipp <anne.philipp@univie.ac.at>
Date: Mon, 8 Oct 2018 17:44:58 +0200
Subject: [PATCH] implemented extraction possibility of EA5 und CERA

---
 run/control/CONTROL_EA5                | 40 +++++++++++++++++
 run/control/CONTROL_EA5.highres        | 37 ++++++++++++++++
 run/control/CONTROL_EA5.public         | 40 +++++++++++++++++
 run/run.sh                             |  8 ++--
 source/python/classes/ControlFile.py   | 61 +++++++++++++++++++++++---
 source/python/classes/EcFlexpart.py    | 35 +++++++++------
 source/python/classes/MarsRetrieval.py | 23 +++++-----
 7 files changed, 209 insertions(+), 35 deletions(-)
 create mode 100644 run/control/CONTROL_EA5
 create mode 100644 run/control/CONTROL_EA5.highres
 create mode 100644 run/control/CONTROL_EA5.public

diff --git a/run/control/CONTROL_EA5 b/run/control/CONTROL_EA5
new file mode 100644
index 0000000..4d8f41f
--- /dev/null
+++ b/run/control/CONTROL_EA5
@@ -0,0 +1,40 @@
+DAY1
+DAY2
+DTIME 1
+TYPE AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN
+TIME 00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23
+STEP 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
+ACCTYPE FC
+ACCTIME 06/18
+ACCMAXSTEP 12
+CLASS EA
+STREAM OPER
+NUMBER OFF
+EXPVER 1
+GRID 1000  
+LEFT -25000
+LOWER 10000
+UPPER 75000
+RIGHT 60000
+LEVEL 137
+LEVELIST 130/to/137
+RESOL 159
+GAUSS 0
+ACCURACY 16
+OMEGA 0
+OMEGADIFF 0
+ETA 1
+ETADIFF 0
+DPDETA 1
+SMOOTH 0
+FORMAT GRIB1
+ADDPAR 186/187/188/235/139/39
+PREFIX EA
+ECSTORAGE 0
+ECTRANS 1
+ECFSDIR ectmp:/${USER}/econdemand/
+MAILFAIL ${USER} 
+MAILOPS ${USER}
+GRIB2FLEXPART 0
+EOF
+
diff --git a/run/control/CONTROL_EA5.highres b/run/control/CONTROL_EA5.highres
new file mode 100644
index 0000000..4ea88fa
--- /dev/null
+++ b/run/control/CONTROL_EA5.highres
@@ -0,0 +1,37 @@
+DAY1
+DAY2
+DTIME 1
+TYPE AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN
+TIME 00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23
+STEP 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
+ACCTYPE FC
+ACCTIME 06/18
+ACCMAXSTEP 12
+M_CLASS EA
+M_STREAM OPER
+M_NUMBER OFF
+M_EXPVER 1
+M_GRID 250  
+M_LEFT -100000
+M_LOWER 00000
+M_UPPER 60000
+M_RIGHT 60000
+M_LEVELIST 60/to/137
+M_RESOL 799
+M_GAUSS 0
+M_ACCURACY 24
+M_OMEGA 0
+M_OMEGADIFF 0
+M_ETA 1
+M_ETADIFF 0
+M_DPDETA 1
+M_SMOOTH 0
+M_FORMAT GRIB2
+M_ADDPAR /186/187/188/235/139/39
+PREFIX EA
+ECSTORAGE 0
+ECTRANS 1
+ECFSDIR ectmp:/${USER}/econdemand/
+MAILOPS ${USER}
+MAILFAIL ${USER}
+EOF
diff --git a/run/control/CONTROL_EA5.public b/run/control/CONTROL_EA5.public
new file mode 100644
index 0000000..40469e9
--- /dev/null
+++ b/run/control/CONTROL_EA5.public
@@ -0,0 +1,40 @@
+DAY1 
+DAY2 
+DTIME 1
+TYPE AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN AN
+TIME 00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23
+STEP 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00
+ACCTYPE FC
+ACCTIME 06/18
+ACCMAXSTEP 12
+CLASS EA
+DATASET ERA5
+STREAM OPER
+NUMBER OFF
+EXPVER 1
+GRID 1000 
+LEFT -15000
+LOWER 30000
+UPPER 75000
+RIGHT 45000
+LEVELIST 1/to/137
+RESOL 213
+ACCURACY 24
+GAUSS 0
+OMEGA 0
+OMEGADIFF 0
+ETA 1
+ETADIFF 0
+DPDETA 1
+SMOOTH 0
+FORMAT GRIB1
+ADDPAR 186/187/188/235/139/39
+PREFIX EApub  
+ECSTORAGE 0
+ECTRANS 1
+ECFSDIR ectmp:/${USER}/econdemand/
+MAILFAIL ${USER} 
+MAILOPS ${USER}
+GRIB2FLEXPART 0
+EOF
+
diff --git a/run/run.sh b/run/run.sh
index 91ee484..65acf26 100755
--- a/run/run.sh
+++ b/run/run.sh
@@ -22,15 +22,15 @@ BASETIME=None
 STEP=None
 LEVELIST=None
 AREA=None
-INPUTDIR='/nas/tmc/Anne/Interpolation/flexextract/flex_extract_v7.1/run/workspace/test'
+INPUTDIR='/nas/tmc/Anne/Interpolation/flexextract/flex_extract_v7.1/run/workspace/test3'
 OUTPUTDIR=None
 FLEXPART_ROOT_SCRIPTS=None 
 PP_ID=None
 JOB_TEMPLATE='job.temp' 
-CONTROLFILE='CONTROL_EI.public' 
+CONTROLFILE='CONTROL_EA5' 
 DEBUG=1 
-REQUEST=1
-PUBLIC=1
+REQUEST=2
+PUBLIC=0
 
 # -----------------------------------------------------------------
 #
diff --git a/source/python/classes/ControlFile.py b/source/python/classes/ControlFile.py
index 25ea969..5f9be59 100644
--- a/source/python/classes/ControlFile.py
+++ b/source/python/classes/ControlFile.py
@@ -112,11 +112,15 @@ class ControlFile(object):
         self.type = None
         self.time = None
         self.step = None
+        self.acctype = None
+        self.acctime = None
+        self.accmaxstep = None
         self.marsclass = None
         self.dataset = None
         self.stream = None
         self.number = 'OFF'
         self.expver = '1'
+        self.gaussian = ''
         self.grid = None
         self.area = ''
         self.left = None
@@ -345,7 +349,7 @@ class ControlFile(object):
 
         # check for having at least a starting date
         # otherwise program is not allowed to run
-        if self.start_date is None:
+        if not self.start_date:
             print('start_date specified neither in command line nor \
                    in CONTROL file ' +  self.controlfile)
             print('Try "' + sys.argv[0].split('/')[-1] +
@@ -353,7 +357,7 @@ class ControlFile(object):
             sys.exit(1)
 
         # retrieve just one day if end_date isn't set
-        if self.end_date is None:
+        if not self.end_date:
             self.end_date = self.start_date
 
         # basetime has only two possible values
@@ -364,18 +368,29 @@ class ControlFile(object):
                 sys.exit(1)
 
         # assure consistency of levelist and level
-        if self.levelist is None and self.level is None:
+        # up-to-date available maximum level numbers at ECMWF, 05.10.2018
+        max_level_list = [16, 19, 31, 40, 50, 60, 62, 91, 137]
+        if not self.levelist and not self.level:
             print('Warning: neither levelist nor level \
                                specified in CONTROL file')
             sys.exit(1)
-        elif self.levelist is None and self.level:
+        elif not self.levelist and self.level:
             self.levelist = '1/to/' + self.level
-        elif (self.levelist and self.level is None) or \
+        elif (self.levelist and not self.level) or \
              (self.levelist[-1] != self.level[-1]):
             self.level = self.levelist.split('/')[-1]
         else:
             pass
 
+        # check if max level is a valid level
+        if int(self.level) not in max_level_list:
+            print('ERROR: ')
+            print('LEVEL must be the maximum level of a specified '
+                  'level list from ECMWF, e.g.')
+            print('[16, 19, 31, 40, 50, 60, 62, 91 or 137]')
+            print('Check parameter "LEVEL" or the max level of "LEVELIST"!')
+            sys.exit(1)
+
         # if area was provided (only from commandline)
         # decompose area into its 4 components
         if self.area:
@@ -403,7 +418,7 @@ class ControlFile(object):
 
         # if maxstep wasn't provided
         # search for it in the "step" parameter
-        if self.maxstep is None:
+        if not self.maxstep:
             self.maxstep = 0
             for s in self.step:
                 if int(s) > self.maxstep:
@@ -461,6 +476,40 @@ class ControlFile(object):
                   'the "dataset"-parameter has to be set in the control file!')
             sys.exit(1)
 
+        if not isinstance(self.type, list):
+            self.type = [self.type]
+
+        for i, val in enumerate(self.type):
+            if self.type[i] == 'AN' and int(self.step[i]) != 0:
+                print('Analysis retrievals must have STEP = 0 (is set to 0)')
+                self.type[i] = 0
+
+        if not isinstance(self.time, list):
+            self.time = [self.time]
+
+        if not isinstance(self.step, list):
+            self.step = [self.step]
+
+        if not self.acctype:
+            print('... Control paramter ACCTYPE was not defined.')
+            try:
+                if len(self.type) > 1 and self.type[1] != 'AN':
+                    print('Use old setting by using TYPE[1] for flux forecast!')
+                    self.acctype = self.type[1]
+            except:
+                print('Use default value "FC" for flux forecast!')
+                self.acctype='FC'
+
+        if not self.acctime:
+            print('... Control paramter ACCTIME was not defined.')
+            print('Use default value "00/12" for flux forecast!')
+            self.acctime='00/12'
+
+        if not self.accmaxstep:
+            print('... Control paramter ACCMAXSTEP was not defined.')
+            print('Use default value "12" for flux forecast!')
+            self.accmaxstep='12'
+
         return
 
     def check_install_conditions(self):
diff --git a/source/python/classes/EcFlexpart.py b/source/python/classes/EcFlexpart.py
index 8ce3f35..eece37d 100644
--- a/source/python/classes/EcFlexpart.py
+++ b/source/python/classes/EcFlexpart.py
@@ -134,7 +134,7 @@ class EcFlexpart(object):
         self.types = dict()
 
         if c.maxstep > len(c.type):    # Pure forecast mode
-            c.type = [c.type[1]]
+            c.type = [c.type[0]]
             c.step = ['{:0>3}'.format(int(c.step[0]))]
             c.time = [c.time[0]]
             for i in range(1, c.maxstep + 1):
@@ -154,8 +154,9 @@ class EcFlexpart(object):
             # (but without 00/12 fields since these are
             # the initialisation times of the flux fields
             # and therefore are zero all the time)
-            self.types[c.type[1]] = {'times': '00/12', 'steps':
-                                     '{}/to/12/by/{}'.format(c.dtime, c.dtime)}
+            self.types[str(c.acctype)] = {'times': str(c.acctime),
+                                          'steps': '{}/to/{}/by/{}'.format(
+                                              c.dtime, c.accmaxstep, c.dtime)}
         else:
             for ty, st, ti in zip(c.type, c.step, c.time):
                 btlist = range(24)
@@ -164,7 +165,12 @@ class EcFlexpart(object):
                 if c.basetime == '00':
                     btlist = [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0]
 
-                if i % int(c.dtime) == 0 and (i in btlist or c.maxstep > 24):
+                if ((ty.upper() == 'AN' and
+                     int(c.time[i]) % int(c.dtime) ==0) or
+                    (ty.upper() != 'AN' and
+                     int(c.step[i]) % int(c.dtime) == 0 and
+                     int(c.step[i]) % int(c.dtime) == 0) ) and \
+                    (int(c.time[i]) in btlist or c.maxstep > 24):
 
                     if ty not in self.types.keys():
                         self.types[ty] = {'times': '', 'steps': ''}
@@ -190,11 +196,7 @@ class EcFlexpart(object):
         self.levelist = c.levelist
         # for gaussian grid retrieval
         self.glevelist = '1/to/' + c.level
-
-        if hasattr(c, 'gaussian') and c.gaussian:
-            self.gaussian = c.gaussian
-        else:
-            self.gaussian = ''
+        self.gaussian = c.gaussian
 
         if 'N' in c.grid:  # Gaussian output grid
             self.grid = c.grid
@@ -238,8 +240,12 @@ class EcFlexpart(object):
                     c.addpar = c.addpar[1:]
                 self.params['OG__SL'][0] += '/' + '/'.join(c.addpar)
 
-            self.params['OG_OROLSM__SL'] = ["160/27/28/173", \
-                                            'SFC', '1', self.grid]
+            if c.marsclass.upper() == 'EA' or c.marsclass.upper() == 'EP':
+                self.params['OG_OROLSM__SL'] = ["160/27/28/244",
+                                                'SFC', '1', self.grid]
+            else:
+                self.params['OG_OROLSM__SL'] = ["160/27/28/173", \
+                                                'SFC', '1', self.grid]
 
             self.params['OG__ML'] = ['T/Q', 'ML', self.levelist, self.grid]
 
@@ -542,7 +548,9 @@ class EcFlexpart(object):
                     pass
                 if pk == 'OG_OROLSM__SL' and not oro:
                     oro = True
-                    retr_param_dict['stream'] = 'OPER'
+                    # in CERA20C (class EP) there is no stream "OPER"!
+                    if self.marsclass.upper() != 'EP':
+                        retr_param_dict['stream'] = 'OPER'
                     retr_param_dict['type'] = 'AN'
                     retr_param_dict['time'] = '00'
                     retr_param_dict['step'] = '000'
@@ -833,7 +841,8 @@ class EcFlexpart(object):
 
                     values = (np.reshape(values, (nj, ni))).flatten() / fak
                     vdp.append(values[:])  # save the accumulated values
-                    if step <= int(c.dtime):
+                    if c.marsclass.upper() == 'EA' or \
+                       step <= int(c.dtime):
                         svdp.append(values[:] / int(c.dtime))
                     else:  # deaccumulate values
                         svdp.append((vdp[-1] - vdp[-2]) / int(c.dtime))
diff --git a/source/python/classes/MarsRetrieval.py b/source/python/classes/MarsRetrieval.py
index 640ad8c..5e5d55c 100644
--- a/source/python/classes/MarsRetrieval.py
+++ b/source/python/classes/MarsRetrieval.py
@@ -462,7 +462,7 @@ class MarsRetrieval(object):
 
         # find all keys without a value and convert all other values to strings
         empty_keys = []
-        for key, value in attrs.itteritems():
+        for key, value in attrs.iteritems():
             if value == '':
                 empty_keys.append(str(key))
             else:
@@ -485,15 +485,16 @@ class MarsRetrieval(object):
                 e = sys.exc_info()[0]
                 print("ERROR: ", e)
                 print('MARS Request failed!')
-            if not self.public and os.stat(target).st_size == 0:
-                print('MARS Request returned no data - please check request')
-                raise IOError
-            elif self.public and os.stat(target).st_size == 0:
-                print('Public MARS Request returned no data - '
-                      'please check request')
-                raise IOError
-            else:
-                raise IOError
+                if not self.public and os.stat(target).st_size == 0:
+                    print('MARS Request returned no data - '
+                          'please check request')
+                    raise IOError
+                elif self.public and os.stat(target).st_size == 0:
+                    print('Public MARS Request returned no data - '
+                          'please check request')
+                    raise IOError
+                else:
+                    raise IOError
         # MARS request via extra process in shell
         else:
             request_str = 'ret'
@@ -514,7 +515,5 @@ class MarsRetrieval(object):
             elif os.stat(target).st_size == 0:
                 print('MARS Request returned no data - please check request')
                 raise IOError
-            else:
-                raise
 
         return
-- 
GitLab