#include <misc.h>#include <params.h>subroutine scanslt (ztodt   ,lat     ,dtr     ,iter    ,pmap    , &                    kdpmpf  ,kdpmph  ,lam     ,phi     ,dphi    , &                    lbasdy  ,lbasdz  ,lbasiy  ,lbasiz  ,lbassi  , &                    detam   ,detai   ,dlam    ,cwava   ,etamid  , &                    etaint  ,grfu    ,grfv    ,grlps1  ,grlps2  , &                    grt1    ,grt2    ,grq1    ,grq2    ,grfu1   , &                    grfu2   ,grfv1   ,grfv2   ,ps      ,u3      , &                    v3      ,t3      ,q3      ,lnpssld ,prhssld , &                    tarrsld ,parrsld ,n3      ,n3m1    ,u3sld   , &                    v3sld   ,etadot  ,nlon    )!-----------------------------------------------------------------------!! Purpose:! Interpolate terms for semi-lagrangian transport and SLD dynamics.! One latitude slice only!! Author:  J. Olson!!-----------------------------------------------------------------------!! $Id: scanslt.F90,v 1.14.2.1 2002/04/22 19:09:51 erik Exp $! $Author: erik $!!-----------------------------------------------------------------------  use precision  use pmgrid,      only: plon, plond, plev, plevp, plat, platd, beglat, endlat, beglatex, &                         endlatex, plndlv, i1, j1  use constituents,only: pcnst, pnats  use comslt,      only: qfcst, gamma, hw1lat  use rgrid,       only: nmmax  use pspect,      only: pmmax  use commap,      only: w, clat, t0  use prognostics, only: ptimelevels  use physconst,   only: cappa  use dynconst,    only: ra  implicit none#include <comctl.h>#include <comfft.h>#include <comhyb.h>!------------------------------Arguments--------------------------------!  real(r8), intent(in)   :: ztodt                 ! twice the time step unless nstep = 0  integer , intent(in)   :: lat                   ! latitude index  real(r8), intent(in)   :: dtr                   ! 1/dt  integer , intent(in)   :: iter                  ! number of iterations for trajectory  integer , intent(in)   :: pmap                  ! dimension of artificial array!                                                 ! used to locate vertical interval!                                                 ! in which departure point falls  integer , intent(in)   :: kdpmpf  (pmap)        ! mapping from artificial array to!                                                 ! model levels  integer , intent(in)   :: kdpmph  (pmap)        ! mapping from artificial array to!                                                 ! model interfaces  real(r8), intent(in)   :: lam     (plond,platd) ! longitude coordinates of model grid  real(r8), intent(in)   :: phi     (platd)       ! latitude  coordinates of model grid  real(r8), intent(in)   :: dphi    (platd)       ! latitudinal grid increments  real(r8), intent(in)   :: lbasdy  (4,2,platd)   ! basis functions for lat deriv est.  real(r8), intent(in)   :: lbasdz  (4,2,plev)    ! basis functions for vert deriv est.  real(r8), intent(in)   :: lbasiy  (4,2,platd)   ! basis functions for Lagrange interp  real(r8), intent(in)   :: lbasiz  (4,2,plev)    ! Lagrange cubic interp wghts (vert)  real(r8), intent(in)   :: lbassi  (4,2,plevp)   ! Lagrange cubic interp wghts (vert)  real(r8), intent(in)   :: detam   (plev)        ! delta eta at levels  real(r8), intent(in)   :: detai   (plevp)       ! delta eta at interfaces  real(r8), intent(in)   :: dlam    (platd)       ! longitudinal grid increment  real(r8), intent(in)   :: cwava   (plat)        ! weight for global water vapor int.  real(r8), intent(in)   :: etamid  (plev)        ! eta at levels  real(r8), intent(in)   :: etaint  (plevp)       ! eta at interfaces  real(r8), intent(in)   :: grfu    (plond,plev,beglat:endlat) ! nonlinear term - u momentum eqn  real(r8), intent(in)   :: grfv    (plond,plev,beglat:endlat) ! nonlinear term - v momentum eqn#if ( defined PVP )  real(r8), intent(out)  :: grlps1(2*pmmax     ,plat/2) ! ------------------------------  real(r8), intent(out)  :: grlps2(2*pmmax     ,plat/2) ! |  real(r8), intent(out)  :: grt1  (2*pmmax,plev,plat/2) ! |  real(r8), intent(out)  :: grt2  (2*pmmax,plev,plat/2) ! |  real(r8), intent(out)  :: grq1  (2*pmmax,plev,plat/2) ! |- see quad for definitions  real(r8), intent(out)  :: grq2  (2*pmmax,plev,plat/2) ! |  real(r8), intent(out)  :: grfu1 (2*pmmax,plev,plat/2) ! |  real(r8), intent(out)  :: grfu2 (2*pmmax,plev,plat/2) ! |  real(r8), intent(out)  :: grfv1 (2*pmmax,plev,plat/2) ! |  real(r8), intent(out)  :: grfv2 (2*pmmax,plev,plat/2) ! ------------------------------#else  real(r8), intent(out)  :: grlps1(     2*pmmax,plat/2) ! ------------------------------  real(r8), intent(out)  :: grlps2(     2*pmmax,plat/2) ! |  real(r8), intent(out)  :: grt1  (plev,2*pmmax,plat/2) ! |  real(r8), intent(out)  :: grt2  (plev,2*pmmax,plat/2) ! |  real(r8), intent(out)  :: grq1  (plev,2*pmmax,plat/2) ! |  real(r8), intent(out)  :: grq2  (plev,2*pmmax,plat/2) ! |- see quad for definitions  real(r8), intent(out)  :: grfu1 (plev,2*pmmax,plat/2) ! |  real(r8), intent(out)  :: grfu2 (plev,2*pmmax,plat/2) ! |  real(r8), intent(out)  :: grfv1 (plev,2*pmmax,plat/2) ! |  real(r8), intent(out)  :: grfv2 (plev,2*pmmax,plat/2) ! ------------------------------#endif  real(r8), intent(in)   :: ps      (plond,beglat:endlat,ptimelevels)  real(r8), intent(in)   :: u3      (plond,plev,beglatex:endlatex,ptimelevels) ! u-wind com  real(r8), intent(in)   :: v3      (plond,plev,beglatex:endlatex,ptimelevels) ! v-wind comp  real(r8), intent(in)   :: t3      (plond,plev,beglatex:endlatex,ptimelevels) ! temperature  real(r8), intent(in)   :: q3      (plond,plev,pcnst+pnats,beglatex:endlatex,ptimelevels)!                                                                    ! q and const  real(r8), intent(in)   :: lnpssld (plond,plev,beglatex:endlatex)   ! RHS Ps term for SLD  real(r8), intent(in)   :: prhssld (plond,plev,beglatex:endlatex)   ! RHS Ps term for SLD  real(r8), intent(in)   :: tarrsld (plond,plev,beglatex:endlatex)   ! T  at arr. pt.(SLD)  real(r8), intent(inout):: parrsld (plond,plev,beglatex:endlatex)   ! Ps at arr. pt.(SLD)  integer , intent(in)   :: n3                                       ! time index  integer , intent(in)   :: n3m1                                     ! time index  real(r8), intent(in)   :: u3sld   (plond,plev ,beglatex:endlatex)  ! u3 inpt for SLD int  real(r8), intent(in)   :: v3sld   (plond,plev ,beglatex:endlatex)  ! v3 inpt for SLD int  real(r8), intent(in)   :: etadot  (plond,plevp,beglatex:endlatex,ptimelevels)! Vertical motion  integer , intent(in)   :: nlon                                     ! # of longitudes!!---------------------------Local workspace-----------------------------!  integer i                    ! index  integer k                    ! index  integer l                    ! index  integer m                    ! constituent index  integer inc                  ! increment for fft991  integer ntr                  ! number of fft's to perform  integer isign                ! flag indicates fft transform directn  integer irow                 ! N/S latitude pair index  integer jcen                 ! lat index (extended grid)!                              ! of forecast  real(r8) fdp  (plon,plev,2)  ! interpolant  real(r8) pmid (plond,plev)   ! pressure at model levels  real(r8) pint (plond,plevp)  ! pressure at interfaces  real(r8) pdel (plond,plev)   ! pressure difference between  real(r8) lamdp(plon,plev)    ! x-coord of dep pt  real(r8) phidp(plon,plev)    ! y-coord of dep pt  real(r8) sigdp(plon,plev)    ! z-coord of dep pt  integer idp   (plon,plev,4)  ! zonal      dep point index  integer jdp   (plon,plev)    ! meridional dep point index  integer kdp   (plon,plev)    ! vertical   dep point index  integer kkdp  (plon,plev)    ! index of z-coordinate of dep pt (alt)  real(r8) xl   (plon,plev,4)  ! weight for x-interpolants (left)  real(r8) xr   (plon,plev,4)  ! weight for x-interpolants (right)  real(r8) wgt1x(plon,plev,4)  ! weight for x-interpolants (Lag Cubic)  real(r8) wgt2x(plon,plev,4)  ! weight for x-interpolants (Lag Cubic)  real(r8) wgt3x(plon,plev,4)  ! weight for x-interpolants (Lag Cubic)  real(r8) wgt4x(plon,plev,4)  ! weight for x-interpolants (Lag Cubic)  real(r8) hl   (plon,plev,4)  ! weight for x-interpolants (Hermite)  real(r8) hr   (plon,plev,4)  ! weight for x-interpolants (Hermite)  real(r8) dhl  (plon,plev,4)  ! weight for x-interpolants (Hermite)  real(r8) dhr  (plon,plev,4)  ! weight for x-interpolants (Hermite)  real(r8) ys   (plon,plev)    ! weight for y-interpolants (south)  real(r8) yn   (plon,plev)    ! weight for y-interpolants (north)  real(r8) wgt1y(plon,plev)    ! weight for y-interpolants (Lag Cubic)  real(r8) wgt2y(plon,plev)    ! weight for y-interpolants (Lag Cubic)  real(r8) wgt3y(plon,plev)    ! weight for y-interpolants (Lag Cubic)  real(r8) wgt4y(plon,plev)    ! weight for y-interpolants (Lag Cubic)  real(r8) hs   (plon,plev)    ! weight for y-interpolants (Hermite)  real(r8) hn   (plon,plev)    ! weight for y-interpolants (Hermite)  real(r8) dhs  (plon,plev)    ! weight for y-interpolants (Hermite)  real(r8) dhn  (plon,plev)    ! weight for y-interpolants (Hermite)  real(r8) rdphi(plon,plev)    ! reciprocal of y-interval  real(r8) wgt1z(plon,plev)    ! weight for z-interpolants (Lag Cubic)  real(r8) wgt2z(plon,plev)    ! weight for z-interpolants (Lag Cubic)  real(r8) wgt3z(plon,plev)    ! weight for z-interpolants (Lag Cubic)  real(r8) wgt4z(plon,plev)    ! weight for z-interpolants (Lag Cubic)  real(r8) hb   (plon,plev)    ! weight for z-interpolants (Hermite)  real(r8) ht   (plon,plev)    ! weight for z-interpolants (Hermite)  real(r8) dhb  (plon,plev)    ! weight for z-interpolants (Hermite)  real(r8) dht  (plon,plev)    ! weight for z-interpolants (Hermite)  real(r8) rdz  (plon,plev)    ! reciprocal of z-interval  real(r8) zt   (plon,plev)    ! top vertical interpolation weight   real(r8) zb   (plon,plev)    ! bot vertical interpolation weight   real(r8) lampr(plon,plev)    ! trajectory increment (x-direction)  real(r8) phipr(plon,plev)    ! trajectory increment (y-direction)  real(r8) upr  (plon,plev)    ! interpolated u field (local geodesic)  real(r8) vpr  (plon,plev)    ! interpolated v field (local geodesic)#if ( ! defined USEFFTLIB )  real(r8) work((plon+1)*5*plev)   ! workspace array for fft991#else  real(r8) work((plon+1)*pcray)    ! workspace array for fft991#endif  real(r8) xnlin(plndlv*4 + plond) ! non-linear terms (equivalence !                                  ! region for following arrays to !                                  ! optimize fft performance)  real(r8) grfulat(plond,plev)     ! non-linear terms for u-momentum   real(r8) grfvlat(plond,plev)     ! non-linear terms for u-momentum   real(r8) grtlat (plond,plev)     ! RHS of T-eqn  real(r8) grqlat (plond,plev)     ! q  real(r8) grpslat(plond)          ! RHS of Ps-eqn!! The following equivalences are for optimal fft performance!#if ( defined SUNOS )!! Stupid Sun compiler hoop-jumping again!  save grfulat, grfvlat, grtlat, grqlat, grpslat, xnlin#endif  equivalence   (grfulat,xnlin(1))  equivalence   (grfvlat,xnlin(1+1*plndlv))  equivalence   (grtlat ,xnlin(1+2*plndlv))  equivalence   (grqlat ,xnlin(1+3*plndlv))  equivalence   (grpslat,xnlin(1+4*plndlv))!  real(r8) pd   (plond)            ! RHS term for Ps and (1/ps)etadot(dp/deta)   real(r8) pdsum(plond)            ! RHS term for Ps and (1/ps)etadot(dp/deta)   real(r8) pd1  (plond)            ! RHS term for Ps and (1/ps)etadot(dp/deta)   real(r8) pdsm1(plond)            ! RHS term for Ps and (1/ps)etadot(dp/deta)   real(r8) pa   (plond)            ! RHS term for Ps and (1/ps)etadot(dp/deta)   real(r8) pasum(plond)            ! RHS term for Ps and (1/ps)etadot(dp/deta)   real(r8) coslat                  ! cos(latitude)  real(r8) tmp1                    ! temp space!  logical limdrh                   ! horizontal derivative limiter flag  logical limdrv                   ! vertical   derivative limiter flag  logical lhrzint                  ! horizontal interp flag  logical lvrtint                  ! vertical   interp flag  logical lhrzwgt                  ! flag to compute horizontal weights  logical lvrtwgt                  ! flag to compute vertical   weights!!-----------------------------------------------------------------------!  if(lat.le.plat/2) then     irow = lat  else     irow = plat + 1 - lat  end if  jcen = j1 - 1 + lat  coslat = cos(clat(lat))!! Initial guess for trajectory midpoints in spherical coords.! Use arrival points as initial guess for trajectory midpoints.!  do k=1,plev     do i=1,nlon        phidp(i,k) = clat(lat)        sigdp(i,k) = etamid(k)     end do  end do!        ! Offset bottom level departure point first guess by epsilon!  do i = 1,nlon     sigdp(i,plev) = sigdp(i,plev)*(1. - 10.*epsilon(sigdp))  end do!! Loop through latitudes producing departure point calculation!  call slttraj(pmap    ,jcen    ,lat     ,ztodt   ,ra      , &               iter    ,lam     ,phi     ,dphi    ,etamid  , &               etaint  ,detam   ,detai   ,lbasiy  ,lbasiz  , &               lbassi  ,kdpmpf  ,kdpmph  ,idp     ,jdp     , &               kdp     ,kkdp    ,xl      ,xr      ,wgt1x   , &               wgt2x   ,wgt3x   ,wgt4x   ,hl      ,hr      , &               dhl     ,dhr     ,ys      ,yn      ,wgt1y   , &               wgt2y   ,wgt3y   ,wgt4y   ,hs      ,hn      , &               dhs     ,dhn     ,rdphi   ,wgt1z   ,wgt2z   , &               wgt3z   ,wgt4z   ,hb      ,ht      ,dhb     , &               dht     ,rdz     ,lampr   ,phipr   ,upr     , &               vpr     ,lamdp   ,phidp   ,sigdp   ,u3      , &               v3      ,u3sld   ,v3sld   ,etadot  ,n3      , &               n3m1    ,dlam    ,nlon    )!! Calculate mass of moisture in field being advected by slt.!  call plevs0(nlon    ,plond   ,plev    ,ps(1,lat,n3),pint    ,pmid    ,pdel)  call qmassa(cwava(lat)  ,w(irow) ,q3(i1,1,1,jcen,n3),pdel    ,hw1lat(1,lat), &              nlon)!! Compute constituent forecast!  lhrzwgt = .true.  lvrtwgt = .true.  lhrzint = .true.  lvrtint = .true.  limdrh  = .true.  limdrv  = .true.  call bandij (dlam    ,phi     ,lamdp   ,phidp   ,idp     , &               jdp     ,nlon    )  call kdpfnd (plev    ,pmap    ,etamid  ,sigdp   ,kdpmpf  , &               kdp     ,nlon    )  call sltwgts(limdrh  ,limdrv  ,lhrzwgt ,lvrtwgt ,plev    , &               idp     ,jdp     ,kdp     ,lam     ,phi     , &               etamid  ,dphi    ,detam   ,lamdp   ,phidp   , &               sigdp   ,lbasiy  ,lbasiz  ,kkdp    ,xl      , &               xr      ,wgt1x   ,wgt2x   ,wgt3x   ,wgt4x   , &               hl      ,hr      ,dhl     ,dhr     ,ys      , &               yn      ,wgt1y   ,wgt2y   ,wgt3y   ,wgt4y   , &               hs      ,hn      ,dhs     ,dhn     ,rdphi   , &               wgt1z   ,wgt2z   ,wgt3z   ,wgt4z   ,hb      , &               ht      ,dhb     ,dht     ,rdz     ,zt      , &