linemsdyn.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

! Change to Courant limiter for non-triangular truncations.
!
         wind = utfac*u3m2(i,k)**2 + vtfac*v3m2(i,k)**2
         vmaxt(k) = max(wind,vmaxt(k))
      end do
   end do
!
! Variables needed in tphysac
!
   coslat = cos(clat(lat))
   rcoslat = 1./coslat
!
! Set current time pressure arrays for model levels etc.
!
   call plevs0(nlon,plond,plev,psm1,pint,pmid,pdel)
   do k=1,plev
      do i=1,nlon
         rpmid(i,k) = 1./pmid(i,k)
         rpdel(i,k) = 1./pdel(i,k)
      end do
   end do
!
! Accumulate statistics for diagnostic print
!
   call stats(lat,     pint,    pdel,      psm1,   &
              vortm1,  divm1,   t3m1,      q3m1(:,:,1), nlon  )
!
! Compute log(surface pressure) for use by grmult and when adding tendency.
!
   do i=1,nlon
      logpsm1(i) = log(psm1(i))
      logpsm2(i) = log(psm2(i))
   end do
!
! Calculate non-linear terms in tendencies
!
   if (adiabatic) t2(:,:) = 0.
   call grmult(rcoslat ,divm1     ,q3m1(1,1,1),t3m1   ,u3m1    , &
               v3m1    ,vortm1    ,t3m2    ,phis    ,dpsl    , &
               dpsm    ,omga    ,pdel    ,pint(1,plevp),logpsm2, &
               logpsm1 ,rpmid   ,rpdel   ,fu      ,fv      , &
               t2      ,ut      ,vt      ,drhs    ,pmid    , &
               etadot  ,etamid  ,engy    ,ddpn    ,vpdsn   , &
               dpslon  ,dpslat  ,nlon    )
!
! Add tendencies to previous timestep values of surface pressure,
! temperature, and (if spectral transport) moisture.  Store *log* surface
! pressure in bpstr array for transform to spectral space.
!
   do i=1,nlon
      bpstr(i) = logpsm2(i) - ztodt*(vpdsn(i)+ddpn(i))/psm1(i)
   end do

   rhypi = 1./hypi(plevp)
   do k=1,plev
      do i=1,nlon
         ddivdt(i,k) = ztodt*(0.5*divm2(i,k) - divm1(i,k))
         bpstr(i) = bpstr(i) - ddivdt(i,k)*hypd(k)*rhypi
         tdyn(i,k) = t3m2(i,k) + ztodt*t2(i,k)
      end do
   end do

   do k=1,plev
      do kk=1,plev
         do i=1,nlon
#ifdef HADVTEST
!
!jr Remove semi-implicit contribution for advection test
!
!jr              tdyn(i,k) = tdyn(i,k) - ddivdt(i,kk)*tau(kk,k)
#else
            tdyn(i,k) = tdyn(i,k) - ddivdt(i,kk)*tau(kk,k)
#endif
         end do
      end do
   end do

!
! Compute maximum vertical Courant number this latitude.
!
   dtime = get_step_size()
   do k=2,plev
      dtdz = dtime/detam(k-1)
      vcour(k) = 0.
      do i=1,nlon
         lvcour = abs(etadot(i,k))*dtdz
         vcour(k) = max(lvcour,vcour(k))
      end do
   end do

   call outfld('ETADOT  ',etadot,plon,lat)

   inc = 1
   isign = -1
!
! FFT of undifferentiated terms
!
!   if(lat==2)write(6,*)'lat,t3 into fft=',lat,tdyn(13,17)
   ntr = plev
   call fft991(tdyn    ,work    ,trig(1,lat),ifax(1,lat),inc     , &
               plond   ,nlon    ,ntr     ,isign   )
!
! FFT of ps,vort,div
!
   ntr = 1
   call fft991(bpstr   ,work    ,trig(1,lat),ifax(1,lat),inc     , &
               plond   ,nlon    ,ntr     ,isign)
   ntr = plev
   vortdyn(:,:) = vortm2(:,:)
   call fft991(vortdyn ,work    ,trig(1,lat),ifax(1,lat),inc     ,&
               plond   ,nlon    ,ntr     ,isign)
   ntr = plev
   divdyn(:,:) = divm2(:,:)
   call fft991(divdyn  ,work    ,trig(1,lat),ifax(1,lat),inc     , &
               plond   ,nlon    ,ntr     ,isign)
!
! Apply cos(lat) to momentum terms before fft
!
   do k=1,plev
      do i=1,nlon
         fu(i,k) = coslat*fu(i,k)
         fv(i,k) = coslat*fv(i,k)
         ut(i,k) = coslat*ut(i,k)
         vt(i,k) = coslat*vt(i,k)
      end do
   end do
!
! FFT of longitudinally and latitudinally differentiated terms
!
   ntr = plev
   call fft991(fu      ,work    ,trig(1,lat),ifax(1,lat),inc     ,&
               plond   ,nlon    ,ntr     ,isign   )
   call fft991(fv      ,work    ,trig(1,lat),ifax(1,lat),inc     ,&
               plond   ,nlon    ,ntr     ,isign   )
   call fft991(ut      ,work    ,trig(1,lat),ifax(1,lat),inc     ,&
               plond   ,nlon    ,ntr     ,isign   )
   call fft991(vt      ,work    ,trig(1,lat),ifax(1,lat),inc     ,&
               plond   ,nlon    ,ntr     ,isign   )
   call fft991(drhs    ,work    ,trig(1,lat),ifax(1,lat),inc     ,&
               plond   ,nlon    ,ntr     ,isign   )

   if (lat.gt.plat/2) then
!
! First of latitude pair (N. hemisphere), save Fourier coefficients
!
      do i=1,2*nmmax
         grlps2(i) = bpstr(i)
      end do

      do k=1,plev
         do i=1,2*nmmax
#if ( defined PVP )
            grt2(i,k) = tdyn(i,k)
            grz2(i,k) = vortdyn(i,k)
            grd2(i,k) = divdyn(i,k)
            grfu2(i,k) = fu(i,k)
            grfv2(i,k) = fv(i,k)
            grut2(i,k) = ut(i,k)
            grvt2(i,k) = vt(i,k)
            grrh2(i,k) = drhs(i,k)
#else
            grt2(k,i) = tdyn(i,k)
            grz2(k,i) = vortdyn(i,k)
            grd2(k,i) = divdyn(i,k)
            grfu2(k,i) = fu(i,k)
            grfv2(k,i) = fv(i,k)
            grut2(k,i) = ut(i,k)
            grvt2(k,i) = vt(i,k)
            grrh2(k,i) = drhs(i,k)
#endif
         end do
      end do

   else
!
! Second of latitude pair (S. hemisphere), save Fourier coefficients
!
      do i=1,2*nmmax
         grlps1(i) = bpstr(i)
      end do

      do k=1,plev
         do i=1,2*nmmax
#if ( defined PVP )
            grt1(i,k) = tdyn(i,k)
            grz1(i,k) = vortdyn(i,k)
            grd1(i,k) = divdyn(i,k)
            grfu1(i,k) = fu(i,k)
            grfv1(i,k) = fv(i,k)
            grut1(i,k) = ut(i,k)
            grvt1(i,k) = vt(i,k)
            grrh1(i,k) = drhs(i,k)
#else
            grt1(k,i) = tdyn(i,k)
            grz1(k,i) = vortdyn(i,k)
            grd1(k,i) = divdyn(i,k)
            grfu1(k,i) = fu(i,k)
            grfv1(k,i) = fv(i,k)
            grut1(k,i) = ut(i,k)
            grvt1(k,i) = vt(i,k)
            grrh1(k,i) = drhs(i,k)
#endif
         end do
      end do
   end if

   return
end subroutine linemsdyn