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
0004
0005
0006
0007 SUBROUTINE MOM_VISC_QGL_STRETCH(
0008 I bi, bj, k,
0009 O stretching, Nsquare,
0010 I myTime, myIter, myThid )
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025 IMPLICIT NONE
0026 #include "SIZE.h"
0027 #include "EEPARAMS.h"
0028 #include "PARAMS.h"
0029 #include "GRID.h"
0030 #include "DYNVARS.h"
0031
0032
0033
0034
0035
0036
0037
0038
0039 INTEGER bi,bj, k
0040 _RL myTime
0041 INTEGER myIter
0042 INTEGER myThid
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
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
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
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
0096
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
0108
0109
0110
0111
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
0120
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
0125
0126
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
0140
0141 Nsquare(i,j) = halfRL * ( Nsquare(i,j) + Nsquarep1 )
0142
0143
0144
0145
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
0151 ELSEIF (k.EQ.kSurfC(i,j,bi,bj) .AND.
0152 & k.EQ.kLowC(i,j,bi,bj)) THEN
0153
0154
0155 stretching(i,j) = 0. _d 0
0156
0157
0158 ELSEIF (k.EQ.kSurfC(i,j,bi,bj) .AND.
0159 & k.LT.kLowC(i,j,bi,bj)) THEN
0160
0161
0162
0163
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
0170
0171
0172
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
0184
0185
0186
0187 Nsquare(i,j) = Nsquarep1
0188
0189
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
0195 ELSEIF (k.GT.kSurfC(i,j,bi,bj) .AND.
0196 & k.EQ.kLowC(i,j,bi,bj)) THEN
0197
0198
0199
0200
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
0207
0208
0209
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
0221
0222
0223
0224
0225
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
0231 ELSE
0232
0233
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