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

Annotation 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.2 - (hide annotations) (download)
Mon Nov 4 21:33:35 2013 UTC (11 years, 3 months ago) by jmc
Branch: MAIN
Changes since 1.1: +24 -16 lines
This is really Louis-Philippe's version (scaling was missing in previous version)

1 jmc 1.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 jmc 1.2 C at current level k (for viscAh3d_00 & viscAh3d_12)
21     C and at level k+1 (for viscAh3d_13 & viscAh3d_23)
22 jmc 1.1
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 jmc 1.2 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 jmc 1.1 _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 jmc 1.2 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
83    
84     twoThird = 2. _d 0 / 3. _d 0
85    
86 jmc 1.1 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 jmc 1.2 tmpFac = smag3D_coeff*twoRL*SQRT(twoRL)*drF(k)**twoThird
132 jmc 1.1
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 jmc 1.2 viscAh3d_00( i , j , k ) =
138     & smag3D_hLsC( i , j , bi , bj )*tmpFac*SQRT(
139 jmc 1.1 & 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 jmc 1.2 & )
149 jmc 1.1 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 jmc 1.2 viscAh3d_12( i , j , k ) =
157     & smag3D_hLsZ( i , j , bi , bj )*tmpFac*SQRT(
158 jmc 1.1 & 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 jmc 1.2 & )
170 jmc 1.1 ENDDO
171     ENDDO
172    
173     C-- ------------------------------------------------------------------
174     C-- calculate next level (k+1) viscosity coeff (uz,vz)
175     C-- ------------------------------------------------------------------
176    
177     C Next level k+1 --> n=3
178     n = 3
179 jmc 1.2 tmpFac = smag3D_coeff*twoRL*SQRT(twoRL)*drC(k+1)**twoThird
180 jmc 1.1
181     C viscAh3d_13 = sqrt( S11+S22+S33+2*(S12+S13+S23) ) @ above uVel
182    
183     DO j=1-OLy,sNy+OLy-1
184     DO i=2-OLx,sNx+OLx
185 jmc 1.2 viscAh3d_13( i , j ,k+1) =
186     & smag3D_hLsW( i , j , bi , bj )*tmpFac*SQRT(
187 jmc 1.1 & 0.25*( S11(i-1, j ,n-1)+S11( i , j ,n-1)
188     & +S11(i-1, j , n )+S11( i , j , n ) )
189     & + 0.25*( S22(i-1, j ,n-1)+S22( i , j ,n-1)
190     & +S22(i-1, j , n )+S22( i , j , n ) )
191     & + 0.25*( S33(i-1, j ,n-1)+S33( i , j ,n-1)
192     & +S33(i-1, j , n )+S33( i , j , n ) )
193     & + 0.5 *( S12( i , j ,n-1)+S12( i ,j+1,n-1)
194     & +S12( i , j , n )+S12( i ,j+1, n ) )
195     & + 2.0 * S13( i , j , n )
196     & + 0.5 *( S23(i-1,j+1, n )+S23( i ,j+1, n )
197     & +S23(i-1, j , n )+S23( i , j , n ) )
198 jmc 1.2 & )
199 jmc 1.1 ENDDO
200     ENDDO
201    
202     C viscAh3d_23 = sqrt( S11+S22+S33+2*(S12+S13+S23) ) @ above vVel
203    
204     DO j=2-OLy,sNy+OLy
205     DO i=1-OLx,sNx+OLx-1
206 jmc 1.2 viscAh3d_23( i , j ,k+1) =
207     & smag3D_hLsS( i , j , bi , bj )*tmpFac*SQRT(
208 jmc 1.1 & 0.25*( S11( i ,j-1,n-1)+S11( i , j ,n-1)
209     & +S11( i ,j-1, n )+S11( i , j , n ) )
210     & + 0.25*( S22( i ,j-1,n-1)+S22( i , j ,n-1)
211     & +S22( i ,j-1, n )+S22( i , j , n ) )
212     & + 0.25*( S33( i ,j-1,n-1)+S33( i , j ,n-1)
213     & +S33( i ,j-1, n )+S33( i , j , n ) )
214     & + 0.5 *( S12( i , j ,n-1)+S12(i+1, j ,n-1)
215     & +S12( i , j , n )+S12(i+1, j , n ) )
216     & + 0.5 *( S13( i , j , n )+S13(i+1, j , n )
217     & +S13( i ,j-1, n )+S13(i+1,j-1, n ) )
218     & + 2.0 * S23( i , j , n )
219 jmc 1.2 & )
220 jmc 1.1 ENDDO
221     ENDDO
222    
223     #ifdef ALLOW_DIAGNOSTICS
224     c IF (useDiagnostics) THEN
225     c CALL DIAGNOSTICS_FILL(viscAh_D,'VISCAHD ',k,1,2,bi,bj,myThid)
226     c ENDIF
227     #endif
228    
229     RETURN
230     END

  ViewVC Help
Powered by ViewVC 1.1.22