C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/dgoldberg/streamice/streamice_advect_thickness.F,v 1.1 2012/03/29 15:59:21 heimbach Exp $
C $Name:  $

#include "STREAMICE_OPTIONS.h"

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


CBOP
      SUBROUTINE STREAMICE_ADVECT_THICKNESS ( 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
      _RL thick_bd
      _RL SLOPE_LIMITER
      _RL sec_per_year, time_step_loc
      external SLOPE_LIMITER

      sec_per_year = 365.*86400.

      time_step_loc = time_step / sec_per_year

      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) = 0. _d 0
          hflux_y_SI (i,j,bi,bj) = 0. _d 0
          hflux_x_SI2 (i,j,bi,bj) = 0. _d 0
          hflux_y_SI2 (i,j,bi,bj) = 0. _d 0
          IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
           h_after_uflux_SI (i,j,bi,bj) = 
     &      H_streamice (i,j,bi,bj)
          ENDIF

          thick_bd = h_bdry_values_SI (i,j,bi,bj)
          IF (thick_bd .ne. 0. _d 0) THEN
            h_after_uflux_SI (i,j,bi,bj) = thick_bd
          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO

!       PRINT *, "H in last row ", H_streamice(81,20,1,1)

      CALL STREAMICE_ADVECT_THICKNESS_X ( myThid,
     O   hflux_x_SI,
     O   h_after_uflux_SI,
     I   time_step_loc )

!       PRINT *, "H in last row ", H_streamice(81,20,1,1)

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          h_after_vflux_SI (i,j,bi,bj) = 
     &     h_after_uflux_SI (i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      CALL STREAMICE_ADVECT_THICKNESS_Y ( myThid,
     O   hflux_y_SI,
     O   h_after_vflux_SI,
     I   time_step_loc )

!       PRINT *, "H in last row ", H_streamice(81,20,1,1)

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
           H_streamice (i,j,bi,bj) = 
     &      h_after_vflux_SI (i,j,bi,bj)
          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO

!       PRINT *, "H in last row ", H_streamice(81,20,1,1)

      CALL STREAMICE_ADV_FRONT ( myThid, time_step_loc )

!       PRINT *, "H in last row ", H_streamice(81,20,1,1)

      _EXCH_XY_RL( H_streamice, myThid )
      _EXCH_XY_RL( area_shelf_streamice, myThid )
      _EXCH_XY_RL( STREAMICE_hmask, myThid )
      
      PRINT *, "END STREAMICE_ADVECT_THICKNESS"

#endif
      RETURN
      END SUBROUTINE STREAMICE_ADVECT_THICKNESS 

!     NEED TO ADD SOME SORT OF CHECK THAT THE HALOS ARE LARGE ENOUGH

      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  ! 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-3,sNy+3
         Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
         IF ((Gj .ge. 1) .and. (Gj .le. Ny)) THEN
          DO i=1-2,sNx+3
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=.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) 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=.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) 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

      SUBROUTINE STREAMICE_ADVECT_THICKNESS_Y ( myThid ,     
     O   hflux_y ,
     O   h , 
     I   time_step )

      IMPLICIT NONE

C     O   hflux_y ! 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_y          (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 vface, phi
      _RL stencil (-1:1)
      LOGICAL H0_valid  ! there are valid cells to calculate a 
                        ! slope-limited 2nd order flux
      _RL SLOPE_LIMITER
      external SLOPE_LIMITER

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-1,sNy+2
         Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
         DO i=1-2,sNx+2
          Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
          IF ((Gi.ge.1) .and. (Gi.le.Nx)) THEN
C         THESE ARRAY BOUNDS INSURE THAT AFTER THIS STEP,
C         VALUES WILL BE RELIABLE 1 GRID CELLS OUT IN THE
C         Y DIRECTION 
           IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .or.
     &         ((STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) .and.
     &          (STREAMICE_hmask(i,j,bi,bj).ne.1.0))) THEN

            IF (STREAMICE_vfacemask(i,j,bi,bj).eq.4.0) THEN
             hflux_y (i,j,bi,bj) = v_flux_bdry_SI (i,j,bi,bj)
            ELSE
           
             vface = .5 * 
     &        (V_streamice(i,j,bi,bj)+V_streamice(i+1,j,bi,bj))

             IF (vface .gt. 0. _d 0) THEN
              DO k=-1,1
               stencil (k) = h(i,j+k-1,bi,bj)
              ENDDO
              IF ((STREAMICE_hmask(i,j,bi,bj).eq.1.0) .and.
     &            (STREAMICE_hmask(i,j-2,bi,bj).eq.1.0)) H0_valid=.true. 
             
              IF ((Gj.eq.1).and.(STREAMICE_hmask(i,j-1,bi,bj).eq.3.0)) 
     &         THEN  ! we are at western bdry and there is a thick. bdry cond
               hflux_y (i,j,bi,bj) = h(i,j-1,bi,bj) * vface
              ELSEIF (H0_valid) THEN
               phi = SLOPE_LIMITER ( 
     &          stencil(0)-stencil(-1),
     &          stencil(1)-stencil(0))
               hflux_y (i,j,bi,bj) = vface *
     &          (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_y (i,j,bi,bj) = vface * stencil(0)
              ENDIF
             ELSE ! uface <= 0
              DO k=-1,1
               stencil (k) = h(i,j-k,bi,bj)
              ENDDO
              IF ((STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) .and.
     &            (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0)) H0_valid=.true. 

              IF ((Gj.eq.Ny).and.(STREAMICE_hmask(i,j+1,bi,bj).eq.3.0)) 
     &         THEN  ! we are at western bdry and there is a thick. bdry cond
               hflux_y (i,j,bi,bj) = h(i,j+1,bi,bj) * vface
              ELSEIF (H0_valid) THEN
               phi = SLOPE_LIMITER ( 
     &          stencil(0)-stencil(-1),
     &          stencil(1)-stencil(0))
               hflux_y (i,j,bi,bj) = vface *
     &          (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_y (i,j,bi,bj) = vface * stencil(0)
              ENDIF

             ENDIF ! uface  0

            ENDIF 
           ENDIF
          ENDIF
         ENDDO
        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-1,sNy+1
         DO i=1-2,sNx+2
          Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
          IF ((Gi .ge. 1) .and. (Gi .le. Nx)) THEN
           IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
            h(i,j,bi,bj) = h(i,j,bi,bj) - time_step * 
     &       (hflux_y(i,j+1,bi,bj)*dxG(i,j+1,bi,bj) -
     &        hflux_y(i,j,bi,bj)*dxG(i,j,bi,bj)) *
     &       recip_rA (i,j,bi,bj)
           ENDIF
          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO

!       CALL WRITE_FLD_XY_RL ("h_after_yflux","",
!      & h, 0, myThid)

#endif
      RETURN
      END SUBROUTINE STREAMICE_ADVECT_THICKNESS_Y