/[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.10 - (show annotations) (download)
Wed Jun 12 20:48:08 2013 UTC (12 years, 1 month ago) by dgoldberg
Branch: MAIN
Changes since 1.9: +25 -9 lines
take care of errors found by debug

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

  ViewVC Help
Powered by ViewVC 1.1.22