/[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.11 - (hide 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
updating contrib streamice repo with latest files, and separated out convergence checks; and parameterised maximum iteration counts and interface w shelfice for coupling

1 dgoldberg 1.11 C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_advect_thickness.F,v 1.7 2014/04/24 12:01:50 dgoldberg Exp $
2 heimbach 1.1 C $Name: $
3    
4     #include "STREAMICE_OPTIONS.h"
5    
6     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7    
8     CBOP
9 dgoldberg 1.8 SUBROUTINE STREAMICE_ADVECT_THICKNESS ( myThid,myIter,time_step )
10 heimbach 1.1
11 dgoldberg 1.11 C *============================================================*
12     C | SUBROUTINE |
13 heimbach 1.1 C | o |
14 dgoldberg 1.11 C *============================================================*
15 heimbach 1.1 C | |
16 dgoldberg 1.11 C *============================================================*
17 heimbach 1.1 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 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
27     # include "tamc.h"
28     #endif
29 dgoldberg 1.11 #ifdef ALLOW_SHELFICE
30     # include "SHELFICE.h"
31     #endif
32 heimbach 1.1
33 dgoldberg 1.8 INTEGER myThid, myIter
34 heimbach 1.1 _RL time_step
35    
36     #ifdef ALLOW_STREAMICE
37    
38 dgoldberg 1.10 INTEGER i, j, bi, bj, Gi, Gj
39 dgoldberg 1.9 _RL thick_bd, uflux, vflux, max_icfl, loc_icfl
40     _RL time_step_full, time_step_rem
41 dgoldberg 1.11 _RL sec_per_year, time_step_loc, MR, SMB, TMB, irho
42 dgoldberg 1.8 _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 dgoldberg 1.11
53 dgoldberg 1.3 CHARACTER*(MAX_LEN_MBUF) msgBuf
54 heimbach 1.1
55     sec_per_year = 365.*86400.
56    
57     time_step_loc = time_step / sec_per_year
58 dgoldberg 1.9 time_step_full = time_step_loc
59     time_step_rem = time_step_loc
60 dgoldberg 1.7 PRINT *, "time_step_loc ", time_step_loc
61    
62 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
63     CADJ STORE streamice_hmask = comlev1, key=ikey_dynamics
64     #endif
65    
66 heimbach 1.1 DO bj=myByLo(myThid),myByHi(myThid)
67     DO bi=myBxLo(myThid),myBxHi(myThid)
68 dgoldberg 1.10 DO j=1,sNy+1
69     DO i=1,sNx+1
70 dgoldberg 1.8
71 dgoldberg 1.11 H_streamice_prev(i,j,bi,bj) =
72 dgoldberg 1.5 & H_streamice(i,j,bi,bj)
73 dgoldberg 1.8
74 heimbach 1.1 hflux_x_SI (i,j,bi,bj) = 0. _d 0
75     hflux_y_SI (i,j,bi,bj) = 0. _d 0
76 dgoldberg 1.11 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 dgoldberg 1.8
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 dgoldberg 1.11 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 dgoldberg 1.8 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 heimbach 1.1 ENDIF
109    
110 dgoldberg 1.8 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 dgoldberg 1.11 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 dgoldberg 1.8 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 heimbach 1.1 ENDIF
143 dgoldberg 1.8
144    
145 heimbach 1.1 ENDDO
146     ENDDO
147     ENDDO
148     ENDDO
149    
150 dgoldberg 1.8 _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 dgoldberg 1.11 _EXCH_XY_RL(h_after_uflux_SI,myThid)
157     _EXCH_XY_RL(h_after_vflux_SI,myThid)
158 dgoldberg 1.3
159 dgoldberg 1.9 #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 dgoldberg 1.10
185 dgoldberg 1.9
186 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
187     CADJ STORE streamice_hmask = comlev1, key=ikey_dynamics
188 dgoldberg 1.8 CADJ STORE H_streamice = comlev1, key=ikey_dynamics
189 heimbach 1.2 #endif
190    
191 dgoldberg 1.9 #ifndef ALLOW_AUTODIFF_TAMC
192 dgoldberg 1.11 do while (time_step_rem .gt. 1.e-15)
193 dgoldberg 1.9 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 dgoldberg 1.8
201     CALL STREAMICE_ADV_FLUX_FL_X ( myThid ,
202     I utrans ,
203     I H_streamice ,
204     I BCMASKX,
205     I BCVALX,
206 heimbach 1.1 O hflux_x_SI,
207     I time_step_loc )
208    
209 dgoldberg 1.3
210 heimbach 1.1
211     DO bj=myByLo(myThid),myByHi(myThid)
212     DO bi=myBxLo(myThid),myBxHi(myThid)
213 dgoldberg 1.8 DO j=1-3,sNy+3
214 dgoldberg 1.10 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 dgoldberg 1.8 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 dgoldberg 1.11 & hflux_x_SI(i,j,bi,bj)*dyG(i,j,bi,bj))
224 dgoldberg 1.8 & * 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 dgoldberg 1.10
231     ENDIF
232 heimbach 1.1 ENDDO
233     ENDDO
234     ENDDO
235     ENDDO
236    
237 dgoldberg 1.8
238 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
239     CADJ STORE streamice_hmask = comlev1, key=ikey_dynamics
240     #endif
241    
242 dgoldberg 1.8 ! 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 heimbach 1.1 O hflux_y_SI,
253     I time_step_loc )
254    
255 dgoldberg 1.3
256 dgoldberg 1.8 DO bj=myByLo(myThid),myByHi(myThid)
257     DO bi=myBxLo(myThid),myBxHi(myThid)
258 dgoldberg 1.10 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 dgoldberg 1.8 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 dgoldberg 1.10 ENDIF
277    
278 dgoldberg 1.8 ENDDO
279     ENDDO
280     ENDDO
281     ENDDO
282 heimbach 1.1
283     DO bj=myByLo(myThid),myByHi(myThid)
284     DO bi=myBxLo(myThid),myBxHi(myThid)
285 dgoldberg 1.8 DO j=1,sNy
286     DO i=1,sNx
287 heimbach 1.1 IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
288 dgoldberg 1.11 H_streamice (i,j,bi,bj) =
289 heimbach 1.1 & h_after_vflux_SI (i,j,bi,bj)
290     ENDIF
291     ENDDO
292     ENDDO
293     ENDDO
294     ENDDO
295    
296 dgoldberg 1.10 ! NOTE: AT THIS POINT H IS NOT VALID ON OVERLAP!!!
297    
298 dgoldberg 1.11 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 dgoldberg 1.3
305 dgoldberg 1.8 #ifdef ALLOW_STREAMICE_2DTRACER
306     CALL STREAMICE_ADVECT_2DTRACER(
307     & myThid,
308     & myIter,
309     & time_step,
310     & uTrans,
311     & vTrans,
312 dgoldberg 1.11 & BCMASKx,
313     & BCMASKy )
314 dgoldberg 1.8 #endif
315 heimbach 1.1
316 dgoldberg 1.9 #ifndef ALLOW_AUTODIFF_TAMC
317 dgoldberg 1.11 time_step_rem = time_step_rem - time_step_loc
318 dgoldberg 1.9 enddo
319     #endif
320    
321    
322 dgoldberg 1.3
323 dgoldberg 1.4 ! NOW WE APPLY MELT RATES !!
324     ! THIS MAY BE MOVED TO A SEPARATE SUBROUTINE
325 dgoldberg 1.11 irho = 1./streamice_density
326 heimbach 1.1
327 dgoldberg 1.4 DO bj=myByLo(myThid),myByHi(myThid)
328     DO bi=myBxLo(myThid),myBxHi(myThid)
329 dgoldberg 1.8 DO j=1,sNy
330     DO i=1,sNx
331 dgoldberg 1.11 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 dgoldberg 1.4 & STREAMICE_hmask(i,j,bi,bj).eq.2.0) THEN
335 dgoldberg 1.11
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 dgoldberg 1.6 & (MR * time_step_loc .gt.
355     & H_streamice (i,j,bi,bj))) THEN
356 dgoldberg 1.4 H_streamice (i,j,bi,bj) = 0. _d 0
357 dgoldberg 1.11 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 dgoldberg 1.7 & H_streamice (i,j,bi,bj) + TMB * time_step_loc
362 dgoldberg 1.4 ENDIF
363 dgoldberg 1.11
364    
365     ENDIF
366 dgoldberg 1.4 ENDDO
367     ENDDO
368     ENDDO
369 dgoldberg 1.11 ENDDO
370    
371     _EXCH_XY_RL(H_streamice,myThid)
372 dgoldberg 1.8
373 dgoldberg 1.10
374 dgoldberg 1.3 WRITE(msgBuf,'(A)') 'END STREAMICE_ADVECT_THICKNESS'
375     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
376     & SQUEEZE_RIGHT , 1)
377 dgoldberg 1.11
378 heimbach 1.1 #endif
379 heimbach 1.2 END

  ViewVC Help
Powered by ViewVC 1.1.22