module get_pseudo_ensperts_mod
use abstract_get_pseudo_ensperts_mod
type, extends(abstract_get_pseudo_ensperts_class) :: get_pseudo_ensperts_class
contains
procedure, pass(this) :: get_pseudo_ensperts => get_pseudo_ensperts_wrf
procedure, nopass :: read_wrf_nmm_tclib
procedure, nopass :: get_bgtc_center
procedure, nopass :: create_pseudo_enpert_blend
procedure, nopass :: pseudo_ens_e2a
end type get_pseudo_ensperts_class
contains
subroutine get_pseudo_ensperts_wrf(this,en_perts,nelen)
!$$$ subprogram documentation block
! . . . .
! subprogram: get_pseudo_ensperts
! prgmmr: mtong org: np22 date: 2011-09-06
!
! abstract: read pseudo ensemble members from TC library, and interpolate to analysis
! grid. Calculate pseudo ensemble perturbations and replace global ensemble
! perturbations in vortex region (300 km from storm center)
!
! program history log:
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use gridmod, only: half_grid,filled_grid,regional
use constants,only: zero,one,two
use mpimod, only: mpi_comm_world,ierror,mype,mpi_rtype
use kinds, only: r_kind,i_kind,r_single
use gsi_4dvar, only: nhr_assimilation
use hybrid_ensemble_parameters, only: n_ens,grd_ens,uv_hyb_ens, &
grid_ratio_ens,write_ens_sprd
use general_sub2grid_mod, only: sub2grid_info,general_sub2grid,general_grid2sub, &
general_sub2grid_create_info
use gsi_io, only: lendian_in
use hybrid_ensemble_parameters, only: n_ens,grd_ens
use control_vectors, only: cvars2d,cvars3d,nc2d,nc3d
use gsi_bundlemod, only: gsi_bundlecreate
use gsi_bundlemod, only: gsi_grid
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use gsi_bundlemod, only: gsi_bundledestroy
use gsi_bundlemod, only: gsi_gridcreate
use get_wrf_nmm_ensperts_mod, only: get_wrf_nmm_ensperts_class
use get_wrf_mass_ensperts_mod, only: get_wrf_mass_ensperts_class
implicit none
! Declare passed variables
class(get_pseudo_ensperts_class), intent(inout) :: this
type(gsi_bundle),allocatable, intent(in ) :: en_perts(:,:)
integer(i_kind), intent(in ) :: nelen
type(sub2grid_info) grd_lib,grd_tmp
real(r_kind),allocatable,dimension(:,:,:) :: u,v,tv,rh
real(r_kind),allocatable,dimension(:,:) :: ps
real(r_kind) bar_norm,sig_norm
integer(i_kind):: ic2,ic3
integer(i_kind) ipc2d(nc2d),ipc3d(nc3d)
real(r_single),pointer,dimension(:,:,:):: w3
real(r_single),pointer,dimension(:,:):: w2
real(r_kind),pointer,dimension(:,:,:):: x3
real(r_kind),pointer,dimension(:,:):: x2
type(gsi_bundle):: en_bar
type(gsi_bundle),allocatable :: lib_perts(:)
type(gsi_grid):: grid_ens
type(get_wrf_nmm_ensperts_class):: nmm_enspert
type(get_wrf_mass_ensperts_class):: mass_enspert
integer(i_kind) :: nord_blend,nn_ens1,nn_ens2
real(r_kind),allocatable,dimension(:,:) :: blend,wgt
real(r_single),allocatable,dimension(:,:) :: outwork
real(r_kind),allocatable::temp(:)
integer(i_kind) i,j,k,n,istatus,mm1
character(120) filename,infofile
! character(24) fileout
logical(4) fexist
integer(i_kind) regional_time(6)
integer(i_kind) nlon_regional,nlat_regional
real(r_single) dlmd,dphd,dlmd_ens,dphd_ens,dlmd_lib,dphd_lib
real(r_kind) dlmd2,dphd2
integer(i_kind) nlon_regional_lib,nlat_regional_lib,nsig,nlon_lib,nlat_lib
integer(i_kind) inner_vars,num_fields
logical,allocatable::vector(:)
integer(i_kind) ku,kv,kt,kq,kps
real(r_kind),allocatable,dimension(:,:,:,:) :: work,work_sub,work_reg
real(r_kind) bc_lon,bc_lat,lclon,lclat,lc_lon,lc_lat,lc_lonm,lc_latm
integer(i_kind) lon_bc,lat_bc,lon_lc,lat_lc,ratio_lon,ratio_lat
logical :: test
associate( this => this ) ! eliminates warning for unused dummy argument needed for binding
end associate
! create ensemble perturbation and ensemble mean bundle
allocate(lib_perts(n_ens))
call gsi_gridcreate(grid_ens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)
call gsi_bundlecreate(en_bar,grid_ens,'ensemble mean',istatus,names2d=cvars2d, &
names3d=cvars3d,bundle_kind=r_kind)
if(istatus/=0) then
write(6,*)' get_pseudo_ensperts: trouble creating en_bar bundle'
call stop2(999)
endif
do n=1,n_ens
call gsi_bundlecreate(lib_perts(n),grid_ens,'library perts',istatus, &
names2d=cvars2d,names3d=cvars3d,bundle_kind=r_single)
if(istatus/=0) then
write(6,*)'get_pseudo_enperts: trouble creating lib_perts bundle'
call stop2(999)
endif
end do
!
! INITIALIZE ENSEMBLE MEAN ACCUMULATORS
en_bar%values=zero
do n=1,n_ens
lib_perts(n)%valuesr4=zero
end do
! read dlmd and dphd from guess
write(filename,'("sigf",i2.2)')nhr_assimilation
open(lendian_in,file=filename,form='unformatted')
rewind lendian_in
read(lendian_in) regional_time,nlon_regional,nlat_regional,nsig, &
dlmd,dphd
close(lendian_in)
open(100,file='filelistp',form='formatted',err=30)
open(200,file='infolist',form='formatted',err=30)
rewind (100)
rewind (200)
do n=1,200
read(100,'(a120)',err=20,end=40)filename
enddo
40 nn_ens1=n-1
do n=1,200
read(200,'(a120)',err=20,end=50)infofile
enddo
50 nn_ens2=n-1
if(nn_ens1 < n_ens .or. nn_ens2 < n_ens)then
print *,'missing ensemble members or infofile'
call stop2(401)
end if
rewind (100)
read(100,'(a120)',err=20,end=60)filename
open(lendian_in,file=trim(filename),form='unformatted')
rewind lendian_in
read(lendian_in) nlon_regional_lib,nlat_regional_lib,nsig,dlmd_lib,dphd_lib
! if(mype == 0)print *,'nlon_lib, nlat_lib, nsig=', nlon_regional_lib,nlat_regional_lib,nsig
! if(mype == 0)print *,'dlmd_lib,dphd_lib=', dlmd_lib,dphd_lib
if(filled_grid) then
nlon_lib=2*nlon_regional_lib-1
nlat_lib=nlat_regional_lib
endif
if(half_grid) then
nlon_lib=nlon_regional_lib
nlat_lib=1+nlat_regional_lib/2
endif
close(lendian_in)
! if(mype == 0)print *,'nlon_lib, nlat_lib, nsig=', nlon_lib, nlat_lib, nsig
num_fields=1+4*nsig
! if(mype == 0)print *,'num_fields=', num_fields
inner_vars=1
allocate(vector(num_fields))
vector=.false.
vector(1:2*nsig)=uv_hyb_ens
call general_sub2grid_create_info(grd_lib,inner_vars,nlat_lib,nlon_lib,nsig,num_fields,regional,vector)
call general_sub2grid_create_info(grd_tmp,inner_vars,grd_ens%nlat,grd_ens%nlon,nsig,num_fields,regional,vector)
deallocate(vector)
! read background TC center longitude and latide and convert to analysis
! grid index
call get_bgtc_center(bc_lon,bc_lat)
! if(mype == 0) print *,'bc_lon,bc_lat=', bc_lon,bc_lat
if(filled_grid)then
dlmd2=dlmd
dphd2=dphd
end if
if(half_grid)then
dlmd2=two*dlmd
dphd2=two*dphd
end if
dlmd2=dlmd2*grid_ratio_ens
dphd2=dlmd2*grid_ratio_ens
dlmd_ens=dlmd*grid_ratio_ens
dphd_ens=dphd*grid_ratio_ens
allocate(blend(grd_ens%nlat,grd_ens%nlon))
blend=zero
nord_blend=4
call create_pseudo_enpert_blend(bc_lon,bc_lat,dlmd2,dphd2,nord_blend,blend)
if(mype==0)then
allocate(outwork(grd_ens%nlon,grd_ens%nlat))
outwork=zero
do j=1,grd_ens%nlon
do i=1,grd_ens%nlat
outwork(j,i)=blend(i,j)
end do
end do
call outgrads1(outwork,grd_ens%nlon,grd_ens%nlat,'blend')
deallocate(outwork)
end if
! Reorganize for eventual distribution to local domains
allocate(temp(grd_ens%itotsub))
if (mype==0) then
do k=1,grd_ens%itotsub
i=grd_ens%ltosi_s(k) ; j=grd_ens%ltosj_s(k)
temp(k)=blend(i,j)
end do
end if
allocate(wgt(grd_ens%lat2,grd_ens%lon2))
wgt=zero
mm1=mype+1
call mpi_scatterv(temp,grd_ens%ijn_s,grd_ens%displs_s,mpi_rtype,&
wgt,grd_ens%ijn_s(mm1),mpi_rtype,0,mpi_comm_world,ierror)
call nmm_enspert%grads3d(grd_ens,wgt,1,mype,'wgt')
deallocate(temp)
rewind(100)
rewind(200)
do n=1,n_ens
read(100,'(a120)',err=20,end=60)filename
filename=trim(filename)
if(mype == 0)print *,'get_pseudo_ensperts: filename=', filename
inquire(file=filename,exist=fexist)
if(.not. fexist) call stop2(401)
read(200,'(a120)',err=20,end=60)infofile
infofile=trim(infofile)
if(mype == 0)print *,'get_pseudo_ensperts: infofile=', infofile
inquire(file=infofile,exist=fexist)
if(.not. fexist) call stop2(401)
! read library tc center index from info file and find the index of the
! point that match the background TC center
open(30,file=infofile,form='formatted')
read(30,*)lclon,lclat,lc_lon,lc_lat
close(30)
! if(mype == 0)print *,'n,lc_lon,lc_lat=', n,lc_lon, lc_lat
if(filled_grid)then
lc_lon=(2*lc_lon)-1
end if
if(half_grid)then
lc_lat=1+lc_lat/2
end if
! if(mype == 0)print *,'n,lc_lon,lc_lat=', n,lc_lon, lc_lat
lc_lonm=lc_lon-(bc_lon-NINT(bc_lon))*dlmd_ens/dlmd_lib
lc_latm=lc_lat-(bc_lat-NINT(bc_lat))*dphd_ens/dphd_lib
! if(mype == 0)print *,'n,lc_lonm,lc_latm=', n,lc_lonm, lc_latm
lon_bc=NINT(bc_lon)
lat_bc=NINT(bc_lat)
lon_lc=NINT(lc_lonm)
lat_lc=NINT(lc_latm)
ratio_lon=INT(dlmd_ens/dlmd_lib)
ratio_lat=INT(dphd_ens/dphd_lib)
! test
! ratio_lon=1
! ratio_lat=1
! test
! if(mype == 0)print *,'lon_bc,lat_bc=',lon_bc,lat_bc
! if(mype == 0)print *,'lon_lc,lat_lc=',lon_lc,lat_lc
! if(mype == 0)print *,'ratio_lon,ratio_lat=',ratio_lon,ratio_lat
allocate( u(grd_lib%lat2,grd_lib%lon2,grd_lib%nsig))
allocate( v(grd_lib%lat2,grd_lib%lon2,grd_lib%nsig))
allocate( tv(grd_lib%lat2,grd_lib%lon2,grd_lib%nsig))
allocate( rh(grd_lib%lat2,grd_lib%lon2,grd_lib%nsig))
allocate( ps(grd_lib%lat2,grd_lib%lon2))
u=zero ; v=zero ; tv=zero ; rh=zero ; ps=zero
! read fields from TC library
call read_wrf_nmm_tclib(grd_lib,filename,mype,ps,u,v,tv,rh)
! write(fileout,'("tclib",i3.3)') n
! call grads3a(grd_lib,u,v,tv,rh,ps,grd_lib%nsig,mype,fileout)
allocate(work_sub(grd_lib%inner_vars,grd_lib%lat2,grd_lib%lon2,num_fields))
work_sub=zero
do k=1,grd_lib%nsig
ku=k ; kv=k+grd_lib%nsig ; kt=k+2*grd_lib%nsig ; kq=k+3*grd_lib%nsig
do j=1,grd_lib%lon2
do i=1,grd_lib%lat2
work_sub(1,i,j,ku)=u(i,j,k)
work_sub(1,i,j,kv)=v(i,j,k)
work_sub(1,i,j,kt)=tv(i,j,k)
work_sub(1,i,j,kq)=rh(i,j,k)
end do
end do
end do
deallocate(u,v,tv,rh)
kps=num_fields
do j=1,grd_lib%lon2
do i=1,grd_lib%lat2
work_sub(1,i,j,kps)=ps(i,j)
end do
end do
deallocate(ps)
allocate(work(grd_lib%inner_vars,grd_lib%nlat,grd_lib%nlon,grd_lib%kbegin_loc:grd_lib%kend_alloc))
work=zero
call general_sub2grid(grd_lib,work_sub,work)
deallocate(work_sub)
! then match to regional analysis grid
allocate(work_reg(grd_tmp%inner_vars,grd_tmp%nlat,grd_tmp%nlon,grd_tmp%kbegin_loc:grd_tmp%kend_alloc))
work_reg=zero
do k=grd_lib%kbegin_loc,grd_lib%kend_loc
call pseudo_ens_e2a(lon_bc,lat_bc,lon_lc,lat_lc,grd_lib%nlon,grd_lib%nlat,work(1,1,1,k), &
blend,ratio_lon,ratio_lat,grd_ens%nlon,grd_ens%nlat,work_reg(1,1,1,k))
end do
deallocate(work)
! next general_grid2sub to go to regional grid subdomains.
allocate(work_sub(grd_tmp%inner_vars,grd_tmp%lat2,grd_tmp%lon2,grd_tmp%num_fields))
work_sub=zero
call general_grid2sub(grd_tmp,work_reg,work_sub)
deallocate(work_reg)
allocate( u(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
allocate( v(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
allocate( tv(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
allocate( rh(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
allocate( ps(grd_ens%lat2,grd_ens%lon2))
do k=1,grd_ens%nsig
ku=k ; kv=ku+grd_ens%nsig ; kt=kv+grd_ens%nsig ; kq=kt+grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
u(i,j,k)=work_sub(1,i,j,ku)
v(i,j,k)=work_sub(1,i,j,kv)
tv(i,j,k)=work_sub(1,i,j,kt) ! now pot virtual temp
rh(i,j,k)=work_sub(1,i,j,kq) ! now rh
end do
end do
end do
kps=num_fields
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
ps(i,j)=work_sub(1,i,j,kps)
end do
end do
deallocate(work_sub)
! write(fileout,'("tclibe2a",i3.3)') n
! call grads3a(grd_ens,u,v,tv,rh,ps,grd_ens%nsig,mype,fileout)
!
! SAVE ENSEMBLE MEMBER DATA IN COLUMN VECTOR
do ic3=1,nc3d
call gsi_bundlegetpointer(lib_perts(n),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)
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'
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('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')
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
end select
end do
deallocate(u,v,tv,rh)
do ic2=1,nc2d
call gsi_bundlegetpointer(lib_perts(n),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)
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'
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
end select
end do
deallocate(ps)
end do ! end do over ensemble
close(100)
close(200)
deallocate(blend)
! Convert to mean
bar_norm = one/float(n_ens)
en_bar%values=en_bar%values*bar_norm
if(write_ens_sprd)then
call mpi_barrier(mpi_comm_world,ierror)
call mass_enspert%ens_spread_dualres_regional(mype,en_perts,nelen,en_bar)
call mpi_barrier(mpi_comm_world,ierror)
end if
test=.false.
if(test)then
allocate(u(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
allocate(v(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
allocate(tv(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
allocate(rh(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
allocate(ps(grd_ens%lat2,grd_ens%lon2))
u=zero; v=zero; tv=zero; rh=zero; ps=zero
call gsi_bundlegetpointer(en_bar,cvars3d,ipc3d,istatus)
if(istatus/=0) then
write(6,*) 'mtong: cannot find 3d pointers'
call stop2(999)
endif
call gsi_bundlegetpointer(en_bar,cvars2d,ipc2d,istatus)
if(istatus/=0) then
write(6,*) 'mtong: cannot find 2d pointers'
call stop2(999)
endif
do ic3=1,nc3d
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
u(i,j,k) = en_bar%r3(ipc3d(ic3))%q(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
v(i,j,k) = en_bar%r3(ipc3d(ic3))%q(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
tv(i,j,k) = en_bar%r3(ipc3d(ic3))%q(i,j,k)
end do
end do
end do
case('q','Q')
do k=1,grd_ens%nsig
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
rh(i,j,k) = en_bar%r3(ipc3d(ic3))%q(i,j,k)
end do
end do
end do
end select
end do
do ic2=1,nc2d
select case (trim(cvars2d(ic2)))
case('ps','PS')
do j=1,grd_ens%lon2
do i=1,grd_ens%lat2
ps(i,j) = en_bar%r2(ipc2d(ic2))%q(i,j)
end do
end do
end select
end do
call nmm_enspert%grads3a(grd_ens,u,v,tv,rh,ps,grd_ens%nsig,mype,'en_bar')
deallocate(u,v,tv,rh,ps)
end if
! Convert ensemble members to perturbations
sig_norm=sqrt(one/max(one,n_ens-one))
do n=1,n_ens
do i=1,nelen
lib_perts(n)%valuesr4(i)=(lib_perts(n)%valuesr4(i)-en_bar%values(i))*sig_norm
end do
end do
call mpi_barrier(mpi_comm_world,ierror)
call gsi_bundledestroy(en_bar,istatus)
if(istatus/=0) then
write(6,*)' in get_pseudo_ensperts: trouble destroying en_bar bundle'
call stop2(999)
endif
! merge pseudo ensemble perturbation to global ensemble perturbation
! 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)
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)
if(istatus/=0) then
write(6,*) 'get_pseudo_enperts: cannot find 2d pointers'
call stop2(999)
endif
do n=1,n_ens
do ic3=1,nc3d
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) * &
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)
end do
end do
end do
end do
do ic2=1,nc2d
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) * &
lib_perts(n)%r2(ipc2d(ic2))%qr4(i,j) + &
(one-wgt(i,j))*en_perts(n,1)%r2(ipc2d(ic2))%qr4(i,j)
end do
end do
end do
end do
do n=1,n_ens
call gsi_bundledestroy(lib_perts(n),istatus)
if(istatus/=0) then
write(6,*)' in get_pseudo_ensperts: trouble destroying lib_perts bundle'
call stop2(999)
endif
end do
deallocate(wgt,lib_perts)
! output perturbations
test=.false.
if(test)then
allocate(u(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
allocate(v(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
allocate(tv(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
allocate(rh(grd_ens%lat2,grd_ens%lon2,grd_ens%nsig))
allocate(ps(grd_ens%lat2,grd_ens%lon2))
u=zero; v=zero; tv=zero; rh=zero; ps=zero
do n=1,n_ens
write(filename,"('enpert',i3.3)") n
call gsi_bundlegetpointer(en_perts(n,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)
if(istatus/=0) then
write(6,*) 'mtong: cannot find 2d pointers'
call stop2(999)
endif
do ic3=1,nc3d
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
u(i,j,k) = en_perts(n,1)%r3(ipc3d(ic3))%qr4(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
v(i,j,k) = en_perts(n,1)%r3(ipc3d(ic3))%qr4(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
tv(i,j,k) = en_perts(n,1)%r3(ipc3d(ic3))%qr4(i,j,k)
end do
end do
end do
case('q','Q')
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)
end do
end do
end do
end select
end do
do ic2=1,nc2d
select case (trim(cvars2d(ic2)))
case('ps','PS')
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)
end do
end do
end select
end do
call nmm_enspert%grads3a(grd_ens,u,v,tv,rh,ps,grd_ens%nsig,mype,filename)
end do
deallocate(u,v,tv,rh,ps)
end if ! if(test)then
return
30 write(6,*) 'GET_PSEUDO_ENSPERTS open filelist failed '
call stop2(400)
20 write(6,*) 'GET_PSEUDO_NMM_ENSPERTS read file list failed ',n
call stop2(400)
60 write(6,*) 'GET_PSEUDO_NMM_ENSPERTS reach end of file list ',n
call stop2(400)
end subroutine get_pseudo_ensperts_wrf
subroutine read_wrf_nmm_tclib(grd,filename,mype,ps,u,v,tv,rh)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_wrf_nmm_tclib read wrf_nmm TC library
! prgmmr: tong org: np22 date: 2011-09-06
!
! program history log:
!
! input argument list:
! grd - structure variable containing information about grid
! (initialized by general_sub2grid_create_info, located in general_sub2grid_mod.f90)
! filename - input sigma file name
! mype - mpi task id
!
! output argument list:
! * - library fields
!
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,r_single,i_kind
use mpimod, only: ierror,mpi_real4,mpi_comm_world,npe
use gridmod, only: half_grid,filled_grid
use constants, only: one,ten,one_tenth,half,fv
use gsi_io, only: lendian_in
use general_sub2grid_mod, only: sub2grid_info
use hybrid_ensemble_parameters, only: q_hyb_ens
use get_wrf_nmm_ensperts_mod, only: get_wrf_nmm_ensperts_class
implicit none
! Declare passed variables here
type(sub2grid_info) ,intent(in ) :: grd
character(120) ,intent(in ) :: filename
integer(i_kind) ,intent(in ) :: mype
real(r_kind),dimension(grd%lat2,grd%lon2) ,intent( out) :: ps
real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent( out) :: u,v,tv,rh
! Declare local parameters
real(r_kind),parameter:: r0_01=0.01_r_kind
real(r_kind),allocatable,dimension(:,:,:) :: tsen,q,prsl
real(r_kind),allocatable,dimension(:,:) :: g_pd
type(get_wrf_nmm_ensperts_class) :: wrf_nmm
! Declare local variables
integer(i_kind) kt,kq,ku,kv
! NMM variable names stuck in here
! other internal variables
integer(i_kind) nlon_regional,nlat_regional,nsig
real(r_single) dlmd,dphd
real(r_single)pt,pdtop
real(r_single),allocatable:: aeta1(:),aeta2(:)
real(r_kind) pdtop_ll,pt_ll
real(r_kind),allocatable:: aeta1_ll(:),aeta2_ll(:)
real(r_single),allocatable::tempa(:)
real(r_single),allocatable::temp1(:,:)
real(r_single),allocatable::all_loc(:,:,:)
integer(i_kind),allocatable::igtype(:),jsiskip(:)
integer(i_kind) irc_s_reg(npe),ird_s_reg(npe)
integer(i_kind) ifld,im,jm,lm,num_nmm_fields
integer(i_kind) num_all_fields,num_loc_groups,num_all_pad
integer(i_kind) i,icount,icount_prev,j,k
integer(i_kind) i_0,i_pd,i_t,i_q,i_u,i_v
integer(i_kind) iderivative
real(r_kind) pd,psfc_this
integer(i_kind) ireturn
logical ice
! character(24) fileout
lm=grd%nsig
num_nmm_fields=1+4*lm
num_all_fields=num_nmm_fields*1
num_loc_groups=num_all_fields/npe
! if(mype == 0) write(6,'(" at 1 in read_wrf_nmm_tclib, lm =",i6)')lm
! if(mype == 0) write(6,'(" at 1 in read_wrf_nmm_tclib, num_nmm_fields=",i6)')num_nmm_fields
! if(mype == 0) write(6,'(" at 1 in read_wrf_nmm_tclib, num_all_fields=",i6)')num_all_fields
! if(mype == 0) write(6,'(" at 1 in read_wrf_nmm_tclib, npe =",i6)')npe
! if(mype == 0) write(6,'(" at 1 in read_wrf_nmm_tclib, num_loc_groups=",i6)')num_loc_groups
do
num_all_pad=num_loc_groups*npe
if(num_all_pad >= num_all_fields) exit
num_loc_groups=num_loc_groups+1
end do
if(mype == 0) write(6,'(" at 1 in read_wrf_nmm_tclib, num_all_pad =",i6)')num_all_pad
if(mype == 0) write(6,'(" at 1 in read_wrf_nmm_tclib, num_loc_groups=",i6)')num_loc_groups
allocate(all_loc(grd%lat2,grd%lon2,num_all_pad))
allocate(jsiskip(num_nmm_fields))
allocate(igtype(num_nmm_fields))
i=0
i=i+1 ; i_pd=i ! pd
jsiskip(i)=2
igtype(i)=1
i_t=i+1
do k=1,lm
i=i+1 ! t(k)
jsiskip(i)=0
igtype(i)=1
end do
i_q=i+1
do k=1,lm
i=i+1 ! q(k)
jsiskip(i)=0 ; igtype(i)=1
end do
i_u=i+1
do k=1,lm
i=i+1 ! u(k)
jsiskip(i)=0 ; igtype(i)=2
end do
i_v=i+1
do k=1,lm
i=i+1 ! v(k)
jsiskip(i)=0 ; igtype(i)=2
end do
do i=1,npe
irc_s_reg(i)=grd%ijn_s(mype+1)
end do
ird_s_reg(1)=0
do i=1,npe
if(i /= 1) ird_s_reg(i)=ird_s_reg(i-1)+irc_s_reg(i-1)
end do
icount=0
icount_prev=1
open(lendian_in,file=filename,form='unformatted') ; rewind lendian_in
read(lendian_in) nlon_regional,nlat_regional,nsig,dlmd,dphd,pt,pdtop
! if(mype == 0)print *,'nlon_regional,nlat_regional,nsig,dlmd,dphd,pt,pdtop'
! if(mype == 0)print *,nlon_regional,nlat_regional,nsig,dlmd,dphd,pt,pdtop
im=nlon_regional
jm=nlat_regional
allocate(aeta1(nsig),aeta2(nsig))
allocate(aeta1_ll(nsig),aeta2_ll(nsig))
allocate(temp1(im,jm))
allocate(tempa(grd%itotsub))
! Read, interpolate, and distribute NMM restart fields
do ifld=1,num_nmm_fields
icount=icount+1
if(jsiskip(ifld) > 0) then
do i=1,jsiskip(ifld)
if(i == 1)then
read(lendian_in)aeta1
else if(i == 2)then
read(lendian_in)aeta2
else
read(lendian_in)
end if
end do
end if
if(mype == mod(icount-1,npe)) then
read(lendian_in)((temp1(i,j),i=1,im),j=1,jm)
if(filled_grid) call wrf_nmm%general_fill_nmm_grid2(grd,temp1,im,jm,tempa,abs(igtype(ifld)),1,ireturn)
if(half_grid) call wrf_nmm%general_half_nmm_grid2(grd,temp1,im,jm,tempa,abs(igtype(ifld)),1,ireturn)
if(ireturn == 1)call stop2(500)
else
read(lendian_in)
end if
! Distribute to local domains everytime we have npe fields
if(mod(icount,npe) == 0.or.icount == num_all_fields) then
call mpi_alltoallv(tempa,grd%ijn_s,grd%displs_s,mpi_real4, &
all_loc(1,1,icount_prev),irc_s_reg,ird_s_reg,mpi_real4,mpi_comm_world,ierror)
icount_prev=icount+1
end if
end do
close(lendian_in)
allocate(tsen(grd%lat2,grd%lon2,grd%nsig),q(grd%lat2,grd%lon2,grd%nsig))
allocate(prsl(grd%lat2,grd%lon2,grd%nsig))
allocate(g_pd(grd%lat2,grd%lon2))
aeta1_ll=aeta1
aeta2_ll=aeta2
pdtop_ll=r0_01*pdtop
pt_ll=r0_01*pt
i_0=0
kt=i_0+i_t-1
kq=i_0+i_q-1
ku=i_0+i_u-1
kv=i_0+i_v-1
do k=1,grd%nsig
kt=kt+1
kq=kq+1
ku=ku+1
kv=kv+1
do i=1,grd%lon2
do j=1,grd%lat2
u(j,i,k) = real(all_loc(j,i,ku),r_kind)
v(j,i,k) = real(all_loc(j,i,kv),r_kind)
q(j,i,k) = real(all_loc(j,i,kq),r_kind)
tsen(j,i,k) = real(all_loc(j,i,kt),r_kind)! actually holds sensible temperature
end do
end do
end do
do i=1,grd%lon2
do j=1,grd%lat2
g_pd(j,i)=real(all_loc(j,i,i_pd),r_kind)
end do
end do
! call grads3d(grd,u,nsig,mype,'u')
! convert wrf nmm pd variable to psfc in mb, and then to log(psfc) in cb
do i=1,grd%lon2
do j=1,grd%lat2
pd=r0_01*real(all_loc(j,i,i_0+i_pd),r_kind)
psfc_this=pd+pdtop_ll+pt_ll
ps(j,i)=one_tenth*psfc_this ! convert from mb to cb
end do
end do
do k=1,grd%nsig
do i=1,grd%lon2
do j=1,grd%lat2
tv(j,i,k) = tsen(j,i,k) * (one+fv*q(j,i,k))
prsl(j,i,k)=one_tenth* &
(aeta1_ll(k)*pdtop_ll + &
aeta2_ll(k)*(ten*ps(j,i)-pdtop_ll-pt_ll) + &
pt_ll)
end do
end do
end do
if (.not.q_hyb_ens) then
ice=.true.
iderivative=0
call genqsat(rh,tsen,prsl,grd%lat2,grd%lon2,grd%nsig,ice,iderivative)
do k=1,grd%nsig
do i=1,grd%lon2
do j=1,grd%lat2
rh(j,i,k)=q(j,i,k)/rh(j,i,k)
end do
end do
end do
else
do k=1,grd%nsig
do i=1,grd%lon2
do j=1,grd%lat2
rh(j,i,k)=q(j,i,k)
end do
end do
end do
end if
deallocate(aeta1,aeta2,aeta1_ll,aeta2_ll)
deallocate(all_loc,jsiskip,igtype)
deallocate(temp1,tempa)
deallocate(tsen,q,prsl,g_pd)
end subroutine read_wrf_nmm_tclib
subroutine get_bgtc_center(bc_lon,bc_lat)
!$$$ subprogram documentation block
! . . . .
! subprogram: get_bgtc_center
!
! prgmmr: mtong org: np22 date: 2011-09-08
!
! abstract: read background TC center longitude and latitude and convert to
! analysis grid index
!
! program history log:
!
! input argument list:
!
! output argument list:
! bc_lon - TC center grid index in x direction
! bc_lat - TC center grid index in y direction
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,r_single,i_kind
use constants, only: zero,one,deg2rad
use gridmod, only: tll2xy
implicit none
real(r_kind) ,intent(inout) :: bc_lon,bc_lat
integer(i_kind) iclat,iclon
character*1 sn,ew
real(r_kind),parameter:: r360=360.0_r_kind
real(r_kind) clat,clon
logical outside
clon=zero
clat=zero
open(32,file='trak.hwrf06',form='formatted')
read(32,'(33x,I3,A1,1x,I4,A1)') iclat,sn,iclon,ew
close(32)
if(sn == 'S')then
clat=- one * float(iclat)/10.0_r_kind
else
clat=float(iclat)/10.0_r_kind
endif
if(ew == 'W')then
clon=360._r_kind - float(iclon)/10.0_r_kind
else
clon=float(iclon)/10.0_r_kind
endif
if (clon >= r360) clon=clon-r360
if (clon < zero) clon=clon+r360
clat = clat * deg2rad
clon = clon * deg2rad
call tll2xy(clon,clat,bc_lon,bc_lat,outside)
if (outside) then
print *,'tc center is outside domain'
call stop2(402)
end if
end subroutine get_bgtc_center
subroutine create_pseudo_enpert_blend(bc_lon,bc_lat,dlmd,dphd,nord_blend,blndmsk)
!$$$ subprogram documentation block
! . . . .
! subprogram: get_pseudo_enpert_blend
! prgmmr: mtong org: np22 date: 2011-09-12
!
! abstract: create blend mask for pseudo ensemble perturbations
! The blending function gradually change from 1 to 0
! witin the blending zone 150 km to 300 km from the
! TC center.
!
! program history log:
! 2011-09-12 mtong, initial documentation
!
! input argument list:
! bc_lon - background TC center x index
! bc_lat - background TC center y index
! dlmd - grid resolution in x direction
! dphd - grid resolution in y direction
! nord_blend - variable used to create variable blend mask
!
! output argument list:
! blndmsk - blend mask
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use constants, only: zero,one
use kinds, only: r_kind,i_kind,r_single
use hybrid_ensemble_parameters, only: grd_ens
use blendmod, only: blend
implicit none
real(r_kind) ,intent(in ) :: bc_lon,bc_lat
real(r_kind ) ,intent(in ) :: dlmd,dphd
integer(i_kind) ,intent(in ) :: nord_blend
real(r_kind) ,intent( out) :: blndmsk(grd_ens%nlat,grd_ens%nlon)
real(r_kind),parameter:: r300=300_r_kind
real(r_kind),parameter:: r150=150_r_kind
real(r_kind),parameter:: r111=111.12_r_kind
real(r_kind) dx,dy,dr,x,y
integer(i_kind) i,j,k,imin,jmin,imax,jmax
integer(i_kind) mm
integer(i_kind),dimension(0:40):: iblend
blndmsk=zero
dx=r111*dlmd
dy=r111*dphd
imin=max(1,INT(bc_lat-r300/dy))
jmin=max(1,INT(bc_lon-r300/dx))
imax=min(grd_ens%nlat,NINT(bc_lat+r300/dy))
jmax=min(grd_ens%nlon,NINT(bc_lon+r300/dx))
! set up blend function
mm=nord_blend
call blend(mm,iblend)
do j=jmin,jmax
do i=imin,imax
dr=sqrt((dy*(i-bc_lat))*(dy*(i-bc_lat))+(dx*(j-bc_lon))*(dx*(j-bc_lon)))
if(dr <= r150)then
blndmsk(i,j)=one
else if(dr <= r300)then
x=abs(dr-r300)/r150
y=iblend(mm)
do k=mm-1,0,-1
y=x*y+iblend(k)
end do
y=y*x**(mm+1)
blndmsk(i,j)=y
end if
end do
end do
end subroutine create_pseudo_enpert_blend
subroutine pseudo_ens_e2a(lon_bc,lat_bc,lon_lc,lat_lc,nlone,nlate,e,blend,ratio_lon, &
ratio_lat,nlona,nlata,a)
!$$$ subprogram documentation block
! . . . .
! subprogram: pseudo_ens_e2a move TC library ensemble grid to analysis grid
! prgmmr: mtong org: np22 date: 2011-09-12
!
! abstract: Move TC library ensemble grid to analysis grid, so that the
! ensemble TC centers match the background TC center
!
! program history log:
! 2011-09-12 mtong, initial documentation
!
! input argument list:
! lon_bc - background TC center x index
! lat_bc - background TC center y index
! lon_lc - ensemble member TC center x index
! lat_lc - ensemble member TC center y index
! nlone - ensemble grid x dimension
! nlate - ensemble grid y dimension
! e - ensemble member values on full ensemble global domain
! blend - blending mask
! ratio_lon - ensemble to analysis grid ratio in x direction
! ratio_lat - ensemble to analysis grid ratio in y direction
! nlona - analysis grid x dimension
! nlata - analysis grid y dimension
!
! output argument list:
! a - ensemble member values on analysis grid
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind,r_single
use constants, only: zero
implicit none
integer(i_kind) ,intent(in ) :: lon_bc,lat_bc,lon_lc,lat_lc
integer(i_kind) ,intent(in ) :: nlone,nlate,nlona,nlata
integer(i_kind) ,intent(in ) :: ratio_lon,ratio_lat
real(r_kind) ,intent(in ) :: e(nlate,nlone)
real(r_kind) ,intent(in ) :: blend(nlata,nlona)
real(r_kind) ,intent( out) :: a(nlata,nlona)
integer(i_kind) i,j,ii,jj
do j=1,nlona
jj=lon_lc+ratio_lon*(j-lon_bc)
if( jj > 0 .and. jj <= nlone )then
do i=1,nlata
ii=lat_lc+ratio_lat*(i-lat_bc)
if(blend(i,j) > zero .and. ii > 0 .and. ii <= nlate)then
a(i,j)=e(ii,jj)
end if
end do
end if
end do
end subroutine pseudo_ens_e2a
end module get_pseudo_ensperts_mod