diff --git a/model/src/calc_gw.F b/model/src/calc_gw.F index 269abc0ac6..b055a84381 100644 --- a/model/src/calc_gw.F +++ b/model/src/calc_gw.F @@ -122,7 +122,7 @@ SUBROUTINE CALC_GW( PARAMETER( jMin = 1 , jMax = sNy ) CEOP #ifdef ALLOW_DIAGNOSTICS - LOGICAL diagDiss, diagAdvec, diag_AB + LOGICAL diagDiss, diagAdvec, diagMetric, diag_AB LOGICAL DIAGNOSTICS_IS_ON EXTERNAL DIAGNOSTICS_IS_ON #endif /* ALLOW_DIAGNOSTICS */ @@ -133,10 +133,12 @@ SUBROUTINE CALC_GW( IF ( useDiagnostics ) THEN diagDiss = DIAGNOSTICS_IS_ON( 'Wm_Diss ', myThid ) diagAdvec = DIAGNOSTICS_IS_ON( 'Wm_Advec', myThid ) + diagMetric= DIAGNOSTICS_IS_ON( 'Wm_Metr ', myThid ) diag_AB = DIAGNOSTICS_IS_ON( 'AB_gW ', myThid ) ELSE diagDiss = .FALSE. diagAdvec = .FALSE. + diagMetric= .FALSE. diag_AB = .FALSE. ENDIF #endif /* ALLOW_DIAGNOSTICS */ @@ -604,7 +606,7 @@ SUBROUTINE CALC_GW( IF ( useNHMTerms .AND. k.GT.1 ) THEN CALL MOM_W_METRIC_NH( - I bi,bj,k, + I bi, bj, k, I uVel, vVel, O gwAdd, I myThid ) @@ -613,11 +615,20 @@ SUBROUTINE CALC_GW( gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)+gwAdd(i,j) ENDDO ENDDO +#ifdef ALLOW_DIAGNOSTICS + IF ( diagMetric ) THEN + CALL DIAGNOSTICS_FILL( gwAdd, 'Wm_Metr ', + I k,1,2, bi,bj, myThid ) +C- note: need to explicitly increment the counter since DIAGNOSTICS_FILL +C does it only if k=1 (never the case here) + IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT('Wm_Metr ',bi,bj,myThid) + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ ENDIF IF ( select3dCoriScheme.GE.1 .AND. k.GT.1 ) THEN CALL MOM_W_CORIOLIS_NH( - I bi,bj,k, - I uVel, vVel, + I bi, bj, k, + I uVel, vVel, recip_rThickC, O gwAdd, I myThid ) DO j=jMin,jMax @@ -633,8 +644,6 @@ SUBROUTINE CALC_GW( 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 diff --git a/model/src/set_parms.F b/model/src/set_parms.F index efe92d7af6..420546f4c7 100644 --- a/model/src/set_parms.F +++ b/model/src/set_parms.F @@ -63,8 +63,13 @@ SUBROUTINE SET_PARMS( myThid ) & SQUEEZE_RIGHT, myThid ) ENDIF IF ( usingCurvilinearGrid ) THEN +C- Horizontal metric terms not implemented for Curvilinear-Grid selectMetricTerms = 0 ENDIF + IF ( usingCylindricalGrid ) THEN +C- No current alternative metric terms formulation for Cylindrical-Grid + selectMetricTerms = MIN( selectMetricTerms, 1 ) + ENDIF IF ( selectCoriMap.EQ.-1 ) THEN IF ( usingCartesianGrid.OR.usingCylindricalGrid ) THEN C default is to use Beta-Plane Coriolis @@ -92,6 +97,8 @@ SUBROUTINE SET_PARMS( myThid ) momPressureForcing= momStepping .AND. momPressureForcing implicitIntGravWave=momPressureForcing .AND. implicitIntGravWave momImplVertAdv = momAdvection .AND. momImplVertAdv + useNHMTerms = momAdvection .AND. useNHMTerms + IF ( .NOT.momAdvection ) selectMetricTerms = 0 implicitViscosity= momViscosity .AND. implicitViscosity useSmag3D = momViscosity .AND. useSmag3D use3Dsolver = nonHydrostatic.OR. implicitIntGravWave diff --git a/pkg/mom_common/MOM_COMMON_OPTIONS.h b/pkg/mom_common/MOM_COMMON_OPTIONS.h index d295357b8f..51c25aff0e 100644 --- a/pkg/mom_common/MOM_COMMON_OPTIONS.h +++ b/pkg/mom_common/MOM_COMMON_OPTIONS.h @@ -26,5 +26,10 @@ C as a function of grid cell thickness and roughness length C zRoughBot (order 0.01m), assuming a von Karman constant = 0.4. #undef ALLOW_BOTTOMDRAG_ROUGHNESS +C-- Deep-Model: use the original vertical advection (with all NH metric terms) +C instead of updated version that advects the product deepFac x (u,v) +C thus removing the need for NH-metric terms: w.(u,v)/r +#undef MOM_USE_OLD_DEEP_VERT_ADV + #endif /* ALLOW_MOM_COMMON */ #endif /* MOM_COMMON_OPTIONS_H */ diff --git a/pkg/mom_common/MOM_VISC.h b/pkg/mom_common/MOM_VISC.h index 3a7905c95f..ee0cb7991f 100644 --- a/pkg/mom_common/MOM_VISC.h +++ b/pkg/mom_common/MOM_VISC.h @@ -20,6 +20,12 @@ C useVariableVisc :: variable (in space or time) viscosity is used _RL L4rdt_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL L4rdt_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) +C-- COMMON /MOM_FLUXFORM_GRID/ just to hold a copy of deepFacC +C deepFacAdv :: copy of deepFacC or just Nr-length vector of 1 when +C MOM_USE_OLD_DEEP_VERT_ADV is defined or useNHMTerms=F + COMMON /MOM_GRID_COPY/ deepFacAdv + _RL deepFacAdv(Nr) + #ifdef ALLOW_SMAG_3D C smag3D_hLsC :: horiz. grid length scale (power 2/3) at grid cell center C smag3D_hLsW :: horiz. grid length scale (power 2/3) at western edge @@ -64,7 +70,7 @@ C viscA4_W :: Horizontal biharmonic viscosity for vertical momentum #ifdef ALLOW_BOTTOMDRAG_ROUGHNESS C-- bottom drag coefficents as a function of grid cell thickness C and roughness length - COMMON /GRID_DRAGCOEFFS_RS/ + COMMON /MOM_DRAGCOEFFS_RS/ & bottomDragCoeffW, bottomDragCoeffS _RS bottomDragCoeffW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RS bottomDragCoeffS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) diff --git a/pkg/mom_common/mom_diagnostics_init.F b/pkg/mom_common/mom_diagnostics_init.F index c87f664d47..90e3fa2f09 100644 --- a/pkg/mom_common/mom_diagnostics_init.F +++ b/pkg/mom_common/mom_diagnostics_init.F @@ -413,6 +413,20 @@ SUBROUTINE MOM_DIAGNOSTICS_INIT( myThid ) #endif /* ALLOW_MOM_VECINV */ +C- diagnostics of metric-terms tendencies: + diagName = 'Um_Metr ' + diagTitle = 'U momentum tendency from Metric terms' + diagCode = 'UUR MR ' + diagMate = diagNum + 2 + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + diagName = 'Vm_Metr ' + diagTitle = 'V momentum tendency from Metric terms' + diagCode = 'VVR MR ' + diagMate = diagNum + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, diagMate, myThid ) + #ifdef ALLOW_NONHYDROSTATIC C- vertical momentum tendencies IF ( usingPCoords ) THEN @@ -420,6 +434,11 @@ SUBROUTINE MOM_DIAGNOSTICS_INIT( myThid ) ELSE diagUnits = 'm/s^2 ' ENDIF + diagName = 'WSidDrag' + diagTitle = 'Vertical momentum tendency from Side Drag' + diagCode = 'WMr LR ' + CALL DIAGNOSTICS_ADDTOLIST( diagNum, + I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) diagName = 'Wm_Diss ' diagTitle = 'W momentum tendency from Dissipation' diagCode = 'WMr LR ' @@ -431,9 +450,8 @@ SUBROUTINE MOM_DIAGNOSTICS_INIT( myThid ) diagCode = 'WMr LR ' CALL DIAGNOSTICS_ADDTOLIST( diagNum, I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) - - diagName = 'WSidDrag' - diagTitle = 'Vertical momentum tendency from Side Drag' + diagName = 'Wm_Metr ' + diagTitle = 'W momentum tendency from Metric terms' diagCode = 'WMr LR ' CALL DIAGNOSTICS_ADDTOLIST( diagNum, I diagName, diagCode, diagUnits, diagTitle, 0, myThid ) diff --git a/pkg/mom_common/mom_init_fixed.F b/pkg/mom_common/mom_init_fixed.F index 9b3000b581..61e4dafb9d 100644 --- a/pkg/mom_common/mom_init_fixed.F +++ b/pkg/mom_common/mom_init_fixed.F @@ -38,11 +38,24 @@ SUBROUTINE MOM_INIT_FIXED( myThid ) CHARACTER*(MAX_LEN_MBUF) msgBuf #endif - k = 1 twoThird = 2. _d 0 / 3. _d 0 recip_dt = 1. _d 0 IF ( deltaTMom.NE.0. ) recip_dt = 1. _d 0/deltaTMom + _BEGIN_MASTER(myThid) + IF ( useNHMTerms ) THEN + DO k=1,Nr + deepFacAdv(k) = deepFacC(k) + ENDDO +#ifndef MOM_USE_OLD_DEEP_VERT_ADV + ELSE +#endif + DO k=1,Nr + deepFacAdv(k) = 1. _d 0 + ENDDO + ENDIF + _END_MASTER(myThid) + DO bj=myByLo(myThid), myByHi(myThid) DO bi=myBxLo(myThid), myBxHi(myThid) diff --git a/pkg/mom_common/mom_u_coriolis_nh.F b/pkg/mom_common/mom_u_coriolis_nh.F index 1c244fe13f..692e44a1b4 100644 --- a/pkg/mom_common/mom_u_coriolis_nh.F +++ b/pkg/mom_common/mom_u_coriolis_nh.F @@ -5,7 +5,7 @@ C !INTERFACE: ========================================================== SUBROUTINE MOM_U_CORIOLIS_NH( - I bi,bj,k,wFld, + I bi, bj, k, wFld, O uCoriolisTerm, I myThid ) @@ -24,11 +24,11 @@ SUBROUTINE MOM_U_CORIOLIS_NH( #include "GRID.h" C !INPUT PARAMETERS: =================================================== -C bi,bj :: tile indices +C bi, bj :: tile indices C k :: vertical level C wFld :: vertical flow C myThid :: thread number - INTEGER bi,bj,k + INTEGER bi, bj, k _RL wFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) INTEGER myThid @@ -37,29 +37,48 @@ SUBROUTINE MOM_U_CORIOLIS_NH( _RL uCoriolisTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy) C !LOCAL VARIABLES: ==================================================== -C i,j :: loop indices - INTEGER i,j,kp1 +C i, j :: loop indices + INTEGER i, j, kp1 _RL wMsk CEOP + kp1 = MIN(k+1,Nr) + wMsk = 1. + IF (k.EQ.Nr) wMsk = 0. - kp1=min(k+1,Nr) - wMsk=1. - IF (k.EQ.Nr) wMsk=0. - -C Energy conserving discretization of 2*Omega*cos(phi)*w - DO j=1-Oly,sNy+Oly - DO i=2-Olx,sNx+Olx - uCoriolisTerm(i,j) = - & 0.5*( fCoriCos( i ,j,bi,bj)*angleCosC( i ,j,bi,bj) - & *0.5*( wFld( i ,j, k ,bi,bj)*rVel2wUnit( k ) - & +wFld( i ,j,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk ) - & + fCoriCos(i-1,j,bi,bj)*angleCosC(i-1,j,bi,bj) - & *0.5*( wFld(i-1,j, k ,bi,bj)*rVel2wUnit( k ) - & +wFld(i-1,j,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk ) - & )*gravitySign + IF ( select3dCoriScheme.EQ.1 ) THEN +C- Original discretization of 2*Omega*cos(phi)*w +C documented as Energy conserving + DO j=1-OLy,sNy+OLy + DO i=2-OLx,sNx+OLx + uCoriolisTerm(i,j) = gravitySign*halfRL + & *( fCoriCos( i ,j,bi,bj)*angleCosC( i ,j,bi,bj)*halfRL + & *( wFld( i ,j, k ,bi,bj)*rVel2wUnit( k ) + & + wFld( i ,j,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk ) + & + fCoriCos(i-1,j,bi,bj)*angleCosC(i-1,j,bi,bj)*halfRL + & *( wFld(i-1,j, k ,bi,bj)*rVel2wUnit( k ) + & + wFld(i-1,j,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk ) + & ) + ENDDO + ENDDO + ELSE +C- Using averaged transport: + DO j=1-OLy,sNy+OLy + DO i=2-OLx,sNx+OLx + uCoriolisTerm(i,j) = gravitySign*halfRL + & *( fCoriCos( i ,j,bi,bj)*angleCosC( i ,j,bi,bj) + & *( wFld( i ,j, k ,bi,bj)*rVel2wUnit( k )*deepFac2F( k ) + & + wFld( i ,j,kp1,bi,bj)*rVel2wUnit(kp1)*deepFac2F(kp1)*wMsk + & )*rA( i ,j,bi,bj)*halfRL + & + fCoriCos(i-1,j,bi,bj)*angleCosC(i-1,j,bi,bj) + & *( wFld(i-1,j, k ,bi,bj)*rVel2wUnit( k )*deepFac2F( k ) + & + wFld(i-1,j,kp1,bi,bj)*rVel2wUnit(kp1)*deepFac2F(kp1)*wMsk + & )*rA(i-1,j,bi,bj)*halfRL + & )*recip_rAw(i,j,bi,bj)*recip_deepFac2C(k) + & *recip_hFacW(i,j,k,bi,bj) + ENDDO ENDDO - ENDDO + ENDIF RETURN END diff --git a/pkg/mom_common/mom_v_coriolis_nh.F b/pkg/mom_common/mom_v_coriolis_nh.F index caca1c55a4..3019177ab2 100644 --- a/pkg/mom_common/mom_v_coriolis_nh.F +++ b/pkg/mom_common/mom_v_coriolis_nh.F @@ -5,7 +5,7 @@ C !INTERFACE: ========================================================== SUBROUTINE MOM_V_CORIOLIS_NH( - I bi,bj,k,wFld, + I bi, bj, k, wFld, O vCoriolisTerm, I myThid ) @@ -24,11 +24,11 @@ SUBROUTINE MOM_V_CORIOLIS_NH( #include "GRID.h" C !INPUT PARAMETERS: =================================================== -C bi,bj :: tile indices +C bi, bj :: tile indices C k :: vertical level C wFld :: vertical flow C myThid :: thread number - INTEGER bi,bj,k + INTEGER bi, bj, k _RL wFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) INTEGER myThid @@ -37,29 +37,48 @@ SUBROUTINE MOM_V_CORIOLIS_NH( _RL vCoriolisTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy) C !LOCAL VARIABLES: ==================================================== -C i,j :: loop indices - INTEGER i,j,kp1 +C i, j :: loop indices + INTEGER i, j, kp1 _RL wMsk CEOP + kp1 = MIN(k+1,Nr) + wMsk = 1. + IF (k.EQ.Nr) wMsk = 0. - kp1=min(k+1,Nr) - wMsk=1. - IF (k.EQ.Nr) wMsk=0. - -C Energy conserving discretization of 2*Omega*cos(phi)*w - DO j=2-Oly,sNy+Oly - DO i=1-Olx,sNx+Olx - vCoriolisTerm(i,j) = - & -0.5*( fCoriCos( i ,j,bi,bj)*angleSinC(i, j ,bi,bj) - & *0.5*( wFld( i ,j, k ,bi,bj)*rVel2wUnit( k ) - & +wFld( i ,j,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk ) - & + fCoriCos(i,j-1,bi,bj)*angleSinC(i,j-1,bi,bj) - & *0.5*( wFld(i,j-1, k ,bi,bj)*rVel2wUnit( k ) - & +wFld(i,j-1,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk ) - & )*gravitySign + IF ( select3dCoriScheme.EQ.1 ) THEN +C- Original discretization of 2*Omega*cos(phi)*w +C documented as Energy conserving + DO j=2-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + vCoriolisTerm(i,j) = -gravitySign*halfRL + & *( fCoriCos(i, j ,bi,bj)*angleSinC(i, j ,bi,bj)*halfRL + & *( wFld(i, j , k ,bi,bj)*rVel2wUnit( k ) + & + wFld(i, j ,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk ) + & + fCoriCos(i,j-1,bi,bj)*angleSinC(i,j-1,bi,bj)*halfRL + & *( wFld(i,j-1, k ,bi,bj)*rVel2wUnit( k ) + & + wFld(i,j-1,kp1,bi,bj)*rVel2wUnit(kp1)*wMsk ) + & ) + ENDDO + ENDDO + ELSE +C- Using averaged transport: + DO j=2-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + vCoriolisTerm(i,j) = -gravitySign*halfRL + & *( fCoriCos(i, j ,bi,bj)*angleSinC(i, j ,bi,bj) + & *( wFld(i, j , k ,bi,bj)*rVel2wUnit( k )*deepFac2F( k ) + & + wFld(i, j ,kp1,bi,bj)*rVel2wUnit(kp1)*deepFac2F(kp1)*wMsk + & )*rA(i, j ,bi,bj)*halfRL + & + fCoriCos(i,j-1,bi,bj)*angleSinC(i,j-1,bi,bj) + & *( wFld(i,j-1, k ,bi,bj)*rVel2wUnit( k )*deepFac2F( k ) + & + wFld(i,j-1,kp1,bi,bj)*rVel2wUnit(kp1)*deepFac2F(kp1)*wMsk + & )*rA(i,j-1,bi,bj)*halfRL + & )*recip_rAs(i,j,bi,bj)*recip_deepFac2C(k) + & *recip_hFacS(i,j,k,bi,bj) + ENDDO ENDDO - ENDDO + ENDIF RETURN END diff --git a/pkg/mom_common/mom_w_coriolis_nh.F b/pkg/mom_common/mom_w_coriolis_nh.F index ef8a2aa443..994fc1f8ee 100644 --- a/pkg/mom_common/mom_w_coriolis_nh.F +++ b/pkg/mom_common/mom_w_coriolis_nh.F @@ -5,8 +5,8 @@ C !INTERFACE: ========================================================== SUBROUTINE MOM_W_CORIOLIS_NH( - I bi,bj,k, - I uFld, vFld, + I bi, bj, k, + I uFld, vFld, recip_rThickC, U wCoriolisTerm, I myThid ) @@ -24,14 +24,17 @@ SUBROUTINE MOM_W_CORIOLIS_NH( #include "GRID.h" C !INPUT PARAMETERS: =================================================== -C bi,bj :: tile indices +C bi, bj :: tile indices C k :: vertical level C uFld :: horizontal flow, u component C vFld :: horizontal flow, v component +C rThickC :: thickness of W-cell (r units) +C recip_rThickC :: reciprocal of W-cell thickness C myThid :: my Thread Id number - INTEGER bi,bj,k + INTEGER bi, bj, k _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) + _RL recip_rThickC(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER myThid C !OUTPUT PARAMETERS: ================================================== @@ -40,15 +43,17 @@ SUBROUTINE MOM_W_CORIOLIS_NH( #ifdef ALLOW_NONHYDROSTATIC C !LOCAL VARIABLES: ==================================================== -C i,j :: loop indices +C i, j :: loop indices INTEGER i,j CEOP - -C Energy conserving discretization of 2*Omega*cos(phi)*u_eastward IF ( k.GT.1 .AND. k.LE.Nr ) THEN - DO j=1-Oly,sNy+Oly-1 - DO i=1-Olx,sNx+Olx-1 + + IF ( select3dCoriScheme.EQ.1 ) THEN +C- Original discretization of 2*Omega*cos(phi)*u_eastward +C documented as Energy conserving + DO j=1-OLy,sNy+OLy-1 + DO i=1-OLx,sNx+OLx-1 wCoriolisTerm(i,j) = & -gravitySign*fCoriCos(i,j,bi,bj)* & ( angleCosC(i,j,bi,bj)*( @@ -62,9 +67,63 @@ SUBROUTINE MOM_W_CORIOLIS_NH( & )*wUnit2rVel(k) ENDDO ENDDO + ELSEIF ( select3dCoriScheme.EQ.2 ) THEN +C- Using thickness-averaged transport (but without hFac): + DO j=1-OLy,sNy+OLy-1 + DO i=1-OLx,sNx+OLx-1 + wCoriolisTerm(i,j) = -gravitySign + & *fCoriCos(i,j,bi,bj) + & *( angleCosC(i,j,bi,bj) + & *( ( uFld(i,j,k-1,bi,bj) + uFld(i+1,j,k-1,bi,bj) ) + & *drF(k-1) + & + ( uFld(i,j, k ,bi,bj) + uFld(i+1,j, k ,bi,bj) ) + & *drF( k ) + & )*0.25 _d 0 + & -angleSinC(i,j,bi,bj) + & *( ( vFld(i,j,k-1,bi,bj) + vFld(i,j+1,k-1,bi,bj) ) + & *drF(k-1) + & + ( vFld(i,j, k ,bi,bj) + vFld(i,j+1, k ,bi,bj) ) + & *drF( k ) + & )*0.25 _d 0 + & )*recip_drC(k)*wUnit2rVel(k) + ENDDO + ENDDO + ELSE +C- Using thickness-averaged transport: +C for now, without dyG*deepFacC weight and without recip_hFacI factor: + DO j=1-OLy,sNy+OLy-1 + DO i=1-OLx,sNx+OLx-1 + wCoriolisTerm(i,j) = -gravitySign + & *fCoriCos(i,j,bi,bj) + & *( angleCosC(i,j,bi,bj) + & *( ( uFld( i ,j,k-1,bi,bj)*hFacW( i ,j,k-1,bi,bj) + & + uFld(i+1,j,k-1,bi,bj)*hFacW(i+1,j,k-1,bi,bj) + & )*drF(k-1) +c & *deepFacC(k-1) + & + ( uFld( i ,j, k ,bi,bj)*hFacW( i ,j, k ,bi,bj) + & + uFld(i+1,j, k ,bi,bj)*hFacW(i+1,j, k ,bi,bj) + & )*drF( k ) +c & *deepFacC( k ) + & )*0.25 _d 0 + & -angleSinC(i,j,bi,bj) + & *( ( vFld(i, j ,k-1,bi,bj)*hFacS(i, j ,k-1,bi,bj) + & + vFld(i,j+1,k-1,bi,bj)*hFacS(i,j+1,k-1,bi,bj) + & )*drF(k-1) +c & *deepFacC(k-1) + & + ( vFld(i, j , k ,bi,bj)*hFacS(i, j , k ,bi,bj) + & + vFld(i,j+1, k ,bi,bj)*hFacS(i,j+1, k ,bi,bj) + & )*drF( k ) +c & *deepFacC( k ) + & )*0.25 _d 0 + & )*recip_rThickC(i,j)*wUnit2rVel(k) +c & )*recip_rThickC(i,j)*recip_deepFacF(k)*wUnit2rVel(k) + ENDDO + ENDDO + ENDIF + ELSE - DO j=1-Oly,sNy+Oly-1 - DO i=1-Olx,sNx+Olx-1 + DO j=1-OLy,sNy+OLy-1 + DO i=1-OLx,sNx+OLx-1 wCoriolisTerm(i,j) = 0. _d 0 ENDDO ENDDO diff --git a/pkg/mom_fluxform/mom_fluxform.F b/pkg/mom_fluxform/mom_fluxform.F index ab5f4b05ff..8a0bd51ee9 100644 --- a/pkg/mom_fluxform/mom_fluxform.F +++ b/pkg/mom_fluxform/mom_fluxform.F @@ -108,10 +108,10 @@ SUBROUTINE MOM_FLUXFORM( C !LOCAL VARIABLES: ==================================================== C i,j :: loop indices -C vF :: viscous flux -C v4F :: bi-harmonic viscous flux +C gAdd :: individual momentum tendency term to add +C del2 :: horizontal second derivative (Laplacian) C uCf,vCf :: Coriolis acceleration -C mT :: Metric terms +C vMT :: vertical/Non-Hyd. metric terms C fZon :: zonal fluxes C fMer :: meridional fluxes C fVrUp,fVrDw :: vertical viscous fluxes at interface k & k+1 @@ -120,11 +120,11 @@ SUBROUTINE MOM_FLUXFORM( C kkey :: tape key (depends on level and tile indices) INTEGER kkey #endif - _RL vF(1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL uCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL vCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy) - _RL mT(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL gAdd(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL del2(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL uCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL vCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL vMT(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL fZon(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL fMer(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL fVrUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy) @@ -171,7 +171,7 @@ SUBROUTINE MOM_FLUXFORM( _RL AhDudxFac _RL vDudyFac _RL AhDudyFac - _RL rVelDudrFac + _RL wDudrFac _RL ArDudrFac _RL fuFac _RL mtFacU @@ -180,13 +180,14 @@ SUBROUTINE MOM_FLUXFORM( _RL AhDvdxFac _RL vDvdyFac _RL AhDvdyFac - _RL rVelDvdrFac + _RL wDvdrFac _RL ArDvdrFac _RL fvFac _RL mtFacV _RL mtNHFacV _RL sideMaskFac - LOGICAL bottomDragTerms, metricTerms + _RL rAdvDeepFac + LOGICAL bottomDragTerms, metricTerms, diagFlag CEOP #ifdef MOM_BOUNDARY_CONSERVE COMMON / MOM_FLUXFORM_LOCAL / uBnd, vBnd @@ -202,11 +203,11 @@ SUBROUTINE MOM_FLUXFORM( C Initialise intermediate terms DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx - vF(i,j) = 0. - v4F(i,j) = 0. + gAdd(i,j) = 0. + del2(i,j) = 0. uCf(i,j) = 0. vCf(i,j) = 0. - mT(i,j) = 0. + vMT(i,j) = 0. fZon(i,j) = 0. fMer(i,j) = 0. fVrUp(i,j)= 0. @@ -237,7 +238,7 @@ SUBROUTINE MOM_FLUXFORM( AhDudxFac = vfFacMom*1. vDudyFac = afFacMom*1. AhDudyFac = vfFacMom*1. - rVelDudrFac = afFacMom*1. + wDudrFac = afFacMom*1. ArDudrFac = vfFacMom*1. mtFacU = mtFacMom*1. mtNHFacU = 1. @@ -247,7 +248,7 @@ SUBROUTINE MOM_FLUXFORM( AhDvdxFac = vfFacMom*1. vDvdyFac = afFacMom*1. AhDvdyFac = vfFacMom*1. - rVelDvdrFac = afFacMom*1. + wDvdrFac = afFacMom*1. ArDvdrFac = vfFacMom*1. mtFacV = mtFacMom*1. mtNHFacV = 1. @@ -407,10 +408,12 @@ SUBROUTINE MOM_FLUXFORM( I myTime, myIter, myThid ) C- Free surface correction term (flux at k=1) - CALL MOM_U_ADV_WU( bi,bj,k,uVel,wVel,rTransU, + CALL MOM_U_ADV_WU( bi, bj, k, deepFacAdv, + I uVel, wVel, rTransU, O fVerUkm, myThid ) - CALL MOM_V_ADV_WV( bi,bj,k,vVel,wVel,rTransV, + CALL MOM_V_ADV_WV( bi, bj, k, deepFacAdv, + I vVel, wVel, rTransV, O fVerVkm, myThid ) C--- endif momAdvection & k=1 @@ -474,13 +477,17 @@ SUBROUTINE MOM_FLUXFORM( IF (momAdvection) THEN C--- Calculate mean fluxes (advection) between cells for zonal flow. +C-- For Deep-Model (deepAtmosphere=T), no need for NHMTerm w*u/r since +C vertical advective flux (in MOM_U_ADV_WU) is changed to W*(u*deepFacC) +C with rAdvDeepFac factor in flux divergence + #ifdef MOM_BOUNDARY_CONSERVE CALL MOM_U_ADV_UU( bi,bj,k,uTrans,uBnd(1-OLx,1-OLy,k,bi,bj), O fZon,myThid ) CALL MOM_U_ADV_VU( bi,bj,k,vTrans,uBnd(1-OLx,1-OLy,k,bi,bj), O fMer,myThid ) - CALL MOM_U_ADV_WU( - I bi,bj,k+1,uBnd,wVel,rTransU, + CALL MOM_U_ADV_WU( bi, bj, k+1, deepFacAdv, + I uBnd, wVel, rTransU, O fVerUkp, myThid ) #else /* MOM_BOUNDARY_CONSERVE */ C-- Zonal flux (fZon is at east face of "u" cell) @@ -493,28 +500,48 @@ SUBROUTINE MOM_FLUXFORM( C-- Vertical flux (fVer is at upper face of "u" cell) C Mean flow component of vertical flux (at k+1) -> fVer - CALL MOM_U_ADV_WU( - I bi,bj,k+1,uVel,wVel,rTransU, + CALL MOM_U_ADV_WU( bi, bj, k+1, deepFacAdv, + I uVel, wVel, rTransU, O fVerUkp, myThid ) #endif /* MOM_BOUNDARY_CONSERVE */ C-- Tendency is minus divergence of the fluxes + coriolis + pressure term - DO j=jMin,jMax - DO i=iMin,iMax - gU(i,j,k,bi,bj) = + rAdvDeepFac = wDudrFac*rkSign +#ifndef MOM_USE_OLD_DEEP_VERT_ADV + IF ( useNHMTerms ) rAdvDeepFac = rAdvDeepFac*recip_deepFacC(k) +#endif + IF ( selectMetricTerms.EQ.3 ) THEN +C- In this case fMer contains the advective flux vTrans*(u*dxC) + DO j=jMin,jMax + DO i=iMin,iMax + gU(i,j,k,bi,bj) = + & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k) + & *recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k) + & *( ( fZon( i, j ) - fZon(i-1,j ) )*uDudxFac + & +( fMer( i,j+1) - fMer( i, j ) )*vDudyFac + & *recip_dxC(i,j,bi,bj) + & +( fVerUkp(i,j) - fVerUkm(i,j) )*rAdvDeepFac + & ) + ENDDO + ENDDO + ELSE + DO j=jMin,jMax + DO i=iMin,iMax + gU(i,j,k,bi,bj) = #ifdef OLD_UV_GEOM - & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/ + & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k)/ & ( 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj)) ) #else - & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k) - & *recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k) + & -_recip_hFacW(i,j,k,bi,bj)*recip_drF(k) + & *recip_rAw(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k) #endif - & *( ( fZon(i,j ) - fZon(i-1,j) )*uDudxFac - & +( fMer(i,j+1) - fMer(i, j) )*vDudyFac - & +( fVerUkp(i,j) - fVerUkm(i,j) )*rkSign*rVelDudrFac - & ) + & *( ( fZon( i, j ) - fZon(i-1,j ) )*uDudxFac + & +( fMer( i,j+1) - fMer( i, j ) )*vDudyFac + & +( fVerUkp(i,j) - fVerUkm(i,j) )*rAdvDeepFac + & ) + ENDDO ENDDO - ENDDO + ENDIF #ifdef ALLOW_DIAGNOSTICS IF ( useDiagnostics ) THEN @@ -575,20 +602,20 @@ SUBROUTINE MOM_FLUXFORM( IF (momViscosity) THEN C--- Calculate eddy fluxes (dissipation) between cells for zonal flow. -C Bi-harmonic term del^2 U -> v4F +C Bi-harmonic term del^2 U -> del2 IF ( useBiharmonicVisc ) & CALL MOM_U_DEL2U( bi, bj, k, uFld, hFacZ, h0FacZ, - O v4f, myThid ) + O del2, myThid ) #ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE v4f(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte +CADJ STORE del2(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte #endif C Laplacian and bi-harmonic terms, Zonal Fluxes -> fZon - CALL MOM_U_XVISCFLUX( bi,bj,k,uFld,v4F,fZon, + CALL MOM_U_XVISCFLUX( bi,bj,k,uFld,del2,fZon, I viscAh_D,viscA4_D,myThid ) C Laplacian and bi-harmonic termis, Merid Fluxes -> fMer - CALL MOM_U_YVISCFLUX( bi,bj,k,uFld,v4F,hFacZ,fMer, + CALL MOM_U_YVISCFLUX( bi,bj,k,uFld,del2,hFacZ,fMer, I viscAh_Z,viscA4_Z,myThid ) C Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw @@ -630,14 +657,14 @@ SUBROUTINE MOM_FLUXFORM( IF (no_slip_sides) THEN C- No-slip BCs impose a drag at walls... CALL MOM_U_SIDEDRAG( bi, bj, k, - I uFld, v4f, h0FacZ, + I uFld, del2, h0FacZ, I viscAh_Z, viscA4_Z, I useHarmonicVisc, useBiharmonicVisc, useVariableVisc, - O vF, + O gAdd, I myThid ) DO j=jMin,jMax DO i=iMin,iMax - guDiss(i,j) = guDiss(i,j) + vF(i,j) + guDiss(i,j) = guDiss(i,j) + gAdd(i,j) ENDDO ENDDO ENDIF @@ -702,32 +729,58 @@ SUBROUTINE MOM_FLUXFORM( c I myTime,myThid) C-- Metric terms for curvilinear grid systems - IF (useNHMTerms) THEN -C o Non-Hydrostatic (spherical) metric terms - CALL MOM_U_METRIC_NH( bi,bj,k,uFld,wVel,mT,myThid ) + diagFlag = .FALSE. +#ifdef MOM_USE_OLD_DEEP_VERT_ADV + IF ( useNHMTerms ) THEN +#else + IF ( useNHMTerms .AND. .NOT.deepAtmosphere ) THEN +#endif +C o Non-Hydrostatic/vertical (spherical) metric terms + CALL MOM_U_METRIC_NH( bi,bj,k,uFld,wVel,vMT,myThid ) + diagFlag = .TRUE. DO j=jMin,jMax DO i=iMin,iMax - gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtNHFacU*mT(i,j) + gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtNHFacU*vMT(i,j) ENDDO ENDDO +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics .AND. + & ( selectMetricTerms.EQ.0 .OR. selectMetricTerms.EQ.3 ) ) THEN + CALL DIAGNOSTICS_FILL( vMT, 'Um_Metr ', k,1,2, bi,bj, myThid ) + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ ENDIF - IF ( usingSphericalPolarGrid .AND. metricTerms ) THEN + IF ( selectMetricTerms.GE.1 .AND. selectMetricTerms.NE.3 ) THEN + IF ( usingSphericalPolarGrid ) THEN C o Spherical polar grid metric terms - CALL MOM_U_METRIC_SPHERE( bi,bj,k,uFld,vFld,mT,myThid ) - DO j=jMin,jMax - DO i=iMin,iMax - gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*mT(i,j) + CALL MOM_U_METRIC_SPHERE( bi,bj,k,uFld,vFld,gAdd,myThid ) + DO j=jMin,jMax + DO i=iMin,iMax + gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*gAdd(i,j) + ENDDO ENDDO - ENDDO - ENDIF - IF ( usingCylindricalGrid .AND. metricTerms ) THEN + ENDIF + IF ( usingCylindricalGrid ) THEN C o Cylindrical grid metric terms - CALL MOM_U_METRIC_CYLINDER( bi,bj,k,uFld,vFld,mT,myThid ) - DO j=jMin,jMax - DO i=iMin,iMax - gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*mT(i,j) + CALL MOM_U_METRIC_CYLINDER( bi,bj,k,uFld,vFld,gAdd,myThid ) + DO j=jMin,jMax + DO i=iMin,iMax + gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mtFacU*gAdd(i,j) + ENDDO ENDDO - ENDDO + ENDIF +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN + IF ( diagFlag ) THEN + DO j=jMin,jMax + DO i=iMin,iMax + gAdd(i,j) = gAdd(i,j) + vMT(i,j) + ENDDO + ENDDO + ENDIF + CALL DIAGNOSTICS_FILL( gAdd, 'Um_Metr ', k,1,2, bi,bj, myThid ) + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| @@ -736,12 +789,17 @@ SUBROUTINE MOM_FLUXFORM( IF (momAdvection) THEN +C-- For Deep-Model (deepAtmosphere=T), no need for NHMTerm w*v/r since +C vertical advective flux (in MOM_V_ADV_WV) is changed to W*(v*deepFacC) +C with rAdvDeepFac factor in flux divergence + #ifdef MOM_BOUNDARY_CONSERVE CALL MOM_V_ADV_UV( bi,bj,k,uTrans,vBnd(1-OLx,1-OLy,k,bi,bj), O fZon,myThid ) CALL MOM_V_ADV_VV( bi,bj,k,vTrans,vBnd(1-OLx,1-OLy,k,bi,bj), O fMer,myThid ) - CALL MOM_V_ADV_WV( bi,bj,k+1,vBnd,wVel,rTransV, + CALL MOM_V_ADV_WV( bi, bj, k+1, deepFacAdv, + I vBnd, wVel, rTransV, O fVerVkp, myThid ) #else /* MOM_BOUNDARY_CONSERVE */ C--- Calculate mean fluxes (advection) between cells for meridional flow. @@ -754,11 +812,16 @@ SUBROUTINE MOM_FLUXFORM( C-- Vertical flux (fVer is at upper face of "v" cell) C Mean flow component of vertical flux (at k+1) -> fVerV - CALL MOM_V_ADV_WV( bi,bj,k+1,vVel,wVel,rTransV, + CALL MOM_V_ADV_WV( bi, bj, k+1, deepFacAdv, + I vVel, wVel, rTransV, O fVerVkp, myThid ) #endif /* MOM_BOUNDARY_CONSERVE */ C-- Tendency is minus divergence of the fluxes + coriolis + pressure term + rAdvDeepFac = wDvdrFac*rkSign +#ifndef MOM_USE_OLD_DEEP_VERT_ADV + IF ( useNHMTerms ) rAdvDeepFac = rAdvDeepFac*recip_deepFacC(k) +#endif DO j=jMin,jMax DO i=iMin,iMax gV(i,j,k,bi,bj) = @@ -771,7 +834,7 @@ SUBROUTINE MOM_FLUXFORM( #endif & *( ( fZon(i+1,j) - fZon(i,j ) )*uDvdxFac & +( fMer(i, j) - fMer(i,j-1) )*vDvdyFac - & +( fVerVkp(i,j) - fVerVkm(i,j) )*rkSign*rVelDvdrFac + & +( fVerVkp(i,j) - fVerVkm(i,j) )*rAdvDeepFac & ) ENDDO ENDDO @@ -834,20 +897,20 @@ SUBROUTINE MOM_FLUXFORM( IF (momViscosity) THEN C--- Calculate eddy fluxes (dissipation) between cells for meridional flow. -C Bi-harmonic term del^2 V -> v4F +C Bi-harmonic term del^2 V -> del2 IF ( useBiharmonicVisc ) & CALL MOM_V_DEL2V( bi, bj, k, vFld, hFacZ, h0FacZ, - O v4f, myThid ) + O del2, myThid ) #ifdef ALLOW_AUTODIFF_TAMC -CADJ STORE v4f(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte +CADJ STORE del2(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte #endif C Laplacian and bi-harmonic terms, Zonal Fluxes -> fZon - CALL MOM_V_XVISCFLUX( bi,bj,k,vFld,v4f,hFacZ,fZon, + CALL MOM_V_XVISCFLUX( bi,bj,k,vFld,del2,hFacZ,fZon, I viscAh_Z,viscA4_Z,myThid ) C Laplacian and bi-harmonic termis, Merid Fluxes -> fMer - CALL MOM_V_YVISCFLUX( bi,bj,k,vFld,v4f,fMer, + CALL MOM_V_YVISCFLUX( bi,bj,k,vFld,del2,fMer, I viscAh_D,viscA4_D,myThid ) C Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw @@ -889,14 +952,14 @@ SUBROUTINE MOM_FLUXFORM( IF (no_slip_sides) THEN C- No-slip BCs impose a drag at walls... CALL MOM_V_SIDEDRAG( bi, bj, k, - I vFld, v4f, h0FacZ, + I vFld, del2, h0FacZ, I viscAh_Z, viscA4_Z, I useHarmonicVisc, useBiharmonicVisc, useVariableVisc, - O vF, + O gAdd, I myThid ) DO j=jMin,jMax DO i=iMin,iMax - gvDiss(i,j) = gvDiss(i,j) + vF(i,j) + gvDiss(i,j) = gvDiss(i,j) + gAdd(i,j) ENDDO ENDDO ENDIF @@ -961,32 +1024,57 @@ SUBROUTINE MOM_FLUXFORM( c I myTime,myThid) C-- Metric terms for curvilinear grid systems - IF (useNHMTerms) THEN -C o Non-Hydrostatic (spherical) metric terms - CALL MOM_V_METRIC_NH( bi,bj,k,vFld,wVel,mT,myThid ) + diagFlag = .FALSE. +#ifdef MOM_USE_OLD_DEEP_VERT_ADV + IF ( useNHMTerms ) THEN +#else + IF ( useNHMTerms .AND. .NOT.deepAtmosphere ) THEN +#endif +C o Non-Hydrostatic/vertical (spherical) metric terms + CALL MOM_V_METRIC_NH( bi,bj,k,vFld,wVel,vMT,myThid ) + diagFlag = .TRUE. DO j=jMin,jMax DO i=iMin,iMax - gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtNHFacV*mT(i,j) + gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtNHFacV*vMT(i,j) ENDDO ENDDO +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics .AND. selectMetricTerms.EQ.0 ) THEN + CALL DIAGNOSTICS_FILL( vMT, 'Vm_Metr ', k,1,2, bi,bj, myThid ) + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ ENDIF - IF ( usingSphericalPolarGrid .AND. metricTerms ) THEN + IF ( selectMetricTerms.GE.1 ) THEN + IF ( usingSphericalPolarGrid ) THEN C o Spherical polar grid metric terms - CALL MOM_V_METRIC_SPHERE( bi,bj,k,uFld,mT,myThid ) - DO j=jMin,jMax - DO i=iMin,iMax - gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*mT(i,j) + CALL MOM_V_METRIC_SPHERE( bi,bj,k,uFld,gAdd,myThid ) + DO j=jMin,jMax + DO i=iMin,iMax + gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*gAdd(i,j) + ENDDO ENDDO - ENDDO - ENDIF - IF ( usingCylindricalGrid .AND. metricTerms ) THEN + ENDIF + IF ( usingCylindricalGrid ) THEN C o Cylindrical grid metric terms - CALL MOM_V_METRIC_CYLINDER( bi,bj,k,uFld,vFld,mT,myThid ) - DO j=jMin,jMax - DO i=iMin,iMax - gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*mT(i,j) + CALL MOM_V_METRIC_CYLINDER( bi,bj,k,uFld,vFld,gAdd,myThid ) + DO j=jMin,jMax + DO i=iMin,iMax + gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mtFacV*gAdd(i,j) + ENDDO ENDDO - ENDDO + ENDIF +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN + IF ( diagFlag ) THEN + DO j=jMin,jMax + DO i=iMin,iMax + gAdd(i,j) = gAdd(i,j) + vMT(i,j) + ENDDO + ENDDO + ENDIF + CALL DIAGNOSTICS_FILL( gAdd, 'Vm_Metr ', k,1,2, bi,bj, myThid ) + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| diff --git a/pkg/mom_fluxform/mom_u_adv_vu.F b/pkg/mom_fluxform/mom_u_adv_vu.F index 75b13ad3dc..62d85837a2 100644 --- a/pkg/mom_fluxform/mom_u_adv_vu.F +++ b/pkg/mom_fluxform/mom_u_adv_vu.F @@ -5,10 +5,10 @@ C !INTERFACE: ========================================================== SUBROUTINE MOM_U_ADV_VU( - I bi,bj,k, + I bi, bj, k, I vTrans, uFld, O AdvectFluxVU, - I myThid) + I myThid ) C !DESCRIPTION: C Calculates the meridional advective flux of zonal momentum: @@ -24,12 +24,12 @@ SUBROUTINE MOM_U_ADV_VU( #include "GRID.h" C !INPUT PARAMETERS: =================================================== -C bi,bj :: tile indices +C bi, bj :: tile indices C k :: vertical level C vTrans :: meridional transport -C uFld :: zonal flow -C myThid :: thread number - INTEGER bi,bj,k +C uFld :: zonal velocity +C myThid :: my Thread Id number + INTEGER bi, bj, k _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER myThid @@ -43,22 +43,39 @@ SUBROUTINE MOM_U_ADV_VU( INTEGER i,j CEOP - DO j=1-Oly+1,sNy+Oly - DO i=1-Olx+1,sNx+Olx - AdvectFluxVU(i,j) = - & 0.25*( vTrans(i,j) + vTrans(i-1,j) ) + IF ( selectMetricTerms.NE.3 ) THEN + DO j=1-OLy+1,sNy+OLy + DO i=1-OLx+1,sNx+OLx + AdvectFluxVU(i,j) = 0.25 _d 0 + & *( vTrans(i,j) + vTrans(i-1,j) ) #ifdef MOM_BOUNDARY_CONSERVE - & *( uFld(i,j)*_maskW(i,j-1,k,bi,bj) + & *( uFld(i,j)*_maskW(i,j-1,k,bi,bj) & + uFld(i,j-1)*_maskW(i,j,k,bi,bj) ) #else - & *( uFld(i,j) + uFld(i,j-1) ) + & *( uFld(i,j) + uFld(i,j-1) ) #endif #ifdef OLD_ADV_BCS - & *_maskW(i,j,k,bi,bj) - & *_maskW(i,j-1,k,bi,bj) + & *_maskW(i,j,k,bi,bj) + & *_maskW(i,j-1,k,bi,bj) #endif /* OLD_ADV_BCS */ - ENDDO - ENDDO + ENDDO + ENDDO + ELSE +C- Advect u*dxC (--> account for metric term u*v*tanPhi/R) + DO j=1-OLy+1,sNy+OLy + DO i=1-OLx+1,sNx+OLx + AdvectFluxVU(i,j) = 0.25 _d 0 + & *( vTrans(i,j) + vTrans(i-1,j) ) +#ifdef MOM_BOUNDARY_CONSERVE + & *( uFld(i, j )*dxC(i, j, bi,bj)*_maskW(i,j-1,k,bi,bj) + & + uFld(i,j-1)*dxC(i,j-1,bi,bj)*_maskW(i, j, k,bi,bj) ) +#else + & *( uFld(i, j )*dxC(i, j, bi,bj) + & + uFld(i,j-1)*dxC(i,j-1,bi,bj) ) +#endif + ENDDO + ENDDO + ENDIF RETURN END diff --git a/pkg/mom_fluxform/mom_u_adv_wu.F b/pkg/mom_fluxform/mom_u_adv_wu.F index ed604a1643..6fa4ad44a3 100644 --- a/pkg/mom_fluxform/mom_u_adv_wu.F +++ b/pkg/mom_fluxform/mom_u_adv_wu.F @@ -5,8 +5,8 @@ C !INTERFACE: ========================================================== SUBROUTINE MOM_U_ADV_WU( - I bi,bj,k, - I uFld,wFld,rTrans, + I bi, bj, k, deepFacA, + I uFld, wFld, rTrans, O advectiveFluxWU, I myThid ) @@ -24,13 +24,15 @@ SUBROUTINE MOM_U_ADV_WU( #include "GRID.h" C !INPUT PARAMETERS: =================================================== -C bi,bj :: tile indices +C bi, bj :: tile indices C k :: vertical level +C deepFacA :: deep-model grid factor at level center C uFld :: zonal velocity C wFld :: vertical velocity C rTrans :: vertical transport (above U point) -C myThid :: thread number - INTEGER bi,bj,k +C myThid :: my Thread Id number + INTEGER bi, bj, k + _RL deepFacA(Nr) _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL wFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy) @@ -41,69 +43,73 @@ SUBROUTINE MOM_U_ADV_WU( _RL advectiveFluxWU(1-OLx:sNx+OLx,1-OLy:sNy+OLy) C !LOCAL VARIABLES: ==================================================== -C i,j :: loop indices - INTEGER i,j +C i, j :: loop indices + INTEGER i, j CEOP IF ( k.EQ.Nr+1 .AND. & useRealFreshWaterFlux .AND. usingPCoords ) THEN - DO j=1-Oly+1,sNy+Oly - DO i=1-Olx+1,sNx+Olx + DO j=1-OLy+1,sNy+OLy + DO i=1-OLx+1,sNx+OLx advectiveFluxWU(i,j) = rTrans(i,j)*uFld(i,j,k-1,bi,bj) + & *deepFacA(k-1) ENDDO ENDDO ELSEIF ( k.GT.Nr .OR. (k.EQ.1.AND.rigidLid) ) THEN C Advective flux = 0 at k=Nr+1 ; = 0 at k=1 if rigid-lid - DO j=1-Oly,sNy+Oly - DO i=1-Olx,sNx+Olx - advectiveFluxWU(i,j) = 0. + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + advectiveFluxWU(i,j) = 0. + ENDDO ENDDO - ENDDO ELSEIF (k.EQ.1) THEN C (linear) Free-surface correction at k=1 - DO j=1-Oly+1,sNy+Oly - DO i=1-Olx+1,sNx+Olx - advectiveFluxWU(i,j) = rTrans(i,j)*uFld(i,j,k,bi,bj) + DO j=1-OLy+1,sNy+OLy + DO i=1-OLx+1,sNx+OLx + advectiveFluxWU(i,j) = rTrans(i,j)*uFld(i,j,k,bi,bj) + & *deepFacA(k) + ENDDO ENDDO - ENDDO ELSE C Vertical advection - interior ; assume uFld & wFld are masked - DO j=1-Oly+1,sNy+Oly - DO i=1-Olx+1,sNx+Olx - advectiveFluxWU(i,j) = rTrans(i,j)* + DO j=1-OLy+1,sNy+OLy + DO i=1-OLx+1,sNx+OLx + advectiveFluxWU(i,j) = rTrans(i,j)*0.5 _d 0 #ifdef MOM_BOUNDARY_CONSERVE - & 0.5 _d 0*( uFld(i,j,k,bi,bj)*_maskW(i,j,k-1,bi,bj) - & +uFld(i,j,k-1,bi,bj)*_maskW(i,j,k,bi,bj) ) + & *( uFld(i,j, k ,bi,bj)*deepFacA( k )*_maskW(i,j,k-1,bi,bj) + & + uFld(i,j,k-1,bi,bj)*deepFacA(k-1)*_maskW(i,j, k ,bi,bj) ) #else - & 0.5 _d 0*( uFld(i,j,k,bi,bj)+uFld(i,j,k-1,bi,bj) ) + & *( uFld(i,j, k ,bi,bj)*deepFacA( k ) + & + uFld(i,j,k-1,bi,bj)*deepFacA(k-1) ) #endif + ENDDO ENDDO - ENDDO - IF ( select_rStar.EQ.0 .AND. .NOT.rigidLid ) THEN -c & .AND. usingPCoords ) THEN + IF ( select_rStar.EQ.0 .AND. .NOT.rigidLid ) THEN +c & .AND. usingPCoords ) THEN C (linear) Free-surface correction at k>1 - DO j=1-Oly+1,sNy+Oly - DO i=1-Olx+1,sNx+Olx + DO j=1-OLy+1,sNy+OLy + DO i=1-OLx+1,sNx+OLx advectiveFluxWU(i,j) = advectiveFluxWU(i,j) - & +0.25*( - & wFld(i, j ,k,bi,bj)*rA(i, j ,bi,bj)* - & (maskC( i ,j,k,bi,bj)-maskC( i ,j,k-1,bi,bj)) - & +wFld(i-1,j,k,bi,bj)*rA(i-1,j,bi,bj)* - & (maskC(i-1,j,k,bi,bj)-maskC(i-1,j,k-1,bi,bj)) - & )*deepFac2F(k)*rhoFacF(k) - & *uFld(i,j,k,bi,bj) + & +0.25 _d 0*( + & wFld( i ,j,k,bi,bj)*rA( i ,j,bi,bj) + & *(maskC( i ,j,k,bi,bj)-maskC( i ,j,k-1,bi,bj)) + & +wFld(i-1,j,k,bi,bj)*rA(i-1,j,bi,bj) + & *(maskC(i-1,j,k,bi,bj)-maskC(i-1,j,k-1,bi,bj)) + & )*deepFac2F(k)*rhoFacF(k) + & *uFld(i,j,k,bi,bj)*deepFacA(k) ENDDO ENDDO -C- endif NOT rigidLid - ENDIF +C- endif select_rStar=0 and NOT rigidLid + ENDIF +C- endif special k cases ENDIF RETURN diff --git a/pkg/mom_fluxform/mom_u_coriolis.F b/pkg/mom_fluxform/mom_u_coriolis.F index 9fc1e030ee..588d1d6a33 100644 --- a/pkg/mom_fluxform/mom_u_coriolis.F +++ b/pkg/mom_fluxform/mom_u_coriolis.F @@ -5,9 +5,9 @@ C !INTERFACE: ========================================================== SUBROUTINE MOM_U_CORIOLIS( - I bi,bj,k,vFld, + I bi, bj, k, vFld, U uCoriolisTerm, - I myThid) + I myThid ) C !DESCRIPTION: C Calculates the horizontal Coriolis term in the zonal equation: @@ -24,11 +24,11 @@ SUBROUTINE MOM_U_CORIOLIS( #include "SURFACE.h" C !INPUT PARAMETERS: =================================================== -C bi,bj :: tile indices +C bi, bj :: tile indices C k :: vertical level C vFld :: meridional flow C myThid :: thread number - INTEGER bi,bj,k + INTEGER bi, bj, k _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER myThid @@ -38,45 +38,54 @@ SUBROUTINE MOM_U_CORIOLIS( C !LOCAL VARIABLES: ==================================================== C i,j :: loop indices - INTEGER i,j - _RS one - PARAMETER( one = 1. _d 0 ) + INTEGER i, j CEOP - IF ( selectCoriScheme.GE.2 ) THEN -C Energy conserving discretization + IF ( selectCoriScheme.LE.1 ) THEN +C- Original discretization DO j=1-OLy,sNy+OLy-1 DO i=1-OLx+1,sNx+OLx - uCoriolisTerm(i,j) = - & 0.5*( _fCori( i ,j,bi,bj) + uCoriolisTerm(i,j) = 0.5 _d 0 + & *( _fCori( i ,j,bi,bj) + _fCori(i-1,j,bi,bj) ) + & *0.25*( vFld( i ,j)+vFld( i ,j+1) + & +vFld(i-1,j)+vFld(i-1,j+1) + & ) + ENDDO + ENDDO + ELSEIF ( selectCoriScheme.LE.3 ) THEN +C- Energy conserving discretization + DO j=1-OLy,sNy+OLy-1 + DO i=1-OLx+1,sNx+OLx + uCoriolisTerm(i,j) = 0.5 _d 0 + & *( _fCori( i ,j,bi,bj) & *0.5*( vFld( i ,j)+vFld( i ,j+1) ) & + _fCori(i-1,j,bi,bj) & *0.5*( vFld(i-1,j)+vFld(i-1,j+1) ) ) ENDDO ENDDO ELSE -C Original discretization +C- Using averaged transport: DO j=1-OLy,sNy+OLy-1 DO i=1-OLx+1,sNx+OLx - uCoriolisTerm(i,j) = - & 0.5*( _fCori( i ,j,bi,bj) + - & _fCori(i-1,j,bi,bj) ) - & *0.25*( - & vFld( i ,j)+vFld( i ,j+1) - & +vFld(i-1,j)+vFld(i-1,j+1) - & ) + uCoriolisTerm(i,j) = 0.5 _d 0 + & *( _fCori( i ,j,bi,bj) + _fCori(i-1,j,bi,bj) ) + & *( vFld( i , j )*dxG( i , j ,bi,bj)*hFacS( i , j ,k,bi,bj) + & + vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*hFacS( i ,j+1,k,bi,bj) + & + vFld(i-1, j )*dxG(i-1, j ,bi,bj)*hFacS(i-1, j ,k,bi,bj) + & + vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*hFacS(i-1,j+1,k,bi,bj) + & )*0.25 _d 0*recip_dxC(i,j,bi,bj)*recip_hFacW(i,j,k,bi,bj) ENDDO ENDDO ENDIF IF ( selectCoriScheme.EQ.1 .OR. selectCoriScheme.EQ.3 ) THEN -C Scale term so that only "wet" points are used -C Due to: Jamart and Ozer, 1986, JGR 91 (C9), 10,621-10,631 +C- Scale term so that only "wet" points are used +C Due to: Jamart and Ozer, 1986, JGR 91 (C9), 10,621-10,631 C "Numerical Boundary Layers and Spurious Residual Flows" DO j=1-OLy,sNy+OLy-1 DO i=1-OLx+1,sNx+OLx uCoriolisTerm(i,j) = uCoriolisTerm(i,j) - & *4. _d 0/MAX( one, + & *4. _d 0/MAX( oneRS, & maskS( i , j ,k,bi,bj)+maskS( i ,j+1,k,bi,bj) & +maskS(i-1, j ,k,bi,bj)+maskS(i-1,j+1,k,bi,bj) ) ENDDO diff --git a/pkg/mom_fluxform/mom_u_metric_sphere.F b/pkg/mom_fluxform/mom_u_metric_sphere.F index b3a4e41cf3..e173a67cfa 100644 --- a/pkg/mom_fluxform/mom_u_metric_sphere.F +++ b/pkg/mom_fluxform/mom_u_metric_sphere.F @@ -5,7 +5,7 @@ C !INTERFACE: ========================================================== SUBROUTINE MOM_U_METRIC_SPHERE( - I bi,bj,k, + I bi, bj, k, I uFld, vFld, O uMetricTerms, I myThid ) @@ -24,12 +24,12 @@ SUBROUTINE MOM_U_METRIC_SPHERE( #include "GRID.h" C !INPUT PARAMETERS: =================================================== -C bi,bj :: tile indices +C bi, bj :: tile indices C k :: vertical level -C uFld :: zonal flow -C vFld :: meridional flow -C myThid :: thread number - INTEGER bi,bj,k +C uFld :: zonal velocity +C vFld :: meridional velocity +C myThid :: my Thread Id number + INTEGER bi, bj, k _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER myThid @@ -39,19 +39,41 @@ SUBROUTINE MOM_U_METRIC_SPHERE( _RL uMetricTerms(1-OLx:sNx+OLx,1-OLy:sNy+OLy) C !LOCAL VARIABLES: ==================================================== -C i,j :: loop indices - INTEGER i,j +C i, j :: loop indices + INTEGER i, j CEOP - DO j=1-Olx,sNy+Oly-1 - DO i=1-Olx+1,sNx+Olx - uMetricTerms(i,j) = - & uFld(i,j)*recip_rSphere*recip_deepFacC(k) - & *0.25*( vFld(i,j )+vFld(i-1,j ) - & +vFld(i,j+1)+vFld(i-1,j+1) - & )*_tanPhiAtU(i,j,bi,bj) - ENDDO - ENDDO + IF ( selectMetricTerms.EQ.1 ) THEN +C- Using analytical expression for tan(Phi) (stored in: tanPhiAtU) + DO j=1-OLx,sNy+OLy-1 + DO i=1-OLx+1,sNx+OLx + uMetricTerms(i,j) = + & uFld(i,j)*recip_rSphere*recip_deepFacC(k) + & *0.25*( vFld(i,j )+vFld(i-1,j ) + & +vFld(i,j+1)+vFld(i-1,j+1) + & )*_tanPhiAtU(i,j,bi,bj) + ENDDO + ENDDO + ELSE +C- Using grid-spacing gradient: estimates tan(Phi)/rSphere as -del^j(dxC)/rAz + DO j=2-OLx,sNy+OLy-1 + DO i=2-OLx,sNx+OLx + uMetricTerms(i,j) = -0.5 _d 0* + & ( ( uFld(i,j-1) + uFld(i,j) )*halfRL + & *( dxC(i,j,bi,bj)-dxC(i,j-1,bi,bj) )*recip_rAz(i, j, bi,bj) + & *( vFld(i-1, j )*dxG(i-1, j, bi,bj)*hFacS(i-1, j, k,bi,bj) + & + vFld( i, j )*dxG( i, j, bi,bj)*hFacS( i, j, k,bi,bj) + & )*halfRL + & + ( uFld(i,j) + uFld(i,j+1) )*halfRL + & *( dxC(i,j+1,bi,bj)-dxC(i,j,bi,bj) )*recip_rAz(i,j+1,bi,bj) + & *( vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*hFacS(i-1,j+1,k,bi,bj) + & + vFld( i, j+1)*dxG( i, j+1,bi,bj)*hFacS( i, j+1,k,bi,bj) + & )*halfRL + & )*recip_hFacW(i,j,k,bi,bj) + & *recip_dxC(i,j,bi,bj)*recip_deepFacC(k) + ENDDO + ENDDO + ENDIF RETURN END diff --git a/pkg/mom_fluxform/mom_v_adv_wv.F b/pkg/mom_fluxform/mom_v_adv_wv.F index 18290a5010..4802e71bc3 100644 --- a/pkg/mom_fluxform/mom_v_adv_wv.F +++ b/pkg/mom_fluxform/mom_v_adv_wv.F @@ -1,15 +1,14 @@ #include "MOM_FLUXFORM_OPTIONS.h" -#undef OLD_MOM_ADV_W CBOP C !ROUTINE: MOM_V_ADV_WV C !INTERFACE: ========================================================== SUBROUTINE MOM_V_ADV_WV( - I bi,bj,k, - I vFld,wFld,rTrans, + I bi, bj, k, deepFacA, + I vFld, wFld, rTrans, O advectiveFluxWV, - I myThid) + I myThid ) C !DESCRIPTION: C Calculates the vertical advective flux of meridional momentum: @@ -25,13 +24,15 @@ SUBROUTINE MOM_V_ADV_WV( #include "GRID.h" C !INPUT PARAMETERS: =================================================== -C bi,bj :: tile indices +C bi, bj :: tile indices C k :: vertical level +C deepFacA :: deep-model grid factor at level center C vFld :: meridional velocity C wFld :: vertical velocity C rTrans :: vertical transport (above V point) -C myThid :: thread number - INTEGER bi,bj,k +C myThid :: my Thread Id number + INTEGER bi, bj, k + _RL deepFacA(Nr) _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL wFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy) @@ -42,69 +43,73 @@ SUBROUTINE MOM_V_ADV_WV( _RL advectiveFluxWV(1-OLx:sNx+OLx,1-OLy:sNy+OLy) C !LOCAL VARIABLES: ==================================================== -C i,j :: loop indices - INTEGER i,j +C i, j :: loop indices + INTEGER i, j CEOP IF ( k.EQ.Nr+1 .AND. & useRealFreshWaterFlux .AND. usingPCoords ) THEN - DO j=1-Oly+1,sNy+Oly - DO i=1-Olx+1,sNx+Olx + DO j=1-OLy+1,sNy+OLy + DO i=1-OLx+1,sNx+OLx advectiveFluxWV(i,j) = rTrans(i,j)*vFld(i,j,k-1,bi,bj) + & *deepFacA(k-1) ENDDO ENDDO ELSEIF ( k.GT.Nr .OR. (k.EQ.1.AND.rigidLid) ) THEN C Advective flux = 0 at k=Nr+1 ; = 0 at k=1 if rigid-lid - DO j=1-Oly,sNy+Oly - DO i=1-Olx,sNx+Olx - advectiveFluxWV(i,j) = 0. + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + advectiveFluxWV(i,j) = 0. + ENDDO ENDDO - ENDDO ELSEIF (k.EQ.1) THEN C (linear) Free-surface correction at k=1 - DO j=1-Oly+1,sNy+Oly - DO i=1-Olx+1,sNx+Olx - advectiveFluxWV(i,j) = rTrans(i,j)*vFld(i,j,k,bi,bj) + DO j=1-OLy+1,sNy+OLy + DO i=1-OLx+1,sNx+OLx + advectiveFluxWV(i,j) = rTrans(i,j)*vFld(i,j,k,bi,bj) + & *deepFacA(k) + ENDDO ENDDO - ENDDO ELSE C Vertical advection - interior ; assume vFld & wFld are masked - DO j=1-Oly+1,sNy+Oly - DO i=1-Olx+1,sNx+Olx - advectiveFluxWV(i,j) = rTrans(i,j)* + DO j=1-OLy+1,sNy+OLy + DO i=1-OLx+1,sNx+OLx + advectiveFluxWV(i,j) = rTrans(i,j)*0.5 _d 0 #ifdef MOM_BOUNDARY_CONSERVE - & 0.5 _d 0*( vFld(i,j,k,bi,bj)*_maskS(i,j,k-1,bi,bj) - & +vFld(i,j,k-1,bi,bj)*_maskS(i,j,k,bi,bj) ) + & *( vFld(i,j, k ,bi,bj)*deepFacA( k )*_maskS(i,j,k-1,bi,bj) + & + vFld(i,j,k-1,bi,bj)*deepFacA(k-1)*_maskS(i,j, k ,bi,bj) ) #else - & 0.5 _d 0*( vFld(i,j,k,bi,bj)+vFld(i,j,k-1,bi,bj) ) + & *( vFld(i,j, k ,bi,bj)*deepFacA( k ) + & + vFld(i,j,k-1,bi,bj)*deepFacA(k-1) ) #endif + ENDDO ENDDO - ENDDO - IF ( select_rStar.EQ.0 .AND. .NOT.rigidLid ) THEN -c & .AND. usingPCoords ) THEN + IF ( select_rStar.EQ.0 .AND. .NOT.rigidLid ) THEN +c & .AND. usingPCoords ) THEN C (linear) Free-surface correction at k>1 - DO j=1-Oly+1,sNy+Oly - DO i=1-Olx+1,sNx+Olx + DO j=1-OLy+1,sNy+OLy + DO i=1-OLx+1,sNx+OLx advectiveFluxWV(i,j) = advectiveFluxWV(i,j) - & +0.25*( - & wFld(i, j ,k,bi,bj)*rA(i, j ,bi,bj)* - & (maskC(i, j ,k,bi,bj)-maskC(i, j ,k-1,bi,bj)) - & +wFld(i,j-1,k,bi,bj)*rA(i,j-1,bi,bj)* - & (maskC(i,j-1,k,bi,bj)-maskC(i,j-1,k-1,bi,bj)) - & )*deepFac2F(k)*rhoFacF(k) - & *vFld(i,j,k,bi,bj) + & +0.25 _d 0*( + & wFld(i, j ,k,bi,bj)*rA(i, j ,bi,bj) + & *(maskC(i, j ,k,bi,bj)-maskC(i, j ,k-1,bi,bj)) + & +wFld(i,j-1,k,bi,bj)*rA(i,j-1,bi,bj) + & *(maskC(i,j-1,k,bi,bj)-maskC(i,j-1,k-1,bi,bj)) + & )*deepFac2F(k)*rhoFacF(k) + & *vFld(i,j,k,bi,bj)*deepFacA(k) ENDDO ENDDO -C- endif NOT rigidLid - ENDIF +C- endif select_rStar=0 and NOT rigidLid + ENDIF +C- endif special k cases ENDIF RETURN diff --git a/pkg/mom_fluxform/mom_v_coriolis.F b/pkg/mom_fluxform/mom_v_coriolis.F index 4dd4892cda..8fc036c5ac 100644 --- a/pkg/mom_fluxform/mom_v_coriolis.F +++ b/pkg/mom_fluxform/mom_v_coriolis.F @@ -5,9 +5,9 @@ C !INTERFACE: ========================================================== SUBROUTINE MOM_V_CORIOLIS( - I bi,bj,k,uFld, + I bi, bj, k, uFld, U vCoriolisTerm, - I myThid) + I myThid ) C !DESCRIPTION: C Calculates the horizontal Coriolis term in the meridional equation: @@ -24,11 +24,11 @@ SUBROUTINE MOM_V_CORIOLIS( #include "SURFACE.h" C !INPUT PARAMETERS: =================================================== -C bi,bj :: tile indices +C bi, bj :: tile indices C k :: vertical level C uFld :: zonal flow C myThid :: thread number - INTEGER bi,bj,k + INTEGER bi, bj, k _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER myThid @@ -38,43 +38,54 @@ SUBROUTINE MOM_V_CORIOLIS( C !LOCAL VARIABLES: ==================================================== C i,j :: loop indices - INTEGER i,j - _RS one - PARAMETER( one = 1. _d 0 ) + INTEGER i, j CEOP - IF ( selectCoriScheme.GE.2 ) THEN -C Energy conserving discretization + IF ( selectCoriScheme.LE.1 ) THEN +C- Original discretization DO j=1-OLy+1,sNy+OLy DO i=1-OLx,sNx+OLx-1 - vCoriolisTerm(i,j) = - & -0.5*( _fCori(i, j ,bi,bj) + vCoriolisTerm(i,j) = -0.5 _d 0 + & *( _fCori(i, j ,bi,bj) + _fCori(i,j-1,bi,bj) ) + & *0.25*( uFld(i, j )+uFld(i+1, j ) + & +uFld(i,j-1)+uFld(i+1,j-1) + & ) + ENDDO + ENDDO + ELSEIF ( selectCoriScheme.LE.3 ) THEN +C- Energy conserving discretization + DO j=1-OLy+1,sNy+OLy + DO i=1-OLx,sNx+OLx-1 + vCoriolisTerm(i,j) = -0.5 _d 0 + & *( _fCori(i, j ,bi,bj) & *0.5*( uFld( i , j )+uFld(i+1, j ) ) & + _fCori(i,j-1,bi,bj) & *0.5*( uFld( i ,j-1)+uFld(i+1,j-1) ) ) ENDDO ENDDO ELSE -C Original discretization +C- Using averaged transport: DO j=1-OLy+1,sNy+OLy DO i=1-OLx,sNx+OLx-1 - vCoriolisTerm(i,j) = - & -0.5*(_fCori(i, j ,bi,bj)+_fCori(i,j-1,bi,bj)) - & *0.25*( uFld(i, j )+uFld(i+1, j ) - & +uFld(i,j-1)+uFld(i+1,j-1) - & ) + vCoriolisTerm(i,j) = -0.5 _d 0 + & *( _fCori(i, j ,bi,bj) + _fCori(i,j-1,bi,bj) ) + & *( uFld( i , j )*dyG( i , j ,bi,bj)*hFacW( i , j ,k,bi,bj) + & + uFld(i+1, j )*dyG(i+1, j ,bi,bj)*hFacW(i+1, j ,k,bi,bj) + & + uFld( i ,j-1)*dyG( i ,j-1,bi,bj)*hFacW( i ,j-1,k,bi,bj) + & + uFld(i+1,j-1)*dyG(i+1,j-1,bi,bj)*hFacW(i+1,j-1,k,bi,bj) + & )*0.25 _d 0*recip_dyC(i,j,bi,bj)*recip_hFacS(i,j,k,bi,bj) ENDDO ENDDO ENDIF IF ( selectCoriScheme.EQ.1 .OR. selectCoriScheme.EQ.3 ) THEN -C Scale term so that only "wet" points are used -C Due to: Jamart and Ozer, 1986, JGR 91 (C9), 10,621-10,631 +C- Scale term so that only "wet" points are used +C Due to: Jamart and Ozer, 1986, JGR 91 (C9), 10,621-10,631 C "Numerical Boundary Layers and Spurious Residual Flows" DO j=1-OLy+1,sNy+OLy DO i=1-OLx,sNx+OLx-1 vCoriolisTerm(i,j) = vCoriolisTerm(i,j) - & *4. _d 0/MAX( one, + & *4. _d 0/MAX( oneRS, & maskW( i , j ,k,bi,bj)+maskW(i+1, j ,k,bi,bj) & +maskW( i ,j-1,k,bi,bj)+maskW(i+1,j-1,k,bi,bj) ) ENDDO diff --git a/pkg/mom_fluxform/mom_v_metric_sphere.F b/pkg/mom_fluxform/mom_v_metric_sphere.F index 63b8438649..756b9f73dc 100644 --- a/pkg/mom_fluxform/mom_v_metric_sphere.F +++ b/pkg/mom_fluxform/mom_v_metric_sphere.F @@ -5,7 +5,7 @@ C !INTERFACE: ========================================================== SUBROUTINE MOM_V_METRIC_SPHERE( - I bi,bj,k, + I bi, bj, k, I uFld, O vMetricTerms, I myThid ) @@ -24,11 +24,11 @@ SUBROUTINE MOM_V_METRIC_SPHERE( #include "GRID.h" C !INPUT PARAMETERS: =================================================== -C bi,bj :: tile indices +C bi, bj :: tile indices C k :: vertical level -C uFld :: zonal flow -C myThid :: thread number - INTEGER bi,bj,k +C uFld :: zonal velocity +C myThid :: my Thread Id number + INTEGER bi, bj, k _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER myThid @@ -37,22 +37,46 @@ SUBROUTINE MOM_V_METRIC_SPHERE( _RL vMetricTerms(1-OLx:sNx+OLx,1-OLy:sNy+OLy) C !LOCAL VARIABLES: ==================================================== -C i,j :: loop indices - INTEGER i,j +C i, j :: loop indices + INTEGER i, j CEOP - DO j=1-Olx+1,sNy+Oly - DO i=1-Olx,sNx+Olx-1 - vMetricTerms(i,j) = -recip_rSphere*recip_deepFacC(k) - & *0.25*( uFld(i,j )+uFld(i+1,j ) - & +uFld(i,j-1)+uFld(i+1,j-1) - & ) - & *0.25*( uFld(i,j )+uFld(i+1,j ) - & +uFld(i,j-1)+uFld(i+1,j-1) - & ) - & *_tanPhiAtV(i,j,bi,bj) - ENDDO - ENDDO + IF ( selectMetricTerms.EQ.1 ) THEN +C- Use analytical expression for tan(Phi) (stored in: tanPhiAtV) + DO j=1-OLx+1,sNy+OLy + DO i=1-OLx,sNx+OLx-1 + vMetricTerms(i,j) = -recip_rSphere*recip_deepFacC(k) + & *0.25*( uFld(i,j )+uFld(i+1,j ) + & +uFld(i,j-1)+uFld(i+1,j-1) + & ) + & *0.25*( uFld(i,j )+uFld(i+1,j ) + & +uFld(i,j-1)+uFld(i+1,j-1) + & ) + & *_tanPhiAtV(i,j,bi,bj) + ENDDO + ENDDO + ELSE +C- Using grid-spacing gradient: estimates tan(Phi)/rSphere as -del^j(dxC)/rAz + DO j=2-OLx,sNy+OLy + DO i=1-OLx,sNx+OLx-1 + vMetricTerms(i,j) = +0.5 _d 0* + & ( ( uFld(i,j-1) + uFld(i,j) )*halfRL + & *( dxC( i, j, bi,bj) + & - dxC( i, j-1,bi,bj) )*recip_rAz( i, j,bi,bj) + & *( uFld( i, j-1)*dyG( i, j-1,bi,bj)*hFacW( i, j-1,k,bi,bj) + & + uFld( i, j )*dyG( i, j, bi,bj)*hFacW( i, j, k,bi,bj) + & )*halfRL + & + ( uFld(i+1,j-1) + uFld(i+1,j) )*halfRL + & *( dxC(i+1, j, bi,bj) + & - dxC(i+1,j-1,bi,bj) )*recip_rAz(i+1,j,bi,bj) + & *( uFld(i+1,j-1)*dyG(i+1,j-1,bi,bj)*hFacW(i+1,j-1,k,bi,bj) + & + uFld(i+1, j )*dyG(i+1, j, bi,bj)*hFacW(i+1, j, k,bi,bj) + & )*halfRL + & )*recip_hFacS(i,j,k,bi,bj) + & *recip_dyC(i,j,bi,bj)*recip_deepFacC(k) + ENDDO + ENDDO + ENDIF RETURN END diff --git a/pkg/mom_vecinv/mom_vecinv.F b/pkg/mom_vecinv/mom_vecinv.F index 13d7d9b63e..a79e94e9c2 100644 --- a/pkg/mom_vecinv/mom_vecinv.F +++ b/pkg/mom_vecinv/mom_vecinv.F @@ -785,13 +785,15 @@ SUBROUTINE MOM_VECINV( C-- Vertical shear terms (-w*du/dr & -w*dv/dr) IF ( .NOT. momImplVertAdv ) THEN - CALL MOM_VI_U_VERTSHEAR(bi,bj,k,uVel,wVel,uCf,myThid) + CALL MOM_VI_U_VERTSHEAR( bi, bj, k, deepFacAdv, + & uVel, wVel, uCf, myThid ) DO j=jMin,jMax DO i=iMin,iMax gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j) ENDDO ENDDO - CALL MOM_VI_V_VERTSHEAR(bi,bj,k,vVel,wVel,vCf,myThid) + CALL MOM_VI_V_VERTSHEAR( bi, bj, k, deepFacAdv, + & vVel, wVel, vCf, myThid ) DO j=jMin,jMax DO i=iMin,iMax gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j) @@ -856,19 +858,25 @@ SUBROUTINE MOM_VECINV( ENDIF C-- Non-Hydrostatic (spherical) metric terms +#ifdef MOM_USE_OLD_DEEP_VERT_ADV IF ( useNHMTerms ) THEN +#else + IF ( useNHMTerms .AND. .NOT.deepAtmosphere ) THEN +#endif CALL MOM_U_METRIC_NH(bi,bj,k,uFld,wVel,uCf,myThid) - DO j=jMin,jMax - DO i=iMin,iMax - gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+uCf(i,j) - ENDDO - ENDDO CALL MOM_V_METRIC_NH(bi,bj,k,vFld,wVel,vCf,myThid) DO j=jMin,jMax DO i=iMin,iMax - gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+vCf(i,j) + gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + uCf(i,j) + gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + vCf(i,j) ENDDO ENDDO +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN + CALL DIAGNOSTICS_FILL( uCf, 'Um_Metr ', k,1,2, bi,bj, myThid ) + CALL DIAGNOSTICS_FILL( vCf, 'Vm_Metr ', k,1,2, bi,bj, myThid ) + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ ENDIF C-- Set du/dt & dv/dt on boundaries to zero diff --git a/pkg/mom_vecinv/mom_vi_u_coriolis.F b/pkg/mom_vecinv/mom_vi_u_coriolis.F index 644536af92..846b1cf4bd 100644 --- a/pkg/mom_vecinv/mom_vi_u_coriolis.F +++ b/pkg/mom_vecinv/mom_vi_u_coriolis.F @@ -4,7 +4,7 @@ C !ROUTINE: MOM_VI_U_CORIOLIS C !INTERFACE: SUBROUTINE MOM_VI_U_CORIOLIS( - I bi, bj, k, + I bi, bj, k, I selectVortScheme, useJamartMomAdv, I vFld, omega3, hFacZ, r_hFacZ, O uCoriolisTerm, @@ -40,16 +40,14 @@ SUBROUTINE MOM_VI_U_CORIOLIS( CEOP C == Local variables == -C msgBuf :: Informational/error meesage buffer +C msgBuf :: Informational/error message buffer CHARACTER*(MAX_LEN_MBUF) msgBuf - LOGICAL upwindVort3 INTEGER i, j _RL vBarXY, vBarXm, vBarXp _RL vort3u _RL vort3mj, vort3ij, vort3mp, vort3ip _RL oneThird, tmpFac _RS epsil - PARAMETER( upwindVort3 = .FALSE. ) epsil = 1. _d -9 tmpFac = 1. _d 0 @@ -59,78 +57,54 @@ SUBROUTINE MOM_VI_U_CORIOLIS( IF ( selectVortScheme.EQ.0 ) THEN C-- using enstrophy conserving scheme (Shallow-Water Eq.) by Sadourny, JAS 75 - DO j=1-Oly,sNy+Oly-1 - DO i=2-Olx,sNx+Olx - vBarXY=0.25*( - & (vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj) - & +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj)) - & +(vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj) - & +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj)) - & ) - IF (upwindVort3) THEN - IF (vBarXY.GT.0.) THEN - vort3u=omega3(i,j)*r_hFacZ(i,j) - ELSE - vort3u=omega3(i,j+1)*r_hFacZ(i,j+1) - ENDIF - ELSE - vort3u=0.5*(omega3(i,j)*r_hFacZ(i,j) - & +omega3(i,j+1)*r_hFacZ(i,j+1)) - ENDIF - uCoriolisTerm(i,j)= +vort3u*vBarXY*recip_dxC(i,j,bi,bj) - & * _maskW(i,j,k,bi,bj) + DO j=1-OLy,sNy+OLy-1 + DO i=2-OLx,sNx+OLx + vBarXY = 0.25 _d 0*( + & ( vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj) + & + vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) ) + & +( vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj) + & + vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) ) + & ) + vort3u = 0.5 _d 0*( omega3(i,j)*r_hFacZ(i,j) + & + omega3(i,j+1)*r_hFacZ(i,j+1) ) + uCoriolisTerm(i,j) = +vort3u*vBarXY*recip_dxC(i,j,bi,bj) + & * _maskW(i,j,k,bi,bj) ENDDO ENDDO ELSEIF ( selectVortScheme.EQ.1 ) THEN C-- same as above, with different formulation (relatively to hFac) - DO j=1-Oly,sNy+Oly-1 - DO i=2-Olx,sNx+Olx - vBarXY= 0.5*( - & (vFld( i , j )*dxG( i , j ,bi,bj)*hFacZ(i, j ) - & +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*hFacZ(i, j )) - & +(vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*hFacZ(i,j+1) - & +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*hFacZ(i,j+1)) - & )/MAX( epsil, hFacZ(i,j)+hFacZ(i,j+1) ) - IF (upwindVort3) THEN - IF (vBarXY.GT.0.) THEN - vort3u=omega3(i,j) - ELSE - vort3u=omega3(i,j+1) - ENDIF - ELSE - vort3u=0.5*(omega3(i,j)+omega3(i,j+1)) - ENDIF - uCoriolisTerm(i,j)= +vort3u*vBarXY*recip_dxC(i,j,bi,bj) - & * _maskW(i,j,k,bi,bj) + DO j=1-OLy,sNy+OLy-1 + DO i=2-OLx,sNx+OLx + vBarXY = 0.5 _d 0*( + & ( vFld( i , j )*dxG( i , j ,bi,bj)*hFacZ(i, j ) + & + vFld(i-1, j )*dxG(i-1, j ,bi,bj)*hFacZ(i, j ) ) + & +( vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*hFacZ(i,j+1) + & + vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*hFacZ(i,j+1) ) + & )/MAX( epsil, hFacZ(i,j)+hFacZ(i,j+1) ) + vort3u = 0.5 _d 0*( omega3(i,j) + omega3(i,j+1) ) + uCoriolisTerm(i,j) = +vort3u*vBarXY*recip_dxC(i,j,bi,bj) + & * _maskW(i,j,k,bi,bj) ENDDO ENDDO ELSEIF ( selectVortScheme.EQ.2 ) THEN C-- using energy conserving scheme (used by Sadourny in JAS 75 paper) - DO j=1-Oly,sNy+Oly-1 - DO i=2-Olx,sNx+Olx - vBarXm=0.5*( + DO j=1-OLy,sNy+OLy-1 + DO i=2-OLx,sNx+OLx + vBarXm = 0.5 _d 0*( & vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj) - & +vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) ) - vBarXp=0.5*( + & + vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) ) + vBarXp = 0.5 _d 0*( & vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj) - & +vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) ) - IF (upwindVort3) THEN - IF ( (vBarXm+vBarXp) .GT.0.) THEN - vort3u=vBarXm*r_hFacZ(i, j )*omega3(i, j ) - ELSE - vort3u=vBarXp*r_hFacZ(i,j+1)*omega3(i,j+1) - ENDIF - ELSE - vort3u = ( vBarXm*r_hFacZ(i, j )*omega3(i, j ) - & +vBarXp*r_hFacZ(i,j+1)*omega3(i,j+1) - & )*0.5 _d 0 - ENDIF - uCoriolisTerm(i,j)= +vort3u*recip_dxC(i,j,bi,bj) - & * _maskW(i,j,k,bi,bj) + & + vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) ) + vort3u = ( vBarXm*r_hFacZ(i, j )*omega3(i, j ) + & + vBarXp*r_hFacZ(i,j+1)*omega3(i,j+1) + & )*0.5 _d 0 + uCoriolisTerm(i,j) = +vort3u*recip_dxC(i,j,bi,bj) + & * _maskW(i,j,k,bi,bj) ENDDO ENDDO @@ -141,36 +115,55 @@ SUBROUTINE MOM_VI_U_CORIOLIS( C domain where uCoriolisTerm is valid : C [ 3-Olx : sNx+Olx-1 ] x [ 2-Oly : sNy+Oly-1 ] C (=> might need overlap of 3 if using CD-scheme) - DO j=1-Oly,sNy+Oly-1 - DO i=2-Olx,sNx+Olx-1 - vort3mj= ( r_hFacZ(i, j )*omega3(i, j ) - & +(r_hFacZ(i,j+1)*omega3(i,j+1) - & +r_hFacZ(i-1,j)*omega3(i-1,j) - & ))*oneThird -c & )*tmpFac)*oneThird + DO j=1-OLy,sNy+OLy-1 + DO i=2-OLx,sNx+OLx-1 + vort3mj = ( r_hFacZ(i, j )*omega3(i, j ) + & +(r_hFacZ(i,j+1)*omega3(i,j+1) + & +r_hFacZ(i-1,j)*omega3(i-1,j) + & ))*oneThird +c & )*tmpFac)*oneThird & *vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) - vort3ij= ( r_hFacZ(i, j )*omega3(i, j ) - & +(r_hFacZ(i,j+1)*omega3(i,j+1) - & +r_hFacZ(i+1,j)*omega3(i+1,j) - & ))*oneThird -c & )*tmpFac)*oneThird + vort3ij = ( r_hFacZ(i, j )*omega3(i, j ) + & +(r_hFacZ(i,j+1)*omega3(i,j+1) + & +r_hFacZ(i+1,j)*omega3(i+1,j) + & ))*oneThird +c & )*tmpFac)*oneThird & *vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj) - vort3mp= ( r_hFacZ(i,j+1)*omega3(i,j+1) - & +(r_hFacZ(i, j )*omega3(i, j ) - & +r_hFacZ(i-1,j+1)*omega3(i-1,j+1) - & ))*oneThird -c & )*tmpFac)*oneThird + vort3mp = ( r_hFacZ(i,j+1)*omega3(i,j+1) + & +(r_hFacZ(i, j )*omega3(i, j ) + & +r_hFacZ(i-1,j+1)*omega3(i-1,j+1) + & ))*oneThird +c & )*tmpFac)*oneThird & *vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) - vort3ip= ( r_hFacZ(i,j+1)*omega3(i,j+1) - & +(r_hFacZ(i, j )*omega3(i, j ) - & +r_hFacZ(i+1,j+1)*omega3(i+1,j+1) - & ))*oneThird -c & )*tmpFac)*oneThird + vort3ip = ( r_hFacZ(i,j+1)*omega3(i,j+1) + & +(r_hFacZ(i, j )*omega3(i, j ) + & +r_hFacZ(i+1,j+1)*omega3(i+1,j+1) + & ))*oneThird +c & )*tmpFac)*oneThird & *vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj) C--- - uCoriolisTerm(i,j)= +( (vort3mj+vort3ij)+(vort3mp+vort3ip) ) - & *0.25 _d 0 *recip_dxC(i,j,bi,bj) - & * _maskW(i,j,k,bi,bj) + uCoriolisTerm(i,j) = +( (vort3mj+vort3ij)+(vort3mp+vort3ip) ) + & *0.25 _d 0 *recip_dxC(i,j,bi,bj) + & * _maskW(i,j,k,bi,bj) + ENDDO + ENDDO + + ELSEIF ( selectVortScheme.EQ.4 ) THEN +C-- using energy conserving scheme without r_hFacZ factor + + DO j=1-OLy,sNy+OLy-1 + DO i=2-OLx,sNx+OLx + vBarXm = 0.5 _d 0*( + & vFld( i , j )*dxG( i , j ,bi,bj)*_hFacS( i , j ,k,bi,bj) + & + vFld(i-1, j )*dxG(i-1, j ,bi,bj)*_hFacS(i-1, j ,k,bi,bj) ) + vBarXp = 0.5 _d 0*( + & vFld( i ,j+1)*dxG( i ,j+1,bi,bj)*_hFacS( i ,j+1,k,bi,bj) + & + vFld(i-1,j+1)*dxG(i-1,j+1,bi,bj)*_hFacS(i-1,j+1,k,bi,bj) ) + vort3u = ( vBarXm*omega3(i, j ) + & + vBarXp*omega3(i,j+1) + & )*0.5 _d 0 + uCoriolisTerm(i,j) = +vort3u*recip_dxC(i,j,bi,bj) + & *recip_hFacW(i,j,k,bi,bj) ENDDO ENDDO @@ -184,8 +177,8 @@ SUBROUTINE MOM_VI_U_CORIOLIS( ENDIF IF ( useJamartMomAdv ) THEN - DO j=1-Oly,sNy+Oly-1 - DO i=2-Olx,sNx+Olx-1 + DO j=1-OLy,sNy+OLy-1 + DO i=2-OLx,sNx+OLx-1 uCoriolisTerm(i,j) = uCoriolisTerm(i,j) & * 4. _d 0 * _hFacW(i,j,k,bi,bj) & / MAX( epsil, diff --git a/pkg/mom_vecinv/mom_vi_u_vertshear.F b/pkg/mom_vecinv/mom_vi_u_vertshear.F index 7a2ca206de..7fb6270309 100644 --- a/pkg/mom_vecinv/mom_vi_u_vertshear.F +++ b/pkg/mom_vecinv/mom_vi_u_vertshear.F @@ -1,57 +1,63 @@ #include "MOM_VECINV_OPTIONS.h" +CBOP +C !ROUTINE: MOM_VI_U_VERTSHEAR + +C !INTERFACE: SUBROUTINE MOM_VI_U_VERTSHEAR( - I bi,bj,K, - I uFld,wFld, + I bi, bj, k, deepFacA, + I uFld, wFld, U uShearTerm, - I myThid) - IMPLICIT NONE + I myThid ) + +C !DESCRIPTION: C *==========================================================* C | S/R MOM_U_VERTSHEAR C *==========================================================* -C *==========================================================* +C !USES: + IMPLICIT NONE C == Global variables == #include "SIZE.h" #include "EEPARAMS.h" #include "GRID.h" #include "PARAMS.h" -C == Routine arguments == - INTEGER bi,bj,K +C !INPUT/OUTPUT PARAMETERS: +C deepFacA :: deep-model grid factor at level center + INTEGER bi, bj, k + _RL deepFacA(Nr) _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL wFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL uShearTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER myThid -C == Local variables == - INTEGER I,J,Kp1,Km1 - _RL mask_Kp1,mask_Km1,wBarXm,wBarXp - _RL uZm,uZp - LOGICAL rAdvAreaWeight -c _RL umask_Kp1,umask_K,umask_Km1 -c LOGICAL freeslipK,noslipK -c PARAMETER(freeslipK=.TRUE.) -c PARAMETER(noslipK=.NOT.freeslipK) -c LOGICAL freeslip1,noslip1 -c PARAMETER(freeslip1=.TRUE.) -c PARAMETER(noslip1=.NOT.freeslip1) -c1 _RL wBarXZ,uZbarZ +C !LOCAL VARIABLES: + INTEGER i, j, kp1, km1 + _RL mask_Kp1, mask_Km1, wBarXm, wBarXp + _RL uZm, uZp, recip_drDeepRho + LOGICAL rAdvAreaWeight +c _RL umask_Kp1,umask_K,umask_Km1 +c1 _RL wBarXZ,uZbarZ +CEOP rAdvAreaWeight =.TRUE. C- Area-weighted average either in KE or in vert. advection: IF ( selectKEscheme.EQ.1 .OR. selectKEscheme.EQ.3 ) & rAdvAreaWeight =.FALSE. - Kp1=min(K+1,Nr) - mask_Kp1=1. - IF (K.EQ.Nr) mask_Kp1=0. - Km1=max(K-1,1) - mask_Km1=1. - IF (K.EQ.1) mask_Km1=0. + kp1 = MIN(k+1,Nr) + mask_Kp1 = oneRL + IF (k.EQ.Nr) mask_Kp1 = zeroRL + km1 = MAX(k-1,1) + mask_Km1 = oneRL + IF (k.EQ.1) mask_Km1 = zeroRL + + recip_drDeepRho = recip_drF(k)/deepFacA(k) + & * recip_deepFac2C(k)*recip_rhoFacC(k) - DO J=1-OLy,sNy+OLy - DO I=2-OLx,sNx+OLx + DO j=1-OLy,sNy+OLy + DO i=2-OLx,sNx+OLx c umask_K=_maskW(i,j,k,bi,bj) @@ -61,43 +67,45 @@ SUBROUTINE MOM_VI_U_VERTSHEAR( c & *mask_Kp1 IF ( rAdvAreaWeight ) THEN -C Transport at interface k : Area weighted average - wBarXm=0.5*( - & wFld(I,J,K,bi,bj)*rA(i,j,bi,bj)*maskC(I,J,Km1,bi,bj) - & +wFld(I-1,J,K,bi,bj)*rA(i-1,j,bi,bj)*maskC(I-1,J,Km1,bi,bj) - & )*mask_Km1*deepFac2F(K)*rhoFacF(K) - & *recip_rAw(i,j,bi,bj) - -C Transport at interface k+1 (here wFld is already masked) - wBarXp=0.5*( - & wFld(I,J,Kp1,bi,bj)*rA(i,j,bi,bj) - & +wFld(I-1,J,Kp1,bi,bj)*rA(i-1,j,bi,bj) - & )*mask_Kp1*deepFac2F(Kp1)*rhoFacF(Kp1) - & *recip_rAw(i,j,bi,bj) +C Transport at interface k : Area weighted average + wBarXm = halfRL*( + & wFld(i,j,k,bi,bj)*rA(i,j,bi,bj)*maskC(i,j,km1,bi,bj) + & + wFld(i-1,j,k,bi,bj)*rA(i-1,j,bi,bj)*maskC(i-1,j,km1,bi,bj) + & )*mask_Km1*deepFac2F(k)*rhoFacF(k) + & *recip_rAw(i,j,bi,bj) + +C Transport at interface k+1 (here wFld is already masked) + wBarXp = halfRL*( + & wFld(i,j,kp1,bi,bj)*rA(i,j,bi,bj) + & + wFld(i-1,j,kp1,bi,bj)*rA(i-1,j,bi,bj) + & )*mask_Kp1*deepFac2F(kp1)*rhoFacF(kp1) + & *recip_rAw(i,j,bi,bj) ELSE -C Transport at interface k : simple average - wBarXm=0.5*( - & wFld(I,J,K,bi,bj)*maskC(I,J,Km1,bi,bj) - & +wFld(I-1,J,K,bi,bj)*maskC(I-1,J,Km1,bi,bj) - & )*mask_Km1*deepFac2F(K)*rhoFacF(K) - -C Transport at interface k+1 (here wFld is already masked) - wBarXp=0.5*( - & wFld(I,J,Kp1,bi,bj) - & +wFld(I-1,J,Kp1,bi,bj) - & )*mask_Kp1*deepFac2F(Kp1)*rhoFacF(Kp1) +C Transport at interface k : simple average + wBarXm = halfRL*( + & wFld(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj) + & + wFld(i-1,j,k,bi,bj)*maskC(i-1,j,km1,bi,bj) + & )*mask_Km1*deepFac2F(k)*rhoFacF(k) + +C Transport at interface k+1 (here wFld is already masked) + wBarXp = halfRL*( + & wFld(i,j,kp1,bi,bj) + & + wFld(i-1,j,kp1,bi,bj) + & )*mask_Kp1*deepFac2F(kp1)*rhoFacF(kp1) ENDIF -C delta_Z( U ) @ interface k +C- delta_Z( U*deepFac ) @ interface k c umask_Km1=mask_Km1*maskW(i,j,Km1,bi,bj) - uZm=(uFld(I,J,K,bi,bj)-mask_Km1*uFld(I,J,Km1,bi,bj))*rkSign + uZm = ( uFld(i,j, k ,bi,bj)*deepFacA( k ) + & - uFld(i,j,km1,bi,bj)*deepFacA(km1)*mask_Km1 )*rkSign c2 & *recip_dRC(K) c IF (freeslip1) uZm=uZm*umask_Km1 c IF (noslip1.AND.umask_Km1.EQ.0.) uZm=uZm*2. -C delta_Z( U ) @ interface k+1 +C- delta_Z( U*deepFac ) @ interface k+1 c umask_Kp1=mask_Kp1*maskW(i,j,Kp1,bi,bj) - uZp=(mask_Kp1*uFld(I,J,Kp1,bi,bj)-uFld(I,J,K,bi,bj))*rkSign + uZp = ( uFld(i,j,kp1,bi,bj)*deepFacA(kp1)*mask_Kp1 + & - uFld(i,j, k ,bi,bj)*deepFacA( k ) )*rkSign c2 & *recip_dRC(Kp1) c IF (freeslipK) uZp=uZp*umask_Kp1 c IF (noslipK.AND.umask_Kp1.EQ.0.) uZp=uZp*2. @@ -117,17 +125,13 @@ SUBROUTINE MOM_VI_U_VERTSHEAR( c2 uShearTerm(I,J)=-0.5*(wBarXp*uZp+wBarXm*uZm) c2 & *_maskW(I,J,K,bi,bj) IF (upwindShear) THEN - uShearTerm(I,J)=-0.5* - & ( (wBarXp*uZp+wBarXm*uZm) - & +(ABS(wBarXp)*uZp-ABS(wBarXm)*uZm) - & )*_recip_hFacW(i,j,k,bi,bj) - & * recip_drF(K) - & * recip_deepFac2C(K)*recip_rhoFacC(K) + uShearTerm(i,j) = -halfRL* + & ( ( wBarXp *uZp + wBarXm *uZm ) + & + ( ABS(wBarXp)*uZp - ABS(wBarXm)*uZm ) + & )*_recip_hFacW(i,j,k,bi,bj)*recip_drDeepRho ELSE - uShearTerm(I,J)=-0.5*(wBarXp*uZp+wBarXm*uZm) - & *_recip_hFacW(i,j,k,bi,bj) - & * recip_drF(K) - & * recip_deepFac2C(K)*recip_rhoFacC(K) + uShearTerm(i,j) = -halfRL*( wBarXp*uZp + wBarXm*uZm ) + & *_recip_hFacW(i,j,k,bi,bj)*recip_drDeepRho ENDIF ENDDO ENDDO diff --git a/pkg/mom_vecinv/mom_vi_v_coriolis.F b/pkg/mom_vecinv/mom_vi_v_coriolis.F index e4629485d5..fa06f5ab5c 100644 --- a/pkg/mom_vecinv/mom_vi_v_coriolis.F +++ b/pkg/mom_vecinv/mom_vi_v_coriolis.F @@ -40,16 +40,14 @@ SUBROUTINE MOM_VI_V_CORIOLIS( CEOP C == Local variables == -C msgBuf :: Informational/error meesage buffer +C msgBuf :: Informational/error message buffer CHARACTER*(MAX_LEN_MBUF) msgBuf - LOGICAL upwindVort3 INTEGER i, j _RL uBarXY, uBarYm, uBarYp _RL vort3v _RL vort3im, vort3ij, vort3pm, vort3pj _RL oneThird, tmpFac _RS epsil - PARAMETER( upwindVort3 = .FALSE. ) epsil = 1. _d -9 tmpFac = 1. _d 0 @@ -59,78 +57,54 @@ SUBROUTINE MOM_VI_V_CORIOLIS( IF ( selectVortScheme.EQ.0 ) THEN C-- using enstrophy conserving scheme (Shallow-Water Eq.) by Sadourny, JAS 75 - DO j=2-Oly,sNy+Oly - DO i=1-Olx,sNx+Olx-1 - uBarXY=0.25*( - & (uFld( i , j )*dyG( i , j ,bi,bj)*_hFacW( i , j ,k,bi,bj) - & +uFld( i ,j-1)*dyG( i ,j-1,bi,bj)*_hFacW( i ,j-1,k,bi,bj)) - & +(uFld(i+1, j )*dyG(i+1, j ,bi,bj)*_hFacW(i+1, j ,k,bi,bj) - & +uFld(i+1,j-1)*dyG(i+1,j-1,bi,bj)*_hFacW(i+1,j-1,k,bi,bj)) - & ) - IF (upwindVort3) THEN - IF (uBarXY.GT.0.) THEN - vort3v=omega3(i,j)*r_hFacZ(i,j) - ELSE - vort3v=omega3(i+1,j)*r_hFacZ(i+1,j) - ENDIF - ELSE - vort3v=0.5*(omega3(i,j)*r_hFacZ(i,j) - & +omega3(i+1,j)*r_hFacZ(i+1,j)) - ENDIF - vCoriolisTerm(i,j)= -vort3v*uBarXY*recip_dyC(i,j,bi,bj) - & * _maskS(i,j,k,bi,bj) + DO j=2-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx-1 + uBarXY = 0.25 _d 0 *( + & ( uFld( i , j )*dyG( i , j ,bi,bj)*_hFacW( i , j ,k,bi,bj) + & + uFld( i ,j-1)*dyG( i ,j-1,bi,bj)*_hFacW( i ,j-1,k,bi,bj) ) + & +( uFld(i+1, j )*dyG(i+1, j ,bi,bj)*_hFacW(i+1, j ,k,bi,bj) + & + uFld(i+1,j-1)*dyG(i+1,j-1,bi,bj)*_hFacW(i+1,j-1,k,bi,bj) ) + & ) + vort3v = 0.5 _d 0*( omega3(i,j)*r_hFacZ(i,j) + & + omega3(i+1,j)*r_hFacZ(i+1,j) ) + vCoriolisTerm(i,j) = -vort3v*uBarXY*recip_dyC(i,j,bi,bj) + & * _maskS(i,j,k,bi,bj) ENDDO ENDDO ELSEIF ( selectVortScheme.EQ.1 ) THEN C-- same as above, with different formulation (relatively to hFac) - DO j=2-Oly,sNy+Oly - DO i=1-Olx,sNx+Olx-1 - uBarXY= 0.5*( - & (uFld( i , j )*dyG( i , j ,bi,bj)*hFacZ( i ,j) - & +uFld( i ,j-1)*dyG( i ,j-1,bi,bj)*hFacZ( i ,j)) - & +(uFld(i+1, j )*dyG(i+1, j ,bi,bj)*hFacZ(i+1,j) - & +uFld(i+1,j-1)*dyG(i+1,j-1,bi,bj)*hFacZ(i+1,j)) - & )/MAX( epsil, hFacZ(i,j)+hFacZ(i+1,j) ) - IF (upwindVort3) THEN - IF (uBarXY.GT.0.) THEN - vort3v=omega3(i,j) - ELSE - vort3v=omega3(i+1,j) - ENDIF - ELSE - vort3v=0.5*(omega3(i,j)+omega3(i+1,j)) - ENDIF - vCoriolisTerm(i,j)= -vort3v*uBarXY*recip_dyC(i,j,bi,bj) - & * _maskS(i,j,k,bi,bj) + DO j=2-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx-1 + uBarXY = 0.5 _d 0*( + & ( uFld( i , j )*dyG( i , j ,bi,bj)*hFacZ( i ,j) + & + uFld( i ,j-1)*dyG( i ,j-1,bi,bj)*hFacZ( i ,j) ) + & +( uFld(i+1, j )*dyG(i+1, j ,bi,bj)*hFacZ(i+1,j) + & + uFld(i+1,j-1)*dyG(i+1,j-1,bi,bj)*hFacZ(i+1,j) ) + & )/MAX( epsil, hFacZ(i,j)+hFacZ(i+1,j) ) + vort3v = 0.5 _d 0*( omega3(i,j) + omega3(i+1,j) ) + vCoriolisTerm(i,j) = -vort3v*uBarXY*recip_dyC(i,j,bi,bj) + & * _maskS(i,j,k,bi,bj) ENDDO ENDDO ELSEIF ( selectVortScheme.EQ.2 ) THEN C-- using energy conserving scheme (used by Sadourny in JAS 75 paper) - DO j=2-Oly,sNy+Oly - DO i=1-Olx,sNx+Olx-1 - uBarYm=0.5*( + DO j=2-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx-1 + uBarYm = 0.5 _d 0*( & uFld( i , j )*dyG( i , j ,bi,bj)*_hFacW( i , j ,k,bi,bj) - & +uFld( i ,j-1)*dyG( i ,j-1,bi,bj)*_hFacW( i ,j-1,k,bi,bj) ) - uBarYp=0.5*( + & + uFld( i ,j-1)*dyG( i ,j-1,bi,bj)*_hFacW( i ,j-1,k,bi,bj) ) + uBarYp = 0.5 _d 0*( & uFld(i+1, j )*dyG(i+1, j ,bi,bj)*_hFacW(i+1, j ,k,bi,bj) - & +uFld(i+1,j-1)*dyG(i+1,j-1,bi,bj)*_hFacW(i+1,j-1,k,bi,bj) ) - IF (upwindVort3) THEN - IF ( (uBarYm+uBarYp) .GT.0.) THEN - vort3v=uBarYm*r_hFacZ( i ,j)*omega3( i ,j) - ELSE - vort3v=uBarYp*r_hFacZ(i+1,j)*omega3(i+1,j) - ENDIF - ELSE - vort3v = ( uBarYm*r_hFacZ( i ,j)*omega3( i ,j) - & +uBarYp*r_hFacZ(i+1,j)*omega3(i+1,j) - & )*0.5 _d 0 - ENDIF - vCoriolisTerm(i,j)= -vort3v*recip_dyC(i,j,bi,bj) - & * _maskS(i,j,k,bi,bj) + & + uFld(i+1,j-1)*dyG(i+1,j-1,bi,bj)*_hFacW(i+1,j-1,k,bi,bj) ) + vort3v = ( uBarYm*r_hFacZ( i ,j)*omega3( i ,j) + & + uBarYp*r_hFacZ(i+1,j)*omega3(i+1,j) + & )*0.5 _d 0 + vCoriolisTerm(i,j) = -vort3v*recip_dyC(i,j,bi,bj) + & * _maskS(i,j,k,bi,bj) ENDDO ENDDO @@ -141,36 +115,55 @@ SUBROUTINE MOM_VI_V_CORIOLIS( C domain where vCoriolisTerm is valid : C [ 2-Olx : sNx+Olx-1 ] x [ 3-Oly : sNy+Oly-1 ] C (=> might need overlap of 3 if using CD-scheme) - DO j=2-Oly,sNy+Oly-1 - DO i=1-Olx,sNx+Olx-1 - vort3im= ( r_hFacZ(i, j )*omega3(i, j ) - & +(r_hFacZ(i+1,j)*omega3(i+1,j) - & +r_hFacZ(i,j-1)*omega3(i,j-1) - & ))*oneThird -c & )*tmpFac)*oneThird + DO j=2-OLy,sNy+OLy-1 + DO i=1-OLx,sNx+OLx-1 + vort3im = ( r_hFacZ(i, j )*omega3(i, j ) + & +(r_hFacZ(i+1,j)*omega3(i+1,j) + & +r_hFacZ(i,j-1)*omega3(i,j-1) + & ))*oneThird +c & )*tmpFac)*oneThird & *uFld( i ,j-1)*dyG( i ,j-1,bi,bj)*_hFacW( i ,j-1,k,bi,bj) - vort3ij= ( r_hFacZ(i, j )*omega3(i, j ) - & +(r_hFacZ(i+1,j)*omega3(i+1,j) - & +r_hFacZ(i,j+1)*omega3(i,j+1) - & ))*oneThird -c & )*tmpFac)*oneThird + vort3ij = ( r_hFacZ(i, j )*omega3(i, j ) + & +(r_hFacZ(i+1,j)*omega3(i+1,j) + & +r_hFacZ(i,j+1)*omega3(i,j+1) + & ))*oneThird +c & )*tmpFac)*oneThird & *uFld( i , j )*dyG( i , j ,bi,bj)*_hFacW( i , j ,k,bi,bj) - vort3pm= ( r_hFacZ(i+1,j)*omega3(i+1,j) - & +(r_hFacZ(i, j )*omega3(i, j ) - & +r_hFacZ(i+1,j-1)*omega3(i+1,j-1) - & ))*oneThird -c & )*tmpFac)*oneThird + vort3pm = ( r_hFacZ(i+1,j)*omega3(i+1,j) + & +(r_hFacZ(i, j )*omega3(i, j ) + & +r_hFacZ(i+1,j-1)*omega3(i+1,j-1) + & ))*oneThird +c & )*tmpFac)*oneThird & *uFld(i+1,j-1)*dyG(i+1,j-1,bi,bj)*_hFacW(i+1,j-1,k,bi,bj) - vort3pj= ( r_hFacZ(i+1,j)*omega3(i+1,j) - & +(r_hFacZ(i, j )*omega3(i, j ) - & +r_hFacZ(i+1,j+1)*omega3(i+1,j+1) - & ))*oneThird -c & )*tmpFac)*oneThird + vort3pj = ( r_hFacZ(i+1,j)*omega3(i+1,j) + & +(r_hFacZ(i, j )*omega3(i, j ) + & +r_hFacZ(i+1,j+1)*omega3(i+1,j+1) + & ))*oneThird +c & )*tmpFac)*oneThird & *uFld(i+1, j )*dyG(i+1, j ,bi,bj)*_hFacW(i+1, j ,k,bi,bj) C--- - vCoriolisTerm(i,j)= -( (vort3im+vort3ij)+(vort3pm+vort3pj) ) - & *0.25 _d 0 *recip_dyC(i,j,bi,bj) - & * _maskS(i,j,k,bi,bj) + vCoriolisTerm(i,j) = -( (vort3im+vort3ij)+(vort3pm+vort3pj) ) + & *0.25 _d 0 *recip_dyC(i,j,bi,bj) + & * _maskS(i,j,k,bi,bj) + ENDDO + ENDDO + + ELSEIF ( selectVortScheme.EQ.4 ) THEN +C-- using energy conserving scheme without r_hFacZ factor + + DO j=2-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx-1 + uBarYm = 0.5 _d 0*( + & uFld( i , j )*dyG( i , j ,bi,bj)*_hFacW( i , j ,k,bi,bj) + & + uFld( i ,j-1)*dyG( i ,j-1,bi,bj)*_hFacW( i ,j-1,k,bi,bj) ) + uBarYp = 0.5 _d 0*( + & uFld(i+1, j )*dyG(i+1, j ,bi,bj)*_hFacW(i+1, j ,k,bi,bj) + & + uFld(i+1,j-1)*dyG(i+1,j-1,bi,bj)*_hFacW(i+1,j-1,k,bi,bj) ) + vort3v = ( uBarYm*omega3( i ,j) + & + uBarYp*omega3(i+1,j) + & )*0.5 _d 0 + vCoriolisTerm(i,j) = -vort3v*recip_dyC(i,j,bi,bj) + & *recip_hFacS(i,j,k,bi,bj) ENDDO ENDDO @@ -184,8 +177,8 @@ SUBROUTINE MOM_VI_V_CORIOLIS( ENDIF IF ( useJamartMomAdv ) THEN - DO j=2-Oly,sNy+Oly-1 - DO i=1-Olx,sNx+Olx-1 + DO j=2-OLy,sNy+OLy-1 + DO i=1-OLx,sNx+OLx-1 vCoriolisTerm(i,j) = vCoriolisTerm(i,j) & * 4. _d 0 * _hFacS(i,j,k,bi,bj) & / MAX( epsil, diff --git a/pkg/mom_vecinv/mom_vi_v_vertshear.F b/pkg/mom_vecinv/mom_vi_v_vertshear.F index 9748c6fcc5..faa151cc4a 100644 --- a/pkg/mom_vecinv/mom_vi_v_vertshear.F +++ b/pkg/mom_vecinv/mom_vi_v_vertshear.F @@ -1,57 +1,63 @@ #include "MOM_VECINV_OPTIONS.h" +CBOP +C !ROUTINE: MOM_VI_V_VERTSHEAR + +C !INTERFACE: SUBROUTINE MOM_VI_V_VERTSHEAR( - I bi,bj,K, - I vFld,wFld, + I bi, bj, k, deepFacA, + I vFld, wFld, U vShearTerm, - I myThid) - IMPLICIT NONE + I myThid ) + +C !DESCRIPTION: C *==========================================================* C | S/R MOM_V_VERTSHEAR C *==========================================================* -C *==========================================================* +C !USES: + IMPLICIT NONE C == Global variables == #include "SIZE.h" #include "EEPARAMS.h" #include "GRID.h" #include "PARAMS.h" -C == Routine arguments == - INTEGER bi,bj,K +C !INPUT/OUTPUT PARAMETERS: +C deepFacA :: deep-model grid factor at level center + INTEGER bi, bj, k + _RL deepFacA(Nr) _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL wFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy) _RL vShearTerm(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER myThid -C == Local variables == - INTEGER I,J,Kp1,Km1 - _RL mask_Kp1,mask_Km1,wBarYm,wBarYp - _RL vZm,vZp - LOGICAL rAdvAreaWeight +C !LOCAL VARIABLES: + INTEGER i, j, kp1, km1 + _RL mask_Kp1, mask_Km1, wBarYm, wBarYp + _RL vZm, vZp, recip_drDeepRho + LOGICAL rAdvAreaWeight c _RL vmask_Kp1,vmask_K,vmask_Km1 -c LOGICAL freeslipK,noslipK -c PARAMETER(freeslipK=.TRUE.) -c PARAMETER(noslipK=.NOT.freeslipK) -c LOGICAL freeslip1,noslip1 -c PARAMETER(freeslip1=.TRUE.) -c PARAMETER(noslip1=.NOT.freeslip1) c1 _RL wBarYZ,vZbarZ +CEOP rAdvAreaWeight =.TRUE. C- Area-weighted average either in KE or in vert. advection: IF ( selectKEscheme.EQ.1 .OR. selectKEscheme.EQ.3 ) & rAdvAreaWeight =.FALSE. - Kp1=min(K+1,Nr) - mask_Kp1=1. - IF (K.EQ.Nr) mask_Kp1=0. - Km1=max(K-1,1) - mask_Km1=1. - IF (K.EQ.1) mask_Km1=0. + kp1 = MIN(k+1,Nr) + mask_Kp1 = oneRL + IF (k.EQ.Nr) mask_Kp1 = zeroRL + km1 = MAX(k-1,1) + mask_Km1 = oneRL + IF (k.EQ.1) mask_Km1 = zeroRL + + recip_drDeepRho = recip_drF(k)/deepFacA(k) + & * recip_deepFac2C(k)*recip_rhoFacC(k) - DO J=2-OLy,sNy+OLy - DO I=1-OLx,sNx+OLx + DO j=2-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx c vmask_K=_maskS(i,j,k,bi,bj) @@ -61,43 +67,45 @@ SUBROUTINE MOM_VI_V_VERTSHEAR( c & *mask_Kp1 IF ( rAdvAreaWeight ) THEN -C Transport at interface k : Area weighted average - wBarYm=0.5*( - & wFld(I,J,K,bi,bj)*rA(i,j,bi,bj)*maskC(i,j,Km1,bi,bj) - & +wFld(I,J-1,K,bi,bj)*rA(i,j-1,bi,bj)*maskC(i,j-1,Km1,bi,bj) - & )*mask_Km1*deepFac2F(K)*rhoFacF(K) - & *recip_rAs(i,j,bi,bj) - -C Transport at interface k+1 (here wFld is already masked) - wBarYp=0.5*( - & wFld(I,J,Kp1,bi,bj)*rA(i,j,bi,bj) - & +wFld(I,J-1,Kp1,bi,bj)*rA(i,j-1,bi,bj) - & )*mask_Kp1*deepFac2F(Kp1)*rhoFacF(Kp1) - & *recip_rAs(i,j,bi,bj) +C Transport at interface k : Area weighted average + wBarYm = halfRL*( + & wFld(i,j,k,bi,bj)*rA(i,j,bi,bj)*maskC(i,j,km1,bi,bj) + & + wFld(i,j-1,k,bi,bj)*rA(i,j-1,bi,bj)*maskC(i,j-1,km1,bi,bj) + & )*mask_Km1*deepFac2F(k)*rhoFacF(k) + & *recip_rAs(i,j,bi,bj) + +C Transport at interface k+1 (here wFld is already masked) + wBarYp = halfRL*( + & wFld(i,j,kp1,bi,bj)*rA(i,j,bi,bj) + & + wFld(i,j-1,kp1,bi,bj)*rA(i,j-1,bi,bj) + & )*mask_Kp1*deepFac2F(kp1)*rhoFacF(kp1) + & *recip_rAs(i,j,bi,bj) ELSE -C Transport at interface k : simple average - wBarYm=0.5*( - & wFld(I,J,K,bi,bj)*maskC(i,j,Km1,bi,bj) - & +wFld(I,J-1,K,bi,bj)*maskC(i,j-1,Km1,bi,bj) - & )*mask_Km1*deepFac2F(K)*rhoFacF(K) - -C Transport at interface k+1 (here wFld is already masked) - wBarYp=0.5*( - & wFld(I,J,Kp1,bi,bj) - & +wFld(I,J-1,Kp1,bi,bj) - & )*mask_Kp1*deepFac2F(Kp1)*rhoFacF(Kp1) +C Transport at interface k : simple average + wBarYm = halfRL*( + & wFld(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj) + & + wFld(i,j-1,k,bi,bj)*maskC(i,j-1,km1,bi,bj) + & )*mask_Km1*deepFac2F(k)*rhoFacF(k) + +C Transport at interface k+1 (here wFld is already masked) + wBarYp = halfRL*( + & wFld(i,j,kp1,bi,bj) + & + wFld(i,j-1,kp1,bi,bj) + & )*mask_Kp1*deepFac2F(kp1)*rhoFacF(kp1) ENDIF -C delta_Z( V ) @ interface k +C- delta_Z( V*deepFac ) @ interface k c vmask_Km1=mask_Km1*maskS(i,j,Km1,bi,bj) - vZm=(vFld(I,J,K,bi,bj)-mask_Km1*vFld(I,J,Km1,bi,bj))*rkSign + vZm = ( vFld(i,j, k ,bi,bj)*deepFacA( k ) + & - vFld(i,j,km1,bi,bj)*deepFacA(km1)*mask_Km1 )*rkSign c2 & *recip_dRC(K) c IF (freeslip1) vZm=vZm*vmask_Km1 c IF (noslip1.AND.vmask_Km1.EQ.0.) vZm=vZm*2. -C delta_Z( V ) @ interface k+1 +C- delta_Z( V*deepFac ) @ interface k+1 c vmask_Kp1=mask_Kp1*maskS(i,j,Kp1,bi,bj) - vZp=(mask_Kp1*vFld(I,J,Kp1,bi,bj)-vFld(I,J,K,bi,bj))*rkSign + vZp = ( vFld(i,j,kp1,bi,bj)*deepFacA(kp1)*mask_Kp1 + & - vFld(i,j, k ,bi,bj)*deepFacA( k ) )*rkSign c2 & *recip_dRC(Kp1) c IF (freeslipK) vZp=vZp*vmask_Kp1 c IF (noslipK.AND.vmask_Kp1.EQ.0.) vZp=vZp*2. @@ -117,17 +125,13 @@ SUBROUTINE MOM_VI_V_VERTSHEAR( c2 vShearTerm(I,J)=-0.5*(wBarYp*vZp+wBarYm*vZm) c2 & *_maskS(I,J,K,bi,bj) IF (upwindShear) THEN - vShearTerm(I,J)=-0.5* - & ( (wBarYp*vZp+wBarYm*vZm) - & +(ABS(wBarYp)*vZp-ABS(wBarYm)*vZm) - & )*_recip_hFacS(i,j,k,bi,bj) - & * recip_drF(K) - & * recip_deepFac2C(K)*recip_rhoFacC(K) + vShearTerm(i,j) = -halfRL* + & ( ( wBarYp *vZp + wBarYm *vZm ) + & + ( ABS(wBarYp)*vZp - ABS(wBarYm)*vZm ) + & )*_recip_hFacS(i,j,k,bi,bj)*recip_drDeepRho ELSE - vShearTerm(I,J)=-0.5*(wBarYp*vZp+wBarYm*vZm) - & *_recip_hFacS(i,j,k,bi,bj) - & * recip_drF(K) - & * recip_deepFac2C(K)*recip_rhoFacC(K) + vShearTerm(i,j) = -halfRL*( wBarYp*vZp + wBarYm*vZm ) + & *_recip_hFacS(i,j,k,bi,bj)*recip_drDeepRho ENDIF ENDDO ENDDO