C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_itd_remap.F,v 1.2 2016/11/29 22:04:08 jmc Exp $
C $Name:  $

C     contains:
C     S/R SEAICE_ITD_REMAP
C     S/R SEAICE_ITD_REMAP_LINEAR
C     S/R SEAICE_ITD_REMAP_CHECK_BOUNDS

#include "SEAICE_OPTIONS.h"

CBOP
C !ROUTINE: SEAICE_ITD_REMAP

C !INTERFACE: ==========================================================
      SUBROUTINE SEAICE_ITD_REMAP(
     I     heffitdpre, areaitdpre,
     I     bi, bj, myTime, myIter, myThid )

C !DESCRIPTION: \bv
C     *===========================================================*
C     | SUBROUTINE SEAICE_ITD_REMAP
C     | o checks if absolute ice thickness in any category
C     |   exceeds its category limits
C     | o remaps sea ice area and volume
C     |   and associated ice properties in thickness space
C     |   following the remapping scheme of Lipscomb (2001), JGR
C     |
C     | Martin Losch, started in May 2014, Martin.Losch@awi.de
C     | with many fixes by Mischa Ungermann (MU)
C     *===========================================================*
C \ev

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

C     === Global variables to be checked and remapped ===
C     AREAITD   :: sea ice area      by category
C     HEFFITD   :: sea ice thickness by category
C
C     === Global variables to be remappped ===
C     HSNOWITD  :: snow thickness    by category
C     enthalpy ?
C     temperature ?
C     salinity ?
C     age ?
C
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"

C !INPUT PARAMETERS: ===================================================
C     === Routine arguments ===
C     bi, bj    :: outer loop counters
C     myTime    :: current time
C     myIter    :: iteration number
C     myThid    :: Thread no. that called this routine.
      _RL myTime
      INTEGER bi,bj
      INTEGER myIter
      INTEGER myThid
      _RL heffitdPre  (1:sNx,1:sNy,1:nITD)
      _RL areaitdPre  (1:sNx,1:sNy,1:nITD)
CEndOfInterface

#ifdef SEAICE_ITD

C !LOCAL VARIABLES: ====================================================
C     === Local variables ===
C     i,j,k       :: inner loop counters
C
      INTEGER i, j, k
      INTEGER kDonor, kRecvr
      _RL slope, area_reg_sq, hice_reg_sq
      _RL etaMin, etaMax, etam, etap, eta2
      _RL dh0, da0, daMax
CMU      _RL oneMinusEps
      _RL third
      PARAMETER ( third = 0.333333333333333333333333333 _d 0 )
C
      _RL dhActual    (1:sNx,1:sNy,1:nITD)
      _RL hActual     (1:sNx,1:sNy,1:nITD)
      _RL hActualPre  (1:sNx,1:sNy,1:nITD)
      _RL dheff, darea, dhsnw
C
      _RL hLimitNew   (1:sNx,1:sNy,0:nITD)
C     coefficients for represent g(h)
C     g0 :: constant coefficient in g(h)
C     g1 :: linear  coefficient in g(h)
C     hL :: left end of range over which g(h) > 0
C     hL :: right end of range over which g(h) > 0
      _RL g0 (1:sNx,1:sNy,0:nITD)
      _RL g1 (1:sNx,1:sNy,0:nITD)
      _RL hL (1:sNx,1:sNy,0:nITD)
      _RL hR (1:sNx,1:sNy,0:nITD)
C     local copy of AREAITD
      _RL aLoc(1:sNx,1:sNy)
C
      LOGICAL doRemapping (1:sNx,1:sNy)
CEOP
C---+-|--1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C     constants
      area_reg_sq = SEAICE_area_reg**2
      hice_reg_sq = SEAICE_hice_reg**2
CMU      oneMinusEps = 1. _d 0 - SEAICE_eps
C     initialisation
      DO j=1,sNy
       DO i=1,sNx
        doRemapping(i,j) = .FALSE.
        IF ( HEFFM(i,j,bi,bj) .NE. 0. _d 0 ) doRemapping(i,j) = .TRUE.
       ENDDO
      ENDDO
C     do not compute regularized hActual as in seaice_growth, because
C     with regularization, hActual deviates too much from the actual
C     category boundaries and the boundary computation fails too often.
      DO k=1,nITD
       DO j=1,sNy
        DO i=1,sNx
         hActualPre (i,j,k) = 0. _d 0
         hActual (i,j,k) = 0. _d 0
         dhActual(i,j,k) = 0. _d 0
         IF (.FALSE.) THEN
          IF ( areaitdPre(i,j,k) .GT. 0. _d 0 ) THEN
           hActualPre(i,j,k) = heffitdPre(i,j,k)
     &         /SQRT( areaitdPre(i,j,k)**2 + area_reg_sq )
CML           hActualPre(i,j,k) = SQRT( hActualPre(i,j,k)**2 + hice_reg_sq )
          ENDIF
          IF ( AREAITD(i,j,k,bi,bj) .GT. 0. _d 0 ) THEN
           hActual(i,j,k) = HEFFITD(i,j,k,bi,bj)
     &         /SQRT( AREAITD(i,j,k,bi,bj)**2 + area_reg_sq )
CML           hActual(i,j,k) = SQRT( hActual(i,j,k)**2 + hice_reg_sq )
          ENDIF
          dhActual(i,j,k) = hActual(i,j,k) - hActualPre(i,j,k)
         ELSE
          IF ( areaitdPre(i,j,k) .GT. SEAICE_area_reg ) THEN
           hActualPre(i,j,k) = heffitdPre(i,j,k)/areaitdPre(i,j,k)
          ENDIF
          IF ( AREAITD(i,j,k,bi,bj) .GT. SEAICE_area_reg ) THEN
           hActual(i,j,k) = HEFFITD(i,j,k,bi,bj)/AREAITD(i,j,k,bi,bj)
          ENDIF
          dhActual(i,j,k) = hActual(i,j,k) - hActualPre(i,j,k)
         ENDIF
        ENDDO
       ENDDO
      ENDDO
C
C     compute new category boundaries
C
      DO j=1,sNy
       DO i=1,sNx
        hLimitNew(i,j,0) = hLimit(0)
       ENDDO
      ENDDO
      DO k=1,nITD-1
       DO j=1,sNy
        DO i=1,sNx
         IF ( hActualPre(i,j,k)  .GT.SEAICE_eps .AND.
     &        hActualPre(i,j,k+1).GT.SEAICE_eps ) THEN
          slope = ( dhActual(i,j,k+1) - dhActual(i,j,k) )
     &         /( hActualPre(i,j,k+1) - hActualPre(i,j,k) )
          hLimitNew(i,j,k) =   hLimit(k) + dhActual(i,j,k)
     &         +     slope * ( hLimit(k) - hActualPre(i,j,k) )
         ELSEIF ( hActualPre(i,j,k)  .GT.SEAICE_eps ) THEN
          hLimitNew(i,j,k) = hLimit(k) + dhActual(i,j,k)
         ELSEIF ( hActualPre(i,j,k+1).GT.SEAICE_eps ) THEN
          hLimitNew(i,j,k) = hLimit(k) + dhActual(i,j,k+1)
         ELSE
          hLimitNew(i,j,k) = hLimit(k)
         ENDIF
C     After computing the new boundary, check
C     (1) if it is between two adjacent thicknesses
         IF ( ( AREAITD(i,j,k,bi,bj).GT.SEAICE_area_reg .AND.
     &          hActual(i,j,k) .GE. hLimitNew(i,j,k) ) .OR.
     &        ( AREAITD(i,j,k+1,bi,bj).GT.SEAICE_area_reg .AND.
     &          hActual(i,j,k+1) .LE. hLimitNew(i,j,k) ) )
     &        doRemapping(i,j) = .FALSE.
C     (2) that it is been the old boudnaries k-1 and k+1
C     (Note from CICE: we could allow this, but would make the code
C     more complicated)
         IF ( ( hLimitNew(i,j,k) .GT. hLimit(k+1) ) .OR.
     &        ( hLimitNew(i,j,k) .LT. hLimit(k-1) ) )
     &        doRemapping(i,j) = .FALSE.
        ENDDO
       ENDDO
      ENDDO
C     Report problems, if there are any. Because this breaks optimization
C     do not do it by default.
C     Where doRemapping is false, the rebinning of seaice_itd_redist
C     (called at the end) will take care of shifting the ice.
      IF ( debugLevel.GE.debLevA )
     &     CALL SEAICE_ITD_REMAP_CHECK_BOUNDS(
     I     AREAITD, hActual, hActualPre, hLimitNew, doRemapping,
     I     bi, bj, myTime, myIter, myThid )
C     computing the upper limit of the thickest category does not require
C     any checks and can be computed now
      k = nITD
      DO j=1,sNy
       DO i=1,sNx
        hLimitNew(i,j,k) = hLimit(k)
        IF ( AREAITD(i,j,k,bi,bj).GT.SEAICE_area_reg )
     &       hLimitNew(i,j,k) = MAX( 3. _d 0*hActual(i,j,k)
     &       - 2. _d 0 * hLimitNew(i,j,k-1), hLimit(k-1) )
       ENDDO
      ENDDO
C
C     end of limit computation, now compute the coefficients of the
C     linear approximations of g(h) => g(eta) = g0 + g1*eta
C
C     CICE does something specical for the first category.
C     compute coefficients for 1st category
      k = 1
      DO j=1,sNy
       DO i=1,sNx
C     initialisation
        aLoc(i,j) = AREAITD(i,j,k,bi,bj)
C     initialise hL and hR
C     this single line is different from the code that follows below
C     for all categories
        hL(i,j,k) = hLimitNew(i,j,k-1)
        hR(i,j,k) = hLimit(k)
       ENDDO
      ENDDO
      CALL SEAICE_ITD_REMAP_LINEAR(
     O     g0(1,1,k), g1(1,1,k),
     U     hL(1,1,k), hR(1,1,k),
     I     hActual(1,1,k), aLoc,
     I     SEAICE_area_reg, SEAICE_eps, doRemapping,
     I     myTime, myIter, myThid )
C
C     Find area lost due to melting of thin (category 1) ice
C
      DO j=1,sNy
       DO i=1,sNx
        IF ( doRemapping(i,j) .AND.
     &       AREAITD(i,j,k,bi,bj) .GT. SEAICE_area_reg ) THEN
CMU if melting of ice in category 1
         IF ( dhActual(i,j,k) .LT. 0. _d 0 ) THEN
C     integrate g(1) from zero to abs(dhActual)
CMU dh0 is max thickness of ice in first category that is melted
          dh0    = MIN(-dhActual(i,j,k),hLimit(k))
          etaMax = MIN(dh0,hR(i,j,k)) - hL(i,j,k)
          IF ( etaMax  0. _d 0 ) THEN
CMU da0 is /int_0^dh0 g dh
           da0 = g0(i,j,k)*etaMax + g1(i,j,k)*etaMax*etaMax*0.5 _d 0
           daMax = AREAITD(i,j,k,bi,bj)
     &          * ( 1. _d 0 - hActual(i,j,k)/hActualPre(i,j,k))
           da0 = MIN( da0, daMax )
CMU adjust thickness to conserve volume
           IF ( (AREAITD(i,j,k,bi,bj)-da0) .GT. SEAICE_area_reg ) THEN
             hActual(i,j,k) = hActual(i,j,k)
     &            * AREAITD(i,j,k,bi,bj)/( AREAITD(i,j,k,bi,bj) - da0 )
           ELSE
             hActual(i,j,k) = ZERO
             da0 = AREAITD(i,j,k,bi,bj)
           ENDIF
CMU increase open water fraction
           AREAITD(i,j,k,bi,bj) = AREAITD(i,j,k,bi,bj) - da0
          ENDIF
         ELSE
CMU H_0* = F_0 * dT
          hLimitNew(i,j,k-1) = MIN( dhActual(i,j,k), hLimit(k) )
         ENDIF
        ENDIF
       ENDDO
      ENDDO
C
C     compute all coefficients
C
      DO k=1,nITD
       DO j=1,sNy
        DO i=1,sNx
C     initialisation
         aLoc(i,j) = AREAITD(i,j,k,bi,bj)
C     initialise hL and hR
         hL(i,j,k) = hLimitNew(i,j,k-1)
         hR(i,j,k) = hLimitNew(i,j,k)
        ENDDO
       ENDDO
       CALL SEAICE_ITD_REMAP_LINEAR(
     O      g0(1,1,k), g1(1,1,k),
     U      hL(1,1,k), hR(1,1,k),
     I      hActual(1,1,k), aLoc,
     I      SEAICE_area_reg, SEAICE_eps, doRemapping,
     I      myTime, myIter, myThid )
      ENDDO
C
      DO k=1,nITD-1
       DO j=1,sNy
        DO i=1,sNx
         dheff = 0. _d 0
         darea = 0. _d 0
         IF ( doRemapping(i,j) ) THEN
C     compute integration limits in eta space
          IF ( hLimitNew(i,j,k) .GT. hLimit(k) ) THEN
           etaMin = MAX(       hLimit(k), hL(i,j,k)) - hL(i,j,k)
           etaMax = MIN(hLimitNew(i,j,k), hR(i,j,k)) - hL(i,j,k)
           kDonor = k
           kRecvr = k+1
          ELSE
           etaMin = 0. _d 0
           etaMax = MIN(hLimit(k), hR(i,j,k+1)) - hL(i,j,k+1)
           kDonor = k+1
           kRecvr = k
          ENDIF
C     compute the area and volume to be moved
          IF ( etaMax .GT. etaMin ) THEN
           etam  = etaMax-etaMin
           etap  = etaMax+etaMin
           eta2  = 0.5*etam*etap
           darea = g0(i,j,kDonor)*etam + g1(i,j,kDonor)*eta2
CML           dheff = g0(i,j,kDonor)*eta2
CML     &          +  g1(i,j,kDonor)*etam*(etap*etap-etaMax*etaMin)*third
CML     &          +  darea*hL(i,j,kDonor)
           dheff = g0(i,j,kDonor)*eta2
     &          +  g1(i,j,kDonor)*(etaMax**3-etaMin**3)*third
     &          +  darea*hL(i,j,kDonor)
          ENDIF
C     ... or shift entire category, if nearly all ice is to be shifted.
CMU          IF ( (darea .GT.AREAITD(i,j,kDonor,bi,bj)*oneMinusEps).OR.
CMU     &         (dheff .GT.HEFFITD(i,j,kDonor,bi,bj)*oneMinusEps) ) THEN
          IF ( (darea .GT.AREAITD(i,j,kDonor,bi,bj)-SEAICE_eps).OR.
     &         (dheff .GT.HEFFITD(i,j,kDonor,bi,bj)-SEAICE_eps) ) THEN
           darea = AREAITD(i,j,kDonor,bi,bj)
           dheff = HEFFITD(i,j,kDonor,bi,bj)
          ENDIF
C     regularize: reset to zero, if there is too little ice to be shifted ...
CMU          IF ( (darea .LT. AREAITD(i,j,kDonor,bi,bj)*SEAICE_eps).OR.
CMU     &         (dheff .LT. HEFFITD(i,j,kDonor,bi,bj)*SEAICE_eps) ) THEN
          IF ( (darea .LT. SEAICE_eps).OR.
     &         (dheff .LT. SEAICE_eps) ) THEN
           darea  = 0. _d 0
           dheff  = 0. _d 0
          ENDIF
C     snow scaled by area
          IF ( AREAITD(i,j,kDonor,bi,bj) .GT. SEAICE_area_reg ) THEN
C     snow scaled by area (why not volume?), CICE also does it in this way
           dhsnw = darea/AREAITD(i,j,kDonor,bi,bj)
     &          * HSNOWITD(i,j,kDonor,bi,bj)
CMU          IF ( HEFFITD(i,j,kDonor,bi,bj) .GT. SEAICE_hice_reg ) THEN
CMU           dhsnw = dheff/HEFFITD(i,j,kDonor,bi,bj)
CMU     &         * HSNOWITD(i,j,kDonor,bi,bj)
          ELSE
           dhsnw = HSNOWITD(i,j,kDonor,bi,bj)
          ENDIF
C     apply increments
          HEFFITD(i,j,kRecvr,bi,bj) = HEFFITD(i,j,kRecvr,bi,bj) + dheff
          HEFFITD(i,j,kDonor,bi,bj) = HEFFITD(i,j,kDonor,bi,bj) - dheff
          AREAITD(i,j,kRecvr,bi,bj) = AREAITD(i,j,kRecvr,bi,bj) + darea
          AREAITD(i,j,kDonor,bi,bj) = AREAITD(i,j,kDonor,bi,bj) - darea
          HSNOWITD(i,j,kRecvr,bi,bj)=HSNOWITD(i,j,kRecvr,bi,bj) + dhsnw
          HSNOWITD(i,j,kDonor,bi,bj)=HSNOWITD(i,j,kDonor,bi,bj) - dhsnw
C     end if doRemapping
         ENDIF
        ENDDO
       ENDDO
      ENDDO

      RETURN
      END


C---+-|--1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: SEAICE_ITD_REMAP_LINEAR C !INTERFACE: ========================================================== SUBROUTINE SEAICE_ITD_REMAP_LINEAR( O g0, g1, U hL, hR, I hActual, area, I SEAICE_area_reg, SEAICE_eps, doRemapping, I myTime, myIter, myThid ) C !DESCRIPTION: \bv C *===========================================================* C | SUBROUTINE SEAICE_ITD_REMAP_LINEAR C | o compute coefficients g0, g1 for piece-wise linear fit C | g(h) = g0 + g1*h C | o compute range boundaries hL, hR for this linear fit C | C | Martin Losch, May 2014, Martin.Losch@awi.de C *===========================================================* C \ev C !USES: =============================================================== IMPLICIT NONE #include "SIZE.h" C !INPUT PARAMETERS: =================================================== C === Routine arguments === C myTime :: current time C myIter :: iteration number C myThid :: Thread no. that called this routine. _RL myTime INTEGER myIter INTEGER myThid C C OUTPUT: coefficients for representing g(h) C g0 :: constant coefficient in g(h) C g1 :: linear coefficient in g(h) C hL :: left end of range over which g(h) > 0 C hL :: right end of range over which g(h) > 0 _RL g0 (1:sNx,1:sNy) _RL g1 (1:sNx,1:sNy) _RL hL (1:sNx,1:sNy) _RL hR (1:sNx,1:sNy) C INPUT: C hActual :: ice thickness of current category C area :: ice concentration of current category _RL hActual (1:sNx,1:sNy) _RL area (1:sNx,1:sNy) C regularization constants _RL SEAICE_area_reg _RL SEAICE_eps C doRemapping :: mask where can be done, excludes points where C new category limits are outside certain bounds LOGICAL doRemapping (1:sNx,1:sNy) CEndOfInterface C !LOCAL VARIABLES: ==================================================== C === Local variables === C i,j :: inner loop counters C INTEGER i, j C auxCoeff :: helper variable C recip_etaR :: reciprocal of range interval in eta space C etaNoR :: ratio of distance to lower limit over etaR _RL auxCoeff _RL recip_etaR, etaNoR _RL third, sixth PARAMETER ( third = 0.333333333333333333333333333 _d 0 ) PARAMETER ( sixth = 0.666666666666666666666666666 _d 0 ) CEOP C C initialisation of hL, hR is done outside this routine C DO j=1,sNy DO i=1,sNx g0(i,j) = 0. _d 0 g1(i,j) = 0. _d 0 IF ( doRemapping(i,j) .AND. & area(i,j) .GT. SEAICE_area_reg .AND. & hR(i,j) - hL(i,j) .GT. SEAICE_eps ) THEN C change hL and hR if hActual falls outside the central third of the range IF ( hActual(i,j) .LT. (2. _d 0*hL(i,j) + hR(i,j))*third ) THEN hR(i,j) = 3. _d 0 * hActual(i,j) - 2. _d 0 * hL(i,j) ELSEIF ( hActual(i,j).GT.(hL(i,j)+2. _d 0*hR(i,j))*third ) THEN hL(i,j) = 3. _d 0 * hActual(i,j) - 2. _d 0 * hR(i,j) ENDIF C calculate new etaR = hR - hL; C catch the case of hR=hL, which can happen when hActual=hR or hL C before entering this routine; in this case g0=g1=0. recip_etaR = 0. _d 0 CMU IF ( hR(i,j) .GT. hL(i,j) ) ! crucial change; lets the model explode IF ( hR(i,j) - hL(i,j) .GT. SEAICE_eps ) & recip_etaR = 1. _d 0 / (hR(i,j) - hL(i,j)) C some abbreviations to avoid computing the same thing multiple times etaNoR = (hActual(i,j) - hL(i,j))*recip_etaR auxCoeff = 6. _d 0 * area(i,j)*recip_etaR C equations (14) of Lipscomb (2001), JGR g0(i,j) = auxCoeff*( sixth - etaNoR ) g1(i,j) = 2. _d 0 * auxCoeff*recip_etaR*( etaNoR - 0.5 _d 0 ) ELSE C not doRemapping C reset hL and hR hL(i,j) = 0. _d 0 hR(i,j) = 0. _d 0 ENDIF ENDDO ENDDO RETURN END


C---+-|--1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CEOP C !ROUTINE: SEAICE_ITD_REMAP_CHECK_BOUNDS C !INTERFACE: ========================================================== SUBROUTINE SEAICE_ITD_REMAP_CHECK_BOUNDS( I AREAITD, hActual, hActualPre, hLimitNew, doRemapping, I bi, bj, myTime, myIter, myThid ) C !DESCRIPTION: \bv C *===========================================================* C | SUBROUTINE SEAICE_ITD_REMAP_CHECK_BOUNDS C | o where doRemapping = .FALSE. print a warning C | C | Martin Losch, May 2014, Martin.Losch@awi.de C *===========================================================* C \ev C !USES: =============================================================== IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "SEAICE_SIZE.h" #include "SEAICE_PARAMS.h" C !INPUT PARAMETERS: =================================================== C === Routine arguments === C bi, bj :: outer loop counters C myTime :: current time C myIter :: iteration number C myThid :: Thread no. that called this routine. _RL myTime INTEGER bi,bj INTEGER myIter INTEGER myThid C hActual :: ice thickness of current category _RL hActual (1:sNx,1:sNy,1:nITD) _RL hActualPre(1:sNx,1:sNy,1:nITD) C hLimitNew :: new "advected" category boundaries after seaice_growth _RL hLimitNew (1:sNx,1:sNy,0:nITD) C AREAITD :: ice concentration of current category _RL AREAITD (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nITD,nSx,nSy) C doRemapping :: mask where can be done, excludes points where C new category limits are outside certain bounds LOGICAL doRemapping (1:sNx,1:sNy) CEndOfInterface C !LOCAL VARIABLES: ==================================================== C === Local variables === C i,j,k :: inner loop counters C INTEGER i, j, k CHARACTER*(MAX_LEN_MBUF) msgBuf CHARACTER*(39) tmpBuf CEOP DO j=1,sNy DO i=1,sNx IF (.NOT.doRemapping(i,j) ) THEN DO k=1,nITD-1 WRITE(tmpBuf,'(A,2I5,A,I10)') & ' at (', i, j, ') in timestep ', myIter IF ( AREAITD(i,j,k,bi,bj).GT.SEAICE_area_reg .AND. & hActual(i,j,k) .GE. hLimitNew(i,j,k) ) THEN WRITE(msgBuf,'(A,I3,A)') & 'SEAICE_ITD_REMAP: hActual(k) >= hLimitNew(k) '// & 'for category ', k, tmpBuf CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) CML PRINT *, hActual(i,j,k), CML & hLimitNew(i,j,k), hLimit(k) ENDIF IF ( AREAITD(i,j,k+1,bi,bj).GT.SEAICE_area_reg .AND. & hActual(i,j,k+1) .LE. hLimitNew(i,j,k) ) THEN WRITE(msgBuf,'(A,I3,A)') & 'SEAICE_ITD_REMAP: hActual(k+1) <= hLimitNew(k) '// & 'for category ', k, tmpBuf CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) PRINT '(8(1X,E10.4))', & AREAITD(i,j,k+1,bi,bj), hActual(i,j,k+1), & hActualPre(i,j,k+1), & AREAITD(i,j,k,bi,bj), hActual(i,j,k), & hActualPre(i,j,k), & hLimitNew(i,j,k), hLimit(k) ENDIF IF ( hLimitNew(i,j,k) .GT. hLimit(k+1) ) THEN WRITE(msgBuf,'(A,I3,A)') & 'SEAICE_ITD_REMAP: hLimitNew(k) > hLimitNew(k+1) '// & 'for category ', k, tmpBuf CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) ENDIF IF ( hLimitNew(i,j,k) .LT. hLimit(k-1) ) THEN WRITE(msgBuf,'(A,I3,A)') & 'SEAICE_ITD_REMAP: hLimitNew(k) < hLimitNew(k-1) '// & 'for category ', k, tmpBuf CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) ENDIF ENDDO ENDIF ENDDO ENDDO #endif /* SEAICE_ITD */ RETURN END