dgoldberg/streamice/streamice_adv_front.F

C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_init_varia.F,v 1.6 2011/06/29 16:24:10 dng Exp $
C $Name:  $

#include "STREAMICE_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP
      SUBROUTINE STREAMICE_ADV_FRONT ( myThid, time_step )

C     /============================================================\
C     | SUBROUTINE                                                 |   
C     | o                                                          |
C     |============================================================|
C     |                                                            |
C     \============================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "STREAMICE.h"
#include "STREAMICE_ADV.h"

      INTEGER myThid
      _RL time_step

#ifdef ALLOW_STREAMICE

      INTEGER i, j, bi, bj, k, n_flux, iter_count, iter_flag
      INTEGER Gi, Gj
      INTEGER new_partial(4)
      _RL href, rho, partial_vol, tot_flux, hpot

      rho = streamice_density
      iter_count = 0
      iter_flag = 1

      DO WHILE (iter_flag .eq. 1)
      
       iter_flag = 0

       IF (iter_count .gt. 0) then
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            hflux_x_SI(i,j,bi,bj)=hflux_x_SI2(i,j,bi,bj)
            hflux_y_SI(i,j,bi,bj)=hflux_y_SI2(i,j,bi,bj)
            hflux_x_SI2(i,j,bi,bj) = 0. _d 0
            hflux_y_SI2(i,j,bi,bj) = 0. _d 0
           ENDDO
          ENDDO
         ENDDO
        ENDDO        
       ENDIF

       iter_count = iter_count + 1

       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-1,sNy+1
          Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
          IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN
           DO i=1-1,sNx+1
            Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
            IF ((Gi .ge. 1) .and. (Gi .le. Nx) .and. 
     &          (STREAMICE_Hmask(i,j,bi,bj).eq.0.0 .or.
     &           STREAMICE_Hmask(i,j,bi,bj).eq.2.0)) THEN
             n_flux = 0
             href = 0. _d 0
             tot_flux = 0. _d 0

             IF (hflux_x_SI(i,j,bi,bj).gt. 0. _d 0) THEN
              n_flux = n_flux + 1
              href = href + H_streamice(i-1,j,bi,bj)
              tot_flux = tot_flux + hflux_x_SI(i,j,bi,bj) * 
     &         dxG(i,j,bi,bj) * time_step
              hflux_x_SI(i,j,bi,bj) = 0. _d 0
             ENDIF

             IF (hflux_x_SI(i+1,j,bi,bj).lt. 0. _d 0) THEN
              n_flux = n_flux + 1
              href = href + H_streamice(i+1,j,bi,bj)
              tot_flux = tot_flux - hflux_x_SI(i+1,j,bi,bj) * 
     &         dxG(i+1,j,bi,bj) * time_step
              hflux_x_SI(i+1,j,bi,bj) = 0. _d 0
             ENDIF

             IF (hflux_y_SI(i,j,bi,bj).gt. 0. _d 0) THEN
              n_flux = n_flux + 1
              href = href + H_streamice(i,j-1,bi,bj)
              tot_flux = tot_flux + hflux_y_SI(i,j,bi,bj) * 
     &         dyG(i,j,bi,bj) * time_step
              hflux_y_SI(i,j,bi,bj) = 0. _d 0
             ENDIF

             IF (hflux_y_SI(i,j+1,bi,bj).lt. 0. _d 0) THEN
              n_flux = n_flux + 1
              href = href + H_streamice(i,j+1,bi,bj)
              tot_flux = tot_flux - hflux_y_SI(i,j+1,bi,bj) * 
     &         dyG(i,j+1,bi,bj) * time_step
              hflux_y_SI(i,j+1,bi,bj) = 0. _d 0
             ENDIF
 
             IF (n_flux .gt. 0) THEN

              href = href / real(n_flux)
              partial_vol = H_streamice (i,j,bi,bj) * 
     &         area_shelf_streamice (i,j,bi,bj) + tot_flux
              hpot = partial_vol * recip_rA(i,j,bi,bj)
              
              IF (hpot .eq. href) THEN ! cell is exactly covered, no overflow
                STREAMICE_hmask (i,j,bi,bj) = 1.0
                H_streamice (i,j,bi,bj) = href
                area_shelf_streamice(i,j,bi,bj) = 
     &           rA(i,j,bi,bj)
              ELSEIF (hpot .lt. href) THEN ! cell still unfilled

!                PRINT *, "PARTIAL CELL INCREASED", tot_flux, i,
!      &         area_shelf_streamice (i,j,bi,bj),
!      &         H_streamice (i,j,bi,bj)

               STREAMICE_hmask (i,j,bi,bj) = 2.0
               area_shelf_streamice (i,j,bi,bj) = partial_vol / href
               H_streamice (i,j,bi,bj) = href
              ELSE ! cell is filled - do overflow

!                PRINT *, "CELL FILLED"

               STREAMICE_hmask (i,j,bi,bj) = 1.0
               area_shelf_streamice(i,j,bi,bj) = 
     &           rA(i,j,bi,bj)

               
               partial_vol = partial_vol - href * rA(i,j,bi,bj)

               iter_flag  = 1

               n_flux = 0 ; 
               DO k=1,4
                new_partial (:) = 0
               ENDDO

               DO k=1,2
                IF (STREAMICE_ufacemask(i-1+k,j,bi,bj).eq.2.0) THEN  ! at a permanent calving boundary - no advance allowed
                   n_flux = n_flux + 1
                ELSEIF (STREAMICE_hmask(i+2*k-3,j,bi,bj).eq.0 _d 0) THEN ! adjacent cell is completely ice free
                   n_flux = n_flux + 1
                   new_partial (k) = 1
                ENDIF
               ENDDO
               DO k=1,2
                IF (STREAMICE_vfacemask (i,j-1+k,bi,bj).eq.2.0) THEN
                    n_flux = n_flux + 1
                ELSEIF (STREAMICE_hmask(i,j+2*k-3,bi,bj).eq.0 _d 0) THEN
                    n_flux = n_flux + 1
                    new_partial (k+2) = 1
                ENDIF
               ENDDO

               IF (n_flux .eq. 0) THEN ! there is nowhere to put the extra ice!
                H_streamice(i,j,bi,bj) = href + partial_vol *
     &             recip_rA(i,j,bi,bj)
               ELSE
                H_streamice(i,j,bi,bj) = href

                DO k=1,2
                 IF (new_partial(k) .eq. 1) THEN
                  hflux_x_SI2(i-1+k,j,bi,bj) = 
     &             partial_vol/time_step/real(n_flux)/
     &               dxG(i-1+k,j,bi,bj)
                 ENDIF
                ENDDO

                DO k=1,2
                 IF (new_partial(k+2) .eq. 1) THEN
                  hflux_y_SI2(i,j-1+k,bi,bj) = 
     &             partial_vol/time_step/real(n_flux)/
     &               dxG(i,j-1+k,bi,bj)
                 ENDIF
                ENDDO

               ENDIF
              ENDIF
             ENDIF
            ENDIF
           ENDDO
          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      IF (iter_count.gt.1) THEN
       PRINT *, "FRONT ADVANCE: ", iter_count, " ITERATIONS"
      ENDIF


#endif
      RETURN
      END
1	heimbach	1.1	C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_init_varia.F,v 1.6 2011/06/29 16:24:10 dng Exp $
2			C $Name: $
3
4			#include "STREAMICE_OPTIONS.h"
5
6			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
7
8			CBOP
9			SUBROUTINE STREAMICE_ADV_FRONT ( myThid, time_step )
10
11			C /============================================================\
12			C \| SUBROUTINE \|
13			C \| o \|
14			C \|============================================================\|
15			C \| \|
16			C \============================================================/
17			IMPLICIT NONE
18
19			C === Global variables ===
20			#include "SIZE.h"
21			#include "GRID.h"
22			#include "EEPARAMS.h"
23			#include "PARAMS.h"
24			#include "STREAMICE.h"
25			#include "STREAMICE_ADV.h"
26
27			INTEGER myThid
28			_RL time_step
29
30			#ifdef ALLOW_STREAMICE
31
32			INTEGER i, j, bi, bj, k, n_flux, iter_count, iter_flag
33			INTEGER Gi, Gj
34			INTEGER new_partial(4)
35			_RL href, rho, partial_vol, tot_flux, hpot
36
37			rho = streamice_density
38			iter_count = 0
39			iter_flag = 1
40
41			DO WHILE (iter_flag .eq. 1)
42
43			iter_flag = 0
44
45			IF (iter_count .gt. 0) then
46			DO bj=myByLo(myThid),myByHi(myThid)
47			DO bi=myBxLo(myThid),myBxHi(myThid)
48			DO j=1-OLy,sNy+OLy
49			DO i=1-OLx,sNx+OLx
50			hflux_x_SI(i,j,bi,bj)=hflux_x_SI2(i,j,bi,bj)
51			hflux_y_SI(i,j,bi,bj)=hflux_y_SI2(i,j,bi,bj)
52			hflux_x_SI2(i,j,bi,bj) = 0. _d 0
53			hflux_y_SI2(i,j,bi,bj) = 0. _d 0
54			ENDDO
55			ENDDO
56			ENDDO
57			ENDDO
58			ENDIF
59
60			iter_count = iter_count + 1
61
62			DO bj=myByLo(myThid),myByHi(myThid)
63			DO bi=myBxLo(myThid),myBxHi(myThid)
64			DO j=1-1,sNy+1
65			Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
66			IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN
67			DO i=1-1,sNx+1
68			Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
69			IF ((Gi .ge. 1) .and. (Gi .le. Nx) .and.
70			& (STREAMICE_Hmask(i,j,bi,bj).eq.0.0 .or.
71			& STREAMICE_Hmask(i,j,bi,bj).eq.2.0)) THEN
72			n_flux = 0
73			href = 0. _d 0
74			tot_flux = 0. _d 0
75
76			IF (hflux_x_SI(i,j,bi,bj).gt. 0. _d 0) THEN
77			n_flux = n_flux + 1
78			href = href + H_streamice(i-1,j,bi,bj)
79			tot_flux = tot_flux + hflux_x_SI(i,j,bi,bj) *
80			& dxG(i,j,bi,bj) * time_step
81			hflux_x_SI(i,j,bi,bj) = 0. _d 0
82			ENDIF
83
84			IF (hflux_x_SI(i+1,j,bi,bj).lt. 0. _d 0) THEN
85			n_flux = n_flux + 1
86			href = href + H_streamice(i+1,j,bi,bj)
87			tot_flux = tot_flux - hflux_x_SI(i+1,j,bi,bj) *
88			& dxG(i+1,j,bi,bj) * time_step
89			hflux_x_SI(i+1,j,bi,bj) = 0. _d 0
90			ENDIF
91
92			IF (hflux_y_SI(i,j,bi,bj).gt. 0. _d 0) THEN
93			n_flux = n_flux + 1
94			href = href + H_streamice(i,j-1,bi,bj)
95			tot_flux = tot_flux + hflux_y_SI(i,j,bi,bj) *
96			& dyG(i,j,bi,bj) * time_step
97			hflux_y_SI(i,j,bi,bj) = 0. _d 0
98			ENDIF
99
100			IF (hflux_y_SI(i,j+1,bi,bj).lt. 0. _d 0) THEN
101			n_flux = n_flux + 1
102			href = href + H_streamice(i,j+1,bi,bj)
103			tot_flux = tot_flux - hflux_y_SI(i,j+1,bi,bj) *
104			& dyG(i,j+1,bi,bj) * time_step
105			hflux_y_SI(i,j+1,bi,bj) = 0. _d 0
106			ENDIF
107
108			IF (n_flux .gt. 0) THEN
109
110			href = href / real(n_flux)
111			partial_vol = H_streamice (i,j,bi,bj) *
112			& area_shelf_streamice (i,j,bi,bj) + tot_flux
113			hpot = partial_vol * recip_rA(i,j,bi,bj)
114
115			IF (hpot .eq. href) THEN ! cell is exactly covered, no overflow
116			STREAMICE_hmask (i,j,bi,bj) = 1.0
117			H_streamice (i,j,bi,bj) = href
118			area_shelf_streamice(i,j,bi,bj) =
119			& rA(i,j,bi,bj)
120			ELSEIF (hpot .lt. href) THEN ! cell still unfilled
121
122			! PRINT *, "PARTIAL CELL INCREASED", tot_flux, i,
123			! & area_shelf_streamice (i,j,bi,bj),
124			! & H_streamice (i,j,bi,bj)
125
126			STREAMICE_hmask (i,j,bi,bj) = 2.0
127			area_shelf_streamice (i,j,bi,bj) = partial_vol / href
128			H_streamice (i,j,bi,bj) = href
129			ELSE ! cell is filled - do overflow
130
131			! PRINT *, "CELL FILLED"
132
133			STREAMICE_hmask (i,j,bi,bj) = 1.0
134			area_shelf_streamice(i,j,bi,bj) =
135			& rA(i,j,bi,bj)
136
137
138			partial_vol = partial_vol - href * rA(i,j,bi,bj)
139
140			iter_flag = 1
141
142			n_flux = 0 ;
143			DO k=1,4
144			new_partial (:) = 0
145			ENDDO
146
147			DO k=1,2
148			IF (STREAMICE_ufacemask(i-1+k,j,bi,bj).eq.2.0) THEN ! at a permanent calving boundary - no advance allowed
149			n_flux = n_flux + 1
150			ELSEIF (STREAMICE_hmask(i+2*k-3,j,bi,bj).eq.0 _d 0) THEN ! adjacent cell is completely ice free
151			n_flux = n_flux + 1
152			new_partial (k) = 1
153			ENDIF
154			ENDDO
155			DO k=1,2
156			IF (STREAMICE_vfacemask (i,j-1+k,bi,bj).eq.2.0) THEN
157			n_flux = n_flux + 1
158			ELSEIF (STREAMICE_hmask(i,j+2*k-3,bi,bj).eq.0 _d 0) THEN
159			n_flux = n_flux + 1
160			new_partial (k+2) = 1
161			ENDIF
162			ENDDO
163
164			IF (n_flux .eq. 0) THEN ! there is nowhere to put the extra ice!
165			H_streamice(i,j,bi,bj) = href + partial_vol *
166			& recip_rA(i,j,bi,bj)
167			ELSE
168			H_streamice(i,j,bi,bj) = href
169
170			DO k=1,2
171			IF (new_partial(k) .eq. 1) THEN
172			hflux_x_SI2(i-1+k,j,bi,bj) =
173			& partial_vol/time_step/real(n_flux)/
174			& dxG(i-1+k,j,bi,bj)
175			ENDIF
176			ENDDO
177
178			DO k=1,2
179			IF (new_partial(k+2) .eq. 1) THEN
180			hflux_y_SI2(i,j-1+k,bi,bj) =
181			& partial_vol/time_step/real(n_flux)/
182			& dxG(i,j-1+k,bi,bj)
183			ENDIF
184			ENDDO
185
186			ENDIF
187			ENDIF
188			ENDIF
189			ENDIF
190			ENDDO
191			ENDIF
192			ENDDO
193			ENDDO
194			ENDDO
195			ENDDO
196
197			IF (iter_count.gt.1) THEN
198			PRINT *, "FRONT ADVANCE: ", iter_count, " ITERATIONS"
199			ENDIF
200
201
202
203			#endif
204			RETURN
205			END