C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_pqm_fun.F,v 1.1 2016/03/13 01:44:03 jmc Exp $
C $Name: $
# include "GAD_OPTIONS.h"
C-- File gad_pqm_fun.F: Routines to form PQM grid-cell polynomial.
C-- Contents
C-- o QUADROOT
C-- o GAD_PPM_FUN_NULL
C-- o GAD_PPM_FUN_MONO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
LOGICAL FUNCTION QUADROOT(aa,bb,cc,xx)
C |================================================================|
C | QUADROOT: find roots of quadratic ax**2 + bx + c = 0. |
C |================================================================|
implicit none
C ====================================================== arguments
_RL aa,bb,cc
_RL xx(1:2)
C ====================================================== variables
_RL sq,a0,b0
a0 = abs(aa)
b0 = abs(bb)
sq = bb * bb - 4. _d 0 * aa * cc
if (a0 .gt. 0. _d 0) then
if (sq .ge. 0. _d 0) then
QUADROOT = .TRUE.
sq = sqrt(sq)
xx(1) = - bb + sq
xx(2) = - bb - sq
aa = 0.5 _d 0 / aa
xx(1) = xx(1) * aa
xx(2) = xx(2) * aa
else
QUADROOT = .FALSE.
end
if
else
if (b0 .gt. 0. _d 0) then
QUADROOT = .TRUE.
xx(1) = - cc / bb
xx(2) = - cc / bb
else
QUADROOT = .FALSE.
end
if
end
if
return
c end function QUADROOT
end
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE GAD_PQM_FUN_NULL(
I ff00,
I fell,ferr,
I dell,derr,
O fhat,mono)
C |================================================================|
C | PQM_FUN_NULL: form PQM grid-cell polynomial. |
C | Piecewise Quartic Method (PQM) - unlimited variant. |
C |================================================================|
implicit none
C ====================================================== arguments
_RL ff00
_RL fell,ferr
_RL dell,derr
_RL fhat(+1:+5)
integer mono
mono = +0
C ============================================== unlimited profile
fhat(1) =
& + (30. _d 0 / 16. _d 0) * ff00
& - ( 7. _d 0 / 16. _d 0) *(ferr+fell)
& + ( 1. _d 0 / 16. _d 0) *(derr-dell)
fhat(2) =
& + ( 3. _d 0 / 4. _d 0) *(ferr-fell)
& - ( 1. _d 0 / 4. _d 0) *(derr+dell)
fhat(3) =
& - (30. _d 0 / 8. _d 0) * ff00
& + (15. _d 0 / 8. _d 0) *(ferr+fell)
& - ( 3. _d 0 / 8. _d 0) *(derr-dell)
fhat(4) =
& - ( 1. _d 0 / 4. _d 0) *(ferr-fell-derr-dell)
fhat(5) =
& + (30. _d 0 / 16. _d 0) * ff00
& - (15. _d 0 / 16. _d 0) *(ferr+fell)
& + ( 5. _d 0 / 16. _d 0) *(derr-dell)
return
c end subroutine GAD_PQM_FUN_NULL
end
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
SUBROUTINE GAD_PQM_FUN_MONO(
I ff00,
I ffll,ffrr,
I fell,ferr,
I dell,derr,
I dfds,
O fhat,mono)
C |================================================================|
C | PQM_FUN_MONO: form PQM grid-cell polynomial. |
C | Piecewise Quartic Method (PQM) - monotonic variant. |
C |================================================================|
implicit none
C ====================================================== arguments
_RL ff00,ffll,ffrr
_RL fell,ferr
_RL dell,derr
_RL dfds(-1:+1)
_RL fhat(+1:+5)
integer mono
C ====================================================== functions
logical QUADROOT
C ====================================================== variables
integer bind
_RL aval,bval,cval
_RL iflx(+1:+2)
_RL dflx(+1:+2)
mono = +0
C ============================================== "flatten" extrema
if((ffrr-ff00) *
& (ff00-ffll) .le. 0. _d 0) then
mono = +1
fhat(1) = ff00
fhat(2) = 0. _d 0
fhat(3) = 0. _d 0
fhat(4) = 0. _d 0
fhat(5) = 0. _d 0
return
end
if
C ============================================== limit edge values
if((ffll - fell) *
& (fell - ff00) .le. 0. _d 0) then
mono = +1
fell = ff00 - dfds(0)
end
if
if((ffrr - ferr) *
& (ferr - ff00) .le. 0. _d 0) then
mono = +1
ferr = ff00 + dfds(0)
end
if
C ============================================== limit edge slopes
if((dell * dfds(-1)) .lt. 0. _d 0) then
mono = +1
dell = dfds(-1)
end
if
if((derr * dfds(+1)) .lt. 0. _d 0) then
mono = +1
derr = dfds(+1)
end
if
C ============================================== limit cell values
fhat(1) =
& + (30. _d 0 / 16. _d 0) * ff00
& - ( 7. _d 0 / 16. _d 0) *(ferr+fell)
& + ( 1. _d 0 / 16. _d 0) *(derr-dell)
fhat(2) =
& + ( 3. _d 0 / 4. _d 0) *(ferr-fell)
& - ( 1. _d 0 / 4. _d 0) *(derr+dell)
fhat(3) =
& - (30. _d 0 / 8. _d 0) * ff00
& + (15. _d 0 / 8. _d 0) *(ferr+fell)
& - ( 3. _d 0 / 8. _d 0) *(derr-dell)
fhat(4) =
& - ( 1. _d 0 / 4. _d 0) *(ferr-fell-derr-dell)
fhat(5) =
& + (30. _d 0 / 16. _d 0) * ff00
& - (15. _d 0 / 16. _d 0) *(ferr+fell)
& + ( 5. _d 0 / 16. _d 0) *(derr-dell)
C ============================= calc. inflexion via 2nd-derivative
aval = 12. _d 0 * fhat(5)
bval = 6. _d 0 * fhat(4)
cval = 2. _d 0 * fhat(3)
if ( QUADROOT(aval,bval,cval,iflx) ) then
bind = +0
if ( ( iflx(1) .gt. -1. _d 0 )
& .and. ( iflx(1) .lt. +1. _d 0 ) ) then
C ============================= check for non-monotonic inflection
dflx(1) = fhat(2)
& + iflx(1) * fhat(3) * 2. _d 0
& +(iflx(1) ** 2) * fhat(4) * 3. _d 0
& +(iflx(1) ** 3) * fhat(5) * 4. _d 0
if (dflx(1)*dfds(+0) .lt. 0. _d 0) then
if (abs(dell)
& .lt. abs(derr) ) then
bind = -1
else
bind = +1
end
if
end
if
end
if
if ( ( iflx(2) .gt. -1. _d 0 )
& .and. ( iflx(2) .lt. +1. _d 0 ) ) then
C ============================= check for non-monotonic inflection
dflx(2) = fhat(2)
& + iflx(2) * fhat(3) * 2. _d 0
& +(iflx(2) ** 2) * fhat(4) * 3. _d 0
& +(iflx(2) ** 3) * fhat(5) * 4. _d 0
if (dflx(2)*dfds(+0) .lt. 0. _d 0) then
if (abs(dell)
& .lt. abs(derr) ) then
bind = -1
else
bind = +1
end
if
end
if
end
if
C ============================= pop non-monotone inflexion to edge
if (bind .eq. -1) then
C ============================= pop inflection points onto -1 edge
mono = +2
derr =
& -( 5. _d 0 / 1. _d 0) * ff00
& +( 3. _d 0 / 1. _d 0) * ferr
& +( 2. _d 0 / 1. _d 0) * fell
dell =
& +( 5. _d 0 / 3. _d 0) * ff00
& -( 1. _d 0 / 3. _d 0) * ferr
& -( 4. _d 0 / 3. _d 0) * fell
if (dell*dfds(-1) .lt. 0. _d 0) then
dell = 0. _d 0
ferr =
& +( 5. _d 0 / 1. _d 0) * ff00
& -( 4. _d 0 / 1. _d 0) * fell
derr =
& +(10. _d 0 / 1. _d 0) * ff00
& -(10. _d 0 / 1. _d 0) * fell
end
if
if (derr*dfds(+1) .lt. 0. _d 0) then
derr = 0. _d 0
fell =
& +( 5. _d 0 / 2. _d 0) * ff00
& -( 3. _d 0 / 2. _d 0) * ferr
dell =
& -( 5. _d 0 / 3. _d 0) * ff00
& +( 5. _d 0 / 3. _d 0) * ferr
end
if
end
if
if (bind .eq. +1) then
C ============================= pop inflection points onto +1 edge
mono = +2
derr =
& -( 5. _d 0 / 3. _d 0) * ff00
& +( 4. _d 0 / 3. _d 0) * ferr
& +( 1. _d 0 / 3. _d 0) * fell
dell =
& +( 5. _d 0 / 1. _d 0) * ff00
& -( 2. _d 0 / 1. _d 0) * ferr
& -( 3. _d 0 / 1. _d 0) * fell
if (dell*dfds(-1) .lt. 0. _d 0) then
dell = 0. _d 0
ferr =
& +( 5. _d 0 / 2. _d 0) * ff00
& -( 3. _d 0 / 2. _d 0) * fell
derr =
& +( 5. _d 0 / 3. _d 0) * ff00
& -( 5. _d 0 / 3. _d 0) * fell
end
if
if (derr*dfds(+1) .lt. 0. _d 0) then
derr = 0. _d 0
fell =
& +( 5. _d 0 / 1. _d 0) * ff00
& -( 4. _d 0 / 1. _d 0) * ferr
dell =
& -(10. _d 0 / 1. _d 0) * ff00
& +(10. _d 0 / 1. _d 0) * ferr
end
if
end
if
end
if
C ============================= re-assemble coefficients on demand
if (mono .eq. +2) then
fhat(1) =
& + (30. _d 0 / 16. _d 0) * ff00
& - ( 7. _d 0 / 16. _d 0) *(ferr+fell)
& + ( 1. _d 0 / 16. _d 0) *(derr-dell)
fhat(2) =
& + ( 3. _d 0 / 4. _d 0) *(ferr-fell)
& - ( 1. _d 0 / 4. _d 0) *(derr+dell)
fhat(3) =
& - (30. _d 0 / 8. _d 0) * ff00
& + (15. _d 0 / 8. _d 0) *(ferr+fell)
& - ( 3. _d 0 / 8. _d 0) *(derr-dell)
fhat(4) =
& - ( 1. _d 0 / 4. _d 0) *(ferr-fell-derr-dell)
fhat(5) =
& + (30. _d 0 / 16. _d 0) * ff00
& - (15. _d 0 / 16. _d 0) *(ferr+fell)
& + ( 5. _d 0 / 16. _d 0) *(derr-dell)
end
if
return
c end subroutine GAD_PQM_FUN_MONO
end