tphysidl.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

!
      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