C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_reg_ridge.F,v 1.5 2014/10/20 03:20:58 gforget Exp $
C $Name:  $

#include "SEAICE_OPTIONS.h"
#ifdef ALLOW_EXF
# include "EXF_OPTIONS.h"
#endif
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif

CBOP
C     !ROUTINE: SEAICE_REG_RIDGE
C     !INTERFACE:
      SUBROUTINE SEAICE_REG_RIDGE( myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *=================================================================*
C     | SUBROUTINE seaice_reg_ridge
C     | o this routine has two purposes:
C     |   (1) clean up after advection (undershoots etc.);
C     |       after advection, the sea ice variables may have unphysical
C     |       values, e.g. < 0 or very thin ice, that are regularized
C     |       here.
C     |   (2) driver for ice ridging;
C     |       concentration as a special case may be > 1 in convergent
C     |       motion and a ridging algorithm redistributes the ice to
C     |       limit the concentration to 1.
C     | o called after S/R seaice_advdiff
C     *=================================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"
#include "SEAICE_TRACER.h"
#ifdef ALLOW_EXF
# include "EXF_FIELDS.h"
#endif
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif /* ALLOW_AUTODIFF_TAMC */

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     myTime :: Simulation time
C     myIter :: Simulation timestep number
C     myThid :: Thread no. that called this routine.
      _RL myTime
      INTEGER myIter, myThid

#ifdef ALLOW_SEAICE
C     !LOCAL VARIABLES:
C     === Local variables ===
C     i,j,bi,bj :: Loop counters
      INTEGER i, j, bi, bj
C     number of surface interface layer
C     IT :: ice thickness category index (ITD and SEAICE_multDim code)
      INTEGER IT
C     reciprocal of time step
      _RL recip_deltaTtherm
C     temporary variables available for the various computations
      _RL tmpscal1, tmpscal2
#ifdef SEAICE_ITD
      _RL tmpscal1itd(1:sNx,1:sNy), tmpscal2itd(1:sNx,1:sNy)
      _RL tmpscal3itd(1:sNx,1:sNy)
C     reciprocal number of ice classes nITD
      _RL recip_nitd
#endif /* SEAICE_ITD */
#ifdef ALLOW_DIAGNOSTICS
C     Helper variables for diagnostics
      _RL DIAGarrayA    (1:sNx,1:sNy)
#endif /* ALLOW_DIAGNOSTICS */
CEOP

C
C     === Routine body ===
C
      recip_deltaTtherm = ONE / SEAICE_deltaTtherm

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)

#ifdef ALLOW_AUTODIFF_TAMC
        act1 = bi - myBxLo(myThid)
        max1 = myBxHi(myThid) - myBxLo(myThid) + 1
        act2 = bj - myByLo(myThid)
        max2 = myByHi(myThid) - myByLo(myThid) + 1
        act3 = myThid - 1
        max3 = nTx*nTy
        act4 = ikey_dynamics - 1
        iicekey = (act1 + 1) + act2*max1
     &                       + act3*max1*max2
     &                       + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */

        DO J=1-OLy,sNy+OLy
         DO I=1-OLx,sNx+OLx
          d_HEFFbyNEG(I,J,bi,bj)    = 0.0 _d 0
          d_HSNWbyNEG(I,J,bi,bj)    = 0.0 _d 0
#ifdef EXF_SEAICE_FRACTION
          d_AREAbyRLX(I,J,bi,bj)    = 0.0 _d 0
          d_HEFFbyRLX(I,J,bi,bj)    = 0.0 _d 0
#endif /* EXF_SEAICE_FRACTION */
#ifdef SEAICE_VARIABLE_SALINITY
          saltFluxAdjust(I,J,bi,bj) = 0.0 _d 0
#endif /* SEAICE_VARIABLE_SALINITY */
         ENDDO
        ENDDO

C =====================================================================
C ========== PART 1: treat pathological cases (post advdiff) ==========
C =====================================================================

#if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
Cgf no dependency through pathological cases treatment
        IF ( SEAICEadjMODE.EQ.0 ) THEN
#endif

#ifdef EXF_SEAICE_FRACTION
CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
C--   (0) relax sea ice concentration towards observation
        IF ( SEAICE_tauAreaObsRelax .GT. zeroRL ) THEN
         DO J=1,sNy
          DO I=1,sNx
           IF ( exf_iceFraction(I,J,bi,bj).GT.AREA(I,J,bi,bj) ) THEN
            d_AREAbyRLX(i,j,bi,bj) =
     &       SEAICE_deltaTtherm/SEAICE_tauAreaObsRelax
     &       * (exf_iceFraction(I,J,bi,bj) - AREA(I,J,bi,bj))
           ENDIF
           IF ( exf_iceFraction(I,J,bi,bj).GT.zeroRS .AND.
     &          AREA(I,J,bi,bj).EQ.0. _d 0) THEN
C           d_HEFFbyRLX(i,j,bi,bj) = 1. _d 1 * siEps * d_AREAbyRLX(i,j,bi,bj)
            d_HEFFbyRLX(i,j,bi,bj) = 1. _d 1 * siEps
           ENDIF
#ifdef SEAICE_ITD
           AREAITD(I,J,1,bi,bj) = AREAITD(I,J,1,bi,bj)
     &                          +  d_AREAbyRLX(i,j,bi,bj)
           HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj)
     &                          +  d_HEFFbyRLX(i,j,bi,bj)
#endif /* SEAICE_ITD */
           AREA(I,J,bi,bj) = AREA(I,J,bi,bj) +  d_AREAbyRLX(i,j,bi,bj)
           HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) +  d_HEFFbyRLX(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDIF
#endif /* EXF_SEAICE_FRACTION */

C--   (1) treat the case of negative values:

#ifdef SEAICE_ITD
        DO IT=1,SEAICE_multDim
         DO J=1,sNy
          DO I=1,sNx
           tmpscal1=0. _d 0
           tmpscal2=0. _d 0
           tmpscal1=MAX(-HEFFITD(I,J,IT,bi,bj),0. _d 0)
           HEFFITD(I,J,IT,bi,bj)=HEFFITD(I,J,IT,bi,bj)+tmpscal1
           d_HEFFbyNEG(I,J,bi,bj)=d_HEFFbyNEG(I,J,bi,bj)+tmpscal1
           tmpscal2=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)
           HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal2
           d_HSNWbyNEG(I,J,bi,bj)=d_HSNWbyNEG(I,J,bi,bj)+tmpscal2
           AREAITD(I,J,IT,bi,bj)=MAX(AREAITD(I,J,IT,bi,bj),0. _d 0)
C     AREA, HEFF, and HSNOW will be updated at end of PART 1
C     by calling SEAICE_ITD_SUM
          ENDDO
         ENDDO
        ENDDO
C     update mean thicknesses HEFF and HSNOW and total ice
C     concentration AREA to match single category values
        CALL SEAICE_ITD_SUM   ( bi, bj, myTime, myIter, myThid )
#else /* ndef SEAICE_ITD */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
        DO J=1,sNy
         DO I=1,sNx
          d_HEFFbyNEG(I,J,bi,bj)=MAX(-HEFF(I,J,bi,bj),0. _d 0)
          HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J,bi,bj)
          d_HSNWbyNEG(I,J,bi,bj)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)
          HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J,bi,bj)
          AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)
         ENDDO
        ENDDO
#endif /* SEAICE_ITD */

C--   (2) treat the case of very thin ice:

#ifdef SEAICE_ITD
C     Here we risk that even though HEFF may be larger than siEps (=1e-5)
C     HEFFITD can have classes with very small (< siEps) non-zero ice volume.
C     We avoid applying the correction to each class because that leads to
C     funny structures in the net heat and freshwater flux into the ocean.
C     Let us keep our fingers crossed, that the model will be benign!
        DO IT=1,SEAICE_multDim
         DO J=1,sNy
          DO I=1,sNx
           IF (HEFF(I,J,bi,bj).LE.siEps) THEN
            HEFFITD(I,J,IT,bi,bj) = 0. _d 0
            HSNOWITD(I,J,IT,bi,bj) = 0. _d 0
           ENDIF
          ENDDO
         ENDDO
        ENDDO
#endif /* SEAICE_ITD */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
        DO J=1,sNy
         DO I=1,sNx
          tmpscal1=0. _d 0
          tmpscal2=0. _d 0
          IF (HEFF(I,J,bi,bj).LE.siEps) THEN
           tmpscal1=-HEFF(I,J,bi,bj)
           tmpscal2=-HSNOW(I,J,bi,bj)
           DO IT=1,SEAICE_multDim
            TICES(I,J,IT,bi,bj)=celsius2K
           ENDDO
          ENDIF
          HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+tmpscal1
          HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal2
          d_HEFFbyNEG(I,J,bi,bj)=d_HEFFbyNEG(I,J,bi,bj)+tmpscal1
          d_HSNWbyNEG(I,J,bi,bj)=d_HSNWbyNEG(I,J,bi,bj)+tmpscal2
         ENDDO
        ENDDO

C--   (3) treat the case of area but no ice/snow:

#ifdef SEAICE_ITD
        DO IT=1,SEAICE_multDim
         DO J=1,sNy
          DO I=1,sNx
           IF ( (HEFFITD(I,J,IT,bi,bj) .EQ.0. _d 0).AND.
     &          (HSNOWITD(I,J,IT,bi,bj).EQ.0. _d 0))
     &          AREAITD(I,J,IT,bi,bj)=0. _d 0
          ENDDO
         ENDDO
        ENDDO
#else /* ndef SEAICE_ITD */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
         DO J=1,sNy
          DO I=1,sNx
           IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
     &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
         ENDDO
        ENDDO
#endif /* SEAICE_ITD */

C--   (4) treat the case of very small area:

#ifndef DISABLE_AREA_FLOOR
#ifdef SEAICE_ITD
        recip_nitd = 1. _d 0 / float(SEAICE_multDim)
        DO IT=1,SEAICE_multDim
         DO J=1,sNy
          DO I=1,sNx
           IF ((HEFFITD(I,J,IT,bi,bj).GT.0).OR.
     &          (HSNOWITD(I,J,IT,bi,bj).GT.0)) THEN
C     SEAICE_area_floor*SEAICE_multDim cannot be allowed to exceed 1
C     hence use SEAICE_area_floor devided by SEAICE_multDim
C     (or install a warning in e.g. seaice_readparms.F)
            AREAITD(I,J,IT,bi,bj)=
     &        MAX(AREAITD(I,J,IT,bi,bj),SEAICE_area_floor*recip_nitd)
          ENDIF
         ENDDO
        ENDDO
       ENDDO
#else /* ndef SEAICE_ITD */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
       DO J=1,sNy
        DO I=1,sNx
         IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN
          AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),SEAICE_area_floor)
         ENDIF
        ENDDO
       ENDDO
#endif /* SEAICE_ITD */
#endif /* DISABLE_AREA_FLOOR */

C     (5) treat sea ice salinity pathological cases
#ifdef SEAICE_VARIABLE_SALINITY
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hsalt(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
        DO J=1,sNy
         DO I=1,sNx
          IF ( (HSALT(I,J,bi,bj) .LT. 0.0).OR.
     &         (HEFF(I,J,bi,bj) .EQ. 0.0)  ) THEN
             saltFluxAdjust(I,J,bi,bj) = - HEFFM(I,J,bi,bj) *
     &            HSALT(I,J,bi,bj) * recip_deltaTtherm
             HSALT(I,J,bi,bj) = 0.0 _d 0
          ENDIF
         ENDDO
        ENDDO
#endif /* SEAICE_VARIABLE_SALINITY */

C =====================================================================
C ========== PART 2: ridging algorithm  ===============================
C =====================================================================

C     treat case of excessive ice cover, e.g., due to ridging:

#ifdef SEAICE_ITD

C     catch up with item (2) that involves category sums AREA and HEFF
        DO J=1,sNy
         DO I=1,sNx
          tmpscal1itd(i,j) = 0. _d 0
          tmpscal2itd(i,j) = 0. _d 0
          tmpscal3itd(i,j) = 0. _d 0
         ENDDO
        ENDDO
        DO IT=1,SEAICE_multDim
         DO J=1,sNy
          DO I=1,sNx
C     TICES was changed above (item 2), now update TICE as ice volume
C     weighted average of TICES
           tmpscal1itd(i,j)=tmpscal1itd(i,j)
     &          +        TICES(I,J,IT,bi,bj) * HEFFITD(I,J,IT,bi,bj)
           tmpscal2itd(i,j)=tmpscal2itd(i,j) + HEFFITD(I,J,IT,bi,bj)
C     also compute total of AREAITD for diagnostics and SItrArea
           tmpscal3itd(i,j)=tmpscal3itd(i,j) + AREAITD(I,J,IT,bi,bj)
          ENDDO
         ENDDO
        ENDDO
        DO J=1,sNy
         DO I=1,sNx
C     save pre-ridging ice concentration for diagnostics:
C     these lines are executed before "ridging" is applied to AREA
C     hence we execute them here before SEAICE_ITD_REDIST is called
C     although this means that AREA has not been completely regularized
#ifdef ALLOW_DIAGNOSTICS
          DIAGarrayA(I,J) = tmpscal3itd(i,j)
#endif
#ifdef ALLOW_SITRACER
          SItrAREA(I,J,bi,bj,1)=tmpscal3itd(i,j)
#endif
         ENDDO
        ENDDO
C     ridge ice according to Lipscomb et al. (2007), Bitz et al. (2001)
C     Thorndyke et al. (1975), Hibler (1980)
        CALL SEAICE_DO_RIDGING( bi, bj, myTime, myIter, myThid )
C     check that all ice thickness categories meet their limits
C     (includes Hibler-type ridging)
        CALL SEAICE_ITD_REDIST( bi, bj, myTime, myIter, myThid )
C     update mean thicknesses HEFF and HSNOW and total ice
C     concentration AREA to match single category values
        CALL SEAICE_ITD_SUM   ( bi, bj, myTime, myIter, myThid )

#else /* ifndef SEAICE_ITD */

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
        DO J=1,sNy
         DO I=1,sNx
C     save pre-ridging ice concentration for diagnostics
#ifdef ALLOW_DIAGNOSTICS
          DIAGarrayA(I,J) = AREA(I,J,bi,bj)
#endif /*  ALLOW_DIAGNOSTICS */
#ifdef ALLOW_SITRACER
          SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)
#endif /*  ALLOW_SITRACER */
C     this is the simple Hibler (1979)-type ridging (capping of
C     concentrations > 1) for the non-ITD sea ice model
          AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),SEAICE_area_max)
         ENDDO
        ENDDO

#endif /* SEAICE_ITD */

#if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
C        end SEAICEadjMODE.EQ.0 statement:
        ENDIF
#endif

#ifdef ALLOW_DIAGNOSTICS
        IF ( useDiagnostics ) THEN
         CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIareaPR',0,1,3,bi,bj,myThid)
        ENDIF
#endif /* ALLOW_DIAGNOSTICS */

C close bi,bj loops
       ENDDO
      ENDDO

#endif /* ALLOW_SEAICE */
      RETURN
      END