subroutine bkgvar_rewgt(sfvar,vpvar,tvar,psvar,mype)
!$$$ subprogram documentation block
! . . . .
! subprogram: bkgvar_rewgt add flow dependence to variances
! prgmmr: kleist org: np22 date: 2007-07-03
!
! abstract: perform background error variance reweighting based on
! flow dependence
!
! program history log:
! 2007-07-03 kleist
! 2008-06-05 safford - rm unused uses
! 2008-09-05 lueken - merged ed's changes into q1fy09 code
! 2009-04-15 wgu - added fpsproj option
! 2009-04-21 wgu - bug fix in routine smooth2d
! 2010-03-31 treadon - replace specmod components with sp_a structure
!
! input argument list:
! sfvar - stream function variance
! vpvar - unbalanced velocity potential variance
! tvar - unbalanced temperature variance
! psvar - unbalanced surface pressure variance
! mype - integer mpi task id
!
! output argument list:
! sfvar - reweighted stream function variance
! vpvar - reweighted unbalanced velocity potential variance
! tvar - reweighted unbalanced temperature variance
! psvar - reweighted unbalanced surface pressure variance
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind,r_quad
use constants, only: ione,one,zero,two,zero_quad,tiny_r_kind
use gridmod, only: nlat,nlon,nsig,lat2,lon2
use guess_grids, only: ges_div,ges_vor,ges_tv,ges_ps,nfldsig
use mpimod, only: npe,mpi_comm_world,ierror,mpi_sum,mpi_rtype,mpi_max
use balmod, only: agvz,wgvz,bvz
use berror, only: bkgv_rewgtfct,bkgv_write,fpsproj
implicit none
! Declare passed variables
real(r_kind),dimension(lat2,lon2,nsig),intent(inout) :: sfvar,vpvar,tvar
real(r_kind),dimension(lat2,lon2) ,intent(inout) :: psvar
integer(i_kind) ,intent(in ) :: mype
! Declare local variables
real(r_kind),dimension(lat2,lon2,nsig):: bald,balt
real(r_kind),dimension(lat2,lon2,nsig):: delpsi,delchi,deltv
real(r_kind),dimension(lat2,lon2):: delps,psresc,balps
real(r_quad),dimension(nsig):: mean_dz,mean_dd,mean_dt
real(r_quad) mean_dps,count
real(r_kind),dimension(nsig,2,npe):: mean_dz0,mean_dd0,mean_dz1,mean_dd1,&
mean_dt0,mean_dt1
real(r_kind),dimension(2,npe):: mean_dps0,mean_dps1
real(r_kind),dimension(nsig):: mean_dzout,mean_ddout,mean_dtout
real(r_kind) mean_dpsout
real(r_kind),dimension(nsig):: max_dz,max_dd,max_dt,max_dz0,max_dd0,&
max_dt0,rmax_dz0,rmax_dd0,rmax_dt0
real(r_kind) max_dps,max_dps0,rmax_dps0
integer(i_kind) i,j,k,l,nsmth,mm1
! Initialize local arrays
psresc=zero
mean_dpsout=zero ; mean_dzout=zero ; mean_ddout=zero ; mean_dtout=zero
mean_dps0 =zero ; mean_dz0 =zero ; mean_dd0 =zero ; mean_dt0 =zero
mean_dps1 =zero ; mean_dz1 =zero ; mean_dd1 =zero ; mean_dt1 =zero
max_dps=zero ; max_dps0=zero
max_dz =zero ; max_dd =zero ; max_dt=zero ; max_dz0=zero
max_dd0=zero ; max_dt0 =zero
balt =zero ; bald =zero ; balps =zero
! Set count to number of global grid points in quad precision
count = float(nlat)*float(nlon)
! Set parameter for communication
mm1=mype+ione
! NOTES:
! This subroutine will only work if FGAT (more than one guess time level)
! is used. For the current operational GDAS with a 6 hour time window,
! the default is to use sigf09-sigf03 to compute the delta variables
!
! Because of the FGAT component and vorticity/divergence issues, this is
! currently set up for global only
!
! No reweighting of ozone, cloud water, moisture yet
! Get stream function and velocity potential from guess vorticity and divergence
call getpsichi(ges_vor(1,1,1,1),ges_vor(1,1,1,nfldsig),delpsi)
call getpsichi(ges_div(1,1,1,1),ges_div(1,1,1,nfldsig),delchi)
! Get delta variables
do k=1,nsig
do j=1,lon2
do i=1,lat2
deltv(i,j,k) =ges_tv(i,j,k,nfldsig)-ges_tv(i,j,k,1)
end do
end do
end do
do j=1,lon2
do i=1,lat2
delps(i,j) = ges_ps(i,j,nfldsig)-ges_ps(i,j,1)
end do
end do
! Balanced surface pressure and velocity potential from delta stream function
do k=1,nsig
do j=1,lon2
do i=1,lat2
bald(i,j,k)=bvz(i,k)*delpsi(i,j,k)
end do
end do
end do
if(fpsproj)then
do k=1,nsig
do j=1,lon2
do i=1,lat2
balps(i,j)=balps(i,j)+wgvz(i,k)*delpsi(i,j,k)
end do
end do
end do
else
do j=1,lon2
do i=1,lat2
do k=1,nsig-ione
balps(i,j)=balps(i,j)+wgvz(i,k)*delpsi(i,j,k)
end do
balps(i,j)=balps(i,j)+wgvz(i,nsig)*(delchi(i,j,1)-bald(i,j,1))
end do
end do
endif
! Balanced temperature from delta stream function
do k=1,nsig
do l=1,nsig
do j=1,lon2
do i=1,lat2
balt(i,j,l)=balt(i,j,l)+agvz(i,l,k)*delpsi(i,j,k)
end do
end do
end do
end do
! Subtract off balanced parts
do k=1,nsig
do j=1,lon2
do i=1,lat2
deltv (i,j,k) = deltv (i,j,k) - balt(i,j,k)
delchi(i,j,k) = delchi(i,j,k) - bald(i,j,k)
end do
end do
end do
do j=1,lon2
do i=1,lat2
delps(i,j) = delps(i,j) - balps(i,j)
end do
end do
! Convert to root mean square
do k=1,nsig
do j=1,lon2
do i=1,lat2
delpsi(i,j,k)=sqrt( delpsi(i,j,k)**two )
delchi(i,j,k)=sqrt( delchi(i,j,k)**two )
deltv (i,j,k)=sqrt( deltv (i,j,k)**two )
end do
end do
end do
do j=1,lon2
do i=1,lat2
delps(i,j)=sqrt( delps(i,j)**two )
end do
end do
! Smooth the delta fields before computing variance reweighting
nsmth=8_i_kind
call smooth2d(delpsi,nsig,nsmth,mype)
call smooth2d(delchi,nsig,nsmth,mype)
call smooth2d(deltv ,nsig,nsmth,mype)
call smooth2d(delps ,ione,nsmth,mype)
! Get global maximum and mean of each of the delta fields; while accounting for
! reproducibility of this kind of calculation
mean_dz =zero_quad ; mean_dd=zero_quad ; mean_dt=zero_quad
mean_dps=zero_quad
do k=1,nsig
do j=2,lon2-ione
do i=2,lat2-ione
mean_dz(k) = mean_dz(k) + delpsi(i,j,k)
mean_dd(k) = mean_dd(k) + delchi(i,j,k)
mean_dt(k) = mean_dt(k) + deltv (i,j,k)
max_dz(k)=max(max_dz(k),delpsi(i,j,k))
max_dd(k)=max(max_dd(k),delchi(i,j,k))
max_dt(k)=max(max_dt(k),deltv (i,j,k))
end do
end do
end do
do j=2,lon2-ione
do i=2,lat2-ione
mean_dps = mean_dps + delps(i,j)
max_dps = max(max_dps,delps(i,j))
end do
end do
! Break quadruple precision into two double precision arrays
do k=1,nsig
mean_dz0(k,1,mm1) = mean_dz(k)
mean_dz0(k,2,mm1) = mean_dz(k) - mean_dz0(k,1,mm1)
mean_dd0(k,1,mm1) = mean_dd(k)
mean_dd0(k,2,mm1) = mean_dd(k) - mean_dd0(k,1,mm1)
mean_dt0(k,1,mm1) = mean_dt(k)
mean_dt0(k,2,mm1) = mean_dt(k) - mean_dt0(k,1,mm1)
end do
mean_dps0(1,mm1) = mean_dps
mean_dps0(2,mm1) = mean_dps - mean_dps0(1,mm1)
! Get task specific max and mean to every task
call mpi_allreduce(mean_dz0,mean_dz1,nsig*2*npe,mpi_rtype,mpi_sum,mpi_comm_world,ierror)
call mpi_allreduce(mean_dd0,mean_dd1,nsig*2*npe,mpi_rtype,mpi_sum,mpi_comm_world,ierror)
call mpi_allreduce(mean_dt0,mean_dt1,nsig*2*npe,mpi_rtype,mpi_sum,mpi_comm_world,ierror)
call mpi_allreduce(mean_dps0,mean_dps1,2*npe,mpi_rtype,mpi_sum,mpi_comm_world,ierror)
call mpi_allreduce(max_dz,max_dz0,nsig,mpi_rtype,mpi_max,mpi_comm_world,ierror)
call mpi_allreduce(max_dd,max_dd0,nsig,mpi_rtype,mpi_max,mpi_comm_world,ierror)
call mpi_allreduce(max_dt,max_dt0,nsig,mpi_rtype,mpi_max,mpi_comm_world,ierror)
call mpi_allreduce(max_dps,max_dps0,ione,mpi_rtype,mpi_max,mpi_comm_world,ierror)
! Reintegrate quad precision number and sum over all mpi tasks
mean_dz =zero_quad ; mean_dd=zero_quad ; mean_dt=zero_quad
mean_dps=zero_quad
do i=1,npe
do k=1,nsig
mean_dz(k) = mean_dz(k) + mean_dz1(k,1,i) + mean_dz1(k,2,i)
mean_dd(k) = mean_dd(k) + mean_dd1(k,1,i) + mean_dd1(k,2,i)
mean_dt(k) = mean_dt(k) + mean_dt1(k,1,i) + mean_dt1(k,2,i)
end do
mean_dps = mean_dps + mean_dps1(1,i) + mean_dps1(2,i)
end do
! Divide by number of grid points to get the mean
do k=1,nsig
mean_dz(k)=mean_dz(k)/count
mean_dd(k)=mean_dd(k)/count
mean_dt(k)=mean_dt(k)/count
end do
mean_dps = mean_dps/count
! Load quad precision array back into double precision array for use
do k=1,nsig
mean_dzout(k)=mean_dz(k)
mean_ddout(k)=mean_dd(k)
mean_dtout(k)=mean_dt(k)
end do
mean_dpsout = mean_dps
! Take reciprocal of max values. Check for tiny values
do k=1,nsig
rmax_dz0(k)=zero
rmax_dd0(k)=zero
rmax_dt0(k)=zero
if (abs(max_dz0(k))>tiny_r_kind) rmax_dz0(k)=one/max_dz0(k)
if (abs(max_dd0(k))>tiny_r_kind) rmax_dd0(k)=one/max_dd0(k)
if (abs(max_dt0(k))>tiny_r_kind) rmax_dt0(k)=one/max_dt0(k)
end do
rmax_dps0=zero
if (abs(max_dps0)>tiny_r_kind) rmax_dps0=one/max_dps0
! Get rescaling factor for each of the variables based on factor, mean, and max
do k=1,nsig
do j=1,lon2
do i=1,lat2
sfvar(i,j,k)=sfvar(i,j,k)* &
(one + bkgv_rewgtfct*(delpsi(i,j,k)-mean_dzout(k))*rmax_dz0(k))
vpvar(i,j,k)=vpvar(i,j,k)* &
(one + bkgv_rewgtfct*(delchi(i,j,k)-mean_ddout(k))*rmax_dd0(k))
tvar(i,j,k)=tvar(i,j,k)* &
(one + bkgv_rewgtfct*(deltv (i,j,k)-mean_dtout(k))*rmax_dt0(k))
end do
end do
end do
do j=1,lon2
do i=1,lat2
psvar(i,j)=psvar(i,j)* &
(one + bkgv_rewgtfct*(delps(i,j)-mean_dpsout)*rmax_dps0)
end do
end do
! Smooth background error variances and write out grid
nsmth=8_i_kind
call smooth2d(sfvar,nsig,nsmth,mype)
call smooth2d(vpvar,nsig,nsmth,mype)
call smooth2d(tvar,nsig,nsmth,mype)
call smooth2d(psvar,ione,nsmth,mype)
if (bkgv_write) call write_bkgvars_grid(sfvar,vpvar,tvar,psvar,mype)
return
end subroutine bkgvar_rewgt
subroutine getpsichi(vordiv1,vordiv2,dpsichi)
!$$$ subprogram documentation block
! . . . .
! subprogram: getpsichi compute psi and chi on subdomains
! prgmmr: kleist org: np22 date: 2007-07-03
!
! abstract: get stream function and velocity potential from vorticity
! and divergence
!
! program history log:
! 2007-07-03 kleist
! 2008-06-05 safford - rm unused uses
! 2008-09-05 lueken - merged ed's changes into q1fy09 code
! 2010-04-01 treadon - move strip,reorder,reorder2 to gridmod
!
! input argument list:
! vor - vorticity on subdomains
! div - divergence on subdomains
!
! output argument list:
! psi - stream function on subdomains
! chi - velocity potential on subdomains
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: ione,zero
use gridmod, only: lat2,nsig,iglobal,lon1,itotsub,lon2,lat1,&
nlat,nlon,ltosi,ltosj,ltosi_s,ltosj_s,sp_a,grd_a,&
strip,reorder,reorder2
use mpimod, only: iscuv_s,ierror,mpi_comm_world,irduv_s,ircuv_s,&
isduv_s,isduv_g,iscuv_g,nnnuvlevs,nuvlevs,irduv_g,ircuv_g,mpi_rtype
implicit none
! Declare passed variables
real(r_kind),dimension(lat2,lon2,nsig),intent(in ) :: vordiv1,vordiv2
real(r_kind),dimension(lat2,lon2,nsig),intent( out) :: dpsichi
! Declare local variables
integer(i_kind) i,j,k,kk,ni1,ni2
real(r_kind),dimension(lat1,lon1,nsig):: vordivsm1,vordivsm2
real(r_kind),dimension(itotsub,nuvlevs):: work1,work2
real(r_kind),dimension(nlat,nlon):: work3
real(r_kind),dimension(sp_a%nc):: spc1
dpsichi=zero
! Zero out work arrays
do k=1,nuvlevs
do j=1,itotsub
work1(j,k)=zero
work2(j,k)=zero
end do
end do
! Strip off endpoints arrays on subdomains
call strip(vordiv1,vordivsm1,nsig)
call strip(vordiv2,vordivsm2,nsig)
vordivsm1=vordivsm2-vordivsm1
call mpi_alltoallv(vordivsm1,iscuv_g,isduv_g,&
mpi_rtype,work2,ircuv_g,irduv_g,mpi_rtype,&
mpi_comm_world,ierror)
! Reorder work arrays
call reorder(work2,nuvlevs,nnnuvlevs)
! Perform scalar g2s on work array
do k=1,nnnuvlevs
spc1=zero
work3=zero
do kk=1,iglobal
ni1=ltosi(kk); ni2=ltosj(kk)
work3(ni1,ni2)=work2(kk,k)
end do
call general_g2s0(grd_a,sp_a,spc1,work3)
! Inverse laplacian
do i=2,sp_a%ncd2
spc1(2*i-ione)=spc1(2*i-ione)/(-sp_a%enn1(i))
spc1(2*i)=spc1(2*i)/(-sp_a%enn1(i))
end do
spc1(1)=zero
spc1(2)=zero
work3=zero
call general_s2g0(grd_a,sp_a,spc1,work3)
do kk=1,itotsub
ni1=ltosi_s(kk); ni2=ltosj_s(kk)
work1(kk,k)=work3(ni1,ni2)
end do
end do !end do nuvlevs
! Reorder the work array for the mpi communication
call reorder2(work1,nuvlevs,nnnuvlevs)
! Get psi/chi back on subdomains
call mpi_alltoallv(work1,iscuv_s,isduv_s,&
mpi_rtype,dpsichi,ircuv_s,irduv_s,mpi_rtype,&
mpi_comm_world,ierror)
return
end subroutine getpsichi
subroutine smooth2d(subd,nlevs,nsmooth,mype)
!$$$ subprogram documentation block
! . . . .
! subprogram: smooth2d perform nine point smoother to interior
! prgmmr: kleist org: np22 date: 2007-07-03
!
! abstract: perform communication necessary, then apply simple nine
! point smoother to interior of domain
!
! program history log:
! 2007-07-03 kleist
!
! input argument list:
! subd - field to be smoothed on subdomains
! nlevs - number of levels to perform smoothing
! nsmooth - number of times to perform smoother
! mype - mpi integer task ik
!
! output argument list:
! subd - smoothed field
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use gridmod,only: nlat,nlon,lat2,lon2
use kinds,only: r_kind,i_kind
use constants, only: izero,ione,zero,half,one,two,three,four
implicit none
integer(i_kind),intent(in ) :: nlevs,mype,nsmooth
real(r_kind) ,intent(inout) :: subd(lat2,lon2,nlevs)
real(r_kind),dimension(nlat,nlon):: grd
real(r_kind),dimension(nlat,0:nlon+ione):: grd2
real(r_kind) corn,cent,side,temp,c2,c3,c4,rterm,sums,sumn
integer i,j,k,n,workpe
workpe=izero
! Get subdomains on the grid
! Set weights for nine point smoother
corn=0.3_r_kind
cent=one
side=half
c4=four
c3=three
c2=two
rterm = one/(cent + c4*side + c4*corn)
do k=1,nlevs
call gather_stuff2(subd(1,1,k),grd,mype,workpe)
if (mype==workpe) then
! Do nsmooth number of passes over the 9pt smoother
do n=1,nsmooth
! Load grd2 which is used in computing the smoothed fields
do j=1,nlon
do i=1,nlat
grd2(i,j)=grd(i,j)
end do
end do
! Load longitude wrapper rows
do i=1,nlat
grd2(i,0)=grd(i,nlon)
grd2(i,nlon+ione)=grd(i,1)
end do
! special treatment for near-poles
sumn=zero
sums=zero
do i=1,nlon
sumn=sumn+grd2(nlat-ione,i)
sums=sums+grd2(2,i)
end do
sumn=sumn/(real(nlon,r_kind))
sums=sums/(real(nlon,r_kind))
do i=0,nlon+ione
grd2(nlat,i)=sumn
grd2(1,i)=sums
end do
! Perform smoother on interior
do j=1,nlon
do i=2,nlat-ione
temp = cent*grd2(i,j) + side*(grd2(i+ione,j) + &
grd2(i-ione,j) + grd2(i,j+ione) + grd2(i,j-ione)) + &
corn*(grd2(i+ione,j+ione) + grd2(i+ione,j-ione) + grd2(i-ione,j-ione) + &
grd2(i-ione,j+ione))
grd(i,j) = temp*rterm
end do
end do
end do ! n smooth number of passes
end if ! mype
call scatter_stuff2(grd,subd(1,1,k),mype,workpe)
end do ! k levs
! Get back on subdomains
return
end subroutine smooth2d
subroutine gather_stuff2(f,g,mype,outpe)
!$$$ subprogram documentation block
! . . . .
! subprogram: gather_stuff2 gather subdomains onto global slabs
! prgmmr: kleist org: np22 date: 2007-07-03
!
! abstract: perform communication necessary to gather subdomains to
! global slabs
!
! program history log:
! 2007-07-03 kleist
! 2010-04-01 treadon - move strip and reorder to gridmod
!
! input argument list:
! f - field on subdomains
! mype - mpi integer task id
! outpe - task to output global slab onto
!
! output argument list:
! g - global slab on task outpe
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: ione
use gridmod, only: iglobal,itotsub,nlat,nlon,lat1,lon1,lat2,lon2,ijn,displs_g,&
ltosi,ltosj,strip,reorder
use mpimod, only: mpi_rtype,mpi_comm_world,ierror
implicit none
integer(i_kind),intent(in ) :: mype,outpe
real(r_kind) ,intent(in ) :: f(lat2,lon2)
real(r_kind) ,intent( out) :: g(nlat,nlon)
real(r_kind) fsm(lat1,lon1)
real(r_kind),allocatable:: tempa(:)
integer(i_kind) i,ii,isize,j
isize=max(iglobal,itotsub)
allocate(tempa(isize))
! Strip off endpoints of input array on subdomains
call strip(f,fsm,ione)
call mpi_gatherv(fsm,ijn(mype+ione),mpi_rtype, &
tempa,ijn,displs_g,mpi_rtype,outpe,mpi_comm_world,ierror)
call reorder(tempa,ione,ione)
do ii=1,iglobal
i=ltosi(ii)
j=ltosj(ii)
g(i,j)=tempa(ii)
end do
deallocate(tempa)
end subroutine gather_stuff2
subroutine scatter_stuff2(g,f,mype,inpe)
!$$$ subprogram documentation block
! . . . .
! subprogram: scatter_stuff2 scatter global slabs to subdomains
! prgmmr: kleist org: np22 date: 2007-07-03
!
! abstract: perform communication necessary to scatter global slabs
! onto the subdomains
!
! program history log:
! 2007-07-03 kleist
! 2008-06-05 safford - rm unused uses
! 2008-09-05 lueken - merged ed's changes into q1fy09 code
!
! input argument list:
! g - field on global slabs
! mype - mpi integer task id
! inpe - task which contains slab
!
! output argument list:
! f - field on subdomains
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: ione
use gridmod, only: iglobal,itotsub,nlat,nlon,lat2,lon2,ijn_s,displs_s,&
ltosi_s,ltosj_s
use mpimod, only: mpi_rtype,mpi_comm_world,ierror
implicit none
integer(i_kind),intent(in ) :: mype,inpe
real(r_kind) ,intent(in ) :: g(nlat,nlon)
real(r_kind) ,intent( out) :: f(lat2,lon2)
real(r_kind),allocatable:: tempa(:)
integer(i_kind) i,ii,isize,j,mm1
isize=max(iglobal,itotsub)
allocate(tempa(isize))
mm1=mype+ione
if (mype==inpe) then
do ii=1,itotsub
i=ltosi_s(ii) ; j=ltosj_s(ii)
tempa(ii)=g(i,j)
end do
end if
call mpi_scatterv(tempa,ijn_s,displs_s,mpi_rtype,&
f,ijn_s(mm1),mpi_rtype,inpe,mpi_comm_world,ierror)
deallocate(tempa)
return
end subroutine scatter_stuff2