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