module tfilt_massfix

contains

subroutine tfilt_massfixrun (ztodt,         lat,    u3m1,   u3,     &
                          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,   hadv ,  &
                          nlon,   pdeldry, pdelm1dry, pdelm2dry)
!----------------------------------------------------------------------- 
! 
! Purpose: 
! Time filter (second half of filter for vorticity and divergence only)
! 
! Method: 
! 
! Author: 
! 
!-----------------------------------------------------------------------
   use shr_kind_mod, only: r8 => shr_kind_r8
   use cam_history,  only: outfld
   use eul_control_mod, only : fixmas,eps
   use pmgrid,       only: plon, plev, plevp, plat
   use pspect
   use commap
   use constituents, only: pcnst, qmin, cnst_cam_outfld, &
                           tottnam, tendnam, cnst_get_type_byind, fixcnam, &
                           hadvnam, vadvnam
   use time_manager, only: get_nstep
   use physconst,    only: cpair, gravit
   use phys_control, only: phys_getopts
   
   implicit none

!
! Input arguments
!
   real(r8), intent(in) :: ztodt                  ! two delta t (unless nstep<2)

   real(r8), intent(inout) :: qfcst(plon,plev,pcnst)! slt moisture forecast
   real(r8), intent(in) :: vort(plon,plev)
   real(r8), intent(in) :: div(plon,plev)
   real(r8), intent(inout) :: vortm2(plon,plev)
   real(r8), intent(inout) :: divm2(plon,plev)
   real(r8), intent(in) :: qminus(plon,plev,pcnst)
   real(r8), intent(inout) :: psm2(plon)
   real(r8), intent(inout) :: um2(plon,plev)
   real(r8), intent(inout) :: vm2(plon,plev)
   real(r8), intent(inout) :: tm2(plon,plev)
   real(r8), intent(inout) :: qm2(plon,plev,pcnst)
   real(r8), intent(inout) :: omga(plon,plev)
   real(r8), intent(in) :: dpsl(plon)
   real(r8), intent(in) :: dpsm(plon)
   real(r8), intent(in) :: beta                   ! energy fixer coefficient
   real(r8), intent(in) :: hadv(plon,plev,pcnst)  ! horizonal q advection tendency
   real(r8), intent(in) :: alpha(pcnst)
   real(r8), intent(in) :: etamid(plev)           ! vertical coords at midpoints 
   real(r8), intent(in) :: u3(plon,plev)
   real(r8), intent(in) :: v3(plon,plev)
   real(r8), intent(inout) :: t3(plon,plev)
   real(r8), intent(inout) :: pdeldry(:,:)      ! dry pressure difference at time n3
   real(r8), intent(inout) :: pdelm1dry(:,:)    ! dry pressure difference at time n3m1
   real(r8), intent(in) :: pdelm2dry(:,:)       ! dry pressure difference at time n3m2


   integer, intent(in) :: lat
   integer, intent(in) :: nlon

! Input/Output arguments

   real(r8), intent(inout) :: q3(plon,plev,pcnst)
   real(r8), intent(inout) :: ps(plon)
   real(r8), intent(inout) :: vortm1(plon,plev)
   real(r8), intent(inout) :: psm1(plon)
   real(r8), intent(inout) :: u3m1(plon,plev)
   real(r8), intent(inout) :: v3m1(plon,plev)
   real(r8), intent(inout) :: t3m1(plon,plev)
   real(r8), intent(inout) :: divm1(plon,plev)
   real(r8), intent(inout) :: q3m1(plon,plev,pcnst)
!
! Local workspace
!
   integer ifcnt                   ! Counter
   integer :: nstep                ! current timestep number
   integer :: timefiltstep         ! 

   real(r8) tfix    (plon)        ! T correction
   real(r8) engycorr(plon,plev)   ! energy equivalent to T correction
   real(r8) rpmid(plon,plev)      ! 1./pmid
   real(r8) pdel(plon,plev)       ! pdel(k)   = pint  (k+1)-pint  (k)
   real(r8) pint(plon,plevp)      ! pressure at model interfaces (n  )
   real(r8) pmid(plon,plev)       ! pressure at model levels (time n)
   real(r8) utend(plon,plev)      ! du/dt
   real(r8) vtend(plon,plev)      ! dv/dt
   real(r8) ttend(plon,plev)      ! dT/dt
   real(r8) qtend(plon,plev,pcnst)! dq/dt
   real(r8) pstend(plon)          ! d(ps)/dt
   real(r8) vadv  (plon,plev,pcnst) ! vertical q advection tendency
   real(r8) pintm1(plon,plevp)    ! pressure at model interfaces (n-1)
   real(r8) pmidm1(plon,plev)     ! pressure at model levels (time n-1)
   real(r8) pdelm1(plon,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(plon,plev,pcnst)   ! total advection of constituents
   real(r8) dqfx3(plon,plev,pcnst)! q tendency due to mass adjustment
   real(r8) coslat                 ! cosine(latitude)
   real(r8) rcoslat(plon)         ! 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,ixnumliq,ixnumice
   real(r8) fixmas_plon(plon)
   real(r8) beta_plon(plon) 
   real(r8) clat_plon(plon)
   real(r8) alpha_plon(plon)

!-----------------------------------------------------------------------
   nstep = get_nstep()

   coslat = cos(clat(lat))
   do i=1,nlon
     rcoslat(i) = 1._r8/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    ,plon   ,plev    ,ps      ,pint    ,pmid    ,pdel)
!
!$OMP PARALLEL DO PRIVATE (K, I)
   do k=1,plev
      do i=1,nlon
         rpmid(i,k) = 1._r8/pmid(i,k)
      enddo
   enddo
!
! Add temperature correction for energy conservation
!
!$OMP PARALLEL DO PRIVATE (K, I)
   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
   do i=1,nlon
      tfix(i) = beta/ztodt
   end do
!
! Output Energy correction term
!
! using do loop and select in order to enable functional parallelism with OpenMP
!$OMP PARALLEL DO PRIVATE (I)
   do i=1,2
      select case (i)
      case (1)
         call outfld ('ENGYCORR',engycorr ,plon   ,lat     )
      case (2)
         call outfld ('TFIX    ',tfix     ,plon   ,lat     )
      end select
   end do

!
! 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_r8
      else
         corm    = 0.1_r8
      end if

!$OMP PARALLEL DO PRIVATE (K, I, IFCNT, WORST, WM, ABSF)
      do k=1,plev
         do i=1,nlon
            dqfx3(i,k,m) = alpha(m)*etamid(k)*abs(qfcst(i,k,m) - qminus(i,k,m))
         end do
         if (lfixlim) then
            ifcnt = 0
            worst = 0._r8
            wm    = 0._r8
            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/real(ifcnt,r8)

! TBH:  Commented out as of CAM CRB meeting on 6/20/03
!              write(iulog,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
            q3   (i,k,m) = qfcst(i,k,m) + ztodt*dqfx3(i,k,m)
            ta   (i,k,m) = (q3    (i,k,m) - qminus(i,k,m))/ztodt
            vadv (i,k,m) = (qfcst(i,k,m) - qminus(i,k,m))/ztodt - hadv(i,k,m)
         end do
      end do
   end do

!$OMP PARALLEL DO PRIVATE (K, I)
   do k=1,plev
      do i=1,nlon
         pdeldry(i,k) = pdel(i,k)*(1._r8-q3(i,k,1))
      end do ! i
   end do ! k

!
! Check for and correct invalid constituents
!
   call qneg3 ('TFILT_MASSFIX',lat   ,nlon    ,plon   ,plev    , &
               1, pcnst, qmin ,q3(1,1,1),.True.)
!
! Send slt tendencies to the history tape
!
!$OMP PARALLEL DO PRIVATE (M)
   do m=1,pcnst
      if ( cnst_cam_outfld(m) ) then
         call outfld(tottnam(m),ta(1,1,m),plon   ,lat     )
      end if
   end do
!
! Calculate vertical motion field
!
   call omcalc (rcoslat ,div     ,u3      ,v3      ,dpsl    ,  &
               dpsm    ,pmid    ,pdel    ,rpmid   ,pint(1,plevp), &
               omga    ,nlon    )

!   write(iulog,*)'tfilt: lat=',lat
!   write(iulog,*)'omga=',omga
!
! Time filter (second half of filter for vorticity and divergence only)
!
!   if(lat.eq.2) then
!      write(iulog,*)'tfilt: ps=',psm2(13),psm1(13),ps(13)
!      write(iulog,*)'tfilt: u=',um2(13,18),u3m1(13,18),u3(13,18)
!      write(iulog,*)'tfilt: t=',tm2(13,18),t3m1(13,18),t3(13,18)
!      write(iulog,*)'tfilt: water=',qm2(13,18,1),q3m1(13,18,1),q3(13,18,1)
!      write(iulog,*)'tfilt: cwat=',qm2(13,18,2),q3m1(13,18,2),q3(13,18,2)
!      write(iulog,*)'tfilt: vort=',vortm2(13,18),vortm1(13,18),vort(13,18)
!      write(iulog,*)'tfilt: div=',divm2(13,18),divm1(13,18),div(13,18)
!   end if

   om2eps = 1._r8 - 2._r8*eps

   if (nstep.ge.2) then
!$OMP PARALLEL DO PRIVATE (K, I, M)
      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
            if (cnst_get_type_byind(m) .eq. 'wet') then 
               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
            endif 
         end do
         do m=2,pcnst
            if (cnst_get_type_byind(m) .eq. 'dry') then 
               do i=1,nlon  ! calculate numerator (timefiltered mass * pdeldry)
                  q3m1(i,k,m) = (om2eps*pdelm1dry(i,k)*q3m1(i,k,m) + &
                       eps*pdelm2dry(i,k)*qm2(i,k,m) + &
                       eps*pdeldry(i,k)*q3(i,k,m))
               end do !i
            endif  !dry
         end do !m
         do i=1,nlon     ! calculate time filtered value of pdeldry
               pdelm1dry(i,k) = om2eps*pdelm1dry(i,k) + &
                    eps*pdelm2dry(i,k) + eps*pdeldry(i,k)
         end do !i
         ! divide time filtered mass*pdeldry by timefiltered pdeldry
         do m=2,pcnst
            if (cnst_get_type_byind(m) == 'dry') then
               do i=1,nlon
                  q3m1(i,k,m) = q3m1(i,k,m)/pdelm1dry(i,k)
               end do !i
            endif ! dry
         end do !m

      end do
      do i=1,nlon
         psm1(i) = om2eps*psm1(i) + eps*psm2(i) + eps*ps(i)
      end do
   end if

   call plevs0 (nlon    ,plon   ,plev    ,psm1    ,pintm1  ,pmidm1  ,pdelm1)
!
! Compute time tendencies:comment out since currently not on h-t
!
!$OMP PARALLEL DO PRIVATE (K, I)
   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

!$OMP PARALLEL DO PRIVATE (M, K, I)
   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

!$OMP PARALLEL DO PRIVATE (M)
   do m=1,pcnst
      if ( cnst_cam_outfld(m) ) then
         call outfld (tendnam(m),qtend(1,1,m),plon,lat)
         call outfld (fixcnam(m),dqfx3(1,1,m),plon,lat)
         call outfld (hadvnam(m),hadv (1,1,m),plon,lat)
         call outfld (vadvnam(m),vadv (1,1,m),plon,lat)
      end if
   end do

! using do loop and select in order to enable functional parallelism with OpenMP
!$OMP PARALLEL DO PRIVATE (I)
   do i=1,4
      select case (i)
      case (1)
         call outfld ('UTEND   ',utend,plon,lat)
      case (2)
         call outfld ('VTEND   ',vtend,plon,lat)
      case (3)
         call outfld ('TTEND   ',ttend,plon,lat)
      case (4)
         call outfld ('LPSTEN  ',pstend,plon,lat)
      end select
   end do

   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_massfixrun

end module tfilt_massfix
