dgoldberg/streamice/streamice_advect_thickness_x.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-|--+----|

      SUBROUTINE STREAMICE_ADVECT_THICKNESS_X ( myThid ,     
     O   hflux_x ,
     O   h , 
     I   time_step )

      IMPLICIT NONE

C     O   hflux_x ! flux per unit width across face
C     O   h
C     I   time_step

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

      INTEGER myThid
      _RL hflux_x          (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

C     LOCAL VARIABLES

      INTEGER i, j, bi, bj, Gi, Gj, k
      _RL uface, phi
      _RL stencil (-1:1)
      LOGICAL H0_valid(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
                        ! there are valid cells to calculate a 
                        ! slope-limited 2nd order flux
      _RL SLOPE_LIMITER
      _RL total_vol_out
      external SLOPE_LIMITER

      total_vol_out = 0.0

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-oly,sNy+oly
         DO i=1-olx,sNx+olx
          H0_valid(i,j,bi,bj)=.false.
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-2,sNy+2
         Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
         IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN
          DO i=1-1,sNx+2
C        THESE ARRAY BOUNDS INSURE THAT AFTER THIS STEP,
C        VALUES WILL BE RELIABLE 2 GRID CELLS OUT IN THE
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(i,j,bi,bj)=.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(i,j,bi,bj)) 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(i,j,bi,bj)=.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(i,j,bi,bj)) 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

           ENDIF
          ENDDO
         ENDIF
        ENDDO
       ENDDO
      ENDDO

C     X-FLUXES AT CELL BOUNDARIES CALCULATED; NOW TAKE FLUX DIVERGENCE TO INCREMENT THICKNESS
      
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-3,sNy+3
         Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
         IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN
          DO i=1-2,sNx+2
           IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
            h(i,j,bi,bj) = h(i,j,bi,bj) - time_step * 
     &       (hflux_x(i+1,j,bi,bj)*dyG(i+1,j,bi,bj) -
     &        hflux_x(i,j,bi,bj)*dyG(i,j,bi,bj)) *
     &       recip_rA (i,j,bi,bj)
           ENDIF
          ENDDO
         ENDIF
        ENDDO
       ENDDO
      ENDDO

!       PRINT *, "TOTAL VOLUME OUT: ", total_vol_out
       
#endif
      RETURN
      END SUBROUTINE STREAMICE_ADVECT_THICKNESS_X

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			SUBROUTINE STREAMICE_ADVECT_THICKNESS_X ( myThid ,
9			O hflux_x ,
10			O h ,
11			I time_step )
12
13			IMPLICIT NONE
14
15			C O hflux_x ! flux per unit width across face
16			C O h
17			C I time_step
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
26			INTEGER myThid
27			_RL hflux_x (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
28			_RL h (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
29			_RL time_step
30
31			#ifdef ALLOW_STREAMICE
32
33			C LOCAL VARIABLES
34
35			INTEGER i, j, bi, bj, Gi, Gj, k
36			_RL uface, phi
37			_RL stencil (-1:1)
38	heimbach	1.2	LOGICAL H0_valid(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
39			! there are valid cells to calculate a
40	heimbach	1.1	! slope-limited 2nd order flux
41			_RL SLOPE_LIMITER
42			_RL total_vol_out
43			external SLOPE_LIMITER
44
45	heimbach	1.2	total_vol_out = 0.0
46	heimbach	1.1
47			DO bj=myByLo(myThid),myByHi(myThid)
48			DO bi=myBxLo(myThid),myBxHi(myThid)
49	heimbach	1.2	DO j=1-oly,sNy+oly
50			DO i=1-olx,sNx+olx
51			H0_valid(i,j,bi,bj)=.false.
52			ENDDO
53			ENDDO
54			ENDDO
55			ENDDO
56
57			DO bj=myByLo(myThid),myByHi(myThid)
58			DO bi=myBxLo(myThid),myBxHi(myThid)
59			DO j=1-2,sNy+2
60	heimbach	1.1	Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
61			IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN
62	heimbach	1.2	DO i=1-1,sNx+2
63	heimbach	1.1	C THESE ARRAY BOUNDS INSURE THAT AFTER THIS STEP,
64			C VALUES WILL BE RELIABLE 2 GRID CELLS OUT IN THE
65			C X DIRECTION AND 3 CELLS OUT IN THE Y DIR
66			IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .or.
67			& ((STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) .and.
68			& (STREAMICE_hmask(i,j,bi,bj).ne.1.0))) THEN
69
70			Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
71
72			IF (STREAMICE_ufacemask(i,j,bi,bj).eq.4.0) THEN
73			hflux_x (i,j,bi,bj) = u_flux_bdry_SI (i,j,bi,bj)
74			ELSE
75
76			uface = .5 *
77			& (U_streamice(i,j,bi,bj)+U_streamice(i,j+1,bi,bj))
78
79
80			IF (uface .gt. 0. _d 0) THEN
81			DO k=-1,1
82			stencil (k) = h(i+k-1,j,bi,bj)
83			ENDDO
84			IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .and.
85	heimbach	1.2	& (STREAMICE_hmask(i-2,j,bi,bj).eq.1.0))
86			& H0_valid(i,j,bi,bj)=.true.
87	heimbach	1.1
88			IF ((Gi.eq.1).and.(STREAMICE_hmask(i-1,j,bi,bj).eq.3.0))
89			& THEN ! we are at western bdry and there is a thick. bdry cond
90			hflux_x (i,j,bi,bj) = h(i-1,j,bi,bj) * uface
91	heimbach	1.2	ELSEIF (H0_valid(i,j,bi,bj)) THEN
92	heimbach	1.1	phi = SLOPE_LIMITER (
93			& stencil(0)-stencil(-1),
94			& stencil(1)-stencil(0))
95			hflux_x (i,j,bi,bj) = uface *
96			& (stencil(0) - phi * .5 * (stencil(0)-stencil(1)))
97			ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme
98			hflux_x (i,j,bi,bj) = uface * stencil(0)
99			ENDIF
100
101			ELSE ! uface <= 0
102
103			DO k=-1,1
104			stencil (k) = h(i-k,j,bi,bj)
105			ENDDO
106			IF ((STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) .and.
107	heimbach	1.2	& (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0))
108			& H0_valid(i,j,bi,bj)=.true.
109	heimbach	1.1
110			IF ((Gi.eq.Nx).and.(STREAMICE_hmask(i+1,j,bi,bj).eq.3.0))
111			& THEN ! we are at western bdry and there is a thick. bdry cond
112			hflux_x (i,j,bi,bj) = h(i+1,j,bi,bj) * uface
113	heimbach	1.2	ELSEIF (H0_valid(i,j,bi,bj)) THEN
114	heimbach	1.1	phi = SLOPE_LIMITER (
115			& stencil(0)-stencil(-1),
116			& stencil(1)-stencil(0))
117			hflux_x (i,j,bi,bj) = uface *
118			& (stencil(0) - phi * .5 * (stencil(0)-stencil(1)))
119			ELSE ! one of the two cells needed for a HO scheme is missing, use FO scheme
120			hflux_x (i,j,bi,bj) = uface * stencil(0)
121			ENDIF
122
123			ENDIF
124
125			ENDIF
126
127			if (streamice_ufacemask(i,j,bi,bj).eq.2.0) THEN
128			total_vol_out = total_vol_out +
129			& hflux_x (i,j,bi,bj) * time_step
130			ENDIF
131
132			ENDIF
133			ENDDO
134			ENDIF
135			ENDDO
136			ENDDO
137			ENDDO
138
139			C X-FLUXES AT CELL BOUNDARIES CALCULATED; NOW TAKE FLUX DIVERGENCE TO INCREMENT THICKNESS
140
141			DO bj=myByLo(myThid),myByHi(myThid)
142			DO bi=myBxLo(myThid),myBxHi(myThid)
143			DO j=1-3,sNy+3
144			Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
145			IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN
146			DO i=1-2,sNx+2
147			IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
148			h(i,j,bi,bj) = h(i,j,bi,bj) - time_step *
149			& (hflux_x(i+1,j,bi,bj)*dyG(i+1,j,bi,bj) -
150			& hflux_x(i,j,bi,bj)dyG(i,j,bi,bj))
151			& recip_rA (i,j,bi,bj)
152			ENDIF
153			ENDDO
154			ENDIF
155			ENDDO
156			ENDDO
157			ENDDO
158
159			! PRINT *, "TOTAL VOLUME OUT: ", total_vol_out
160
161			#endif
162			RETURN
163			END SUBROUTINE STREAMICE_ADVECT_THICKNESS_X
164