subroutine getvvel(t,t_thor,prsth,prdif,what) !$$$ subprogram documentation block ! . . . . ! subprogram: get_vvel compute vertical velocity ! prgmmr: kleist org: np20 date: 2007-05-08 ! ! abstract: get model level vertical velocity for generalized coordinate ! ! program history log: ! 2007-05-08 kleist ! ! usage: ! input argument list: ! t - temperature ! t_thor - horizontal component to temperature tendency (advection) ! prsth - horizontal component to pressure tendency (advection) ! prdif - delta pressure ! ! output argument list: ! what - model level vertical velocity ! ! notes: ! see NCEP Office Note 445 (Juang 2005) ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds,only: r_kind,i_kind use constants, only: ione,zero,one,two,rd,cp,half use gridmod, only: lat2,lon2,nsig,bk5,ck5,tref5 use tendsmod, only: adiag9,bdiag9,cdiag9,factk9,wint9,wint9_f,& r_bdiag9 implicit none ! Declare passed variables: real(r_kind),dimension(lat2,lon2,nsig) ,intent(in ) :: t,t_thor,prdif real(r_kind),dimension(lat2,lon2,nsig+ione),intent(in ) :: prsth real(r_kind),dimension(lat2,lon2,nsig+ione),intent( out) :: what ! Declare local variables: real(r_kind),dimension(lat2,lon2,nsig):: r_prdif,tsum,t_tsum,tdiff real(r_kind) kapr integer(i_kind):: i,j,k ! IN: ! t: temperature ! t_t: horizontal part of temperature tendency ! prsth: horizontal part of pressure tendency ! prdif: P(k)-P(k+1) ! OUT: ! what: model level vertical velocity ! ***NOTES*** on vertical velocity for general coordinates ! general formula for pressure (on 445, eq 7.1) ! p(i,j,k) = ak5(k) + bk5(k)*psfc(i,j)+ck5(k)*(t(i,j,k)/tref5(k))**(cp/rd) ! following taken from on445, eq 7.8a ! solve adiag(k)*w(k-1) + bdiag(k)*w(k) + cdiag(k)*w(k+1) = rhs(k) ! ! factor = (ck5(k)/(t(i,j,k-1)+t(i,j,k))) * (cp/(2*rd)) * & ! ( (t(i,j,k-1)+t(i,j,k))/(tref(k-1)+tref(k)) )**(cp/rd) ! Multiplies w(k-1) : ! adiag(k)=factor*(t(i,j,k-2)-t(i,j,k-1))*r_prdif9(i,j,k-1) ! Multiplies w(k) : ! bdiag(k)=factor*(t(i,j,k-1)-t(i,j,k))*(r_prdif9(i,j,k-1)+r_prdif9(i,j,k)) - one ! Multiplies w(k+1) : ! cdiag(k)=factor*(t(i,j,k)-t(i,j,k+1))*r_prdif9(i,j,k) ! Right hand forcing: ! rhs(k)=bk5*prsth9(i,j,1)-prsth9(i,j,k)+2*factor*(t_t(i,j,k-1)+t_t(i,j,k)) ! ! to solve above tridiagonal matrix, do following: ! **** do forward elimination in a loop from k=2 to nsig, ! do k=3,nsig ! rhs(k)=rhs(k)-adiag(k)*rhs(k-1)/bdiag(k-1) ! bdiag(k)=bdiag(k)-adiag(k)*cdiag(k-1)/bdiag(k-1) ! end do ! ! **** then backward substitution in loop from k=nsig to k=2. ! w(nsig)=rhs(nsig)/bdiag(nsig) ! do k=nsig-1,2,-1 ! w(k)=(rhs(k)-cdiag(k)*w(k+1))/bdiag(k) ! end do ! kapr=cp/rd do k=1,nsig+ione do j=1,lon2 do i=1,lat2 what(i,j,k)=zero wint9(i,j,k)=zero wint9_f(i,j,k)=zero end do end do end do do k=1,nsig do j=1,lon2 do i=1,lat2 adiag9(i,j,k)=zero bdiag9(i,j,k)=zero cdiag9(i,j,k)=zero factk9(i,j,k)=zero r_bdiag9(i,j,k)=zero tsum(i,j,k)=zero t_tsum(i,j,k)=zero tdiff(i,j,k)=zero end do end do end do do k=1,nsig-ione do j=1,lon2 do i=1,lat2 tsum(i,j,k)=t(i,j,k)+t(i,j,k+ione) t_tsum(i,j,k)=t_thor(i,j,k)+t_thor(i,j,k+ione) tdiff(i,j,k)=t(i,j,k)-t(i,j,k+ione) end do end do end do do k=1,nsig do j=1,lon2 do i=1,lat2 r_prdif(i,j,k)=one/prdif(i,j,k) end do end do end do ! forward elimination: k=2_i_kind do j=1,lon2 do i=1,lat2 factk9(i,j,k) = ck5(k)*half*kapr/tsum(i,j,k-ione) * & ( half*tsum(i,j,k-ione)/tref5(k-ione) )**kapr bdiag9(i,j,k)=factk9(i,j,k)*tdiff(i,j,k-ione)*(r_prdif(i,j,k-ione) + & r_prdif(i,j,k)) - one cdiag9(i,j,k)=factk9(i,j,k)*tdiff(i,j,k)*r_prdif(i,j,k) wint9(i,j,k)=bk5(k)*prsth(i,j,1)-prsth(i,j,k)+two*factk9(i,j,k)*t_tsum(i,j,k-ione) wint9_f(i,j,k)=wint9(i,j,k) r_bdiag9(i,j,k)=one/bdiag9(i,j,k) end do end do do k=3,nsig do j=1,lon2 do i=1,lat2 factk9(i,j,k) = ck5(k)*half*kapr/tsum(i,j,k-ione) * & ( half*tsum(i,j,k-ione)/tref5(k-ione) )**kapr adiag9(i,j,k)=factk9(i,j,k)*tdiff(i,j,k-2_i_kind)*r_prdif(i,j,k-ione) bdiag9(i,j,k)=factk9(i,j,k)*tdiff(i,j,k-ione)*(r_prdif(i,j,k-ione)+ & r_prdif(i,j,k)) - one cdiag9(i,j,k)=factk9(i,j,k)*tdiff(i,j,k)*r_prdif(i,j,k) wint9 (i,j,k)=bk5(k)*prsth(i,j,1)-prsth(i,j,k)+two*factk9(i,j,k)*t_tsum(i,j,k-ione) wint9 (i,j,k)=wint9 (i,j,k)-adiag9(i,j,k)*wint9 (i,j,k-ione)*r_bdiag9(i,j,k-ione) bdiag9(i,j,k)=bdiag9(i,j,k)-adiag9(i,j,k)*cdiag9(i,j,k-ione)*r_bdiag9(i,j,k-ione) wint9_f(i,j,k)=wint9(i,j,k) r_bdiag9(i,j,k)=one/bdiag9(i,j,k) end do end do end do ! back substitution k=nsig do j=1,lon2 do i=1,lat2 wint9(i,j,k)=wint9(i,j,k)*r_bdiag9(i,j,k) end do end do do k=nsig-ione,2,-1 do j=1,lon2 do i=1,lat2 wint9(i,j,k)=(wint9(i,j,k)-cdiag9(i,j,k)*wint9(i,j,k+ione))*r_bdiag9(i,j,k) end do end do end do do k=1,nsig+ione do j=1,lon2 do i=1,lat2 what(i,j,k)=wint9(i,j,k) end do end do end do return end subroutine getvvel subroutine getvvel_tl(t,t_thor,t_thor9,prsth,prdif,what) !$$$ subprogram documentation block ! . . . . ! subprogram: get_vvel_tl tlm of getvvel ! prgmmr: kleist org: np20 date: 2007-05-08 ! ! abstract: tangent linear calculation of vertical velocity ! ! program history log: ! 2007-05-08 kleist ! 2008-06-04 safford - rm unused uses ! ! usage: ! input argument list: ! t - temperature ! t_thor - horizontal component to temperature tendency (advection) ! t_thor9 - basic state horizontal component to temperature tendency (advection) ! prsth - horizontal component to pressure tendency (advection) ! prdif - delta pressure ! ! output argument list: ! what - perturbation model level vertical velocity ! ! notes: ! see NCEP Office Note 445 (Juang 2005) ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds,only: r_kind,i_kind use constants, only: ione,zero,one,two,rd,cp,half use gridmod, only: lat2,lon2,nsig,bk5,ck5,tref5 use guess_grids, only: ges_tv,ntguessig use tendsmod, only: prdif9,r_prdif9,adiag9,bdiag9,cdiag9,factk9,& r_bdiag9,wint9,wint9_f implicit none ! Declare passed variables: real(r_kind),dimension(lat2,lon2,nsig) ,intent(in ) :: t,t_thor,t_thor9,prdif real(r_kind),dimension(lat2,lon2,nsig+ione),intent(in ) :: prsth real(r_kind),dimension(lat2,lon2,nsig+ione),intent( out) :: what ! Declare local variables: real(r_kind),dimension(lat2,lon2,nsig):: tsum,t_tsum,tdiff,factk real(r_kind),dimension(lat2,lon2,nsig):: tsum9,t_tsum9,tdiff9 real(r_kind),dimension(lat2,lon2,nsig):: adiag,bdiag,cdiag real(r_kind) kapr,kaprm1,terma,termb real(r_kind) c1,c2,fprime integer(i_kind):: i,j,k,it ! Constants/Parameters: kapr=cp/rd kaprm1=kapr-one it=ntguessig do k=1,nsig+ione do j=1,lon2 do i=1,lat2 what(i,j,k)=zero end do end do end do do k=1,nsig do j=1,lon2 do i=1,lat2 adiag(i,j,k)=zero bdiag(i,j,k)=zero cdiag(i,j,k)=zero factk(i,j,k)=zero tsum(i,j,k)=zero t_tsum(i,j,k)=zero tdiff(i,j,k)=zero tsum9(i,j,k)=zero t_tsum9(i,j,k)=zero tdiff9(i,j,k)=zero end do end do end do ! do k=1,nsig-ione do j=1,lon2 do i=1,lat2 tsum(i,j,k)=t(i,j,k)+t(i,j,k+ione) t_tsum(i,j,k)=t_thor(i,j,k)+t_thor(i,j,k+ione) tdiff(i,j,k)=t(i,j,k)-t(i,j,k+ione) tsum9(i,j,k)=ges_tv(i,j,k,it)+ges_tv(i,j,k+ione,it) tdiff9(i,j,k)=ges_tv(i,j,k,it)-ges_tv(i,j,k+ione,it) t_tsum9(i,j,k)=t_thor9(i,j,k)+t_thor9(i,j,k+ione) end do end do end do ! forward elimination: k=2_i_kind do j=1,lon2 do i=1,lat2 c1=ck5(k)*half*kapr c2=half/tref5(k-ione) fprime=kapr*c2*tsum(i,j,k-ione)*((c2*tsum9(i,j,k-ione))**kaprm1) factk(i,j,k) = c1* ( (tsum9(i,j,k-ione)*fprime) - (tsum(i,j,k-ione)* & ((c2*tsum9(i,j,k-ione))**kapr)) ) / (tsum9(i,j,k-1)**two) terma = ( prdif9(i,j,k-ione)*(factk(i,j,k)*tdiff9(i,j,k-ione) + tdiff(i,j,k-ione)* & factk9(i,j,k)) - (factk9(i,j,k)*tdiff9(i,j,k-ione)*prdif(i,j,k-ione)) ) * & (r_prdif9(i,j,k-ione)**two) termb= ( prdif9(i,j,k)*(factk(i,j,k)*tdiff9(i,j,k-ione) + tdiff(i,j,k-ione)* & factk9(i,j,k)) - (factk9(i,j,k)*tdiff9(i,j,k-ione)*prdif(i,j,k)) ) * & (r_prdif9(i,j,k)**two) bdiag(i,j,k) = terma+termb cdiag(i,j,k) = ( prdif9(i,j,k)*(factk(i,j,k)*tdiff9(i,j,k) + & tdiff(i,j,k)*factk9(i,j,k)) - factk9(i,j,k)*tdiff9(i,j,k)* & prdif(i,j,k) ) * (r_prdif9(i,j,k)**two) what(i,j,k)=bk5(k)*prsth(i,j,1)-prsth(i,j,k) + two* & (factk9(i,j,k)*t_tsum(i,j,k-ione) + factk(i,j,k)*t_tsum9(i,j,k-ione)) end do end do do k=3,nsig do j=1,lon2 do i=1,lat2 c1=ck5(k)*half*kapr c2=half/tref5(k-ione) fprime=kapr*c2*tsum(i,j,k-ione)*((c2*tsum9(i,j,k-ione))**kaprm1) factk(i,j,k) = c1* ( (tsum9(i,j,k-ione)*fprime) - (tsum(i,j,k-ione)* & ((c2*tsum9(i,j,k-ione))**kapr)) ) / (tsum9(i,j,k-ione)**two) adiag(i,j,k) = ( prdif9(i,j,k-ione)*(factk(i,j,k)*tdiff9(i,j,k-2_i_kind) + & tdiff(i,j,k-2_i_kind)*factk9(i,j,k)) - factk9(i,j,k)*tdiff9(i,j,k-2_i_kind)* & prdif(i,j,k-ione) ) * (r_prdif9(i,j,k-ione)**two) terma = ( prdif9(i,j,k-ione)*(factk(i,j,k)*tdiff9(i,j,k-ione) + & tdiff(i,j,k-ione)*factk9(i,j,k)) - factk9(i,j,k)*tdiff9(i,j,k-ione)* & prdif(i,j,k-ione) ) * (r_prdif9(i,j,k-ione)**two) termb = ( prdif9(i,j,k)*(factk(i,j,k)*tdiff9(i,j,k-ione) + & tdiff(i,j,k-ione)*factk9(i,j,k)) - factk9(i,j,k)*tdiff9(i,j,k-ione)* & prdif(i,j,k) ) * (r_prdif9(i,j,k)**two) bdiag(i,j,k) = terma+termb cdiag(i,j,k) = ( prdif9(i,j,k)*(factk(i,j,k)*tdiff9(i,j,k) + & tdiff(i,j,k)*factk9(i,j,k)) - factk9(i,j,k)*tdiff9(i,j,k)* & prdif(i,j,k) ) * (r_prdif9(i,j,k)**two) what(i,j,k)=bk5(k)*prsth(i,j,1)-prsth(i,j,k) + two* & (factk9(i,j,k)*t_tsum(i,j,k-ione) + factk(i,j,k)*t_tsum9(i,j,k-ione)) what(i,j,k) = what(i,j,k) - ( (bdiag9(i,j,k-ione)*(adiag(i,j,k)*wint9_f(i,j,k-ione) + & what(i,j,k-ione)*adiag9(i,j,k)) - adiag9(i,j,k)*wint9_f(i,j,k-ione)* & bdiag(i,j,k-ione)) * (r_bdiag9(i,j,k-ione)**two) ) bdiag(i,j,k)=bdiag(i,j,k) - ( (bdiag9(i,j,k-ione)*(adiag(i,j,k)*cdiag9(i,j,k-ione) + & cdiag(i,j,k-ione)*adiag9(i,j,k)) - adiag9(i,j,k)*cdiag9(i,j,k-ione)* & bdiag(i,j,k-ione)) * (r_bdiag9(i,j,k-ione)**two) ) end do end do end do ! back substitution k=nsig do j=1,lon2 do i=1,lat2 what(i,j,k) = (bdiag9(i,j,k)*what(i,j,k) - wint9_f(i,j,k)*bdiag(i,j,k)) * & (r_bdiag9(i,j,k)**two) end do end do do k=nsig-ione,2,-1 do j=1,lon2 do i=1,lat2 what(i,j,k) = ( bdiag9(i,j,k)*(what(i,j,k)-(cdiag9(i,j,k)*what(i,j,k+ione) + & cdiag(i,j,k)*wint9(i,j,k+ione)) ) - bdiag(i,j,k)*(wint9_f(i,j,k)- & cdiag9(i,j,k)*wint9(i,j,k+ione)) ) * (r_bdiag9(i,j,k)**two) end do end do end do return end subroutine getvvel_tl subroutine getvvel_ad(t,t_thor,t_thor9,prsth,prdif,whatin) !$$$ subprogram documentation block ! . . . . ! subprogram: get_vvel_ad adjoint of getvvel ! prgmmr: kleist org: np20 date: 2007-05-08 ! ! abstract: adjoint of calculation of vertical velocity ! ! program history log: ! 2007-05-08 kleist ! 2008-06-04 safford - rm unused uses ! ! usage: ! input argument list: ! t - temperature ! t_thor - horizontal component to temperature tendency (advection) ! t_thor9 - basic state horizontal component to temperature tendency (advection) ! prsth - horizontal component to pressure tendency (advection) ! prdif - delta pressure ! what - model level vertical velocity ! ! output argument list: ! t - temperature ! t_thor - horizontal component to temperature tendency (advection) ! prsth - horizontal component to pressure tendency (advection) ! prdif - delta pressure ! ! notes: ! see NCEP Office Note 445 (Juang 2005) ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds,only: r_kind,i_kind use constants, only: ione,zero,one,two,rd,cp,half use gridmod, only: lat2,lon2,nsig,bk5,ck5,tref5 use guess_grids, only: ges_tv,ntguessig use tendsmod, only: prdif9,r_prdif9,adiag9,bdiag9,cdiag9,factk9,& r_bdiag9,wint9,wint9_f implicit none ! Declare passed variables: real(r_kind),dimension(lat2,lon2,nsig+ione),intent(in ) :: whatin real(r_kind),dimension(lat2,lon2,nsig) ,intent(in ) :: t_thor9 real(r_kind),dimension(lat2,lon2,nsig) ,intent(inout) :: t,t_thor,prdif real(r_kind),dimension(lat2,lon2,nsig+ione),intent(inout) :: prsth ! Declare local variables: real(r_kind),dimension(lat2,lon2,nsig+ione):: what real(r_kind),dimension(lat2,lon2,nsig):: tsum,t_tsum,tdiff,factk real(r_kind),dimension(lat2,lon2,nsig):: tsum9,t_tsum9,tdiff9 real(r_kind),dimension(lat2,lon2,nsig):: adiag,bdiag,cdiag real(r_kind) kapr,kaprm1,terma,termb real(r_kind) c1,c2,fprime,tmp1,tmp2 integer(i_kind):: i,j,k,it ! Constants/Parameters: kapr=cp/rd kaprm1=kapr-one it=ntguessig fprime=zero do k=1,nsig do j=1,lon2 do i=1,lat2 adiag(i,j,k)=zero bdiag(i,j,k)=zero cdiag(i,j,k)=zero factk(i,j,k)=zero tsum(i,j,k)=zero t_tsum(i,j,k)=zero tdiff(i,j,k)=zero tsum9(i,j,k)=zero t_tsum9(i,j,k)=zero tdiff9(i,j,k)=zero end do end do end do do k=1,nsig-ione do j=1,lon2 do i=1,lat2 tsum9 (i,j,k)=ges_tv (i,j,k,it)+ges_tv (i,j,k+ione,it) tdiff9 (i,j,k)=ges_tv (i,j,k,it)-ges_tv (i,j,k+ione,it) t_tsum9(i,j,k)=t_thor9(i,j,k )+t_thor9(i,j,k+ione ) end do end do end do do k=1,nsig+ione do j=1,lon2 do i=1,lat2 what(i,j,k)=whatin(i,j,k) end do end do end do ! back substitution do k=2,nsig-ione do j=1,lon2 do i=1,lat2 ! Use the _f (forwared/fixed) value of wint9 for the k-level here, but use the update ! real value of wint9 for the k+1 level tmp1=what(i,j,k)*(r_bdiag9(i,j,k)**two) what(i,j,k+ione) = what(i,j,k+ione) - tmp1*bdiag9(i,j,k)*cdiag9(i,j,k) cdiag(i,j,k) = cdiag(i,j,k) - tmp1*bdiag9(i,j,k)*wint9(i,j,k+ione) bdiag(i,j,k) = bdiag(i,j,k) - tmp1*(wint9_f(i,j,k)-cdiag9(i,j,k)*wint9(i,j,k+ione)) ! Update what(i,j,k) last what(i,j,k) = bdiag9(i,j,k)*tmp1 end do end do end do k=nsig do j=1,lon2 do i=1,lat2 ! use the _f (forward sub / fixed value) here at the k level tmp1=what(i,j,k)*(r_bdiag9(i,j,k)**two) bdiag(i,j,k) = bdiag(i,j,k) - wint9_f(i,j,k)*tmp1 what(i,j,k) = bdiag9(i,j,k)*tmp1 end do end do do k=nsig,3,-1 do j=1,lon2 do i=1,lat2 c1=ck5(k)*half*kapr c2=half/tref5(k-ione) tmp1=bdiag(i,j,k)*(r_bdiag9(i,j,k-ione)**two) adiag(i,j,k ) = adiag(i,j,k ) - tmp1*bdiag9(i,j,k-ione)*cdiag9(i,j,k-ione) cdiag(i,j,k-ione) = cdiag(i,j,k-ione) - tmp1*bdiag9(i,j,k-ione)*adiag9(i,j,k ) bdiag(i,j,k-ione) = bdiag(i,j,k-ione) + tmp1*adiag9(i,j,k )*cdiag9(i,j,k-ione) tmp2=what(i,j,k)*(r_bdiag9(i,j,k-ione)**two) adiag(i,j,k ) = adiag(i,j,k ) - tmp2*bdiag9(i,j,k-ione)*wint9_f(i,j,k-ione) bdiag(i,j,k-ione) = bdiag(i,j,k-ione) + tmp2*adiag9(i,j,k )*wint9_f(i,j,k-ione) what (i,j,k-ione) = what (i,j,k-ione) - tmp2*bdiag9(i,j,k-ione)*adiag9(i,j,k ) prsth(i,j,1) = prsth(i,j,1) + bk5(k)*what(i,j,k) prsth(i,j,k) = prsth(i,j,k) - what(i,j,k) t_tsum(i,j,k-ione) = t_tsum(i,j,k-ione) + two*factk9(i,j,k)*what(i,j,k) factk(i,j,k) = factk(i,j,k) + two*t_tsum9(i,j,k-ione)*what(i,j,k) tmp1=cdiag(i,j,k)*(r_prdif9(i,j,k)**two) factk(i,j,k) = factk(i,j,k) + tmp1*prdif9(i,j,k)*tdiff9(i,j,k) tdiff(i,j,k) = tdiff(i,j,k) + tmp1*prdif9(i,j,k)*factk9(i,j,k) prdif(i,j,k) = prdif(i,j,k) - tmp1*factk9(i,j,k)*tdiff9(i,j,k) tmp1=bdiag(i,j,k)*(r_prdif9(i,j,k)**two) factk(i,j,k ) = factk(i,j,k ) + tmp1*prdif9(i,j,k)*tdiff9(i,j,k-ione) tdiff(i,j,k-ione) = tdiff(i,j,k-ione) + tmp1*prdif9(i,j,k)*factk9(i,j,k ) prdif(i,j,k ) = prdif(i,j,k ) - tmp1*factk9(i,j,k)*tdiff9(i,j,k-ione) tmp2=bdiag(i,j,k)*(r_prdif9(i,j,k-ione)**two) factk(i,j,k ) = factk(i,j,k ) + tmp2*prdif9(i,j,k-ione)*tdiff9(i,j,k-ione) tdiff(i,j,k-ione) = tdiff(i,j,k-ione) + tmp2*prdif9(i,j,k-ione)*factk9(i,j,k ) prdif(i,j,k-ione) = prdif(i,j,k-ione) - tmp2*factk9(i,j,k )*tdiff9(i,j,k-ione) tmp1=adiag(i,j,k)*(r_prdif9(i,j,k-1)**two) factk(i,j,k ) = factk(i,j,k ) + tmp1*prdif9(i,j,k-ione)*tdiff9(i,j,k-2_i_kind) tdiff(i,j,k-2_i_kind) = tdiff(i,j,k-2_i_kind) + tmp1*prdif9(i,j,k-ione)*factk9(i,j,k ) prdif(i,j,k-ione ) = prdif(i,j,k-ione ) - tmp1*factk9(i,j,k )*tdiff9(i,j,k-2_i_kind) tmp2=c1*factk(i,j,k)/(tsum9(i,j,k-ione)**two) tsum(i,j,k-ione) = tsum(i,j,k-ione) - tmp2*((c2*tsum9(i,j,k-ione))**kapr) fprime = tmp2*tsum9(i,j,k-ione) tsum(i,j,k-ione) = tsum(i,j,k-ione) + kapr*c2*fprime*((c2*tsum9(i,j,k-ione))**kaprm1) end do end do end do k=2_i_kind do j=1,lon2 do i=1,lat2 c1=ck5(k)*half*kapr c2=half/tref5(k-ione) prsth(i,j,1) = prsth(i,j,1) + bk5(k)*what(i,j,k) prsth(i,j,k) = prsth(i,j,k) - what(i,j,k) t_tsum(i,j,k-ione) = t_tsum(i,j,k-ione) + two*factk9(i,j,k)*what(i,j,k) factk(i,j,k) = factk(i,j,k) + two*t_tsum9(i,j,k-ione)*what(i,j,k) tmp1=cdiag(i,j,k)*(r_prdif9(i,j,k)**two) factk(i,j,k) = factk(i,j,k) + tmp1*prdif9(i,j,k)*tdiff9(i,j,k) tdiff(i,j,k) = tdiff(i,j,k) + tmp1*prdif9(i,j,k)*factk9(i,j,k) prdif(i,j,k) = prdif(i,j,k) - tmp1*factk9(i,j,k)*tdiff9(i,j,k) terma = bdiag(i,j,k) ; termb=bdiag(i,j,k) tmp1=bdiag(i,j,k)*(r_prdif9(i,j,k)**two) factk(i,j,k) = factk(i,j,k) + tmp1*prdif9(i,j,k)*tdiff9(i,j,k-ione) tdiff(i,j,k-ione) = tdiff(i,j,k-ione) + tmp1*prdif9(i,j,k)*factk9(i,j,k) prdif(i,j,k) = prdif(i,j,k) - tmp1*factk9(i,j,k)*tdiff9(i,j,k-ione) tmp2=bdiag(i,j,k)*(r_prdif9(i,j,k-ione)**two) factk(i,j,k) = factk(i,j,k) + tmp2*prdif9(i,j,k-ione)*tdiff9(i,j,k-ione) tdiff(i,j,k-ione) = tdiff(i,j,k-ione) + tmp2*prdif9(i,j,k-ione)*factk9(i,j,k) prdif(i,j,k-ione) = prdif(i,j,k-ione) - tmp2*factk9(i,j,k)*tdiff9(i,j,k-ione) tmp1=c1*factk(i,j,k)/(tsum9(i,j,k-ione)**two) tsum(i,j,k-ione) = tsum(i,j,k-ione) - tmp1*((c2*tsum9(i,j,k-ione))**kapr) fprime = tsum9(i,j,k-ione)*tmp1 tsum(i,j,k-ione) = tsum(i,j,k-ione) + kapr*c2*fprime*((c2*tsum9(i,j,k-ione))**kaprm1) end do end do do k=1,nsig-ione do j=1,lon2 do i=1,lat2 t(i,j,k) = t(i,j,k) + tsum(i,j,k) + tdiff(i,j,k) t(i,j,k+ione) = t(i,j,k+ione) + tsum(i,j,k) - tdiff(i,j,k) t_thor(i,j,k) = t_thor(i,j,k) + t_tsum(i,j,k) t_thor(i,j,k+ione) = t_thor(i,j,k+ione) + t_tsum(i,j,k) end do end do end do return end subroutine getvvel_ad