C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_ppm_flx_y.F,v 1.1 2016/03/13 01:44:03 jmc Exp $ C $Name: $ # include "GAD_OPTIONS.h" SUBROUTINE GAD_PPM_FLX_Y(bi,bj,kk,ix, & calc_CFL,delT,vvel, & vfac,fhat,flux,myThid ) C |================================================================| C | PPM_FLX_Y: evaluate PPM 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 ix :: x-index. C calc_CFL :: TRUE to calc. CFL from vel. C delT :: time-step. C vvel :: vel.-comp in y-direction. C vfac :: vel.-flux in y-direction. C fhat :: row of poly. coeff. C flux :: adv.-flux in y-direction. C myThid :: thread number. C ================================================================ integer bi,bj,kk,ix logical calc_CFL _RL delT _RL vvel(1-OLx:sNx+OLx, & 1-OLy:sNy+OLy) _RL vfac(1-OLx:sNx+OLx, & 1-OLy:sNy+OLy) _RL fhat(1:3, & 1-OLy:sNy+OLy) _RL flux(1-OLx:sNx+OLx, & 1-OLy:sNy+OLy) integer myThid C ================================================================ C iy :: y-indexing. C vCFL :: CFL number. C intF :: upwind tracer edge-value. C ss11,ss22 :: int. endpoints. C ivec :: int. basis vec. C ================================================================ integer iy _RL vCFL,intF _RL ss11,ss22 _RL ivec(1:3) 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 iy = 1-OLy+3, sNy+OLy-2 if (vvel(ix,iy) .eq. 0. _d 0) then flux(ix,iy) = 0. _d 0 else if (vvel(ix,iy) .gt. 0. _d 0) then C ==================== integrate PPM profile over upwind cell IY-1 if ( calc_CFL ) then vCFL = vvel(ix,iy) * delT & * recip_dyF(ix,iy-1,bi,bj) & * recip_deepFacC(kk) else vCFL = vvel(ix,iy) end
if ss11 = +1. _d 0 - 2. _d 0 * vCFL 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) intF = ivec(1) * fhat(1,iy-1) & + ivec(2) * fhat(2,iy-1) & + ivec(3) * fhat(3,iy-1) intF = intF / (ss22 - ss11) else C ==================== integrate PPM profile over upwind cell IY+0 if ( calc_CFL ) then vCFL = vvel(ix,iy) * delT & * recip_dyF(ix,iy-0,bi,bj) & * recip_deepFacC(kk) else vCFL = vvel(ix,iy) end
if ss11 = -1. _d 0 - 2. _d 0 * vCFL 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) intF = ivec(1) * fhat(1,iy-0) & + ivec(2) * fhat(2,iy-0) & + ivec(3) * fhat(3,iy-0) intF = intF / (ss22 - ss11) end
if C ==================== calc. flux = upwind tracer * face-transport flux(ix,iy) = + vfac(ix,iy) * intF end
if end
do return c end subroutine GAD_PPM_FLX_Y end