!      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      )   returnend subroutine tphysidl