subroutine pbl_tl(u,v,t,ps,jstart,jstop) !$$$ subprogram documentation block ! . . . ! subprogram: pbl_tl ! ! prgrmmr: m. rancic ! ! abstract: tangent linear of pbl ! ! program history log: ! 2008-04-02 safford -- add subprogram doc block, rm unused uses ! ! input argument list: ! u - ! v - ! ps - ! t - ! ! output 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: uges0,vges0,oges0,pges0,tges0,uges1,vges1,oges1 use pblmod, only: dudz,dvdz,dodz,zi,rdzi,rdzl,eps_m,nsig_pbl use pblmod, only: lmbd,karm,karm0,alph,beta,epxilon use tends4pertmod, only: time_step_half implicit none ! Declare passed variables real(r_kind),dimension(lat2,lon2) ,intent(in ):: ps real(r_kind),dimension(lat2,lon2,nsig),intent(inout):: u,v,t 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_tl,dzi_tl real(r_kind),dimension(nsig_pbl):: rdzl_tl,dudz_tl,dvdz_tl,dodz_tl,km_tl 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,rho_bg real(r_kind),dimension(2:nsig_pbl):: lmix_tl,ri_tl real(r_kind),dimension(nsig_pbl):: u_tl,v_tl,t_tl,zl_tl,o_tl real(r_kind),dimension(nsig_pbl+1):: p_tl,zi_tl,ck_tl,ck_bg real(r_kind),dimension(nsig_pbl):: u_out,v_out,t_out real(r_kind) rssq_bg,ssq_tl real(r_kind) zmix_tl real(r_kind):: acoef,bcoef,ccoef real(r_kind):: ax1,ax2 real(r_kind):: aux_tl,ckplus_bg,ckmnus_bg,ckplus_tl,ckmnus_tl real(r_kind) rho_tl real(r_kind),dimension(nsig_pbl):: aux_bg,ckdif_bg,cksum_bg 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_tl,b_tl,c_tl 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 call getprs_bck_tl(ps,t,prs) ua(1)=zero; va(1)=zero; oa(1)=zero uc(nsig_pbl)=zero; vc(nsig_pbl)=zero; oc(nsig_pbl)=zero dudz_tl(1)=zero; dvdz_tl(1)=zero; dodz_tl(1)=zero zi_tl(1)=zero; rdzl_tl(1)=zero; km_tl(1)=zero do j=jstart,jstop do i=1,lat2 u_out=zero; v_out=zero; t_out=zero ! 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) o_bg(k)=oges0(i,j,k) p_bg(k)=pges0(i,j,k) u1_bg(k)=uges1(i,j,k) v1_bg(k)=vges1(i,j,k) o1_bg(k)=oges1(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 do k=2,nsig_pbl ck_bg(k)=km_bg(k)*rdzl_bg(k) ! << NL model repeat >> end do ck_bg(1)=zero ; ck_bg(nsig_pbl+1)=zero ! Preliminaries 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) do k=1,nsig_pbl 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 ! Start tangent linear ! Perturbation fields do k=1,nsig_pbl t_tl(k)=t(i,j,k) u_tl(k)=u(i,j,k) v_tl(k)=v(i,j,k) p_tl(k)=prs(i,j,k) end do p_tl(nsig_pbl+1)=prs(i,j,nsig_pbl+1) !(1) Perturbation of potential temperature do k=1,nsig_pbl ax1=o_bg(k)/t_bg(k) ax2=o_bg(k)/(p_bg(k)+p_bg(k+1))*rd_over_cp o_tl(k)=ax1*t_tl(k)-ax2*(p_tl(k)+p_tl(k+1)) end do !(2) Perturbation of heights do k=1,nsig_pbl dzi_tl(k)=(zi_bg(k+1)-zi_bg(k))/t_bg(k)*t_tl(k) & -four*rd_over_g*t_bg(k)/(p_bg(k)+p_bg(k+1))**2 & *(p_bg(k)*p_tl(k+1)-p_bg(k+1)*p_tl(k)) end do do k=1,nsig_pbl zi_tl(k+1)=zi_tl(k)+dzi_tl(k) end do do k=1,nsig_pbl zl_tl(k)=half*(zi_tl(k+1)+zi_tl(k)) rdzi_tl(k)=-rdzi_bg(k)**2 *(zi_tl(k+1)-zi_tl(k)) end do do k=2,nsig_pbl rdzl_tl(k)=-rdzl_bg(k)**2 *(zl_tl(k)-zl_tl(k-1)) end do !(3) Perturbation of vertical gradients do k=2,nsig_pbl dodz_tl(k)=(o_tl(k)-o_tl(k-1))*rdzl_bg(k) & +(o_bg(k)-o_bg(k-1))*rdzl_tl(k) dudz_tl(k)=(u_tl(k)-u_tl(k-1))*rdzl_bg(k) & +(u_bg(k)-u_bg(k-1))*rdzl_tl(k) dvdz_tl(k)=(v_tl(k)-v_tl(k-1))*rdzl_bg(k) & +(v_bg(k)-v_bg(k-1))*rdzl_tl(k) end do !(4) Perturbation of mixing coefficients do k=2,nsig_pbl ax1=one/(o_bg(k)+o_bg(k-1)) zmix_tl=zi_tl(k) lmix_tl(k)=karm*zmix_tl/(one+karm0*zmix_bg(k))**2 ssq_tl=dudz_bg(k)*dudz_tl(k)+dvdz_bg(k)*dvdz_tl(k) if(ssq_bg(k)< eps_m) then ri_tl(k)=ri_bg(k)*( dodz_tl(k)/dodz_bg(k)-& ax1*(o_tl(k)+o_tl(k-1))) else rssq_bg=one/ssq_bg(k) ri_tl(k)=ri_bg(k)*( dodz_tl(k)/dodz_bg(k)-& ax1*(o_tl(k)+o_tl(k-1))-& rssq_bg*two*ssq_tl ) end if if(ri_bg(k)<zero)then rho_tl=-r10/rho_bg(k) *ri_tl(k) else rho_tl=-r10*(sqrt(rho_bg(k)))**3*ri_tl(k) end if if(ssq_bg(k)< eps_m) then km_tl(k)=km_bg(k)*( two*lmix_tl(k)/lmix_bg(k) & +rho_tl/rho_bg(k) ) else km_tl(k)=km_bg(k)*( two*lmix_tl(k)/lmix_bg(k) & +ssq_tl*rssq_bg+rho_tl/rho_bg(k) ) end if end do do k=2,nsig_pbl ck_tl(k)=ck_bg(k)*(km_tl(k)/km_bg(k)+rdzl_tl(k)/rdzl_bg(k)) end do ck_tl(1)=zero ; ck_tl(nsig_pbl+1)=zero !(5) Perturbation of trapezoidal scheme ! Start tangent linear !(5.1) do k=1,nsig_pbl aux_tl=time_step_half*rdzi_tl(k) ckplus_tl=ck_tl(k+1)+ck_tl(k) ckmnus_tl=(ck_tl(k+1)-ck_tl(k))*epxilon a_tl(k)=aux_bg(k)*(ckplus_tl-ckmnus_tl)+ & aux_tl*ckdif_bg(k) c_tl(k)=aux_bg(k)*(ckplus_tl+ckmnus_tl)+ & aux_tl*cksum_bg(k) b_tl(k)=-a_tl(k)-c_tl(k) end do do k=1,nsig_pbl pu(k)=ua(k)*a_tl(k)+ub(k)*b_tl(k)+uc(k)*c_tl(k) pv(k)=va(k)*a_tl(k)+vb(k)*b_tl(k)+vc(k)*c_tl(k) po(k)=oa(k)*a_tl(k)+ob(k)*b_tl(k)+oc(k)*c_tl(k) end do !(5.2) call multi_tridiag(acoef0,bcoef0,ccoef0,u_tl,qu,nsig_pbl) call multi_tridiag(acoef0,bcoef0,ccoef0,v_tl,qv,nsig_pbl) call multi_tridiag(acoef0,bcoef0,ccoef0,o_tl,qo,nsig_pbl) fu(1:nsig_pbl)=qu(1:nsig_pbl)-pu(1:nsig_pbl) fv(1:nsig_pbl)=qv(1:nsig_pbl)-pv(1:nsig_pbl) fo(1:nsig_pbl)=qo(1:nsig_pbl)-po(1:nsig_pbl) !(5.3) call tridiag(acoef1,bcoef1,ccoef1,fu,nsig_pbl) call tridiag(acoef1,bcoef1,ccoef1,fv,nsig_pbl) call tridiag(acoef1,bcoef1,ccoef1,fo,nsig_pbl) ! Update perturbations of wind and potential temperature u_out(1:nsig_pbl)=fu(1:nsig_pbl) v_out(1:nsig_pbl)=fv(1:nsig_pbl) o_tl(1:nsig_pbl)=fo(1:nsig_pbl) !(5.4) ! Update perturbation of temperature do k=1,nsig_pbl ax1=t_bg(k)/o_bg(k) ax2=t_bg(k)/(p_bg(k)+p_bg(k+1))*rd_over_cp t_out(k)=ax1*o_tl(k)+ax2*(p_tl(k)+p_tl(k+1)) end do do k=1,nsig_pbl u(i,j,k)=u_out(k) v(i,j,k)=v_out(k) t(i,j,k)=t_out(k) end do end do end do return end subroutine pbl_tl subroutine tridiag(a,b,c,f,jmx) ! Solves 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)=-c(1)/b(1) f(1)= f(1)/b(1) do j=2,jmx p=one/(b(j)+a(j)*q(j-1)) q(j)=-c(j)*p f(j)=( f(j)-a(j)*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 subroutine multi_tridiag(a,b,c,u,f,jmx) ! multiply a tridiagonal matrix (a,b,c) with a vector u ! 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)*a(j)+u(j)*b(j)+u(j+1)*c(j) end do f(1)= u(1)*b(1)+u(1+1)*c(1) f(jmx)=u(jmx-1)*a(jmx)+u(jmx)*b(jmx) end subroutine multi_tridiag