/[MITgcm]/MITgcm/pkg/mom_common/mom_calc_smag_3d.F
ViewVC logotype

Contents of /MITgcm/pkg/mom_common/mom_calc_smag_3d.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Tue Nov 5 13:31:50 2013 UTC (10 years, 7 months ago) by jmc
Branch: MAIN
- move here (previously in tutorial_deep_convection code) 2nd version of isotropic
  3-D Smagorinsky code interface: strain and viscosity are locally declared in
  dynmics.F and pass as argument to CALC_GW; ensure that all field value that are
  used are set.

1 C $Header: /u/gcmpack/MITgcm/verification/tutorial_deep_convection/code/mom_calc_smag_3d.F,v 1.2 2013/11/04 21:33:35 jmc Exp $
2 C $Name: $
3
4 #include "MOM_COMMON_OPTIONS.h"
5
6 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7 CBOP
8 C !ROUTINE: MOM_CALC_SMAG_3D
9
10 C !INTERFACE:
11 SUBROUTINE MOM_CALC_SMAG_3D(
12 I str11, str22, str33, str12, str13, str23,
13 O viscAh3d_00, viscAh3d_12,
14 O viscAh3d_13, viscAh3d_23,
15 I smag3D_hLsC, smag3D_hLsW, smag3D_hLsS, smag3D_hLsZ,
16 I k, bi, bj, myThid )
17
18 C !DESCRIPTION:
19 C Calculate Smagorinsky 3-D (harmonic) viscosities
20 C at current level k (for viscAh3d_00 & viscAh3d_12)
21 C and at level k+1 (for viscAh3d_13 & viscAh3d_23)
22
23 C !USES:
24 IMPLICIT NONE
25
26 C == Global variables ==
27 #include "SIZE.h"
28 #include "EEPARAMS.h"
29 #include "PARAMS.h"
30 #include "GRID.h"
31 c#include "MOM_VISC.h"
32
33 C !INPUT PARAMETERS:
34 C str11 :: strain component Vxx @ grid-cell center
35 C str22 :: strain component Vyy @ grid-cell center
36 C str33 :: strain component Vzz @ grid-cell center
37 C str12 :: strain component Vxy @ grid-cell corner
38 C str13 :: strain component Vxz @ above uVel
39 C str23 :: strain component Vyz @ above vVel
40 C smag3D_hLsC :: horiz. grid length scale (power 2/3) at grid cell center
41 C smag3D_hLsW :: horiz. grid length scale (power 2/3) at western edge
42 C smag3D_hLsS :: horiz. grid length scale (power 2/3) at southern egde
43 C smag3D_hLsZ :: horiz. grid length scale (power 2/3) at grid cell corner
44 C k :: current level index
45 C bi, bj :: tile indices
46 C myThid :: my Thread Id number
47 _RL str11(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
48 _RL str22(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
49 _RL str33(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
50 _RL str12(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
51 _RL str13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
52 _RL str23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
53 _RS smag3D_hLsC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54 _RS smag3D_hLsW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
55 _RS smag3D_hLsS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
56 _RS smag3D_hLsZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
57 INTEGER k, bi, bj
58 INTEGER myThid
59
60 C !OUTPUT PARAMETERS:
61 C viscAh3d_00 :: Smagorinsky viscosity @ grid-cell center
62 C viscAh3d_12 :: Smagorinsky viscosity @ grid-cell corner
63 C viscAh3d_13 :: Smagorinsky viscosity @ above uVel
64 C viscAh3d_23 :: Smagorinsky viscosity @ above vVel
65 _RL viscAh3d_00(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
66 _RL viscAh3d_12(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
67 _RL viscAh3d_13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
68 _RL viscAh3d_23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
69 CEOP
70
71 C !LOCAL VARIABLES:
72 INTEGER i, j
73 INTEGER kp, n
74 _RL twoThird, tmpFac
75 _RL S11(1-OLx:sNx+OLx,1-OLy:sNy+OLy,3)
76 _RL S22(1-OLx:sNx+OLx,1-OLy:sNy+OLy,3)
77 _RL S33(1-OLx:sNx+OLx,1-OLy:sNy+OLy,3)
78 _RL S12(1-OLx:sNx+OLx,1-OLy:sNy+OLy,3)
79 _RL S13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,3)
80 _RL S23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,3)
81
82 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
83
84 twoThird = 2. _d 0 / 3. _d 0
85
86 n = 1
87 DO kp=k-1,k+1
88 IF ( kp.GE.1 .AND. kp.LE.Nr ) THEN
89 DO j=1-OLy,sNy+OLy
90 DO i=1-OLx,sNx+OLx
91 S11(i,j,n) = str11(i,j,kp)**2
92 S22(i,j,n) = str22(i,j,kp)**2
93 S33(i,j,n) = str33(i,j,kp)**2
94 S12(i,j,n) = str12(i,j,kp)**2
95 S13(i,j,n) = str13(i,j,kp)**2
96 S23(i,j,n) = str23(i,j,kp)**2
97 ENDDO
98 ENDDO
99 ELSEIF ( kp.GE.1 ) THEN
100 DO j=1-OLy,sNy+OLy
101 DO i=1-OLx,sNx+OLx
102 S11(i,j,n) = 0. _d 0
103 S22(i,j,n) = 0. _d 0
104 S33(i,j,n) = 0. _d 0
105 S12(i,j,n) = 0. _d 0
106 S13(i,j,n) = str13(i,j,kp)**2
107 S23(i,j,n) = str23(i,j,kp)**2
108 ENDDO
109 ENDDO
110 ELSE
111 DO j=1-OLy,sNy+OLy
112 DO i=1-OLx,sNx+OLx
113 S11(i,j,n) = 0. _d 0
114 S22(i,j,n) = 0. _d 0
115 S33(i,j,n) = 0. _d 0
116 S12(i,j,n) = 0. _d 0
117 S13(i,j,n) = 0. _d 0
118 S23(i,j,n) = 0. _d 0
119 ENDDO
120 ENDDO
121 ENDIF
122 n=n+1
123 ENDDO
124
125 C-- ------------------------------------------------------------------
126 C-- calculate current level Smag viscosity coeff
127 C-- ------------------------------------------------------------------
128
129 C Current level k --> n=2
130 n = 2
131 tmpFac = smag3D_coeff*twoRL*SQRT(twoRL)*drF(k)**twoThird
132
133 C viscAh3d_00 = sqrt( S11+S22+S33+2*(S12+S13+S23) ) @ grid-cell center
134
135 DO j=1-OLy,sNy+OLy-1
136 DO i=1-OLx,sNx+OLx-1
137 viscAh3d_00(i,j,k) =
138 & smag3D_hLsC(i,j,bi,bj)*tmpFac*SQRT(
139 & S11( i , j , n )
140 & + S22( i , j , n )
141 & + S33( i , j , n )
142 & + 0.5*( S12( i ,j+1, n )+S12(i+1,j+1, n )
143 & +S12( i , j , n )+S12(i+1, j , n ) )
144 & + 0.5*( S13( i , j , n )+S13(i+1, j , n )
145 & +S13( i , j ,n+1)+S13(i+1, j ,n+1) )
146 & + 0.5*( S23( i , j , n )+S23( i ,j+1, n )
147 & +S23( i , j ,n+1)+S23( i ,j+1,n+1) )
148 & )
149 ENDDO
150 ENDDO
151
152 C viscAh3d_12 = sqrt( S11+S22+S33+2*(S12+S13+S23) ) @ grid-cell corner
153
154 DO j=2-OLy,sNy+OLy
155 DO i=2-OLx,sNx+OLx
156 viscAh3d_12(i,j,k) =
157 & smag3D_hLsZ(i,j,bi,bj)*tmpFac*SQRT(
158 & 0.25*( S11(i-1, j , n )+S11( i , j , n )
159 & +S11(i-1,j-1, n )+S11( i ,j-1, n ) )
160 & + 0.25*( S22(i-1, j , n )+S22( i , j , n )
161 & +S22(i-1,j-1, n )+S22( i ,j-1, n ) )
162 & + 0.25*( S33(i-1, j , n )+S33( i , j , n )
163 & +S33(i-1,j-1, n )+S33( i ,j-1, n ) )
164 & + 2.0 * S12( i , j , n )
165 & + 0.5 *( S13( i ,j-1, n )+S13( i , j , n )
166 & +S13( i ,j-1,n+1)+S13( i , j ,n+1) )
167 & + 0.5 *( S23(i-1, j , n )+S23( i , j , n )
168 & +S23(i-1, j ,n+1)+S23( i , j ,n+1) )
169 & )
170 ENDDO
171 ENDDO
172
173 C-- ------------------------------------------------------------------
174 C-- calculate next level (k+1) viscosity coeff (uz,vz)
175 C-- ------------------------------------------------------------------
176
177 IF ( k.EQ.1 ) THEN
178 DO j=1-OLy,sNy+OLy
179 DO i=1-OLx,sNx+OLx
180 viscAh3d_13(i,j,k) = 0. _d 0
181 viscAh3d_23(i,j,k) = 0. _d 0
182 ENDDO
183 ENDDO
184 ENDIF
185
186 C Next level k+1 --> n=3
187 n = 3
188 tmpFac = smag3D_coeff*twoRL*SQRT(twoRL)*drC(k+1)**twoThird
189
190 C viscAh3d_13 = sqrt( S11+S22+S33+2*(S12+S13+S23) ) @ above uVel
191
192 DO j=1-OLy,sNy+OLy-1
193 DO i=2-OLx,sNx+OLx
194 viscAh3d_13(i,j,k+1) =
195 & smag3D_hLsW(i,j,bi,bj)*tmpFac*SQRT(
196 & 0.25*( S11(i-1, j ,n-1)+S11( i , j ,n-1)
197 & +S11(i-1, j , n )+S11( i , j , n ) )
198 & + 0.25*( S22(i-1, j ,n-1)+S22( i , j ,n-1)
199 & +S22(i-1, j , n )+S22( i , j , n ) )
200 & + 0.25*( S33(i-1, j ,n-1)+S33( i , j ,n-1)
201 & +S33(i-1, j , n )+S33( i , j , n ) )
202 & + 0.5 *( S12( i , j ,n-1)+S12( i ,j+1,n-1)
203 & +S12( i , j , n )+S12( i ,j+1, n ) )
204 & + 2.0 * S13( i , j , n )
205 & + 0.5 *( S23(i-1,j+1, n )+S23( i ,j+1, n )
206 & +S23(i-1, j , n )+S23( i , j , n ) )
207 & )
208 ENDDO
209 ENDDO
210
211 C viscAh3d_23 = sqrt( S11+S22+S33+2*(S12+S13+S23) ) @ above vVel
212
213 DO j=2-OLy,sNy+OLy
214 DO i=1-OLx,sNx+OLx-1
215 viscAh3d_23(i,j,k+1) =
216 & smag3D_hLsS(i,j,bi,bj)*tmpFac*SQRT(
217 & 0.25*( S11( i ,j-1,n-1)+S11( i , j ,n-1)
218 & +S11( i ,j-1, n )+S11( i , j , n ) )
219 & + 0.25*( S22( i ,j-1,n-1)+S22( i , j ,n-1)
220 & +S22( i ,j-1, n )+S22( i , j , n ) )
221 & + 0.25*( S33( i ,j-1,n-1)+S33( i , j ,n-1)
222 & +S33( i ,j-1, n )+S33( i , j , n ) )
223 & + 0.5 *( S12( i , j ,n-1)+S12(i+1, j ,n-1)
224 & +S12( i , j , n )+S12(i+1, j , n ) )
225 & + 0.5 *( S13( i , j , n )+S13(i+1, j , n )
226 & +S13( i ,j-1, n )+S13(i+1,j-1, n ) )
227 & + 2.0 * S23( i , j , n )
228 & )
229 ENDDO
230 ENDDO
231
232 #ifdef ALLOW_DIAGNOSTICS
233 c IF (useDiagnostics) THEN
234 c CALL DIAGNOSTICS_FILL(viscAh_D,'VISCAHD ',k,1,2,bi,bj,myThid)
235 c ENDIF
236 #endif
237
238 RETURN
239 END

  ViewVC Help
Powered by ViewVC 1.1.22