tfilt_massfix.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

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

subroutine tfilt_massfix (ztodt   ,lat     ,u3m1    ,v3m1    ,t3m1    , &
                          q3      ,q3m1    ,ps      ,cwava   ,alpha   , &
                          etamid  ,qfcst   ,div     ,phis    ,omga    , &
                          dpsl    ,dpsm    ,nlon    ,t3      ,beta )
!-----------------------------------------------------------------------
!
! Purpose:
! Atmosphere and constituent mass fixer
!
! Author:  J. Olson
!
!-----------------------------------------------------------------------

  use precision
  use pmgrid
  use pspect
  use commap
  use history, only: outfld
  use constituents, only: pcnst, pnats, qmin
  use tracers, only: tottnam

  implicit none

#include <comctl.h>
#include <comlun.h>
#include <comqfl.h>
#include <comtfc.h>
!
!------------------------------Arguments--------------------------------
!
  real(r8), intent(in)   :: ztodt                          ! timestep
  integer , intent(in)   :: lat                            ! latitude index
  real(r8), intent(in)   :: u3m1    (plond,plev)             ! u
  real(r8), intent(in)   :: v3m1    (plond,plev)             ! v
  real(r8), intent(inout):: t3m1    (plond,plev)             ! T
  real(r8), intent(in)   :: q3      (plond,plev,pcnst+pnats) ! q + consts (time level n-1)
  real(r8), intent(inout):: q3m1    (plond,plev,pcnst+pnats) ! q + consts (time level n  )
  real(r8), intent(inout):: ps    (plond)                  ! Ps
  real(r8), intent(in)   :: cwava                          ! weight for global integrals
  real(r8), intent(in)   :: alpha (pcnst)                  ! slt fixer coefficient
  real(r8), intent(in)   :: etamid(plev)                   ! vertical coords at midpoints 
  real(r8), intent(in)   :: qfcst (plond,plev,pcnst)       ! slt moisture forecast
  real(r8), intent(in)   :: div   (plond,plev)             ! divergence
  real(r8), intent(in)   :: phis  (plond)                  ! Geopotential field
  real(r8), intent(out)  :: omga  (plond,plev)             ! vertical motion
  real(r8), intent(in)   :: dpsl  (plond)                  ! long comp of grad ln(ps)
  real(r8), intent(in)   :: dpsm  (plond)                  ! lat comp  of grad ln(ps)
  integer , intent(in)   :: nlon                           ! number of longitudes
  real(r8), intent(in)   :: t3(plond,plev) ! temperature
  real(r8), intent(in)   :: beta                           ! energy fixer coefficient
!
!---------------------------Local workspace-----------------------------
!
  integer ifcnt                     ! Counter
  real(r8) rpmid (plond,plev)       ! 1./pmid
  real(r8) pdel  (plond,plev)       ! pdel(k)   = pint  (k+1)-pint  (k)
  real(r8) pint  (plond,plevp)      ! pressure at model interfaces (n  )
  real(r8) pmid  (plond,plev)       ! pressure at model levels (time n)
! real(r8) utend (plond,plev)       ! du/dt
! real(r8) vtend (plond,plev)       ! dv/dt
! real(r8) ttend (plond,plev)       ! dT/dt
! real(r8) qtend (plond,plev,pcnst) ! dq/dt
! real(r8) pstend(plond)            ! d(ps)/dt
  real(r8) psl   (plond)            ! sea level pressure
  real(r8) corm                     ! fixer limit
  real(r8) wm                       ! accumulator 
  real(r8) absf                     ! absolute value of fixer
  real(r8) worst                    ! largest fixer contribution at each model level
  logical lfixlim                   ! flag to turn on fixer limiter

  real(r8) ta    (plond,plev,pcnst) ! total advection of constituents
  real(r8) dqfx3 (plond,plev,pcnst) ! q tendency due to mass adjustment
  real(r8) coslat                   ! cosine(latitude)
  real(r8) rcoslat(plond)           ! 1./cosine(latitude)
! real(r8) engt                     ! Thermal   energy integral
! real(r8) engk                     ! Kinetic   energy integral
! real(r8) engp                     ! Potential energy integral

  integer i, k, m                   ! indices
!
!-----------------------------------------------------------------------
!
  coslat  = cos(clat(lat))
  do i=1,nlon
     rcoslat(i) = 1./coslat
  enddo
  lfixlim = .true.
  corm    = 0.1
!
! Fix temperature to counter-act energy imbalance created by dynamics only
!
  do k=1,plev
    do i=1,nlon
      t3m1(i,k) = t3m1(i,k) + beta*abs(t3m1(i,k) - t3(i,k))
    end do
  end do
!
! Set average dry mass to specified constant preserving horizontal
! gradients of ln(ps). Proportionality factor was calculated in STEPON
! for nstep=0 or SCAN2 otherwise from integrals calculated in INIDAT
! and SCAN2 respectively.
! Set p*.
!
  do i=1,nlon
     ps(i) = ps(i)*fixmas
  end do
!
! Set current time pressure arrays for model levels etc.
!
  call plevs0(nlon    ,plond   ,plev    ,ps      ,pint    ,pmid    ,pdel)
!
  rpmid(:nlon,:plev) = 1./pmid(:nlon,:plev)
!
! Compute q tendency due to mass adjustment
! If LFIXLIM = .T., then:
! Check to see if fixer is exceeding a desired fractional limit of the
! constituent mixing ratio ("corm").  If so, then limit the fixer to
! that specified limit.
!
  do m=1,pcnst
     do k=1,plev
        do i=1,nlon
           dqfx3(i,k,m) = alpha(m)*etamid(k)*abs(qfcst(i,k,m) - q3(i,k,m))
        end do

        if (lfixlim) then
           ifcnt = 0
           worst = 0.
           wm    = 0.
           do i = 1,nlon
              absf = abs(dqfx3(i,k,m))
              if (absf.gt.corm) then
                 ifcnt = ifcnt + 1
                 worst = max(absf,worst)
                 wm = wm + absf
                 dqfx3(i,k,m) = sign(corm,dqfx3(i,k,m))
              endif
           end do
           if (ifcnt.gt.0) then
              wm = wm/float(ifcnt)
              write (6,1000) m,corm,ifcnt,k,lat,wm,worst
           endif
        endif

        do i=1,nlon
           dqfx3(i,k,m) = qfcst(i,k,m)*dqfx3(i,k,m)/ztodt
#ifdef HADVTEST
           q3m1(i,k,m) = qfcst(i,k,m)
#else
           q3m1(i,k,m) = qfcst(i,k,m) + ztodt*dqfx3(i,k,m)
#endif
           ta(i,k,m) = (q3m1(i,k,m) - q3(i,k,m))/ztodt
        end do
     end do
  end do

!
! Check for and correct invalid constituents
!
  call qneg3('TFILT_MASSFIX',lat   ,nlon    ,plond   ,plev    , &
             pcnst+pnats,qmin ,q3m1(1,1,1))
!
! Send slt tendencies to the history tape
!
  do m=1,pcnst
     call outfld(tottnam(m),ta(1,1,m)   ,plond   ,lat     )
  end do
!
! Calculate vertical motion field
!
  call omcalc(rcoslat ,div     ,u3m1    ,v3m1    ,dpsl    , &
              dpsm    ,pmid    ,pdel    ,rpmid   ,pint(1,plevp), &
              omga    ,nlon    )

  call plevs0(nlon    ,plond   ,plev    ,ps      ,pint    ,pmid    ,pdel)
!
! Compute time tendencies:comment out since currently not on h-t
!
!      do k=1,plev
!        do i=1,nlon
!          ttend(i,k) = (t3m1(i,k)-tm2(i,k))/ztodt
!          utend(i,k) = (u3m1(i,k)-um2(i,k))/ztodt
!          vtend(i,k) = (v3m1(i,k)-vm2(i,k))/ztodt
!        end do  
!      end do  
!      do m=1,pcnst
!        do k=1,plev
!          do i=1,nlon
!            qtend(i,k,m) = (q3m1(i,k,m) - qm2(i,k,m))/ztodt
!          end do
!        end do
!      end do
!      do i=1,nlon
!        pstend(i) = (ps(i) - psm2(i))/ztodt
!      end do

!
! do m=1,pcnst
!    call outfld(tendnam(m),qtend(1,1,m),plond   ,lat     )
! end do
! call outfld('UTEND   ',utend   ,plond   ,lat     )
! call outfld('VTEND   ',vtend   ,plond   ,lat     )
! call outfld('TTEND   ',ttend   ,plond   ,lat     )
! call outfld('LPSTEN  ',pstend  ,plond   ,lat     )

  call plevs0(nlon    ,plond   ,plev    ,ps      ,pint    ,pmid    ,pdel)

  return
1000 format(' TIMEFILTER: WARNING: fixer for tracer ',i3,' exceeded ', &
            f8.5,' for ',i5,' points at k,lat = ',2i4,' Avg/Worst = ',1p2e10.2)
end subroutine tfilt_massfix