C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_pqm_flx_x.F,v 1.1 2016/03/13 01:44:03 jmc Exp $
C $Name: $
# include "GAD_OPTIONS.h"
SUBROUTINE GAD_PQM_FLX_X(bi,bj,kk,iy,
& calc_CFL,delT,uvel,
& ufac,fhat,flux,myThid)
C |================================================================|
C | PQM_FLX_X: evaluate PQM flux on grid-cell edges. |
C |================================================================|
implicit none
C =============================================== global variables
# include "SIZE.h"
# include "GRID.h"
# include "GAD.h"
C ================================================================
C bi,bj :: tile indexing.
C kk :: r-index.
C iy :: y-index.
C calc_CFL :: TRUE to calc. CFL from vel.
C delT :: time-step.
C uvel :: vel.-comp in x-direction.
C ufac :: vel.-flux in x-direction.
C fhat :: row of poly. coeff.
C flux :: adv.-flux in x-direction.
C myThid :: thread number.
C ================================================================
integer bi,bj,kk,iy
logical calc_CFL
_RL delT
_RL uvel(1-OLx:sNx+OLx,
& 1-OLy:sNy+OLy)
_RL ufac(1-OLx:sNx+OLx,
& 1-OLy:sNy+OLy)
_RL fhat(1:5,
& 1-OLx:sNx+OLx)
_RL flux(1-OLx:sNx+OLx,
& 1-OLy:sNy+OLy)
integer myThid
C ================================================================
C ix :: x-indexing.
C uCFL :: CFL number.
C intF :: upwind tracer edge-value.
C ss11,ss22 :: int. endpoints.
C ivec :: int. basis vec.
C ================================================================
integer ix
_RL uCFL,intF
_RL ss11,ss22
_RL ivec(1:5)
C ================================================================
C (1): calc. "departure-points" for each grid-cell edge by int-
C egrating edge position backward in time over one single
C time-step. This is a "single-cell" implementation: requ-
C ires CFL <= 1.0.
C (2): calc. flux as the integral of the upwind grid-cell poly-
C nomial over the deformation interval found in (1).
C ================================================================
do ix = 1-OLx+4, sNx+OLx-3
if (uvel(ix,iy) .eq. 0. _d 0) then
flux(ix,iy) = 0. _d 0
else
if (uvel(ix,iy) .gt. 0. _d 0) then
C ==================== integrate PQM profile over upwind cell IX-1
if ( calc_CFL ) then
uCFL = uvel(ix,iy) * delT
& * recip_dxF(ix-1,iy,bi,bj)
& * recip_deepFacC(kk)
else
uCFL = uvel(ix,iy)
end
if
ss11 = +1. _d 0 - 2. _d 0 * uCFL
ss22 = +1. _d 0
C ==================== integrate profile over region swept by face
ivec(1) = ss22 - ss11
ivec(2) =(ss22 ** 2
& - ss11 ** 2)*(1. _d 0 / 2. _d 0)
ivec(3) =(ss22 ** 3
& - ss11 ** 3)*(1. _d 0 / 3. _d 0)
ivec(4) =(ss22 ** 4
& - ss11 ** 4)*(1. _d 0 / 4. _d 0)
ivec(5) =(ss22 ** 5
& - ss11 ** 5)*(1. _d 0 / 5. _d 0)
intF = ivec(1) * fhat(1,ix-1)
& + ivec(2) * fhat(2,ix-1)
& + ivec(3) * fhat(3,ix-1)
& + ivec(4) * fhat(4,ix-1)
& + ivec(5) * fhat(5,ix-1)
intF = intF / (ss22 - ss11)
else
C ==================== integrate PQM profile over upwind cell IX+0
if ( calc_CFL ) then
uCFL = uvel(ix,iy) * delT
& * recip_dxF(ix-0,iy,bi,bj)
& * recip_deepFacC(kk)
else
uCFL = uvel(ix,iy)
end
if
ss11 = -1. _d 0 - 2. _d 0 * uCFL
ss22 = -1. _d 0
C ==================== integrate profile over region swept by face
ivec(1) = ss22 - ss11
ivec(2) =(ss22 ** 2
& - ss11 ** 2)*(1. _d 0 / 2. _d 0)
ivec(3) =(ss22 ** 3
& - ss11 ** 3)*(1. _d 0 / 3. _d 0)
ivec(4) =(ss22 ** 4
& - ss11 ** 4)*(1. _d 0 / 4. _d 0)
ivec(5) =(ss22 ** 5
& - ss11 ** 5)*(1. _d 0 / 5. _d 0)
intF = ivec(1) * fhat(1,ix-0)
& + ivec(2) * fhat(2,ix-0)
& + ivec(3) * fhat(3,ix-0)
& + ivec(4) * fhat(4,ix-0)
& + ivec(5) * fhat(5,ix-0)
intF = intF / (ss22 - ss11)
end
if
C ==================== calc. flux = upwind tracer * face-transport
flux(ix,iy) = + ufac(ix,iy) * intF
end
if
end
do
return
c end subroutine GAD_PQM_FLX_X
end