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 20 3DEnVar members.
First, We need to create the working directory for EDA test and move to that directory:
$ mkdir -p /glade/scratch/$USER/mpas_jedi_tutorial/eda
$ cd /glade/scratch/$USER/mpas_jedi_tutorial/eda
Then, we need to prepare necessary files for EDA test.
Link or copy MPAS physics-related data and table, MPAS graph, MPAS streams, and MPAS namelist files:
$ ln -sf ../MPAS_namelist_stream_physics_files/*TBL .
$ ln -sf ../MPAS_namelist_stream_physics_files/*DATA .
$ ln -sf ../MPAS_namelist_stream_physics_files/*DBL .
$ ln -sf ../MPAS_namelist_stream_physics_files/VERSION .
$ ln -sf ../MPAS_namelist_stream_physics_files/x1.10242.graph.info.part.36 .
$ 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 ../background/2018041418/x1.10242.init.2018-04-14_18.00.00.nc ./static.10242.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 ../mpas_bundle_v2_SP/code/mpas-jedi/test/testinput/namelists/geovars.yaml ./
$ ln -sf ../mpas_bundle_v2_SP/code/mpas-jedi/test/testinput/namelists/keptvars.yaml ./
$ ln -sf ../mpas_bundle_v2_SP/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 11 12 13 14 15 16 17 18 19 20)
mkdir -p bg/mem0${imem}
cp ../ensemble/2018041418/${imem}/mpasout.2018-04-15_00.00.00.nc bg/mem0${imem}
end
cp -r bg an
rename mpasout bg bg/mem*/*.nc
rename mpasout an an/mem*/*.nc
cp -r an/mem001 an/mem000
After running prep_an_bg_directories.csh, the an and bg directories and files are prepared.
$ csh ./prep_an_bg_directories.csh
Next step, we create directory dbIn and link the own-generated observation IODA files to it. And also, we need to create directory dbOut to store the output observation files of all EDA members. In case the IODA were not genereated successfully in obs_ioda, we may link the IODA files from the obs_ioda_pregenerated.
$ mkdir dbIn
$ cd dbIn
$ ln -sf ../../obs_ioda/2018041500/aircraft_obs_2018041500.h5 ./
$ ln -sf ../../obs_ioda/2018041500/gnssro_obs_2018041500.h5 ./gnssrorefncep_obs_2018041500.h5
$ ln -sf ../../obs_ioda/2018041500/sfc_obs_2018041500.h5 ./
$ ln -sf ../../obs_ioda/2018041500/satwind_obs_2018041500.h5 ./
$ ln -sf ../../obs_ioda/2018041500/sondes_obs_2018041500.h5 ./
$ cd /glade/scratch/$USER/mpas_jedi_tutorial/eda
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 11 12 13 14 15 16 17 18 19 20)
mkdir -p dbOut/mem0${imem}
end
By running the prep_dbout_directory.csh with the command:
$ 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 ./bump_loc
$ ln -sf ../B_Matrix/* ./
Last, copy related yaml files and scripts for EDA running:
$ cp ../MPAS_JEDI_yamls_scripts/eda_3denvar.yaml ./
$ cp ../MPAS_JEDI_yamls_scripts/eda_hybrid.yaml ./
$ cp ../MPAS_JEDI_yamls_scripts/eda_rtpp.yaml ./
$ cp ../MPAS_JEDI_yamls_scripts/run_eda_members.csh ./
$ cp ../MPAS_JEDI_yamls_scripts/run_eda_rtpp.csh ./
The eda_3denvar.yaml is a template yaml file for all 3DEnVar EDA members. The different seetings among all EDA members can be found by looking for {{imem}} and {{IMEM}} in the yaml file.
To run the EDA test, make sure the main_dir in run_eda_members.csh is set to /glade/scratch/$USER/mpas_jedi_tutorial/eda. And we need also specify the yaml_file as eda_3denvar.yaml.
If we run all 20 EDA members serialy in one job, it will take about 45 minutes. In run_eda_members.csh, the command line foreach i_ens (1 2 3 4 5) indicates taht EDA members from 1 to 5 will be done in a single job. We can change it to foreach i_ens (6 7 8 9 10) to do the EDA members from 6 to 10, and simialr thing can be done for members from 11 to 15, and from 16 to 20. In that case, finishing 20-member EDA analysis needs four pareall jobs, and each job will take about 12 minutes. To save time, we can only do a 5-member EDA test.
Sumbit the job script run_eda_members.csh to do the listed EDA members with the command:
$ 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 new observation files generated in dbOut/mem*, namely obsout_da_*.h5. In each observation file, we will find new groups are generated, like hofx0, hofx1, ombg, oman, EffectiveQC0, EffectiveQC1, EffectiveError0, EffectiveError1.
Different like LETKF that does posterior inflation itself. A posterior RTPP inflation is done using an external executable for EDA. We can submit the script run_eda_rtpp.csh:
$ 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, and also output analysis mean as an/mem000/an.2018-04-15_00.00.00.nc.
In previous EDA test, we employ the 3DEnVar for each EDA member. We can also use other variational methods for EDA, like 3DVar, hybrid-3DEnVar, and also 4DEnVar. Try to modify the eda_3denvar.yaml to let it work with hybrid-3DEnVar capability. A yaml file named eda_hybrid.yaml is given for reference. After modifying eda_3denvar.yaml, we can test 20-member hybrid-3DEnVar EDA members. If we want to keep the results of 3DEnVar EDA, backup the an and dbOut directories before submitting the jobs for the new test.
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 libray is available:
$ which ncl
If the system fails to find the executable ncl, it may be because the spack-stack modules for GNU are in use. In that case, we can use the following commands to load the NCL module:
$ module reset
$ module load intel/19.0.5
$ module load ncl/6.6.2
Try which ncl again and the system will find the ncl executable in /glade/u/apps/ch/opt/ncl/6.6.2/intel/19.0.5/bin/ncl.
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 = 20
We can draw the increments with the command:
$ ncl plot_eda_ensmean_anainc.ncl
and then, we will get the figure 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
and after running we will get 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 the script 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 = 20
After running the script:
$ ncl plot_eda_vert_profile_err.ncl.
we can get the figure file eda_vert_profile_sondes_windEastward.png:
In this subsection, we will further run the single-time LETKF analysis. First, we need to create the working directory for LETKF test and move to that directory:
$ mkdir -p /glade/scratch/$USER/mpas_jedi_tutorial/letkf
$ cd /glade/scratch/$USER/mpas_jedi_tutorial/letkf
Similar to the EDA run, we first prepare necessary files for LETKF running. The following files for LETKF are exactly the same as that in EDA, including MPAS physics-related data and table, MPAS graph, MPAS streams, MPAS namelist, MPAS 2-stream files, static yaml files for MPAS-JEDI, as well as the background and analysis directory and files. Just follow the instructions in section 6.1 before creating the observations directories.
There are several differences for the observation directory, yaml file, and scripts in LETKF test.
For the observation directory, we need to create dbIn, dbOut, and dbAna directories, and copy/link IODA observation files to dbIn directory:
$ mkdir dbIn dbOut dbAna
$ cd dbIn
$ ln -sf ../../obs_ioda/2018041500/aircraft_obs_2018041500.h5 ./
$ ln -sf ../../obs_ioda/2018041500/gnssro_obs_2018041500.h5 ./gnssrorefncep_obs_2018041500.h5
$ ln -sf ../../obs_ioda/2018041500/sfc_obs_2018041500.h5 ./
$ ln -sf ../../obs_ioda/2018041500/satwind_obs_2018041500.h5 ./
$ ln -sf ../../obs_ioda/2018041500/sondes_obs_2018041500.h5 ./
$ cd ..
copy the yaml files and scripts for LETKF:
$ cp ../MPAS_JEDI_yamls_scripts/letkf_observer.yaml ./
$ cp ../MPAS_JEDI_yamls_scripts/letkf_solver.yaml ./
$ cp ../MPAS_JEDI_yamls_scripts/letkf_solver_vloc.yaml ./
$ cp ../MPAS_JEDI_yamls_scripts/letkf_diagoma.yaml ./
$ cp ../MPAS_JEDI_yamls_scripts/run_letkf_observer.csh ./
$ cp ../MPAS_JEDI_yamls_scripts/run_letkf_solver.csh ./
$ cp ../MPAS_JEDI_yamls_scripts/run_letkf_diagoma.csh ./
The letkf_observer.yaml, letkf_solver.yaml, and letkf_diagoma.yaml are used for the three steps of LETKF. The letkf_observer.yaml is to used to calculate the ensemble HofX and do quality control based on the ensemble mean. The letkf_solver.yaml is for the LETKF solver step, which updates both the ensemble mean and ensemble members. The lketkf_diagoma.yaml is used to calculate the ensemble HofX of analysis for OMA diagnosis.
To run LETKF test, we first submit the script run_letkf_observer.csh:
$ qsub run_letkf_observer.csh
If this step runs successfully, we will get obsout_da_*.h5 files in dbOut. In each hdf5 file, we can see new groups of hofx_0_*, hofx_y_mean_xb0, ombg, EffectiveQC0, EffectiveError0 are generated.
Then, we sumbit the script run_letkf_solver.csh:
$ qsub run_letkf_solver.csh
If this step runs successfull, we will get updated 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, we submit the script run_letkf_diagoma.csh:
$ qsub run_letkf_diagoma.csh
The final step will generate obsout_*.h5 files in dbAna, similar to the files in dbOut but for HofX of ensemble analysis.
The previous run only test the horizontal localization. If we want to also employ the vertical localization, just chage
_height-based obs localizations: &heightAndHorizObsLoc
<<: *2DLETKFheightAndHorizObsLoc
to
_height-based obs localizations: &heightAndHorizObsLoc
<<: *3DLETKFheightAndHorizObsLoc
in letkf_solver.yaml, and then resumbit run_letkf_solver.csh and run_letkf_diagoma.csh. Backup an and dbAna directories If we want to store the results of LETKF with horizontal localization only.
We also prepared some NCL scripts to plot graphics to see the impacts of DA. Copy them to the working directory. They are plot_letkf_ensmean_anainc.ncl, plot_letkf_ens_spread.ncl, and plot_letkf_vert_profile_err.ncl.
$ cp ../ncl_scripts/plot_letkf_*.ncl .
Be sure that the NCL module is loaded.
$ module reset
$ module load intel/19.0.5
$ module load ncl/6.6.2
We can plot the analysis increment of ensemble mean with the script 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
Plot the ensemble mean increments with the command:
$ ncl plot_letkf_ensmean_anainc.ncl
and then, we will get the figure letkf_ens_mean_inc_surface_pressure_lev0.png.
With the plot_letkf_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_letkf_ens_spread.ncl:
$ ncl plot_letkf_ens_spread.ncl
and after running we will get figure letkf_ens_spread_surface_pressure_lev0.png.
We can also get the vertical profiles of RMSE, BIAS, and total spread of the Innovations (Obs minus Bak) and Residuals (Obs minus Ana) with the script plot_letkf_vert_profile_err.ncl. We need to specify the diag_obs and diag_var in the script:
diag_obs="sondes"
diag_var="windEastward"
After running the script:
$ ncl plot_letkf_vert_profile_err.ncl
we can get the figure file letkf_vert_profile_sondes_windEastward.png.
