C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_diffusion.F,v 1.2 2010/01/03 21:09:28 jmc Exp $
C $Name:  $

#include "THSICE_OPTIONS.h"

CBOP
C !ROUTINE: THSICE_DIFFUSION

C !INTERFACE: ==========================================================
      SUBROUTINE THSICE_DIFFUSION(
     I                  maskOc,
     U                  uIce, vIce,
     I                  bi, bj, myTime, myIter, myThid )

C !DESCRIPTION: \bv
C     *===========================================================*
C     | SUBROUTINE THSICE_DIFFUSION
C     | o Account for total (ice+snow) thickness diffusion by
C     |   modifying ice-velocity:
C     |   If no velocity in the first place, and if using 1rst Order
C     |   upwind adv.scheme, this is equivalent to a diffusion of
C     |   ice+snow thichness.
C     *===========================================================*
C \ev

C !USES: ===============================================================
      IMPLICIT NONE

C     === Global variables ===

#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "THSICE_SIZE.h"
#include "THSICE_PARAMS.h"
#include "THSICE_VARS.h"

C !INPUT PARAMETERS: ===================================================
C     === Routine arguments ===
C     maskOc    :: ocean surface mask (0=land ; 1=ocean)
C     uIce/vIce :: current ice velocity on C-grid [m/s]
C     bi,bj     :: Tile indices
C     myTime    :: Current time in simulation (s)
C     myIter    :: Current iteration number
C     myThid    :: My Thread Id number
      _RS     maskOc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     uIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     vIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER bi,bj
      _RL     myTime
      INTEGER myIter
      INTEGER myThid

#ifdef ALLOW_THSICE
C !LOCAL VARIABLES: ====================================================
C     === Local variables ===
C     i,j,      :: Loop counters
C     iceFld    :: sea-ice + snow mass density
C     msgBuf    :: Informational/error message buffer
      INTEGER i, j
c     CHARACTER*(MAX_LEN_MBUF) msgBuf
      _RL     iceFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     tmpFld, hIceEpsil
      LOGICAL dBugFlag
c#include "THSICE_DEBUG.h"
CEOP

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

      dBugFlag = debugLevel.GE.debLevB
      hIceEpsil = 1. _d -10

      IF ( thSIce_diffK .GT. 0. ) THEN
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          iceFld(i,j)  = ( rhos*snowHeight(i,j,bi,bj)
     &                    +rhoi*iceHeight(i,j,bi,bj) )
c                        *iceMask(i,j,bi,bj)
         ENDDO
        ENDDO

        DO j=1-OLy,sNy+OLy
         DO i=1-OLx+1,sNx+OLx
          tmpFld = MAX( iceFld(i-1,j),iceFld(i,j) )
     &                * maskOc(i-1,j)*maskOc(i,j)
          IF ( tmpFld.GT.hIceEpsil )
     &    uIce(i,j) = uIce(i,j)
     &              + thSIce_diffK*( iceFld(i-1,j)-iceFld(i,j) )
     &                            *recip_dxC(i,j,bi,bj)/tmpFld
         ENDDO
        ENDDO

        DO j=1-OLy+1,sNy+OLy
         DO i=1-OLx,sNx+OLx
          tmpFld = MAX( iceFld(i,j-1),iceFld(i,j) )
     &                 *maskOc(i,j-1)*maskOc(i,j)
          IF ( tmpFld.GT.hIceEpsil )
     &    vIce(i,j) = vIce(i,j)
     &              + thSIce_diffK*( iceFld(i,j-1)-iceFld(i,j) )
     &                            *recip_dyC(i,j,bi,bj)/tmpFld
         ENDDO
        ENDDO

      ENDIF

#endif /* ALLOW_THSICE */

      RETURN
      END