pkg/mom_common/mom_quasihydrostatic.F

af53f2701c Jean*0001 #include "MOM_COMMON_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: MOM_QUASIHYDROSTATIC
                0005
                0006 C !INTERFACE: ==========================================================
                0007       SUBROUTINE MOM_QUASIHYDROSTATIC(
ee19c4a296 Jean*0008      I                bi, bj, k,
af53f2701c Jean*0009      I                uFld, vFld,
ee19c4a296 Jean*0010      U                effectiveBuoy,
                0011      I                myTime, myIter, myThid )
af53f2701c Jean*0012
                0013 C !DESCRIPTION:
                0014 C     *==========================================================*
                0015 C     | o SUBROUTINE MOM_QUASIHYDROSTATIC
                0016 C     |   Add Quasi-Hydrostatic Terms to buoyancy
                0017 C     *==========================================================*
                0018
                0019 C !USES: ===============================================================
                0020       IMPLICIT NONE
                0021
                0022 C--   == Global data ==
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
                0026 #include "GRID.h"
ee19c4a296 Jean*0027 #ifdef ALLOW_QHYD_STAGGER_TS
                0028 # include "RESTART.h"
                0029 # include "NH_VARS.h"
                0030 #endif
af53f2701c Jean*0031
                0032 C--   == Routine arguments ==
                0033 C !INPUT VARIABLES: ====================================================
ee19c4a296 Jean*0034 C  bi, bj        :: tile indices
af53f2701c Jean*0035 C  k             :: vertical level
                0036 C  uFld          :: zonal flow
                0037 C  vFld          :: meridional flow
ee19c4a296 Jean*0038 C  myTime        :: current time in simulation
                0039 C  myIter        :: current iteration number
af53f2701c Jean*0040 C  myThid        :: my Thread Id number
ee19c4a296 Jean*0041       INTEGER bi, bj, k
af53f2701c Jean*0042       _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0043       _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
ee19c4a296 Jean*0044       _RL myTime
                0045       INTEGER myIter
af53f2701c Jean*0046       INTEGER myThid
                0047
                0048 C !OUTPUT/MODIFIED VARIABLES: ==========================================
c178b9393a Jean*0049 C  effectiveBuoy :: Density (z-coord) / specific volume (p-coord) anomaly
                0050       _RL effectiveBuoy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
af53f2701c Jean*0051
                0052 C !LOCAL VARIABLES: ====================================================
ee19c4a296 Jean*0053 C  i, j          :: loop indices
                0054 C  gWinBuoy      :: vertical acceleraton term to add to buoyancy
                0055 C  gw_AB         :: tendency increment from Adams-Bashforth
af53f2701c Jean*0056 C  scalingFactor :: scaling factor (from acceleration to density)
43262d8cdd Jean*0057       INTEGER iMin,iMax,jMin,jMax
                0058       PARAMETER( iMin = 0 , iMax = sNx+1 )
                0059       PARAMETER( jMin = 0 , jMax = sNy+1 )
ee19c4a296 Jean*0060       INTEGER i, j
                0061       _RL gWinBuoy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0062 #ifdef ALLOW_QHYD_STAGGER_TS
                0063       _RL gw_AB   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0064 #endif
af53f2701c Jean*0065       _RL scalingFactor
                0066 CEOP
                0067
427e24e121 Jean*0068       IF ( select3dCoriScheme.GE.1 .OR. useNHMTerms ) THEN
ee19c4a296 Jean*0069
                0070         IF ( usingZCoords ) THEN
c178b9393a Jean*0071 C--   Z-coordinate case: Input is density anomaly
                0072
ee19c4a296 Jean*0073           scalingFactor = rhoConst*gravitySign
                0074      &                            *recip_gravity*recip_gravFacC(k)
af53f2701c Jean*0075 C-    to reproduce old (wrong) results:
ee19c4a296 Jean*0076 c         scalingFactor=gravitySign*recip_gravity
c178b9393a Jean*0077
ee19c4a296 Jean*0078         ELSEIF ( fluidIsWater ) THEN
c178b9393a Jean*0079 C--   P-coordinate, oceanic case: Input is specific-volume anomaly
                0080
5b172de0d2 Jean*0081 c         scalingFactor = recip_rhoConst*recip_gravity
                0082 C     better to use rhoRef = rho(tRef,sRef,p) and consistent with
                0083 C-    omega <-> w-velocity conversion:
                0084           scalingFactor = ( oneRL / rhoRef(k) )*recip_gravity
c178b9393a Jean*0085
ee19c4a296 Jean*0086         ELSE
c178b9393a Jean*0087 C--   P-coord., Ideal-Gas case: Input is virtual potential temp. anomaly
                0088 C     (see White & Bromley, QJRMS 1995)
ee19c4a296 Jean*0089           scalingFactor = tRef(k)*recip_gravity
c178b9393a Jean*0090
ee19c4a296 Jean*0091         ENDIF
af53f2701c Jean*0092
ee19c4a296 Jean*0093         DO j=1-OLy,sNy+OLy
                0094          DO i=1-OLx,sNx+OLx
                0095            gWinBuoy(i,j) = 0. _d 0
                0096          ENDDO
9496c6c9ef Jean*0097         ENDDO
af53f2701c Jean*0098
427e24e121 Jean*0099         IF ( select3dCoriScheme.GE.1 ) THEN
ee19c4a296 Jean*0100          DO j=jMin,jMax
                0101           DO i=iMin,iMax
                0102            gWinBuoy(i,j) = fCoriCos(i,j,bi,bj)*
                0103      &       ( angleCosC(i,j,bi,bj)*halfRL
                0104      &             *( uFld(i,j,k,bi,bj) + uFld(i+1,j,k,bi,bj) )
                0105      &        -angleSinC(i,j,bi,bj)*halfRL
                0106      &             *( vFld(i,j,k,bi,bj) + vFld(i,j+1,k,bi,bj) )
                0107      &       )
                0108           ENDDO
                0109          ENDDO
                0110         ENDIF
                0111
                0112         IF ( useNHMTerms ) THEN
                0113          DO j=jMin,jMax
                0114           DO i=iMin,iMax
                0115            gWinBuoy(i,j) = gWinBuoy(i,j)
                0116      &       + (  ( uFld( i ,j,k,bi,bj)*uFld( i ,j,k,bi,bj)
                0117      &            + uFld(i+1,j,k,bi,bj)*uFld(i+1,j,k,bi,bj) )
                0118      &          + ( vFld(i, j ,k,bi,bj)*vFld(i, j ,k,bi,bj)
                0119      &            + vFld(i,j+1,k,bi,bj)*vFld(i,j+1,k,bi,bj) )
                0120      &         )* halfRL*recip_rSphere*recip_deepFacC(k)
                0121           ENDDO
                0122          ENDDO
                0123         ENDIF
                0124
                0125 #ifdef ALLOW_QHYD_STAGGER_TS
                0126         IF ( staggerTimeStep ) THEN
                0127 # ifdef ALLOW_ADAMSBASHFORTH_3
                0128           CALL ADAMS_BASHFORTH3( bi, bj, k, 1,
                0129      U                           gWinBuoy,
                0130      U                           QHydGwNm,
                0131      O                           gw_AB,
                0132      I                           qHydStartAB, myIter, myThid )
                0133 # else /* ALLOW_ADAMSBASHFORTH_3 */
                0134           CALL ADAMS_BASHFORTH2( bi, bj, 1, 1,
                0135      U                           gWinBuoy,
                0136      U                           QHydGwNm(1-OLx,1-OLy,k,bi,bj),
                0137      O                           gw_AB,
                0138      I                           qHydStartAB, myIter, myThid )
                0139 # endif /* ALLOW_ADAMSBASHFORTH_3 */
                0140         ENDIF
                0141 #endif /* ALLOW_QHYD_STAGGER_TS */
                0142
                0143         DO j=jMin,jMax
                0144          DO i=iMin,iMax
                0145            effectiveBuoy(i,j) = effectiveBuoy(i,j)
                0146      &                        + scalingFactor*gWinBuoy(i,j)
                0147          ENDDO
af53f2701c Jean*0148         ENDDO
ee19c4a296 Jean*0149
af53f2701c Jean*0150       ENDIF
                0151
                0152       RETURN
                0153       END