Merge pull request #451 from TingLei-daprediction/feature/gsi_fv3lam_…

…mix_ens

  GitHub Issue #450 Using global and regional ensembles simultaneously for fv3-lam GSI EnVar (add l_both_fv3sar_gfs_ens option)

src/gsi/cplr_get_fv3_regional_ensperts.f90

@@ -48,6 +48,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
      use constants, only: zero,one,half,zero_single,rd_over_cp,one_tenth
      use mpimod, only: mpi_comm_world,ierror,mype
      use hybrid_ensemble_parameters, only: n_ens,grd_ens
+     use hybrid_ensemble_parameters, only: l_both_fv3sar_gfs_ens, n_ens_gfs,n_ens_fv3sar
      use hybrid_ensemble_parameters, only: ntlevs_ens,ensemble_path
      use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d
      use gsi_bundlemod, only: gsi_bundlecreate
@@ -104,6 +105,18 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
 
      character(255) ensfilenam_str
      type(type_fv3regfilenameg)::fv3_filename 
+     integer(i_kind):: imem_start,n_fv3sar
+
+     if(n_ens/=(n_ens_gfs+n_ens_fv3sar)) then
+        write(6,*)'wrong, the sum of  n_ens_gfs and n_ens_fv3sar not equal n_ens, stop'
+        write(6,*)"n_ens, n_ens_gfs and n_ens_fv3sar are",n_ens, n_ens_gfs , n_ens_fv3sar
+        call stop2(222)
+     endif
+     if(l_both_fv3sar_gfs_ens) then
+        imem_start=n_ens_gfs+1
+     else
+        imem_start=1
+     endif
 
      call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)
      ! Allocate bundle to hold mean of ensemble members
@@ -205,6 +218,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
           call stop2(9991)
        endif
     enddo ! for m 
+
 
      ! print info message for dirZDA
     if(mype==0)then
@@ -233,7 +247,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
  ! INITIALIZE ENSEMBLE MEAN ACCUMULATORS
         en_bar(m)%values=zero
 
-        do n=1,n_ens
+        do n=imem_start,n_ens
            en_perts(n,m)%valuesr4 = zero
         enddo
 
@@ -242,8 +256,9 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
         kapr=one/rd_over_cp
  !
  ! LOOP OVER ENSEMBLE MEMBERS 
-        do n=1,n_ens
-           write(ensfilenam_str,22) trim(adjustl(ensemble_path)),ens_fhrlevs(m),n
+        do n_fv3sar=1,n_ens_fv3sar
+           n=n_ens_gfs+n_fv3sar
+           write(ensfilenam_str,22) trim(adjustl(ensemble_path)),ens_fhrlevs(m),n_fv3sar
 22  format(a,'fv3SAR',i2.2,'_ens_mem',i3.3)
  ! DEFINE INPUT FILE NAME
            fv3_filename%grid_spec=trim(ensfilenam_str)//'-fv3_grid_spec' !exmaple thinktobe
@@ -493,7 +508,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
         enddo 
  !
  ! CALCULATE ENSEMBLE MEAN
-        bar_norm = one/float(n_ens)
+        bar_norm = one/float(n_ens_fv3sar)
         en_bar(m)%values=en_bar(m)%values*bar_norm
 
  ! Copy pbar to module array.  ps_bar may be needed for vertical localization
@@ -521,9 +536,9 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
  !
  !
  ! CONVERT ENSEMBLE MEMBERS TO ENSEMBLE PERTURBATIONS
-        sig_norm=sqrt(one/max(one,n_ens-one))
+        sig_norm=sqrt(one/max(one,n_ens_fv3sar-one))
 
-        do n=1,n_ens
+        do n=imem_start,n_ens
            do i=1,nelen
               en_perts(n,m)%valuesr4(i)=(en_perts(n,m)%valuesr4(i)-en_bar(m)%values(i))*sig_norm
            end do

src/gsi/get_gefs_for_regional.f90

@@ -41,10 +41,9 @@ subroutine get_gefs_for_regional
                      fv3_regional
   use hybrid_ensemble_parameters, only: region_lat_ens,region_lon_ens
   use hybrid_ensemble_parameters, only: en_perts,ps_bar,nelen
-  use hybrid_ensemble_parameters, only: n_ens,grd_ens,grd_a1,grd_e1,p_e2a,uv_hyb_ens,dual_res
+  use hybrid_ensemble_parameters, only: n_ens_gfs,grd_ens,grd_a1,grd_e1,p_e2a,uv_hyb_ens,dual_res
   use hybrid_ensemble_parameters, only: full_ensemble,q_hyb_ens,l_ens_in_diff_time,write_ens_sprd
   use hybrid_ensemble_parameters, only: ntlevs_ens,ensemble_path,jcap_ens
- !use hybrid_ensemble_parameters, only: add_bias_perturbation
   use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d
   use gsi_bundlemod, only: gsi_bundlecreate
   use gsi_bundlemod, only: gsi_bundle
@@ -229,14 +228,20 @@ subroutine get_gefs_for_regional
   do n=1,200
      read(10,'(a)',err=20,end=40)filename 
   enddo
-40 n_ens=n-1
+40 n_ens_temp=n-1
+  if(n_ens_gfs/=n_ens_temp) then
+     n_ens_gfs=n_ens_temp
+     if(mype == 0) then
+         write(6,*)'the n_ens_gfs is adjusted to the actual number of ensemble members ',n_ens_temp
+     endif
+  endif
 
 !    set n_ens_temp depending on if we want to add bias perturbation to the ensemble
 
   if(add_bias_perturbation) then
-     n_ens_temp=n_ens+1
+     n_ens_temp=n_ens_gfs+1
   else
-     n_ens_temp=n_ens
+     n_ens_temp=n_ens_gfs
   end if
 
   rewind (10) 
@@ -585,13 +590,13 @@ subroutine get_gefs_for_regional
                     grd_gfs%nlat,sp_gfs%rlats,grd_gfs%nlon,sp_gfs%rlons,nord_g2r,p_g2r)
 
 !  allocate mix ensemble space--horizontal on regional domain, vertical still gefs 
-  allocate(st_eg(grd_mix%lat2,grd_mix%lon2,grd_mix%nsig,n_ens))
-  allocate(vp_eg(grd_mix%lat2,grd_mix%lon2,grd_mix%nsig,n_ens))
-  allocate( t_eg(grd_mix%lat2,grd_mix%lon2,grd_mix%nsig,n_ens))
-  allocate(rh_eg(grd_mix%lat2,grd_mix%lon2,grd_mix%nsig,n_ens))
-  allocate(oz_eg(grd_mix%lat2,grd_mix%lon2,grd_mix%nsig,n_ens))
-  allocate(cw_eg(grd_mix%lat2,grd_mix%lon2,grd_mix%nsig,n_ens))
-  allocate( p_eg_nmmb(grd_mix%lat2,grd_mix%lon2,n_ens))
+  allocate(st_eg(grd_mix%lat2,grd_mix%lon2,grd_mix%nsig,n_ens_gfs))
+  allocate(vp_eg(grd_mix%lat2,grd_mix%lon2,grd_mix%nsig,n_ens_gfs))
+  allocate( t_eg(grd_mix%lat2,grd_mix%lon2,grd_mix%nsig,n_ens_gfs))
+  allocate(rh_eg(grd_mix%lat2,grd_mix%lon2,grd_mix%nsig,n_ens_gfs))
+  allocate(oz_eg(grd_mix%lat2,grd_mix%lon2,grd_mix%nsig,n_ens_gfs))
+  allocate(cw_eg(grd_mix%lat2,grd_mix%lon2,grd_mix%nsig,n_ens_gfs))
+  allocate( p_eg_nmmb(grd_mix%lat2,grd_mix%lon2,n_ens_gfs))
   st_eg=zero ; vp_eg=zero ; t_eg=zero ; rh_eg=zero ; oz_eg=zero ; cw_eg=zero 
   p_eg_nmmb=zero
 
@@ -616,7 +621,7 @@ subroutine get_gefs_for_regional
 
   rewind(10)
   inithead=.true.
-  do n=1,n_ens
+  do n=1,n_ens_gfs
      read(10,'(a)',err=20,end=20)filename 
      filename=trim(ensemble_path) // trim(filename)
 !     write(filename,100) n
@@ -966,7 +971,7 @@ subroutine get_gefs_for_regional
 !   compute mean state
   stbar=zero ; vpbar=zero ; tbar=zero ; rhbar=zero ; ozbar=zero ; cwbar=zero 
   pbar_nmmb=zero
-  do n=1,n_ens
+  do n=1,n_ens_gfs
      do k=1,grd_mix%nsig
         do j=1,grd_mix%lon2
            do i=1,grd_mix%lat2
@@ -987,7 +992,7 @@ subroutine get_gefs_for_regional
   end do
 
 ! Convert to mean
-  bar_norm = one/float(n_ens)
+  bar_norm = one/float(n_ens_gfs)
   do k=1,grd_mix%nsig
      do j=1,grd_mix%lon2
         do i=1,grd_mix%lat2
@@ -1020,7 +1025,7 @@ subroutine get_gefs_for_regional
 !www  ensemble perturbation for all but the first member if full_ensemble
   if(full_ensemble)n1=2
 
-  do n=n1,n_ens
+  do n=n1,n_ens_gfs
      do k=1,grd_mix%nsig
         do j=1,grd_mix%lon2
            do i=1,grd_mix%lat2
@@ -1116,7 +1121,7 @@ subroutine get_gefs_for_regional
   allocate(rht(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
   allocate(ozt(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
   allocate(cwt(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
-  do n=1,n_ens
+  do n=1,n_ens_gfs
      do j=1,grd_ens%lon2
         do i=1,grd_ens%lat2
            do k=1,grd_mix%nsig

src/gsi/gsimod.F90

@@ -150,6 +150,7 @@ module gsimod
                          readin_localization,write_ens_sprd,eqspace_ensgrid,grid_ratio_ens,&
                          readin_beta,use_localization_grid,use_gfs_ens,q_hyb_ens,i_en_perts_io, &
                          l_ens_in_diff_time,ensemble_path,ens_fast_read,sst_staticB,limqens
+  use hybrid_ensemble_parameters,only : l_both_fv3sar_gfs_ens,n_ens_gfs,n_ens_fv3sar
   use rapidrefresh_cldsurf_mod, only: init_rapidrefresh_cldsurf, &
                             dfi_radar_latent_heat_time_period,metar_impact_radius,&
                             metar_impact_radius_lowcloud,l_gsd_terrain_match_surftobs, &
@@ -1350,7 +1351,7 @@ module gsimod
 !     sst_staticB - use only static background error covariance for SST statistic
 !              
 !                         
-  namelist/hybrid_ensemble/l_hyb_ens,uv_hyb_ens,q_hyb_ens,aniso_a_en,generate_ens,n_ens,nlon_ens,nlat_ens,jcap_ens,&
+  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,nlon_ens,nlat_ens,jcap_ens,&
                 pseudo_hybens,merge_two_grid_ensperts,regional_ensemble_option,fv3sar_bg_opt,fv3sar_ensemble_opt,full_ensemble,pwgtflg,&
                 jcap_ens_test,beta_s0,beta_e0,s_ens_h,s_ens_v,readin_localization,eqspace_ensgrid,readin_beta,&
                 grid_ratio_ens, &
@@ -1745,6 +1746,23 @@ subroutine gsimain_initialize
         c_varqc=c_varqc_new
      end if
   end if
+  if(l_both_fv3sar_gfs_ens) then
+    if(n_ens /= n_ens_gfs + n_ens_fv3sar .or. regional_ensemble_option /= 5 ) then 
+       write(6,*)'the set up for l_both_fv3sar_gfs_ens=.true. is wrong,stop'
+       call stop2(137)
+    endif
+  else
+    if (regional_ensemble_option==5) then 
+       n_ens_gfs=0
+       n_ens_fv3sar=n_ens
+    elseif (regional_ensemble_option==1) then 
+       n_ens_gfs=n_ens
+       n_ens_fv3sar=0
+    else 
+       write(6,*)'n_ens_gfs and n_ens_fv3sar won"t be used if not regional_ensemble_option==5' 
+    endif
+
+  endif
   if(ltlint) then
      if(vqc .or. njqc .or. nvqc)then
        vqc = .false.

src/gsi/hybrid_ensemble_isotropic.F90

@@ -1169,13 +1169,14 @@ subroutine load_ensemble
                                           pseudo_hybens,regional_ensemble_option,&
                                           i_en_perts_io
     use hybrid_ensemble_parameters, only: nelen,en_perts,ps_bar
+    use hybrid_ensemble_parameters, only: l_both_fv3sar_gfs_ens 
     use gsi_enscouplermod, only: gsi_enscoupler_put_gsi_ens
     use mpimod, only: mype
     use get_pseudo_ensperts_mod, only: get_pseudo_ensperts_class
     use get_wrf_mass_ensperts_mod, only: get_wrf_mass_ensperts_class
     use get_fv3_regional_ensperts_mod, only: get_fv3_regional_ensperts_class
     use get_wrf_nmm_ensperts_mod, only: get_wrf_nmm_ensperts_class
-  use hybrid_ensemble_parameters, only: region_lat_ens,region_lon_ens
+    use hybrid_ensemble_parameters, only: region_lat_ens,region_lon_ens
     use mpimod, only: mpi_comm_world
 
     implicit none
@@ -1358,6 +1359,10 @@ subroutine load_ensemble
 
                 call get_nmmb_ensperts
              case(5)
+                if (l_both_fv3sar_gfs_ens) then ! first read in gfs ensembles for regional 
+                   call get_gefs_for_regional
+                endif
+
 !     regional_ensemble_option = 5: ensembles are fv3 regional.
                 call fv3_regional_enspert%get_fv3_regional_ensperts(en_perts,nelen,ps_bar)
 

src/gsi/hybrid_ensemble_parameters.f90

@@ -258,6 +258,7 @@ module hybrid_ensemble_parameters
 ! set passed variables to public
   public :: generate_ens,n_ens,nlon_ens,nlat_ens,jcap_ens,jcap_ens_test,l_hyb_ens,&
        s_ens_h,oz_univ_static,vvlocal
+  public :: n_ens_gfs,n_ens_fv3sar
   public :: uv_hyb_ens,q_hyb_ens,s_ens_v,beta_s0,beta_e0,aniso_a_en,s_ens_hv,s_ens_vv
   public :: readin_beta,beta_s,beta_e
   public :: readin_localization
@@ -292,6 +293,7 @@ module hybrid_ensemble_parameters
   public :: region_lat_ens,region_lon_ens
   public :: region_dx_ens,region_dy_ens
   public :: ens_fast_read
+  public :: l_both_fv3sar_gfs_ens 
   public :: sst_staticB
   public :: limqens
 
@@ -311,8 +313,10 @@ module hybrid_ensemble_parameters
   logical vvlocal
   logical l_ens_in_diff_time
   logical ens_fast_read
+  logical l_both_fv3sar_gfs_ens
   integer(i_kind) i_en_perts_io
   integer(i_kind) n_ens,nlon_ens,nlat_ens,jcap_ens,jcap_ens_test
+  integer(i_kind) n_ens_gfs,n_ens_fv3sar
   real(r_kind) beta_s0,beta_e0,s_ens_h,s_ens_v,grid_ratio_ens
   type(sub2grid_info),save :: grd_ens,grd_loc,grd_sploc,grd_anl,grd_e1,grd_a1
   type(spec_vars),save :: sp_ens,sp_loc
@@ -420,7 +424,9 @@ subroutine init_hybrid_ensemble_parameters
   ensemble_path = './'       ! default for path to ensemble members
   ens_fast_read=.false.
   limqens=1.0_r_single       ! default for limiting ensemble RH (+/-)
-
+  l_both_fv3sar_gfs_ens=.false.
+  n_ens_gfs=0 
+  n_ens_fv3sar=0
 
 end subroutine init_hybrid_ensemble_parameters