#include <misc.h>
#include <params.h>
subroutine tfilt_massfix (ztodt, lat, u3m1, u3, & 1,39
v3m1, v3, t3m1, t3, q3m1, &
q3, psm1, ps, alpha, &
etamid, qfcst, vort, div, vortm2, &
divm2, qminus, psm2, um2, &
vm2, tm2, qm2, vortm1, divm1, &
omga, dpsl, dpsm, beta, nlon)
!-----------------------------------------------------------------------
!
! Purpose:
! Time filter (second half of filter for vorticity and divergence only)
!
! Method:
!
! Author:
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use history, only: outfld
use pmgrid
use pspect
use commap
use constituents, only: pcnst, pnats, qmin
use constituents, only: cnst_get_type_byind
use time_manager, only: get_nstep
use physconst, only: cpair, gravit
#if ( defined SCAM)
use scamMod
#endif
#if ( defined BFB_CAM_SCAM_IOP )
use iop
use constituents, only: ppcnst,cnst_get_ind,cnst_name
#endif
implicit none
#include <comctl.h>
#include <comlun.h>
#include <comqfl.h>
#include <comtfc.h>
!
! Input arguments
!
real(r8), intent(in) :: ztodt ! two delta t (unless nstep<2)
real(r8), intent(in) :: qfcst(plond,plev,pcnst)! slt moisture forecast
real(r8), intent(in) :: vort(plond,plev)
real(r8), intent(in) :: div(plond,plev)
real(r8), intent(inout) :: vortm2(plond,plev)
real(r8), intent(inout) :: divm2(plond,plev)
real(r8), intent(in) :: qminus(plond,plev,pcnst)
real(r8), intent(inout) :: psm2(plond)
real(r8), intent(inout) :: um2(plond,plev)
real(r8), intent(inout) :: vm2(plond,plev)
real(r8), intent(inout) :: tm2(plond,plev)
real(r8), intent(inout) :: qm2(plond,plev,pcnst+pnats)
real(r8), intent(inout) :: omga(plond,plev)
real(r8), intent(in) :: dpsl(plond)
real(r8), intent(in) :: dpsm(plond)
real(r8), intent(in) :: beta ! energy fixer coefficient
real(r8), intent(in) :: alpha(pcnst)
real(r8), intent(in) :: etamid(plev) ! vertical coords at midpoints
real(r8), intent(in) :: u3(plond,plev)
real(r8), intent(in) :: v3(plond,plev)
real(r8), intent(inout) :: t3(plond,plev)
integer, intent(in) :: lat
integer, intent(in) :: nlon
! Input/Output arguments
real(r8), intent(inout) :: q3(plond,plev,pcnst+pnats)
real(r8), intent(inout) :: ps(plond)
real(r8), intent(inout) :: vortm1(plond,plev)
real(r8), intent(inout) :: psm1(plond)
real(r8), intent(inout) :: u3m1(plond,plev)
real(r8), intent(inout) :: v3m1(plond,plev)
real(r8), intent(inout) :: t3m1(plond,plev)
real(r8), intent(inout) :: divm1(plond,plev)
real(r8), intent(inout) :: q3m1(plond,plev,pcnst+pnats)
!
! Local workspace
!
integer ifcnt ! Counter
integer :: nstep ! current timestep number
real(r8) qfcst1(plond,plev,pcnst) ! Workspace for qfcst
real(r8) engycorr(plond,plev) ! energy equivalent to T correction
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) pintm1(plond,plevp) ! pressure at model interfaces (n-1)
real(r8) pmidm1(plond,plev) ! pressure at model levels (time n-1)
real(r8) pdelm1(plond,plev) ! pdelm1(k) = pintm1(k+1)-pintm1(k)
real(r8) om2eps
real(r8) corm
real(r8) wm
real(r8) absf
real(r8) worst
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,j,ixcldliq,ixcldice
#if ( defined BFB_CAM_SCAM_IOP )
real(r8) :: t3forecast(plond,plev),delta_t3(plond,plev)
real(r8) :: q3forecast(plond,plev,ppcnst),delta_q3(plond,plev,ppcnst)
real(r8) divt3d(plond,plev) ! 1./pmid
real(r8) divq3d(plond,plev,ppcnst) ! 1./pmid
integer j,ixcldliq,ixcldice
#endif
!-----------------------------------------------------------------------
qfcst1(1:nlon,:,:) = qfcst(i1:nlon+i1-1,:,:)
nstep = get_nstep()
#if ( defined BFB_CAM_SCAM_IOP )
!
! Calculate 3d dynamics term
!
do k=1,plev
do i=1,nlon
divt3d(i,k)=(t3(i,k)-tm2(i,k))/ztodt -t2sav(i,k,lat)
t3forecast(i,k)=tm2(i,k)+ztodt*t2sav(i,k,lat)+ztodt*divt3d(i,k)
end do
end do
do i=1,nlon
do m=1,ppcnst
do k=1,plev
divq3d(i,k,m)= (qfcst1(i,k,m)-qminus(i,k,m))/ztodt
q3forecast(i,k,m)=qminus(i,k,m)+divq3d(i,k,m)*ztodt
end do
end do
end do
q3(:nlon,:,:)=q3forecast(:nlon,:,:)
t3(:nlon,:)=t3forecast(:nlon,:)
qfcst1(:nlon,:,:)=q3(:nlon,:,:)
!
! outflds for iop history tape - to get bit for bit with scam
! the n-1 values are put out. After the fields are written out
! the current time level of info will be buffered for output next
! timestep
!
call outfld('t',t3sav(1,1,lat) ,plond ,lat )
call outfld('q',q3sav(1,1,1,lat),plond ,lat )
call outfld('Ps',pssav(1,lat) ,plond ,lat )
call outfld('u',u3sav(1,1,lat) ,plond ,lat )
call outfld('v',v3sav(1,1,lat) ,plond ,lat )
!
! read single values into plon arrays for output to history tape
! it would be nice if history tape supported 1 dimensional array variables
!
fixmas_plond(:)=fixmassav(lat)
beta_plond(:)=betasav(lat)
clat_plond=clat(lat)
call outfld('fixmas',fixmas_plond,plond ,lat )
call outfld('beta',beta_plond ,plond ,lat )
call outfld('CLAT ',clat_plond ,plond ,lat )
call outfld('divT3d',divt3dsav(1,1,lat) ,plond ,lat )
call outfld('divq3d',divq3dsav(1,1,1,lat) ,plond ,lat )
call cnst_get_ind('CLDLIQ', ixcldliq)
call cnst_get_ind('CLDICE', ixcldice)
if (nstep.eq.0) then
call outfld('dcldliq',q3sav(1,1,ixcldliq,lat),plond ,lat )
call outfld('dcldice',q3sav(1,1,ixcldice,lat),plond ,lat )
else
call outfld('dcldliq',divq3dsav(1,1,ixcldliq,lat),plond ,lat )
call outfld('dcldice',divq3dsav(1,1,ixcldice,lat),plond ,lat )
endif
!
! buffer fields for output at next timestep
!
divq3dsav(:nlon,:,:,lat)= divq3d(:nlon,:,:)
divt3dsav(:nlon,:,lat)= divt3d(:nlon,:)
t3sav(:nlon,:,lat) = t3(:nlon,:)
u3sav(:nlon,:,lat) = u3(:nlon,:)
v3sav(:nlon,:,lat) = v3(:nlon,:)
q3sav(:nlon,:,:,lat) = q3(:nlon,:,:)
pssav(:nlon,lat) = ps(:nlon)
betasav(lat)=beta
fixmassav(lat)=fixmas
#endif
coslat = cos(clat(lat))
do i=1,nlon
rcoslat(i) = 1./coslat
enddo
lfixlim = .true.
!
! 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)
!
! Add temperature correction for energy conservation
!
do k=1,plev
do i=1,nlon
engycorr(i,k) = (cpair/gravit)*beta*pdel(i,k)/ztodt
t3 (i,k) = t3(i,k) + beta
end do
end do
!
! Output Energy correction term
!
call outfld('ENGYCORR',engycorr ,plond ,lat )
!
! 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
if (cnst_get_type_byind(m).eq.'dry' ) then
corm = 1.e36
else
corm = 0.1
end if
do k=1,plev
do i=1,nlon
#if ( defined SCAM )
dqfx3(i,k,m) = dqfxcam(i,k,m)
#else
dqfx3(i,k,m) = alpha(m)*etamid(k)*abs(qfcst1(i,k,m) - qminus(i,k,m))
#if ( defined BFB_CAM_SCAM_IOP )
dqfx3sav(i,k,m) = dqfx3(i,k,m)
#endif
#endif
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)
#ifndef HADVTEST
! TBH: Commented out as of CAM CRB meeting on 6/20/03
! write (6,1000) m,corm,ifcnt,k,lat,wm,worst
#endif
endif
endif
do i=1,nlon
dqfx3(i,k,m) = qfcst1(i,k,m)*dqfx3(i,k,m)/ztodt
#ifdef HADVTEST
q3(i,k,m) = qfcst1(i,k,m)
#else
q3(i,k,m) = qfcst1(i,k,m) + ztodt*dqfx3(i,k,m)
#endif
ta(i,k,m) = (q3(i,k,m) - qminus(i,k,m))/ztodt
end do
end do
end do
#if ( !defined SCAM )
#if ( defined BFB_CAM_SCAM_IOP )
do m=1,pcnst
alpha_plond(:,m)= alphasav(m,lat)
call outfld(alphanam(m),alpha_plond(1,m) ,plond ,lat )
call outfld(dqfxnam(m),dqfx3savm1(1,1,m,lat) ,plond ,lat )
alphasav(m,lat)= alpha(m)
dqfx3savm1(:,:,m,lat)= dqfx3sav(:,:,m)
end do
#endif
#endif
!
! Check for and correct invalid constituents
!
call qneg3 ('TFILT_MASSFIX',lat ,nlon ,plond ,plev , &
pcnst+pnats,qmin ,q3(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
#if ( !defined SCAM )
!
! Calculate vertical motion field
!
call omcalc (rcoslat ,div ,u3 ,v3 ,dpsl , &
dpsm ,pmid ,pdel ,rpmid ,pint(1,plevp), &
omga ,nlon )
#endif
! write(6,*)'tfilt: lat=',lat
! write(6,*)'omga=',omga
!
! Time filter (second half of filter for vorticity and divergence only)
!
! if(lat.eq.2) then
! write(6,*)'tfilt: ps=',psm2(13),psm1(13),ps(13)
! write(6,*)'tfilt: u=',um2(13,18),u3m1(13,18),u3(13,18)
! write(6,*)'tfilt: t=',tm2(13,18),t3m1(13,18),t3(13,18)
! write(6,*)'tfilt: water=',qm2(13,18,1),q3m1(13,18,1),q3(13,18,1)
! write(6,*)'tfilt: cwat=',qm2(13,18,2),q3m1(13,18,2),q3(13,18,2)
! write(6,*)'tfilt: vort=',vortm2(13,18),vortm1(13,18),vort(13,18)
! write(6,*)'tfilt: div=',divm2(13,18),divm1(13,18),div(13,18)
! end if
om2eps = 1. - 2.*eps
#if ( defined SCAM )
!
! for scam I moved the time advance and readiop before dynamics to get forcasted
! values for divt3d and divq3d
!
if (nstep.ge.3) then
#else
if (nstep.ge.2) then
#endif
do k=1,plev
do i=1,nlon
u3m1(i,k) = om2eps*u3m1(i,k) + eps*um2(i,k) + eps*u3(i,k)
v3m1(i,k) = om2eps*v3m1(i,k) + eps*vm2(i,k) + eps*v3(i,k)
t3m1(i,k) = om2eps*t3m1(i,k) + eps*tm2(i,k) + eps*t3(i,k)
q3m1(i,k,1) = om2eps*q3m1(i,k,1) + eps*qm2(i,k,1) + eps*q3(i,k,1)
vortm1(i,k) = om2eps*vortm1(i,k) + eps*vortm2(i,k) + eps*vort(i,k)
divm1(i,k) = om2eps*divm1(i,k) + eps*divm2(i,k) + eps*div(i,k)
end do
do m=2,pcnst+pnats
do i=1,nlon
q3m1(i,k,m) = om2eps*q3m1(i,k,m) + eps*qm2(i,k,m) + eps*q3(i,k,m)
end do
end do
end do
do i=1,nlon
psm1(i) = om2eps*psm1(i) + eps*psm2(i) + eps*ps(i)
end do
end if
call plevs0 (nlon ,plond ,plev ,psm1 ,pintm1 ,pmidm1 ,pdelm1)
!
! Compute time tendencies:comment out since currently not on h-t
!
do k=1,plev
do i=1,nlon
ttend(i,k) = (t3(i,k)-tm2(i,k))/ztodt
utend(i,k) = (u3(i,k)-um2(i,k))/ztodt
vtend(i,k) = (v3(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) = (q3(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)
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