/[MITgcm]/MITgcm/pkg/streamice/streamice_driving_stress.F
ViewVC logotype

Contents of /MITgcm/pkg/streamice/streamice_driving_stress.F

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


Revision 1.10 - (show annotations) (download)
Sun Oct 18 14:37:49 2015 UTC (8 years, 6 months ago) by dgoldberg
Branch: MAIN
CVS Tags: checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, HEAD
Changes since 1.9: +13 -13 lines
fix error for anisotropic grid cell width

1 C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_driving_stress.F,v 1.9 2015/07/04 10:54:24 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( 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 grd_below_sl
43 _RL sx, sy, diffx, diffy, geom_fac
44 _RL i_rhow
45 _RL avg_density
46 & (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
47 #ifdef STREAMICE_FIRN_CORRECTION
48 _RL firn_depth, h
49 #endif
50
51
52 IF (myXGlobalLo.eq.1) at_west_bdry = .true.
53 IF (myYGlobalLo.eq.1) at_south_bdry = .true.
54 IF (myXGlobalLo-1+sNx*nSx.eq.Nx)
55 & at_east_bdry = .false.
56 IF (myYGlobalLo-1+sNy*nSy.eq.Ny)
57 & at_north_bdry = .false.
58 #ifdef STREAMICE_FIRN_CORRECTION
59 firn_depth = streamice_density *
60 & streamice_firn_correction
61 & / (streamice_density-streamice_density_firn)
62 #endif
63 i_rhow = 1./streamice_density_ocean_avg
64
65 DO bj = myByLo(myThid), myByHi(myThid)
66 DO bi = myBxLo(myThid), myBxHi(myThid)
67 DO j=1-OLy+1,sNy+OLy-1
68 DO i=1-OLx+1,sNx+OLx-1
69 taudx_SI(i,j,bi,bj) = 0. _d 0
70 taudy_SI(i,j,bi,bj) = 0. _d 0
71 #ifdef STREAMICE_FIRN_CORRECTION
72 if (STREAMICE_apply_firn_correction) then
73 if (streamice_hmask(i,j,bi,bj).eq.1) then
74 h = h_streamice(i,j,bi,bj)
75 if (h.lt.firn_depth) then
76 avg_density(i,j,bi,bj) = streamice_density_firn
77 else
78 avg_density(i,j,bi,bj) = streamice_density *
79 & (h - streamice_firn_correction) / h
80 endif
81 endif
82 else
83 #endif
84 avg_density(i,j,bi,bj) = streamice_density
85 #ifdef STREAMICE_FIRN_CORRECTION
86 endif
87 #endif
88 ENDDO
89 ENDDO
90 ENDDO
91 ENDDO
92
93 DO bj = myByLo(myThid), myByHi(myThid)
94 DO bi = myBxLo(myThid), myBxHi(myThid)
95
96 DO i=1,sNx
97 DO j=1,sNy
98
99 diffx = 0. _d 0
100 diffy = 0. _d 0
101 sx = 0. _d 0
102 sy = 0. _d 0
103
104 Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
105 Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
106
107 IF (streamice_umask(i,j,bi,bj).eq.1.0) THEN
108
109 IF (streamice_hmask(i-1,j,bi,bj).eq.1.0.AND.
110 & streamice_hmask(i,j,bi,bj).eq.1.0) THEN
111
112 ! geom_fac = sqrt(rA(i-1,j,bi,bj)*recip_rA(i,j,bi,bj)*
113 ! & dxF(i,j,bi,bj)*recip_dxF(i-1,j,bi,bj))
114
115
116 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
117 & 0.25 * dyG(i,j,bi,bj) *
118 & gravity *
119 & (H_streamice(i,j,bi,bj)*avg_density(i,j,bi,bj)+
120 & H_streamice(i-1,j,bi,bj)*avg_density(i-1,j,bi,bj)) *
121 & (surf_el_streamice(i,j,bi,bj)-
122 & surf_el_streamice(i-1,j,bi,bj))
123
124 CCC
125 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
126 & streamice_density * gravity *
127 & streamice_bg_surf_slope_x * .25 * rA(i,j,bi,bj) *
128 & (H_streamice(i-1,j,bi,bj) + H_streamice(i,j,bi,bj))
129 CCC
130
131 ELSE IF (streamice_hmask(i-1,j,bi,bj).eq.1.0) THEN
132
133 IF (float_frac_streamice(i-1,j,bi,bj) .eq. 1.0) THEN
134
135 #ifdef USE_ALT_RLOW
136 IF (R_low_si(i-1,j,bi,bj) .lt. 0. _d 0) then
137 #else
138 IF (R_low(i-1,j,bi,bj) .lt. 0. _d 0) then
139 #endif
140 grd_below_sl = 1. _d 0
141 else
142 grd_below_sl = 0. _d 0
143 endif
144
145 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
146 & 0.25 * dyG(i,j,bi,bj) *
147 & gravity *
148 & (avg_density(i-1,j,bi,bj) * H_streamice(i-1,j,bi,bj)**2-
149 #ifdef USE_ALT_RLOW
150 & streamice_density_ocean_avg * grd_below_sl *
151 & R_low_si(i-1,j,bi,bj)**2)
152 #else
153 & streamice_density_ocean_avg * grd_below_sl *
154 & R_low(i-1,j,bi,bj)**2)
155 #endif
156
157 ELSE
158
159 #ifdef STREAMICE_FIRN_CORRECTION
160
161 if (STREAMICE_apply_firn_correction) then
162 if (H_streamice(i-1,j,bi,bj).lt.firn_depth) then
163 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
164 & 0.25 * dyG(i,j,bi,bj) *
165 & streamice_density_firn * gravity *
166 & (1-streamice_density_firn*i_rhow) *
167 & H_streamice(i-1,j,bi,bj)**2
168 else
169 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
170 & 0.25 * dyG(i,j,bi,bj) * gravity * (
171 & streamice_density_firn * firn_depth**2 +
172 & (h_streamice(i-1,j,bi,bj)-firn_depth) *
173 & (streamice_density_firn*firn_depth+streamice_density*
174 & (h_streamice(i-1,j,bi,bj)-streamice_firn_correction)) -
175 & streamice_density**2*i_rhow*
176 & (h_streamice(i-1,j,bi,bj)-streamice_firn_correction)**2
177 & )
178 endif
179 else
180
181 #endif
182 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
183 & 0.25 * dyG(i,j,bi,bj) *
184 & streamice_density * gravity *
185 & (1-streamice_density*i_rhow) *
186 & H_streamice(i-1,j,bi,bj)**2
187
188 #ifdef STREAMICE_FIRN_CORRECTION
189 endif
190 #endif
191
192 ENDIF
193
194 ELSE IF (streamice_hmask(i,j,bi,bj).eq.1.0) THEN
195
196 IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) THEN
197
198 #ifdef USE_ALT_RLOW
199 IF (R_low_si(i,j,bi,bj) .lt. 0. _d 0) then
200 #else
201 IF (R_low(i,j,bi,bj) .lt. 0. _d 0) then
202 #endif
203 grd_below_sl = 1. _d 0
204 else
205 grd_below_sl = 0. _d 0
206 endif
207
208 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
209 & 0.25 * dyG(i,j,bi,bj) *
210 & gravity *
211 & (avg_density(i,j,bi,bj) * H_streamice(i,j,bi,bj)**2 -
212 #ifdef USE_ALT_RLOW
213 & streamice_density_ocean_avg * grd_below_sl *
214 & R_low_si(i,j,bi,bj)**2)
215 #else
216 & streamice_density_ocean_avg * grd_below_sl *
217 & R_low(i,j,bi,bj)**2)
218 #endif
219
220 ELSE
221
222 #ifdef STREAMICE_FIRN_CORRECTION
223
224 if (STREAMICE_apply_firn_correction) then
225 if (H_streamice(i,j,bi,bj).lt.firn_depth) then
226 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
227 & 0.25 * dyG(i,j,bi,bj) *
228 & streamice_density_firn * gravity *
229 & (1-streamice_density_firn*i_rhow) *
230 & H_streamice(i,j,bi,bj)**2
231 else
232 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
233 & 0.25 * dyG(i,j,bi,bj) * gravity * (
234 & streamice_density_firn * firn_depth**2 +
235 & (h_streamice(i,j,bi,bj)-firn_depth) *
236 & (streamice_density_firn*firn_depth+streamice_density*
237 & (h_streamice(i,j,bi,bj)-streamice_firn_correction)) -
238 & streamice_density**2*i_rhow*
239 & (h_streamice(i,j,bi,bj)-streamice_firn_correction)**2
240 & )
241 endif
242 else
243
244 #endif
245 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
246 & 0.25 * dyG(i,j,bi,bj) *
247 & streamice_density * gravity *
248 & (1-streamice_density*i_rhow) *
249 & H_streamice(i,j,bi,bj)**2
250 #ifdef STREAMICE_FIRN_CORRECTION
251 endif
252 #endif
253
254
255 END IF
256 END IF
257
258 ! cells below
259
260 IF (streamice_hmask(i-1,j-1,bi,bj).eq.1.0.AND.
261 & streamice_hmask(i,j-1,bi,bj).eq.1.0) THEN
262
263 ! geom_fac = sqrt(rA(i-1,j-1,bi,bj)*recip_rA(i,j-1,bi,bj)*
264 ! & dxF(i,j-1,bi,bj)*recip_dxF(i-1,j-1,bi,bj))
265
266 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
267 & 0.25 * dyg(i,j-1,bi,bj) *
268 & gravity *
269 & (H_streamice(i,j-1,bi,bj)*avg_density(i,j-1,bi,bj)+
270 & H_streamice(i-1,j-1,bi,bj)*avg_density(i-1,j-1,bi,bj)) *
271 & (surf_el_streamice(i,j-1,bi,bj)-
272 & surf_el_streamice(i-1,j-1,bi,bj))
273
274 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
275 & streamice_density * gravity *
276 & streamice_bg_surf_slope_x * .25 * rA(i,j-1,bi,bj) *
277 & (H_streamice(i-1,j-1,bi,bj) + H_streamice(i,j-1,bi,bj))
278
279 ELSE IF (streamice_hmask(i-1,j-1,bi,bj).eq.1.0) THEN
280
281 IF (float_frac_streamice(i-1,j-1,bi,bj) .eq. 1.0) THEN
282
283 #ifdef USE_ALT_RLOW
284 IF (R_low_si(i-1,j-1,bi,bj) .lt. 0. _d 0) then
285 #else
286 IF (R_low(i-1,j-1,bi,bj) .lt. 0. _d 0) then
287 #endif
288 grd_below_sl = 1. _d 0
289 else
290 grd_below_sl = 0. _d 0
291 endif
292
293 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
294 & 0.25 * dyg(i,j-1,bi,bj) *
295 & gravity *
296 & (avg_density(i-1,j-1,bi,bj)*H_streamice(i-1,j-1,bi,bj)**2-
297 #ifdef USE_ALT_RLOW
298 & streamice_density_ocean_avg*grd_below_sl *
299 & R_low_si(i-1,j-1,bi,bj)**2)
300 #else
301 & streamice_density_ocean_avg * grd_below_sl *
302 & R_low(i-1,j-1,bi,bj)**2)
303 #endif
304
305 ELSE
306
307 #ifdef STREAMICE_FIRN_CORRECTION
308
309 if (STREAMICE_apply_firn_correction) then
310 if (H_streamice(i-1,j-1,bi,bj).lt.firn_depth) then
311 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
312 & 0.25 * dyG(i,j-1,bi,bj) *
313 & streamice_density_firn * gravity *
314 & (1-streamice_density_firn*i_rhow) *
315 & H_streamice(i-1,j-1,bi,bj)**2
316 else
317 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
318 & 0.25 * dyG(i,j-1,bi,bj) * gravity * (
319 & streamice_density_firn * firn_depth**2 +
320 & (h_streamice(i-1,j-1,bi,bj)-firn_depth) *
321 & (streamice_density_firn*firn_depth+streamice_density*
322 & (h_streamice(i-1,j-1,bi,bj)-streamice_firn_correction))-
323 & streamice_density**2*i_rhow*
324 & (h_streamice(i-1,j-1,bi,bj)-streamice_firn_correction)**2
325 & )
326 endif
327 else
328
329 #endif
330 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
331 & 0.25 * dyG(i,j-1,bi,bj) *
332 & streamice_density * gravity *
333 & (1-streamice_density*i_rhow) *
334 & H_streamice(i-1,j-1,bi,bj)**2
335 #ifdef STREAMICE_FIRN_CORRECTION
336 endif
337 #endif
338
339
340
341 ENDIF
342
343 ELSE IF (streamice_hmask(i,j-1,bi,bj).eq.1.0) THEN
344
345 IF (float_frac_streamice(i,j-1,bi,bj) .eq. 1.0) THEN
346
347 #ifdef USE_ALT_RLOW
348 IF (R_low_si(i,j-1,bi,bj) .lt. 0. _d 0) then
349 #else
350 IF (R_low(i,j-1,bi,bj) .lt. 0. _d 0) then
351 #endif
352 grd_below_sl = 1. _d 0
353 else
354 grd_below_sl = 0. _d 0
355 endif
356
357 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
358 & 0.25 * dyg(i,j-1,bi,bj) *
359 & gravity *
360 & (avg_density(i,j-1,bi,bj) * H_streamice(i,j-1,bi,bj)**2 -
361 #ifdef USE_ALT_RLOW
362 & streamice_density_ocean_avg * grd_below_sl *
363 & R_low_si(i,j-1,bi,bj)**2)
364 #else
365 & streamice_density_ocean_avg * grd_below_sl *
366 & R_low(i,j-1,bi,bj)**2)
367 #endif
368
369 ELSE
370
371 #ifdef STREAMICE_FIRN_CORRECTION
372 if (STREAMICE_apply_firn_correction) then
373 if (H_streamice(i,j-1,bi,bj).lt.firn_depth) then
374 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
375 & 0.25 * dyG(i,j-1,bi,bj) *
376 & streamice_density_firn * gravity *
377 & (1-streamice_density_firn*i_rhow) *
378 & H_streamice(i,j-1,bi,bj)**2
379 else
380 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
381 & 0.25 * dyG(i,j-1,bi,bj) * gravity * (
382 & streamice_density_firn * firn_depth**2 +
383 & (h_streamice(i,j-1,bi,bj)-firn_depth) *
384 & (streamice_density_firn*firn_depth+streamice_density*
385 & (h_streamice(i,j-1,bi,bj)-streamice_firn_correction))-
386 & streamice_density**2*i_rhow*
387 & (h_streamice(i,j-1,bi,bj)-streamice_firn_correction)**2
388 & )
389 endif
390 else
391
392 #endif
393 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
394 & 0.25 * dyG(i,j-1,bi,bj) *
395 & streamice_density * gravity *
396 & (1-streamice_density*i_rhow) *
397 & H_streamice(i,j-1,bi,bj)**2
398 #ifdef STREAMICE_FIRN_CORRECTION
399 endif
400 #endif
401
402
403 END IF
404 END IF
405 END IF ! if umask==1
406
407 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
408
409 IF (streamice_vmask(i,j,bi,bj).eq.1.0) THEN
410
411 IF (streamice_hmask(i,j-1,bi,bj).eq.1.0.AND.
412 & streamice_hmask(i,j,bi,bj).eq.1.0) THEN
413
414 ! geom_fac = sqrt(rA(i,j-1,bi,bj)*recip_rA(i,j,bi,bj)*
415 ! & dxF(i,j,bi,bj)*recip_dyF(i,j-1,bi,bj))
416
417 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
418 & 0.25 * dxG(i,j,bi,bj) *
419 & gravity *
420 & (H_streamice(i,j,bi,bj)*avg_density(i,j,bi,bj)+
421 & H_streamice(i,j-1,bi,bj)*avg_density(i,j-1,bi,bj))*
422 & (surf_el_streamice(i,j,bi,bj)-
423 & surf_el_streamice(i,j-1,bi,bj))
424
425 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
426 & streamice_density * gravity *
427 & streamice_bg_surf_slope_y * .25 * rA(i,j,bi,bj) *
428 & (H_streamice(i,j-1,bi,bj) + H_streamice(i,j,bi,bj))
429
430 ELSE IF (streamice_hmask(i,j-1,bi,bj).eq.1.0) THEN
431
432 IF (float_frac_streamice(i,j-1,bi,bj) .eq. 1.0) THEN
433
434 #ifdef USE_ALT_RLOW
435 IF (R_low_si(i,j-1,bi,bj) .lt. 0. _d 0) then
436 #else
437 IF (R_low(i,j-1,bi,bj) .lt. 0. _d 0) then
438 #endif
439 grd_below_sl = 1. _d 0
440 else
441 grd_below_sl = 0. _d 0
442 endif
443
444 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
445 & 0.25 * dxG(i,j,bi,bj) *
446 & gravity *
447 & (avg_density(i,j-1,bi,bj) * H_streamice(i,j-1,bi,bj)**2 -
448 #ifdef USE_ALT_RLOW
449 & streamice_density_ocean_avg * grd_below_sl *
450 & R_low_si(i,j-1,bi,bj)**2)
451 #else
452 & streamice_density_ocean_avg * grd_below_sl *
453 & R_low(i,j-1,bi,bj)**2)
454 #endif
455
456 ELSE
457
458 #ifdef STREAMICE_FIRN_CORRECTION
459 if (STREAMICE_apply_firn_correction) then
460 if (H_streamice(i,j-1,bi,bj).lt.firn_depth) then
461 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
462 & 0.25 * dxG(i,j,bi,bj) *
463 & streamice_density_firn * gravity *
464 & (1-streamice_density_firn*i_rhow) *
465 & H_streamice(i,j-1,bi,bj)**2
466 else
467 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
468 & 0.25 * dxG(i,j,bi,bj) * gravity * (
469 & streamice_density_firn * firn_depth**2 +
470 & (h_streamice(i,j-1,bi,bj)-firn_depth) *
471 & (streamice_density_firn*firn_depth+streamice_density*
472 & (h_streamice(i,j-1,bi,bj)-streamice_firn_correction))-
473 & streamice_density**2*i_rhow*
474 & (h_streamice(i,j-1,bi,bj)-streamice_firn_correction)**2
475 & )
476 endif
477 else
478
479 #endif
480 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
481 & 0.25 * dxG(i,j,bi,bj) *
482 & streamice_density * gravity *
483 & (1-streamice_density*i_rhow) *
484 & H_streamice(i,j-1,bi,bj)**2
485 #ifdef STREAMICE_FIRN_CORRECTION
486 endif
487 #endif
488
489
490 ENDIF
491
492 ELSE IF (streamice_hmask(i,j,bi,bj).eq.1.0) THEN
493
494 IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) THEN
495
496 #ifdef USE_ALT_RLOW
497 IF (R_low_si(i,j,bi,bj) .lt. 0. _d 0) then
498 #else
499 IF (R_low(i,j,bi,bj) .lt. 0. _d 0) then
500 #endif
501 grd_below_sl = 1. _d 0
502 else
503 grd_below_sl = 0. _d 0
504 endif
505
506 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
507 & 0.25 * dxG(i,j,bi,bj) *
508 & gravity *
509 & (avg_density(i,j,bi,bj) * H_streamice(i,j,bi,bj)**2 -
510 #ifdef USE_ALT_RLOW
511 & streamice_density_ocean_avg * grd_below_sl *
512 & R_low_si(i,j,bi,bj)**2)
513 #else
514 & streamice_density_ocean_avg * grd_below_sl *
515 & R_low(i,j,bi,bj)**2)
516 #endif
517
518 ELSE
519
520
521 #ifdef STREAMICE_FIRN_CORRECTION
522 if (STREAMICE_apply_firn_correction) then
523 if (H_streamice(i,j,bi,bj).lt.firn_depth) then
524 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
525 & 0.25 * dxG(i,j,bi,bj) *
526 & streamice_density_firn * gravity *
527 & (1-streamice_density_firn*i_rhow) *
528 & H_streamice(i,j,bi,bj)**2
529 else
530 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
531 & 0.25 * dxG(i,j,bi,bj) * gravity * (
532 & streamice_density_firn * firn_depth**2 +
533 & (h_streamice(i,j,bi,bj)-firn_depth) *
534 & (streamice_density_firn*firn_depth+streamice_density*
535 & (h_streamice(i,j,bi,bj)-streamice_firn_correction))-
536 & streamice_density**2*i_rhow*
537 & (h_streamice(i,j,bi,bj)-streamice_firn_correction)**2
538 & )
539 endif
540 else
541
542 #endif
543 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
544 & 0.25 * dxG(i,j,bi,bj) *
545 & streamice_density * gravity *
546 & (1-streamice_density*i_rhow) *
547 & H_streamice(i,j,bi,bj)**2
548 #ifdef STREAMICE_FIRN_CORRECTION
549 endif
550 #endif
551
552
553 END IF
554 END IF
555
556 ! cells to left
557
558 IF (streamice_hmask(i-1,j-1,bi,bj).eq.1.0.AND.
559 & streamice_hmask(i-1,j,bi,bj).eq.1.0) THEN
560
561 ! geom_fac = sqrt(rA(i-1,j-1,bi,bj)*recip_rA(i-1,j,bi,bj)*
562 ! & dxF(i-1,j,bi,bj)*recip_dxF(i-1,j-1,bi,bj))
563
564 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
565 & 0.25 * dxG(i-1,j,bi,bj) *
566 & gravity *
567 & (H_streamice(i-1,j,bi,bj)*avg_density(i-1,j,bi,bj)+
568 & H_streamice(i-1,j-1,bi,bj)*avg_density(i-1,j-1,bi,bj))*
569 & (surf_el_streamice(i-1,j,bi,bj)-
570 & surf_el_streamice(i-1,j-1,bi,bj))
571
572 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
573 & streamice_density * gravity *
574 & streamice_bg_surf_slope_y * .25 * rA(i-1,j,bi,bj) *
575 & (H_streamice(i-1,j-1,bi,bj) + H_streamice(i-1,j,bi,bj))
576
577
578 ELSE IF (streamice_hmask(i-1,j-1,bi,bj).eq.1.0) THEN
579
580 IF (float_frac_streamice(i-1,j-1,bi,bj) .eq. 1.0) THEN
581
582 #ifdef USE_ALT_RLOW
583 IF (R_low_si(i-1,j-1,bi,bj) .lt. 0. _d 0) then
584 #else
585 IF (R_low(i-1,j-1,bi,bj) .lt. 0. _d 0) then
586 #endif
587 grd_below_sl = 1. _d 0
588 else
589 grd_below_sl = 0. _d 0
590 endif
591
592 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
593 & 0.25 * dxG(i-1,j,bi,bj) *
594 & gravity *
595 & (avg_density(i-1,j-1,bi,bj)*H_streamice(i-1,j-1,bi,bj)**2-
596 #ifdef USE_ALT_RLOW
597 & streamice_density_ocean_avg*grd_below_sl *
598 & R_low_si(i-1,j-1,bi,bj)**2)
599 #else
600 & streamice_density_ocean_avg * grd_below_sl *
601 & R_low(i-1,j-1,bi,bj)**2)
602 #endif
603
604 ELSE
605
606 #ifdef STREAMICE_FIRN_CORRECTION
607 if (STREAMICE_apply_firn_correction) then
608 if (H_streamice(i-1,j-1,bi,bj).lt.firn_depth) then
609 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
610 & 0.25 * dxG(i-1,j,bi,bj) *
611 & streamice_density_firn * gravity *
612 & (1-streamice_density_firn*i_rhow) *
613 & H_streamice(i-1,j-1,bi,bj)**2
614 else
615 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
616 & 0.25 * dxG(i-1,j,bi,bj) * gravity * (
617 & streamice_density_firn * firn_depth**2 +
618 & (h_streamice(i-1,j-1,bi,bj)-firn_depth) *
619 & (streamice_density_firn*firn_depth+streamice_density*
620 & (h_streamice(i-1,j-1,bi,bj)-streamice_firn_correction))-
621 & streamice_density**2*i_rhow*
622 & (h_streamice(i-1,j-1,bi,bj)-streamice_firn_correction)**2
623 & )
624 endif
625 else
626 #endif
627 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
628 & 0.25 * dxG(i-1,j,bi,bj) *
629 & streamice_density * gravity *
630 & (1-streamice_density*i_rhow) *
631 & H_streamice(i-1,j-1,bi,bj)**2
632 #ifdef STREAMICE_FIRN_CORRECTION
633 endif
634 #endif
635
636
637
638 ENDIF
639
640 ELSE IF (streamice_hmask(i-1,j,bi,bj).eq.1.0) THEN
641
642 IF (float_frac_streamice(i-1,j,bi,bj) .eq. 1.0) THEN
643
644 #ifdef USE_ALT_RLOW
645 IF (R_low_si(i-1,j,bi,bj) .lt. 0. _d 0) then
646 #else
647 IF (R_low(i-1,j,bi,bj) .lt. 0. _d 0) then
648 #endif
649 grd_below_sl = 1. _d 0
650 else
651 grd_below_sl = 0. _d 0
652 endif
653
654 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
655 & 0.25 * dxG(i-1,j,bi,bj) *
656 & gravity *
657 & (avg_density(i-1,j,bi,bj) * H_streamice(i-1,j,bi,bj)**2 -
658 #ifdef USE_ALT_RLOW
659 & streamice_density_ocean_avg * grd_below_sl *
660 & R_low_si(i-1,j,bi,bj)**2)
661 #else
662 & streamice_density_ocean_avg * grd_below_sl *
663 & R_low(i-1,j,bi,bj)**2)
664 #endif
665
666 ELSE
667
668 #ifdef STREAMICE_FIRN_CORRECTION
669 if (STREAMICE_apply_firn_correction) then
670 if (H_streamice(i-1,j,bi,bj).lt.firn_depth) then
671 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
672 & 0.25 * dxG(i-1,j,bi,bj) *
673 & streamice_density_firn * gravity *
674 & (1-streamice_density_firn*i_rhow) *
675 & H_streamice(i-1,j,bi,bj)**2
676 else
677 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
678 & 0.25 * dxG(i-1,j,bi,bj) * gravity * (
679 & streamice_density_firn * firn_depth**2 +
680 & (h_streamice(i-1,j,bi,bj)-firn_depth) *
681 & (streamice_density_firn*firn_depth+streamice_density*
682 & (h_streamice(i-1,j,bi,bj)-streamice_firn_correction))-
683 & streamice_density**2*i_rhow*
684 & (h_streamice(i-1,j,bi,bj)-streamice_firn_correction)**2
685 & )
686 endif
687 else
688 #endif
689 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
690 & 0.25 * dxG(i-1,j,bi,bj) *
691 & streamice_density * gravity *
692 & (1-streamice_density*i_rhow) *
693 & H_streamice(i-1,j,bi,bj)**2
694 #ifdef STREAMICE_FIRN_CORRECTION
695 endif
696 #endif
697
698
699 END IF
700 END IF
701 END IF ! if vmask ==1
702
703 ENDDO
704 ENDDO
705 ENDDO
706 ENDDO
707
708 ! taudx_SI (1,1,1,1) = taudx_SI (1,1,1,1) +
709 ! & streamice_v_normal_pert (1,1,1,1)
710 ! call write_fld_xy_rl("TAUDX_SI","",taudx_si,0,mythid)
711 ! call write_fld_xy_rl("TAUDY_SI","",taudy_si,0,mythid)
712
713
714
715 #endif
716 RETURN
717 END
718
719

  ViewVC Help
Powered by ViewVC 1.1.22