module mapz_module 1

 public map1_ppm, mapn_ppm

CONTAINS

!----------------------------------------------------------------------- 
!BOP
! !ROUTINE:  map1_ppm --- Piecewise parabolic mapping, variant 1
!
! !INTERFACE:

  subroutine map1_ppm( km,   pe1,   q1,  kn,   pe2,   q2,                & 3,2
                       ng_s, ng_n, itot, i1, i2,                         &
                       j, jfirst, jlast, iv, kord)

! !USES:
  use shr_kind_mod, only: r8 => shr_kind_r8
      implicit none

! !INPUT PARAMETERS:
      integer, intent(in) :: i1                ! Starting longitude
      integer, intent(in) :: i2                ! Finishing longitude
      integer, intent(in) :: itot              ! Total latitudes
      integer, intent(in) :: iv                ! Mode: 0 ==  constituents  1 == ???
      integer, intent(in) :: kord              ! Method order
      integer, intent(in) :: j                 ! Current latitude
      integer, intent(in) :: jfirst            ! Starting latitude
      integer, intent(in) :: jlast             ! Finishing latitude
      integer, intent(in) :: ng_s              ! Ghosted latitudes south
      integer, intent(in) :: ng_n              ! Ghosted latitudes north
      integer, intent(in) :: km                ! Original vertical dimension
      integer, intent(in) :: kn                ! Target vertical dimension

      real(r8), intent(in) ::  pe1(itot,km+1)  ! pressure at layer edges 
                                               ! (from model top to bottom surface)
                                               ! in the original vertical coordinate
      real(r8), intent(in) ::  pe2(itot,kn+1)  ! pressure at layer edges 
                                               ! (from model top to bottom surface)
                                               ! in the new vertical coordinate
      real(r8), intent(in) ::  q1(itot,jfirst-ng_s:jlast+ng_n,km) ! Field input

! !INPUT/OUTPUT PARAMETERS:
      real(r8), intent(inout)::  q2(itot,jfirst-ng_s:jlast+ng_n,kn) ! Field output

! !DESCRIPTION:
!
!     Perform piecewise parabolic method on a given latitude    
! IV = 0: constituents
! pe1: pressure at layer edges (from model top to bottom surface)
!      in the original vertical coordinate
! pe2: pressure at layer edges (from model top to bottom surface)
!      in the new vertical coordinate
!
! !REVISION HISTORY: 
!    00.04.24   Lin       Last modification
!    01.03.26   Sawyer    Added ProTeX documentation
!    02.04.04   Sawyer    incorporated latest FVGCM version
!    02.06.20   Sawyer    made Q2 inout since the args for Q1/Q2 same
!    03.07.22   Parks     Cleaned main loop, removed gotos
!
!EOP
!-----------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
      real(r8)       r3, r23
      parameter (r3 = 1./3., r23 = 2./3.)
      real(r8)   dp1(i1:i2,km)
      real(r8)  q4(4,i1:i2,km)

      integer i, k, l, ll, k0
      real(r8)    pl, pr, qsum, delp, esl

      do k=1,km
         do i=i1,i2
             dp1(i,k) = pe1(i,k+1) - pe1(i,k)
            q4(1,i,k) = q1(i,j,k)
         enddo
      enddo

! Compute vertical subgrid distribution
      call ppm2m( q4, dp1, km, i1, i2, iv, kord )

! Mapping
      do i=i1,i2

        k0 = 1
        do k=1,kn

          do l=k0,km

! locate the top edge: pe2(i,k)
            if(pe2(i,k) >= pe1(i,l) .and. pe2(i,k) <= pe1(i,l+1)) then
              pl = (pe2(i,k)-pe1(i,l)) / dp1(i,l)
              if(pe2(i,k+1) <= pe1(i,l+1)) then
! entire new grid is within the original grid
                pr = (pe2(i,k+1)-pe1(i,l)) / dp1(i,l)
                q2(i,j,k) = q4(2,i,l) + 0.5*(q4(4,i,l)+q4(3,i,l)-q4(2,i,l))  &
                       *(pr+pl)-q4(4,i,l)*r3*(pr*(pr+pl)+pl**2)
                k0 = l
                exit

              else

! Fractional area...
                qsum = (pe1(i,l+1)-pe2(i,k))*(q4(2,i,l)+0.5*(q4(4,i,l)+   &
                    q4(3,i,l)-q4(2,i,l))*(1.+pl)-q4(4,i,l)*           &
                     (r3*(1.+pl*(1.+pl))))

                do ll=l+1,km
! locate the bottom edge: pe2(i,k+1)
                  if(pe2(i,k+1) > pe1(i,ll+1) ) then
! Whole layer..
                    qsum = qsum + dp1(i,ll)*q4(1,i,ll)
                  else
                    delp = pe2(i,k+1)-pe1(i,ll)
                    esl = delp / dp1(i,ll)
                    qsum = qsum + delp*(q4(2,i,ll)+0.5*esl*               &
                         (q4(3,i,ll)-q4(2,i,ll)+q4(4,i,ll)*(1.-r23*esl)))
                    k0 = ll
                    exit
                  endif
                enddo

                q2(i,j,k) = qsum / ( pe2(i,k+1) - pe2(i,k) )
                exit
              endif
            endif
          end do                ! do l = k0, km

        end do                  ! do k = 1, kn

      end do                    ! do i = i1, i2

      return
!EOC
 end subroutine map1_ppm
!----------------------------------------------------------------------- 

!-----------------------------------------------------------------------
!BOP
! !ROUTINE:  mapn_ppm --- Piecewise parabolic mapping, multiple tracers
!
! !INTERFACE:

 subroutine mapn_ppm(km,   pe1,   q1, nq,                  & 1,2
                     kn,   pe2,   q2, ng_s, ng_n,          &
                     itot, i1, i2, j,                      &
                     jfirst, jlast, iv, kord)

! !USES:
  use shr_kind_mod, only: r8 => shr_kind_r8
      implicit none

! !INPUT PARAMETERS:
      integer, intent(in) :: i1                ! Starting longitude
      integer, intent(in) :: i2                ! Finishing longitude
      integer, intent(in) :: itot              ! Total latitudes
      integer, intent(in) :: iv                ! Mode: 0 ==  constituents  1 == ???
      integer, intent(in) :: kord              ! Method order
      integer, intent(in) :: j                 ! Current latitude
      integer, intent(in) :: jfirst            ! Starting latitude
      integer, intent(in) :: jlast             ! Finishing latitude
      integer, intent(in) :: ng_s              ! Ghosted latitudes south
      integer, intent(in) :: ng_n              ! Ghosted latitudes north
      integer, intent(in) :: km                ! Original vertical dimension
      integer, intent(in) :: kn                ! Target vertical dimension
      integer, intent(in) :: nq                ! Number of tracers

      real(r8), intent(in) :: pe1(itot,km+1)   ! pressure at layer edges 
                                               ! (from model top to bottom surface)
                                               ! in the original vertical coordinate
      real(r8), intent(in) :: pe2(itot,kn+1)   ! pressure at layer edges 
                                               ! (from model top to bottom surface)
                                               ! in the new vertical coordinate
      real(r8), intent(in) ::  q1(itot,jfirst-ng_s:jlast+ng_n,km,nq) ! Field input
! !INPUT/OUTPUT PARAMETERS:
      real(r8), intent(inout)::  q2(itot,jfirst-ng_s:jlast+ng_n,kn,nq) ! Field output

! !DESCRIPTION:
!
!     Perform piecewise parabolic method on a given latitude    
! IV = 0: constituents
! pe1: pressure at layer edges (from model top to bottom surface)
!      in the original vertical coordinate
! pe2: pressure at layer edges (from model top to bottom surface)
!      in the new vertical coordinate
!
! !REVISION HISTORY: 
!    02.04.04   Sawyer    incorporated latest FVGCM version, ProTeX
!    02.06.20   Sawyer    made Q2 inout since the args for Q1/Q2 same
!    03.07.22   Sawyer    Cleaned main loop, removed gotos (from Parks)
!
!EOP
!-----------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:

      real(r8)       r3, r23
      parameter (r3 = 1./3., r23 = 2./3.)

      real(r8)   dp1(i1:i2,km)
      real(r8)  q4(4,i1:i2,km)

      integer i, k, l, ll, k0, iq
      real(r8)    pl, pr, qsum, delp, esl

      do k=1,km
         do i=i1,i2
             dp1(i,k) = pe1(i,k+1) - pe1(i,k)
         enddo
      enddo

    do iq=1,nq
      do k=1,km
         do i=i1,i2
            q4(1,i,k) = q1(i,j,k,iq)
         enddo
      enddo

! Compute vertical subgrid distribution
      call ppm2m( q4, dp1, km, i1, i2, iv, kord )

! Mapping
      do i=i1,i2
        k0 = 1
        do k=1,kn
          do l=k0,km
! locate the top edge: pe2(i,k)
            if(pe2(i,k) >= pe1(i,l) .and. pe2(i,k) <= pe1(i,l+1)) then
              pl = (pe2(i,k)-pe1(i,l)) / dp1(i,l)
              if(pe2(i,k+1) <= pe1(i,l+1)) then
! entire new grid is within the original grid
                pr = (pe2(i,k+1)-pe1(i,l)) / dp1(i,l)
                q2(i,j,k,iq) = q4(2,i,l) + 0.5*(q4(4,i,l)+q4(3,i,l)-q4(2,i,l))  &
                          *(pr+pl)-q4(4,i,l)*r3*(pr*(pr+pl)+pl**2)
                k0 = l
                exit
              else
! Fractional area...
                qsum = (pe1(i,l+1)-pe2(i,k))*(q4(2,i,l)+0.5*(q4(4,i,l)+     &
                    q4(3,i,l)-q4(2,i,l))*(1.+pl)-q4(4,i,l)*             &
                     (r3*(1.+pl*(1.+pl))))
                do ll=l+1,km
! locate the bottom edge: pe2(i,k+1)
                  if(pe2(i,k+1) > pe1(i,ll+1) ) then
! Whole layer..
                    qsum = qsum + dp1(i,ll)*q4(1,i,ll)
                  else
                    delp = pe2(i,k+1)-pe1(i,ll)
                    esl = delp / dp1(i,ll)
                    qsum = qsum + delp*(q4(2,i,ll)+0.5*esl*               &
                       (q4(3,i,ll)-q4(2,i,ll)+q4(4,i,ll)*(1.-r23*esl)))
                    k0 = ll
                    exit
                  endif
                enddo

                q2(i,j,k,iq) = qsum / ( pe2(i,k+1) - pe2(i,k) )
                exit
              endif
            endif

          enddo   ! do l=k0,km

        enddo   ! do k=1,kn

      enddo   ! do i=i1,i2

    enddo   ! do iq=1,nq

    return
!EOC
 end subroutine mapn_ppm
!----------------------------------------------------------------------- 


!----------------------------------------------------------------------- 
!BOP
! !ROUTINE:  ppm2m --- Piecewise parabolic method for fields
!
! !INTERFACE:

 subroutine ppm2m(a4, delp, km, i1, i2, iv, kord) 2,6

! !USES:
 use shr_kind_mod, only: r8 => shr_kind_r8
 implicit none

! !INPUT PARAMETERS:
 integer, intent(in):: iv      ! iv =-1: winds
                               ! iv = 0: positive definite scalars
                               ! iv = 1: others
 integer, intent(in):: i1      ! Starting longitude
 integer, intent(in):: i2      ! Finishing longitude
 integer, intent(in):: km      ! vertical dimension
 integer, intent(in):: kord    ! Order (or more accurately method no.):
                               ! 
 real (r8), intent(in):: delp(i1:i2,km)     ! layer pressure thickness

! !INPUT/OUTPUT PARAMETERS:
 real (r8), intent(inout):: a4(4,i1:i2,km)  ! Interpolated values

! !DESCRIPTION:
!
!   Perform the piecewise parabolic method 
! 
! !REVISION HISTORY: 
!   ??.??.??    Lin        Creation
!   02.04.04    Sawyer     Newest release from FVGCM
!   02.04.23    Sawyer     Incorporated minor algorithmic change to 
!                          maintain CAM zero diffs (see comments inline)
! 
!EOP
!-----------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
! local arrays:
      real(r8)  dc(i1:i2,km)
      real(r8)  h2(i1:i2,km)
      real(r8) delq(i1:i2,km)
      real(r8) df2(i1:i2,km)
      real(r8) d4(i1:i2,km)

! local scalars:
      real(r8) fac
      real(r8) a1, a2, c1, c2, c3, d1, d2
      real(r8) qmax, qmin, cmax, cmin
      real(r8) qm, dq, tmp

      integer i, k, km1, lmt
      real(r8) qmp, pmp
      real(r8) lac
      integer it

      km1 = km - 1
       it = i2 - i1 + 1

      do k=2,km
         do i=i1,i2
            delq(i,k-1) =   a4(1,i,k) - a4(1,i,k-1)
              d4(i,k  ) = delp(i,k-1) + delp(i,k)
         enddo
      enddo

      do k=2,km1
         do i=i1,i2
            c1  = (delp(i,k-1)+0.5*delp(i,k))/d4(i,k+1)
            c2  = (delp(i,k+1)+0.5*delp(i,k))/d4(i,k)
            tmp = delp(i,k)*(c1*delq(i,k) + c2*delq(i,k-1)) /      &
                                    (d4(i,k)+delp(i,k+1))
            qmax = max(a4(1,i,k-1),a4(1,i,k),a4(1,i,k+1)) - a4(1,i,k)
            qmin = a4(1,i,k) - min(a4(1,i,k-1),a4(1,i,k),a4(1,i,k+1))
             dc(i,k) = sign(min(abs(tmp),qmax,qmin), tmp)
            df2(i,k) = tmp
         enddo
      enddo

!****6***0*********0*********0*********0*********0*********0**********72
! 4th order interpolation of the provisional cell edge value
!****6***0*********0*********0*********0*********0*********0**********72

      do k=3,km1
      do i=i1,i2
      c1 = delq(i,k-1)*delp(i,k-1) / d4(i,k)
      a1 = d4(i,k-1) / (d4(i,k) + delp(i,k-1))
      a2 = d4(i,k+1) / (d4(i,k) + delp(i,k))
      a4(2,i,k) = a4(1,i,k-1) + c1 + 2./(d4(i,k-1)+d4(i,k+1)) *    &
                ( delp(i,k)*(c1*(a1 - a2)+a2*dc(i,k-1)) -          &
                                delp(i,k-1)*a1*dc(i,k  ) )
      enddo
      enddo

      call steepz(i1, i2, km, a4, df2, dc, delq, delp, d4)

! Area preserving cubic with 2nd deriv. = 0 at the boundaries
! Top
      do i=i1,i2
      d1 = delp(i,1)
      d2 = delp(i,2)
      qm = (d2*a4(1,i,1)+d1*a4(1,i,2)) / (d1+d2)
      dq = 2.*(a4(1,i,2)-a4(1,i,1)) / (d1+d2)
      c1 = 4.*(a4(2,i,3)-qm-d2*dq) / ( d2*(2.*d2*d2+d1*(d2+3.*d1)) )
      c3 = dq - 0.5*c1*(d2*(5.*d1+d2)-3.*d1**2)
      a4(2,i,2) = qm - 0.25*c1*d1*d2*(d2+3.*d1)
      a4(2,i,1) = d1*(2.*c1*d1**2-c3) + a4(2,i,2)
      dc(i,1) =  a4(1,i,1) - a4(2,i,1)
! No over- and undershoot condition
      cmax = max(a4(1,i,1), a4(1,i,2))
      cmin = min(a4(1,i,1), a4(1,i,2))
      a4(2,i,2) = max(cmin,a4(2,i,2))
      a4(2,i,2) = min(cmax,a4(2,i,2))
      enddo

      if( iv == 0 ) then
         do i=i1,i2
!
! WS: 02.04.23  Algorithmic difference with FVGCM.  FVGCM does this:
!
!!!            a4(2,i,1) = a4(1,i,1)
!!!            a4(3,i,1) = a4(1,i,1)
!
!     CAM does this:
!
            a4(2,i,1) = max(0.,a4(2,i,1))
            a4(2,i,2) = max(0.,a4(2,i,2))
         enddo
      elseif ( iv == -1 ) then
! Winds:
        if( km > 32 ) then
          do i=i1,i2
! More dampping: top layer as the sponge
             a4(2,i,1) = a4(1,i,1)
             a4(3,i,1) = a4(1,i,1)
          enddo
        else
          do i=i1,i2
             if( a4(1,i,1)*a4(2,i,1) <=  0. ) then
                 a4(2,i,1) = 0.
             else
                 a4(2,i,1) = sign(min(abs(a4(1,i,1)),    &
                                      abs(a4(2,i,1))),   &
                                          a4(1,i,1)  )
            endif
          enddo
        endif
      endif

! Bottom
! Area preserving cubic with 2nd deriv. = 0 at the surface
      do i=i1,i2
         d1 = delp(i,km)
         d2 = delp(i,km1)
         qm = (d2*a4(1,i,km)+d1*a4(1,i,km1)) / (d1+d2)
         dq = 2.*(a4(1,i,km1)-a4(1,i,km)) / (d1+d2)
         c1 = (a4(2,i,km1)-qm-d2*dq) / (d2*(2.*d2*d2+d1*(d2+3.*d1)))
         c3 = dq - 2.0*c1*(d2*(5.*d1+d2)-3.*d1**2)
         a4(2,i,km) = qm - c1*d1*d2*(d2+3.*d1)
         a4(3,i,km) = d1*(8.*c1*d1**2-c3) + a4(2,i,km)
         dc(i,km) = a4(3,i,km) -  a4(1,i,km)
! No over- and under-shoot condition
         cmax = max(a4(1,i,km), a4(1,i,km1))
         cmin = min(a4(1,i,km), a4(1,i,km1))
         a4(2,i,km) = max(cmin,a4(2,i,km))
         a4(2,i,km) = min(cmax,a4(2,i,km))
      enddo

! Enforce constraint at the surface

      if ( iv == 0 ) then
! Positive definite scalars:
           do i=i1,i2
              a4(3,i,km) = max(0., a4(3,i,km))
           enddo
      elseif ( iv == -1 ) then
! Winds:
           do i=i1,i2
              if( a4(1,i,km)*a4(3,i,km) <=  0. ) then
                  a4(3,i,km) = 0.
              else
                  a4(3,i,km) = sign( min(abs(a4(1,i,km)),   &
                                         abs(a4(3,i,km))),  &
                                             a4(1,i,km)  )
              endif
           enddo
      endif

      do k=1,km1
         do i=i1,i2
            a4(3,i,k) = a4(2,i,k+1)
         enddo
      enddo
 
! f(s) = AL + s*[(AR-AL) + A6*(1-s)]         ( 0 <= s  <= 1 )
 
! Top 2 and bottom 2 layers always use monotonic mapping
      do k=1,2
         do i=i1,i2
            a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
         enddo
         call kmppm(dc(i1,k), a4(1,i1,k), it, 0)
      enddo

      if(kord >= 7) then
!****6***0*********0*********0*********0*********0*********0**********72
! Huynh's 2nd constraint
!****6***0*********0*********0*********0*********0*********0**********72
      do k=2, km1
         do i=i1,i2
! Method#1
!           h2(i,k) = delq(i,k) - delq(i,k-1)
! Method#2
!           h2(i,k) = 2.*(dc(i,k+1)/delp(i,k+1) - dc(i,k-1)/delp(i,k-1))
!    &               / ( delp(i,k)+0.5*(delp(i,k-1)+delp(i,k+1)) )
!    &               * delp(i,k)**2
! Method#3
            h2(i,k) = dc(i,k+1) - dc(i,k-1)
         enddo
      enddo

      if( kord == 7 ) then
         fac = 1.5           ! original quasi-monotone
      else
         fac = 0.125         ! full monotone
      endif

      do k=3, km-2
        do i=i1,i2
! Right edges
!        qmp   = a4(1,i,k) + 2.0*delq(i,k-1)
!        lac   = a4(1,i,k) + fac*h2(i,k-1) + 0.5*delq(i,k-1)
!
         pmp   = 2.*dc(i,k)
         qmp   = a4(1,i,k) + pmp
         lac   = a4(1,i,k) + fac*h2(i,k-1) + dc(i,k)
         qmin  = min(a4(1,i,k), qmp, lac)
         qmax  = max(a4(1,i,k), qmp, lac)
         a4(3,i,k) = min(max(a4(3,i,k), qmin), qmax)
! Left  edges
!        qmp   = a4(1,i,k) - 2.0*delq(i,k)
!        lac   = a4(1,i,k) + fac*h2(i,k+1) - 0.5*delq(i,k)
!
         qmp   = a4(1,i,k) - pmp
         lac   = a4(1,i,k) + fac*h2(i,k+1) - dc(i,k)
         qmin  = min(a4(1,i,k), qmp, lac)
         qmax  = max(a4(1,i,k), qmp, lac)
         a4(2,i,k) = min(max(a4(2,i,k), qmin), qmax)
! Recompute A6
         a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
        enddo
! Additional constraint to prevent negatives when kord=7
         if (iv == 0 .and. kord == 7) then
             call kmppm(dc(i1,k), a4(1,i1,k), it, 2)
         endif
      enddo

      else
 
         lmt = kord - 3
         lmt = max(0, lmt)
         if (iv == 0) lmt = min(2, lmt)

      do k=3, km-2
      if( kord /= 4) then
         do i=i1,i2
            a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
         enddo
      endif
         call kmppm(dc(i1,k), a4(1,i1,k), it, lmt)
      enddo
      endif

      do k=km1,km
         do i=i1,i2
            a4(4,i,k) = 3.*(2.*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
         enddo
         call kmppm(dc(i1,k), a4(1,i1,k), it, 0)
      enddo

      return
!EOC
 end subroutine ppm2m
!-----------------------------------------------------------------------

!----------------------------------------------------------------------- 
!BOP
! !ROUTINE:  kmppm --- Perform piecewise parabolic method in vertical
!
! !INTERFACE:

 subroutine kmppm(dm, a4, itot, lmt) 4,1

! !USES:
 use shr_kind_mod, only: r8 => shr_kind_r8
 implicit none

! !INPUT PARAMETERS:
      real(r8), intent(in)::     dm(*)     ! ??????
      integer, intent(in) ::     itot      ! Total Longitudes
      integer, intent(in) ::     lmt       ! 0: Standard PPM constraint
                                           ! 1: Improved full monotonicity constraint (Lin)
                                           ! 2: Positive definite constraint
                                           ! 3: do nothing (return immediately)

! !INPUT/OUTPUT PARAMETERS:
      real(r8), intent(inout) :: a4(4,*)   ! ???????
                                           ! AA <-- a4(1,i)
                                           ! AL <-- a4(2,i)
                                           ! AR <-- a4(3,i)
                                           ! A6 <-- a4(4,i)

! !DESCRIPTION:
!
!    Writes a standard set of data to the history buffer. 
!
! !REVISION HISTORY: 
!    00.04.24   Lin       Last modification
!    01.03.26   Sawyer    Added ProTeX documentation
!    02.04.04   Sawyer    Incorporated newest FVGCM version
!
!EOP
!-----------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:

      real(r8)       r12
      parameter (r12 = 1./12.)

      real(r8) qmp
      integer i
      real(r8) da1, da2, a6da
      real(r8) fmin

! Developer: S.-J. Lin, NASA-GSFC
! Last modified: Apr 24, 2000

      if ( lmt == 3 ) return

      if(lmt == 0) then
! Standard PPM constraint
      do i=1,itot
      if(dm(i) == 0.) then
         a4(2,i) = a4(1,i)
         a4(3,i) = a4(1,i)
         a4(4,i) = 0.
      else
         da1  = a4(3,i) - a4(2,i)
         da2  = da1**2
         a6da = a4(4,i)*da1
         if(a6da < -da2) then
            a4(4,i) = 3.*(a4(2,i)-a4(1,i))
            a4(3,i) = a4(2,i) - a4(4,i)
         elseif(a6da > da2) then
            a4(4,i) = 3.*(a4(3,i)-a4(1,i))
            a4(2,i) = a4(3,i) - a4(4,i)
         endif
      endif
      enddo

      elseif (lmt == 1) then

! Improved full monotonicity constraint (Lin)
! Note: no need to provide first guess of A6 <-- a4(4,i)
      do i=1, itot
           qmp = 2.*dm(i)
         a4(2,i) = a4(1,i)-sign(min(abs(qmp),abs(a4(2,i)-a4(1,i))), qmp)
         a4(3,i) = a4(1,i)+sign(min(abs(qmp),abs(a4(3,i)-a4(1,i))), qmp)
         a4(4,i) = 3.*( 2.*a4(1,i) - (a4(2,i)+a4(3,i)) )
      enddo

      elseif (lmt == 2) then

! Positive definite constraint
      do i=1,itot
      if( abs(a4(3,i)-a4(2,i)) < -a4(4,i) ) then
      fmin = a4(1,i)+0.25*(a4(3,i)-a4(2,i))**2/a4(4,i)+a4(4,i)*r12
         if( fmin < 0. ) then
         if(a4(1,i)<a4(3,i) .and. a4(1,i)<a4(2,i)) then
            a4(3,i) = a4(1,i)
            a4(2,i) = a4(1,i)
            a4(4,i) = 0.
         elseif(a4(3,i) > a4(2,i)) then
            a4(4,i) = 3.*(a4(2,i)-a4(1,i))
            a4(3,i) = a4(2,i) - a4(4,i)
         else
            a4(4,i) = 3.*(a4(3,i)-a4(1,i))
            a4(2,i) = a4(3,i) - a4(4,i)
         endif
         endif
      endif
      enddo

      endif

      return
!EOC
 end subroutine kmppm
!-----------------------------------------------------------------------

!----------------------------------------------------------------------- 
!BOP
! !ROUTINE:  steepz --- Calculate attributes for PPM
!
! !INTERFACE:

 subroutine steepz(i1, i2, km, a4, df2, dm, dq, dp, d4) 1,1

! !USES:
   use shr_kind_mod, only: r8 => shr_kind_r8
   implicit none

! !INPUT PARAMETERS:
      integer, intent(in) :: km                   ! Total levels
      integer, intent(in) :: i1                   ! Starting longitude
      integer, intent(in) :: i2                   ! Finishing longitude
      real(r8), intent(in) ::  dp(i1:i2,km)       ! grid size
      real(r8), intent(in) ::  dq(i1:i2,km)       ! backward diff of q
      real(r8), intent(in) ::  d4(i1:i2,km)       ! backward sum:  dp(k)+ dp(k-1) 
      real(r8), intent(in) :: df2(i1:i2,km)       ! first guess mismatch
      real(r8), intent(in) ::  dm(i1:i2,km)       ! monotonic mismatch

! !INPUT/OUTPUT PARAMETERS:
      real(r8), intent(inout) ::  a4(4,i1:i2,km)  ! first guess/steepened

!
! !DESCRIPTION:
!   This is complicated stuff related to the Piecewise Parabolic Method
!   and I need to read the Collela/Woodward paper before documenting
!   thoroughly.
!
! !REVISION HISTORY: 
!   ??.??.??    Lin?       Creation
!   01.03.26    Sawyer     Added ProTeX documentation
!
!EOP
!-----------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
      integer i, k
      real(r8) alfa(i1:i2,km)
      real(r8)    f(i1:i2,km)
      real(r8)  rat(i1:i2,km)
      real(r8)  dg2

! Compute ratio of dq/dp
      do k=2,km
         do i=i1,i2
            rat(i,k) = dq(i,k-1) / d4(i,k)
         enddo
      enddo

! Compute F
      do k=2,km-1
         do i=i1,i2
            f(i,k) = (rat(i,k+1) - rat(i,k))                             &
                     / ( dp(i,k-1)+dp(i,k)+dp(i,k+1) )
         enddo
      enddo

      do k=3,km-2
         do i=i1,i2
         if(f(i,k+1)*f(i,k-1)<0. .and. df2(i,k)/=0.) then
            dg2 = (f(i,k+1)-f(i,k-1))*((dp(i,k+1)-dp(i,k-1))**2          &
                   + d4(i,k)*d4(i,k+1) )
            alfa(i,k) = max(0., min(0.5, -0.1875*dg2/df2(i,k))) 
         else
            alfa(i,k) = 0.
         endif
         enddo
      enddo

      do k=4,km-2
         do i=i1,i2
            a4(2,i,k) = (1.-alfa(i,k-1)-alfa(i,k)) * a4(2,i,k) +         &
                        alfa(i,k-1)*(a4(1,i,k)-dm(i,k))    +             &
                        alfa(i,k)*(a4(1,i,k-1)+dm(i,k-1))
         enddo
      enddo

      return
!EOC
 end subroutine steepz
!----------------------------------------------------------------------- 

end module mapz_module