/[MITgcm]/MITgcm/pkg/seaice/seaice_calc_viscosities.F
ViewVC logotype

Annotation of /MITgcm/pkg/seaice/seaice_calc_viscosities.F

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

Revision 1.3 - (hide annotations) (download)
Mon Mar 20 21:36:12 2006 UTC (18 years, 1 month ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint58l_post, checkpoint58e_post, checkpoint58u_post, checkpoint58w_post, checkpoint58r_post, checkpoint58n_post, checkpoint58x_post, checkpoint58t_post, checkpoint58h_post, checkpoint58q_post, checkpoint58j_post, checkpoint59, checkpoint58f_post, checkpoint58d_post, checkpoint58c_post, checkpoint58i_post, checkpoint58g_post, checkpoint58o_post, checkpoint58y_post, checkpoint58k_post, checkpoint58v_post, checkpoint58s_post, checkpoint58p_post, checkpoint58m_post
Changes since 1.2: +51 -3 lines 
  seaice: add an EVP solver following Hunke and Dukowicz (1997) and the
  documentation of CICE. Turn on by defining SEAICE_ALLOW_EVP in
  SEAICE_OPTIONS.h and SEAICEuseEVP=.true. in data.seaice. Works only
  with SEAICE_CGRID defined.
  Use at own risk.
1 mlosch 1.3 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_calc_viscosities.F,v 1.2 2006/03/18 00:32:43 mlosch Exp $
2 mlosch 1.1 C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5    
6     CStartOfInterface
7     SUBROUTINE SEAICE_CALC_VISCOSITIES(
8     I uFld, vFld, zMin, zMax, hEffM, press0,
9     O eta, zeta, press,
10 mlosch 1.3 #ifdef SEAICE_ALLOW_EVP
11     O seaice_div, seaice_tension, seaice_shear,
12     #endif /* SEAICE_ALLOW_EVP */
13 mlosch 1.1 I myThid )
14     C /==========================================================\
15     C | SUBROUTINE SEAICE_CALC_VISCOSITIES |
16     C | o compute shear and bulk viscositites eta, zeta and the |
17     C | ice strength P |
18     C | (see Zhang and Hibler, JGR, 102, 8691-8702, 1997) |
19     C |==========================================================|
20     C | written by Martin Losch, Mar 2006 |
21     C \==========================================================/
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     #include "SEAICE_PARAMS.h"
30    
31     #ifdef ALLOW_AUTODIFF_TAMC
32     # include "tamc.h"
33     #endif
34    
35     C === Routine arguments ===
36     C myThid - Thread no. that called this routine.
37     INTEGER myThid
38     _RL uFld (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
39     _RL vFld (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
40     _RL zMin (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
41     _RL zMax (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
42     _RL hEffM (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
43     _RL press0(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
44     _RL press (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
45     _RL eta (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
46     _RL zeta (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
47 mlosch 1.3 #ifdef SEAICE_ALLOW_EVP
48     _RL seaice_div (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
49     _RL seaice_tension(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
50     _RL seaice_shear (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
51     #endif /* SEAICE_ALLOW_EVP */
52 mlosch 1.1 CEndOfInterface
53    
54     #ifdef SEAICE_CGRID
55     #ifdef SEAICE_ALLOW_DYNAMICS
56     C === Local variables ===
57     C i,j,bi,bj - Loop counters
58     C e11, e12, e22 - components of strain rate tensor
59     C ecm2 - inverse of square of eccentricity of yield curve
60     INTEGER i, j, bi, bj
61     _RL e11, e12, e22
62     _RL ECM2, DELT1, DELT2
63 mlosch 1.3 #ifdef SEAICE_ALLOW_EVP
64     _RL delta (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65     #endif
66 mlosch 1.1
67     C-- FIRST SET UP BASIC CONSTANTS
68     ecm2=0. _d 0
69     IF ( SEAICE_eccen .NE. 0. _d 0 ) ecm2=ONE/(SEAICE_eccen**2)
70     C
71     DO bj=myByLo(myThid),myByHi(myThid)
72     DO bi=myBxLo(myThid),myBxHi(myThid)
73 mlosch 1.3 DO j=1-Oly+1,sNy+Oly-1
74     DO i=1-Olx+1,sNx+Olx-1
75 mlosch 1.1 C NOW EVALUATE STRAIN RATES
76     e11= _recip_dxF(I,J,bi,bj) *
77 mlosch 1.2 & (uFld(I+1,J,bi,bj)-uFld(I,J,bi,bj))
78 mlosch 1.1 & -HALF*
79 mlosch 1.2 & (vFld(I,J,bi,bj)+vFld(I,J+1,bi,bj))
80 mlosch 1.1 & * _tanPhiAtU(I,J,bi,bj)*recip_rSphere
81     e22= _recip_dyF(I,J,bi,bj) *
82     & (vFld(I,J+1,bi,bj)-vFld(I,J,bi,bj))
83 mlosch 1.3 C one metric term is missing
84 mlosch 1.1 e12=HALF*(
85     & (uFld(I,J+1,bi,bj)+uFld(I+1,J+1,bi,bj)
86     & -uFld(I,J-1,bi,bj)-uFld(I+1,J-1,bi,bj))
87     & * 1. _d 0 / (dyC(I,J,bi,bj) + dyC(I,J-1,bi,bj))
88     & +
89     & (vFld(I+1,J+1,bi,bj)+vFld(I+1,J,bi,bj)
90     & -vFld(I-1,J+1,bi,bj)-vFld(I-1,J,bi,bj))
91     & * 1. _d 0 / (dxC(I,J,bi,bj) + dxC(I-1,J,bi,bj))
92     & +HALF*
93     & (uFld(I, J, bi,bj)+uFld(I+1,J, bi,bj))
94     & * _tanPhiAtU(I,J,bi,bj)*recip_rSphere)
95     C NOW EVALUATE VISCOSITIES
96     DELT1=(e11**2+e22**2)*(ONE+ECM2)
97     & +4. _d 0*ECM2*e12**2
98     & +TWO*e11*e22*(ONE-ECM2)
99     IF ( DELT1 .LE. SEAICE_EPS_SQ ) THEN
100     DELT2=SEAICE_EPS
101     ELSE
102     DELT2=SQRT(DELT1)
103     ENDIF
104     zeta(I,J,bi,bj)=HALF*PRESS0(I,J,bi,bj)/DELT2
105     C NOW PUT MIN AND MAX VISCOSITIES IN
106     zeta(I,J,bi,bj)=MIN(ZMAX(I,J,bi,bj),zeta(I,J,bi,bj))
107     zeta(I,J,bi,bj)=MAX(ZMIN(I,J,bi,bj),zeta(I,J,bi,bj))
108     C NOW SET VISCOSITIES TO ZERO AT HEFFMFLOW PTS
109     zeta (I,J,bi,bj) = zeta(I,J,bi,bj)*HEFFM(I,J,bi,bj)
110     eta (I,J,bi,bj) = ECM2*zeta(I,J,bi,bj)
111     press(I,J,bi,bj) = TWO*zeta(I,J,bi,bj)*DELT2
112 mlosch 1.3 #ifdef SEAICE_ALLOW_EVP
113     delta(I,J) = DELT2
114     seaice_div (I,J,bi,bj) = (e11+e22)/DELT2*hEffM(I,J,bi,bj)
115     seaice_tension(I,J,bi,bj) = (e11-e22)/DELT2*hEffM(I,J,bi,bj)
116     C recompute e12 on Z point so that seaice_shear and seaice_sigma12
117     C are also defined on Z points (one half and two drop out)
118     seaice_shear (I,J,bi,bj) = (
119     & (uFld(I ,J ,bi,bj) * _dxC(I ,J ,bi,bj)
120     & -uFld(I ,J-1,bi,bj) * _dxC(I ,J-1,bi,bj)
121     & +vFld(I ,J ,bi,bj) * _dyC(I ,J ,bi,bj)
122     & -vFld(I-1,J ,bi,bj) * _dyC(I-1,J ,bi,bj))
123     & * recip_rAz(I,J,bi,bj)
124     & +
125     & 0.25 _d 0 * (uFld(I,J,bi,bj)+uFld(I ,J-1,bi,bj))
126     & * ( _tanPhiAtU(I,J,bi,bj) + _tanPhiAtU(I,J-1,bi,bj) )
127     & *recip_rSphere
128     & )
129     C one metric term is missing
130     #endif /* SEAICE_ALLOW_EVP */
131     ENDDO
132     ENDDO
133     #ifdef SEAICE_ALLOW_EVP
134     DO j=1,sNy
135     DO i=1,sNx
136     DELT2 = 0.25 * ( delta(I,J ) + delta(I-1,J )
137     & + delta(I,J-1) + delta(I-1,J-1) )
138     seaice_shear(I,J,bi,bj) =
139     & seaice_shear(I,J,bi,bj)/DELT2
140     & *hEffM(I ,J ,bi,bj)*hEffM(I-1,J ,bi,bj)
141     & *hEffM(I ,J-1,bi,bj)*hEffM(I-1,J-1,bi,bj)
142 mlosch 1.1 ENDDO
143     ENDDO
144 mlosch 1.3 #endif /* SEAICE_ALLOW_EVP */
145 mlosch 1.1 ENDDO
146     ENDDO
147    
148     C-- Update overlap regions
149     _EXCH_XY_R8( eta, myThid)
150     _EXCH_XY_R8( zeta, myThid)
151     _EXCH_XY_R8(press, myThid)
152 mlosch 1.3 #ifdef SEAICE_ALLOW_EVP
153     _EXCH_XY_R8(seaice_div , myThid)
154     _EXCH_XY_R8(seaice_tension, myThid)
155     _EXCH_XY_R8(seaice_shear , myThid)
156     #endif /* SEAICE_ALLOW_EVP */
157 mlosch 1.1 #endif /* SEAICE_ALLOW_DYNAMICS */
158     #endif /* SEAICE_CGRID */
159     RETURN
160     END
  ViewVC Help
Powered by ViewVC 1.1.22