C**** QUSEM12 E001M12 SOMTQ QUSB261AM12
C**** QUSBM9=QUSB140M9 with correction to second order moment calc.
C**** Changes for constant pressure above LS1 + REAL*8
C**** QUS   is Russell quadratic upstream scheme for temperature
C**** and water vapor advection, with limits applied to water vapor.
C**** Changes for constant pressure above LS1
C**** FQU,FQV for additional diagnostics
C**** Routines included: AADVT, AADVTX, AADVTY, AADVTZ


      MODULE QUSCOM 13,1
!@sum  QUSCOM contains gcm-specific advection parameters/workspace
!@auth Maxwell Kelley
      USE QUSDEF
      IMPLICIT NONE
      SAVE
      INTEGER :: IM,JM,LM
      INTEGER :: XSTRIDE,YSTRIDE,ZSTRIDE
      REAL*8 :: BYIM
C**** AIR MASS FLUXES
      REAL*8, DIMENSION(:,:,:), ALLOCATABLE :: MFLX
C**** WORKSPACE FOR AADVTX
cc    REAL*8, DIMENSION(:), ALLOCATABLE :: AM,F_I
cc    REAL*8, DIMENSION(:,:), ALLOCATABLE :: FMOM_I
C**** WORKSPACE FOR AADVTY
cc    REAL*8, DIMENSION(:), ALLOCATABLE :: BM,F_J
cc    REAL*8, DIMENSION(:,:), ALLOCATABLE :: FMOM_J
C**** WORKSPACE FOR AADVTZ
cc    REAL*8, DIMENSION(:), ALLOCATABLE :: CM,F_L
cc    REAL*8, DIMENSION(:,:), ALLOCATABLE :: FMOM_L
      END MODULE QUSCOM


      SUBROUTINE init_QUS(IM_GCM,JM_GCM,LM_GCM) 1,2
!@sum  init_QUS sets gcm-specific advection parameters/workspace
!@auth Maxwell Kelley
      use QUSCOM
      USE PARAM
      INTEGER, INTENT(IN) :: IM_GCM,JM_GCM,LM_GCM
C**** SET RESOLUTION
      IM = IM_GCM
      JM = JM_GCM
      LM = LM_GCM
      BYIM = 1.D0/REAL(IM,KIND=8)
      XSTRIDE = 1
      YSTRIDE = IM
      ZSTRIDE = IM*JM
C**** ALLOCATE SPACE FOR AIR MASS FLUXES
      ALLOCATE(MFLX(IM,JM,LM))
C**** ALLOCATE WORKSPACE FOR AADVTX
cc    ALLOCATE(AM(IM),F_I(IM),FMOM_I(NMOM,IM))
C**** ALLOCATE WORKSPACE FOR AADVTY
cc    ALLOCATE(BM(JM),F_J(JM),FMOM_J(NMOM,JM))
C**** ALLOCATE WORKSPACE FOR AADVTZ
cc    ALLOCATE(CM(LM),F_L(LM),FMOM_L(NMOM,LM))

      call sync_param("prather_limits",prather_limits)

      RETURN
      END SUBROUTINE init_QUS


      SUBROUTINE AADVT (MA,RM,RMOM,SD,PU,PV,DT,QLIMIT,FQU,FQV) 1,6
!@sum  AADVT advection driver
!@auth G. Russell, modified by Maxwell Kelley
c****
c**** AADVT advects tracers using the Quadradic Upstream Scheme.
c****
c**** input:
c****  pu,pv,sd (kg/s) = east-west,north-south,vertical mass fluxes
c****      qlimit = whether moment limitations should be used
C****         DT (s) = time step
c****
c**** input/output:
c****     rm = tracer concentration
c****   rmom = moments of tracer concentration
c****     ma (kg) = fluid mass
c****
      USE QUSDEF
      USE QUSCOM, ONLY : IM,JM,LM, MFLX
      IMPLICIT NONE

      REAL*8, dimension(im,jm,lm) :: rm,ma
      REAL*8, dimension(NMOM,IM,JM,LM) :: rmom

      REAL*8, INTENT(IN) :: DT
      REAL*8, dimension(im,jm,lm), intent(in) :: pu,pv
      REAL*8, dimension(im,jm,lm-1), intent(in) :: sd
      LOGICAL, INTENT(IN) :: QLIMIT

      REAL*8, dimension(im,jm), intent(inout) :: fqu,fqv

      INTEGER :: I,J,L,N
      REAL*8 :: BYMA

C**** Initialise diagnostics
      FQU=0.  ; FQV=0.

C**** Fill in values at the poles
!$OMP  PARALLEL DO PRIVATE(I,L,N)
      DO L=1,LM
         DO I=2,IM
           RM(I,1 ,L) =   RM(1,1 ,L)
           RM(I,JM,L) =   RM(1,JM,L)
           DO N=1,NMOM
             RMOM(N,I,1 ,L) =  RMOM(N,1,1 ,L)
             RMOM(N,I,JM,L) =  RMOM(N,1,JM,L)
         enddo
         enddo
      enddo
!$OMP  END PARALLEL DO
C****
C**** convert from concentration to mass units
C****
!$OMP  PARALLEL DO PRIVATE(I,J,L)
      DO L=1,LM
      DO J=1,JM
      DO I=1,IM
         RM(I,J,L)=RM(I,J,L)*MA(I,J,L)
         RMOM(:,I,J,L)=RMOM(:,I,J,L)*MA(I,J,L)
      enddo
      enddo
      enddo
!$OMP  END PARALLEL DO
C****
C**** Advect the tracer using the quadratic upstream scheme
C****
CC    mflx(:,:,:)=pu(:,:,:)*(.5*dt)
!$OMP  PARALLEL DO PRIVATE(L)
       DO L=1,LM
          mflx(:,:,l)=pu(:,:,l)*(.5*dt)
       ENDDO
!$OMP  END PARALLEL DO
      CALL AADVTX (RM,RMOM,MA,MFLX,QLIMIT,FQU)
CC    mflx(:,1:jm-1,:)=pv(:,2:jm,:)*dt
CC    mflx(:,jm,:)=0.
!$OMP  PARALLEL DO PRIVATE(L)
       DO L=1,LM
          mflx(:,1:jm-1,l)=pv(:,2:jm,l)*dt
          mflx(:,jm,l)=0.
       ENDDO
!$OMP  END PARALLEL DO
      CALL AADVTY (RM,RMOM,MA,MFLX,QLIMIT,FQV)
CC    mflx(:,:,1:lm-1)=sd(:,:,1:lm-1)*(-dt)
CC    mflx(:,:,lm)=0.
!$OMP  PARALLEL DO PRIVATE(L)
      DO L=1,LM
         IF(L.NE.LM)  THEN
            MFLX(:,:,L)=SD(:,:,L)*(-DT)
         ELSE
            MFLX(:,:,L)=0.
         END IF
      ENDDO
!$OMP  END PARALLEL DO
      CALL AADVTZ (RM,RMOM,MA,MFLX,QLIMIT)
CC    mflx(:,:,:)=pu(:,:,:)*(.5*dt)
!$OMP  PARALLEL DO PRIVATE(L)
       DO L=1,LM
          mflx(:,:,l)=pu(:,:,l)*(.5*dt)
       ENDDO
!$OMP  END PARALLEL DO
      CALL AADVTX (RM,RMOM,MA,MFLX,QLIMIT,FQU)
C****
C**** convert from mass to concentration units
C****
!$OMP  PARALLEL DO PRIVATE(I,J,L,BYMA)
      DO L=1,LM
      DO J=1,JM
      DO I=1,IM
         BYMA = 1.D0/MA(I,J,L)
         RM(I,J,L)=RM(I,J,L)*BYMA
         RMOM(:,I,J,L)=RMOM(:,I,J,L)*BYMA
      enddo
      enddo
      enddo
!$OMP  END PARALLEL DO
      RETURN
      END


      subroutine aadvtx(rm,rmom,mass,mu,qlimit,fqu) 2,4
!@sum  AADVTX advection driver for x-direction
!@auth Maxwell Kelley
c****
c**** aadvtx advects tracers in the west to east direction using the
c**** quadratic upstream scheme.  if qlimit is true, the moments are
c**** limited to prevent the mean tracer from becoming negative.
c****
c**** input:
c****     mu (kg) = west-east mass flux, positive eastward
c****      qlimit = whether moment limitations should be used
c****
c**** input/output:
c****     rm (kg) = tracer mass
c****   rmom (kg) = moments of tracer mass
c****   mass (kg) = fluid mass
c****
      use QUSDEF
ccc   use QUSCOM, only : im,jm,lm, xstride,am,f_i,fmom_i
      use QUSCOM, only : im,jm,lm, xstride
      implicit none
      REAL*8, dimension(im,jm,lm) :: rm,mass,mu,hfqu
      REAL*8, dimension(NMOM,IM,JM,LM) :: rmom
      logical ::  qlimit
      REAL*8, INTENT(OUT), DIMENSION(IM,JM) :: FQU
      REAL*8  AM(IM), F_I(IM), FMOM_I(NMOM,IM)
     &     ,MASS_I(IM), COURMAX, BYNSTEP
      integer :: i,ip1,j,l,ierr,nerr,ICKERR,ns,nstep
c**** loop over layers and latitudes
      ICKERR=0
!$OMP  PARALLEL DO PRIVATE(J,L,AM,F_I,FMOM_I,IERR,NERR,
!$OMP*             I,IP1,NS,NSTEP,BYNSTEP,COURMAX,MASS_I)
!$OMP* SHARED(IM,QLIMIT,XSTRIDE)
!$OMP* REDUCTION(+:ICKERR)
      do l=1,lm
      do j=2,jm-1
c****
c**** decide how many timesteps to take
c****
      nstep=0
      courmax = 2.
      do while(courmax.gt.1. .and. nstep.lt.20)
        nstep = nstep+1
        bynstep = 1d0/real(nstep,kind=8)
        am(:) = mu(:,j,l)*bynstep
        mass_i(:)  = mass(:,j,l)
        courmax = 0.
        do ns=1,nstep
          i = im
          do ip1=1,im
            if(am(i).gt.0.) then
               courmax = max(courmax,+am(i)/mass_i(i))
            else
               courmax = max(courmax,-am(i)/mass_i(ip1))
            endif
            i = ip1
          enddo
          if(ns.lt.nstep) then
             i = im
             do ip1=1,im
                mass_i(ip1) = mass_i(ip1) + (am(i)-am(ip1))
                i = ip1
             enddo
          endif
        enddo
      enddo
      if(courmax.gt.1.) then
         write(6,*) 'aadvtx: j,l,courmax=',j,l,courmax
         ICKERR=ICKERR+1
      endif

c      am(:) = mu(:,j,l)*bynstep ! am already set
      hfqu(:,j,l)  = 0.
c****
c**** loop over timesteps
c****
      do ns=1,nstep
c****
c**** call 1-d advection routine
c****
      call adv1d(rm(1,j,l),rmom(1,1,j,l), f_i,fmom_i, mass(1,j,l),
     &        am, im, qlimit,xstride,xdir,ierr,nerr)
      if (ierr.gt.0) then
        write(6,*) "Error in aadvtx: i,j,l=",nerr,j,l
        if (ierr.eq.2) then
          write(0,*) "Error in qlimit: abs(a) > 1"
CCC       call stop_model('Error in qlimit: abs(a) > 1',11)
          ICKERR=ICKERR+1
        end if
      end if
c****
c**** store tracer flux in fqu array
c****
CCC   fqu(:,j)  = fqu(:,j) + f_i(:)
      hfqu(:,j,l)  = hfqu(:,j,l) + f_i(:)
      enddo ! ns
      enddo ! j
      enddo ! l
!$OMP  END PARALLEL DO
c
c     now sum into fqu
c
!$OMP  PARALLEL DO PRIVATE(J,L)
      do j=2,jm-1
      do l=1,lm
         fqu(:,j)  = fqu(:,j) + hfqu(:,j,l)
      enddo ! j
      enddo ! l
!$OMP  END PARALLEL DO
C
      IF(ICKERR.GT.0)  CALL stop_model('Stopped in aadvtx',11)
C
      return
c****
      end subroutine aadvtx


      subroutine aadvty(rm,rmom,mass,mv,qlimit,fqv) 1,4
!@sum  AADVTY advection driver for y-direction
!@auth Maxwell Kelley
c****
c**** aadvty advects tracers in the south to north direction using the
c**** quadratic upstream scheme.  if qlimit is true, the moments are
c**** limited to prevent the mean tracer from becoming negative.
c****
c**** input:
c****     mv (kg) = north-south mass flux, positive northward
c****      qlimit = whether moment limitations should be used
c****
c**** input/output:
c****     rm (kg) = tracer mass
c****   rmom (kg) = moments of tracer mass
c****   mass (kg) = fluid mass
c****
      use QUSDEF
ccc   use QUSCOM, only : im,jm,lm, ystride,bm,f_j,fmom_j, byim
      use QUSCOM, only : im,jm,lm, ystride,               byim
      implicit none
      REAL*8, dimension(im,jm,lm) :: rm,mass,mv
      REAL*8, dimension(NMOM,IM,JM,LM) :: rmom
      logical ::  qlimit
      REAL*8, intent(out), dimension(im,jm) :: fqv
      REAL*8  HFQV(IM,JM,LM),BM(JM),F_J(JM),FMOM_J(NMOM,JM)
      integer :: i,j,l,ierr,nerr,ICKERR
      REAL*8 ::
     &     m_sp,m_np,rm_sp,rm_np,rzm_sp,rzm_np,rzzm_sp,rzzm_np
c**** loop over layers
      ICKERR=0
!$OMP  PARALLEL DO PRIVATE(I,L,M_SP,M_NP,RM_SP,RM_NP,RZM_SP,RZZM_SP,
!$OMP*             F_J,FMOM_J,RZM_NP,RZZM_NP,BM,IERR,NERR)
!$OMP* SHARED(JM,QLIMIT,YSTRIDE)
!$OMP* REDUCTION(+:ICKERR)
      do l=1,lm
c**** scale polar boxes to their full extent
      mass(:,1:jm:jm-1,l)=mass(:,1:jm:jm-1,l)*im
      m_sp = mass(1,1 ,l)
      m_np = mass(1,jm,l)
      rm(:,1:jm:jm-1,l)=rm(:,1:jm:jm-1,l)*im
      rm_sp = rm(1,1 ,l)
      rm_np = rm(1,jm,l)
      do i=1,im
         rmom(:,i,1 ,l)=rmom(:,i,1 ,l)*im
         rmom(:,i,jm,l)=rmom(:,i,jm,l)*im
      enddo
      rzm_sp  = rmom(mz ,1,1 ,l)
      rzzm_sp = rmom(mzz,1,1 ,l)
      rzm_np  = rmom(mz ,1,jm,l)
      rzzm_np = rmom(mzz,1,jm,l)
c**** loop over longitudes
      do i=1,im
c****
c**** load 1-dimensional arrays
c****
      bm   (:) = mv(i,:,l) !/nstep
      bm(jm)= 0.
      rmom(ihmoms,i,1 ,l) = 0.! horizontal moments are zero at pole
      rmom(ihmoms,i,jm,l) = 0.
c****
c**** call 1-d advection routine
c****
      call adv1d(rm(i,1,l),rmom(1,i,1,l), f_j,fmom_j, mass(i,1,l),
     &     bm, jm,qlimit,ystride,ydir,ierr,nerr)
      if (ierr.gt.0) then
        write(6,*) "Error in aadvty: i,j,l=",i,nerr,l
        if (ierr.eq.2) then
          write(0,*) "Error in qlimit: abs(b) > 1"
ccc       call stop_model('Error in qlimit: abs(b) > 1',11)
          ICKERR=ICKERR+1
        endif
      end if
c**** store tracer flux in fqv array
ccc   fqv(i,:) = fqv(i,:) + f_j(:)
ccc   fqv(i,jm) = 0.   ! play it safe
      hfqv(i,:,l) = f_j(:)
      rmom(ihmoms,i,1 ,l) = 0.! horizontal moments are zero at pole
      rmom(ihmoms,i,jm,l) = 0.
      enddo ! end loop over longitudes
c**** average and unscale polar boxes
      mass(:,1 ,l) = (m_sp + sum(mass(:,1 ,l)-m_sp))*byim
      mass(:,jm,l) = (m_np + sum(mass(:,jm,l)-m_np))*byim
      rm(:,1 ,l) = (rm_sp + sum(rm(:,1 ,l)-rm_sp))*byim
      rm(:,jm,l) = (rm_np + sum(rm(:,jm,l)-rm_np))*byim
      rmom(mz ,:,1 ,l) = (rzm_sp  + sum(rmom(mz ,:,1 ,l)-rzm_sp ))*byim
      rmom(mzz,:,1 ,l) = (rzzm_sp + sum(rmom(mzz,:,1 ,l)-rzzm_sp))*byim
      rmom(mz ,:,jm,l) = (rzm_np  + sum(rmom(mz ,:,jm,l)-rzm_np ))*byim
      rmom(mzz,:,jm,l) = (rzzm_np + sum(rmom(mzz,:,jm,l)-rzzm_np))*byim
      enddo ! end loop over levels
!$OMP  END PARALLEL DO
c
c     sum into fqv
c
!$OMP  PARALLEL DO PRIVATE(J,L)
      do j=1,jm
      do l=1,lm
         fqv(:,j)  = fqv(:,j) + hfqv(:,j,l)
      enddo ! j
      enddo ! l
!$OMP  END PARALLEL DO
      fqv(:,jm) = 0. ! not really needed
C
      IF(ICKERR.GT.0)  CALL stop_model('Stopped in aadvty',11)
C
      return
c****
      end subroutine aadvty



      subroutine aadvtz(rm,rmom,mass,mw,qlimit) 1,4
!@sum  AADVTZ advection driver for z-direction
!@auth Maxwell Kelley
c****
c**** aadvtz advects tracers in the upward vertical direction using the
c**** quadratic upstream scheme.  if qlimit is true, the moments are
c**** limited to prevent the mean tracer from becoming negative.
c****
c**** input:
c****     mw (kg) = vertical mass flux, positive upward
c****      qlimit = whether moment limitations should be used
c****
c**** input/output:
c****     rm (kg) = tracer mass
c****   rmom (kg) = moments of tracer mass
c****   mass (kg) = fluid mass
c****
      use QUSDEF
ccc   use QUSCOM, only : im,jm,lm, zstride,cm,f_l,fmom_l
      use QUSCOM, only : im,jm,lm, zstride
      implicit none
      REAL*8, dimension(im,jm,lm) :: rm,mass,mw
      REAL*8, dimension(NMOM,IM,JM,LM) :: rmom
      logical ::  qlimit
      REAL*8  CM(LM),F_L(LM),FMOM_L(NMOM,LM),MASS_L(LM)
      real*8 bynstep,courmax
      integer :: i,j,l,ierr,nerr,ICKERR,ns,nstep
c**** loop over latitudes and longitudes
      ICKERR=0
!$OMP  PARALLEL DO PRIVATE(I,J,CM,F_L,FMOM_L,IERR,NERR,
!$OMP* NSTEP,COURMAX,BYNSTEP,MASS_L,NS,L)
!$OMP* SHARED(LM,QLIMIT,ZSTRIDE)
!$OMP* REDUCTION(+:ICKERR)
      do j=1,jm
      do i=1,im
c****
c**** decide how many timesteps to take
c****
      nstep=0
      courmax = 2.
      do while(courmax.gt.1. .and. nstep.lt.20)
        nstep = nstep+1
        bynstep = 1d0/real(nstep,kind=8)
        cm(:) = mw(i,j,:)*bynstep
        cm(lm) = 0.
        mass_l(:)  = mass(i,j,:)
        courmax = 0.
        do ns=1,nstep
          do l=1,lm-1
            if(cm(l).gt.0.) then
               courmax = max(courmax,+cm(l)/mass_l(l))
            else
               courmax = max(courmax,-cm(l)/mass_l(l+1))
            endif
          enddo
          if(ns.lt.nstep) then
            do l=1,lm-1
               mass_l(l  ) = mass_l(l  ) - cm(l)
               mass_l(l+1) = mass_l(l+1) + cm(l)
            enddo
          endif
        enddo
      enddo
      if(courmax.gt.1.) then
         write(6,*) 'aadvtz: i,j,courmax=',i,j,courmax
         ICKERR=ICKERR+1
      endif

c      cm(:) = mw(i,j,:)*bynstep ! cm already set
c      cm(lm)= 0.

c****
c**** loop over timesteps
c****
      do ns=1,nstep
c****
c**** call 1-d advection routine
c****
      call adv1d(rm(i,j,1),rmom(1,i,j,1),f_l,fmom_l,mass(i,j,1),
     &        cm,lm,qlimit,zstride,zdir,ierr,nerr)
      if (ierr.gt.0) then
        write(6,*) "Error in aadvtz: i,j,l=",i,j,nerr
        if (ierr.eq.2) then
          write(0,*) "Error in qlimit: abs(c) > 1"
ccc       call stop_model('Error in qlimit: abs(c) > 1',11)
          ICKERR=ICKERR+1
        endif
      end if
      enddo ! ns
      enddo ! i
      enddo ! j
!$OMP  END PARALLEL DO
C
      IF(ICKERR.GT.0) call stop_model('Stopped in aadvtz',11)
      return
c****
      end subroutine aadvtz