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