C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_implicit_r.F,v 1.27 2016/10/26 21:36:03 jmc Exp $
C $Name: $
#include "GAD_OPTIONS.h"
CBOP
C !ROUTINE: GAD_IMPLICIT_R
C !INTERFACE:
SUBROUTINE GAD_IMPLICIT_R(
I implicitAdvection, advectionScheme, tracerIdentity,
I deltaTLev,
I kappaRX, recip_hFac, wFld, tracer,
U gTracer,
I bi, bj, myTime, myIter, myThid )
C !DESCRIPTION:
C Solve implicitly vertical advection and diffusion for one tracer.
C !USES:
IMPLICIT NONE
C == Global data ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "GAD.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine Arguments ==
C implicitAdvection :: if True, treat vertical advection implicitly
C advectionScheme :: advection scheme to use
C tracerIdentity :: Identifier for the tracer
C kappaRX :: 3-D array for vertical diffusion coefficient
C recip_hFac :: inverse of cell open-depth factor
C wFld :: Advection velocity field, vertical component
C tracer :: tracer field at current time step
C gTracer :: future tracer field
C bi,bj :: tile indices
C myTime :: current time
C myIter :: current iteration number
C myThid :: thread number
LOGICAL implicitAdvection
INTEGER advectionScheme
INTEGER tracerIdentity
_RL deltaTLev(Nr)
_RL kappaRX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RS recip_hFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL tracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL gTracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
INTEGER bi, bj
_RL myTime
INTEGER myIter, myThid
#ifdef ALLOW_DIAGNOSTICS
C !FUNCTIONS:
CHARACTER*4 GAD_DIAG_SUFX
EXTERNAL
LOGICAL DIAGNOSTICS_IS_ON
EXTERNAL
#endif
C !LOCAL VARIABLES:
C == Local variables ==
C iMin,iMax,jMin,jMax :: computational domain
C i,j,k :: loop indices
C a5d :: 2nd lower diagonal of the pentadiagonal matrix
C b5d :: 1rst lower diagonal of the pentadiagonal matrix
C c5d :: main diagonal of the pentadiagonal matrix
C d5d :: 1rst upper diagonal of the pentadiagonal matrix
C e5d :: 2nd upper diagonal of the pentadiagonal matrix
C rTrans :: vertical volume transport at interface k
C localTr :: local copy of tracer (for Non-Lin Adv.Scheme)
C diagonalNumber :: number of non-zero diagonals in the matrix
C errCode :: > 0 if singular matrix
C msgBuf :: Informational/error message buffer
INTEGER iMin,iMax,jMin,jMax
PARAMETER( iMin = 1, iMax = sNx )
PARAMETER( jMin = 1, jMax = sNy )
INTEGER i,j,k
INTEGER diagonalNumber, errCode
_RL a5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL b5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL c5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL d5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL e5d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL localTr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
#ifdef ALLOW_DIAGNOSTICS
CHARACTER*8 diagName
CHARACTER*4 diagSufx
LOGICAL diagDif, diagAdv
INTEGER km1, km2, kp1, kp2
_RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL div(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL flx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# ifdef SOLVE_DIAGONAL_LOWMEMORY
CHARACTER*(MAX_LEN_MBUF) msgBuf
# endif /* SOLVE_DIAGONAL_LOWMEMORY */
#endif /* ALLOW_DIAGNOSTICS */
CEOP
C-- no need to solve anything with only 1 level:
IF (Nr.GT.1) THEN
C-- Initialise
DO k=1,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
a5d(i,j,k) = 0. _d 0
b5d(i,j,k) = 0. _d 0
c5d(i,j,k) = 1. _d 0
d5d(i,j,k) = 0. _d 0
e5d(i,j,k) = 0. _d 0
ENDDO
ENDDO
ENDDO
diagonalNumber = 1
IF (implicitDiffusion) THEN
C-- set the tri-diagonal matrix to solve the implicit diffusion problem
diagonalNumber = 3
C- 1rst lower diagonal :
DO k=2,Nr
DO j=jMin,jMax
DO i=iMin,iMax
b5d(i,j,k) = -deltaTLev(k)*maskC(i,j,k-1,bi,bj)
& *recip_hFac(i,j,k)*recip_drF(k)
& *recip_deepFac2C(k)*recip_rhoFacC(k)
& *kappaRX(i,j, k )*recip_drC( k )
& *deepFac2F( k )*rhoFacF( k )
ENDDO
ENDDO
ENDDO
C- 1rst upper diagonal :
DO k=1,Nr-1
DO j=jMin,jMax
DO i=iMin,iMax
d5d(i,j,k) = -deltaTLev(k)*maskC(i,j,k+1,bi,bj)
& *recip_hFac(i,j,k)*recip_drF(k)
& *recip_deepFac2C(k)*recip_rhoFacC(k)
& *kappaRX(i,j,k+1)*recip_drC(k+1)
& *deepFac2F(k+1)*rhoFacF(k+1)
ENDDO
ENDDO
ENDDO
C- Main diagonal :
DO k=1,Nr
DO j=jMin,jMax
DO i=iMin,iMax
c5d(i,j,k) = 1. _d 0 - ( b5d(i,j,k) + d5d(i,j,k) )
C to recover older (prior to 2016-10-05) results:
c c5d(i,j,k) = 1. _d 0 - b5d(i,j,k) - d5d(i,j,k)
ENDDO
ENDDO
ENDDO
C-- end if implicitDiffusion
ENDIF
IF (implicitAdvection) THEN
C-- Non-Linear Advection scheme: keep a local copy of tracer field
IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT .OR.
& advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
IF ( multiDimAdvection ) THEN
DO k=1,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
localTr(i,j,k) = gTracer(i,j,k)
ENDDO
ENDDO
ENDDO
ELSE
DO k=1,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
localTr(i,j,k) = tracer(i,j,k)
ENDDO
ENDDO
ENDDO
ENDIF
ENDIF
DO k=Nr,1,-1
C-- Compute transport
IF (k.EQ.1) THEN
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
rTrans(i,j) = 0. _d 0
ENDDO
ENDDO
ELSE
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
rTrans(i,j) = wFld(i,j,k)*rA(i,j,bi,bj)
& *deepFac2F(k)*rhoFacF(k)
& *maskC(i,j,k-1,bi,bj)
ENDDO
ENDDO
ENDIF
#ifdef ALLOW_AIM
C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
IF ( k.GE.2 .AND.
& (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.k.LT.Nr)
& ) THEN
#else
IF ( k.GE.2 ) THEN
#endif
IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
diagonalNumber = 3
CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
I deltaTLev, rTrans, recip_hFac,
U b5d, c5d, d5d,
I myThid )
ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
& .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
diagonalNumber = 3
CALL GAD_DST2U1_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
I advectionScheme, deltaTLev,
I rTrans, recip_hFac,
U b5d, c5d, d5d,
I myThid )
ELSEIF ( advectionScheme.EQ.ENUM_FLUX_LIMIT ) THEN
diagonalNumber = 3
CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
I deltaTLev, rTrans, recip_hFac, localTr,
U b5d, c5d, d5d,
I myThid )
ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_3RD
& .OR. advectionScheme.EQ.ENUM_CENTERED_4TH
& .OR. advectionScheme.EQ.ENUM_DST3 ) THEN
diagonalNumber = 5
CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
I advectionScheme, deltaTLev,
I rTrans, recip_hFac,
U a5d, b5d, c5d, d5d, e5d,
I myThid )
ELSEIF ( advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
diagonalNumber = 5
CALL GAD_DST3FL_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
I deltaTLev, rTrans, recip_hFac, localTr,
U a5d, b5d, c5d, d5d, e5d,
I myThid )
ELSE
STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'
ENDIF
ENDIF
C-- end k loop
ENDDO
C-- end if implicitAdvection
ENDIF
IF ( diagonalNumber .EQ. 3 ) THEN
C-- Solve tri-diagonal system :
errCode = -1
CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
I b5d, c5d, d5d,
U gTracer,
O errCode,
I bi, bj, myThid )
IF (errCode.GE.1) THEN
STOP 'GAD_IMPLICIT_R: error when solving 3-Diag problem'
ENDIF
ELSEIF ( diagonalNumber .EQ. 5 ) THEN
C-- Solve penta-diagonal system :
errCode = -1
CALL SOLVE_PENTADIAGONAL( iMin,iMax, jMin,jMax,
I a5d, b5d, c5d, d5d, e5d,
U gTracer,
O errCode,
I bi, bj, myThid )
IF (errCode.GE.1) THEN
STOP 'GAD_IMPLICIT_R: error when solving 5-Diag problem'
ENDIF
ELSEIF ( diagonalNumber .NE. 1 ) THEN
STOP 'GAD_IMPLICIT_R: no solver available'
ENDIF
#ifdef ALLOW_DIAGNOSTICS
C-- Set diagnostic suffix for the current tracer
IF ( useDiagnostics ) THEN
diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
diagName = 'DFrI'//diagSufx
diagDif = implicitDiffusion
IF ( diagDif ) diagDif = DIAGNOSTICS_IS_ON(diagName,myThid)
diagName = 'ADVr'//diagSufx
diagAdv = implicitAdvection
IF ( diagAdv ) diagAdv = DIAGNOSTICS_IS_ON(diagName,myThid)
IF ( diagDif .OR. diagAdv ) THEN
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
flx(i,j) = 0. _d 0
ENDDO
ENDDO
C-- start diagnostics k loop
DO k= Nr,1,-1
IF ( implicitDiffusion .AND. k.GE.2 ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
df(i,j) =
cc#ifdef ALLOW_AUTODIFF_OPENAD
cc & -rA(i,j,bi,bj)%v
cc#else
& -rA(i,j,bi,bj)*deepFac2F(k)*rhoFacF(k)
cc#endif
& * kappaRX(i,j,k)*recip_drC(k)*rkSign
& * (gTracer(i,j,k) - gTracer(i,j,k-1))
& * maskC(i,j,k,bi,bj)
& * maskC(i,j,k-1,bi,bj)
ENDDO
ENDDO
ELSE
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
df(i,j) = 0. _d 0
ENDDO
ENDDO
ENDIF
C- Note: Needs to explicitly increment counter (call DIAGNOSTICS_COUNT)
C since skipping k=1 DIAGNOSTICS_FILL call.
IF ( diagDif .AND. k.GE.2 ) THEN
diagName = 'DFrI'//diagSufx
CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
#ifdef ALLOW_LAYERS
IF ( useLayers ) THEN
CALL LAYERS_FILL( df, tracerIdentity, 'DFR',
& k, 1, 2,bi,bj, myThid )
ENDIF
#endif /* ALLOW_LAYERS */
ENDIF
IF ( diagAdv ) THEN
#ifdef SOLVE_DIAGONAL_LOWMEMORY
diagName = 'ADVr'//diagSufx
WRITE(msgBuf,'(4A)') 'GAD_IMPLICIT_R: ',
& 'unable to compute Diagnostic "', diagName, '" with'
CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
& SQUEEZE_RIGHT, myThid )
WRITE(msgBuf,'(4A)') 'GAD_IMPLICIT_R: ',
& '#define SOLVE_DIAGONAL_LOWMEMORY (in CPP_OPTIONS.h)'
CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
& SQUEEZE_RIGHT, myThid )
STOP 'ABNORMAL END: S/R GAD_IMPLICIT_R'
#endif /* SOLVE_DIAGONAL_LOWMEMORY */
km1=MAX(1,k-1)
km2=MAX(1,k-2)
kp1=MIN(Nr,k+1)
kp2=MIN(Nr,k+2)
C-- Flux_divergence*deltaT = Tr^n - Tr^n+1 = [A-I](Tr^n+1)
C = deltaT*rkSign*[ Flx_k+1 - Flx_k ]/dz
DO j=jMin,jMax
DO i=iMin,iMax
div(i,j) = gTracer(i,j,k)*( c5d(i,j,k) - 1. _d 0 )
& + gTracer(i,j,km1)*b5d(i,j,k)
& + gTracer(i,j,kp1)*d5d(i,j,k)
ENDDO
ENDDO
IF ( diagonalNumber .EQ. 5 ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
div(i,j) = div(i,j)
& + gTracer(i,j,km2)*a5d(i,j,k)
& + gTracer(i,j,kp2)*e5d(i,j,k)
ENDDO
ENDDO
ENDIF
#ifdef NONLIN_FRSURF
IF ( nonlinFreeSurf.GT.0 ) THEN
C-- use future hFac to stay consistent with solver matrix
IF ( select_rStar.GT.0 ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
div(i,j) = div(i,j)*h0FacC(i,j,k,bi,bj)*drF(k)
& *rStarFacC(i,j,bi,bj)
ENDDO
ENDDO
ELSEIF ( selectSigmaCoord.NE.0 ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
div(i,j) = div(i,j)*(
& _hFacC(i,j,k,bi,bj)*drF(k)
& + dBHybSigF(k)*dEtaHdt(i,j,bi,bj)*deltaTFreeSurf
& )
ENDDO
ENDDO
ELSE
DO j=jMin,jMax
DO i=iMin,iMax
IF ( k.EQ.kSurfC(i,j,bi,bj) ) THEN
div(i,j) = div(i,j)*hFac_surfC(i,j,bi,bj)*drF(k)
ELSE
div(i,j) = div(i,j)*_hFacC(i,j,k,bi,bj)*drF(k)
ENDIF
ENDDO
ENDDO
ENDIF
ELSE
#else /* NONLIN_FRSURF */
IF ( .TRUE. ) THEN
#endif /* NONLIN_FRSURF */
C-- use current hFac (consistent with solver matrix)
DO j=jMin,jMax
DO i=iMin,iMax
div(i,j) = div(i,j)*_hFacC(i,j,k,bi,bj)*drF(k)
ENDDO
ENDDO
ENDIF
DO j=jMin,jMax
DO i=iMin,iMax
flx(i,j) = flx(i,j)
cc#ifdef ALLOW_AUTODIFF_OPENAD
cc & - rkSign*div(i,j)*rA(i,j,bi,bj)%v/deltaTLev(k)
cc#else
& - rkSign*div(i,j)*rA(i,j,bi,bj)
& *deepFac2C(k)*rhoFacC(k)/deltaTLev(k)
cc#endif
af(i,j) = flx(i,j) - df(i,j)
ENDDO
ENDDO
diagName = 'ADVr'//diagSufx
CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
#ifdef ALLOW_LAYERS
IF ( useLayers ) THEN
CALL LAYERS_FILL(af,tracerIdentity,'AFR',
& k,1,2,bi,bj,myThid)
ENDIF
#endif /* ALLOW_LAYERS */
ENDIF
C-- end diagnostics k loop
ENDDO
ENDIF
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
C-- end if Nr > 1
ENDIF
RETURN
END