From 237ba4a28189c066cd88ab0cb21cc579067ada1f Mon Sep 17 00:00:00 2001 From: "yongming.wang" Date: Sun, 28 Aug 2022 19:20:06 -0500 Subject: [PATCH] GitHub Issue NOAA-EMC/GSI#468. Add capability for multiscale DA using scale- and variable-dependent localization --- fix | 1 - src/gsi/apply_scaledepwgts.f90 | 189 +++ 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 | 101 +- src/gsi/cplr_get_fv3_regional_ensperts.f90 | 370 ++++- 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 | 12 +- src/gsi/ens_spread_mod.f90 | 4 +- src/gsi/ensctl2model.f90 | 44 +- src/gsi/ensctl2model_ad.f90 | 54 +- src/gsi/ensctl2state.f90 | 4 +- src/gsi/ensctl2state_ad.f90 | 4 +- src/gsi/general_specmod.f90 | 41 + src/gsi/get_gefs_ensperts_dualres.f90 | 70 +- src/gsi/get_gefs_for_regional.f90 | 8 +- src/gsi/get_nmmb_ensperts.f90 | 8 +- src/gsi/gsi_files.cmake | 1 + src/gsi/gsi_rfv3io_mod.f90 | 152 +- src/gsi/gsimod.F90 | 8 +- src/gsi/hybrid_ensemble_isotropic.F90 | 1561 ++++++++++++------- src/gsi/hybrid_ensemble_parameters.f90 | 43 +- src/gsi/jfunc.f90 | 6 +- src/gsi/obsmod.F90 | 2 +- src/gsi/pcgsoi.f90 | 1 - src/gsi/read_dbz_nc.f90 | 8 +- src/gsi/setupdbz.f90 | 6 + 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 +- 35 files changed, 2037 insertions(+), 739 deletions(-) delete mode 160000 fix create mode 100644 src/gsi/apply_scaledepwgts.f90 diff --git a/fix b/fix deleted file mode 160000 index 023f81d4b..000000000 --- a/fix +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 023f81d4ba631d235156172e04a529a4bf273617 diff --git a/src/gsi/apply_scaledepwgts.f90 b/src/gsi/apply_scaledepwgts.f90 new file mode 100644 index 000000000..c7a40c169 --- /dev/null +++ b/src/gsi/apply_scaledepwgts.f90 @@ -0,0 +1,189 @@ +!$$$ subprogram documentation block +!28-03-2018 T. Lei and D. Kleist consoliated and added codes +! for the scale dependent scale localization scheme + function fwgtofwvlen (rvlft,rvrgt,rcons,rlen,rinput) + use kinds, only: r_kind,i_kind,r_single + implicit none + real (r_kind):: fwgtofwvlen + real(r_kind):: rvlft,rvrgt,rcons,rlen,rinput + real (r_kind):: rlen1,rtem1,rconshalf +! write(6,*)'rinput and ,rvlft rvlft is ',rinput,rvlft,rvrgt + rlen1=rlen/10.0 ! rlen corresponds to a (-5,5) region + rconshalf=0.5*rcons + if(rinput.gt.rvlft.and.rinput.lt.rvrgt) then + fwgtofwvlen=rcons + else + rtem1=min(abs(rinput-rvlft),abs(rinput-rvrgt)) + fwgtofwvlen=rconshalf*(1.0+tanh(5.0-rtem1/rlen1)) + endif + + + end function fwgtofwvlen +! . . . . + subroutine init_mult_spc_wgts(jcap_in) + use kinds, only: r_kind,i_kind,r_single + use hybrid_ensemble_parameters,only: s_ens_hv,sp_loc,grd_ens,grd_loc,sp_ens,n_ens,p_sploc2ens,grd_sploc + use hybrid_ensemble_parameters,only: use_localization_grid + use gridmod,only: use_sp_eqspace + use general_specmod, only: general_init_spec_vars + use constants, only: zero,half,one,two,three,rearth,pi + use constants, only: rad2deg + use mpimod, only: mype + use general_sub2grid_mod, only: general_sub2grid_create_info + use egrid2agrid_mod,only: g_create_egrid2agrid + use general_sub2grid_mod, only: sub2grid_info + use gsi_io, only: verbose + use hybrid_ensemble_parameters, only: nsclgrp + use hybrid_ensemble_parameters, only: spc_multwgt,spcwgt_params,l_sum_spc_weights + implicit none + + integer(i_kind),intent(in ) :: jcap_in + real(r_kind),allocatable::totwvlength(:) + + integer(i_kind) i,ii,j,k,l,n,kk,nsigend + integer ig + real(r_kind) rwv0,rtem1,rtem2 + real (r_kind):: fwgtofwvlen + + + allocate(totwvlength(jcap_in)) + + rwv0=2*pi*rearth*.001_r_kind + do i=1,jcap_in + totwvlength(i)= rwv0/i + enddo + do i=1,jcap_in + rtem1=0 + do ig=1,nsclgrp + if(ig.ne.2) then + spc_multwgt(i,ig)=fwgtofwvlen(spcwgt_params(1,ig),spcwgt_params(2,ig),& + spcwgt_params(3,ig),spcwgt_params(4,ig),totwvlength(i)) !fwv2wgt(twvlength) + if(l_sum_spc_weights.eq.0 ) then + rtem1=rtem1+spc_multwgt(i,ig) + else + rtem1=rtem1+spc_multwgt(i,ig)**2 + endif + endif + enddo + rtem2=1-rtem1 + if(abs(rtem2).ge.zero) then + + if(l_sum_spc_weights.eq.0 ) then + spc_multwgt(i,2)=rtem2 + else + spc_multwgt(i,2)=sqrt(rtem2) + endif + endif + enddo + spc_multwgt=max(spc_multwgt,0.0) +! if(mype.eq.0) then +! open(121,file="test-spcwght.dat",form="formatted") +! do i=1,jcap_in +! write(121,111)(totwvlength(i)),((spc_multwgt(i,j)),j=1,nsclgrp) +! enddo +! close (121) +! endif +! 111 format(g15.6,3(g11.4,1x)) + + deallocate(totwvlength) + return + end subroutine init_mult_spc_wgts + + + + subroutine apply_scaledepwgts(grd_in,sp_in,wbundle,spwgts,wbundle2,igrp,nensid) +! 2017-03-30 J. Kay, X. Wang - copied from Kleist's apply_scaledepwgts and +! add the calculation of scale-dependent weighting for mixed resolution ensemble +! POC: xuguang.wang@ou.edu +! + use constants, only: one + use control_vectors, only: nrf_var,cvars2d,cvars3d,control_vector + use kinds, only: r_kind,i_kind + use kinds, only: r_single + use mpimod, only: mype,nvar_id,levs_id + use hybrid_ensemble_parameters, only: oz_univ_static + use general_specmod, only: general_spec_multwgt + use gsi_bundlemod, only: gsi_bundle + use general_sub2grid_mod, only: general_sub2grid,general_grid2sub + use general_specmod, only: spec_vars + use general_sub2grid_mod, only: sub2grid_info + use mpimod, only: mpi_comm_world,mype,npe,ierror + use file_utility, only : get_lun + implicit none + +! Declare passed variables + type(gsi_bundle),intent(in) :: wbundle + type(gsi_bundle),intent(inout) :: wbundle2 + type(spec_vars),intent (in):: sp_in + type(sub2grid_info),intent(in)::grd_in + real(r_kind),dimension(0:sp_in%jcap),intent(in):: spwgts + integer(i_kind), intent(in):: igrp !group index + integer(i_kind), intent(in):: nensid !ensemble member index + +! Declare local variables + integer(i_kind) ii,iflg,kk + integer(i_kind) i,j,lunit + + real(r_kind),dimension(grd_in%lat2,grd_in%lon2):: slndt,sicet,sst + real(r_kind),dimension(grd_in%nlat*grd_in%nlon*grd_in%nlevs_alloc) :: hwork + real(r_kind),dimension(grd_in%nlat,grd_in%nlon,grd_in%nlevs_alloc) :: work + real(r_kind),dimension(grd_in%nlat,grd_in%nlon,grd_in%nlevs_alloc) :: work1 + real(r_kind),dimension(sp_in%nc):: spc1 + real(r_single) outwork(grd_in%nlon,grd_in%nlat) + character*64 :: fname1 + character*5:: varname1 + + iflg=1 +! Beta1 first +! Get from subdomains to + call general_sub2grid(grd_in,wbundle%values,hwork) + work=reshape(hwork,(/grd_in%nlat,grd_in%nlon,grd_in%nlevs_alloc/)) + work1=work +! Transfer work to spectral space + do kk=1,grd_in%nlevs_alloc +! nvar_id is not supposed to be used for ensemble-related structures fr !being now +! Shut off this bit if univariate ozone or surface temperature +! if ( (nrf_var(nvar_id(kk))=='oz').or.(nrf_var(nvar_id(kk))=='OZ').and.oz_univ_static ) then +! cycle +! elseif ( (nrf_var(nvar_id(kk))=='sst').or.(nrf_var(nvar_id(kk))=='SST') ) then +! cycle +! elseif ( (nrf_var(nvar_id(kk))=='stl').or.(nrf_var(nvar_id(kk))=='STL') ) then +! cycle +! elseif ( (nrf_var(nvar_id(kk))=='sti').or.(nrf_var(nvar_id(kk))=='STI') ) then +! cycle +! end if +! if(mype==0) then + do j=1,grd_in%nlon + do i=1,grd_in%nlat + outwork(j,i)=work(i,j,kk) + enddo + enddo +! i=5-len_trim( nrf_var(nvar_id(kk))) +! varname1=repeat("_",i) +! varname1=trim(varname1)//trim(nrf_var(nvar_id(kk))) +! write(fname1,'("grp",i2.2,"mem",i2.2,"-out_vname",A5,"lev-",i3.3)')igrp,nensid, varname1,levs_id(kk) + + +! Transform from physical space to spectral space + call general_g2s0(grd_in,sp_in,spc1,work(:,:,kk)) + +! Apply spectral weights + call general_spec_multwgt(sp_in,spc1,spwgts) +! Transform back to physical space + call general_s2g0(grd_in,sp_in,spc1,work(:,:,kk)) + work1(:,:,kk)=work(:,:,kk)-work1(:,:,kk) +! if(mype==0) then +! do j=1,grd_in%nlon +! do i=1,grd_in%nlat +! outwork(j,i)=work(i,j,kk) +! enddo +! enddo + + end do + +! Transfer work back to subdomains + hwork=reshape(work,(/grd_in%nlat*grd_in%nlon*grd_in%nlevs_alloc/)) + call general_grid2sub(grd_in,hwork,wbundle2%values) + + return +end subroutine apply_scaledepwgts diff --git a/src/gsi/class_get_fv3_regional_ensperts.f90 b/src/gsi/class_get_fv3_regional_ensperts.f90 index 6c51a4fd8..bda4a47cd 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 703a62df3..beee90f5e 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 902b14ad2..b392bcd17 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 19170033e..9f496c7b0 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 a1c0a61ed..50fa81fae 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() @@ -445,10 +445,10 @@ subroutine allocate_cv(ycv) ! !$$$ end documentation block - use hybrid_ensemble_parameters, only: grd_ens + use hybrid_ensemble_parameters, only: grd_ens,nsclgrp implicit none type(control_vector), intent( out) :: ycv - integer(i_kind) :: ii,jj,nn,ndim,ierror,n_step,n_aens + integer(i_kind) :: ii,jj,nn,ndim,ierror,n_step,n_aens,ig character(len=256)::bname character(len=max_varname_length)::ltmp(1) type(gsi_grid) :: grid_motley @@ -481,7 +481,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,nsclgrp,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 @@ -522,21 +522,21 @@ subroutine allocate_cv(ycv) ! Set ensemble-based part of control vector if (l_hyb_ens) then - + do ig=1,nsclgrp ltmp(1)='a_en' do nn=1,n_ens - ycv%aens(jj,nn)%values => ycv%values(ii+1:ii+n_aens) + 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,nn),ycv%grid_aens,bname,ierror,names3d=ltmp,bundle_kind=r_kind) + 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,nn)%ndim + ndim=ndim+ycv%aens(jj,ig,nn)%ndim ii=ii+n_aens enddo - + end do endif enddo @@ -620,16 +620,20 @@ subroutine deallocate_cv(ycv) ! !$$$ end documentation block + use hybrid_ensemble_parameters, only: nsclgrp + implicit none type(control_vector), intent(inout) :: ycv - integer(i_kind) :: ii,nn,ierror + integer(i_kind) :: ii,nn,ierror,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) - enddo + do ig=1,nsclgrp + do nn=n_ens,1,-1 + call GSI_BundleUnset(ycv%aens(ii,ig,nn),ierror) + enddo + end do endif ! if (beta_s0>tiny_r_kind) then if(mvars>0) then @@ -845,9 +849,11 @@ real(r_quad) function qdot_prod_sub(xcv,ycv) ! !$$$ end documentation block + use hybrid_ensemble_parameters, only: nsclgrp + implicit none type(control_vector), intent(in ) :: xcv, ycv - integer(i_kind) :: ii,nn,m3d,m2d,i,j,itot + integer(i_kind) :: ii,nn,m3d,m2d,i,j,itot,ig,nigtmp real(r_quad),allocatable,dimension(:) :: partsum qdot_prod_sub=zero_quad @@ -858,18 +864,20 @@ 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(:) ) - end do - end do + do ig=1,nsclgrp + 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 else do ii=1,nsubwin 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+n_ens*nsclgrp allocate(partsum(itot)) !$omp parallel do schedule(dynamic,1) private(i) do i = 1,m3d @@ -880,10 +888,13 @@ 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 -!$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) - end do +!$omp parallel do schedule(dynamic,1) private(i,ig,nigtmp) + do ig=1,nsclgrp + nigtmp=n_ens*(ig-1) + do i = 1,n_ens + partsum(m3d+m2d+i+nigtmp) = 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 qdot_prod_sub = qdot_prod_sub + partsum(i) @@ -933,13 +944,15 @@ subroutine qdot_prod_vars_eb(xcv,ycv,prods,eb) ! !$$$ end documentation block + use hybrid_ensemble_parameters, only: nsclgrp + implicit none type(control_vector), intent(in ) :: xcv, ycv character(len=*) , intent(in ) :: eb real(r_quad) , intent( out) :: prods(nsubwin+1) real(r_quad) :: zz(nsubwin) - integer(i_kind) :: ii,i,nn,m3d,m2d + integer(i_kind) :: ii,i,nn,m3d,m2d,ig,ngtmp,nn0 real(r_quad),allocatable,dimension(:) :: partsum prods(:)=zero_quad @@ -953,11 +966,13 @@ 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(:) ) - end do - end do + do ig=1,nsclgrp + 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 else if(trim(eb) == 'cost_b') then @@ -981,20 +996,24 @@ 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*nsclgrp)) + do ig=1,nsclgrp + 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 + nn0=nn+ngtmp + partsum(nn0) = zero_quad + m3d=xcv%aens(ii,ig,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) + 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,nn)%n2d + m2d=xcv%aens(ii,ig,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) + 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 - do nn=1,n_ens + end do + do nn=1,n_ens*nsclgrp zz(ii)=zz(ii)+partsum(nn) end do deallocate(partsum) @@ -1348,11 +1367,13 @@ subroutine random_cv(ycv,kseed) ! !$$$ end documentation block +use hybrid_ensemble_parameters, only: nsclgrp + implicit none type(control_vector) , intent(inout) :: ycv integer(i_kind), optional, intent(in ) :: kseed -integer(i_kind):: ii,jj,nn,iseed +integer(i_kind):: ii,jj,nn,iseed,ig integer, allocatable :: nseed(:) ! Intentionaly default integer real(r_kind), allocatable :: zz(:) @@ -1381,12 +1402,14 @@ subroutine random_cv(ycv,kseed) if (nval_lenz_en>0) then allocate(zz(nval_lenz_en)) do nn=1,n_ens + do ig=1,nsclgrp do jj=1,nsubwin call random_number(zz) do ii=1,nval_lenz_en - ycv%aens(jj,nn)%values(ii) = two*zz(ii)-one + ycv%aens(jj,ig,nn)%values(ii) = two*zz(ii)-one enddo enddo + end do enddo deallocate(zz) endif diff --git a/src/gsi/cplr_get_fv3_regional_ensperts.f90 b/src/gsi/cplr_get_fv3_regional_ensperts.f90 index d64f26602..ebca5e52f 100644 --- a/src/gsi/cplr_get_fv3_regional_ensperts.f90 +++ b/src/gsi/cplr_get_fv3_regional_ensperts.f90 @@ -8,12 +8,15 @@ module get_fv3_regional_ensperts_mod procedure, pass(this) :: get_fv3_regional_ensperts => get_fv3_regional_ensperts_run procedure, pass(this) :: ens_spread_dualres_regional => ens_spread_dualres_regional_fv3_regional procedure, pass(this) :: general_read_fv3_regional + procedure, nopass :: GFilter + procedure, nopass :: scale_decomposition end type get_fv3_regional_ensperts_class - type(sub2grid_info):: grd_fv3lam_ens_dynvar_io_nouv,grd_fv3lam_ens_tracer_io_nouv,grd_fv3lam_ens_uv + type(sub2grid_info):: grd_fv3lam_ens_dynvar_io_nouv,grd_fv3lam_ens_tracer_io_nouv,grd_fv3lam_ens_phyvar_io_nouv,grd_fv3lam_ens_uv character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_dynmetvars3d_nouv ! copy of cvars3d excluding uv 3-d fields character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_tracermetvars3d_nouv ! copy of cvars3d excluding uv 3-d fields + character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_phymetvars3d_nouv character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_dynmetvars2d_nouv character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_ens_io_tracermetvars2d_nouv contains @@ -34,6 +37,11 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) ! ! 2021-08-10 lei - modify for fv3-lam ensemble spread output ! 2021-11-01 lei - modify for fv3-lam parallel IO + ! + ! 2022-08-29 Y. Wang, X. Wang - add scale decomposition (Gaussian smoothing) for multiscale DA using + ! scale- and variable-dependent localization, + ! Wang and Wang (2022ab, MWR) + ! poc: xuguang.wang@ou.edu ! input argument list: ! ! output argument list: @@ -65,27 +73,34 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) use directDA_radaruse_mod, only: l_cvpnr use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info use gridmod,only: regional - use gsi_rfv3io_mod, only: fv3lam_io_dynmetvars3d_nouv,fv3lam_io_tracermetvars3d_nouv + use gsi_rfv3io_mod, only: fv3lam_io_dynmetvars3d_nouv,fv3lam_io_tracermetvars3d_nouv, fv3lam_io_phymetvars3d_nouv use gsi_rfv3io_mod, only: fv3lam_io_dynmetvars2d_nouv,fv3lam_io_tracermetvars2d_nouv use netcdf , only: nf90_open, nf90_close,nf90_nowrite,nf90_inquire,nf90_format_netcdf4 use netcdf_mod , only: nc_check + use hybrid_ensemble_parameters, only: nsclgrp + use obsmod,only: if_model_dbz 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 real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig):: u,v,tv,oz,rh real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2):: ps - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)::w,ql,qi,qr,qg,qs,qnr + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)::w,ql,qi,qr,qg,qs,qnr,dbz real(r_single),pointer,dimension(:,:,:):: w3 =>NULL() real(r_single),pointer,dimension(:,:):: w2 =>NULL() real(r_kind),pointer,dimension(:,:,:):: x3 =>NULL() real(r_kind),pointer,dimension(:,:):: x2 =>NULL() type(gsi_bundle),allocatable,dimension(:):: en_bar + + type(gsi_bundle),allocatable :: en_pertstmp2(:) + type(gsi_bundle) :: en_pertstmp1 + integer(i_kind) :: ig0,ig + type(gsi_grid):: grid_ens real(r_kind):: bar_norm,sig_norm,kapr,kap1 @@ -96,7 +111,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) integer(i_kind):: i,j,k,n,mm1,istatus integer(i_kind):: ndynvario2d,ntracerio2d - integer(r_kind):: ndynvario3d,ntracerio3d + integer(r_kind):: ndynvario3d,ntracerio3d,nphyvario3d integer(i_kind):: inner_vars,numfields integer(i_kind):: ilev,ic2,ic3 integer(i_kind):: m @@ -119,6 +134,16 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) endif call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig) + + if( nsclgrp > 1 )then + allocate(en_pertstmp2(nsclgrp)) + + call gsi_bundlecreate(en_pertstmp1,grid_ens,'ensemble1',istatus,names2d=cvars2d,names3d=cvars3d,bundle_kind=r_single) + do ig=1,nsclgrp + call gsi_bundlecreate(en_pertstmp2(ig),grid_ens,'ensemble2',istatus,names2d=cvars2d,names3d=cvars3d,bundle_kind=r_single) + end do + end if + ! Allocate bundle to hold mean of ensemble members allocate(en_bar(ntlevs_ens)) @@ -127,6 +152,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) ntracerio2d=0 ndynvario3d=size(fv3lam_io_dynmetvars3d_nouv) ntracerio3d=size(fv3lam_io_tracermetvars3d_nouv) + if(if_model_dbz) nphyvario3d=size(fv3lam_io_phymetvars3d_nouv) if (allocated(fv3lam_io_dynmetvars2d_nouv))then ndynvario2d=size(fv3lam_io_dynmetvars2d_nouv) endif @@ -137,13 +163,14 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) allocate(fv3lam_ens_io_tracermetvars3d_nouv(ndynvario3d)) fv3lam_ens_io_dynmetvars3d_nouv=fv3lam_io_dynmetvars3d_nouv fv3lam_ens_io_tracermetvars3d_nouv=fv3lam_io_tracermetvars3d_nouv + if(if_model_dbz) fv3lam_ens_io_phymetvars3d_nouv=fv3lam_io_phymetvars3d_nouv if (ndynvario2d > 0 ) then allocate(fv3lam_ens_io_dynmetvars2d_nouv(ndynvario2d)) fv3lam_ens_io_dynmetvars2d_nouv=fv3lam_io_dynmetvars2d_nouv endif if (ntracerio2d > 0 ) then allocate(fv3lam_ens_io_tracermetvars2d_nouv(ntracerio2d)) - fv3lam_ens_io_tracermetvars2d_nouv =fv3lam_io_tracermetvars3d_nouv + fv3lam_ens_io_tracermetvars2d_nouv =fv3lam_io_tracermetvars2d_nouv endif @@ -191,7 +218,23 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) grd_ens%nlon,grd_ens%nsig,numfields,regional,names=names,lnames=lnames) + if( if_model_dbz )then + inner_vars=1 + numfields=inner_vars*(nphyvario3d*grd_ens%nsig+ntracerio2d) + deallocate(lnames,names) + allocate(lnames(1,numfields),names(1,numfields)) + ilev=1 + do i=1,nphyvario3d + do k=1,grd_ens%nsig + lnames(1,ilev)=k + names(1,ilev)=fv3lam_ens_io_phymetvars3d_nouv(i) + ilev=ilev+1 + enddo + enddo + call general_sub2grid_create_info(grd_fv3lam_ens_phyvar_io_nouv,inner_vars,grd_ens%nlat,& + grd_ens%nlon,grd_ens%nsig,numfields,regional,names=names,lnames=lnames) + end if numfields=grd_ens%nsig @@ -239,6 +282,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) end if + ig0 = 1 do m=1,ntlevs_ens @@ -248,7 +292,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,ig0,m)%valuesr4 = zero enddo mm1=mype+1 @@ -265,6 +309,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) fv3_filename%ak_bk=trim(ensfilenam_str)//'-fv3_akbk' fv3_filename%dynvars=trim(ensfilenam_str)//'-fv3_dynvars' fv3_filename%tracers=trim(ensfilenam_str)//"-fv3_tracer" + fv3_filename%phyvars=trim(ensfilenam_str)//'-fv3_phyvars' fv3_filename%sfcdata=trim(ensfilenam_str)//"-fv3_sfcdata" fv3_filename%couplerres=trim(ensfilenam_str)//"-coupler.res" @@ -301,17 +346,22 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) ! READ ENEMBLE MEMBERS DATA if (mype == 0) write(6,'(a,a)') & 'CALL READ_FV3_REGIONAL_ENSPERTS FOR ENS DATA with the filename str : ',trim(ensfilenam_str) - if (.not.l_use_dbz_directDA ) then ! Read additional hydrometers and w for dirZDA + if ( (.not.l_use_dbz_directDA) .and. (.not. if_model_dbz) ) then ! Read additional hydrometers and w for dirZDA call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz) else - call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz, & - g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_qnr=qnr,g_w=w) + if( l_use_dbz_directDA ) then + call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz, & + g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_qnr=qnr,g_w=w) + else if( if_model_dbz )then + call this%general_read_fv3_regional(fv3_filename,ps,u,v,tv,rh,oz, & + g_ql=ql,g_qi=qi,g_qr=qr,g_qs=qs,g_qg=qg,g_qnr=qnr,g_w=w,g_dbz=dbz) + end if end if ! 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,ig0,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) @@ -463,6 +513,16 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) end do end do end do + + case('dbz','DBZ') + do k=1,grd_ens%nsig + do i=1,grd_ens%lon2 + do j=1,grd_ens%lat2 + w3(j,i,k) = dbz(j,i,k) + x3(j,i,k)=x3(j,i,k)+dbz(j,i,k) + end do + end do + end do end select @@ -471,7 +531,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,ig0,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 +600,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,ig0,m)%valuesr4(i)=(en_perts(n,ig0,m)%valuesr4(i)-en_bar(m)%values(i))*sig_norm end do end do @@ -561,6 +621,27 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) deallocate(en_bar) ! + + + if( nsclgrp > 1 )then + do m=1,ntlevs_ens + do n=1,n_ens + if(mype==0) print*,"Diffusion operator for mem ",n + en_pertstmp1%valuesr4 = en_perts(n,ig0,m)%valuesr4 + call scale_decomposition(grd_ens,en_pertstmp1,en_pertstmp2,nsclgrp) + do ig=1,nsclgrp + en_perts(n,ig,m)%valuesr4 = en_pertstmp2(ig)%valuesr4 + end do + end do + end do + endif + + call gsi_bundledestroy(en_pertstmp1,istatus) + if( nsclgrp > 1 )then + do ig=1,nsclgrp + call gsi_bundledestroy(en_pertstmp2(ig),istatus) + end do + end if return @@ -572,7 +653,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar) end subroutine get_fv3_regional_ensperts_run subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g_rh,g_oz, & - g_ql,g_qi,g_qr,g_qs,g_qg,g_qnr,g_w) + g_ql,g_qi,g_qr,g_qs,g_qg,g_qnr,g_w,g_dbz) !$$$ subprogram documentation block ! first compied from general_read_arw_regional . . . . ! subprogram: general_read_fv3_regional read fv3sar model ensemble members @@ -623,6 +704,7 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g use hybrid_ensemble_parameters, only: grd_ens use directDA_radaruse_mod, only: l_use_cvpqx, cvpqx_pval, cld_nt_updt use directDA_radaruse_mod, only: l_cvpnr, cvpnr_pval + use obsmod, only:if_model_dbz @@ -632,7 +714,7 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g class(get_fv3_regional_ensperts_class), intent(inout) :: this type (type_fv3regfilenameg) , intent (in) :: fv3_filenameginput real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),intent(out)::g_u,g_v,g_tv,g_rh,g_oz - real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),optional,intent(out)::g_ql,g_qi,g_qr + real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),optional,intent(out)::g_ql,g_qi,g_qr,g_dbz real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig),optional,intent(out)::g_qs,g_qg,g_qnr,g_w real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2),intent(out):: g_ps @@ -656,11 +738,13 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g character(len=24),parameter :: myname_ = 'general_read_fv3_regional' type(gsi_bundle) :: gsibundle_fv3lam_ens_dynvar_nouv type(gsi_bundle) :: gsibundle_fv3lam_ens_tracer_nouv + type(gsi_bundle) :: gsibundle_fv3lam_ens_phyvar_nouv type(gsi_grid):: grid_ens character(len=:),allocatable :: grid_spec !='fv3_grid_spec' character(len=:),allocatable :: ak_bk !='fv3_akbk' character(len=:),allocatable :: dynvars !='fv3_dynvars' + character(len=:),allocatable :: phyvars !='fv3_phyvars' character(len=:),allocatable :: tracers !='fv3_tracer' character(len=:),allocatable :: sfcdata !='fv3_sfcdata' character(len=:),allocatable :: couplerres!='coupler.res' @@ -677,6 +761,7 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g ak_bk=fv3_filenameginput%ak_bk dynvars=fv3_filenameginput%dynvars tracers=fv3_filenameginput%tracers + if(if_model_dbz) phyvars=fv3_filenameginput%phyvars sfcdata=fv3_filenameginput%sfcdata couplerres=fv3_filenameginput%couplerres @@ -701,6 +786,10 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g endif + if( if_model_dbz )then + call gsi_bundlecreate(gsibundle_fv3lam_ens_phyvar_nouv,grid_ens,'gsibundle_fv3lam_ens_phyvar_nouv',istatus, & + names3d=fv3lam_ens_io_phymetvars3d_nouv) + end if @@ -714,6 +803,9 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g fv3_filenameginput%dynvars,fv3_filenameginput) call gsi_fv3ncdf_read(grd_fv3lam_ens_tracer_io_nouv,gsibundle_fv3lam_ens_tracer_nouv,& fv3_filenameginput%tracers,fv3_filenameginput) + if( if_model_dbz )& + call gsi_fv3ncdf_read(grd_fv3lam_ens_phyvar_io_nouv,gsibundle_fv3lam_ens_phyvar_nouv,& + fv3_filenameginput%phyvars,fv3_filenameginput) else call gsi_fv3ncdf_read_v1(grd_fv3lam_ens_dynvar_io_nouv,gsibundle_fv3lam_ens_dynvar_nouv,& fv3_filenameginput%dynvars,fv3_filenameginput) @@ -724,14 +816,17 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'tsen' ,g_tsen ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'q' ,g_q ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'oz' ,g_oz ,istatus );ier=ier+istatus - if (l_use_dbz_directDA) then + if (l_use_dbz_directDA .or. if_model_dbz ) then call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'ql' ,g_ql ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qi' ,g_qi ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qr' ,g_qr ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qs' ,g_qs ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qg' ,g_qg ,istatus );ier=ier+istatus - call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qnr',g_qnr ,istatus );ier=ier+istatus + if (l_use_dbz_directDA) & + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_tracer_nouv, 'qnr',g_qnr ,istatus );ier=ier+istatus call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_dynvar_nouv, 'w' , g_w ,istatus );ier=ier+istatus + if( if_model_dbz )& + call GSI_Bundlegetvar ( gsibundle_fv3lam_ens_phyvar_nouv, 'dbz' , g_dbz ,istatus );ier=ier+istatus end if @@ -834,6 +929,8 @@ subroutine general_read_fv3_regional(this,fv3_filenameginput,g_ps,g_u,g_v,g_tv,g enddo call gsi_bundledestroy(gsibundle_fv3lam_ens_dynvar_nouv) call gsi_bundledestroy(gsibundle_fv3lam_ens_tracer_nouv) + + if(if_model_dbz)call gsi_bundledestroy(gsibundle_fv3lam_ens_phyvar_nouv) return @@ -880,14 +977,14 @@ 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 type(gsi_grid):: grid_ens,grid_anl real(r_kind) sp_norm - integer(i_kind) i,n + integer(i_kind) i,n,ig0 integer(i_kind) istatus character(255) enspreadfilename @@ -917,12 +1014,13 @@ subroutine ens_spread_dualres_regional_fv3_regional(this,mype,en_perts,nelen) sp_norm=(one/float(n_ens)) sube%values=zero + ig0=1 ! 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,ig0,1)%valuesr4(i))*(en_perts(n,ig0,1)%valuesr4(i)) end do end do @@ -934,4 +1032,236 @@ subroutine ens_spread_dualres_regional_fv3_regional(this,mype,en_perts,nelen) return end subroutine ens_spread_dualres_regional_fv3_regional + + subroutine scale_decomposition( grd_in, en_p, en_out, nscale ) + ! + ! 2022-08-29 Y. Wang, X. Wang - Scale decomposition for multiscale DA using SDL and VDL + ! poc: xuguang.wang@ou.edu + ! + use hybrid_ensemble_parameters, only: smooth_scales,vdl_scale,small_loc_varlist, smooth_scales_num + use kinds, only: r_kind,i_kind,r_single + use general_sub2grid_mod, only: general_sub2grid,general_grid2sub,general_sub2grid_create_info, & + sub2grid_info,general_sub2grid_destroy_info + use gsi_bundlemod, only: gsi_bundlegetpointer,gsi_bundlemerge,gsi_bundle,gsi_bundledup,gsi_bundledestroy + use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d,nrf_var + use gridmod,only: nsig + use mpimod, only: mype,nvar_id + + implicit none + ! + integer(i_kind) :: scale_num,k,i,j,niter,n_ens, nscale + integer(i_kind) :: lon2,lat2,nvert,kend_loc,ic3,istatus,iscl,ic2 + integer(i_kind) :: nlon, nlat, kbegin_loc, kend_alloc,ii,ie,is + + type(sub2grid_info),intent(in) :: grd_in + type(gsi_bundle), intent(in) :: en_p + type(gsi_bundle), intent(inout) :: en_out(nscale) + + real(r_single),dimension(grd_in%nlat*grd_in%nlon*grd_in%nlevs_alloc) :: hwork + real(r_single),dimension(grd_in%nlat,grd_in%nlon,grd_in%nlevs_alloc) :: work + real(r_single),dimension(grd_in%nlat,grd_in%nlon,grd_in%nlevs_alloc) :: tmpwork + real(r_single),dimension(grd_in%nlat,grd_in%nlon,grd_in%nlevs_alloc) :: work_out + real(r_single),dimension(2,grd_in%nlat,grd_in%nlon,grd_in%nlevs_alloc) :: work_vdl + real(r_single),dimension(2,grd_in%nlat,grd_in%nlon,grd_in%nlevs_alloc) :: work_vdl_last + real(r_single),allocatable,dimension(:,:) :: kernal + + logical :: maskflg + integer(i_kind) :: nn, ig, ivdl, vdl_group + real(r_single) :: sigma + + hwork=0.0_r_single + work=0.0_r_single + work_out=0.0_r_single + tmpwork = 0.0_r_single + + nlat = grd_in%nlat + nlon = grd_in%nlon + + call general_sub2grid(grd_in,en_p%valuesr4,hwork) + work=reshape(hwork,(/grd_in%nlat,grd_in%nlon,grd_in%nlevs_alloc/)) + + ig = 1 + + do iscl = 1, smooth_scales_num + do k = 1,grd_in%nlevs_alloc + + if( any( trim(nrf_var(nvar_id(k))) == small_loc_varlist ) ) then + maskflg = .true. + else + maskflg = .false. + end if + + nn = 2 * int(smooth_scales(iscl)) + 1 + sigma = real( smooth_scales(iscl) / 3.0 ) + if(allocated(kernal))deallocate(kernal) + allocate(kernal(nn,nn)) + call Gaussian_kernal(nn, sigma, kernal) + call GFilter( work(:,:,k), work_out(:,:,k), kernal, nlat,nlon, int(smooth_scales(iscl)) ) + tmpwork(:,:,k) = work(:,:,k) - work_out(:,:,k) + + if( vdl_scale(ig) == 1 )then + if( maskflg ) then + work_vdl(vdl_scale(ig),:,:,k) = tmpwork(:,:,k) + work_vdl(vdl_scale(ig+1),:,:,k) = 0.0_r_kind + else + work_vdl(vdl_scale(ig),:,:,k) = 0.0_r_kind + work_vdl(vdl_scale(ig+1),:,:,k) = tmpwork(:,:,k) + end if + + else if ( vdl_scale(ig) == 0 )then + work_vdl(1,:,:,k) = tmpwork(:,:,k) + end if + + + + work(:,:,k) = work_out(:,:,k) + end do + + if( vdl_scale(ig) == 1 )then + vdl_group = 2 + else if ( vdl_scale(ig) == 0 )then + vdl_group = 1 + end if + + do ivdl = 1, vdl_group ! we separate variables into 2 groups in this version for variable-dependent localization + hwork=reshape(work_vdl(ivdl,:,:,:),(/grd_in%nlat*grd_in%nlon*grd_in%nlevs_alloc/)) + call general_grid2sub(grd_in,hwork,en_out(ig)%valuesr4) + ig = ig + 1 + end do + + if( iscl == smooth_scales_num )then + do k = 1,grd_in%nlevs_alloc + + if( any( trim(nrf_var(nvar_id(k))) == small_loc_varlist ) ) then + maskflg = .true. + else + maskflg = .false. + end if + + if( vdl_scale(ig) == 1 )then + if( maskflg ) then + work_vdl(vdl_scale(ig),:,:,k) = work_out(:,:,k) + work_vdl(vdl_scale(ig+1),:,:,k) = 0.0_r_kind + else + work_vdl(vdl_scale(ig),:,:,k) = 0.0_r_kind + work_vdl(vdl_scale(ig+1),:,:,k) = work_out(:,:,k) + end if + + else if ( vdl_scale(ig) == 0 )then + work_vdl(1,:,:,k) = work_out(:,:,k) + end if + end do + end if + + + end do ! iscl + + if( vdl_scale(ig) == 1 )then + vdl_group = 2 + else if ( vdl_scale(ig) == 0 )then + vdl_group = 1 + end if + + do ivdl = 1, vdl_group + hwork=reshape(work_vdl(ivdl,:,:,:),(/grd_in%nlat*grd_in%nlon*grd_in%nlevs_alloc/)) + call general_grid2sub(grd_in,hwork,en_out(ig)%valuesr4) + ig = ig + 1 + end do + + if( mype == 0 )then + if( ig - 1 /= nscale ) print*,'ig does not match the number of decomposed scales',ig,nscale + end if + + end subroutine scale_decomposition + + SUBROUTINE GFilter( field, fieldout, kernal, nx, ny, nn ) + ! + ! 2022-08-29 Y. Wang, X. Wang - Gaussian filter for scale decomposition + ! poc: xuguang.wang@ou.edu + + implicit none + ! + + integer, intent(in) :: nx,ny,nn + + real,dimension(nx,ny), intent(in) :: field + real,dimension(2*nn+1,2*nn+1), intent(in) :: kernal + real,dimension(nx,ny), intent(out) :: fieldout + + ! Local declare + real,parameter :: pi = 3.14159 + integer :: ix,iy,ii,jj + real :: summ + real,dimension(nx+2*nn,ny+2*nn) :: field_ext + + ! Fill the center + field_ext(nn+1:nx+nn,nn+1:ny+nn) = field + + ! Extend bottom side + do ix = 1, nn + field_ext(ix,nn+1:ny+nn) = field_ext(2*nn+1-ix,nn+1:ny+nn) + end do + + ! Extend top side + do ix = 1, nn + field_ext(nx+nn+ix,nn+1:ny+nn) = field_ext(nx+nn+1-ix,nn+1:ny+nn) + end do + + ! Extend left + do iy = 1, nn + field_ext(:,iy) = field_ext(:,2*nn+1-iy) + end do + + ! Extend right + do iy = 1, nn + field_ext(:,ny+nn+iy) = field_ext(:,ny+nn+1-iy) + end do + + do ix = nn + 1, nx + nn + do iy = nn + 1, ny + nn + summ = 0.0 + do ii = 1, 2*nn + 1 + do jj = 1, 2*nn + 1 + summ = summ + field_ext(ix-nn+ii-1,iy-nn+jj-1)*kernal(ii,jj) + end do + end do + fieldout(ix-nn,iy-nn) = summ + end do + end do + + END SUBROUTINE GFilter + + SUBROUTINE Gaussian_kernal(nn, sigma, kernal) + ! + ! 2022-08-29 Y. Wang, X. Wang - Calculate kernal for Gaussian filter + ! poc: xuguang.wang@ou.edu + + + implicit none + integer, intent(in) :: nn + real, intent(in) :: sigma + real,dimension(nn, nn), intent(out) :: kernal + real,parameter :: pi = 3.141593 + + ! Local declare + integer :: i, j, k + real :: norm, summ + + k = int(nn/2) + 1 + + summ = 0.0 + do i = 1, nn + do j = 1, nn + norm = real(i-k)**2.0 + real(j-k)**2.0 + kernal(i,j) = exp(-1.0*norm/(2*sigma**2.0))/sqrt(2.0*pi*sigma**2.0) + summ = summ + kernal(i,j) + end do + end do + + kernal = kernal / summ + + END SUBROUTINE Gaussian_kernal + + + end module get_fv3_regional_ensperts_mod diff --git a/src/gsi/cplr_get_pseudo_ensperts.f90 b/src/gsi/cplr_get_pseudo_ensperts.f90 index 8c59cd0be..9d1216540 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 eb94c6c0f..775c9c419 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 9de32ea38..3e3082263 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 582574b52..5c9934919 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 1a9dc58a5..b9d7db36e 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) @@ -312,7 +312,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 +324,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 afab9737a..da70d31ab 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 525bf33af..35cc03444 100644 --- a/src/gsi/ensctl2model.f90 +++ b/src/gsi/ensctl2model.f90 @@ -26,7 +26,7 @@ subroutine ensctl2model(xhat,mval,eval) use control_vectors, only: control_vector,cvars3d use gsi_4dvar, only: ibin_anl use hybrid_ensemble_parameters, only: uv_hyb_ens,dual_res,ntlevs_ens -use hybrid_ensemble_parameters, only: n_ens,q_hyb_ens +use hybrid_ensemble_parameters, only: n_ens,q_hyb_ens,nsclgrp use hybrid_ensemble_isotropic, only: ensemble_forward_model,ensemble_forward_model_dual_res use hybrid_ensemble_isotropic, only: sqrt_beta_s_mult,sqrt_beta_e_mult, & ckgcov_a_en_new_factorization @@ -56,9 +56,9 @@ subroutine ensctl2model(xhat,mval,eval) integer(i_kind) :: jj,ic,id,istatus,nclouds,nn integer(i_kind), parameter :: ncvars = 5 -integer(i_kind) :: icps(ncvars) +integer(i_kind) :: icps(ncvars),ig 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 '/) @@ -117,17 +117,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 - enddo + allocate(ebundle(nsclgrp,n_ens)) + do ig=1,nsclgrp + 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 + end do ! Apply square-root of ensemble error covariance - call ckgcov_a_en_new_factorization(xhat%aens(jj,1)%values(:),ebundle) + do ig=1,nsclgrp + call ckgcov_a_en_new_factorization(ig,xhat%aens(jj,ig,1)%values(:),ebundle(ig,:)) + end do call sqrt_beta_e_mult(ebundle) ! Initialize ensemble contribution to zero @@ -230,13 +234,15 @@ 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 - enddo + do ig=1,nsclgrp + 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 + end do deallocate(ebundle) end do ! ntlevs diff --git a/src/gsi/ensctl2model_ad.f90 b/src/gsi/ensctl2model_ad.f90 index 8c491b6de..13caa8bb4 100644 --- a/src/gsi/ensctl2model_ad.f90 +++ b/src/gsi/ensctl2model_ad.f90 @@ -23,7 +23,7 @@ subroutine ensctl2model_ad(eval,mval,grad) use kinds, only: r_kind,i_kind use control_vectors, only: control_vector,cvars3d use gsi_4dvar, only: ibin_anl -use hybrid_ensemble_parameters, only: uv_hyb_ens,dual_res,nval_lenz_en,ntlevs_ens,n_ens,q_hyb_ens +use hybrid_ensemble_parameters, only: uv_hyb_ens,dual_res,nval_lenz_en,ntlevs_ens,n_ens,q_hyb_ens,nsclgrp use hybrid_ensemble_isotropic, only: ensemble_forward_model_ad use hybrid_ensemble_isotropic, only: ckgcov_a_en_new_factorization_ad use hybrid_ensemble_isotropic, only: ensemble_forward_model_ad_dual_res @@ -55,9 +55,9 @@ subroutine ensctl2model_ad(eval,mval,grad) integer(i_kind) :: ii,jj,ic,id,istatus,nclouds,nn integer(i_kind), parameter :: ncvars = 5 -integer(i_kind) :: icps(ncvars) +integer(i_kind) :: icps(ncvars),ig 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 ', & @@ -122,13 +122,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(nsclgrp,n_ens)) + do ig=1,nsclgrp + 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 +212,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,nsclgrp + 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 +226,9 @@ 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,nsclgrp + call ckgcov_a_en_new_factorization_ad(ig,grade,ebundle(ig,:)) + end do call gsi_bundledestroy(wbundle_c,istatus) if (istatus/=0) then @@ -230,19 +236,23 @@ 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 ig=1,nsclgrp + do ii=1,n_ens + grad%aens(jj,ig,ii)%values=grad%aens(jj,ig,ii)%values + grade + enddo + end do ! 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 + do ig=1,nsclgrp + 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', nn, istatus + call stop2(999) + endif + enddo enddo deallocate(ebundle) ! deallocate(grade) diff --git a/src/gsi/ensctl2state.f90 b/src/gsi/ensctl2state.f90 index 4214b01d8..0d6d3042c 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 212d0312c..4c038c8c6 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/general_specmod.f90 b/src/gsi/general_specmod.f90 index 439e26e43..6bf9f212b 100644 --- a/src/gsi/general_specmod.f90 +++ b/src/gsi/general_specmod.f90 @@ -64,6 +64,7 @@ module general_specmod ! set subroutines to public public :: general_init_spec_vars public :: general_destroy_spec_vars + public :: general_spec_multwgt ! set passed variables to public public :: spec_vars public :: spec_cut @@ -306,6 +307,46 @@ subroutine general_init_spec_vars(sp,jcap,jcap_test,nlat_a,nlon_a,eqspace) return end subroutine general_init_spec_vars + subroutine general_spec_multwgt(sp,spcwrk,spcwgt) +! 2017-03-30 J. Kay, X. Wang - add general_spec_multwgt for scale-dependent +! weighting of mixed resolution ensemble, +! POC: xuguang.wang@ou.edu +! + implicit none + type(spec_vars),intent(in) :: sp + real(r_kind),dimension(sp%nc),intent(inout) :: spcwrk + real(r_kind),dimension(0:sp%jcap),intent(in) :: spcwgt + + integer(i_kind) ii1,l,m,jmax,ks,n + +! ii1=0 +! do l=0,sp%jcap +! do m=0,sp%jcap-l +! ii1=ii1+2 +! spcwrk(ii1-1)=spcwrk(ii1-1)*spcwgt(l) +! spcwrk(ii1)=spcwrk(ii1)*spcwgt(l) +! end do +! end do +!! Code borrowed from spvar in splib + JMAX=sp%jcap + + L=0 + DO N=0,JMAX + KS=2*N + spcwrk(KS+1)=spcwrk(KS+1)*spcwgt(N) + ENDDO + DO N=0,JMAX + DO L=MAX(1,N-JMAX),MIN(N,JMAX) + KS=L*(2*JMAX+(-1)*(L-1))+2*N + spcwrk(KS+1) = spcwrk(KS+1)*spcwgt(N) + spcwrk(KS+2) = spcwrk(KS+2)*spcwgt(N) + ENDDO + ENDDO + + + return + end subroutine general_spec_multwgt + subroutine general_destroy_spec_vars(sp) !$$$ subprogram documentation block ! . . . . diff --git a/src/gsi/get_gefs_ensperts_dualres.f90 b/src/gsi/get_gefs_ensperts_dualres.f90 index c547a9f06..af54758a0 100644 --- a/src/gsi/get_gefs_ensperts_dualres.f90 +++ b/src/gsi/get_gefs_ensperts_dualres.f90 @@ -68,6 +68,7 @@ subroutine get_gefs_ensperts_dualres use gsi_enscouplermod, only: gsi_enscoupler_create_sub2grid_info use gsi_enscouplermod, only: gsi_enscoupler_destroy_sub2grid_info use general_sub2grid_mod, only: sub2grid_info,general_sub2grid_create_info,general_sub2grid_destroy_info + use hybrid_ensemble_parameters, only: nsclgrp,sp_ens,spc_multwgt implicit none real(r_kind),pointer,dimension(:,:) :: ps @@ -77,7 +78,9 @@ subroutine get_gefs_ensperts_dualres real(r_kind),pointer,dimension(:,:,:):: p3 real(r_kind),pointer,dimension(:,:):: x2 type(gsi_bundle),allocatable,dimension(:) :: en_read - type(gsi_bundle):: en_bar + type(gsi_bundle) :: en_bar + type(gsi_bundle) :: en_pertstmp1 + type(gsi_bundle) :: en_pertstmp2 ! type(gsi_grid) :: grid_ens real(r_kind) bar_norm,sig_norm,kapr,kap1 ! real(r_kind),allocatable,dimension(:,:):: z,sst2 @@ -91,34 +94,35 @@ subroutine get_gefs_ensperts_dualres ! integer(i_kind) il,jl logical ice,hydrometeor type(sub2grid_info) :: grd_tmp + integer(i_kind) :: ig0,ig ! 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 +130,24 @@ 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 temporary usage + if( nsclgrp > 1 )then + call gsi_bundlecreate(en_pertstmp1,en_perts(1,1,1)%grid,'aux-ens-read',istatus,names2d=cvars2d,names3d=cvars3d) + call gsi_bundlecreate(en_pertstmp2,en_perts(1,1,1)%grid,'aux-ens-read',istatus,names2d=cvars2d,names3d=cvars3d) + if(istatus/=0) then + write(6,*)' get_gefs_ensperts_dualres: trouble creating en_read like tempbundle' + call stop2(999) + endif + end if + ! 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 @@ -144,10 +158,12 @@ subroutine get_gefs_ensperts_dualres ! sst2=zero ! for now, sst not used in ensemble perturbations, so if sst array is called for ! then sst part of en_perts will be zero when sst2=zero + ig0 = 1 + !$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,ig0,m)%valuesr4=zero_single end do end do @@ -254,7 +270,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,ig0,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 +311,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,ig0,m)%valuesr4(i)=en_perts(n,ig0,m)%valuesr4(i)-en_bar%values(i) end do if(.not. q_hyb_ens) then do ic3=1,nc3d @@ -304,8 +320,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,ig0,m)%r3(ipic)%qr4(i,j,k) = min(en_perts(n,ig0,m)%r3(ipic)%qr4(i,j,k),limqens) + en_perts(n,ig0,m)%r3(ipic)%qr4(i,j,k) = max(en_perts(n,ig0,m)%r3(ipic)%qr4(i,j,k),-limqens) end do end do end do @@ -313,11 +329,25 @@ 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,ig0,m)%valuesr4(i)=en_perts(n,ig0,m)%valuesr4(i)*sig_norm end do end do end do ntlevs_ens_loop !end do over bins + if(nsclgrp > 1) then + do m=1,ntlevs_ens + do n=1,n_ens + en_pertstmp1%values=en_perts(n,ig0,m)%valuesr4 + do ig=1,nsclgrp + call apply_scaledepwgts(grd_ens,sp_ens,en_pertstmp1,spc_multwgt(:,ig),en_pertstmp2,ig,n) + en_perts(n,ig,m)%valuesr4=en_pertstmp2%values + enddo + enddo + + enddo + + endif + do n=n_ens,1,-1 call gsi_bundledestroy(en_read(n),istatus) if ( istatus /= 0 ) & @@ -429,6 +459,7 @@ subroutine ens_spread_dualres(en_bar,ibin) logical regional integer(i_kind) num_fields,inner_vars,istatus logical,allocatable::vector(:) + integer(i_kind) :: ig0 ! create simple regular grid call gsi_gridcreate(grid_anl,grd_anl%lat2,grd_anl%lon2,grd_anl%nsig) @@ -450,12 +481,13 @@ subroutine ens_spread_dualres(en_bar,ibin) endif sp_norm=(one/float(n_ens)) + ig0 = 1 sube%values=zero 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,ig0,ibin)%valuesr4(i)-en_bar%values(i))*(en_perts(n,ig0,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 958849a17..a076f0ccf 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 93d23c837..ece1780c0 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/gsi_files.cmake b/src/gsi/gsi_files.cmake index 461b49ddf..1658c83bf 100644 --- a/src/gsi/gsi_files.cmake +++ b/src/gsi/gsi_files.cmake @@ -87,6 +87,7 @@ anisofilter_glb.f90 antcorr_application.f90 antest_maps0.f90 antest_maps0_glb.f90 +apply_scaledepwgts.f90 atms_spatial_average_mod.f90 balmod.f90 berror.f90 diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index 6f50d1386..34436c6c2 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -65,6 +65,7 @@ module gsi_rfv3io_mod character(len=:),allocatable :: ak_bk !='fv3_akbk' character(len=:),allocatable :: dynvars !='fv3_dynvars' character(len=:),allocatable :: tracers !='fv3_tracer' + character(len=:),allocatable :: phyvars !='fv3_phyvars' character(len=:),allocatable :: sfcdata !='fv3_sfcdata' character(len=:),allocatable :: couplerres!='coupler.res' contains @@ -87,8 +88,9 @@ module gsi_rfv3io_mod type(sub2grid_info) :: grd_fv3lam_dynvar_ionouv type(sub2grid_info) :: grd_fv3lam_tracer_ionouv type(sub2grid_info) :: grd_fv3lam_tracerchem_ionouv + type(sub2grid_info) :: grd_fv3lam_phyvar_ionouv type(sub2grid_info) :: grd_fv3lam_uv - integer(i_kind) ,parameter:: ndynvarslist=13, ntracerslist=8 + integer(i_kind) ,parameter:: ndynvarslist=13, ntracerslist=8, nphyvarslist=1 character(len=max_varname_length), dimension(ndynvarslist), parameter :: & vardynvars = [character(len=max_varname_length) :: & @@ -96,12 +98,14 @@ module gsi_rfv3io_mod character(len=max_varname_length), dimension(ntracerslist+naero_cmaq_fv3+7), parameter :: & vartracers = [character(len=max_varname_length) :: & 'q','oz','ql','qi','qr','qs','qg','qnr',aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk'] - character(len=max_varname_length), dimension(15+naero_cmaq_fv3+7), parameter :: & + character(len=max_varname_length), dimension(nphyvarslist), parameter :: & + varphyvars = [character(len=max_varname_length) :: 'dbz'] + character(len=max_varname_length), dimension(16+naero_cmaq_fv3+7), parameter :: & varfv3name = [character(len=max_varname_length) :: & - 'u','v','W','T','delp','sphum','o3mr','liq_wat','ice_wat','rainwat','snowwat','graupel','rain_nc','ps','DZ', & + 'u','v','W','T','delp','sphum','o3mr','liq_wat','ice_wat','rainwat','snowwat','graupel','rain_nc','ref_f3d','ps','DZ', & aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk'], & vgsiname = [character(len=max_varname_length) :: & - 'u','v','w','tsen','delp','q','oz','ql','qi','qr','qs','qg','qnr','ps','delzinc', & + 'u','v','w','tsen','delp','q','oz','ql','qi','qr','qs','qg','qnr','dbz','ps','delzinc', & aeronames_cmaq_fv3,'pm25at','pm25ac','pm25co','pm2_5','amassi','amassj','amassk'] character(len=max_varname_length),dimension(:),allocatable:: name_metvars2d character(len=max_varname_length),dimension(:),allocatable:: name_metvars3d @@ -126,6 +130,7 @@ module gsi_rfv3io_mod public :: k_snwdph,k_f10m,mype_2d,n2d,k_orog,k_psfc,k_t2m,k_q2m public :: ijns,ijns2d,displss,displss2d,ijnz,displsz_g public :: fv3lam_io_dynmetvars3d_nouv,fv3lam_io_tracermetvars3d_nouv,fv3lam_io_tracerchemvars3d_nouv + public :: fv3lam_io_phymetvars3d_nouv public :: fv3lam_io_dynmetvars2d_nouv,fv3lam_io_tracermetvars2d_nouv integer(i_kind) mype_u,mype_v,mype_t,mype_q,mype_p,mype_delz,mype_oz,mype_ql @@ -153,6 +158,7 @@ module gsi_rfv3io_mod ! copy of cvars3d excluding uv 3-d fields character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_tracermetvars3d_nouv ! copy of cvars3d excluding uv 3-d fields + character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_phymetvars3d_nouv character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_tracerchemvars3d_nouv ! copy of cvars3d excluding uv 3-d fields character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_io_dynmetvars2d_nouv @@ -163,8 +169,10 @@ module gsi_rfv3io_mod !to define names in gsibundle character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_names_gsibundle_tracer_nouv !to define names in gsibundle + character(len=max_varname_length),allocatable,dimension(:) :: fv3lam_names_gsibundle_phyvar_nouv type(gsi_bundle):: gsibundle_fv3lam_dynvar_nouv type(gsi_bundle):: gsibundle_fv3lam_tracer_nouv + type(gsi_bundle):: gsibundle_fv3lam_phyvar_nouv type(gsi_bundle):: gsibundle_fv3lam_tracerchem_nouv contains @@ -197,6 +205,12 @@ subroutine fv3regfilename_init(this,it) write(filename,"(A11,I2.2)") 'fv3_tracer_',ifilesig(it) this%tracers=trim(filename) endif + if (it == ntguessig) then + this%phyvars='fv3_phyvars' + else + write(filename,"(A12,I2.2)") 'fv3_phyvars_',ifilesig(it) + this%phyvars=trim(filename) + endif if (it == ntguessig) then this%sfcdata='fv3_sfcdata' else @@ -761,6 +775,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) use gsi_metguess_mod, only: gsi_metguess_get use netcdf, only:nf90_open,nf90_close,nf90_inquire,nf90_nowrite, nf90_format_netcdf4 use gsi_chemguess_mod, only: gsi_chemguess_get + use obsmod, only: if_model_dbz implicit none @@ -790,7 +805,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) real(r_kind),dimension(:,:,:),pointer::ges_qg=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_qnr=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_w=>NULL() - + real(r_kind),dimension(:,:,:),pointer::ges_dbz=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_aalj=>NULL() real(r_kind),dimension(:,:,:),pointer::ges_acaj=>NULL() @@ -880,8 +895,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) integer(i_kind),dimension(:,:),allocatable:: lnames integer(i_kind),dimension(:,:),allocatable:: uvlnames integer(i_kind):: inner_vars,numfields - integer(i_kind):: ndynvario2d,ntracerio2d,ilev,jdynvar,jtracer - integer(r_kind):: iuv,ndynvario3d,ntracerio3d,ntracerchemio3d + integer(i_kind):: ndynvario2d,ntracerio2d,ilev,jdynvar,jtracer,jphyvar + integer(r_kind):: iuv,ndynvario3d,ntracerio3d,nphyvario3d,ntracerchemio3d integer(i_kind):: loc_id,ncfmt !clt this block is still maintained for they would be needed for a certain 2d fields IO @@ -946,6 +961,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) iuv=0 ndynvario3d=0 ntracerio3d=0 + nphyvario3d=0 do i=1,size(name_metvars3d) vartem=trim(name_metvars3d(i)) if(trim(vartem)=='u'.or.trim(vartem)=='v') then @@ -956,6 +972,8 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) ndynvario3d=ndynvario3d+1 else if (ifindstrloc(vartracers,trim(vartem))> 0) then ntracerio3d=ntracerio3d+1 + else if (ifindstrloc(varphyvars,trim(vartem))> 0) then + nphyvario3d=nphyvario3d+1 else write(6,*)'the metvarname1 ',trim(vartem),' has not been considered yet, stop' call stop2(333) @@ -967,6 +985,12 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) write(6,*)"the set up for met variable is not as expected, abort" call stop2(222) endif + if ( if_model_dbz ) then + if( nphyvario3d<=0 ) then + write(6,*)"the set up for met variable (phyvar) is not as expected, abort" + call stop2(223) + end if + endif if (fv3sar_bg_opt == 0.and.ifindstrloc(name_metvars3d,'delp') <= 0)then ndynvario3d=ndynvario3d+1 ! for delp endif @@ -976,6 +1000,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) if (l_reg_update_hydro_delz.and.fv3sar_bg_opt==0) ndynvario3d=ndynvario3d+1 ! for delzinc allocate(fv3lam_io_dynmetvars3d_nouv(ndynvario3d)) allocate(fv3lam_io_tracermetvars3d_nouv(ntracerio3d)) + allocate(fv3lam_io_phymetvars3d_nouv(nphyvario3d)) if (laeroana_fv3cmaq) then allocate(fv3lam_io_tracerchemvars3d_nouv(naero_cmaq_fv3+7)) @@ -983,6 +1008,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) jdynvar=0 jtracer=0 + jphyvar=0 do i=1,size(name_metvars3d) vartem=trim(name_metvars3d(i)) if(.not.(trim(vartem)=='u'.or.trim(vartem)=='v'.or.trim(vartem)=='iqr')) then @@ -997,6 +1023,9 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) else if (ifindstrloc(vartracers,trim(vartem)) > 0) then jtracer=jtracer+1 fv3lam_io_tracermetvars3d_nouv(jtracer)=trim(vartem) + else if (ifindstrloc(varphyvars,trim(vartem)) > 0) then + jphyvar=jphyvar+1 + fv3lam_io_phymetvars3d_nouv(jphyvar)=trim(vartem) else write(6,*)'the metvarname ',vartem,' is not expected, stop' call flush(6) @@ -1012,7 +1041,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) jdynvar=jdynvar+1 fv3lam_io_dynmetvars3d_nouv(jdynvar)="delzinc" endif - if(jdynvar /= ndynvario3d.or.jtracer /= ntracerio3d ) then + if(jdynvar /= ndynvario3d.or.jtracer /= ntracerio3d.or.jphyvar /= nphyvario3d ) then write(6,*)'ndynvario3d is not as expected, stop' call flush(6) call stop2(333) @@ -1020,6 +1049,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) if(mype == 0) then write(6,*) ' fv3lam_io_dynmetvars3d_nouv is ',(trim(fv3lam_io_dynmetvars3d_nouv(i)),i=1,ndynvario3d) write(6,*) ' fv3lam_io_tracermevars3d_nouv is ',(trim(fv3lam_io_tracermetvars3d_nouv(i)),i=1,ntracerio3d) + write(6,*) ' fv3lam_io_phymetvars3d_nouv is ',(trim(fv3lam_io_phymetvars3d_nouv(i)),i=1,nphyvario3d) endif ndynvario2d=0 @@ -1032,7 +1062,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) else if (ifindstrloc(vartracers,trim(vartem)) > 0) then ntracerio2d=ntracerio2d+1 else if(trim(vartem)=='z') then - write(6,*)'the metvarname ',trim(vartem),' will be dealt separately' + if(mype == 0) write(6,*)'the metvarname ',trim(vartem),' will be dealt separately' else if(trim(vartem)=='t2m') then else if(trim(vartem)=='q2m') then else @@ -1120,6 +1150,11 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) ntracerio2d=0 endif + if( if_model_dbz )then + call gsi_bundlecreate(gsibundle_fv3lam_phyvar_nouv,GSI_MetGuess_Bundle(it)%grid,'gsibundle_fv3lam_phyvar_nouv',istatus, & + names3d=fv3lam_io_phymetvars3d_nouv) + end if + if (laeroana_fv3cmaq) then if (allocated(fv3lam_io_tracerchemvars3d_nouv) ) then call gsi_bundlecreate(gsibundle_fv3lam_tracerchem_nouv,GSI_ChemGuess_Bundle(it)%grid,'gsibundle_fv3lam_tracerchem_nouv',istatus, & @@ -1183,6 +1218,23 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) endif + if( if_model_dbz )then + inner_vars=1 + numfields=inner_vars*(nphyvario3d*grd_a%nsig) + deallocate(lnames,names) + allocate(lnames(1,numfields),names(1,numfields)) + ilev=1 + do i=1,nphyvario3d + do k=1,grd_a%nsig + lnames(1,ilev)=k + names(1,ilev)=trim(fv3lam_io_phymetvars3d_nouv(i)) + ilev=ilev+1 + enddo + enddo + call general_sub2grid_create_info(grd_fv3lam_phyvar_ionouv,inner_vars,grd_a%nlat,& + grd_a%nlon,grd_a%nsig,numfields,regional,names=names,lnames=lnames) + end if + inner_vars=2 numfields=grd_a%nsig allocate(uvlnames(inner_vars,numfields),uvnames(inner_vars,numfields)) @@ -1207,15 +1259,19 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'tv' ,ges_tv ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'q' ,ges_q ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'oz' ,ges_oz ,istatus );ier=ier+istatus - if (l_use_dbz_directDA) then + if (l_use_dbz_directDA .or. if_model_dbz) then call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ql' ,ges_ql ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qi' ,ges_qi ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qr' ,ges_qr ,istatus );ier=ier+istatus - call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'iqr' ,ges_iqr ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qs' ,ges_qs ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qg' ,ges_qg ,istatus );ier=ier+istatus - call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnr',ges_qnr ,istatus );ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'w' , ges_w ,istatus );ier=ier+istatus + if (l_use_dbz_directDA) then + call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'iqr' ,ges_iqr ,istatus );ier=ier+istatus + call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnr',ges_qnr ,istatus );ier=ier+istatus + end if + if(if_model_dbz) & + call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'dbz' , ges_dbz ,istatus );ier=ier+istatus end if if (ier/=0) call die(trim(myname),'cannot get pointers for fv3 met-fields, ier =',ier) @@ -1265,6 +1321,10 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) & ,fv3filenamegin(it)%dynvars,fv3filenamegin(it)) call gsi_fv3ncdf_read(grd_fv3lam_tracer_ionouv,gsibundle_fv3lam_tracer_nouv & & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) + if( if_model_dbz )then + call gsi_fv3ncdf_read(grd_fv3lam_phyvar_ionouv,gsibundle_fv3lam_phyvar_nouv & + & ,fv3filenamegin(it)%phyvars,fv3filenamegin(it)) + end if if (laeroana_fv3cmaq) then call gsi_fv3ncdf_read(grd_fv3lam_tracerchem_ionouv,gsibundle_fv3lam_tracerchem_nouv & & ,fv3filenamegin(it)%tracers,fv3filenamegin(it)) @@ -1330,6 +1390,7 @@ subroutine read_fv3_netcdf_guess(fv3filenamegin) endif call gsi_copy_bundle(gsibundle_fv3lam_dynvar_nouv,GSI_MetGuess_Bundle(it)) call gsi_copy_bundle(gsibundle_fv3lam_tracer_nouv,GSI_MetGuess_Bundle(it)) + if(if_model_dbz) call gsi_copy_bundle(gsibundle_fv3lam_phyvar_nouv,GSI_MetGuess_Bundle(it)) call GSI_BundleGetPointer ( gsibundle_fv3lam_dynvar_nouv, 'tsen' ,ges_tsen_readin ,istatus );ier=ier+istatus !! tsen2tv !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do k=1,nsig @@ -2042,7 +2103,8 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) character(len=max_varname_length) :: varname,vgsiname character(len=max_varname_length) :: name character(len=max_varname_length) :: filenamein2 - + real(r_kind),allocatable,dimension(:,:):: uu2d_tmp + integer(i_kind) :: countloc_tmp(3),startloc_tmp(3) integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3) integer(i_kind) ilev,ilevtot,inative @@ -2103,6 +2165,13 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) inative=nzp1-ilev startloc=(/1,1,inative/) countloc=(/nxcase,nycase,1/) + ! Variable ref_f3d in phy_data.nc has a smaller domain size than + ! dynvariables and tracers as well as a reversed order in vertical + if ( trim(adjustl(varname)) == 'ref_f3d' )then + allocate(uu2d_tmp(nxcase,nycase)) + startloc_tmp=(/1,1,ilev/) + countloc_tmp=(/nxcase,nycase,1/) + end if if(fv3_io_layout_y > 1) then do nio=0,fv3_io_layout_y-1 @@ -2115,7 +2184,18 @@ subroutine gsi_fv3ncdf_read(grd_ionouv,cstate_nouv,filenamein,fv3filenamegin) enddo else iret=nf90_inq_varid(gfile_loc,trim(adjustl(varname)),var_id) - iret=nf90_get_var(gfile_loc,var_id,uu2d,start=startloc,count=countloc) + if ( trim(adjustl(varname)) == 'ref_f3d' )then + uu2d = 0.0_r_kind + iret=nf90_get_var(gfile_loc,var_id,uu2d_tmp,start=startloc_tmp,count=countloc_tmp) + where(uu2d_tmp < 0.0_r_kind) + uu2d_tmp = 0.0_r_kind + endwhere + + uu2d(1:nxcase,1:nycase) = uu2d_tmp + deallocate(uu2d_tmp) + else + iret=nf90_get_var(gfile_loc,var_id,uu2d,start=startloc,count=countloc) + end if endif call fv3_h_to_ll(uu2d,hwork(1,:,:,ilevtot),nxcase,nycase,nloncase,nlatcase,grid_reverse_flag) @@ -2571,6 +2651,7 @@ subroutine wrfv3_netcdf(fv3filenamegin) use directDA_radaruse_mod, only: l_cvpnr, cvpnr_pval use gridmod, only: eta1_ll,eta2_ll use constants, only: one + use obsmod, only: if_model_dbz implicit none @@ -2596,6 +2677,7 @@ subroutine wrfv3_netcdf(fv3filenamegin) real(r_kind),pointer,dimension(:,:,:):: ges_qg =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_qnr =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_w =>NULL() + real(r_kind),pointer,dimension(:,:,:):: ges_dbz =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_delzinc =>NULL() real(r_kind),pointer,dimension(:,:,:):: ges_delp =>NULL() real(r_kind),dimension(:,: ),allocatable:: ges_ps_write @@ -2685,14 +2767,17 @@ subroutine wrfv3_netcdf(fv3filenamegin) call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'u' , ges_u ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'v' , ges_v ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'q' ,ges_q ,istatus);ier=ier+istatus - if (l_use_dbz_directDA) then + if (l_use_dbz_directDA .or. if_model_dbz) then call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'ql' ,ges_ql ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qi' ,ges_qi ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qr' ,ges_qr ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qs' ,ges_qs ,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qg' ,ges_qg ,istatus);ier=ier+istatus + if (l_use_dbz_directDA) & call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'qnr',ges_qnr,istatus);ier=ier+istatus call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'w' , ges_w ,istatus);ier=ier+istatus + if( if_model_dbz )& + call GSI_BundleGetPointer ( GSI_MetGuess_Bundle(it), 'dbz' , ges_dbz ,istatus);ier=ier+istatus end if if(i_use_2mq4b > 0 .and. i_use_2mt4b > 0 ) then call GSI_BundleGetPointer (GSI_MetGuess_Bundle(it),'q2m',ges_q2m,istatus); ier=ier+istatus @@ -2822,6 +2907,7 @@ subroutine wrfv3_netcdf(fv3filenamegin) call gsi_copy_bundle(GSI_MetGuess_Bundle(it),gsibundle_fv3lam_dynvar_nouv) call gsi_copy_bundle(GSI_MetGuess_Bundle(it),gsibundle_fv3lam_tracer_nouv) + if( if_model_dbz ) call gsi_copy_bundle(GSI_MetGuess_Bundle(it),gsibundle_fv3lam_phyvar_nouv) if (laeroana_fv3cmaq) then call gsi_copy_bundle(GSI_ChemGuess_Bundle(it),gsibundle_fv3lam_tracerchem_nouv) end if @@ -2866,6 +2952,10 @@ subroutine wrfv3_netcdf(fv3filenamegin) add_saved,fv3filenamegin%dynvars,fv3filenamegin) call gsi_fv3ncdf_write(grd_fv3lam_tracer_ionouv,gsibundle_fv3lam_tracer_nouv, & add_saved,fv3filenamegin%tracers,fv3filenamegin) + if( if_model_dbz ) then + call gsi_fv3ncdf_write(grd_fv3lam_phyvar_ionouv,gsibundle_fv3lam_phyvar_nouv,& + add_saved,fv3filenamegin%phyvars,fv3filenamegin) + end if call gsi_fv3ncdf_writeuv(grd_fv3lam_uv,ges_u,ges_v,add_saved,fv3filenamegin) if (laeroana_fv3cmaq) then call gsi_fv3ncdf_write(grd_fv3lam_tracerchem_ionouv,gsibundle_fv3lam_tracerchem_nouv, & @@ -3490,6 +3580,7 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file character(len=max_varname_length) :: varname,vgsiname integer(i_kind) nlatcase,nloncase,nxcase,nycase,countloc(3),startloc(3) + integer(i_kind) countloc_tmp(3),startloc_tmp(3) integer(i_kind) kbgn,kend integer(i_kind) inative,ilev,ilevtot integer(i_kind) :: VarId,gfile_loc @@ -3497,6 +3588,7 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file real(r_kind),allocatable,dimension(:,:):: work_a real(r_kind),allocatable,dimension(:,:):: work_b real(r_kind),allocatable,dimension(:,:):: workb2,worka2 + real(r_kind),allocatable,dimension(:,:):: work_b_tmp ! for io_layout > 1 real(r_kind),allocatable,dimension(:,:):: work_b_layout @@ -3551,7 +3643,11 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file work_a=hwork(1,:,:,ilevtot) - + if( trim(varname) == 'ref_f3d' )then + allocate(work_b_tmp(nlon_regional,nlat_regional)) + countloc_tmp=(/nxcase,nycase,1/) + startloc_tmp=(/1,1,ilev/) + end if call check( nf90_inq_varid(gfile_loc,trim(varname),VarId) ) @@ -3581,7 +3677,16 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file deallocate(work_b_layout) enddo else - call check( nf90_get_var(gfile_loc,VarId,work_b,start = startloc, count = countloc) ) + if( trim(varname) == 'ref_f3d' )then + work_b = 0.0_r_kind + call check( nf90_get_var(gfile_loc,VarId,work_b_tmp,start = startloc_tmp, count = countloc_tmp) ) + where(work_b_tmp < 0.0_r_kind) + work_b_tmp = 0.0_r_kind + end where + work_b(1:nxcase,1:nycase) = work_b_tmp + else + call check( nf90_get_var(gfile_loc,VarId,work_b,start = startloc, count = countloc) ) + end if endif call fv3_h_to_ll(work_b(:,:),worka2,nlon_regional,nlat_regional,nloncase,nlatcase,grid_reverse_flag) !!!!!!!! analysis_inc: work_a !!!!!!!!!!!!!!!! @@ -3601,7 +3706,16 @@ subroutine gsi_fv3ncdf_write(grd_ionouv,cstate_nouv,add_saved,filenamein,fv3file deallocate(work_b_layout) enddo else - call check( nf90_put_var(gfile_loc,VarId,work_b, start = startloc, count = countloc) ) + if( trim(varname) == 'ref_f3d' )then + work_b_tmp = work_b(1:nxcase,1:nycase) + where(work_b_tmp < 0.0_r_kind) + work_b_tmp = 0.0_r_kind + end where + call check( nf90_put_var(gfile_loc,VarId,work_b_tmp, start = startloc_tmp, count = countloc_tmp) ) + deallocate(work_b_tmp) + else + call check( nf90_put_var(gfile_loc,VarId,work_b, start = startloc, count = countloc) ) + end if endif enddo !ilevtotl loop @@ -4321,6 +4435,8 @@ subroutine getfv3lamfilevname(vgsinamein,fv3filenamegref,filenameout,vname) filenameout=fv3filenamegref%dynvars else if(ifindstrloc(vartracers,vgsinamein)> 0 ) then filenameout=fv3filenamegref%tracers + else if(ifindstrloc(varphyvars,vgsinamein)> 0) then + filenameout=fv3filenamegref%phyvars else write(6,*)'the filename corresponding to var ',trim(vgsinamein),' is not found, stop ' call stop2(333) diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index ff636bd70..aa61ed45f 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,& + nsclgrp,smooth_scales,cross_correlation_reset,para_covwithsclgrp,l_sum_spc_weights,& + vdl_scale,small_loc_varlist,smooth_scales_num 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,& @@ -1356,7 +1358,9 @@ module gsimod 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,smooth_scales,cross_correlation_reset,para_covwithsclgrp,l_sum_spc_weights,& + vdl_scale,small_loc_varlist,smooth_scales_num ! rapidrefresh_cldsurf (options for cloud analysis and surface ! enhancement for RR appilcation ): diff --git a/src/gsi/hybrid_ensemble_isotropic.F90 b/src/gsi/hybrid_ensemble_isotropic.F90 index 717815751..07fb115d2 100644 --- a/src/gsi/hybrid_ensemble_isotropic.F90 +++ b/src/gsi/hybrid_ensemble_isotropic.F90 @@ -48,6 +48,10 @@ module hybrid_ensemble_isotropic ! 2015-04-07 carley - bug fix to allow grd_loc%nlat=grd_loc%nlon ! 2016-05-13 parrish - remove beta12mult ! 2018-02-15 wu - add code for fv3_regional option +! 2022-08-29 Y. Wang, X. Wang - add capabilties of multiscale DA by using +! scale- and variable-dependent localization, +! Wang and Wang (2022ab, MWR) +! poc: xuguang.wang@ou.edu ! ! subroutines included: ! sub init_rf_z - initialize localization recursive filter (z direction) @@ -133,6 +137,8 @@ module hybrid_ensemble_isotropic public :: bkerror_a_en public :: ckgcov_a_en_new_factorization public :: ckgcov_a_en_new_factorization_ad + public :: ckgcov_a_en_new_factorization_sdl + public :: ckgcov_a_en_new_factorization_ad_sdl public :: hybens_grid_setup public :: hybens_localization_setup public :: convert_km_to_grid_units @@ -160,17 +166,17 @@ module hybrid_ensemble_isotropic real(r_kind),allocatable:: fmatz(:,:,:,:) real(r_kind),allocatable:: fmat0z(:,:,:) - real(r_kind),allocatable:: fmatx(:,:,:,:,:) - real(r_kind),allocatable:: fmat0x(:,:,:,:) - real(r_kind),allocatable:: fmaty(:,:,:,:) - real(r_kind),allocatable:: fmat0y(:,:,:) + real(r_kind),allocatable:: fmatx(:,:,:,:,:,:) + real(r_kind),allocatable:: fmat0x(:,:,:,:,:) + real(r_kind),allocatable:: fmaty(:,:,:,:,:) + real(r_kind),allocatable:: fmat0y(:,:,:,:) real(r_kind),allocatable:: znorm_new(:,:) - real(r_kind),allocatable:: xnorm_new(:,:,:) - real(r_kind),allocatable:: ynorm_new(:,:) + real(r_kind),allocatable:: xnorm_new(:,:,:,:) + real(r_kind),allocatable:: ynorm_new(:,:,:) real(r_kind),allocatable:: psbar(:) ! Other local variables for horizontal/spectral localization - real(r_kind),allocatable,dimension(:,:) :: spectral_filter,sqrt_spectral_filter + real(r_kind),allocatable,dimension(:,:,:) :: spectral_filter,sqrt_spectral_filter integer(i_kind),allocatable,dimension(:) :: k_index ! following is for special subdomain to slab variables used when internally generating ensemble members @@ -374,37 +380,39 @@ subroutine init_rf_x(x_len,kl) ! !$$$ - use hybrid_ensemble_parameters, only: grd_loc + use hybrid_ensemble_parameters, only: grd_loc,nsclgrp use hybrid_ensemble_parameters, only: region_dx_ens,region_dy_ens use constants, only: half implicit none - real(r_kind),intent(in ) :: x_len(kl) + real(r_kind),intent(in ) :: x_len(nsclgrp,kl) integer(i_kind),intent(in ) :: kl - integer(i_kind) i,j,k,l,kk + integer(i_kind) i,j,k,l,kk,ig real(r_kind) aspect(grd_loc%nlon) real(r_kind) fmatc(2,grd_loc%nlon,2),fmat0c(grd_loc%nlon,2) ! use new factorization: if(allocated(fmatx)) deallocate(fmatx) if(allocated(fmat0x)) deallocate(fmat0x) - allocate(fmatx(grd_loc%nlat,2,grd_loc%nlon,2,kl),fmat0x(grd_loc%nlat,grd_loc%nlon,2,kl)) + allocate(fmatx(nsclgrp,grd_loc%nlat,2,grd_loc%nlon,2,kl),fmat0x(nsclgrp,grd_loc%nlat,grd_loc%nlon,2,kl)) do k=1,kl do i=1,grd_loc%nlat - do j=1,grd_loc%nlon - aspect(j)=(x_len(k)*region_dy_ens(grd_loc%nlat/2,grd_loc%nlon/2)/region_dx_ens(i,j))**2 ! only works for rotated lat-lon grids - enddo - call get_new_alpha_beta(aspect,grd_loc%nlon,fmatc,fmat0c) - do kk=1,2 - do j=1,grd_loc%nlon - do l=1,2 - fmatx(i,l,j,kk,k)=fmatc(l,j,kk) - enddo - fmat0x(i,j,kk,k)=fmat0c(j,kk) - enddo - enddo + do ig=1,nsclgrp + do j=1,grd_loc%nlon + aspect(j)=(x_len(ig,k)*region_dy_ens(grd_loc%nlat/2,grd_loc%nlon/2)/region_dx_ens(i,j))**2 ! only works for rotated lat-lon grids + enddo + call get_new_alpha_beta(aspect,grd_loc%nlon,fmatc,fmat0c) + do kk=1,2 + do j=1,grd_loc%nlon + do l=1,2 + fmatx(ig,i,l,j,kk,k)=fmatc(l,j,kk) + enddo + fmat0x(ig,i,j,kk,k)=fmat0c(j,kk) + enddo + enddo + end do enddo enddo return @@ -437,26 +445,28 @@ subroutine init_rf_y(y_len,kl) ! !$$$ - use hybrid_ensemble_parameters, only: grd_loc + use hybrid_ensemble_parameters, only: grd_loc,nsclgrp use constants, only: half implicit none - real(r_kind),intent(in ) :: y_len(kl) + real(r_kind),intent(in ) :: y_len(nsclgrp,kl) integer(i_kind),intent(in ) :: kl real(r_kind) aspect(grd_loc%nlat) - integer(i_kind) i,k + integer(i_kind) i,k,ig ! use new factorization: if(allocated(fmaty)) deallocate(fmaty) if(allocated(fmat0y)) deallocate(fmat0y) - allocate(fmaty(2,grd_loc%nlat,2,kl),fmat0y(grd_loc%nlat,2,kl)) + allocate(fmaty(nsclgrp,2,grd_loc%nlat,2,kl),fmat0y(nsclgrp,grd_loc%nlat,2,kl)) do k=1,kl - do i=1,grd_loc%nlat - aspect(i)=y_len(k)**2 - enddo - call get_new_alpha_beta(aspect,grd_loc%nlat,fmaty(1,1,1,k),fmat0y(1,1,k)) + do ig=1,nsclgrp + do i=1,grd_loc%nlat + aspect(i)=y_len(ig,k)**2 + enddo + call get_new_alpha_beta(aspect,grd_loc%nlat,fmaty(ig,:,:,:,k),fmat0y(ig,:,:,k)) + end do enddo return @@ -537,7 +547,7 @@ subroutine new_factorization_rf_z(f,iadvance,iback) end subroutine new_factorization_rf_z -subroutine new_factorization_rf_x(f,iadvance,iback,nlevs) +subroutine new_factorization_rf_x(iscale,f,iadvance,iback,nlevs) !$$$ subprogram documentation block ! . . . ! subprogram: new_factorization_rf_x @@ -568,7 +578,7 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs) use hybrid_ensemble_parameters, only: grd_loc,vvlocal implicit none - integer(i_kind),intent(in ) :: iadvance,iback,nlevs + integer(i_kind),intent(in ) :: iadvance,iback,nlevs,iscale real(r_kind) ,intent(inout) :: f(grd_loc%nlat,grd_loc%nlon,max(nlevs,1)) integer(i_kind) i,j,k,l,ny,nx,nz @@ -582,7 +592,7 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs) if(iadvance == 1) then do j=1,nx do i=1,ny - f(i,j,k)=xnorm_new(i,j,k)*f(i,j,k) + f(i,j,k)=xnorm_new(iscale,i,j,k)*f(i,j,k) enddo enddo end if @@ -590,29 +600,29 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs) do j=1,nx do l=1,min(2,j-1) do i=1,ny - f(i,j,k)=f(i,j,k)-fmatx(i,l,j,iadvance,k)*f(i,j-l,k) + f(i,j,k)=f(i,j,k)-fmatx(iscale,i,l,j,iadvance,k)*f(i,j-l,k) enddo enddo do i=1,ny - f(i,j,k)=fmat0x(i,j,iadvance,k)*f(i,j,k) + f(i,j,k)=fmat0x(iscale,i,j,iadvance,k)*f(i,j,k) enddo enddo do j=nx,1,-1 do l=1,min(2,nx-j) do i=1,ny - f(i,j,k)=f(i,j,k)-fmatx(i,l,j+l,iback,k)*f(i,j+l,k) + f(i,j,k)=f(i,j,k)-fmatx(iscale,i,l,j+l,iback,k)*f(i,j+l,k) enddo enddo do i=1,ny - f(i,j,k)=fmat0x(i,j,iback,k)*f(i,j,k) + f(i,j,k)=fmat0x(iscale,i,j,iback,k)*f(i,j,k) enddo enddo if(iadvance == 2) then do j=1,nx do i=1,ny - f(i,j,k)=xnorm_new(i,j,k)*f(i,j,k) + f(i,j,k)=xnorm_new(iscale,i,j,k)*f(i,j,k) enddo enddo end if @@ -625,7 +635,7 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs) if(iadvance == 1) then do j=1,nx do i=1,ny - f(i,j,k)=xnorm_new(i,j,1)*f(i,j,k) + f(i,j,k)=xnorm_new(iscale,i,j,1)*f(i,j,k) enddo enddo end if @@ -633,29 +643,29 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs) do j=1,nx do l=1,min(2,j-1) do i=1,ny - f(i,j,k)=f(i,j,k)-fmatx(i,l,j,iadvance,1)*f(i,j-l,k) + f(i,j,k)=f(i,j,k)-fmatx(iscale,i,l,j,iadvance,1)*f(i,j-l,k) enddo enddo do i=1,ny - f(i,j,k)=fmat0x(i,j,iadvance,1)*f(i,j,k) + f(i,j,k)=fmat0x(iscale,i,j,iadvance,1)*f(i,j,k) enddo enddo do j=nx,1,-1 do l=1,min(2,nx-j) do i=1,ny - f(i,j,k)=f(i,j,k)-fmatx(i,l,j+l,iback,1)*f(i,j+l,k) + f(i,j,k)=f(i,j,k)-fmatx(iscale,i,l,j+l,iback,1)*f(i,j+l,k) enddo enddo do i=1,ny - f(i,j,k)=fmat0x(i,j,iback,1)*f(i,j,k) + f(i,j,k)=fmat0x(iscale,i,j,iback,1)*f(i,j,k) enddo enddo if(iadvance == 2) then do j=1,nx do i=1,ny - f(i,j,k)=xnorm_new(i,j,1)*f(i,j,k) + f(i,j,k)=xnorm_new(iscale,i,j,1)*f(i,j,k) enddo enddo end if @@ -665,7 +675,7 @@ subroutine new_factorization_rf_x(f,iadvance,iback,nlevs) return end subroutine new_factorization_rf_x -subroutine new_factorization_rf_y(f,iadvance,iback,nlevs) +subroutine new_factorization_rf_y(iscale,f,iadvance,iback,nlevs) !$$$ subprogram documentation block ! . . . ! subprogram: new_factorization_rf_y @@ -697,7 +707,7 @@ subroutine new_factorization_rf_y(f,iadvance,iback,nlevs) ! use mpimod, only: mype implicit none - integer(i_kind),intent(in ) :: iadvance,iback,nlevs + integer(i_kind),intent(in ) :: iadvance,iback,nlevs,iscale real(r_kind) ,intent(inout) :: f(grd_loc%nlat,grd_loc%nlon,max(nlevs,1)) integer(i_kind) i,j,k,l,nx,ny,nz @@ -710,27 +720,27 @@ subroutine new_factorization_rf_y(f,iadvance,iback,nlevs) if(iadvance == 1) then do i=1,ny - f(i,j,k)=ynorm_new(i,k)*f(i,j,k) + f(i,j,k)=ynorm_new(iscale,i,k)*f(i,j,k) enddo end if do i=1,ny do l=1,min(2,i-1) - f(i,j,k)=f(i,j,k)-fmaty(l,i,iadvance,k)*f(i-l,j,k) + f(i,j,k)=f(i,j,k)-fmaty(iscale,l,i,iadvance,k)*f(i-l,j,k) enddo - f(i,j,k)=fmat0y(i,iadvance,k)*f(i,j,k) + f(i,j,k)=fmat0y(iscale,i,iadvance,k)*f(i,j,k) enddo do i=ny,1,-1 do l=1,min(2,ny-i) - f(i,j,k)=f(i,j,k)-fmaty(l,i+l,iback,k)*f(i+l,j,k) + f(i,j,k)=f(i,j,k)-fmaty(iscale,l,i+l,iback,k)*f(i+l,j,k) enddo - f(i,j,k)=fmat0y(i,iback,k)*f(i,j,k) + f(i,j,k)=fmat0y(iscale,i,iback,k)*f(i,j,k) enddo if(iadvance == 2) then do i=1,ny - f(i,j,k)=ynorm_new(i,k)*f(i,j,k) + f(i,j,k)=ynorm_new(iscale,i,k)*f(i,j,k) enddo end if @@ -742,27 +752,27 @@ subroutine new_factorization_rf_y(f,iadvance,iback,nlevs) if(iadvance == 1) then do i=1,ny - f(i,j,k)=ynorm_new(i,1)*f(i,j,k) + f(i,j,k)=ynorm_new(iscale,i,1)*f(i,j,k) enddo end if do i=1,ny do l=1,min(2,i-1) - f(i,j,k)=f(i,j,k)-fmaty(l,i,iadvance,1)*f(i-l,j,k) + f(i,j,k)=f(i,j,k)-fmaty(iscale,l,i,iadvance,1)*f(i-l,j,k) enddo - f(i,j,k)=fmat0y(i,iadvance,1)*f(i,j,k) + f(i,j,k)=fmat0y(iscale,i,iadvance,1)*f(i,j,k) enddo do i=ny,1,-1 do l=1,min(2,ny-i) - f(i,j,k)=f(i,j,k)-fmaty(l,i+l,iback,1)*f(i+l,j,k) + f(i,j,k)=f(i,j,k)-fmaty(iscale,l,i+l,iback,1)*f(i+l,j,k) enddo - f(i,j,k)=fmat0y(i,iback,1)*f(i,j,k) + f(i,j,k)=fmat0y(iscale,i,iback,1)*f(i,j,k) enddo if(iadvance == 2) then do i=1,ny - f(i,j,k)=ynorm_new(i,1)*f(i,j,k) + f(i,j,k)=ynorm_new(iscale,i,1)*f(i,j,k) enddo end if @@ -875,12 +885,12 @@ subroutine normal_new_factorization_rf_x !$$$ end documentation block use kinds, only: r_kind,i_kind - use hybrid_ensemble_parameters, only: grd_loc,vvlocal + use hybrid_ensemble_parameters, only: grd_loc,vvlocal,nsclgrp use constants, only: zero,one implicit none - integer(i_kind) i,j,k,iadvance,iback,kl + integer(i_kind) i,j,k,iadvance,iback,kl,ig real(r_kind) f(grd_loc%nlat,grd_loc%nlon,grd_loc%kend_alloc+1-grd_loc%kbegin_loc) real(r_kind),allocatable:: diag(:,:,:) ! real(r_kind) diag(grd_loc%nlat,grd_loc%nlon) @@ -896,55 +906,59 @@ subroutine normal_new_factorization_rf_x kl=1 endif if(allocated(xnorm_new)) deallocate(xnorm_new) - allocate(xnorm_new(grd_loc%nlat,grd_loc%nlon,kl)) + allocate(xnorm_new(nsclgrp,grd_loc%nlat,grd_loc%nlon,kl)) if(allocated(diag)) deallocate(diag) allocate(diag(grd_loc%nlat,grd_loc%nlon,kl)) xnorm_new=one - do j=1,grd_loc%nlon - f=zero - do k=1,kl - do i=1,grd_loc%nlat - f(i,j,k)=one - enddo - enddo - iadvance=1 ; iback=2 - call new_factorization_rf_x(f,iadvance,iback,kl) - iadvance=2 ; iback=1 - call new_factorization_rf_x(f,iadvance,iback,kl) - do k=1,kl - do i=1,grd_loc%nlat - diag(i,j,k)=sqrt(one/f(i,j,k)) - enddo - enddo - enddo - do k=1,kl - do j=1,grd_loc%nlon - do i=1,grd_loc%nlat - xnorm_new(i,j,k)=diag(i,j,k) - enddo - enddo - enddo + do ig=1,nsclgrp + do j=1,grd_loc%nlon + f=zero + do k=1,kl + do i=1,grd_loc%nlat + f(i,j,k)=one + enddo + enddo + iadvance=1 ; iback=2 + call new_factorization_rf_x(ig,f,iadvance,iback,kl) + iadvance=2 ; iback=1 + call new_factorization_rf_x(ig,f,iadvance,iback,kl) + do k=1,kl + do i=1,grd_loc%nlat + diag(i,j,k)=sqrt(one/f(i,j,k)) + enddo + enddo + enddo + do k=1,kl + do j=1,grd_loc%nlon + do i=1,grd_loc%nlat + xnorm_new(ig,i,j,k)=diag(i,j,k) + enddo + enddo + enddo + end do ! check accuracy of xnorm if(debug) then - do j=1,grd_loc%nlon - f=zero - do k=1,kl - do i=1,grd_loc%nlat - f(i,j,k)=one - enddo - enddo - iadvance=1 ; iback=2 - call new_factorization_rf_x(f,iadvance,iback,kl) - iadvance=2 ; iback=1 - call new_factorization_rf_x(f,iadvance,iback,kl) - do k=1,kl - do i=1,grd_loc%nlat - diag(i,j,k)=f(i,j,k) - enddo - enddo - enddo - write(6,*)' in normal_new_factorization_rf_x,min,max(diag)=',minval(diag),maxval(diag) + do ig=1,nsclgrp + do j=1,grd_loc%nlon + f=zero + do k=1,kl + do i=1,grd_loc%nlat + f(i,j,k)=one + enddo + enddo + iadvance=1 ; iback=2 + call new_factorization_rf_x(ig,f,iadvance,iback,kl) + iadvance=2 ; iback=1 + call new_factorization_rf_x(ig,f,iadvance,iback,kl) + do k=1,kl + do i=1,grd_loc%nlat + diag(i,j,k)=f(i,j,k) + enddo + enddo + enddo + write(6,*)' in normal_new_factorization_rf_x,min,max(diag)=',minval(diag),maxval(diag) + end do endif return @@ -976,11 +990,11 @@ subroutine normal_new_factorization_rf_y !$$$ end documentation block use kinds, only: r_kind,i_kind - use hybrid_ensemble_parameters, only: grd_loc,vvlocal + use hybrid_ensemble_parameters, only: grd_loc,vvlocal,nsclgrp use constants, only: zero,one implicit none - integer(i_kind) i,k,lend,lcount,iadvance,iback,kl,loop,ll,iend + integer(i_kind) i,k,lend,lcount,iadvance,iback,kl,loop,ll,iend,ig ! real(r_kind) f(grd_loc%nlat,grd_loc%nlon*(grd_loc%kend_alloc+1-grd_loc%kbegin_loc)),diag(grd_loc%nlat) real(r_kind) f(grd_loc%nlat,grd_loc%nlon,grd_loc%kend_alloc+1-grd_loc%kbegin_loc) real(r_kind),allocatable:: diag(:,:) @@ -997,7 +1011,7 @@ subroutine normal_new_factorization_rf_y endif if(allocated(ynorm_new)) deallocate(ynorm_new) - allocate(ynorm_new(grd_loc%nlat,kl)) + allocate(ynorm_new(nsclgrp,grd_loc%nlat,kl)) if(allocated(diag)) deallocate(diag) allocate(diag(grd_loc%nlat,kl)) @@ -1012,58 +1026,62 @@ subroutine normal_new_factorization_rf_y iend=grd_loc%nlon endif - do loop=1,lend - ll=(loop-1)*iend - f=zero - do k=1,kl - do i=1,iend - lcount=ll+i - f(lcount,i,k)=one - if(lcount == grd_loc%nlat) exit - enddo - enddo - - iadvance=1 ; iback=2 - call new_factorization_rf_y(f,iadvance,iback,kl) - iadvance=2 ; iback=1 - call new_factorization_rf_y(f,iadvance,iback,kl) - - do k=1,kl - do i=1,iend - lcount=ll+i - diag(lcount,k)=sqrt(one/f(lcount,i,k)) - ynorm_new(lcount,k)=diag(lcount,k) - if(lcount == grd_loc%nlat) exit - enddo - enddo - enddo + do ig=1,nsclgrp + do loop=1,lend + ll=(loop-1)*iend + f=zero + do k=1,kl + do i=1,iend + lcount=ll+i + f(lcount,i,k)=one + if(lcount == grd_loc%nlat) exit + enddo + enddo + + iadvance=1 ; iback=2 + call new_factorization_rf_y(ig,f,iadvance,iback,kl) + iadvance=2 ; iback=1 + call new_factorization_rf_y(ig,f,iadvance,iback,kl) + + do k=1,kl + do i=1,iend + lcount=ll+i + diag(lcount,k)=sqrt(one/f(lcount,i,k)) + ynorm_new(ig,lcount,k)=diag(lcount,k) + if(lcount == grd_loc%nlat) exit + enddo + enddo + enddo + end do ! check that ynorm is corect if(debug) then - do loop=1,lend - ll=(loop-1)*iend - f=zero - do k=1,kl - do i=1,iend - lcount=ll+i - f(lcount,i,k)=one - if(lcount == grd_loc%nlat) exit - enddo - enddo + do ig=1,nsclgrp + do loop=1,lend + ll=(loop-1)*iend + f=zero + do k=1,kl + do i=1,iend + lcount=ll+i + f(lcount,i,k)=one + if(lcount == grd_loc%nlat) exit + enddo + enddo - iadvance=1 ; iback=2 - call new_factorization_rf_y(f,iadvance,iback,kl) - iadvance=2 ; iback=1 - call new_factorization_rf_y(f,iadvance,iback,kl) + iadvance=1 ; iback=2 + call new_factorization_rf_y(ig,f,iadvance,iback,kl) + iadvance=2 ; iback=1 + call new_factorization_rf_y(ig,f,iadvance,iback,kl) - do k=1,kl - do i=1,iend - lcount=ll+i - diag(lcount,k)=sqrt(one/f(lcount,i,k)) - if(lcount == grd_loc%nlat) exit - enddo - enddo - enddo - write(6,*)' in normal_new_factorization_rf_y, min,max(diag)=',minval(diag),maxval(diag) + do k=1,kl + do i=1,iend + lcount=ll+i + diag(lcount,k)=sqrt(one/f(lcount,i,k)) + if(lcount == grd_loc%nlat) exit + enddo + enddo + enddo + write(6,*)' in normal_new_factorization_rf_y, min,max(diag)=',minval(diag),maxval(diag) + end do endif return end subroutine normal_new_factorization_rf_y @@ -1093,30 +1111,32 @@ subroutine create_ensemble ! !$$$ use hybrid_ensemble_parameters, only: n_ens,grd_ens,ntlevs_ens - use hybrid_ensemble_parameters, only: nelen,en_perts,ps_bar + use hybrid_ensemble_parameters, only: nelen,en_perts,ps_bar,nsclgrp implicit none type(gsi_grid) :: grid_ens - integer(i_kind) n,istatus,m + integer(i_kind) n,istatus,m,ig character(len=*),parameter::myname_=trim(myname)//'*create_ensemble' nelen=grd_ens%latlon11*(max(0,nc3d)*grd_ens%nsig+max(0,nc2d)) ! create ensemble perturbations bundles (using newly added r_single capability - allocate(en_perts(n_ens,ntlevs_ens)) + allocate(en_perts(n_ens,nsclgrp,ntlevs_ens)) call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig) do m=1,ntlevs_ens - do n=1,n_ens - call gsi_bundlecreate(en_perts(n,m),grid_ens,'ensemble perts',istatus, & - names2d=cvars2d,names3d=cvars3d,bundle_kind=r_single) - if(istatus/=0) then - write(6,*)trim(myname_),': trouble creating en_perts bundle' - call stop2(999) - endif - enddo + do ig =1,nsclgrp + do n=1,n_ens + call gsi_bundlecreate(en_perts(n,ig,m),grid_ens,'ensemble perts',istatus, & + names2d=cvars2d,names3d=cvars3d,bundle_kind=r_single) + if(istatus/=0) then + write(6,*)trim(myname_),': trouble creating en_perts bundle' + call stop2(999) + endif + enddo + end do enddo @@ -1176,7 +1196,7 @@ subroutine load_ensemble 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,nsclgrp use mpimod, only: mpi_comm_world implicit none @@ -1188,7 +1208,7 @@ subroutine load_ensemble type(gsi_bundle),allocatable:: en_bar(:) type(gsi_bundle):: bundle_anl,bundle_ens type(gsi_grid) :: grid_anl,grid_ens - integer(i_kind) i,j,n,ii,m + integer(i_kind) i,j,n,ii,m,ig integer(i_kind) istatus real(r_kind),allocatable:: seed(:,:) real(r_kind),pointer,dimension(:,:) :: cv_ps=>NULL() @@ -1245,7 +1265,7 @@ subroutine load_ensemble do n=1,n_ens call generate_one_ensemble_perturbation(bundle_anl,bundle_ens,seed) do ii=1,nelen - en_perts(n,m)%valuesr4(ii)=bundle_ens%values(ii) + en_perts(n,1,m)%valuesr4(ii)=bundle_ens%values(ii) en_bar(m)%values(ii)=en_bar(m)%values(ii)+bundle_ens%values(ii) enddo enddo @@ -1276,9 +1296,9 @@ subroutine load_ensemble do m=1,ntlevs_ens do n=1,n_ens do ii=1,nelen - en_perts(n,m)%valuesr4(ii)=(en_perts(n,m)%valuesr4(ii)-en_bar(m)%values(ii)*bar_norm)*sig_norm + en_perts(n,1,m)%valuesr4(ii)=(en_perts(n,1,m)%valuesr4(ii)-en_bar(m)%values(ii)*bar_norm)*sig_norm enddo - call gsi_enscoupler_put_gsi_ens(grd_ens,n,m,en_perts(n,m),istatus) + call gsi_enscoupler_put_gsi_ens(grd_ens,n,m,en_perts(n,1,m),istatus) if(istatus/=0) then write(6,*)trim(myname_),': trouble writing perts' call stop2(999) @@ -1633,16 +1653,16 @@ subroutine rescale_ensemble_rh_perturbations use kinds, only: r_kind,i_kind use gridmod, only: regional use hybrid_ensemble_parameters, only: n_ens,grd_ens,grd_anl,grd_a1,grd_e1,p_e2a,ntlevs_ens - use hybrid_ensemble_parameters, only: en_perts + use hybrid_ensemble_parameters, only: en_perts,nsclgrp use general_sub2grid_mod, only: general_suba2sube use berror, only: qvar3d implicit none - integer(i_kind) i,j,k,n,istatus,m + integer(i_kind) i,j,k,n,istatus,m,ig real(r_kind) qvar3d_ens(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig,1) real(r_single),pointer,dimension(:,:,:):: w3=>NULL() - call gsi_bundlegetpointer(en_perts(1,1),'q',w3,istatus) + call gsi_bundlegetpointer(en_perts(1,1,1),'q',w3,istatus) if(istatus/=0) then write(6,*)' rh variable not available, skip subroutine rescale_ensemble_rh_perturbations' return @@ -1656,20 +1676,22 @@ subroutine rescale_ensemble_rh_perturbations end if do m=1,ntlevs_ens !$omp parallel do schedule(dynamic,1) private(n,i,j,k,w3,istatus) - do n=1,n_ens - call gsi_bundlegetpointer(en_perts(n,m),'q',w3,istatus) - if(istatus/=0) then - write(6,*)' error retrieving pointer to rh variable for ensemble number ',n - call stop2(999) - end if - do k=1,grd_ens%nsig - do j=1,grd_ens%lon2 - do i=1,grd_ens%lat2 - w3(i,j,k)=qvar3d_ens(i,j,k,1)*w3(i,j,k) - enddo - enddo - enddo - enddo + do ig=1,nsclgrp + do n=1,n_ens + call gsi_bundlegetpointer(en_perts(n,ig,m),'q',w3,istatus) + if(istatus/=0) then + write(6,*)' error retrieving pointer to rh variable for ensemble number ',n + call stop2(999) + end if + do k=1,grd_ens%nsig + do j=1,grd_ens%lon2 + do i=1,grd_ens%lat2 + w3(i,j,k)=qvar3d_ens(i,j,k,1)*w3(i,j,k) + enddo + enddo + enddo + enddo + end do enddo return @@ -1698,20 +1720,22 @@ subroutine destroy_ensemble ! !$$$ use hybrid_ensemble_parameters, only: l_hyb_ens,n_ens,ntlevs_ens - use hybrid_ensemble_parameters, only: en_perts,ps_bar + use hybrid_ensemble_parameters, only: en_perts,ps_bar,nsclgrp implicit none - integer(i_kind) istatus,n,m + integer(i_kind) istatus,n,m,ig if(l_hyb_ens) then do m=1,ntlevs_ens - do n=1,n_ens - call gsi_bundleunset(en_perts(n,m),istatus) - if(istatus/=0) then - write(6,*)'in destroy_ensemble: trouble destroying en_perts bundle' - call stop2(999) - endif - enddo + do ig=1,nsclgrp + do n=1,n_ens + call gsi_bundleunset(en_perts(n,ig,m),istatus) + if(istatus/=0) then + write(6,*)'in destroy_ensemble: trouble destroying en_perts bundle' + call stop2(999) + endif + enddo + end do enddo deallocate(ps_bar) deallocate(en_perts) @@ -1758,17 +1782,17 @@ subroutine ensemble_forward_model(cvec,a_en,ibin) ! !$$$ use hybrid_ensemble_parameters, only: n_ens,pwgtflg,pwgt - use hybrid_ensemble_parameters, only: en_perts + use hybrid_ensemble_parameters, only: en_perts,nsclgrp use constants, only: zero implicit none type(gsi_bundle),intent(inout) :: cvec - type(gsi_bundle),intent(in) :: a_en(n_ens) + type(gsi_bundle),intent(in) :: a_en(nsclgrp,n_ens) integer,intent(in) :: ibin character(len=*),parameter :: myname_=trim(myname)//'*ensemble_forward_model' logical :: nogood - integer(i_kind) :: i,j,k,n,im,jm,km,ic2,ic3,ipic,ipx,km_tmp + integer(i_kind) :: i,j,k,n,im,jm,km,ic2,ic3,ipic,ipx,km_tmp,ig integer(i_kind) :: ipc3d(nc3d),ipc2d(nc2d),istatus im=cvec%grid%im @@ -1776,7 +1800,7 @@ subroutine ensemble_forward_model(cvec,a_en,ibin) km=cvec%grid%km ! Check resolution consistency between static and ensemble components - nogood=im/=a_en(1)%grid%im.or.jm/=a_en(1)%grid%jm.or.km/=a_en(1)%grid%km + nogood=im/=a_en(1,1)%grid%im.or.jm/=a_en(1,1)%grid%jm.or.km/=a_en(1,1)%grid%km if (nogood) then write(6,*) myname_,': static&ensemble vectors have inconsistent dims' call stop2(999) @@ -1798,7 +1822,7 @@ subroutine ensemble_forward_model(cvec,a_en,ibin) ipx=1 -!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ipic) +!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ipic,ig) do k=1,km do ic3=1,nc3d ipic=ipc3d(ic3) @@ -1808,17 +1832,19 @@ subroutine ensemble_forward_model(cvec,a_en,ibin) enddo enddo do n=1,n_ens - do j=1,jm - do i=1,im - cvec%r3(ipic)%q(i,j,k)=cvec%r3(ipic)%q(i,j,k) & - +a_en(n)%r3(ipx)%q(i,j,k)*en_perts(n,ibin)%r3(ipic)%qr4(i,j,k) - enddo - enddo + do ig=1,nsclgrp + do j=1,jm + do i=1,im + cvec%r3(ipic)%q(i,j,k)=cvec%r3(ipic)%q(i,j,k) & + +a_en(ig,n)%r3(ipx)%q(i,j,k)*en_perts(n,ig,ibin)%r3(ipic)%qr4(i,j,k) + enddo + enddo + end do enddo enddo enddo -!$omp parallel do schedule(dynamic,1) private(j,n,k,i,ic2,ipic) +!$omp parallel do schedule(dynamic,1) private(j,n,k,i,ic2,ipic,ig) do ic2=1,nc2d ipic=ipc2d(ic2) do j=1,jm @@ -1839,12 +1865,14 @@ subroutine ensemble_forward_model(cvec,a_en,ibin) do n=1,n_ens do j=1,jm - do k=1,km_tmp - do i=1,im - cvec%r2(ipic)%q(i,j)=cvec%r2(ipic)%q(i,j) & - +a_en(n)%r3(ipx)%q(i,j,k)*en_perts(n,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k) - enddo - enddo + do ig=1,nsclgrp + do k=1,km_tmp + do i=1,im + cvec%r2(ipic)%q(i,j)=cvec%r2(ipic)%q(i,j) & + +a_en(ig,n)%r3(ipx)%q(i,j,k)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k) + enddo + enddo + end do enddo enddo ! enddo n_ens @@ -1853,8 +1881,10 @@ subroutine ensemble_forward_model(cvec,a_en,ibin) do n=1,n_ens do j=1,jm do i=1,im - cvec%r2(ipic)%q(i,j)=cvec%r2(ipic)%q(i,j) & - +a_en(n)%r3(ipx)%q(i,j,1)*en_perts(n,ibin)%r2(ipic)%qr4(i,j) + do ig=1,nsclgrp + cvec%r2(ipic)%q(i,j)=cvec%r2(ipic)%q(i,j) & + +a_en(ig,n)%r3(ipx)%q(i,j,1)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j) + end do enddo enddo enddo ! enddo n_ens @@ -1908,20 +1938,20 @@ subroutine ensemble_forward_model_dual_res(cvec,a_en,ibin) !$$$ use hybrid_ensemble_parameters, only: n_ens,pwgtflg,pwgt use hybrid_ensemble_parameters, only: grd_ens,grd_anl,p_e2a - use hybrid_ensemble_parameters, only: en_perts + use hybrid_ensemble_parameters, only: en_perts,nsclgrp use general_sub2grid_mod, only: general_sube2suba use gridmod,only: regional use constants, only: zero implicit none type(gsi_bundle),intent(inout) :: cvec - type(gsi_bundle),intent(in) :: a_en(n_ens) + type(gsi_bundle),intent(in) :: a_en(nsclgrp,n_ens) integer,intent(in) :: ibin character(len=*),parameter::myname_=trim(myname)//'*ensemble_forward_model_dual_res' type(gsi_grid) :: grid_ens,grid_anl type(gsi_bundle) :: work_ens,work_anl - integer(i_kind) :: i,j,k,n,im,jm,km,ic2,ic3,ipic,ipx,km_tmp + integer(i_kind) :: i,j,k,n,im,jm,km,ic2,ic3,ipic,ipx,km_tmp,ig integer(i_kind) :: ipc2d(nc2d),ipc3d(nc3d),istatus ! Request ensemble-corresponding fields from control vector @@ -1958,7 +1988,7 @@ subroutine ensemble_forward_model_dual_res(cvec,a_en,ibin) im=work_ens%grid%im jm=work_ens%grid%jm km=work_ens%grid%km -!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ipic) +!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ipic,ig) do k=1,km do ic3=1,nc3d ipic=ipc3d(ic3) @@ -1968,16 +1998,18 @@ subroutine ensemble_forward_model_dual_res(cvec,a_en,ibin) enddo enddo do n=1,n_ens - do j=1,jm - do i=1,im - work_ens%r3(ipic)%q(i,j,k)=work_ens%r3(ipic)%q(i,j,k) & - +a_en(n)%r3(ipx)%q(i,j,k)*en_perts(n,ibin)%r3(ipic)%qr4(i,j,k) - enddo + do ig=1,nsclgrp + do j=1,jm + do i=1,im + work_ens%r3(ipic)%q(i,j,k)=work_ens%r3(ipic)%q(i,j,k) & + +a_en(ig,n)%r3(ipx)%q(i,j,k)*en_perts(n,ig,ibin)%r3(ipic)%qr4(i,j,k) + enddo + enddo enddo enddo enddo enddo -!$omp parallel do schedule(dynamic,1) private(j,n,k,i,ic2,ipic) +!$omp parallel do schedule(dynamic,1) private(j,n,k,i,ic2,ipic,ig) do ic2=1,nc2d ipic=ipc2d(ic2) do j=1,jm @@ -1997,24 +2029,28 @@ subroutine ensemble_forward_model_dual_res(cvec,a_en,ibin) endif do n=1,n_ens - do k=1,km_tmp - do j=1,jm - do i=1,im - work_ens%r2(ipic)%q(i,j)=work_ens%r2(ipic)%q(i,j) & - +a_en(n)%r3(ipx)%q(i,j,k)*en_perts(n,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k) - enddo - enddo + do ig=1,nsclgrp + do k=1,km_tmp + do j=1,jm + do i=1,im + work_ens%r2(ipic)%q(i,j)=work_ens%r2(ipic)%q(i,j) & + +a_en(ig,n)%r3(ipx)%q(i,j,k)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k) + enddo + enddo + enddo enddo enddo ! enddo n_ens case('SST') do n=1,n_ens - do j=1,jm - do i=1,im - work_ens%r2(ipic)%q(i,j)=work_ens%r2(ipic)%q(i,j) & - +a_en(n)%r3(ipx)%q(i,j,1)*en_perts(n,ibin)%r2(ipic)%qr4(i,j) - enddo + do ig=1,nsclgrp + do j=1,jm + do i=1,im + work_ens%r2(ipic)%q(i,j)=work_ens%r2(ipic)%q(i,j) & + +a_en(ig,n)%r3(ipx)%q(i,j,1)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j) + enddo + enddo enddo enddo ! enddo n_ens @@ -2081,23 +2117,23 @@ subroutine ensemble_forward_model_ad(cvec,a_en,ibin) !$$$ use hybrid_ensemble_parameters, only: n_ens,pwgtflg,pwgt - use hybrid_ensemble_parameters, only: en_perts + use hybrid_ensemble_parameters, only: en_perts,nsclgrp implicit none type(gsi_bundle),intent(inout) :: cvec - type(gsi_bundle),intent(inout) :: a_en(n_ens) + type(gsi_bundle),intent(inout) :: a_en(nsclgrp,n_ens) integer,intent(in) :: ibin character(len=*),parameter :: myname_=trim(myname)//'*ensemble_forward_model_ad' logical :: nogood - integer(i_kind) :: i,j,k,n,im,jm,km,ic2,ic3,ipx,ipic,km_tmp + integer(i_kind) :: i,j,k,n,im,jm,km,ic2,ic3,ipx,ipic,km_tmp,ig integer(i_kind) :: ipc3d(nc3d),ipc2d(nc2d),istatus im=cvec%grid%im jm=cvec%grid%jm km=cvec%grid%km ! Check resolution consistency between static and ensemble components - nogood=im/=a_en(1)%grid%im.or.jm/=a_en(1)%grid%jm.or.km/=a_en(1)%grid%km + nogood=im/=a_en(1,1)%grid%im.or.jm/=a_en(1,1)%grid%jm.or.km/=a_en(1,1)%grid%km if (nogood) then write(6,*) myname_,': static/ensemble vectors have inconsistent dims' call stop2(999) @@ -2118,53 +2154,55 @@ subroutine ensemble_forward_model_ad(cvec,a_en,ibin) endif ipx=1 -!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ic2,ipic) - do n=1,n_ens - do ic3=1,nc3d - ipic=ipc3d(ic3) - do k=1,km - do j=1,jm - do i=1,im - a_en(n)%r3(ipx)%q(i,j,k)=a_en(n)%r3(ipx)%q(i,j,k) & - +cvec%r3(ipic)%q(i,j,k)*en_perts(n,ibin)%r3(ipic)%qr4(i,j,k) - enddo - enddo - enddo - enddo - do ic2=1,nc2d - - ipic=ipc2d(ic2) - select case ( trim(StrUpCase(cvars2d(ic2))) ) - - case('PS') - - if ( pwgtflg ) then - km_tmp = km - else - km_tmp = 1 - endif - - do k=1,km_tmp - do j=1,jm - do i=1,im - a_en(n)%r3(ipx)%q(i,j,k)=a_en(n)%r3(ipx)%q(i,j,k) & - +cvec%r2(ipic)%q(i,j)*en_perts(n,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k) - enddo - enddo - enddo +!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ic2,ipic,ig) + do ig=1,nsclgrp + do n=1,n_ens + do ic3=1,nc3d + ipic=ipc3d(ic3) + do k=1,km + do j=1,jm + do i=1,im + a_en(ig,n)%r3(ipx)%q(i,j,k)=a_en(ig,n)%r3(ipx)%q(i,j,k) & + +cvec%r3(ipic)%q(i,j,k)*en_perts(n,ig,ibin)%r3(ipic)%qr4(i,j,k) + enddo + enddo + enddo + enddo + do ic2=1,nc2d - case('SST') + ipic=ipc2d(ic2) + select case ( trim(StrUpCase(cvars2d(ic2))) ) + + case('PS') + + if ( pwgtflg ) then + km_tmp = km + else + km_tmp = 1 + endif - do j=1,jm - do i=1,im - a_en(n)%r3(ipx)%q(i,j,1)=a_en(n)%r3(ipx)%q(i,j,1) & - +cvec%r2(ipic)%q(i,j)*en_perts(n,ibin)%r2(ipic)%qr4(i,j) - enddo - enddo - - end select - enddo - enddo ! enddo n_ens + do k=1,km_tmp + do j=1,jm + do i=1,im + a_en(ig,n)%r3(ipx)%q(i,j,k)=a_en(ig,n)%r3(ipx)%q(i,j,k) & + +cvec%r2(ipic)%q(i,j)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k) + enddo + enddo + enddo + + case('SST') + + do j=1,jm + do i=1,im + a_en(ig,n)%r3(ipx)%q(i,j,1)=a_en(ig,n)%r3(ipx)%q(i,j,1) & + +cvec%r2(ipic)%q(i,j)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j) + enddo + enddo + + end select + enddo + enddo ! enddo n_ens + end do return end subroutine ensemble_forward_model_ad @@ -2211,20 +2249,20 @@ subroutine ensemble_forward_model_ad_dual_res(cvec,a_en,ibin) !$$$ use hybrid_ensemble_parameters, only: n_ens,pwgtflg,pwgt use hybrid_ensemble_parameters, only: n_ens,grd_ens,grd_anl,p_e2a - use hybrid_ensemble_parameters, only: en_perts + use hybrid_ensemble_parameters, only: en_perts,nsclgrp use general_sub2grid_mod, only: general_sube2suba_ad use gridmod,only: regional use constants, only: zero implicit none type(gsi_bundle),intent(inout) :: cvec - type(gsi_bundle),intent(inout) :: a_en(n_ens) + type(gsi_bundle),intent(inout) :: a_en(nsclgrp,n_ens) integer,intent(in) :: ibin character(len=*),parameter::myname_=trim(myname)//'*ensemble_forward_model_ad_dual_res' type(gsi_grid) :: grid_ens,grid_anl type(gsi_bundle) :: work_ens,work_anl - integer(i_kind) :: i,j,k,n,im,jm,km,ic2,ic3,ipx,ipic,km_tmp + integer(i_kind) :: i,j,k,n,im,jm,km,ic2,ic3,ipx,ipic,km_tmp,ig integer(i_kind) :: ipc2d(nc2d),ipc3d(nc3d),istatus ! Request ensemble-corresponding fields from control vector @@ -2271,56 +2309,58 @@ subroutine ensemble_forward_model_ad_dual_res(cvec,a_en,ibin) endif ipx=1 - im=a_en(1)%grid%im - jm=a_en(1)%grid%jm - km=a_en(1)%grid%km -!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ic2,ipic) - do n=1,n_ens - do ic3=1,nc3d - ipic=ipc3d(ic3) - do k=1,km - do j=1,jm - do i=1,im - a_en(n)%r3(ipx)%q(i,j,k)=a_en(n)%r3(ipx)%q(i,j,k) & - +work_ens%r3(ipic)%q(i,j,k)*en_perts(n,ibin)%r3(ipic)%qr4(i,j,k) - enddo - enddo - enddo - enddo - do ic2=1,nc2d - - ipic=ipc2d(ic2) - select case ( trim(StrUpCase(cvars2d(ic2))) ) - - case('PS') - - if ( pwgtflg ) then - km_tmp = km - else - km_tmp = 1 - endif - - do k=1,km_tmp - do j=1,jm - do i=1,im - a_en(n)%r3(ipx)%q(i,j,k)=a_en(n)%r3(ipx)%q(i,j,k) & - +work_ens%r2(ipic)%q(i,j)*en_perts(n,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k) - enddo - enddo - enddo - - case('SST') - - do j=1,jm - do i=1,im - a_en(n)%r3(ipx)%q(i,j,1)=a_en(n)%r3(ipx)%q(i,j,1) & - +work_ens%r2(ipic)%q(i,j)*en_perts(n,ibin)%r2(ipic)%qr4(i,j) - enddo - enddo - - end select - enddo - enddo ! enddo n_ens + im=a_en(1,1)%grid%im + jm=a_en(1,1)%grid%jm + km=a_en(1,1)%grid%km +!$omp parallel do schedule(dynamic,1) private(j,n,ic3,k,i,ic2,ipic,ig) + do ig=1,nsclgrp + do n=1,n_ens + do ic3=1,nc3d + ipic=ipc3d(ic3) + do k=1,km + do j=1,jm + do i=1,im + a_en(ig,n)%r3(ipx)%q(i,j,k)=a_en(ig,n)%r3(ipx)%q(i,j,k) & + +work_ens%r3(ipic)%q(i,j,k)*en_perts(n,ig,ibin)%r3(ipic)%qr4(i,j,k) + enddo + enddo + enddo + enddo + do ic2=1,nc2d + + ipic=ipc2d(ic2) + select case ( trim(StrUpCase(cvars2d(ic2))) ) + + case('PS') + + if ( pwgtflg ) then + km_tmp = km + else + km_tmp = 1 + endif + + do k=1,km_tmp + do j=1,jm + do i=1,im + a_en(ig,n)%r3(ipx)%q(i,j,k)=a_en(ig,n)%r3(ipx)%q(i,j,k) & + +work_ens%r2(ipic)%q(i,j)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j)*pwgt(i,j,k) + enddo + enddo + enddo + + case('SST') + + do j=1,jm + do i=1,im + a_en(ig,n)%r3(ipx)%q(i,j,1)=a_en(ig,n)%r3(ipx)%q(i,j,1) & + +work_ens%r2(ipic)%q(i,j)*en_perts(n,ig,ibin)%r2(ipic)%qr4(i,j) + enddo + enddo + + end select + enddo + enddo ! enddo n_ens + end do call gsi_bundledestroy(work_ens,istatus) if(istatus/=0) then write(6,*)trim(myname_),': trouble destroying work ens bundle' @@ -2758,6 +2798,7 @@ subroutine sqrt_beta_e_mult_cvec(grady) use kinds, only: r_kind,i_kind use gsi_4dvar, only: nsubwin use hybrid_ensemble_parameters, only: n_ens,sqrt_beta_e,grd_ens + use hybrid_ensemble_parameters, only: nsclgrp use control_vectors,only: control_vector use timermod, only: timer_ini,timer_fnl @@ -2770,7 +2811,7 @@ subroutine sqrt_beta_e_mult_cvec(grady) ! Declare local variables character(len=*),parameter::myname_=myname//'*sqrt_beta_e_mult' - integer(i_kind) :: i,j,k,ii,nn + integer(i_kind) :: i,j,k,ii,nn,ig ! Initialize timer call timer_ini('sqrt_beta_e_mult') @@ -2779,12 +2820,14 @@ subroutine sqrt_beta_e_mult_cvec(grady) !$omp parallel do schedule(dynamic,1) private(nn,k,j,i,ii) do j=1,grd_ens%lon2 do ii=1,nsubwin - do nn=1,n_ens - do k=1,nsig - do i=1,grd_ens%lat2 - grady%aens(ii,nn)%r3(1)%q(i,j,k) = sqrt_beta_e(k)*grady%aens(ii,nn)%r3(1)%q(i,j,k) - enddo - enddo + do ig=1,nsclgrp + do nn=1,n_ens + do k=1,nsig + do i=1,grd_ens%lat2 + grady%aens(ii,ig,nn)%r3(1)%q(i,j,k) = sqrt_beta_e(k)*grady%aens(ii,ig,nn)%r3(1)%q(i,j,k) + enddo + enddo + enddo enddo enddo enddo @@ -2823,6 +2866,7 @@ subroutine sqrt_beta_e_mult_bundle(aens) !$$$ end documentation block use kinds, only: r_kind,i_kind use hybrid_ensemble_parameters, only: n_ens,sqrt_beta_e,grd_ens + use hybrid_ensemble_parameters, only: nsclgrp use gsi_bundlemod, only: gsi_bundle use timermod, only: timer_ini,timer_fnl use gridmod, only: nsig @@ -2830,11 +2874,11 @@ subroutine sqrt_beta_e_mult_bundle(aens) implicit none ! Declare passed variables - type(gsi_bundle),intent(inout) :: aens(n_ens) + type(gsi_bundle),intent(inout) :: aens(nsclgrp,n_ens) ! Declare local variables character(len=*),parameter::myname_=myname//'*sqrt_beta_e_mult' - integer(i_kind) :: i,j,k,nn + integer(i_kind) :: i,j,k,nn,ig ! Initialize timer call timer_ini('sqrt_beta_e_mult') @@ -2842,12 +2886,14 @@ subroutine sqrt_beta_e_mult_bundle(aens) ! multiply by sqrt_beta_e !$omp parallel do schedule(dynamic,1) private(nn,k,j,i) do j=1,grd_ens%lon2 - do nn=1,n_ens - do k=1,nsig - do i=1,grd_ens%lat2 - aens(nn)%r3(1)%q(i,j,k) = sqrt_beta_e(k)*aens(nn)%r3(1)%q(i,j,k) - enddo - enddo + do ig=1,nsclgrp + do nn=1,n_ens + do k=1,nsig + do i=1,grd_ens%lat2 + aens(ig,nn)%r3(1)%q(i,j,k) = sqrt_beta_e(k)*aens(ig,nn)%r3(1)%q(i,j,k) + enddo + enddo + enddo enddo enddo @@ -2884,7 +2930,7 @@ subroutine init_sf_xy(jcap_in) use kinds, only: r_kind,i_kind,r_single use hybrid_ensemble_parameters,only: s_ens_hv,sp_loc,grd_ens,grd_loc,sp_ens,n_ens,p_sploc2ens,grd_sploc - use hybrid_ensemble_parameters,only: use_localization_grid + use hybrid_ensemble_parameters,only: use_localization_grid,nsclgrp use gridmod,only: use_sp_eqspace use general_specmod, only: general_init_spec_vars use constants, only: zero,half,one,two,three,rearth,pi @@ -2899,7 +2945,7 @@ subroutine init_sf_xy(jcap_in) integer(i_kind),intent(in ) :: jcap_in - integer(i_kind) i,ii,j,k,l,n,jcap,kk,nsigend + integer(i_kind) i,ii,j,k,l,n,jcap,kk,nsigend,ig real(r_kind),allocatable::g(:),gsave(:) real(r_kind) factor real(r_kind),allocatable::rkm(:),f(:,:),f0(:,:) @@ -2924,17 +2970,19 @@ subroutine init_sf_xy(jcap_in) ! make sure s_ens_hv is within allowable range ( pi*rearth*.001/jcap_in <= s_ens_hv <= 5500 ) s_ens_h_min=pi*rearth*.001_r_kind/jcap_in - do k=1,grd_ens%nsig - if(s_ens_hv(k) < s_ens_h_min) then - if(mype == 0) write(6,*)' s_ens_hv(',k,') = ',s_ens_hv(k),' km--too small, min value = ', & - s_ens_h_min,' km.' - if(mype == 0) write(6,*)' s_ens_hv(',k,') reset to min value' - s_ens_hv(k)=s_ens_h_min - else if(s_ens_hv(k) > 5500._r_kind) then - if(mype == 0) write(6,*)' s_ens_hv(',k,') = ',s_ens_hv(k),' km--too large, max value = 5500 km.' - if(mype == 0) write(6,*)' s_ens_hv(',k,') reset to max value' - s_ens_hv(k)=5500._r_kind - end if + do ig=1,nsclgrp + do k=1,grd_ens%nsig + if(s_ens_hv(k,ig) < s_ens_h_min) then + if(mype == 0) write(6,*)' s_ens_hv(',k,') = ',s_ens_hv(k,ig),' km--too small, min value = ', & + s_ens_h_min,' km.' + if(mype == 0) write(6,*)' s_ens_hv(',k,') reset to min value' + s_ens_hv(k,ig)=s_ens_h_min + else if(s_ens_hv(k,ig) > 5500._r_kind) then + if(mype == 0) write(6,*)' s_ens_hv(',k,') = ',s_ens_hv(k,ig),' km--too large, max value = 5500 km.' + if(mype == 0) write(6,*)' s_ens_hv(',k,') reset to max value' + s_ens_hv(k,ig)=5500._r_kind + end if + enddo enddo @@ -3059,98 +3107,102 @@ subroutine init_sf_xy(jcap_in) rkm(1+(grd_sploc%nlat-2)/2), & -rkm(grd_sploc%nlat-(grd_sploc%nlat-2)/2)+rkm(1+(grd_sploc%nlat-2)/2),' km' - if(.not.allocated(spectral_filter)) allocate(spectral_filter(sp_loc%nc,grd_sploc%nsig)) - if(.not.allocated(sqrt_spectral_filter)) allocate(sqrt_spectral_filter(sp_loc%nc,grd_sploc%nsig)) + if(.not.allocated(spectral_filter)) allocate(spectral_filter(nsclgrp,sp_loc%nc,grd_sploc%nsig)) + if(.not.allocated(sqrt_spectral_filter)) allocate(sqrt_spectral_filter(nsclgrp,sp_loc%nc,grd_sploc%nsig)) allocate(g(sp_loc%nc),gsave(sp_loc%nc)) allocate(pn0_npole(0:sp_loc%jcap)) allocate(ksame(grd_sploc%nsig)) - ksame=.false. - do k=2,grd_sploc%nsig - if(s_ens_hv(k) == s_ens_hv(k-1))ksame(k)=.true. - enddo - spectral_filter=zero - do k=1,grd_sploc%nsig - if(ksame(k))then - spectral_filter(:,k)=spectral_filter(:,k-1) - else - do i=1,grd_sploc%nlat - f0(i,1)=exp(-half*(rkm(i)/s_ens_hv(k))**2) - enddo - - do j=2,grd_sploc%nlon - do i=1,grd_sploc%nlat - f0(i,j)=f0(i,1) - enddo - enddo - - call general_g2s0(grd_sploc,sp_loc,g,f0) - - call general_s2g0(grd_sploc,sp_loc,g,f) - -! adjust so value at np = 1 - f=f/f(grd_sploc%nlat,1) - f0=f - call general_g2s0(grd_sploc,sp_loc,g,f) - call general_s2g0(grd_sploc,sp_loc,g,f) - if(mype == 0)then - nsigend=k - do kk=k+1,grd_sploc%nsig - if(s_ens_hv(kk) /= s_ens_hv(k))exit - nsigend=nsigend+1 - enddo - write(6,900)k,nsigend,sp_loc%jcap,s_ens_hv(k),maxval(abs(f0-f)) - 900 format(' in init_sf_xy, jcap,s_ens_hv(',i5,1x,'-',i5,'), max diff(f0-f)=', & - i10,f10.2,e20.10) - end if - -! correct spectrum by dividing by pn0_npole - gsave=g - -! obtain pn0_npole - do n=0,sp_loc%jcap - g=zero - g(2*n+1)=one - call general_s2g0(grd_sploc,sp_loc,g,f) - pn0_npole(n)=f(grd_sploc%nlat,1) - enddo - - g=zero - do n=0,sp_loc%jcap - g(2*n+1)=gsave(2*n+1)/pn0_npole(n) - enddo - -! obtain spectral_filter - - ii=0 - do l=0,sp_loc%jcap - factor=one - if(l > 0) factor=half - do n=l,sp_loc%jcap - ii=ii+1 - if(sp_loc%factsml(ii)) then - spectral_filter(ii,k)=zero - else - spectral_filter(ii,k)=factor*g(2*n+1) - end if - ii=ii+1 - if(l == 0 .or. sp_loc%factsml(ii)) then - spectral_filter(ii,k)=zero - else - spectral_filter(ii,k)=factor*g(2*n+1) - end if - enddo - enddo - end if - enddo + do ig = 1,nsclgrp + ksame=.false. + do k=2,grd_sploc%nsig + if(s_ens_hv(k,ig) == s_ens_hv(k-1,ig))ksame(k)=.true. + enddo + spectral_filter=zero + do k=1,grd_sploc%nsig + if(ksame(k))then + spectral_filter(ig,:,k)=spectral_filter(ig,:,k-1) + else + do i=1,grd_sploc%nlat + f0(i,1)=exp(-half*(rkm(i)/s_ens_hv(k,1))**2) + enddo + + do j=2,grd_sploc%nlon + do i=1,grd_sploc%nlat + f0(i,j)=f0(i,1) + enddo + enddo + + call general_g2s0(grd_sploc,sp_loc,g,f0) + + call general_s2g0(grd_sploc,sp_loc,g,f) + + ! adjust so value at np = 1 + f=f/f(grd_sploc%nlat,1) + f0=f + call general_g2s0(grd_sploc,sp_loc,g,f) + call general_s2g0(grd_sploc,sp_loc,g,f) + if(mype == 0)then + nsigend=k + do kk=k+1,grd_sploc%nsig + if(s_ens_hv(kk,ig) /= s_ens_hv(k,ig))exit + nsigend=nsigend+1 + enddo + write(6,900)k,nsigend,sp_loc%jcap,s_ens_hv(k,ig),maxval(abs(f0-f)) + 900 format(' in init_sf_xy, jcap,s_ens_hv(',i5,1x,'-',i5,'), max diff(f0-f)=', & + i10,f10.2,e20.10) + end if + + ! correct spectrum by dividing by pn0_npole + gsave=g + + ! obtain pn0_npole + do n=0,sp_loc%jcap + g=zero + g(2*n+1)=one + call general_s2g0(grd_sploc,sp_loc,g,f) + pn0_npole(n)=f(grd_sploc%nlat,1) + enddo + + g=zero + do n=0,sp_loc%jcap + g(2*n+1)=gsave(2*n+1)/pn0_npole(n) + enddo + + ! obtain spectral_filter + + ii=0 + do l=0,sp_loc%jcap + factor=one + if(l > 0) factor=half + do n=l,sp_loc%jcap + ii=ii+1 + if(sp_loc%factsml(ii)) then + spectral_filter(ig,ii,k)=zero + else + spectral_filter(ig,ii,k)=factor*g(2*n+1) + end if + ii=ii+1 + if(l == 0 .or. sp_loc%factsml(ii)) then + spectral_filter(ig,ii,k)=zero + else + spectral_filter(ig,ii,k)=factor*g(2*n+1) + end if + enddo + enddo + end if + enddo + end do deallocate(g,gsave,pn0_npole,ksame) ! Compute sqrt(spectral_filter). Ensure spectral_filter >=0 zero !$omp parallel do schedule(dynamic,1) private(k,i) - do k=1,grd_sploc%nsig - do i=1,sp_loc%nc - if (spectral_filter(i,k) < zero) spectral_filter(i,k)=zero - sqrt_spectral_filter(i,k) = sqrt(spectral_filter(i,k)) - end do + do ig=1,nsclgrp + do k=1,grd_sploc%nsig + do i=1,sp_loc%nc + if (spectral_filter(ig,i,k) < zero) spectral_filter(ig,i,k)=zero + sqrt_spectral_filter(ig,i,k) = sqrt(spectral_filter(ig,i,k)) + end do + end do end do ! assign array k_index for each processor, based on grd_loc%kbegin_loc,grd_loc%kend_loc @@ -3167,7 +3219,7 @@ subroutine init_sf_xy(jcap_in) do k=grd_loc%kbegin_loc,grd_loc%kend_loc ftest(grd_ens%nlat/2,grd_ens%nlon/2,k)=one enddo - call sf_xy(ftest,grd_loc%kbegin_loc,grd_loc%kend_loc) + call sf_xy(1,ftest,grd_loc%kbegin_loc,grd_loc%kend_loc) if(mype==0) then do j=1,grd_ens%nlon do i=1,grd_ens%nlat @@ -3183,7 +3235,7 @@ subroutine init_sf_xy(jcap_in) end subroutine init_sf_xy -subroutine sf_xy(f,k_start,k_end) +subroutine sf_xy(iscale,f,k_start,k_end) !$$$ subprogram documentation block ! . . . ! subprogram: sf_xy spectral isotropic localization for global domain @@ -3219,6 +3271,7 @@ subroutine sf_xy(f,k_start,k_end) use egrid2agrid_mod,only: g_egrid2agrid,g_egrid2agrid_ad implicit none + integer(i_kind),intent(in ) :: iscale integer(i_kind),intent(in ) :: k_start,k_end real(r_kind) ,intent(inout) :: f(grd_ens%nlat,grd_ens%nlon,k_start:max(k_start,k_end)) @@ -3230,7 +3283,7 @@ subroutine sf_xy(f,k_start,k_end) !$omp parallel do schedule(dynamic,1) private(k) do k=k_start,k_end - call sfilter(grd_ens,sp_loc,spectral_filter(1,k_index(k)),f(1,1,k)) + call sfilter(grd_ens,sp_loc,spectral_filter(iscale,:,k_index(k)),f(1,1,k)) enddo else @@ -3239,7 +3292,7 @@ subroutine sf_xy(f,k_start,k_end) !$omp parallel do schedule(dynamic,1) private(k,work) do k=k_start,k_end call g_egrid2agrid_ad(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k)) - call sfilter(grd_ens,sp_loc,spectral_filter(:,k_index(k)),f(1,1,k)) + call sfilter(grd_ens,sp_loc,spectral_filter(iscale,:,k_index(k)),f(1,1,k)) call g_egrid2agrid(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k)) enddo @@ -3248,7 +3301,7 @@ subroutine sf_xy(f,k_start,k_end) end subroutine sf_xy -subroutine sqrt_sf_xy(z,f,k_start,k_end) +subroutine sqrt_sf_xy(iscale,z,f,k_start,k_end) !$$$ subprogram documentation block ! . . . ! subprogram: sqrt_sf_xy sqrt(sf_xy) @@ -3282,6 +3335,7 @@ subroutine sqrt_sf_xy(z,f,k_start,k_end) use egrid2agrid_mod,only: g_egrid2agrid implicit none + integer(i_kind),intent(in ) :: iscale integer(i_kind),intent(in ) :: k_start,k_end real(r_kind) ,intent(in ) :: z(sp_loc%nc,k_start:max(k_start,k_end)) real(r_kind) ,intent( out) :: f(grd_ens%nlat,grd_ens%nlon,k_start:max(k_start,k_end)) @@ -3294,7 +3348,7 @@ subroutine sqrt_sf_xy(z,f,k_start,k_end) if(.not.use_localization_grid) then do k=k_start,k_end - g(:)=z(:,k)*sqrt_spectral_filter(:,k_index(k)) + g(:)=z(:,k)*sqrt_spectral_filter(iscale,:,k_index(k)) call general_s2g0(grd_ens,sp_loc,g,f(:,:,k)) enddo @@ -3302,7 +3356,7 @@ subroutine sqrt_sf_xy(z,f,k_start,k_end) vector=.false. do k=k_start,k_end - g(:)=z(:,k)*sqrt_spectral_filter(:,k_index(k)) + g(:)=z(:,k)*sqrt_spectral_filter(iscale,:,k_index(k)) call general_s2g0(grd_sploc,sp_loc,g,work) call g_egrid2agrid(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k)) enddo @@ -3312,7 +3366,7 @@ subroutine sqrt_sf_xy(z,f,k_start,k_end) end subroutine sqrt_sf_xy -subroutine sqrt_sf_xy_ad(z,f,k_start,k_end) +subroutine sqrt_sf_xy_ad(iscale, z,f,k_start,k_end) !$$$ subprogram documentation block ! . . . ! subprogram: sqrt_sf_xy_ad adjoint of sqrt_sf_xy @@ -3347,6 +3401,7 @@ subroutine sqrt_sf_xy_ad(z,f,k_start,k_end) use egrid2agrid_mod,only: g_egrid2agrid_ad implicit none + integer(i_kind),intent(in ) :: iscale integer(i_kind),intent(in ) :: k_start,k_end real(r_kind) ,intent(inout) :: z(sp_loc%nc,k_start:max(k_start,k_end)) real(r_kind) ,intent(inout) :: f(grd_ens%nlat,grd_ens%nlon,k_start:max(k_start,k_end)) @@ -3360,7 +3415,7 @@ subroutine sqrt_sf_xy_ad(z,f,k_start,k_end) do k=k_start,k_end call general_s2g0_ad(grd_ens,sp_loc,g,f(:,:,k)) - z(:,k)=g(:)*sqrt_spectral_filter(:,k_index(k)) + z(:,k)=g(:)*sqrt_spectral_filter(iscale,:,k_index(k)) enddo else @@ -3369,7 +3424,7 @@ subroutine sqrt_sf_xy_ad(z,f,k_start,k_end) do k=k_start,k_end call g_egrid2agrid_ad(p_sploc2ens,work,f(:,:,k:k),k,k,vector(k:k)) call general_s2g0_ad(grd_sploc,sp_loc,g,work) - z(:,k)=g(:)*sqrt_spectral_filter(:,k_index(k)) + z(:,k)=g(:)*sqrt_spectral_filter(iscale,:,k_index(k)) enddo end if @@ -3463,15 +3518,18 @@ subroutine bkerror_a_en(grady) use gsi_4dvar, only: nsubwin, lsqrtb use control_vectors, only: control_vector use timermod, only: timer_ini,timer_fnl - use hybrid_ensemble_parameters, only: n_ens + use hybrid_ensemble_parameters, only: n_ens,nsclgrp,alphacvarsclgrpmat + use hybrid_ensemble_parameters, only: cross_correlation_reset + use hybrid_ensemble_parameters, only: l_sum_spc_weights use gsi_bundlemod,only: gsi_bundlegetpointer implicit none ! Declare passed variables type(control_vector),intent(inout) :: grady + type(gsi_bundle),allocatable :: ebundle(:,:) ! Declare local variables - integer(i_kind) ii,ip,istatus + integer(i_kind) ii,ip,istatus,ig,ilevel,jlevel,nn if (lsqrtb) then write(6,*)'bkerror_a_en: not for use with lsqrtb' @@ -3482,7 +3540,7 @@ subroutine bkerror_a_en(grady) call timer_ini('bkerror_a_en') ! Put things in grady first since operations change input variables - call gsi_bundlegetpointer ( grady%aens(1,1),'a_en',ip,istatus) + call gsi_bundlegetpointer ( grady%aens(1,1,1),'a_en',ip,istatus) if(istatus/=0) then write(6,*)'bkerror_a_en: trouble getting pointer to ensemble CV' call stop2(317) @@ -3491,6 +3549,19 @@ subroutine bkerror_a_en(grady) ! multiply by sqrt_beta_e_mult call sqrt_beta_e_mult(grady) + if( nsclgrp > 1 .and. l_sum_spc_weights == 0 )then + allocate(ebundle(nsclgrp,n_ens)) + do ig=1,nsclgrp + do nn=1,n_ens + call gsi_bundlecreate (ebundle(ig,nn),grady%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 + end if + ! Apply variances, as well as vertical & horizontal parts of background error do ii=1,nsubwin !if(test_sqrt_localization) then @@ -3501,12 +3572,48 @@ subroutine bkerror_a_en(grady) ! deallocate(z) !else ! write(6,*)' using bkgcov_a_en_new_factorization' - call bkgcov_a_en_new_factorization(grady%aens(ii,1:n_ens)) + if( nsclgrp == 1 )then + call bkgcov_a_en_new_factorization(1,grady%aens(ii,1,1:n_ens)) + else + if (l_sum_spc_weights /= 0) then + do ig=1,nsclgrp + call bkgcov_a_en_new_factorization(ig,grady%aens(ii,ig,1:n_ens)) + enddo + else + do ig=1,nsclgrp + call ckgcov_a_en_new_factorization_ad_sdl(ig,ebundle(ig,:),grady%aens(ii,ig,1:n_ens)) + end do + + do ilevel = 1, nsclgrp + do jlevel = 1, nsclgrp + alphacvarsclgrpmat(ilevel,jlevel) = cross_correlation_reset(ilevel,jlevel) + end do + end do + + call covsclgrp_a_en_multmatrix(ebundle,alphacvarsclgrpmat) + + do ig=1,nsclgrp + call ckgcov_a_en_new_factorization_sdl(ig,ebundle(ig,:),grady%aens(ii,ig,1:n_ens)) + end do + end if + endif !end if enddo ! multiply by sqrt_beta_e_mult call sqrt_beta_e_mult(grady) + if( nsclgrp > 1 .and. l_sum_spc_weights == 0 )then + do ig=1,nsclgrp + 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) + end if ! Finalize timer call timer_fnl('bkerror_a_en') @@ -3514,7 +3621,70 @@ subroutine bkerror_a_en(grady) return end subroutine bkerror_a_en -subroutine bkgcov_a_en_new_factorization(a_en) +subroutine covsclgrp_a_en_multmatrix(a_en,cmatrix) +!$$$ subprogram documentation block +! . . . . +! subprogram: left multiplie a matrix of (ngvarloc,ngnvarloc) +! +! program history log: +! input argument list: +! a_en - +! +! output argument list: +! a_en - +! +! remarks: +! need to reconcile grid in gsi_bundle w/ grid_ens/grid_anl +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ + use hybrid_ensemble_parameters,only: n_ens,grd_ens,grd_anl + use hybrid_ensemble_parameters,only: nsclgrp + implicit none + + type(gsi_bundle),intent(inout) :: a_en(nsclgrp,n_ens) + + character(len=*),parameter::myname_='csclgrp_a_en_multmatrix' + logical nogood + integer(i_kind) i,j,k,n,im,jm,km,ic2,ic3 + integer(i_kind) :: ia_en2,ia_en + integer(i_kind) ipc3d(nc3d),ipc2d(nc2d),ipe(1),istatus + real(r_kind) cmatrix(nsclgrp,nsclgrp) + real(r_kind) rtem1(nsclgrp) + +! Check resolution consistency between static and ensemble components + +! Request ensemble-corresponding fields from control vector +! NOTE: because ensemble perturbation bundle structure is same as control +! vector, use same ipc3d and + im=a_en(1,1)%grid%im + jm=a_en(1,1)%grid%jm + km=a_en(1,1)%grid%km + ipe(1)=1 + do n=1,n_ens + do k=1,km + do j=1,jm + do i=1,im + rtem1=0.0_r_kind + do ia_en=1,nsclgrp + do ia_en2=1,nsclgrp + rtem1(ia_en)=rtem1(ia_en)+a_en(ia_en2,n)%r3(ipe(1))%q(i,j,k)*cmatrix(ia_en,ia_en2) + enddo + enddo + do ia_en=1,nsclgrp + a_en(ia_en,n)%r3(ipe(1))%q(i,j,k)=rtem1(ia_en) + enddo + end do ! i + end do !j + end do !k + end do !n_ens + +end subroutine covsclgrp_a_en_multmatrix + +subroutine bkgcov_a_en_new_factorization(iscale,a_en) !$$$ subprogram documentation block ! . . . . ! subprogram: bkgcov_a_en copy of bkgcov for hybrid ens var a_en @@ -3553,6 +3723,7 @@ subroutine bkgcov_a_en_new_factorization(a_en) ! Passed Variables ! real(r_kind),dimension(grd_loc%latlon1n,n_ens),intent(inout) :: a_en + integer(i_kind),intent(in ) :: iscale type(gsi_bundle),intent(inout) :: a_en(n_ens) ! Local Variables @@ -3592,13 +3763,13 @@ subroutine bkgcov_a_en_new_factorization(a_en) ! Apply horizontal smoother for number of horizontal scales if(regional) then iadvance=1 ; iback=2 - call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) - call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) + call new_factorization_rf_x(1,hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) + call new_factorization_rf_y(1,hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) iadvance=2 ; iback=1 - call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) - call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) + call new_factorization_rf_y(1,hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) + call new_factorization_rf_x(1,hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) else - call sf_xy(hwork,grd_loc%kbegin_loc,grd_loc%kend_loc) + call sf_xy(iscale,hwork,grd_loc%kbegin_loc,grd_loc%kend_loc) end if ! Put back onto subdomains @@ -3620,10 +3791,11 @@ subroutine bkgcov_a_en_new_factorization(a_en) return end subroutine bkgcov_a_en_new_factorization -subroutine ckgcov_a_en_new_factorization(z,a_en) +subroutine ckgcov_a_en_new_factorization(iscale,z,a_en) !$$$ subprogram documentation block ! . . . . -! subprogram: ckgcov_a_en_new_factorization sqrt(bkgcov_a_en_new_factorization) +! subprogram: ckgcov_a_en_new_factorization +! sqrt(bkgcov_a_en_new_factorization) ! prgmmr: parrish org: np22 date: 2011-06-27 ! ! abstract: make a copy of bkgcov_a_en_new_factorization and form sqrt. @@ -3632,10 +3804,12 @@ subroutine ckgcov_a_en_new_factorization(z,a_en) ! 2011-06-27 parrish, initial documentation ! ! input argument list: -! z - long vector containing sqrt control vector for ensemble extended control variable +! z - long vector containing sqrt control vector for ensemble +! extended control variable ! ! output argument list: -! a_en - bundle containing intermediate control variable after multiplication by sqrt(S), the +! a_en - bundle containing intermediate control variable after +! multiplication by sqrt(S), the ! ensemble localization correlation. ! ! attributes: @@ -3654,6 +3828,7 @@ subroutine ckgcov_a_en_new_factorization(z,a_en) implicit none ! Passed Variables + integer(i_kind),intent(in ) :: iscale type(gsi_bundle),intent(inout) :: a_en(n_ens) real(r_kind),dimension(nval_lenz_en),intent(in ) :: z @@ -3664,8 +3839,10 @@ subroutine ckgcov_a_en_new_factorization(z,a_en) ! and nhoriz = grd_loc%nlat*grd_loc%nlon for regional, ! nhoriz = (sp_loc%jcap+1)*(sp_loc%jcap+2) for global ! but internal array hwork always has -! dimension grd_loc%nlat*grd_loc%nlon*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1) -! which just happens to match up with nval_lenz_en for regional case, but not global. +! dimension +! grd_loc%nlat*grd_loc%nlon*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1) +! which just happens to match up with nval_lenz_en for regional case, but not +! global. real(r_kind),allocatable,dimension(:):: a_en_work call gsi_bundlegetpointer(a_en(1),'a_en',ipnt,istatus) @@ -3676,8 +3853,10 @@ subroutine ckgcov_a_en_new_factorization(z,a_en) if(grd_loc%kend_loc+1-grd_loc%kbegin_loc==0) then -! no work to be done on this processor, but hwork still has allocated space, since -! grd_loc%kend_alloc = grd_loc%kbegin_loc in this case, so set to zero. +! no work to be done on this processor, but hwork still has allocated space, +! since +! grd_loc%kend_alloc = grd_loc%kbegin_loc in this case, so +! set to zero. hwork=zero else ! Apply horizontal smoother for number of horizontal scales @@ -3685,11 +3864,11 @@ subroutine ckgcov_a_en_new_factorization(z,a_en) ! Make a copy of input variable z to hwork hwork=z iadvance=2 ; iback=1 - call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) - call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) + call new_factorization_rf_y(iscale,hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) + call new_factorization_rf_x(iscale,hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) else #ifdef LATER - call sqrt_sf_xy(z,hwork,grd_loc%kbegin_loc,grd_loc%kend_loc) + call sqrt_sf_xy(iscale,z,hwork,grd_loc%kbegin_loc,grd_loc%kend_loc) #else write(6,*) ' problem with ibm compiler with "use hybrid_ensemble_isotropic, only: sqrt_sf_xy"' #endif /*LATER*/ @@ -3725,10 +3904,262 @@ subroutine ckgcov_a_en_new_factorization(z,a_en) return end subroutine ckgcov_a_en_new_factorization -subroutine ckgcov_a_en_new_factorization_ad(z,a_en) +subroutine ckgcov_a_en_new_factorization_ad(iscale,z,a_en) +!$$$ subprogram documentation block +! . . . . +! subprogram: ckgcov_a_en_new_factorization_ad adjoint of +! ckgcov_a_en_new_factorization +! prgmmr: parrish org: np22 date: 2011-06-27 +! +! abstract: adjoint of ckgcov_a_en_new_factorization. Calling +! ckgcov_a_en_new_factorization_ad, +! followed by ckgcov_a_en_new_factorization is the equivalent of one +! call to +! subroutine bkgcov_a_en_new_factorization. +! +! program history log: +! 2011-06-27 parrish, initial documentation +! +! input argument list: +! z - long vector containing sqrt control vector for ensemble +! extended control variable +! a_en - bundle containing intermediate control variable after +! multiplication by sqrt(S), the +! ensemble localization correlation. +! +! output argument list: +! z - long vector containing sqrt control vector for ensemble +! extended control variable +! a_en - bundle containing intermediate control variable after +! multiplication by sqrt(S), the +! ensemble localization correlation. +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +!$$$ + use kinds, only: r_kind,i_kind + use constants, only: zero + use gridmod, only: regional + use hybrid_ensemble_parameters, only: n_ens,grd_loc + use hybrid_ensemble_parameters, only: nval_lenz_en + use general_sub2grid_mod, only: general_sub2grid + use gsi_bundlemod, only: gsi_bundle + use gsi_bundlemod, only: gsi_bundlegetpointer + + implicit none + +! Passed Variables + integer(i_kind),intent(in ) :: iscale + type(gsi_bundle),intent(inout) :: a_en(n_ens) + real(r_kind),dimension(nval_lenz_en),intent(inout) :: z + +! Local Variables + integer(i_kind) ii,k,iadvance,iback,is,ie,ipnt,istatus + real(r_kind) hwork(grd_loc%nlat*grd_loc%nlon*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1)) +!NOTE: nval_lenz_en = nhoriz*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1) +! and nhoriz = grd_loc%nlat*grd_loc%nlon for regional, +! nhoriz = (sp_loc%jcap+1)*(sp_loc%jcap+2) for global +! but internal array hwork always has +! dimension +! grd_loc%nlat*grd_loc%nlon*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1) +! which just happens to match up with nval_lenz_en for regional case, but not +! global. + real(r_kind),allocatable,dimension(:):: a_en_work + + call gsi_bundlegetpointer(a_en(1),'a_en',ipnt,istatus) + if(istatus/=0) then + write(6,*)'ckgcov_a_en_new_factorization_ad: trouble getting pointer to ensemble CV' + call stop2(999) + endif + +! Apply vertical smoother on each ensemble member + do k=1,n_ens + + iadvance=1 ; iback=2 + call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback) + + enddo + +! To avoid my having to touch the general sub2grid and grid2sub, +! get copy for ensemble components to work array + allocate(a_en_work(n_ens*a_en(1)%ndim),stat=istatus) + if(istatus/=0) then + write(6,*)'ckgcov_a_en_new_factorization_ad: trouble in alloc(a_en_work)' + call stop2(999) + endif + ii=0 + do k=1,n_ens + is=ii+1 + ie=ii+a_en(1)%ndim + a_en_work(is:ie)=a_en(k)%values(1:a_en(k)%ndim) + ii=ii+a_en(1)%ndim + enddo + +! Convert from subdomain to full horizontal field distributed among processors + call general_sub2grid(grd_loc,a_en_work,hwork) + deallocate(a_en_work) + + if(grd_loc%kend_loc+1-grd_loc%kbegin_loc==0) then +! no work to be done on this processor, but z still has allocated space, +! since +! grd_loc%kend_alloc = grd_loc%kbegin_loc in this case, so +! set to zero. + z=zero + else +! Apply horizontal smoother for number of horizontal scales + if(regional) then + iadvance=1 ; iback=2 + call new_factorization_rf_x(iscale,hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) + call new_factorization_rf_y(iscale,hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) + z=hwork + else + call sqrt_sf_xy_ad(iscale,z,hwork,grd_loc%kbegin_loc,grd_loc%kend_loc) + end if + end if + + return +end subroutine ckgcov_a_en_new_factorization_ad + + +subroutine ckgcov_a_en_new_factorization_sdl(iscale,a_en_in,a_en) !$$$ subprogram documentation block ! . . . . -! subprogram: ckgcov_a_en_new_factorization_ad adjoint of ckgcov_a_en_new_factorization +! subprogram: ckgcov_a_en_new_factorization_sdl sqrt(bkgcov_a_en_new_factorization) +! prgmmr: parrish org: np22 date: 2011-06-27 +! +! abstract: make a copy of bkgcov_a_en_new_factorization and form sqrt. +! +! program history log: +! 2011-06-27 parrish, initial documentation +! +! input argument list: +! z - long vector containing sqrt control vector for ensemble extended control variable +! +! output argument list: +! a_en - bundle containing intermediate control variable after multiplication by sqrt(S), the +! ensemble localization correlation. +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +!$$$ + use kinds, only: r_kind,i_kind + use constants, only: zero + use gridmod, only: regional + use hybrid_ensemble_parameters, only: n_ens,grd_loc,sp_loc + use hybrid_ensemble_parameters, only: nval_lenz_en + use general_sub2grid_mod, only: general_sub2grid,general_grid2sub + use gsi_bundlemod, only: gsi_bundle + use gsi_bundlemod, only: gsi_bundlegetpointer + + implicit none + +! Passed Variables + integer(i_kind),intent(in ) :: iscale + type(gsi_bundle),intent(inout) :: a_en(n_ens),a_en_in(n_ens) + +! Local Variables + integer(i_kind) ii,k,iadvance,iback,is,ie,ipnt,istatus + real(r_kind) hwork(grd_loc%nlat*grd_loc%nlon*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1)) + integer(i_kind) ilatlon,kk,i,j + real(r_kind) rwork(grd_loc%nlat,grd_loc%nlon) +!NOTE: nval_lenz_en = nhoriz*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1) +! and nhoriz = grd_loc%nlat*grd_loc%nlon for regional, +! nhoriz = (sp_loc%jcap+1)*(sp_loc%jcap+2) for global +! but internal array hwork always has +! dimension grd_loc%nlat*grd_loc%nlon*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1) +! which just happens to match up with nval_lenz_en for regional case, but not global. + real(r_kind),allocatable,dimension(:):: a_en_work + real(r_kind) :: z(nval_lenz_en) + real(r_kind) :: zsp(sp_loc%nc,grd_loc%kbegin_loc:max(grd_loc%kbegin_loc,grd_loc%kend_alloc)) + + call gsi_bundlegetpointer(a_en(1),'a_en',ipnt,istatus) + if(istatus/=0) then + write(6,*)'ckgcov_a_en_new_factorization_sdl: trouble getting pointer to ensemble CV' + call stop2(999) + endif + +! Convert from subdomain to full horizontal field distributed among processors + allocate(a_en_work(n_ens*a_en(1)%ndim),stat=istatus) + if(istatus/=0) then + write(6,*)'ckgcov_a_en_new_factorization_sdl: trouble in alloc(a_en_work)' + call stop2(999) + endif + + ii=0 + do k=1,n_ens + is=ii+1 + ie=ii+a_en(1)%ndim + a_en_work(is:ie)=a_en_in(k)%values(1:a_en_in(k)%ndim) + ii=ii+a_en(1)%ndim + enddo + + call general_sub2grid(grd_loc,a_en_work,hwork) + + ! global application requires spectral space + if( .not. regional )then + ilatlon=grd_loc%nlat*grd_loc%nlon + do kk=grd_loc%kbegin_loc,grd_loc%kend_alloc + k=kk-grd_loc%kbegin_loc + do j=1,grd_loc%nlon + do i=1,grd_loc%nlat + rwork(i,j)=hwork(k*ilatlon+(j-1)*grd_loc%nlat+i) + enddo + enddo + call general_g2s0(grd_loc,sp_loc,zsp(:,kk),rwork(:,:)) + enddo + z=reshape(zsp,(/nval_lenz_en/)) + end if + + if(grd_loc%kend_loc+1-grd_loc%kbegin_loc==0) then +! no work to be done on this processor, but hwork still has allocated space, since +! grd_loc%kend_alloc = grd_loc%kbegin_loc in this case, so set to zero. + hwork=zero + else +! Apply horizontal smoother for number of horizontal scales + if(regional) then +! Make a copy of input variable z to hwork + iadvance=2 ; iback=1 + call new_factorization_rf_y(iscale,hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) + call new_factorization_rf_x(iscale,hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) + else +#ifdef LATER + call sqrt_sf_xy(iscale,z,hwork,grd_loc%kbegin_loc,grd_loc%kend_loc) +#else + write(6,*) ' problem with ibm compiler with "use hybrid_ensemble_isotropic, only: sqrt_sf_xy"' +#endif /*LATER*/ + end if + end if + +! Put back onto subdomains + call general_grid2sub(grd_loc,hwork,a_en_work) + +! Retrieve ensemble components from long vector + ii=0 + do k=1,n_ens + is=ii+1 + ie=ii+a_en(1)%ndim + a_en(k)%values(1:a_en(k)%ndim)=a_en_work(is:ie) + ii=ii+a_en(1)%ndim + enddo + deallocate(a_en_work) + +! Apply vertical smoother on each ensemble member + do k=1,n_ens + + iadvance=2 ; iback=1 + call new_factorization_rf_z(a_en(k)%r3(ipnt)%q,iadvance,iback) + + enddo + + return +end subroutine ckgcov_a_en_new_factorization_sdl + +subroutine ckgcov_a_en_new_factorization_ad_sdl(iscale,a_en_out,a_en) +!$$$ subprogram documentation block +! . . . . +! subprogram: ckgcov_a_en_new_factorization_ad_sdl adjoint of ckgcov_a_en_new_factorization ! prgmmr: parrish org: np22 date: 2011-06-27 ! ! abstract: adjoint of ckgcov_a_en_new_factorization. Calling ckgcov_a_en_new_factorization_ad, @@ -3757,19 +4188,25 @@ subroutine ckgcov_a_en_new_factorization_ad(z,a_en) use gridmod, only: regional use hybrid_ensemble_parameters, only: n_ens,grd_loc use hybrid_ensemble_parameters, only: nval_lenz_en - use general_sub2grid_mod, only: general_sub2grid + use general_sub2grid_mod, only: general_sub2grid,general_grid2sub use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer + use egrid2agrid_mod,only: g_egrid2agrid + use hybrid_ensemble_parameters, only: sp_loc,grd_sploc implicit none ! Passed Variables - type(gsi_bundle),intent(inout) :: a_en(n_ens) - real(r_kind),dimension(nval_lenz_en),intent(inout) :: z + integer(i_kind),intent(in ) :: iscale + type(gsi_bundle),intent(inout) :: a_en(n_ens),a_en_out(n_ens) + real(r_kind),dimension(nval_lenz_en) :: z + real(r_kind),dimension(nval_lenz_en) :: z2 ! Local Variables - integer(i_kind) ii,k,iadvance,iback,is,ie,ipnt,istatus + integer(i_kind) ii,k,kk,iadvance,iback,is,ie,ipnt,istatus real(r_kind) hwork(grd_loc%nlat*grd_loc%nlon*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1)) + real(r_kind):: zsp(sp_loc%nc,grd_loc%kbegin_loc:max(grd_loc%kbegin_loc,grd_loc%kend_alloc)) + !NOTE: nval_lenz_en = nhoriz*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1) ! and nhoriz = grd_loc%nlat*grd_loc%nlon for regional, ! nhoriz = (sp_loc%jcap+1)*(sp_loc%jcap+2) for global @@ -3777,10 +4214,15 @@ subroutine ckgcov_a_en_new_factorization_ad(z,a_en) ! dimension grd_loc%nlat*grd_loc%nlon*(grd_loc%kend_alloc-grd_loc%kbegin_loc+1) ! which just happens to match up with nval_lenz_en for regional case, but not global. real(r_kind),allocatable,dimension(:):: a_en_work + real(r_kind) rwork(grd_loc%nlat,grd_loc%nlon) + real(r_kind) rwork1(grd_loc%nlat,grd_loc%nlon,1) + integer(i_kind):: i,j,ilatlon + real(r_kind) work_sp(grd_sploc%nlat,grd_sploc%nlon,1) + logical vector(grd_loc%kbegin_loc:max(grd_loc%kend_alloc,grd_loc%kbegin_loc)) call gsi_bundlegetpointer(a_en(1),'a_en',ipnt,istatus) if(istatus/=0) then - write(6,*)'ckgcov_a_en_new_factorization_ad: trouble getting pointer to ensemble CV' + write(6,*)'ckgcov_a_en_new_factorization_ad_sdl: trouble getting pointer to ensemble CV' call stop2(999) endif @@ -3796,7 +4238,7 @@ subroutine ckgcov_a_en_new_factorization_ad(z,a_en) ! get copy for ensemble components to work array allocate(a_en_work(n_ens*a_en(1)%ndim),stat=istatus) if(istatus/=0) then - write(6,*)'ckgcov_a_en_new_factorization_ad: trouble in alloc(a_en_work)' + write(6,*)'ckgcov_a_en_new_factorization_ad_sdl: trouble in alloc(a_en_work)' call stop2(999) endif ii=0 @@ -3809,7 +4251,6 @@ subroutine ckgcov_a_en_new_factorization_ad(z,a_en) ! Convert from subdomain to full horizontal field distributed among processors call general_sub2grid(grd_loc,a_en_work,hwork) - deallocate(a_en_work) if(grd_loc%kend_loc+1-grd_loc%kbegin_loc==0) then ! no work to be done on this processor, but z still has allocated space, since @@ -3819,16 +4260,43 @@ subroutine ckgcov_a_en_new_factorization_ad(z,a_en) ! Apply horizontal smoother for number of horizontal scales if(regional) then iadvance=1 ; iback=2 - call new_factorization_rf_x(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) - call new_factorization_rf_y(hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) - z=hwork + call new_factorization_rf_x(iscale,hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) + call new_factorization_rf_y(iscale,hwork,iadvance,iback,grd_loc%kend_loc+1-grd_loc%kbegin_loc) else - call sqrt_sf_xy_ad(z,hwork,grd_loc%kbegin_loc,grd_loc%kend_loc) + call sqrt_sf_xy_ad(iscale,z,hwork,grd_loc%kbegin_loc,grd_loc%kend_loc) end if end if + if( .not. regional )then + zsp=reshape(z,(/sp_loc%nc,max(grd_loc%kbegin_loc,grd_loc%kend_alloc)-grd_loc%kbegin_loc+1/)) + z2=reshape(zsp,(/nval_lenz_en/)) + ilatlon=grd_loc%nlat*grd_loc%nlon + do kk=grd_loc%kbegin_loc,grd_loc%kend_loc + + call general_s2g0(grd_loc,sp_loc,zsp(:,kk),rwork(:,:)) + k=kk-grd_loc%kbegin_loc + do j=1,grd_loc%nlon + do i=1,grd_loc%nlat + hwork(k*ilatlon+(j-1)*grd_loc%nlat+i)=rwork(i,j) + enddo + enddo + enddo + end if + + ! Put back to subdomain + call general_grid2sub(grd_loc,hwork,a_en_work) + + do k=1,n_ens + ii=(k-1)*a_en_out(1)%ndim + is=ii+1 + ie=ii+a_en_out(1)%ndim + a_en_out(k)%values(1:a_en_out(k)%ndim)=a_en_work(is:ie) + enddo + + deallocate(a_en_work) + return -end subroutine ckgcov_a_en_new_factorization_ad +end subroutine ckgcov_a_en_new_factorization_ad_sdl ! ------------------------------------------------------------------------------ ! ------------------------------------------------------------------------------ @@ -3871,6 +4339,7 @@ subroutine hybens_grid_setup use constants, only: zero,one use control_vectors, only: cvars3d,nc2d,nc3d use gridmod, only: region_lat,region_lon,region_dx,region_dy + use hybrid_ensemble_parameters, only:nsclgrp, spc_multwgt,spcwgt_params implicit none @@ -3971,6 +4440,30 @@ subroutine hybens_grid_setup end if end if + allocate(spc_multwgt(0:jcap_ens,nsclgrp)) + allocate(spcwgt_params(4,nsclgrp)) + spc_multwgt=1.0 + if(nsclgrp > 1) then + spcwgt_params(1,1)=1.E+4 + spcwgt_params(2,1)=1.0E+8 + spcwgt_params(3,1)=1.0 + spcwgt_params(4,1)=8000.0 + + spcwgt_params(1,2)=2000.0 + spcwgt_params(2,2)=2000.0 + spcwgt_params(3,2)=1.0 + spcwgt_params(4,2)=6000.0 + + spcwgt_params(1,3)=0.0 + spcwgt_params(2,3)=500.0 + spcwgt_params(3,3)=1.0 + spcwgt_params(4,3)=500.0 + endif + + if(nsclgrp > 1) then + call init_mult_spc_wgts(jcap_ens) + endif + return end subroutine hybens_grid_setup @@ -4008,16 +4501,16 @@ subroutine hybens_localization_setup nval_lenz_en,regional_ensemble_option use hybrid_ensemble_parameters, only: readin_beta,beta_s,beta_e,beta_s0,beta_e0,sqrt_beta_s,sqrt_beta_e use hybrid_ensemble_parameters, only: readin_localization,create_hybens_localization_parameters, & - vvlocal,s_ens_h,s_ens_hv,s_ens_v,s_ens_vv + vvlocal,s_ens_h,s_ens_hv,s_ens_v,s_ens_vv,nsclgrp use gsi_io, only: verbose implicit none integer(i_kind),parameter :: lunin = 47 character(len=40),parameter :: fname = 'hybens_info' - integer(i_kind) :: k,msig,istat,nz,kl + integer(i_kind) :: k,msig,istat,nz,kl,ig logical :: lexist,print_verbose - real(r_kind),allocatable:: s_ens_h_gu_x(:),s_ens_h_gu_y(:) + real(r_kind),allocatable:: s_ens_h_gu_x(:,:),s_ens_h_gu_y(:,:) print_verbose=.false. .and. mype == 0 if(verbose .and. mype == 0)print_verbose=.true. @@ -4044,8 +4537,8 @@ 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), beta_s(k), beta_e(k) + if(mype==0) write(6,101) s_ens_hv(k,1), s_ens_vv(k), beta_s(k), beta_e(k) enddo close(lunin) @@ -4060,8 +4553,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(nsclgrp,grd_loc%nsig*n_ens)) + if(.not.allocated(s_ens_h_gu_y)) allocate(s_ens_h_gu_y(nsclgrp,grd_loc%nsig*n_ens)) endif endif ! if ( readin_localization .or. readin_beta ) @@ -4091,9 +4584,11 @@ 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 + if(.not.allocated(s_ens_h_gu_x)) allocate(s_ens_h_gu_x(nsclgrp,1)) + if(.not.allocated(s_ens_h_gu_y)) allocate(s_ens_h_gu_y(nsclgrp,1)) + do ig=1,nsclgrp + s_ens_hv(:,ig) = s_ens_h(ig) + end do s_ens_vv = s_ens_v endif @@ -4103,10 +4598,10 @@ 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) + call convert_km_to_grid_units(s_ens_h_gu_x,s_ens_h_gu_y,nsclgrp,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 init_rf_x(s_ens_h_gu_x(nsclgrp,grd_loc%kbegin_loc:grd_loc%kend_alloc),kl) + call init_rf_y(s_ens_h_gu_y(nsclgrp,grd_loc%kbegin_loc:grd_loc%kend_alloc),kl) else call init_rf_x(s_ens_h_gu_x,kl) call init_rf_y(s_ens_h_gu_y,kl) @@ -4139,8 +4634,10 @@ 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' - do k=1,grd_ens%nsig - write(6,101) s_ens_hv(k), s_ens_vv(k), beta_s(k), beta_e(k) + do ig=1,nsclgrp + do k=1,grd_ens%nsig + write(6,101) s_ens_hv(k,ig), s_ens_vv(k), beta_s(k), beta_e(k) + enddo enddo endif @@ -4148,7 +4645,7 @@ subroutine hybens_localization_setup end subroutine hybens_localization_setup -subroutine convert_km_to_grid_units(s_ens_h_gu_x,s_ens_h_gu_y,nz) +subroutine convert_km_to_grid_units(s_ens_h_gu_x,s_ens_h_gu_y,nscale,nz) !$$$ subprogram documentation block ! . . . . ! subprogram: convert_km_to_grid_units @@ -4174,16 +4671,16 @@ 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: grd_loc,n_ens,s_ens_hv,nsclgrp use hybrid_ensemble_parameters, only: region_dx_ens,region_dy_ens 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) + integer(i_kind) ,intent(in ) ::nz,nscale + real(r_kind),intent( out) ::s_ens_h_gu_x(nscale,nz*n_ens),s_ens_h_gu_y(nscale,nz*n_ens) logical :: print_verbose real(r_kind) dxmax,dymax - integer(i_kind) k,n,nk + integer(i_kind) k,n,nk,ig print_verbose=.false. if(verbose) print_verbose=.true. @@ -4196,21 +4693,25 @@ subroutine convert_km_to_grid_units(s_ens_h_gu_x,s_ens_h_gu_y,nz) .001_r_kind*minval(region_dy_ens),.001_r_kind*dymax 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) + do ig=1,nsclgrp + do k=1,nz + s_ens_h_gu_x(ig,k)=s_ens_hv(k,ig)/(.001_r_kind*dxmax) + s_ens_h_gu_y(ig,k)=s_ens_hv(k,ig)/(.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,ig),s_ens_h_gu_x(ig,k),s_ens_h_gu_y(ig,k) + enddo 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 + do ig=1,nsclgrp + do n=2,n_ens + nk=(n-1)*grd_loc%nsig + do k=1,grd_loc%nsig + s_ens_h_gu_x(ig,nk+k)=s_ens_h_gu_x(k,ig) + s_ens_h_gu_y(ig,nk+k)=s_ens_h_gu_y(k,ig) + enddo + enddo enddo endif return diff --git a/src/gsi/hybrid_ensemble_parameters.f90 b/src/gsi/hybrid_ensemble_parameters.f90 index 9ca7770c4..0ed43713f 100644 --- a/src/gsi/hybrid_ensemble_parameters.f90 +++ b/src/gsi/hybrid_ensemble_parameters.f90 @@ -151,6 +151,9 @@ 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-08-29 Y. Wang, X. Wang - add parameters for multiscale DA by using scale- and variable-dependent localization +! Wang and Wang (2022ab, MWR) +! poc: xuguang.wang@ou.edu ! ! subroutines included: @@ -297,6 +300,16 @@ module hybrid_ensemble_parameters public :: sst_staticB public :: limqens + public :: nsclgrp + public :: alphacvarsclgrpmat + public :: para_covwithsclgrp + public :: spc_multwgt + public :: spcwgt_params + public :: l_sum_spc_weights + public :: smooth_scales + public :: cross_correlation_reset + public :: vdl_scale,small_loc_varlist,smooth_scales_num + logical l_hyb_ens,uv_hyb_ens,q_hyb_ens,oz_univ_static,sst_staticB logical aniso_a_en logical full_ensemble,pwgtflg @@ -316,12 +329,20 @@ module hybrid_ensemble_parameters 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),parameter::max_aens=10 + real(r_kind) s_ens_h(max_aens) + real(r_kind) smooth_scales(max_aens) + character(len=3) small_loc_varlist(100) + integer(i_kind) vdl_scale(max_aens) + real(r_kind), dimension(max_aens,max_aens) :: cross_correlation_reset + integer(i_kind) smooth_scales_num 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,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 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_hv + real(r_kind),allocatable,dimension(:) :: s_ens_vv 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 @@ -337,13 +358,20 @@ module hybrid_ensemble_parameters integer(i_kind) fv3sar_ensemble_opt character(len=512),save :: ensemble_path + real(r_kind),allocatable,dimension(:,:) :: alphacvarsclgrpmat + real(r_kind),allocatable,dimension(:,:) :: spc_multwgt + real(r_kind),allocatable,dimension(:,:) :: spcwgt_params + real (r_kind) :: para_covwithsclgrp=1 + integer(i_kind):: nsclgrp=1 + integer(i_kind):: l_sum_spc_weights=0 + ! following is for storage of ensemble perturbations: ! def en_perts - array of ensemble perturbations ! 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 @@ -412,6 +440,11 @@ subroutine init_hybrid_ensemble_parameters beta_s0=one beta_e0=-one grid_ratio_ens=one + cross_correlation_reset=1.0 + smooth_scales=-999.0 + smooth_scales_num = -999 + vdl_scale = 0 + small_loc_varlist = 'aaa' s_ens_h = 2828._r_kind ! km (this was optimal value in ! Wang, X.,D. M. Barker, C. Snyder, and T. M. Hamill, 2008: A hybrid ! ETKF.3DVAR data assimilation scheme for the WRF Model. Part II: @@ -435,15 +468,17 @@ 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,nsclgrp),s_ens_vv(grd_ens%nsig) ) 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(nsclgrp,nsclgrp)) beta_s =one beta_e =zero sqrt_beta_s=one sqrt_beta_e=zero pwgt=zero + alphacvarsclgrpmat=one end subroutine create_hybens_localization_parameters diff --git a/src/gsi/jfunc.f90 b/src/gsi/jfunc.f90 index d65b4c8fd..cfddf1f7d 100644 --- a/src/gsi/jfunc.f90 +++ b/src/gsi/jfunc.f90 @@ -737,7 +737,7 @@ subroutine set_pointer use constants, only : max_varname_length 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: l_hyb_ens,n_ens,generate_ens,grd_ens,nval_lenz_en,nsclgrp implicit none integer(i_kind) n_ensz,nval_lenz_tot,nval_lenz_enz @@ -748,7 +748,7 @@ 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_len=nval_len+nsclgrp*n_ens*nsig*grd_ens%latlon11 nval_levs_ens=nval_levs+n_ens*nsig end if nsclen=npred*jpch_rad @@ -776,7 +776,7 @@ subroutine set_pointer if(lsqrtb.and.l_hyb_ens) then n_ensz=n_ens nval_lenz_enz=nval_lenz_en - nval_lenz_tot=nval_lenz+n_ensz*nval_lenz_enz + nval_lenz_tot=nval_lenz+nsclgrp*n_ensz*nval_lenz_enz endif nclenz=nsubwin*nval_lenz_tot+nsclen+npclen+ntclen else diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90 index 1fb03d094..3dd936d94 100644 --- a/src/gsi/obsmod.F90 +++ b/src/gsi/obsmod.F90 @@ -747,7 +747,7 @@ subroutine init_obsmod_dflts if_vterminal=.false. l2rwthin =.false. if_vrobs_raw=.false. - if_model_dbz=.true. + if_model_dbz=.false. inflate_obserr=.false. whichradar="KKKK" diff --git a/src/gsi/pcgsoi.f90 b/src/gsi/pcgsoi.f90 index a4ae2431b..4d1de421b 100644 --- a/src/gsi/pcgsoi.f90 +++ b/src/gsi/pcgsoi.f90 @@ -334,7 +334,6 @@ subroutine pcgsoi() end if end if - ! 3. Calculate new norm of gradients and factors going into b calculation dprod(1) = qdot_prod_sub(gradx,grady) if(iter > 0 .and. .not. lanlerr)then diff --git a/src/gsi/read_dbz_nc.f90 b/src/gsi/read_dbz_nc.f90 index 89eebde8b..813b5d7c4 100644 --- a/src/gsi/read_dbz_nc.f90 +++ b/src/gsi/read_dbz_nc.f90 @@ -71,7 +71,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no one_tenth,r1000,r60,r60inv,r100,r400,grav_equator, & eccentricity,somigliana,grav_ratio,grav,semi_major_axis,flattening use gridmod, only: tll2xy,nsig,nlat,nlon - use obsmod, only: iadate,doradaroneob, & + use obsmod, only: iadate,doradaroneob,oneoblat,oneoblon,oneobheight,oneobradid,& mintiltdbz,maxtiltdbz,minobrangedbz,maxobrangedbz,& static_gsi_nopcp_dbz,rmesh_dbz,zmesh_dbz use hybrid_ensemble_parameters,only : l_hyb_ens @@ -380,6 +380,12 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no thislon = lon(i,j) thislat = lat(i,j) + + if(doradaroneob) then + thislat=oneoblat + thislon=oneoblon + thishgt=oneobheight + endif !-Check format of longitude and correct if necessary diff --git a/src/gsi/setupdbz.f90 b/src/gsi/setupdbz.f90 index 9bbf5ed34..f30f27f3f 100644 --- a/src/gsi/setupdbz.f90 +++ b/src/gsi/setupdbz.f90 @@ -1451,6 +1451,11 @@ subroutine setupdbz(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,radardbz_d ! Write information to diagnostic file if(radardbz_diagsave .and. ii>0 )then + if( .not. l_use_dbz_directDA )then + write(7)'dbz',nchar,nreal,ii,mype,ioff0 + write(7)cdiagbuf(1:ii),rdiagbuf(:,1:ii) + deallocate(cdiagbuf,rdiagbuf) + else write(string,600) jiter 600 format('radardbz_',i2.2) diag_file=trim(dirname) // trim(string) @@ -1476,6 +1481,7 @@ subroutine setupdbz(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,radardbz_d write(lu_diag)cdiagbuf(1:ii),rdiagbuf(:,1:ii) deallocate(cdiagbuf,rdiagbuf) close(lu_diag) + end if end if write(6,*)'mype, irefsmlobs,irejrefsmlobs are ',mype,' ',irefsmlobs, ' ',irejrefsmlobs ! close(52) !simulated obs diff --git a/src/gsi/stub_get_pseudo_ensperts.f90 b/src/gsi/stub_get_pseudo_ensperts.f90 index 7bababc0c..ca39b4b6e 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 3089cb933..bf98f3106 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 14ab869ff..db27125d7 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