/[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.1 by heimbach, Thu Mar 29 15:59:21 2012 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 24  C     === Global variables === Line 23  C     === Global variables ===
23  #include "PARAMS.h"  #include "PARAMS.h"
24  #include "STREAMICE.h"  #include "STREAMICE.h"
25  #include "STREAMICE_ADV.h"  #include "STREAMICE_ADV.h"
26    #ifdef ALLOW_AUTODIFF_TAMC
27    # include "tamc.h"
28    #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        _RL sec_per_year, time_step_loc, MR, SMB, TMB, irho
42        external SLOPE_LIMITER        _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
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
63    CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics
64    #endif
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    
71              H_streamice_prev(i,j,bi,bj) =
72         &     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  
            h_after_uflux_SI (i,j,bi,bj) =  
      &      H_streamice (i,j,bi,bj)  
           ENDIF  
79    
80            thick_bd = h_bdry_values_SI (i,j,bi,bj)            IF (STREAMICE_ufacemask(i,j,bi,bj).eq.3.0) THEN
81            IF (thick_bd .ne. 0. _d 0) THEN             BCMASKX(i,j,bi,bj) = 3.0
82              h_after_uflux_SI (i,j,bi,bj) = thick_bd             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
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
109    
110  !       PRINT *, "H in last row ", H_streamice(81,20,1,1)            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               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    
       CALL STREAMICE_ADVECT_THICKNESS_X ( myThid,  
      O   hflux_x_SI,  
      O   h_after_uflux_SI,  
      I   time_step_loc )  
128    
 !       PRINT *, "H in last row ", H_streamice(81,20,1,1)  
129    
130        DO bj=myByLo(myThid),myByHi(myThid)             vtrans(i,j,bi,bj) = 1.0
131         DO bi=myBxLo(myThid),myBxHi(myThid)            ELSEIF (.not.(
132          DO j=1-OLy,sNy+OLy       &     STREAMICE_hmask(i,j,bi,bj).eq.1.0.OR.
133           DO i=1-OLx,sNx+OLx       &     STREAMICE_hmask(i,j-1,bi,bj).eq.1.0)) THEN
134            h_after_vflux_SI (i,j,bi,bj) =             BCMASKY(i,j,bi,bj) = 0.0
135       &     h_after_uflux_SI (i,j,bi,bj)             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
143    
144    
145           ENDDO           ENDDO
146          ENDDO          ENDDO
147         ENDDO         ENDDO
148        ENDDO        ENDDO
149    
150        CALL STREAMICE_ADVECT_THICKNESS_Y ( myThid,        _EXCH_XY_RL(utrans,myThid)
151       O   hflux_y_SI,        _EXCH_XY_RL(vtrans,myThid)
152       O   h_after_vflux_SI,        _EXCH_XY_RS(BCMASKx,myThid)
153       I   time_step_loc )        _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  !       PRINT *, "H in last row ", H_streamice(81,20,1,1)        max_icfl = 1.e-20
162    
163        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
164         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
165          DO j=1-OLy,sNy+OLy          DO j=1,sNy
166           DO i=1-OLx,sNx+OLx           DO i=1,sNx
167            IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN            IF (streamice_hmask(i,j,bi,bj).eq.1.0) THEN
168             H_streamice (i,j,bi,bj) =             loc_icfl=max(abs(utrans(i,j,bi,bj)),
169       &      h_after_vflux_SI (i,j,bi,bj)       &                  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            ENDIF
176           ENDDO           ENDDO
177          ENDDO          ENDDO
178         ENDDO         ENDDO
179        ENDDO        ENDDO
180    
181  !       PRINT *, "H in last row ", H_streamice(81,20,1,1)        CALL GLOBAL_MAX_R8 (max_icfl, myThid)
182    
183        CALL STREAMICE_ADV_FRONT ( myThid, time_step_loc )  #endif
184    
 !       PRINT *, "H in last row ", H_streamice(81,20,1,1)  
185    
186        _EXCH_XY_RL( H_streamice, myThid )  #ifdef ALLOW_AUTODIFF_TAMC
187        _EXCH_XY_RL( area_shelf_streamice, myThid )  CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics
188        _EXCH_XY_RL( STREAMICE_hmask, myThid )  CADJ STORE H_streamice  = comlev1, key=ikey_dynamics
189          #endif
       PRINT *, "END STREAMICE_ADVECT_THICKNESS"  
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
       RETURN  
       END SUBROUTINE STREAMICE_ADVECT_THICKNESS  
200    
201  !     NEED TO ADD SOME SORT OF CHECK THAT THE HALOS ARE LARGE ENOUGH        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    
       SUBROUTINE STREAMICE_ADVECT_THICKNESS_X ( myThid ,      
      O   hflux_x ,  
      O   h ,  
      I   time_step )  
209    
       IMPLICIT NONE  
210    
211  C     O   hflux_x ! flux per unit width across face        DO bj=myByLo(myThid),myByHi(myThid)
212  C     O   h         DO bi=myBxLo(myThid),myBxHi(myThid)
213  C     I   time_step          DO j=1-3,sNy+3
214             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  C     === Global variables ===            IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
221  #include "SIZE.h"             h_after_uflux_SI (i,j,bi,bj) = H_streamice(i,j,bi,bj) -
222  #include "GRID.h"       &      (hflux_x_SI(i+1,j,bi,bj)*dyG(i+1,j,bi,bj) -
223  #include "EEPARAMS.h"       &        hflux_x_SI(i,j,bi,bj)*dyG(i,j,bi,bj))
224  #include "PARAMS.h"       &        * recip_rA (i,j,bi,bj) * time_step_loc
225  #include "STREAMICE.h"             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        INTEGER myThid            ENDIF
232        _RL hflux_x          (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)           ENDDO
233        _RL h                (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)            ENDDO
234        _RL time_step         ENDDO
235                ENDDO
 #ifdef ALLOW_STREAMICE  
236    
 C     LOCAL VARIABLES  
237    
238        INTEGER i, j, bi, bj, Gi, Gj, k  #ifdef ALLOW_AUTODIFF_TAMC
239        _RL uface, phi  CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics
240        _RL stencil (-1:1)  #endif
241        LOGICAL H0_valid  ! there are valid cells to calculate a  
242                          ! slope-limited 2nd order flux  !      CALL STREAMICE_ADVECT_THICKNESS_Y ( myThid,
243        _RL SLOPE_LIMITER  !     O   hflux_y_SI,
244        _RL total_vol_out  !     O   h_after_vflux_SI,
245        external SLOPE_LIMITER  !     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,
253         I   time_step_loc )
254    
        total_vol_out = 0.0  
255    
256        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
257         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
258          DO j=1-3,sNy+3          DO j=1,sNy
259           Gj = (myYGlobalLo-1)+(bj-1)*sNy+j           DO i=1,sNx
260           IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN            Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
261            DO i=1-2,sNx+3            Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
262  C        THESE ARRAY BOUNDS INSURE THAT AFTER THIS STEP,            IF (((Gj .ge. 1) .and. (Gj .le. Ny))
263  C        VALUES WILL BE RELIABLE 2 GRID CELLS OUT IN THE       &       .or.STREAMICE_EW_PERIODIC) THEN
 C        X DIRECTION AND 3 CELLS OUT IN THE Y DIR  
            IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .or.  
      &         ((STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) .and.  
      &          (STREAMICE_hmask(i,j,bi,bj).ne.1.0))) THEN  
   
             Gi = (myXGlobalLo-1)+(bi-1)*sNx+i  
               
             IF (STREAMICE_ufacemask(i,j,bi,bj).eq.4.0) THEN  
              hflux_x (i,j,bi,bj) = u_flux_bdry_SI (i,j,bi,bj)  
             ELSE  
             
              uface = .5 *  
      &        (U_streamice(i,j,bi,bj)+U_streamice(i,j+1,bi,bj))  
   
   
              IF (uface .gt. 0. _d 0) THEN  
               DO k=-1,1  
                stencil (k) = h(i+k-1,j,bi,bj)  
               ENDDO  
               IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .and.  
      &            (STREAMICE_hmask(i-2,j,bi,bj).eq.1.0)) H0_valid=.true.  
               
               IF ((Gi.eq.1).and.(STREAMICE_hmask(i-1,j,bi,bj).eq.3.0))  
      &         THEN  ! we are at western bdry and there is a thick. bdry cond  
                hflux_x (i,j,bi,bj) = h(i-1,j,bi,bj) * uface  
               ELSEIF (H0_valid) THEN  
                phi = SLOPE_LIMITER (  
      &          stencil(0)-stencil(-1),  
      &          stencil(1)-stencil(0))  
                hflux_x (i,j,bi,bj) = uface *  
      &          (stencil(0) - phi * .5 * (stencil(0)-stencil(1)))  
               ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme  
                hflux_x (i,j,bi,bj) = uface * stencil(0)  
               ENDIF  
               
              ELSE ! uface <= 0  
   
               DO k=-1,1  
                stencil (k) = h(i-k,j,bi,bj)  
               ENDDO  
               IF ((STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) .and.  
      &            (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0)) H0_valid=.true.  
   
               IF ((Gi.eq.Nx).and.(STREAMICE_hmask(i+1,j,bi,bj).eq.3.0))  
      &         THEN  ! we are at western bdry and there is a thick. bdry cond  
                hflux_x (i,j,bi,bj) = h(i+1,j,bi,bj) * uface  
               ELSEIF (H0_valid) THEN  
                phi = SLOPE_LIMITER (  
      &          stencil(0)-stencil(-1),  
      &          stencil(1)-stencil(0))  
                hflux_x (i,j,bi,bj) = uface *  
      &          (stencil(0) - phi * .5 * (stencil(0)-stencil(1)))  
               ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme  
                hflux_x (i,j,bi,bj) = uface * stencil(0)  
               ENDIF  
   
              ENDIF  
   
             ENDIF  
   
             if (streamice_ufacemask(i,j,bi,bj).eq.2.0) THEN  
              total_vol_out = total_vol_out +  
      &         hflux_x (i,j,bi,bj) * time_step  
             ENDIF  
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             ENDIF
274            ENDDO  
275              ENDIF
276           ENDIF           ENDIF
277    
278             ENDDO
279          ENDDO          ENDDO
280         ENDDO         ENDDO
281        ENDDO        ENDDO
282    
 C     X-FLUXES AT CELL BOUNDARIES CALCULATED; NOW TAKE FLUX DIVERGENCE TO INCREMENT THICKNESS  
         
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-3,sNy+3          DO j=1,sNy
286           Gj = (myYGlobalLo-1)+(bj-1)*sNy+j           DO i=1,sNx
287           IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN            IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
288            DO i=1-2,sNx+2             H_streamice (i,j,bi,bj) =
289             IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN       &      h_after_vflux_SI (i,j,bi,bj)
290              h(i,j,bi,bj) = h(i,j,bi,bj) - time_step *            ENDIF
291       &       (hflux_x(i+1,j,bi,bj)*dyG(i+1,j,bi,bj) -           ENDDO
      &        hflux_x(i,j,bi,bj)*dyG(i,j,bi,bj)) *  
      &       recip_rA (i,j,bi,bj)  
            ENDIF  
           ENDDO  
          ENDIF  
292          ENDDO          ENDDO
293         ENDDO         ENDDO
294        ENDDO        ENDDO
295    
296  !       PRINT *, "TOTAL VOLUME OUT: ", total_vol_out  ! NOTE: AT THIS POINT H IS NOT VALID ON OVERLAP!!!
         
 #endif  
       RETURN  
       END SUBROUTINE STREAMICE_ADVECT_THICKNESS_X  
   
       SUBROUTINE STREAMICE_ADVECT_THICKNESS_Y ( myThid ,      
      O   hflux_y ,  
      O   h ,  
      I   time_step )  
   
       IMPLICIT NONE  
297    
298  C     O   hflux_y ! flux per unit width across face        if (streamice_move_front) then
299  C     O   h        CALL STREAMICE_ADV_FRONT (
300  C     I   time_step       &  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  C     === Global variables ===  #ifndef ALLOW_AUTODIFF_TAMC
317  #include "SIZE.h"        time_step_rem = time_step_rem - time_step_loc
318  #include "GRID.h"        enddo
319  #include "EEPARAMS.h"  #endif
 #include "PARAMS.h"  
 #include "STREAMICE.h"  
320    
       INTEGER myThid  
       _RL hflux_y          (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)  
       _RL h                (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)    
       _RL time_step  
         
 #ifdef ALLOW_STREAMICE  
321    
 C     LOCAL VARIABLES  
322    
323        INTEGER i, j, bi, bj, Gi, Gj, k  ! NOW WE APPLY MELT RATES !!
324        _RL vface, phi  ! THIS MAY BE MOVED TO A SEPARATE SUBROUTINE
325        _RL stencil (-1:1)        irho = 1./streamice_density
       LOGICAL H0_valid  ! there are valid cells to calculate a  
                         ! slope-limited 2nd order flux  
       _RL SLOPE_LIMITER  
       external SLOPE_LIMITER  
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-1,sNy+2          DO j=1,sNy
330           Gj = (myYGlobalLo-1)+(bj-1)*sNy+j           DO i=1,sNx
          DO i=1-2,sNx+2  
331            Gi = (myXGlobalLo-1)+(bi-1)*sNx+i            Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
332            IF ((Gi.ge.1) .and. (Gi.le.Nx)) THEN            Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
333  C         THESE ARRAY BOUNDS INSURE THAT AFTER THIS STEP,            IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0 .or.
334  C         VALUES WILL BE RELIABLE 1 GRID CELLS OUT IN THE       &       STREAMICE_hmask(i,j,bi,bj).eq.2.0) THEN
335  C         Y DIRECTION  
336             IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .or.             IF (STREAMICE_allow_cpl) THEN
337       &         ((STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) .and.  #ifdef ALLOW_SHELFICE
338       &          (STREAMICE_hmask(i,j,bi,bj).ne.1.0))) THEN              MR = -1. * (1.-float_frac_streamice(i,j,bi,bj)) *
339         &        shelfIceFreshWaterFlux(I,J,bi,bj) * irho *
340              IF (STREAMICE_vfacemask(i,j,bi,bj).eq.4.0) THEN       &        sec_per_year
341               hflux_y (i,j,bi,bj) = v_flux_bdry_SI (i,j,bi,bj)  
342              ELSE  #else
343                          STOP 'SHELFICE IS NOT ENABLED'
344               vface = .5 *  #endif      
345       &        (V_streamice(i,j,bi,bj)+V_streamice(i+1,j,bi,bj))             ELSE
346                  MR = (1.-float_frac_streamice(i,j,bi,bj)) *
347               IF (vface .gt. 0. _d 0) THEN       &             (BDOT_STREAMICE(i,j,bi,bj) +
348                DO k=-1,1       &              BDOT_pert(i,j,bi,bj))
349                 stencil (k) = h(i,j+k-1,bi,bj)             ENDIF
               ENDDO  
               IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .and.  
      &            (STREAMICE_hmask(i,j-2,bi,bj).eq.1.0)) H0_valid=.true.  
               
               IF ((Gj.eq.1).and.(STREAMICE_hmask(i,j-1,bi,bj).eq.3.0))  
      &         THEN  ! we are at western bdry and there is a thick. bdry cond  
                hflux_y (i,j,bi,bj) = h(i,j-1,bi,bj) * vface  
               ELSEIF (H0_valid) THEN  
                phi = SLOPE_LIMITER (  
      &          stencil(0)-stencil(-1),  
      &          stencil(1)-stencil(0))  
                hflux_y (i,j,bi,bj) = vface *  
      &          (stencil(0) - phi * .5 * (stencil(0)-stencil(1)))  
               ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme  
                hflux_y (i,j,bi,bj) = vface * stencil(0)  
               ENDIF  
              ELSE ! uface <= 0  
               DO k=-1,1  
                stencil (k) = h(i,j-k,bi,bj)  
               ENDDO  
               IF ((STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) .and.  
      &            (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0)) H0_valid=.true.  
   
               IF ((Gj.eq.Ny).and.(STREAMICE_hmask(i,j+1,bi,bj).eq.3.0))  
      &         THEN  ! we are at western bdry and there is a thick. bdry cond  
                hflux_y (i,j,bi,bj) = h(i,j+1,bi,bj) * vface  
               ELSEIF (H0_valid) THEN  
                phi = SLOPE_LIMITER (  
      &          stencil(0)-stencil(-1),  
      &          stencil(1)-stencil(0))  
                hflux_y (i,j,bi,bj) = vface *  
      &          (stencil(0) - phi * .5 * (stencil(0)-stencil(1)))  
               ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme  
                hflux_y (i,j,bi,bj) = vface * stencil(0)  
               ENDIF  
   
              ENDIF ! uface  0  
350    
351              ENDIF             SMB = ADOT_STREAMICE(i,j,bi,bj)
352               TMB = SMB - MR
353               IF ((TMB.lt.0.0) .and.
354         &         (MR * time_step_loc .gt.
355         &          H_streamice (i,j,bi,bj))) THEN
356                    H_streamice (i,j,bi,bj) = 0. _d 0
357                    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         &          H_streamice (i,j,bi,bj) + TMB * time_step_loc
362             ENDIF             ENDIF
           ENDIF  
          ENDDO  
         ENDDO  
        ENDDO  
       ENDDO  
363    
 C     X-FLUXES AT CELL BOUNDARIES CALCULATED; NOW TAKE FLUX DIVERGENCE TO INCREMENT THICKNESS  
364    
         
         
       DO bj=myByLo(myThid),myByHi(myThid)  
        DO bi=myBxLo(myThid),myBxHi(myThid)  
         DO j=1-1,sNy+1  
          DO i=1-2,sNx+2  
           Gi = (myXGlobalLo-1)+(bi-1)*sNx+i  
           IF ((Gi .ge. 1) .and. (Gi .le. Nx)) THEN  
            IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN  
             h(i,j,bi,bj) = h(i,j,bi,bj) - time_step *  
      &       (hflux_y(i,j+1,bi,bj)*dxG(i,j+1,bi,bj) -  
      &        hflux_y(i,j,bi,bj)*dxG(i,j,bi,bj)) *  
      &       recip_rA (i,j,bi,bj)  
            ENDIF  
365            ENDIF            ENDIF
366           ENDDO           ENDDO
367          ENDDO          ENDDO
368         ENDDO         ENDDO
369        ENDDO        ENDDO
370    
371  !       CALL WRITE_FLD_XY_RL ("h_after_yflux","",        _EXCH_XY_RL(H_streamice,myThid)
 !      & h, 0, myThid)  
372    
 #endif  
       RETURN  
       END SUBROUTINE STREAMICE_ADVECT_THICKNESS_Y  
373    
374          WRITE(msgBuf,'(A)') 'END STREAMICE_ADVECT_THICKNESS'
375           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
376         &                     SQUEEZE_RIGHT , 1)
377    
378    #endif
379          END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22