module stpjcmod
!$$$ module documentation block
! . . . .
! module: stpjcmod module for weak constraint stp routines
! pgrmmr: kleist
!
! abstract: module for Jc step routines
!
! program history log:
! 2012-01-21 kleist - consolidation of Jc step routines into single module
! 2014-03-19 pondeca - add stepzise calculation for wspd10m weak constraint term
! 2014-05-07 pondeca - add stepzise calculation for howv weak constraint term
! 2014-06-17 carley/zhu - add stepzise calculation for lcbas weak constraint term
! 2015-07-10 pondeca - add stepzise calculation for cldch weak constraint term
! 2019-03-05 martin - update stplimq to weight factqmin/max by latitude
! 2019-03-14 eliu - add stplimqc to constraint negative hydrometeors
!
! subroutines included:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind,r_quad
use constants, only: zero,two,one,half,zero_quad,one_quad,two_quad
use gsi_bundlemod, only: gsi_bundle,gsi_bundlegetpointer
use gsi_metguess_mod, only: gsi_metguess_bundle
implicit none
PRIVATE
PUBLIC stplimqc,stplimq,stplimg,stplimp,stplimv,stplimw10m,stplimhowv,stplimcldch,stpliml,stpjcdfi,stpjcpdry
contains
subroutine stplimq(rval,sval,sges,outmin,outmax,nstep,itbin)
!$$$ subprogram documentation block
! . . . .
! subprogram: stplimq calculate penalty and stepsize for limit of q
! prgmmr: derber org: np23 date: 1996-11-19
!
! abstract: calculate stepsize contribution and penalty for limiting q
!
! program history log:
! 1996-11-19 derber
! 1998-07-10 weiyu yang
! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version
! 2004-03-15 kleist, d., derber, j., treadon, r., use negative q only
! 2004-06-02 kleist, add penalty for excess moisture
! 2004-07-29 treadon - add only to module use, add intent in/out
! 2004-11-22 derber - modify for openMP
! 2006-09-18 derber - modify output b1 and b3
! 2007-06-04 derber - use quad precision to get reproducability over number of processors
! 2007-06-04 derber - use quad precision to get reproducability over number of processors
! 2008-08-14 derber - optimize
! 2010-05-13 todling - update to use gsi_bundle
! 2010-07-10 todling - merge w/ r8741 (trunk); qx(:)->qx (who made the change?)
! 2011-12-27 kleist - add bins for 4d capability (4densvar option)
! 2019-03-05 martin - update to weight factqmin/max by latitude
!
! input argument list:
! rq - search direction
! sq - increment in grid space
! sges - step size estimates (4)
! nstep - number of step size estimates if == 0 then just do outer loop
! itbin - observation bin number (time level)
!
! output argument list:
! outmin(1:nstep) - current penalty for negative q sges(1:nstep)
! outmax(1:nstep) - current penalty for excess q sges(1:nstep)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use gridmod, only: lat1,lon1,nsig,istart,wgtfactlats
use jfunc, only: factqmin,factqmax
use guess_grids, only: ges_qsat
use mpimod, only: mype
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: nstep,itbin
real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges
real(r_quad),dimension(max(1,nstep)),intent( out) :: outmin,outmax
type(gsi_bundle) ,intent(in ) :: rval,sval
! Declare local variables
integer(i_kind) i,j,k,kk,ier,istatus,ii,mm1
real(r_kind) q,qx
real(r_kind),pointer,dimension(:,:,:) :: rq,sq
real(r_kind),pointer,dimension(:,:,:) :: ges_q_it=>NULL()
outmin=zero_quad; outmax=zero_quad
if (factqmin==zero .and. factqmax==zero) return
mm1=mype+1
! Retrieve pointers
! Simply return if any pointer not found
ier=0
call gsi_bundlegetpointer(sval,'q',sq,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'q',rq,istatus);ier=istatus+ier
if(ier/=0)return
call gsi_bundlegetpointer(gsi_metguess_bundle(itbin),'q',ges_q_it,ier)
if(ier/=0)return
! Loop over interior of subdomain
if(nstep > 0)then
do k = 1,nsig
do j = 2,lon1+1
do i = 2,lat1+1
ii=istart(mm1)+i-2
! Values for q using stepsizes
q = ges_q_it(i,j,k) + sq(i,j,k)
do kk=1,nstep
qx = q + sges(kk)*rq(i,j,k)
if(qx < zero)then
outmin(kk)=outmin(kk)+(factqmin*wgtfactlats(ii))*qx*qx &
/(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin))
else
if(qx > ges_qsat(i,j,k,itbin))then
outmax(kk)=outmax(kk)+(factqmax*wgtfactlats(ii))*(qx-ges_qsat(i,j,k,itbin))* &
(qx-ges_qsat(i,j,k,itbin))/(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin))
end if
end if
end do
end do
end do
end do
else
do k = 1,nsig
do j = 2,lon1+1
do i = 2,lat1+1
ii=istart(mm1)+i-2
! Values for q using stepsizes
q = ges_q_it(i,j,k)
if(q < zero)then
outmin(1)=outmin(1)+(factqmin*wgtfactlats(ii))*q*q/(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin))
else
if(q > ges_qsat(i,j,k,itbin))then
outmax(1)=outmax(1)+(factqmax*wgtfactlats(ii))*(q-ges_qsat(i,j,k,itbin))*&
(q-ges_qsat(i,j,k,itbin))/(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin))
end if
end if
end do
end do
end do
end if
do kk=2,nstep
outmin(kk)=outmin(kk)-outmin(1)
outmax(kk)=outmax(kk)-outmax(1)
end do
return
end subroutine stplimq
subroutine stplimqc(rval,sval,sges,out,nstep,itbin,cldtype)
!$$$ subprogram documentation block
! . . . .
! subprogram: stplimqc calculate penalty and stepsize for limit of qc
! prgmmr: eliu org: np23 date: 2018-05-30
!
! abstract: calculate stepsize contribution and penalty for limiting q
!
! program history log:
! 2018-05-30 eliu - based on stplimq
!
! input argument list:
! rqc - search direction
! sqc - increment in grid space
! sges - step size estimates (4)
! nstep - number of step size estimates if == 0 then just do outer loop
! itbin - observation bin number (time level)
!
! output argument list:
! outmin(1:nstep) - current penalty for negative q sges(1:nstep)
! outmax(1:nstep) - current penalty for excess q sges(1:nstep)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use mpimod, only: mype
use gridmod, only: lat1,lon1,nsig
use jfunc, only: factql,factqi,factqr,factqs,factqg
use guess_grids, only: ges_qsat
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: nstep,itbin
real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges
real(r_quad),dimension(max(1,nstep)),intent( out) :: out
type(gsi_bundle) ,intent(in ) :: rval,sval
character(2) ,intent(in ) :: cldtype
! Declare local variables
integer(i_kind) i,j,k,kk,ier,ier1,istatus
real(r_kind) qc,qx
real(r_kind) factqc
real(r_kind),pointer,dimension(:,:,:) :: rqc,sqc
real(r_kind),pointer,dimension(:,:,:) :: ges_qc_it=>NULL()
out=zero_quad
! Retrieve pointers
! Simply return if any pointer not found
ier=0; ier1=0; istatus=0
if (cldtype == 'ql') then
factqc =factql
call gsi_bundlegetpointer(sval,'ql',sqc,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'ql',rqc,istatus);ier=istatus+ier
call gsi_bundlegetpointer(gsi_metguess_bundle(itbin),'ql',ges_qc_it,ier1)
endif
if (cldtype == 'qi') then
factqc =factqi
call gsi_bundlegetpointer(sval,'qi',sqc,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'qi',rqc,istatus);ier=istatus+ier
call gsi_bundlegetpointer(gsi_metguess_bundle(itbin),'qi',ges_qc_it,ier1)
endif
if (cldtype == 'qr') then
factqc =factqr
call gsi_bundlegetpointer(sval,'qr',sqc,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'qr',rqc,istatus);ier=istatus+ier
call gsi_bundlegetpointer(gsi_metguess_bundle(itbin),'qr',ges_qc_it,ier1)
endif
if (cldtype == 'qs') then
factqc =factqs
call gsi_bundlegetpointer(sval,'qs',sqc,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'qs',rqc,istatus);ier=istatus+ier
call gsi_bundlegetpointer(gsi_metguess_bundle(itbin),'qs',ges_qc_it,ier1)
endif
if (cldtype == 'qg') then
factqc =factqg
call gsi_bundlegetpointer(sval,'qg',sqc,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'qg',rqc,istatus);ier=istatus+ier
call gsi_bundlegetpointer(gsi_metguess_bundle(itbin),'qg',ges_qc_it,ier1)
endif
if (mype==0) write(6,*)'stplimqc: factqc = ', factqc
if (mype==0) write(6,*)'stplimqc: ier ier1 = ', ier, ier1
if ( factqc==0 ) return
if ( ier/=0 .or. ier1/=0 ) return
! Loop over interior of subdomain
if(nstep > 0)then
do k = 1,nsig
do j = 2,lon1+1
do i = 2,lat1+1
! Values for q using stepsizes
qc = ges_qc_it(i,j,k) + sqc(i,j,k)
do kk=1,nstep
qx = qc + sges(kk)*rqc(i,j,k)
if(qx < zero)then
out(kk)=out(kk)+factqc*qx*qx/(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin))
end if
end do
end do
end do
end do
else
do k = 1,nsig
do j = 2,lon1+1
do i = 2,lat1+1
! Values for q using stepsizes
qc = ges_qc_it(i,j,k)
if(qc < zero)then
out(1)=out(1)+factqc*qc*qc/(ges_qsat(i,j,k,itbin)*ges_qsat(i,j,k,itbin))
end if
end do
end do
end do
end if
do kk=2,nstep
out(kk)=out(kk)-out(1)
end do
return
end subroutine stplimqc
subroutine stplimg(rval,sval,sges,out,nstep)
!$$$ subprogram documentation block
! . . . .
! subprogram: stplimg calculate penalty and stepsize for limit of q
! prgmmr: derber org: np23 date: 1996-11-19
!
! abstract: calculate stepsize contribution and penalty for limiting q
!
! program history log:
! 2011-02-23 zhu
!
! input argument list:
! rg - search direction
! sg - increment in grid space
! sges - step size estimates (4)
! nstep - number of step size estimates if == 0 then just do outer loop
!
! output argument list:
! out(1:nstep) - current penalty for negative gust sges(1:nstep)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use gridmod, only: lat1,lon1
use derivsmod, only: ggues
use jfunc, only: factg
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: nstep
real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges
real(r_quad),dimension(max(1,nstep)),intent( out) :: out
type(gsi_bundle) ,intent(in ) :: rval,sval
! Declare local variables
integer(i_kind) i,j,kk,ier,istatus
real(r_kind) gust,gx
real(r_kind),pointer,dimension(:,:) :: rg,sg
out=zero_quad
if (factg==zero) return
! Retrieve pointers
! Simply return if any pointer not found
ier=0
call gsi_bundlegetpointer(sval,'gust',sg,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'gust',rg,istatus);ier=istatus+ier
if(ier/=0)return
! Loop over interior of subdomain
if(nstep > 0)then
do j = 2,lon1+1
do i = 2,lat1+1
! Values for gust using stepsizes
gust = ggues(i,j) + sg(i,j)
do kk=1,nstep
gx = gust + sges(kk)*rg(i,j)
if(gx < zero)then
out(kk)=out(kk)+factg*gx*gx/(ggues(i,j)*ggues(i,j))
end if
end do
end do
end do
end if
do kk=2,nstep
out(kk)=out(kk)-out(1)
end do
return
end subroutine stplimg
subroutine stplimp(rval,sval,sges,out,nstep)
!$$$ subprogram documentation block
! . . . .
! subprogram: stplimp calculate penalty and stepsize for limit of q
! prgmmr: derber org: np23 date: 1996-11-19
!
! abstract: calculate stepsize contribution and penalty for limiting q
!
! program history log:
! 2011-02-23 zhu
!
! input argument list:
! rp - search direction
! sp - increment in grid space
! sges - step size estimates (4)
! nstep - number of step size estimates if == 0 then just do outer loop
!
! output argument list:
! out(1:nstep) - current penalty for negative pblh sges(1:nstep)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use gridmod, only: lat1,lon1
use derivsmod, only: pgues
use jfunc, only: factp
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: nstep
real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges
real(r_quad),dimension(max(1,nstep)),intent( out) :: out
type(gsi_bundle) ,intent(in ) :: rval,sval
! Declare local variables
integer(i_kind) i,j,kk,ier,istatus
real(r_kind) pblh,px
real(r_kind),pointer,dimension(:,:) :: rp,sp
out=zero_quad
if (factp==zero) return
! Retrieve pointers
! Simply return if any pointer not found
ier=0
call gsi_bundlegetpointer(sval,'pblh',sp,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'pblh',rp,istatus);ier=istatus+ier
if(ier/=0)return
! Loop over interior of subdomain
if(nstep > 0)then
do j = 2,lon1+1
do i = 2,lat1+1
! Values for pblh using stepsizes
pblh = pgues(i,j) + sp(i,j)
do kk=1,nstep
px = pblh + sges(kk)*rp(i,j)
if(px < zero)then
out(kk)=out(kk)+factp*px*px/(pgues(i,j)*pgues(i,j))
end if
end do
end do
end do
end if
do kk=2,nstep
out(kk)=out(kk)-out(1)
end do
return
end subroutine stplimp
subroutine stplimv(rval,sval,sges,out,nstep)
!$$$ subprogram documentation block
! . . . .
! subprogram: stplimv calculate penalty and stepsize for limit of q
! prgmmr: derber org: np23 date: 1996-11-19
!
! abstract: calculate stepsize contribution and penalty for limiting q
!
! program history log:
! 2011-02-23 zhu
!
! input argument list:
! rg - search direction
! sg - increment in grid space
! sges - step size estimates (4)
! nstep - number of step size estimates if == 0 then just do outer loop
!
! output argument list:
! out(1:nstep) - current penalty for negative vis sges(1:nstep)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use gridmod, only: lat1,lon1
use jfunc, only: factv
use derivsmod, only: vgues
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: nstep
real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges
real(r_quad),dimension(max(1,nstep)),intent( out) :: out
type(gsi_bundle) ,intent(in ) :: rval,sval
! Declare local variables
integer(i_kind) i,j,kk,ier,istatus
real(r_kind) vis,vx
real(r_kind),pointer,dimension(:,:) :: rg,sg
out=zero_quad
if (factv==zero) return
! Retrieve pointers
! Simply return if any pointer not found
ier=0
call gsi_bundlegetpointer(sval,'vis',sg,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'vis',rg,istatus);ier=istatus+ier
if(ier/=0)return
! Loop over interior of subdomain
if(nstep > 0)then
do j = 2,lon1+1
do i = 2,lat1+1
! Values for vis using stepsizes
vis = vgues(i,j) + sg(i,j)
do kk=1,nstep
vx = vis + sges(kk)*rg(i,j)
if(vx < zero)then
out(kk)=out(kk)+factv*vx*vx/(vgues(i,j)*vgues(i,j))
end if
end do
end do
end do
end if
do kk=2,nstep
out(kk)=out(kk)-out(1)
end do
return
end subroutine stplimv
subroutine stplimw10m(rval,sval,sges,out,nstep)
!$$$ subprogram documentation block
! . . . .
! subprogram: stplimw10m calculate penalty and stepsize for limit of q
! prgmmr: derber org: np23 date: 1996-11-19
!
! abstract: calculate stepsize contribution and penalty for limiting wspd10m
!
! program history log:
! 2014-03-19 pondeca
!
! input argument list:
! rg - search direction
! sg - increment in grid space
! sges - step size estimates (4)
! nstep - number of step size estimates if == 0 then just do outer loop
!
! output argument list:
! out(1:nstep) - current penalty for negative wspd10m sges(1:nstep)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use gridmod, only: lat1,lon1
use jfunc, only: factw10m
use derivsmod, only: w10mgues
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: nstep
real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges
real(r_quad),dimension(max(1,nstep)),intent( out) :: out
type(gsi_bundle) ,intent(in ) :: rval,sval
! Declare local variables
integer(i_kind) i,j,kk,ier,istatus
real(r_kind) wspd10m,gx
real(r_kind),pointer,dimension(:,:) :: rg,sg
out=zero_quad
if (factw10m==zero) return
! Retrieve pointers
! Simply return if any pointer not found
ier=0
call gsi_bundlegetpointer(sval,'wspd10m',sg,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'wspd10m',rg,istatus);ier=istatus+ier
if(ier/=0)return
! Loop over interior of subdomain
if(nstep > 0)then
do j = 2,lon1+1
do i = 2,lat1+1
! Values for wspd10m using stepsizes
wspd10m = w10mgues(i,j) + sg(i,j)
do kk=1,nstep
gx = wspd10m + sges(kk)*rg(i,j)
if(gx < zero)then
out(kk)=out(kk)+factw10m*gx*gx/(w10mgues(i,j)*w10mgues(i,j))
end if
end do
end do
end do
end if
do kk=2,nstep
out(kk)=out(kk)-out(1)
end do
return
end subroutine stplimw10m
subroutine stplimhowv(rval,sval,sges,out,nstep)
!$$$ subprogram documentation block
! . . . .
! subprogram: stplimhowv calculate penalty and stepsize for limit of howv
! prgmmr: pondeca org: np23 date: 2014-05-07
!
! abstract: calculate stepsize contribution and penalty for limiting howv
!
! program history log:
! 2014-05-07 pondeca
!
! input argument list:
! rg - search direction
! sg - increment in grid space
! sges - step size estimates (4)
! nstep - number of step size estimates if == 0 then just do outer loop
!
! output argument list:
! out(1:nstep) - current penalty for negative howv sges(1:nstep)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use gridmod, only: lat1,lon1
use jfunc, only: facthowv
use derivsmod, only: howvgues
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: nstep
real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges
real(r_quad),dimension(max(1,nstep)),intent( out) :: out
type(gsi_bundle) ,intent(in ) :: rval,sval
! Declare local variables
integer(i_kind) i,j,kk,ier,istatus
real(r_kind) howv,gx
real(r_kind),pointer,dimension(:,:) :: rg,sg
out=zero_quad
if (facthowv==zero) return
! Retrieve pointers
! Simply return if any pointer not found
ier=0
call gsi_bundlegetpointer(sval,'howv',sg,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'howv',rg,istatus);ier=istatus+ier
if(ier/=0)return
! Loop over interior of subdomain
if(nstep > 0)then
do j = 2,lon1+1
do i = 2,lat1+1
! Values for howv using stepsizes
howv = howvgues(i,j) + sg(i,j)
do kk=1,nstep
gx = howv + sges(kk)*rg(i,j)
if(gx < zero)then
out(kk)=out(kk)+facthowv*gx*gx/(howvgues(i,j)*howvgues(i,j))
end if
end do
end do
end do
end if
do kk=2,nstep
out(kk)=out(kk)-out(1)
end do
return
end subroutine stplimhowv
subroutine stplimcldch(rval,sval,sges,out,nstep)
!$$$ subprogram documentation block
! . . . .
! subprogram: stplimcldch calculate penalty and stepsize for limit of cldch
! prgmmr: derber org: np23 date: 1996-11-19
!
! abstract: calculate stepsize contribution and penalty for limiting cldch
!
! program history log:
! 2015-07-10 pondeca
!
! input argument list:
! rg - search direction
! sg - increment in grid space
! sges - step size estimates (4)
! nstep - number of step size estimates if == 0 then just do outer loop
!
! output argument list:
! out(1:nstep) - current penalty for negative cldch sges(1:nstep)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use gridmod, only: lat1,lon1
use jfunc, only: factcldch
use derivsmod, only: cldchgues
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: nstep
real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges
real(r_quad),dimension(max(1,nstep)),intent( out) :: out
type(gsi_bundle) ,intent(in ) :: rval,sval
! Declare local variables
integer(i_kind) i,j,kk,ier,istatus
real(r_kind) cldch,cx
real(r_kind),pointer,dimension(:,:) :: rg,sg
out=zero_quad
if (factcldch==zero) return
! Retrieve pointers
! Simply return if any pointer not found
ier=0
call gsi_bundlegetpointer(sval,'cldch',sg,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'cldch',rg,istatus);ier=istatus+ier
if(ier/=0)return
! Loop over interior of subdomain
if(nstep > 0)then
do j = 2,lon1+1
do i = 2,lat1+1
! Values for cldch using stepsizes
cldch = cldchgues(i,j) + sg(i,j)
do kk=1,nstep
cx = cldch + sges(kk)*rg(i,j)
if(cx < zero)then
out(kk)=out(kk)+factcldch*cx*cx/(cldchgues(i,j)*cldchgues(i,j))
end if
end do
end do
end do
end if
do kk=2,nstep
out(kk)=out(kk)-out(1)
end do
return
end subroutine stplimcldch
subroutine stpliml(rval,sval,sges,out,nstep)
!$$$ subprogram documentation block
! . . . .
! subprogram: stpliml calculate penalty and stepsize for limit of q
! prgmmr: derber org: np23 date: 1996-11-19
!
! abstract: calculate stepsize contribution and penalty for limiting q
!
! program history log:
! 2012-04-23 zhu
!
! input argument list:
! rg - search direction
! sg - increment in grid space
! sges - step size estimates (4)
! nstep - number of step size estimates if == 0 then just do outer loop
!
! output argument list:
! out(1:nstep) - current penalty for negative lcbas sges(1:nstep)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use gridmod, only: lat1,lon1
use jfunc, only: factl
use derivsmod, only: lgues
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: nstep
real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges
real(r_quad),dimension(max(1,nstep)),intent( out) :: out
type(gsi_bundle) ,intent(in ) :: rval,sval
! Declare local variables
integer(i_kind) i,j,kk,ier,istatus
real(r_kind) lcbas,vx
real(r_kind),pointer,dimension(:,:) :: rg,sg
out=zero_quad
if (factl==zero) return
! Retrieve pointers
! Simply return if any pointer not found
ier=0
call gsi_bundlegetpointer(sval,'lcbas',sg,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'lcbas',rg,istatus);ier=istatus+ier
if(ier/=0)return
! Loop over interior of subdomain
if(nstep > 0)then
do j = 2,lon1+1
do i = 2,lat1+1
! Values for lcbas using stepsizes
lcbas = lgues(i,j) + sg(i,j)
do kk=1,nstep
vx = lcbas + sges(kk)*rg(i,j)
if(vx < zero)then
out(kk)=out(kk)+factl*vx*vx/(lgues(i,j)*lgues(i,j))
end if
end do
end do
end do
end if
do kk=2,nstep
out(kk)=out(kk)-out(1)
end do
return
end subroutine stpliml
subroutine stpjcpdry(rval,sval,pen,b,c,nbins)
!$$$ subprogram documentation block
! . . . .
! subprogram: stpjcpdry penalty and stp size for mean dry ps conservation
! prgmmr: kleist org: np23 date: 2009-07-07
!
! abstract: calculate stepsize contribution and penalty for limiting changes to global
! mean dry ps increment
!
! program history log:
! 2009-07-07 kleist
! 2010-05-13 todling - update to use gsi_bundle
! 2010-05-25 derber - modify to decrease number of communications
! 2010-08-18 hu - add qpvals= to mpl_allreduce call
! 2011-11-01 eliu - add handling for ql & qi increments and search directions
! 2013-05-05 todling - separate dry mass from the rest (zero-diff change)
!
! input argument list:
! rq - q search direction
! rc - cloud water search direction
! rp - surface pressure search direction
! sq - q increment
! sc - cloud water increment
! sp - increment in grid space
!
! output argument list:
! pen - current penalty for mean dry pressure constraint
! b - contribution to numerator
! c - contribution to denomenator
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use mpimod, only: mype
use gridmod, only: lat2,lon2,nsig,wgtlats,nlon,istart
use guess_grids, only: ges_prsi,ntguessig
use mpl_allreducemod, only: mpl_reduce
use jcmod, only: bamp_jcpdry
use gsi_bundlemod, only: assignment(=)
use gsi_metguess_mod, only: gsi_metguess_get
implicit none
! Declare passed variables
type(gsi_bundle),dimension(nbins),intent(in ) :: sval
type(gsi_bundle),dimension(nbins),intent(in ) :: rval
real(r_quad) ,intent( out) :: pen,b,c
integer(i_kind) ,intent(in ) :: nbins
! Declare local variables
real(r_quad),dimension(2*nbins):: dmass
real(r_quad) :: rcon,con
integer(i_kind) i,j,k,it,mm1,ii,ier,icw,iql,iqi,iqr,iqs,iqg,istatus,n
real(r_kind),pointer,dimension(:,:,:) :: rq,sq,rc,sc,rql,rqi,sql,sqi
real(r_kind),pointer,dimension(:,:,:) :: rqr,rqs,rqg,sqr,sqs,sqg
real(r_kind),pointer,dimension(:,:) :: rp,sp
logical return_now
real(r_quad) :: dmn, dmn2
pen=zero_quad ; b=zero_quad ; c=zero_quad
it=ntguessig
dmass=zero_quad
rcon=one_quad/(two_quad*float(nlon))
mm1=mype+1
return_now = .false.
do n=1,nbins
! Retrieve pointers
! Simply return if any pointer not found
ier=0; icw=0; iql=0; iqi=0; iqr=0; iqs=0; iqg=0; istatus=0
call gsi_bundlegetpointer(sval(n),'q' ,sq, istatus) ; ier=istatus+ier
call gsi_bundlegetpointer(sval(n),'cw',sc, istatus) ; icw=istatus+icw
call gsi_bundlegetpointer(sval(n),'ql',sql,istatus) ; iql=istatus+iql
call gsi_bundlegetpointer(sval(n),'qi',sqi,istatus) ; iqi=istatus+iqi
call gsi_bundlegetpointer(sval(n),'qr',sqr,istatus) ; iqr=istatus+iqr
call gsi_bundlegetpointer(sval(n),'qs',sqs,istatus) ; iqs=istatus+iqs
call gsi_bundlegetpointer(sval(n),'qg',sqg,istatus) ; iqg=istatus+iqg
call gsi_bundlegetpointer(sval(n),'ps',sp, istatus) ; ier=istatus+ier
call gsi_bundlegetpointer(rval(n),'q' ,rq, istatus) ; ier=istatus+ier
call gsi_bundlegetpointer(rval(n),'cw',rc, istatus) ; icw=istatus+icw
call gsi_bundlegetpointer(rval(n),'ql',rql,istatus) ; iql=istatus+iql
call gsi_bundlegetpointer(rval(n),'qi',rqi,istatus) ; iqi=istatus+iqi
call gsi_bundlegetpointer(rval(n),'qr',rqr,istatus) ; iqr=istatus+iqr
call gsi_bundlegetpointer(rval(n),'qs',rqs,istatus) ; iqs=istatus+iqs
call gsi_bundlegetpointer(rval(n),'qg',rqg,istatus) ; iqg=istatus+iqg
call gsi_bundlegetpointer(rval(n),'ps',rp, istatus) ; ier=istatus+ier
if( ier/=0 .or. ((iql+iqi)/=0 .and. icw/=0) ) then
if (mype==0) write(6,*) 'stpjcpdry: warning - missing some required variables'
if (mype==0) write(6,*) 'stpjcpdry: dry mass constraint not performed'
return
end if
! Calculate mean surface pressure contribution in subdomain
dmn = dmass(n)
dmn2 = dmass(n + nbins)
!$omp parallel do private(i,j,ii) firstprivate(rcon,con) reduction(+:dmn) reduction(+:dmn2) collapse(2)
do j=2,lon2-1
do i=2,lat2-1
ii=istart(mm1)+i-2
con=wgtlats(ii)*rcon
dmn=dmn+sp(i,j)*con
dmn2=dmn2+rp(i,j)*con
end do
end do
!$omp end parallel do
dmass(n) = dmn
dmass(n+nbins) = dmn2
! Remove water to get incremental dry ps
!$omp parallel do private(k,j,i,ii) firstprivate(rcon,con) reduction(-:dmn) reduction(-:dmn2) collapse(2)
do k=1,nsig
do j=2,lon2-1
do i=2,lat2-1
ii=istart(mm1)+i-2
con=(ges_prsi(i,j,k,it)-ges_prsi(i,j,k+1,it))*wgtlats(ii)*rcon
dmn=dmn - sq(i,j,k)*con
dmn2=dmn2 - rq(i,j,k)*con
if(icw==0)then
dmn=dmn - sc(i,j,k)*con
dmn2=dmn2 - rc(i,j,k)*con
else
dmn=dmn - (sql(i,j,k)+sqi(i,j,k))*con
dmn2=dmn2 - (rql(i,j,k)+rqi(i,j,k))*con
if (iqr==0) then
dmn = dmn - sqr(i,j,k)*con
dmn2= dmn2- rqr(i,j,k)*con
endif
if (iqs==0) then
dmn = dmn - sqs(i,j,k)*con
dmn2= dmn2- rqs(i,j,k)*con
endif
if (iqg==0) then
dmn = dmn - sqg(i,j,k)*con
dmn2= dmn2- rqg(i,j,k)*con
endif
endif
end do
end do
end do
!$omp end parallel do
dmass(n) = dmn
dmass(n+nbins) = dmn2
end do
call mpl_reduce(2*nbins,0,qpvals=dmass)
! Now penalize non-zero global mean dry ps increment
! Notice there will only be a contribution from PE=0
if(mype == 0)then
do n=1,nbins
pen = pen + bamp_jcpdry*dmass(n)*dmass(n)
b = b - bamp_jcpdry*dmass(n+nbins)*dmass(n)
c = c + bamp_jcpdry*dmass(n+nbins)*dmass(n+nbins)
end do
end if
return
end subroutine stpjcpdry
subroutine stpjcdfi(rval,sval,pen,b,c)
!$$$ subprogram documentation block
! . . . .
! subprogram: stpjcdfi penalty and stp size for weak constraint DFI
! prgmmr: kleist org: np23 date: 2012-01-19
!
! abstract: calculate stepsize contribution and penalty for Jc DFI
!
! program history log:
! 2012-01-19 kleist - adaptation of evaljcdfi
! 2013-05-15 todling - penalty should be non-zero on root only
!
! input argument list:
! rval -
! sval -
!
! output argument list:
! pen - current penalty for Jc DFI
! b - contribution to numerator
! c - contribution to denomenator
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use jcmod, only: wgtdfi,alphajc
use gsi_bundlemod, only : self_add,self_mul,assignment(=)
use gsi_4dvar, only: nobs_bins
use mpimod, only: mype
use state_vectors, only : allocate_state,deallocate_state
use state_vectors, only : dot_product
implicit none
! Declare passed variables
type(gsi_bundle),dimension(nobs_bins),intent(in) :: sval
type(gsi_bundle),dimension(nobs_bins),intent(in) :: rval
real(r_quad),intent(out) :: pen,b,c
! Declare local variables
integer(i_kind) :: jj,idfi
real(r_quad) :: pjc,rjc
type(gsi_bundle) :: sfilter,afilter
type(gsi_bundle) :: rfilter,bfilter
pen=zero_quad ; b=zero_quad ; c=zero_quad
!************************************************************************************
idfi = (nobs_bins-1)/2+1
call allocate_state(sfilter)
call allocate_state(afilter)
call allocate_state(rfilter)
call allocate_state(bfilter)
! Compute filtered state
sfilter=zero
rfilter=zero
do jj=1,nobs_bins
call self_add(sfilter,wgtdfi(jj),sval(jj))
call self_add(rfilter,wgtdfi(jj),rval(jj))
enddo
! Compute difference from filtered state
call self_add(sfilter,-one,sval(idfi))
call self_add(rfilter,-one,rval(idfi))
! Apply Jc multiplicative factor
call self_mul(sfilter,alphajc)
call self_mul(rfilter,alphajc)
! Convert to energy norm/apply grid factors
call enorm_state_red(sfilter,pjc,afilter)
call enorm_state_red(rfilter,rjc,bfilter)
! Penalty, b, c
if (mype==0) then
pen = dot_product(afilter,afilter)
b = -dot_product(afilter,bfilter)
c = dot_product(bfilter,bfilter)
write(6,*)'Jc DFI penalty = ',pen
endif
call deallocate_state(sfilter)
call deallocate_state(afilter)
call deallocate_state(rfilter)
call deallocate_state(bfilter)
return
end subroutine stpjcdfi
end module stpjcmod