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

Contents of /MITgcm_contrib/dgoldberg/streamice/streamice_advect_thickness.F

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


Revision 1.11 - (show annotations) (download)
Wed Aug 27 19:29:13 2014 UTC (10 years, 10 months ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
Changes since 1.10: +86 -41 lines
Error occurred while calculating annotation data.
updating contrib streamice repo with latest files, and separated out convergence checks; and parameterised maximum iteration counts and interface w shelfice for coupling

1 C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_advect_thickness.F,v 1.7 2014/04/24 12:01:50 dgoldberg 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_ADVECT_THICKNESS ( myThid,myIter,time_step )
10
11 C *============================================================*
12 C | SUBROUTINE |
13 C | o |
14 C *============================================================*
15 C | |
16 C *============================================================*
17 IMPLICIT NONE
18
19 C === Global variables ===
20 #include "SIZE.h"
21 #include "GRID.h"
22 #include "EEPARAMS.h"
23 #include "PARAMS.h"
24 #include "STREAMICE.h"
25 #include "STREAMICE_ADV.h"
26 #ifdef ALLOW_AUTODIFF_TAMC
27 # include "tamc.h"
28 #endif
29 #ifdef ALLOW_SHELFICE
30 # include "SHELFICE.h"
31 #endif
32
33 INTEGER myThid, myIter
34 _RL time_step
35
36 #ifdef ALLOW_STREAMICE
37
38 INTEGER i, j, bi, bj, Gi, Gj
39 _RL thick_bd, uflux, vflux, max_icfl, loc_icfl
40 _RL time_step_full, time_step_rem
41 _RL sec_per_year, time_step_loc, MR, SMB, TMB, irho
42 _RL BCVALX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
43 _RL BCVALY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
44 _RS BCMASKX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
45 _RS BCMASKY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
46 _RL utrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
47 _RL vtrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
48 _RL h_after_uflux_SI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
49 _RL h_after_vflux_SI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
50 _RL hflux_x_SI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
51 _RL hflux_y_SI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
52
53 CHARACTER*(MAX_LEN_MBUF) msgBuf
54
55 sec_per_year = 365.*86400.
56
57 time_step_loc = time_step / sec_per_year
58 time_step_full = time_step_loc
59 time_step_rem = time_step_loc
60 PRINT *, "time_step_loc ", time_step_loc
61
62 #ifdef ALLOW_AUTODIFF_TAMC
63 CADJ STORE streamice_hmask = comlev1, key=ikey_dynamics
64 #endif
65
66 DO bj=myByLo(myThid),myByHi(myThid)
67 DO bi=myBxLo(myThid),myBxHi(myThid)
68 DO j=1,sNy+1
69 DO i=1,sNx+1
70
71 H_streamice_prev(i,j,bi,bj) =
72 & H_streamice(i,j,bi,bj)
73
74 hflux_x_SI (i,j,bi,bj) = 0. _d 0
75 hflux_y_SI (i,j,bi,bj) = 0. _d 0
76 h_after_uflux_SI(i,j,bi,bj) = H_streamice(i,j,bi,bj)
77 h_after_vflux_SI(i,j,bi,bj) = H_streamice(i,j,bi,bj)
78
79
80 IF (STREAMICE_ufacemask(i,j,bi,bj).eq.3.0) THEN
81 BCMASKX(i,j,bi,bj) = 3.0
82 BCVALX(i,j,bi,bj) = h_ubdry_values_SI(i,j,bi,bj)
83 utrans(i,j,bi,bj) = .5 * (
84 & u_streamice(i,j,bi,bj)+u_streamice(i,j+1,bi,bj))
85 ELSEIF (STREAMICE_ufacemask(i,j,bi,bj).eq.4.0) THEN
86 IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
87 uflux = u_flux_bdry_SI(i,j,bi,bj)
88 BCMASKX(i,j,bi,bj) = 3.0
89 BCVALX(i,j,bi,bj) = uflux
90 utrans(i,j,bi,bj) = 1.0
91 ELSEIF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0) THEN
92 uflux = u_flux_bdry_SI(i,j,bi,bj)
93 BCMASKX(i,j,bi,bj) = 3.0
94 BCVALX(i,j,bi,bj) = uflux
95 utrans(i,j,bi,bj) = -1.0
96 ENDIF
97 ELSEIF (.not.(
98 & STREAMICE_hmask(i,j,bi,bj).eq.1.0.OR.
99 & STREAMICE_hmask(i-1,j,bi,bj).eq.1.0)) THEN
100 BCMASKX(i,j,bi,bj) = 0.0
101 BCVALX(i,j,bi,bj) = 0. _d 0
102 utrans(i,j,bi,bj) = 0. _d 0
103 ELSE
104 BCMASKX(i,j,bi,bj) = 0.0
105 BCVALX(i,j,bi,bj) = 0. _d 0
106 utrans(i,j,bi,bj) = .5 * (
107 & u_streamice(i,j,bi,bj)+u_streamice(i,j+1,bi,bj))
108 ENDIF
109
110 IF (STREAMICE_vfacemask(i,j,bi,bj).eq.3.0) THEN
111 BCMASKy(i,j,bi,bj) = 3.0
112 BCVALy(i,j,bi,bj) = h_vbdry_values_SI(i,j,bi,bj)
113 vtrans(i,j,bi,bj) = .5 * (
114 & v_streamice(i,j,bi,bj)+v_streamice(i+1,j,bi,bj))
115 ELSEIF (STREAMICE_vfacemask(i,j,bi,bj).eq.4.0) THEN
116 IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
117 vflux = v_flux_bdry_SI(i,j,bi,bj)
118 BCMASKY(i,j,bi,bj) = 3.0
119 BCVALY(i,j,bi,bj) = vflux
120 vtrans(i,j,bi,bj) = 1.0
121 ELSEIF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0) THEN
122 vflux = v_flux_bdry_SI(i,j,bi,bj)
123 BCMASKY(i,j,bi,bj) = 3.0
124 BCVALY(i,j,bi,bj) = vflux
125 vtrans(i,j,bi,bj) = -1.0
126 ENDIF
127
128
129
130 vtrans(i,j,bi,bj) = 1.0
131 ELSEIF (.not.(
132 & STREAMICE_hmask(i,j,bi,bj).eq.1.0.OR.
133 & STREAMICE_hmask(i,j-1,bi,bj).eq.1.0)) THEN
134 BCMASKY(i,j,bi,bj) = 0.0
135 BCVALY(i,j,bi,bj) = 0. _d 0
136 vtrans(i,j,bi,bj) = 0. _d 0
137 ELSE
138 BCMASKy(i,j,bi,bj) = 0.0
139 BCVALy(i,j,bi,bj) = 0. _d 0
140 vtrans(i,j,bi,bj) = .5 * (
141 & v_streamice(i,j,bi,bj)+v_streamice(i+1,j,bi,bj))
142 ENDIF
143
144
145 ENDDO
146 ENDDO
147 ENDDO
148 ENDDO
149
150 _EXCH_XY_RL(utrans,myThid)
151 _EXCH_XY_RL(vtrans,myThid)
152 _EXCH_XY_RS(BCMASKx,myThid)
153 _EXCH_XY_RS(BCMASKy,myThid)
154 _EXCH_XY_RL(BCVALX,myThid)
155 _EXCH_XY_RL(BCVALY,myThid)
156 _EXCH_XY_RL(h_after_uflux_SI,myThid)
157 _EXCH_XY_RL(h_after_vflux_SI,myThid)
158
159 #ifndef ALLOW_AUTODIFF_TAMC
160
161 max_icfl = 1.e-20
162
163 DO bj=myByLo(myThid),myByHi(myThid)
164 DO bi=myBxLo(myThid),myBxHi(myThid)
165 DO j=1,sNy
166 DO i=1,sNx
167 IF (streamice_hmask(i,j,bi,bj).eq.1.0) THEN
168 loc_icfl=max(abs(utrans(i,j,bi,bj)),
169 & abs(utrans(i+1,j,bi,bj))) / dxF(i,j,bi,bj)
170 loc_icfl=max(loc_icfl,max(abs(vtrans(i,j,bi,bj)),
171 & abs(vtrans(i,j+1,bi,bj))) / dyF(i,j,bi,bj))
172 if (loc_icfl.gt.max_icfl) then
173 max_icfl = loc_icfl
174 ENDIF
175 ENDIF
176 ENDDO
177 ENDDO
178 ENDDO
179 ENDDO
180
181 CALL GLOBAL_MAX_R8 (max_icfl, myThid)
182
183 #endif
184
185
186 #ifdef ALLOW_AUTODIFF_TAMC
187 CADJ STORE streamice_hmask = comlev1, key=ikey_dynamics
188 CADJ STORE H_streamice = comlev1, key=ikey_dynamics
189 #endif
190
191 #ifndef ALLOW_AUTODIFF_TAMC
192 do while (time_step_rem .gt. 1.e-15)
193 time_step_loc = min (
194 & streamice_cfl_factor / max_icfl,
195 & time_step_rem )
196 if (time_step_loc .lt. time_step_full) then
197 PRINT *, "TAKING PARTIAL TIME STEP", time_step_loc
198 endif
199 #endif
200
201 CALL STREAMICE_ADV_FLUX_FL_X ( myThid ,
202 I utrans ,
203 I H_streamice ,
204 I BCMASKX,
205 I BCVALX,
206 O hflux_x_SI,
207 I time_step_loc )
208
209
210
211 DO bj=myByLo(myThid),myByHi(myThid)
212 DO bi=myBxLo(myThid),myBxHi(myThid)
213 DO j=1-3,sNy+3
214 DO i=1,sNx
215 Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
216 Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
217 IF (((Gj .ge. 1) .and. (Gj .le. Ny))
218 & .or.STREAMICE_NS_PERIODIC) THEN
219
220 IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
221 h_after_uflux_SI (i,j,bi,bj) = H_streamice(i,j,bi,bj) -
222 & (hflux_x_SI(i+1,j,bi,bj)*dyG(i+1,j,bi,bj) -
223 & hflux_x_SI(i,j,bi,bj)*dyG(i,j,bi,bj))
224 & * recip_rA (i,j,bi,bj) * time_step_loc
225 IF ( h_after_uflux_SI (i,j,bi,bj).le.0.0) THEN
226 PRINT *, "h neg after x", i,j,hflux_x_SI(i+1,j,bi,bj),
227 & hflux_x_SI(i,j,bi,bj)
228 ENDIF
229 ENDIF
230
231 ENDIF
232 ENDDO
233 ENDDO
234 ENDDO
235 ENDDO
236
237
238 #ifdef ALLOW_AUTODIFF_TAMC
239 CADJ STORE streamice_hmask = comlev1, key=ikey_dynamics
240 #endif
241
242 ! CALL STREAMICE_ADVECT_THICKNESS_Y ( myThid,
243 ! O hflux_y_SI,
244 ! O h_after_vflux_SI,
245 ! I time_step_loc )
246
247 CALL STREAMICE_ADV_FLUX_FL_Y ( myThid ,
248 I vtrans ,
249 I h_after_uflux_si ,
250 I BCMASKY,
251 I BCVALY,
252 O hflux_y_SI,
253 I time_step_loc )
254
255
256 DO bj=myByLo(myThid),myByHi(myThid)
257 DO bi=myBxLo(myThid),myBxHi(myThid)
258 DO j=1,sNy
259 DO i=1,sNx
260 Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
261 Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
262 IF (((Gj .ge. 1) .and. (Gj .le. Ny))
263 & .or.STREAMICE_EW_PERIODIC) THEN
264
265 IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
266 h_after_vflux_SI (i,j,bi,bj) = h_after_uflux_SI(i,j,bi,bj) -
267 & (hflux_y_SI(i,j+1,bi,bj)*dxG(i,j+1,bi,bj) -
268 & hflux_y_SI(i,j,bi,bj)*dxG(i,j,bi,bj)) *
269 & recip_rA (i,j,bi,bj) * time_step_loc
270 IF ( h_after_vflux_SI (i,j,bi,bj).le.0.0) THEN
271 PRINT *, "h neg after y", i,j,hflux_y_SI(i,j+1,bi,bj),
272 & hflux_y_SI(i,j,bi,bj)
273 ENDIF
274
275 ENDIF
276 ENDIF
277
278 ENDDO
279 ENDDO
280 ENDDO
281 ENDDO
282
283 DO bj=myByLo(myThid),myByHi(myThid)
284 DO bi=myBxLo(myThid),myBxHi(myThid)
285 DO j=1,sNy
286 DO i=1,sNx
287 IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
288 H_streamice (i,j,bi,bj) =
289 & h_after_vflux_SI (i,j,bi,bj)
290 ENDIF
291 ENDDO
292 ENDDO
293 ENDDO
294 ENDDO
295
296 ! NOTE: AT THIS POINT H IS NOT VALID ON OVERLAP!!!
297
298 if (streamice_move_front) then
299 CALL STREAMICE_ADV_FRONT (
300 & myThid, time_step_loc,
301 & hflux_x_SI, hflux_y_SI )
302 endif
303
304
305 #ifdef ALLOW_STREAMICE_2DTRACER
306 CALL STREAMICE_ADVECT_2DTRACER(
307 & myThid,
308 & myIter,
309 & time_step,
310 & uTrans,
311 & vTrans,
312 & BCMASKx,
313 & BCMASKy )
314 #endif
315
316 #ifndef ALLOW_AUTODIFF_TAMC
317 time_step_rem = time_step_rem - time_step_loc
318 enddo
319 #endif
320
321
322
323 ! NOW WE APPLY MELT RATES !!
324 ! THIS MAY BE MOVED TO A SEPARATE SUBROUTINE
325 irho = 1./streamice_density
326
327 DO bj=myByLo(myThid),myByHi(myThid)
328 DO bi=myBxLo(myThid),myBxHi(myThid)
329 DO j=1,sNy
330 DO i=1,sNx
331 Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
332 Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
333 IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0 .or.
334 & STREAMICE_hmask(i,j,bi,bj).eq.2.0) THEN
335
336 IF (STREAMICE_allow_cpl) THEN
337 #ifdef ALLOW_SHELFICE
338 MR = -1. * (1.-float_frac_streamice(i,j,bi,bj)) *
339 & shelfIceFreshWaterFlux(I,J,bi,bj) * irho *
340 & sec_per_year
341
342 #else
343 STOP 'SHELFICE IS NOT ENABLED'
344 #endif
345 ELSE
346 MR = (1.-float_frac_streamice(i,j,bi,bj)) *
347 & (BDOT_STREAMICE(i,j,bi,bj) +
348 & BDOT_pert(i,j,bi,bj))
349 ENDIF
350
351 SMB = ADOT_STREAMICE(i,j,bi,bj)
352 TMB = SMB - MR
353 IF ((TMB.lt.0.0) .and.
354 & (MR * time_step_loc .gt.
355 & H_streamice (i,j,bi,bj))) THEN
356 H_streamice (i,j,bi,bj) = 0. _d 0
357 STREAMICE_hmask(i,j,bi,bj) = 0.0
358 PRINT *, "GOT HERE melted away! ", i,j
359 ELSE
360 H_streamice (i,j,bi,bj) =
361 & H_streamice (i,j,bi,bj) + TMB * time_step_loc
362 ENDIF
363
364
365 ENDIF
366 ENDDO
367 ENDDO
368 ENDDO
369 ENDDO
370
371 _EXCH_XY_RL(H_streamice,myThid)
372
373
374 WRITE(msgBuf,'(A)') 'END STREAMICE_ADVECT_THICKNESS'
375 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
376 & SQUEEZE_RIGHT , 1)
377
378 #endif
379 END

  ViewVC Help
Powered by ViewVC 1.1.22