map3_ppm.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

!----------------------------------------------------------------------- 
!BOP
! !ROUTINE:  map3_ppm --- Piecewise parabolic mapping, variant 3
!
! !INTERFACE:
    subroutine map3_ppm( km,   pe1,      q1,          &
                         kn,   pe2,      q2,   dp2,   &
                         itot, i1, i2, iv, kord )

! !USES:
    use precision
    implicit none

! !INPUT PARAMETERS:
    integer i1                        ! Starting longitude
    integer i2                        ! Finishing longitude
    integer itot                      ! Total latitudes
    integer iv                        ! Mode: 0 ==  constituents  1 == ???
    integer kord                      ! Method order
    integer km                        ! Original vertical dimension
    integer kn                        ! Target vertical dimension

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

! !OUTPUT PARAMETERS:
    real(r8) q2(itot,kn)              ! Field output

! !DESCRIPTION:
!
!     Perform piecewise parabolic method on a given latitude.  
!
! !REVISION HISTORY: 
!    00.04.24   Lin       Last modification
!    01.03.26   Sawyer    Added ProTeX documentation
!
!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,k)
       enddo
    enddo

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

! Mapping
    do 1000 i=i1,i2
         k0 = 1
      do 555 k=1,kn
      do 100 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,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
               goto 555
          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
                 goto 123
                 endif
              enddo
              goto 123
           endif
      endif
100   continue
123   q2(i,k) = qsum / dp2(i,k)
555   continue
1000  continue

      return
!EOC
      end subroutine map3_ppm
!-----------------------------------------------------------------------