C $Header: /u/gcmpack/MITgcm/model/src/calc_gw.F,v 1.55 2014/12/24 19:09:33 jmc Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" #ifdef ALLOW_MOM_COMMON # include "MOM_COMMON_OPTIONS.h" #endif #define CALC_GW_NEW_THICK CBOP C !ROUTINE: CALC_GW C !INTERFACE: SUBROUTINE CALC_GW( I bi, bj, kappaRU, kappaRV, I str13, str23, str33, I viscAh3d_00, viscAh3d_13, viscAh3d_23, I myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R CALC_GW C | o Calculate vertical velocity tendency terms C | ( Non-Hydrostatic only ) C *==========================================================* C | In NH, the vertical momentum tendency must be C | calculated explicitly and included as a source term C | for a 3d pressure eqn. Calculate that term here. C | This routine is not used in HYD calculations. C *==========================================================* C \ev C !USES: IMPLICIT NONE C == Global variables == #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "RESTART.h" #include "SURFACE.h" #include "DYNVARS.h" #include "NH_VARS.h" #ifdef ALLOW_ADDFLUID # include "FFIELDS.h" #endif #ifdef ALLOW_MOM_COMMON # include "MOM_VISC.h" #endif C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == C bi,bj :: current tile indices C kappaRU :: vertical viscosity at U points C kappaRV :: vertical viscosity at V points C myTime :: Current time in simulation C myIter :: Current iteration number in simulation C myThid :: Thread number for this instance of the routine. INTEGER bi,bj _RL kappaRU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1) _RL kappaRV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1) #ifdef ALLOW_SMAG_3D _RL str13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1) _RL str23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1) _RL str33(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL viscAh3d_00(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL viscAh3d_13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1) _RL viscAh3d_23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1) #else /* ALLOW_SMAG_3D */ _RL str13(1), str23(1), str33(1) _RL viscAh3d_00(1), viscAh3d_13(1), viscAh3d_23(1) #endif /* ALLOW_SMAG_3D */ _RL myTime INTEGER myIter INTEGER myThid #ifdef ALLOW_NONHYDROSTATIC #ifdef ALLOW_MOM_COMMON C !LOCAL VARIABLES: C == Local variables == C biharmonicVisc:: use horizontal biharmonic viscosity for vertical momentum C iMin, iMax :: Ranges and sub-block indices on which calculations C jMin, jMax are applied. C xA :: W-Cell face area normal to X C yA :: W-Cell face area normal to Y C rThickC_W :: thickness (in r-units) of W-Cell at Western Edge C rThickC_S :: thickness (in r-units) of W-Cell at Southern Edge C rThickC_C :: thickness (in r-units) of W-Cell (centered on W pt) C recip_rThickC :: reciprol thickness of W-Cell (centered on W-point) C flx_NS :: vertical momentum flux, meridional direction C flx_EW :: vertical momentum flux, zonal direction C flxAdvUp :: vertical mom. advective flux, vertical direction (@ level k-1) C flxDisUp :: vertical mom. dissipation flux, vertical direction (@ level k-1) C flx_Dn :: vertical momentum flux, vertical direction (@ level k) C gwDiss :: vertical momentum dissipation tendency C gwAdd :: other tendencies (Coriolis, Metric-terms) C gw_AB :: tendency increment from Adams-Bashforth C del2w :: laplacian of wVel C wFld :: local copy of wVel C i,j,k :: Loop counters LOGICAL biharmonicVisc INTEGER iMin,iMax,jMin,jMax _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL rThickC_W (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL rThickC_S (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL rThickC_C (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL recip_rThickC(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL flxAdvUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL flxDisUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL gwDiss(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL gwAdd (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL gw_AB (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL del2w (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER i,j,k, km1, kp1 _RL mskM1, mskP1 _RL tmp_WbarZ _RL uTrans, vTrans, rTrans _RL viscLoc PARAMETER( iMin = 1 , iMax = sNx ) PARAMETER( jMin = 1 , jMax = sNy ) CEOP #ifdef ALLOW_DIAGNOSTICS LOGICAL diagDiss, diagAdvec, diag_AB LOGICAL DIAGNOSTICS_IS_ON EXTERNAL #endif /* ALLOW_DIAGNOSTICS */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef ALLOW_DIAGNOSTICS IF ( useDiagnostics ) THEN diagDiss = DIAGNOSTICS_IS_ON( 'Wm_Diss ', myThid ) diagAdvec = DIAGNOSTICS_IS_ON( 'Wm_Advec', myThid ) diag_AB = DIAGNOSTICS_IS_ON( 'AB_gW ', myThid ) ELSE diagDiss = .FALSE. diagAdvec = .FALSE. diag_AB = .FALSE. ENDIF #endif /* ALLOW_DIAGNOSTICS */ biharmonicVisc = viscA4W.NE.zeroRL & .OR. ( useVariableVisc .AND. useBiharmonicVisc ) C-- Initialise gW to zero DO k=1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx gW(i,j,k,bi,bj) = 0. ENDDO ENDDO ENDDO C- Initialise gwDiss to zero DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx gwDiss(i,j) = 0. ENDDO ENDDO IF (momViscosity) THEN C- Initialize del2w to zero: DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx del2w(i,j) = 0. _d 0 ENDDO ENDDO ENDIF C-- Boundaries condition at top (vertical advection of vertical momentum): DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx flxAdvUp(i,j) = 0. c flxDisUp(i,j) = 0. ENDDO ENDDO C--- Sweep down column DO k=1,Nr km1 = MAX( k-1, 1 ) kp1 = MIN( k+1,Nr ) mskM1 = 1. mskP1 = 1. IF ( k.EQ. 1 ) mskM1 = 0. IF ( k.EQ.Nr ) mskP1 = 0. IF ( k.GT.1 ) THEN C-- Compute grid factor arround a W-point: #ifdef CALC_GW_NEW_THICK DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx IF ( maskC(i,j,k-1,bi,bj).EQ.0. .OR. & maskC(i,j, k ,bi,bj).EQ.0. ) THEN recip_rThickC(i,j) = 0. ELSE C- valid in z & p coord.; also accurate if Interface @ middle between 2 centers recip_rThickC(i,j) = 1. _d 0 / & ( MIN( Ro_surf(i,j,bi,bj),rC(k-1) ) & - MAX( R_low(i,j,bi,bj), rC(k) ) & ) ENDIF ENDDO ENDDO IF (momViscosity) THEN DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx rThickC_C(i,j) = MAX( zeroRS, & MIN( Ro_surf(i,j,bi,bj), rC(k-1) ) & -MAX( R_low(i,j,bi,bj), rC(k) ) & ) ENDDO ENDDO DO j=1-OLy,sNy+OLy DO i=1-OLx+1,sNx+OLx rThickC_W(i,j) = MAX( zeroRS, & MIN( rSurfW(i,j,bi,bj), rC(k-1) ) & -MAX( rLowW(i,j,bi,bj), rC(k) ) & ) C W-Cell Western face area: xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j) c & *deepFacF(k) ENDDO ENDDO DO j=1-OLy+1,sNy+OLy DO i=1-OLx,sNx+OLx rThickC_S(i,j) = MAX( zeroRS, & MIN( rSurfS(i,j,bi,bj), rC(k-1) ) & -MAX( rLowS(i,j,bi,bj), rC(k) ) & ) C W-Cell Southern face area: yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j) c & *deepFacF(k) C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC. C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted) ENDDO ENDDO ENDIF #else /* CALC_GW_NEW_THICK */ DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx C- note: assume fluid @ smaller k than bottom: does not work in p-coordinate ! IF ( maskC(i,j,k,bi,bj).EQ.0. ) THEN recip_rThickC(i,j) = 0. ELSE recip_rThickC(i,j) = 1. _d 0 / & ( drF(k-1)*halfRS & + drF( k )*MIN( _hFacC(i,j, k ,bi,bj), halfRS ) & ) ENDIF c IF (momViscosity) THEN #ifdef NONLIN_FRSURF rThickC_C(i,j) = & drF(k-1)*MAX( h0FacC(i,j,k-1,bi,bj)-halfRS, zeroRS ) & + drF( k )*MIN( h0FacC(i,j,k ,bi,bj), halfRS ) #else rThickC_C(i,j) = & drF(k-1)*MAX( _hFacC(i,j,k-1,bi,bj)-halfRS, zeroRS ) & + drF( k )*MIN( _hFacC(i,j,k ,bi,bj), halfRS ) #endif rThickC_W(i,j) = & drF(k-1)*MAX( _hFacW(i,j,k-1,bi,bj)-halfRS, zeroRS ) & + drF( k )*MIN( _hFacW(i,j,k ,bi,bj), halfRS ) rThickC_S(i,j) = & drF(k-1)*MAX( _hFacS(i,j,k-1,bi,bj)-halfRS, zeroRS ) & + drF( k )*MIN( _hFacS(i,j, k ,bi,bj), halfRS ) C W-Cell Western face area: xA(i,j) = _dyG(i,j,bi,bj)*rThickC_W(i,j) c & *deepFacF(k) C W-Cell Southern face area: yA(i,j) = _dxG(i,j,bi,bj)*rThickC_S(i,j) c & *deepFacF(k) C deep-model: xA,yA is only used for viscous flux, in terms like: xA/dxC,yA/dyC. C this gives deepFacF*recip_deepFacF => cancel each other (and therefore omitted) c ENDIF ENDDO ENDDO #endif /* CALC_GW_NEW_THICK */ ELSEIF ( selectNHfreeSurf.GE.1 ) THEN DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx recip_rThickC(i,j) = recip_drC(k) c rThickC_C(i,j) = drC(k) c rThickC_W(i,j) = drC(k) c rThickC_S(i,j) = drC(k) c xA(i,j) = _dyG(i,j,bi,bj)*drC(k) c yA(i,j) = _dxG(i,j,bi,bj)*drC(k) ENDDO ENDDO ENDIF C-- local copy of wVel: DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx wFld(i,j) = wVel(i,j,k,bi,bj) ENDDO ENDDO C-- horizontal bi-harmonic dissipation IF ( momViscosity .AND. k.GT.1 .AND. biharmonicVisc ) THEN C- calculate the horizontal Laplacian of vertical flow C Zonal flux d/dx W IF ( useCubedSphereExchange ) THEN C to compute d/dx(W), fill corners with appropriate values: CALL FILL_CS_CORNER_TR_RL( 1, .FALSE., & wFld, bi,bj, myThid ) ENDIF DO j=1-OLy,sNy+OLy flx_EW(1-OLx,j)=0. DO i=1-OLx+1,sNx+OLx flx_EW(i,j) = & ( wFld(i,j) - wFld(i-1,j) ) & *_recip_dxC(i,j,bi,bj)*xA(i,j) #ifdef COSINEMETH_III & *sqCosFacU(j,bi,bj) #endif #ifdef ALLOW_OBCS & *maskInW(i,j,bi,bj) #endif ENDDO ENDDO C Meridional flux d/dy W IF ( useCubedSphereExchange ) THEN C to compute d/dy(W), fill corners with appropriate values: CALL FILL_CS_CORNER_TR_RL( 2, .FALSE., & wFld, bi,bj, myThid ) ENDIF DO i=1-OLx,sNx+OLx flx_NS(i,1-OLy)=0. ENDDO DO j=1-OLy+1,sNy+OLy DO i=1-OLx,sNx+OLx flx_NS(i,j) = & ( wFld(i,j) - wFld(i,j-1) ) & *_recip_dyC(i,j,bi,bj)*yA(i,j) #ifdef ISOTROPIC_COS_SCALING #ifdef COSINEMETH_III & *sqCosFacV(j,bi,bj) #endif #endif #ifdef ALLOW_OBCS & *maskInS(i,j,bi,bj) #endif ENDDO ENDDO C del^2 W C Divergence of horizontal fluxes DO j=1-OLy,sNy+OLy-1 DO i=1-OLx,sNx+OLx-1 del2w(i,j) = ( ( flx_EW(i+1,j)-flx_EW(i,j) ) & +( flx_NS(i,j+1)-flx_NS(i,j) ) & )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j) & *recip_deepFac2F(k) ENDDO ENDDO C end if biharmonic viscosity ENDIF IF ( momViscosity .AND. k.GT.1 ) THEN C Viscous Flux on Western face DO j=jMin,jMax DO i=iMin,iMax+1 flx_EW(i,j)= & - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i-1,j,k,bi,bj))*halfRL & *(wVel(i,j,k,bi,bj)-wVel(i-1,j,k,bi,bj)) & *_recip_dxC(i,j,bi,bj)*xA(i,j) & *cosFacU(j,bi,bj) & + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i-1,j,k,bi,bj))*halfRL & *(del2w(i,j)-del2w(i-1,j)) & *_recip_dxC(i,j,bi,bj)*xA(i,j) #ifdef COSINEMETH_III & *sqCosFacU(j,bi,bj) #else & *cosFacU(j,bi,bj) #endif ENDDO ENDDO C Viscous Flux on Southern face DO j=jMin,jMax+1 DO i=iMin,iMax flx_NS(i,j)= & - (viscAh_W(i,j,k,bi,bj)+viscAh_W(i,j-1,k,bi,bj))*halfRL & *(wVel(i,j,k,bi,bj)-wVel(i,j-1,k,bi,bj)) & *_recip_dyC(i,j,bi,bj)*yA(i,j) #ifdef ISOTROPIC_COS_SCALING & *cosFacV(j,bi,bj) #endif & + (viscA4_W(i,j,k,bi,bj)+viscA4_W(i,j-1,k,bi,bj))*halfRL & *(del2w(i,j)-del2w(i,j-1)) & *_recip_dyC(i,j,bi,bj)*yA(i,j) #ifdef ISOTROPIC_COS_SCALING #ifdef COSINEMETH_III & *sqCosFacV(j,bi,bj) #else & *cosFacV(j,bi,bj) #endif #endif ENDDO ENDDO C Viscous Flux on Lower face of W-Cell (= at tracer-cell center, level k) DO j=jMin,jMax DO i=iMin,iMax C Interpolate vert viscosity to center of tracer-cell (level k): viscLoc = ( kappaRU(i,j,k) +kappaRU(i+1,j,k) & +kappaRU(i,j,k+1)+kappaRU(i+1,j,k+1) & +kappaRV(i,j,k) +kappaRV(i,j+1,k) & +kappaRV(i,j,k+1)+kappaRV(i,j+1,k+1) & )*0.125 _d 0 flx_Dn(i,j) = & - viscLoc*( wVel(i,j,kp1,bi,bj)*mskP1 & -wVel(i,j, k ,bi,bj) )*rkSign & *recip_drF(k)*rA(i,j,bi,bj) & *deepFac2C(k)*rhoFacC(k) ENDDO ENDDO IF ( k.EQ.2 ) THEN C Viscous Flux on Upper face of W-Cell (= at tracer-cell center, level k-1) DO j=jMin,jMax DO i=iMin,iMax C Interpolate horizontally (but not vertically) vert viscosity to center: C Although background visc. might be defined at k=1, this is not C generally true when using variable visc. (from vertical mixing scheme). C Therefore, no vert. interp. and only horizontal interpolation. viscLoc = ( kappaRU(i,j,k) + kappaRU(i+1,j,k) & +kappaRV(i,j,k) + kappaRV(i,j+1,k) & )*0.25 _d 0 flxDisUp(i,j) = & - viscLoc*( wVel(i,j, k ,bi,bj) & -wVel(i,j,k-1,bi,bj) )*rkSign & *recip_drF(k-1)*rA(i,j,bi,bj) & *deepFac2C(k-1)*rhoFacC(k-1) C to recover old (before 2009/11/30) results (since flxDisUp(k=2) was zero) c flxDisUp(i,j) = 0. ENDDO ENDDO ENDIF C Tendency is minus divergence of viscous fluxes: C anelastic: vert.visc.flx is scaled by rhoFac but hor.visc.fluxes are not DO j=jMin,jMax DO i=iMin,iMax gwDiss(i,j) = & -( ( flx_EW(i+1,j)-flx_EW(i,j) ) & + ( flx_NS(i,j+1)-flx_NS(i,j) ) & + ( flx_Dn(i,j)-flxDisUp(i,j) )*rkSign & *recip_rhoFacF(k) & )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j) & *recip_deepFac2F(k) C-- prepare for next level (k+1) flxDisUp(i,j)=flx_Dn(i,j) ENDDO ENDDO ENDIF IF ( momViscosity .AND. k.GT.1 .AND. no_slip_sides ) THEN C- No-slip BCs impose a drag at walls... CALL MOM_W_SIDEDRAG( I bi,bj,k, I wVel, del2w, I rThickC_C, recip_rThickC, I viscAh_W, viscA4_W, O gwAdd, I myThid ) DO j=jMin,jMax DO i=iMin,iMax gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j) ENDDO ENDDO ENDIF #ifdef ALLOW_SMAG_3D IF ( useSmag3D .AND. k.GT.1 ) THEN CALL MOM_W_SMAG_3D( I str13, str23, str33, I viscAh3d_00, viscAh3d_13, viscAh3d_23, I rThickC_W, rThickC_S, rThickC_C, recip_rThickC, O gwAdd, I iMin,iMax,jMin,jMax, k, bi, bj, myThid ) DO j = jMin,jMax DO i = iMin,iMax gwDiss(i,j) = gwDiss(i,j) + gwAdd(i,j) ENDDO ENDDO ENDIF #endif /* ALLOW_SMAG_3D */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| IF ( momAdvection ) THEN IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN C Advective Flux on Western face DO j=jMin,jMax DO i=iMin,iMax+1 C transport through Western face area: uTrans = ( & drF(km1)*_hFacW(i,j,km1,bi,bj)*uVel(i,j,km1,bi,bj) & *rhoFacC(km1)*mskM1 & + drF( k )*_hFacW(i,j, k ,bi,bj)*uVel(i,j, k ,bi,bj) & *rhoFacC(k) & )*halfRL*_dyG(i,j,bi,bj)*deepFacF(k) flx_EW(i,j) = uTrans*(wFld(i,j)+wFld(i-1,j))*halfRL c flx_EW(i,j)= c & uTrans*(wVel(i,j,k,bi,bj)+wVel(i-1,j,k,bi,bj))*halfRL ENDDO ENDDO C Advective Flux on Southern face DO j=jMin,jMax+1 DO i=iMin,iMax C transport through Southern face area: vTrans = ( & drF(km1)*_hFacS(i,j,km1,bi,bj)*vVel(i,j,km1,bi,bj) & *rhoFacC(km1)*mskM1 & +drF( k )*_hFacS(i,j, k ,bi,bj)*vVel(i,j, k ,bi,bj) & *rhoFacC(k) & )*halfRL*_dxG(i,j,bi,bj)*deepFacF(k) flx_NS(i,j) = vTrans*(wFld(i,j)+wFld(i,j-1))*halfRL c flx_NS(i,j)= c & vTrans*(wVel(i,j,k,bi,bj)+wVel(i,j-1,k,bi,bj))*halfRL ENDDO ENDDO ENDIF C Advective Flux on Lower face of W-Cell (= at tracer-cell center, level k) c IF (.TRUE.) THEN DO j=jMin,jMax DO i=iMin,iMax C NH in p-coord.: advect wSpeed [m/s] with rTrans tmp_WbarZ = halfRL* & ( wVel(i,j, k ,bi,bj)*rVel2wUnit( k ) & +wVel(i,j,kp1,bi,bj)*rVel2wUnit(kp1)*mskP1 ) C transport through Lower face area: rTrans = halfRL* & ( wVel(i,j, k ,bi,bj)*deepFac2F( k )*rhoFacF( k ) & +wVel(i,j,kp1,bi,bj)*deepFac2F(kp1)*rhoFacF(kp1) & *mskP1 & )*rA(i,j,bi,bj) flx_Dn(i,j) = rTrans*tmp_WbarZ ENDDO ENDDO c ENDIF IF ( k.EQ.1 .AND. selectNHfreeSurf.GE.1 ) THEN C Advective Flux on Upper face of W-Cell (= at surface) DO j=jMin,jMax DO i=iMin,iMax tmp_WbarZ = wVel(i,j,k,bi,bj)*rVel2wUnit(k) rTrans = wVel(i,j,k,bi,bj)*deepFac2F(k)*rhoFacF(k) & *rA(i,j,bi,bj) flxAdvUp(i,j) = rTrans*tmp_WbarZ c flxAdvUp(i,j) = 0. ENDDO ENDDO ENDIF IF ( k.GT.1 .OR. selectNHfreeSurf.GE.1 ) THEN C Tendency is minus divergence of advective fluxes: C anelastic: all transports & advect. fluxes are scaled by rhoFac DO j=jMin,jMax DO i=iMin,iMax C to recover old (before 2009/11/30) results (since flxAdvUp(k=2) was zero) c IF (k.EQ.2) flxAdvUp(i,j) = 0. gW(i,j,k,bi,bj) = & -( ( flx_EW(i+1,j)-flx_EW(i,j) ) & + ( flx_NS(i,j+1)-flx_NS(i,j) ) & + ( flx_Dn(i,j)-flxAdvUp(i,j) )*rkSign*wUnit2rVel(k) & )*recip_rA(i,j,bi,bj)*recip_rThickC(i,j) & *recip_deepFac2F(k)*recip_rhoFacF(k) ENDDO ENDDO #ifdef ALLOW_ADDFLUID IF ( selectAddFluid.GE.1 ) THEN DO j=jMin,jMax DO i=iMin,iMax gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj) & + wVel(i,j,k,bi,bj)*mass2rUnit*0.5 _d 0 & *( addMass(i,j,k,bi,bj) & +addMass(i,j,km1,bi,bj)*mskM1 ) & *recip_rA(i,j,bi,bj)*recip_rThickC(i,j) & *recip_deepFac2F(k)*recip_rhoFacF(k) ENDDO ENDDO ENDIF #endif /* ALLOW_ADDFLUID */ ENDIF DO j=jMin,jMax DO i=iMin,iMax C-- prepare for next level (k+1) flxAdvUp(i,j)=flx_Dn(i,j) ENDDO ENDDO c ELSE C- if momAdvection / else c DO j=1-OLy,sNy+OLy c DO i=1-OLx,sNx+OLx c gW(i,j,k,bi,bj) = 0. _d 0 c ENDDO c ENDDO C- endif momAdvection. ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| IF ( useNHMTerms .AND. k.GT.1 ) THEN CALL MOM_W_METRIC_NH( I bi,bj,k, I uVel, vVel, O gwAdd, I myThid ) DO j=jMin,jMax DO i=iMin,iMax gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j) ENDDO ENDDO ENDIF IF ( use3dCoriolis .AND. k.GT.1 ) THEN CALL MOM_W_CORIOLIS_NH( I bi,bj,k, I uVel, vVel, O gwAdd, I myThid ) DO j=jMin,jMax DO i=iMin,iMax gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j) ENDDO ENDDO ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #ifdef ALLOW_DIAGNOSTICS IF ( diagDiss ) THEN CALL DIAGNOSTICS_FILL( gwDiss, 'Wm_Diss ', & k, 1, 2, bi,bj, myThid ) C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL C does it only if k=1 (never the case here) c IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Diss ',bi,bj,myThid) ENDIF IF ( diagAdvec ) THEN CALL DIAGNOSTICS_FILL( gW, 'Wm_Advec', & k,Nr, 1, bi,bj, myThid ) c IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Advec',bi,bj,myThid) ENDIF #endif /* ALLOW_DIAGNOSTICS */ C-- Dissipation term inside the Adams-Bashforth: IF ( momViscosity .AND. momDissip_In_AB) THEN DO j=jMin,jMax DO i=iMin,iMax gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j) ENDDO ENDDO ENDIF C- Compute effective gW_[n+1/2] terms (including Adams-Bashforth weights) C and save gW_[n] into gwNm1 for the next time step. #ifdef ALLOW_ADAMSBASHFORTH_3 CALL ADAMS_BASHFORTH3( I bi, bj, k, Nr, U gW(1-OLx,1-OLy,1,bi,bj), gwNm, O gw_AB, I nHydStartAB, myIter, myThid ) #else /* ALLOW_ADAMSBASHFORTH_3 */ CALL ADAMS_BASHFORTH2( I bi, bj, k, Nr, U gW(1-OLx,1-OLy,1,bi,bj), U gwNm1(1-OLx,1-OLy,1,bi,bj), O gw_AB, I nHydStartAB, myIter, myThid ) #endif /* ALLOW_ADAMSBASHFORTH_3 */ #ifdef ALLOW_DIAGNOSTICS IF ( diag_AB ) THEN CALL DIAGNOSTICS_FILL(gw_AB,'AB_gW ',k,1,2,bi,bj,myThid) ENDIF #endif /* ALLOW_DIAGNOSTICS */ C-- Dissipation term outside the Adams-Bashforth: IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN DO j=jMin,jMax DO i=iMin,iMax gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwDiss(i,j) ENDDO ENDDO ENDIF C- end of the k loop ENDDO #ifdef ALLOW_DIAGNOSTICS IF (useDiagnostics) THEN CALL DIAGNOSTICS_FILL(viscAh_W,'VISCAHW ',0,Nr,1,bi,bj,myThid) CALL DIAGNOSTICS_FILL(viscA4_W,'VISCA4W ',0,Nr,1,bi,bj,myThid) ENDIF #endif /* ALLOW_DIAGNOSTICS */ #endif /* ALLOW_MOM_COMMON */ #endif /* ALLOW_NONHYDROSTATIC */ RETURN END