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

Contents 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 - (show 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 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_driving_stress.F,v 1.1 2012/03/29 15:59:21 heimbach 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 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 IF (Gi.eq.1.AND..NOT.STREAMICE_EW_periodic) THEN
78 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 ELSEIF (Gi.eq.Nx.AND..NOT.STREAMICE_EW_periodic) THEN
85 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 IF (Gj.eq.1.AND..NOT.STREAMICE_NS_periodic) THEN
114 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 ELSEIF (Gj.eq.Ny.AND..NOT.STREAMICE_NS_periodic) THEN
121 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 & 0.25 * streamice_density * gravity *
153 & (streamice_bg_surf_slope_x+sx) *
154 & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
155 ! & (streamice_bg_surf_slope_x) *
156 ! & 1000. * rA(i,j,bi,bj)
157 taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
158 & 0.25 * streamice_density * gravity *
159 & (streamice_bg_surf_slope_y+sy) *
160 & 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