#include <misc.h>
#include <params.h>

subroutine linemsbc(lat     ,tdt     ,ps      ,psm1    ,&
                    u3      ,u3m1    ,v3      ,v3m1    ,t3      ,&
                    t3m1    ,q3      ,q3m1    ,&
                    dpsl    ,phis    ,&
                    asdir   ,asdif   ,aldir   ,aldif   ,dpsm    ,&
                    oro     ,ts      ,snowh   ,srfrad  ,&
                    precc   ,precl   ,precsc  ,precsl  ,sols    ,&
                    soll    ,solsd   ,solld   ,flwds   ,lwup    ,&
                    vort    ,div     ,&
                    vorttl  ,divtl   ,psm2    ,um2     ,vm2     ,&
                    tm2     ,qm2     ,pblht   ,qpert   ,tpert   ,&
                    fsns    ,qrs     ,qrl     ,         &
                    vortm1  ,divm1   ,omga    ,psomc   ,thm     ,&
                    zm      ,pmln    ,pilnm1  ,sgh     ,cnt     ,&
                    cnb     ,cmfmc   )

!-----------------------------------------------------------------------
!
! Purpose:
! Control adjustment physics, time filter, history tape write-up,
! tendency physics, non-linear dynamical terms, FFT and combine terms
! in preparation for Fourier -> spectral quadrature.
!
! Method:
! The naming convention is as follows:
!  - prefix gr contains grid point values before FFT and Fourier
!     coefficients after
!  - t, q, d, z and ps refer to temperature, specific humidity,
!     divergence, vorticity and surface pressure
!  - "1" suffix to an array => symmetric component current latitude pair
!  - "2" suffix to an array => antisymmetric component.
!
! Author:
! Original version:  CCM1
! Standardized:      J. Rosinski, June 1992
! Reviewed:          D. Williamson, J. Hack, August 1992
! Reviewed:          D. Williamson, March 1996
!
!-----------------------------------------------------------------------
!
! $Id: linemsbc.F90,v 1.1.6.2 2000/09/07 16:44:54 rosinski Exp $
! $Author: rosinski $

   use precision
   use history
   use pmgrid
   use pspect
!-----------------------------------------------------------------------
   implicit none
!------------------------------Commons----------------------------------
#include <comcon.h>
!-----------------------------------------------------------------------
#include <comctl.h>
!-----------------------------------------------------------------------
#include <comindx.h>
!-----------------------------------------------------------------------
#include <comlun.h>
!-----------------------------------------------------------------------
#include <commap.h>
!-----------------------------------------------------------------------
#include <comqfl.h>
!-----------------------------------------------------------------------
#include <comtfc.h>
!-----------------------------------------------------------------------
#include <comtim.h>
!-----------------------------------------------------------------------
#include <comtrcnm.h>
!-----------------------------------------------------------------------
!
! Input arguments
!
   integer, intent(in) :: lat                 ! latitude index for S->N storage

   real(r8), intent(in) :: tdt                  ! two delta t (unless nstep<2)
   real(r8), intent(in) :: dpsl(plond)
   real(r8), intent(in) :: phis(plond)
   real(r8), intent(in) :: asdir(plond)           ! albedo short wave direct
   real(r8), intent(in) :: asdif(plond)           ! albedo short wave diffuse
   real(r8), intent(in) :: aldir(plond)           ! albedo long wave direct
   real(r8), intent(in) :: aldif(plond)           ! albedo long wave diffuse
   real(r8), intent(in) :: dpsm(plond)           ! albedo long wave diffuse
   real(r8), intent(in) :: oro(plond)             ! land/ocean/sea ice flag
   real(r8), intent(in) :: snowh(plond)           ! equivalent liquid water snow depth
   real(r8), intent(inout) :: srfrad(plond)
   real(r8), intent(inout) :: precc(plond)           ! convective precip rate
   real(r8), intent(inout) :: precl(plond)           ! large scale stable precip rate
   real(r8), intent(inout) :: precsl(plond)          ! snow component large scale precip
   real(r8), intent(inout) :: precsc(plond)          ! snow component convective precip
   real(r8), intent(inout) :: sols(plond)            ! direct beam solar rad onto srf(sw)
   real(r8), intent(inout) :: soll(plond)            ! direct beam solar rad onto srf(lw)
   real(r8), intent(inout) :: solsd(plond)           ! diffuse solar radiation onto srf(sw)
   real(r8), intent(inout) :: solld(plond)           ! diffuse solar radiation onto srf(lw)
   real(r8), intent(inout) :: flwds(plond)           ! surface longwave down flux
   real(r8), intent(in) :: lwup(plond)            ! surface longwave up flux
   real(r8), intent(in) :: vort(plond,plev)

! Input/output variables

   real(r8), intent(inout) :: ts(plond)              ! surface temperature
   real(r8), intent(inout) :: ps(plond)              ! surface pressure (time n)
   real(r8), intent(inout) :: psm1(plond)            ! surface pressure (time n-1)
   real(r8), intent(inout) :: u3(plond,plev)         ! u-wind (time n)
   real(r8), intent(inout) :: u3m1(plond,plev)       ! u-wind (time n-1)
   real(r8), intent(inout) :: v3(plond,plev)         ! v-wind (time n)
   real(r8), intent(inout) :: v3m1(plond,plev)       ! v-wind (time n-1)
   real(r8), intent(inout) :: t3(plond,plev)         ! temperature (time n)
   real(r8), intent(inout) :: t3m1(plond,plev)       ! temperature (time n-1)
   real(r8), intent(inout) :: q3(plond,plev,pcnst+pnats)   ! constituent conc(time n: h2o first)
   real(r8), intent(inout) :: q3m1(plond,plev,pcnst+pnats) ! constituent conc(time n-1:h2o first)
   real(r8) cnt(plond)             ! top of convective layer
   real(r8) cnb(plond)             ! bottom of convective layer
   real(r8) cmfmc(plond,plev)      ! convective mass flux--m sub c
!
!---------------------------Local workspace-----------------------------
!
   real(r8) qc(plond,plev)         ! q tendency due to total rainout
   real(r8) pmid(plond,plev)       ! pressure at model levels (time n)
   real(r8) pmidm1(plond,plev)     ! pressure at model levels (time n-1)
   real(r8) rpmid(plond,plev)      ! 1./pmid
   real(r8) pint(plond,plevp)      ! pressure at model interfaces (n  )
   real(r8) pintm1(plond,plevp)    ! pressure at model interfaces (n-1)
   real(r8) piln(plond,plevp)      ! log(pint)
   real(r8) pdel(plond,plev)       ! pdel(k)   = pint  (k+1)-pint  (k)
   real(r8) pdelm1(plond,plev)     ! pdelm1(k) = pintm1(k+1)-pintm1(k)
   real(r8) rpdel(plond,plev)      ! 1./pdel
   real(r8) div(plond,plev)
   real(r8) vorttl(plond,plev)
   real(r8) divtl(plond,plev)
   real(r8) psm2(plond)
   real(r8) um2(plond,plev)
   real(r8) vm2(plond,plev)
   real(r8) tm2(plond,plev)
   real(r8) qm2(plond,plev,pcnst+pnats)
   real(r8) pblht(plond)
   real(r8) qpert(plond,pcnst+pnats)
   real(r8) tpert(plond)
   real(r8) fsns(plond)
   real(r8) qrs(plond,plev)
   real(r8) qrl(plond,plev)
   real(r8) vortm1(plond,plev)
   real(r8) divm1(plond,plev)
   real(r8) omga(plond,plev)
   real(r8) psomc(plond,plev)
   real(r8) thm(plond,plev)
   real(r8) zm(plond,plev)
   real(r8) pmln(plond,plev)
   real(r8) pilnm1(plond,plevp)
   real(r8) sgh(plond)             ! std dev of orography of outfld

   real(r8) dtcond(plond,plev)     ! dT/dt due to convective adjustment
   real(r8) dqcond(plond,plev,pcnst+pnats)! dq/dt due to convective adjustment
   real(r8) cmfdt(plond,plev)      ! hack convective temperature tendency
   real(r8) cmfdq(plond,plev)      ! hack convective moisture tendency
   real(r8) zmdt(plond,plev)       ! zhang convective temp tendency
   real(r8) zmdq(plond,plev)       ! zhang convective moisture tendency
   real(r8) cmfdqr(plond,plev)     ! dq/dt due to convective rainout
   real(r8) cmfsl(plond,plev)      ! liquid water static energy flux
   real(r8) cmflq(plond,plev)      ! total water flux
   real(r8) pmlnm1(plond,plev)     ! log(pmidm1)
!
   real(r8) om2eps                 ! 1. - 2.*eps for time filter
!
   integer i,k                     ! longitude,level,constituent indices
!
   real(r8) coslat                 ! cosine(latitude)
   real(r8) rcoslat                ! 1./cosine(latitude)
   real(r8) hr(plond,plev)            ! heating rate for d(theta)/dt calc.
   real(r8) precst(plond)          ! total snow precip
   real(r8) tvm(plond,plev)        ! virtual temperature
!
!-----------------------------------------------------------------------
!
   coslat = cos(clat(lat))
   rcoslat = 1./coslat

   call outfld('SGH     ',sgh,plond   ,lat     )
!
! Set average dry mass to specified constant preserving horizontal
! gradients of ln(ps). Proportionality factor was calculated in STEPON
! for nstep=0.
! Mass fixer only needs to be applied at nstep=0 since its applied in
! timefilter.F for other timesteps
!
   if (nstep == 0) then
      do i=1,plon
         ps(i) = ps(i)*fixmas
      end do
   end if
!
! Set current time pressure arrays for model levels etc.
!
   call plevs0(ps      ,pint    ,pmid    ,pdel    )
!
   do k=1,plev
      do i=1,plon
         pmln(i,k) = log(pmid(i,k))
         rpmid(i,k) = 1./pmid(i,k)
         rpdel(i,k) = 1./pdel(i,k)
         piln(i,k) = log(pint(i,k))
      end do
   end do
   do i=1,plon
      piln(i,plevp) = log(pint(i,plevp))
   end do

   if (nstep.eq.0) then
!
! Copy past time pressure arrays for model levels for first time step
!
      psm1(:plon) = ps(:plon)
!
! Give qm2 meaningful data on startup
!
      qm2(:plon,:plev,:pcnst+pnats) = q3m1(:plon,:plev,:pcnst+pnats)
   endif
!
! Calculate vertical motion field
!
   call omcalc(rcoslat ,div     ,u3      ,v3      ,dpsl    ,  &
               dpsm    ,pmid    ,pdel    ,rpmid   ,pint(1,plevp), &
               omga    )

   if (nstep == 0 .and. .not. adiabatic .and. .not. idealdyn) then
      if (pcnst .gt. 1) call mr2mf( q3 )
!
! Do adjustment physics on initial data (dry convective adjustment,
! convection, and stable precip)
!
      call aphys(nstep   ,lat     ,tdt     ,phis    ,pblht   , &
                 tpert   ,pmid    ,pint    ,pdel    ,rpdel   , &
                 pmln    ,piln    ,t3      ,q3      ,ts      , &
                 qpert   ,precl   ,precc   ,precsl  ,precsc  , &
                 dtcond  ,dqcond  ,cnt     ,cnb     ,qc      , &
                 cmfdt   ,cmfdq   ,cmfdqr  ,cmfmc   ,cmfsl   , &
                 cmflq   ,zmdt    ,zmdq    )
!
! Convert mass fractions of non-water tracers back to mixing ratios.
! (Overwrite non-water portions of q3).
!
      if (pcnst .gt. 1) call mf2mr( q3 )
   end if

   if (nstep.eq.0) then
!
! Copy adjusted time n values to time n-1 first time step only
!
      q3m1(:plon,:plev,:pcnst+pnats) = q3(:plon,:plev,:pcnst+pnats)
      t3m1(:plon,:plev) = t3(:plon,:plev)
   endif
!
! Time filter (second half of filter for vorticity and divergence only)
!
   if (nstep < 3) then
!
! Copy vorticity, divergence when not time filtered
!
      vortm1(:plon,:plev) = vorttl(:plon,:plev)
      divm1(:plon,:plev) = divtl(:plon,:plev)
   end if

   if (nstep.ge.2) then
!
! Time filter (first half for vorticity and divergence only)
!
      om2eps = 1. - 2.*eps
      vorttl(:plon,:) = om2eps*vort(:plon,:) + eps*vortm1(:plon,:)
      divtl(:plon,:) = om2eps*div(:plon,:) + eps*divm1(:plon,:)
   else
!
! Copy time filter fields for first steps when not filtered
!
      vorttl(:plon,:plev) = vort(:plon,:plev)
      divtl(:plon,:plev) = div(:plon,:plev)
   end if
!
! Copy n-1 fields to n-2 time level
!
   um2(:plon,:) = u3m1(:plon,:)
   vm2(:plon,:) = v3m1(:plon,:)
   tm2(:plon,:) = t3m1(:plon,:)
   qm2(:plon,:,:pcnst+pnats) = q3m1(:plon,:,:pcnst+pnats)
   psm2(:plon) = psm1(:plon)

   call plevs0(psm1    ,pintm1  ,pmidm1  ,pdelm1  )
!
! log(pressure) arrays neede by theta, z calculation and in next scan
!
   pmlnm1(:plon,:)  = log(pmidm1(:plon,:))
   pilnm1(:plon,:)  = log(pintm1(:plon,:))
!
! Virtual temperatures, potential temperatures, and model level heights
! for use in radiation, surface, and diffusion calculations
!
   call virtem(t3m1    ,q3m1    ,zvir    ,tvm     )
   call pottem(pintm1  ,pmidm1  ,t3m1    ,cappa   ,psomc   , &
               thm     )
   call zmid(pilnm1(1,plevp),pmlnm1  ,rair    ,gravit  ,tvm     , &
             zm      )
!
! Convert mixing ratio of non-water tracers to mass fraction of total
! atmospheric mass (Overwrite non-water portions of q3m1).
!
   if (.not. adiabatic .and. .not. idealdyn) then
      if (pcnst .gt. 1) call mr2mf( q3m1 )
!
! Calculate total snow precipitation rate:
!
      do i=1,plon
         precst(i) = precsc(i) + precsl(i)
      end do
!
! Do tendency physics: clouds, radiation, surface temp and fluxes,
! pbl and vertical diffusion,rayleigh friction and gravity wave drag
!
      call tphysbc(lat     ,clat(lat),clon   ,coslat  , &
                   oro     ,pmidm1  ,pintm1  ,pdelm1  ,omga    , &
                   u3m1    ,v3m1    ,t3m1    ,q3m1    ,psm1    , &
                   pmlnm1  ,pilnm1  ,zm      ,asdir   ,asdif   , &
                   aldir   ,aldif   ,snowh   ,cnt     ,cnb     , &
                   qrs     ,qrl     ,flwds   ,fsns    ,tvm     , &
                   thm     ,cmfmc   ,lwup    ,precst  ,srfrad  , &
                   sols    ,soll    ,solsd   ,solld   )
!
! Convert mass fractions of non-water tracers back to mixing ratios.
! (Overwrite non-water portions of q3m1).
!
      if (pcnst+pnats .gt. 1) call mf2mr (q3m1)
!
! Surface temperature driver, flux saver and
! heating rate needed for d(theta)/dt computation
!
      call outfld('QRS     ',qrs     ,plond   ,lat     )
      call outfld('QRL     ',qrl     ,plond   ,lat     )
      call outfld('PRECL   ',precl   ,plond   ,lat     )
      call outfld('PRECC   ',precc   ,plond   ,lat     )
!ku --------------------------------------------------------------------
      call outfld('ALDIF     ',aldif     ,plond   ,lat     )
      call outfld('ASDIF     ',asdif     ,plond   ,lat     )
      call outfld('ALDIR     ',aldir     ,plond   ,lat     )
      call outfld('ASDIR     ',asdir     ,plond   ,lat     )
!ku --------------------------------------------------------------------
      do k=1,plev
         do i=1,plon
            hr(i,k) = (qrs(i,k)+ qrl(i,k))*(1.e5/pmid(i,k))**cappa
         end do
      end do
      call outfld('HR      ',hr      ,plond   ,lat     )
   end if

   return
end subroutine linemsbc