#include <misc.h>#include <params.h>subroutine depinc(jcen    ,dt      ,ra      ,locgeo  ,lvert   , &                  lam     ,phi     ,ump     ,vmp     ,wmp     , &                  upr     ,vpr     ,lamdp   ,phidp   ,lampr   , &                  phipr   ,sigpr   ,nlon    )!!-----------------------------------------------------------------------!! Purpose:! Compute departure point increments based upon interpolated wind components.!! Author:  J. Olson!!-----------------------------------------------------------------------!! $Id: depinc.F90,v 1.2.42.1 2002/04/22 19:09:49 erik Exp $! $Author: erik $!!-----------------------------------------------------------------------  use precision  use pmgrid  implicit none!------------------------------Arguments--------------------------------!  integer , intent(in)   :: jcen              ! index of lat slice(extend)  real(r8), intent(in)   :: dt                ! time step (seconds)  real(r8), intent(in)   :: ra                ! 1./(radius of earth)  logical , intent(in)   :: locgeo            ! local geodesic flag  logical , intent(in)   :: lvert             ! flag to compute vertical trajectory  real(r8), intent(in)   :: lam  (plond)      ! long. coord of model grid  real(r8), intent(in)   :: phi  (platd)      ! lat.  coord of model grid  real(r8), intent(inout):: ump  (plon,plev)  ! interpolated u field  real(r8), intent(inout):: vmp  (plon,plev)  ! interpolated v field  real(r8), intent(in)   :: wmp  (plon,plev)  ! interpolated w field  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(in)   :: lamdp(plon,plev)  ! zonal      dep pt. coord.  real(r8), intent(in)   :: phidp(plon,plev)  ! meridional dep pt. coord.  real(r8), intent(out)  :: lampr(plon,plev)  ! trajectory increment (x-direction)  real(r8), intent(inout):: phipr(plon,plev)  ! trajectory increment (y-direction)  real(r8), intent(out)  :: sigpr(plon,plev)  ! trajectory increment (vertical)  integer , intent(in)   :: nlon              ! number of longitudes for this lat!!---------------------------Local workspace-----------------------------!  real(r8) phi0    ! Current latitude (radians)             real(r8) cphi0   ! cosine of latitude  real(r8) sphi0   ! sine   of latitude  real(r8) cdlam   ! |  real(r8) clamgc  ! |  real(r8) cphid   ! |  real(r8) cphigc  ! |  real(r8) dlamgc  ! | -- temporary variables  real(r8) sdlam   ! |  real(r8) slamgc  ! |  real(r8) sphid   ! |  real(r8) sphigc  ! |  integer i        ! |  integer ii       ! | -- indices  integer k        ! |!!-----------------------------------------------------------------------!  phi0    = phi(jcen)  cphi0   = cos(phi0)  sphi0   = sin(phi0)!! Compute trajectory increment!  do k = 1,plev!! Place u/v on unit sphere!     do i = 1,nlon        ump(i,k) = ump(i,k)*ra        vmp(i,k) = vmp(i,k)*ra        if (lvert) sigpr(i,k) = -dt*wmp(i,k)     end do!! If near the pole, convert to g.c.!     if (locgeo) then!! Compute winds at departure points in g.c.!        do i = 1,nlon           ii = i + i1 - 1           dlamgc = lam(ii) - lamdp(i,k)           sdlam  = sin( dlamgc )           cdlam  = cos( dlamgc )           sphid  = sin( phidp(i,k) )           cphid  = cos( phidp(i,k) )           sphigc = sphid*cphi0 - cphid*sphi0*cdlam           cphigc = cos( asin( sphigc ) )           slamgc = -sdlam*cphid/cphigc           clamgc = cos( asin( slamgc ) )           vpr(i,k) = (vmp(i,k)*(cphid*cphi0 + sphid*sphi0*cdlam) &                - ump(i,k)*sphi0*sdlam)/cphigc           upr(i,k) = (ump(i,k)*cdlam + vmp(i,k)*sphid*sdlam + &                vpr(i,k)*slamgc*sphigc)/clamgc!! Compute trajectory increment between arrival and departure points.!           lampr(i,k) = -dt*upr(i,k)/cos(phipr(i,k))           phipr(i,k) = -dt*vpr(i,k)        end do     else        do i = 1,nlon           lampr(i,k) = -dt*ump(i,k)/cos(phidp(i,k))           phipr(i,k) = -dt*vmp(i,k)        end do     end if  end do!  returnend subroutine depinc