slttraj.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
#include <params.h>
subroutine 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    )
!
!-----------------------------------------------------------------------
!
! Purpose:
! Determine trajectory departure point for each arrival point.
! Assume that the first guess for the dep point coordinates are the
! arrival point coordinates
!
! Author:  J. Olson
!
!-----------------------------------------------------------------------
!
! $Id: slttraj.F90,v 1.4.8.1 2002/04/22 19:09:52 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

  use precision
  use pmgrid
  use comslt
  use prognostics, only: ptimelevels
#if (!defined CRAY)
  use srchutil
#endif
  implicit none

!------------------------------Arguments--------------------------------
!
  integer , intent(in)   :: pmap                ! artificial vert grid dim.
  integer , intent(in)   :: jcen                ! index of lat slice(extend)
  integer , intent(in)   :: lat                 ! index of lat slice (model)
  real(r8), intent(in)   :: ztodt               ! time step (seconds)
  real(r8), intent(in)   :: ra                  ! 1./(radius of earth)
  integer , intent(in)   :: iter                ! iteration count
  real(r8), intent(in)   :: lam   (plond,platd) ! long. coord of model grid
  real(r8), intent(in)   :: phi   (platd)       ! lat.  coord of model grid
  real(r8), intent(in)   :: dphi  (platd)       ! increment between lats.
  real(r8), intent(in)   :: etamid(plev)        ! vertical full levels
  real(r8), intent(in)   :: etaint(plevp)       ! vertical half levels
  real(r8), intent(in)   :: detam (plev)        ! increment between full levs
  real(r8), intent(in)   :: detai (plevp)       ! increment between half levs
  real(r8), intent(in)   :: lbasiy(4,2,platd)   ! lat interp wts(lagrng)
  real(r8), intent(in)   :: lbasiz(4,2,plev)    ! vert interp wghts (lagrng)
  real(r8), intent(in)   :: lbassi(4,2,plevp)   ! vert interp wghts (lagrng)

  integer , intent(in)   :: kdpmpf(pmap)        ! artificial vert grid indices
  integer , intent(in)   :: kdpmph(pmap)        ! artificial vert grid indices
  integer , intent(out)  :: idp   (plon,plev,4) ! zonal      dep point index
  integer , intent(out)  :: jdp   (plon,plev)   ! meridional dep point index
  integer , intent(out)  :: kdp   (plon,plev)   ! vertical   dep point index
  integer , intent(in)   :: kkdp  (plon,plev)   ! index of z-coordinate of dep pt (alt)

  real(r8), intent(out)  :: xl    (plon,plev,4) ! weight for x-interpolants (left)
  real(r8), intent(out)  :: xr    (plon,plev,4) ! weight for x-interpolants (right)
  real(r8), intent(in)   :: wgt1x (plon,plev,4) ! weight for x-interpolants (Lag Cubic)
  real(r8), intent(in)   :: wgt2x (plon,plev,4) ! weight for x-interpolants (Lag Cubic)
  real(r8), intent(in)   :: wgt3x (plon,plev,4) ! weight for x-interpolants (Lag Cubic)
  real(r8), intent(in)   :: wgt4x (plon,plev,4) ! weight for x-interpolants (Lag Cubic)
  real(r8), intent(in)   :: hl    (plon,plev,4) ! weight for x-interpolants (Hermite)
  real(r8), intent(in)   :: hr    (plon,plev,4) ! weight for x-interpolants (Hermite)
  real(r8), intent(in)   :: dhl   (plon,plev,4) ! weight for x-interpolants (Hermite)
  real(r8), intent(in)   :: dhr   (plon,plev,4) ! weight for x-interpolants (Hermite)

  real(r8), intent(out)  :: ys    (plon,plev)   ! weight for y-interpolants (south)
  real(r8), intent(out)  :: yn    (plon,plev)   ! weight for y-interpolants (north)
  real(r8), intent(in)   :: wgt1y (plon,plev)   ! weight for y-interpolants (Lag Cubic)
  real(r8), intent(in)   :: wgt2y (plon,plev)   ! weight for y-interpolants (Lag Cubic)
  real(r8), intent(in)   :: wgt3y (plon,plev)   ! weight for y-interpolants (Lag Cubic)
  real(r8), intent(in)   :: wgt4y (plon,plev)   ! weight for y-interpolants (Lag Cubic)
  real(r8), intent(in)   :: hs    (plon,plev)   ! weight for y-interpolants (Hermite)
  real(r8), intent(in)   :: hn    (plon,plev)   ! weight for y-interpolants (Hermite)
  real(r8), intent(in)   :: dhs   (plon,plev)   ! weight for y-interpolants (Hermite)
  real(r8), intent(in)   :: dhn   (plon,plev)   ! weight for y-interpolants (Hermite)
  real(r8), intent(in)   :: rdphi (plon,plev)   ! reciprocal of y-interval

  real(r8), intent(in)   :: wgt1z (plon,plev)   ! weight for z-interpolants (Lag Cubic)
  real(r8), intent(in)   :: wgt2z (plon,plev)   ! weight for z-interpolants (Lag Cubic)
  real(r8), intent(in)   :: wgt3z (plon,plev)   ! weight for z-interpolants (Lag Cubic)
  real(r8), intent(in)   :: wgt4z (plon,plev)   ! weight for z-interpolants (Lag Cubic)
  real(r8), intent(in)   :: hb    (plon,plev)   ! weight for z-interpolants (Hermite)
  real(r8), intent(in)   :: ht    (plon,plev)   ! weight for z-interpolants (Hermite)
  real(r8), intent(in)   :: dhb   (plon,plev)   ! weight for z-interpolants (Hermite)
  real(r8), intent(in)   :: dht   (plon,plev)   ! weight for z-interpolants (Hermite)
  real(r8), intent(in)   :: rdz   (plon,plev)   ! reciprocal of z-interval

  real(r8), intent(out)  :: lampr (plon,plev)   ! trajectory increment (x-direction)
  real(r8), intent(out)  :: phipr (plon,plev)   ! trajectory increment (y-direction)
  real(r8), intent(out)  :: upr   (plon,plev)   ! interpolated u field (local geodesic)
  real(r8), intent(out)  :: vpr   (plon,plev)   ! interpolated v field (local geodesic)
  real(r8), intent(out)  :: lamdp (plon,plev)   ! zonal      departure pt. coord.
  real(r8), intent(inout):: phidp (plon,plev)   ! meridional departure pt. coord.
  real(r8), intent(inout):: sigdp (plon,plev)   ! vertical   departure pt. coord.
  real(r8), intent(in)   :: u3    (plond, plev, beglatex:endlatex,ptimelevels) ! u wind component
  real(r8), intent(in)   :: v3    (plond, plev, beglatex:endlatex,ptimelevels) ! v wind component
  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)   :: n3                                       ! time indicies
  integer , intent(in)   :: n3m1                                     ! time indicies
  real(r8), intent(in)   :: dlam  (platd)                            ! d-long for each lat
  integer , intent(in)   :: nlon                                     ! # of longitudes
#if (defined CRAY)
  integer,external :: ismin
#endif
!
!---------------------------Local workspace-----------------------------
!
  real(r8) fac                  ! 1-eps
  real(r8) tmp                  ! temp variable
  real(r8) ump    (plon,plev)   ! interpolated u field
  real(r8) vmp    (plon,plev)   ! interpolated v field
  real(r8) wmp    (plon,plev)   ! interpolated w field
  real(r8) wmpa   (plon,plev)   ! interpolated w field (arrival level)
  real(r8) sigpr  (plon,plev)   ! trajectory increment (vertical)
  real(r8) etamin               ! minimum eta
  real(r8) etamax               ! maximum eta
  real(r8) delsig (plev)        ! dist between dep point and model levs
  real(r8) zt     (plon,plev)   ! top vertical interpolation weight 
  real(r8) zb     (plon,plev)   ! bot vertical interpolation weight 
  real(r8) zttmp                ! top vertical interpolation weight 
  real(r8) zbtmp                ! bot vertical interpolation weight 
  real(r8) rdx    (platd)       ! reciprocal of del-x

  logical locgeo                ! local geodesic flag
  logical lvert                 ! flag to compute vertical trajectory
  logical larrival(plon)        ! flag to indicate whether or not to
!                               ! place the alternative dep point at the
!                               ! arrival point.

  integer kount                 ! index counter
  integer kkmin                 ! index of minimum "delsig"
  integer i                     ! |
  integer ii                    ! |
  integer j                     ! |
  integer jj2                   ! | - indices
  integer jj3                   ! |
  integer k                     ! |
  integer kk                    ! |
  integer n                     ! |
!
!-----------------------------------------------------------------------
!
  fac     = 1. - 10.*epsilon (fac)
  locgeo  = .false.
  if(abs(phi(jcen)) .ge. phigs) locgeo = .true.
!
! Interpolate arrival and departure vertical velocities
!
  do k = 1,plev
     do i = 1,nlon
        zttmp = ( etaint(k+1) - sigdp(i,k) )/detai(k)
        zbtmp = 1. - zttmp
        wmpa(i,k) = etadot(i1+i-1,k  ,jcen,n3  )*zttmp &
                  + etadot(i1+i-1,k+1,jcen,n3  )*zbtmp
        wmp (i,k) = etadot(i1+i-1,k  ,jcen,n3m1)*zttmp &
                  + etadot(i1+i-1,k+1,jcen,n3m1)*zbtmp
     end do
  end do
!
! Set up computation of trajectory
!
  do k = 1,plev
     do i = 1,nlon
        ii = i + i1 - 1
!
! Place u/v on unit sphere
!
        ump(i,k) = 0.5*(u3sld(ii,k,jcen) + u3(ii,k,jcen,n3))*ra
        vmp(i,k) = 0.5*(v3sld(ii,k,jcen) + v3(ii,k,jcen,n3))*ra
        wmp(i,k) = 0.5*(wmp(i,k)         + wmpa(i,k)       )
!
! Estimate departure point of parcel trajectory.
!
        lampr(i,k) = -ztodt*ump(i,k)
        phipr(i,k) = -ztodt*vmp(i,k)
        sigpr(i,k) = -ztodt*wmp(i,k)
        if (.not. locgeo) lampr(i,k) = lampr(i,k)/cos( phidp(i,k) )
!
! Initialize winds for use in g.c. calculations.
!
        if (locgeo) then
           upr(i,k) = ump(i,k)
           vpr(i,k) = vmp(i,k)
        end if
     end do
  end do
!
! Estimate initial trajectory departure points
!
  lvert = .true.
  call trajdp(lat     ,jcen    ,ztodt   ,locgeo  ,lvert   , &
              lam(1,jcen),phi  ,etamid  ,upr     ,vpr     , &
              lampr   ,phipr   ,sigpr   ,lamdp   ,phidp   , &
              sigdp   ,nlon    )
!
! Loop over departure point iterates.
!
  do n = 1,iter
!
! Determine departure point indicies and interpolate u,v,w
!
     call bandij (dlam    ,phi     ,lamdp   ,phidp   ,idp     , &
                  jdp     ,nlon    )
     call kdpfnd (plev    ,pmap    ,etamid  ,sigdp   ,kdpmpf  , &
                  kdp     ,nlon    )
!
! Compute weights for x,y,z dimensions
!
     do j = 1,platd
        rdx(j) = 1./(lam(nxpt+2,j) - lam(nxpt+1,j))
     end do
     do k = 1,plev
        do i = 1,nlon
           jj2 = jdp(i,k)
           jj3 = jdp(i,k) + 1
           xl(i,k,2) = (lam(idp(i,k,2)+1,jj2) - lamdp(i,k))*rdx(jj2)
           xl(i,k,3) = (lam(idp(i,k,3)+1,jj3) - lamdp(i,k))*rdx(jj3)
           xr(i,k,2) = 1. - xl(i,k,2)
           xr(i,k,3) = 1. - xl(i,k,3)
           ys(i,k)   = ( phi(jdp(i,k)+1) - phidp(i,k) )/dphi(jdp(i,k))
           yn(i,k)   = 1. - ys(i,k)
           zt(i,k)   = (etamid(kdp(i,k)+1)-sigdp(i,k))/detam(kdp(i,k))
           zb(i,k)   = 1. - zt(i,k)
        end do
     end do
!
     call sltlinint(plev    ,1       ,u3sld(1,1,beglatex),xl      ,xr      , &
                    ys      ,yn      ,zt                 ,zb      ,idp     , &
                    jdp     ,kdp     ,ump                ,nlon    )
     call sltlinint(plev    ,1       ,v3sld(1,1,beglatex),xl      ,xr      , &
                    ys      ,yn      ,zt                 ,zb      ,idp     , &
                    jdp     ,kdp     ,vmp                ,nlon    )
     call kdpfnd (plevp   ,pmap    ,etaint  ,sigdp   ,kdpmph  , &
                  kdp     ,nlon    )
     do k = 1,plev
        do i = 1,nlon
           zt(i,k)   = (etaint(kdp(i,k)+1)-sigdp(i,k))/detai(kdp(i,k))
           zb(i,k)   = 1. - zt(i,k)
        end do
     end do
     call sltlinint(plevp   ,1       ,etadot(1,1,beglatex,n3m1),xl      ,xr      , &
                    ys      ,yn      ,zt                       ,zb      ,idp     , &
                    jdp     ,kdp     ,wmp                      ,nlon    )
!
! Add arrival point velocities (n) to interpolated velocities
!
     do k = 1,plev
        do i = 1,nlon
           ii = i + i1 - 1
           ump(i,k) = 0.5*(ump(i,k) + u3(ii,k,jcen,n3))
           vmp(i,k) = 0.5*(vmp(i,k) + v3(ii,k,jcen,n3))
           wmp(i,k) = 0.5*(wmp(i,k) + wmpa(i,k)       )
        end do
     end do
!
! Compute new trajectory departure points
!
     lvert = .true.
     call depinc (jcen    ,ztodt   ,ra      ,locgeo  ,lvert   , &
                  lam(1,jcen),phi  ,ump     ,vmp     ,wmp     , &
                  upr     ,vpr     ,lamdp   ,phidp   ,lampr   , &
                  phipr   ,sigpr   ,nlon    )
     call trajdp (lat     ,jcen    ,ztodt   ,locgeo  ,lvert   , &
                  lam(1,jcen),phi  ,etamid  ,upr     ,vpr     , &
                  lampr   ,phipr   ,sigpr   ,lamdp   ,phidp   , &
                  sigdp   ,nlon    )
  end do
!
! Compile departure point binning statistics
!
  do k = 1,plev
     if(k .eq. 1) then
        etamin = etamid(k  )
        etamax = etamid(k+1)
     elseif(k .eq. plev) then
        etamin = etamid(k-1)
        etamax = etamid(k  )
     else
        etamin = etamid(k-1)
        etamax = etamid(k+1)
     end if
     do i = 1,nlon
        larrival(i) = sigdp(i,k) .ge. etamin .and. sigdp(i,k) .le. etamax
     end do
!
! If:  Departure point is within one grid interval from the arrival
! point, bin it in the ARRIVAL point bin
!
     tmp = 0.
     do i = 1,nlon
        if (larrival(i)) then
           tmp           = tmp + 1.
        end if
     end do
     levkntl(k,k,lat) = levkntl(k,k,lat) + tmp
!
! Else:  find departure point bin
!
     kount = 0
     do i=1,nlon
        if (larrival(i)) kount = kount + 1
     end do
     if (kount .ne. nlon) then
        do i = 1,nlon
           if (.not. larrival(i)) then
              delsig(1) = 1.e+35
              do kk = 2,plev
                 delsig(kk) = abs( sigdp(i,k) - etaint(kk) )
              end do
              kkmin = ismin( plev,delsig,1 )
              if(kkmin .gt. k) kkmin = kkmin-1
              levkntl(kkmin,k,lat) = levkntl(kkmin,k,lat) + 1.
           end if
        end do
     end if
  end do
!
  return
end subroutine slttraj