#define MODHS 1
subroutine tphysidl(ztodt, state, tend) 1,18
!-----------------------------------------------------------------------
!
! 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 ppgrid, only: pcols, pver
use phys_grid, only: get_rlat_all_p
use physics_types, only: physics_state, physics_tend, physics_ptend, &
physics_ptend_init, physics_update
use physconst, only: gravit, cappa, rair, cpair
use abortutils, only: endrun
use hycoef, only: ps0, etamid
use cam_history, only: outfld
use cam_logfile, only: iulog
use time_manager, only: get_nstep
use check_energy, only: check_energy_chng
implicit none
!
! Input arguments
!
real(r8), intent(in) :: ztodt ! Two times model timestep (2 delta-t)
type(physics_state), intent(inout) :: state
!
! Output arguments
!
type(physics_tend ), intent(inout) :: tend
!
!---------------------------Local workspace-----------------------------
!
type(physics_ptend) :: ptend
integer :: lchnk ! chunk identifier
integer :: ncol ! number of atmospheric columns
integer :: nstep ! timestep number
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
real(r8) :: zero(pcols)
!
!-----------------------------------------------------------------------
!
idlflag = 1
lchnk = state%lchnk
ncol = state%ncol
nstep = get_nstep()
zero = 0._r8
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
! initialize individual parameterization tendencies
call physics_ptend_init(ptend)
ptend%name = 'tphysidl'
ptend%ls = .true.
ptend%lu = .true.
ptend%lv = .true.
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.
!
!-----------------------------------------------------------------------
!
! Compute idealized radiative heating rates (as dry static energy)
!
efoldf = 1._r8
efolda = 40._r8
efolds = 4._r8
sigmab = 0.7_r8
t00 = 200._r8
!
onemsig = 1._r8 - sigmab
!
ka = 1._r8/(86400._r8*efolda)
ks = 1._r8/(86400._r8*efolds)
!
do k=1,pver
if (etamid(k) > sigmab) then
do i=1,ncol
kt = ka + (ks - ka)*cossqsq(i)*(etamid(k) - sigmab)/onemsig
#ifdef MODHS
tmp = kt/(1._r8+ ztodt*kt)
#else
tmp = kt
#endif
trefc = 315._r8 - 60._r8*sinsq(i)
trefa = (trefc - 10._r8*cossq(i)*log((pmid(i,k)/ps0)))*(pmid(i,k)/ps0)**cappa
trefa = max(t00,trefa)
ptend%s(i,k) = (trefa - state%t(i,k))*tmp*cpair
end do
else
#ifdef MODHS
tmp = ka/(1._r8+ ztodt*ka)
#else
tmp = ka
#endif
do i=1,ncol
trefc = 315._r8 - 60._r8*sinsq(i)
trefa = (trefc - 10._r8*cossq(i)*log((pmid(i,k)/ps0)))*(pmid(i,k)/ps0)**cappa
trefa = max(t00,trefa)
ptend%s(i,k) = (trefa - state%t(i,k))*tmp*cpair
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._r8
ptend%v(i,k) = 0._r8
end do
end do
kf = 1._r8/(86400._r8*efoldf)
!
do k=1,pver
if (etamid(k) > sigmab) then
kv = kf*(etamid(k) - sigmab)/onemsig
#ifdef MODHS
tmp = -kv/(1._r8+ ztodt*kv)
#else
tmp = -kv
#endif
do i=1,ncol
ptend%u(i,k) = tmp*state%u(i,k)
ptend%v(i,k) = tmp*state%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.
!
!-----------------------------------------------------------------------
!
! Compute idealized radiative heating rates (as dry static energy)
!
efoldf = 1._r8
efolda = 40._r8
efoldaa = 40._r8
efolds = 4._r8
sigmab = 0.7_r8
t00 = 200._r8
!
onemsig = 1._r8 - sigmab
!
ka = 1._r8/(86400._r8*efolda)
kaa = 1._r8/(86400._r8*efoldaa)
ks = 1._r8/(86400._r8*efolds)
!
pi = 4._r8*atan(1._r8)
phi0 = 60._r8*pi/180._r8
dphi0 = 15._r8*pi/180._r8
a0 = 2.65_r8/dphi0
aeq = 10000._r8
apole = 200._r8
lapsew = -3.345e-03_r8
constw = rair*lapsew/gravit
lapsec = 2.00e-03_r8
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_r8*(1._r8 + tanh(a0*(acoslat - phi0)))
tmp = kt/(1._r8+ ztodt*kt)
trefc = 315._r8 - 60._r8*sinsq(i)
trefa = (trefc - 10._r8*cossq(i)*log((pmid(i,k)/ps0)))*(pmid(i,k)/ps0)** &
cappa
trefa = max(t00,trefa)
if (pmid(i,k) < 10000._r8) then
trefa = t00*((pmid(i,k)/10000._r8))**constc
tmp = kaa/(1._r8+ ztodt*kaa)
endif
if (pmid(i,k) < p0strat) then
trefa = trefa + t00*( ((pmid(i,k)/p0strat))**constw - 1._r8 )
tmp = kaa/(1._r8+ ztodt*kaa)
endif
ptend%s(i,k) = (trefa - state%t(i,k))*tmp*cpair
end do
else
do i=1,ncol
acoslat = abs(acos(coslat(i)))
p0strat = aeq - (aeq - apole)*0.5_r8*(1._r8 + tanh(a0*(acoslat - phi0)))
tmp = ka/(1._r8+ ztodt*ka)
trefc = 315._r8 - 60._r8*sinsq(i)
trefa = (trefc - 10._r8*cossq(i)*log((pmid(i,k)/ps0)))*(pmid(i,k)/ps0)** &
cappa
trefa = max(t00,trefa)
if (pmid(i,k) < 10000._r8) then
trefa = t00*((pmid(i,k)/10000._r8))**constc
tmp = kaa/(1._r8+ ztodt*kaa)
endif
if (pmid(i,k) < p0strat) then
trefa = trefa + t00*( ((pmid(i,k)/p0strat))**constw - 1._r8 )
tmp = kaa/(1._r8+ ztodt*kaa)
endif
ptend%s(i,k) = (trefa - state%t(i,k))*tmp*cpair
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._r8
ptend%v(i,k) = 0._r8
end do
end do
kf = 1._r8/(86400._r8*efoldf)
!
do k=1,pver
if (etamid(k) > sigmab) then
kv = kf*(etamid(k) - sigmab)/onemsig
tmp = -kv/(1._r8+ ztodt*kv)
do i=1,ncol
ptend%u(i,k) = tmp*state%u(i,k)
ptend%v(i,k) = tmp*state%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.
!
!-----------------------------------------------------------------------
!
! Compute idealized radiative heating rates (as dry static energy)
!
efoldf = 1._r8
efolda = 40._r8
efolds = 4._r8
efold_strat = 40._r8
efold_meso = 10._r8
efoldv = 0.5_r8
sigmab = 0.7_r8
lapse = 0.00225_r8
h0 = 7000._r8
t00 = 200._r8
p_infint = 0.01_r8
!
onemsig = 1._r8 - sigmab
!
ka = 1._r8/(86400._r8*efolda)
ks = 1._r8/(86400._r8*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._r8+ ztodt*kt)
trefc = 315._r8 - 60._r8*sinsq(i)
trefa = (trefc - 10._r8*cossq(i)*log((pmid(i,k)/ps0)))*(pmid(i,k)/ps0)**cappa
trefa = max(t00,trefa)
ptend%s(i,k) = (trefa - state%t(i,k))*tmp*cpair
end do
else
do i=1,ncol
tmp = ka/(1._r8+ ztodt*ka)
pressmb = pmid(i,k)*0.01_r8
trefc = 315._r8 - 60._r8*sinsq(i)
trefa = (trefc - 10._r8*cossq(i)*log((pmid(i,k)/ps0)))*(pmid(i,k)/ps0)** &
cappa
trefa = max(t00,trefa)
if (pressmb <= 100._r8 .and. pressmb > 1._r8) then
trefa = t00 + lapse*h0*coslat(i)*log(100._r8/pressmb)
tmp = (efold_strat-efold_meso)*log(pressmb)/log(100._r8)
tmp = efold_meso + tmp
tmp = 1._r8/(86400._r8*tmp)
tmp = tmp/(1._r8+ ztodt*tmp)
endif
if (pressmb <= 1._r8 .and. pressmb > 0.01_r8) then
trefa = t00 + lapse*h0*coslat(i)*log(100._r8*pressmb)
tmp = 1._r8/(86400._r8*efold_meso)
tmp = tmp/(1._r8+ ztodt*tmp)
endif
if (pressmb <= 0.01_r8) then
tmp = 1._r8/(86400._r8*efold_meso)
tmp = tmp/(1._r8+ ztodt*tmp)
endif
ptend%s(i,k) = (trefa - state%t(i,k))*tmp*cpair
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._r8
ptend%v(i,k) = 0._r8
end do
end do
kf = 1._r8/(86400._r8*efoldf)
!
do k=1,pver
if (etamid(k) > sigmab) then
kv = kf*(etamid(k) - sigmab)/onemsig
tmp = -kv/(1._r8+ ztodt*kv)
do i=1,ncol
ptend%u(i,k) = tmp*state%u(i,k)
ptend%v(i,k) = tmp*state%v(i,k)
end do
else
do i=1,ncol
pressmb = pmid(i,k)*0.01_r8
if (pressmb <= 100._r8) then
kv = 1._r8/(86400._r8*efoldv)
tmp = 1._r8 + tanh(1.5_r8*log10(p_infint/pressmb))
kv = kv*tmp
tmp = -kv/(1._r8+ ztodt*kv)
ptend%u(i,k) = tmp*state%u(i,k)
ptend%v(i,k) = tmp*state%v(i,k)
endif
end do
endif
end do
else
write(iulog,*) 'TPHYSIDL: flag for choosing desired type of idealized ', &
'physics ("idlflag") is set incorrectly.'
write(iulog,*) 'The valid options are 1, 2, or 3.'
write(iulog,*) 'idlflag is currently set to: ',idlflag
call endrun('tphysidl: invalid option')
endif
! update the state and total physics tendency
call physics_update(state, tend, ptend, ztodt)
! Can't turn on conservation error messages unless the appropriate heat
! surface flux is computed and supplied as an argument to
! check_energy_chng to account for how the ideal physics forcings are
! changing the total energy.
call check_energy_chng(state, tend, "tphysidl", nstep, ztodt, zero, zero, zero, zero)
call outfld('QRS', tend%dtdt, pcols, lchnk)
end subroutine tphysidl