Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

GitHub Issue NOAA-EMC/GSI#459 Add Scale/Variable/Time-Dependent Localization for EnVar #460

Merged
merged 1 commit into from
Nov 22, 2022
Merged

GitHub Issue NOAA-EMC/GSI#459 Add Scale/Variable/Time-Dependent Localization for EnVar #460

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/gsi/class_get_fv3_regional_ensperts.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
2 changes: 1 addition & 1 deletion src/gsi/class_get_pseudo_ensperts.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions src/gsi/class_get_wrf_mass_ensperts.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/gsi/class_get_wrf_nmm_ensperts.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
122 changes: 81 additions & 41 deletions src/gsi/control_vectors.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down Expand Up @@ -296,6 +296,7 @@ subroutine init_anacv
! language: f90
! machine: ibm rs/6000 sp
!
use hybrid_ensemble_parameters,only:idaen3d,idaen2d
implicit none
!character(len=*),parameter:: rcname='anavinfo.txt'
character(len=*),parameter:: rcname='anavinfo' ! filename should have extension
Expand Down Expand Up @@ -345,6 +346,7 @@ subroutine init_anacv
allocate(cvarsmd(mvars))
allocate(atsfc_sdv(mvars))
allocate(an_amp0(nvars))
allocate(idaen3d(nc3d),idaen2d(nc2d))
allocate(incvars_to_zero(nvars))
allocate(incvars_zero_strat(nvars))
incvars_to_zero(:) = 'NONE'
Expand All @@ -371,12 +373,22 @@ subroutine init_anacv
cvars2d(nc2d)=trim(adjustl(var))
nrf2_loc(nc2d)=ii ! rid of soon
as2d(nc2d)=aas
if(itracer>10) then
idaen2d(nc2d)=2
else
idaen2d(nc2d)=1
endif
else
nc3d=nc3d+1
cvars3d(nc3d)=trim(adjustl(var))
nrf3_loc(nc3d)=ii ! rid of soon
nrf_3d(ii)=.true.
as3d(nc3d)=aas
if(itracer>10) then
idaen3d(nc3d)=2
else
idaen3d(nc3d)=1
endif
endif
endif
nrf_var(ii)=trim(adjustl(var))
Expand Down Expand Up @@ -446,9 +458,11 @@ subroutine allocate_cv(ycv)
!$$$ end documentation block

use hybrid_ensemble_parameters, only: grd_ens
use hybrid_ensemble_parameters, only: naensgrp
implicit none
type(control_vector), intent( out) :: ycv
integer(i_kind) :: ii,jj,nn,ndim,ierror,n_step,n_aens
integer(i_kind) :: ig
character(len=256)::bname
character(len=max_varname_length)::ltmp(1)
type(gsi_grid) :: grid_motley
Expand Down Expand Up @@ -481,7 +495,7 @@ subroutine allocate_cv(ycv)
! If so, define grid of ensemble control vector
n_aens=0
if (l_hyb_ens) then
ALLOCATE(ycv%aens(nsubwin,n_ens))
ALLOCATE(ycv%aens(nsubwin,naensgrp,n_ens))
call GSI_GridCreate(ycv%grid_aens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)
if (lsqrtb) then
n_aens=nval_lenz_en
Expand Down Expand Up @@ -524,17 +538,19 @@ subroutine allocate_cv(ycv)
if (l_hyb_ens) then

ltmp(1)='a_en'
do nn=1,n_ens
ycv%aens(jj,nn)%values => ycv%values(ii+1:ii+n_aens)
write(bname,'(a,i3.3,a,i4.4)') 'Ensemble Control Bundle subwin-',jj,' and member-',nn
call GSI_BundleSet(ycv%aens(jj,nn),ycv%grid_aens,bname,ierror,names3d=ltmp,bundle_kind=r_kind)
if (ierror/=0) then
write(6,*)'allocate_cv: error alloc(ensemble bundle)'
call stop2(109)
endif
ndim=ndim+ycv%aens(jj,nn)%ndim

ii=ii+n_aens
do ig=1,naensgrp
do nn=1,n_ens
ycv%aens(jj,ig,nn)%values => ycv%values(ii+1:ii+n_aens)
write(bname,'(a,i3.3,a,i4.4)') 'Ensemble Control Bundle subwin-',jj,' and member-',nn
call GSI_BundleSet(ycv%aens(jj,ig,nn),ycv%grid_aens,bname,ierror,names3d=ltmp,bundle_kind=r_kind)
if (ierror/=0) then
write(6,*)'allocate_cv: error alloc(ensemble bundle)'
call stop2(109)
endif
ndim=ndim+ycv%aens(jj,ig,nn)%ndim

ii=ii+n_aens
enddo
enddo

endif
Expand Down Expand Up @@ -620,15 +636,19 @@ subroutine deallocate_cv(ycv)
!
!$$$ end documentation block

use hybrid_ensemble_parameters, only: naensgrp
implicit none
type(control_vector), intent(inout) :: ycv
integer(i_kind) :: ii,nn,ierror
integer(i_kind) :: ig

if (ycv%lallocated) then
do ii=1,nsubwin
if (l_hyb_ens) then
do nn=n_ens,1,-1
call GSI_BundleUnset(ycv%aens(ii,nn),ierror)
do ig=1,naensgrp
do nn=n_ens,1,-1
call GSI_BundleUnset(ycv%aens(ii,ig,nn),ierror)
enddo
enddo
endif
! if (beta_s0>tiny_r_kind) then
Expand Down Expand Up @@ -845,9 +865,11 @@ real(r_quad) function qdot_prod_sub(xcv,ycv)
!
!$$$ end documentation block

use hybrid_ensemble_parameters, only: naensgrp
implicit none
type(control_vector), intent(in ) :: xcv, ycv
integer(i_kind) :: ii,nn,m3d,m2d,i,j,itot
integer(i_kind) :: ig,nigtmp
real(r_quad),allocatable,dimension(:) :: partsum

qdot_prod_sub=zero_quad
Expand All @@ -858,9 +880,11 @@ real(r_quad) function qdot_prod_sub(xcv,ycv)
qdot_prod_sub=qdot_prod_sub+qdot_product( xcv%step(ii)%values(:) ,ycv%step(ii)%values(:) )
end do
if(l_hyb_ens) then
do nn=1,n_ens
do ii=1,nsubwin
qdot_prod_sub=qdot_prod_sub+qdot_product( xcv%aens(ii,nn)%values(:) ,ycv%aens(ii,nn)%values(:) )
do ig=1,naensgrp
do nn=1,n_ens
do ii=1,nsubwin
qdot_prod_sub=qdot_prod_sub+qdot_product( xcv%aens(ii,ig,nn)%values(:) ,ycv%aens(ii,ig,nn)%values(:) )
end do
end do
end do
endif
Expand All @@ -869,7 +893,7 @@ real(r_quad) function qdot_prod_sub(xcv,ycv)
m3d=xcv%step(ii)%n3d
m2d=xcv%step(ii)%n2d
itot=max(m3d,0)+max(m2d,0)
if(l_hyb_ens)itot=itot+n_ens
if(l_hyb_ens)itot=itot+naensgrp*n_ens
allocate(partsum(itot))
!$omp parallel do schedule(dynamic,1) private(i)
do i = 1,m3d
Expand All @@ -880,9 +904,12 @@ real(r_quad) function qdot_prod_sub(xcv,ycv)
partsum(m3d+i) = dplevs(xcv%step(ii)%r2(i)%q,ycv%step(ii)%r2(i)%q,ihalo=1)
enddo
if(l_hyb_ens) then
do ig=1,naensgrp
nigtmp=n_ens*(ig-1)
!$omp parallel do schedule(dynamic,1) private(i)
do i = 1,n_ens
partsum(m3d+m2d+i) = dplevs(xcv%aens(ii,i)%r3(1)%q,ycv%aens(ii,i)%r3(1)%q,ihalo=1)
do i = 1,n_ens
partsum(m3d+m2d+nigtmp+i) = dplevs(xcv%aens(ii,ig,i)%r3(1)%q,ycv%aens(ii,ig,i)%r3(1)%q,ihalo=1)
end do
end do
end if
do i=1,itot
Expand Down Expand Up @@ -933,6 +960,7 @@ subroutine qdot_prod_vars_eb(xcv,ycv,prods,eb)
!
!$$$ end documentation block

use hybrid_ensemble_parameters, only: naensgrp
implicit none
type(control_vector), intent(in ) :: xcv, ycv
character(len=*) , intent(in ) :: eb
Expand All @@ -941,6 +969,8 @@ subroutine qdot_prod_vars_eb(xcv,ycv,prods,eb)
real(r_quad) :: zz(nsubwin)
integer(i_kind) :: ii,i,nn,m3d,m2d
real(r_quad),allocatable,dimension(:) :: partsum
integer(i_kind) :: ig
integer(i_kind) ::ngtmp,nn0

prods(:)=zero_quad
zz(:)=zero_quad
Expand All @@ -953,9 +983,11 @@ subroutine qdot_prod_vars_eb(xcv,ycv,prods,eb)
end do
endif
if(trim(eb) == 'cost_e') then
do nn=1,n_ens
do ii=1,nsubwin
zz(ii)=zz(ii)+qdot_product( xcv%aens(ii,nn)%values(:) ,ycv%aens(ii,nn)%values(:) )
do ig=1,naensgrp
do nn=1,n_ens
do ii=1,nsubwin
zz(ii)=zz(ii)+qdot_product( xcv%aens(ii,ig,nn)%values(:) ,ycv%aens(ii,ig,nn)%values(:) )
end do
end do
end do
endif
Expand All @@ -981,21 +1013,25 @@ subroutine qdot_prod_vars_eb(xcv,ycv,prods,eb)
end if
if(trim(eb) == 'cost_e') then
do ii=1,nsubwin ! RTod: somebody could work in opt/zing this ...
allocate(partsum(n_ens))
allocate(partsum(n_ens*naensgrp))
do ig=1,naensgrp
ngtmp=(ig-1)*n_ens
!$omp parallel do schedule(dynamic,1) private(nn,m3d,m2d)
do nn=1,n_ens
partsum(nn) = zero_quad
m3d=xcv%aens(ii,nn)%n3d
do i = 1,m3d
partsum(nn)= partsum(nn) + dplevs(xcv%aens(ii,nn)%r3(i)%q,ycv%aens(ii,nn)%r3(i)%q,ihalo=1)
do nn=1,n_ens
nn0=nn+ngtmp
partsum(nn0) = zero_quad
m3d=xcv%aens(ii,ig,nn)%n3d
do i = 1,m3d
partsum(nn0)= partsum(nn0) + dplevs(xcv%aens(ii,ig,nn)%r3(i)%q,ycv%aens(ii,ig,nn)%r3(i)%q,ihalo=1)
enddo
m2d=xcv%aens(ii,ig,nn)%n2d
do i = 1,m2d
partsum(nn0)= partsum(nn0) + dplevs(xcv%aens(ii,ig,nn)%r2(i)%q,ycv%aens(ii,ig,nn)%r2(i)%q,ihalo=1)
enddo
enddo
m2d=xcv%aens(ii,nn)%n2d
do i = 1,m2d
partsum(nn)= partsum(nn) + dplevs(xcv%aens(ii,nn)%r2(i)%q,ycv%aens(ii,nn)%r2(i)%q,ihalo=1)
enddo
enddo
do nn=1,n_ens
zz(ii)=zz(ii)+partsum(nn)
end do
do nn=1,n_ens*naensgrp
zz(ii)=zz(ii)+partsum(nn)
end do
deallocate(partsum)
end do
Expand Down Expand Up @@ -1348,13 +1384,15 @@ subroutine random_cv(ycv,kseed)
!
!$$$ end documentation block

use hybrid_ensemble_parameters, only: naensgrp
implicit none
type(control_vector) , intent(inout) :: ycv
integer(i_kind), optional, intent(in ) :: kseed

integer(i_kind):: ii,jj,nn,iseed
integer, allocatable :: nseed(:) ! Intentionaly default integer
real(r_kind), allocatable :: zz(:)
integer(i_kind) :: ig

iseed=iadatebgn
if (present(kseed)) iseed=iseed+kseed
Expand All @@ -1381,10 +1419,12 @@ subroutine random_cv(ycv,kseed)
if (nval_lenz_en>0) then
allocate(zz(nval_lenz_en))
do nn=1,n_ens
do jj=1,nsubwin
call random_number(zz)
do ii=1,nval_lenz_en
ycv%aens(jj,nn)%values(ii) = two*zz(ii)-one
do ig=1,naensgrp
do jj=1,nsubwin
call random_number(zz)
do ii=1,nval_lenz_en
ycv%aens(jj,ig,nn)%values(ii) = two*zz(ii)-one
enddo
enddo
enddo
enddo
Expand Down
14 changes: 7 additions & 7 deletions src/gsi/cplr_get_fv3_regional_ensperts.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)

implicit none
class(get_fv3_regional_ensperts_class), intent(inout) :: this
type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:)
type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:,:)
integer(i_kind), intent(in ):: nelen
real(r_single),dimension(:,:,:),allocatable,intent(inout):: ps_bar

Expand Down Expand Up @@ -248,7 +248,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
en_bar(m)%values=zero

do n=imem_start,n_ens
en_perts(n,m)%valuesr4 = zero
en_perts(n,1,m)%valuesr4 = zero
enddo

mm1=mype+1
Expand Down Expand Up @@ -311,7 +311,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)
! SAVE ENSEMBLE MEMBER DATA IN COLUMN VECTOR
do ic3=1,nc3d

call gsi_bundlegetpointer(en_perts(n,m),trim(cvars3d(ic3)),w3,istatus)
call gsi_bundlegetpointer(en_perts(n,1,m),trim(cvars3d(ic3)),w3,istatus)
if(istatus/=0) then
write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for ensemble member ',n
call stop2(9992)
Expand Down Expand Up @@ -471,7 +471,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)

do ic2=1,nc2d

call gsi_bundlegetpointer(en_perts(n,m),trim(cvars2d(ic2)),w2,istatus)
call gsi_bundlegetpointer(en_perts(n,1,m),trim(cvars2d(ic2)),w2,istatus)
if(istatus/=0) then
write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for ensemble member ',n
call stop2(9994)
Expand Down Expand Up @@ -540,7 +540,7 @@ subroutine get_fv3_regional_ensperts_run(this,en_perts,nelen,ps_bar)

do n=imem_start,n_ens
do i=1,nelen
en_perts(n,m)%valuesr4(i)=(en_perts(n,m)%valuesr4(i)-en_bar(m)%values(i))*sig_norm
en_perts(n,1,m)%valuesr4(i)=(en_perts(n,1,m)%valuesr4(i)-en_bar(m)%values(i))*sig_norm
end do
end do

Expand Down Expand Up @@ -880,7 +880,7 @@ subroutine ens_spread_dualres_regional_fv3_regional(this,mype,en_perts,nelen)

class(get_fv3_regional_ensperts_class), intent(inout) :: this
integer(i_kind),intent(in):: mype
type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:)
type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:,:)
integer(i_kind), intent(in ):: nelen

type(gsi_bundle):: sube,suba
Expand Down Expand Up @@ -922,7 +922,7 @@ subroutine ens_spread_dualres_regional_fv3_regional(this,mype,en_perts,nelen)
do n=1,n_ens
do i=1,nelen
sube%values(i)=sube%values(i) &
+(en_perts(n,1)%valuesr4(i))*(en_perts(n,1)%valuesr4(i))
+(en_perts(n,1,1)%valuesr4(i))*(en_perts(n,1,1)%valuesr4(i))
end do
end do

Expand Down
Loading