pkg/generic_advdiff/gad_plm_fun.F

8e4c181d69 Jean*0001 #     include "GAD_OPTIONS.h"
                0002
                0003 C--  File gad_plm_fun.F: Routines for monotone piecewise linear method.
                0004 C--   Contents
                0005 C--   o GAD_PLM_FUN_U
                0006 C--   o GAD_PLM_FUN_V
                0007
                0008 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0009
                0010       SUBROUTINE GAD_PLM_FUN_U(
                0011      I           ffll,ff00,ffrr,
                0012      O           dfds)
                0013 C     |================================================================|
                0014 C     | PLM_FUN_U: monotone piecewise linear method.                   |
                0015 C     |     - uniform grid-spacing variant.                            |
                0016 C     |================================================================|
                0017
                0018           implicit none
                0019
                0020 C     ====================================================== arguments
                0021           _RL ffll,ff00,ffrr
                0022           _RL dfds(-1:+1)
                0023
                0024 C     ====================================================== variables
                0025           _RL fell,ferr,scal
                0026           _RL epsil
                0027           PARAMETER( epsil = 1. _d -16 )
                0028
                0029           dfds(-1) = ff00 - ffll
                0030           dfds(+1) = ffrr - ff00
                0031
                0032           if (dfds(-1) * dfds(+1) .gt. 0.0 _d 0) then
                0033
                0034 C     ======================================= calc. ll//rr edge values
                0035               fell = 0.5 _d 0 * (ffll + ff00)
                0036               ferr = 0.5 _d 0 * (ff00 + ffrr)
                0037
                0038 C     ======================================= calc. centred derivative
                0039               dfds(+0) =
                0040      &               0.5 _d 0 * (ferr - fell)
                0041
                0042 C     ======================================= monotonic slope-limiting
                0043               scal = min(abs(dfds(-1)),
                0044      &                   abs(dfds(+1)))
                0045      &             / max(abs(dfds(+0)), epsil)
                0046 c    &             / max(abs(dfds(+0)), epsilon(ff00))
                0047               scal = min(scal, 1.0 _d 0)
                0048
                0049               dfds(+0) = scal * dfds(+0)
                0050
                0051           else
                0052
                0053 C     ======================================= flatten if local extrema
                0054               dfds(+0) = 0.0 _d 0
                0055
                0056           end if
                0057
                0058           dfds(-1) = 0.5 _d 0 * dfds(-1)
                0059           dfds(+1) = 0.5 _d 0 * dfds(+1)
                0060
                0061           return
                0062
                0063 c     end subroutine GAD_PLM_FUN_U
                0064       end
                0065
                0066 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0067
                0068       SUBROUTINE GAD_PLM_FUN_V(
                0069      I           ffll,ff00,ffrr,
                0070      I           ddll,dd00,ddrr,
                0071      O           dfds)
                0072 C     |================================================================|
                0073 C     | PLM_FUN_V: monotone piecewise linear method.                   |
                0074 C     |     - variable grid-spacing variant.                           |
                0075 C     |================================================================|
                0076
                0077           implicit none
                0078
                0079 C     ====================================================== arguments
                0080           _RL ffll,ff00,ffrr
                0081           _RL ddll,dd00,ddrr
                0082           _RL dfds(-1:+1)
                0083
                0084 C     ====================================================== variables
                0085           _RL fell,ferr,scal
                0086           _RL epsil
                0087           PARAMETER( epsil = 1. _d -16 )
                0088
                0089           dfds(-1) = ff00 - ffll
                0090           dfds(+1) = ffrr - ff00
                0091
                0092           if (dfds(-1) * dfds(+1) .gt. 0.0 _d 0) then
                0093
                0094 C     ======================================= calc. ll//rr edge values
                0095               fell = (dd00 * ffll + ddll * ff00)
                0096      &             / (ddll + dd00)
                0097               ferr = (ddrr * ff00 + dd00 * ffrr)
                0098      &             / (dd00 + ddrr)
                0099
                0100 C     ======================================= calc. centred derivative
                0101               dfds(+0) =
                0102      &               0.5 _d 0 * (ferr - fell)
                0103
                0104 C     ======================================= monotonic slope-limiting
                0105               scal = min(abs(dfds(-1)),
                0106      &                   abs(dfds(+1)))
                0107      &             / max(abs(dfds(+0)), epsil)
                0108 c    &             / max(abs(dfds(+0)), epsilon(ff00))
                0109               scal = min(scal, 1.0 _d 0)
                0110
                0111               dfds(+0) = scal * dfds(+0)
                0112
                0113           else
                0114
                0115 C     ======================================= flatten if local extrema
                0116               dfds(+0) = 0.0 _d 0
                0117
                0118           end if
                0119
                0120 C     == !! check this
                0121           dfds(-1) = dfds(-1) / (ddll + dd00) * dd00
                0122           dfds(+1) = dfds(+1) / (dd00 + ddrr) * dd00
                0123
                0124           return
                0125
                0126 c     end subroutine GAD_PLM_FUN_V
                0127       end