From 3b9705748bd59be02428916b172535759c7f4a88 Mon Sep 17 00:00:00 2001 From: Sho Yokota Date: Tue, 26 Jul 2022 16:06:39 -0500 Subject: [PATCH] GitHub Issue NOAA-EMC/GSI#459 Add Scale/Variable/Time-Dependent Localization for EnVar --- src/gsi/class_get_fv3_regional_ensperts.f90 | 2 +- src/gsi/class_get_pseudo_ensperts.f90 | 2 +- src/gsi/class_get_wrf_mass_ensperts.f90 | 4 +- src/gsi/class_get_wrf_nmm_ensperts.f90 | 2 +- src/gsi/control_vectors.f90 | 122 +- src/gsi/cplr_get_fv3_regional_ensperts.f90 | 14 +- src/gsi/cplr_get_pseudo_ensperts.f90 | 28 +- src/gsi/cplr_get_wrf_mass_ensperts.f90 | 16 +- src/gsi/cplr_get_wrf_nmm_ensperts.f90 | 14 +- src/gsi/cplr_gfs_ensmod.f90 | 2 +- src/gsi/en_perts_io.f90 | 17 +- src/gsi/ens_spread_mod.f90 | 4 +- src/gsi/ensctl2model.f90 | 38 +- src/gsi/ensctl2model_ad.f90 | 59 +- src/gsi/ensctl2state.f90 | 4 +- src/gsi/ensctl2state_ad.f90 | 4 +- src/gsi/get_gefs_ensperts_dualres.f90 | 36 +- src/gsi/get_gefs_for_regional.f90 | 8 +- src/gsi/get_nmmb_ensperts.f90 | 8 +- src/gsi/gsimod.F90 | 50 +- src/gsi/hybrid_ensemble_isotropic.F90 | 1347 ++++++++++++------- src/gsi/hybrid_ensemble_parameters.f90 | 79 +- src/gsi/jfunc.f90 | 5 +- src/gsi/stub_get_pseudo_ensperts.f90 | 2 +- src/gsi/stub_get_wrf_mass_ensperts.f90 | 4 +- src/gsi/stub_get_wrf_nmm_ensperts.f90 | 2 +- 26 files changed, 1203 insertions(+), 670 deletions(-) diff --git a/src/gsi/class_get_fv3_regional_ensperts.f90 b/src/gsi/class_get_fv3_regional_ensperts.f90 index 6c51a4fd88..bda4a47cd3 100644 --- a/src/gsi/class_get_fv3_regional_ensperts.f90 +++ b/src/gsi/class_get_fv3_regional_ensperts.f90 @@ -31,7 +31,7 @@ subroutine get_fv3_regional_ensperts(this,en_perts,nelen,ps_bar) import abstract_get_fv3_regional_ensperts_class implicit none class(abstract_get_fv3_regional_ensperts_class),intent(inout) :: this - type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:) integer(i_kind), intent(in ):: nelen real(r_single),dimension(:,:,:),allocatable, intent(inout):: ps_bar diff --git a/src/gsi/class_get_pseudo_ensperts.f90 b/src/gsi/class_get_pseudo_ensperts.f90 index 703a62df3f..beee90f5e9 100644 --- a/src/gsi/class_get_pseudo_ensperts.f90 +++ b/src/gsi/class_get_pseudo_ensperts.f90 @@ -12,7 +12,7 @@ subroutine get_pseudo_ensperts(this,en_perts,nelen) import abstract_get_pseudo_ensperts_class implicit none class(abstract_get_pseudo_ensperts_class), intent(inout) :: this - type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:,:) integer(i_kind), intent(in ) :: nelen end subroutine get_pseudo_ensperts end interface diff --git a/src/gsi/class_get_wrf_mass_ensperts.f90 b/src/gsi/class_get_wrf_mass_ensperts.f90 index 902b14ad21..b392bcd174 100644 --- a/src/gsi/class_get_wrf_mass_ensperts.f90 +++ b/src/gsi/class_get_wrf_mass_ensperts.f90 @@ -12,7 +12,7 @@ subroutine get_wrf_mass_ensperts(this,en_perts,nelen,ps_bar) import abstract_get_wrf_mass_ensperts_class implicit none class(abstract_get_wrf_mass_ensperts_class), intent(inout) :: this - type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:) integer(i_kind), intent(in ):: nelen real(r_single),dimension(:,:,:),allocatable:: ps_bar end subroutine get_wrf_mass_ensperts @@ -25,7 +25,7 @@ subroutine ens_spread_dualres_regional(this,mype,en_perts,nelen,en_bar) implicit none class(abstract_get_wrf_mass_ensperts_class), intent(inout) :: this integer(i_kind),intent(in):: mype - type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:,:) integer(i_kind), intent(in ):: nelen type(gsi_bundle),OPTIONAL,intent(in):: en_bar end subroutine ens_spread_dualres_regional diff --git a/src/gsi/class_get_wrf_nmm_ensperts.f90 b/src/gsi/class_get_wrf_nmm_ensperts.f90 index 19170033e2..9f496c7b03 100644 --- a/src/gsi/class_get_wrf_nmm_ensperts.f90 +++ b/src/gsi/class_get_wrf_nmm_ensperts.f90 @@ -13,7 +13,7 @@ subroutine get_wrf_nmm_ensperts(this,en_perts,nelen,region_lat_ens,region_lon_en import abstract_get_wrf_nmm_ensperts_class implicit none class(abstract_get_wrf_nmm_ensperts_class),intent(inout) :: this - type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:) integer(i_kind), intent(in ):: nelen real(r_kind),allocatable, intent(inout):: region_lat_ens(:,:),region_lon_ens(:,:) real(r_single),dimension(:,:,:),allocatable, intent(inout):: ps_bar diff --git a/src/gsi/control_vectors.f90 b/src/gsi/control_vectors.f90 index a1c0a61ed9..0847257777 100644 --- a/src/gsi/control_vectors.f90 +++ b/src/gsi/control_vectors.f90 @@ -144,7 +144,7 @@ module control_vectors type(GSI_Bundle), pointer :: step(:) type(GSI_Bundle), pointer :: motley(:) type(GSI_Grid) :: grid_aens - type(GSI_Bundle), pointer :: aens(:,:) + type(GSI_Bundle), pointer :: aens(:,:,:) real(r_kind), pointer :: predr(:) => NULL() real(r_kind), pointer :: predp(:) => NULL() real(r_kind), pointer :: predt(:) => NULL() @@ -296,6 +296,7 @@ subroutine init_anacv ! language: f90 ! machine: ibm rs/6000 sp ! +use hybrid_ensemble_parameters,only:idaen3d,idaen2d implicit none !character(len=*),parameter:: rcname='anavinfo.txt' character(len=*),parameter:: rcname='anavinfo' ! filename should have extension @@ -345,6 +346,7 @@ subroutine init_anacv allocate(cvarsmd(mvars)) allocate(atsfc_sdv(mvars)) allocate(an_amp0(nvars)) +allocate(idaen3d(nc3d),idaen2d(nc2d)) allocate(incvars_to_zero(nvars)) allocate(incvars_zero_strat(nvars)) incvars_to_zero(:) = 'NONE' @@ -371,12 +373,22 @@ subroutine init_anacv cvars2d(nc2d)=trim(adjustl(var)) nrf2_loc(nc2d)=ii ! rid of soon as2d(nc2d)=aas + if(itracer>10) then + idaen2d(nc2d)=2 + else + idaen2d(nc2d)=1 + endif else nc3d=nc3d+1 cvars3d(nc3d)=trim(adjustl(var)) nrf3_loc(nc3d)=ii ! rid of soon nrf_3d(ii)=.true. as3d(nc3d)=aas + if(itracer>10) then + idaen3d(nc3d)=2 + else + idaen3d(nc3d)=1 + endif endif endif nrf_var(ii)=trim(adjustl(var)) @@ -446,9 +458,11 @@ subroutine allocate_cv(ycv) !$$$ end documentation block use hybrid_ensemble_parameters, only: grd_ens + use hybrid_ensemble_parameters, only: naensgrp implicit none type(control_vector), intent( out) :: ycv integer(i_kind) :: ii,jj,nn,ndim,ierror,n_step,n_aens + integer(i_kind) :: ig character(len=256)::bname character(len=max_varname_length)::ltmp(1) type(gsi_grid) :: grid_motley @@ -481,7 +495,7 @@ subroutine allocate_cv(ycv) ! If so, define grid of ensemble control vector n_aens=0 if (l_hyb_ens) then - ALLOCATE(ycv%aens(nsubwin,n_ens)) + ALLOCATE(ycv%aens(nsubwin,naensgrp,n_ens)) call GSI_GridCreate(ycv%grid_aens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig) if (lsqrtb) then n_aens=nval_lenz_en @@ -524,17 +538,19 @@ subroutine allocate_cv(ycv) if (l_hyb_ens) then ltmp(1)='a_en' - do nn=1,n_ens - ycv%aens(jj,nn)%values => ycv%values(ii+1:ii+n_aens) - write(bname,'(a,i3.3,a,i4.4)') 'Ensemble Control Bundle subwin-',jj,' and member-',nn - call GSI_BundleSet(ycv%aens(jj,nn),ycv%grid_aens,bname,ierror,names3d=ltmp,bundle_kind=r_kind) - if (ierror/=0) then - write(6,*)'allocate_cv: error alloc(ensemble bundle)' - call stop2(109) - endif - ndim=ndim+ycv%aens(jj,nn)%ndim - - ii=ii+n_aens + do ig=1,naensgrp + do nn=1,n_ens + ycv%aens(jj,ig,nn)%values => ycv%values(ii+1:ii+n_aens) + write(bname,'(a,i3.3,a,i4.4)') 'Ensemble Control Bundle subwin-',jj,' and member-',nn + call GSI_BundleSet(ycv%aens(jj,ig,nn),ycv%grid_aens,bname,ierror,names3d=ltmp,bundle_kind=r_kind) + if (ierror/=0) then + write(6,*)'allocate_cv: error alloc(ensemble bundle)' + call stop2(109) + endif + ndim=ndim+ycv%aens(jj,ig,nn)%ndim + + ii=ii+n_aens + enddo enddo endif @@ -620,15 +636,19 @@ subroutine deallocate_cv(ycv) ! !$$$ end documentation block + use hybrid_ensemble_parameters, only: naensgrp implicit none type(control_vector), intent(inout) :: ycv integer(i_kind) :: ii,nn,ierror + integer(i_kind) :: ig if (ycv%lallocated) then do ii=1,nsubwin if (l_hyb_ens) then - do nn=n_ens,1,-1 - call GSI_BundleUnset(ycv%aens(ii,nn),ierror) + do ig=1,naensgrp + do nn=n_ens,1,-1 + call GSI_BundleUnset(ycv%aens(ii,ig,nn),ierror) + enddo enddo endif ! if (beta_s0>tiny_r_kind) then @@ -845,9 +865,11 @@ real(r_quad) function qdot_prod_sub(xcv,ycv) ! !$$$ end documentation block + use hybrid_ensemble_parameters, only: naensgrp implicit none type(control_vector), intent(in ) :: xcv, ycv integer(i_kind) :: ii,nn,m3d,m2d,i,j,itot + integer(i_kind) :: ig,nigtmp real(r_quad),allocatable,dimension(:) :: partsum qdot_prod_sub=zero_quad @@ -858,9 +880,11 @@ real(r_quad) function qdot_prod_sub(xcv,ycv) qdot_prod_sub=qdot_prod_sub+qdot_product( xcv%step(ii)%values(:) ,ycv%step(ii)%values(:) ) end do if(l_hyb_ens) then - do nn=1,n_ens - do ii=1,nsubwin - qdot_prod_sub=qdot_prod_sub+qdot_product( xcv%aens(ii,nn)%values(:) ,ycv%aens(ii,nn)%values(:) ) + do ig=1,naensgrp + do nn=1,n_ens + do ii=1,nsubwin + qdot_prod_sub=qdot_prod_sub+qdot_product( xcv%aens(ii,ig,nn)%values(:) ,ycv%aens(ii,ig,nn)%values(:) ) + end do end do end do endif @@ -869,7 +893,7 @@ real(r_quad) function qdot_prod_sub(xcv,ycv) m3d=xcv%step(ii)%n3d m2d=xcv%step(ii)%n2d itot=max(m3d,0)+max(m2d,0) - if(l_hyb_ens)itot=itot+n_ens + if(l_hyb_ens)itot=itot+naensgrp*n_ens allocate(partsum(itot)) !$omp parallel do schedule(dynamic,1) private(i) do i = 1,m3d @@ -880,9 +904,12 @@ real(r_quad) function qdot_prod_sub(xcv,ycv) partsum(m3d+i) = dplevs(xcv%step(ii)%r2(i)%q,ycv%step(ii)%r2(i)%q,ihalo=1) enddo if(l_hyb_ens) then + do ig=1,naensgrp + nigtmp=n_ens*(ig-1) !$omp parallel do schedule(dynamic,1) private(i) - do i = 1,n_ens - partsum(m3d+m2d+i) = dplevs(xcv%aens(ii,i)%r3(1)%q,ycv%aens(ii,i)%r3(1)%q,ihalo=1) + do i = 1,n_ens + partsum(m3d+m2d+nigtmp+i) = dplevs(xcv%aens(ii,ig,i)%r3(1)%q,ycv%aens(ii,ig,i)%r3(1)%q,ihalo=1) + end do end do end if do i=1,itot @@ -933,6 +960,7 @@ subroutine qdot_prod_vars_eb(xcv,ycv,prods,eb) ! !$$$ end documentation block + use hybrid_ensemble_parameters, only: naensgrp implicit none type(control_vector), intent(in ) :: xcv, ycv character(len=*) , intent(in ) :: eb @@ -941,6 +969,8 @@ subroutine qdot_prod_vars_eb(xcv,ycv,prods,eb) real(r_quad) :: zz(nsubwin) integer(i_kind) :: ii,i,nn,m3d,m2d real(r_quad),allocatable,dimension(:) :: partsum + integer(i_kind) :: ig + integer(i_kind) ::ngtmp,nn0 prods(:)=zero_quad zz(:)=zero_quad @@ -953,9 +983,11 @@ subroutine qdot_prod_vars_eb(xcv,ycv,prods,eb) end do endif if(trim(eb) == 'cost_e') then - do nn=1,n_ens - do ii=1,nsubwin - zz(ii)=zz(ii)+qdot_product( xcv%aens(ii,nn)%values(:) ,ycv%aens(ii,nn)%values(:) ) + do ig=1,naensgrp + do nn=1,n_ens + do ii=1,nsubwin + zz(ii)=zz(ii)+qdot_product( xcv%aens(ii,ig,nn)%values(:) ,ycv%aens(ii,ig,nn)%values(:) ) + end do end do end do endif @@ -981,21 +1013,25 @@ subroutine qdot_prod_vars_eb(xcv,ycv,prods,eb) end if if(trim(eb) == 'cost_e') then do ii=1,nsubwin ! RTod: somebody could work in opt/zing this ... - allocate(partsum(n_ens)) + allocate(partsum(n_ens*naensgrp)) + do ig=1,naensgrp + ngtmp=(ig-1)*n_ens !$omp parallel do schedule(dynamic,1) private(nn,m3d,m2d) - do nn=1,n_ens - partsum(nn) = zero_quad - m3d=xcv%aens(ii,nn)%n3d - do i = 1,m3d - partsum(nn)= partsum(nn) + dplevs(xcv%aens(ii,nn)%r3(i)%q,ycv%aens(ii,nn)%r3(i)%q,ihalo=1) + do nn=1,n_ens + nn0=nn+ngtmp + partsum(nn0) = zero_quad + m3d=xcv%aens(ii,ig,nn)%n3d + do i = 1,m3d + partsum(nn0)= partsum(nn0) + dplevs(xcv%aens(ii,ig,nn)%r3(i)%q,ycv%aens(ii,ig,nn)%r3(i)%q,ihalo=1) + enddo + m2d=xcv%aens(ii,ig,nn)%n2d + do i = 1,m2d + partsum(nn0)= partsum(nn0) + dplevs(xcv%aens(ii,ig,nn)%r2(i)%q,ycv%aens(ii,ig,nn)%r2(i)%q,ihalo=1) + enddo enddo - m2d=xcv%aens(ii,nn)%n2d - do i = 1,m2d - partsum(nn)= partsum(nn) + dplevs(xcv%aens(ii,nn)%r2(i)%q,ycv%aens(ii,nn)%r2(i)%q,ihalo=1) - enddo - enddo - do nn=1,n_ens - zz(ii)=zz(ii)+partsum(nn) + end do + do nn=1,n_ens*naensgrp + zz(ii)=zz(ii)+partsum(nn) end do deallocate(partsum) end do @@ -1348,6 +1384,7 @@ subroutine random_cv(ycv,kseed) ! !$$$ end documentation block +use hybrid_ensemble_parameters, only: naensgrp implicit none type(control_vector) , intent(inout) :: ycv integer(i_kind), optional, intent(in ) :: kseed @@ -1355,6 +1392,7 @@ subroutine random_cv(ycv,kseed) integer(i_kind):: ii,jj,nn,iseed integer, allocatable :: nseed(:) ! Intentionaly default integer real(r_kind), allocatable :: zz(:) +integer(i_kind) :: ig iseed=iadatebgn if (present(kseed)) iseed=iseed+kseed @@ -1381,10 +1419,12 @@ subroutine random_cv(ycv,kseed) if (nval_lenz_en>0) then allocate(zz(nval_lenz_en)) do nn=1,n_ens - do jj=1,nsubwin - call random_number(zz) - do ii=1,nval_lenz_en - ycv%aens(jj,nn)%values(ii) = two*zz(ii)-one + do ig=1,naensgrp + do jj=1,nsubwin + call random_number(zz) + do ii=1,nval_lenz_en + ycv%aens(jj,ig,nn)%values(ii) = two*zz(ii)-one + enddo enddo enddo enddo diff --git a/src/gsi/cplr_get_fv3_regional_ensperts.f90 b/src/gsi/cplr_get_fv3_regional_ensperts.f90 index d64f266026..f99ca39790 100644 --- a/src/gsi/cplr_get_fv3_regional_ensperts.f90 +++ b/src/gsi/cplr_get_fv3_regional_ensperts.f90 @@ -73,7 +73,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) implicit none class(get_fv3_regional_ensperts_class), intent(inout) :: this - type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:) integer(i_kind), intent(in ):: nelen real(r_single),dimension(:,:,:),allocatable,intent(inout):: ps_bar @@ -248,7 +248,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) en_bar(m)%values=zero do n=imem_start,n_ens - en_perts(n,m)%valuesr4 = zero + en_perts(n,1,m)%valuesr4 = zero enddo mm1=mype+1 @@ -311,7 +311,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) ! SAVE ENSEMBLE MEMBER DATA IN COLUMN VECTOR do ic3=1,nc3d - call gsi_bundlegetpointer(en_perts(n,m),trim(cvars3d(ic3)),w3,istatus) + call gsi_bundlegetpointer(en_perts(n,1,m),trim(cvars3d(ic3)),w3,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n call stop2(9992) @@ -471,7 +471,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) do ic2=1,nc2d - call gsi_bundlegetpointer(en_perts(n,m),trim(cvars2d(ic2)),w2,istatus) + call gsi_bundlegetpointer(en_perts(n,1,m),trim(cvars2d(ic2)),w2,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n call stop2(9994) @@ -540,7 +540,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) 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 + en_perts(n,1,m)%valuesr4(i)=(en_perts(n,1,m)%valuesr4(i)-en_bar(m)%values(i))*sig_norm end do end do @@ -880,7 +880,7 @@ subroutine ens_spread_dualres_regional_fv3_regional(this,mype,en_perts,nelen) class(get_fv3_regional_ensperts_class), intent(inout) :: this integer(i_kind),intent(in):: mype - type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:,:) integer(i_kind), intent(in ):: nelen type(gsi_bundle):: sube,suba @@ -922,7 +922,7 @@ subroutine ens_spread_dualres_regional_fv3_regional(this,mype,en_perts,nelen) do n=1,n_ens do i=1,nelen sube%values(i)=sube%values(i) & - +(en_perts(n,1)%valuesr4(i))*(en_perts(n,1)%valuesr4(i)) + +(en_perts(n,1,1)%valuesr4(i))*(en_perts(n,1,1)%valuesr4(i)) end do end do diff --git a/src/gsi/cplr_get_pseudo_ensperts.f90 b/src/gsi/cplr_get_pseudo_ensperts.f90 index 8c59cd0be8..9d12165409 100644 --- a/src/gsi/cplr_get_pseudo_ensperts.f90 +++ b/src/gsi/cplr_get_pseudo_ensperts.f90 @@ -57,7 +57,7 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen) implicit none ! Declare passed variables class(get_pseudo_ensperts_class), intent(inout) :: this - type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:,:) integer(i_kind), intent(in ) :: nelen @@ -600,12 +600,12 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen) ! NOTE: because library perturbation bundle structure is same as ensemble perturbation, ! use same ipc3d and ipc2d indices for lib_perts and en_perts bundles. - call gsi_bundlegetpointer (en_perts(1,1),cvars3d,ipc3d,istatus) + call gsi_bundlegetpointer (en_perts(1,1,1),cvars3d,ipc3d,istatus) if(istatus/=0) then write(6,*) 'get_pseudo_enperts: cannot find 3d pointers' call stop2(999) endif - call gsi_bundlegetpointer (en_perts(1,1),cvars2d,ipc2d,istatus) + call gsi_bundlegetpointer (en_perts(1,1,1),cvars2d,ipc2d,istatus) if(istatus/=0) then write(6,*) 'get_pseudo_enperts: cannot find 2d pointers' call stop2(999) @@ -618,9 +618,9 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen) do k=1,grd_ens%nsig do j=1,grd_ens%lon2 do i=1,grd_ens%lat2 - en_perts(n,1)%r3(ipc3d(ic3))%qr4(i,j,k) = wgt(i,j) * & + en_perts(n,1,1)%r3(ipc3d(ic3))%qr4(i,j,k) = wgt(i,j) * & lib_perts(n)%r3(ipc3d(ic3))%qr4(i,j,k) + & - (one-wgt(i,j))*en_perts(n,1)%r3(ipc3d(ic3))%qr4(i,j,k) + (one-wgt(i,j))*en_perts(n,1,1)%r3(ipc3d(ic3))%qr4(i,j,k) end do end do end do @@ -631,9 +631,9 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen) do j=1,grd_ens%lon2 do i=1,grd_ens%lat2 - en_perts(n,1)%r2(ipc2d(ic2))%qr4(i,j) = wgt(i,j) * & + en_perts(n,1,1)%r2(ipc2d(ic2))%qr4(i,j) = wgt(i,j) * & lib_perts(n)%r2(ipc2d(ic2))%qr4(i,j) + & - (one-wgt(i,j))*en_perts(n,1)%r2(ipc2d(ic2))%qr4(i,j) + (one-wgt(i,j))*en_perts(n,1,1)%r2(ipc2d(ic2))%qr4(i,j) end do end do @@ -667,12 +667,12 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen) write(filename,"('enpert',i3.3)") n - call gsi_bundlegetpointer(en_perts(n,1),cvars3d,ipc3d,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),cvars3d,ipc3d,istatus) if(istatus/=0) then write(6,*) 'mtong: cannot find 3d pointers' call stop2(999) endif - call gsi_bundlegetpointer(en_perts(n,1),cvars2d,ipc2d,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),cvars2d,ipc2d,istatus) if(istatus/=0) then write(6,*) 'mtong: cannot find 2d pointers' call stop2(999) @@ -686,7 +686,7 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen) do k=1,grd_ens%nsig do j=1,grd_ens%lon2 do i=1,grd_ens%lat2 - u(i,j,k) = en_perts(n,1)%r3(ipc3d(ic3))%qr4(i,j,k) + u(i,j,k) = en_perts(n,1,1)%r3(ipc3d(ic3))%qr4(i,j,k) end do end do end do @@ -696,7 +696,7 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen) do k=1,grd_ens%nsig do j=1,grd_ens%lon2 do i=1,grd_ens%lat2 - v(i,j,k) = en_perts(n,1)%r3(ipc3d(ic3))%qr4(i,j,k) + v(i,j,k) = en_perts(n,1,1)%r3(ipc3d(ic3))%qr4(i,j,k) end do end do end do @@ -706,7 +706,7 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen) do k=1,grd_ens%nsig do j=1,grd_ens%lon2 do i=1,grd_ens%lat2 - tv(i,j,k) = en_perts(n,1)%r3(ipc3d(ic3))%qr4(i,j,k) + tv(i,j,k) = en_perts(n,1,1)%r3(ipc3d(ic3))%qr4(i,j,k) end do end do end do @@ -716,7 +716,7 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen) do k=1,grd_ens%nsig do j=1,grd_ens%lon2 do i=1,grd_ens%lat2 - rh(i,j,k) = en_perts(n,1)%r3(ipc3d(ic3))%qr4(i,j,k) + rh(i,j,k) = en_perts(n,1,1)%r3(ipc3d(ic3))%qr4(i,j,k) end do end do end do @@ -731,7 +731,7 @@ subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen) do j=1,grd_ens%lon2 do i=1,grd_ens%lat2 - ps(i,j) = en_perts(n,1)%r2(ipc2d(ic2))%qr4(i,j) + ps(i,j) = en_perts(n,1,1)%r2(ipc2d(ic2))%qr4(i,j) end do end do diff --git a/src/gsi/cplr_get_wrf_mass_ensperts.f90 b/src/gsi/cplr_get_wrf_mass_ensperts.f90 index eb94c6c0f3..775c9c4195 100644 --- a/src/gsi/cplr_get_wrf_mass_ensperts.f90 +++ b/src/gsi/cplr_get_wrf_mass_ensperts.f90 @@ -63,7 +63,7 @@ subroutine get_wrf_mass_ensperts_wrf(this,en_perts,nelen,ps_bar) implicit none class(get_wrf_mass_ensperts_class), intent(inout) :: this - type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:) integer(i_kind), intent(in ):: nelen real(r_single),dimension(:,:,:),allocatable:: ps_bar @@ -148,7 +148,7 @@ subroutine get_wrf_mass_ensperts_wrf(this,en_perts,nelen,ps_bar) en_bar%values=zero do n=1,n_ens - en_perts(n,it)%valuesr4 = zero + en_perts(n,1,it)%valuesr4 = zero enddo if ( .not. l_use_dbz_directDA ) then !direct reflectivity DA option does not employ this @@ -253,7 +253,7 @@ subroutine get_wrf_mass_ensperts_wrf(this,en_perts,nelen,ps_bar) ! SAVE ENSEMBLE MEMBER DATA IN COLUMN VECTOR member_data_loop: do ic3=1,nc3d - call gsi_bundlegetpointer(en_perts(n,it),trim(cvars3d(ic3)),w3,istatus) + call gsi_bundlegetpointer(en_perts(n,1,it),trim(cvars3d(ic3)),w3,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n call stop2(999) @@ -442,7 +442,7 @@ subroutine get_wrf_mass_ensperts_wrf(this,en_perts,nelen,ps_bar) member_mass_loop: do ic2=1,nc2d - call gsi_bundlegetpointer(en_perts(n,it),trim(cvars2d(ic2)),w2,istatus) + call gsi_bundlegetpointer(en_perts(n,1,it),trim(cvars2d(ic2)),w2,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n call stop2(999) @@ -515,7 +515,7 @@ subroutine get_wrf_mass_ensperts_wrf(this,en_perts,nelen,ps_bar) do n=1,n_ens do i=1,nelen - en_perts(n,it)%valuesr4(i)=(en_perts(n,it)%valuesr4(i)-en_bar%values(i))*sig_norm + en_perts(n,1,it)%valuesr4(i)=(en_perts(n,1,it)%valuesr4(i)-en_bar%values(i))*sig_norm end do end do @@ -2147,7 +2147,7 @@ subroutine ens_spread_dualres_regional_wrf(this,mype,en_perts,nelen,en_bar) class(get_wrf_mass_ensperts_class), intent(inout) :: this type(gsi_bundle),OPTIONAL,intent(in):: en_bar integer(i_kind),intent(in):: mype - type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:,:) integer(i_kind), intent(in ):: nelen type(gsi_bundle):: sube,suba @@ -2201,7 +2201,7 @@ subroutine ens_spread_dualres_regional_wrf(this,mype,en_perts,nelen,en_bar) do n=1,n_ens do i=1,nelen sube%values(i)=sube%values(i) & - +en_perts(n,1)%valuesr4(i)*en_perts(n,1)%valuesr4(i) + +en_perts(n,1,1)%valuesr4(i)*en_perts(n,1,1)%valuesr4(i) end do end do @@ -2212,7 +2212,7 @@ subroutine ens_spread_dualres_regional_wrf(this,mype,en_perts,nelen,en_bar) do n=1,n_ens do i=1,nelen sube%values(i)=sube%values(i) & - +(en_perts(n,1)%valuesr4(i)-en_bar%values(i))*(en_perts(n,1)%valuesr4(i)-en_bar%values(i)) + +(en_perts(n,1,1)%valuesr4(i)-en_bar%values(i))*(en_perts(n,1,1)%valuesr4(i)-en_bar%values(i)) end do end do diff --git a/src/gsi/cplr_get_wrf_nmm_ensperts.f90 b/src/gsi/cplr_get_wrf_nmm_ensperts.f90 index 9de32ea384..3e30822639 100644 --- a/src/gsi/cplr_get_wrf_nmm_ensperts.f90 +++ b/src/gsi/cplr_get_wrf_nmm_ensperts.f90 @@ -72,7 +72,7 @@ subroutine get_wrf_nmm_ensperts_wrf(this,en_perts,nelen,region_lat_ens,region_lo implicit none class(get_wrf_nmm_ensperts_class), intent(inout) :: this - type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:) integer(i_kind), intent(in ):: nelen real(r_kind),allocatable, intent(inout):: region_lat_ens(:,:),region_lon_ens(:,:) real(r_single),dimension(:,:,:),allocatable, intent(inout):: ps_bar @@ -689,7 +689,7 @@ subroutine get_wrf_nmm_ensperts_wrf(this,en_perts,nelen,region_lat_ens,region_lo ! SAVE ENSEMBLE MEMBER DATA IN COLUMN VECTOR do ic3=1,nc3d - call gsi_bundlegetpointer(en_perts(n,1),trim(cvars3d(ic3)),w3,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars3d(ic3)),w3,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n call stop2(999) @@ -791,7 +791,7 @@ subroutine get_wrf_nmm_ensperts_wrf(this,en_perts,nelen,region_lat_ens,region_lo do ic2=1,nc2d - call gsi_bundlegetpointer(en_perts(n,1),trim(cvars2d(ic2)),w2,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars2d(ic2)),w2,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n call stop2(999) @@ -828,7 +828,7 @@ subroutine get_wrf_nmm_ensperts_wrf(this,en_perts,nelen,region_lat_ens,region_lo if (use_gfs_stratosphere)then do n=1,n_ens do ic3=1,nc3d - call gsi_bundlegetpointer(en_perts(n,1),trim(cvars3d(ic3)),w3,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars3d(ic3)),w3,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' & for ensemble member ',n @@ -893,7 +893,7 @@ subroutine get_wrf_nmm_ensperts_wrf(this,en_perts,nelen,region_lat_ens,region_lo do n=1,n_ens do i=1,nelen - en_perts(n,1)%valuesr4(i)=(en_perts(n,1)%valuesr4(i)-en_bar%values(i))*sig_norm + en_perts(n,1,1)%valuesr4(i)=(en_perts(n,1,1)%valuesr4(i)-en_bar%values(i))*sig_norm end do end do @@ -2917,7 +2917,7 @@ subroutine ens_member_mean_dualres_regional(this,en_bar,mype,en_perts,nelen) class(get_wrf_nmm_ensperts_class), intent(inout) :: this type(gsi_bundle),intent(in):: en_bar integer(i_kind),intent(in):: mype - type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:) integer(i_kind), intent(in ):: nelen type(gsi_bundle):: sube,suba @@ -2961,7 +2961,7 @@ subroutine ens_member_mean_dualres_regional(this,en_bar,mype,en_perts,nelen) do i=1,nelen if(n <= n_ens)then - sube%values(i)=en_perts(n,1)%valuesr4(i)*sig_norm_inv+en_bar%values(i) + sube%values(i)=en_perts(n,1,1)%valuesr4(i)*sig_norm_inv+en_bar%values(i) else sube%values(i)=en_bar%values(i) end if diff --git a/src/gsi/cplr_gfs_ensmod.f90 b/src/gsi/cplr_gfs_ensmod.f90 index 582574b52f..5c9934919b 100644 --- a/src/gsi/cplr_gfs_ensmod.f90 +++ b/src/gsi/cplr_gfs_ensmod.f90 @@ -440,7 +440,7 @@ subroutine move2bundle_(grd3d,en_loc3,atm_bundle,m_cvars2d,m_cvars3d,iret) ! if(trim(cvars2d(m))=='sst') sst=en_loc3(:,:,m_cvars2d(m)) !no sst for now enddo - km = en_perts(1,1)%grid%km + km = en_perts(1,1,1)%grid%km !$omp parallel do schedule(dynamic,1) private(m) do m=1,nc3d if(trim(cvars3d(m))=='sf')then diff --git a/src/gsi/en_perts_io.f90 b/src/gsi/en_perts_io.f90 index 1a9dc58a50..f873a8520f 100644 --- a/src/gsi/en_perts_io.f90 +++ b/src/gsi/en_perts_io.f90 @@ -124,7 +124,7 @@ subroutine en_perts_get_from_save_fulldomain do ic3=1,nc3d - call gsi_bundlegetpointer(en_perts(n,1),trim(cvars3d(ic3)),w3,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars3d(ic3)),w3,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n call stop2(999) @@ -142,7 +142,7 @@ subroutine en_perts_get_from_save_fulldomain do ic2=1,nc2d - call gsi_bundlegetpointer(en_perts(n,1),trim(cvars2d(ic2)),w2,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars2d(ic2)),w2,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n call stop2(999) @@ -221,7 +221,7 @@ subroutine en_perts_get_from_save ! do ic3=1,nc3d - call gsi_bundlegetpointer(en_perts(n,1),trim(cvars3d(ic3)),w3,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars3d(ic3)),w3,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n call stop2(999) @@ -238,7 +238,7 @@ subroutine en_perts_get_from_save do ic2=1,nc2d - call gsi_bundlegetpointer(en_perts(n,1),trim(cvars2d(ic2)),w2,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars2d(ic2)),w2,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n call stop2(999) @@ -289,6 +289,7 @@ subroutine en_perts_save use kinds, only: r_kind,i_kind,r_single use gsi_bundlemod, only: GSI_BundleGetPointer use mpeu_util, only: die + use hybrid_ensemble_parameters, only: nsclgrp implicit none real(r_single),pointer,dimension(:,:,:):: w3 @@ -300,6 +301,10 @@ subroutine en_perts_save integer(i_kind) ic3,ic2 integer(i_kind) iunit + if(nsclgrp>1) then + write(6,*)"nsclgrp >1 is not considerred in this part, stop" + stop + endif iunit=20 write(filename,'(a,I4.4)') 'saved_en_perts.pe',mype @@ -312,7 +317,7 @@ subroutine en_perts_save ! do ic3=1,nc3d - call gsi_bundlegetpointer(en_perts(n,1),trim(cvars3d(ic3)),w3,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars3d(ic3)),w3,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n call stop2(999) @@ -324,7 +329,7 @@ subroutine en_perts_save end do do ic2=1,nc2d - call gsi_bundlegetpointer(en_perts(n,1),trim(cvars2d(ic2)),w2,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars2d(ic2)),w2,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n call stop2(999) diff --git a/src/gsi/ens_spread_mod.f90 b/src/gsi/ens_spread_mod.f90 index afab9737ac..da70d31ab0 100644 --- a/src/gsi/ens_spread_mod.f90 +++ b/src/gsi/ens_spread_mod.f90 @@ -90,7 +90,7 @@ subroutine ens_spread_dualres_regional(en_bar) do n=1,n_ens do i=1,nelen sube%values(i)=sube%values(i) & - +en_perts(n,1)%valuesr4(i)*en_perts(n,1)%valuesr4(i) + +en_perts(n,1,1)%valuesr4(i)*en_perts(n,1,1)%valuesr4(i) end do end do @@ -101,7 +101,7 @@ subroutine ens_spread_dualres_regional(en_bar) do n=1,n_ens do i=1,nelen sube%values(i)=sube%values(i) & - +(en_perts(n,1)%valuesr4(i)-en_bar%values(i))*(en_perts(n,1)%valuesr4(i)-en_bar%values(i)) + +(en_perts(n,1,1)%valuesr4(i)-en_bar%values(i))*(en_perts(n,1,1)%valuesr4(i)-en_bar%values(i)) end do end do diff --git a/src/gsi/ensctl2model.f90 b/src/gsi/ensctl2model.f90 index 525bf33af6..12e1fe374e 100644 --- a/src/gsi/ensctl2model.f90 +++ b/src/gsi/ensctl2model.f90 @@ -43,6 +43,7 @@ subroutine ensctl2model(xhat,mval,eval) use gsi_metguess_mod, only: gsi_metguess_get use mod_strong, only: tlnmc_option use timermod, only: timer_ini,timer_fnl +use hybrid_ensemble_parameters,only: naensgrp implicit none ! Declare passed variables @@ -58,7 +59,7 @@ subroutine ensctl2model(xhat,mval,eval) integer(i_kind), parameter :: ncvars = 5 integer(i_kind) :: icps(ncvars) type(gsi_bundle):: wbundle_c ! work bundle -type(gsi_bundle),allocatable :: ebundle(:) +type(gsi_bundle),allocatable :: ebundle(:,:) character(len=3), parameter :: mycvars(ncvars) = (/ & ! vars from CV needed here 'sf ', 'vp ', 'ps ', 't ', & 'q '/) @@ -77,6 +78,7 @@ subroutine ensctl2model(xhat,mval,eval) logical :: do_getprs_tl,do_normal_rh_to_q,do_tv_to_tsen,do_getuv,lstrong_bk_vars logical :: do_tlnmc,do_q_copy +integer(i_kind) :: ig ! **************************************************************************** ! Initialize timer @@ -117,17 +119,21 @@ subroutine ensctl2model(xhat,mval,eval) do_tlnmc = lstrong_bk_vars .and. ( (tlnmc_option==3) .or. & (jj==ibin_anl .and. tlnmc_option==2) ) - allocate(ebundle(n_ens)) - do nn=1,n_ens - call gsi_bundlecreate (ebundle(nn),xhat%aens(1,1),'c2m ensemble work',istatus) - if(istatus/=0) then - write(6,*) trim(myname), ': trouble creating work ens-bundle' - call stop2(999) - endif + allocate(ebundle(naensgrp,n_ens)) + do ig=1,naensgrp + do nn=1,n_ens + call gsi_bundlecreate (ebundle(ig,nn),xhat%aens(1,1,1),'c2m ensemble work',istatus) + if(istatus/=0) then + write(6,*) trim(myname), ': trouble creating work ens-bundle' + call stop2(999) + endif + enddo enddo ! Apply square-root of ensemble error covariance - call ckgcov_a_en_new_factorization(xhat%aens(jj,1)%values(:),ebundle) + do ig=1,naensgrp + call ckgcov_a_en_new_factorization(ig,xhat%aens(jj,ig,1)%values(:),ebundle(ig,:)) + enddo call sqrt_beta_e_mult(ebundle) ! Initialize ensemble contribution to zero @@ -230,12 +236,14 @@ subroutine ensctl2model(xhat,mval,eval) call stop2(999) endif - do nn=n_ens,1,-1 ! first in; last out - call gsi_bundledestroy(ebundle(nn),istatus) - if(istatus/=0) then - write(6,*) trim(myname), ': trouble destroying work ens bundle, ', istatus - call stop2(999) - endif + do ig=1,naensgrp + do nn=n_ens,1,-1 ! first in; last out + call gsi_bundledestroy(ebundle(ig,nn),istatus) + if(istatus/=0) then + write(6,*) trim(myname), ': trouble destroying work ens bundle, ', istatus + call stop2(999) + endif + enddo enddo deallocate(ebundle) diff --git a/src/gsi/ensctl2model_ad.f90 b/src/gsi/ensctl2model_ad.f90 index 8c491b6de8..706dafc59c 100644 --- a/src/gsi/ensctl2model_ad.f90 +++ b/src/gsi/ensctl2model_ad.f90 @@ -42,6 +42,7 @@ subroutine ensctl2model_ad(eval,mval,grad) use gsi_metguess_mod, only: gsi_metguess_get use mod_strong, only: tlnmc_option use timermod, only: timer_ini,timer_fnl +use hybrid_ensemble_parameters,only: naensgrp implicit none ! Declare passed variables @@ -57,7 +58,7 @@ subroutine ensctl2model_ad(eval,mval,grad) integer(i_kind), parameter :: ncvars = 5 integer(i_kind) :: icps(ncvars) type(gsi_bundle):: wbundle_c ! work bundle -type(gsi_bundle),allocatable :: ebundle(:) +type(gsi_bundle),allocatable :: ebundle(:,:) real(r_kind) :: grade(nval_lenz_en) character(len=3), parameter :: mycvars(ncvars) = (/ & ! vars from CV needed here 'sf ', 'vp ', 'ps ', 't ', & @@ -77,7 +78,7 @@ subroutine ensctl2model_ad(eval,mval,grad) logical :: do_getuv,do_tv_to_tsen_ad,do_normal_rh_to_q_ad,do_getprs_ad logical :: do_tlnmc,lstrong_bk_vars,do_q_copy - +integer(i_kind) :: ig !**************************************************************************** ! Initialize timer @@ -122,13 +123,15 @@ subroutine ensctl2model_ad(eval,mval,grad) (jj==ibin_anl .and. tlnmc_option==2) ) !allocate(grade(nval_lenz_en)) - allocate(ebundle(n_ens)) - do nn=1,n_ens - call gsi_bundlecreate (ebundle(nn),grad%aens(1,1),'m2c ensemble work',istatus) - if(istatus/=0) then - write(6,*) trim(myname), ': trouble creating work ens-bundle' - call stop2(999) - endif + allocate(ebundle(naensgrp,n_ens)) + do ig=1,naensgrp + do nn=1,n_ens + call gsi_bundlecreate (ebundle(ig,nn),grad%aens(1,1,1),'m2c ensemble work',istatus) + if(istatus/=0) then + write(6,*) trim(myname), ': trouble creating work ens-bundle' + call stop2(999) + endif + enddo enddo ! Create a temporary bundle similar to grad, and copy contents of grad into it @@ -210,8 +213,10 @@ subroutine ensctl2model_ad(eval,mval,grad) if(do_getprs_ad) call getprs_ad(cv_ps,cv_tv,rv_prse) end if - do nn=1,n_ens - ebundle(nn)%values=grad%aens(jj,nn)%values + do ig=1,naensgrp + do nn=1,n_ens + ebundle(ig,nn)%values=grad%aens(jj,ig,nn)%values + enddo enddo if(dual_res) then call ensemble_forward_model_ad_dual_res(wbundle_c,ebundle,jj) @@ -222,7 +227,23 @@ subroutine ensctl2model_ad(eval,mval,grad) ! Apply square-root of ensemble error covariance call sqrt_beta_e_mult(ebundle) - call ckgcov_a_en_new_factorization_ad(grade,ebundle) + do ig=1,naensgrp + call ckgcov_a_en_new_factorization_ad(ig,grade,ebundle(ig,:)) + + do ii=1,n_ens + grad%aens(jj,ig,ii)%values=grad%aens(jj,ig,ii)%values + grade + enddo + enddo + + do ig=1,naensgrp + do nn=n_ens,1,-1 ! first in; last out + call gsi_bundledestroy(ebundle(ig,nn),istatus) + if(istatus/=0) then + write(6,*) trim(myname), ': trouble destroying work ens bundle', ig, nn, istatus + call stop2(999) + endif + enddo + enddo call gsi_bundledestroy(wbundle_c,istatus) if (istatus/=0) then @@ -230,20 +251,6 @@ subroutine ensctl2model_ad(eval,mval,grad) call stop2(999) endif - do ii=1,n_ens - grad%aens(jj,ii)%values=grad%aens(jj,ii)%values + grade - enddo -! do ii=1,nval_lenz_en -! grad%aens(jj,1)%values(ii)=grad%aens(jj,1)%values(ii)+grade(ii) -! enddo - - do nn=n_ens,1,-1 ! first in; last out - call gsi_bundledestroy(ebundle(nn),istatus) - if(istatus/=0) then - write(6,*) trim(myname), ': trouble destroying work ens bundle', nn, istatus - call stop2(999) - endif - enddo deallocate(ebundle) ! deallocate(grade) diff --git a/src/gsi/ensctl2state.f90 b/src/gsi/ensctl2state.f90 index 4214b01d85..0d6d3042c5 100644 --- a/src/gsi/ensctl2state.f90 +++ b/src/gsi/ensctl2state.f90 @@ -165,9 +165,9 @@ subroutine ensctl2state(xhat,mval,eval) ! For 4densvar, this is the "3D/Time-invariant contribution from static B" if(dual_res) then - call ensemble_forward_model_dual_res(wbundle_c,xhat%aens(1,:),jj) + call ensemble_forward_model_dual_res(wbundle_c,xhat%aens(1,:,:),jj) else - call ensemble_forward_model(wbundle_c,xhat%aens(1,:),jj) + call ensemble_forward_model(wbundle_c,xhat%aens(1,:,:),jj) end if ! Get pointers to required state variables diff --git a/src/gsi/ensctl2state_ad.f90 b/src/gsi/ensctl2state_ad.f90 index 212d0312c9..4c038c8c6e 100644 --- a/src/gsi/ensctl2state_ad.f90 +++ b/src/gsi/ensctl2state_ad.f90 @@ -260,9 +260,9 @@ subroutine ensctl2state_ad(eval,mval,grad) !$omp end parallel sections if(dual_res) then - call ensemble_forward_model_ad_dual_res(wbundle_c,grad%aens(1,:),jj) + call ensemble_forward_model_ad_dual_res(wbundle_c,grad%aens(1,:,:),jj) else - call ensemble_forward_model_ad(wbundle_c,grad%aens(1,:),jj) + call ensemble_forward_model_ad(wbundle_c,grad%aens(1,:,:),jj) end if end do diff --git a/src/gsi/get_gefs_ensperts_dualres.f90 b/src/gsi/get_gefs_ensperts_dualres.f90 index c547a9f062..fa3d0ecbdd 100644 --- a/src/gsi/get_gefs_ensperts_dualres.f90 +++ b/src/gsi/get_gefs_ensperts_dualres.f90 @@ -93,32 +93,32 @@ subroutine get_gefs_ensperts_dualres type(sub2grid_info) :: grd_tmp ! Create perturbations grid and get variable names from perturbations - if(en_perts(1,1)%grid%im/=grd_ens%lat2.or. & - en_perts(1,1)%grid%jm/=grd_ens%lon2.or. & - en_perts(1,1)%grid%km/=grd_ens%nsig ) then + if(en_perts(1,1,1)%grid%im/=grd_ens%lat2.or. & + en_perts(1,1,1)%grid%jm/=grd_ens%lon2.or. & + en_perts(1,1,1)%grid%km/=grd_ens%nsig ) then if (mype==0) then write(6,*) 'get_gefs_ensperts_dualres: grd_ens ', grd_ens%lat2,grd_ens%lon2,grd_ens%nsig - write(6,*) 'get_gefs_ensperts_dualres: pertgrid', en_perts(1,1)%grid%im, en_perts(1,1)%grid%jm, en_perts(1,1)%grid%km + write(6,*) 'get_gefs_ensperts_dualres: pertgrid', en_perts(1,1,1)%grid%im, en_perts(1,1,1)%grid%jm, en_perts(1,1,1)%grid%km write(6,*) 'get_gefs_ensperts_dualres: inconsistent dims, aborting ...' endif call stop2(999) endif - call gsi_bundlegetpointer (en_perts(1,1),cvars3d,ipc3d,istatus) + call gsi_bundlegetpointer (en_perts(1,1,1),cvars3d,ipc3d,istatus) if(istatus/=0) then write(6,*) ' get_gefs_ensperts_dualres',': cannot find 3d pointers' call stop2(999) endif - call gsi_bundlegetpointer (en_perts(1,1),cvars2d,ipc2d,istatus) + call gsi_bundlegetpointer (en_perts(1,1,1),cvars2d,ipc2d,istatus) if(istatus/=0) then write(6,*) ' get_gefs_ensperts_dualres',': cannot find 2d pointers' call stop2(999) endif - im=en_perts(1,1)%grid%im - jm=en_perts(1,1)%grid%jm - km=en_perts(1,1)%grid%km + im=en_perts(1,1,1)%grid%im + jm=en_perts(1,1,1)%grid%jm + km=en_perts(1,1,1)%grid%km bar_norm = one/float(n_ens) sig_norm=sqrt(one/max(one,n_ens-one)) @@ -126,14 +126,14 @@ subroutine get_gefs_ensperts_dualres call gsi_enscoupler_create_sub2grid_info(grd_tmp,km,npe,grd_ens) ! Allocate bundle to hold mean of ensemble members - call gsi_bundlecreate(en_bar,en_perts(1,1)%grid,'ensemble',istatus,names2d=cvars2d,names3d=cvars3d) + call gsi_bundlecreate(en_bar,en_perts(1,1,1)%grid,'ensemble',istatus,names2d=cvars2d,names3d=cvars3d) if ( istatus /= 0 ) & call die('get_gefs_ensperts_dualres',': trouble creating en_bar bundle, istatus =',istatus) ! Allocate bundle used for reading members allocate(en_read(n_ens)) do n=1,n_ens - call gsi_bundlecreate(en_read(n),en_perts(1,1)%grid,'ensemble member',istatus,names2d=cvars2d,names3d=cvars3d) + call gsi_bundlecreate(en_read(n),en_perts(1,1,1)%grid,'ensemble member',istatus,names2d=cvars2d,names3d=cvars3d) if ( istatus /= 0 ) & call die('get_gefs_ensperts_dualres',': trouble creating en_read bundle, istatus =',istatus) end do @@ -147,7 +147,7 @@ subroutine get_gefs_ensperts_dualres !$omp parallel do schedule(dynamic,1) private(m,n) do m=1,ntlevs_ens do n=1,n_ens - en_perts(n,m)%valuesr4=zero_single + en_perts(n,1,m)%valuesr4=zero_single end do end do @@ -254,7 +254,7 @@ subroutine get_gefs_ensperts_dualres end do !c3d do i=1,nelen - en_perts(n,m)%valuesr4(i)=en_read(n)%values(i) + en_perts(n,1,m)%valuesr4(i)=en_read(n)%values(i) en_bar%values(i)=en_bar%values(i)+en_read(n)%values(i) end do @@ -295,7 +295,7 @@ subroutine get_gefs_ensperts_dualres !$omp parallel do schedule(dynamic,1) private(n,i,ic3,ipic,k,j) do n=1,n_ens do i=1,nelen - en_perts(n,m)%valuesr4(i)=en_perts(n,m)%valuesr4(i)-en_bar%values(i) + en_perts(n,1,m)%valuesr4(i)=en_perts(n,1,m)%valuesr4(i)-en_bar%values(i) end do if(.not. q_hyb_ens) then do ic3=1,nc3d @@ -304,8 +304,8 @@ subroutine get_gefs_ensperts_dualres do k=1,km do j=1,jm do i=1,im - en_perts(n,m)%r3(ipic)%qr4(i,j,k) = min(en_perts(n,m)%r3(ipic)%qr4(i,j,k),limqens) - en_perts(n,m)%r3(ipic)%qr4(i,j,k) = max(en_perts(n,m)%r3(ipic)%qr4(i,j,k),-limqens) + en_perts(n,1,m)%r3(ipic)%qr4(i,j,k) = min(en_perts(n,1,m)%r3(ipic)%qr4(i,j,k),limqens) + en_perts(n,1,m)%r3(ipic)%qr4(i,j,k) = max(en_perts(n,1,m)%r3(ipic)%qr4(i,j,k),-limqens) end do end do end do @@ -313,7 +313,7 @@ subroutine get_gefs_ensperts_dualres end do end if do i=1,nelen - en_perts(n,m)%valuesr4(i)=en_perts(n,m)%valuesr4(i)*sig_norm + en_perts(n,1,m)%valuesr4(i)=en_perts(n,1,m)%valuesr4(i)*sig_norm end do end do end do ntlevs_ens_loop !end do over bins @@ -455,7 +455,7 @@ subroutine ens_spread_dualres(en_bar,ibin) do n=1,n_ens do i=1,nelen sube%values(i)=sube%values(i) & - +(en_perts(n,ibin)%valuesr4(i)-en_bar%values(i))*(en_perts(n,ibin)%valuesr4(i)-en_bar%values(i)) + +(en_perts(n,1,ibin)%valuesr4(i)-en_bar%values(i))*(en_perts(n,1,ibin)%valuesr4(i)-en_bar%values(i)) end do end do diff --git a/src/gsi/get_gefs_for_regional.f90 b/src/gsi/get_gefs_for_regional.f90 index 958849a177..a076f0ccfd 100644 --- a/src/gsi/get_gefs_for_regional.f90 +++ b/src/gsi/get_gefs_for_regional.f90 @@ -1333,9 +1333,9 @@ subroutine get_gefs_for_regional do ic3=1,nc3d if(ntlevs_ens > 1) then - call gsi_bundlegetpointer(en_perts(n,it),trim(cvars3d(ic3)),w3,istatus) + call gsi_bundlegetpointer(en_perts(n,1,it),trim(cvars3d(ic3)),w3,istatus) else - call gsi_bundlegetpointer(en_perts(n,1),trim(cvars3d(ic3)),w3,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars3d(ic3)),w3,istatus) endif if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n @@ -1422,9 +1422,9 @@ subroutine get_gefs_for_regional do ic2=1,nc2d if(ntlevs_ens > 1) then - call gsi_bundlegetpointer(en_perts(n,it),trim(cvars2d(ic2)),w2,istatus) + call gsi_bundlegetpointer(en_perts(n,1,it),trim(cvars2d(ic2)),w2,istatus) else - call gsi_bundlegetpointer(en_perts(n,1),trim(cvars2d(ic2)),w2,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars2d(ic2)),w2,istatus) endif if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n diff --git a/src/gsi/get_nmmb_ensperts.f90 b/src/gsi/get_nmmb_ensperts.f90 index 93d23c837d..ece1780c03 100644 --- a/src/gsi/get_nmmb_ensperts.f90 +++ b/src/gsi/get_nmmb_ensperts.f90 @@ -59,7 +59,7 @@ subroutine get_nmmb_ensperts endif do n=1,n_ens - en_perts(n,1)%valuesr4=zero + en_perts(n,1,1)%valuesr4=zero end do en_bar%values=zero @@ -116,7 +116,7 @@ subroutine get_nmmb_ensperts do ic3=1,nc3d - call gsi_bundlegetpointer(en_perts(n,1),trim(cvars3d(ic3)),w3,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars3d(ic3)),w3,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n,' in get_nmmb_ensperts' call stop2(999) @@ -277,7 +277,7 @@ subroutine get_nmmb_ensperts do ic2=1,nc2d - call gsi_bundlegetpointer(en_perts(n,1),trim(cvars2d(ic2)),w2,istatus) + call gsi_bundlegetpointer(en_perts(n,1,1),trim(cvars2d(ic2)),w2,istatus) if(istatus/=0) then write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n, ' in get_nmmb_ensperts' call stop2(999) @@ -344,7 +344,7 @@ subroutine get_nmmb_ensperts do n=1,n_ens do i=1,nelen - en_perts(n,1)%valuesr4(i)=(en_perts(n,1)%valuesr4(i)-en_bar%values(i))*sig_norm + en_perts(n,1,1)%valuesr4(i)=(en_perts(n,1,1)%valuesr4(i)-en_bar%values(i))*sig_norm end do end do diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index ff636bd70f..859d4cedb1 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -149,7 +149,9 @@ module gsimod beta_s0,beta_e0,s_ens_h,s_ens_v,init_hybrid_ensemble_parameters,& 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 + l_ens_in_diff_time,ensemble_path,ens_fast_read,sst_staticB,limqens, & + ntotensgrp,nsclgrp,naensgrp,ngvarloc,ntlevs_ens,naensloc, & + i_ensloccov4tim,i_ensloccov4var,i_ensloccov4scl,l_timloc_opt 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,& @@ -497,6 +499,7 @@ module gsimod ! 2. fv3_regional = .true. ! 3. fv3_cmaq_regional = .true. ! 4. berror_fv3_cmaq_regional = .true. +! 09-15-2022 yokota - add scale/variable/time-dependent localization ! !EOP !------------------------------------------------------------------------- @@ -1309,9 +1312,19 @@ module gsimod ! beta_e0 - default weight given to ensemble background error covariance ! (if .not. readin_beta). if beta_e0<0, then it is set to ! 1.-beta_s0 (this is the default) -! s_ens_h - homogeneous isotropic horizontal ensemble localization scale (km) -! s_ens_v - vertical localization scale (grid units for now) -! s_ens_h, s_ens_v, and beta_s0 are tunable parameters. +! s_ens_h - horizontal localization correlation length of Gaussian exp(-0.5*(r/L)**2) +! (units of km), default = 2828.0 +! s_ens_v - vertical localization correlation length of Gaussian exp(-0.5*(r/L)**2) +! (grid units if s_ens_v>=0, or units of ln(p) if s_ens_v<0), default = 30.0 +! in scale/variable/time-dependent localization (SDL/VDL/TDL), +! localization length for i-th scale, j-th variable, and k-th time is +! s_ens_[hv]( i + nsclgrp*(j-1) + nsclgrp*ngvarloc*(k-1) ) +! in SDL(nsclgrp>1), i = 1(largest scale) .. nsclgrp(smallest scale) +! in VDL(ngvarloc=2), j = 1(itracer<=10) .. 2(itracer>=11) +! in TDL(l_timloc_opt=true), k = 1(first time bin) .. ntlevs_ens(last time bin) +! in SDL, scale separation length for i-th scale is also set here as +! s_ens_[hv]( naensgrp+i ) - naensgrp is the total number of localization lengths for SDL/VDL/TDL +! in applying SDL only horizontally, set s_ens_v(naensgrp+i)=0.0 ! use_gfs_ens - controls use of global ensemble: .t. use GFS (default); .f. uses user-defined ens ! readin_localization - flag to read (.true.)external localization information file ! readin_beta - flag to read (.true.) the vertically varying beta parameters beta_s and beta_e @@ -1349,14 +1362,27 @@ module gsimod ! ensemble_path - path to ensemble members; default './' ! ens_fast_read - read ensemble in parallel; default '.false.' ! sst_staticB - use only static background error covariance for SST statistic -! -! +! nsclgrp - number of scale-dependent localization lengths +! l_timloc_opt - if true, then turn on time-dependent localization +! ngvarloc - number of variable-dependent localization lengths +! naensloc - total number of spatial localization lengths and scale separation lengths (should be naensgrp+nsclgrp-1) +! i_ensloccov4tim - flag of cross-temporal localization +! =0: cross-temporal covariance is retained +! =1: cross-temporal covariance is zero +! i_ensloccov4var - flag of cross-variable localization +! =0: cross-variable covariance is retained +! =1: cross-variable covariance is zero +! i_ensloccov4scl - flag of cross-scale localization +! =0: cross-scale covariance is retained +! =1: cross-scale covariance is zero +! 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, & oz_univ_static,write_ens_sprd,use_localization_grid,use_gfs_ens, & - i_en_perts_io,l_ens_in_diff_time,ensemble_path,ens_fast_read,sst_staticB,limqens + i_en_perts_io,l_ens_in_diff_time,ensemble_path,ens_fast_read,sst_staticB,limqens, & + nsclgrp,l_timloc_opt,ngvarloc,naensloc,i_ensloccov4tim,i_ensloccov4var,i_ensloccov4scl ! rapidrefresh_cldsurf (options for cloud analysis and surface ! enhancement for RR appilcation ): @@ -1807,6 +1833,16 @@ subroutine gsimain_initialize call stop2(137) endif + ntotensgrp=nsclgrp*ngvarloc + if(l_timloc_opt) then + naensgrp=ntotensgrp*ntlevs_ens + else + naensgrp=ntotensgrp + endif + if(naensloc NULL() + character(len=*),parameter::myname_=myname//'*hybens_localization_setup' + + l_read_success=.false. print_verbose=.false. .and. mype == 0 if(verbose .and. mype == 0)print_verbose=.true. @@ -4044,10 +4214,25 @@ subroutine hybens_localization_setup endif if(mype==0) write(6,'(" LOCALIZATION, BETA_S, BETA_E VERTICAL PROFILES FOLLOW")') do k = 1,grd_ens%nsig - read(lunin,101) s_ens_hv(k), s_ens_vv(k), beta_s(k), beta_e(k) - if(mype==0) write(6,101) s_ens_hv(k), s_ens_vv(k), beta_s(k), beta_e(k) + read(lunin,101) s_ens_hv(k,1), s_ens_vv(k,1), beta_s(k), beta_e(k) + if(mype==0) write(6,101) s_ens_hv(k,1), s_ens_vv(k,1), beta_s(k), beta_e(k) enddo + do ig=2,naensloc + do k = 1,grd_ens%nsig + read(lunin,101,end=300) s_ens_hv(k,ig),s_ens_vv(k,ig) + enddo + enddo + l_read_success=.true. close(lunin) +300 continue + if(.not.l_read_success) then + do ig=2,naensloc + do k = 1,grd_ens%nsig + s_ens_hv(k,ig)=s_ens_hv(k,1) + s_ens_vv(k,ig)=s_ens_vv(k,1) + enddo + enddo + endif else @@ -4060,8 +4245,8 @@ subroutine hybens_localization_setup vvlocal = .true. nz = msig kl = grd_loc%kend_alloc-grd_loc%kbegin_loc+1 - if(.not.allocated(s_ens_h_gu_x)) allocate(s_ens_h_gu_x(grd_loc%nsig*n_ens)) - if(.not.allocated(s_ens_h_gu_y)) allocate(s_ens_h_gu_y(grd_loc%nsig*n_ens)) + if(.not.allocated(s_ens_h_gu_x)) allocate(s_ens_h_gu_x(grd_loc%nsig*n_ens,naensloc)) + if(.not.allocated(s_ens_h_gu_y)) allocate(s_ens_h_gu_y(grd_loc%nsig*n_ens,naensloc)) endif endif ! if ( readin_localization .or. readin_beta ) @@ -4091,10 +4276,12 @@ subroutine hybens_localization_setup if ( .not. readin_localization ) then ! assign all levels to same value, s_ens_h, s_ens_v nz = 1 kl = 1 - if(.not.allocated(s_ens_h_gu_x)) allocate(s_ens_h_gu_x(1)) - if(.not.allocated(s_ens_h_gu_y)) allocate(s_ens_h_gu_y(1)) - s_ens_hv = s_ens_h - s_ens_vv = s_ens_v + if(.not.allocated(s_ens_h_gu_x)) allocate(s_ens_h_gu_x(1,naensloc)) + if(.not.allocated(s_ens_h_gu_y)) allocate(s_ens_h_gu_y(1,naensloc)) + do ig=1,naensloc + s_ens_hv(:,ig) = s_ens_h(ig) + s_ens_vv(:,ig) = s_ens_v(ig) + enddo endif ! Set up localization filters @@ -4103,11 +4290,19 @@ subroutine hybens_localization_setup call normal_new_factorization_rf_z if ( regional ) then ! convert s_ens_h from km to grid units. - call convert_km_to_grid_units(s_ens_h_gu_x,s_ens_h_gu_y,nz) if ( vvlocal ) then - call init_rf_x(s_ens_h_gu_x(grd_loc%kbegin_loc:grd_loc%kend_alloc),kl) - call init_rf_y(s_ens_h_gu_y(grd_loc%kbegin_loc:grd_loc%kend_alloc),kl) + call convert_km_to_grid_units(s_ens_h_gu_x(1:nz,:),s_ens_h_gu_y(1:nz,:),nz) + do n=2,n_ens + nk=(n-1)*nz + do k=1,nz + s_ens_h_gu_x(nk+k,:)=s_ens_h_gu_x(k,:) + s_ens_h_gu_y(nk+k,:)=s_ens_h_gu_y(k,:) + enddo + enddo + call init_rf_x(s_ens_h_gu_x(grd_loc%kbegin_loc:grd_loc%kend_alloc,:),kl) + call init_rf_y(s_ens_h_gu_y(grd_loc%kbegin_loc:grd_loc%kend_alloc,:),kl) else + call convert_km_to_grid_units(s_ens_h_gu_x,s_ens_h_gu_y,nz) call init_rf_x(s_ens_h_gu_x,kl) call init_rf_y(s_ens_h_gu_y,kl) endif @@ -4117,6 +4312,80 @@ subroutine hybens_localization_setup call init_sf_xy(jcap_ens) endif + if(ntotensgrp>1) then + call gsi_bundlegetpointer(en_perts(1,1,1),cvars3d,ipc3d,istatus) + if(istatus/=0) then + write(6,*) myname_,': cannot find 3d pointers' + call stop2(999) + endif + call gsi_bundlegetpointer(en_perts(1,1,1),cvars2d,ipc2d,istatus) + if(istatus/=0) then + write(6,*) myname_,': cannot find 2d pointers' + call stop2(999) + endif + if(nsclgrp>1) then + call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig) + allocate(values(grd_ens%latlon11*grd_ens%nsig*n_ens)) + do ig=1,nsclgrp-1 + ii=0 + do n=1,n_ens + a_en(n)%values => values(ii+1:ii+grd_ens%latlon11*grd_ens%nsig) + call gsi_bundleset(a_en(n),grid_ens,'Ensemble Bundle',istatus,names3d=(/'a_en'/),bundle_kind=r_kind) + if (istatus/=0) then + write(6,*) myname_,': error alloc(ensemble bundle)' + call stop2(999) + endif + ii=ii+grd_ens%latlon11*grd_ens%nsig + enddo + do m=1,ntlevs_ens + do n=1,n_ens + en_perts(n,ig+1,m)%valuesr4=en_perts(n,ig,m)%valuesr4 + enddo + do ic3=1,nc3d + ipic=ipc3d(ic3) + do n=1,n_ens + do k=1,grd_ens%nsig + a_en(n)%r3(1)%q(:,:,k)=en_perts(n,ig,m)%r3(ipic)%qr4(:,:,k) + enddo + enddo + call bkgcov_a_en_new_factorization(naensgrp+ig,a_en) + do n=1,n_ens + do k=1,grd_ens%nsig + en_perts(n,ig,m)%r3(ipic)%qr4(:,:,k)=a_en(n)%r3(1)%q(:,:,k) + enddo + enddo + enddo + do ic2=1,nc2d + ipic=ipc2d(ic2) + do n=1,n_ens + do k=1,grd_ens%nsig + a_en(n)%r3(1)%q(:,:,k)=en_perts(n,ig,m)%r2(ipic)%qr4(:,:) + enddo + enddo + call bkgcov_a_en_new_factorization(naensgrp+ig,a_en) + do n=1,n_ens + en_perts(n,ig,m)%r2(ipic)%qr4(:,:)=a_en(n)%r3(1)%q(:,:,1) + enddo + enddo + do n=1,n_ens + en_perts(n,ig+1,m)%valuesr4=en_perts(n,ig+1,m)%valuesr4-en_perts(n,ig,m)%valuesr4 + enddo + enddo + do n=1,n_ens + call gsi_bundleunset(a_en(n),istatus) + enddo + enddo + deallocate(values) + endif + do ig=nsclgrp+1,ntotensgrp + do m=1,ntlevs_ens + do n=1,n_ens + en_perts(n,ig,m)%valuesr4=en_perts(n,ig-nsclgrp,m)%valuesr4 + enddo + enddo + enddo + endif + !!!!!!!! setup beta_s, beta_e!!!!!!!!!!!! ! vertical variation of static and ensemble weights @@ -4138,12 +4407,14 @@ subroutine hybens_localization_setup ! write out final values for s_ens_hv, s_ens_vv, beta_s, beta_e if ( print_verbose ) then - write(6,*) 'HYBENS_LOCALIZATION_SETUP: s_ens_hv,s_ens_vv,beta_s,beta_e' + write(6,*) 'HYBENS_LOCALIZATION_SETUP: s_ens_hv(:,1),s_ens_vv(:,1),beta_s,beta_e' do k=1,grd_ens%nsig - write(6,101) s_ens_hv(k), s_ens_vv(k), beta_s(k), beta_e(k) + write(6,101) s_ens_hv(k,1), s_ens_vv(k,1), beta_s(k), beta_e(k) enddo endif + call setup_ensgrp2aensgrp + return end subroutine hybens_localization_setup @@ -4174,16 +4445,17 @@ subroutine convert_km_to_grid_units(s_ens_h_gu_x,s_ens_h_gu_y,nz) !$$$ use kinds, only: r_kind,i_kind - use hybrid_ensemble_parameters, only: grd_loc,n_ens,s_ens_hv + use hybrid_ensemble_parameters, only: s_ens_hv use hybrid_ensemble_parameters, only: region_dx_ens,region_dy_ens + use hybrid_ensemble_parameters, only: naensloc use gsi_io, only: verbose implicit none integer(i_kind) ,intent(in ) ::nz - real(r_kind),intent( out) ::s_ens_h_gu_x(nz*n_ens),s_ens_h_gu_y(nz*n_ens) + real(r_kind),intent( out) ::s_ens_h_gu_x(nz,naensloc),s_ens_h_gu_y(nz,naensloc) logical :: print_verbose real(r_kind) dxmax,dymax - integer(i_kind) k,n,nk + integer(i_kind) k print_verbose=.false. if(verbose) print_verbose=.true. @@ -4197,25 +4469,14 @@ subroutine convert_km_to_grid_units(s_ens_h_gu_x,s_ens_h_gu_y,nz) end if do k=1,nz - s_ens_h_gu_x(k)=s_ens_hv(k)/(.001_r_kind*dxmax) - s_ens_h_gu_y(k)=s_ens_hv(k)/(.001_r_kind*dymax) - if(print_verbose) write(6,*)' in convert_km_to_grid_units,s_ens_h,s_ens_h_gu_x,y=', & - s_ens_hv(k),s_ens_h_gu_x(k),s_ens_h_gu_y(k) - + s_ens_h_gu_x(k,:)=s_ens_hv(k,:)/(.001_r_kind*dxmax) + s_ens_h_gu_y(k,:)=s_ens_hv(k,:)/(.001_r_kind*dymax) + if(print_verbose) write(6,*)' in convert_km_to_grid_units,s_ens_h(k,1),s_ens_h_gu_x,y(k,1)=', & + s_ens_hv(k,1),s_ens_h_gu_x(k,1),s_ens_h_gu_y(k,1) enddo - if(nz>1)then - do n=2,n_ens - nk=(n-1)*grd_loc%nsig - do k=1,grd_loc%nsig - s_ens_h_gu_x(nk+k)=s_ens_h_gu_x(k) - s_ens_h_gu_y(nk+k)=s_ens_h_gu_y(k) - enddo - enddo - endif return - end subroutine convert_km_to_grid_units subroutine grads1(f,nvert,mype,fname) @@ -5146,4 +5407,116 @@ subroutine setup_pwgt end subroutine setup_pwgt +subroutine setup_ensgrp2aensgrp +!$$$ subprogram documentation block +! . . . . +! subprogram: set a matrix of (naensgrp,naensgrp) +! +! program history log: +! 2022-09-15 yokota - add scale/variable/time-dependent localization +! +! input argument list: +! +! output argument list: +! +! remarks: +! need to reconcile grid in gsi_bundle w/ grid_ens/grid_anl +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ end documentation block + use constants, only: zero,one + use hybrid_ensemble_parameters, only: l_timloc_opt,i_ensloccov4tim,i_ensloccov4var,i_ensloccov4scl + use hybrid_ensemble_parameters, only: ensloccov4tim,ensloccov4var,ensloccov4scl + use hybrid_ensemble_parameters, only: ntotensgrp,naensgrp,ntlevs_ens,nsclgrp,ngvarloc + use hybrid_ensemble_parameters, only: ensgrp2aensgrp + use hybrid_ensemble_parameters, only: idaen2d,idaen3d + use hybrid_ensemble_parameters, only: alphacvarsclgrpmat + implicit none + + integer (i_kind):: i,j + integer (i_kind):: ig,ibin,ic,iscl,itim1,itim2,igvar1,igvar2,iscl1,iscl2,ivargrp + integer (i_kind):: ntimloc,interval4aens + + if (l_timloc_opt) then + ntimloc=ntlevs_ens + interval4aens=ntotensgrp + else + ntimloc=1 + interval4aens=0 + endif + if(naensgrp/=ntimloc*ngvarloc*nsclgrp) then + write(6,*)'setup_ensgrp2aensgrp: wrong naensgrp' + call stop2(666) + endif + if(ntotensgrp/=ngvarloc*nsclgrp) then + write(6,*)'setup_ensgrp2aensgrp: wrong ntotensgrp' + call stop2(666) + endif + ensgrp2aensgrp=-999 + do ibin=1,ntlevs_ens + do ic=1,nc3d+nc2d + if(ngvarloc>1) then + if(ic<=nc3d) ivargrp=idaen3d(ic) + if(ic> nc3d) ivargrp=idaen2d(ic-nc3d) + else + ivargrp=1 + endif + do iscl=1,nsclgrp + ig=(ivargrp-1)*nsclgrp+iscl + ensgrp2aensgrp(ig,ic,ibin)=(ibin-1)*interval4aens+(ivargrp-1)*nsclgrp+iscl + enddo + enddo + enddo + + if (i_ensloccov4tim==0) then + ensloccov4tim=one + elseif (i_ensloccov4tim==1)then + ensloccov4tim=zero + ensloccov4tim(1)=one + else + write(6,*)'setup_ensgrp2aensgrp: wrong i_ensloccov4tim' + call stop2(666) + endif + if (i_ensloccov4var==0) then + ensloccov4var=one + elseif (i_ensloccov4var==1)then + ensloccov4var=zero + ensloccov4var(1)=one + else + write(6,*)'setup_ensgrp2aensgrp: wrong i_ensloccov4var' + call stop2(666) + endif + if (i_ensloccov4scl==0) then + ensloccov4scl=one + elseif (i_ensloccov4scl==1)then + ensloccov4scl=zero + ensloccov4scl(1)=one + else + write(6,*)'setup_ensgrp2aensgrp: wrong i_ensloccov4scl' + call stop2(666) + endif + + do itim2=1,ntimloc + do itim1=1,ntimloc + do igvar1=1,ngvarloc + do igvar2=1,ngvarloc + do iscl1=1,nsclgrp + do iscl2=1,nsclgrp + i=(itim1-1)*interval4aens+(igvar1-1)*nsclgrp+iscl1 + j=(itim2-1)*interval4aens+(igvar2-1)*nsclgrp+iscl2 + alphacvarsclgrpmat(i,j)=ensloccov4tim(abs(itim1-itim2)+1) & + *ensloccov4var(abs(igvar1-igvar2)+1) & + *ensloccov4scl(abs(iscl1-iscl2)+1) !first for ttime covariance + enddo + enddo + enddo + enddo + enddo + enddo + +end subroutine setup_ensgrp2aensgrp + end module hybrid_ensemble_isotropic diff --git a/src/gsi/hybrid_ensemble_parameters.f90 b/src/gsi/hybrid_ensemble_parameters.f90 index 9ca7770c41..7b1c963764 100644 --- a/src/gsi/hybrid_ensemble_parameters.f90 +++ b/src/gsi/hybrid_ensemble_parameters.f90 @@ -90,8 +90,19 @@ module hybrid_ensemble_parameters ! beta_e0 - default weight given to ensemble background error covariance ! (if .not. readin_beta). if beta_e0<0, then it is set to ! 1.-beta_s0 (this is the default) -! s_ens_h: horizontal localization correlation length (units of km), default = 2828.0 -! s_ens_v: vertical localization correlation length (grid units), default = 30.0 +! s_ens_h: horizontal localization correlation length of Gaussian exp(-0.5*(r/L)**2) +! (units of km), default = 2828.0 +! s_ens_v: vertical localization correlation length of Gaussian exp(-0.5*(r/L)**2) +! (grid units if s_ens_v>=0, or units of ln(p) if s_ens_v<0), default = 30.0 +! in scale/variable/time-dependent localization (SDL/VDL/TDL), +! localization length for i-th scale, j-th variable, and k-th time is +! s_ens_[hv]( i + nsclgrp*(j-1) + nsclgrp*ngvarloc*(k-1) ) +! in SDL(nsclgrp>1), i = 1(largest scale) .. nsclgrp(smallest scale) +! in VDL(ngvarloc=2), j = 1(itracer<=10) .. 2(itracer>=11) +! in TDL(l_timloc_opt=true), k = 1(first time bin) .. ntlevs_ens(last time bin) +! in SDL, scale separation length for i-th scale is also set here as +! s_ens_[hv]( naensgrp+i ) - naensgrp is the total number of localization lengths for SDL/VDL/TDL +! in applying SDL only horizontally, set s_ens_v(naensgrp+i)=0.0 ! generate_ens: if .true., generate ensemble perturbations internally as random samples of background B. ! (used primarily for testing/debugging) ! if .false., read external ensemble perturbations @@ -118,6 +129,19 @@ module hybrid_ensemble_parameters ! ensemble_path: path to ensemble members; default './' ! ens_fast_read: read ensemble in parallel; default '.false.' ! sst_staticB: if .true. (default) uses only static part of B error covariance for SST +! nsclgrp: number of scale-dependent localization lengths +! l_timloc_opt: if true, then turn on time-dependent localization +! ngvarloc: number of variable-dependent localization lengths +! naensloc: total number of spatial localization lengths and scale separation lengths (should be naensgrp+nsclgrp-1) +! i_ensloccov4tim: flag of cross-temporal localization +! =0: cross-temporal covariance is retained +! =1: cross-temporal covariance is zero +! i_ensloccov4var: flag of cross-variable localization +! =0: cross-variable covariance is retained +! =1: cross-variable covariance is zero +! i_ensloccov4scl: flag of cross-scale localization +! =0: cross-scale covariance is retained +! =1: cross-scale covariance is zero !===================================================================================================== ! ! @@ -151,6 +175,7 @@ module hybrid_ensemble_parameters ! 2015-01-22 Hu - add flag i_en_perts_io to control reading ensemble perturbation. ! 2015-02-11 Hu - add flag l_ens_in_diff_time to force GSI hybrid use ensembles not available at analysis time ! 2015-09-18 todling - add sst_staticB to control use of ensemble SST error covariance +! 2022-09-15 yokota - add scale/variable/time-dependent localization ! ! subroutines included: @@ -174,8 +199,10 @@ module hybrid_ensemble_parameters ! beta_s(:) = beta_s0 , vertically varying weights given to B ; ! beta_e(:) = 1 - beta_s0 , vertically varying weights given to A. ! If (readin_beta) then beta_s and beta_e are read from a file and beta_s0 is not used. -! def s_ens_h - homogeneous isotropic horizontal ensemble localization scale (km) -! def s_ens_v - vertical localization scale (grid units for now) +! def s_ens_h - horizontal localization correlation length of Gaussian exp(-0.5*(r/L)**2) +! (units of km), default = 2828.0 +! def s_ens_v - vertical localization correlation length of Gaussian exp(-0.5*(r/L)**2) +! (grid units if s_ens_v>=0, or units of ln(p) if s_ens_v<0), default = 30.0 ! def readin_localization - flag to read (.true.)external localization information file ! def eqspace_ensgrid - if .true., then ensemble grid is equal spaced, staggered 1/2 grid unit off ! poles. if .false., then Gaussian grid assumed for ensemble (global only) @@ -292,12 +319,20 @@ module hybrid_ensemble_parameters public :: en_perts,ps_bar public :: region_lat_ens,region_lon_ens public :: region_dx_ens,region_dy_ens + public :: naensgrp,ntotensgrp,nsclgrp,naensloc,ngvarloc + public :: ensgrp2aensgrp + public :: ensloccov4tim,ensloccov4var,ensloccov4scl + public :: alphacvarsclgrpmat + public :: l_timloc_opt + public :: i_ensloccov4tim,i_ensloccov4var,i_ensloccov4scl + public :: idaen3d,idaen2d public :: ens_fast_read public :: l_both_fv3sar_gfs_ens public :: sst_staticB public :: limqens logical l_hyb_ens,uv_hyb_ens,q_hyb_ens,oz_univ_static,sst_staticB + logical l_timloc_opt logical aniso_a_en logical full_ensemble,pwgtflg logical generate_ens @@ -317,11 +352,16 @@ module hybrid_ensemble_parameters 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 + real(r_kind) beta_s0,beta_e0,grid_ratio_ens + integer(i_kind),parameter::max_naensloc=20 + integer(i_kind),parameter::max_nvars=100 + real(r_kind) s_ens_h(max_naensloc) + real(r_kind) s_ens_v(max_naensloc) type(sub2grid_info),save :: grd_ens,grd_loc,grd_sploc,grd_anl,grd_e1,grd_a1 type(spec_vars),save :: sp_ens,sp_loc type(egrid2agrid_parm),save :: p_e2a,p_sploc2ens - real(r_kind),allocatable,dimension(:) :: s_ens_hv,s_ens_vv + real(r_kind),allocatable,dimension(:,:) :: s_ens_vv + real(r_kind),allocatable,dimension(:,:) :: s_ens_hv real(r_kind),allocatable,dimension(:) :: sqrt_beta_s,sqrt_beta_e real(r_kind),allocatable,dimension(:) :: beta_s,beta_e real(r_kind),allocatable,dimension(:,:,:) :: pwgt @@ -336,6 +376,18 @@ module hybrid_ensemble_parameters integer(i_kind) regional_ensemble_option integer(i_kind) fv3sar_ensemble_opt character(len=512),save :: ensemble_path + real(r_kind),allocatable,dimension(:,:) :: alphacvarsclgrpmat + integer(i_kind),allocatable,dimension(:,:,:) :: ensgrp2aensgrp + real(r_kind),allocatable,dimension(:) :: ensloccov4tim,ensloccov4var,ensloccov4scl + integer(i_kind) :: nsclgrp=1 + integer(i_kind) :: naensgrp=1 + integer(i_kind) :: ntotensgrp=1 + integer(i_kind) :: naensloc=1 + integer(i_kind) :: ngvarloc=1 + integer(i_kind) :: i_ensloccov4tim=0 + integer(i_kind) :: i_ensloccov4var=0 + integer(i_kind) :: i_ensloccov4scl=0 + integer(i_kind),allocatable,dimension(:) :: idaen3d,idaen2d ! following is for storage of ensemble perturbations: @@ -343,7 +395,7 @@ module hybrid_ensemble_parameters ! def nelen - length of one ensemble perturbation vector integer(i_kind) nelen - type(gsi_bundle),save,allocatable :: en_perts(:,:) + type(gsi_bundle),save,allocatable :: en_perts(:,:,:) real(r_single),dimension(:,:,:),allocatable:: ps_bar real(r_single):: limqens @@ -384,6 +436,7 @@ subroutine init_hybrid_ensemble_parameters implicit none l_hyb_ens=.false. + l_timloc_opt=.false. full_ensemble=.false. pwgtflg=.false. uv_hyb_ens=.false. @@ -435,15 +488,23 @@ subroutine create_hybens_localization_parameters use constants, only: zero implicit none - allocate( s_ens_hv(grd_ens%nsig),s_ens_vv(grd_ens%nsig) ) + allocate( s_ens_hv(grd_ens%nsig,naensloc),s_ens_vv(grd_ens%nsig,naensloc) ) allocate( beta_s(grd_ens%nsig),beta_e(grd_ens%nsig)) allocate( sqrt_beta_s(grd_ens%nsig),sqrt_beta_e(grd_ens%nsig) ) allocate( pwgt(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig) ) + allocate( alphacvarsclgrpmat(naensgrp,naensgrp) ) + allocate( ensgrp2aensgrp(ntotensgrp,max_nvars,ntlevs_ens) ) + allocate( ensloccov4tim(ntlevs_ens),ensloccov4var(ngvarloc),ensloccov4scl(nsclgrp) ) beta_s =one beta_e =zero sqrt_beta_s=one sqrt_beta_e=zero pwgt=zero + alphacvarsclgrpmat=one + ensgrp2aensgrp=1 + ensloccov4tim=one + ensloccov4var=one + ensloccov4scl=one end subroutine create_hybens_localization_parameters @@ -453,6 +514,8 @@ subroutine destroy_hybens_localization_parameters deallocate(s_ens_vv,s_ens_hv) deallocate(beta_s,beta_e) deallocate(sqrt_beta_s,sqrt_beta_e,pwgt) + deallocate(alphacvarsclgrpmat) + deallocate(ensgrp2aensgrp) end subroutine destroy_hybens_localization_parameters diff --git a/src/gsi/jfunc.f90 b/src/gsi/jfunc.f90 index d65b4c8fd3..1b92ad2e94 100644 --- a/src/gsi/jfunc.f90 +++ b/src/gsi/jfunc.f90 @@ -738,6 +738,7 @@ subroutine set_pointer use gsi_4dvar, only: nsubwin, lsqrtb use bias_predictors, only: setup_predictors use hybrid_ensemble_parameters, only: l_hyb_ens,n_ens,generate_ens,grd_ens,nval_lenz_en + use hybrid_ensemble_parameters, only: naensgrp implicit none integer(i_kind) n_ensz,nval_lenz_tot,nval_lenz_enz @@ -748,8 +749,8 @@ subroutine set_pointer nval_levs=max(0,nc3d)*nsig+max(0,nc2d) nval_len=nval_levs*latlon11 if(l_hyb_ens) then - nval_len=nval_len+n_ens*nsig*grd_ens%latlon11 - nval_levs_ens=nval_levs+n_ens*nsig + nval_len=nval_len+naensgrp*n_ens*nsig*grd_ens%latlon11 + nval_levs_ens=nval_levs+naensgrp*n_ens*nsig end if nsclen=npred*jpch_rad npclen=npredp*npcptype diff --git a/src/gsi/stub_get_pseudo_ensperts.f90 b/src/gsi/stub_get_pseudo_ensperts.f90 index 7bababc0c3..ca39b4b6e2 100644 --- a/src/gsi/stub_get_pseudo_ensperts.f90 +++ b/src/gsi/stub_get_pseudo_ensperts.f90 @@ -32,7 +32,7 @@ subroutine get_pseudo_ensperts_dummy(this,en_perts,nelen) use kinds, only: i_kind implicit none class(get_pseudo_ensperts_class), intent(inout) :: this - type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:,:) integer(i_kind), intent(in ) :: nelen write(6,*)'get_pseudo_ensperts: ***WARNING*** dummy call ... does nothing!' diff --git a/src/gsi/stub_get_wrf_mass_ensperts.f90 b/src/gsi/stub_get_wrf_mass_ensperts.f90 index 3089cb9336..bf98f3106e 100644 --- a/src/gsi/stub_get_wrf_mass_ensperts.f90 +++ b/src/gsi/stub_get_wrf_mass_ensperts.f90 @@ -38,7 +38,7 @@ subroutine ens_spread_dualres_regional_dummy(this,mype,en_perts,nelen,en_bar) implicit none class(get_wrf_mass_ensperts_class), intent(inout) :: this integer(i_kind),intent(in):: mype - type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:,:) integer(i_kind), intent(in ):: nelen type(gsi_bundle),OPTIONAL,intent(in):: en_bar write(6,*)'get_wrf_mass_ensperts: ***WARNING*** dummy call ... does nothing!' @@ -50,7 +50,7 @@ subroutine get_wrf_mass_ensperts_dummy(this,en_perts,nelen,ps_bar) use gsi_bundlemod, only: gsi_bundle implicit none class(get_wrf_mass_ensperts_class), intent(inout) :: this - type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:) integer(i_kind), intent(in ):: nelen real(r_single),dimension(:,:,:),allocatable:: ps_bar write(6,*)'get_wrf_mass_ensperts: ***WARNING*** dummy call ... does nothing!' diff --git a/src/gsi/stub_get_wrf_nmm_ensperts.f90 b/src/gsi/stub_get_wrf_nmm_ensperts.f90 index 14ab869ffd..db27125d70 100644 --- a/src/gsi/stub_get_wrf_nmm_ensperts.f90 +++ b/src/gsi/stub_get_wrf_nmm_ensperts.f90 @@ -13,7 +13,7 @@ subroutine get_wrf_nmm_ensperts_dummy(this,en_perts,nelen,region_lat_ens,region_ implicit none class(get_wrf_nmm_ensperts_class), intent(inout) :: this - type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:) + type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:) integer(i_kind), intent(in ):: nelen real(r_kind),allocatable, intent(inout):: region_lat_ens(:,:),region_lon_ens(:,:) real(r_single),dimension(:,:,:),allocatable, intent(inout):: ps_bar