C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_advect_thickness.F,v 1.11 2015/04/20 14:26:38 dgoldberg Exp $
C $Name:  $

#include "STREAMICE_OPTIONS.h"
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif

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

CBOP
      SUBROUTINE STREAMICE_ADVECT_THICKNESS ( myThid,myIter,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"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif
#ifdef ALLOW_SHELFICE
# include "SHELFICE.h"
#endif

      INTEGER myThid, myIter
      _RL time_step

#ifdef ALLOW_STREAMICE

      INTEGER i, j, bi, bj, Gi, Gj
      _RL thick_bd, uflux, vflux, max_icfl, loc_icfl
      _RL time_step_full, time_step_rem
      _RL sec_per_year, time_step_loc, MR, SMB, TMB, irho
      _RL BCVALX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL BCVALY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS BCMASKX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS BCMASKY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL utrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vtrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL h_after_uflux_SI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL h_after_vflux_SI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL hflux_x_SI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL hflux_y_SI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)

      CHARACTER*(MAX_LEN_MBUF) msgBuf

      CALL TIMER_START ('STREAMICE_ADVECT_THICKNESS',myThid)

      sec_per_year = 365.*86400.

      time_step_loc = time_step / sec_per_year
      time_step_full = time_step_loc
      time_step_rem = time_step_loc
      PRINT *, "time_step_loc ", time_step_loc

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics
#endif

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1,sNy+1
         DO i=1,sNx+1

          H_streamice_prev(i,j,bi,bj) =
     &     H_streamice(i,j,bi,bj)

          hflux_x_SI (i,j,bi,bj) = 0. _d 0
          hflux_y_SI (i,j,bi,bj) = 0. _d 0
          h_after_uflux_SI(i,j,bi,bj) = H_streamice(i,j,bi,bj)
          h_after_vflux_SI(i,j,bi,bj) = H_streamice(i,j,bi,bj)

          IF (STREAMICE_ufacemask(i,j,bi,bj).eq.3.0) THEN
           BCMASKX(i,j,bi,bj) = 3.0
           BCVALX(i,j,bi,bj) = h_ubdry_values_SI(i,j,bi,bj)
           utrans(i,j,bi,bj) = .5 * (
     &      u_streamice(i,j,bi,bj)+u_streamice(i,j+1,bi,bj))
          ELSEIF (STREAMICE_ufacemask(i,j,bi,bj).eq.4.0) THEN
           IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
            uflux = u_flux_bdry_SI(i,j,bi,bj)
            BCMASKX(i,j,bi,bj) = 3.0
            BCVALX(i,j,bi,bj) = uflux
            utrans(i,j,bi,bj) = 1.0
           ELSEIF (STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN
            uflux = u_flux_bdry_SI(i,j,bi,bj)
            BCMASKX(i,j,bi,bj) = 3.0
            BCVALX(i,j,bi,bj) = uflux
            utrans(i,j,bi,bj) = -1.0
           ENDIF
          ELSEIF (.not.(
     &     STREAMICE_hmask(i,j,bi,bj).eq.1.0.OR.
     &     STREAMICE_hmask(i-1,j,bi,bj).eq.1.0)) THEN
           BCMASKX(i,j,bi,bj) = 0.0
           BCVALX(i,j,bi,bj) = 0. _d 0
           utrans(i,j,bi,bj) = 0. _d 0
          ELSE
           BCMASKX(i,j,bi,bj) = 0.0
           BCVALX(i,j,bi,bj) = 0. _d 0
           utrans(i,j,bi,bj) = .5 * (
     &      u_streamice(i,j,bi,bj)+u_streamice(i,j+1,bi,bj))
          ENDIF

          IF (STREAMICE_vfacemask(i,j,bi,bj).eq.3.0) THEN
           BCMASKy(i,j,bi,bj) = 3.0
           BCVALy(i,j,bi,bj) = h_vbdry_values_SI(i,j,bi,bj)
           vtrans(i,j,bi,bj) = .5 * (
     &      v_streamice(i,j,bi,bj)+v_streamice(i+1,j,bi,bj))
          ELSEIF (STREAMICE_vfacemask(i,j,bi,bj).eq.4.0) THEN
           IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
            vflux = v_flux_bdry_SI(i,j,bi,bj)
            BCMASKY(i,j,bi,bj) = 3.0
            BCVALY(i,j,bi,bj) = vflux
            vtrans(i,j,bi,bj) = 1.0
           ELSEIF (STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN
            vflux = v_flux_bdry_SI(i,j,bi,bj)
            BCMASKY(i,j,bi,bj) = 3.0
            BCVALY(i,j,bi,bj) = vflux
            vtrans(i,j,bi,bj) = -1.0
           ENDIF

           vtrans(i,j,bi,bj) = 1.0
          ELSEIF (.not.(
     &     STREAMICE_hmask(i,j,bi,bj).eq.1.0.OR.
     &     STREAMICE_hmask(i,j-1,bi,bj).eq.1.0)) THEN
           BCMASKY(i,j,bi,bj) = 0.0
           BCVALY(i,j,bi,bj) = 0. _d 0
           vtrans(i,j,bi,bj) = 0. _d 0
          ELSE
           BCMASKy(i,j,bi,bj) = 0.0
           BCVALy(i,j,bi,bj) = 0. _d 0
           vtrans(i,j,bi,bj) =  .5 * (
     &      v_streamice(i,j,bi,bj)+v_streamice(i+1,j,bi,bj))
          ENDIF

         ENDDO
        ENDDO
       ENDDO
      ENDDO

      _EXCH_XY_RL(utrans,myThid)
      _EXCH_XY_RL(vtrans,myThid)
      _EXCH_XY_RS(BCMASKx,myThid)
      _EXCH_XY_RS(BCMASKy,myThid)
      _EXCH_XY_RL(BCVALX,myThid)
      _EXCH_XY_RL(BCVALY,myThid)
      _EXCH_XY_RL(h_after_uflux_SI,myThid)
      _EXCH_XY_RL(h_after_vflux_SI,myThid)

#ifndef ALLOW_AUTODIFF

      max_icfl = 1.e-20

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1,sNy
         DO i=1,sNx
          IF (streamice_hmask(i,j,bi,bj).eq.1.0) THEN
           loc_icfl=max(abs(utrans(i,j,bi,bj)),
     &                  abs(utrans(i+1,j,bi,bj))) / dxF(i,j,bi,bj)
           loc_icfl=max(loc_icfl,max(abs(vtrans(i,j,bi,bj)),
     &                  abs(vtrans(i,j+1,bi,bj))) / dyF(i,j,bi,bj))
           if (loc_icfl.gt.max_icfl) then
            max_icfl = loc_icfl
           ENDIF
          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      CALL GLOBAL_MAX_R8 (max_icfl, myThid)

#endif /* ALLOW_AUTODIFF */

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics
CADJ STORE H_streamice  = comlev1, key=ikey_dynamics
#endif

#ifndef ALLOW_AUTODIFF
      do while (time_step_rem .gt. 1.e-15)
       time_step_loc = min (
     &  streamice_cfl_factor / max_icfl,
     &  time_step_rem )
       if (time_step_loc .lt. time_step_full) then
        PRINT *, "TAKING PARTIAL TIME STEP", time_step_loc
       endif
#endif /* ALLOW_AUTODIFF */

      CALL STREAMICE_ADV_FLUX_FL_X ( myThid ,
     I   utrans ,
     I   H_streamice ,
     I   BCMASKX,
     I   BCVALX,
     O   hflux_x_SI,
     I   time_step_loc )

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-3,sNy+3
         DO i=1,sNx
          Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
          Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
          IF (((Gj .ge. 1) .and. (Gj .le. Ny))
     &       .or.STREAMICE_NS_PERIODIC) THEN

          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) -
     &      (hflux_x_SI(i+1,j,bi,bj)*dyG(i+1,j,bi,bj) -
     &        hflux_x_SI(i,j,bi,bj)*dyG(i,j,bi,bj))
     &        * recip_rA (i,j,bi,bj) * time_step_loc
           IF ( h_after_uflux_SI (i,j,bi,bj).le.0.0) THEN
            PRINT *, "h neg after x", i,j,hflux_x_SI(i+1,j,bi,bj),
     &        hflux_x_SI(i,j,bi,bj)
           ENDIF
          ENDIF

          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE streamice_hmask  = comlev1, key=ikey_dynamics
#endif

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

      CALL STREAMICE_ADV_FLUX_FL_Y ( myThid ,
     I   vtrans ,
     I   h_after_uflux_si ,
     I   BCMASKY,
     I   BCVALY,
     O   hflux_y_SI,
     I   time_step_loc )

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1,sNy
         DO i=1,sNx
          Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
          Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
          IF (((Gj .ge. 1) .and. (Gj .le. Ny))
     &       .or.STREAMICE_EW_PERIODIC) THEN

          IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
           h_after_vflux_SI (i,j,bi,bj) = h_after_uflux_SI(i,j,bi,bj) -
     &      (hflux_y_SI(i,j+1,bi,bj)*dxG(i,j+1,bi,bj) -
     &        hflux_y_SI(i,j,bi,bj)*dxG(i,j,bi,bj)) *
     &        recip_rA (i,j,bi,bj) * time_step_loc
           IF ( h_after_vflux_SI (i,j,bi,bj).le.0.0) THEN
            PRINT *, "h neg after y", i,j,hflux_y_SI(i,j+1,bi,bj),
     &        hflux_y_SI(i,j,bi,bj)
           ENDIF

          ENDIF
         ENDIF

         ENDDO
        ENDDO
       ENDDO
      ENDDO

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1,sNy
         DO i=1,sNx
          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

! NOTE: AT THIS POINT H IS NOT VALID ON OVERLAP!!!

      if (streamice_move_front) then
      CALL STREAMICE_ADV_FRONT (
     &  myThid, time_step_loc,
     &  hflux_x_SI, hflux_y_SI )
      endif

#ifdef ALLOW_STREAMICE_2DTRACER
      CALL STREAMICE_ADVECT_2DTRACER(
     &  myThid,
     &  myIter,
     &  time_step,
     &  uTrans,
     &  vTrans,
     &  BCMASKx,
     &  BCMASKy )
#endif

#ifndef ALLOW_AUTODIFF
      time_step_rem = time_step_rem - time_step_loc
      enddo
#endif /* ALLOW_AUTODIFF */

! NOW WE APPLY MELT RATES !!
! THIS MAY BE MOVED TO A SEPARATE SUBROUTINE
      irho = 1./streamice_density

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1,sNy
         DO i=1,sNx
          Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
          Gj = (myYGlobalLo-1)+(bj-1)*sNy+j
          IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0 .or.
     &       STREAMICE_hmask(i,j,bi,bj).eq.2.0) THEN

           IF (STREAMICE_allow_cpl) THEN
#ifdef ALLOW_SHELFICE
            MR = -1. * (1.-float_frac_streamice(i,j,bi,bj)) *
     &        shelfIceFreshWaterFlux(I,J,bi,bj) * irho *
     &        sec_per_year

#else
            STOP 'SHELFICE IS NOT ENABLED'
#endif
           ELSE
              MR = (1.-float_frac_streamice(i,j,bi,bj)) *
     &             (BDOT_STREAMICE(i,j,bi,bj) +
     &              BDOT_pert(i,j,bi,bj))
           ENDIF

           SMB = ADOT_STREAMICE(i,j,bi,bj)
           TMB = SMB - MR
           IF ((TMB.lt.0.0) .and.
     &         (MR * time_step_loc .gt.
     &          H_streamice (i,j,bi,bj))) THEN
                H_streamice (i,j,bi,bj) = 0. _d 0
                STREAMICE_hmask(i,j,bi,bj) = 0.0
               PRINT *, "GOT HERE melted away! ", i,j
           ELSE
               H_streamice (i,j,bi,bj) =
     &          H_streamice (i,j,bi,bj) + TMB * time_step_loc
           ENDIF

          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      _EXCH_XY_RL(H_streamice,myThid)

      WRITE(msgBuf,'(A)') 'END STREAMICE_ADVECT_THICKNESS'
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                     SQUEEZE_RIGHT , 1)
      CALL TIMER_STOP ('STREAMICE_ADVECT_THICKNESS',myThid)

#endif
      RETURN
      END