ppm2m.f90

来自「CCSM Research Tools: Community Atmospher」· F90 代码 · 共 271 行

!----------------------------------------------------------------------- 
!BOP
! !ROUTINE:  ppm2m --- Piecewise parabolic method for fields
!
! !INTERFACE:
 subroutine ppm2m(a4, delp, km, i1, i2, iv, kord)

! !USES:
 use precision
 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
! 
!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)

 real (r8) fac
 real (r8) a1, a2, c1, c2, c3, d1, d2
 real (r8) qmax, qmin, cmax, cmin
 real (r8) qm, dq, tmp
 real (r8) qmp, pmp
 real (r8) lac
 integer lmt
 integer i, k, km1
 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

!------------------------------------------------------------
! 4th order interpolation of the provisional cell edge value
!------------------------------------------------------------

    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
            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:
        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

! 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 .ge. 7) then

!----------------------------------------
! Huynh's 2nd constraint
!----------------------------------------

      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 .eq. 0 .and. kord .eq. 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 .ne. 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
!-----------------------------------------------------------------------