/[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.1 - (show annotations) (download)
Thu Mar 29 15:59:21 2012 UTC (13 years, 4 months ago) by heimbach
Branch: MAIN
Initial check-in of Dan Goldberg's streamice package

1 C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_init_varia.F,v 1.6 2011/06/29 16:24:10 dng 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) 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) 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) 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) 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 * sx *
153 & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
154 taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
155 & 0.25 * streamice_density * gravity * sy *
156 & H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
157
158 ENDIF
159 ENDDO
160 ENDDO
161
162 IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) then
163 neu_val = .5 * gravity *
164 & (streamice_density * H_streamice (i,j,bi,bj) ** 2 -
165 & streamice_density_ocean_avg * R_low(i,j,bi,bj) ** 2)
166 ELSE
167 neu_val = .5 * gravity *
168 & (1-streamice_density/streamice_density_ocean_avg) *
169 & streamice_density * H_streamice(i,j,bi,bj) ** 2
170 ENDIF
171
172 IF ((STREAMICE_ufacemask(i,j,bi,bj) .eq. 2)
173 & .OR. (STREAMICE_hmask(i-1,j,bi,bj) .eq. 0)
174 & .OR. (STREAMICE_hmask(i-1,j,bi,bj) .eq. 2) ) THEN ! left face of the cell is at a stress boundary
175 ! the depth-integrated longitudinal stress is equal to the difference of depth-integrated pressure on either side of the face
176 ! on the ice side, it is rho g h^2 / 2
177 ! on the ocean side, it is rhow g (delta OD)^2 / 2
178 ! 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
179 ! ice in the current cell
180
181 taudx_SI(i,j,bi,bj) = taudx_SI(i,j,bi,bj) -
182 & .5 * dyG(i,j,bi,bj) * neu_val ! note negative sign is due to direction of normal vector
183 taudx_SI(i,j+1,bi,bj) = taudx_SI(i,j+1,bi,bj) -
184 & .5 * dyG(i,j,bi,bj) * neu_val
185 ENDIF
186
187 IF ((STREAMICE_ufacemask(i+1,j,bi,bj) .eq. 2)
188 & .OR. (STREAMICE_hmask(i+1,j,bi,bj) .eq. 0)
189 & .OR. (STREAMICE_hmask(i+1,j,bi,bj) .eq. 2) ) THEN
190
191 taudx_SI(i+1,j,bi,bj) = taudx_SI(i+1,j,bi,bj) +
192 & .5 * dyG(i+1,j,bi,bj) * neu_val ! note negative sign is due to direction of normal vector
193 taudx_SI(i+1,j+1,bi,bj) = taudx_SI(i+1,j+1,bi,bj) +
194 & .5 * dyG(i+1,j,bi,bj) * neu_val
195 ENDIF
196
197 IF ((STREAMICE_vfacemask(i,j,bi,bj) .eq. 2)
198 & .OR. (STREAMICE_hmask(i,j-1,bi,bj) .eq. 0)
199 & .OR. (STREAMICE_hmask(i,j-1,bi,bj) .eq. 2) ) THEN
200
201 taudy_SI(i,j,bi,bj) = taudy_SI(i,j,bi,bj) -
202 & .5 * dxG(i,j,bi,bj) * neu_val ! note negative sign is due to direction of normal vector
203 taudy_SI(i+1,j,bi,bj) = taudy_SI(i+1,j,bi,bj) -
204 & .5 * dxG(i,j,bi,bj) * neu_val
205 ENDIF
206
207 IF ((STREAMICE_vfacemask(i,j+1,bi,bj) .eq. 2)
208 & .OR. (STREAMICE_hmask(i,j+1,bi,bj) .eq. 0)
209 & .OR. (STREAMICE_hmask(i,j+1,bi,bj) .eq. 2) ) THEN
210
211 taudy_SI(i,j+1,bi,bj) = taudy_SI(i,j+1,bi,bj) +
212 & .5 * dxG(i,j+1,bi,bj) * neu_val ! note negative sign is due to direction of normal vector
213 taudy_SI(i+1,j+1,bi,bj) = taudy_SI(i+1,j+1,bi,bj) +
214 & .5 * dxG(i,j+1,bi,bj) * neu_val
215 ENDIF
216
217 ENDIF
218 ENDDO
219 ENDDO
220 ENDDO
221 ENDDO
222
223
224
225 #endif
226 RETURN
227 END
228
229

  ViewVC Help
Powered by ViewVC 1.1.22