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

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

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

revision 1.6 by dgoldberg, Wed Jan 9 21:56:18 2013 UTC revision 1.11 by dgoldberg, Wed Aug 27 19:29:13 2014 UTC
# Line 5  C $Name$ Line 5  C $Name$
5    
6  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7    
   
8  CBOP  CBOP
9        SUBROUTINE STREAMICE_ADVECT_THICKNESS ( myThid, time_step )        SUBROUTINE STREAMICE_ADVECT_THICKNESS ( myThid,myIter,time_step )
10    
11  C     /============================================================\  C     *============================================================*
12  C     | SUBROUTINE                                                 |    C     | SUBROUTINE                                                 |
13  C     | o                                                          |  C     | o                                                          |
14  C     |============================================================|  C     *============================================================*
15  C     |                                                            |  C     |                                                            |
16  C     \============================================================/  C     *============================================================*
17        IMPLICIT NONE        IMPLICIT NONE
18    
19  C     === Global variables ===  C     === Global variables ===
# Line 27  C     === Global variables === Line 26  C     === Global variables ===
26  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
27  # include "tamc.h"  # include "tamc.h"
28  #endif  #endif
29    #ifdef ALLOW_SHELFICE
30    # include "SHELFICE.h"
31    #endif
32    
33        INTEGER myThid        INTEGER myThid, myIter
34        _RL time_step        _RL time_step
35    
36  #ifdef ALLOW_STREAMICE  #ifdef ALLOW_STREAMICE
37    
38        INTEGER i, j, bi, bj        INTEGER i, j, bi, bj, Gi, Gj
39        _RL thick_bd        _RL thick_bd, uflux, vflux, max_icfl, loc_icfl
40        _RL SLOPE_LIMITER        _RL time_step_full, time_step_rem
41        _RL sec_per_year, time_step_loc, MR, SMB, TMB        _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        CHARACTER*(MAX_LEN_MBUF) msgBuf
       external SLOPE_LIMITER  
54    
55        sec_per_year = 365.*86400.        sec_per_year = 365.*86400.
56    
57        time_step_loc = time_step / sec_per_year        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  #ifdef ALLOW_AUTODIFF_TAMC
63  CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics  CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics
# Line 50  CADJ STORE streamice_hmask  = comlev1, k Line 65  CADJ STORE streamice_hmask  = comlev1, k
65    
66        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
67         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
68          DO j=1-OLy,sNy+OLy          DO j=1,sNy+1
69           DO i=1-OLx,sNx+OLx           DO i=1,sNx+1
70            H_streamice_prev(i,j,bi,bj) =  
71              H_streamice_prev(i,j,bi,bj) =
72       &     H_streamice(i,j,bi,bj)       &     H_streamice(i,j,bi,bj)
73    
74            hflux_x_SI (i,j,bi,bj) = 0. _d 0            hflux_x_SI (i,j,bi,bj) = 0. _d 0
75            hflux_y_SI (i,j,bi,bj) = 0. _d 0            hflux_y_SI (i,j,bi,bj) = 0. _d 0
76            hflux_x_SI2 (i,j,bi,bj) = 0. _d 0            h_after_uflux_SI(i,j,bi,bj) = H_streamice(i,j,bi,bj)
77            hflux_y_SI2 (i,j,bi,bj) = 0. _d 0            h_after_vflux_SI(i,j,bi,bj) = H_streamice(i,j,bi,bj)
78            IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN  
79             h_after_uflux_SI (i,j,bi,bj) =  
80       &      H_streamice (i,j,bi,bj)            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            ENDIF
109    
110            thick_bd = h_bdry_values_SI (i,j,bi,bj)            IF (STREAMICE_vfacemask(i,j,bi,bj).eq.3.0) THEN
111            IF (thick_bd .ne. 0. _d 0) THEN             BCMASKy(i,j,bi,bj) = 3.0
112              h_after_uflux_SI (i,j,bi,bj) = thick_bd             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            ENDIF
143    
144    
145           ENDDO           ENDDO
146          ENDDO          ENDDO
147         ENDDO         ENDDO
148        ENDDO        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  #ifdef ALLOW_AUTODIFF_TAMC
 CADJ STORE h_after_uflux_si = comlev1, key=ikey_dynamics  
187  CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics  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  #endif
200    
201        CALL STREAMICE_ADVECT_THICKNESS_X ( myThid,        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,       O   hflux_x_SI,
      O   h_after_uflux_SI,  
207       I   time_step_loc )       I   time_step_loc )
208    
209    
210    
211        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
212         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
213          DO j=1-OLy,sNy+OLy          DO j=1-3,sNy+3
214           DO i=1-OLx,sNx+OLx           DO i=1,sNx
215            h_after_vflux_SI (i,j,bi,bj) =            Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
216       &     h_after_uflux_SI (i,j,bi,bj)            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           ENDDO
233          ENDDO          ENDDO
234         ENDDO         ENDDO
235        ENDDO        ENDDO
236    
237    
238  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 CADJ STORE h_after_vflux_si = comlev1, key=ikey_dynamics  
239  CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics  CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics
240  #endif  #endif
241    
242        CALL STREAMICE_ADVECT_THICKNESS_Y ( myThid,  !      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,       O   hflux_y_SI,
      O   h_after_vflux_SI,  
253       I   time_step_loc )       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)        DO bj=myByLo(myThid),myByHi(myThid)
284         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
285          DO j=1-OLy,sNy+OLy          DO j=1,sNy
286           DO i=1-OLx,sNx+OLx           DO i=1,sNx
287            IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN            IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
288             H_streamice (i,j,bi,bj) =             H_streamice (i,j,bi,bj) =
289       &      h_after_vflux_SI (i,j,bi,bj)       &      h_after_vflux_SI (i,j,bi,bj)
290            ENDIF            ENDIF
291           ENDDO           ENDDO
# Line 121  CADJ STORE streamice_hmask  = comlev1, k Line 293  CADJ STORE streamice_hmask  = comlev1, k
293         ENDDO         ENDDO
294        ENDDO        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    
       CALL STREAMICE_ADV_FRONT ( myThid, time_step_loc )  
321    
322    
323  ! NOW WE APPLY MELT RATES !!  ! NOW WE APPLY MELT RATES !!
324  ! THIS MAY BE MOVED TO A SEPARATE SUBROUTINE  ! THIS MAY BE MOVED TO A SEPARATE SUBROUTINE
325          irho = 1./streamice_density
326    
327        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
328         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
329          DO j=1-OLy,sNy+OLy          DO j=1,sNy
330           DO i=1-OLx,sNx+OLx           DO i=1,sNx
331             IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0 .or.            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       &       STREAMICE_hmask(i,j,bi,bj).eq.2.0) THEN
335               MR = (1.-float_frac_streamice(i,j,bi,bj)) *  
336       &             BDOT_STREAMICE(i,j,bi,bj)             IF (STREAMICE_allow_cpl) THEN
337               SMB = ADOT_STREAMICE(i,j,bi,bj)  #ifdef ALLOW_SHELFICE
338               TMB = SMB - MR              MR = -1. * (1.-float_frac_streamice(i,j,bi,bj)) *
339               IF ((TMB.lt.0.0) .and.       &        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.       &         (MR * time_step_loc .gt.
355       &          H_streamice (i,j,bi,bj))) THEN       &          H_streamice (i,j,bi,bj))) THEN
356                  H_streamice (i,j,bi,bj) = 0. _d 0                  H_streamice (i,j,bi,bj) = 0. _d 0
357                  STREAMICE_hmask(i,j,bi,bj) = 0.                  STREAMICE_hmask(i,j,bi,bj) = 0.0
358               ELSE                 PRINT *, "GOT HERE melted away! ", i,j
359                 H_streamice (i,j,bi,bj) =             ELSE
360       &          H_streamice (i,j,bi,bj) + TMB                 H_streamice (i,j,bi,bj) =
361               ENDIF       &          H_streamice (i,j,bi,bj) + TMB * time_step_loc
362             ENDIF             ENDIF
363    
364    
365              ENDIF
366           ENDDO           ENDDO
367          ENDDO          ENDDO
368         ENDDO         ENDDO
369        ENDDO            ENDDO
370    
371                _EXCH_XY_RL(H_streamice,myThid)
372    
373    
374        WRITE(msgBuf,'(A)') 'END STREAMICE_ADVECT_THICKNESS'        WRITE(msgBuf,'(A)') 'END STREAMICE_ADVECT_THICKNESS'
375         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
376       &                     SQUEEZE_RIGHT , 1)       &                     SQUEEZE_RIGHT , 1)
377          
378  #endif  #endif
379        END        END
   
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22