/[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.3 by dgoldberg, Thu Jul 26 16:13:18 2012 UTC revision 1.10 by dgoldberg, Wed Jun 12 20:48:08 2013 UTC
# Line 7  C---+----1----+----2----+----3----+----4 Line 7  C---+----1----+----2----+----3----+----4
7    
8    
9  CBOP  CBOP
10        SUBROUTINE STREAMICE_ADVECT_THICKNESS ( myThid, time_step )        SUBROUTINE STREAMICE_ADVECT_THICKNESS ( myThid,myIter,time_step )
11    
12  C     /============================================================\  C     /============================================================\
13  C     | SUBROUTINE                                                 |    C     | SUBROUTINE                                                 |  
# Line 28  C     === Global variables === Line 28  C     === Global variables ===
28  # include "tamc.h"  # include "tamc.h"
29  #endif  #endif
30    
31        INTEGER myThid        INTEGER myThid, myIter
32        _RL time_step        _RL time_step
33    
34  #ifdef ALLOW_STREAMICE  #ifdef ALLOW_STREAMICE
35    
36        INTEGER i, j, bi, bj        INTEGER i, j, bi, bj, Gi, Gj
37        _RL thick_bd        _RL thick_bd, uflux, vflux, max_icfl, loc_icfl
38        _RL SLOPE_LIMITER        _RL time_step_full, time_step_rem
39        _RL sec_per_year, time_step_loc        _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        CHARACTER*(MAX_LEN_MBUF) msgBuf
       external SLOPE_LIMITER  
52    
53        sec_per_year = 365.*86400.        sec_per_year = 365.*86400.
54    
55        time_step_loc = time_step / sec_per_year        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  #ifdef ALLOW_AUTODIFF_TAMC
61  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 63  CADJ STORE streamice_hmask  = comlev1, k
63    
64        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
65         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
66          DO j=1-OLy,sNy+OLy          DO j=1,sNy+1
67           DO i=1-OLx,sNx+OLx           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            hflux_x_SI (i,j,bi,bj) = 0. _d 0
73            hflux_y_SI (i,j,bi,bj) = 0. _d 0            hflux_y_SI (i,j,bi,bj) = 0. _d 0
74            hflux_x_SI2 (i,j,bi,bj) = 0. _d 0  
75            hflux_y_SI2 (i,j,bi,bj) = 0. _d 0  
76            IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN            IF (STREAMICE_ufacemask(i,j,bi,bj).eq.3.0) THEN
77             h_after_uflux_SI (i,j,bi,bj) =             BCMASKX(i,j,bi,bj) = 3.0
78       &      H_streamice (i,j,bi,bj)             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            ENDIF
98    
99            thick_bd = h_bdry_values_SI (i,j,bi,bj)            IF (STREAMICE_vfacemask(i,j,bi,bj).eq.3.0) THEN
100            IF (thick_bd .ne. 0. _d 0) THEN             BCMASKy(i,j,bi,bj) = 3.0
101              h_after_uflux_SI (i,j,bi,bj) = thick_bd             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            ENDIF
121    
122    
123           ENDDO           ENDDO
124          ENDDO          ENDDO
125         ENDDO         ENDDO
126        ENDDO        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  #ifdef ALLOW_AUTODIFF_TAMC
 CADJ STORE h_after_uflux_si = comlev1, key=ikey_dynamics  
163  CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics  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  #endif
176    
177        CALL STREAMICE_ADVECT_THICKNESS_X ( myThid,        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,       O   hflux_x_SI,
      O   h_after_uflux_SI,  
183       I   time_step_loc )       I   time_step_loc )
184    
185    
186    
187        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
188         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
189          DO j=1-OLy,sNy+OLy          DO j=1-3,sNy+3
190           DO i=1-OLx,sNx+OLx           DO i=1,sNx
191            h_after_vflux_SI (i,j,bi,bj) =            Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
192       &     h_after_uflux_SI (i,j,bi,bj)            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           ENDDO
209          ENDDO          ENDDO
210         ENDDO         ENDDO
211        ENDDO        ENDDO
212    
213    
214  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 CADJ STORE h_after_vflux_si = comlev1, key=ikey_dynamics  
215  CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics  CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics
216  #endif  #endif
217    
218        CALL STREAMICE_ADVECT_THICKNESS_Y ( myThid,  !      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,       O   hflux_y_SI,
      O   h_after_vflux_SI,  
229       I   time_step_loc )       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)        DO bj=myByLo(myThid),myByHi(myThid)
260         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
261          DO j=1-OLy,sNy+OLy          DO j=1,sNy
262           DO i=1-OLx,sNx+OLx           DO i=1,sNx
263            IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN            IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
264             H_streamice (i,j,bi,bj) =             H_streamice (i,j,bi,bj) =
265       &      h_after_vflux_SI (i,j,bi,bj)       &      h_after_vflux_SI (i,j,bi,bj)
# Line 120  CADJ STORE streamice_hmask  = comlev1, k Line 269  CADJ STORE streamice_hmask  = comlev1, k
269         ENDDO         ENDDO
270        ENDDO        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        CALL STREAMICE_ADV_FRONT ( myThid, time_step_loc )  ! 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    
       _EXCH_XY_RL( H_streamice, myThid )  
       _EXCH_XY_RL( area_shelf_streamice, myThid )  
       _EXCH_XY_RL( STREAMICE_hmask, myThid )  
327                
328        WRITE(msgBuf,'(A)') 'END STREAMICE_ADVECT_THICKNESS'        WRITE(msgBuf,'(A)') 'END STREAMICE_ADVECT_THICKNESS'
329         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22