map1_ppm.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

!----------------------------------------------------------------------- 
!BOP
! !ROUTINE:  map1_ppm --- Piecewise parabolic mapping, variant 1
!
! !INTERFACE:
      subroutine map1_ppm( km,     pe1,   q1,    kn,   pe2,               &
                           q2,     itot,  i1,    i2,   j,                 &
                           jfirst, jlast, ngs,   ngn,  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 j                              ! Current latitude
      integer jfirst                         ! Starting latitude
      integer jlast                          ! Finishing latitude
      integer ngs                            ! Ghosted latitudes south
      integer ngn                            ! Ghosted latitudes north
      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)  q1(itot,jfirst-ngs:jlast+ngn,km)     ! Field input

! !OUTPUT PARAMETERS:
      real(r8)  q2(itot,jfirst-ngs:jlast+ngn,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,j,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) .ge. pe1(i,l) .and. pe2(i,k) .le. pe1(i,l+1)) then
         pl = (pe2(i,k)-pe1(i,l)) / dp1(i,l)
         if(pe2(i,k+1) .le. 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
               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) .gt. 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,j,k) = qsum / ( pe2(i,k+1) - pe2(i,k) )
555   continue
1000  continue

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