Skip to content

Commit

Permalink
Add the Multigrid Beta Filter (MGBF) for ensemble localization (#699) (
Browse files Browse the repository at this point in the history
…#700)

**DUE DATE for merger of this PR into `develop` is 3/29/2024 (six weeks
after PR creation).**

**Description**

Resolves #699

This PR is to add the option to apply Multigrid Beta Filter (MGBF;
[Purser et al. 2022](https://doi.org/10.1175/MWR-D-20-0405.1)) for
ensemble localization instead of Recursive Filter (RF). This work
includes to add an initial version of the MGBF as a subdirectory in GSI.

To apply the MGBF, set "l_mgbf_loc=true" in the namelist and
additionally input "mgbf_loc01.nml". (In Scale/Variable-Dependent
Localization, input also "mgbf_locXX.nml" (XX=02,03,...) with the same
number of grid points.)

<details>

**<summary>How to set MGBF parameters in mgbf_locXX.nml</summary>**

An example of mgbf_locXX.nml is as follows:

```
&PARAMETERS_MGBETA
  mg_ampl01=1.125,            ! length of vertical beta filter (standard deviation; filter grid unit)
  mg_ampl02=2.4,              ! length of horizontal beta filter (standard deviation; filter grid unit)
  mg_ampl03=0.85,             ! length of 3D beta filter (standard deviation; filter grid unit)
  mg_weig1=0.,                ! weight of generation 1
  mg_weig2=0.,                ! weight of generation 2
  mg_weig3=0.,                ! weight of generation 3
  mg_weig4=1.,                ! weight of generation 4
  hx=5,                       ! number of halo grid points in x-direction
  hy=5,                       ! number of halo grid points in y-direction
  hz=3,                       ! number of halo grid points in z-direction
  p=2,                        ! beta filter exponent
  mgbf_line=.false.,          ! set false except for mgbf_proc=2,4,7
  mgbf_proc=8,                ! 1-2: 3D filter; 3-5: 2D filter for static B; 6-8: 2D filter for localization (1,3,6: radial filter; 2,4,7: line filter; 5,8: isotropic line filter)
  lm_a=65,                    ! number of vertical layers in analysis grid
  lm=33,                      ! number of vertical layers in filter grid
  km2=0,                      ! number of 2D variables (set 0 for localization)
  km3=1,                      ! number of 3D variables (set 1 for localization)
  n_ens=30,                   ! ensemble size
  l_loc=.true.,               ! set true in localization
  l_filt_g1=.false.,          ! set false in skipping generation 1
  l_lin_vertical=.true.,      ! set true in applying linear vertical interpolation for analysis-filter mapping
  l_lin_horizontal=.true.,    ! set true in applying linear horizontal interpolation for analysis-filter mapping
  l_quad_horizontal=.false.,  ! set true in applying quadratic horizontal interpolation for analysis-filter mapping
  l_new_map=.true.,           ! set true in applying efficient vertical interpolation for analysis-filter mapping
  l_vertical_filter=.true.,   ! set true in applying vertical beta filter outside 2D filter
  ldelta=.false.,             ! (not used)
  lquart=.false.,             ! set true in applying quadratic horizontal interpolation for up/down-sending
  lhelm=.false.,              ! set true in applying Helmholtz differential operator for weighting
  nm0=1975,                   ! number of analysis grid points in x-direction
  mm0=1350,                   ! number of analysis grid points in y-direction
  gm_max=4,                   ! highest generation (max: 4)
  nxPE=79,                    ! number of MPI processors in x-direction
  nyPE=54,                    ! number of MPI processors in y-direction
  im_filt=8,                  ! number of filter grid points in each MPI processor in x-direction
  jm_filt=8,                  ! number of filter grid points in each MPI processor in y-direction
 /
```

Here, to make the result of MGBF-based localization similar to RF-based
one, we can set the beta filter length ( mg_ampl0[12] ) from the
recursive filter length in the GSI namelist ( s_ens_[vh] ) as:

- $\text{mg\\_ampl01} = \left[\text{s\\_ens\\_v (grid unit)} *
\frac{1}{\sqrt{2}} * \frac{\text{lm}-1}{\text{lm\\_a}-1} \right]^2$
- $\text{mg\\_ampl02} = \left[\frac{\text{s\\_ens\\_h
(km)}}{\text{analysis grid interval (km)}} * \frac{1}{\sqrt{2}} *
\frac{\text{im\\_filt} * \text{nxPE}}{\text{nm0}} * \frac{1}{2} *
\frac{1}{2} * \frac{1}{2} \right]^2$ (in case mg_weig[1-4]=[0,0,0,1])

Please note there are some limitations for the other MGBF parameters
such as:

- The number of MPI processors input in GSI should be nxPE x nyPE 
- (nm0, mm0, lm_a) should be the same as the GSI analysis grid
- nm0 should be divisible by nxPE
- mm0 should be divisible by nyPE
- nm0 / nxPE = mm0 / nyPE

</details>

<details>

**<summary>How to run RRFS regression tests with MGBF-based
localization</summary>**

Change settings in regression/ as follows, and run Test#3
(rrfs_3denvar_glbens)

```diff
diff --git a/regression/regression_namelists.sh b/regression/regression_namelists.sh
index 7ca183ef3..671d028ff 100755
--- a/regression/regression_namelists.sh
+++ b/regression/regression_namelists.sh
@@ -457,7 +457,7 @@ OBS_INPUT::
    beta_s0=0.15,s_ens_h=110,s_ens_v=3,
    regional_ensemble_option=1,
    pseudo_hybens = .false.,
-   grid_ratio_ens = 3,
+   grid_ratio_ens = 5.1,
    l_ens_in_diff_time=.true.,
    ensemble_path='',
    i_en_perts_io=1,
@@ -465,6 +465,7 @@ OBS_INPUT::
    fv3sar_bg_opt=0,
    readin_localization=.true.,
    ens_fast_read=.false.,
+   l_mgbf_loc=.true.,
  /
  &RAPIDREFRESH_CLDSURF
    dfi_radar_latent_heat_time_period=20.0,
```

```diff
diff --git a/regression/regression_namelists_db.sh b/regression/regression_namelists_db.sh
index e03917e88..36b8b6a22 100755
--- a/regression/regression_namelists_db.sh
+++ b/regression/regression_namelists_db.sh
@@ -438,7 +438,7 @@ OBS_INPUT::
    beta_s0=0.15,s_ens_h=110,s_ens_v=3,
    regional_ensemble_option=1,
    pseudo_hybens = .false.,
-   grid_ratio_ens = 3,
+   grid_ratio_ens = 5.1,
    l_ens_in_diff_time=.true.,
    ensemble_path='',
    i_en_perts_io=1,
@@ -446,6 +446,7 @@ OBS_INPUT::
    fv3sar_bg_opt=0,
    readin_localization=.true.,
    ens_fast_read=.false.,
+   l_mgbf_loc=.true.,
  /
  &RAPIDREFRESH_CLDSURF
    dfi_radar_latent_heat_time_period=20.0,
```

```diff
diff --git a/regression/regression_param.sh b/regression/regression_param.sh
index 2ac615fc4..6186acdbb 100755
--- a/regression/regression_param.sh
+++ b/regression/regression_param.sh
@@ -87,23 +87,23 @@ case $regtest in
     rrfs_3denvar_glbens)

         if [[ "$machine" = "Hera" ]]; then
-           topts[1]="0:15:00" ; popts[1]="20/1/"  ; ropts[1]="/1"
-           topts[2]="0:15:00" ; popts[2]="20/2/"  ; ropts[2]="/1"
+           topts[1]="0:15:00" ; popts[1]="30/2/"  ; ropts[1]="/1"
+           topts[2]="0:15:00" ; popts[2]="15/4/"  ; ropts[2]="/1"
         elif [[ "$machine" = "Orion" ]]; then
-           topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1"
-           topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/2"
+           topts[1]="0:15:00" ; popts[1]="30/2/"  ; ropts[1]="/1"
+           topts[2]="0:15:00" ; popts[2]="15/4/"  ; ropts[2]="/2"
         elif [[ "$machine" = "Hercules" ]]; then
-           topts[1]="0:15:00" ; popts[1]="20/1/" ; ropts[1]="/1"
-           topts[2]="0:15:00" ; popts[2]="20/2/" ; ropts[2]="/2"
+           topts[1]="0:15:00" ; popts[1]="30/2/"  ; ropts[1]="/1"
+           topts[2]="0:15:00" ; popts[2]="15/4/"  ; ropts[2]="/2"
         elif [[ "$machine" = "Jet" ]]; then
-           topts[1]="0:15:00" ; popts[1]="20/1/"  ; ropts[1]="/1"
-           topts[2]="0:15:00" ; popts[2]="20/2/"  ; ropts[2]="/1"
+           topts[1]="0:15:00" ; popts[1]="15/4/"  ; ropts[1]="/1"
+           topts[2]="0:15:00" ; popts[2]="10/6/"  ; ropts[2]="/1"
         elif [[ "$machine" = "Gaea" ]]; then
-           topts[1]="0:15:00" ; popts[1]="18/1/"  ; ropts[1]="/1"
-           topts[2]="0:15:00" ; popts[2]="18/2/"  ; ropts[2]="/1"
+           topts[1]="0:15:00" ; popts[1]="15/4/"  ; ropts[1]="/1"
+           topts[2]="0:15:00" ; popts[2]="10/6/"  ; ropts[2]="/1"
         elif [[ "$machine" = "wcoss2" ]]; then
-           topts[1]="0:15:00" ; popts[1]="64/1/" ; ropts[1]="/1"
-           topts[2]="0:15:00" ; popts[2]="128/2/" ; ropts[2]="/1"
+           topts[1]="0:15:00" ; popts[1]="60/1/"  ; ropts[1]="/1"
+           topts[2]="0:15:00" ; popts[2]="30/2/"  ; ropts[2]="/1"
         fi

         if [ "$debug" = ".true." ] ; then
```

```diff
diff --git a/regression/rrfs_3denvar_glbens.sh b/regression/rrfs_3denvar_glbens.sh
index af5da5117..04fd73d57 100755
--- a/regression/rrfs_3denvar_glbens.sh
+++ b/regression/rrfs_3denvar_glbens.sh
@@ -272,6 +272,46 @@ $gsi_namelist

 EOF

+cat << EOF > mgbf_loc01.nml
+&PARAMETERS_MGBETA
+  mg_ampl01=1.125,            ! length of vertical beta filter (standard deviation; filter grid unit)
+  mg_ampl02=1.615,            ! length of horizontal beta filter (standard deviation; filter grid unit)
+  mg_ampl03=0.85,             ! length of 3D beta filter (standard deviation; filter grid unit)
+  mg_weig1=0.,                ! weight of generation 1
+  mg_weig2=1.,                ! weight of generation 2
+  mg_weig3=0.,                ! weight of generation 3
+  mg_weig4=0.,                ! weight of generation 4
+  hx=4,                       ! number of halo grid points in x-direction
+  hy=4,                       ! number of halo grid points in y-direction
+  hz=3,                       ! number of halo grid points in z-direction
+  p=2,                        ! beta filter exponent
+  mgbf_line=.false.,          ! set false except for mgbf_proc=2,4,7
+  mgbf_proc=8,                ! 1-2: 3D filter; 3-5: 2D filter for static B; 6-8: 2D filter for localization (1,3,6: radial filter; 2,4,7: line filter; 5,8: isotropic line filter)
+  lm_a=65,                    ! number of vertical layers in analysis grid
+  lm=33,                      ! number of vertical layers in filter grid
+  km2=0,                      ! number of 2D variables (set 0 for localization)
+  km3=1,                      ! number of 3D variables (set 1 for localization)
+  n_ens=10,                   ! ensemble size
+  l_loc=.true.,               ! set true in localization
+  l_filt_g1=.false.,          ! set false in skipping generation 1
+  l_lin_vertical=.true.,      ! set true in applying linear vertical interpolation for analysis-filter mapping
+  l_lin_horizontal=.true.,    ! set true in applying linear horizontal interpolation for analysis-filter mapping
+  l_quad_horizontal=.false.,  ! set true in applying quadratic horizontal interpolation for analysis-filter mapping
+  l_new_map=.true.,           ! set true in applying efficient vertical interpolation for analysis-filter mapping
+  l_vertical_filter=.true.,   ! set true in applying vertical beta filter outside 2D filter
+  ldelta=.false.,             ! (not used)
+  lquart=.false.,             ! set true in applying quadratic horizontal interpolation for up/down-sending
+  lhelm=.false.,              ! set true in applying Helmholtz differential operator for weighting
+  nm0=40,                     ! number of analysis grid points in x-direction
+  mm0=24,                     ! number of analysis grid points in y-direction
+  gm_max=2,                   ! highest generation (max: 4)
+  nxPE=10,                    ! number of MPI processors in x-direction
+  nyPE=6,                     ! number of MPI processors in y-direction
+  im_filt=4,                  ! number of filter grid points in each MPI processor in x-direction
+  jm_filt=4,                  ! number of filter grid points in each MPI processor in y-direction
+ /
+EOF
+
 # Copy executable and fixed files to $tmpdir
 if [[ $exp == *"updat"* ]]; then
    $ncp $gsiexec_updat  ./gsi.x
```

</details>

**Type of change**

Please delete options that are not relevant.

- [ ] Bug fix (non-breaking change which fixes an issue)
- [x] New feature (non-breaking change which adds functionality)
- [ ] Breaking change (fix or feature that would cause existing
functionality to not work as expected)
- [ ] This change requires a documentation update

**How Has This Been Tested?**

EnVar for NA-domain RRFS was tested with "mgbf_locXX.nml" (XX=01) shown
above on Orion. The resulting analysis increment was similar to the
original and the computation time for localization became short.
  
**Checklist**

- [x] My code follows the style guidelines of this project
- [x] I have performed a self-review of my own code
- [x] I have commented my code, particularly in hard-to-understand areas
- [x] New and existing tests pass with my changes
- [x] Any dependent changes have been merged and published

Co-authored-by: Sho Yokota <syokota@Orion-login-1.HPC.MsState.Edu>
  • Loading branch information
shoyokota and Sho Yokota committed Mar 26, 2024
1 parent f7e93ab commit 1cc934e
Show file tree
Hide file tree
Showing 36 changed files with 26,816 additions and 141 deletions.
2 changes: 2 additions & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ endif()
option(OPENMP "Enable OpenMP Threading" OFF)
option(ENABLE_MKL "Use MKL for LAPACK implementation (if available)" ON)
option(BUILD_GSDCLOUD "Build GSD Cloud Analysis Library" OFF)
option(BUILD_MGBF "Build MGBF Library" ON)
option(BUILD_GSI "Build GSI" ON)
option(BUILD_ENKF "Build EnKF" ON)
option(BUILD_REG_TESTING "Build the Regression Testing Suite" OFF)
Expand All @@ -37,6 +38,7 @@ option(BUILD_REG_TESTING "Build the Regression Testing Suite" OFF)
message(STATUS "OPENMP ................. ${OPENMP}")
message(STATUS "ENABLE_MKL ............. ${ENABLE_MKL}")
message(STATUS "BUILD_GSDCLOUD ......... ${BUILD_GSDCLOUD}")
message(STATUS "BUILD_MGBF ............. ${BUILD_MGBF}")
message(STATUS "BUILD_GSI .............. ${BUILD_GSI}")
message(STATUS "BUILD_ENKF ............. ${BUILD_ENKF}")
message(STATUS "BUILD_REG_TESTING ...... ${BUILD_REG_TESTING}")
Expand Down
1 change: 1 addition & 0 deletions INSTALL.md
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,7 @@ CMake allows for various options that can be specified on the command line via `
| `OPENMP` | Enable OpenMP Threading (`OFF`) |
| `ENABLE_MKL` | Use MKL (`ON`), If not found use LAPACK |
| `BUILD_GSDCLOUD` | Build GSD Cloud Library (`OFF`) |
| `BUILD_MGBF` | Build MGBF Library (`ON`) |
| `BUILD_GSI` | Build GSI library and executable (`ON`) |
| `BUILD_ENKF` | Build EnKF library and executable (`ON`) |
| `BUILD_REG_TESTING` | Enable Regression Testing (`ON`) |
Expand Down
5 changes: 5 additions & 0 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,11 @@ if(BUILD_GSDCLOUD)
add_subdirectory(GSD)
endif()

if(BUILD_MGBF)
message(STATUS "Building MGBF library")
add_subdirectory(mgbf)
endif()

if(BUILD_GSI)
message(STATUS "Building GSI")
add_subdirectory(gsi)
Expand Down
15 changes: 15 additions & 0 deletions src/gsi/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ endif()
option(OPENMP "Enable OpenMP Threading" OFF)
option(ENABLE_MKL "Use MKL for LAPACK implementation (if available)" ON)
option(USE_GSDCLOUD "Use GSD Cloud Analysis library" OFF)
option(USE_MGBF "Use MGBF library" ON)

set(GSI_VALID_MODES "GFS" "Regional")
set(GSI_MODE "GFS" CACHE STRING "Choose the GSI Application.")
Expand All @@ -43,6 +44,7 @@ endif()
message(STATUS "GSI: OPENMP ................. ${OPENMP}")
message(STATUS "GSI: ENABLE_MKL ............. ${ENABLE_MKL}")
message(STATUS "GSI: USE_GSDCLOUD ........... ${USE_GSDCLOUD}")
message(STATUS "GSI: USE_MGBF ............... ${USE_MGBF}")
message(STATUS "GSI: GSI_MODE ............... ${GSI_MODE}")

# Dependencies
Expand Down Expand Up @@ -87,6 +89,13 @@ if(USE_GSDCLOUD)
endif()
endif()

# MGBF library dependency
if(USE_MGBF)
if(NOT TARGET mgbf)
find_package(mgbf REQUIRED)
endif()
endif()

# Get compiler flags for the GSI application
include(gsiapp_compiler_flags)

Expand Down Expand Up @@ -158,6 +167,12 @@ if(USE_GSDCLOUD)
endif()
target_link_libraries(gsi_fortran_obj PUBLIC gsdcloud::gsdcloud)
endif()
if(USE_MGBF)
if(TARGET mgbf)
add_dependencies(gsi_fortran_obj mgbf)
endif()
target_link_libraries(gsi_fortran_obj PUBLIC mgbf::mgbf)
endif()
if(OpenMP_Fortran_FOUND)
target_link_libraries(gsi_fortran_obj PRIVATE OpenMP::OpenMP_Fortran)
endif()
Expand Down
18 changes: 16 additions & 2 deletions src/gsi/gsimod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -161,7 +161,7 @@ module gsimod
ntotensgrp,nsclgrp,naensgrp,ngvarloc,ntlevs_ens,naensloc, &
r_ensloccov4tim,r_ensloccov4var,r_ensloccov4scl,l_timloc_opt,&
vdl_scale,vloc_varlist,&
global_spectral_filter_sd,assign_vdl_nml,parallelization_over_ensmembers
global_spectral_filter_sd,assign_vdl_nml,parallelization_over_ensmembers,l_mgbf_loc
use hybrid_ensemble_parameters,only : l_both_fv3sar_gfs_ens,n_ens_gfs,n_ens_fv3sar,weight_ens_gfs,weight_ens_fv3sar
use rapidrefresh_cldsurf_mod, only: init_rapidrefresh_cldsurf, &
dfi_radar_latent_heat_time_period,metar_impact_radius,&
Expand Down Expand Up @@ -529,6 +529,7 @@ module gsimod
! - innov_use_model_fed=.true. : Use FED from BG to calculate innovation.
! this requires if_model_fed=.true.
! it works either an EnVar DA run or a GSI observer run.
! 02-20-2024 yokota - add MGBF-based localization
!
!EOP
!-------------------------------------------------------------------------
Expand Down Expand Up @@ -1452,6 +1453,7 @@ module gsimod
! ^ ^ ^ ^ ^
! s_ens_h = v1L1 v2L1 v3L1 v1L2 v2L2
! Then localization lengths will be assigned as above.
! l_mgbf_loc - if true, multi-grid beta filter is used for localization instead of recursive filter
!
namelist/hybrid_ensemble/l_hyb_ens,uv_hyb_ens,q_hyb_ens,aniso_a_en,generate_ens,n_ens,&
l_both_fv3sar_gfs_ens,n_ens_gfs,n_ens_fv3sar,weight_ens_gfs,weight_ens_fv3sar,nlon_ens,nlat_ens,jcap_ens,&
Expand All @@ -1462,7 +1464,7 @@ module gsimod
i_en_perts_io,l_ens_in_diff_time,ensemble_path,ens_fast_read,sst_staticB,limqens, &
nsclgrp,l_timloc_opt,ngvarloc,naensloc,r_ensloccov4tim,r_ensloccov4var,r_ensloccov4scl,&
vdl_scale,vloc_varlist,&
global_spectral_filter_sd,assign_vdl_nml,parallelization_over_ensmembers
global_spectral_filter_sd,assign_vdl_nml,parallelization_over_ensmembers,l_mgbf_loc

! rapidrefresh_cldsurf (options for cloud analysis and surface
! enhancement for RR appilcation ):
Expand Down Expand Up @@ -1985,6 +1987,18 @@ subroutine gsimain_initialize
regional=wrf_nmm_regional.or.wrf_mass_regional.or.twodvar_regional.or.nems_nmmb_regional .or. cmaq_regional
regional=regional.or.fv3_regional.or.fv3_cmaq_regional

! Force turn off MGBF-based localization except for regional application
if(.not.regional.and.l_mgbf_loc) then
l_mgbf_loc=.false.
if(mype==0) write(6,*)'GSIMOD: for global app, l_mgbf_loc is not applicable, reset l_mgbf_loc=',l_mgbf_loc
end if

! Force turn off MGBF-based localization for lsqrtb=.true.
if(lsqrtb.and.l_mgbf_loc) then
l_mgbf_loc=.false.
if(mype==0) write(6,*)'GSIMOD: for lsqrtb=.true., l_mgbf_loc is not applicable, reset l_mgbf_loc=',l_mgbf_loc
end if

! Currently only able to have use_gfs_stratosphere=.true. for nems_nmmb_regional=.true.
use_gfs_stratosphere=use_gfs_stratosphere.and.(nems_nmmb_regional.or.wrf_nmm_regional)
if(mype==0) write(6,*) 'in gsimod: use_gfs_stratosphere,nems_nmmb_regional,wrf_nmm_regional= ', &
Expand Down
Loading

0 comments on commit 1cc934e

Please sign in to comment.