pkg/generic_advdiff/gad_pqm_fun.F

8e4c181d69 Jean*0001 #     include "GAD_OPTIONS.h"
                0002
                0003 C--  File gad_pqm_fun.F: Routines to form PQM grid-cell polynomial.
                0004 C--   Contents
                0005 C--   o QUADROOT
                0006 C--   o GAD_PPM_FUN_NULL
                0007 C--   o GAD_PPM_FUN_MONO
                0008
                0009 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0010
                0011       LOGICAL FUNCTION QUADROOT(aa,bb,cc,xx)
                0012 C     |================================================================|
                0013 C     | QUADROOT: find roots of quadratic ax**2 + bx + c = 0.          |
                0014 C     |================================================================|
                0015
                0016           implicit none
                0017
                0018 C     ====================================================== arguments
                0019           _RL aa,bb,cc
                0020           _RL xx(1:2)
                0021
                0022 C     ====================================================== variables
                0023           _RL sq,a0,b0
                0024
                0025           a0 = abs(aa)
                0026           b0 = abs(bb)
                0027
                0028           sq = bb * bb - 4. _d 0 * aa * cc
                0029
                0030           if (a0 .gt. 0. _d 0) then
                0031
                0032           if (sq .ge. 0. _d 0) then
                0033
                0034               QUADROOT =  .TRUE.
                0035
                0036               sq = sqrt(sq)
                0037
                0038               xx(1) =  - bb + sq
                0039               xx(2) =  - bb - sq
                0040
                0041               aa = 0.5 _d 0 / aa
                0042
                0043               xx(1) = xx(1) * aa
                0044               xx(2) = xx(2) * aa
                0045
                0046           else
                0047
                0048               QUADROOT = .FALSE.
                0049
                0050           end if
                0051
                0052           else
                0053
                0054           if (b0 .gt. 0. _d 0) then
                0055
                0056               QUADROOT =  .TRUE.
                0057
                0058               xx(1) =  - cc / bb
                0059               xx(2) =  - cc / bb
                0060
                0061           else
                0062
                0063               QUADROOT = .FALSE.
                0064
                0065           end if
                0066
                0067           end if
                0068
                0069           return
                0070
                0071 c     end function QUADROOT
                0072       end
                0073
                0074 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0075
                0076       SUBROUTINE GAD_PQM_FUN_NULL(
                0077      I           ff00,
                0078      I           fell,ferr,
                0079      I           dell,derr,
                0080      O           fhat,mono)
                0081 C     |================================================================|
                0082 C     | PQM_FUN_NULL: form PQM grid-cell polynomial.                   |
                0083 C     | Piecewise Quartic Method (PQM) - unlimited variant.            |
                0084 C     |================================================================|
                0085
                0086           implicit none
                0087
                0088 C     ====================================================== arguments
                0089           _RL ff00
                0090           _RL fell,ferr
                0091           _RL dell,derr
                0092           _RL fhat(+1:+5)
                0093           integer mono
                0094
                0095           mono = +0
                0096
                0097 C     ============================================== unlimited profile
                0098           fhat(1) =
                0099      & + (30. _d 0 / 16. _d 0) * ff00
                0100      & - ( 7. _d 0 / 16. _d 0) *(ferr+fell)
                0101      & + ( 1. _d 0 / 16. _d 0) *(derr-dell)
                0102           fhat(2) =
                0103      & + ( 3. _d 0 /  4. _d 0) *(ferr-fell)
                0104      & - ( 1. _d 0 /  4. _d 0) *(derr+dell)
                0105           fhat(3) =
                0106      & - (30. _d 0 /  8. _d 0) * ff00
                0107      & + (15. _d 0 /  8. _d 0) *(ferr+fell)
                0108      & - ( 3. _d 0 /  8. _d 0) *(derr-dell)
                0109           fhat(4) =
                0110      & - ( 1. _d 0 /  4. _d 0) *(ferr-fell-derr-dell)
                0111           fhat(5) =
                0112      & + (30. _d 0 / 16. _d 0) * ff00
                0113      & - (15. _d 0 / 16. _d 0) *(ferr+fell)
                0114      & + ( 5. _d 0 / 16. _d 0) *(derr-dell)
                0115
                0116           return
                0117
                0118 c     end subroutine GAD_PQM_FUN_NULL
                0119       end
                0120
                0121 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0122
                0123       SUBROUTINE GAD_PQM_FUN_MONO(
                0124      I           ff00,
                0125      I           ffll,ffrr,
                0126      I           fell,ferr,
                0127      I           dell,derr,
                0128      I           dfds,
                0129      O           fhat,mono)
                0130 C     |================================================================|
                0131 C     | PQM_FUN_MONO: form PQM grid-cell polynomial.                   |
                0132 C     | Piecewise Quartic Method (PQM) - monotonic variant.            |
                0133 C     |================================================================|
                0134
                0135           implicit none
                0136
                0137 C     ====================================================== arguments
                0138           _RL ff00,ffll,ffrr
                0139           _RL fell,ferr
                0140           _RL dell,derr
                0141           _RL dfds(-1:+1)
                0142           _RL fhat(+1:+5)
                0143           integer mono
                0144
                0145 C     ====================================================== functions
                0146           logical QUADROOT
                0147
                0148 C     ====================================================== variables
                0149           integer bind
                0150           _RL aval,bval,cval
                0151           _RL iflx(+1:+2)
                0152           _RL dflx(+1:+2)
                0153
                0154           mono = +0
                0155
                0156 C     ============================================== "flatten" extrema
                0157           if((ffrr-ff00) *
                0158      &       (ff00-ffll) .le. 0. _d 0) then
                0159
                0160               mono = +1
                0161
                0162               fhat(1) = ff00
                0163
                0164               fhat(2) = 0. _d 0
                0165               fhat(3) = 0. _d 0
                0166               fhat(4) = 0. _d 0
                0167               fhat(5) = 0. _d 0
                0168
                0169               return
                0170
                0171           end if
                0172
                0173 C     ============================================== limit edge values
                0174           if((ffll - fell) *
                0175      &       (fell - ff00) .le. 0. _d 0) then
                0176
                0177               mono = +1
                0178
                0179               fell = ff00 - dfds(0)
                0180
                0181           end if
                0182
                0183           if((ffrr - ferr) *
                0184      &       (ferr - ff00) .le. 0. _d 0) then
                0185
                0186               mono = +1
                0187
                0188               ferr = ff00 + dfds(0)
                0189
                0190           end if
                0191
                0192 C     ============================================== limit edge slopes
                0193           if((dell * dfds(-1)) .lt. 0. _d 0) then
                0194
                0195               mono = +1
                0196
                0197               dell = dfds(-1)
                0198
                0199           end if
                0200
                0201           if((derr * dfds(+1)) .lt. 0. _d 0) then
                0202
                0203               mono = +1
                0204
                0205               derr = dfds(+1)
                0206
                0207           end if
                0208
                0209 C     ============================================== limit cell values
                0210           fhat(1) =
                0211      & + (30. _d 0 / 16. _d 0) * ff00
                0212      & - ( 7. _d 0 / 16. _d 0) *(ferr+fell)
                0213      & + ( 1. _d 0 / 16. _d 0) *(derr-dell)
                0214           fhat(2) =
                0215      & + ( 3. _d 0 /  4. _d 0) *(ferr-fell)
                0216      & - ( 1. _d 0 /  4. _d 0) *(derr+dell)
                0217           fhat(3) =
                0218      & - (30. _d 0 /  8. _d 0) * ff00
                0219      & + (15. _d 0 /  8. _d 0) *(ferr+fell)
                0220      & - ( 3. _d 0 /  8. _d 0) *(derr-dell)
                0221           fhat(4) =
                0222      & - ( 1. _d 0 /  4. _d 0) *(ferr-fell-derr-dell)
                0223           fhat(5) =
                0224      & + (30. _d 0 / 16. _d 0) * ff00
                0225      & - (15. _d 0 / 16. _d 0) *(ferr+fell)
                0226      & + ( 5. _d 0 / 16. _d 0) *(derr-dell)
                0227
                0228 C     ============================= calc. inflexion via 2nd-derivative
                0229           aval = 12. _d 0 * fhat(5)
                0230           bval =  6. _d 0 * fhat(4)
                0231           cval =  2. _d 0 * fhat(3)
                0232
                0233           if ( QUADROOT(aval,bval,cval,iflx) ) then
                0234
                0235               bind = +0
                0236
                0237               if ( ( iflx(1) .gt. -1. _d 0 )
                0238      &       .and. ( iflx(1) .lt. +1. _d 0 ) ) then
                0239
                0240 C     ============================= check for non-monotonic inflection
                0241               dflx(1) =       fhat(2)
                0242      &      + iflx(1)       * fhat(3) * 2. _d 0
                0243      &      +(iflx(1) ** 2) * fhat(4) * 3. _d 0
                0244      &      +(iflx(1) ** 3) * fhat(5) * 4. _d 0
                0245
                0246               if (dflx(1)*dfds(+0) .lt. 0. _d 0) then
                0247
                0248                   if (abs(dell)
                0249      &           .lt. abs(derr) ) then
                0250
                0251                       bind = -1
                0252
                0253                   else
                0254
                0255                       bind = +1
                0256
                0257                   end if
                0258
                0259               end if
                0260
                0261               end if
                0262
                0263               if ( ( iflx(2) .gt. -1. _d 0 )
                0264      &       .and. ( iflx(2) .lt. +1. _d 0 ) ) then
                0265
                0266 C     ============================= check for non-monotonic inflection
                0267               dflx(2) =       fhat(2)
                0268      &      + iflx(2)       * fhat(3) * 2. _d 0
                0269      &      +(iflx(2) ** 2) * fhat(4) * 3. _d 0
                0270      &      +(iflx(2) ** 3) * fhat(5) * 4. _d 0
                0271
                0272               if (dflx(2)*dfds(+0) .lt. 0. _d 0) then
                0273
                0274                   if (abs(dell)
                0275      &           .lt. abs(derr) ) then
                0276
                0277                       bind = -1
                0278
                0279                   else
                0280
                0281                       bind = +1
                0282
                0283                   end if
                0284
                0285               end if
                0286
                0287               end if
                0288
                0289 C     ============================= pop non-monotone inflexion to edge
                0290
                0291               if (bind .eq. -1) then
                0292
                0293 C     ============================= pop inflection points onto -1 edge
                0294                   mono = +2
                0295
                0296                   derr =
                0297      &          -( 5. _d 0 / 1. _d 0) * ff00
                0298      &          +( 3. _d 0 / 1. _d 0) * ferr
                0299      &          +( 2. _d 0 / 1. _d 0) * fell
                0300                   dell =
                0301      &          +( 5. _d 0 / 3. _d 0) * ff00
                0302      &          -( 1. _d 0 / 3. _d 0) * ferr
                0303      &          -( 4. _d 0 / 3. _d 0) * fell
                0304
                0305                   if (dell*dfds(-1) .lt. 0. _d 0) then
                0306
                0307                       dell = 0. _d 0
                0308
                0309                       ferr =
                0310      &              +( 5. _d 0 / 1. _d 0) * ff00
                0311      &              -( 4. _d 0 / 1. _d 0) * fell
                0312                       derr =
                0313      &              +(10. _d 0 / 1. _d 0) * ff00
                0314      &              -(10. _d 0 / 1. _d 0) * fell
                0315
                0316                   end if
                0317
                0318                   if (derr*dfds(+1) .lt. 0. _d 0) then
                0319
                0320                       derr = 0. _d 0
                0321
                0322                       fell =
                0323      &              +( 5. _d 0 / 2. _d 0) * ff00
                0324      &              -( 3. _d 0 / 2. _d 0) * ferr
                0325                       dell =
                0326      &              -( 5. _d 0 / 3. _d 0) * ff00
                0327      &              +( 5. _d 0 / 3. _d 0) * ferr
                0328
                0329                   end if
                0330
                0331               end if
                0332
                0333               if (bind .eq. +1) then
                0334
                0335 C     ============================= pop inflection points onto +1 edge
                0336                   mono = +2
                0337
                0338                   derr =
                0339      &          -( 5. _d 0 / 3. _d 0) * ff00
                0340      &          +( 4. _d 0 / 3. _d 0) * ferr
                0341      &          +( 1. _d 0 / 3. _d 0) * fell
                0342                   dell =
                0343      &          +( 5. _d 0 / 1. _d 0) * ff00
                0344      &          -( 2. _d 0 / 1. _d 0) * ferr
                0345      &          -( 3. _d 0 / 1. _d 0) * fell
                0346
                0347                   if (dell*dfds(-1) .lt. 0. _d 0) then
                0348
                0349                       dell = 0. _d 0
                0350
                0351                       ferr =
                0352      &              +( 5. _d 0 / 2. _d 0) * ff00
                0353      &              -( 3. _d 0 / 2. _d 0) * fell
                0354                       derr =
                0355      &              +( 5. _d 0 / 3. _d 0) * ff00
                0356      &              -( 5. _d 0 / 3. _d 0) * fell
                0357
                0358                   end if
                0359
                0360                   if (derr*dfds(+1) .lt. 0. _d 0) then
                0361
                0362                       derr = 0. _d 0
                0363
                0364                       fell =
                0365      &              +( 5. _d 0 / 1. _d 0) * ff00
                0366      &              -( 4. _d 0 / 1. _d 0) * ferr
                0367                       dell =
                0368      &              -(10. _d 0 / 1. _d 0) * ff00
                0369      &              +(10. _d 0 / 1. _d 0) * ferr
                0370
                0371                   end if
                0372
                0373               end if
                0374
                0375           end if
                0376
                0377 C     ============================= re-assemble coefficients on demand
                0378           if (mono .eq. +2) then
                0379
                0380           fhat(1) =
                0381      & + (30. _d 0 / 16. _d 0) * ff00
                0382      & - ( 7. _d 0 / 16. _d 0) *(ferr+fell)
                0383      & + ( 1. _d 0 / 16. _d 0) *(derr-dell)
                0384           fhat(2) =
                0385      & + ( 3. _d 0 /  4. _d 0) *(ferr-fell)
                0386      & - ( 1. _d 0 /  4. _d 0) *(derr+dell)
                0387           fhat(3) =
                0388      & - (30. _d 0 /  8. _d 0) * ff00
                0389      & + (15. _d 0 /  8. _d 0) *(ferr+fell)
                0390      & - ( 3. _d 0 /  8. _d 0) *(derr-dell)
                0391           fhat(4) =
                0392      & - ( 1. _d 0 /  4. _d 0) *(ferr-fell-derr-dell)
                0393           fhat(5) =
                0394      & + (30. _d 0 / 16. _d 0) * ff00
                0395      & - (15. _d 0 / 16. _d 0) *(ferr+fell)
                0396      & + ( 5. _d 0 / 16. _d 0) *(derr-dell)
                0397
                0398           end if
                0399
                0400           return
                0401
                0402 c     end subroutine GAD_PQM_FUN_MONO
                0403       end