cycling directly from wrfout

config/cfg.py

@@ -9,7 +9,7 @@ class ExperimentConfiguration(object):
 
 
 exp = ExperimentConfiguration()
-exp.expname = "exp_v1.12_IR-WV-VIS"
+exp.expname = "exp_v1.16_P0-4_IR-VIS-WV"
 exp.model_dx = 2000
 exp.timestep = 10
 exp.n_ens = 40
@@ -44,7 +44,7 @@ psfc = dict(plotname='SYNOP Pressure', plotunits='[dBz]',
             cov_loc_radius_km=32, cov_loc_vert_km=5)
 
 
-exp.observations = [ir108, wv73, vis]
+exp.observations = [ir108, vis, wv73]  # ir108, wv73, vis]
 
 # directory paths depend on the name of the experiment
 cluster.expname = exp.expname

config/clusters.py

@@ -7,7 +7,7 @@ class ClusterConfig(object):
         pass
 
     def archivedir(self):
-        return '/gpfs/data/fs71386/lkugler/sim_archive/'+self.expname  # #return '/raid61/scratch/lkugler/VSC/'+self.expname
+        return self.archive_base+'/'+self.expname
 
     def wrf_rundir(self, iens):
         return self.userdir+'/run_WRF/'+self.expname+'/'+str(iens)
@@ -21,15 +21,17 @@ vsc.python = '/home/fs71386/lkugler/miniconda3/envs/DART/bin/python'
 vsc.ncks = '/home/fs71386/lkugler/miniconda3/envs/DART/bin/ncks'
 vsc.userdir = '/home/fs71386/lkugler'
 vsc.srcdir = '/home/fs71386/lkugler/compile/WRF/WRF-4.2.1/run'
+vsc.archive_base = '/gpfs/data/fs71386/lkugler/sim_archive/'
 vsc.dart_srcdir = '/home/fs71386/lkugler/DART/DART_WRF_RTTOV_early_access/models/wrf/work'
 vsc.dartrundir = '/home/fs71386/lkugler/run_DART'
 vsc.scriptsdir = '/home/fs71386/lkugler/DART-WRF/scripts'
 
-vsc.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.12_LMU_nature/2008-07-30_06:00/2/wrfout_d01_%Y-%m-%d_%H:%M:%S'
+vsc.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.13_P0_nature/2008-07-30_06:00/4/wrfout_d01_%Y-%m-%d_%H:%M:%S'
+#vsc.input_profile = '/home/fs71386/lkugler/wrf_sounding/data/wrf/ens/from_LMU/raso.nat.<iens>.wrfprof'
 vsc.input_profile = '/home/fs71386/lkugler/wrf_sounding/data/wrf/ens/from_LMU/raso.raso.<iens>.wrfprof'
 
-vsc.ideal = vsc.userdir+'/compile/bin/ideal-v4.2.1_v1.11.exe'
-vsc.wrfexe = vsc.userdir+'/compile/bin/wrf-v4.2.1_v1.12.exe'
+vsc.ideal = vsc.userdir+'/compile/bin/ideal-v4.2.2_v1.16.exe'
+vsc.wrfexe = vsc.userdir+'/compile/bin/wrf-v4.2.2_v1.16.exe'
 vsc.namelist = vsc.scriptsdir+'/../templates/namelist.input'
 vsc.run_WRF = '/gpfs/data/fs71386/lkugler/DART-WRF/scripts/osse/run_ens.vsc.sh'
 
@@ -44,6 +46,7 @@ jet.python = '/jetfs/home/lkugler/miniconda3/bin/python'
 jet.ncks = 'ncks'
 jet.userdir = '/jetfs/home/lkugler'
 jet.srcdir = '/jetfs/home/lkugler/compile/WRF/WRF-4.1.5/run'
+jet.archive_base = '/jetfs/home/lkugler/data_jetfs/sim_archive/'
 jet.dartrundir = '/jetfs/home/lkugler/DART-WRF/rundir'
 jet.scriptsdir = '/jetfs/home/lkugler/DART-WRF/scripts/osse'
 jet.nature_wrfout = '/raid61/scratch/lkugler/VSC/sim_archive/OSSE_v1.10_LMU+shear/2/single/wrfout_d01_%Y-%m-%d_%H:%M:%S'

scheduler.py

@@ -7,19 +7,20 @@ import os, sys, shutil
 import datetime as dt
 from slurmpy import Slurm
 
-# necessary to find modules in folder, since SLURM runs the script elsewhere
-sys.path.append(os.getcwd())
-
 from config.cfg import exp, cluster
 from scripts.utils import script_to_str, symlink
 
-# allow scripts to access the configuration
-symlink(cluster.scriptsdir+'/../config', cluster.scriptsdir+'/config')
+if __name__ == "__main__":
+    # necessary to find modules in folder, since SLURM runs the script elsewhere
+    sys.path.append(os.getcwd())
+
+    # allow scripts to access the configuration
+    symlink(cluster.scriptsdir+'/../config', cluster.scriptsdir+'/config')
 
-log_dir = cluster.archivedir()+'/logs/'
-slurm_scripts_dir = cluster.archivedir()+'/slurm-scripts/'
-print('logging to', log_dir)
-print('scripts, which are submitted to SLURM:', slurm_scripts_dir)
+    log_dir = cluster.archivedir()+'/logs/'
+    slurm_scripts_dir = cluster.archivedir()+'/slurm-scripts/'
+    print('logging to', log_dir)
+    print('scripts, which are submitted to SLURM:', slurm_scripts_dir)
 
 def my_Slurm(*args, cfg_update=dict(), **kwargs):
     """Shortcut to slurmpy's class; keep certain default kwargs
@@ -58,7 +59,8 @@ def backup_scripts():
 
     # reproducibility
     for f in ['scheduler.py', 'config/clusters.py', 'config/cfg.py']:
-        func(current+'/../'+f,  main_a+f, shutil.copy)
+        fname = os.path.basename(f)
+        func(current+'/../'+f,  main_a+'/'+fname, shutil.copy)
         
     for f in os.listdir(current):
         func(os.path.join(current, f), main_a+'/', shutil.copy)
@@ -103,35 +105,36 @@ def update_wrfinput_from_archive(valid_time, background_init_time, exppath, depe
                 +IC_path, depends_on=[depends_on])
     return id
 
-def run_ENS(begin, end, depends_on=None):
+def run_ENS(begin, end, depends_on=None, first_minute=True):
     """Run forecast for 1 minute, save output. 
     Then run whole timespan with 5 minutes interval.
     """
     id = depends_on
 
-    # first minute forecast (needed for validating an assimilation)
-    hist_interval = 1
-    radt = 1  # calc CFRAC also in first minute
-    begin_plus1 = begin+dt.timedelta(minutes=1)
-    s = my_Slurm("preWRF1", cfg_update=dict(time="2"))
-    id = s.run(' '.join([cluster.python,
-               cluster.scriptsdir+'/prepare_namelist.py',
-               begin.strftime('%Y-%m-%d_%H:%M'),
-               begin_plus1.strftime('%Y-%m-%d_%H:%M'),
-               str(hist_interval), str(radt),]), 
-             depends_on=[id])
-
-    s = my_Slurm("runWRF1", cfg_update={"nodes": "1", "array": "1-"+str(exp.n_nodes),
-                "time": "2", "mem-per-cpu": "2G"})
-    cmd = script_to_str(cluster.run_WRF).replace('<expname>', exp.expname)
-    id = s.run(cmd, depends_on=[id])
-
-    # apply forward operator (DART filter without assimilation)
-    s = my_Slurm("fwOP-1m", cfg_update=dict(time="10", ntasks=48))
-    id = s.run(cluster.python+' '+cluster.scriptsdir+'/apply_obs_op_dart.py '
-               + begin.strftime('%Y-%m-%d_%H:%M')+' '
-               + begin_plus1.strftime('%Y-%m-%d_%H:%M'),
-               depends_on=[id])
+    if first_minute:
+        # first minute forecast (needed for validating an assimilation)
+        hist_interval = 1
+        radt = 1  # calc CFRAC also in first minute
+        begin_plus1 = begin+dt.timedelta(minutes=1)
+        s = my_Slurm("preWRF1", cfg_update=dict(time="2"))
+        id = s.run(' '.join([cluster.python,
+                   cluster.scriptsdir+'/prepare_namelist.py',
+                   begin.strftime('%Y-%m-%d_%H:%M'),
+                   begin_plus1.strftime('%Y-%m-%d_%H:%M'),
+                   str(hist_interval), str(radt),]), 
+                 depends_on=[id])
+
+        s = my_Slurm("runWRF1", cfg_update={"nodes": "1", "array": "1-"+str(exp.n_nodes),
+                    "time": "2", "mem-per-cpu": "2G"})
+        cmd = script_to_str(cluster.run_WRF).replace('<expname>', exp.expname)
+        id = s.run(cmd, depends_on=[id])
+
+        # apply forward operator (DART filter without assimilation)
+        s = my_Slurm("fwOP-1m", cfg_update=dict(time="10", ntasks=48))
+        id = s.run(cluster.python+' '+cluster.scriptsdir+'/apply_obs_op_dart.py '
+                   + begin.strftime('%Y-%m-%d_%H:%M')+' '
+                   + begin_plus1.strftime('%Y-%m-%d_%H:%M'),
+                   depends_on=[id])
 
     # whole forecast timespan
     hist_interval = 5
@@ -210,9 +213,10 @@ def assimilate(assim_time, prior_init_time,
     return id
 
 
-def create_satimages(depends_on=None):
-    s = my_Slurm("pRTTOV", cfg_update={"ntasks": "48", "time": "40"})
-    s.run(cluster.python+' /home/fs71386/lkugler/RTTOV-WRF/loop.py '+exp.expname,
+def create_satimages(init_time, depends_on=None):
+    s = my_Slurm("pRTTOV", cfg_update={"ntasks": "48", "time": "30"})
+    s.run(cluster.python+' /home/fs71386/lkugler/RTTOV-WRF/run_init.py '+cluster.archivedir()
+          +init_time.strftime('/%Y-%m-%d_%H:%M/'),
           depends_on=[depends_on])
 
 def mailme(depends_on=None):
@@ -222,66 +226,66 @@ def mailme(depends_on=None):
 
 
 ################################
-print('starting osse')
+if __name__ == "__main__":
+    print('starting osse')
 
-timedelta_integrate = dt.timedelta(minutes=45)
-timedelta_btw_assim = dt.timedelta(minutes=30)
+    timedelta_integrate = dt.timedelta(minutes=75)
+    timedelta_btw_assim = dt.timedelta(minutes=60)
 
-backup_scripts()
-id = None
+    backup_scripts()
+    id = None
 
-start_from_existing_state = False
-is_new_run = not start_from_existing_state
+    start_from_existing_state = False
+    is_new_run = not start_from_existing_state
 
-if is_new_run:
-    id = prepare_wrfinput()  # create initial conditions
+    if is_new_run:
+        id = prepare_wrfinput()  # create initial conditions
 
-    # spin up the ensemble
-    background_init_time = dt.datetime(2008, 7, 30, 6, 0)
-    integration_end_time = dt.datetime(2008, 7, 30, 10, 0)
-    id = run_ENS(begin=background_init_time,
-                end=integration_end_time,
-                depends_on=id)
-    time = integration_end_time
-    first_guess = False
-    
-elif start_from_existing_state:
-    id = prepare_wrfinput()  # create initial conditions
-    
-    # get initial conditions from archive
-    init_time = dt.datetime(2008, 7, 30, 6)
-    time = dt.datetime(2008, 7, 30, 10)
-    exppath_arch = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.11_LMU_filter'
-    id = update_wrfinput_from_archive(time, init_time, exppath_arch, depends_on=id)
-
-# values for assimilation
-assim_time = time
-prior_init_time = init_time
-prior_path_exp = exppath_arch
-
-while time <= dt.datetime(2008, 7, 30, 17):
-
-    id = assimilate(assim_time,
-                    prior_init_time,
-                    prior_path_exp=prior_path_exp,
+        # spin up the ensemble
+        init_time = dt.datetime(2008, 7, 30, 6, 0)
+        integration_end_time = dt.datetime(2008, 7, 30, 9, 0)
+        id = run_ENS(begin=init_time,
+                    end=integration_end_time,
+                    first_minute=False,
                     depends_on=id)
-    prior_path_exp = False  # use own exp path
-
-    # integration
-    this_forecast_init = assim_time  # start integration from here
-    this_forecast_end = assim_time + timedelta_integrate
-
-    id = run_ENS(begin=this_forecast_init,
-                 end=this_forecast_end,
-                 depends_on=id)
+        time = integration_end_time
+        
+    elif start_from_existing_state:
+        
+        # get initial conditions from archive
+        init_time = dt.datetime(2008, 7, 30, 6)
+        integration_end_time = dt.datetime(2008, 7, 30, 10)
+        exppath_arch = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.13_P0_ps-t2'
+        id = update_wrfinput_from_archive(integration_end_time, init_time, exppath_arch, depends_on=id)
+
+    # values for assimilation
+    assim_time = integration_end_time
+    prior_init_time = init_time
+    prior_path_exp = False #exppath_arch
+
+    while time <= dt.datetime(2008, 7, 30, 18):
+
+        id = assimilate(assim_time,
+                        prior_init_time,
+                        prior_path_exp=prior_path_exp,
+                        depends_on=id)
+        prior_path_exp = False  # use own exp path
+
+        # integration
+        this_forecast_init = assim_time  # start integration from here
+        this_forecast_end = assim_time + timedelta_integrate
+
+        id = run_ENS(begin=this_forecast_init,
+                    end=this_forecast_end,
+                    depends_on=id)
+        
+        create_satimages(this_forecast_init, depends_on=id)
 
-    # increment time
-    time += timedelta_btw_assim
+        # increment time
+        time += timedelta_btw_assim
 
-    # values for next iteration
-    assim_time = time
-    prior_init_time = assim_time - timedelta_btw_assim
-    
-    create_satimages(depends_on=id)
+        # values for next iteration
+        assim_time = time
+        prior_init_time = time - timedelta_btw_assim
 
-mailme(id)
+    mailme(id)

scripts/update_wrfinput_from_filteroutput.py

@@ -9,22 +9,11 @@ time = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')
 background_init_time = dt.datetime.strptime(sys.argv[2], '%Y-%m-%d_%H:%M')
 exppath_firstguess = str(sys.argv[3])
 
-#if cluster.name != 'srvx8':
-#    copy = copy_scp_srvx8
-#    mkdir = mkdir_srvx8
-cycle_vars = ['U', 'V', 'P', 'PH', 'T', 'MU', 'QVAPOR', 'QCLOUD', 'QRAIN', 'QICE', 'QSNOW',
-              'QGRAUP', 'QNICE', 'QNRAIN', 'U10', 'V10', 'T2', 'Q2', 'PSFC', 'TSLB',
-              'SMOIS', 'TSK']
-
-update_vars = ['Times', 'U', 'V', 'T', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'PSFC', 'TSK', 'CLDFRA']
-# note: CLDFRA is just diagnostic, but maybe it helps that RTTOV can simulate clouds at initialization time
-
-# variables which are updated need not to be cycled
-for var in update_vars:
-    if var in cycle_vars:
-        cycle_vars.remove(var)
-
-cycles = ','.join(cycle_vars)
+"""
+# assumes T = THM (dry potential temperature as prognostic variable)
+
+"""
+update_vars = ['Times', 'U', 'V', 'T', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'TSK', 'CLDFRA']
 updates = ','.join(update_vars)
 
 print('move output to WRF dir as new initial conditions')
@@ -36,13 +25,14 @@ for iens in range(1, exp.n_ens+1):
     wrf_ic = cluster.wrf_rundir(iens) + '/wrfinput_d01'
 
     # cycles variables from wrfout (prior state)
-    print('cycling', cycles, 'into', wrf_ic, 'from', prior_wrf)
-    os.system(cluster.ncks+' -A -v '+cycles+' '+prior_wrf+' '+wrf_ic)
+    print('copy prior', prior_wrf, 'to wrfinput', wrf_ic)
+    # os.system(cluster.ncks+' -A -v '+cycles+' '+prior_wrf+' '+wrf_ic)
+    copy(prior_wrf, wrf_ic)
 
     print('updating', updates, 'in', wrf_ic, 'from', filter_out)
     os.system(cluster.ncks+' -A -v '+updates+' '+filter_out+' '+wrf_ic)
 
-    print('writing T into THM of wrfinput')
+    print('writing T into THM of wrfinput')  # assumes T = THM (dry potential temperature as prognostic variable)
     thm_in = nc.Dataset(filter_out, 'r').variables['T'][:]
     dsout = nc.Dataset(wrf_ic, 'r+')
     dsout.variables['THM'][:] = thm_in

scripts/update_wrfinput_from_wrfout.py

@@ -9,8 +9,9 @@ exppath = str(sys.argv[1])
 # time = dt.datetime.strptime(sys.argv[1], '%Y-%m-%d_%H:%M')
 # background_init_time = dt.datetime.strptime(sys.argv[2], '%Y-%m-%d_%H:%M')
 
-update_vars = ['Times', 'U', 'V', 'PH', 'T', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'PSFC', 'TSK']
-vars = ','.join(update_vars)
+#update_vars = ['Times', 'U', 'V', 'PH', 'T', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'PSFC', 'TSK']
+# vars = 'Times,U,V,T,THM,CLDFRA,P,PH,PHB,MU,QVAPOR,QCLOUD,QRAIN,QICE,QSNOW,QGRAUP,QNICE,QNRAIN,U10,V10,T2,Q2,PSFC,TSLB,SMOIS,TSK'  # ','.join(update_vars)
+# vars = 'Times,LU_INDEX,ZNU,ZNW,ZS,DZS,VAR_SSO,U,V,W,PH,PHB,T,THM,MU,MUB,P,PB,FNM,FNP,RDNW,RDN,DNW,DN,CFN,CFN1,THIS_IS_AN_IDEAL_RUN,P_HYD,Q2,T2,TH2,PSFC,U10,V10,RDX,RDY,RESM,ZETATOP,CF1,CF2,CF3,QVAPOR,QCLOUD,QRAIN,QICE,QSNOW,QGRAUP,QNICE,QNRAIN,SHDMAX,SHDMIN,SNOALB,TSLB,SMOIS,SH2O,SMCREL,SEAICE,IVGTYP,ISLTYP,VEGFRA,SNOW,SNOWH,CANWAT,LAI,QKE,VAR,F,E,SINALPHA,COSALPHA,HGT,TSK,P_TOP,T00,P00,TLP,TISO,TLP_STRAT,P_STRAT,CLDFRA,CLAT,ALBBCK,TMN,XLAND,SNOWC,SR,SAVE_TOPO_FROM_REAL,C1H,C2H,C1F,C2F,C3H,C4H,C3F,C4F,PCB,PC,LANDMASK,LAKEMASK,SST'  # these are all variables, that are in wrfout and wrfinput
 
 for iens in range(1, exp.n_ens+1):
     print('update state in wrfinput  wrfout file to DART background file')
@@ -21,4 +22,5 @@ for iens in range(1, exp.n_ens+1):
     assert os.path.isfile(wrfout), wrfout
 
     # overwrite variables in wrfinput file
-    os.system(cluster.ncks+' -A -v '+vars+' '+wrfout+' '+wrfin)
+    # os.system(cluster.ncks+' -A -v '+vars+' '+wrfout+' '+wrfin)
+    copy(wrfout, wrfin)