/[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.3 - (show annotations) (download)
Thu Nov 7 23:09:45 2013 UTC (10 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint65, HEAD
Changes since 1.2: +50 -71 lines
simplify + reduce local storage for squared strain component

1 C $Header: /u/gcmpack/MITgcm/pkg/mom_common/mom_calc_smag_3d.F,v 1.1 2013/11/05 13:31:50 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 #ifdef ALLOW_SMAG_3D
72 C !LOCAL VARIABLES:
73 C S66 :: Sum of squared of the 3 strain component @ grid-cell center
74 C S12 :: Squared of strain component Vxy @ grid-cell corner
75 C S13 :: Squared of strain component Vxz @ above uVel
76 C S23 :: Squared of strain component Vyz @ above vVel
77 INTEGER i, j
78 INTEGER kl, n
79 _RL twoThird, tmpFac
80 _RL S66(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
81 _RL S12(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
82 _RL S13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
83 _RL S23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
84
85 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
86
87 twoThird = 2. _d 0 / 3. _d 0
88
89 DO n=1,2
90 kl = k + n - 1
91 IF ( kl.LE.Nr ) THEN
92 DO j=1-OLy,sNy+OLy
93 DO i=1-OLx,sNx+OLx
94 S66(i,j,n) = str11(i,j,kl)**2
95 & + str22(i,j,kl)**2
96 & + str33(i,j,kl)**2
97 S12(i,j,n) = str12(i,j,kl)**2
98 ENDDO
99 ENDDO
100 ELSE
101 DO j=1-OLy,sNy+OLy
102 DO i=1-OLx,sNx+OLx
103 S66(i,j,n) = 0. _d 0
104 S12(i,j,n) = 0. _d 0
105 ENDDO
106 ENDDO
107 ENDIF
108 DO j=1-OLy,sNy+OLy
109 DO i=1-OLx,sNx+OLx
110 S13(i,j,n) = str13(i,j,kl)**2
111 S23(i,j,n) = str23(i,j,kl)**2
112 ENDDO
113 ENDDO
114 ENDDO
115
116 C-- ------------------------------------------------------------------
117 C-- calculate current level Smag viscosity coeff
118 C-- ------------------------------------------------------------------
119
120 C Current level k --> n=1
121 kl = k
122 n = 1
123 tmpFac = smag3D_coeff*twoRL*SQRT(twoRL)*drF(kl)**twoThird
124
125 C viscAh3d_00 = sqrt( S11+S22+S33+2*(S12+S13+S23) ) @ grid-cell center
126
127 DO j=1-OLy,sNy+OLy-1
128 DO i=1-OLx,sNx+OLx-1
129 viscAh3d_00(i,j,kl) =
130 & smag3D_hLsC(i,j,bi,bj)*tmpFac*SQRT(
131 & S66( i , j , n )
132 & + 0.5*( S12( i ,j+1, n )+S12(i+1,j+1, n )
133 & +S12( i , j , n )+S12(i+1, j , n ) )
134 & + 0.5*( S13( i , j , n )+S13(i+1, j , n )
135 & +S13( i , j ,n+1)+S13(i+1, j ,n+1) )
136 & + 0.5*( S23( i , j , n )+S23( i ,j+1, n )
137 & +S23( i , j ,n+1)+S23( i ,j+1,n+1) )
138 & )
139 ENDDO
140 ENDDO
141
142 C viscAh3d_12 = sqrt( S11+S22+S33+2*(S12+S13+S23) ) @ grid-cell corner
143
144 DO j=2-OLy,sNy+OLy
145 DO i=2-OLx,sNx+OLx
146 viscAh3d_12(i,j,kl) =
147 & smag3D_hLsZ(i,j,bi,bj)*tmpFac*SQRT(
148 & 0.25*( S66(i-1, j , n )+S66( i , j , n )
149 & +S66(i-1,j-1, n )+S66( i ,j-1, n ) )
150 & + 2.0 * S12( i , j , n )
151 & + 0.5 *( S13( i ,j-1, n )+S13( i , j , n )
152 & +S13( i ,j-1,n+1)+S13( i , j ,n+1) )
153 & + 0.5 *( S23(i-1, j , n )+S23( i , j , n )
154 & +S23(i-1, j ,n+1)+S23( i , j ,n+1) )
155 & )
156 ENDDO
157 ENDDO
158
159 C-- ------------------------------------------------------------------
160 C-- calculate next level (k+1) viscosity coeff (uz,vz)
161 C-- ------------------------------------------------------------------
162
163 IF ( k.EQ.1 ) THEN
164 DO j=1-OLy,sNy+OLy
165 DO i=1-OLx,sNx+OLx
166 viscAh3d_13(i,j,k) = 0. _d 0
167 viscAh3d_23(i,j,k) = 0. _d 0
168 ENDDO
169 ENDDO
170 ENDIF
171
172 C Next level k+1 --> n=2
173 kl = k+1
174 n = 2
175 tmpFac = smag3D_coeff*twoRL*SQRT(twoRL)*drC(kl)**twoThird
176
177 C viscAh3d_13 = sqrt( S11+S22+S33+2*(S12+S13+S23) ) @ above uVel
178
179 DO j=1-OLy,sNy+OLy-1
180 DO i=2-OLx,sNx+OLx
181 viscAh3d_13(i,j,kl) =
182 & smag3D_hLsW(i,j,bi,bj)*tmpFac*SQRT(
183 & 0.25*( S66(i-1, j ,n-1)+S66( i , j ,n-1)
184 & +S66(i-1, j , n )+S66( i , j , n ) )
185 & + 0.5 *( S12( i , j ,n-1)+S12( i ,j+1,n-1)
186 & +S12( i , j , n )+S12( i ,j+1, n ) )
187 & + 2.0 * S13( i , j , n )
188 & + 0.5 *( S23(i-1,j+1, n )+S23( i ,j+1, n )
189 & +S23(i-1, j , n )+S23( i , j , n ) )
190 & )
191 ENDDO
192 ENDDO
193
194 C viscAh3d_23 = sqrt( S11+S22+S33+2*(S12+S13+S23) ) @ above vVel
195
196 DO j=2-OLy,sNy+OLy
197 DO i=1-OLx,sNx+OLx-1
198 viscAh3d_23(i,j,kl) =
199 & smag3D_hLsS(i,j,bi,bj)*tmpFac*SQRT(
200 & 0.25*( S66( i ,j-1,n-1)+S66( i , j ,n-1)
201 & +S66( i ,j-1, n )+S66( i , j , n ) )
202 & + 0.5 *( S12( i , j ,n-1)+S12(i+1, j ,n-1)
203 & +S12( i , j , n )+S12(i+1, j , n ) )
204 & + 0.5 *( S13( i , j , n )+S13(i+1, j , n )
205 & +S13( i ,j-1, n )+S13(i+1,j-1, n ) )
206 & + 2.0 * S23( i , j , n )
207 & )
208 ENDDO
209 ENDDO
210
211 #ifdef ALLOW_DIAGNOSTICS
212 c IF ( useDiagnostics.AND. k.EQ.Nr ) THEN
213 c CALL DIAGNOSTICS_FILL( viscAh3d_00, 'Smag3D_C',
214 c & 0, Nr, 2, bi, bj, myThid )
215 c ENDIF
216 #endif
217
218 #endif /* ALLOW_SMAG_3D */
219 RETURN
220 END

  ViewVC Help
Powered by ViewVC 1.1.22