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

Annotation 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.2 - (hide annotations) (download)
Fri Sep 28 02:55:40 2012 UTC (12 years, 9 months ago) by heimbach
Branch: MAIN
Changes since 1.1: +2 -2 lines
Bug fixes

1 heimbach 1.2 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_driving_stress_ppm.F,v 1.1 2012/09/27 23:37:39 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_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 heimbach 1.2 & streamice_density_ocean_avg * R_low(i,j,bi,bj) ** 2)
354 dgoldberg 1.1 #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 ! note negative sign is due to direction of normal vector
372     taudx_SI(i,j+1,bi,bj) = taudx_SI(i,j+1,bi,bj) -
373     & .5 * dyG(i,j,bi,bj) * neu_val
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 ! 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
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 ! note negative sign is due to direction of normal vector
392     taudy_SI(i+1,j,bi,bj) = taudy_SI(i+1,j,bi,bj) -
393     & .5 * dxG(i,j,bi,bj) * neu_val
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 ! note negative sign is due to direction of normal vector
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
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