works-ish with mask value changed

NOAA-EMC · Jul 11, 2023 · 2a0dfdc · 2a0dfdc
1 parent d765863
commit 2a0dfdc
Show file tree

Hide file tree

Showing 11 changed files with 164 additions and 87 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -94,7 +94,7 @@ if(BUILD_GDASBUNDLE)
   endif()
   ecbuild_bundle( PROJECT gsw GIT "https://github.com/jcsda-internal/GSW-Fortran.git" BRANCH develop )
   ecbuild_bundle( PROJECT mom6 GIT "https://github.com/jcsda-internal/MOM6.git" BRANCH main-ecbuild RECURSIVE )
-  ecbuild_bundle( PROJECT soca GIT "https://github.com/jcsda-internal/soca.git" BRANCH develop )
+  ecbuild_bundle( PROJECT soca GIT "https://github.com/jcsda-internal/soca.git" BRANCH develop )# feature/remap_uv )
 
 # Build IODA converters
   option(BUILD_IODA_CONVERTERS "Build IODA Converters" ON)

diff --git a/parm/soca/berror/saber_blocks.yaml b/parm/soca/berror/saber_blocks.yaml
@@ -45,27 +45,23 @@ components:
         standard deviation: true
 
       - linear variable change name: BalanceSOCA
-#        kst:
-#          dsdtmax: 0.1
-#          dsdzmin: 3.0e-6
-#          dtdzmin: 1.0e-6
-#          nlayers: 20
-
         ksshts:
           nlayers: 10
   weight:
-    value: 0.5
+    value: 1.0
 - covariance:
     covariance model: ensemble
     members from template:
       template:
         read_from_file: 1
         date: '{{ATM_WINDOW_MIDDLE}}'
         basename: ./static_ens/
-        ocn_filename: ocn.%mem%.nc
-        ice_filename: ice.%mem%.nc
-        remap_filename: ./INPUT/MOM.res.nc
-        state variables: [tocn, socn, ssh, uocn, vocn, hocn, cicen, hicen, hsnon, mld, layer_depth]
+        ocn_filename: 'ocn.filtered.%mem%.incr.{{ATM_WINDOW_BEGIN}}.nc'
+        ice_filename: 'ice.filtered.%mem%.incr.{{ATM_WINDOW_BEGIN}}.nc'
+#        ocn_filename: 'ocn.%mem%.nc'
+#        ice_filename: 'ice.%mem%.nc'
+#        remap_filename: ./INPUT/MOM.res.nc
+        state variables: [tocn, socn, ssh, uocn, vocn, cicen, hicen, hsnon]
       pattern: '%mem%'
       nmembers: ${CLIM_ENS_SIZE}
     localization:
@@ -80,15 +76,16 @@ components:
             read local nicas: true
           model:
             do not cross mask boundaries: false
-#            groups:
-#            - group name: group 1
-#              variables: [tocn, socn, uocn, vocn, ssh]
-#            - group name: group 2
-#              variables: [cicen, hicen, hsnon]
-
+##            groups:
+##            - group name: group 1
+##              variables: [tocn, socn, uocn, vocn, ssh]
+##            - group name: group 2
+##              variables: [cicen, hicen, hsnon]
+##
   weight:
-    read_from_file: 3
-    basename: ./
-    ocn_filename: 'ocn.ens_weights.incr.{{ATM_WINDOW_MIDDLE}}.nc'
-    ice_filename: 'ice.ens_weights.incr.{{ATM_WINDOW_MIDDLE}}.nc'
-    date: '{{ATM_WINDOW_MIDDLE}}'
+    value: 1.0
+#    read_from_file: 3
+#    basename: ./
+#    ocn_filename: 'ocn.ens_weights.incr.{{ATM_WINDOW_MIDDLE}}.nc'
+#    ice_filename: 'ice.ens_weights.incr.{{ATM_WINDOW_MIDDLE}}.nc'
+#    date: '{{ATM_WINDOW_MIDDLE}}'
diff --git a/parm/soca/berror/soca_apply_bm1.yaml b/parm/soca/berror/soca_apply_bm1.yaml
@@ -11,17 +11,17 @@ output geometry:
 linear variable change:
   input variables: [tocn, socn, ssh, uocn, vocn, hocn]
   output variables: [tocn, socn, ssh, uocn, vocn, hocn]
-  do inverse: true
+  #do inverse: true
   linear variable changes:
   - linear variable change name: BalanceSOCA
-    kst:
-      dsdtmax: 0.1
-      dsdzmin: 3.0e-6
-      dtdzmin: 1.0e-6
-      nlayers: 20
+#    kst:
+#      dsdtmax: 0.1
+#      dsdzmin: 3.0e-6
+#      dtdzmin: 1.0e-6
+#      nlayers: 20
 
-    ksshts:
-      nlayers: 10
+#    ksshts:
+#      nlayers: 10
 
 increments:
 - date: '{{ATM_WINDOW_BEGIN}}'

diff --git a/parm/soca/berror/soca_apply_steric.yaml b/parm/soca/berror/soca_apply_steric.yaml
@@ -0,0 +1,43 @@
+input geometry:
+  geom_grid_file: soca_gridspec.nc
+  mom6_input_nml: mom_input.nml
+  fields metadata: fields_metadata.yaml
+
+output geometry:
+  geom_grid_file: soca_gridspec.nc
+  mom6_input_nml: mom_input.nml
+  fields metadata: fields_metadata.yaml
+
+linear variable change:
+  input variables: [tocn, socn, ssh, uocn, vocn, hocn, cicen, hicen, hsnon]
+  output variables: [tocn, socn, ssh, uocn, vocn, hocn, cicen, hicen, hsnon]
+  do inverse: false
+  linear variable changes:
+  - linear variable change name: BkgErrFILT
+    ocean_depth_min: 500  # zero where ocean is shallower than 500m
+    rescale_bkgerr: 0.3   # rescale perturbation
+    efold_z: 1500.0       # Apply exponential decay
+  - linear variable change name: BalanceSOCA  # linear steric height from (T,S) perturbation
+
+increments:
+- date: '{{ATM_WINDOW_BEGIN}}'
+  input variables: [tocn, socn, ssh, uocn, vocn, hocn, cicen, hicen, hsnon]
+  input:
+     read_from_file: 1
+     basename: ./static_ens/
+     ocn_filename: 'ocn.bal.ens.MEMNUM.{{ATM_WINDOW_BEGIN}}.PT0S.nc'
+     ice_filename: 'ice.bal.ens.MEMNUM.{{ATM_WINDOW_BEGIN}}.PT0S.nc'
+     date: '{{ATM_WINDOW_BEGIN}}'
+     state variables: [ssh, tocn, socn, uocn, vocn, cicen, hicen, hsnon]
+  trajectory:
+    read_from_file: 1
+    basename: ./INPUT/
+    ocn_filename: MOM.res.nc
+    ice_filename: cice.res.nc
+    date: '{{ATM_WINDOW_BEGIN}}'
+    state variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh, hocn, mld, layer_depth]
+  output:
+     datadir: ./static_ens
+     exp: filtered.MEMNUM
+     type: incr
+     date: '{{ATM_WINDOW_BEGIN}}'
diff --git a/parm/soca/berror/soca_clim_ens_perts.yaml b/parm/soca/berror/soca_clim_ens_perts.yaml
@@ -3,7 +3,7 @@ geometry:
   mom6_input_nml: mom_input.nml
   fields metadata: fields_metadata.yaml
 
-recenter variables:  [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh]
+recenter variables:  [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh, hocn, mld, layer_depth]
 
 zero center: True
 
@@ -13,7 +13,7 @@ center:
   ocn_filename: MOM.res.nc
   ice_filename: cice.res.nc
   date: '{{ATM_WINDOW_BEGIN}}'
-  state variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh, hocn]
+  state variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh, hocn, mld, layer_depth]
 
 ensemble:
   members from template:
@@ -24,7 +24,7 @@ ensemble:
       ocn_filename: ocn.%mem%.nc
       ice_filename: ice.%mem%.nc
       remap_filename: ./INPUT/MOM.res.nc
-      state variables: [tocn, socn, ssh, uocn, vocn, hocn, cicen, hicen, hsnon]
+      state variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh, hocn, mld, layer_depth]
     pattern: '%mem%'
     nmembers: ${CLIM_ENS_SIZE}
 

diff --git a/parm/soca/berror/soca_ensweights.yaml b/parm/soca/berror/soca_ensweights.yaml
@@ -16,8 +16,9 @@ background:
   read_from_file: 1
 
 weights:
-  ice: 0.000001
-  ocean: 0.5
+  # Need to provide weights^2 when reading from file
+  ice: 0.0025   #  5% of original variance
+  ocean: 0.01   # 10%      "       "
 
 output:
   datadir: ./

diff --git a/parm/soca/variational/3dvarfgat.yaml b/parm/soca/variational/3dvarfgat.yaml
@@ -3,14 +3,6 @@ cost function:
   window begin: '{{ATM_WINDOW_BEGIN}}'
   window length: $(ATM_WINDOW_LENGTH)
   analysis variables: &soca_ana_vars  [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh]
-#  variable change:
-#    variable change name: Ana2Model
-#    rotate:
-#      u: []
-#      v: []
-#    log:
-#      var: []
-#    output variables: *soca_ana_vars
 
   geometry:
     mom6_input_nml: mom_input.nml
@@ -53,7 +45,7 @@ output:
   datadir: Data
   exp: 3dvarfgat_pseudo
   type: an
-  frequency: PT1H # can only writes out the analysis at the middle of the window
+  frequency: PT3H # can only writes out the analysis at the middle of the window
                   # when using 3D-FGAT as cost function
 
 final:

diff --git a/scripts/exgdas_global_marine_analysis_bmat.sh b/scripts/exgdas_global_marine_analysis_bmat.sh
@@ -112,18 +112,17 @@ if [ $err -gt 0  ]; then
 fi
 
 ################################################################################
-# Compute the ens std. dev of B^-1 X
+# Compute the ens std. dev, ignore ssh variance
 # TODO (G): Implement what's below into one single oops application
 # 0 - Diagnostic: Compute moments of original ensemble
 # 1 - Ensemble perturbations:
-#       o Vertically Interpolate into the deterministic layers
+#       o Vertically Interpolate to the deterministic layers
 #       o Recenter around 0 to create an ensemble of perturbations
-# 2 - Apply the inverse of the balance to the ensemble members
-
+# 2 - Filter members and apply the linear steric height balance to each members
 # 3 - Copy h from deterministic to unbalanced perturbations
-# 3 - Compute moments of converted ensemble perturbations
+# 4 - Compute moments of converted ensemble perturbations
 
-# Compute ensemble moments (not really needed but a useful diagnostic)
+# Compute ensemble moments
 clean_yaml soca_clim_ens_moments.yaml
 $APRUN_OCNANAL $JEDI_BIN/soca_ensmeanandvariance.x soca_clim_ens_moments.yaml
 export err=$?; err_chk
@@ -132,38 +131,82 @@ if [ $err -gt 0  ]; then
 fi
 
 # Compute perturbations
-#clean_yaml soca_clim_ens_perts.yaml
-#$APRUN_OCNANAL $JEDI_BIN/soca_ensrecenter.x soca_clim_ens_perts.yaml
-#export err=$?; err_chk
-#if [ $err -gt 0  ]; then
-#    exit $err
-#fi
-
-# Compute (B^-1 X)
-#clean_yaml soca_apply_bm1.yaml
-#$APRUN_OCNANAL $JEDI_BIN/soca_convertincrement.x soca_apply_bm1.yaml
-#export err=$?; err_chk
-#if [ $err -gt 0  ]; then
-#    exit $err
-#fi
+clean_yaml soca_clim_ens_perts.yaml
+$APRUN_OCNANAL $JEDI_BIN/soca_ensrecenter.x soca_clim_ens_perts.yaml
+export err=$?; err_chk
+if [ $err -gt 0  ]; then
+    exit $err
+fi
 
-# Recenter around deterministic
-#clean_yaml soca_recenter_det.yaml
-#$APRUN_OCNANAL $JEDI_BIN/soca_ensrecenter.x soca_recenter_det.yaml
-#export err=$?; err_chk
-#if [ $err -gt 0  ]; then
-#    exit $err
-#fi
+# Filter perturbations
+clean_yaml soca_apply_steric.yaml
+$APRUN_OCNANAL $JEDI_BIN/soca_convertincrement.x soca_apply_steric.yaml
+export err=$?; err_chk
+if [ $err -gt 0  ]; then
+    exit $err
+fi
 
 # Append hocn from deterministic to members
-#unbal_pert_files=$(ls ./static_ens/ocn.unbal.*.nc)
-#ncks -O -v h ./INPUT/MOM.res.nc ./static_ens/hocn.nc
-#ncrename -d zl,zaxis_1 -d yh,yaxis_1 -d xh,xaxis_1 -d time,Time ./static_ens/hocn.nc
-#for pert_file in ${unbal_pert_files}; do
-#    echo "ooooooooooo ${pert_file}"
-#    ncks -O -x -v h ${pert_file} ${pert_file}
-#    ncks -A -v h ./static_ens/hocn.nc ${pert_file}
-#done
+filt_pert_files=$(ls ./static_ens/ocn.filtered.*.nc)
+ncks -O -v h ./INPUT/MOM.res.nc ./static_ens/hocn.nc
+ncrename -d zl,zaxis_1 -d yh,yaxis_1 -d xh,xaxis_1 -d time,Time ./static_ens/hocn.nc
+for pert_file in ${filt_pert_files}; do
+    echo "ooooooooooo ${pert_file}"
+    ncks -O -x -v h ${pert_file} ${pert_file}
+    ncks -A -v h ./static_ens/hocn.nc ${pert_file}
+done
+
+## Replace undef with zero's
+## TODO: this step should not be necessary ...
+flist=$(ls ./static_ens/ocn.filtered.*.nc)
+for f in $flist; do
+  fout=$(basename $f)
+  ncatted -O -a _FillValue,,d,, -a missing_value,,d,, $f $fout
+  ncap2 -O -s "where(ave_ssh!=ave_ssh) ave_ssh = 0" \
+           -s "where(abs(Temp) > 999) Temp = 0" \
+           -s "where(abs(Salt) > 999) Salt = 0" \
+           -s "where(abs(u) > 999) u = 0" \
+           -s "where(abs(v) > 999) v = 0" \
+           -s "where(abs(h) > 999) h = 0" \
+           $fout $f
+done
+
+flist=$(ls ./static_ens/ice.filtered.*.nc)
+for f in $flist; do
+  echo $f
+  fout=$(basename $f)
+  ncatted -O -a _FillValue,,d,, -a missing_value,,d,, $f $fout
+  ncap2 -O -s "where(abs(aicen) > 999) aicen = 0" \
+           -s "where(abs(hicen) > 999) hicen = 0" \
+           -s "where(abs(hsnon) > 999) hsnon = 0" \
+           $fout $f
+done
+
+# Fix background
+flist=$(ls ./bkg/*.ocnf00*.nc)
+for f in $flist; do
+  fout=$(basename $f)
+  ncatted -O -a _FillValue,,d,, -a missing_value,,d,, $f $fout
+  ncap2 -O -s "where(abs(ave_ssh) > 999) ave_ssh = 0" \
+           -s "where(abs(Temp) > 999) Temp = 0" \
+           -s "where(abs(Salt) > 999) Salt = 0" \
+           -s "where(abs(u) > 999) u = 0" \
+           -s "where(abs(v) > 999) v = 0" \
+           -s "where(abs(h) > 999) h = 0" \
+           -s "where(abs(MLD) > 999) MLD = 0" \
+           $fout $f
+done
+
+flist=$(ls ls ./static_ens/ice.filtered.*.nc)
+flist=$(ls ./bkg/*.agg_icef00*.nc)
+for f in $flist; do
+  fout=$(basename $f)
+  ncatted -O -a _FillValue,,d,, -a missing_value,,d,, $f $fout
+  ncap2 -O -s "where(abs(aicen) > 999) aicen = 0" \
+           -s "where(abs(hicen) > 999) hicen = 0" \
+           -s "where(abs(hsnon) > 999) hsnon = 0" \
+           $fout $f
+done
 
 # Compute std. dev. (B^-1 X)
 #clean_yaml soca_clim_ens_stddev_b.yaml

diff --git a/scripts/exgdas_global_marine_analysis_bmat_vrfy.sh b/scripts/exgdas_global_marine_analysis_bmat_vrfy.sh
@@ -41,7 +41,7 @@ export DIRAC="${HOMEgfs}/sorc/gdas.cd/ush/ufsda/dirac_yaml.py"
 
 arr=("tocn" "socn" "ssh" "cicen" "hicen")
 level=1
-ndiracs=100
+ndiracs=5
 
 for i in "${!arr[@]}"
 do
@@ -54,7 +54,7 @@ exp: dirac_${var}_${level}
 type: an
 EOL
 
-    ${DIRAC} --varyaml 'var.yaml' \
+    ${DIRAC} --varyaml 'var_original.yaml' \
              --fields 'soca_gridspec.nc' \
              --dim1 'xaxis_1' \
              --dim2 'yaxis_1' \

diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py
@@ -372,9 +372,9 @@ def find_clim_ens(input_date):
 config = Template.substitute_structure(config, TemplateConstants.DOLLAR_PARENTHESES, envconfig.get)
 config.save(berr_yaml)
 
-logging.info(f"---------------- generate soca_apply_bm1.yaml")
-berr_yaml = os.path.join(anl_dir, 'soca_apply_bm1.yaml')
-berr_yaml_template = os.path.join(berror_yaml_dir, 'soca_apply_bm1.yaml')
+logging.info(f"---------------- generate soca_apply_steric.yaml")
+berr_yaml = os.path.join(anl_dir, 'soca_apply_steric.yaml')
+berr_yaml_template = os.path.join(berror_yaml_dir, 'soca_apply_steric.yaml')
 config = YAMLFile(path=berr_yaml_template)
 config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get)
 config = Template.substitute_structure(config, TemplateConstants.DOLLAR_PARENTHESES, envconfig.get)

diff --git a/ush/ufsda/dirac_yaml.py b/ush/ufsda/dirac_yaml.py
@@ -32,15 +32,16 @@ def var2dirac(args):
     diracconfig['background error'] = varconfig['cost function']['background error']
 
     # Overwrite variables
-    state_vars = ['cicen', 'hicen', 'hsnon', 'socn', 'tocn', 'ssh', 'hocn', 'mld', 'layer_depth']
+    state_vars = ['cicen', 'hicen', 'hsnon', 'socn', 'tocn', 'uocn', 'vocn', 'ssh', 'hocn', 'mld', 'layer_depth']
     diracconfig['initial condition']['state variables'] = state_vars
 
     # str to integer
     nmem=int(diracconfig['background error']['components'][1]['covariance']['members from template']['nmembers'] )
     diracconfig['background error']['components'][1]['covariance']['members from template']['nmembers'] = nmem
 
     # Set date to the middle of the window
-    window_middle = diracconfig['background error']['components'][1]['weight']['date']
+    #window_middle = varconfig['background error']['components'][1]['weight']['date']
+    window_middle = varconfig['cost function']['model']['states'][0]['date']
     diracconfig['initial condition']['date'] = window_middle
 
     # Remove static B