C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_som_advect.F,v 1.14 2015/06/03 13:39:22 rpa Exp $ C $Name: $ #include "GAD_OPTIONS.h" C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: GAD_SOM_ADVECT C !INTERFACE: ========================================================== SUBROUTINE GAD_SOM_ADVECT( I implicitAdvection, advectionScheme, vertAdvecScheme, I tracerIdentity, deltaTLev, I uFld, vFld, wFld, tracer, U smTr, O gTracer, I bi,bj, myTime,myIter,myThid) C !DESCRIPTION: C Calculates the tendency of a tracer due to advection. C It uses the 2nd-Order moment advection scheme with multi-dimensional method C see Prather, 1986, JGR, v.91, D-6, pp.6671-6681. C C The tendency (output) is over-written by this routine. C !USES: =============================================================== IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "GAD.h" #ifdef ALLOW_EXCH2 #include "W2_EXCH2_SIZE.h" #include "W2_EXCH2_TOPOLOGY.h" #endif /* ALLOW_EXCH2 */ C !INPUT PARAMETERS: =================================================== C implicitAdvection :: implicit vertical advection (later on) C advectionScheme :: advection scheme to use (Horizontal plane) C vertAdvecScheme :: advection scheme to use (vertical direction) C tracerIdentity :: tracer identifier (required only for OBCS) C uFld :: Advection velocity field, zonal component C vFld :: Advection velocity field, meridional component C wFld :: Advection velocity field, vertical component C tracer :: tracer field C bi,bj :: tile indices C myTime :: current time C myIter :: iteration number C myThid :: thread number LOGICAL implicitAdvection INTEGER advectionScheme, vertAdvecScheme INTEGER tracerIdentity _RL deltaTLev(Nr) _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL vFld (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,nSx,nSy) INTEGER bi,bj _RL myTime INTEGER myIter INTEGER myThid C !OUTPUT PARAMETERS: ================================================== C smTr :: tracer 1rst & 2nd Order moments C gTracer :: tendency array _RL smTr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy,nSOM) _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) C !LOCAL VARIABLES: ==================================================== C maskUp :: 2-D array mask for W points C i,j,k :: loop indices C kUp :: index into 2 1/2D array, toggles between 1 and 2 C kDown :: index into 2 1/2D array, toggles between 2 and 1 C xA,yA :: areas of X and Y face of tracer cells C uTrans,vTrans :: 2-D arrays of volume transports at U,V points C rTrans :: 2-D arrays of volume transports at W points C afx :: 2-D array for horizontal advective flux, x direction C afy :: 2-D array for horizontal advective flux, y direction C afr :: 2-D array for vertical advective flux C fVerT :: 2 1/2D arrays for vertical advective flux C calc_fluxes_X :: logical to indicate to calculate fluxes in X dir C calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir C interiorOnly :: only update the interior of myTile, but not the edges C overlapOnly :: only update the edges of myTile, but not the interior C npass :: number of passes in multi-dimensional method C ipass :: number of the current pass being made C myTile :: variables used to determine which cube face C nCFace :: owns a tile for cube grid runs using C :: multi-dim advection. C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube C msgBuf :: Informational/error message buffer _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER i,j,k,km1,kUp,kDown _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL afx (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL afr (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL smVol (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL smTr0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL alp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL aln (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fp_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fn_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fp_o (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fn_o (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fp_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fn_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fp_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fn_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fp_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fn_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fp_xx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fn_xx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fp_yy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fn_yy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fp_zz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fn_zz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fp_xy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fn_xy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fp_xz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fn_xz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fp_yz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL fn_yz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2) _RL smCorners(OLx,OLy,4,-1:nSOM) c _RL localTr LOGICAL calc_fluxes_X, calc_fluxes_Y LOGICAL interiorOnly, overlapOnly INTEGER limiter INTEGER npass, ipass INTEGER nCFace, n LOGICAL N_edge, S_edge, E_edge, W_edge LOGICAL noFlowAcrossSurf CHARACTER*(MAX_LEN_MBUF) msgBuf #ifdef ALLOW_EXCH2 INTEGER myTile #endif #ifdef ALLOW_DIAGNOSTICS CHARACTER*8 diagName CHARACTER*4 diagSufx LOGICAL doDiagAdvX, doDiagAdvY, doDiagAdvR C- Functions: CHARACTER*4 GAD_DIAG_SUFX EXTERNAL GAD_DIAG_SUFX LOGICAL DIAGNOSTICS_IS_ON EXTERNAL DIAGNOSTICS_IS_ON #endif CEOP #ifdef ALLOW_DIAGNOSTICS C-- Set diagnostics flags and suffix for the current tracer doDiagAdvX = .FALSE. doDiagAdvY = .FALSE. doDiagAdvR = .FALSE. IF ( useDiagnostics ) THEN diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid ) diagName = 'ADVx'//diagSufx doDiagAdvX = DIAGNOSTICS_IS_ON( diagName, myThid ) diagName = 'ADVy'//diagSufx doDiagAdvY = DIAGNOSTICS_IS_ON( diagName, myThid ) diagName = 'ADVr'//diagSufx doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid ) ENDIF #endif C-- Set up work arrays with valid (i.e. not NaN) values C These inital values do not alter the numerical results. They C just ensure that all memory references are to valid floating C point numbers. This prevents spurious hardware signals due to C uninitialised but inert locations. DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx afx(i,j) = 0. afy(i,j) = 0. C- xA,yA,uTrans,vTrans are set over the full domain C => no need for extra initialisation c xA(i,j) = 0. _d 0 c yA(i,j) = 0. _d 0 c uTrans(i,j) = 0. _d 0 c vTrans(i,j) = 0. _d 0 C- rTrans is set over the full domain: no need for extra initialisation c rTrans(i,j) = 0. _d 0 ENDDO ENDDO DO n=-1,nSOM DO k=1,4 DO j=1,OLy DO i=1,OLx smCorners(i,j,k,n) = 0. ENDDO ENDDO ENDDO ENDDO IF ( implicitAdvection ) THEN WRITE(msgBuf,'(2A)') 'S/R GAD_SOM_ADVECT: ', & 'not coded for implicit-vertical Advection' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R GAD_SOM_ADVECT' ENDIF IF ( vertAdvecScheme .NE. advectionScheme ) THEN WRITE(msgBuf,'(2A)') 'S/R GAD_SOM_ADVECT: ', & 'not coded for different vertAdvecScheme' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R GAD_SOM_ADVECT' ENDIF C-- Set tile-specific parameters for horizontal fluxes IF (useCubedSphereExchange) THEN npass = 3 #ifdef ALLOW_EXCH2 myTile = W2_myTileList(bi,bj) nCFace = exch2_myFace(myTile) N_edge = exch2_isNedge(myTile).EQ.1 S_edge = exch2_isSedge(myTile).EQ.1 E_edge = exch2_isEedge(myTile).EQ.1 W_edge = exch2_isWedge(myTile).EQ.1 #else nCFace = bi N_edge = .TRUE. S_edge = .TRUE. E_edge = .TRUE. W_edge = .TRUE. #endif ELSE npass = 2 nCFace = 0 N_edge = .FALSE. S_edge = .FALSE. E_edge = .FALSE. W_edge = .FALSE. ENDIF limiter = MOD(advectionScheme, 10) C-- Start of k loop for horizontal fluxes DO k=1,Nr C-- Get temporary terms used by tendency routines DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k) & *drF(k)*_hFacW(i,j,k,bi,bj) yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k) & *drF(k)*_hFacS(i,j,k,bi,bj) ENDDO ENDDO C-- Calculate "volume transports" through tracer cell faces. C anelastic: scaled by rhoFacC (~ mass transport) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx uTrans(i,j) = uFld(i,j,k)*xA(i,j)*rhoFacC(k) vTrans(i,j) = vFld(i,j,k)*yA(i,j)*rhoFacC(k) ENDDO ENDDO C-- grid-box volume and tracer content (zero order moment) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx smVol(i,j,k) = rA(i,j,bi,bj)*deepFac2C(k) & *drF(k)*hFacC(i,j,k,bi,bj) & *rhoFacC(k) smTr0(i,j,k) = tracer(i,j,k,bi,bj)*smVol(i,j,k) C- fill empty grid-box: smVol(i,j,k) = smVol(i,j,k) & + (1. _d 0 - maskC(i,j,k,bi,bj)) ENDDO ENDDO C-- Multiple passes for different directions on different tiles C-- For cube need one pass for each of red, green and blue axes. DO ipass=1,npass interiorOnly = .FALSE. overlapOnly = .FALSE. IF (useCubedSphereExchange) THEN C- CubedSphere : pass 3 times, with partial update of local tracer field IF (ipass.EQ.1) THEN overlapOnly = MOD(nCFace,3).EQ.0 interiorOnly = MOD(nCFace,3).NE.0 calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2 calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5 ELSEIF (ipass.EQ.2) THEN overlapOnly = MOD(nCFace,3).EQ.2 interiorOnly = MOD(nCFace,3).EQ.1 calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4 calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1 ELSE interiorOnly = .TRUE. calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6 calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3 ENDIF ELSE C- not CubedSphere calc_fluxes_X = MOD(ipass,2).EQ.1 calc_fluxes_Y = .NOT.calc_fluxes_X ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- X direction C- Do not compute fluxes if C a) needed in overlap only C and b) the overlap of myTile are not cube-face Edges IF ( calc_fluxes_X .AND. & (.NOT.overlapOnly .OR. N_edge .OR. S_edge) & ) THEN C- Internal exchange for calculations in X IF ( useCubedSphereExchange .AND. .NOT.interiorOnly ) THEN CALL GAD_SOM_PREP_CS_CORNER( U smVol, smTr0, smTr, smCorners, I .TRUE., overlapOnly, interiorOnly, I N_edge, S_edge, E_edge, W_edge, I ipass, k, Nr, bi, bj, myThid ) ENDIF C- Solve advection in X and update moments IF ( advectionScheme.EQ.ENUM_SOM_PRATHER & .OR. advectionScheme.EQ.ENUM_SOM_LIMITER ) THEN CALL GAD_SOM_ADV_X( I bi,bj,k, limiter, I overlapOnly, interiorOnly, I N_edge, S_edge, E_edge, W_edge, I deltaTLev(k), uTrans, I maskInC(1-OLx,1-OLy,bi,bj), U smVol(1-OLx,1-OLy,k), U smTr0(1-OLx,1-OLy,k), U smTr(1-OLx,1-OLy,k,bi,bj,1), U smTr(1-OLx,1-OLy,k,bi,bj,2), U smTr(1-OLx,1-OLy,k,bi,bj,3), U smTr(1-OLx,1-OLy,k,bi,bj,4), U smTr(1-OLx,1-OLy,k,bi,bj,5), U smTr(1-OLx,1-OLy,k,bi,bj,6), U smTr(1-OLx,1-OLy,k,bi,bj,7), U smTr(1-OLx,1-OLy,k,bi,bj,8), U smTr(1-OLx,1-OLy,k,bi,bj,9), O afx, myThid ) ELSE STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM' ENDIF C-- End of X direction ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- Y direction C- Do not compute fluxes if C a) needed in overlap only C and b) the overlap of myTile are not cube-face edges IF ( calc_fluxes_Y .AND. & (.NOT.overlapOnly .OR. E_edge .OR. W_edge) & ) THEN C- Internal exchange for calculations in Y IF ( useCubedSphereExchange .AND. .NOT.interiorOnly ) THEN CALL GAD_SOM_PREP_CS_CORNER( U smVol, smTr0, smTr, smCorners, I .FALSE., overlapOnly, interiorOnly, I N_edge, S_edge, E_edge, W_edge, I iPass, k, Nr, bi, bj, myThid ) ENDIF C- Solve advection in Y and update moments IF ( advectionScheme.EQ.ENUM_SOM_PRATHER & .OR. advectionScheme.EQ.ENUM_SOM_LIMITER ) THEN CALL GAD_SOM_ADV_Y( I bi,bj,k, limiter, I overlapOnly, interiorOnly, I N_edge, S_edge, E_edge, W_edge, I deltaTLev(k), vTrans, I maskInC(1-OLx,1-OLy,bi,bj), U smVol(1-OLx,1-OLy,k), U smTr0(1-OLx,1-OLy,k), U smTr(1-OLx,1-OLy,k,bi,bj,1), U smTr(1-OLx,1-OLy,k,bi,bj,2), U smTr(1-OLx,1-OLy,k,bi,bj,3), U smTr(1-OLx,1-OLy,k,bi,bj,4), U smTr(1-OLx,1-OLy,k,bi,bj,5), U smTr(1-OLx,1-OLy,k,bi,bj,6), U smTr(1-OLx,1-OLy,k,bi,bj,7), U smTr(1-OLx,1-OLy,k,bi,bj,8), U smTr(1-OLx,1-OLy,k,bi,bj,9), O afy, myThid ) ELSE STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM' ENDIF C-- End of Y direction ENDIF C-- End of ipass loop ENDDO IF ( implicitAdvection ) THEN C- explicit advection is done ; store tendency in gTracer: DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx C-- without rescaling of tendencies: c localTr = smTr0(i,j,k)/smVol(i,j,k) c gTracer(i,j,k) = ( localTr - tracer(i,j,k,bi,bj) ) c & / deltaTLev(k) C-- consistent with rescaling of tendencies (in FREESURF_RESCALE_G): gTracer(i,j,k) = & ( smTr0(i,j,k) - tracer(i,j,k,bi,bj)*smVol(i,j,k) ) & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj) & *recip_rhoFacC(k) & /deltaTLev(k) ENDDO ENDDO ENDIF #ifdef ALLOW_DIAGNOSTICS IF ( doDiagAdvX ) THEN diagName = 'ADVx'//diagSufx CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid ) ENDIF IF ( doDiagAdvY ) THEN diagName = 'ADVy'//diagSufx CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid ) ENDIF #ifdef ALLOW_LAYERS IF ( useLayers ) THEN CALL LAYERS_FILL(afx,tracerIdentity,'AFX',k,1,2,bi,bj,myThid) CALL LAYERS_FILL(afy,tracerIdentity,'AFY',k,1,2,bi,bj,myThid) ENDIF #endif /* ALLOW_LAYERS */ #endif #ifdef ALLOW_DEBUG IF ( debugLevel .GE. debLevC & .AND. tracerIdentity.EQ.GAD_TEMPERATURE & .AND. k.LE.3 .AND. myIter.EQ.1+nIter0 & .AND. nPx.EQ.1 .AND. nPy.EQ.1 & .AND. useCubedSphereExchange ) THEN CALL DEBUG_CS_CORNER_UV( ' afx,afy from GAD_SOM_ADVECT', & afx,afy, k, standardMessageUnit,bi,bj,myThid ) ENDIF #endif /* ALLOW_DEBUG */ C-- End of K loop for horizontal fluxes ENDDO C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| noFlowAcrossSurf = rigidLid .OR. nonlinFreeSurf.GE.1 & .OR. select_rStar.NE.0 IF ( .NOT.implicitAdvection ) THEN C-- Apply limiter (if any): CALL GAD_SOM_LIM_R( bi,bj, limiter, U smVol, U smTr0, U smTr(1-OLx,1-OLy,1,bi,bj,1), U smTr(1-OLx,1-OLy,1,bi,bj,2), U smTr(1-OLx,1-OLy,1,bi,bj,3), U smTr(1-OLx,1-OLy,1,bi,bj,4), U smTr(1-OLx,1-OLy,1,bi,bj,5), U smTr(1-OLx,1-OLy,1,bi,bj,6), U smTr(1-OLx,1-OLy,1,bi,bj,7), U smTr(1-OLx,1-OLy,1,bi,bj,8), U smTr(1-OLx,1-OLy,1,bi,bj,9), I myThid ) C-- Start of k loop for vertical flux DO k=Nr,1,-1 C-- kUp Cycles through 1,2 to point to w-layer above C-- kDown Cycles through 2,1 to point to w-layer below kUp = 1+MOD(Nr-k,2) kDown= 1+MOD(Nr-k+1,2) IF (k.EQ.Nr) THEN C-- Set advective fluxes at the very bottom: DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx alp (i,j,kDown) = 0. _d 0 aln (i,j,kDown) = 0. _d 0 fp_v (i,j,kDown) = 0. _d 0 fn_v (i,j,kDown) = 0. _d 0 fp_o (i,j,kDown) = 0. _d 0 fn_o (i,j,kDown) = 0. _d 0 fp_x (i,j,kDown) = 0. _d 0 fn_x (i,j,kDown) = 0. _d 0 fp_y (i,j,kDown) = 0. _d 0 fn_y (i,j,kDown) = 0. _d 0 fp_z (i,j,kDown) = 0. _d 0 fn_z (i,j,kDown) = 0. _d 0 fp_xx(i,j,kDown) = 0. _d 0 fn_xx(i,j,kDown) = 0. _d 0 fp_yy(i,j,kDown) = 0. _d 0 fn_yy(i,j,kDown) = 0. _d 0 fp_zz(i,j,kDown) = 0. _d 0 fn_zz(i,j,kDown) = 0. _d 0 fp_xy(i,j,kDown) = 0. _d 0 fn_xy(i,j,kDown) = 0. _d 0 fp_xz(i,j,kDown) = 0. _d 0 fn_xz(i,j,kDown) = 0. _d 0 fp_yz(i,j,kDown) = 0. _d 0 fn_yz(i,j,kDown) = 0. _d 0 ENDDO ENDDO ENDIF C-- Compute Vertical transport #ifdef ALLOW_AIM C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr c IF ( k.EQ.1 .OR. c & (useAIM .AND. tracerIdentity.EQ.GAD_SALINITY .AND. k.EQ.Nr) c & ) THEN #else c IF ( k.EQ.1 ) THEN #endif IF ( noFlowAcrossSurf .AND. k.EQ.1 ) THEN C- Surface interface : DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx rTrans(i,j) = 0. maskUp(i,j) = 0. ENDDO ENDDO ELSEIF ( noFlowAcrossSurf ) THEN C- Interior interface : 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) maskUp(i,j) = 1. ENDDO ENDDO ELSE C- Linear Free-Surface: do not mask rTrans : km1= MAX(k-1,1) 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) maskUp(i,j) = maskC(i,j,km1,bi,bj)*maskC(i,j,k,bi,bj) ENDDO ENDDO C- end Surface/Interior if bloc ENDIF C- Compute vertical advective flux in the interior: IF ( vertAdvecScheme.EQ.ENUM_SOM_PRATHER & .OR. vertAdvecScheme.EQ.ENUM_SOM_LIMITER ) THEN CALL GAD_SOM_ADV_R( I bi,bj,k, kUp, kDown, I deltaTLev(k), rTrans, maskUp, I maskInC(1-OLx,1-OLy,bi,bj), U smVol, U smTr0, U smTr(1-OLx,1-OLy,1,bi,bj,1), U smTr(1-OLx,1-OLy,1,bi,bj,2), U smTr(1-OLx,1-OLy,1,bi,bj,3), U smTr(1-OLx,1-OLy,1,bi,bj,4), U smTr(1-OLx,1-OLy,1,bi,bj,5), U smTr(1-OLx,1-OLy,1,bi,bj,6), U smTr(1-OLx,1-OLy,1,bi,bj,7), U smTr(1-OLx,1-OLy,1,bi,bj,8), U smTr(1-OLx,1-OLy,1,bi,bj,9), U alp, aln, fp_v, fn_v, fp_o, fn_o, U fp_x, fn_x, fp_y, fn_y, fp_z, fn_z, U fp_xx, fn_xx, fp_yy, fn_yy, fp_zz, fn_zz, U fp_xy, fn_xy, fp_xz, fn_xz, fp_yz, fn_yz, O afr, myThid ) ELSE STOP 'GAD_SOM_ADVECT: adv. scheme incompatibale with SOM' ENDIF C-- Compute new tracer value and store tracer tendency DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx C-- without rescaling of tendencies: c localTr = smTr0(i,j,k)/smVol(i,j,k) c gTracer(i,j,k) = ( localTr - tracer(i,j,k,bi,bj) ) c & / deltaTLev(k) C-- Non-Lin Free-Surf: consistent with rescaling of tendencies C (in FREESURF_RESCALE_G) and RealFreshFlux/addMass. C Also valid for linear Free-Surf (r & r* coords) except that surf tracer C loss/gain is computed (in GAD_SOM_ADV_R) from partially updated tracer C (instead of from Tr^n as fresh-water dilution effect) resulting in C inaccurate linFSConserveTr and "surfExpan_" monitor. gTracer(i,j,k) = & ( smTr0(i,j,k) - tracer(i,j,k,bi,bj)*smVol(i,j,k) ) & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k) & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj) & *recip_rhoFacC(k) & /deltaTLev(k) ENDDO ENDDO #ifdef ALLOW_DIAGNOSTICS IF ( doDiagAdvR ) THEN diagName = 'ADVr'//diagSufx CALL DIAGNOSTICS_FILL( afr, & diagName, k,1, 2,bi,bj, myThid ) ENDIF #ifdef ALLOW_LAYERS IF ( useLayers ) THEN CALL LAYERS_FILL(afr,tracerIdentity,'AFR', & k,1,2,bi,bj,myThid) ENDIF #endif /* ALLOW_LAYERS */ #endif C-- End of k loop for vertical flux ENDDO C-- end of if not.implicitAdvection block ENDIF RETURN END