subroutine pbl_ad(u,v,t,ps,jstart,jstop)
!$$$ subprogram documentation block
! . . .
! subprogram: pbl_ad
!
! prgrmmr: m. rancic
!
! abstract: adjoint of pbl_tl
!
! program history log:
! 2008-04-02 safford -- add subprogram doc block, rm unused uses
!
! output argument list:
! ps -
! u -
! v -
! t -
!
! input argument list:
! u -
! v -
! t -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds,only: r_kind,i_kind
use constants,only: one,zero,two,four,five,half,rd_over_g,rd_over_cp,grav
use gridmod,only: lat2,lon2,nsig
use pblmod, only: dudz,dvdz,dodz,zi,rdzi,rdzl,eps_m
use pblmod, only: lmbd,karm,karm0,alph,beta,epxilon,nsig_pbl
use pblmod, only: uges0,vges0,oges0,pges0,tges0,uges1,vges1,oges1
use tends4pertmod, only: time_step_half
implicit none
! Declare passed variables
real(r_kind),dimension(lat2,lon2,nsig),intent(inout):: u,v,t
real(r_kind),dimension(lat2,lon2) ,intent( out):: ps
integer(i_kind) ,intent(in ):: jstart,jstop
! Declare local parameters
real(r_kind),parameter:: r10 = 10.0_r_kind
real(r_kind),parameter:: r20 = 20.0_r_kind
! Declare local variables
real(r_kind),dimension(lat2,lon2,nsig+1):: prs
real(r_kind),dimension(nsig_pbl):: u1_bg,v1_bg,o1_bg
real(r_kind),dimension(nsig_pbl):: u_bg,v_bg,t_bg,o_bg,rdzi_bg,rdzi_ad
real(r_kind),dimension(nsig_pbl):: rdzl_ad,dudz_ad,dvdz_ad,dodz_ad,km_ad
real(r_kind),dimension(nsig_pbl+1):: p_bg,zi_bg
real(r_kind),dimension(2:nsig_pbl):: rdzl_bg
real(r_kind),dimension(2:nsig_pbl):: lmix_bg,ri_bg,km_bg,zmix_bg,ssq_bg
real(r_kind),dimension(2:nsig_pbl):: dudz_bg,dvdz_bg,dodz_bg
real(r_kind),dimension(2:nsig_pbl):: lmix_ad,ri_ad,rho_bg
real(r_kind),dimension(nsig_pbl):: u_ad,v_ad,t_ad,zl_ad,o_ad
real(r_kind),dimension(nsig_pbl+1):: p_ad,zi_ad,ck_ad,ck_bg
real(r_kind),dimension(nsig_pbl+1):: dzi_ad
real(r_kind),dimension(nsig_pbl):: aux_bg,ckdif_bg,cksum_bg
real(r_kind) ckplus_bg,ckmnus_bg,ckplus_ad,ckmnus_ad
real(r_kind) rssq_bg,ssq_ad
real(r_kind) zmix_ad
real(r_kind):: acoef,bcoef,ccoef
real(r_kind):: ax1,ax2
real(r_kind):: aux_ad
real(r_kind) rho_ad
real(r_kind),dimension(nsig_pbl):: acoef0,bcoef0,ccoef0
real(r_kind),dimension(nsig_pbl):: acoef1,bcoef1,ccoef1
real(r_kind),dimension(nsig_pbl):: a_ad,b_ad,c_ad
real(r_kind),dimension(nsig_pbl):: wu,ua,ub,uc,qu,pu
real(r_kind),dimension(nsig_pbl):: wv,va,vb,vc,qv,pv
real(r_kind),dimension(nsig_pbl):: wo,oa,ob,oc,qo,po
real(r_kind),dimension(nsig_pbl):: fu,fv,fo
integer(i_kind) i,j,k
! Preliminaries
ua(1)=zero; va(1)=zero; oa(1)=zero
uc(nsig_pbl)=zero; vc(nsig_pbl)=zero; oc(nsig_pbl)=zero
do j=jstart,jstop
do i=1,lat2
! Background fields saved from the nonlinear model
do k=1,nsig_pbl
t_bg(k)=tges0(i,j,k)
u_bg(k)=uges0(i,j,k)
v_bg(k)=vges0(i,j,k)
o1_bg(k)=oges1(i,j,k)
u1_bg(k)=uges1(i,j,k)
v1_bg(k)=vges1(i,j,k)
p_bg(k)=pges0(i,j,k)
o_bg(k)=oges0(i,j,k)
zi_bg(k)=zi(i,j,k)
rdzi_bg(k)=rdzi(i,j,k)
end do
p_bg(nsig_pbl+1)=pges0(i,j,nsig_pbl+1)
zi_bg(nsig_pbl+1)=zi(i,j,nsig_pbl+1)
do k=2,nsig_pbl
rdzl_bg(k)=rdzl(i,j,k)
dodz_bg(k)=dodz(i,j,k)
dudz_bg(k)=dudz(i,j,k)
dvdz_bg(k)=dvdz(i,j,k)
end do
do k=2,nsig_pbl
zmix_bg(k)=zi(i,j,k)-zi(i,j,1)
lmix_bg(k)=karm*zmix_bg(k)/(one+karm0*zmix_bg(k))
ssq_bg(k)=dudz_bg(k)**2+dvdz_bg(k)**2
if(ssq_bg(k) < eps_m) then
ri_bg(k)=two*grav*dodz(i,j,k)/((o_bg(k)+o_bg(k-1))*eps_m)
else
ri_bg(k)=two*grav*dodz(i,j,k)/((o_bg(k)+o_bg(k-1))*ssq_bg(k))
end if
if( ri_bg(k) < zero) then
rho_bg(k)=sqrt(one-r20*ri_bg(k))
else
rho_bg(k)=one/(one+five*ri_bg(k))**2
end if
if(ssq_bg(k) < eps_m) then
km_bg(k)=lmix_bg(k)**2 *sqrt(eps_m )*rho_bg(k)
else
km_bg(k)=lmix_bg(k)**2 *sqrt(ssq_bg(k))*rho_bg(k)
end if
end do
ck_bg(1)=zero
do k=2,nsig_pbl
ck_bg(k)=km_bg(k)*rdzl_bg(k) ! << NL model repeat >>
end do
ck_bg(nsig_pbl+1)=zero
! Preliminaries
do k=nsig_pbl,1,-1 ! << NL model repeat >>
aux_bg(k)=time_step_half*rdzi_bg(k)
ckplus_bg=ck_bg(k+1)+ck_bg(k)
ckmnus_bg=(ck_bg(k+1)-ck_bg(k))*epxilon
ckdif_bg(k)=ckplus_bg-ckmnus_bg
cksum_bg(k)=ckplus_bg+ckmnus_bg
acoef=aux_bg(k)*ckdif_bg(k)
ccoef=aux_bg(k)*cksum_bg(k)
bcoef=-acoef-ccoef
acoef0(k)=-acoef*beta
bcoef0(k)=-bcoef*beta-one
ccoef0(k)=-ccoef*beta
acoef1(k)=acoef*alph
bcoef1(k)=bcoef*alph-one
ccoef1(k)=ccoef*alph
end do
wu(1:nsig_pbl)=alph*u1_bg(1:nsig_pbl)+beta*u_bg(1:nsig_pbl)
wv(1:nsig_pbl)=alph*v1_bg(1:nsig_pbl)+beta*v_bg(1:nsig_pbl)
wo(1:nsig_pbl)=alph*o1_bg(1:nsig_pbl)+beta*o_bg(1:nsig_pbl)
ua(2:nsig_pbl)=wu(1:nsig_pbl-1)
ub(1:nsig_pbl)=wu(1:nsig_pbl)
uc(1:nsig_pbl-1)=wu(2:nsig_pbl)
va(2:nsig_pbl)=wv(1:nsig_pbl-1)
vb(1:nsig_pbl)=wv(1:nsig_pbl)
vc(1:nsig_pbl-1)=wv(2:nsig_pbl)
oa(2:nsig_pbl)=wo(1:nsig_pbl-1)
ob(1:nsig_pbl)=wo(1:nsig_pbl)
oc(1:nsig_pbl-1)=wo(2:nsig_pbl)
! Perturbation fields
do k=1,nsig_pbl
t_ad(k)=t(i,j,k)
u_ad(k)=u(i,j,k)
v_ad(k)=v(i,j,k)
end do
!! p_ad(nsig_pbl+1)=zero
p_ad=zero
rdzl_ad(1)=zero
zi_ad(1)=zero
!(5) Adjoint of trapezoidal scheme
!5.4) Adjoint of perturbation of temperature
do k=nsig_pbl,1,-1
ax1=t_bg(k)/o_bg(k)
ax2=t_bg(k)/(p_bg(k)+p_bg(k+1))*rd_over_cp
o_ad(k)=ax1*t_ad(k)
p_ad(k)=ax2*t_ad(k)
p_ad(k+1)=ax2*t_ad(k)+p_ad(k+1)
end do
!5.3) Adjoint of tridiagonal system
fo(1:nsig_pbl)=o_ad(1:nsig_pbl)
fu(1:nsig_pbl)=u_ad(1:nsig_pbl)
fv(1:nsig_pbl)=v_ad(1:nsig_pbl)
call tridiag_ad(acoef1,bcoef1,ccoef1,fo,nsig_pbl)
call tridiag_ad(acoef1,bcoef1,ccoef1,fu,nsig_pbl)
call tridiag_ad(acoef1,bcoef1,ccoef1,fv,nsig_pbl)
!5.2) Adjoint of u,v,o
call multi_tridiag_ad(acoef0,bcoef0,ccoef0,fo,o_ad,nsig_pbl)
call multi_tridiag_ad(acoef0,bcoef0,ccoef0,fu,u_ad,nsig_pbl)
call multi_tridiag_ad(acoef0,bcoef0,ccoef0,fv,v_ad,nsig_pbl)
!5.1) Adjoint of a_ad,b_ad,c_ad
do k=nsig_pbl,1,-1
a_ad(k)=-(ua(k)*fu(k)+va(k)*fv(k)+oa(k)*fo(k))
b_ad(k)=-(ub(k)*fu(k)+vb(k)*fv(k)+ob(k)*fo(k))
c_ad(k)=-(uc(k)*fu(k)+vc(k)*fv(k)+oc(k)*fo(k))
end do
ck_ad(nsig_pbl+1)=zero
do k=nsig_pbl,1,-1
a_ad(k)=a_ad(k)-b_ad(k)
c_ad(k)=c_ad(k)-b_ad(k)
ckplus_ad=aux_bg(k)*( a_ad(k)+c_ad(k) )
ckmnus_ad=aux_bg(k)*(-a_ad(k)+c_ad(k) )
aux_ad =ckdif_bg(k)*a_ad(k)+cksum_bg(k)*c_ad(k)
ck_ad(k) = ckplus_ad-epxilon*ckmnus_ad
ck_ad(k+1) = ckplus_ad+epxilon*ckmnus_ad + ck_ad(k+1)
rdzi_ad(k)=aux_ad*time_step_half
end do
do k=nsig_pbl,2,-1
km_ad(k)=ck_ad(k)*rdzl_bg(k)
rdzl_ad(k)=ck_ad(k)*km_bg(k)
end do
km_ad(1)=zero
rdzl_ad(1)=zero
!(4) Adjoint of perturbation of mixing coefficients
do k=nsig_pbl,2,-1
!Preliminaries
ax1=one/(o_bg(k)+o_bg(k-1))
if(ssq_bg(k)<eps_m) then
!Eq(46)
lmix_ad(k)=two*km_bg(k)/lmix_bg(k)*km_ad(k)
rho_ad = km_bg(k)/rho_bg(k) *km_ad(k)
ssq_ad = zero !! test
else
!Eq(44)
rssq_bg=one/ssq_bg(k)
lmix_ad(k)=two*km_bg(k)/lmix_bg(k)*km_ad(k)
rho_ad = km_bg(k)/rho_bg(k) *km_ad(k)
ssq_ad = km_bg(k)*rssq_bg *km_ad(k)
end if
!Eq(43)
if(ri_bg(k)<zero)then
ri_ad(k)=-r10/rho_bg(k) *rho_ad
else
ri_ad(k)=-r10*(sqrt(rho_bg(k)))**3*rho_ad
end if
if(ssq_bg(k)<eps_m) then
!Eq(45)
dodz_ad(k)=(ri_bg(k)/dodz_bg(k))*ri_ad(k)
o_ad(k )=-ri_bg(k)*ax1*ri_ad(k)+o_ad(k)
o_ad(k-1)=-ri_bg(k)*ax1*ri_ad(k)+o_ad(k-1) !<<test
else
!Eq(42)
dodz_ad(k)=(ri_bg(k)/dodz_bg(k))*ri_ad(k)
o_ad(k )=-ri_bg(k)*ax1*ri_ad(k)+o_ad(k)
o_ad(k-1)=-ri_bg(k)*ax1*ri_ad(k)+o_ad(k-1) !<<test
ssq_ad =-ri_bg(k)*rssq_bg*two*ri_ad(k)+ ssq_ad
end if
!Eq(41)
zmix_ad=karm*lmix_ad(k)/(one+karm0*zmix_bg(k))**2
!Eq(40)
dudz_ad(k)=dudz_bg(k)*ssq_ad
dvdz_ad(k)=dvdz_bg(k)*ssq_ad
!Eq(39)
zi_ad(k)=zmix_ad
end do
!(3) Adjoint of perturbation of vertical gradiens
do k=nsig_pbl,2,-1
u_ad(k-1) = -dudz_ad(k)*rdzl_bg(k)+u_ad(k-1) !!<test
u_ad(k ) = dudz_ad(k)*rdzl_bg(k)+u_ad(k)
v_ad(k-1) = -dvdz_ad(k)*rdzl_bg(k)+v_ad(k-1) !!<test
v_ad(k ) = dvdz_ad(k)*rdzl_bg(k)+v_ad(k)
o_ad(k-1) = -dodz_ad(k)*rdzl_bg(k)+o_ad(k-1)
o_ad(k ) = dodz_ad(k)*rdzl_bg(k)+o_ad(k)
rdzl_ad(k)= &
+(u_bg(k)-u_bg(k-1))*dudz_ad(k) &
+(v_bg(k)-v_bg(k-1))*dvdz_ad(k) &
+(o_bg(k)-o_bg(k-1))*dodz_ad(k) &
+rdzl_ad(k) ! <<test
end do
!(2) Adjoint of perturbation of heights
zl_ad(nsig_pbl)=zero
do k=nsig_pbl,2,-1
zl_ad(k-1)= rdzl_bg(k)**2*rdzl_ad(k)
zl_ad(k )=-rdzl_bg(k)**2*rdzl_ad(k)+zl_ad(k)
end do
zi_ad(nsig_pbl+1)=zero
zi_ad(1)=zero
do k=nsig_pbl,1,-1
zi_ad(k )=half*zl_ad(k)+rdzi_bg(k)**2*rdzi_ad(k) + zi_ad(k)
zi_ad(k+1)=half*zl_ad(k)-rdzi_bg(k)**2*rdzi_ad(k) + zi_ad(k+1)
end do
dzi_ad(nsig_pbl+1)=zero
do k=nsig_pbl,1,-1
dzi_ad(k)=dzi_ad(k+1)+zi_ad(k+1)
end do
do k=nsig_pbl,1,-1
ax1 = four*rd_over_g*t_bg(k)/(p_bg(k)+p_bg(k+1))**2
t_ad(k)=(zi_bg(k+1)-zi_bg(k))/t_bg(k)*dzi_ad(k)
p_ad(k )= ax1 *p_bg(k+1)*dzi_ad(k)+p_ad(k)
p_ad(k+1)=-ax1 *p_bg(k)*dzi_ad(k)+p_ad(k+1)
end do
!(1) Adjoint of perturbation of potential temperature
do k=nsig_pbl,1,-1
ax1=o_bg(k)/t_bg(k)
ax2=o_bg(k)/(p_bg(k)+p_bg(k+1))*rd_over_cp
t_ad(k)= ax1*o_ad(k) + t_ad(k)
p_ad(k )=-ax2*o_ad(k)+p_ad(k )
p_ad(k+1)=-ax2*o_ad(k)+p_ad(k+1)
end do
!(0) Copy back to 3d output arrays
do k=nsig_pbl,1,-1
u(i,j,k)=u_ad(k)
v(i,j,k)=v_ad(k)
t(i,j,k)=t_ad(k)
prs(i,j,k)=p_ad(k)
end do
prs(i,j,nsig_pbl+1)=p_ad(nsig_pbl+1)
end do
end do
ps=zero
call getprs_bck_ad(ps,t,prs)
return
end subroutine pbl_ad
subroutine tridiag_ad(a,b,c,f,jmx)
! Adjoint of solver for a standard tridiagonal system
! - 'Thomas' tridiagonal algorithm -
! (Adapted from Durran (1999), p. 440)
! a - sub (lower) diagonal
! b - center diagonal
! c - super (upper) diagonal
! f - right hand side and solution
!
! ( a(1) and c(jmx) need not to be initialized )
use kinds,only: r_kind,i_kind,r_quad
use constants,only: one,zero
implicit none
integer(i_kind), intent(in):: jmx
real(r_kind), dimension(jmx), intent(in):: a,b,c
real(r_kind), dimension(jmx), intent(inout):: f
real(r_kind), dimension(jmx):: q
real(r_kind) p
integer(i_kind) j
! Forward elimination sweep
q(1)=-a(2)/b(1)
f(1)= f(1)/b(1)
do j=2,jmx
p=one/(b(j)+c(j-1)*q(j-1))
q(j)=-a(j+1)*p
f(j)=( f(j)-c(j-1)*f(j-1) )*p
end do
! Backward pass
do j=jmx-1,1,-1
f(j)=f(j)+q(j)*f(j+1)
end do
end subroutine tridiag_ad
subroutine multi_tridiag_ad(a,b,c,u,f,jmx)
! multimply adjoint of tridiagonal matrix (a,b,c) with
! vector u and store result in vector f
use kinds,only: r_kind,i_kind
implicit none
integer(i_kind), intent(in):: jmx
real(r_kind), dimension(jmx), intent(in):: a,b,c,u
real(r_kind), dimension(jmx), intent(out):: f
integer(i_kind) j
do j=2,jmx-1
f(j)=u(j-1)*c(j-1)+u(j)*b(j)+u(j+1)*a(j+1)
end do
f(1)= u(1)*b(1)+u(1+1)*a(1+1)
f(jmx)=u(jmx-1)*c(jmx-1)+u(jmx)*b(jmx)
end subroutine multi_tridiag_ad