/[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.2 - (hide annotations) (download)
Wed Oct 31 16:34:43 2012 UTC (12 years, 9 months ago) by torge
Branch: MAIN
Changes since 1.1: +12 -6 lines 
merge with changes on main branch
1 torge 1.2 C $Header: /u/gcmpack/MITgcm_contrib/torge/itd/code/seaice_calc_viscosities.F,v 1.1 2012/10/24 21:48:53 torge Exp $
2 torge 1.1 C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5    
6     CStartOfInterface
7     SUBROUTINE SEAICE_CALC_VISCOSITIES(
8     I e11, e22, e12, zMin, zMax, hEffM, press0,
9     O eta, zeta, press,
10     I iStep, myTime, myIter, myThid )
11     C /==========================================================\
12     C | SUBROUTINE SEAICE_CALC_VISCOSITIES |
13     C | o compute shear and bulk viscositites eta, zeta and the |
14     C | corrected ice strength P |
15     C | (see Zhang and Hibler, JGR, 102, 8691-8702, 1997) |
16     C |==========================================================|
17     C | written by Martin Losch, Mar 2006 |
18     C \==========================================================/
19     IMPLICIT NONE
20    
21     C === Global variables ===
22     #include "SIZE.h"
23     #include "EEPARAMS.h"
24     #include "PARAMS.h"
25     #include "GRID.h"
26     #include "SEAICE_SIZE.h"
27     #include "SEAICE_PARAMS.h"
28    
29     #ifdef ALLOW_AUTODIFF_TAMC
30     # include "tamc.h"
31     #endif
32    
33     C === Routine arguments ===
34     C iStep :: Sub-time-step number
35     C myTime :: Simulation time
36     C myIter :: Simulation timestep number
37     C myThid :: My Thread Id. number
38     INTEGER iStep
39     _RL myTime
40     INTEGER myIter
41     INTEGER myThid
42     C strain rate tensor
43     _RL e11 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
44     _RL e22 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
45     _RL e12 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
46     C
47     _RL zMin (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
48     _RL zMax (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
49     _RL hEffM (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
50     C
51     _RL press0(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
52     _RL press (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
53     C bulk viscosity
54     _RL eta (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
55     C shear viscosity
56     _RL zeta (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
57     CEndOfInterface
58    
59     #if ( defined (SEAICE_CGRID) && defined (SEAICE_ALLOW_DYNAMICS) )
60     C === Local variables ===
61     C i,j,bi,bj - Loop counters
62     C e11, e12, e22 - components of strain rate tensor
63     C ecm2 - inverse of square of eccentricity of yield curve
64     INTEGER i, j, bi, bj
65     _RL ECM2, deltaC, deltaCreg
66     _RL e12C (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
67     #ifdef SEAICE_ALLOW_TEM
68     _RL etaMax, etaDen
69     #endif /* SEAICE_ALLOW_TEM */
70 torge 1.2 #ifdef SEAICE_ZETA_SMOOTHREG
71     _RL argTmp
72     #endif /* SEAICE_ZETA_SMOOTHREG */
73 torge 1.1
74     C-- FIRST SET UP BASIC CONSTANTS
75     ecm2=0. _d 0
76     IF ( SEAICE_eccen .NE. 0. _d 0 ) ecm2=ONE/(SEAICE_eccen**2)
77     C
78     DO bj=myByLo(myThid),myByHi(myThid)
79     DO bi=myBxLo(myThid),myBxHi(myThid)
80     C need to do this beforehand for easier vectorization after
81     C TAFization
82     DO j=1-Oly+1,sNy+Oly-1
83     DO i=1-Olx+1,sNx+Olx-1
84     e12C(i,j) = 0.25 * ( e12(I,J ,bi,bj) + e12(I+1,J ,bi,bj)
85     & + e12(I,J+1,bi,bj) + e12(I+1,J+1,bi,bj) )
86     ENDDO
87     ENDDO
88     DO j=1-Oly+1,sNy+Oly-1
89     DO i=1-Olx+1,sNx+Olx-1
90     deltaC = (e11(i,j,bi,bj)**2+e22(i,j,bi,bj)**2)*(ONE+ecm2)
91     & + 4. _d 0*ecm2*e12C(i,j)**2
92     & + 2. _d 0*e11(i,j,bi,bj)*e22(i,j,bi,bj)*(ONE-ecm2)
93     #ifdef ALLOW_AUTODIFF_TAMC
94     C avoid sqrt of 0
95     IF ( deltaC .GT. 0. _d 0 ) deltaC = SQRT(deltaC)
96     #else
97     deltaC = SQRT(deltaC)
98     #endif /* ALLOW_AUTODIFF_TAMC */
99     deltaCreg = MAX(deltaC,SEAICE_EPS)
100     C "replacement pressure"
101     zeta (I,J,bi,bj) = HALF*press0(I,J,bi,bj)/deltaCreg
102     C put min and max viscosities in
103     #ifdef SEAICE_ZETA_SMOOTHREG
104     C regularize zeta to zmax with a smooth tanh-function instead
105     C of a min(zeta,zmax). This improves convergence of iterative
106     C solvers (Lemieux and Tremblay 2009, JGR). No effect on EVP
107 torge 1.2 IF ( ZMAX(I,J,bi,bj) .GT. 0. _d 0 ) THEN
108     argTmp = exp(-2.*zeta(I,J,bi,bj)/ZMAX(I,J,bi,bj))
109     zeta (I,J,bi,bj) = ZMAX(I,J,bi,bj)
110     & *(1. _d 0 - argTmp)/(1. _d 0 + argTmp)
111     ELSE
112     zeta (I,J,bi,bj) = ZMAX(I,J,bi,bj)
113     ENDIF
114 torge 1.1 #else
115     zeta (I,J,bi,bj) = MIN(ZMAX(I,J,bi,bj),zeta(I,J,bi,bj))
116     #endif /* SEAICE_ZETA_SMOOTHREG */
117 torge 1.2 zeta (I,J,bi,bj) = MAX(ZMIN(I,J,bi,bj),zeta(I,J,bi,bj))
118 torge 1.1 C set viscosities to zero at hEffM flow pts
119     zeta (I,J,bi,bj) = zeta(I,J,bi,bj)*HEFFM(I,J,bi,bj)
120     eta (I,J,bi,bj) = ECM2*zeta(I,J,bi,bj)
121     press(I,J,bi,bj) = TWO *zeta(I,J,bi,bj)*deltaC
122     ENDDO
123     ENDDO
124     #ifdef SEAICE_ALLOW_TEM
125     IF ( SEAICEuseTEM ) THEN
126     DO j=1-Oly+1,sNy+Oly-1
127     DO i=1-Olx+1,sNx+Olx-1
128     etaDen = (e11(I,J,bi,bj)-e22(I,J,bi,bj))**2
129     & + 4. _d 0*e12C(i,j)**2
130     etaDen = SQRT(MAX(SEAICE_EPS_SQ,etaDen))
131     etaMax = ( 0.5 _d 0*press(I,J,bi,bj)-zeta(I,J,bi,bj)
132     & *( e11(I,J,bi,bj)+e22(I,J,bi,bj) )
133     & )/etaDen
134     eta(I,J,bi,bj) = MIN(eta(I,J,bi,bj),etaMax)
135     ENDDO
136     ENDDO
137     ENDIF
138     #endif /* SEAICE_ALLOW_TEM */
139     ENDDO
140     ENDDO
141    
142     #endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID */
143     RETURN
144     END
  ViewVC Help
Powered by ViewVC 1.1.22