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

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

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

Revision 1.11 - (hide annotations) (download)
Wed Mar 18 12:58:17 2009 UTC (15 years, 2 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint61n, checkpoint61o, checkpoint61l, checkpoint61m
Changes since 1.10: +62 -1 lines 
change the discretization of strain rates to be consistent with
the overall finite-volume discretization. This change alone changes
all seaice verification experiments, therefore the old variant is
keep for now with a CPP-flag SEAICE_OLD_AND_BAD_DISCRETIZATION
1 mlosch 1.10 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_calc_strainrates.F,v 1.9 2007/11/14 15:55:48 mlosch Exp $
2 mlosch 1.1 C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5    
6     CStartOfInterface
7 jmc 1.8 SUBROUTINE SEAICE_CALC_STRAINRATES(
8 mlosch 1.1 I uFld, vFld,
9     O e11, e22, e12,
10 mlosch 1.9 I kSize, iStep, myTime, myIter, myThid )
11 mlosch 1.1 C /==========================================================\
12     C | SUBROUTINE SEAICE_CALC_STRAINRATES |
13     C | o compute strain rates from ice velocities |
14     C |==========================================================|
15     C | written by Martin Losch, Apr 2007 |
16     C \==========================================================/
17     IMPLICIT NONE
18    
19     C === Global variables ===
20     #include "SIZE.h"
21     #include "EEPARAMS.h"
22     #include "PARAMS.h"
23     #include "GRID.h"
24     #include "SEAICE_PARAMS.h"
25 mlosch 1.11 #include "SEAICE.h"
26 mlosch 1.1
27     #ifdef ALLOW_AUTODIFF_TAMC
28     # include "tamc.h"
29     #endif
30    
31     C === Routine arguments ===
32 jmc 1.8 C iStep :: Sub-time-step number
33     C myTime :: Simulation time
34     C myIter :: Simulation timestep number
35     C myThid :: My Thread Id. number
36 mlosch 1.9 C kSize :: length of 3rd dimension of velocity variables
37 jmc 1.8 INTEGER iStep
38     _RL myTime
39     INTEGER myIter
40 mlosch 1.1 INTEGER myThid
41 mlosch 1.9 INTEGER kSize
42 mlosch 1.1 C ice velocities
43 mlosch 1.9 _RL uFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,kSize,nSx,nSy)
44     _RL vFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,kSize,nSx,nSy)
45 mlosch 1.1 C strain rate tensor
46     _RL e11 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
47     _RL e22 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
48     _RL e12 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
49     CEndOfInterface
50    
51     #ifdef SEAICE_CGRID
52     #ifdef SEAICE_ALLOW_DYNAMICS
53     C === Local variables ===
54     C i,j,bi,bj - Loop counters
55     INTEGER i, j, bi, bj
56 mlosch 1.2 C hFacU, hFacV - determine the no-slip boundary condition
57     INTEGER k
58 mlosch 1.11 _RS hFacU, hFacV, noSlipFac
59 mlosch 1.2
60 mlosch 1.4 k = 1
61 mlosch 1.11 noSlipFac = 0. _d 0
62     IF ( SEAICE_no_slip ) noSlipFac = 1. _d 0
63 mlosch 1.1 C
64 mlosch 1.11 #ifndef SEAICE_OLD_AND_BAD_DISCRETIZATION
65     DO bj=myByLo(myThid),myByHi(myThid)
66     DO bi=myBxLo(myThid),myBxHi(myThid)
67     DO j=1-Oly,sNy+Oly-1
68     DO i=1-Olx,sNx+Olx-1
69     C evaluate strain rates
70     e11(I,J,bi,bj) = _recip_dxF(I,J,bi,bj) *
71     & (uFld(I+1,J,1,bi,bj)-uFld(I,J,1,bi,bj))
72     & +HALF*
73     & (vFld(I,J,1,bi,bj)+vFld(I,J+1,1,bi,bj))
74     & * k2AtC(I,J,bi,bj)
75     e22(I,J,bi,bj) = _recip_dyF(I,J,bi,bj) *
76     & (vFld(I,J+1,1,bi,bj)-vFld(I,J,1,bi,bj))
77     & +HALF*
78     & (uFld(I,J,1,bi,bj)+uFld(I+1,J,1,bi,bj))
79     & * k1AtC(I,J,bi,bj)
80     C one metric term is missing
81     ENDDO
82     ENDDO
83     DO j=1-Oly+1,sNy+Oly
84     DO i=1-Olx+1,sNx+Olx
85     hFacU = _maskW(i,j,k,bi,bj) - _maskW(i,j-1,k,bi,bj)
86     hFacV = _maskS(i,j,k,bi,bj) - _maskS(i-1,j,k,bi,bj)
87     e12(I,J,bi,bj) = HALF*(
88     & ( uFld(I,J,1,bi,bj) - uFld(I ,J-1,1,bi,bj) )
89     & * _recip_dyU(I,J,bi,bj)
90     & + ( vFld(I,J,1,bi,bj) - vFld(I-1,J ,1,bi,bj) )
91     & * _recip_dxV(I,J,bi,bj)
92     & - k1AtZ(I,J,bi,bj)
93     & * 0.5 _d 0 * (vFld(I,J,1,bi,bj)+vFld(I-1,J ,1,bi,bj))
94     & - k2AtZ(I,J,bi,bj)
95     & * 0.5 _d 0 * (uFld(I,J,1,bi,bj)+uFld(I ,J-1,1,bi,bj))
96     & )
97     & *maskC(I ,J ,k,bi,bj)*maskC(I-1,J ,k,bi,bj)
98     & *maskC(I ,J-1,k,bi,bj)*maskC(I-1,J-1,k,bi,bj)
99     & + 2.0 _d 0 * noSlipFac * (
100     & ( uFld(I,J,1,bi,bj) + uFld(I ,J-1,1,bi,bj) )
101     & * _recip_dyU(I,J,bi,bj) * hFacU
102     & + ( vFld(I,J,1,bi,bj) + vFld(I-1,J ,1,bi,bj) )
103     & * _recip_dxV(I,J,bi,bj) * hFacV
104     & )
105     C no slip at the boundary implies u(j)+u(j-1)=0 and v(i)+v(i-1)=0
106     C accross the boundary; this is already accomplished by masking so
107     C that the following lines are not necessary
108     c$$$ & - hFacV * k1AtZ(I,J,bi,bj)
109     c$$$ & * 0.5 _d 0 * (vFld(I,J,1,bi,bj)+vFld(I-1,J ,1,bi,bj))
110     c$$$ & - hFacU * k2AtZ(I,J,bi,bj)
111     c$$$ & * 0.5 _d 0 * (uFld(I,J,1,bi,bj)+uFld(I ,J-1,1,bi,bj))
112     ENDDO
113     ENDDO
114    
115     c$$$ ENDIF
116     ENDDO
117     ENDDO
118     #else
119     C this the old and incomplete discretization, here I also erroneously
120     C used finite-volumes to discretize the strain rates
121 mlosch 1.1 DO bj=myByLo(myThid),myByHi(myThid)
122     DO bi=myBxLo(myThid),myBxHi(myThid)
123 mlosch 1.5 DO j=1-Oly,sNy+Oly-1
124     DO i=1-Olx,sNx+Olx-1
125     C evaluate strain rates
126 mlosch 1.6 e11(I,J,bi,bj) = _recip_dxF(I,J,bi,bj) *
127 mlosch 1.9 & (uFld(I+1,J,1,bi,bj)-uFld(I,J,1,bi,bj))
128 mlosch 1.1 & -HALF*
129 mlosch 1.9 & (vFld(I,J,1,bi,bj)+vFld(I,J+1,1,bi,bj))
130 mlosch 1.1 & * _tanPhiAtU(I,J,bi,bj)*recip_rSphere
131 mlosch 1.6 e22(I,J,bi,bj) = _recip_dyF(I,J,bi,bj) *
132 mlosch 1.9 & (vFld(I,J+1,1,bi,bj)-vFld(I,J,1,bi,bj))
133 mlosch 1.1 C one metric term is missing
134 mlosch 1.5 ENDDO
135     ENDDO
136     DO j=1-Oly+1,sNy+Oly
137     DO i=1-Olx+1,sNx+Olx
138 mlosch 1.6 e12(I,J,bi,bj) = HALF*(
139 mlosch 1.9 & (uFld(I ,J ,1,bi,bj) * _dxC(I ,J ,bi,bj)
140     & -uFld(I ,J-1,1,bi,bj) * _dxC(I ,J-1,bi,bj)
141     & +vFld(I ,J ,1,bi,bj) * _dyC(I ,J ,bi,bj)
142     & -vFld(I-1,J ,1,bi,bj) * _dyC(I-1,J ,bi,bj))
143 mlosch 1.1 & * recip_rAz(I,J,bi,bj)
144     & +
145 mlosch 1.9 & 0.25 _d 0 * (uFld(I,J,1,bi,bj)+uFld(I ,J-1,1,bi,bj))
146 mlosch 1.1 & * ( _tanPhiAtU(I,J,bi,bj) + _tanPhiAtU(I,J-1,bi,bj) )
147     & *recip_rSphere
148     & )
149 mlosch 1.4 & *maskC(I ,J ,k,bi,bj)*maskC(I-1,J ,k,bi,bj)
150     & *maskC(I ,J-1,k,bi,bj)*maskC(I-1,J-1,k,bi,bj)
151 mlosch 1.1 C one metric term is missing
152     ENDDO
153     ENDDO
154 mlosch 1.2 IF ( SEAICE_no_slip ) THEN
155 mlosch 1.3 C no slip boundary conditions apply only to e12
156 mlosch 1.5 DO j=1-Oly+1,sNy+Oly
157     DO i=1-Olx+1,sNx+Olx
158 mlosch 1.2 hFacU = _maskW(i,j,k,bi,bj) - _maskW(i,j-1,k,bi,bj)
159     hFacV = _maskS(i,j,k,bi,bj) - _maskS(i-1,j,k,bi,bj)
160    
161 mlosch 1.6 e12(I,J,bi,bj) = e12(I,J,bi,bj)
162 mlosch 1.10 & + recip_rAz(i,j,bi,bj) * 2. _d 0 *
163 mlosch 1.9 & ( hFacU * ( _dxC(i,j-1,bi,bj)*uFld(i,j ,1,bi,bj)
164     & + _dxC(i,j, bi,bj)*uFld(i,j-1,1,bi,bj) )
165     & + hFacV * ( _dyC(i-1,j,bi,bj)*vFld(i ,j,1,bi,bj)
166     & + _dyC(i, j,bi,bj)*vFld(i-1,j,1,bi,bj) ) )
167 jmc 1.8 & - hFacU
168 mlosch 1.9 & * 0.25 _d 0 * (uFld(I,J,1,bi,bj)+uFld(I ,J-1,1,bi,bj))
169 mlosch 1.2 & * ( _tanPhiAtU(I,J,bi,bj) + _tanPhiAtU(I,J-1,bi,bj) )
170     & *recip_rSphere
171     C one metric term is missing
172     ENDDO
173     ENDDO
174    
175     ENDIF
176 mlosch 1.1 ENDDO
177     ENDDO
178 mlosch 1.11 #endif /* SEAICE_OLD_AND_BAD_DISCRETIZATION */
179 mlosch 1.1 #endif /* SEAICE_ALLOW_DYNAMICS */
180     #endif /* SEAICE_CGRID */
181     RETURN
182     END
  ViewVC Help
Powered by ViewVC 1.1.22