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

Annotation of /MITgcm_contrib/dgoldberg/streamice/streamice_advect_thickness.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    
9     CBOP
10     SUBROUTINE STREAMICE_ADVECT_THICKNESS ( myThid, time_step )
11    
12     C /============================================================\
13     C | SUBROUTINE |
14     C | o |
15     C |============================================================|
16     C | |
17     C \============================================================/
18     IMPLICIT NONE
19    
20     C === Global variables ===
21     #include "SIZE.h"
22     #include "GRID.h"
23     #include "EEPARAMS.h"
24     #include "PARAMS.h"
25     #include "STREAMICE.h"
26     #include "STREAMICE_ADV.h"
27    
28     INTEGER myThid
29     _RL time_step
30    
31     #ifdef ALLOW_STREAMICE
32    
33     INTEGER i, j, bi, bj
34     _RL thick_bd
35     _RL SLOPE_LIMITER
36     _RL sec_per_year, time_step_loc
37     external SLOPE_LIMITER
38    
39     sec_per_year = 365.*86400.
40    
41     time_step_loc = time_step / sec_per_year
42    
43     DO bj=myByLo(myThid),myByHi(myThid)
44     DO bi=myBxLo(myThid),myBxHi(myThid)
45     DO j=1-OLy,sNy+OLy
46     DO i=1-OLx,sNx+OLx
47     hflux_x_SI (i,j,bi,bj) = 0. _d 0
48     hflux_y_SI (i,j,bi,bj) = 0. _d 0
49     hflux_x_SI2 (i,j,bi,bj) = 0. _d 0
50     hflux_y_SI2 (i,j,bi,bj) = 0. _d 0
51     IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
52     h_after_uflux_SI (i,j,bi,bj) =
53     & H_streamice (i,j,bi,bj)
54     ENDIF
55    
56     thick_bd = h_bdry_values_SI (i,j,bi,bj)
57     IF (thick_bd .ne. 0. _d 0) THEN
58     h_after_uflux_SI (i,j,bi,bj) = thick_bd
59     ENDIF
60     ENDDO
61     ENDDO
62     ENDDO
63     ENDDO
64    
65     ! PRINT *, "H in last row ", H_streamice(81,20,1,1)
66    
67     CALL STREAMICE_ADVECT_THICKNESS_X ( myThid,
68     O hflux_x_SI,
69     O h_after_uflux_SI,
70     I time_step_loc )
71    
72     ! PRINT *, "H in last row ", H_streamice(81,20,1,1)
73    
74     DO bj=myByLo(myThid),myByHi(myThid)
75     DO bi=myBxLo(myThid),myBxHi(myThid)
76     DO j=1-OLy,sNy+OLy
77     DO i=1-OLx,sNx+OLx
78     h_after_vflux_SI (i,j,bi,bj) =
79     & h_after_uflux_SI (i,j,bi,bj)
80     ENDDO
81     ENDDO
82     ENDDO
83     ENDDO
84    
85     CALL STREAMICE_ADVECT_THICKNESS_Y ( myThid,
86     O hflux_y_SI,
87     O h_after_vflux_SI,
88     I time_step_loc )
89    
90     ! PRINT *, "H in last row ", H_streamice(81,20,1,1)
91    
92     DO bj=myByLo(myThid),myByHi(myThid)
93     DO bi=myBxLo(myThid),myBxHi(myThid)
94     DO j=1-OLy,sNy+OLy
95     DO i=1-OLx,sNx+OLx
96     IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
97     H_streamice (i,j,bi,bj) =
98     & h_after_vflux_SI (i,j,bi,bj)
99     ENDIF
100     ENDDO
101     ENDDO
102     ENDDO
103     ENDDO
104    
105     ! PRINT *, "H in last row ", H_streamice(81,20,1,1)
106    
107     CALL STREAMICE_ADV_FRONT ( myThid, time_step_loc )
108    
109     ! PRINT *, "H in last row ", H_streamice(81,20,1,1)
110    
111     _EXCH_XY_RL( H_streamice, myThid )
112     _EXCH_XY_RL( area_shelf_streamice, myThid )
113     _EXCH_XY_RL( STREAMICE_hmask, myThid )
114    
115     PRINT *, "END STREAMICE_ADVECT_THICKNESS"
116    
117     #endif
118     RETURN
119     END SUBROUTINE STREAMICE_ADVECT_THICKNESS
120    
121     ! NEED TO ADD SOME SORT OF CHECK THAT THE HALOS ARE LARGE ENOUGH
122    
123     SUBROUTINE STREAMICE_ADVECT_THICKNESS_X ( myThid ,
124     O hflux_x ,
125     O h ,
126     I time_step )
127    
128     IMPLICIT NONE
129    
130     C O hflux_x ! flux per unit width across face
131     C O h
132     C I time_step
133    
134     C === Global variables ===
135     #include "SIZE.h"
136     #include "GRID.h"
137     #include "EEPARAMS.h"
138     #include "PARAMS.h"
139     #include "STREAMICE.h"
140    
141     INTEGER myThid
142     _RL hflux_x (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
143     _RL h (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
144     _RL time_step
145    
146     #ifdef ALLOW_STREAMICE
147    
148     C LOCAL VARIABLES
149    
150     INTEGER i, j, bi, bj, Gi, Gj, k
151     _RL uface, phi
152     _RL stencil (-1:1)
153     LOGICAL H0_valid ! there are valid cells to calculate a
154     ! slope-limited 2nd order flux
155     _RL SLOPE_LIMITER
156     _RL total_vol_out
157     external SLOPE_LIMITER
158    
159     total_vol_out = 0.0
160    
161     DO bj=myByLo(myThid),myByHi(myThid)
162     DO bi=myBxLo(myThid),myBxHi(myThid)
163     DO j=1-3,sNy+3
164     Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
165     IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN
166     DO i=1-2,sNx+3
167     C THESE ARRAY BOUNDS INSURE THAT AFTER THIS STEP,
168     C VALUES WILL BE RELIABLE 2 GRID CELLS OUT IN THE
169     C X DIRECTION AND 3 CELLS OUT IN THE Y DIR
170     IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .or.
171     & ((STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) .and.
172     & (STREAMICE_hmask(i,j,bi,bj).ne.1.0))) THEN
173    
174     Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
175    
176     IF (STREAMICE_ufacemask(i,j,bi,bj).eq.4.0) THEN
177     hflux_x (i,j,bi,bj) = u_flux_bdry_SI (i,j,bi,bj)
178     ELSE
179    
180     uface = .5 *
181     & (U_streamice(i,j,bi,bj)+U_streamice(i,j+1,bi,bj))
182    
183    
184     IF (uface .gt. 0. _d 0) THEN
185     DO k=-1,1
186     stencil (k) = h(i+k-1,j,bi,bj)
187     ENDDO
188     IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .and.
189     & (STREAMICE_hmask(i-2,j,bi,bj).eq.1.0)) H0_valid=.true.
190    
191     IF ((Gi.eq.1).and.(STREAMICE_hmask(i-1,j,bi,bj).eq.3.0))
192     & THEN ! we are at western bdry and there is a thick. bdry cond
193     hflux_x (i,j,bi,bj) = h(i-1,j,bi,bj) * uface
194     ELSEIF (H0_valid) THEN
195     phi = SLOPE_LIMITER (
196     & stencil(0)-stencil(-1),
197     & stencil(1)-stencil(0))
198     hflux_x (i,j,bi,bj) = uface *
199     & (stencil(0) - phi * .5 * (stencil(0)-stencil(1)))
200     ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme
201     hflux_x (i,j,bi,bj) = uface * stencil(0)
202     ENDIF
203    
204     ELSE ! uface <= 0
205    
206     DO k=-1,1
207     stencil (k) = h(i-k,j,bi,bj)
208     ENDDO
209     IF ((STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) .and.
210     & (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0)) H0_valid=.true.
211    
212     IF ((Gi.eq.Nx).and.(STREAMICE_hmask(i+1,j,bi,bj).eq.3.0))
213     & THEN ! we are at western bdry and there is a thick. bdry cond
214     hflux_x (i,j,bi,bj) = h(i+1,j,bi,bj) * uface
215     ELSEIF (H0_valid) THEN
216     phi = SLOPE_LIMITER (
217     & stencil(0)-stencil(-1),
218     & stencil(1)-stencil(0))
219     hflux_x (i,j,bi,bj) = uface *
220     & (stencil(0) - phi * .5 * (stencil(0)-stencil(1)))
221     ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme
222     hflux_x (i,j,bi,bj) = uface * stencil(0)
223     ENDIF
224    
225     ENDIF
226    
227     ENDIF
228    
229     if (streamice_ufacemask(i,j,bi,bj).eq.2.0) THEN
230     total_vol_out = total_vol_out +
231     & hflux_x (i,j,bi,bj) * time_step
232     ENDIF
233    
234     ENDIF
235     ENDDO
236     ENDIF
237     ENDDO
238     ENDDO
239     ENDDO
240    
241     C X-FLUXES AT CELL BOUNDARIES CALCULATED; NOW TAKE FLUX DIVERGENCE TO INCREMENT THICKNESS
242    
243     DO bj=myByLo(myThid),myByHi(myThid)
244     DO bi=myBxLo(myThid),myBxHi(myThid)
245     DO j=1-3,sNy+3
246     Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
247     IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN
248     DO i=1-2,sNx+2
249     IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
250     h(i,j,bi,bj) = h(i,j,bi,bj) - time_step *
251     & (hflux_x(i+1,j,bi,bj)*dyG(i+1,j,bi,bj) -
252     & hflux_x(i,j,bi,bj)*dyG(i,j,bi,bj)) *
253     & recip_rA (i,j,bi,bj)
254     ENDIF
255     ENDDO
256     ENDIF
257     ENDDO
258     ENDDO
259     ENDDO
260    
261     ! PRINT *, "TOTAL VOLUME OUT: ", total_vol_out
262    
263     #endif
264     RETURN
265     END SUBROUTINE STREAMICE_ADVECT_THICKNESS_X
266    
267     SUBROUTINE STREAMICE_ADVECT_THICKNESS_Y ( myThid ,
268     O hflux_y ,
269     O h ,
270     I time_step )
271    
272     IMPLICIT NONE
273    
274     C O hflux_y ! flux per unit width across face
275     C O h
276     C I time_step
277    
278     C === Global variables ===
279     #include "SIZE.h"
280     #include "GRID.h"
281     #include "EEPARAMS.h"
282     #include "PARAMS.h"
283     #include "STREAMICE.h"
284    
285     INTEGER myThid
286     _RL hflux_y (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
287     _RL h (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
288     _RL time_step
289    
290     #ifdef ALLOW_STREAMICE
291    
292     C LOCAL VARIABLES
293    
294     INTEGER i, j, bi, bj, Gi, Gj, k
295     _RL vface, phi
296     _RL stencil (-1:1)
297     LOGICAL H0_valid ! there are valid cells to calculate a
298     ! slope-limited 2nd order flux
299     _RL SLOPE_LIMITER
300     external SLOPE_LIMITER
301    
302     DO bj=myByLo(myThid),myByHi(myThid)
303     DO bi=myBxLo(myThid),myBxHi(myThid)
304     DO j=1-1,sNy+2
305     Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
306     DO i=1-2,sNx+2
307     Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
308     IF ((Gi.ge.1) .and. (Gi.le.Nx)) THEN
309     C THESE ARRAY BOUNDS INSURE THAT AFTER THIS STEP,
310     C VALUES WILL BE RELIABLE 1 GRID CELLS OUT IN THE
311     C Y DIRECTION
312     IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .or.
313     & ((STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) .and.
314     & (STREAMICE_hmask(i,j,bi,bj).ne.1.0))) THEN
315    
316     IF (STREAMICE_vfacemask(i,j,bi,bj).eq.4.0) THEN
317     hflux_y (i,j,bi,bj) = v_flux_bdry_SI (i,j,bi,bj)
318     ELSE
319    
320     vface = .5 *
321     & (V_streamice(i,j,bi,bj)+V_streamice(i+1,j,bi,bj))
322    
323     IF (vface .gt. 0. _d 0) THEN
324     DO k=-1,1
325     stencil (k) = h(i,j+k-1,bi,bj)
326     ENDDO
327     IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .and.
328     & (STREAMICE_hmask(i,j-2,bi,bj).eq.1.0)) H0_valid=.true.
329    
330     IF ((Gj.eq.1).and.(STREAMICE_hmask(i,j-1,bi,bj).eq.3.0))
331     & THEN ! we are at western bdry and there is a thick. bdry cond
332     hflux_y (i,j,bi,bj) = h(i,j-1,bi,bj) * vface
333     ELSEIF (H0_valid) THEN
334     phi = SLOPE_LIMITER (
335     & stencil(0)-stencil(-1),
336     & stencil(1)-stencil(0))
337     hflux_y (i,j,bi,bj) = vface *
338     & (stencil(0) - phi * .5 * (stencil(0)-stencil(1)))
339     ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme
340     hflux_y (i,j,bi,bj) = vface * stencil(0)
341     ENDIF
342     ELSE ! uface <= 0
343     DO k=-1,1
344     stencil (k) = h(i,j-k,bi,bj)
345     ENDDO
346     IF ((STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) .and.
347     & (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0)) H0_valid=.true.
348    
349     IF ((Gj.eq.Ny).and.(STREAMICE_hmask(i,j+1,bi,bj).eq.3.0))
350     & THEN ! we are at western bdry and there is a thick. bdry cond
351     hflux_y (i,j,bi,bj) = h(i,j+1,bi,bj) * vface
352     ELSEIF (H0_valid) THEN
353     phi = SLOPE_LIMITER (
354     & stencil(0)-stencil(-1),
355     & stencil(1)-stencil(0))
356     hflux_y (i,j,bi,bj) = vface *
357     & (stencil(0) - phi * .5 * (stencil(0)-stencil(1)))
358     ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme
359     hflux_y (i,j,bi,bj) = vface * stencil(0)
360     ENDIF
361    
362     ENDIF ! uface 0
363    
364     ENDIF
365     ENDIF
366     ENDIF
367     ENDDO
368     ENDDO
369     ENDDO
370     ENDDO
371    
372     C X-FLUXES AT CELL BOUNDARIES CALCULATED; NOW TAKE FLUX DIVERGENCE TO INCREMENT THICKNESS
373    
374    
375    
376     DO bj=myByLo(myThid),myByHi(myThid)
377     DO bi=myBxLo(myThid),myBxHi(myThid)
378     DO j=1-1,sNy+1
379     DO i=1-2,sNx+2
380     Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
381     IF ((Gi .ge. 1) .and. (Gi .le. Nx)) THEN
382     IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
383     h(i,j,bi,bj) = h(i,j,bi,bj) - time_step *
384     & (hflux_y(i,j+1,bi,bj)*dxG(i,j+1,bi,bj) -
385     & hflux_y(i,j,bi,bj)*dxG(i,j,bi,bj)) *
386     & recip_rA (i,j,bi,bj)
387     ENDIF
388     ENDIF
389     ENDDO
390     ENDDO
391     ENDDO
392     ENDDO
393    
394     ! CALL WRITE_FLD_XY_RL ("h_after_yflux","",
395     ! & h, 0, myThid)
396    
397     #endif
398     RETURN
399     END SUBROUTINE STREAMICE_ADVECT_THICKNESS_Y
400    

  ViewVC Help
Powered by ViewVC 1.1.22