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

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

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


Revision 1.10 - (hide annotations) (download)
Sun Oct 18 14:37:49 2015 UTC (8 years, 7 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 dgoldberg 1.10 C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_driving_stress.F,v 1.9 2015/07/04 10:54:24 dgoldberg Exp $
2 dgoldberg 1.1 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 dgoldberg 1.6 _RL grd_below_sl
43 dgoldberg 1.2 _RL sx, sy, diffx, diffy, geom_fac
44 dgoldberg 1.7 _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 dgoldberg 1.1
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 dgoldberg 1.7 #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 dgoldberg 1.1
65     DO bj = myByLo(myThid), myByHi(myThid)
66     DO bi = myBxLo(myThid), myBxHi(myThid)
67 dgoldberg 1.4 DO j=1-OLy+1,sNy+OLy-1
68     DO i=1-OLx+1,sNx+OLx-1
69 dgoldberg 1.1 taudx_SI(i,j,bi,bj) = 0. _d 0
70     taudy_SI(i,j,bi,bj) = 0. _d 0
71 dgoldberg 1.7 #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 dgoldberg 1.1 ENDDO
89     ENDDO
90     ENDDO
91     ENDDO
92    
93     DO bj = myByLo(myThid), myByHi(myThid)
94     DO bi = myBxLo(myThid), myBxHi(myThid)
95    
96 dgoldberg 1.2 DO i=1,sNx
97     DO j=1,sNy
98 dgoldberg 1.1
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 dgoldberg 1.2 IF (streamice_umask(i,j,bi,bj).eq.1.0) THEN
108 dgoldberg 1.1
109 dgoldberg 1.2 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 dgoldberg 1.5 ! 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 dgoldberg 1.2
115 dgoldberg 1.7
116 dgoldberg 1.2 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
117 dgoldberg 1.5 & 0.25 * dyG(i,j,bi,bj) *
118 dgoldberg 1.7 & 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 dgoldberg 1.2 & (surf_el_streamice(i,j,bi,bj)-
122 dgoldberg 1.5 & surf_el_streamice(i-1,j,bi,bj))
123 dgoldberg 1.2
124 dgoldberg 1.7 CCC
125 dgoldberg 1.3 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 dgoldberg 1.7 CCC
130 dgoldberg 1.3
131 dgoldberg 1.2 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 dgoldberg 1.6 #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 dgoldberg 1.2 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
146 dgoldberg 1.5 & 0.25 * dyG(i,j,bi,bj) *
147 dgoldberg 1.2 & gravity *
148 dgoldberg 1.7 & (avg_density(i-1,j,bi,bj) * H_streamice(i-1,j,bi,bj)**2-
149 dgoldberg 1.2 #ifdef USE_ALT_RLOW
150 dgoldberg 1.6 & streamice_density_ocean_avg * grd_below_sl *
151     & R_low_si(i-1,j,bi,bj)**2)
152 dgoldberg 1.2 #else
153 dgoldberg 1.6 & streamice_density_ocean_avg * grd_below_sl *
154     & R_low(i-1,j,bi,bj)**2)
155 dgoldberg 1.2 #endif
156    
157     ELSE
158 dgoldberg 1.1
159 dgoldberg 1.7 #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 dgoldberg 1.8 & streamice_density_firn * firn_depth**2 +
172 dgoldberg 1.7 & (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 dgoldberg 1.2 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
183 dgoldberg 1.5 & 0.25 * dyG(i,j,bi,bj) *
184 dgoldberg 1.2 & streamice_density * gravity *
185 dgoldberg 1.7 & (1-streamice_density*i_rhow) *
186 dgoldberg 1.2 & H_streamice(i-1,j,bi,bj)**2
187 dgoldberg 1.1
188 dgoldberg 1.7 #ifdef STREAMICE_FIRN_CORRECTION
189     endif
190     #endif
191    
192 dgoldberg 1.2 ENDIF
193 dgoldberg 1.1
194 dgoldberg 1.2 ELSE IF (streamice_hmask(i,j,bi,bj).eq.1.0) THEN
195 dgoldberg 1.1
196 dgoldberg 1.2 IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) THEN
197    
198 dgoldberg 1.6 #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 dgoldberg 1.2 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
209 dgoldberg 1.5 & 0.25 * dyG(i,j,bi,bj) *
210 dgoldberg 1.2 & gravity *
211 dgoldberg 1.7 & (avg_density(i,j,bi,bj) * H_streamice(i,j,bi,bj)**2 -
212 dgoldberg 1.2 #ifdef USE_ALT_RLOW
213 dgoldberg 1.6 & streamice_density_ocean_avg * grd_below_sl *
214     & R_low_si(i,j,bi,bj)**2)
215 dgoldberg 1.2 #else
216 dgoldberg 1.6 & streamice_density_ocean_avg * grd_below_sl *
217     & R_low(i,j,bi,bj)**2)
218 dgoldberg 1.2 #endif
219 dgoldberg 1.1
220     ELSE
221    
222 dgoldberg 1.7 #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 dgoldberg 1.8 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
227 dgoldberg 1.7 & 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 dgoldberg 1.8 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
233 dgoldberg 1.7 & 0.25 * dyG(i,j,bi,bj) * gravity * (
234 dgoldberg 1.8 & streamice_density_firn * firn_depth**2 +
235 dgoldberg 1.7 & (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 dgoldberg 1.8 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
246 dgoldberg 1.5 & 0.25 * dyG(i,j,bi,bj) *
247 dgoldberg 1.2 & streamice_density * gravity *
248 dgoldberg 1.7 & (1-streamice_density*i_rhow) *
249     & H_streamice(i,j,bi,bj)**2
250     #ifdef STREAMICE_FIRN_CORRECTION
251     endif
252     #endif
253    
254 dgoldberg 1.2
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 dgoldberg 1.5 ! 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 dgoldberg 1.2
266     taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
267 dgoldberg 1.5 & 0.25 * dyg(i,j-1,bi,bj) *
268 dgoldberg 1.7 & 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 dgoldberg 1.2 & (surf_el_streamice(i,j-1,bi,bj)-
272 dgoldberg 1.5 & surf_el_streamice(i-1,j-1,bi,bj))
273 dgoldberg 1.2
274 dgoldberg 1.3 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
275     & streamice_density * gravity *
276 dgoldberg 1.5 & streamice_bg_surf_slope_x * .25 * rA(i,j-1,bi,bj) *
277 dgoldberg 1.3 & (H_streamice(i-1,j-1,bi,bj) + H_streamice(i,j-1,bi,bj))
278    
279 dgoldberg 1.2 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 dgoldberg 1.6 #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 dgoldberg 1.2 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
294 dgoldberg 1.5 & 0.25 * dyg(i,j-1,bi,bj) *
295 dgoldberg 1.2 & gravity *
296 dgoldberg 1.7 & (avg_density(i-1,j-1,bi,bj)*H_streamice(i-1,j-1,bi,bj)**2-
297 dgoldberg 1.2 #ifdef USE_ALT_RLOW
298 dgoldberg 1.6 & streamice_density_ocean_avg*grd_below_sl *
299     & R_low_si(i-1,j-1,bi,bj)**2)
300 dgoldberg 1.2 #else
301 dgoldberg 1.6 & streamice_density_ocean_avg * grd_below_sl *
302     & R_low(i-1,j-1,bi,bj)**2)
303 dgoldberg 1.2 #endif
304    
305     ELSE
306    
307 dgoldberg 1.7 #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 dgoldberg 1.8 & streamice_density_firn * firn_depth**2 +
320 dgoldberg 1.7 & (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 dgoldberg 1.2 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) +
331 dgoldberg 1.7 & 0.25 * dyG(i,j-1,bi,bj) *
332 dgoldberg 1.2 & streamice_density * gravity *
333 dgoldberg 1.7 & (1-streamice_density*i_rhow) *
334 dgoldberg 1.2 & H_streamice(i-1,j-1,bi,bj)**2
335 dgoldberg 1.7 #ifdef STREAMICE_FIRN_CORRECTION
336     endif
337     #endif
338    
339    
340 dgoldberg 1.1
341     ENDIF
342    
343 dgoldberg 1.2 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 dgoldberg 1.1
347 dgoldberg 1.6 #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 dgoldberg 1.2 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
358 dgoldberg 1.5 & 0.25 * dyg(i,j-1,bi,bj) *
359 dgoldberg 1.2 & gravity *
360 dgoldberg 1.7 & (avg_density(i,j-1,bi,bj) * H_streamice(i,j-1,bi,bj)**2 -
361 dgoldberg 1.2 #ifdef USE_ALT_RLOW
362 dgoldberg 1.6 & streamice_density_ocean_avg * grd_below_sl *
363     & R_low_si(i,j-1,bi,bj)**2)
364 dgoldberg 1.2 #else
365 dgoldberg 1.6 & streamice_density_ocean_avg * grd_below_sl *
366     & R_low(i,j-1,bi,bj)**2)
367 dgoldberg 1.2 #endif
368 dgoldberg 1.1
369 dgoldberg 1.2 ELSE
370 dgoldberg 1.1
371 dgoldberg 1.7 #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 dgoldberg 1.8 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
375 dgoldberg 1.7 & 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 dgoldberg 1.8 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
381 dgoldberg 1.7 & 0.25 * dyG(i,j-1,bi,bj) * gravity * (
382 dgoldberg 1.8 & streamice_density_firn * firn_depth**2 +
383 dgoldberg 1.7 & (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 dgoldberg 1.8 taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) -
394 dgoldberg 1.7 & 0.25 * dyG(i,j-1,bi,bj) *
395 dgoldberg 1.2 & streamice_density * gravity *
396 dgoldberg 1.7 & (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 dgoldberg 1.2
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 dgoldberg 1.5 ! 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 dgoldberg 1.2
417     taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
418 dgoldberg 1.5 & 0.25 * dxG(i,j,bi,bj) *
419 dgoldberg 1.7 & 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 dgoldberg 1.2 & (surf_el_streamice(i,j,bi,bj)-
423 dgoldberg 1.5 & surf_el_streamice(i,j-1,bi,bj))
424 dgoldberg 1.2
425 dgoldberg 1.3 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 dgoldberg 1.2 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 dgoldberg 1.6 #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 dgoldberg 1.2 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
445 dgoldberg 1.5 & 0.25 * dxG(i,j,bi,bj) *
446 dgoldberg 1.2 & gravity *
447 dgoldberg 1.7 & (avg_density(i,j-1,bi,bj) * H_streamice(i,j-1,bi,bj)**2 -
448 dgoldberg 1.2 #ifdef USE_ALT_RLOW
449 dgoldberg 1.6 & streamice_density_ocean_avg * grd_below_sl *
450     & R_low_si(i,j-1,bi,bj)**2)
451 dgoldberg 1.2 #else
452 dgoldberg 1.6 & streamice_density_ocean_avg * grd_below_sl *
453     & R_low(i,j-1,bi,bj)**2)
454 dgoldberg 1.2 #endif
455 dgoldberg 1.1
456     ELSE
457    
458 dgoldberg 1.7 #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 dgoldberg 1.9 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
462 dgoldberg 1.10 & 0.25 * dxG(i,j,bi,bj) *
463 dgoldberg 1.7 & streamice_density_firn * gravity *
464     & (1-streamice_density_firn*i_rhow) *
465     & H_streamice(i,j-1,bi,bj)**2
466     else
467 dgoldberg 1.9 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
468 dgoldberg 1.10 & 0.25 * dxG(i,j,bi,bj) * gravity * (
469 dgoldberg 1.8 & streamice_density_firn * firn_depth**2 +
470 dgoldberg 1.7 & (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 dgoldberg 1.9 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
481 dgoldberg 1.10 & 0.25 * dxG(i,j,bi,bj) *
482 dgoldberg 1.2 & streamice_density * gravity *
483 dgoldberg 1.7 & (1-streamice_density*i_rhow) *
484 dgoldberg 1.2 & H_streamice(i,j-1,bi,bj)**2
485 dgoldberg 1.7 #ifdef STREAMICE_FIRN_CORRECTION
486     endif
487     #endif
488    
489 dgoldberg 1.1
490     ENDIF
491    
492 dgoldberg 1.2 ELSE IF (streamice_hmask(i,j,bi,bj).eq.1.0) THEN
493 dgoldberg 1.1
494 dgoldberg 1.2 IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) THEN
495 dgoldberg 1.1
496 dgoldberg 1.6 #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 dgoldberg 1.2 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
507 dgoldberg 1.5 & 0.25 * dxG(i,j,bi,bj) *
508 dgoldberg 1.2 & gravity *
509 dgoldberg 1.7 & (avg_density(i,j,bi,bj) * H_streamice(i,j,bi,bj)**2 -
510 dgoldberg 1.2 #ifdef USE_ALT_RLOW
511 dgoldberg 1.6 & streamice_density_ocean_avg * grd_below_sl *
512     & R_low_si(i,j,bi,bj)**2)
513 dgoldberg 1.2 #else
514 dgoldberg 1.6 & streamice_density_ocean_avg * grd_below_sl *
515     & R_low(i,j,bi,bj)**2)
516 dgoldberg 1.2 #endif
517 dgoldberg 1.1
518     ELSE
519    
520 dgoldberg 1.7
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 dgoldberg 1.9 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
525 dgoldberg 1.10 & 0.25 * dxG(i,j,bi,bj) *
526 dgoldberg 1.7 & streamice_density_firn * gravity *
527     & (1-streamice_density_firn*i_rhow) *
528     & H_streamice(i,j,bi,bj)**2
529     else
530 dgoldberg 1.9 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
531 dgoldberg 1.10 & 0.25 * dxG(i,j,bi,bj) * gravity * (
532 dgoldberg 1.8 & streamice_density_firn * firn_depth**2 +
533 dgoldberg 1.7 & (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 dgoldberg 1.9 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
544 dgoldberg 1.10 & 0.25 * dxG(i,j,bi,bj) *
545 dgoldberg 1.2 & streamice_density * gravity *
546 dgoldberg 1.7 & (1-streamice_density*i_rhow) *
547     & H_streamice(i,j,bi,bj)**2
548     #ifdef STREAMICE_FIRN_CORRECTION
549     endif
550     #endif
551    
552 dgoldberg 1.2
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 dgoldberg 1.5 ! 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 dgoldberg 1.2
564     taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
565 dgoldberg 1.5 & 0.25 * dxG(i-1,j,bi,bj) *
566 dgoldberg 1.7 & 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 dgoldberg 1.2 & (surf_el_streamice(i-1,j,bi,bj)-
570 dgoldberg 1.5 & surf_el_streamice(i-1,j-1,bi,bj))
571 dgoldberg 1.2
572 dgoldberg 1.3 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
573     & streamice_density * gravity *
574 dgoldberg 1.5 & streamice_bg_surf_slope_y * .25 * rA(i-1,j,bi,bj) *
575 dgoldberg 1.3 & (H_streamice(i-1,j-1,bi,bj) + H_streamice(i-1,j,bi,bj))
576    
577    
578 dgoldberg 1.2 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 dgoldberg 1.6 #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 dgoldberg 1.2 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
593 dgoldberg 1.5 & 0.25 * dxG(i-1,j,bi,bj) *
594 dgoldberg 1.2 & gravity *
595 dgoldberg 1.7 & (avg_density(i-1,j-1,bi,bj)*H_streamice(i-1,j-1,bi,bj)**2-
596 dgoldberg 1.2 #ifdef USE_ALT_RLOW
597 dgoldberg 1.6 & streamice_density_ocean_avg*grd_below_sl *
598     & R_low_si(i-1,j-1,bi,bj)**2)
599 dgoldberg 1.2 #else
600 dgoldberg 1.6 & streamice_density_ocean_avg * grd_below_sl *
601     & R_low(i-1,j-1,bi,bj)**2)
602 dgoldberg 1.2 #endif
603 dgoldberg 1.1
604     ELSE
605    
606 dgoldberg 1.7 #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 dgoldberg 1.9 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
610 dgoldberg 1.10 & 0.25 * dxG(i-1,j,bi,bj) *
611 dgoldberg 1.7 & streamice_density_firn * gravity *
612     & (1-streamice_density_firn*i_rhow) *
613     & H_streamice(i-1,j-1,bi,bj)**2
614     else
615 dgoldberg 1.9 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
616 dgoldberg 1.10 & 0.25 * dxG(i-1,j,bi,bj) * gravity * (
617 dgoldberg 1.8 & streamice_density_firn * firn_depth**2 +
618 dgoldberg 1.7 & (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 dgoldberg 1.9 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) +
628 dgoldberg 1.10 & 0.25 * dxG(i-1,j,bi,bj) *
629 dgoldberg 1.2 & streamice_density * gravity *
630 dgoldberg 1.7 & (1-streamice_density*i_rhow) *
631 dgoldberg 1.2 & H_streamice(i-1,j-1,bi,bj)**2
632 dgoldberg 1.7 #ifdef STREAMICE_FIRN_CORRECTION
633     endif
634     #endif
635    
636    
637 dgoldberg 1.1
638 dgoldberg 1.2 ENDIF
639 dgoldberg 1.1
640 dgoldberg 1.2 ELSE IF (streamice_hmask(i-1,j,bi,bj).eq.1.0) THEN
641 dgoldberg 1.1
642 dgoldberg 1.2 IF (float_frac_streamice(i-1,j,bi,bj) .eq. 1.0) THEN
643 dgoldberg 1.1
644 dgoldberg 1.6 #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 dgoldberg 1.2 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
655 dgoldberg 1.5 & 0.25 * dxG(i-1,j,bi,bj) *
656 dgoldberg 1.2 & gravity *
657 dgoldberg 1.7 & (avg_density(i-1,j,bi,bj) * H_streamice(i-1,j,bi,bj)**2 -
658 dgoldberg 1.1 #ifdef USE_ALT_RLOW
659 dgoldberg 1.6 & streamice_density_ocean_avg * grd_below_sl *
660     & R_low_si(i-1,j,bi,bj)**2)
661 dgoldberg 1.2 #else
662 dgoldberg 1.6 & streamice_density_ocean_avg * grd_below_sl *
663     & R_low(i-1,j,bi,bj)**2)
664 dgoldberg 1.1 #endif
665    
666 dgoldberg 1.2 ELSE
667    
668 dgoldberg 1.7 #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 dgoldberg 1.9 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
672 dgoldberg 1.10 & 0.25 * dxG(i-1,j,bi,bj) *
673 dgoldberg 1.7 & streamice_density_firn * gravity *
674     & (1-streamice_density_firn*i_rhow) *
675     & H_streamice(i-1,j,bi,bj)**2
676     else
677 dgoldberg 1.9 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
678 dgoldberg 1.10 & 0.25 * dxG(i-1,j,bi,bj) * gravity * (
679 dgoldberg 1.8 & streamice_density_firn * firn_depth**2 +
680 dgoldberg 1.7 & (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 dgoldberg 1.9 taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) -
690 dgoldberg 1.10 & 0.25 * dxG(i-1,j,bi,bj) *
691 dgoldberg 1.2 & streamice_density * gravity *
692 dgoldberg 1.7 & (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 dgoldberg 1.2
699     END IF
700     END IF
701     END IF ! if vmask ==1
702    
703 dgoldberg 1.1 ENDDO
704     ENDDO
705     ENDDO
706     ENDDO
707    
708 dgoldberg 1.4 ! taudx_SI (1,1,1,1) = taudx_SI (1,1,1,1) +
709     ! & streamice_v_normal_pert (1,1,1,1)
710 dgoldberg 1.7 ! call write_fld_xy_rl("TAUDX_SI","",taudx_si,0,mythid)
711     ! call write_fld_xy_rl("TAUDY_SI","",taudy_si,0,mythid)
712 dgoldberg 1.4
713 dgoldberg 1.1
714    
715     #endif
716     RETURN
717     END
718    
719    

  ViewVC Help
Powered by ViewVC 1.1.22