Merge pull request #238 from GEOS-ESM/feature/wjiang/add_regrid_py

adding python version of regrid.pl

CHANGELOG.md

@@ -9,6 +9,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- series of regrid*.py are added to GEOS_Util/post
+
 ### Changed
 
 - Changed restart file geosachem to achem in regrid.pl

GEOS_Util/post/CMakeLists.txt

@@ -50,6 +50,11 @@ install(
    PROGRAMS ${perlscripts} ec2grd.csh
    DESTINATION bin)
 
+file(GLOB pythonscripts CONFIGURE_DEPENDS ./remap_restart/remap*)
+install(
+   PROGRAMS ${pythonscripts} 
+   DESTINATION bin)
+
 set_property(SOURCE rsg3_vinterp.F APPEND_STRING PROPERTY COMPILE_FLAGS "${FREAL8} ${BYTERECLEN} ${EXTENDED_SOURCE}")
 set_property(SOURCE rs_vinterp.F90 APPEND_STRING PROPERTY COMPILE_FLAGS "${FREAL8} ${BYTERECLEN}")
 set_property(SOURCE rs_vinterp_scm.F90 APPEND_STRING PROPERTY COMPILE_FLAGS "${FREAL8} ${BYTERECLEN}")

GEOS_Util/post/remap_restart/remap_analysis.py

@@ -0,0 +1,137 @@
+#!/usr/bin/env python3
+#
+# source install/bin/g5_modules
+#
+# Newer GEOS code should load a module with GEOSpyD Python3 if not run:
+#   module load python/GEOSpyD/Min4.10.3_py3.9
+#
+
+import os
+from datetime import datetime, timedelta
+import subprocess
+import shlex
+import shutil
+import glob
+import fileinput
+import ruamel.yaml
+from remap_base import remap_base
+
+class analysis(remap_base):
+  def __init__(self, **configs):
+     super().__init__(**configs)
+
+  def remap(self):
+     config = self.config
+     bkg = config['output']['analysis']['bkg']
+     if ( not bkg ): return
+
+     analysis_in = self.find_analysis()
+     if len(analysis_in) ==0 :
+       print("\n There are no analysis files. \n")
+       return
+
+     print("\n Remapping or copying analysis files...\n")
+
+     cwdir  = os.getcwd()
+     bindir = os.path.dirname(os.path.realpath(__file__))
+     in_bcsdir  = config['input']['shared']['bcs_dir']
+     out_bcsdir = config['output']['shared']['bcs_dir']
+     out_dir    = config['output']['shared']['out_dir']
+     if not os.path.exists(out_dir) : os.makedirs(out_dir)
+     print( "cd " + out_dir)
+     os.chdir(out_dir)
+
+     tmpdir = out_dir+'/ana_data/'
+     if os.path.exists(tmpdir) : subprocess.call(['rm', '-rf',tmpdir])
+     print ("mkdir " + tmpdir)
+     os.makedirs(tmpdir)
+
+     print( "cd " + tmpdir)
+     os.chdir(tmpdir)
+
+     yyyymmddhh_ = str(config['input']['shared']['yyyymmddhh'])
+     yyyy_ = yyyymmddhh_[0:4]
+     mm_   = yyyymmddhh_[4:6]
+     dd_   = yyyymmddhh_[6:8]
+     hh_   = yyyymmddhh_[8:10]
+     rst_time = datetime(year=int(yyyy_), month=int(mm_), day=int(dd_), hour = int(hh_))
+     expid_in  = config['input']['shared']['expid']
+     expid_out = config['output']['shared']['expid']
+     if (expid_out) :
+        expid_out = expid_out + '.'
+     else:
+        expid_out = ''
+
+     aqua   = config['output']['analysis']['aqua']
+     local_fs=[]
+     for f in analysis_in:
+       print(f)
+       fname    = os.path.basename(f)
+       out_name = fname.replace(expid_in + '.', expid_out)
+       f_tmp = tmpdir+'/'+out_name
+       local_fs.append(f_tmp) 
+       shutil.copy(f,f_tmp)
+       if out_name.find('satbias') != -1 :
+         if (aqua):
+            f_ = open(f_tmp, 'w')
+            for line in fileinput.input(f):
+               f_.write(line.replace('airs281SUBSET_aqua', 'airs281_aqua      '))
+            f_.close() 
+
+     nlevel    = config['output']['air']['nlevel']
+     agrid_out = config['output']['shared']['agrid']
+     flags = "-g5 -res " + self.get_grid_kind(agrid_out.upper()) + " -nlevs " + str(nlevel)
+     bkg_files = glob.glob(tmpdir+'/*.bkg??_eta_rst*')
+     for f in bkg_files:
+        f_orig = f + ".orig"
+        shutil.move(f,f_orig)
+        cmd = bindir + '/dyn2dyn.x ' + flags + ' -o ' + f + ' ' + f_orig
+        print(cmd)
+        subprocess.call(shlex.split(cmd))
+
+     for f in local_fs:
+       fname = os.path.basename(f)
+       shutil.move(f, out_dir+'/'+fname)
+    # write lcv 
+     lcv = config['output']['analysis']['lcv']
+     if lcv :
+       ymd_ = yyyymmddhh_[0:8]
+       hh_  = yyyymmddhh_[8:10]
+       hms_ = hh_+'0000'
+       rstlcvOut = out_dir+'/'+expid_out+'rst.lcv.'+ymd_+'_'+hh_+'z.bin'
+       cmd = bindir+'/mkdrstdate.x ' + ymd_ + ' ' + hms_ +' ' + rstlcvOut
+       print(cmd)
+       subprocess.call(shlex.split(cmd))
+     print( "cd " + cwdir)
+     os.chdir(cwdir)
+
+  def get_grid_kind(this, grid):
+     hgrd = {}
+     hgrd['C12']   = 'a'
+     hgrd['C24']   = 'a'
+     hgrd['C48']   = 'b'
+     hgrd['C90']   = 'c'
+     hgrd['C180']  = 'd'
+     hgrd['C360']  = 'd'
+     hgrd['C500']  = 'd'
+     hgrd['C720']  = 'e'
+     hgrd['C1000'] = 'e'
+     hgrd['C1440'] = 'e'
+     hgrd['C2000'] = 'e'
+     hgrd['C2880'] = 'e'
+     hgrd['C5760'] = 'e'
+     return hgrd[grid]
+
+  def find_analysis(self):
+     analysis_in = []
+     rst_dir = self.config['input']['shared']['rst_dir']
+     bkgs = glob.glob(rst_dir + '/*_eta_rst*')
+     sfcs = glob.glob(rst_dir + '/*_sfc_rst*')
+     anasat = glob.glob(rst_dir + '/*ana_satb*')
+     traks= glob.glob(rst_dir + '/*.trak.GDA.rst*')
+     analysis_in = bkgs + sfcs + traks + anasat
+     return list(dict.fromkeys(analysis_in))
+
+if __name__ == '__main__' :
+   ana = analysis(params_file='remap_params.yaml')
+   ana.remap()

GEOS_Util/post/remap_restart/remap_base.py

@@ -0,0 +1,30 @@
+#!/usr/bin/env python3
+#
+import os
+import ruamel.yaml
+import shutil
+
+class remap_base(object):
+  def __init__(self, **configs):
+     for key, value in configs.items():
+        if (key == 'params_file'):
+          print( "use Config yaml file: " + value)
+          yaml = ruamel.yaml.YAML()
+          stream =''
+          with  open(value, 'r') as f:
+            stream = f.read()
+          self.config = yaml.load(stream)
+          out_dir    = self.config['output']['shared']['out_dir']
+          if not os.path.exists(out_dir) : os.makedirs(out_dir)
+          f = os.path.basename(value)
+          dest = out_dir+'/'+f
+          try:
+            shutil.copy(value, dest)
+          except shutil.SameFileError:
+            pass
+        if (key == 'config_obj'):
+          print( "use Config obj")
+          self.config = value
+          out_dir    = self.config['output']['shared']['out_dir']
+          if not os.path.exists(out_dir) : os.makedirs(out_dir)
+        break

GEOS_Util/post/remap_restart/remap_catchANDcn.py

@@ -0,0 +1,134 @@
+#!/usr/bin/env python3
+#
+import os
+import sys
+import subprocess
+import shutil
+import glob
+import ruamel.yaml
+import shlex
+from remap_base import remap_base
+
+class catchANDcn(remap_base):
+  def __init__(self, **configs):
+     super().__init__(**configs)
+
+  def remap(self):
+     if not self.config['output']['surface']['remap_catch']:
+        return
+
+     config = self.config
+     rst_dir = config['input']['shared']['rst_dir']
+     model = config['input']['surface']['catch_model']
+     in_rstfile = ''
+     yyyymmddhh_ = str(self.config['input']['shared']['yyyymmddhh'])
+     time = yyyymmddhh_[0:8]+'_'+yyyymmddhh_[8:10]
+     in_rstfiles = glob.glob(rst_dir+'/*'+model+'_*'+time+'*')
+     if len(in_rstfiles) == 0:
+        print('\n try catchXX file without time stamp')
+        in_rstfiles = glob.glob(rst_dir+'/*'+model+'_*')
+     if len(in_rstfiles) == 0:
+        return
+     in_rstfile = in_rstfiles[0]  
+
+     print("\nRemapping " + model + ".....\n")
+
+     cwdir  = os.getcwd()
+     bindir = os.path.dirname(os.path.realpath(__file__))
+     in_bcsdir  = config['input']['shared']['bcs_dir']
+     out_bcsdir = config['output']['shared']['bcs_dir']
+     out_dir    = config['output']['shared']['out_dir']
+     expid      = config['output']['shared']['expid']
+     in_wemin   = config['input']['surface']['wemin']
+     out_wemin  = config['output']['surface']['wemin']
+     surflay    = config['output']['surface']['surflay']
+     in_tilefile  = glob.glob(in_bcsdir+ '/*-Pfafstetter.til')[0]
+     out_tilefile = glob.glob(out_bcsdir+ '/*-Pfafstetter.til')[0]
+     account    = config['slurm']['account']
+     # even the input is binary, the output si nc4
+     suffix     = time+'z.nc4'
+
+     if (expid) :
+        expid = expid + '.'
+     else:
+        expid = ''
+     suffix = '_rst.' + suffix
+     out_rstfile = expid + os.path.basename(in_rstfile).split('_rst')[0].split('.')[-1]+suffix
+
+     if not os.path.exists(out_dir) : os.makedirs(out_dir)
+     print( "cd " + out_dir)
+     os.chdir(out_dir)
+
+     InData_dir = out_dir+'/InData/'
+     print ("mkdir -p " + InData_dir)
+     os.makedirs(InData_dir, exist_ok = True)
+
+     f = os.path.basename(in_rstfile)
+     dest = InData_dir+'/'+f
+     # file got copy because the computing node cannot access archive
+     print('\nCopy ' + in_rstfile + ' to ' +dest)
+     shutil.copyfile(in_rstfile,dest)
+     in_rstfile = dest
+
+     log_name = out_dir+'/'+'mk_catchANDcn_log'
+     mk_catch_j_template = """#!/bin/csh -f
+#SBATCH --account={account}
+#SBATCH --ntasks=56
+#SBATCH --time=1:00:00
+#SBATCH --job-name=mk_catchANDcn
+#SBATCH --qos=debug
+#SBATCH --output={log_name}
+#
+
+source {Bin}/g5_modules
+
+limit stacksize unlimited
+
+set esma_mpirun_X = ( {Bin}/esma_mpirun -np 56 )
+set mk_catchANDcnRestarts_X   = ( {Bin}/mk_catchANDcnRestarts.x )
+
+set params = ( -model {model}  -time {time} -in_tilefile {in_tilefile} )
+set params = ( $params -out_bcs {out_bcs} -out_tilefile {out_tilefile} -out_dir {out_dir} )
+set params = ( $params -surflay {surflay} -in_wemin {in_wemin} -out_wemin {out_wemin} ) 
+set params = ( $params -in_rst {in_rstfile} -out_rst {out_rstfile} ) 
+$esma_mpirun_X $mk_catchANDcnRestarts_X $params
+
+"""
+     catch1script =  mk_catch_j_template.format(Bin = bindir, account = account, out_bcs = out_bcsdir, \
+                  model = model, out_dir = out_dir, surflay = surflay, log_name = log_name,  \
+                  in_wemin   = in_wemin, out_wemin = out_wemin, out_tilefile = out_tilefile, in_tilefile = in_tilefile, \
+                  in_rstfile = in_rstfile, out_rstfile = out_rstfile, time = yyyymmddhh_ )
+
+     script_name = './mk_catchANDcn.j'
+
+     catch_scrpt = open(script_name,'wt')
+     catch_scrpt.write(catch1script)
+     catch_scrpt.close()
+
+     interactive = os.getenv('SLURM_JOB_ID', default = None)
+     if ( interactive ) :
+       print('interactive mode\n')
+       ntasks = os.getenv('SLURM_NTASKS', default = None)
+       if ( not ntasks):
+         nnodes = int(os.getenv('SLURM_NNODES', default = '1'))
+         ncpus  = int(os.getenv('SLURM_CPUS_ON_NODE', default = '28'))
+         ntasks = nnodes * ncpus
+       ntasks = int(ntasks)
+       NPE = 56
+       if (ntasks < NPE):
+         print("\nYou should have at least {NPE} cores. Now you only have {ntasks} cores ".format(NPE=NPE, ntasks=ntasks))
+
+       subprocess.call(['chmod', '755', script_name])
+       print(script_name+  '  1>' + log_name  + '  2>&1')
+       os.system(script_name + ' 1>' + log_name+ ' 2>&1')
+
+     else:
+       print("sbatch -W " + script_name +"\n")
+       subprocess.call(['sbatch','-W', script_name])
+
+     print( "cd " + cwdir)
+     os.chdir(cwdir)
+
+if __name__ == '__main__' :
+   catch = catchANDcn(params_file='remap_params.yaml')
+   catch.remap()