/[MITgcm]/MITgcm/verification/tutorial_deep_convection/code/mom_calc_smag_3d.F
ViewVC logotype

Contents of /MITgcm/verification/tutorial_deep_convection/code/mom_calc_smag_3d.F

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


Revision 1.1 - (show annotations) (download)
Mon Sep 30 18:19:41 2013 UTC (11 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64p
add 1rst version of isotropic 3-D Smagorinsky code (from Louis-Philippe),
  for now all in this exp. code dir.

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

  ViewVC Help
Powered by ViewVC 1.1.22