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

Annotation of /MITgcm_contrib/dgoldberg/streamice/streamice_driving_stress.F

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


Revision 1.2 - (hide annotations) (download)
Tue Sep 18 17:06:48 2012 UTC (12 years, 10 months ago) by dgoldberg
Branch: MAIN
Changes since 1.1: +11 -7 lines
changes for periodic boundary conds and hybrid stress balance

1 dgoldberg 1.2 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_driving_stress.F,v 1.1 2012/03/29 15:59:21 heimbach Exp $
2 heimbach 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     _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 dgoldberg 1.2 IF (Gi.eq.1.AND..NOT.STREAMICE_EW_periodic) THEN
78 heimbach 1.1 IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0) THEN
79     sx = (surf_el_streamice(i+1,j,bi,bj)-
80     & surf_el_streamice(i,j,bi,bj))/dxC(i+1,j,bi,bj)
81     ELSE
82     sx = 0. _d 0
83     ENDIF
84 dgoldberg 1.2 ELSEIF (Gi.eq.Nx.AND..NOT.STREAMICE_EW_periodic) THEN
85 heimbach 1.1 IF (STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN
86     sx = (surf_el_streamice(i,j,bi,bj)-
87     & surf_el_streamice(i-1,j,bi,bj))/dxC(i,j,bi,bj)
88     ELSE
89     sx = 0. _d 0
90     ENDIF
91     ELSE
92     IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0) THEN
93     diffx = diffx + dxC(i+1,j,bi,bj)
94     sx = surf_el_streamice(i+1,j,bi,bj)
95     ELSE
96     sx = surf_el_streamice(i,j,bi,bj)
97     ENDIF
98     IF (STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN
99     diffx = diffx + dxC(i,j,bi,bj)
100     sx = sx - surf_el_streamice(i-1,j,bi,bj)
101     ELSE
102     sx = sx - surf_el_streamice(i,j,bi,bj)
103     ENDIF
104     IF (diffx .gt. 0. _d 0) THEN
105     sx = sx / diffx
106     ELSE
107     sx = 0. _d 0
108     ENDIF
109     ENDIF
110    
111    
112    
113 dgoldberg 1.2 IF (Gj.eq.1.AND..NOT.STREAMICE_NS_periodic) THEN
114 heimbach 1.1 IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0) THEN
115     sy = (surf_el_streamice(i,j+1,bi,bj)-
116     & surf_el_streamice(i,j,bi,bj))/dyC(i,j+1,bi,bj)
117     ELSE
118     sy = 0. _d 0
119     ENDIF
120 dgoldberg 1.2 ELSEIF (Gj.eq.Ny.AND..NOT.STREAMICE_NS_periodic) THEN
121 heimbach 1.1 IF (STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN
122     sy = (surf_el_streamice(i,j,bi,bj)-
123     & surf_el_streamice(i,j-1,bi,bj))/dyC(i,j,bi,bj)
124     ELSE
125     sy = 0. _d 0
126     ENDIF
127     ELSE
128     IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0) THEN
129    
130     diffy = diffy + dyC(i,j+1,bi,bj)
131     sy = surf_el_streamice(i,j+1,bi,bj)
132     ELSE
133     sy = surf_el_streamice(i,j,bi,bj)
134     ENDIF
135     IF (STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN
136     diffy = diffy + dyC(i,j,bi,bj)
137     sy = sy - surf_el_streamice(i,j-1,bi,bj)
138     ELSE
139     sy = sy - surf_el_streamice(i,j,bi,bj)
140     ENDIF
141     IF (diffy .gt. 0. _d 0) THEN
142     sy = sy / diffy
143     ELSE
144     sy = 0. _d 0
145     ENDIF
146     ENDIF
147    
148     DO k=0,1
149     DO l=0,1
150     IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
151     taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
152 dgoldberg 1.2 & 0.25 * streamice_density * gravity *
153     & (streamice_bg_surf_slope_x+sx) *
154 heimbach 1.1 & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
155 dgoldberg 1.2 ! & (streamice_bg_surf_slope_x) *
156     ! & 1000. * rA(i,j,bi,bj)
157 heimbach 1.1 taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
158 dgoldberg 1.2 & 0.25 * streamice_density * gravity *
159     & (streamice_bg_surf_slope_y+sy) *
160 heimbach 1.1 & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
161    
162     ENDIF
163     ENDDO
164     ENDDO
165    
166     IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) then
167     neu_val = .5 * gravity *
168     & (streamice_density * H_streamice (i,j,bi,bj) ** 2 -
169     & streamice_density_ocean_avg * R_low(i,j,bi,bj) ** 2)
170     ELSE
171     neu_val = .5 * gravity *
172     & (1-streamice_density/streamice_density_ocean_avg) *
173     & streamice_density * H_streamice(i,j,bi,bj) ** 2
174     ENDIF
175    
176     IF ((STREAMICE_ufacemask(i,j,bi,bj) .eq. 2)
177     & .OR. (STREAMICE_hmask(i-1,j,bi,bj) .eq. 0)
178     & .OR. (STREAMICE_hmask(i-1,j,bi,bj) .eq. 2) ) THEN ! left face of the cell is at a stress boundary
179     ! the depth-integrated longitudinal stress is equal to the difference of depth-integrated pressure on either side of the face
180     ! on the ice side, it is rho g h^2 / 2
181     ! on the ocean side, it is rhow g (delta OD)^2 / 2
182     ! 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
183     ! ice in the current cell
184    
185     taudx_SI(i,j,bi,bj) = taudx_SI(i,j,bi,bj) -
186     & .5 * dyG(i,j,bi,bj) * neu_val ! note negative sign is due to direction of normal vector
187     taudx_SI(i,j+1,bi,bj) = taudx_SI(i,j+1,bi,bj) -
188     & .5 * dyG(i,j,bi,bj) * neu_val
189     ENDIF
190    
191     IF ((STREAMICE_ufacemask(i+1,j,bi,bj) .eq. 2)
192     & .OR. (STREAMICE_hmask(i+1,j,bi,bj) .eq. 0)
193     & .OR. (STREAMICE_hmask(i+1,j,bi,bj) .eq. 2) ) THEN
194    
195     taudx_SI(i+1,j,bi,bj) = taudx_SI(i+1,j,bi,bj) +
196     & .5 * dyG(i+1,j,bi,bj) * neu_val ! note negative sign is due to direction of normal vector
197     taudx_SI(i+1,j+1,bi,bj) = taudx_SI(i+1,j+1,bi,bj) +
198     & .5 * dyG(i+1,j,bi,bj) * neu_val
199     ENDIF
200    
201     IF ((STREAMICE_vfacemask(i,j,bi,bj) .eq. 2)
202     & .OR. (STREAMICE_hmask(i,j-1,bi,bj) .eq. 0)
203     & .OR. (STREAMICE_hmask(i,j-1,bi,bj) .eq. 2) ) THEN
204    
205     taudy_SI(i,j,bi,bj) = taudy_SI(i,j,bi,bj) -
206     & .5 * dxG(i,j,bi,bj) * neu_val ! note negative sign is due to direction of normal vector
207     taudy_SI(i+1,j,bi,bj) = taudy_SI(i+1,j,bi,bj) -
208     & .5 * dxG(i,j,bi,bj) * neu_val
209     ENDIF
210    
211     IF ((STREAMICE_vfacemask(i,j+1,bi,bj) .eq. 2)
212     & .OR. (STREAMICE_hmask(i,j+1,bi,bj) .eq. 0)
213     & .OR. (STREAMICE_hmask(i,j+1,bi,bj) .eq. 2) ) THEN
214    
215     taudy_SI(i,j+1,bi,bj) = taudy_SI(i,j+1,bi,bj) +
216     & .5 * dxG(i,j+1,bi,bj) * neu_val ! note negative sign is due to direction of normal vector
217     taudy_SI(i+1,j+1,bi,bj) = taudy_SI(i+1,j+1,bi,bj) +
218     & .5 * dxG(i,j+1,bi,bj) * neu_val
219     ENDIF
220    
221     ENDIF
222     ENDDO
223     ENDDO
224     ENDDO
225     ENDDO
226    
227    
228    
229     #endif
230     RETURN
231     END
232    
233    

  ViewVC Help
Powered by ViewVC 1.1.22