/[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.1 - (hide 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 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    
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