pkg/generic_advdiff/gad_ppm_fun.F

8e4c181d69 Jean*0001 #     include "GAD_OPTIONS.h"
                0002
                0003 C--  File gad_ppm_fun.F: Routines to form PPM grid-cell polynomial.
                0004 C--   Contents
                0005 C--   o GAD_PPM_FUN_NULL
                0006 C--   o GAD_PPM_FUN_MONO
                0007
                0008 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0009
                0010       SUBROUTINE GAD_PPM_FUN_NULL(
                0011      I           ff00,
                0012      I           fell,ferr,
                0013      O           fhat,mono)
                0014 C     |================================================================|
                0015 C     | PPM_FUN_NULL: form PPM grid-cell polynomial.                   |
                0016 C     | Piecewise Parabolic Method (PPM), unlimited variant.           |
                0017 C     |================================================================|
                0018
                0019           implicit none
                0020
                0021 C     ====================================================== arguments
                0022           _RL ff00
                0023           _RL fell,ferr
                0024           _RL fhat(+1:+3)
                0025           integer  mono
                0026
                0027           mono = +0
                0028
                0029 C     ============================================== unlimited profile
                0030           fhat( 1 ) =
                0031      &  +(3. _d 0 / 2. _d 0) * ff00
                0032      &  -(1. _d 0 / 4. _d 0) *(ferr+fell)
                0033           fhat( 2 ) =
                0034      &  +(1. _d 0 / 2. _d 0) *(ferr-fell)
                0035           fhat( 3 ) =
                0036      &  -(3. _d 0 / 2. _d 0) * ff00
                0037      &  +(3. _d 0 / 4. _d 0) *(ferr+fell)
                0038
                0039           return
                0040
                0041 c     end subroutine GAD_PPM_FUN_NULL
                0042       end
                0043
                0044 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0045
                0046       SUBROUTINE GAD_PPM_FUN_MONO(
                0047      I           ff00,
                0048      I           ffll,ffrr,
                0049      I           fell,ferr,
                0050      I           dfds,
                0051      O           fhat,mono)
                0052 C     |================================================================|
                0053 C     | PPM_FUN_MONO: form PPM grid-cell polynomial.                   |
                0054 C     | Piecewise Parabolic Method (PPM) - monotonic variant.          |
                0055 C     |================================================================|
                0056
                0057           implicit none
                0058
                0059 C     ====================================================== arguments
                0060           _RL ff00
                0061           _RL ffll,ffrr
                0062           _RL fell,ferr
                0063           _RL dfds(-1:+1)
                0064           _RL fhat(+1:+3)
                0065           integer mono
                0066
                0067 C     ====================================================== variables
                0068           _RL turn
                0069
                0070           mono = 0
                0071
                0072 C     ============================================== "flatten" extrema
                0073           if((ffrr-ff00) *
                0074      &       (ff00-ffll) .le. 0. _d 0) then
                0075
                0076               mono = +1
                0077
                0078               fhat(1) = ff00
                0079
                0080               fhat(2) = 0. _d 0
                0081               fhat(3) = 0. _d 0
                0082
                0083               return
                0084
                0085           end if
                0086
                0087 C     ============================================== limit edge values
                0088           if((ffll-fell) *
                0089      &       (fell-ff00) .le. 0. _d 0) then
                0090
                0091               mono = +1
                0092
                0093               fell = ff00 - dfds(0)
                0094
                0095           end if
                0096
                0097           if((ffrr-ferr) *
                0098      &       (ferr-ff00) .le. 0. _d 0) then
                0099
                0100               mono = +1
                0101
                0102               ferr = ff00 + dfds(0)
                0103
                0104           end if
                0105
                0106 C     ============================================== limit cell values
                0107           fhat( 1 ) =
                0108      &  +(3. _d 0 / 2. _d 0) * ff00
                0109      &  -(1. _d 0 / 4. _d 0) *(ferr+fell)
                0110           fhat( 2 ) =
                0111      &  +(1. _d 0 / 2. _d 0) *(ferr-fell)
                0112           fhat( 3 ) =
                0113      &  -(3. _d 0 / 2. _d 0) * ff00
                0114      &  +(3. _d 0 / 4. _d 0) *(ferr+fell)
                0115
                0116           if (abs(fhat(3)) .gt.
                0117      &        abs(fhat(2))*.5 _d 0) then
                0118
                0119           turn = -0.5 _d 0 * fhat(2)
                0120      &                     / fhat(3)
                0121
                0122           if ((turn .ge. -1. _d 0)
                0123      &   .and.(turn .le. +0. _d 0)) then
                0124
                0125               mono = +2
                0126
                0127 C     ====================================== push TURN onto lower edge
                0128               ferr = +3. _d 0 * ff00
                0129      &               -2. _d 0 * fell
                0130
                0131           end if
                0132
                0133           if ((turn .gt. +0. _d 0)
                0134      &   .and.(turn .le. +1. _d 0)) then
                0135
                0136               mono = +2
                0137
                0138 C     ====================================== push TURN onto upper edge
                0139               fell = +3. _d 0 * ff00
                0140      &               -2. _d 0 * ferr
                0141
                0142           end if
                0143
                0144           end if
                0145
                0146           if (mono .gt. +1) then
                0147
                0148 C     ====================================== re-calc. coeff. on demand
                0149           fhat( 1 ) =
                0150      &  +(3. _d 0 / 2. _d 0) * ff00
                0151      &  -(1. _d 0 / 4. _d 0) *(ferr+fell)
                0152           fhat( 2 ) =
                0153      &  +(1. _d 0 / 2. _d 0) *(ferr-fell)
                0154           fhat( 3 ) =
                0155      &  -(3. _d 0 / 2. _d 0) * ff00
                0156      &  +(3. _d 0 / 4. _d 0) *(ferr+fell)
                0157
                0158           end if
                0159
                0160           return
                0161
                0162 c     end subroutine GAD_PPM_FUN_MONO
                0163       end