/[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.9 - (show annotations) (download)
Tue Jun 11 17:42:17 2013 UTC (12 years, 1 month ago) by dgoldberg
Branch: MAIN
Changes since 1.8: +47 -3 lines
allow subcycling for thickness advec

1 C $Header: /u/gcmpack/MITgcm_contrib/dgoldberg/streamice/streamice_advect_thickness.F,v 1.8 2013/06/08 22:15:33 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
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-3,sNy+3
67 DO i=1-3,sNx+3
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-1,sNx+1
191 IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
192 h_after_uflux_SI (i,j,bi,bj) = H_streamice(i,j,bi,bj) -
193 & (hflux_x_SI(i+1,j,bi,bj)*dyG(i+1,j,bi,bj) -
194 & hflux_x_SI(i,j,bi,bj)*dyG(i,j,bi,bj))
195 & * recip_rA (i,j,bi,bj) * time_step_loc
196 IF ( h_after_uflux_SI (i,j,bi,bj).le.0.0) THEN
197 PRINT *, "h neg after x", i,j,hflux_x_SI(i+1,j,bi,bj),
198 & hflux_x_SI(i,j,bi,bj)
199 ENDIF
200 ENDIF
201 ENDDO
202 ENDDO
203 ENDDO
204 ENDDO
205
206
207
208 #ifdef ALLOW_AUTODIFF_TAMC
209 CADJ STORE streamice_hmask = comlev1, key=ikey_dynamics
210 #endif
211
212 ! CALL STREAMICE_ADVECT_THICKNESS_Y ( myThid,
213 ! O hflux_y_SI,
214 ! O h_after_vflux_SI,
215 ! I time_step_loc )
216
217 CALL STREAMICE_ADV_FLUX_FL_Y ( myThid ,
218 I vtrans ,
219 I h_after_uflux_si ,
220 I BCMASKY,
221 I BCVALY,
222 O hflux_y_SI,
223 I time_step_loc )
224
225
226 DO bj=myByLo(myThid),myByHi(myThid)
227 DO bi=myBxLo(myThid),myBxHi(myThid)
228 DO j=1-1,sNy+1
229 DO i=1-1,sNx+1
230 IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
231 h_after_vflux_SI (i,j,bi,bj) = h_after_uflux_SI(i,j,bi,bj) -
232 & (hflux_y_SI(i,j+1,bi,bj)*dxG(i,j+1,bi,bj) -
233 & hflux_y_SI(i,j,bi,bj)*dxG(i,j,bi,bj)) *
234 & recip_rA (i,j,bi,bj) * time_step_loc
235 IF ( h_after_vflux_SI (i,j,bi,bj).le.0.0) THEN
236 PRINT *, "h neg after y", i,j,hflux_y_SI(i,j+1,bi,bj),
237 & hflux_y_SI(i,j,bi,bj)
238 ENDIF
239
240 ENDIF
241 ENDDO
242 ENDDO
243 ENDDO
244 ENDDO
245
246 DO bj=myByLo(myThid),myByHi(myThid)
247 DO bi=myBxLo(myThid),myBxHi(myThid)
248 DO j=1,sNy
249 DO i=1,sNx
250 IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
251 H_streamice (i,j,bi,bj) =
252 & h_after_vflux_SI (i,j,bi,bj)
253 ENDIF
254 ENDDO
255 ENDDO
256 ENDDO
257 ENDDO
258
259 ! CALL STREAMICE_ADV_FRONT (
260 ! & myThid, time_step_loc,
261 ! & hflux_x_SI, hflux_y_SI )
262
263
264 #ifdef ALLOW_STREAMICE_2DTRACER
265 CALL STREAMICE_ADVECT_2DTRACER(
266 & myThid,
267 & myIter,
268 & time_step,
269 & uTrans,
270 & vTrans,
271 & bcMaskx,
272 & bcMasky )
273 #endif
274
275 #ifndef ALLOW_AUTODIFF_TAMC
276 time_step_rem = time_step_rem - time_step_loc
277 enddo
278 #endif
279
280
281
282 ! NOW WE APPLY MELT RATES !!
283 ! THIS MAY BE MOVED TO A SEPARATE SUBROUTINE
284
285 DO bj=myByLo(myThid),myByHi(myThid)
286 DO bi=myBxLo(myThid),myBxHi(myThid)
287 DO j=1,sNy
288 DO i=1,sNx
289 IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0 .or.
290 & STREAMICE_hmask(i,j,bi,bj).eq.2.0) THEN
291 MR = (1.-float_frac_streamice(i,j,bi,bj)) *
292 & BDOT_STREAMICE(i,j,bi,bj)
293 SMB = ADOT_STREAMICE(i,j,bi,bj)
294 TMB = SMB - MR
295 IF ((TMB.lt.0.0) .and.
296 & (MR * time_step_loc .gt.
297 & H_streamice (i,j,bi,bj))) THEN
298 H_streamice (i,j,bi,bj) = 0. _d 0
299 STREAMICE_hmask(i,j,bi,bj) = 0.
300 ELSE
301 H_streamice (i,j,bi,bj) =
302 & H_streamice (i,j,bi,bj) + TMB * time_step_loc
303 ENDIF
304 ENDIF
305 ENDDO
306 ENDDO
307 ENDDO
308 ENDDO
309
310 _EXCH_XY_RL (H_streamice,myThid)
311
312 WRITE(msgBuf,'(A)') 'END STREAMICE_ADVECT_THICKNESS'
313 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
314 & SQUEEZE_RIGHT , 1)
315
316 #endif
317 END
318

  ViewVC Help
Powered by ViewVC 1.1.22