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


subroutine tphysidl (ztodt   ,taux    ,tauy    ,etamid  , state   , & 1,14
                     tend)
!----------------------------------------------------------------------- 
! 
! Purpose: 
!  algorithm 1: Held/Suarez IDEALIZED physics
!  algorithm 2: Held/Suarez IDEALIZED physics (Williamson modified stratosphere
!  algorithm 3: Held/Suarez IDEALIZED physics (Lin/Williamson modified strato/meso-sphere
!  algorithm 4: Boer/Denis  IDEALIZED physics
!
! Author: J. Olson
! 
!-----------------------------------------------------------------------
   use shr_kind_mod, only: r8 => shr_kind_r8
   use pmgrid            , only: plev,plat,plevp
   use ppgrid
   use phys_grid         , only: get_lat_all_p, get_rlat_all_p
   use vertical_diffusion, only: vd_intr
   use physics_types     , only: physics_state, physics_tend, physics_ptend
   use geopotential      , only: geopotential_t
   use history,            only: outfld
   use physconst,          only: gravit, cappa, rair
   use constituents,       only: pcnst, pnats
   use abortutils,         only: endrun

   implicit none

#include <comhyb.h>
!
! Input arguments
!
   real(r8), intent(in) :: ztodt                   ! Two times model timestep (2 delta-t)
!
! Output arguments
!
   real(r8), intent(out) :: taux(pcols)            ! X surface stress (zonal)
   real(r8), intent(out) :: tauy(pcols)            ! Y surface stress (meridional)
   real(r8), intent(in)  :: etamid(pver)           ! midpoint values of eta (a+b)

   type(physics_state), intent(inout) :: state
   type(physics_tend ), intent(inout) :: tend
!
!---------------------------Local workspace-----------------------------
!
   type(physics_ptend)   :: ptend                  ! indivdual parameterization tendencies

   integer :: lchnk                                ! chunk identifier
   integer :: ncol                                 ! number of atmospheric columns

   real(r8) clat(pcols)                        ! latitudes(radians) for columns
   real(r8) pmid(pcols,pver)                   ! mid-point pressure
   integer  i,k                                ! Longitude, level indices

   real(r8) tmp                                ! temporary
   real(r8) kf                                 ! 1./efolding_time for wind dissipation
   real(r8) ka                                 ! 1./efolding_time for temperature diss.
   real(r8) kaa                                ! 1./efolding_time for temperature diss.
   real(r8) ks                                 ! 1./efolding_time for temperature diss.
   real(r8) kv                                 ! 1./efolding_time (normalized) for wind
   real(r8) kt                                 ! 1./efolding_time for temperature diss.
   real(r8) trefa                              ! "radiative equilibrium" T
   real(r8) trefc                              ! used in calc of "radiative equilibrium" T
   real(r8) cossq(pcols)                       ! coslat**2
   real(r8) cossqsq(pcols)                     ! coslat**4
   real(r8) sinsq(pcols)                       ! sinlat**2
   real(r8) onemsig                            ! 1. - sigma_reference
   real(r8) efoldf                             ! efolding time for wind dissipation
   real(r8) efolda                             ! efolding time for T dissipation
   real(r8) efoldaa                            ! efolding time for T dissipation
   real(r8) efolds                             ! efolding time for T dissipation
   real(r8) efold_strat                        ! efolding time for T dissipation in Strat
   real(r8) efold_meso                         ! efolding time for T dissipation in Meso
   real(r8) efoldv                             ! efolding time for wind dissipation
   real(r8) p_infint                           ! effective top of model
   real(r8) constw                             ! constant
   real(r8) lapsew                             ! lapse rate
   real(r8) p0strat                            ! threshold pressure
   real(r8) phi0                               ! threshold latitude
   real(r8) dphi0                              ! del-latitude
   real(r8) a0                                 ! coefficient
   real(r8) aeq                                ! 100 mb
   real(r8) apole                              ! 2   mb
   real(r8) pi                                 ! 3.14159...
   real(r8) coslat(pcols)                      ! cosine(latitude)
   real(r8) acoslat                            ! abs(acos(coslat))
   real(r8) constc                             ! constant
   real(r8) lapsec                             ! lapse rate
   real(r8) lapse                              ! lapse rate
   real(r8) h0                                 ! scale height (7 km)
   real(r8) sigmab                             ! threshold sigma level
   real(r8) pressmb                            ! model pressure in mb
   real(r8) t00                                ! minimum reference temperature
   integer  idlflag                            ! Flag to choose which idealized physics
!
!-----------------------------------------------------------------------
!
   idlflag = 1

   lchnk = state%lchnk
   ncol  = state%ncol
!
! Copy pressures into local array
!
   call get_rlat_all_p(lchnk, ncol, clat)
   do i=1,ncol
      coslat (i) = cos(clat(i))
      sinsq  (i) = sin(clat(i))*sin(clat(i))
      cossq  (i) = coslat(i)*coslat(i)
      cossqsq(i) = cossq (i)*cossq (i)
   end do
   do k=1,pver
      do i=1,ncol
         pmid(i,k) = state%pmid(i,k)
      end do
   end do

   if (idlflag == 1) then
!
!-----------------------------------------------------------------------
!
! Held/Suarez IDEALIZED physics algorithm:
!
!   Held, I. M., and M. J. Suarez, 1994: A proposal for the
!   intercomparison of the dynamical cores of atmospheric general
!   circulation models.
!   Bulletin of the Amer. Meteor. Soc., vol. 75, pp. 1825-1830.
!
!-----------------------------------------------------------------------
!
! Add idealized radiative heating rates to temperature tendency
!
      efoldf =  1.
      efolda = 40.
      efolds =  4.
      sigmab =  0.7
      t00    = 200.
!
      onemsig = 1. - sigmab
!
      ka = 1./(86400.*efolda)
      ks = 1./(86400.*efolds)
!
      do k=1,pver
         if (etamid(k) > sigmab) then
            do i=1,ncol
               kt = ka + (ks - ka)*cossqsq(i)*(etamid(k) - sigmab)/onemsig
               tmp   = kt/(1.+ ztodt*kt)
               trefc   = 315. - 60.*sinsq(i)
               trefa = (trefc - 10.*cossq(i)*log((pmid(i,k)/ps0)))*(pmid(i,k)/ps0)**cappa
               trefa    = max(t00,trefa)
               tend%dtdt (i,k) = (trefa - state%t(i,k))*tmp
            end do
         else
            tmp   = ka/(1.+ ztodt*ka)
            do i=1,ncol
               trefc   = 315. - 60.*sinsq(i)
               trefa = (trefc - 10.*cossq(i)*log((pmid(i,k)/ps0)))*(pmid(i,k)/ps0)**cappa
               trefa    = max(t00,trefa)
               tend%dtdt (i,k) = (trefa - state%t(i,k))*tmp
            end do
         endif
      end do
!
! Add diffusion near the surface for the wind fields
!
      do k=1,pver
         do i=1,pcols
            ptend%u(i,k) = 0.
            ptend%v(i,k) = 0.
         end do
      end do
      do i=1,pcols
         taux(i) = 0.
         tauy(i) = 0.
      end do
!
      kf = 1./(86400.*efoldf)
!
      do k=1,pver
         if (etamid(k) > sigmab) then
            kv  = kf*(etamid(k) - sigmab)/onemsig
            tmp = -kv/(1.+ ztodt*kv)
            do i=1,ncol
               ptend%u(i,k) = tmp*state%u(i,k)
               ptend%v(i,k) = tmp*state%v(i,k)
               tend%dudt(i,k)  = tend%dudt(i,k) + ptend%u(i,k)
               tend%dvdt(i,k)  = tend%dvdt(i,k) + ptend%v(i,k)
            end do
         endif
      end do

   elseif (idlflag == 2) then
!
!-----------------------------------------------------------------------
!
! Modified Held/Suarez IDEALIZED physics algorithm
! (modified with Williamson stratosphere):
!
!   Williamson, D. L., J. G. Olson and B. A. Boville, 1998: A comparison
!   of semi--Lagrangian and Eulerian tropical climate simulations.
!   Mon. Wea. Rev., vol 126, pp. 1001-1012.
!
!-----------------------------------------------------------------------
!
! Add idealized radiative heating rates to temperature tendency
!
      efoldf  =  1.
      efolda  = 40.
      efoldaa = 40.
      efolds  =  4.
      sigmab  =  0.7
      t00     = 200.
!
      onemsig = 1. - sigmab
!
      ka  = 1./(86400.*efolda)
      kaa = 1./(86400.*efoldaa)
      ks  = 1./(86400.*efolds)
!
      pi     = 4.*atan(1.)
      phi0   = 60.*pi/180.
      dphi0  = 15.*pi/180.
      a0     = 2.65/dphi0
      aeq    = 10000.
      apole  = 200.
      lapsew = -3.345e-03
      constw = rair*lapsew/gravit
      lapsec =  2.00e-03
      constc = rair*lapsec/gravit
      do k=1,pver
         if (etamid(k) > sigmab) then
            do i=1,ncol
               kt = ka + (ks - ka)*cossqsq(i)*(etamid(k) - sigmab)/onemsig
               acoslat = abs(acos(coslat(i)))
               p0strat = aeq - (aeq - apole)*0.5*(1. + tanh(a0*(acoslat - phi0)))
               tmp     = kt/(1.+ ztodt*kt)
               trefc   = 315. - 60.*sinsq(i)
               trefa   = (trefc - 10.*cossq(i)*log((pmid(i,k)/ps0)))*(pmid(i,k)/ps0)** &
                                                                                     cappa
               trefa   = max(t00,trefa)
               if (pmid(i,k) < 10000.) then
                  trefa = t00*((pmid(i,k)/10000.))**constc
                  tmp   = kaa/(1.+ ztodt*kaa)
               endif
               if (pmid(i,k) < p0strat) then
                  trefa = trefa + t00*( ((pmid(i,k)/p0strat))**constw - 1. )
                  tmp   = kaa/(1.+ ztodt*kaa)
               endif
               tend%dtdt (i,k) = (trefa - state%t(i,k))*tmp
            end do
         else
            do i=1,ncol
               acoslat = abs(acos(coslat(i)))
               p0strat = aeq - (aeq - apole)*0.5*(1. + tanh(a0*(acoslat - phi0)))
               tmp     = ka/(1.+ ztodt*ka)
               trefc   = 315. - 60.*sinsq(i)
               trefa   = (trefc - 10.*cossq(i)*log((pmid(i,k)/ps0)))*(pmid(i,k)/ps0)** &
                                                                                     cappa
               trefa   = max(t00,trefa)
               if (pmid(i,k) < 10000.) then
                  trefa = t00*((pmid(i,k)/10000.))**constc
                  tmp   = kaa/(1.+ ztodt*kaa)
               endif
               if (pmid(i,k) < p0strat) then
                  trefa = trefa + t00*( ((pmid(i,k)/p0strat))**constw - 1. )
                  tmp   = kaa/(1.+ ztodt*kaa)
               endif
               tend%dtdt (i,k) = (trefa - state%t(i,k))*tmp
            end do
         endif
      end do
!
! Add diffusion near the surface for the wind fields
!
      do k=1,pver
         do i=1,pcols
            ptend%u(i,k) = 0.
            ptend%v(i,k) = 0.
         end do
      end do
      do i=1,pcols
         taux(i) = 0.
         tauy(i) = 0.
      end do
!
      kf = 1./(86400.*efoldf)
!
      do k=1,pver
         if (etamid(k) > sigmab) then
            kv  = kf*(etamid(k) - sigmab)/onemsig
            tmp = -kv/(1.+ ztodt*kv)
            do i=1,ncol
               ptend%u(i,k) = tmp*state%u(i,k)
               ptend%v(i,k) = tmp*state%v(i,k)
               tend%dudt(i,k)  = tend%dudt(i,k) + ptend%u(i,k)
               tend%dvdt(i,k)  = tend%dvdt(i,k) + ptend%v(i,k)
            end do
         endif
      end do

   elseif (idlflag == 3) then
!
!-----------------------------------------------------------------------
!
! Held/Suarez IDEALIZED physics algorithm:
! (modified with Lin/Williamson stratosphere/mesosphere):
!
!   Held, I. M., and M. J. Suarez, 1994: A proposal for the
!   intercomparison of the dynamical cores of atmospheric general
!   circulation models.
!   Bulletin of the Amer. Meteor. Soc., vol. 75, pp. 1825-1830.
!
!-----------------------------------------------------------------------
!
! Add idealized radiative heating rates to temperature tendency
!
      efoldf      =  1.
      efolda      = 40.
      efolds      =  4.
      efold_strat = 40.
      efold_meso  = 10.
      efoldv      = 0.5
      sigmab      = 0.7
      lapse       = 0.00225
      h0          = 7000.
      t00         = 200.
      p_infint    = 0.01
!
      onemsig = 1. - sigmab
!
      ka = 1./(86400.*efolda)
      ks = 1./(86400.*efolds)
!
      do k=1,pver
         if (etamid(k) > sigmab) then
            do i=1,ncol
               kt    = ka + (ks - ka)*cossqsq(i)*(etamid(k) - sigmab)/onemsig
               tmp   = kt/(1.+ ztodt*kt)
               trefc = 315. - 60.*sinsq(i)
               trefa = (trefc - 10.*cossq(i)*log((pmid(i,k)/ps0)))*(pmid(i,k)/ps0)**cappa
               trefa = max(t00,trefa)
               tend%dtdt (i,k) = (trefa - state%t(i,k))*tmp
            end do
         else
            do i=1,ncol
               tmp     = ka/(1.+ ztodt*ka)
               pressmb = pmid(i,k)*0.01
               trefc   = 315. - 60.*sinsq(i)
               trefa   = (trefc - 10.*cossq(i)*log((pmid(i,k)/ps0)))*(pmid(i,k)/ps0)** &
                                                                                     cappa
               trefa   = max(t00,trefa)
               if (pressmb <= 100. .and. pressmb > 1.) then
                  trefa = t00 + lapse*h0*coslat(i)*log(100./pressmb)
                  tmp   = (efold_strat-efold_meso)*log(pressmb)/log(100.)
                  tmp   = efold_meso + tmp
                  tmp   = 1./(86400.*tmp)
                  tmp   = tmp/(1.+ ztodt*tmp)
               endif
               if (pressmb <= 1. .and. pressmb > 0.01) then
                  trefa = t00 + lapse*h0*coslat(i)*log(100.*pressmb)
                  tmp   = 1./(86400.*efold_meso)
                  tmp   = tmp/(1.+ ztodt*tmp)
               endif
               if (pressmb <= 0.01) then
                  tmp   = 1./(86400.*efold_meso)
                  tmp   = tmp/(1.+ ztodt*tmp)
               endif
               tend%dtdt (i,k) = (trefa - state%t(i,k))*tmp
            end do
         endif
      end do
!
! Add diffusion near the surface for the wind fields
!
      do k=1,pver
         do i=1,pcols
            ptend%u(i,k) = 0.
            ptend%v(i,k) = 0.
         end do
      end do
      do i=1,pcols
         taux(i) = 0.
         tauy(i) = 0.
      end do
!
      kf = 1./(86400.*efoldf)
!
      do k=1,pver
         if (etamid(k) > sigmab) then
            kv  = kf*(etamid(k) - sigmab)/onemsig
            tmp = -kv/(1.+ ztodt*kv)
            do i=1,ncol
               ptend%u(i,k) = tmp*state%u(i,k)
               ptend%v(i,k) = tmp*state%v(i,k)
               tend%dudt(i,k)  = tend%dudt(i,k) + ptend%u(i,k)
               tend%dvdt(i,k)  = tend%dvdt(i,k) + ptend%v(i,k)
            end do
         else
            do i=1,ncol
               pressmb  = pmid(i,k)*0.01
               if (pressmb <= 100.) then
                  kv       = 1./(86400.*efoldv)
                  tmp      = 1. + tanh(1.5*log10(p_infint/pressmb))
                  kv       = kv*tmp
                  tmp      = -kv/(1.+ ztodt*kv)
                  ptend%u(i,k) = tmp*state%u(i,k)
                  ptend%v(i,k) = tmp*state%v(i,k)
                  tend%dudt(i,k)  = tend%dudt(i,k) + ptend%u(i,k)
                  tend%dvdt(i,k)  = tend%dvdt(i,k) + ptend%v(i,k)
               endif
            end do
         endif
      end do

   else
      write(6,*) 'TPHYSIDL: flag for choosing desired type of idealized ', &
                 'physics ("idlflag") is set incorrectly.'
      write(6,*) 'The valid options are 1, 2, or 3.'
      write(6,*) 'idlflag is currently set to: ',idlflag
      call endrun
   endif
!
! Archive idealized temperature tendency
!
   call outfld('QRS     ',tend%dtdt      ,pcols   ,lchnk      )

   return
end subroutine tphysidl