subroutine pcp_k(km,dtp,del_in,sl_in,rbs,&
slmask,xkt2,ncloud,frain,rmmhr,&
psexp,dplat,dplon,&
t_in,q_in,u_in,v_in,div_in,cwm_in,&
tsen_ten_in,q_ten_in,p_ten_in,&
tsas_o,qsas_o,rnsas_of,cldwrk,kbcon,ktcon,jmin,kuo,&
tlrg_o,qlrg_o,rnlrg_of,&
t_out,q_out,cwm_out,u_out,v_out,rn_out,&
t_in_ad,q_in_ad,cwm_in_ad,u_in_ad,v_in_ad,div_in_ad,&
t_out_ad,q_out_ad,cwm_out_ad,u_out_ad,v_out_ad,rn_out_ad)
!$$$ subprogram documentation block
! . . . .
! subprogram: pcp_k driver for precipitation forward & adjoint models
! prgmmr: treadon org: np23 date: 2003-12-18
!
! abstract: This subroutine calls GFS precipitation physics routines followed
! by a call to adjoint of these routines. The GFS precipitation
! physics include both convective and grid-scale (explicit) processes.
! Convective precipitation is parameterized using the Simplified
! Arakawa-Schubert scheme originally developed by George Grell.
! Explicit precipitation processes are modeled folliwng the work of
! Q. Zhao.
!
! program history log:
! 2003-12-18 treadon - initial routine
! 2004-06-14 treadon - reformat documenation
! 2006-04-12 treadon - change del and sl from 1d to 2d arrays
! 2006-09-15 treadon - change (k,i) arrays to (k) arrays
! 2006-10-12 treadon - remove virtual temperature
! 2008-04-29 safford - rm unused uses
! 2010-03-31 treadon - replace jcap with sp_a%jcap
! 2013-01-15 parrish - convert gscond_ad.f90,nlmsas_ad.f90,omegas_ad.f90 to modules
! and add interfaces to account for type mismatch
!
! input argument list:
! km - number of levels in vertical profile
! dtp - physics timestep
! del - "sigma" thickness of layers
! sl - "sigma" value at layer midpoints
! rbs - 1/(sin(latitude)**2)
! slmask - sea (=0), land (=1), ice/snow (=2) mask
! xkt2 - random number for cloud top selection in SAS
! ncloud - flag to turn on cloud liquid water detrainment in SAS
! frain - factor to account for difference between physics and model timestep
! rmmhr - conversion factor to get rain rate units of mm/hr
! psexp - surface pressure (cb)
! dplat - partial derivative of ln(ps) with respect to latitude
! dplon - partial derivative of ln(ps) with respect to longitude
! del_in -
! sl_in -
! t_in - sensible temperature
! q_in - specific humidity
! u_in - zonal wind component
! v_in - meridional wind component
! div_in - divergence
! cwm_in - cloud condensate mixing ratio
! t_out_ad - temperature perturbation
! q_out_ad - q perturbation
! cwm_out_ad- cloud condensate mixing ratio perturbation
! u_out_ad - zonal wind perturbation
! v_out_ad - meridional wind perturbation
! rn_out_ad - rain rate perturbation
! tsen_ten_in
! q_ten_in
! p_ten_in
!
! output argument list:
! tsas_o - temperature following call to SAS
! qsas_o - q following call to SAS
! rnsas_of - convective precipitation rate
! cldwrk - cloud work function computed/used by SAS
! kbcon - integer index of model level at which SAS updraft originates
! ktcon - integer index of model level for SAS cloud top
! jmin - integer index of model level at which SAS downdraft originates
! kuo - integer flag for convective activity (0=no convection, 1=convection)
! tlrg_o - temperature following call to explicit precipitation routines
! qlrg_o - q following call to explicit precipitation routines
! rnlrg_of - precipitation rate following call to explicit precipitation routines
! t_out - temperature following call to all GFS precipitation physics
! q_out - q following call to all GFS precipitation physics
! cwm_out - cloud condensate mixing ratio following call to all GFS precipitation physics
! u_out - zonal wind component following call to all GFS precipitation physics
! v_out - meridional wind component following call to all GFS precipitation physics
! rn_out - precipitation rate following call to all GFS precipitation physics
! t_in_ad - partial derivative of temperature with respect to rain rate
! q_in_ad - partial derivative of q with respect to rain rate
! cwm_in_ad - partial derivative of cloud condensate mixing ratio with respect to rain rate
! u_in_ad - partial derivative of zonal wind with respect to rain rate
! v_in_ad - partial derivative of meridional with respect to rain rate
! div_in_ad - partial derivative of divergence with respect to rain rate
!
! remarks:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
!
use kinds, only: r_kind,i_kind
use constants, only: rhcbot,rhctop,dx_inv,dx_min,one,zero
use pcpinfo, only: tiny_obs
use gridmod, only: nlon,sp_a
use gscond_ad_mod, only: gscond_ad
use nlmsas_ad_mod, only: nlmsas_ad
use omegas_ad_mod, only: omegas_ad
implicit none
! Declare passed variables
integer(i_kind) , intent(in ) :: km,ncloud
integer(i_kind) , intent( out) :: kbcon,jmin,ktcon,kuo
real(r_kind) , intent(in ) :: dtp,frain,rmmhr
real(r_kind) , intent(in ) :: rbs,dplat,dplon,slmask,psexp
real(r_kind) , intent( out) :: cldwrk,rn_out
real(r_kind),dimension(km), intent(in ) :: del_in,sl_in
real(r_kind),dimension(km), intent(in ) :: t_in,q_in,cwm_in,u_in,v_in,div_in,&
t_out_ad,q_out_ad,cwm_out_ad,u_out_ad,v_out_ad,&
tsen_ten_in,q_ten_in,p_ten_in
real(r_kind),dimension(km), intent( out) :: t_out,q_out,cwm_out,u_out,v_out,&
t_in_ad,q_in_ad,u_in_ad,v_in_ad,div_in_ad,cwm_in_ad
! Declare local parameters
real(r_kind),parameter:: r0_99=0.99_r_kind
! Declare local arrays
logical:: skipsas,skiplrg
logical adjoint
integer(i_kind):: k,kb,im,ix
real(r_kind):: ps,rcl,rcs,xkt,xkt2,rnlrg_o,&
rnlrg_o_ad,rnsas_o,rnsas_o_ad,rnsas_of,rnlrg_of,rnsas_of_ad,&
rnlrg_of_ad,rn_out_ad,work2,tem,work1
real(r_kind),dimension(km):: rhc,u_i,v_i,div_i,vvel_o,&
u_i_ad,v_i_ad,div_i_ad,vvel_o_ad,&
t_ten_i,q_ten_i,p_ten_i,qgs_i,cwmgs_i,tgs_i,&
qgs_o,cwmgs_o,tgs_o, qlrg_i, cwmlrg_i, tlrg_i,&
qlrg_o, cwmlrg_o, tlrg_o,t_ten_i_ad,q_ten_i_ad,p_ten_i_ad,&
qgs_i_ad,cwmgs_i_ad,tgs_i_ad,qgs_o_ad,cwmgs_o_ad,tgs_o_ad,&
qlrg_i_ad, cwmlrg_i_ad, tlrg_i_ad,&
qlrg_o_ad, cwmlrg_o_ad, tlrg_o_ad,&
tsas_i,qsas_i,cwmsas_i,usas_i,vsas_i,wsas_i,&
tsas_o,qsas_o,cwmsas_o,usas_o,vsas_o,&
tsas_i_ad,qsas_i_ad,cwmsas_i_ad,usas_i_ad,vsas_i_ad,wsas_i_ad,&
tsas_o_ad,qsas_o_ad,cwmsas_o_ad,usas_o_ad,vsas_o_ad,&
save_tlrg,save_qlrg,save_cwmlrg,&
save_tsas,save_qsas,save_cwmsas,save_usas,save_vsas,&
save_wsas,del_i,sl_i
!**************************************************************************
! Initialize output arrays to zero
im=1
ix=1
cldwrk = zero
kb = 0
jmin = 0
kbcon = 0
ktcon = 0
kuo = 0
rn_out = zero
do k = 1,km
t_out(k) = zero
q_out(k) = zero
cwm_out(k) = zero
u_out(k) = zero
v_out(k) = zero
t_in_ad(k) = zero
q_in_ad(k) = zero
u_in_ad(k) = zero
v_in_ad(k) = zero
div_in_ad(k)= zero
cwm_in_ad(k)= zero
end do
! SASCNV needs vertical velocity, so compute it.
rcl = rbs
ps = psexp
do k = 1,km
u_i(k) = u_in(k)
v_i(k) = v_in(k)
div_i(k) = div_in(k)
del_i(k) = del_in(k)
sl_i(k) = sl_in(k)
vvel_o(k) = zero
u_i_ad(k) = zero
v_i_ad(k) = zero
div_i_ad(k) = zero
vvel_o_ad(k)= zero
end do
adjoint = .false.
call omegas_ad(im,ix,km,dplat,dplon,u_i,v_i,div_i,ps,rcl,&
del_i,sl_i,vvel_o,u_i_ad,v_i_ad,div_i_ad,vvel_o_ad,adjoint)
! Call convective parameterization, SASCNV.
rnsas_o = zero
rnlrg_o = zero
rcs = sqrt(rcl)
xkt = xkt2
rnsas_o_ad = zero
do k = 1,km
usas_i(k) = u_in(k)
vsas_i(k) = v_in(k)
wsas_i(k) = vvel_o(k)
tsas_i(k) = t_in(k)
qsas_i(k) = q_in(k)
cwmsas_i(k) = cwm_in(k)
del_i(k) = del_in(k)
sl_i(k) = sl_in(k)
usas_o(k) = zero
vsas_o(k) = zero
tsas_o(k) = zero
qsas_o(k) = zero
cwmsas_o(k) = zero
usas_o_ad(k) = zero
vsas_o_ad(k) = zero
tsas_o_ad(k) = zero
qsas_o_ad(k) = zero
cwmsas_o_ad(k) = zero
tsas_i_ad(k) = zero
qsas_i_ad(k) = zero
cwmsas_i_ad(k) = zero
usas_i_ad(k) = zero
vsas_i_ad(k) = zero
wsas_i_ad(k) = zero
end do
adjoint = .false.
call nlmsas_ad(im,ix,km,sp_a%jcap,dtp,del_i,sl_i,rcs,&
slmask,xkt,ncloud,psexp,&
tsas_i,qsas_i,cwmsas_i,usas_i,vsas_i,wsas_i,&
tsas_o,qsas_o,cwmsas_o,usas_o,vsas_o,rnsas_o,&
cldwrk,kbcon,ktcon,jmin,kuo,kb,&
tsas_i_ad,qsas_i_ad,cwmsas_i_ad,usas_i_ad,vsas_i_ad,wsas_i_ad,&
tsas_o_ad,qsas_o_ad,cwmsas_o_ad,usas_o_ad,vsas_o_ad,rnsas_o_ad,&
adjoint)
! Transfer u and v to output arrays
do k=1,km
u_out(k) = usas_o(k)
v_out(k) = vsas_o(k)
end do
! Compute critical threshold relative humidities
tem = (rhctop-rhcbot)/(km-one)
work1 = (log(one/(rcs*nlon))-dx_min) * dx_inv
work2 = one - work1
do k=1,km
rhc(k) = rhcbot + tem*(k-1)
rhc(k) = r0_99*work1 + rhc(k)*work2
end do
! Call gridscale condensation routine, GSCOND
do k = 1,km
t_ten_i(k) = tsen_ten_in(k)
q_ten_i(k) = q_ten_in(k)
p_ten_i(k) = p_ten_in(k)
sl_i(k) = sl_in(k)
tgs_i(k) = tsas_o(k)
qgs_i(k) = qsas_o(k)
cwmgs_i(k) = cwmsas_o(k)
tgs_o(k) = zero
qgs_o(k) = zero
cwmgs_o(k) = zero
t_ten_i_ad(k) = zero
q_ten_i_ad(k) = zero
p_ten_i_ad(k) = zero
tgs_i_ad(k) = zero
qgs_i_ad(k) = zero
cwmgs_i_ad(k) = zero
tgs_o_ad(k) = zero
qgs_o_ad(k) = zero
cwmgs_o_ad(k) = zero
end do
adjoint = .false.
call gscond_ad(im,ix,km,dtp,sl_i,psexp,rhc,&
t_ten_i,q_ten_i,p_ten_i,&
qgs_i,cwmgs_i,tgs_i,&
qgs_o,cwmgs_o,tgs_o,&
t_ten_i_ad,q_ten_i_ad,p_ten_i_ad,&
qgs_i_ad,cwmgs_i_ad,tgs_i_ad,&
qgs_o_ad,cwmgs_o_ad,tgs_o_ad,&
adjoint)
! Call gridscale precipitation routine, PRECPD
rnlrg_o = zero
rnlrg_o_ad = zero
do k = 1,km
tlrg_i(k) = tgs_o(k)
qlrg_i(k) = qgs_o(k)
cwmlrg_i(k) = cwmgs_o(k)
del_i(k) = del_in(k)
sl_i(k) = sl_in(k)
tlrg_o(k) = zero
qlrg_o(k) = zero
cwmlrg_o(k) = zero
tlrg_i_ad(k) = zero
qlrg_i_ad(k) = zero
cwmlrg_i_ad(k) = zero
tlrg_o_ad(k) = zero
qlrg_o_ad(k) = zero
cwmlrg_o_ad(k) = zero
end do
adjoint = .false.
call precpd_ad(km,dtp,del_i,sl_i,psexp,rhc,&
qlrg_i, cwmlrg_i, tlrg_i,&
qlrg_o, cwmlrg_o, tlrg_o,rnlrg_o,&
qlrg_i_ad, cwmlrg_i_ad, tlrg_i_ad,&
qlrg_o_ad, cwmlrg_o_ad, tlrg_o_ad,rnlrg_o_ad,&
adjoint)
! Combine convective and gridscale precipitation to get the total
rnsas_of = frain*rnsas_o*rmmhr
rnlrg_of = frain*rnlrg_o*rmmhr
rn_out = rnsas_of + rnlrg_of
skipsas = .false.
skiplrg = .false.
! If convective or gridscale precipitation are too small (ie, < tiny_obs),
! then reset precipitation contribution to zero
if (rnsas_of<tiny_obs) then
skipsas = .true.
rnsas_of = zero
endif
if (rnlrg_of<tiny_obs) then
skiplrg = .true.
rnlrg_of = zero
endif
rn_out = rnsas_of + rnlrg_of
! Load output arrays with adjusted temperature and
! moisture profiles. Convert dry temperature to
! virtual temperature.
do k = 1,km
cwm_out(k) = cwmlrg_o(k)
q_out(k) = qlrg_o(k)
t_out(k) = tlrg_o(k)
end do
!==================================================================
! BEGIN ADJOINT
!
! Adjoint of combine convective and gridscale precipitation
! to get the total.
rnsas_of_ad = rn_out_ad
rnlrg_of_ad = rn_out_ad
if (skipsas) rnsas_of_ad = zero
if (skiplrg) rnlrg_of_ad = zero
rnsas_o_ad = frain*rnsas_of_ad*rmmhr
rnlrg_o_ad = frain*rnlrg_of_ad*rmmhr
! Call adjoint of gridscale precipitation routine, PRECPD
rnlrg_o = zero
do k = 1,km
tlrg_i(k) = tgs_o(k)
qlrg_i(k) = qgs_o(k)
cwmlrg_i(k) = cwmgs_o(k)
del_i(k) = del_in(k)
sl_i(k) = sl_in(k)
tlrg_o(k) = zero
qlrg_o(k) = zero
cwmlrg_o(k) = zero
tlrg_i_ad(k) = zero
qlrg_i_ad(k) = zero
cwmlrg_i_ad(k) = zero
tlrg_o_ad(k) = t_out_ad(k)
qlrg_o_ad(k) = q_out_ad(k)
cwmlrg_o_ad(k) = cwm_out_ad(k)
save_tlrg(k) = tlrg_o_ad(k)
save_qlrg(k) = qlrg_o_ad(k)
save_cwmlrg(k) = cwmlrg_o_ad(k)
end do
adjoint = .true.
call precpd_ad(km,dtp,del_i,sl_i,psexp,rhc,&
qlrg_i, cwmlrg_i, tlrg_i,&
qlrg_o, cwmlrg_o, tlrg_o,rnlrg_o,&
qlrg_i_ad, cwmlrg_i_ad, tlrg_i_ad,&
qlrg_o_ad, cwmlrg_o_ad, tlrg_o_ad,rnlrg_o_ad,&
adjoint)
! If gridscale precipitation < tiny_obs, do not include forcing from
! gridscale precipitation in adjoint (ie, gradient).
if (skiplrg) then
do k=1,km
tlrg_i_ad(k) = save_tlrg(k)
qlrg_i_ad(k) = save_qlrg(k)
cwmlrg_i_ad(k) = save_cwmlrg(k)
end do
endif
! Adjoint of call to GSCOND
do k = 1,km
t_ten_i(k) = tsen_ten_in(k)
q_ten_i(k) = q_ten_in(k)
p_ten_i(k) = p_ten_in(k)
sl_i(k) = sl_in(k)
tgs_i(k) = tsas_o(k)
qgs_i(k) = qsas_o(k)
cwmgs_i(k) = cwmsas_o(k)
tgs_o(k) = zero
qgs_o(k) = zero
cwmgs_o(k) = zero
t_ten_i_ad(k) = zero
q_ten_i_ad(k) = zero
p_ten_i_ad(k) = zero
tgs_i_ad(k) = zero
qgs_i_ad(k) = zero
cwmgs_i_ad(k) = zero
tgs_o_ad(k) = tlrg_i_ad(k)
qgs_o_ad(k) = qlrg_i_ad(k)
cwmgs_o_ad(k) = cwmlrg_i_ad(k)
end do
adjoint = .true.
call gscond_ad(im,ix,km,dtp,sl_i,psexp,rhc,&
t_ten_i,q_ten_i,p_ten_i,&
qgs_i,cwmgs_i,tgs_i,&
qgs_o,cwmgs_o,tgs_o,&
t_ten_i_ad,q_ten_i_ad,p_ten_i_ad,&
qgs_i_ad,cwmgs_i_ad,tgs_i_ad,&
qgs_o_ad,cwmgs_o_ad,tgs_o_ad,&
adjoint)
! Adjoint of u and v transfers
do k=1,km
usas_o_ad(k) = u_out_ad(k)
vsas_o_ad(k) = v_out_ad(k)
end do
! Adjoint of call to convective parameterization, SASCNV
rnsas_o = zero
xkt = xkt2
do k = 1,km
usas_i(k) = u_in(k)
vsas_i(k) = v_in(k)
wsas_i(k) = vvel_o(k)
tsas_i(k) = t_in(k)
qsas_i(k) = q_in(k)
cwmsas_i(k) = cwm_in(k)
del_i(k) = del_in(k)
sl_i(k) = sl_in(k)
usas_o(k) = zero
vsas_o(k) = zero
tsas_o(k) = zero
qsas_o(k) = zero
cwmsas_o(k) = zero
tsas_o_ad(k) = tgs_i_ad(k)
qsas_o_ad(k) = qgs_i_ad(k)
cwmsas_o_ad(k) = cwmgs_i_ad(k)
tsas_i_ad(k) = zero
qsas_i_ad(k) = zero
cwmsas_i_ad(k) = zero
usas_i_ad(k) = zero
vsas_i_ad(k) = zero
wsas_i_ad(k) = zero
save_usas(k) = usas_o_ad(k)
save_vsas(k) = vsas_o_ad(k)
save_wsas(k) = zero
save_qsas(k) = qgs_i_ad(k)
save_tsas(k) = tgs_i_ad(k)
save_cwmsas(k) = cwmgs_i_ad(k)
end do
adjoint = .true.
call nlmsas_ad(im,ix,km,sp_a%jcap,dtp,del_i,sl_i,rcs,&
slmask,xkt,ncloud,psexp,&
tsas_i,qsas_i,cwmsas_i,usas_i,vsas_i,wsas_i,&
tsas_o,qsas_o,cwmsas_o,usas_o,vsas_o,rnsas_o,&
cldwrk,kbcon,ktcon,jmin,kuo,kb,&
tsas_i_ad,qsas_i_ad,cwmsas_i_ad,usas_i_ad,vsas_i_ad,wsas_i_ad,&
tsas_o_ad,qsas_o_ad,cwmsas_o_ad,usas_o_ad,vsas_o_ad,rnsas_o_ad,&
adjoint)
! If convective precipitation < tiny_obs, do not include forcing from
! convective precipitation in adjoint (ie, gradient).
if (skipsas) then
do k=1,km
usas_i_ad(k) = save_usas(k)
vsas_i_ad(k) = save_vsas(k)
wsas_i_ad(k) = save_wsas(k)
qsas_i_ad(k) = save_qsas(k)
tsas_i_ad(k) = save_tsas(k)
cwmsas_i_ad(k) = save_cwmsas(k)
end do
endif
! Adjoint of "0" --> "sas0" transfer
do k = 1,km
u_in_ad(k) = usas_i_ad(k)
v_in_ad(k) = vsas_i_ad(k)
vvel_o_ad(k) = wsas_i_ad(k)
q_in_ad(k) = qsas_i_ad(k)
t_in_ad(k) = tsas_i_ad(k)
cwm_in_ad(k) = cwmsas_i_ad(k)
end do
! Adjoint of call to omegas
rcl = rbs
ps = psexp
do k = 1,km
u_i(k) = u_in(k)
v_i(k) = v_in(k)
div_i(k) = div_in(k)
del_i(k) = del_in(k)
sl_i(k) = sl_in(k)
vvel_o(k) = zero
! u_i_ad(k) = u_in_ad(k)
! v_i_ad(k) = v_in_ad(k)
u_i_ad(k) = zero
v_i_ad(k) = zero
div_i_ad(k) = zero
end do
adjoint = .true.
call omegas_ad(im,ix,km,dplat,dplon,u_i,v_i,div_i,ps,rcl,&
del_i,sl_i,vvel_o,u_i_ad,v_i_ad,div_i_ad,vvel_o_ad,adjoint)
! Adjoint of _i --> _in transfer
do k = 1,km
u_in_ad(k) = u_in_ad(k) + u_i_ad(k)
v_in_ad(k) = v_in_ad(k) + v_i_ad(k)
div_in_ad(k) = div_i_ad(k)
end do
! End of routine
return
end subroutine pcp_k