subroutine get_nmmb_ensperts
!$$$ subprogram documentation block
! . . . .
! subprogram: get_nmmb_ensperts adaptation of get_gefs_ensperts_dualres
! prgmmr: kleist org: np22 date: 2010-01-05
!
! abstract: read ensemble members, and construct ensemble perturbations, for use
! with hybrid ensemble option.
!
! program history log:
! 2011-07-01 carley - initial adaptation for NMMB (not yet dual-res compat.)
! 2011-09-19 carley - implement single precision bundle changes
! 2017-05-12 Y. Wang and X. wang - add one option to read hydrometeors and W
! for radar DA, POC: xuguang.wang@ou.edu
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind,r_single
use gridmod, only: pt_ll,pdtop_ll,aeta2_ll,aeta1_ll
use hybrid_ensemble_parameters, only: en_perts,ps_bar,nelen
use constants,only: zero,one,one_tenth,ten
use mpimod, only: mpi_comm_world,ierror,mype
use hybrid_ensemble_parameters, only: n_ens,grd_ens,q_hyb_ens
use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d
use gsi_bundlemod, only: gsi_bundlecreate,gsi_bundleset,gsi_grid,gsi_bundle, &
gsi_bundlegetpointer,gsi_bundledestroy,gsi_gridcreate
use mpeu_util, only: getindex
implicit none
real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig):: u,v,tv,q,oz,qs,rh,tsen,prsl
real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig):: w, qr, qli, ql, dbz, dw, qi
real(r_kind),dimension(grd_ens%lat2,grd_ens%lon2):: z,ps,sst2
real(r_kind),pointer,dimension(:,:,:):: x3
real(r_single),pointer,dimension(:,:,:) :: w3
real(r_kind),pointer,dimension(:,:):: x2
real(r_single),pointer,dimension(:,:):: w2
type(gsi_bundle):: en_bar
type(gsi_grid) :: grid_ens
real(r_kind) bar_norm,sig_norm
integer(i_kind) istatus,i,ic2,ic3,j,k,n,iderivative,i_radar_qr,i_radar_qli
character(70) filename
logical ice
call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)
call gsi_bundlecreate(en_bar,grid_ens,'ensemble',istatus,names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind)
if(istatus/=0) then
write(6,*)' get_nmmb_ensperts: trouble creating en_bar bundle'
call stop2(999)
endif
do n=1,n_ens
en_perts(n,1)%valuesr4=zero
end do
en_bar%values=zero
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
! Determine if qr and qli are control variables for radar data assimilation,
i_radar_qr=0
i_radar_qli=0
i_radar_qr=getindex(cvars3d,'qr')
i_radar_qli=getindex(cvars3d,'qli')
do n=1,n_ens
write(filename,100) n !make the filename
100 format('nmmb_ens_mem',i3.3)
if (mype==0)write(6,*) 'CALL GENERAL_READ_NMMB FOR ENS FILE : ',filename
if( i_radar_qr > 0 .and. i_radar_qli > 0 )then
call general_read_nmmb_radar(grd_ens,filename,mype,z,ps,u,v,w,qr,qli,ql,qi,dbz,dw,tv,tsen,q,oz)
else
call general_read_nmmb(grd_ens,filename,mype,z,ps,u,v,tv,tsen,q,oz)
end if
! For regional application (NMMB) use the the u,v option (i.e. uv_hyb_ens)
! Compute RH
! get 3d pressure at layer midpoints
! using code adapted from subroutine load_prsges for nmmb
! (in guess_grids.F90)
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
prsl(i,j,k)=one_tenth* &
(aeta1_ll(k)*pdtop_ll + &
aeta2_ll(k)*(ten*ps(i,j)-pdtop_ll-pt_ll) + &
pt_ll)
end do
end do
end do
if (.not.q_hyb_ens) then
ice=.true.
iderivative=0
call genqsat(qs,tsen(1,1,1),prsl(1,1,1),grd_ens%lat2,grd_ens%lon2,grd_ens%nsig,ice,iderivative)
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
rh(i,j,k) = q(i,j,k)/qs(i,j,k)
end do
end do
end do
end if
do ic3=1,nc3d
call gsi_bundlegetpointer(en_perts(n,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)
end if
call gsi_bundlegetpointer(en_bar,trim(cvars3d(ic3)),x3,istatus)
if(istatus/=0) then
write(6,*)' error retrieving pointer to ',trim(cvars3d(ic3)),' for en_bar in get_nmmb_ensperts'
call stop2(999)
end if
select case (trim(cvars3d(ic3)))
case('sf','SF')
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w3(i,j,k) = u(i,j,k)
x3(i,j,k)=x3(i,j,k)+u(i,j,k)
end do
end do
end do
case('vp','VP')
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w3(i,j,k) = v(i,j,k)
x3(i,j,k)=x3(i,j,k)+v(i,j,k)
end do
end do
end do
case('w','W')
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w3(i,j,k) = w(i,j,k)
x3(i,j,k)=x3(i,j,k)+w(i,j,k)
end do
end do
end do
case('dw','DW')
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w3(i,j,k) = dw(i,j,k)
x3(i,j,k)=x3(i,j,k)+dw(i,j,k)
end do
end do
end do
case('t','T')
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w3(i,j,k) = tv(i,j,k)
x3(i,j,k)=x3(i,j,k)+tv(i,j,k)
end do
end do
end do
case('q','Q')
if (.not.q_hyb_ens) then ! use RH
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w3(i,j,k) = rh(i,j,k)
x3(i,j,k)=x3(i,j,k)+rh(i,j,k)
end do
end do
end do
else ! use Q
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w3(i,j,k) = q(i,j,k)
x3(i,j,k)=x3(i,j,k)+q(i,j,k)
end do
end do
end do
end if
case('qr','QR')
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w3(i,j,k) = qr(i,j,k)
x3(i,j,k)=x3(i,j,k)+qr(i,j,k)
end do
end do
end do
case('qli','QLI')
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w3(i,j,k) = qli(i,j,k)
x3(i,j,k)=x3(i,j,k)+qli(i,j,k)
end do
end do
end do
case('qi','QI')
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w3(i,j,k) = qi(i,j,k)
x3(i,j,k)=x3(i,j,k)+qi(i,j,k)
end do
end do
end do
case('ql','QL')
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w3(i,j,k) = ql(i,j,k)
x3(i,j,k)=x3(i,j,k)+ql(i,j,k)
end do
end do
end do
case('dbz','DBZ')
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w3(i,j,k) = dbz(i,j,k)
x3(i,j,k)=x3(i,j,k)+dbz(i,j,k)
end do
end do
end do
case('oz','OZ')
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w3(i,j,k) = oz(i,j,k)
x3(i,j,k)=x3(i,j,k)+oz(i,j,k)
end do
end do
end do
end select
end do
do ic2=1,nc2d
call gsi_bundlegetpointer(en_perts(n,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)
end if
call gsi_bundlegetpointer(en_bar,trim(cvars2d(ic2)),x2,istatus)
if(istatus/=0) then
write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for en_bar in get_nmmb_ensperts'
call stop2(999)
end if
select case (trim(cvars2d(ic2)))
case('ps','PS')
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w2(i,j) = ps(i,j)
x2(i,j)=x2(i,j)+ps(i,j)
end do
end do
case('sst','SST')
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
w2(i,j) = sst2(i,j)
x2(i,j)=x2(i,j)+sst2(i,j)
end do
end do
end select
end do
end do ! end do over ensemble
! Convert to mean
bar_norm = one/float(n_ens)
en_bar%values=en_bar%values*bar_norm
! Copy pbar to module array. ps_bar may be needed for vertical localization
! in terms of scale heights/normalized p/p
do ic2=1,nc2d
if(trim(cvars2d(ic2)) == 'ps'.or.trim(cvars2d(ic2)) == 'PS') then
call gsi_bundlegetpointer(en_bar,trim(cvars2d(ic2)),x2,istatus)
if(istatus/=0) then
write(6,*)' error retrieving pointer to ',trim(cvars2d(ic2)),' for en_bar in get_nmmb_ensperts'
call stop2(999)
end if
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
ps_bar(i,j,1)=x2(i,j)
end do
end do
exit
end if
end do
call mpi_barrier(mpi_comm_world,ierror)
! Convert ensemble members to perturbations
sig_norm=sqrt(one/max(one,n_ens-one))
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
end do
end do
call gsi_bundledestroy(en_bar,istatus)
if(istatus/=0) then
write(6,*)' in get_nmmb_ensperts: trouble destroying en_bar bundle in get_nmmb_ensperts'
call stop2(999)
endif
if (mype==0)write(6,*) 'get_nmmb_ensperts DONE'
return
end subroutine get_nmmb_ensperts