ldas_setup

#!/usr/bin/env python3

import os
import sys
import glob
import copy
import linecache
import shutil
import argparse
import fileinput
import time
import resource
import subprocess as sp
import shlex
import tempfile
import netCDF4

from dateutil                     import rrule
from datetime                     import datetime
from datetime                     import timedelta
from collections                  import OrderedDict
from dateutil.relativedelta       import relativedelta
from remap_utils                  import *
from remap_lake_landice_saltwater import *
from remap_catchANDcn             import *
from lenkf_j_template             import *

"""
This script is intended to be run from any installed directory with GEOSldas.x and ldas_setup
(The default setup is ../install/bin)
"""

class LDASsetup:

    def __init__(self, cmdLineArgs):
        """
        """
        # ------
        # Required exe input fields
        # These fields are needed to pre-compute exp dir structure
        # ------
        rqdExeInpKeys = ['EXP_ID', 'EXP_DOMAIN', 'NUM_LDAS_ENSEMBLE',
            'BEG_DATE', 'END_DATE','RESTART_PATH',
            'RESTART_DOMAIN','RESTART_ID','MET_TAG','MET_PATH','FORCE_DTSTEP','BCS_PATH', 'BCS_RESOLUTION']
        rqdExeInpKeys_rst = ['EXP_ID', 'EXP_DOMAIN', 'NUM_LDAS_ENSEMBLE',
            'BEG_DATE', 'END_DATE','MET_TAG','MET_PATH','FORCE_DTSTEP','BCS_PATH', 'BCS_RESOLUTION']

        # These keywords are excluded from LDAS.rc (i.e., only needed in pre- or post-processing)
        self.NoneLDASrcKeys=['EXP_ID', 'EXP_DOMAIN',
            'BEG_DATE', 'END_DATE','RESTART','RESTART_PATH',
            'RESTART_DOMAIN','RESTART_ID','BCS_PATH','TILING_FILE','GRN_FILE','LAI_FILE','LNFM_FILE','NIRDF_FILE',
            'VISDF_FILE','CATCH_DEF_FILE','NDVI_FILE',
            'NML_INPUT_PATH','HISTRC_FILE','RST_FROM_GLOBAL','JOB_SGMT','NUM_SGMT','POSTPROC_HIST',
            'MINLON','MAXLON','MINLAT','MAXLAT','EXCLUDE_FILE','INCLUDE_FILE','MWRTM_PATH','GRIDNAME',
            'ADAS_EXPDIR', 'BCS_RESOLUTION', 'TILE_FILE_FORMAT' ]

        self.GEOS_SITE = "@GEOS_SITE@"

        # ------
        # Required resource manager input fields
        # ------
        rqdRmInpKeys    = ['account', 'walltime', 'ntasks_model']
        # ------
        # Optional resource manager input fields
        # ------
        optSlurmInpKeys = ['job_name', 'qos', 'oserver_nodes', 'writers-per-node', 'constraint']

        
        # ===============================================================================================
        #    
        # ------
        # ./ldas_setup sample ...
        # ------
        #
        # "sample" sub-command:
        #   '--exeinp' and '--batinp' are mutually exclusive command line arguments.
        #   Specifying one will set it to True and set the other one to False.
        #   That is, we can have either: {'exeinp': False, 'batinp': True }
        #                            or: {'exeinp': True,  'batinp': False} 
        
        if 'exeinp' in cmdLineArgs:     # 'exeinp' is always present in "sample" mode.
            
            if   cmdLineArgs['exeinp']:
                _produceExeInput()
            elif cmdLineArgs['batinp']:
                _printRmInputKeys( rqdRmInpKeys,  optSlurmInpKeys)
            else:
                raise Exception('unrecognized option')
            #
            # EXIT after completing "sample" sub-command
            sys.exit(0)
            

        # ===============================================================================================
        #    
        # ------
        # ./ldas_setup setup ...
        # ------
        # Instance variables
        self.exeinpfile         = cmdLineArgs['exeinpfile']
        self.batinpfile         = cmdLineArgs['batinpfile']
        exphome_                = cmdLineArgs['exphome'].rstrip('/')
        assert os.path.isdir(exphome_)                                   # exphome should exist
        self.exphome            = os.path.abspath(exphome_)
        self.verbose            = cmdLineArgs['verbose']

        # command line args for coupled land-atm DAS (see "help" strings in parseCmdLine() for details)
        self.ladas_cpl          = cmdLineArgs['ladas_cpl']
        self.nymdb              = cmdLineArgs['nymdb']
        self.nhmsb              = cmdLineArgs['nhmsb']
        self.agcm_res           = cmdLineArgs['agcm_res']
        self.bcs_version        = cmdLineArgs['bcs_version']
        self.rstloc             = cmdLineArgs['rstloc']
        self.varwindow          = cmdLineArgs['varwindow'] 
        self.nens               = cmdLineArgs['nens']

        # obsolete command line args
        self.runmodel           = cmdLineArgs['runmodel']
        if self.runmodel :
          print('\n The option "--runmodel" is out of date, not necessary anymore. \n')

        self.daysperjob         = cmdLineArgs['daysperjob']
        self.monthsperjob       = cmdLineArgs['monthsperjob']

        self.rqdExeInp          = OrderedDict()
        self.rqdRmInp           = OrderedDict()
        self.optRmInp           = OrderedDict()
        self.rundir             = None
        self.blddir             = None
        self.blddirLn           = None
        self.outdir             = None
        self.out_path           = None
        self.inpdir             = None
        self.exefyl             = None
        self.isZoomIn           = False
        self.catch              = ''
        self.has_mwrtm          = False
        self.has_vegopacity     = False
        self.assim              = False
        self.has_landassim_seed = False
        self.has_geos_pert      = False
        self.nSegments          = 1
        self.perturb            = 0
        self.first_ens_id       = 0
        self.in_rstfile         = None
        self.in_tilefile        = None                 # default string
        self.ens_id_width       = 6                    # _eXXXX
        self.bcs_dir_land       = ''
        self.bcs_dir_geom       = ''
        self.bcs_dir_landshared = ''
        self.tile_types         = ''
        self.with_land          = False
        self.with_landice       = False
        self.adas_expdir        = ''

        # ------
        # Read exe input file which is required to set up the dir
        # ------
        if self.ladas_cpl is None:
           self.ladas_cpl = 0
        else:
           self.ladas_cpl = int(self.ladas_cpl)

        self.rqdExeInp = {}
        if self.ladas_cpl == 0:
           self.rqdExeInp = self._parseInputFile(cmdLineArgs['exeinpfile'])
        else:
           _produceExeInput(out_dict=self.rqdExeInp, ladas_cpl=self.ladas_cpl)

        # verifing the required input
        if 'RESTART' not in self.rqdExeInp :
            self.rqdExeInp['RESTART'] = "1"

        if self.rqdExeInp['RESTART'].isdigit() :
            if int(self.rqdExeInp['RESTART']) ==0 :
                rqdExeInpKeys = rqdExeInpKeys_rst
                self.rqdExeInp['RESTART_ID']     = 'None'
                self.rqdExeInp['RESTART_DOMAIN'] = 'None'
                self.rqdExeInp['RESTART_PATH']   = 'None'
        else:
            if self.rqdExeInp['RESTART'] =='G' :
                rqdExeInpKeys = rqdExeInpKeys_rst
                self.rqdExeInp['RESTART_DOMAIN'] = 'None'
            else:
                self.rqdExeInp['RESTART_ID']     = 'None'
                self.rqdExeInp['RESTART_DOMAIN'] = 'None'
                self.rqdExeInp['RESTART_PATH']   = 'None'
    
        ### check if ldas is coupled to adas; if so, set/overwrite input parameters accordingly
        if self.ladas_cpl > 0 :

            # make sure all necessary command line arguments were supplied
            assert self.nymdb         is not None, "Error. Must have command line arg nymdb       for coupled land-atm DAS.\n"
            assert self.nhmsb         is not None, "Error. Must have command line arg nhmsb       for coupled land-atm DAS.\n"
            assert self.agcm_res      is not None, "Error. Must have command line arg agcm_res    for coupled land-atm DAS.\n"
            assert self.bcs_version   is not None, "Error. Must have command line arg bcs_version for coupled land-atm DAS.\n"
            assert self.rstloc        is not None, "Error. Must have command line arg rstloc      for coupled land-atm DAS.\n"
            assert self.varwindow     is not None, "Error. Must have command line arg varwindow   for coupled land-atm DAS.\n"
            assert self.nens          is not None, "Error. Must have command line arg nens        for coupled land-atmensDAS.\n"
           
            self.rqdExeInp['BEG_DATE']                 = f"{self.nymdb} {self.nhmsb}"
            rstloc_                                    = self.rstloc.rstrip('/')                  # remove trailing '/'
            assert os.path.isdir(rstloc_)                                                         # make sure rstloc_ is a valid directory
            self.rstloc                                = os.path.abspath(rstloc_)
            self.rqdExeInp['RESTART_PATH']             = os.path.dirname( self.rstloc) 
            self.rqdExeInp['RESTART_ID']               = os.path.basename(self.rstloc)
            self.adas_expdir                           = os.path.dirname( self.exphome)
            self.rqdExeInp['ADAS_EXPDIR']              = self.adas_expdir
            self.adas_expid                            = os.path.basename(self.adas_expdir)
            self.rqdExeInp['MET_TAG']                  = self.adas_expid + '__bkg'
            
            if   self.ladas_cpl == 1 :
                # ldas coupled with determistic component of ADAS
                self.rqdExeInp['EXP_ID']               = self.adas_expid  + '_LDAS'
                self.rqdExeInp['MET_PATH']             = self.adas_expdir + '/recycle/holdpredout'
                self.rqdExeInp['ENSEMBLE_FORCING']     = 'NO' 
            elif self.ladas_cpl == 2 :
                # ldas coupled with ensemble component of ADAS
                self.rqdExeInp['EXP_ID']               = self.adas_expid  + '_LDAS4ens'
                self.rqdExeInp['MET_PATH']             = self.adas_expdir + '/atmens/mem'
                self.rqdExeInp['ENSEMBLE_FORCING']     = 'YES' 
            else :
                exit("Error. Unknown value of self.ladas_cpl.\n")
           
            self.rqdExeInp['NUM_LDAS_ENSEMBLE']        = self.nens                                # fvsetup finds Nens by counting restart files
            self.first_ens_id                          = 1                                        # match ADAS convention
            self.rqdExeInp['FIRST_ENS_ID']             = self.first_ens_id

            self.agcm_res                              = 'CF' + self.agcm_res                     # change format to "CFnnnn"
            self.rqdExeInp['EXP_DOMAIN']               = self.agcm_res +'x6C_GLOBAL'
           
            # when coupled to ADAS, "BCS_PATH" EXCLUDE bcs version info  
            # hard-wired BCS_PATH for now 
            self.rqdExeInp['BCS_PATH']                 = "/discover/nobackup/projects/gmao/bcs_shared/fvInput/ExtData/esm/tiles"
            self.rqdExeInp['BCS_PATH']                 = self.rqdExeInp['BCS_PATH'].rstrip('/') + '/' + self.bcs_version 
            if self.bcs_version == "Icarus-NLv3" :
                self.rqdExeInp['BCS_PATH']             = self.rqdExeInp['BCS_PATH']  + '_new_layout'  
            self.rqdExeInp['BCS_RESOLUTION']           = self.agcm_res +'x6C_'       + self.agcm_res +'x6C'
            self.rqdExeInp['RESTART_DOMAIN']           = self.agcm_res +'x6C_GLOBAL'

            # the following are not in default rqdExeInp list; hardwire for now
            self.rqdExeInp['MWRTM_PATH']               = '/discover/nobackup/projects/gmao/smap/LDAS_inputs_for_LADAS/RTM_params/RTMParam_SMAP_L4SM_v006/'
            self.rqdExeInp['LAND_ASSIM']               = "YES"
            self.rqdExeInp['MET_HINTERP']              = 0
            self.landassim_dt                          = 10800                                    # seconds
            # make sure ADAS analysis window [minutes] is multiple of LANDASSIM_DT [seconds]
            if int(self.varwindow) % (self.landassim_dt/60) == 0 : 
                self.rqdExeInp['LANDASSIM_DT']         = self.landassim_dt
            else :
                exit("Error. LANDASSIM_DT is inconsistent with ADAS analysis window.\n")
            self.rqdExeInp['LANDASSIM_T0']             = "013000"                                 # HHMMSS
            jsgmt1                                     = "00000000"
            jsgmt2                                     = hours_to_hhmmss(int(self.varwindow)/60)  # convert minutes to HHMMSS
            self.rqdExeInp['JOB_SGMT']                 = f"{jsgmt1} {jsgmt2}"
            self.rqdExeInp['NUM_SGMT']                 = 1 
            self.rqdExeInp['FORCE_DTSTEP']             = 3600

            # determine END_DATE = BEG_DATE + TIME_STEP_OF_ADAS_CYCLE
            _beg_date                                  = datetime.strptime( self.rqdExeInp['BEG_DATE'], "%Y%m%d %H%M%S")
            _hours                                     = int(self.rqdExeInp['JOB_SGMT'][ 9:11])
            _end_date                                  = _beg_date + timedelta(hours=int(self.varwindow)/60) 
            self.rqdExeInp['END_DATE']                 = _end_date.strftime("%Y%m%d %H%M%S") 
           
            # end if self.ladas_cpl > 0 -----------------------------------------------------------------------------------------

        for key in rqdExeInpKeys :
           assert key in self.rqdExeInp,' "%s" is required in the input file %s' % (key,self.exeinpfile)

        # print rqd exe inputs 
        if self.verbose:       
            print ('\nInputs from exeinp file:\n')
            _printdict(self.rqdExeInp)

        self.tile_types = self.rqdExeInp.get('TILE_TYPES',"100").split()
        if "100" in self.tile_types :
          self.with_land = True
        if  "20" in self.tile_types :
          self.with_landice = True

        # nens is an integer and =1 for model run
        self.nens = int(self.rqdExeInp['NUM_LDAS_ENSEMBLE']) # fail if Nens's val is not int
        assert self.nens>0, 'NUM_LDAS_ENSEMBLE [%d] <= 0' % self.nens
        _mydir = self.exphome + '/' + self.rqdExeInp['EXP_ID']
        assert not os.path.isdir(_mydir), 'Dir [%s] already exists!' % _mydir
        _mydir = None
        self.first_ens_id = int(self.rqdExeInp.get('FIRST_ENS_ID',0))
        
        self.perturb = int(self.rqdExeInp.get('PERTURBATIONS',0))
        if self.nens > 1:
           self.perturb = 1
        self.ensdirs  = ['ens%04d'%iens for iens in range(self.first_ens_id, self.nens  + self.first_ens_id)]
        # if self.ens_id_width = 4, _width = '_e%04d'
        _width = '_e%0{}d'.format(self.ens_id_width-2)
        # self.ensids will be a list of [_e0000, _e0001, ...]
        self.ensids   = [ _width%iens   for iens in range(self.first_ens_id, self.nens  + self.first_ens_id)]
        if (self.nens == 1) :
            self.ensdirs_avg = self.ensdirs
            self.ensids=['']
        else :
            self.ensdirs_avg = self.ensdirs + ['ens_avg']

        ## convert date-time strings to datetime object
        ## start/end_time are converted to lists
        ## ensure end>start

        self.begDates=[]
        self.endDates=[]
        self.begDates.append(
                datetime.strptime(
                    self.rqdExeInp['BEG_DATE'],
                    '%Y%m%d %H%M%S'
                    )
                )
        self.endDates.append(
                datetime.strptime(
                    self.rqdExeInp['END_DATE'],
                    '%Y%m%d %H%M%S'
                    )
                )
        if self.rqdExeInp['RESTART'].isdigit() :
            if int(self.rqdExeInp['RESTART']) == 0 :
                print ("No restart file (cold restart): Forcing start date to January 1, 0z")
                year = self.begDates[0].year
                self.begDates[0]=datetime(year =year,month=1,day =1,hour =0, minute= 0,second= 0)

        assert self.endDates[0]>self.begDates[0], \
            'END_DATE <= BEG_DATE'

        self.job_sgmt = []
        if 'JOB_SGMT' in self.rqdExeInp:
            self.job_sgmt.append("JOB_SGMT: "+self.rqdExeInp['JOB_SGMT'])
        else:
            _datediff = relativedelta(self.endDates[0],self.begDates[0])
            self.rqdExeInp['JOB_SGMT'] = "%04d%02d%02d %02d%02d%02d" %(_datediff.years,
                                               _datediff.months,
                                               _datediff.days,
                                               _datediff.hours,
                                               _datediff.minutes,
                                               _datediff.seconds)
            self.job_sgmt.append("JOB_SGMT: "+self.rqdExeInp['JOB_SGMT'])

        if 'NUM_SGMT' not in self.rqdExeInp:
           self.rqdExeInp['NUM_SGMT'] = 1

        _years  = int(self.rqdExeInp['JOB_SGMT'][ 0: 4])
        _months = int(self.rqdExeInp['JOB_SGMT'][ 4: 6])
        _days   = int(self.rqdExeInp['JOB_SGMT'][ 6: 8])
        assert self.rqdExeInp['JOB_SGMT'][8] == ' ' and self.rqdExeInp['JOB_SGMT'][9] != ' ', "JOB_SGMT format is not right"
        _hours  = int(self.rqdExeInp['JOB_SGMT'][ 9:11])
        _mins   = int(self.rqdExeInp['JOB_SGMT'][11:13])
        _seconds= int(self.rqdExeInp['JOB_SGMT'][13:15])


        _difftime =timedelta(days = _years*365+_months*30+_days,hours = _hours,minutes=_mins,seconds=_seconds)
        _difftime = int(self.rqdExeInp['NUM_SGMT'])*_difftime
        _d = self.begDates[0]
        _endDate = self.endDates[0]
        _d = _d + _difftime
        while _d < _endDate :
            print (_difftime.days)
            self.nSegments +=1
            print (_d.year, _d.month, _d.day)
            self.begDates.append(_d)
            self.endDates.insert(-1,_d)
            _d = _d+ _difftime

        # assemble bcs sub-directories 
        self.bcs_dir_land       = self.rqdExeInp['BCS_PATH']+ '/land/'       + self.rqdExeInp['BCS_RESOLUTION']+'/'
        self.bcs_dir_geom       = self.rqdExeInp['BCS_PATH']+ '/geometry/'   + self.rqdExeInp['BCS_RESOLUTION']+'/'
        self.bcs_dir_landshared = self.rqdExeInp['BCS_PATH']+ '/land/shared/'

        # make sure MET_PATH and RESTART_PATH have trailing '/'
        if self.rqdExeInp['MET_PATH'][-1] != '/':
           self.rqdExeInp['MET_PATH'] = self.rqdExeInp['MET_PATH']+'/'
        if self.rqdExeInp['RESTART_PATH'][-1] != '/':
           self.rqdExeInp['RESTART_PATH'] = self.rqdExeInp['RESTART_PATH']+'/'

        # make sure catchment and vegdyn restart files ( at least one for each) exist
        if 'CATCH_DEF_FILE' not in self.rqdExeInp :
           self.rqdExeInp['CATCH_DEF_FILE']= self.bcs_dir_land + 'clsm/catchment.def'
        if (self.with_land) :
           assert os.path.isfile(self.rqdExeInp['CATCH_DEF_FILE']),"[%s] file does not exist " % self.rqdExeInp['CATCH_DEF_FILE']

        self.rqdExeInp['RST_FROM_GLOBAL'] = 1
        # skip checking. It is users' reponsibility to make it right!
        #if self.rqdExeInp['RESTART'].isdigit() :
        #    if int(self.rqdExeInp['RESTART']) == 1 :
        #        _numg = int(linecache.getline(self.rqdExeInp['CATCH_DEF_FILE'], 1).strip())
        #        _numd = _numg
        #        ldas_domain = self.rqdExeInp['RESTART_PATH']+ \
        #                      self.rqdExeInp['RESTART_ID'] + \
        #                      '/output/'+self.rqdExeInp['RESTART_DOMAIN']+'/rc_out/'+self.rqdExeInp['RESTART_ID']+'.ldas_domain.txt'
        #        if os.path.isfile(ldas_domain) :
        #            _numd = int(linecache.getline(ldas_domain, 1).strip())
        #
        #            if _numg != _numd :
        #                self.rqdExeInp['RST_FROM_GLOBAL'] = 0

        self.rqdExeInp['LNFM_FILE'] = ''
        tile_file_format = self.rqdExeInp.get('TILE_FILE_FORMAT', 'DEFAULT')
        if int(self.rqdExeInp['RST_FROM_GLOBAL']) == 1 :
           txt_tile = glob.glob(self.bcs_dir_geom + '*.til')
           nc4_tile = glob.glob(self.bcs_dir_geom + '*.nc4')
           if tile_file_format.upper() == 'TXT'     : self.rqdExeInp['TILING_FILE'] =  txt_tile[0]
           if tile_file_format.upper() == 'DEFAULT' : self.rqdExeInp['TILING_FILE'] = (txt_tile+nc4_tile)[-1]

           self.rqdExeInp['GRN_FILE']    = glob.glob(self.bcs_dir_land + 'green_clim_*.data')[0]
           self.rqdExeInp['LAI_FILE']    = glob.glob(self.bcs_dir_land + 'lai_clim_*.data'  )[0]
           tmp_ = glob.glob(self.bcs_dir_land + 'lnfm_clim_*.data')
           if (len(tmp_) ==1) :
              self.rqdExeInp['LNFM_FILE'] = tmp_[0]
           self.rqdExeInp['NDVI_FILE']   = glob.glob(self.bcs_dir_land + 'ndvi_clim_*.data' )[0]
           self.rqdExeInp['NIRDF_FILE']  = glob.glob(self.bcs_dir_land + 'nirdf_*.dat'      )[0]
           self.rqdExeInp['VISDF_FILE']  = glob.glob(self.bcs_dir_land + 'visdf_*.dat'      )[0]
        else :
           inpdir=self.rqdExeInp['RESTART_PATH']+self.rqdExeInp['RESTART_ID']+'/input/'
           self.rqdExeInp['TILING_FILE'] = os.path.realpath(glob.glob(inpdir+'*tile.data')[0])
           self.rqdExeInp['GRN_FILE']    = os.path.realpath(glob.glob(inpdir+'green*data')[0])
           self.rqdExeInp['LAI_FILE']    = os.path.realpath(glob.glob(inpdir+'lai*data'  )[0])
           tmp_ = glob.glob(self.bcs_dir_land + 'lnfm_clim_*.data')
           if (len(tmp_) == 1) :
              self.rqdExeInp['LNFM_FILE'] = tmp_[0]
           self.rqdExeInp['NDVI_FILE']   = os.path.realpath(glob.glob(inpdir+'ndvi*data' )[0])
           self.rqdExeInp['NIRDF_FILE']  = os.path.realpath(glob.glob(inpdir+'nirdf*data')[0])
           self.rqdExeInp['VISDF_FILE']  = os.path.realpath(glob.glob(inpdir+'visdf*data')[0])

        if self.rqdExeInp['RESTART'].isdigit() :
            if int(self.rqdExeInp['RESTART']) == 2 :
                self.rqdExeInp['RST_FROM_GLOBAL'] = 1
                ldas_domain = self.rqdExeInp['RESTART_PATH']+ \
                              self.rqdExeInp['RESTART_ID'] + \
                              '/output/'+self.rqdExeInp['RESTART_DOMAIN']+'/rc_out/'+self.rqdExeInp['RESTART_ID']+'.ldas_domain.txt'
                inpdir=self.rqdExeInp['RESTART_PATH']+self.rqdExeInp['RESTART_ID']+'/input/'
                in_tilefiles_ = glob.glob(inpdir+'*tile.data')
                if len(in_tilefiles_) == 0 :
                   inpdir=self.rqdExeInp['RESTART_PATH']+self.rqdExeInp['RESTART_ID']+'/output/'+self.rqdExeInp['RESTART_DOMAIN']+'/rc_out/'
                   in_tilefiles_ = glob.glob(inpdir+'MAPL_*.til')
                   if len(in_tilefiles_) == 0 :
                      in_tilefiles_ = glob.glob(inpdir+'/*.til')
                   if len(in_tilefiles_) == 0 :
                      in_tilefiles_ = glob.glob(inpdir+'/*.nc4')

                self.in_tilefile =os.path.realpath(in_tilefiles_[0])

                if os.path.isfile(ldas_domain):
                    txt_tile = glob.glob(self.bcs_dir_geom + '*.til')
                    nc4_tile = glob.glob(self.bcs_dir_geom + '*.nc4')
                    if tile_file_format.upper() == 'TXT'     : self.rqdExeInp['TILING_FILE'] = txt_tile[0]
                    if tile_file_format.upper() == 'DEFAULT' : self.rqdExeInp['TILING_FILE'] = (txt_tile+nc4_tile)[-1]

                    self.rqdExeInp['GRN_FILE']   = glob.glob(self.bcs_dir_land + 'green_clim_*.data')[0]
                    self.rqdExeInp['LAI_FILE']   = glob.glob(self.bcs_dir_land + 'lai_clim_*.data'  )[0]
                    tmp_ = glob.glob(self.bcs_dir_land + 'lnfm_clim_*.data')
                    if (len(tmp_) == 1) :
                       self.rqdExeInp['LNFM_FILE'] = tmp_[0]
                    self.rqdExeInp['LNFM_FILE']  = glob.glob(self.bcs_dir_land + 'lnfm_clim_*.data' )[0]
                    self.rqdExeInp['NDVI_FILE']  = glob.glob(self.bcs_dir_land + 'ndvi_clim_*.data' )[0]
                    self.rqdExeInp['NIRDF_FILE'] = glob.glob(self.bcs_dir_land + 'nirdf_*.dat'      )[0]
                    self.rqdExeInp['VISDF_FILE'] = glob.glob(self.bcs_dir_land + 'visdf_*.dat'      )[0]

        if 'GRIDNAME' not in self.rqdExeInp :
            tmptile = os.path.realpath(self.rqdExeInp['TILING_FILE'])
            extension = os.path.splitext(tmptile)[1]
            if extension == '.domain':
              extension = os.path.splitext(tmptile)[0]
            gridname_ =''
            if extension == '.til':
               gridname_ = linecache.getline(tmptile, 3).strip()
            else:
               nc_file = netCDF4.Dataset(tmptile,'r')
               gridname_ = nc_file.getncattr('Grid_Name')
            # in case it is an old name: SMAP-EASEvx-Mxx
            gridname_ = gridname_.replace('SMAP-','').replace('-M','_M')
            self.rqdExeInp['GRIDNAME']  = gridname_
               
        if 'LSM_CHOICE' not in self.rqdExeInp:
            self.rqdExeInp['LSM_CHOICE'] = 1

        if  int(self.rqdExeInp['LSM_CHOICE']) == 1 :
            self.catch = 'catch'
        if  int(self.rqdExeInp['LSM_CHOICE']) == 2 :
            self.catch = 'catchcnclm40'

        if self.with_land: 
           assert int(self.rqdExeInp['LSM_CHOICE']) <= 2, "\nLSM_CHOICE=3 (Catchment-CN4.5) is no longer supported. Please set LSM_CHOICE to 1 (Catchment) or 2 (Catchment-CN4.0)"

        if 'POSTPROC_HIST' not in self.rqdExeInp:
            self.rqdExeInp['POSTPROC_HIST'] = 0

        if 'RUN_IRRIG' not in self.rqdExeInp:
            self.rqdExeInp['RUN_IRRIG'] = 0

        if 'AEROSOL_DEPOSITION' not in self.rqdExeInp:
            self.rqdExeInp['AEROSOL_DEPOSITION'] = 0
        # default is global
        _domain_dic=OrderedDict()
        _domain_dic['MINLON']=-180.
        _domain_dic['MAXLON']= 180.
        _domain_dic['MINLAT']= -90.
        _domain_dic['MAXLAT']=  90.
        _domain_dic['EXCLUDE_FILE']= "''"
        _domain_dic['INCLUDE_FILE']= "''"

        for key,val in _domain_dic.items() :
            if key in self.rqdExeInp :
                _domain_dic[key]= self.rqdExeInp[key]
        self.domain_def = tempfile.NamedTemporaryFile(mode='w', delete=False)
        self.domain_def.write('&domain_inputs\n')
        for key,val in _domain_dic.items() :
            keyn=(key+"  = ").ljust(16)
            valn = str(val)
            if '_FILE' in key:
               self.domain_def.write(keyn+ "'"+valn+"'"+'\n')
            else :
               self.domain_def.write(keyn+ valn +'\n')
        self.domain_def.write('/\n')
        self.domain_def.close()

        # make sure bcs files exist
        if self.rqdExeInp['RESTART'].isdigit() and self.with_land :
            if int(self.rqdExeInp['RESTART']) >= 1 :
                y4m2='Y%4d/M%02d' % (self.begDates[0].year, self.begDates[0].month)
                y4m2d2_h2m2='%4d%02d%02d_%02d%02d' % (self.begDates[0].year, self.begDates[0].month,
                                                      self.begDates[0].day,self.begDates[0].hour,self.begDates[0].minute)
                tmpFile=self.rqdExeInp['RESTART_ID']+'.'+self.catch+'_internal_rst.'+y4m2d2_h2m2
                tmpRstDir=self.rqdExeInp['RESTART_PATH']+'/'.join([self.rqdExeInp['RESTART_ID'],'output',
                                                                   self.rqdExeInp['RESTART_DOMAIN'],'rs',self.ensdirs[0],y4m2])
                catchRstFile=tmpRstDir+'/'+tmpFile

                assert os.path.isfile(catchRstFile), self.catch+'_internal_rst file [%s] does not exist!' %(catchRstFile)
                self.in_rstfile = catchRstFile

                if int(self.rqdExeInp['RESTART']) == 1 :
                    tmpFile=self.rqdExeInp['RESTART_ID']+'.vegdyn_internal_rst'
                    tmpRstDir=self.rqdExeInp['RESTART_PATH']+'/'.join([self.rqdExeInp['RESTART_ID'],'output',
                                                                       self.rqdExeInp['RESTART_DOMAIN'],'rs',self.ensdirs[0]])
                    vegdynRstFile=tmpRstDir+'/'+tmpFile
                    if not os.path.isfile(vegdynRstFile):
                        assert int(self.rqdExeInp['RST_FROM_GLOBAL']) == 1, 'restart from LDASsa should be global'

                    tmpFile=self.rqdExeInp['RESTART_ID']+'.landpert_internal_rst.'+y4m2d2_h2m2
                    tmpRstDir=self.rqdExeInp['RESTART_PATH']+'/'.join([self.rqdExeInp['RESTART_ID'],'output',
                                                                       self.rqdExeInp['RESTART_DOMAIN'],'rs',self.ensdirs[0],y4m2])
                    landpertRstFile=tmpRstDir+'/'+tmpFile
                    if ( os.path.isfile(landpertRstFile)) :
                        self.has_geos_pert = True

            elif (int(self.rqdExeInp['RESTART']) == 0) :
              if (self.catch == 'catch'):
                self.in_rstfile = '/discover/nobackup/projects/gmao/ssd/land/l_data/LandRestarts_for_Regridding' \
                                  '/Catch/M09/20170101/catch_internal_rst'
                self.in_tilefile = '/discover/nobackup/projects/gmao/ssd/land/l_data/geos5/bcs/CLSM_params' \
                                   '/mkCatchParam_SMAP_L4SM_v002/SMAP_EASEv2_M09/SMAP_EASEv2_M09_3856x1624.til'
              elif (self.catch == 'catchcnclm40'):
                self.in_rstfile  = '/discover/nobackup/projects/gmao/ssd/land/l_data/LandRestarts_for_Regridding' \
                                   '/CatchCN/M36/20150301_0000/catchcnclm40_internal_dummy'
                self.in_tilefile = '/discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Heracles-NL/SMAP_EASEv2_M36/SMAP_EASEv2_M36_964x406.til'
              elif (self.catch == 'catchcnclm45'):
                self.in_rstfile  = '/discover/nobackup/projects/gmao/ssd/land/l_data/LandRestarts_for_Regridding' \
                                   '/CatchCN/M36/19800101_0000/catchcnclm45_internal_dummy'
                self.in_tilefile = '/discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Icarus-NLv3/Icarus-NLv3_EASE/SMAP_EASEv2_M36/SMAP_EASEv2_M36_964x406.til'
              else:
                sys.exit('need to provide at least dummy files')

        # DEAL WITH mwRTM input from exec
        self.assim = True if self.rqdExeInp.get('LAND_ASSIM', 'NO').upper() == 'YES' and self.with_land else False
        # verify mwrtm file
        if 'MWRTM_PATH' in self.rqdExeInp and self.with_land :
           self.rqdExeInp['MWRTM_PATH'] = self.rqdExeInp['MWRTM_PATH']+'/'+ self.rqdExeInp['BCS_RESOLUTION']+'/'
           mwrtm_param_file_ = self.rqdExeInp['MWRTM_PATH']+'mwRTM_param.nc4'
           vegopacity_file_  = self.rqdExeInp['MWRTM_PATH']+'vegopacity.bin'
           if os.path.isfile(mwrtm_param_file_) :
              self.has_mwrtm = True
              self.mwrtm_file = mwrtm_param_file_
           else :
              assert not mwrtm_param_file_.strip(),  ' MWRTM_PATH: %s should contain mwRTM_param.nc4'% self.rqdExeInp['MWRTM_PATH']
              del self.rqdExeInp['MWRTM_PATH']
           if os.path.isfile(vegopacity_file_) :
              self.has_vegopacity = True
              self.rqdExeInp['VEGOPACITY_FILE'] = vegopacity_file_


        # ------
        # Read rm input file
        # Read (and pop from inpfile) the input required fields in to
        # self.rqdRmInp. Fields left in inpDictFromFile are then
        # read in to self.optRmInp
        # ------
        # re-using inpDictFromFile

        if self.ladas_cpl == 0 :
            inpDictFromFile = self._parseInputFile(cmdLineArgs['batinpfile'])
            # REQUIRED inputs
            for key in rqdRmInpKeys:
                self.rqdRmInp[key] = inpDictFromFile.pop(key)

            # checks on rqd rm inputs
            ## account and walltime should exist
            assert self.rqdRmInp['account']
            assert self.rqdRmInp['walltime']
            ## ntasks_model is a +ve integer
            _ntasks = int(self.rqdRmInp['ntasks_model'])
            assert _ntasks>0
            self.rqdRmInp['ntasks_model'] = _ntasks
            _ntasks = None

            # OPTIONAL inputs
            for key in inpDictFromFile:
                assert key in optSlurmInpKeys, \
                    'unknown resource manager key [%s]' % key
                self.optRmInp[key] = inpDictFromFile[key]
        else :
            self.rqdRmInp['account']         = cmdLineArgs['account']
            self.rqdRmInp['walltime']        = "01:00:00" 
            self.rqdRmInp['ntasks_model']    = 120


        # print rm inputs
        if self.verbose:
            print ('\n\nRequired inputs for resource manager:')
            _printdict(self.rqdRmInp)
            print ('\n\nOptional inputs for resource manager:')
            _printdict(self.optRmInp)
            print ('\n\n')
            
        # ------
        # set top level directories
        #     rundir, inpdir, outdir, blddir
        # executable
        #     exefyl
        # ------

        self.bindir     = os.path.dirname(os.path.realpath(__file__))
        self.blddir     = self.bindir.rsplit('/',1)[0]
        exefyl          = '/bin/GEOSldas.x'
        tmp_execfyl     = self.blddir   + exefyl
        assert os.path.isfile(tmp_execfyl),\
            'Executable [%s] does not exist!' % tmp_execfyl
        self.expdir     = self.exphome  + '/' + self.rqdExeInp['EXP_ID']
        self.rundir     = self.expdir   + '/run'
        self.inpdir     = self.expdir   + '/input'
        self.outdir     = self.expdir   + '/output'
        self.scratchdir = self.expdir   + '/scratch'
        self.blddirLn   = self.expdir   + '/build'
        self.out_path   = self.outdir   + '/'+self.rqdExeInp['EXP_DOMAIN']
        self.bcsdir     = self.outdir   + '/'+self.rqdExeInp['EXP_DOMAIN']+'/rc_out/'
        self.rstdir     = self.outdir   + '/'+self.rqdExeInp['EXP_DOMAIN']+'/rs/'
        self.exefyl     = self.blddirLn + exefyl

        # default is set to 0 ( no output server)
        if 'oserver_nodes' not in self.optRmInp :
           self.optRmInp['oserver_nodes'] = 0

        if (int(self.optRmInp['oserver_nodes']) >=1) :
           self.rqdExeInp['WRITE_RESTART_BY_OSERVER'] = "YES"
           # set default for now
           if 'writers-per-node' not in self.optRmInp:
              self.optRmInp['writers-per-node'] = 5
        else:
           self.optRmInp['writers-per-node'] = 0

           
    # end __init__
    
    # -----------------------------------------------------------------------------------
           
    def _parseInputFile(self, inpfile):
        """
        Private method: parse input file and return a dict of options
        Input: input file
        Output: dict
        """

        inpdict = OrderedDict()
        errstr = "line [%d] of [%s] is not in the form 'key: value'"

        # determine which default values to pick from GEOS_SurfaceGridComp.rc
        if self.ladas_cpl == 0 :
            use_rc_defaults = 'GEOSldas=>'    # use defaults for LDAS
        else :
            use_rc_defaults = 'GEOSagcm=>'    # use defaults for AGCM
        
        fin = open(inpfile, 'r')
        linenum = 0
        for line in fin:
            linenum += 1
            line = line.strip()
            # blank line
            if not line:
                continue
            if '"GEOSagcm=>"' in line:        # echo lines that contain "GEOSagcm=>" (w/ quotation marks) [GEOS_SurfaceGridComp.rc]
                continue
            if '"GEOSldas=>"' in line:        # echo lines that contain "GEOSldas=>" (w/ quotation marks) [GEOS_SurfaceGridComp.rc]
                continue
            # get 'GEOSldas=>' or 'GEOSagcm=>' defaults in GEOS_SurfaceGridComp.rc
            if use_rc_defaults in line:
                line = line.split(use_rc_defaults)[1]
            # handle comments
            position = line.find('#')
            if position==0: # comment line
                continue
            if position>0:  # strip out comment
                line = line[:position]
            # we expect a line to be of the form
            # key : value
            assert ':' in line, errstr % (linenum, inpfile)

            key, val = line.split(':',1)
            key = key.strip()
            val = val.strip()
            if not key or not val:
                print ("WARNING: " + errstr % (linenum, inpfile))
                continue
                #raise Exception(errstr % (linenum, inpfile))
            if key in inpdict:
                raise Exception('Duplicate key [%s] in [%s]' % (key, inpfile))
            inpdict[key] = val.strip()
        fin.close()

        return inpdict

    # -----------------------------------------------------------------------------------
           
    def _mkdir_p(self,path):
        """
        Private method: implement 'mkdir -p' functionality
        """

        if os.path.isdir(path):
            return
        else:
            os.makedirs(path)

    # -----------------------------------------------------------------------------------
           
    def createDirStructure(self):
        """
        Create required dir structure
        """

        status = False

        # shorthands
        _nens = self.nens

        # run/inp/wrk dirs
        self._mkdir_p(self.exphome+'/'+self.rqdExeInp['EXP_ID'])
        self._mkdir_p(self.rundir)
        self._mkdir_p(self.inpdir)
        self._mkdir_p(self.outdir)
        self._mkdir_p(self.scratchdir)

        #-start-shorthand-function-
        def _getDirName(outtyp, ensdir, yyyymm):
            return '/'.join([
                    self.outdir,
                    self.rqdExeInp['EXP_DOMAIN'],
                    outtyp, # ana/cat/rs/rc_out
                    ensdir,
                    yyyymm
                    ])
        #-end-shorthand-function-

        # met forcing dir
        myMetDir = self.inpdir + '/met_forcing'
        self._mkdir_p(myMetDir)

        # ensxxxx directories
        nSegments = self.nSegments
        for iseg in range(nSegments):
            _start = self.begDates[iseg]
            _end = self.endDates[iseg]

            # Yyyyy/Mmm between StartDateTime and EndDateTime
            newDate = _start
            y4m2_list = [('Y%4d/M%02d' % (newDate.year, newDate.month))]
            while newDate<_end:
                newDate += relativedelta(months=1)
                y4m2_list.append('Y%4d/M%02d' % (newDate.year, newDate.month))

            # ExpDomain/ana/, /cat/ directories
            for ensdir in self.ensdirs_avg:
                for y4m2 in y4m2_list:
                    self._mkdir_p(_getDirName('ana', ensdir, y4m2))
                    self._mkdir_p(_getDirName('cat', ensdir, y4m2))

            # ExpDomain/rs/ directories
            for ensdir in self.ensdirs:
                for y4m2 in y4m2_list:
                    self._mkdir_p(_getDirName('rs', ensdir, y4m2))

            # ExpDomain/rc_out/ - only for _start
            self._mkdir_p(_getDirName('rc_out', '', y4m2_list[0]))

        # restart dir
        self._mkdir_p(self.inpdir + '/restart')

        status = True
        return status

    # -----------------------------------------------------------------------------------
           
    # create links to BCs, restarts, met forcing, ...
    def createLnRstBc(self) :
        # link bld dir
        status = False

        _nens = self.nens

        os.symlink(self.blddir, self.blddirLn)

        # met forcing dir
        self.ensemble_forcing = True if self.rqdExeInp.get('ENSEMBLE_FORCING', 'NO').upper() == 'YES' else False

        myMetPath =''
        for _i in range(self.first_ens_id, _nens + self.first_ens_id) :
           str_ens = ''
           if ( _nens != 1 and self.ensemble_forcing):
              str_ens = '%03d'%(_i)
           metpath = self.rqdExeInp['MET_PATH'].rstrip('/')+str_ens
           myMetDir = self.inpdir + '/met_forcing'
           myMetPath = myMetDir + '/' + metpath.split('/')[-1]
           os.symlink(metpath, myMetPath)
           # update 'met_path' to use relative path from outdir
           if ( not self.ensemble_forcing):
              break
        if ( _nens !=1 and self.ensemble_forcing) :
           # replace last three character with '%s"
           self.rqdExeInp['MET_PATH'] = os.path.relpath(myMetPath, self.rundir)[:-3]+'%s'
        else:
           self.rqdExeInp['MET_PATH'] = os.path.relpath(myMetPath, self.rundir)

        # update tile file
        tile= self.rqdExeInp['TILING_FILE']
        short_tile= os.path.basename(self.rqdExeInp['TILING_FILE'])
        newtile = self.bcsdir+'/'+short_tile
        shutil.copy(tile, newtile)
        tile=newtile
        # if three extra lines exist, remove them and save it to inputdir

        print ('\nCorrect the tile file if it is an old EASE tile format... \n')
        EASEtile=self.bcsdir+'/MAPL_'+short_tile
        cmd = self.bindir + '/preprocess_ldas.x correctease  '+ tile + ' '+ EASEtile
        print ("cmd:   " + cmd)

        sp.call(shlex.split(cmd))

        if os.path.isfile(EASEtile) :
            #update tile file name
            short_tile ='MAPL_'+short_tile
            tile=EASEtile
        # setup BC files

        catchment_def = self.rqdExeInp['CATCH_DEF_FILE']
        exp_id = self.rqdExeInp['EXP_ID']

        _start = self.begDates[0]
        _y4m2d2h2m2 ='%4d%02d%02d%02d%02d' % (_start.year, _start.month,_start.day,_start.hour,_start.minute)

        dzsf = '50.0'
        if 'SURFLAY' in self.rqdExeInp :
           dzsf = self.rqdExeInp['SURFLAY']

        # These are dummy values for *cold* restart:
        wemin_in  = '13'  #  WEmin input/output for scale_catch(cn),
        wemin_out = '13'  #
        if 'WEMIN_IN' in self.rqdExeInp :
           wemin_in = self.rqdExeInp['WEMIN_IN']
        if 'WEMIN_OUT' in self.rqdExeInp :
           wemin_out = self.rqdExeInp['WEMIN_OUT']

        tmp_f2g_file = tempfile.NamedTemporaryFile(delete=False) 
        cmd = self.bindir +'/preprocess_ldas.x c_f2g ' + tile + ' ' + self.domain_def.name + ' '+ self.out_path + ' ' + catchment_def + ' ' + exp_id + ' ' + _y4m2d2h2m2 + ' '+ dzsf + ' ' + tmp_f2g_file.name + ' ' + '_'.join(self.tile_types) 

        print ('Creating f2g file if necessary: '+ tmp_f2g_file.name +'....\n')
        print ("cmd:   " + cmd)
        sp.call(shlex.split(cmd)) 
        # check if it is local or global
        if os.path.getsize(tmp_f2g_file.name) !=0 :
           self.isZoomIn= True
        #os.remove(self.domain_def.name)

        # update tile domain
        if self.isZoomIn:
            newZoominTile = tile+'.domain'
            print ("\nCreating local tile file :"+ newZoominTile)
            print ("\nAdding 1000 to type of tiles to be excluded from domain...\n")
            cmd = self.bindir +'/preprocess_ldas.x zoomin_tile ' + tile + ' ' + newZoominTile + ' '+ tmp_f2g_file.name
            print ("cmd:   " + cmd)
            sp.call(shlex.split(cmd))
            short_tile=short_tile +'.domain'
            tile = newZoominTile

        myTile=self.inpdir+'/tile.data'
        os.symlink(tile,myTile)

        if self.with_land:
           bcs=[self.rqdExeInp['GRN_FILE'],
                self.rqdExeInp['LAI_FILE'],
                self.rqdExeInp['NDVI_FILE'],
                self.rqdExeInp['NIRDF_FILE'],
                self.rqdExeInp['VISDF_FILE'] ]
           if (self.rqdExeInp['LNFM_FILE'] != ''):
              bcs += [self.rqdExeInp['LNFM_FILE']]
           if (self.has_vegopacity):
              bcs += [self.rqdExeInp['VEGOPACITY_FILE']]
           bcstmp=[]
           for bcf in bcs :
               shutil.copy(bcf, self.bcsdir+'/')
               bcstmp=bcstmp+[self.bcsdir+'/'+os.path.basename(bcf)]
           bcs=bcstmp

           if self.isZoomIn:
              print ("Creating the boundary files for the simulation domain...\n")
              bcs_tmp=[]
              for bcf in bcs :
                  cmd = self.bindir +'/preprocess_ldas.x zoomin_bc ' + bcf + ' '+ bcf+'.domain' + ' '+ tmp_f2g_file.name
                  print ("cmd:   " + cmd)
                  sp.call(shlex.split(cmd))
                  bcs_tmp=bcs_tmp+[bcf+'.domain']
              bcs=bcs_tmp


           # link BC
           print ("linking bcs...")
           bcnames=['green','lai','ndvi','nirdf','visdf']
           if (self.rqdExeInp['LNFM_FILE'] != ''):
              bcnames += ['lnfm']
           if (self.has_vegopacity):
              bcnames += ['vegopacity']
           for bcln,bc in zip(bcnames,bcs) :
              myBC=self.inpdir+'/'+bcln+'.data'
              os.symlink(bc,myBC)

           if ("catchcn" in self.catch):
              os.symlink(self.bcs_dir_landshared + 'CO2_MonthlyMean_DiurnalCycle.nc4', \
                          self.inpdir+'/CO2_MonthlyMean_DiurnalCycle.nc4')

        # create and link restart
        print ("Creating and linking restart...")
        _start = self.begDates[0]

        y4m2='Y%4d/M%02d'%(_start.year, _start.month)
        y4m2d2_h2m2 ='%4d%02d%02d_%02d%02d' % (_start.year, _start.month,_start.day,_start.hour,_start.minute)

        myRstDir = self.inpdir + '/restart/'

        rstpath = self.rqdExeInp['RESTART_PATH']+ \
                  self.rqdExeInp['RESTART_ID'] + \
                  '/output/'+self.rqdExeInp['RESTART_DOMAIN']+'/rs/'
        rcoutpath = self.rqdExeInp['RESTART_PATH']+ \
                  self.rqdExeInp['RESTART_ID'] + \
                  '/output/'+self.rqdExeInp['RESTART_DOMAIN']+'/rc_out/'

        # pass into remap_config_ldas
        exp_id    = self.rqdExeInp['EXP_ID']
        RESTART_str  = str(self.rqdExeInp['RESTART'])
        YYYYMMDD  = '%4d%02d%02d' % (_start.year, _start.month,_start.day)
        YYYYMMDDHH= '%4d%02d%02d%02d' % (_start.year, _start.month,_start.day, _start.hour)
        rstid     = self.rqdExeInp['RESTART_ID']
        rstdomain = self.rqdExeInp['RESTART_DOMAIN']
        rstpath0  = self.rqdExeInp['RESTART_PATH']

        # just copy the landassim pert seed if it exists
        for iens in range(self.nens) :
           _ensdir = self.ensdirs[iens]
           _ensid  = self.ensids[iens]
           landassim_seeds = rstpath + _ensdir + '/' + y4m2+'/' + rstid + '.landassim_obspertrseed_rst.'+y4m2d2_h2m2
           if os.path.isfile(landassim_seeds) and self.assim :
               _seeds = self.rstdir + _ensdir + '/' + y4m2+'/' + exp_id + '.landassim_obspertrseed_rst.'+y4m2d2_h2m2
               shutil.copy(landassim_seeds, _seeds)
               os.symlink(_seeds, myRstDir+ '/landassim_obspertrseed'+ _ensid +'_rst')
               self.has_landassim_seed = True
        mk_outdir = self.exphome+'/'+exp_id+'/mk_restarts/'

        if (RESTART_str != '1' and (self.with_land or self.with_landice)):
           bcs_path = self.rqdExeInp['BCS_PATH']
           while bcs_path[-1] == '/' : bcs_path = bcs_path[0:-1]
           bc_base    = os.path.dirname(bcs_path)
           bc_version = os.path.basename(bcs_path)

           remap_tpl = os.path.dirname(os.path.realpath(__file__)) + '/remap_params.tpl'
           config = yaml_to_config(remap_tpl)

           config['slurm_pbs']['account'] = self.rqdRmInp['account']
           config['slurm_pbs']['qos'] = 'debug'

           config['input']['surface']['catch_tilefile']  = self.in_tilefile
           config['input']['shared']['expid']            = self.rqdExeInp['RESTART_ID']
           config['input']['shared']['yyyymmddhh']       = YYYYMMDDHH
           if RESTART_str != 'M':
             config['input']['shared']['rst_dir']        = os.path.dirname(self.in_rstfile)+'/'
           config['input']['surface']['wemin']           = wemin_in
           config['input']['surface']['catch_model']     = self.catch

           config['output']['shared']['out_dir']         = mk_outdir
           config['output']['surface']['catch_remap']    = True
           config['output']['surface']['catch_tilefile'] = self.rqdExeInp['TILING_FILE']
           config['output']['shared']['bc_base']         = bc_base
           config['output']['shared']['bc_version']      = bc_version
           config['output']['surface']['EASE_grid']      = self.rqdExeInp['BCS_RESOLUTION']

           config['output']['shared']['expid']           = self.rqdExeInp['EXP_ID']
           config['output']['surface']['surflay']        = dzsf
           config['output']['surface']['wemin']          = wemin_out

           if RESTART_str == "M" : # restart from merra2
             yyyymm = int(YYYYMMDDHH[0:6])
             merra2_expid = "d5124_m2_jan10"
             if yyyymm < 197901 :
                exit("Error. MERRA-2 data < 1979 not available\n")
             elif (yyyymm < 199201):
                merra2_expid = "d5124_m2_jan79"
             elif (yyyymm < 200106):
                merra2_expid = "d5124_m2_jan91"
             elif (yyyymm < 201101):
                merra2_expid = "d5124_m2_jan00"
             elif (yyyymm < 202106):
                merra2_expid = "d5124_m2_jan10"
             # There was a rewind in MERRA2 from Jun 2021 to Sept 2021
             elif (yyyymm < 202110):
                merra2_expid = "d5124_m2_jun21"
             config['input']['shared']['expid']      = merra2_expid
             config['input']['shared']['rst_dir']    = mk_outdir+ '/merra2_tmp_'+ YYYYMMDDHH
             config['input']['surface']['wemin']     = 26
             config['input']['shared']['bc_base']    = '/discover/nobackup/projects/gmao/bcs_shared/fvInput/ExtData/esm/tiles'
             config['input']['shared']['bc_version'] = 'GM4'
             config['input']['shared']['agrid']      = 'C180'
             config['input']['shared']['ogrid']      = '1440x720'
             config['input']['shared']['omodel']     = 'data'
             config['input']['shared']['MERRA-2']    = True
             config['input']['surface']['catch_tilefile'] = '/discover/nobackup/projects/gmao/bcs_shared/fvInput/ExtData/esm/tiles/GM4/geometry/CF0180x6C_DE1440xPE0720/CF0180x6C_DE1440xPE0720-Pfafstetter.til'

           if self.with_land:
             catch_obj = catchANDcn(config_obj = config)
             catch_obj.remap()
           if self.with_landice:
             config['output']['surface']['remap_water'] = True
             config['input']['surface']['zoom'] = '2'
             landice_obj = lake_landice_saltwater(config_obj = config)
             landice_obj.remap()

        #for ens in self.ensdirs :
        catchRstFile0  = ''
        vegdynRstFile0 = ''
        landiceRstFile0 = ''
        for iens in range(self.nens) :
            ensdir   = self.ensdirs[iens]
            ensid    = self.ensids[iens]
            myCatchRst = myRstDir+'/'+self.catch +ensid +'_internal_rst'
            myLandiceRst = myRstDir+'/'+ 'landice' +ensid +'_internal_rst'
            myVegRst   = myRstDir+'/'+'vegdyn'+ensid +'_internal_rst'
            myPertRst  = myRstDir+'/'+ 'landpert' +ensid +'_internal_rst'

            catchRstFile  = ''
            vegdynRstFile = ''
            pertRstFile   = ''
            print ("restart: " + self.rqdExeInp['RESTART'])

            if self.rqdExeInp['RESTART'].isdigit() :

                if int(self.rqdExeInp['RESTART']) == 0 or int(self.rqdExeInp['RESTART']) == 2 :
                    vegdynRstFile = glob.glob(self.bcs_dir_land + 'vegdyn_*.dat')[0]
                    catchRstFile  = glob.glob(self.exphome+'/'+exp_id+'/mk_restarts/*'+self.catch+'_internal_rst.'+YYYYMMDD+'*')[0]
                else : # RESTART == 1
                    catchRstFile = rstpath+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['RESTART_ID']+'.'+self.catch+'_internal_rst.'+y4m2d2_h2m2
                    vegdynRstFile= rstpath+ensdir +'/'+self.rqdExeInp['RESTART_ID']+ '.vegdyn_internal_rst'
                    if not os.path.isfile(vegdynRstFile): # no vegdyn restart from LDASsa
                       if not os.path.isfile(vegdynRstFile0):
                         vegdynRstFile = glob.glob(self.bcs_dir_land + 'vegdyn_*.dat')[0]
            else :
                vegdynRstFile = glob.glob(self.bcs_dir_land + 'vegdyn_*.dat')[0]
                if self.with_land:
                  catchRstFile  = glob.glob(self.exphome+'/'+exp_id+'/mk_restarts/*'+self.catch+'_internal_rst.'+YYYYMMDD+'*')[0]

            # catchment restart file
            if os.path.isfile(catchRstFile) and self.with_land :
                catchLocal = self.rstdir+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.'+self.catch+'_internal_rst.'+y4m2d2_h2m2
                if self.isZoomIn :
                    print( "Creating local catchment restart file... \n")
                    cmd=self.bindir +'/preprocess_ldas.x zoomin_catchrst '+ catchRstFile +' ' + catchLocal + ' '+ tmp_f2g_file.name
                    print ("cmd:  "+cmd)
                    sp.call(shlex.split(cmd))
                else :
                    shutil.copy(catchRstFile,catchLocal)

                catchRstFile = catchLocal

                if '0000' in ensdir :
                    catchRstFile0 = catchRstFile
            else : # re-use 0000 catch file
                catchRstFile = catchRstFile0

            # vegdyn restart file
            if os.path.isfile(vegdynRstFile) and self.with_land :
                vegdynLocal = self.rstdir+ensdir +'/'+self.rqdExeInp['EXP_ID']+'.vegdyn_internal_rst'
                if self.isZoomIn :
                    print ("Creating the local veg restart file... \n")
                    cmd=self.bindir + '/preprocess_ldas.x zoomin_vegrst '+ vegdynRstFile +' ' + vegdynLocal + ' '+ tmp_f2g_file.name
                    print ("cmd:   " + cmd)
                    sp.call(shlex.split(cmd))
                else :
                    shutil.copy(vegdynRstFile,vegdynLocal)

                vegdynRstFile = vegdynLocal

                if '0000' in ensdir :
                    vegdynRstFile0 = vegdynRstFile
            else :
                vegdynRstFile = vegdynRstFile0

            landiceRstFile = ''
            if self.with_landice :
               if self.rqdExeInp['RESTART'].isdigit():
                  if int(self.rqdExeInp['RESTART']) == 0 or int(self.rqdExeInp['RESTART']) == 2 :
                     print("RESTART=0 and RESTART=2 not supported for landice tiles. Please use RESTART=M (MERRA-2).")
                  landiceRstFile = rstpath+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['RESTART_ID']+'.'+'landice_internal_rst.'+y4m2d2_h2m2
               else:
                  landiceRstFile = glob.glob(self.exphome+'/'+exp_id+'/mk_restarts/*'+'landice_internal_rst.'+YYYYMMDD+'*')[0]

            if os.path.isfile(landiceRstFile) :
               landiceLocal = self.rstdir+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.landice_internal_rst.'+y4m2d2_h2m2
               if self.isZoomIn :
                  print ("Creating zoom-in of landice restart file... \n")
                  cmd=self.bindir + '/preprocess_ldas.x zoomin_landicerst '+ landiceRstFile +' ' + landiceLocal + ' '+ tmp_f2g_file.name
                  print ("cmd:   " + cmd)
                  sp.call(shlex.split(cmd))
               else :
                  shutil.copy(landiceRstFile,landiceLocal)

               landiceRstFile = landiceLocal

               if '0000' in ensdir :
                  landiceRstFile0 = landiceRstFile
            else :
                landiceRstFile = landiceRstFile0
            
            if (self.has_geos_pert and self.perturb == 1) :
                pertRstFile = rstpath+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['RESTART_ID']+'.landpert_internal_rst.'+y4m2d2_h2m2
                pertLocal   = self.rstdir+ensdir +'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.landpert_internal_rst.'+y4m2d2_h2m2
                shutil.copy(pertRstFile,pertLocal)
                pertRstFile = pertLocal

            if self.with_land :
               print ('catchRstFile: '  +  catchRstFile)
               print ('vegdynRstFile: ' +  vegdynRstFile)
               os.symlink(catchRstFile,   myCatchRst)
               os.symlink(vegdynRstFile,  myVegRst)
            if self.with_landice :
               print("link landice restart: " + myLandiceRst)
               os.symlink(landiceRstFile, myLandiceRst)
            if ( self.has_geos_pert and  self.perturb == 1 ):
               os.symlink(pertRstFile,    myPertRst)

        # catch_param restar file
        catch_param_file = self.bcsdir+'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.ldas_catparam.'+y4m2d2_h2m2+'z.bin'
        if self.with_land:
           assert os.path.isfile(catch_param_file), "need catch_param file %s" % catch_param_file

        if self.has_mwrtm :
           mwRTMRstFile = self.mwrtm_file
           mwRTMLocal = self.bcsdir+'/'+ y4m2+'/'+self.rqdExeInp['EXP_ID']+'.ldas_mwRTMparam.'+y4m2d2_h2m2+'z.nc4'
           if self.isZoomIn :
              print ("Creating the local mwRTM restart file... \n")
              cmd= self.bindir +'/preprocess_ldas.x zoomin_mwrtmrst '+ mwRTMRstFile +' ' + mwRTMLocal + ' '+ tmp_f2g_file.name

              print ("cmd:   " + cmd)
              sp.call(shlex.split(cmd))
           else :
              shutil.copy(mwRTMRstFile,mwRTMLocal)

           mwRTMRstFile = mwRTMLocal
           mymwRTMRst = myRstDir+'/mwrtm_param_rst'
           os.symlink(mwRTMRstFile,  mymwRTMRst)

        # update 'restart_path' to use relative path from outdir
        print ("Updating restart path...")
        self.rqdExeInp['RESTART_PATH'] = myRstDir
        #if os.path.isfile(tmp_f2g_file.name):
        #   os.remove(tmp_f2g_file.name)
        status = True
        return status

    # -----------------------------------------------------------------------------------
           
    def createRCFiles(self):
        """
        (1) get resource files form DEFAULT rc files from /etc
        (2) update from customed rc files
        (2) write rc files to the run directory
        """

        status = False

        for mydir in [self.blddirLn, self.rundir]:
            assert os.path.isdir(mydir), \
                'dir [%s] does not exist!' % mydir

        if self.ladas_cpl == 0:            
            # copy ldas_setup exeinp and batinp input files to rundir (for the record)
            # if a file w/ the same name already exists at rundir
            # append 1,2,3 etc, to the filename
            ## exe inp file
            exefilename = self.exeinpfile.rstrip('/').split('/')[-1]
            newfilename = exefilename
            _nens = self.nens
            ctr = 0
            while os.path.isfile(self.rundir+'/'+newfilename):
                ctr += 1
                newfilename = exefilename + '.%d' % ctr
            shutil.copy(self.exeinpfile, self.rundir+'/'+newfilename)
            ## bat inp file
            batfilename = self.batinpfile.rstrip('/').split('/')[-1]
            newfilename = batfilename
            ctr = 0
            while os.path.isfile(self.rundir+'/'+newfilename):
                ctr += 1
                newfilename = batfilename + '.%d' % ctr
            shutil.copy(self.batinpfile, self.rundir+'/'+newfilename)

        # -----------------------------------    
            
        etcdir = self.blddirLn + '/etc'

        #defalt nml
        default_nml = glob.glob(etcdir+'/LDASsa_DEFAULT_inputs_*.nml')
        for nmlfile in default_nml:
            shortfile=self.rundir+'/'+nmlfile.split('/')[-1]
            shutil.copy2(nmlfile, shortfile)
        # special nml
        special_nml=[]
        if self.ladas_cpl > 0:
            special_nml= glob.glob(etcdir+'/LDASsa_SPECIAL_inputs_*.nml')
        else : 
            if 'NML_INPUT_PATH' in self.rqdExeInp : 
                special_nml = glob.glob(self.rqdExeInp['NML_INPUT_PATH']+'/LDASsa_SPECIAL_inputs_*.nml')

        for nmlfile in special_nml:
            shortfile=self.rundir+'/'+nmlfile.split('/')[-1]
            shutil.copy2(nmlfile, shortfile)   

        if self.ladas_cpl > 0:
            # edit resolution info in ensupd nml file
            sp.run(['sed', '-i', 's/<CFnnnn>/'+self.agcm_res+'/g', self.rundir+'/LDASsa_SPECIAL_inputs_ensupd.nml'])
 
        # get optimzed NX and IMS
        optimized_distribution_file = tempfile.NamedTemporaryFile(delete=False)
        print  ("Optimizing... decomposition of processes.... \n")
        cmd = self.bindir + '/preprocess_ldas.x optimize '+ self.inpdir+'/tile.data '+ str(self.rqdRmInp['ntasks_model']) + ' ' + optimized_distribution_file.name + ' ' + self.rundir + ' ' + '_'.join(self.tile_types)
        print ("cmd:   " + cmd)
        print ("IMS.rc or JMS.rc would be generated on " + self.rundir)
        sp.call(shlex.split(cmd))
        optinxny=self._parseInputFile(optimized_distribution_file.name)
        if (int(optinxny['NX']) == 1):
            if int(optinxny['NY']) != int(self.rqdRmInp['ntasks_model']):
                self.rqdRmInp['ntasks_model']=optinxny['NY']
                print ('adjust ntasks_model %d for cubed-sphere grid' % int(self.rqdRmInp['ntasks_model']))


        #os.remove(optimized_distribution_file.name)

        # DEFAULT rc files
        default_rc = glob.glob(etcdir+'/GEOSldas_*.rc')
        assert len(default_rc)==6
        print (default_rc)
        for rcfile in default_rc:
            shortfile=rcfile.rsplit('GEOSldas_',1)[1]
            print (shortfile + ' ' +  etcdir + ' ' + self.rundir)
            if shortfile =='HIST.rc':
                tmprcfile=self.rundir+'/HISTORY.rc'
                histrc_file=rcfile

                _file_found = False
                if 'HISTRC_FILE' in self.rqdExeInp :
                    _tmpfile = self.rqdExeInp['HISTRC_FILE'].replace("'",'').replace('"','')
                    if(os.path.isfile(_tmpfile)) :
                        _file_found = True
                    else :
                        assert not _tmpfile.strip(), "HISTRC_FILE: %s is NOT a file. " %_tmpfile

                if _file_found :
                    histrc_file = self.rqdExeInp['HISTRC_FILE']
                    shutil.copy2(histrc_file,tmprcfile)
                else :
                    shutil.copy2(histrc_file,tmprcfile)
                    GRID='EASE ' + self.rqdExeInp['GRIDNAME'] + ' ' +tmprcfile
                    if '-CF' in self.rqdExeInp['GRIDNAME'] :
                        GRID ='CUBE ' + self.rqdExeInp['GRIDNAME'] + ' ' +tmprcfile
                    _assim = '1' if self.assim else '0'
                    cmd =self.bindir +'/process_hist.csh '+ str(self.rqdExeInp['LSM_CHOICE']) + ' ' + str(self.rqdExeInp['AEROSOL_DEPOSITION']) + \
                          ' ' + GRID + ' ' + str(self.rqdExeInp['RUN_IRRIG']) + ' ' + _assim + ' '+ str(self.nens)
                    print(cmd)
                    #os.system(cmd)
                    sp.call(shlex.split(cmd))
                for line in fileinput.input(tmprcfile,inplace=True):
                    print (line.rstrip().replace('GEOSldas_expid',self.rqdExeInp['EXP_ID']))

            # if coupled land-atm DAS, always use either GEOSldas_HISTdet.rc or GEOSldas_HISTens.rc (depending on ladas_cpl)
            if ( shortfile =='HISTdet.rc' and self.ladas_cpl == 1 ) or ( shortfile =='HISTens.rc' and self.ladas_cpl == 2 ):
                tmprcfile=self.rundir+'/HISTORY.rc'
                histrc_file=rcfile
                shutil.copy2(rcfile, tmprcfile)
                for line in fileinput.input(tmprcfile,inplace=True):
                    print (line.rstrip().replace('GEOSldas_expid',self.rqdExeInp['EXP_ID']))
                for line in fileinput.input(tmprcfile,inplace=True):
                    print (line.rstrip().replace('GRIDNAME',self.rqdExeInp['GRIDNAME']))

            # just copy an empty ExtData.rc
            if shortfile=='ExtData.rc' :
                shutil.copy2(rcfile, self.rundir+'/'+shortfile)

            if shortfile == 'CAP.rc':
                tmprcfile = self.rundir+'/CAP.rc'
                shutil.copy2(rcfile,tmprcfile)

                _num_sgmt = int(self.rqdExeInp['NUM_SGMT'])

                for line in fileinput.input(tmprcfile,inplace=True):
                    print (line.rstrip().replace('JOB_SGMT:',self.job_sgmt[0]))
                for line in fileinput.input(tmprcfile,inplace=True):
                    print (line.rstrip().replace('NUM_SGMT:','NUM_SGMT: %d'% _num_sgmt))
                for line in fileinput.input(tmprcfile,inplace=True):
                    print (line.rstrip().replace('BEG_DATE:',self.begDates[ 0].strftime('BEG_DATE: %Y%m%d %H%M%S')))
                for line in fileinput.input(tmprcfile,inplace=True):
                    print (line.rstrip().replace('END_DATE:',self.endDates[-1].strftime('END_DATE: %Y%m%d %H%M%S')))

            if shortfile == 'LDAS.rc' :
                ldasrcInp = OrderedDict()
                # land default
                default_surfrcInp = self._parseInputFile(etcdir+'/GEOS_SurfaceGridComp.rc')
                for key,val in default_surfrcInp.items() :
                    ldasrcInp[key] = val

                # ldas default, may overwrite land default
                default_ldasrcInp = self._parseInputFile(rcfile)
                for key,val in default_ldasrcInp.items() :
                    ldasrcInp[key] = val

                # exeinp, may overwrite ldas default
                for key,val in self.rqdExeInp.items():
                    if key not in self.NoneLDASrcKeys:
                        ldasrcInp[key]= val

                # overide by optimized distribution
                #for key,val in optinxny.items():
                #    ldasrcInp[key]= val

                # create BC in rc file
                tmpl_ = ''
                if self.nens >1 :
                   tmpl_='%s'
                if self.perturb == 1:
                   ldasrcInp['PERTURBATIONS'] ='1'
                rstkey =[]
                rstval =[]
                if self.with_land :
                   bcval=['../input/green','../input/lai','../input/lnfm','../input/ndvi','../input/nirdf','../input/visdf']
                   bckey=['GREEN','LAI','LNFM','NDVI','NIRDF','VISDF']
                   for key, val in zip(bckey,bcval):
                      keyn = key+'_FILE'
                      valn = val+'.data'
                      ldasrcInp[keyn]= valn
                   if('catchcn' in self.catch):
                      ldasrcInp['CO2_MonthlyMean_DiurnalCycle_FILE']= '../input/CO2_MonthlyMean_DiurnalCycle.nc4'
                   else:
                      # remove catchcn-specific entries that do not apply to catch model
                      ldasrcInp.pop('DTCN',None)
                      ldasrcInp.pop('ATM_CO2',None)
                      ldasrcInp.pop('CO2',None)
                      ldasrcInp.pop('CO2_YEAR',None)
                      ldasrcInp.pop('PRESCRIBE_DVG',None)

                # create restart item in RC
                   catch_ = self.catch.upper()

                   if catch_+'_INTERNAL_RESTART_TYPE' in ldasrcInp :
                      # avoid duplicate
                      del ldasrcInp[ catch_ +'_INTERNAL_RESTART_TYPE']
                   if catch_+'_INTERNAL_CHECKPOINT_TYPE' in ldasrcInp :
                      # avoid duplicate
                      del ldasrcInp[ catch_ +'_INTERNAL_CHECKPOINT_TYPE']
                   if 'VEGDYN_INTERNAL_RESTART_TYPE' in ldasrcInp :
                      # avoid duplicate
                      del ldasrcInp['VEGDYN_INTERNAL_RESTART_TYPE']

                   rstkey.append(catch_)
                   rstkey.append('VEGDYN')
                   rstval.append(self.catch)
                   rstval.append('vegdyn')

                if self.with_landice:
                  rstkey.append('LANDICE')
                  rstval.append('landice')

                if self.has_mwrtm :
                   keyn='LANDASSIM_INTERNAL_RESTART_FILE'
                   valn='../input/restart/mwrtm_param_rst'
                   ldasrcInp[keyn]= valn
                if self.has_vegopacity :
                   keyn='VEGOPACITY_FILE'
                   valn='../input/vegopacity.data'
                   ldasrcInp[keyn]= valn

                if self.nens > 1 :
                   keyn='ENS_ID_WIDTH'
                   valn=str(self.ens_id_width)
                   ldasrcInp[keyn]= valn

                if self.has_landassim_seed and self.assim :
                   keyn='LANDASSIM_OBSPERTRSEED_RESTART_FILE'
                   valn='../input/restart/landassim_obspertrseed'+tmpl_+'_rst'
                   ldasrcInp[keyn]= valn

                if self.assim:
                   keyn='LANDASSIM_OBSPERTRSEED_CHECKPOINT_FILE'
                   valn='landassim_obspertrseed'+tmpl_+'_checkpoint'
                   ldasrcInp[keyn]= valn

                for key,val in zip(rstkey,rstval) :
                    keyn = key+ '_INTERNAL_RESTART_FILE'
                    valn = '../input/restart/'+val+tmpl_+'_internal_rst'
                    ldasrcInp[keyn]= valn

                # checkpoint file and its type
                if self.with_land :
                   keyn = catch_ + '_INTERNAL_CHECKPOINT_FILE'
                   valn = self.catch+tmpl_+'_internal_checkpoint'
                   ldasrcInp[keyn]= valn

                if self.with_landice :
                   keyn = 'LANDICE_INTERNAL_CHECKPOINT_FILE'
                   valn = 'landice'+tmpl_+'_internal_checkpoint'
                   ldasrcInp[keyn]= valn
                # specify LANDPERT restart file
                if (self.perturb == 1):
                    keyn = 'LANDPERT_INTERNAL_RESTART_FILE'
                    valn = '../input/restart/landpert'+tmpl_+'_internal_rst'
                    ldasrcInp[keyn]= valn
                    # for lat/lon and EASE tile space, specify LANDPERT checkpoint file here (via MAPL);
                    # for cube-sphere tile space, Landpert GC will set up LANDPERT checkpoint file
                    if ('-CF' not in self.rqdExeInp['GRIDNAME']):
                        keyn = 'LANDPERT_INTERNAL_CHECKPOINT_FILE'
                        valn = 'landpert'+tmpl_+'_internal_checkpoint'
                        ldasrcInp[keyn]= valn

                # add items for stretched grid
                if '-SG' in self.rqdExeInp['BCS_RESOLUTION']:
                   pos_ = self.rqdExeInp['BCS_RESOLUTION'].find('-SG')
                   SG   = self.rqdExeInp['BCS_RESOLUTION'][pos_+1:pos_+6]    # get ID of stretched grid (e.g., SG002)
                   ldasrcInp['STRETCH_FACTOR'] = STRETCH_GRID[SG][0]
                   ldasrcInp['TARGET_LAT']     = STRETCH_GRID[SG][1]
                   ldasrcInp['TARGET_LON']     = STRETCH_GRID[SG][2]

                # write LDAS.rc
                fout =open(self.rundir+'/'+shortfile,'w')
                # ldasrcInp['NUM_LDAS_ENSEMBLE']=ldasrcInp.pop('NUM_ENSEMBLE')
                for key,val in optinxny.items():
                    keyn=(key+":").ljust(36)
                    fout.write(keyn+str(val)+'\n')
                for key,val in ldasrcInp.items() :
                    keyn=(key+":").ljust(36)
                    fout.write(keyn+str(val)+'\n')
                fout.write("OUT_PATH:".ljust(36)+self.out_path+'\n')
                fout.write("EXP_ID:".ljust(36)+self.rqdExeInp['EXP_ID']+'\n')
                fout.write("TILING_FILE:".ljust(36)+"../input/tile.data\n")

                fout.close()

        fout=open(self.rundir+'/'+'cap_restart','w')
        #fout.write(self.rqdExeInp['BEG_DATE'])
        fout.write(self.begDates[0].strftime('%Y%m%d %H%M%S'))
        fout.close()
        status=True
        return status

    # -----------------------------------------------------------------------------------
           
    def createBatchRun(self):
        """
        """

        status = False

        os.chdir(self.rundir)
        fout =open(self.rundir+'/ldas_batchrun.j','w')
        fout.write("#!/bin/bash -f\n")
        jobid = None
        SBATCHQSUB = 'sbatch'
        expid = self.rqdExeInp['EXP_ID']
        if self.GEOS_SITE == 'NAS':
          SBATCHQSUB = 'qsub'
        fout.write("\nsed -i 's/if($capdate<$enddate) "+SBATCHQSUB+"/#if($capdate<$enddate) "+SBATCHQSUB+"/g' lenkf.j\n\n")
        nSegments = self.nSegments
        for iseg in range(nSegments):
            if iseg ==0 :
                fout.write("jobid%d=$(echo $(sbatch lenkf.j) | cut -d' ' -f 4)\n"%(iseg))
                fout.write("echo $jobid%d\n"%iseg )
            else :
                _start = self.begDates[iseg]
                myDateTime = '%04d%02d%02d_%02d%02dz' % \
                   (_start.year, _start.month, _start.day,_start.hour,_start.minute)
                _logfile = os.path.relpath(
                   '/'.join([
                        self.outdir,
                        self.rqdExeInp['EXP_DOMAIN'],
                        'rc_out',
                        'Y%04d' % _start.year,
                        'M%02d' % _start.month,
                        '.'.join([expid, 'ldas_log', myDateTime, 'txt']),
                        ]),
                    self.rundir)
                _errfile = os.path.relpath(
                   '/'.join([
                        self.outdir,
                        self.rqdExeInp['EXP_DOMAIN'],
                        'rc_out',
                        'Y%04d' % _start.year,
                        'M%02d' % _start.month,
                        '.'.join([expid, 'ldas_err', myDateTime, 'txt']),
                        ]),
                    self.rundir)

                #fout.write("jobid%d=$(echo $(sbatch --dependency=afterany:$jobid%d --output=%s --error=%s lenkf.j) | cut -d' ' -f 4)\n"%(iseg,iseg-1,_logfile, _errfile))
                fout.write("jobid%d=$(echo $(sbatch --dependency=afterok:$jobid%d lenkf.j) | cut -d' ' -f 4)\n"%(iseg,iseg-1))
                fout.write("echo $jobid%d\n"%iseg )
        fout.write("\nsed -i 's/#if($capdate<$enddate) "+SBATCHQSUB+"/if($capdate<$enddate) "+SBATCHQSUB+"/g' lenkf.j\n\n")
        fout.close()

        sp.call(['chmod', '755', self.rundir+'/ldas_batchrun.j'])
        status = True
        return status

    # -----------------------------------------------------------------------------------
           
    def createRunScripts(self):
        """
        """

        status = False

        os.chdir(self.rundir)

        my_qos='allnccs'
        if self.GEOS_SITE == 'NAS': my_qos = 'normal'
        if 'qos' in self.optRmInp :
            my_qos = self.optRmInp['qos']

        my_job=self.rqdExeInp['EXP_ID']
        if 'job_name' in self.optRmInp :
            my_job = self.optRmInp['job_name']

        start = self.begDates[0]
        expid = self.rqdExeInp['EXP_ID']
        myDateTime = '%04d%02d%02d_%02d%02dz' % \
                (start.year, start.month, start.day,start.hour,start.minute)
        my_logfile = os.path.relpath(
                '/'.join([
                        self.outdir,
                        self.rqdExeInp['EXP_DOMAIN'],
                        'rc_out',
                        'Y%04d' % start.year,
                        'M%02d' % start.month,
                        '.'.join([expid, 'ldas_log', myDateTime, 'txt']),
                        ]),
                self.rundir)
        my_errfile = os.path.relpath(
                '/'.join([
                        self.outdir,
                        self.rqdExeInp['EXP_DOMAIN'],
                        'rc_out',
                        'Y%04d' % start.year,
                        'M%02d' % start.month,
                        '.'.join([expid, 'ldas_err', myDateTime, 'txt']),
                        ]),
                self.rundir)

        constraint = '"[mil|cas]"'
        if self.GEOS_SITE == "NAS" :
            constraint = 'cas_ait'

        if 'constraint' in self.optRmInp:
            constraint = self.optRmInp['constraint']

        SBATCHQSUB = 'sbatch'
        if self.GEOS_SITE == 'NAS':
            SBATCHQSUB = 'qsub'

        DETECTED_MPI_STACK = "@MPI_STACK@"

        job_head = job_directive[self.GEOS_SITE]
        lenkf_str= (job_head+job_body).format(
                  SBATCHQSUB         = SBATCHQSUB,
                  MY_ACCOUNT         = self.rqdRmInp['account'],
                  MY_WALLTIME        = self.rqdRmInp['walltime'],
                  MY_NTASKS_MODEL    = str(self.rqdRmInp['ntasks_model']),
                  MY_CONSTRAINT      = constraint,
                  MY_OSERVER_NODES   = str(self.optRmInp['oserver_nodes']),
                  MY_WRITERS_NPES    = str(self.optRmInp['writers-per-node']),
                  MY_QOS             = my_qos,
                  MY_JOB             = my_job,
                  MY_EXPID           = self.rqdExeInp['EXP_ID'],
                  MY_EXPDOMAIN       = self.rqdExeInp['EXP_DOMAIN'],
                  MY_LOGFILE         = my_logfile,
                  MY_ERRFILE         = my_errfile,
                  MY_LANDMODEL       = self.catch,
                  MY_POSTPROC_HIST   = str(self.rqdExeInp['POSTPROC_HIST']),
                  MY_FIRST_ENS_ID    = str(self.first_ens_id),
                  MY_LADAS_COUPLING  = str(self.ladas_cpl),
                  MY_ENSEMBLE_FORCING= self.rqdExeInp.get('ENSEMBLE_FORCING', 'NO').upper(),
                  MY_ADAS_EXPDIR     = self.adas_expdir,
                  MY_EXPDIR          = self.expdir,
                  DETECTED_MPI_STACK = DETECTED_MPI_STACK,
                  )

        with open('lenkf.j','wt') as fout :
           fout.write(lenkf_str)
        sp.call(['chmod', '755', 'lenkf.j'])

        print ('\nExperiment directory: %s' % self.expdir)
        print ()
        status = True
        return status

# -----------------------------------------------------------------------------------
           
def _printdict(d):
    """
    Private method: print a 'flat' dictionary
    """

    for key, val in d.items():
        print (key.ljust(23), ':', val)

# -----------------------------------------------------------------------------------
    
def _produceExeInput(out_dict=None,ladas_cpl=0):
    """
    Private method: (1) Print sample exeinp file to stdout for offline GEOSldas setup:
                        _produceExeInput().
                    *or* 
                    (2) Create dictionary w/ default parameters from GEOSldas_LDAS.rc and 
                        GEOS_SurfaceGridComp.rc for coupled land-atm DAS setup:
                        _produceExeInput(out_dict, ladas_cpl=[1,2]).
    """
    
    if ladas_cpl > 0:
      assert out_dict is not None , " Need out_dict to hold the default parameters"

    # stand-alone (offline) LDAS:
    #   sample exeinp file includes placeholders for inputs that the user needs to provide
    if ladas_cpl == 0:
       print ('####################################################################################')
       print ('#                                                                                  #')
       print ('#                             REQUIRED INPUTS                                      #')
       print ('#                                                                                  #')
       print ('#             These inputs are needed to set up output dir structure.              #')
       print ('#                                                                                  #')
       print ('####################################################################################')
       print ()        
       print ('############################################################')
       print ('#                                                          #')
       print ('#                 EXPERIMENT INFO                          #')
       print ('#                                                          #')
       print ('# Format for start/end times is yyyymmdd hhmmss.           #')
       print ('#                                                          #')
       print ('############################################################')
       print ()
       print ('EXP_ID:')
       print ('EXP_DOMAIN:')
       print ('NUM_LDAS_ENSEMBLE:')
       print ('BEG_DATE:')
       print ('END_DATE:')
       print ()
       print ('############################################################')
       print ('#                                                          #')
       print ('#                 RESTART INFO                             #')
       print ('#                                                          #')
       print ('# (i) Select "RESTART" option:                             #')
       print ('#                                                          #')
       print ('# Use one of the following options if you *have* a         #')
       print ('#   GEOSldas restart file:                                 #')
       print ('#                                                          #')
       print ('# RESTART: 1                                               #')
       print ('#   YES, have restart file from GEOSldas                   #')
       print ('#        in SAME tile space (grid) with SAME boundary      #')
       print ('#        conditions and SAME snow model parameter (WEMIN). #')
       print ('#        The restart domain can be for the same or         #')
       print ('#        a larger one.                                     #')
       print ('#                                                          #')
       print ('# RESTART: 2                                               #')
       print ('#   YES, have restart file from GEOSldas but               #')
       print ('#        in a DIFFERENT tile space (grid) or with          #')
       print ('#        DIFFERENT boundary conditions or DIFFERENT snow   #')
       print ('#        model parameter (WEMIN).                          #')
       print ('#        Restart *must* be for the GLOBAL domain.          #')
       print ('#                                                          #')
       print ('# Use one of the following options if you DO NOT have a    #')
       print ('#   GEOSldas restart file                                  #')
       print ('#   (works for global domain ONLY!):                       #')
       print ('#                                                          #')
       print ('# RESTART: 0                                               #')
       print ('#   Cold start from some old restart for Jan 1, 0z.        #')
       print ('#                                                          #')
       print ('# RESTART: M                                               #')
       print ('#   Re-tile from archived MERRA-2 restart file.            #')
       print ('#                                                          #')
       print ('# -------------------------------------------------------- #')
       print ('# IMPORTANT:                                               #')
       print ('#   Except for RESTART=1, SPIN-UP is REQUIRED in almost    #')
       print ('#   all cases.                                             #')
       print ('# -------------------------------------------------------- #')
       print ('#                                                          #')
       print ('#                                                          #')
       print ('# (ii) Specify experiment ID/location of restart file:     #')
       print ('#                                                          #')
       print ('# For RESTART=1 or RESTART=2:                              #')
       print ('#   Specify RESTART_ID, RESTART_PATH, RESTART_DOMAIN with  #')
       print ('#   restarts stored as follows:                            #')
       print ('#     RESTART_PATH/RESTART_ID/output/RESTART_DOMAIN/rs/    #')
       print ('#                                                          #')
       print ('# For RESTART=0 or RESTART=M:                              #')
       print ('#   There is no need to specify RESTART_ID, RESTART_PATH,  #')
       print ('#   and RESTART_DOMAIN.                                    #')
       print ('#                                                          #')
       print ('############################################################')
       print ()
       print ('RESTART:')
       print ('#RESTART_ID:')
       print ('#RESTART_PATH:')
       print ('#RESTART_DOMAIN:')
       print ()
       print ('############################################################')
       print ('#                                                          #')
       print ('#         SURFACE METEOROLOGICAL FORCING                   #')
       print ('#                                                          #')
       print ('#  Surface meteorological forcing time step is in seconds. #')
       print ('#                                                          #')
       print ('#  NOTE:                                                   #')
       print ('#    When forcing is on cube-sphere (CS) grid, must use:   #')
       print ('#    - Model tile space (BCS) derived from same CS grid.   #')
       print ('#    - Nearest-neighbor interpolation (MET_HINTERP: 0).    #')    
       print ('#                                                          #')
       print ('#  For more information, see:                              #')
       print ('#    GEOSldas/doc/README.MetForcing_and_BCS.md             #')
       print ('#                                                          #')
       print ('############################################################')
       print ()
       print ('MET_TAG:')
       print ('MET_PATH:')
       print ('FORCE_DTSTEP:')
       print ()
       print ('############################################################')
       print ('#                                                          #')
       print ('#          LAND BOUNDARY CONDITIONS (BCS)                  #')
       print ('#                                                          #')
       print ('#  Path to and (atmospheric) resolution of BCS.            #')
       print ('#    Path includes BCS_VERSION.                            #')    
       print ('#                                                          #')
       print ('#  For more information, see:                              #')
       print ('#    GEOSldas/doc/README.MetForcing_and_BCS.md             #')
       print ('#    [..]/GEOSsurface_GridComp/Utils/Raster/make_bcs       #')
       print ('#                                                          #')
       print ('############################################################')
       print ()
       print ('BCS_PATH:')
       print ('BCS_RESOLUTION:')
       print ()
       print ('############################################################')

       # end if ladas_cpl==0
       
    # add defaults from rc template files
       
    current_directory = os.path.dirname(__file__)    # path where ldas_setup is

    # rc template files are in [current_directory]/../etc/

    # add defaults from GEOSldas_LDAS.rc
    
    _fn = current_directory+'/../etc/GEOSldas_LDAS.rc'      
    lines = []
    with open(_fn) as _f:
       i_ = 1
       for line in _f:
          if ( i_ < 5 or i_ >10):            # ignore lines 5-10 - may need to change if GEOSldas_LDAS.rc is edited
             if ladas_cpl == 0:
                sys.stdout.write(line)
                sys.stdout.flush()
             else:
                lines.append(line)
          i_ += 1
    print ()
    print ()

    # add land model parameter defaults from GEOS_SurfaceGridComp.rc

    _fn = current_directory+'/../etc/GEOS_SurfaceGridComp.rc'  

    if ladas_cpl == 0 :
        use_rc_defaults = 'GEOSldas=>'       # use defaults for LDAS
    else :
        use_rc_defaults = 'GEOSagcm=>'       # use defaults for AGCM
        
    with open(_fn) as _f :
        i_ = 1
        for line in _f:
            i_ +=1
            # skip over lines 5-21 (content does not apply after this processing) - may need to change if GEOS_SurfaceGridComp.rc is edited
            if ( 5<=i_ and i_<=21 ) :        
                continue
            if ladas_cpl == 0:
                # process lines that contain string "use_rc_defaults"
                if use_rc_defaults in line:                            
                    line0 = line.split(use_rc_defaults)[1]
                    sys.stdout.write(line0)
                    sys.stdout.flush()
                # echo blank lines and comment lines (except if they contain 'GEOSldas=>' or 'GEOSagcm=>')
                if (not 'GEOSldas=>' in line) and (not 'GEOSagcm=>' in line) :  
                    if (not line.strip()) or line.strip().startswith('#') :   
                        sys.stdout.write(line)
                        sys.stdout.flush()
            else:
                # process lines that contain string "use_rc_defaults"
                if use_rc_defaults in line:                            
                    line0 = line.split(use_rc_defaults)[1]
                    lines.append(line0)
    print ()
    print ()

    if ladas_cpl > 0:
       for line in lines:
          line = line.strip()
          if not line:
             continue
          if line.startswith('#'):
             continue
          # strip out inline comment (if present)
          position = line.find('#')
          if position>0:
             line = line[:position]
                    # extract key/value pair and add to dictionary
          key, val = line.split(":",1)
          out_dict[key.strip()] = val.strip()

# -----------------------------------------------------------------------------------
    
def _printRmInputKeys(rqdRmInpKeys, optRmInpKeys):
    """
    Private method: print sample resource manager input
    """

    print ('#')
    print ('# REQUIRED inputs')
    print ('#')
    print ('# NOTE:')
    print ('# - account          = computational project number')
    print ('#                      [At NCCS: Use command "getsponsor" to see available account number(s).]' )
    print ('# - walltime         = walltime requested; format is HH:MM:SS (hours/minutes/seconds)')
    print ('# - ntasks_model     = number of processors requested for the model (typically 126; output server is not included)')
    print ('#')
    for key in rqdRmInpKeys:
        print (key + ':')
    print ()
    print ('#')
    print ('# OPTIONAL inputs')
    print ('#')
    print ('# NOTE:')
    print ('# - job_name         = name of experiment; default is "exp_id"')
    print ('# - qos              = quality-of-service; do not specify by default; specify "debug" for faster but limited service.')
    print ('# - oserver_nodes    = number of nodes for oserver ( default is 0, for future use )')
    print ('# - writers-per-node = tasks per oserver_node for writing ( default is 5, for future use ),')
    print ('#      IMPORTANT REQUIREMENT: total #writers = writers-per-node * oserver_nodes >= 2')
    print ('#      Jobs will hang when oserver_nodes = writers-per-node = 1.')
    print ('# - constraint       = name of chip set(s) (NCCS default is "[mil|cas]", NAS default is "cas_ait").')
    print ('#')
    for key in optRmInpKeys:
        print ('#'+key + ':')

# -----------------------------------------------------------------------------------
                   
def parseCmdLine():
    """
    parse command line arguments and return a dict of options
    """

    #print 'in: parseCmdLine'
    p = argparse.ArgumentParser(
        description= \
            "Script to setup a GEOSldas experiment. "\
            "To set up an *offline* GEOSldas experiment, the script requires "\
            "two (2) input files, one for the experiment inputs ('exeinp') and the "\
            "other for the (SLURM) resource manager ('batinp'). To create sample input "\
            "files use './ldas_setup sample -h'. "\
            "When used to set up a coupled land-atm DAS experiment, the script "\
            "*ignores* the input files, but dummy strings must still be specified. ",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
        )
    p_sub = p.add_subparsers(help='sub-command help')

    
    # subparser: "sample" command
    p_sample = p_sub.add_parser(
        'sample',
        help='write sample input files',
        description='Print sample input files - either for the '\
            'experiment inputs or the resource manager (SLURM)',
        )
    group = p_sample.add_mutually_exclusive_group(required=True)
    group.add_argument(
        '--exeinp',
        help='print sample input file used to generate RC files for GEOSldas App.',
        action='store_true',
        )
    group.add_argument(
        '--batinp',
        help='print sample input file for the resource manager (SLURM)',
        action='store_true',
        )

    
    # subparser: "setup" command
    p_setup = p_sub.add_parser(
        'setup',
        help='setup LDAS experiment',
        description="The 'setup' sub-command is used to set up a GEOSldas " \
            "experiment. The positional argument 'exphome' is used to create " \
            "work_path (=exphome+/output) and run_path (=exphome+/run)."
        )
    p_setup.add_argument(
        '-v',
        '--verbose',
        help='verbose output',
        action='store_true',
        )
    p_setup.add_argument('exphome', help='experiment location')
    p_setup.add_argument(
        'exeinpfile',
        help='input file with arguments used to generate RC files for GEOSldas App',
        )
    p_setup.add_argument(
        'batinpfile',
        help='input file with arguments for SLURM',
        )

    # the following command line arguments, if present, take precedence over what
    #   is specified in the exeinp and batinp files
    ladas_cpl_help='land-atm DAS coupling mode: default=0 (not coupled); LDAS coupled with (1) ADAS deterministic member or (2) ADAS ensemble. Required when exeinp is dummy'
    p_setup.add_argument(
        '--ladas_cpl',
        help=ladas_cpl_help,
        type=str
        )
    p_setup.add_argument(
        '--account',
        help='overwrites computing/sponsor account from batinp file',
        type=str
        )
    p_setup.add_argument(
        '--nymdb',
        help='overwrites date in BEG_DATE from exeinp file: yyyymmdd. Required when exeinp is dummy',
        type=str
        )
    p_setup.add_argument(
        '--nhmsb',
        help='overwrites time in BEG_DATE from exeinp file: hhmmss. Required when exeinp is dummy',
        type=str
        ) 
    p_setup.add_argument(
        '--agcm_res',
        help='AGCM resolution associated with boundary conditions: nnnn (4 digits). Required when exeinp is dummy',
        type=str
        )
    p_setup.add_argument(
        '--bcs_version',
        help='boundary conditions version. Required when exeinp is dummy',
        type=str
        )
    p_setup.add_argument(
        '--rstloc',
        help='location of LDAS restarts (restart_path/restart_id). Required when exeinp is dummy',
        type=str
        )
    p_setup.add_argument(
        '--varwindow',
        help='ADAS analysis window (minutes). Required when exeinp is dummy',
        type=str
        )
    p_setup.add_argument(
        '--nens',
        help='number of ensemble members. Required when exeinp is dummy',          
        type=str
        )  

    # obsolete command line args
    p_setup.add_argument(
        '--runmodel',
        help='Obsolete.',
        action='store_true',
        )
    spltgrp = p_setup.add_mutually_exclusive_group()
    spltgrp.add_argument(
        '--daysperjob',
        type=int,
        metavar='N',
        help='Obsolete.  Use NUM_SGMT and JOB_SGMT in exeinp file.',
        )
    spltgrp.add_argument(
        '--monthsperjob',
        type=int,
        metavar='N',
        help='Obsolete.  Use NUM_SGMT and JOB_SGMT in exeinp file.',
        )

    return p.parse_args()



def hours_to_hhmmss(hours):
    
    # Convert hours to timedelta
    td = timedelta(hours=hours)
    
    # Extract hours, minutes, seconds
    total_seconds = int(td.total_seconds())
    hours, remainder = divmod(total_seconds, 3600)
    minutes, seconds = divmod(remainder, 60)
    
    # Format as HHMMSS
    return f"{hours:02d}{minutes:02d}{seconds:02d}"



if __name__=='__main__':

    resource.setrlimit(resource.RLIMIT_STACK, (resource.RLIM_INFINITY, resource.RLIM_INFINITY))
    #print "reading params...."
    args = vars(parseCmdLine()) # vars converts to dict
    ld = LDASsetup(args)

    print ("creating dir structure")
    status = ld.createDirStructure()
    assert(status)

    print ("creating links to restarts, BCs, met forcing, ...")
    status = ld.createLnRstBc()
    assert(status)

    print ("creating RC Files")
    status = ld.createRCFiles()
    assert(status)

    print ("creating gcm style batch Run scripts lenkf.j")
    status = ld.createRunScripts()
    assert(status)
    
    print ("creating batch Run scripts")
    status = ld.createBatchRun()
    assert (status)

# =================== EOF =======================================================================