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

1 heimbach 1.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