/[MITgcm]/MITgcm_contrib/ksnow/press_release/code_expt/streamice_advect_thickness.F
ViewVC logotype

Annotation of /MITgcm_contrib/ksnow/press_release/code_expt/streamice_advect_thickness.F

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


Revision 1.1 - (hide annotations) (download)
Fri Dec 16 15:25:29 2016 UTC (9 years, 8 months ago) by ksnow
Branch: MAIN
Adding shelfice_remeshing files for experiment

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

  ViewVC Help
Powered by ViewVC 1.1.22