/[MITgcm]/MITgcm_contrib/torge/itd/code/seaice_calc_viscosities.F
ViewVC logotype

Annotation of /MITgcm_contrib/torge/itd/code/seaice_calc_viscosities.F

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


Revision 1.4 - (hide annotations) (download)
Wed Mar 27 18:59:52 2013 UTC (12 years, 4 months ago) by torge
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +16 -5 lines
updating my MITgcm_contrib directory to include latest changes on main branch;
settings are to run a 1D test szenario with ITD code and 7 categories

1 torge 1.4 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_calc_viscosities.F,v 1.22 2013/02/28 16:26:31 mlosch Exp $
2 torge 1.1 C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5    
6 torge 1.3 CBOP
7 torge 1.1 CStartOfInterface
8     SUBROUTINE SEAICE_CALC_VISCOSITIES(
9     I e11, e22, e12, zMin, zMax, hEffM, press0,
10 torge 1.3 O eta, etaZ, zeta, press,
11 torge 1.1 I iStep, myTime, myIter, myThid )
12     C /==========================================================\
13     C | SUBROUTINE SEAICE_CALC_VISCOSITIES |
14     C | o compute shear and bulk viscositites eta, zeta and the |
15     C | corrected ice strength P |
16     C | (see Zhang and Hibler, JGR, 102, 8691-8702, 1997) |
17     C |==========================================================|
18     C | written by Martin Losch, Mar 2006 |
19     C \==========================================================/
20     IMPLICIT NONE
21    
22     C === Global variables ===
23     #include "SIZE.h"
24     #include "EEPARAMS.h"
25     #include "PARAMS.h"
26     #include "GRID.h"
27     #include "SEAICE_SIZE.h"
28     #include "SEAICE_PARAMS.h"
29    
30     #ifdef ALLOW_AUTODIFF_TAMC
31     # include "tamc.h"
32     #endif
33    
34     C === Routine arguments ===
35     C iStep :: Sub-time-step number
36     C myTime :: Simulation time
37     C myIter :: Simulation timestep number
38     C myThid :: My Thread Id. number
39     INTEGER iStep
40     _RL myTime
41     INTEGER myIter
42     INTEGER myThid
43     C strain rate tensor
44     _RL e11 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
45     _RL e22 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
46     _RL e12 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
47     C
48     _RL zMin (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
49     _RL zMax (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
50     _RL hEffM (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
51     C
52     _RL press0(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
53     _RL press (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
54     C bulk viscosity
55     _RL eta (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
56 torge 1.3 _RL etaZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
57 torge 1.1 C shear viscosity
58     _RL zeta (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
59     CEndOfInterface
60    
61     #if ( defined (SEAICE_CGRID) && defined (SEAICE_ALLOW_DYNAMICS) )
62     C === Local variables ===
63     C i,j,bi,bj - Loop counters
64     C e11, e12, e22 - components of strain rate tensor
65     C ecm2 - inverse of square of eccentricity of yield curve
66     INTEGER i, j, bi, bj
67 torge 1.3 _RL ECM2, deltaC, deltaCreg, tmp
68     _RL e12Csqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
69     _RL e11Zsqr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70     _RL e22Zsqr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71     _RL e11e22Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72     _RL pressZ (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
73 torge 1.1 #ifdef SEAICE_ALLOW_TEM
74     _RL etaMax, etaDen
75     #endif /* SEAICE_ALLOW_TEM */
76 torge 1.3 INTEGER k
77     _RL zetaZ, deltaZ, deltaZreg, sumNorm, zetaZmax
78 torge 1.2 #ifdef SEAICE_ZETA_SMOOTHREG
79     _RL argTmp
80     #endif /* SEAICE_ZETA_SMOOTHREG */
81 torge 1.3 CEOP
82 torge 1.1
83     C-- FIRST SET UP BASIC CONSTANTS
84 torge 1.3 k=1
85 torge 1.1 ecm2=0. _d 0
86     IF ( SEAICE_eccen .NE. 0. _d 0 ) ecm2=ONE/(SEAICE_eccen**2)
87     C
88     DO bj=myByLo(myThid),myByHi(myThid)
89     DO bi=myBxLo(myThid),myBxHi(myThid)
90     C need to do this beforehand for easier vectorization after
91     C TAFization
92 torge 1.3 IF ( SEAICEetaZmethod .LT. 2 ) THEN
93     DO j=1-Oly+1,sNy+Oly-1
94     DO i=1-Olx+1,sNx+Olx-1
95     tmp = 0.25 *
96     & ( e12(I,J ,bi,bj) + e12(I+1,J ,bi,bj)
97     & + e12(I,J+1,bi,bj) + e12(I+1,J+1,bi,bj) )
98     e12Csqr(i,j) = tmp*tmp
99     ENDDO
100     ENDDO
101 torge 1.4 ELSEIF ( SEAICEetaZmethod .GE. 2 ) THEN
102 torge 1.3 DO j=1-Oly+1,sNy+Oly-1
103     DO i=1-Olx+1,sNx+Olx-1
104 torge 1.4 e12Csqr(i,j) = 0.25 _d 0 *
105 torge 1.3 & ( e12(I,J ,bi,bj)**2 + e12(I+1,J ,bi,bj)**2
106     & + e12(I,J+1,bi,bj)**2 + e12(I+1,J+1,bi,bj)**2 )
107     ENDDO
108 torge 1.1 ENDDO
109 torge 1.3 ENDIF
110 torge 1.1 DO j=1-Oly+1,sNy+Oly-1
111     DO i=1-Olx+1,sNx+Olx-1
112     deltaC = (e11(i,j,bi,bj)**2+e22(i,j,bi,bj)**2)*(ONE+ecm2)
113 torge 1.3 & + 4. _d 0*ecm2*e12Csqr(i,j)
114 torge 1.1 & + 2. _d 0*e11(i,j,bi,bj)*e22(i,j,bi,bj)*(ONE-ecm2)
115     #ifdef ALLOW_AUTODIFF_TAMC
116     C avoid sqrt of 0
117     IF ( deltaC .GT. 0. _d 0 ) deltaC = SQRT(deltaC)
118     #else
119     deltaC = SQRT(deltaC)
120     #endif /* ALLOW_AUTODIFF_TAMC */
121     deltaCreg = MAX(deltaC,SEAICE_EPS)
122     C "replacement pressure"
123     zeta (I,J,bi,bj) = HALF*press0(I,J,bi,bj)/deltaCreg
124     C put min and max viscosities in
125     #ifdef SEAICE_ZETA_SMOOTHREG
126     C regularize zeta to zmax with a smooth tanh-function instead
127     C of a min(zeta,zmax). This improves convergence of iterative
128     C solvers (Lemieux and Tremblay 2009, JGR). No effect on EVP
129 torge 1.3 argTmp = exp(-1. _d 0/(deltaCreg*SEAICE_zetaMaxFac))
130     zeta (I,J,bi,bj) = ZMAX(I,J,bi,bj)
131 torge 1.2 & *(1. _d 0 - argTmp)/(1. _d 0 + argTmp)
132 torge 1.1 #else
133     zeta (I,J,bi,bj) = MIN(ZMAX(I,J,bi,bj),zeta(I,J,bi,bj))
134     #endif /* SEAICE_ZETA_SMOOTHREG */
135 torge 1.2 zeta (I,J,bi,bj) = MAX(ZMIN(I,J,bi,bj),zeta(I,J,bi,bj))
136 torge 1.1 C set viscosities to zero at hEffM flow pts
137     zeta (I,J,bi,bj) = zeta(I,J,bi,bj)*HEFFM(I,J,bi,bj)
138     eta (I,J,bi,bj) = ECM2*zeta(I,J,bi,bj)
139     press(I,J,bi,bj) = TWO *zeta(I,J,bi,bj)*deltaC
140     ENDDO
141     ENDDO
142     #ifdef SEAICE_ALLOW_TEM
143     IF ( SEAICEuseTEM ) THEN
144     DO j=1-Oly+1,sNy+Oly-1
145     DO i=1-Olx+1,sNx+Olx-1
146     etaDen = (e11(I,J,bi,bj)-e22(I,J,bi,bj))**2
147 torge 1.3 & + 4. _d 0*e12Csqr(i,j)
148 torge 1.1 etaDen = SQRT(MAX(SEAICE_EPS_SQ,etaDen))
149     etaMax = ( 0.5 _d 0*press(I,J,bi,bj)-zeta(I,J,bi,bj)
150     & *( e11(I,J,bi,bj)+e22(I,J,bi,bj) )
151     & )/etaDen
152     eta(I,J,bi,bj) = MIN(eta(I,J,bi,bj),etaMax)
153     ENDDO
154     ENDDO
155     ENDIF
156     #endif /* SEAICE_ALLOW_TEM */
157 torge 1.3 C compute eta at Z-points
158     IF ( SEAICEetaZmethod .eq. 0 ) THEN
159     C This is the old and stupid way
160     DO j=1-Oly+1,sNy+Oly-1
161     DO i=1-Olx+1,sNx+Olx-1
162     etaZ(I,J,bi,bj) =
163     & ( eta (I,J ,bi,bj) + eta (I-1,J ,bi,bj)
164     & + eta (I,J-1,bi,bj) + eta (I-1,J-1,bi,bj) )
165     & / MAX(1.D0,maskC(I,J, k,bi,bj)+maskC(I-1,J, k,bi,bj)
166     & + maskC(I,J-1,k,bi,bj)+maskC(I-1,J-1,k,bi,bj) )
167     ENDDO
168     ENDDO
169 torge 1.4 ELSEIF ( SEAICEetaZmethod .GT. 2 ) THEN
170     DO j=1-Oly+1,sNy+Oly-1
171     DO i=1-Olx+1,sNx+Olx-1
172     sumNorm = maskC(I,J, k,bi,bj)+maskC(I-1,J, k,bi,bj)
173     & + maskC(I,J-1,k,bi,bj)+maskC(I-1,J-1,k,bi,bj)
174     IF ( sumNorm.GT.0. _d 0 ) sumNorm = 1. _d 0 / sumNorm
175     etaZ(I,J,bi,bj) = sumNorm *
176     & ( eta (I,J ,bi,bj) + eta (I-1,J ,bi,bj)
177     & + eta (I,J-1,bi,bj) + eta (I-1,J-1,bi,bj) )
178     ENDDO
179     ENDDO
180 torge 1.3 ELSE
181     IF ( SEAICEetaZmethod .eq. 1 ) THEN
182     DO j=1-Oly+1,sNy+Oly-1
183     DO i=1-Olx+1,sNx+Olx-1
184     sumNorm = maskC(I,J, k,bi,bj)+maskC(I-1,J, k,bi,bj)
185     & + maskC(I,J-1,k,bi,bj)+maskC(I-1,J-1,k,bi,bj)
186     IF ( sumNorm.GT.0. _d 0 ) sumNorm = 1. _d 0 / sumNorm
187     pressZ(i,j) = sumNorm *
188     & ( press0(i,j, bi,bj) + press0(i-1,j, bi,bj)
189     & + press0(i,j-1,bi,bj) + press0(i-1,j-1,bi,bj) )
190     e11Zsqr(i,j) = ( sumNorm *
191     & ( e11(i,j, bi,bj) + e11(i-1,j, bi,bj)
192     & + e11(i,j-1,bi,bj) + e11(i-1,j-1,bi,bj) ) )**2
193     e22Zsqr(i,j) = ( sumNorm *
194     & ( e22(i,j, bi,bj) + e22(i-1,j, bi,bj)
195     & + e22(i,j-1,bi,bj) + e22(i-1,j-1,bi,bj) ) )**2
196     e11e22Z(i,j) = sumNorm * sumNorm *
197     & ( e11(i, j, bi,bj) + e11(i-1,j, bi,bj)
198     & + e11(i, j-1,bi,bj) + e11(i-1,j-1,bi,bj) ) *
199     & ( e22(i, j, bi,bj) + e22(i-1,j, bi,bj)
200     & + e22(i ,j-1,bi,bj) + e22(i-1,j-1,bi,bj) )
201     ENDDO
202     ENDDO
203     ELSEIF ( SEAICEetaZmethod .eq. 2 ) THEN
204     DO j=1-Oly+1,sNy+Oly-1
205     DO i=1-Olx+1,sNx+Olx-1
206     sumNorm = maskC(I,J, k,bi,bj)+maskC(I-1,J, k,bi,bj)
207     & + maskC(I,J-1,k,bi,bj)+maskC(I-1,J-1,k,bi,bj)
208     IF ( sumNorm.GT.0. _d 0 ) sumNorm = 1. _d 0 / sumNorm
209     pressZ(i,j) = sumNorm *
210     & ( press0(i,j, bi,bj) + press0(i-1,j, bi,bj)
211     & + press0(i,j-1,bi,bj) + press0(i-1,j-1,bi,bj) )
212     e11Zsqr(i,j) = sumNorm *
213     & ( e11(i,j, bi,bj)**2 + e11(i-1,j, bi,bj)**2
214     & + e11(i,j-1,bi,bj)**2 + e11(i-1,j-1,bi,bj)**2 )
215     e22Zsqr(i,j) = sumNorm *
216     & ( e22(i,j, bi,bj)**2 + e22(i-1,j, bi,bj)**2
217     & + e22(i,j-1,bi,bj)**2 + e22(i-1,j-1,bi,bj)**2 )
218     e11e22Z(i,j) = sumNorm *
219     & ( e11(i, j, bi,bj)*e22(i, j, bi,bj)
220     & + e11(i-1,j, bi,bj)*e22(i-1,j, bi,bj)
221     & + e11(i, j-1,bi,bj)*e22(i ,j-1,bi,bj)
222     & + e11(i-1,j-1,bi,bj)*e22(i-1,j-1,bi,bj) )
223     ENDDO
224     ENDDO
225     ENDIF
226     DO j=1-Oly+1,sNy+Oly-1
227     DO i=1-Olx+1,sNx+Olx-1
228     deltaZ = (e11Zsqr(i,j)+e22Zsqr(i,j))*(ONE+ecm2)
229     & + 4. _d 0*ecm2*e12(i,j,bi,bj)**2
230     & + 2. _d 0*e11e22Z(i,j)*(ONE-ecm2)
231     #ifdef ALLOW_AUTODIFF_TAMC
232     C avoid sqrt of 0
233     IF ( deltaZ .GT. 0. _d 0 ) deltaZ = SQRT(deltaZ)
234     #else
235     deltaZ = SQRT(deltaZ)
236     #endif /* ALLOW_AUTODIFF_TAMC */
237     deltaZreg = MAX(deltaZ,SEAICE_EPS)
238     C "replacement pressure"
239     zetaZ = HALF*pressZ(i,j)/deltaZreg
240     C put min and max viscosities in
241     zetaZmax = SEAICE_zetaMaxFac*pressZ(i,j)
242     #ifdef SEAICE_ZETA_SMOOTHREG
243     C regularize zeta to zmax with a smooth tanh-function instead
244     C of a min(zeta,zmax). This improves convergence of iterative
245     C solvers (Lemieux and Tremblay 2009, JGR). No effect on EVP
246     argTmp = exp(-1. _d 0/(deltaZreg*SEAICE_zetaMaxFac))
247     zetaZ = zetaZmax*(1. _d 0-argTmp)/(1. _d 0+argTmp)
248     #else
249     zetaZ = MIN(zetaZmax,zetaZ)
250     #endif /* SEAICE_ZETA_SMOOTHREG */
251     zetaZ = MAX(ZMIN(I,J,bi,bj),zetaZ)
252     etaZ(I,J,bi,bj) = ECM2*zetaZ
253     ENDDO
254     ENDDO
255     #ifdef SEAICE_ALLOW_TEM
256     IF ( SEAICEuseTEM ) THEN
257     STOP 'OOPS'
258     CML DO j=1-Oly+1,sNy+Oly-1
259     CML DO i=1-Olx+1,sNx+Olx-1
260     CML etaDen = (e11Zsqr(i,j)+e22Zsqr(i,j)-2.*e11e22Z(i,j))
261     CML & + 4. _d 0*e12(I,J,bi,bj)**2
262     CML etaDen = SQRT(MAX(SEAICE_EPS_SQ,etaDen))
263     CML etaMax = ( 0.5 _d 0*press(I,J,bi,bj)-zeta(I,J,bi,bj)
264     CML & *( e11(I,J,bi,bj)+e22(I,J,bi,bj) )
265     CML & )/etaDen
266     CML etaZ(I,J,bi,bj) = MIN(etaZ(I,J,bi,bj),etaMax)
267     CML ENDDO
268     CML ENDDO
269     ENDIF
270     #endif /* SEAICE_ALLOW_TEM */
271     C SEAICEetaZmethod
272     ENDIF
273     C free-slip means no lateral stress, which is best achieved masking
274     C eta on vorticity(=Z)-points; from now on we only need to worry
275     C about the no-slip boundary conditions
276     IF (.NOT.SEAICE_no_slip) THEN
277 torge 1.4 DO J=1-Oly+1,sNy+Oly-1
278     DO I=1-Olx+1,sNx+Olx-1
279 torge 1.3 etaZ(I,J,bi,bj) = etaZ(I,J,bi,bj)
280     & *maskC(I,J, k,bi,bj)*maskC(I-1,J, k,bi,bj)
281     & *maskC(I,J-1,k,bi,bj)*maskC(I-1,J-1,k,bi,bj)
282     ENDDO
283     ENDDO
284     ENDIF
285 torge 1.1 ENDDO
286     ENDDO
287    
288     #endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID */
289     RETURN
290     END

  ViewVC Help
Powered by ViewVC 1.1.22