Back to home page

MITgcm
 
 

    

File indexing completed on 2019-04-27 05:10:53 UTC

view on githubraw file Latest commit f59d76b0 on 2019-04-27 00:31:52 UTC
f59d76b0dd Ed D*0001 #include "MOM_COMMON_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C !ROUTINE: MOM_VISC_QGL_STRETCH
                0005 
                0006 C !INTERFACE: ==========================================================
                0007       SUBROUTINE MOM_VISC_QGL_STRETCH(
                0008      I                                 bi, bj, k,
                0009      O                                 stretching, Nsquare,
                0010      I                                 myTime, myIter, myThid )
                0011 
                0012 C     !DESCRIPTION:
                0013 C     *==========================================================*
                0014 C     | SUBROUTINE MOM_VISC_QGL_STRETCH
                0015 C     | Calculates the stratification and vortex stretching terms
                0016 C     | for the Quasi-Geostrophic implementation of Leith
                0017 C     | dynamic viscosity.
                0018 C     | At the upper and lower boundaries, the stratification is
                0019 C     | assumed constant. This means that 'stretching' is
                0020 C     | dictated by the buoyancy difference between the
                0021 C     | boundary cell and its vertical neighbour.
                0022 C     *==========================================================*
                0023 
                0024 C !USES: ===============================================================
                0025       IMPLICIT NONE
                0026 #include "SIZE.h"
                0027 #include "EEPARAMS.h"
                0028 #include "PARAMS.h"
                0029 #include "GRID.h"
                0030 #include "DYNVARS.h"
                0031 c#include "MOM_VISC.h"
                0032 
                0033 C !INPUT PARAMETERS: ===================================================
                0034 C  bi,bj                :: tile indices
                0035 C  k                    :: vertical level
                0036 C  myTime               :: current time of simulation ( s )
                0037 C  myIter               :: current iteration number of simulation
                0038 C  myThid               :: my Thread Id number
                0039       INTEGER bi,bj, k
                0040       _RL     myTime
                0041       INTEGER myIter
                0042       INTEGER myThid
                0043 
                0044 C !OUTPUT PARAMETERS: ==================================================
                0045 C  stretching           :: vortex stretching contribution, calculated at
                0046 C                          cell centre, except near top and bottom, when
                0047 C                          it is calculated at cell face, but assumed to
                0048 C                          reside at cell centre.
                0049 C  Nsquare              :: buoyancy frequency, averaged to cell centre,
                0050 C                          except near top and bottom, when constant
                0051 C                          stratification is assumed, and it is copied
                0052 C                          to cell centre.
                0053       _RL stretching(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0054       _RL Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0055 
                0056 #ifdef ALLOW_LEITH_QG
                0057 
                0058 C !LOCAL VARIABLES: ====================================================
                0059       INTEGER i, j
                0060       INTEGER iMin, iMax
                0061       INTEGER jMin, jMax
                0062 
                0063       _RL QGL_epsil
                0064 
                0065       _RL kernel_1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0066       _RL kernel_2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0067 
                0068       _RL buoy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0069       _RL buoy_m1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0070       _RL buoy_p1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0071       _RL buoy_1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0072       _RL buoy_2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0073 
                0074       _RL Nsquarep1
                0075 
                0076 CEOP
                0077 
                0078       DO j=1-OLy,sNy+OLy
                0079        DO i=1-OLx,sNx+OLx
                0080         kernel_1(i,j)   = 0. _d 0
                0081         kernel_2(i,j)   = 0. _d 0
                0082         buoy(i,j)       = 0. _d 0
                0083         buoy_m1(i,j)    = 0. _d 0
                0084         buoy_p1(i,j)    = 0. _d 0
                0085         buoy_1(i,j)     = 0. _d 0
                0086         buoy_2(i,j)     = 0. _d 0
                0087         stretching(i,j) = 0. _d 0
                0088         Nsquare(i,j)    = 0. _d 0
                0089        ENDDO
                0090       ENDDO
                0091 
                0092       QGL_epsil = 1. _d -12
                0093       Nsquarep1 = 0. _d 0
                0094 
                0095 C use sigmaRfield, rhoInSitu from common block DYNVARS.h
                0096 C sigmaRfield is correctly masked, so no need to mask Nsquare
                0097 
                0098       iMin = 1-OLx
                0099       iMax = sNx+OLx
                0100       jMin = 1-OLy
                0101       jMax = sNy+OLy
                0102 
                0103       DO j=jMin,jMax
                0104         DO i=iMin,iMax
                0105           IF (k.GT.kSurfC(i,j,bi,bj) .AND.
                0106      &                 k.LT.kLowC(i,j,bi,bj)) THEN
                0107 C           In the ocean interior. Standard calculation.
                0108 
                0109 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                0110 C              buoyancy
                0111 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                0112             buoy(i,j)    = gravity*gravitySign*recip_rhoConst*
                0113      &                       rhoInSitu(i,j,k,bi,bj)
                0114             buoy_m1(i,j) = gravity*gravitySign*recip_rhoConst*
                0115      &                       rhoInSitu(i,j,k-1,bi,bj)
                0116             buoy_p1(i,j) = gravity*gravitySign*recip_rhoConst*
                0117      &                       rhoInSitu(i,j,k+1,bi,bj)
                0118 
                0119 C           Interpolate buoyancy to upper and lower cell faces.
                0120 C           (same location as Nsquare is calculated)
                0121             buoy_1(i,j)= halfRL * (buoy(i,j) + buoy_m1(i,j))
                0122             buoy_2(i,j)= halfRL * (buoy(i,j) + buoy_p1(i,j))
                0123 
                0124 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                0125 C              (f/N^2) * b
                0126 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                0127             Nsquare(i,j) = gravity*gravitySign*recip_rhoConst
                0128      &                  * sigmaRfield(i,j,k,bi,bj)
                0129             Nsquarep1    = gravity*gravitySign*recip_rhoConst
                0130      &                  * sigmaRfield(i,j,k+1,bi,bj)
                0131 
                0132             kernel_1(i,j) = (fCori(i,j,bi,bj)/
                0133      &                     MAX(Nsquare(i,j),QGL_epsil))*
                0134      &                        buoy_1(i,j)
                0135             kernel_2(i,j)=(fCori(i,j,bi,bj)/
                0136      &                     MAX(Nsquarep1,QGL_epsil))*
                0137      &                        buoy_2(i,j)
                0138 
                0139 C           Average Nsquare to cell centre for use in
                0140 C           MOM_VISC_QGL_LIMIT
                0141             Nsquare(i,j) = halfRL * ( Nsquare(i,j) + Nsquarep1 )
                0142 
                0143 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                0144 C              d/dz of it
                0145 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                0146             stretching(i,j)= maskC(i,j,k,bi,bj)
                0147      &                *recip_drF(k)*rkSign
                0148      &                *(kernel_2(i,j)-kernel_1(i,j))
                0149 
                0150 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                0151           ELSEIF (k.EQ.kSurfC(i,j,bi,bj) .AND.
                0152      &                 k.EQ.kLowC(i,j,bi,bj)) THEN
                0153 C           Ocean only has one level. There is no possibility
                0154 C           for vertical stratification. Fail gracefully.
                0155             stretching(i,j) = 0. _d 0
                0156 
                0157 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                0158           ELSEIF (k.EQ.kSurfC(i,j,bi,bj) .AND.
                0159      &                 k.LT.kLowC(i,j,bi,bj)) THEN
                0160 C           Ocean has at least two levels - currently in the uppermost one.
                0161 C           Use stratification from k+1, and assume it is constant.
                0162 
                0163 C           buoyancy
                0164             buoy(i,j)    = gravity*gravitySign*recip_rhoConst*
                0165      &                       rhoInSitu(i,j,k,bi,bj)
                0166             buoy_p1(i,j) = gravity*gravitySign*recip_rhoConst*
                0167      &                       rhoInSitu(i,j,k+1,bi,bj)
                0168 
                0169 C           Assuming constant stratification, so no need to
                0170 C           interpolate these.
                0171 
                0172 C           (f/N^2) * b
                0173             Nsquarep1    = gravity*gravitySign*recip_rhoConst
                0174      &                  * sigmaRfield(i,j,k+1,bi,bj)
                0175 
                0176             kernel_1(i,j) = (fCori(i,j,bi,bj)/
                0177      &                     MAX(Nsquarep1,QGL_epsil))*
                0178      &                        buoy(i,j)
                0179             kernel_2(i,j)=(fCori(i,j,bi,bj)/
                0180      &                     MAX(Nsquarep1,QGL_epsil))*
                0181      &                        buoy_p1(i,j)
                0182 
                0183 C           Average Nsquare to cell centre for use in
                0184 C           MOM_VISC_QGL_LIMIT
                0185 C           (have assumed constant stratification, so
                0186 C           just assign it)
                0187             Nsquare(i,j) = Nsquarep1
                0188 
                0189 C           d/dz of it
                0190             stretching(i,j)= maskC(i,j,k,bi,bj)
                0191      &                *recip_drC(k+1)*rkSign
                0192      &                *(kernel_2(i,j)-kernel_1(i,j))
                0193 
                0194 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                0195           ELSEIF (k.GT.kSurfC(i,j,bi,bj) .AND.
                0196      &                 k.EQ.kLowC(i,j,bi,bj)) THEN
                0197 C         Ocean has at least two levels - currently in the lowest one.
                0198 C         Use stratification from this level, and assume it is constant.
                0199 
                0200 C           buoyancy
                0201             buoy(i,j)    = gravity*gravitySign*recip_rhoConst*
                0202      &                       rhoInSitu(i,j,k,bi,bj)
                0203             buoy_m1(i,j) = gravity*gravitySign*recip_rhoConst*
                0204      &                       rhoInSitu(i,j,k-1,bi,bj)
                0205 
                0206 C           Assuming constant stratification, so no need to
                0207 C           interpolate these.
                0208 
                0209 C           (f/N^2) * b
                0210             Nsquare(i,j) = gravity*gravitySign*recip_rhoConst
                0211      &                  * sigmaRfield(i,j,k,bi,bj)
                0212 
                0213             kernel_1(i,j) = (fCori(i,j,bi,bj)/
                0214      &                     MAX(Nsquare(i,j),QGL_epsil))*
                0215      &                        buoy_m1(i,j)
                0216             kernel_2(i,j)=(fCori(i,j,bi,bj)/
                0217      &                     MAX(Nsquare(i,j),QGL_epsil))*
                0218      &                        buoy(i,j)
                0219 
                0220 C           Average Nsquare to cell centre for use in
                0221 C           MOM_VISC_QGL_LIMIT
                0222 C           (have assumed constant stratification, so
                0223 C           no need to do anything)
                0224 
                0225 C           d/dz of it
                0226             stretching(i,j)= maskC(i,j,k,bi,bj)
                0227      &                *recip_drC(k)*rkSign
                0228      &                *(kernel_2(i,j)-kernel_1(i,j))
                0229 
                0230 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                0231           ELSE
                0232 C           Not in the ocean (probably in an ice shelf)
                0233 C           Do nothing - stretching should be zero
                0234             stretching(i,j) = 0. _d 0
                0235 
                0236           ENDIF
                0237 
                0238         ENDDO
                0239       ENDDO
                0240 
                0241 #endif /* ALLOW_LEITH_QG */
                0242 
                0243       RETURN
                0244       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!