pkg/thsice/thsice_diffusion.F

a4eca6e929 Jean*0001 #include "THSICE_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: THSICE_DIFFUSION
                0005
                0006 C !INTERFACE: ==========================================================
                0007       SUBROUTINE THSICE_DIFFUSION(
                0008      I                  maskOc,
                0009      U                  uIce, vIce,
                0010      I                  bi, bj, myTime, myIter, myThid )
                0011
                0012 C !DESCRIPTION: \bv
                0013 C     *===========================================================*
                0014 C     | SUBROUTINE THSICE_DIFFUSION
                0015 C     | o Account for total (ice+snow) thickness diffusion by
                0016 C     |   modifying ice-velocity:
                0017 C     |   If no velocity in the first place, and if using 1rst Order
                0018 C     |   upwind adv.scheme, this is equivalent to a diffusion of
                0019 C     |   ice+snow thichness.
                0020 C     *===========================================================*
                0021 C \ev
                0022
                0023 C !USES: ===============================================================
                0024       IMPLICIT NONE
                0025
                0026 C     === Global variables ===
                0027
                0028 #include "SIZE.h"
                0029 #include "EEPARAMS.h"
                0030 #include "PARAMS.h"
                0031 #include "GRID.h"
                0032 #include "THSICE_SIZE.h"
                0033 #include "THSICE_PARAMS.h"
                0034 #include "THSICE_VARS.h"
                0035
                0036 C !INPUT PARAMETERS: ===================================================
                0037 C     === Routine arguments ===
                0038 C     maskOc    :: ocean surface mask (0=land ; 1=ocean)
                0039 C     uIce/vIce :: current ice velocity on C-grid [m/s]
                0040 C     bi,bj     :: Tile indices
                0041 C     myTime    :: Current time in simulation (s)
                0042 C     myIter    :: Current iteration number
                0043 C     myThid    :: My Thread Id number
                0044       _RS     maskOc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0045       _RL     uIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0046       _RL     vIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0047       INTEGER bi,bj
                0048       _RL     myTime
                0049       INTEGER myIter
                0050       INTEGER myThid
                0051
                0052 #ifdef ALLOW_THSICE
                0053 C !LOCAL VARIABLES: ====================================================
                0054 C     === Local variables ===
                0055 C     i,j,      :: Loop counters
                0056 C     iceFld    :: sea-ice + snow mass density
a444616fe4 Jean*0057 C     msgBuf    :: Informational/error message buffer
a4eca6e929 Jean*0058       INTEGER i, j
a444616fe4 Jean*0059 c     CHARACTER*(MAX_LEN_MBUF) msgBuf
a4eca6e929 Jean*0060       _RL     iceFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0061       _RL     tmpFld, hIceEpsil
                0062       LOGICAL dBugFlag
                0063 c#include "THSICE_DEBUG.h"
                0064 CEOP
                0065
                0066 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0067
ae605e558b Jean*0068       dBugFlag = debugLevel.GE.debLevC
a4eca6e929 Jean*0069       hIceEpsil = 1. _d -10
                0070
                0071       IF ( thSIce_diffK .GT. 0. ) THEN
                0072         DO j=1-OLy,sNy+OLy
                0073          DO i=1-OLx,sNx+OLx
                0074           iceFld(i,j)  = ( rhos*snowHeight(i,j,bi,bj)
                0075      &                    +rhoi*iceHeight(i,j,bi,bj) )
                0076 c                        *iceMask(i,j,bi,bj)
                0077          ENDDO
                0078         ENDDO
                0079
                0080         DO j=1-OLy,sNy+OLy
                0081          DO i=1-OLx+1,sNx+OLx
                0082           tmpFld = MAX( iceFld(i-1,j),iceFld(i,j) )
                0083      &                * maskOc(i-1,j)*maskOc(i,j)
                0084           IF ( tmpFld.GT.hIceEpsil )
                0085      &    uIce(i,j) = uIce(i,j)
                0086      &              + thSIce_diffK*( iceFld(i-1,j)-iceFld(i,j) )
                0087      &                            *recip_dxC(i,j,bi,bj)/tmpFld
                0088          ENDDO
                0089         ENDDO
                0090
                0091         DO j=1-OLy+1,sNy+OLy
                0092          DO i=1-OLx,sNx+OLx
                0093           tmpFld = MAX( iceFld(i,j-1),iceFld(i,j) )
                0094      &                 *maskOc(i,j-1)*maskOc(i,j)
                0095           IF ( tmpFld.GT.hIceEpsil )
                0096      &    vIce(i,j) = vIce(i,j)
                0097      &              + thSIce_diffK*( iceFld(i,j-1)-iceFld(i,j) )
                0098      &                            *recip_dyC(i,j,bi,bj)/tmpFld
                0099          ENDDO
                0100         ENDDO
                0101
                0102       ENDIF
                0103
                0104 #endif /* ALLOW_THSICE */
                0105
                0106       RETURN
                0107       END