From 4beb45f0c75c1e2e532987c109b6dc8db881e441 Mon Sep 17 00:00:00 2001
From: lkugler <lukas.kugler@gmail.com>
Date: Fri, 19 May 2023 18:07:32 +0200
Subject: [PATCH] evaluate only if necessary

---
 dartwrf/assim_synth_obs.py | 51 ++++++++++++++++++++++----------------
 1 file changed, 30 insertions(+), 21 deletions(-)

diff --git a/dartwrf/assim_synth_obs.py b/dartwrf/assim_synth_obs.py
index 738ac1b..d1cef07 100755
--- a/dartwrf/assim_synth_obs.py
+++ b/dartwrf/assim_synth_obs.py
@@ -259,26 +259,32 @@ def get_parametrized_error(obscfg, osf_prior):
         NotImplementedError('sat_channel not implemented', obscfg.get("sat_channel"))
 
 
-def set_obserr_assimilate_in_obsseqout(oso, osf_prior, outfile="./obs_seq.out"):
+def set_obserr_assimilate_in_obsseqout(oso, outfile="./obs_seq.out"):
     """"Overwrite existing variance values in obs_seq.out files
     
     Args:
         oso (ObsSeq): python representation of obs_seq.out file, will be modified and written to file
-        osf_prior (ObsSeq): python representation of obs_seq.final (output of filter in evaluate-mode without posterior)
-                            contains prior values; used for parameterized errors
 
     Returns:
         None    (writes to file)
+
+    Variables:
+        osf_prior (ObsSeq): python representation of obs_seq.final (output of filter in evaluate-mode without posterior)
+                        contains prior values; used for parameterized errors
     """
+
     for obscfg in exp.observations:
-        kind_str = obscfg['kind']
-        kind = osq.obs_kind_nrs[kind_str]
+        kind_str = obscfg['kind']  # e.g. 'RADIOSONDE_TEMPERATURE'
+        kind = osq.obs_kind_nrs[kind_str]  # e.g. 263
 
-        # modify each kind separately, one after each other
+        # modify observation error of each kind sequentially
         where_oso_iskind = oso.df.kind == kind
-        where_osf_iskind = osf_prior.df.kind == kind
-
+        
         if obscfg["error_assimilate"] == False:
+            osf_prior = obsseq.ObsSeq(cluster.dartrundir + "/obs_seq.final")  # this file will be generated by `evaluate()`
+            
+            where_osf_iskind = osf_prior.df.kind == kind
+
             assim_err = get_parametrized_error(obscfg, osf_prior.df[where_osf_iskind])
             oso.df.loc[where_oso_iskind, 'variance'] = assim_err**2
             #assert np.allclose(assim_err, oso.df['variance']**2)  # check
@@ -288,7 +294,9 @@ def set_obserr_assimilate_in_obsseqout(oso, osf_prior, outfile="./obs_seq.out"):
 
     oso.to_dart(outfile)
 
-def qc_obs(time, oso, osf_prior):
+def qc_obs(time, oso):
+    osf_prior = obsseq.ObsSeq(cluster.dartrundir + "/obs_seq.final")
+
     # obs should be superobbed already!
     for i, obscfg in enumerate(exp.observations): 
         if i > 0:
@@ -373,9 +381,6 @@ def evaluate(assim_time,
     copy(cluster.dartrundir + "/obs_seq.final", fout)
     print(fout, "saved.")
 
-    osf = obsseq.ObsSeq(cluster.dartrundir + "/obs_seq.final")
-    return osf
-
 
 
 def generate_obsseq_out(time):
@@ -522,9 +527,14 @@ def main(time, prior_init_time, prior_valid_time, prior_path_exp):
         None
     """
     nproc = cluster.max_nproc
+    do_QC = getattr(exp, "reject_smallFGD", False)  # True: triggers additional evaluations of prior & posterior
 
+    # for which observation type do we have a parametrized observation error?
+    error_is_parametrized = [obscfg["error_assimilate"] == False for obscfg in exp.observations]
+
+    # create directory to run DART in
     archive_time = cluster.archivedir + time.strftime("/%Y-%m-%d_%H:%M/")
-    os.makedirs(cluster.dartrundir, exist_ok=True)  # create directory to run DART in
+    os.makedirs(cluster.dartrundir, exist_ok=True)  
     os.chdir(cluster.dartrundir)
 
     # link DART binaries to run_DART
@@ -534,8 +544,6 @@ def main(time, prior_init_time, prior_valid_time, prior_path_exp):
     os.system("rm -f input.nml obs_seq.in obs_seq.out obs_seq.out-orig obs_seq.final")  
     set_DART_nml()
 
-    do_QC_here = getattr(exp, "reject_smallFGD", False)  # True: triggers additional evaluations of prior & posterior
-
     print("prepare nature")
     prepare_nature_dart(time)  # link WRF files to DART directory
 
@@ -545,16 +553,17 @@ def main(time, prior_init_time, prior_valid_time, prior_path_exp):
     print(" get observations with specified obs-error")
     oso = get_obsseq_out(time)
 
-    if do_QC_here:
+    # is any observation error parametrized?
+    if any(error_is_parametrized) or do_QC:
         print(" (optional) evaluate prior for all observations (incl rejected)")
-        osf_prior = evaluate(time, output_format="%Y-%m-%d_%H:%M_obs_seq.final-eval_prior_allobs")
+        evaluate(time, output_format="%Y-%m-%d_%H:%M_obs_seq.final-eval_prior_allobs")
 
     print(" assign observation-errors for assimilation ")
-    set_obserr_assimilate_in_obsseqout(oso, osf_prior, outfile=cluster.dartrundir + "/obs_seq.out")
+    set_obserr_assimilate_in_obsseqout(oso, outfile=cluster.dartrundir + "/obs_seq.out")
 
-    if do_QC_here:
+    if do_QC:
         print(" 2.3) reject observations? ")
-        qc_obs(time, oso, osf_prior)
+        qc_obs(time, oso)
 
     if exp.prior_inflation == 2:
         prepare_inflation_2(time, prior_init_time)
@@ -567,7 +576,7 @@ def main(time, prior_init_time, prior_valid_time, prior_path_exp):
     if exp.prior_inflation == 2:
         archive_inflation_2(time)
 
-    if do_QC_here:
+    if do_QC:
         print(" 4) evaluate posterior observations for all observations (incl rejected)")
         write_list_of_inputfiles_posterior(time)
         
-- 
GitLab