/[MITgcm]/MITgcm_contrib/dgoldberg/streamice/streamice_driving_stress_ppm.F
ViewVC logotype

Contents of /MITgcm_contrib/dgoldberg/streamice/streamice_driving_stress_ppm.F

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


Revision 1.4 - (show annotations) (download)
Wed Aug 27 19:29:13 2014 UTC (10 years, 10 months ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +1 -1 lines
Error occurred while calculating annotation data.
updating contrib streamice repo with latest files, and separated out convergence checks; and parameterised maximum iteration counts and interface w shelfice for coupling

1 C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_driving_stress_ppm.F,v 1.1 2013/06/12 21:30:22 dgoldberg Exp $
2 C $Name: $
3
4 #include "STREAMICE_OPTIONS.h"
5
6 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7
8 CBOP
9 SUBROUTINE STREAMICE_DRIVING_STRESS_PPM( myThid )
10 ! O taudx,
11 ! O taudy )
12
13 C /============================================================\
14 C | SUBROUTINE |
15 C | o |
16 C |============================================================|
17 C | |
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 "STREAMICE.h"
27 #include "STREAMICE_CG.h"
28
29 C !INPUT/OUTPUT ARGUMENTS
30 INTEGER myThid
31 ! _RL taudx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
32 ! _RL taudx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
33
34 #ifdef ALLOW_STREAMICE
35
36
37 C LOCAL VARIABLES
38 INTEGER i, j, bi, bj, k, l,
39 & Gi, Gj
40 LOGICAL at_west_bdry, at_east_bdry,
41 & at_north_bdry, at_south_bdry
42 _RL sx, sy, diffx, diffy, neu_val
43
44 IF (myXGlobalLo.eq.1) at_west_bdry = .true.
45 IF (myYGlobalLo.eq.1) at_south_bdry = .true.
46 IF (myXGlobalLo-1+sNx*nSx.eq.Nx)
47 & at_east_bdry = .false.
48 IF (myYGlobalLo-1+sNy*nSy.eq.Ny)
49 & at_north_bdry = .false.
50
51 DO bj = myByLo(myThid), myByHi(myThid)
52 DO bi = myBxLo(myThid), myBxHi(myThid)
53 DO j=1-OLy,sNy+OLy
54 DO i=1-OLx,sNx+OLx
55 taudx_SI(i,j,bi,bj) = 0. _d 0
56 taudy_SI(i,j,bi,bj) = 0. _d 0
57 ENDDO
58 ENDDO
59 ENDDO
60 ENDDO
61
62 DO bj = myByLo(myThid), myByHi(myThid)
63 DO bi = myBxLo(myThid), myBxHi(myThid)
64
65 DO i=0,sNx+1
66 DO j=0,sNy+1
67
68 diffx = 0. _d 0
69 diffy = 0. _d 0
70 sx = 0. _d 0
71 sy = 0. _d 0
72
73 Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
74 Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
75
76 IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
77
78 ! we are in an "active" cell
79
80 IF (Gi.eq.1.AND..NOT.STREAMICE_EW_periodic) THEN
81
82 ! western boundary - only one sided possible
83
84 IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0) THEN
85
86 ! cell to east is active
87
88 sx = (surf_el_streamice(i+1,j,bi,bj)-
89 & surf_el_streamice(i,j,bi,bj))/dxC(i+1,j,bi,bj)
90 ELSE
91
92 ! cell to east is empty
93
94 sx = 0. _d 0
95 ENDIF
96
97 DO k=0,1
98 DO l=0,1
99 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
100 taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
101 & 0.25 * streamice_density * gravity *
102 & (streamice_bg_surf_slope_x+sx) *
103 & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
104 ENDIF
105 ENDDO
106 ENDDO
107
108 ELSEIF (Gi.eq.Nx.AND..NOT.STREAMICE_EW_periodic) THEN
109
110 ! eastern boundary - only one sided possible
111
112 IF (STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN
113
114 ! cell to west is active
115
116 sx = (surf_el_streamice(i,j,bi,bj)-
117 & surf_el_streamice(i-1,j,bi,bj))/dxC(i,j,bi,bj)
118 ELSE
119
120 ! cell to west is inactive
121
122 sx = 0. _d 0
123 ENDIF
124
125 DO k=0,1
126 DO l=0,1
127 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
128 taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
129 & 0.25 * streamice_density * gravity *
130 & (streamice_bg_surf_slope_x+sx) *
131 & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
132 ENDIF
133 ENDDO
134 ENDDO
135
136 ELSE
137
138 ! interior (west-east) cell
139
140 IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0 .and.
141 & STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN
142
143 k = 0
144 DO l=0,1
145 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
146 taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
147 & streamice_density * gravity * (1./6.) *
148 & (-2.*surf_el_streamice(i-1,j,bi,bj) +
149 & surf_el_streamice(i,j,bi,bj) +
150 & surf_el_streamice(i+1,j,bi,bj) +
151 & 3.*streamice_bg_surf_slope_x * dxF(i,j,bi,bj)) *
152 & H_streamice(i,j,bi,bj) * .5 * dyF(i,j,bi,bj)
153 ENDIF
154 ENDDO
155
156 k = 1
157 DO l=0,1
158 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
159 taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
160 & streamice_density * gravity * (1./6.) *
161 & (-surf_el_streamice(i-1,j,bi,bj) -
162 & surf_el_streamice(i,j,bi,bj) +
163 & 2*surf_el_streamice(i+1,j,bi,bj) +
164 & 3.*streamice_bg_surf_slope_x * dxF(i,j,bi,bj)) *
165 & H_streamice(i,j,bi,bj) * .5 * dyF(i,j,bi,bj)
166 ENDIF
167 ENDDO
168
169
170 ELSE
171
172 IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0) THEN
173
174 sx = (surf_el_streamice(i+1,j,bi,bj)-
175 & surf_el_streamice(i,j,bi,bj))/dxC(i+1,j,bi,bj)
176
177 ELSEIF (STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN
178
179 sx = (surf_el_streamice(i,j,bi,bj)-
180 & surf_el_streamice(i-1,j,bi,bj))/dxC(i,j,bi,bj)
181
182 ELSE
183
184 sx = 0. _d 0
185
186 ENDIF
187
188 DO k=0,1
189 DO l=0,1
190 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
191 taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
192 & 0.25 * streamice_density * gravity *
193 & (streamice_bg_surf_slope_x+sx) *
194 & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
195 ENDIF
196 ENDDO
197 ENDDO
198
199 ENDIF
200
201 ENDIF
202
203 !!!!!!!! DONE WITH X-GRADIENT
204
205 IF (Gj.eq.1.AND..NOT.STREAMICE_NS_periodic) THEN
206
207 ! western boundary - only one sided possible
208
209 IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0) THEN
210
211 ! cell to east is active
212
213 sy = (surf_el_streamice(i,j+1,bi,bj)-
214 & surf_el_streamice(i,j,bi,bj))/dyC(i,j+1,bi,bj)
215 ELSE
216
217 ! cell to east is empty
218
219 sy = 0. _d 0
220 ENDIF
221
222 DO k=0,1
223 DO l=0,1
224 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
225 taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
226 & 0.25 * streamice_density * gravity *
227 & (streamice_bg_surf_slope_y+sy) *
228 & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
229 ENDIF
230 ENDDO
231 ENDDO
232
233 ELSEIF (Gj.eq.Ny.AND..NOT.STREAMICE_NS_periodic) THEN
234
235 ! eastern boundary - only one sided possible
236
237 IF (STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN
238
239 ! cell to west is active
240
241 sy = (surf_el_streamice(i,j,bi,bj)-
242 & surf_el_streamice(i,j-1,bi,bj))/dyC(i,j,bi,bj)
243
244 ELSE
245
246 ! cell to west is inactive
247
248 sy = 0. _d 0
249 ENDIF
250
251 DO k=0,1
252 DO l=0,1
253 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
254 taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
255 & 0.25 * streamice_density * gravity *
256 & (streamice_bg_surf_slope_y+sy) *
257 & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
258 ENDIF
259 ENDDO
260 ENDDO
261
262 ELSE
263
264 ! interior (west-east) cell
265
266 IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0 .and.
267 & STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN
268
269 l = 0
270 DO k=0,1
271 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
272 taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
273 & streamice_density * gravity * (1./6.) *
274 & (-2.*surf_el_streamice(i,j-1,bi,bj) +
275 & surf_el_streamice(i,j,bi,bj) +
276 & surf_el_streamice(i,j+1,bi,bj) +
277 & 3.*streamice_bg_surf_slope_y * dyF(i,j,bi,bj)) *
278 & H_streamice(i,j,bi,bj) * .5 * dxF(i,j,bi,bj)
279 ENDIF
280 ENDDO
281
282 l = 1
283 DO k=0,1
284 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
285 taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
286 & streamice_density * gravity * (1./6.) *
287 & (-surf_el_streamice(i,j-1,bi,bj) -
288 & surf_el_streamice(i,j,bi,bj) +
289 & 2*surf_el_streamice(i,j+1,bi,bj) +
290 & 3.*streamice_bg_surf_slope_y * dyF(i,j,bi,bj)) *
291 & H_streamice(i,j,bi,bj) * .5 * dxF(i,j,bi,bj)
292 ENDIF
293 ENDDO
294
295
296 ELSE
297
298 IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0) THEN
299
300 sy = (surf_el_streamice(i,j+1,bi,bj)-
301 & surf_el_streamice(i,j,bi,bj))/dxC(i,j+1,bi,bj)
302
303 ELSEIF (STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN
304
305 sy = (surf_el_streamice(i,j,bi,bj)-
306 & surf_el_streamice(i,j-1,bi,bj))/dxC(i,j,bi,bj)
307
308 ELSE
309
310 sy = 0. _d 0
311
312 ENDIF
313
314 DO k=0,1
315 DO l=0,1
316 IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
317 taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
318 & 0.25 * streamice_density * gravity *
319 & (streamice_bg_surf_slope_y+sy) *
320 & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
321 ENDIF
322 ENDDO
323 ENDDO
324
325 ENDIF
326
327 ENDIF
328
329 ! DO k=0,1
330 ! DO l=0,1
331 ! IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
332 ! taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
333 ! & 0.25 * streamice_density * gravity *
334 ! & (streamice_bg_surf_slope_x+sx) *
335 ! & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
336 ! taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
337 ! & 0.25 * streamice_density * gravity *
338 ! & (streamice_bg_surf_slope_y+sy) *
339 ! & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
340 !
341 ! ENDIF
342 ! ENDDO
343 ! ENDDO
344
345 IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) then
346 #ifdef USE_ALT_RLOW
347 neu_val = .5 * gravity *
348 & (streamice_density * H_streamice (i,j,bi,bj) ** 2 -
349 & streamice_density_ocean_avg * R_low_si(i,j,bi,bj) ** 2)
350 #else
351 neu_val = .5 * gravity *
352 & (streamice_density * H_streamice (i,j,bi,bj) ** 2 -
353 & streamice_density_ocean_avg * R_low(i,j,bi,bj) ** 2)
354 #endif
355 ELSE
356 neu_val = .5 * gravity *
357 & (1-streamice_density/streamice_density_ocean_avg) *
358 & streamice_density * H_streamice(i,j,bi,bj) ** 2
359 ENDIF
360
361 IF ((STREAMICE_ufacemask(i,j,bi,bj) .eq. 2)
362 & .OR. (STREAMICE_hmask(i-1,j,bi,bj) .eq. 0)
363 & .OR. (STREAMICE_hmask(i-1,j,bi,bj) .eq. 2) ) THEN ! left face of the cell is at a stress boundary
364 ! the depth-integrated longitudinal stress is equal to the difference of depth-integrated pressure on either side of the face
365 ! on the ice side, it is rho g h^2 / 2
366 ! on the ocean side, it is rhow g (delta OD)^2 / 2
367 ! OD can be zero under the ice; but it is ASSUMED on the ice-free side of the face, topography elevation is not above the base of the
368 ! ice in the current cell
369
370 taudx_SI(i,j,bi,bj) = taudx_SI(i,j,bi,bj) -
371 & .5 * dyG(i,j,bi,bj)*(neu_val+streamice_addl_backstress)
372 taudx_SI(i,j+1,bi,bj) = taudx_SI(i,j+1,bi,bj) -
373 & .5 * dyG(i,j,bi,bj)*(neu_val+streamice_addl_backstress)
374 ENDIF
375
376 IF ((STREAMICE_ufacemask(i+1,j,bi,bj) .eq. 2)
377 & .OR. (STREAMICE_hmask(i+1,j,bi,bj) .eq. 0)
378 & .OR. (STREAMICE_hmask(i+1,j,bi,bj) .eq. 2) ) THEN
379
380 taudx_SI(i+1,j,bi,bj) = taudx_SI(i+1,j,bi,bj) +
381 & .5 * dyG(i+1,j,bi,bj)*(neu_val+streamice_addl_backstress) ! note negative sign is due to direction of normal vector
382 taudx_SI(i+1,j+1,bi,bj) = taudx_SI(i+1,j+1,bi,bj) +
383 & .5 * dyG(i+1,j,bi,bj)*(neu_val+streamice_addl_backstress)
384 ENDIF
385
386 IF ((STREAMICE_vfacemask(i,j,bi,bj) .eq. 2)
387 & .OR. (STREAMICE_hmask(i,j-1,bi,bj) .eq. 0)
388 & .OR. (STREAMICE_hmask(i,j-1,bi,bj) .eq. 2) ) THEN
389
390 taudy_SI(i,j,bi,bj) = taudy_SI(i,j,bi,bj) -
391 & .5 * dxG(i,j,bi,bj)*(neu_val+streamice_addl_backstress)
392 taudy_SI(i+1,j,bi,bj) = taudy_SI(i+1,j,bi,bj) -
393 & .5 * dxG(i,j,bi,bj)*(neu_val+streamice_addl_backstress)
394 ENDIF
395
396 IF ((STREAMICE_vfacemask(i,j+1,bi,bj) .eq. 2)
397 & .OR. (STREAMICE_hmask(i,j+1,bi,bj) .eq. 0)
398 & .OR. (STREAMICE_hmask(i,j+1,bi,bj) .eq. 2) ) THEN
399
400 taudy_SI(i,j+1,bi,bj) = taudy_SI(i,j+1,bi,bj) +
401 & .5 * dxG(i,j+1,bi,bj)*(neu_val+streamice_addl_backstress)
402 taudy_SI(i+1,j+1,bi,bj) = taudy_SI(i+1,j+1,bi,bj) +
403 & .5 * dxG(i,j+1,bi,bj)*(neu_val+streamice_addl_backstress)
404 ENDIF
405
406 ENDIF
407 ENDDO
408 ENDDO
409 ENDDO
410 ENDDO
411
412
413
414 #endif
415 RETURN
416 END
417
418

  ViewVC Help
Powered by ViewVC 1.1.22