In this section, we'll get some experience in running single-time EDA and LETKF analysis test.
Here, we will run the single-time EDA analysis, which has a total of 10 3DEnVar members.
First, We need to create the working directory for EDA test and move to that directory:
$ mkdir -p /glade/derecho/scratch/$USER/mpas_jedi_tutorial/eda
$ cd /glade/derecho/scratch/$USER/mpas_jedi_tutorial/eda
Set environment variable for the mpas-bundle directory,
$ export bundle_dir=/glade/derecho/scratch/liuz/mpas_bundle_v3_public_gnuSP
Then, we need to prepare necessary files for EDA test.
Link MPAS physics-related data and table, MPAS graph, MPAS streams, and MPAS namelist files:
$ ln -sf ../MPAS_namelist_stream_physics_files/*BL .
$ ln -sf ../MPAS_namelist_stream_physics_files/*DATA .
$ ln -sf ../MPAS_namelist_stream_physics_files/x1.10242.graph.info.part.64 .
$ ln -sf ../MPAS_namelist_stream_physics_files/streams.atmosphere_240km .
$ ln -sf ../MPAS_namelist_stream_physics_files/stream_list.atmosphere.* .
$ ln -sf ../MPAS_namelist_stream_physics_files/namelist.atmosphere_240km .
Link MPAS 2-stream files:
$ ln -sf ../MPAS_namelist_stream_physics_files/x1.10242.invariant.nc .
$ ln -sf ../background/2018041418/mpasout.2018-04-15_00.00.00.nc ./templateFields.10242.nc
Link static yaml files for MPAS-JEDI:
$ ln -sf ${bundle_dir}/code/mpas-jedi/test/testinput/namelists/geovars.yaml ./
$ ln -sf ${bundle_dir}/code/mpas-jedi/test/testinput/namelists/keptvars.yaml ./
$ ln -sf ${bundle_dir}/code/mpas-jedi/test/testinput/obsop_name_map.yaml ./
Prepare the ensemble background as well as ensemble analysis directories and files for EDA. Since we have 20 members, it will be better that we deal with them by creating a simple cshell script prep_an_bg_directories.csh:
#!/bin/csh
foreach imem (01 02 03 04 05 06 07 08 09 10)
mkdir -p bg/mem0${imem} an/mem0${imem}
set tpath=/glade/derecho/scratch/$USER/mpas_jedi_tutorial
ln -sf ${tpath}/ensemble/2018041418/${imem}/mpasout.2018-04-15_00.00.00.nc bg/mem0${imem}/bg.2018-04-15_00.00.00.nc
cp ../ensemble/2018041418/${imem}/mpasout.2018-04-15_00.00.00.nc an/mem0${imem}/an.2018-04-15_00.00.00.nc
end
After running prep_an_bg_directories.csh, the an and bg directories and files are prepared.
$ chmod 755 ./prep_an_bg_directories.csh
$ ./prep_an_bg_directories.csh
Next step, we create a directory dbIn for observation input, and a directory dbOut to store the observation feedback files of all EDA members.
$ mkdir dbIn
$ cd dbIn
$ ln -sf ../../obs_ioda_pregenerated/2018041500/aircraft_obs_2018041500.h5 ./
$ ln -sf ../../obs_ioda_pregenerated/2018041500/gnssro_obs_2018041500.h5 ./
$ ln -sf ../../obs_ioda_pregenerated/2018041500/sfc_obs_2018041500.h5 ./
$ ln -sf ../../obs_ioda_pregenerated/2018041500/satwind_obs_2018041500.h5 ./
$ ln -sf ../../obs_ioda_pregenerated/2018041500/sondes_obs_2018041500.h5 ./
$ cd ..
For dbOut directory, we need to create subdirectories for all EDA members with a simple cshell script, prep_dbout_directory.csh:
#!/bin/csh
foreach imem (01 02 03 04 05 06 07 08 09 10)
mkdir -p dbOut/mem0${imem}
end
By running the prep_dbout_directory.csh with the command:
$ chmod 755 ./prep_dbout_directory.csh
$ ./prep_dbout_directory.csh
the dbOut directory and its subdirectories are created.
Also, link the localization matrix and static BEC files.
$ ln -sf ../localization_pregenerated ./bumploc
$ ln -sf ../B_Matrix ./
Last, copy related yaml files and scripts for EDA running:
$ ln -sf ../MPAS_JEDI_yamls_scripts/eda_3denvar.yaml ./
$ ln -sf ../MPAS_JEDI_yamls_scripts/eda_rtpp.yaml ./
$ ln -sf ../MPAS_JEDI_yamls_scripts/run_eda_members.csh ./
$ ln -sf ../MPAS_JEDI_yamls_scripts/run_eda_rtpp.csh ./
The eda_3denvar.yaml is a template yaml file for all 3DEnVar EDA members. the eda_3denvar.yaml file is very similar to 3denvar.yaml. You may do "diff 3denvar.yaml eda_3denvar.yaml" to see their difference.
18,19c18,19
< member: 1
< number of members: 1
---
> member: {{imem}}
> number of members: 10
38d37
<
40c39
< filename: ./an.$Y-$M-$D_$h.$m.$s.nc
---
> filename: ../../an/mem{{IMEM}}/an.$Y-$M-$D_$h.$m.$s.nc
50,51c49,50
< - <<: *iterationConfig
< ninner: 50
---
> # - <<: *iterationConfig
> # ninner: 50
69c68
< filename: ./bg.2018-04-15_00.00.00.nc
---
> filename: ../../bg/mem{{IMEM}}/bg.2018-04-15_00.00.00.nc
80c79
< files prefix: ./BUMP_files/bumploc_1200km6km
---
> files prefix: ../../bumploc/bumploc_1200km6km
89c88
< filename: ../ensemble/2018041418/%iMember%/mpasout.2018-04-15_00.00.00.nc
---
> filename: ../../bg/mem%iMember%/bg.2018-04-15_00.00.00.nc
92,93c91,93
< zero padding: 2
< nmembers: 10
---
> zero padding: 3
> nmembers: 9
> except: [{{imem}}]
95c95
< obs perturbations: false
---
> obs perturbations: true
103c103
< obsfile: ./aircraft_obs_2018041500.h5
---
> obsfile: ../../dbIn/aircraft_obs_2018041500.h5
107c107
< obsfile: ./obsout_da_aircraft.h5
---
> obsfile: ../../dbOut/mem{{IMEM}}/obsout_da_aircraft.h5
131c131
< obsfile: ./gnssro_obs_2018041500.h5
---
> obsfile: ../../dbIn/gnssro_obs_2018041500.h5
135c135
< obsfile: ./obsout_da_gnssroref.h5
---
> obsfile: ../../dbOut/mem{{IMEM}}/obsout_da_gnssroref.h5
175c175
< obsfile: ./satwind_obs_2018041500.h5
---
> obsfile: ../../dbIn/satwind_obs_2018041500.h5
179c179
< obsfile: ./obsout_da_satwind.h5
---
> obsfile: ../../dbOut/mem{{IMEM}}/obsout_da_satwind.h5
228c228
< obsfile: ./sfc_obs_2018041500.h5
---
> obsfile: ../../dbIn/sfc_obs_2018041500.h5
232c232
< obsfile: ./obsout_da_sfc.h5
---
> obsfile: ../../dbOut/mem{{IMEM}}/obsout_da_sfc.h5
262c262
< obsfile: ./sondes_obs_2018041500.h5
---
> obsfile: ../../dbIn/sondes_obs_2018041500.h5
266c266
< obsfile: ./obsout_da_sondes.h5
---
> obsfile: ../../dbOut/mem{{IMEM}}/obsout_da_sondes.h5
Even though we can run each EDA member simulataneously, for this tutorial we run all 10 EDA members serially in one job, which will take about 10 minutes (one outer loop setting). In run_eda_members.csh, you may modify the line foreach iens to run less EDA members to save time.
Submit the job script run_eda_members.csh with:
$ qsub run_eda_members.csh
When all EDA members are done, we will see that all an/mem*/an.2018-04-15_00.00.00.nc are updated, and there are observation feedback files generated in dbOut/mem*, namely obsout_da_*.h5.
Different from LETKF, which does posterior inflation internally, an offline posterior RTPP inflation can be run after EDA analyses with:
$ qsub run_eda_rtpp.csh
The posterior RTPP inflation will update all EDA analysis members an/mem*/an.2018-04-15_00.00.00.nc.
This EDA test is ensemble of 3DEnVar, you could also run ensemble of 3DVar, hybrid-3DEnVar, or even hybrid-4DEnVar, with a corresponding yaml setting. A yaml file named eda_3dhyb.yaml is provided for you to play with.
We also prepared some NCL scripts for graphics to see the impacts of DA. Copy them to the working directory.
$ cp ../ncl_scripts/plot_eda_*.ncl .
Check whether the NCL library is available:
$ which ncl
If the system fails to find ncl, you may need to load the NCL module:
$ module reset
$ module load ncl/6.6.2
We can plot the analysis increment of ensemble mean with the script plot_eda_ensmean_anainc.ncl. Specify the diag_var, diag_lev (for 3-D variables), and nens (if run fewer EDA members) in the script,
diag_var="surface_pressure"
diag_lev=10
nens = 10
and then run:
$ ncl plot_eda_ensmean_anainc.ncl
This will produce a figure named eda_ens_mean_inc_surface_pressure_lev0.png.
With the plot_eda_ens_spread.ncl, we can plot the ensemble spread of both background and analysis. Similar to the ensemble mean increments, just change the diag_var and diag_lev in the NCL script, and run plot_eda_ens_spread.ncl:
$ ncl plot_eda_ens_spread.ncl
This will produce a figure eda_ens_spread_surface_pressure_lev0.png.
We can also plot the vertical profiles of RMSE, BIAS, and total spread of the Innovations (Obs minus Bak) and Residuals (Obs minus Ana) with plot_eda_vert_profile_err.ncl. We need to specify the diag_obs, diag_var and nens (if run fewer EDA members) in the script:
diag_obs="sondes"
diag_var="windEastward"
nens = 10
Running:
$ ncl plot_eda_vert_profile_err.ncl
will produce a figure file eda_vert_profile_sondes_windEastward.png:
In this subsection, we will further run a single-time LETKF analysis. First, we need to create the working directory for LETKF and move to that directory:
$ mkdir -p /glade/derecho/scratch/$USER/mpas_jedi_tutorial/letkf
$ cd /glade/derecho/scratch/$USER/mpas_jedi_tutorial/letkf
Set environment variable for the mpas-bundle directory,
$ export bundle_dir=/glade/derecho/scratch/liuz/mpas_bundle_v3_public_gnuSP
Similar to EDA, we first prepare necessary files for LETKF.
$ ln -sf ../MPAS_namelist_stream_physics_files/*BL .
$ ln -sf ../MPAS_namelist_stream_physics_files/*DATA .
$ ln -sf ../MPAS_namelist_stream_physics_files/x1.10242.graph.info.part.64 .
$ ln -sf ../MPAS_namelist_stream_physics_files/streams.atmosphere_240km .
$ ln -sf ../MPAS_namelist_stream_physics_files/stream_list.atmosphere.* .
$ ln -sf ../MPAS_namelist_stream_physics_files/namelist.atmosphere_240km .
$ ln -sf ../MPAS_namelist_stream_physics_files/x1.10242.invariant.nc .
$ ln -sf ../background/2018041418/mpasout.2018-04-15_00.00.00.nc ./templateFields.10242.nc
$ ln -sf ${bundle_dir}/code/mpas-jedi/test/testinput/namelists/geovars.yaml ./
$ ln -sf ${bundle_dir}/code/mpas-jedi/test/testinput/namelists/keptvars.yaml ./
$ ln -sf ${bundle_dir}/code/mpas-jedi/test/testinput/obsop_name_map.yaml ./
Now prepare 'bg' and 'an' folders under letkf.
$ cp ../eda/prep_an_bg_directories.csh .
$ ./prep_an_bg_directories.csh
$ cp -r an/mem001 an/mem000
The reason of creating an/mem000 is that LETKF also writes out the mean of ensemble analyses under an/mem000.
For observation-related directories, we need to create dbIn, dbOut, and dbAna directories, and link IODA observation files to dbIn directory:
$ mkdir dbIn dbOut dbAna
$ cd dbIn
$ ln -sf ../../obs_ioda_pregenerated/2018041500/aircraft_obs_2018041500.h5 ./
$ ln -sf ../../obs_ioda_pregenerated/2018041500/gnssro_obs_2018041500.h5 ./
$ ln -sf ../../obs_ioda_pregenerated/2018041500/sfc_obs_2018041500.h5 ./
$ ln -sf ../../obs_ioda_pregenerated/2018041500/satwind_obs_2018041500.h5 ./
$ ln -sf ../../obs_ioda_pregenerated/2018041500/sondes_obs_2018041500.h5 ./
$ cd ..
link the yaml files and scripts for LETKF:
$ ln -sf ../MPAS_JEDI_yamls_scripts/letkf_observer.yaml ./
$ ln -sf ../MPAS_JEDI_yamls_scripts/letkf_solver.yaml ./
$ ln -sf ../MPAS_JEDI_yamls_scripts/letkf_solver_vloc.yaml ./
$ ln -sf ../MPAS_JEDI_yamls_scripts/letkf_diagoma.yaml ./
$ ln -sf ../MPAS_JEDI_yamls_scripts/run_letkf_observer.csh ./
$ ln -sf ../MPAS_JEDI_yamls_scripts/run_letkf_solver.csh ./
$ ln -sf ../MPAS_JEDI_yamls_scripts/run_letkf_diagoma.csh ./
letkf_observer.yaml, letkf_solver.yaml, and letkf_diagoma.yaml are used for the three steps of LETKF. letkf_observer.yaml is for calculating the ensemble HofX of background and performing quality control based on the ensemble mean. letkf_solver.yaml is for updating the ensemble mean and ensemble members. lketkf_diagoma.yaml is for calculating the ensemble HofX of analysis for OMA diagnosis.
We first submit run_letkf_observer.csh:
$ qsub run_letkf_observer.csh
This observer step will produce obsout_da_*.h5 files in dbOut. In each df5 file, we can see groups of hofx_0_*, hofx_y_mean_xb0, ombg, EffectiveQC0, EffectiveError0.
Then, submit run_letkf_solver.csh:
$ qsub run_letkf_solver.csh
This solver step will update an/mem*/an.2018-04-15_00.00.00.nc files, including the ensemble mean of analysis an/mem000/an.2018-04-15_00.00.00.nc.
Last, submit run_letkf_diagoma.csh:
$ qsub run_letkf_diagoma.csh
The analysis-diagnos step will generate obsout_*.h5 files in dbAna, similar to the files in dbOut but for HofX of ensemble analysis.
The previous run only tests the horizontal localization. If we want to also employ the vertical localization, just change
_height-based obs localizations: &heightAndHorizObsLoc
<<: *2DLETKFheightAndHorizObsLoc
to
_height-based obs localizations: &heightAndHorizObsLoc
<<: *3DLETKFheightAndHorizObsLoc
in letkf_solver.yaml, and then resubmit run_letkf_solver.csh and run_letkf_diagoma.csh. Backup an and dbAna directories before another run.
Similar to the EDA test, we also prepared some NCL scripts to plot LETKF results.
$ cp ../ncl_scripts/plot_letkf_*.ncl .
Plotting the analysis increment of ensemble mean with plot_letkf_ensmean_anainc.ncl. We need to specify the diag_var and diag_lev (for 3-D variables) in the NCL script.
diag_var="surface_pressure"
diag_lev=10
Then do:
$ ncl plot_letkf_ensmean_anainc.ncl
This will produce a figure letkf_ens_mean_inc_surface_pressure_lev0.png.
For plotting ensemble spread of both background and analysis, you could change diag_var and diag_lev in plot_letkf_ens_spread.ncl before running:
$ ncl plot_letkf_ens_spread.ncl
This will produce a figure letkf_ens_spread_surface_pressure_lev0.png.
For plotting the vertical profiles of RMSE, BIAS, and total spread of the Innovations (Obs minus Bak) and Residuals (Obs minus Ana) with plot_letkf_vert_profile_err.ncl, we could specify the diag_obs and diag_var in the script:
diag_obs="sondes"
diag_var="windEastward"
then run:
$ ncl plot_letkf_vert_profile_err.ncl
This will produce a figure letkf_vert_profile_sondes_windEastward.png.
