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

Diff of /MITgcm_contrib/dgoldberg/streamice/streamice_adv_front.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.8 by dgoldberg, Sat Jun 8 22:15:33 2013 UTC
# Line 6  C $Name$ Line 6  C $Name$
6  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7    
8  CBOP  CBOP
9        SUBROUTINE STREAMICE_ADV_FRONT ( myThid, time_step )        SUBROUTINE STREAMICE_ADV_FRONT (
10         & myThid,
11         & time_step,
12         & hflux_x_si,
13         & hflux_y_si )
14    
15  C     /============================================================\  C     /============================================================\
16  C     | SUBROUTINE                                                 |    C     | SUBROUTINE                                                 |  
# Line 23  C     === Global variables === Line 27  C     === Global variables ===
27  #include "PARAMS.h"  #include "PARAMS.h"
28  #include "STREAMICE.h"  #include "STREAMICE.h"
29  #include "STREAMICE_ADV.h"  #include "STREAMICE_ADV.h"
30    #ifdef ALLOW_AUTODIFF_TAMC
31    # include "tamc.h"
32    #endif
33    
34        INTEGER myThid        INTEGER myThid
35        _RL time_step        _RL time_step
36          _RL hflux_x_SI (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
37          _RL hflux_y_SI (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
38    
39  #ifdef ALLOW_STREAMICE  #ifdef ALLOW_STREAMICE
40    
41        INTEGER i, j, bi, bj, k, n_flux, iter_count, iter_flag        INTEGER i, j, bi, bj, k, iter_count, iter_rpt
42        INTEGER Gi, Gj        INTEGER Gi, Gj
43        INTEGER new_partial(4)        INTEGER new_partial(4)
44          INTEGER ikey_front, ikey_1
45          _RL iter_flag
46          _RL n_flux_1, n_flux_2
47        _RL href, rho, partial_vol, tot_flux, hpot        _RL href, rho, partial_vol, tot_flux, hpot
48          CHARACTER*(MAX_LEN_MBUF) msgBuf
49          _RL hflux_x_SI2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
50          _RL hflux_y_SI2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
51    
52    
53        rho = streamice_density        rho = streamice_density
54        iter_count = 0  cph      iter_count = 0
55        iter_flag = 1        iter_flag = 1. _d 0
56          iter_rpt = 0
57    
58            DO bj=myByLo(myThid),myByHi(myThid)
59             DO bi=myBxLo(myThid),myBxHi(myThid)
60              DO j=1-OLy,sNy+OLy
61               DO i=1-OLx,sNx+OLx
62                hflux_x_SI2(i,j,bi,bj) = 0. _d 0
63                hflux_y_SI2(i,j,bi,bj) = 0. _d 0
64               ENDDO
65              ENDDO
66             ENDDO
67            ENDDO
68    
69    
70          DO iter_count = 0, 3
71    
72        DO WHILE (iter_flag .eq. 1)  #ifdef ALLOW_AUTODIFF_TAMC
73                   ikey_front = (ikey_dynamics-1)*4 + iter_count + 1
74         iter_flag = 0  CADJ STORE area_shelf_streamice
75    CADJ &     = comlev1_stream_front, key = ikey_front
76    CADJ STORE h_streamice
77    CADJ &     = comlev1_stream_front, key = ikey_front
78    CADJ STORE hflux_x_si
79    CADJ &     = comlev1_stream_front, key = ikey_front
80    CADJ STORE hflux_x_si2
81    CADJ &     = comlev1_stream_front, key = ikey_front
82    CADJ STORE hflux_y_si
83    CADJ &     = comlev1_stream_front, key = ikey_front
84    CADJ STORE hflux_y_si2
85    CADJ &     = comlev1_stream_front, key = ikey_front
86    CADJ STORE streamice_hmask
87    CADJ &     = comlev1_stream_front, key = ikey_front
88    CADJ STORE iter_flag
89    CADJ &     = comlev1_stream_front, key = ikey_front
90    #endif
91    
92           IF ( iter_flag .GT. 0. ) THEN
93    
94           iter_flag = 0. _d 0
95    
96         IF (iter_count .gt. 0) then         IF (iter_count .gt. 0) then
97          DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
# Line 57  C     === Global variables === Line 108  C     === Global variables ===
108          ENDDO                  ENDDO        
109         ENDIF         ENDIF
110    
111         iter_count = iter_count + 1  !       iter_count = iter_count + 1
112           iter_rpt = iter_rpt + 1
113    
114         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
115          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
116    
117           DO j=1-1,sNy+1           DO j=1-1,sNy+1
118            Gj = (myYGlobalLo-1)+(bj-1)*sNy+j            Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
119            IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN            IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN
120             DO i=1-1,sNx+1             DO i=1-1,sNx+1
121              Gi = (myXGlobalLo-1)+(bi-1)*sNx+i              Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
122              IF ((Gi .ge. 1) .and. (Gi .le. Nx) .and.  
123    #ifdef ALLOW_AUTODIFF_TAMC
124              act1 = bi - myBxLo(myThid)
125              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
126              act2 = bj - myByLo(myThid)
127              max2 = myByHi(myThid) - myByLo(myThid) + 1
128              act3 = myThid - 1
129              max3 = nTx*nTy
130              act4 = ikey_front - 1
131              ikey_1 = i + 1
132         &         + (sNx+2)*(j)
133         &         + (sNx+2)*(sNy+2)*act1
134         &         + (sNx+2)*(sNy+2)*max1*act2
135         &         + (sNx+2)*(sNy+2)*max1*max2*act3
136         &         + (sNx+2)*(sNy+2)*max1*max2*max3*act4
137    CADJ STORE area_shelf_streamice(i,j,bi,bj)
138    CADJ &     = comlev1_stream_ij, key = ikey_1
139    CADJ STORE h_streamice(i,j,bi,bj)
140    CADJ &     = comlev1_stream_ij, key = ikey_1
141    CADJ STORE hflux_x_si(i,j,bi,bj)
142    CADJ &     = comlev1_stream_ij, key = ikey_1
143    CADJ STORE hflux_y_si(i,j,bi,bj)
144    CADJ &     = comlev1_stream_ij, key = ikey_1
145    CADJ STORE streamice_hmask(i,j,bi,bj)
146    CADJ &     = comlev1_stream_ij, key = ikey_1
147    #endif
148    
149                IF ((Gi .ge. 1) .and. (Gi .le. Nx) .and.
150       &          (STREAMICE_Hmask(i,j,bi,bj).eq.0.0 .or.       &          (STREAMICE_Hmask(i,j,bi,bj).eq.0.0 .or.
151       &           STREAMICE_Hmask(i,j,bi,bj).eq.2.0)) THEN       &           STREAMICE_Hmask(i,j,bi,bj).eq.2.0)) THEN
152               n_flux = 0               n_flux_1 = 0. _d 0
153               href = 0. _d 0               href = 0. _d 0
154               tot_flux = 0. _d 0               tot_flux = 0. _d 0
155    
156    #ifdef ALLOW_AUTODIFF_TAMC
157    CADJ STORE hflux_x_SI(i,j,bi,bj)
158    CADJ &     = comlev1_stream_ij, key = ikey_1
159    #endif
160               IF (hflux_x_SI(i,j,bi,bj).gt. 0. _d 0) THEN               IF (hflux_x_SI(i,j,bi,bj).gt. 0. _d 0) THEN
161                n_flux = n_flux + 1                n_flux_1 = n_flux_1 + 1. _d 0
162                href = href + H_streamice(i-1,j,bi,bj)                href = href + H_streamice(i-1,j,bi,bj)
163                tot_flux = tot_flux + hflux_x_SI(i,j,bi,bj) *                tot_flux = tot_flux + hflux_x_SI(i,j,bi,bj) *
164       &         dxG(i,j,bi,bj) * time_step       &         dxG(i,j,bi,bj) * time_step
165                hflux_x_SI(i,j,bi,bj) = 0. _d 0                hflux_x_SI(i,j,bi,bj) = 0. _d 0
166               ENDIF               ENDIF
167    
168    #ifdef ALLOW_AUTODIFF_TAMC
169    CADJ STORE hflux_x_SI(i,j,bi,bj)
170    CADJ &     = comlev1_stream_ij, key = ikey_1
171    #endif
172               IF (hflux_x_SI(i+1,j,bi,bj).lt. 0. _d 0) THEN               IF (hflux_x_SI(i+1,j,bi,bj).lt. 0. _d 0) THEN
173                n_flux = n_flux + 1                n_flux_1 = n_flux_1 + 1. _d 0
174                href = href + H_streamice(i+1,j,bi,bj)                href = href + H_streamice(i+1,j,bi,bj)
175                tot_flux = tot_flux - hflux_x_SI(i+1,j,bi,bj) *                tot_flux = tot_flux - hflux_x_SI(i+1,j,bi,bj) *
176       &         dxG(i+1,j,bi,bj) * time_step       &         dxG(i+1,j,bi,bj) * time_step
177                hflux_x_SI(i+1,j,bi,bj) = 0. _d 0                hflux_x_SI(i+1,j,bi,bj) = 0. _d 0
178               ENDIF               ENDIF
179    
180    #ifdef ALLOW_AUTODIFF_TAMC
181    CADJ STORE hflux_y_SI(i,j,bi,bj)
182    CADJ &     = comlev1_stream_ij, key = ikey_1
183    #endif
184               IF (hflux_y_SI(i,j,bi,bj).gt. 0. _d 0) THEN               IF (hflux_y_SI(i,j,bi,bj).gt. 0. _d 0) THEN
185                n_flux = n_flux + 1                n_flux_1 = n_flux_1 + 1. _d 0
186                href = href + H_streamice(i,j-1,bi,bj)                href = href + H_streamice(i,j-1,bi,bj)
187                tot_flux = tot_flux + hflux_y_SI(i,j,bi,bj) *                tot_flux = tot_flux + hflux_y_SI(i,j,bi,bj) *
188       &         dyG(i,j,bi,bj) * time_step       &         dyG(i,j,bi,bj) * time_step
189                hflux_y_SI(i,j,bi,bj) = 0. _d 0                hflux_y_SI(i,j,bi,bj) = 0. _d 0
190               ENDIF               ENDIF
191    
192    #ifdef ALLOW_AUTODIFF_TAMC
193    CADJ STORE hflux_y_SI(i,j,bi,bj)
194    CADJ &     = comlev1_stream_ij, key = ikey_1
195    #endif
196               IF (hflux_y_SI(i,j+1,bi,bj).lt. 0. _d 0) THEN               IF (hflux_y_SI(i,j+1,bi,bj).lt. 0. _d 0) THEN
197                n_flux = n_flux + 1                n_flux_1 = n_flux_1 + 1. _d 0
198                href = href + H_streamice(i,j+1,bi,bj)                href = href + H_streamice(i,j+1,bi,bj)
199                tot_flux = tot_flux - hflux_y_SI(i,j+1,bi,bj) *                tot_flux = tot_flux - hflux_y_SI(i,j+1,bi,bj) *
200       &         dyG(i,j+1,bi,bj) * time_step       &         dyG(i,j+1,bi,bj) * time_step
201                hflux_y_SI(i,j+1,bi,bj) = 0. _d 0                hflux_y_SI(i,j+1,bi,bj) = 0. _d 0
202               ENDIF               ENDIF
203    
204               IF (n_flux .gt. 0) THEN               IF (n_flux_1 .gt. 0.) THEN
205    
206                href = href / real(n_flux)                href = href / n_flux_1
207                partial_vol = H_streamice (i,j,bi,bj) *                partial_vol = H_streamice (i,j,bi,bj) *
208       &         area_shelf_streamice (i,j,bi,bj) + tot_flux       &         area_shelf_streamice (i,j,bi,bj) + tot_flux
209                hpot = partial_vol * recip_rA(i,j,bi,bj)                hpot = partial_vol * recip_rA(i,j,bi,bj)
# Line 119  C     === Global variables === Line 215  C     === Global variables ===
215       &           rA(i,j,bi,bj)       &           rA(i,j,bi,bj)
216                ELSEIF (hpot .lt. href) THEN ! cell still unfilled                ELSEIF (hpot .lt. href) THEN ! cell still unfilled
217    
218  !                PRINT *, "PARTIAL CELL INCREASED", tot_flux, i,  
 !      &         area_shelf_streamice (i,j,bi,bj),  
 !      &         H_streamice (i,j,bi,bj)  
219    
220                 STREAMICE_hmask (i,j,bi,bj) = 2.0                 STREAMICE_hmask (i,j,bi,bj) = 2.0
221                 area_shelf_streamice (i,j,bi,bj) = partial_vol / href                 area_shelf_streamice (i,j,bi,bj) = partial_vol / href
222                 H_streamice (i,j,bi,bj) = href                 H_streamice (i,j,bi,bj) = href
223                ELSE ! cell is filled - do overflow                ELSE ! cell is filled - do overflow
224    
 !                PRINT *, "CELL FILLED"  
225    
226                 STREAMICE_hmask (i,j,bi,bj) = 1.0                 STREAMICE_hmask (i,j,bi,bj) = 1.0
227                 area_shelf_streamice(i,j,bi,bj) =                 area_shelf_streamice(i,j,bi,bj) =
# Line 137  C     === Global variables === Line 230  C     === Global variables ===
230                                
231                 partial_vol = partial_vol - href * rA(i,j,bi,bj)                 partial_vol = partial_vol - href * rA(i,j,bi,bj)
232    
233                 iter_flag  = 1                 iter_flag  = 1. _d 0
234    
235                 n_flux = 0 ;                 n_flux_2 = 0. _d 0 ;
236                 DO k=1,4                 DO k=1,4
237                  new_partial (:) = 0                  new_partial (:) = 0
238                 ENDDO                 ENDDO
239    
240                 DO k=1,2                 DO k=1,2
241                  IF (STREAMICE_ufacemask(i-1+k,j,bi,bj).eq.2.0) THEN  ! at a permanent calving boundary - no advance allowed                  IF (STREAMICE_ufacemask(i-1+k,j,bi,bj).eq.2.0) THEN  ! at a permanent calving boundary - no advance allowed
242                     n_flux = n_flux + 1                     n_flux_2 = n_flux_2 + 1. _d 0
243                  ELSEIF (STREAMICE_hmask(i+2*k-3,j,bi,bj).eq.0 _d 0) THEN ! adjacent cell is completely ice free                  ELSEIF (STREAMICE_hmask(i+2*k-3,j,bi,bj).eq.0 _d 0) THEN ! adjacent cell is completely ice free
244                     n_flux = n_flux + 1                     n_flux_2 = n_flux_2 + 1. _d 0
245                     new_partial (k) = 1                     new_partial (k) = 1
246                  ENDIF                  ENDIF
247                 ENDDO                 ENDDO
248                 DO k=1,2                 DO k=1,2
249                  IF (STREAMICE_vfacemask (i,j-1+k,bi,bj).eq.2.0) THEN                  IF (STREAMICE_vfacemask (i,j-1+k,bi,bj).eq.2.0) THEN
250                      n_flux = n_flux + 1                      n_flux_2 = n_flux_2 + 1. _d 0
251                  ELSEIF (STREAMICE_hmask(i,j+2*k-3,bi,bj).eq.0 _d 0) THEN                  ELSEIF (STREAMICE_hmask(i,j+2*k-3,bi,bj).eq.0 _d 0) THEN
252                      n_flux = n_flux + 1                      n_flux_2 = n_flux_2 + 1. _d 0
253                      new_partial (k+2) = 1                      new_partial (k+2) = 1
254                  ENDIF                  ENDIF
255                 ENDDO                 ENDDO
256    
257                 IF (n_flux .eq. 0) THEN ! there is nowhere to put the extra ice!                 IF (n_flux_2 .eq. 0.) THEN ! there is nowhere to put the extra ice!
258                  H_streamice(i,j,bi,bj) = href + partial_vol *                  H_streamice(i,j,bi,bj) = href + partial_vol *
259       &             recip_rA(i,j,bi,bj)       &             recip_rA(i,j,bi,bj)
260                 ELSE                 ELSE
# Line 170  C     === Global variables === Line 263  C     === Global variables ===
263                  DO k=1,2                  DO k=1,2
264                   IF (new_partial(k) .eq. 1) THEN                   IF (new_partial(k) .eq. 1) THEN
265                    hflux_x_SI2(i-1+k,j,bi,bj) =                    hflux_x_SI2(i-1+k,j,bi,bj) =
266       &             partial_vol/time_step/real(n_flux)/       &             partial_vol/time_step/n_flux_2/
267       &               dxG(i-1+k,j,bi,bj)       &               dxG(i-1+k,j,bi,bj)
268                   ENDIF                   ENDIF
269                  ENDDO                  ENDDO
# Line 178  C     === Global variables === Line 271  C     === Global variables ===
271                  DO k=1,2                  DO k=1,2
272                   IF (new_partial(k+2) .eq. 1) THEN                   IF (new_partial(k+2) .eq. 1) THEN
273                    hflux_y_SI2(i,j-1+k,bi,bj) =                    hflux_y_SI2(i,j-1+k,bi,bj) =
274       &             partial_vol/time_step/real(n_flux)/       &             partial_vol/time_step/n_flux_2/
275       &               dxG(i,j-1+k,bi,bj)       &               dxG(i,j-1+k,bi,bj)
276                   ENDIF                   ENDIF
277                  ENDDO                  ENDDO
# Line 186  C     === Global variables === Line 279  C     === Global variables ===
279                 ENDIF                 ENDIF
280                ENDIF                ENDIF
281               ENDIF               ENDIF
282    
283              ENDIF              ENDIF
284             ENDDO             ENDDO
285            ENDIF            ENDIF
286           ENDDO           ENDDO
287    c
288          ENDDO          ENDDO
289         ENDDO         ENDDO
290    c
291          ENDIF
292        ENDDO        ENDDO
293    
294        IF (iter_count.gt.1) THEN        IF (iter_rpt.gt.1) THEN
295         PRINT *, "FRONT ADVANCE: ", iter_count, " ITERATIONS"         WRITE(msgBuf,'(A,I5,A)') 'FRONT ADVANCE: ',iter_rpt,
296         &  ' ITERATIONS'
297           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
298         &                     SQUEEZE_RIGHT , 1)      
299        ENDIF        ENDIF
300    
301    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22