/[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.7 by heimbach, Thu Oct 4 03:45:53 2012 UTC
# Line 23  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    
30        INTEGER myThid        INTEGER myThid
31        _RL time_step        _RL time_step
32    
33  #ifdef ALLOW_STREAMICE  #ifdef ALLOW_STREAMICE
34    
35        INTEGER i, j, bi, bj, k, n_flux, iter_count, iter_flag        INTEGER i, j, bi, bj, k, iter_count, iter_rpt
36        INTEGER Gi, Gj        INTEGER Gi, Gj
37        INTEGER new_partial(4)        INTEGER new_partial(4)
38          INTEGER ikey_front, ikey_1
39          _RL iter_flag
40          _RL n_flux_1, n_flux_2
41        _RL href, rho, partial_vol, tot_flux, hpot        _RL href, rho, partial_vol, tot_flux, hpot
42          CHARACTER*(MAX_LEN_MBUF) msgBuf
43    
44        rho = streamice_density        rho = streamice_density
45        iter_count = 0  cph      iter_count = 0
46        iter_flag = 1        iter_flag = 1. _d 0
47          iter_rpt = 0
48    
49          DO iter_count = 0, 3
50    
51    #ifdef ALLOW_AUTODIFF_TAMC
52             ikey_front = (ikey_dynamics-1)*4 + iter_count + 1
53    CADJ STORE area_shelf_streamice
54    CADJ &     = comlev1_stream_front, key = ikey_front
55    CADJ STORE h_streamice
56    CADJ &     = comlev1_stream_front, key = ikey_front
57    CADJ STORE hflux_x_si
58    CADJ &     = comlev1_stream_front, key = ikey_front
59    CADJ STORE hflux_x_si2
60    CADJ &     = comlev1_stream_front, key = ikey_front
61    CADJ STORE hflux_y_si
62    CADJ &     = comlev1_stream_front, key = ikey_front
63    CADJ STORE hflux_y_si2
64    CADJ &     = comlev1_stream_front, key = ikey_front
65    CADJ STORE streamice_hmask
66    CADJ &     = comlev1_stream_front, key = ikey_front
67    CADJ STORE iter_flag
68    CADJ &     = comlev1_stream_front, key = ikey_front
69    #endif
70    
71        DO WHILE (iter_flag .eq. 1)         IF ( iter_flag .GT. 0. ) THEN
72          
73         iter_flag = 0         iter_flag = 0. _d 0
74    
75         IF (iter_count .gt. 0) then         IF (iter_count .gt. 0) then
76          DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
# Line 57  C     === Global variables === Line 87  C     === Global variables ===
87          ENDDO                  ENDDO        
88         ENDIF         ENDIF
89    
90         iter_count = iter_count + 1  !       iter_count = iter_count + 1
91           iter_rpt = iter_rpt + 1
92    
93         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
94          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
95    
96           DO j=1-1,sNy+1           DO j=1-1,sNy+1
97            Gj = (myYGlobalLo-1)+(bj-1)*sNy+j            Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
98            IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN            IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN
99             DO i=1-1,sNx+1             DO i=1-1,sNx+1
100              Gi = (myXGlobalLo-1)+(bi-1)*sNx+i              Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
101              IF ((Gi .ge. 1) .and. (Gi .le. Nx) .and.  
102    #ifdef ALLOW_AUTODIFF_TAMC
103              act1 = bi - myBxLo(myThid)
104              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
105              act2 = bj - myByLo(myThid)
106              max2 = myByHi(myThid) - myByLo(myThid) + 1
107              act3 = myThid - 1
108              max3 = nTx*nTy
109              act4 = ikey_front - 1
110              ikey_1 = i + 1
111         &         + (sNx+2)*(j)
112         &         + (sNx+2)*(sNy+2)*act1
113         &         + (sNx+2)*(sNy+2)*max1*act2
114         &         + (sNx+2)*(sNy+2)*max1*max2*act3
115         &         + (sNx+2)*(sNy+2)*max1*max2*max3*act4
116    CADJ STORE area_shelf_streamice(i,j,bi,bj)
117    CADJ &     = comlev1_stream_ij, key = ikey_1
118    CADJ STORE h_streamice(i,j,bi,bj)
119    CADJ &     = comlev1_stream_ij, key = ikey_1
120    CADJ STORE hflux_x_si(i,j,bi,bj)
121    CADJ &     = comlev1_stream_ij, key = ikey_1
122    CADJ STORE hflux_y_si(i,j,bi,bj)
123    CADJ &     = comlev1_stream_ij, key = ikey_1
124    CADJ STORE streamice_hmask(i,j,bi,bj)
125    CADJ &     = comlev1_stream_ij, key = ikey_1
126    #endif
127    
128                IF ((Gi .ge. 1) .and. (Gi .le. Nx) .and.
129       &          (STREAMICE_Hmask(i,j,bi,bj).eq.0.0 .or.       &          (STREAMICE_Hmask(i,j,bi,bj).eq.0.0 .or.
130       &           STREAMICE_Hmask(i,j,bi,bj).eq.2.0)) THEN       &           STREAMICE_Hmask(i,j,bi,bj).eq.2.0)) THEN
131               n_flux = 0               n_flux_1 = 0. _d 0
132               href = 0. _d 0               href = 0. _d 0
133               tot_flux = 0. _d 0               tot_flux = 0. _d 0
134    
135    #ifdef ALLOW_AUTODIFF_TAMC
136    CADJ STORE hflux_x_SI(i,j,bi,bj)
137    CADJ &     = comlev1_stream_ij, key = ikey_1
138    #endif
139               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
140                n_flux = n_flux + 1                n_flux_1 = n_flux_1 + 1. _d 0
141                href = href + H_streamice(i-1,j,bi,bj)                href = href + H_streamice(i-1,j,bi,bj)
142                tot_flux = tot_flux + hflux_x_SI(i,j,bi,bj) *                tot_flux = tot_flux + hflux_x_SI(i,j,bi,bj) *
143       &         dxG(i,j,bi,bj) * time_step       &         dxG(i,j,bi,bj) * time_step
144                hflux_x_SI(i,j,bi,bj) = 0. _d 0                hflux_x_SI(i,j,bi,bj) = 0. _d 0
145               ENDIF               ENDIF
146    
147    #ifdef ALLOW_AUTODIFF_TAMC
148    CADJ STORE hflux_x_SI(i,j,bi,bj)
149    CADJ &     = comlev1_stream_ij, key = ikey_1
150    #endif
151               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
152                n_flux = n_flux + 1                n_flux_1 = n_flux_1 + 1. _d 0
153                href = href + H_streamice(i+1,j,bi,bj)                href = href + H_streamice(i+1,j,bi,bj)
154                tot_flux = tot_flux - hflux_x_SI(i+1,j,bi,bj) *                tot_flux = tot_flux - hflux_x_SI(i+1,j,bi,bj) *
155       &         dxG(i+1,j,bi,bj) * time_step       &         dxG(i+1,j,bi,bj) * time_step
156                hflux_x_SI(i+1,j,bi,bj) = 0. _d 0                hflux_x_SI(i+1,j,bi,bj) = 0. _d 0
157               ENDIF               ENDIF
158    
159    #ifdef ALLOW_AUTODIFF_TAMC
160    CADJ STORE hflux_y_SI(i,j,bi,bj)
161    CADJ &     = comlev1_stream_ij, key = ikey_1
162    #endif
163               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
164                n_flux = n_flux + 1                n_flux_1 = n_flux_1 + 1. _d 0
165                href = href + H_streamice(i,j-1,bi,bj)                href = href + H_streamice(i,j-1,bi,bj)
166                tot_flux = tot_flux + hflux_y_SI(i,j,bi,bj) *                tot_flux = tot_flux + hflux_y_SI(i,j,bi,bj) *
167       &         dyG(i,j,bi,bj) * time_step       &         dyG(i,j,bi,bj) * time_step
168                hflux_y_SI(i,j,bi,bj) = 0. _d 0                hflux_y_SI(i,j,bi,bj) = 0. _d 0
169               ENDIF               ENDIF
170    
171    #ifdef ALLOW_AUTODIFF_TAMC
172    CADJ STORE hflux_y_SI(i,j,bi,bj)
173    CADJ &     = comlev1_stream_ij, key = ikey_1
174    #endif
175               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
176                n_flux = n_flux + 1                n_flux_1 = n_flux_1 + 1. _d 0
177                href = href + H_streamice(i,j+1,bi,bj)                href = href + H_streamice(i,j+1,bi,bj)
178                tot_flux = tot_flux - hflux_y_SI(i,j+1,bi,bj) *                tot_flux = tot_flux - hflux_y_SI(i,j+1,bi,bj) *
179       &         dyG(i,j+1,bi,bj) * time_step       &         dyG(i,j+1,bi,bj) * time_step
180                hflux_y_SI(i,j+1,bi,bj) = 0. _d 0                hflux_y_SI(i,j+1,bi,bj) = 0. _d 0
181               ENDIF               ENDIF
182    
183               IF (n_flux .gt. 0) THEN               IF (n_flux_1 .gt. 0.) THEN
184    
185                href = href / real(n_flux)                href = href / n_flux_1
186                partial_vol = H_streamice (i,j,bi,bj) *                partial_vol = H_streamice (i,j,bi,bj) *
187       &         area_shelf_streamice (i,j,bi,bj) + tot_flux       &         area_shelf_streamice (i,j,bi,bj) + tot_flux
188                hpot = partial_vol * recip_rA(i,j,bi,bj)                hpot = partial_vol * recip_rA(i,j,bi,bj)
# Line 119  C     === Global variables === Line 194  C     === Global variables ===
194       &           rA(i,j,bi,bj)       &           rA(i,j,bi,bj)
195                ELSEIF (hpot .lt. href) THEN ! cell still unfilled                ELSEIF (hpot .lt. href) THEN ! cell still unfilled
196    
197  !                PRINT *, "PARTIAL CELL INCREASED", tot_flux, i,  
 !      &         area_shelf_streamice (i,j,bi,bj),  
 !      &         H_streamice (i,j,bi,bj)  
198    
199                 STREAMICE_hmask (i,j,bi,bj) = 2.0                 STREAMICE_hmask (i,j,bi,bj) = 2.0
200                 area_shelf_streamice (i,j,bi,bj) = partial_vol / href                 area_shelf_streamice (i,j,bi,bj) = partial_vol / href
201                 H_streamice (i,j,bi,bj) = href                 H_streamice (i,j,bi,bj) = href
202                ELSE ! cell is filled - do overflow                ELSE ! cell is filled - do overflow
203    
 !                PRINT *, "CELL FILLED"  
204    
205                 STREAMICE_hmask (i,j,bi,bj) = 1.0                 STREAMICE_hmask (i,j,bi,bj) = 1.0
206                 area_shelf_streamice(i,j,bi,bj) =                 area_shelf_streamice(i,j,bi,bj) =
# Line 137  C     === Global variables === Line 209  C     === Global variables ===
209                                
210                 partial_vol = partial_vol - href * rA(i,j,bi,bj)                 partial_vol = partial_vol - href * rA(i,j,bi,bj)
211    
212                 iter_flag  = 1                 iter_flag  = 1. _d 0
213    
214                 n_flux = 0 ;                 n_flux_2 = 0. _d 0 ;
215                 DO k=1,4                 DO k=1,4
216                  new_partial (:) = 0                  new_partial (:) = 0
217                 ENDDO                 ENDDO
218    
219                 DO k=1,2                 DO k=1,2
220                  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
221                     n_flux = n_flux + 1                     n_flux_2 = n_flux_2 + 1. _d 0
222                  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
223                     n_flux = n_flux + 1                     n_flux_2 = n_flux_2 + 1. _d 0
224                     new_partial (k) = 1                     new_partial (k) = 1
225                  ENDIF                  ENDIF
226                 ENDDO                 ENDDO
227                 DO k=1,2                 DO k=1,2
228                  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
229                      n_flux = n_flux + 1                      n_flux_2 = n_flux_2 + 1. _d 0
230                  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
231                      n_flux = n_flux + 1                      n_flux_2 = n_flux_2 + 1. _d 0
232                      new_partial (k+2) = 1                      new_partial (k+2) = 1
233                  ENDIF                  ENDIF
234                 ENDDO                 ENDDO
235    
236                 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!
237                  H_streamice(i,j,bi,bj) = href + partial_vol *                  H_streamice(i,j,bi,bj) = href + partial_vol *
238       &             recip_rA(i,j,bi,bj)       &             recip_rA(i,j,bi,bj)
239                 ELSE                 ELSE
# Line 170  C     === Global variables === Line 242  C     === Global variables ===
242                  DO k=1,2                  DO k=1,2
243                   IF (new_partial(k) .eq. 1) THEN                   IF (new_partial(k) .eq. 1) THEN
244                    hflux_x_SI2(i-1+k,j,bi,bj) =                    hflux_x_SI2(i-1+k,j,bi,bj) =
245       &             partial_vol/time_step/real(n_flux)/       &             partial_vol/time_step/n_flux_2/
246       &               dxG(i-1+k,j,bi,bj)       &               dxG(i-1+k,j,bi,bj)
247                   ENDIF                   ENDIF
248                  ENDDO                  ENDDO
# Line 178  C     === Global variables === Line 250  C     === Global variables ===
250                  DO k=1,2                  DO k=1,2
251                   IF (new_partial(k+2) .eq. 1) THEN                   IF (new_partial(k+2) .eq. 1) THEN
252                    hflux_y_SI2(i,j-1+k,bi,bj) =                    hflux_y_SI2(i,j-1+k,bi,bj) =
253       &             partial_vol/time_step/real(n_flux)/       &             partial_vol/time_step/n_flux_2/
254       &               dxG(i,j-1+k,bi,bj)       &               dxG(i,j-1+k,bi,bj)
255                   ENDIF                   ENDIF
256                  ENDDO                  ENDDO
# Line 186  C     === Global variables === Line 258  C     === Global variables ===
258                 ENDIF                 ENDIF
259                ENDIF                ENDIF
260               ENDIF               ENDIF
261    
262              ENDIF              ENDIF
263             ENDDO             ENDDO
264            ENDIF            ENDIF
265           ENDDO           ENDDO
266    c
267          ENDDO          ENDDO
268         ENDDO         ENDDO
269    c
270          ENDIF
271        ENDDO        ENDDO
272    
273        IF (iter_count.gt.1) THEN        IF (iter_rpt.gt.1) THEN
274         PRINT *, "FRONT ADVANCE: ", iter_count, " ITERATIONS"         WRITE(msgBuf,'(A,I5,A)') 'FRONT ADVANCE: ',iter_rpt,
275         &  ' ITERATIONS'
276           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
277         &                     SQUEEZE_RIGHT , 1)      
278        ENDIF        ENDIF
279    
280    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22