** Warning **

Issuing rollback() due to DESTROY without explicit disconnect() of DBD::mysql::db handle dbname=MITgcm at /usr/local/share/lxr/lib/LXR/Common.pm line 1224.

Last-Modified: Fri, 10 Sep 2026 05:09:17 GMT Content-Type: text/html; charset=utf-8 MITgcm/MITgcm/pkg/seaice/seaice_itd_remap.F
Back to home page

MITgcm
 
 

    

File indexing completed on 2026-05-05 05:09:03 UTC

view on githubraw file Latest commit 3f0f10fc on 2026-05-04 14:55:37 UTC
ed2f6fecc4 Mart*0001 C     contains:
                0002 C     S/R SEAICE_ITD_REMAP
                0003 C     S/R SEAICE_ITD_REMAP_LINEAR
                0004 C     S/R SEAICE_ITD_REMAP_CHECK_BOUNDS
                0005 
                0006 #include "SEAICE_OPTIONS.h"
                0007 
                0008 CBOP
                0009 C !ROUTINE: SEAICE_ITD_REMAP
                0010 
                0011 C !INTERFACE: ==========================================================
                0012       SUBROUTINE SEAICE_ITD_REMAP(
                0013      I     heffitdpre, areaitdpre,
                0014      I     bi, bj, myTime, myIter, myThid )
                0015 
                0016 C !DESCRIPTION: \bv
                0017 C     *===========================================================*
                0018 C     | SUBROUTINE SEAICE_ITD_REMAP
                0019 C     | o checks if absolute ice thickness in any category
                0020 C     |   exceeds its category limits
                0021 C     | o remaps sea ice area and volume
                0022 C     |   and associated ice properties in thickness space
                0023 C     |   following the remapping scheme of Lipscomb (2001), JGR
                0024 C     |
                0025 C     | Martin Losch, started in May 2014, Martin.Losch@awi.de
                0026 C     | with many fixes by Mischa Ungermann (MU)
                0027 C     *===========================================================*
                0028 C \ev
                0029 
                0030 C !USES: ===============================================================
                0031       IMPLICIT NONE
                0032 
                0033 C     === Global variables to be checked and remapped ===
                0034 C     AREAITD   :: sea ice area      by category
                0035 C     HEFFITD   :: sea ice thickness by category
                0036 C
                0037 C     === Global variables to be remappped ===
                0038 C     HSNOWITD  :: snow thickness    by category
                0039 C     enthalpy ?
                0040 C     temperature ?
                0041 C     salinity ?
                0042 C     age ?
                0043 C
                0044 #include "SIZE.h"
                0045 #include "EEPARAMS.h"
                0046 #include "PARAMS.h"
                0047 #include "SEAICE_SIZE.h"
                0048 #include "SEAICE_PARAMS.h"
3f0f10fc37 Mart*0049 #include "SEAICE_GRID.h"
ed2f6fecc4 Mart*0050 #include "SEAICE.h"
                0051 
                0052 C !INPUT PARAMETERS: ===================================================
                0053 C     === Routine arguments ===
                0054 C     bi, bj    :: outer loop counters
                0055 C     myTime    :: current time
                0056 C     myIter    :: iteration number
                0057 C     myThid    :: Thread no. that called this routine.
                0058       _RL myTime
                0059       INTEGER bi,bj
                0060       INTEGER myIter
                0061       INTEGER myThid
                0062       _RL heffitdPre  (1:sNx,1:sNy,1:nITD)
                0063       _RL areaitdPre  (1:sNx,1:sNy,1:nITD)
                0064 
                0065 #ifdef SEAICE_ITD
                0066 
                0067 C !LOCAL VARIABLES: ====================================================
                0068 C     === Local variables ===
                0069 C     i,j,k       :: inner loop counters
                0070 C
                0071       INTEGER i, j, k
                0072       INTEGER kDonor, kRecvr
                0073       _RL slope, area_reg_sq, hice_reg_sq
                0074       _RL etaMin, etaMax, etam, etap, eta2
                0075       _RL dh0, da0, daMax
                0076 CMU      _RL oneMinusEps
                0077       _RL third
                0078       PARAMETER ( third = 0.333333333333333333333333333 _d 0 )
dcd6ed0c75 Jean*0079 C
ed2f6fecc4 Mart*0080       _RL dhActual    (1:sNx,1:sNy,1:nITD)
                0081       _RL hActual     (1:sNx,1:sNy,1:nITD)
                0082       _RL hActualPre  (1:sNx,1:sNy,1:nITD)
                0083       _RL dheff, darea, dhsnw
dcd6ed0c75 Jean*0084 C
ed2f6fecc4 Mart*0085       _RL hLimitNew   (1:sNx,1:sNy,0:nITD)
                0086 C     coefficients for represent g(h)
                0087 C     g0 :: constant coefficient in g(h)
                0088 C     g1 :: linear  coefficient in g(h)
                0089 C     hL :: left end of range over which g(h) > 0
                0090 C     hL :: right end of range over which g(h) > 0
                0091       _RL g0 (1:sNx,1:sNy,0:nITD)
                0092       _RL g1 (1:sNx,1:sNy,0:nITD)
                0093       _RL hL (1:sNx,1:sNy,0:nITD)
                0094       _RL hR (1:sNx,1:sNy,0:nITD)
                0095 C     local copy of AREAITD
                0096       _RL aLoc(1:sNx,1:sNy)
                0097       LOGICAL doRemapping (1:sNx,1:sNy)
                0098 CEOP
                0099 C---+-|--1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0100 
                0101 C     constants
                0102       area_reg_sq = SEAICE_area_reg**2
                0103       hice_reg_sq = SEAICE_hice_reg**2
                0104 CMU      oneMinusEps = 1. _d 0 - SEAICE_eps
                0105 C     initialisation
                0106       DO j=1,sNy
                0107        DO i=1,sNx
                0108         doRemapping(i,j) = .FALSE.
                0109         IF ( HEFFM(i,j,bi,bj) .NE. 0. _d 0 ) doRemapping(i,j) = .TRUE.
                0110        ENDDO
                0111       ENDDO
dcd6ed0c75 Jean*0112 C     do not compute regularized hActual as in seaice_growth, because
ed2f6fecc4 Mart*0113 C     with regularization, hActual deviates too much from the actual
                0114 C     category boundaries and the boundary computation fails too often.
                0115       DO k=1,nITD
                0116        DO j=1,sNy
                0117         DO i=1,sNx
                0118          hActualPre (i,j,k) = 0. _d 0
                0119          hActual (i,j,k) = 0. _d 0
                0120          dhActual(i,j,k) = 0. _d 0
                0121          IF (.FALSE.) THEN
                0122           IF ( areaitdPre(i,j,k) .GT. 0. _d 0 ) THEN
                0123            hActualPre(i,j,k) = heffitdPre(i,j,k)
                0124      &         /SQRT( areaitdPre(i,j,k)**2 + area_reg_sq )
dcd6ed0c75 Jean*0125 CML           hActualPre(i,j,k) = SQRT( hActualPre(i,j,k)**2 + hice_reg_sq )
ed2f6fecc4 Mart*0126           ENDIF
                0127           IF ( AREAITD(i,j,k,bi,bj) .GT. 0. _d 0 ) THEN
                0128            hActual(i,j,k) = HEFFITD(i,j,k,bi,bj)
                0129      &         /SQRT( AREAITD(i,j,k,bi,bj)**2 + area_reg_sq )
dcd6ed0c75 Jean*0130 CML           hActual(i,j,k) = SQRT( hActual(i,j,k)**2 + hice_reg_sq )
ed2f6fecc4 Mart*0131           ENDIF
                0132           dhActual(i,j,k) = hActual(i,j,k) - hActualPre(i,j,k)
                0133          ELSE
                0134           IF ( areaitdPre(i,j,k) .GT. SEAICE_area_reg ) THEN
                0135            hActualPre(i,j,k) = heffitdPre(i,j,k)/areaitdPre(i,j,k)
                0136           ENDIF
                0137           IF ( AREAITD(i,j,k,bi,bj) .GT. SEAICE_area_reg ) THEN
                0138            hActual(i,j,k) = HEFFITD(i,j,k,bi,bj)/AREAITD(i,j,k,bi,bj)
                0139           ENDIF
                0140           dhActual(i,j,k) = hActual(i,j,k) - hActualPre(i,j,k)
                0141          ENDIF
                0142         ENDDO
                0143        ENDDO
                0144       ENDDO
dcd6ed0c75 Jean*0145 C
ed2f6fecc4 Mart*0146 C     compute new category boundaries
dcd6ed0c75 Jean*0147 C
ed2f6fecc4 Mart*0148       DO j=1,sNy
                0149        DO i=1,sNx
                0150         hLimitNew(i,j,0) = hLimit(0)
                0151        ENDDO
                0152       ENDDO
                0153       DO k=1,nITD-1
                0154        DO j=1,sNy
                0155         DO i=1,sNx
                0156          IF ( hActualPre(i,j,k)  .GT.SEAICE_eps .AND.
                0157      &        hActualPre(i,j,k+1).GT.SEAICE_eps ) THEN
dcd6ed0c75 Jean*0158           slope = ( dhActual(i,j,k+1) - dhActual(i,j,k) )
ed2f6fecc4 Mart*0159      &         /( hActualPre(i,j,k+1) - hActualPre(i,j,k) )
                0160           hLimitNew(i,j,k) =   hLimit(k) + dhActual(i,j,k)
                0161      &         +     slope * ( hLimit(k) - hActualPre(i,j,k) )
                0162          ELSEIF ( hActualPre(i,j,k)  .GT.SEAICE_eps ) THEN
                0163           hLimitNew(i,j,k) = hLimit(k) + dhActual(i,j,k)
                0164          ELSEIF ( hActualPre(i,j,k+1).GT.SEAICE_eps ) THEN
                0165           hLimitNew(i,j,k) = hLimit(k) + dhActual(i,j,k+1)
                0166          ELSE
                0167           hLimitNew(i,j,k) = hLimit(k)
                0168          ENDIF
                0169 C     After computing the new boundary, check
                0170 C     (1) if it is between two adjacent thicknesses
                0171          IF ( ( AREAITD(i,j,k,bi,bj).GT.SEAICE_area_reg .AND.
                0172      &          hActual(i,j,k) .GE. hLimitNew(i,j,k) ) .OR.
                0173      &        ( AREAITD(i,j,k+1,bi,bj).GT.SEAICE_area_reg .AND.
dcd6ed0c75 Jean*0174      &          hActual(i,j,k+1) .LE. hLimitNew(i,j,k) ) )
ed2f6fecc4 Mart*0175      &        doRemapping(i,j) = .FALSE.
                0176 C     (2) that it is been the old boudnaries k-1 and k+1
                0177 C     (Note from CICE: we could allow this, but would make the code
                0178 C     more complicated)
                0179          IF ( ( hLimitNew(i,j,k) .GT. hLimit(k+1) ) .OR.
                0180      &        ( hLimitNew(i,j,k) .LT. hLimit(k-1) ) )
                0181      &        doRemapping(i,j) = .FALSE.
                0182         ENDDO
                0183        ENDDO
                0184       ENDDO
                0185 C     Report problems, if there are any. Because this breaks optimization
                0186 C     do not do it by default.
                0187 C     Where doRemapping is false, the rebinning of seaice_itd_redist
                0188 C     (called at the end) will take care of shifting the ice.
                0189       IF ( debugLevel.GE.debLevA )
                0190      &     CALL SEAICE_ITD_REMAP_CHECK_BOUNDS(
                0191      I     AREAITD, hActual, hActualPre, hLimitNew, doRemapping,
                0192      I     bi, bj, myTime, myIter, myThid )
                0193 C     computing the upper limit of the thickest category does not require
                0194 C     any checks and can be computed now
                0195       k = nITD
                0196       DO j=1,sNy
                0197        DO i=1,sNx
                0198         hLimitNew(i,j,k) = hLimit(k)
                0199         IF ( AREAITD(i,j,k,bi,bj).GT.SEAICE_area_reg )
                0200      &       hLimitNew(i,j,k) = MAX( 3. _d 0*hActual(i,j,k)
                0201      &       - 2. _d 0 * hLimitNew(i,j,k-1), hLimit(k-1) )
                0202        ENDDO
                0203       ENDDO
dcd6ed0c75 Jean*0204 C
                0205 C     end of limit computation, now compute the coefficients of the
ed2f6fecc4 Mart*0206 C     linear approximations of g(h) => g(eta) = g0 + g1*eta
dcd6ed0c75 Jean*0207 C
                0208 C     CICE does something specical for the first category.
ed2f6fecc4 Mart*0209 C     compute coefficients for 1st category
                0210       k = 1
                0211       DO j=1,sNy
                0212        DO i=1,sNx
                0213 C     initialisation
                0214         aLoc(i,j) = AREAITD(i,j,k,bi,bj)
                0215 C     initialise hL and hR
dcd6ed0c75 Jean*0216 C     this single line is different from the code that follows below
ed2f6fecc4 Mart*0217 C     for all categories
                0218         hL(i,j,k) = hLimitNew(i,j,k-1)
                0219         hR(i,j,k) = hLimit(k)
                0220        ENDDO
                0221       ENDDO
                0222       CALL SEAICE_ITD_REMAP_LINEAR(
dcd6ed0c75 Jean*0223      O     g0(1,1,k), g1(1,1,k),
ed2f6fecc4 Mart*0224      U     hL(1,1,k), hR(1,1,k),
dcd6ed0c75 Jean*0225      I     hActual(1,1,k), aLoc,
ed2f6fecc4 Mart*0226      I     SEAICE_area_reg, SEAICE_eps, doRemapping,
                0227      I     myTime, myIter, myThid )
                0228 C
                0229 C     Find area lost due to melting of thin (category 1) ice
                0230 C
                0231       DO j=1,sNy
                0232        DO i=1,sNx
dcd6ed0c75 Jean*0233         IF ( doRemapping(i,j) .AND.
ed2f6fecc4 Mart*0234      &       AREAITD(i,j,k,bi,bj) .GT. SEAICE_area_reg ) THEN
                0235 CMU if melting of ice in category 1
                0236          IF ( dhActual(i,j,k) .LT. 0. _d 0 ) THEN
                0237 C     integrate g(1) from zero to abs(dhActual)
                0238 CMU dh0 is max thickness of ice in first category that is melted
                0239           dh0    = MIN(-dhActual(i,j,k),hLimit(k))
                0240           etaMax = MIN(dh0,hR(i,j,k)) - hL(i,j,k)
                0241           IF ( etaMax > 0. _d 0 ) THEN
                0242 CMU da0 is /int_0^dh0 g dh
                0243            da0 = g0(i,j,k)*etaMax + g1(i,j,k)*etaMax*etaMax*0.5 _d 0
                0244            daMax = AREAITD(i,j,k,bi,bj)
                0245      &          * ( 1. _d 0 - hActual(i,j,k)/hActualPre(i,j,k))
                0246            da0 = MIN( da0, daMax )
                0247 CMU adjust thickness to conserve volume
                0248            IF ( (AREAITD(i,j,k,bi,bj)-da0) .GT. SEAICE_area_reg ) THEN
                0249              hActual(i,j,k) = hActual(i,j,k)
                0250      &            * AREAITD(i,j,k,bi,bj)/( AREAITD(i,j,k,bi,bj) - da0 )
                0251            ELSE
                0252              hActual(i,j,k) = ZERO
                0253              da0 = AREAITD(i,j,k,bi,bj)
                0254            ENDIF
                0255 CMU increase open water fraction
                0256            AREAITD(i,j,k,bi,bj) = AREAITD(i,j,k,bi,bj) - da0
                0257           ENDIF
                0258          ELSE
                0259 CMU H_0* = F_0 * dT
                0260           hLimitNew(i,j,k-1) = MIN( dhActual(i,j,k), hLimit(k) )
                0261          ENDIF
                0262         ENDIF
                0263        ENDDO
                0264       ENDDO
                0265 C
                0266 C     compute all coefficients
                0267 C
                0268       DO k=1,nITD
                0269        DO j=1,sNy
                0270         DO i=1,sNx
                0271 C     initialisation
                0272          aLoc(i,j) = AREAITD(i,j,k,bi,bj)
                0273 C     initialise hL and hR
                0274          hL(i,j,k) = hLimitNew(i,j,k-1)
                0275          hR(i,j,k) = hLimitNew(i,j,k)
                0276         ENDDO
                0277        ENDDO
                0278        CALL SEAICE_ITD_REMAP_LINEAR(
dcd6ed0c75 Jean*0279      O      g0(1,1,k), g1(1,1,k),
ed2f6fecc4 Mart*0280      U      hL(1,1,k), hR(1,1,k),
                0281      I      hActual(1,1,k), aLoc,
                0282      I      SEAICE_area_reg, SEAICE_eps, doRemapping,
                0283      I      myTime, myIter, myThid )
                0284       ENDDO
dcd6ed0c75 Jean*0285 C
ed2f6fecc4 Mart*0286       DO k=1,nITD-1
                0287        DO j=1,sNy
                0288         DO i=1,sNx
                0289          dheff = 0. _d 0
                0290          darea = 0. _d 0
                0291          IF ( doRemapping(i,j) ) THEN
                0292 C     compute integration limits in eta space
                0293           IF ( hLimitNew(i,j,k) .GT. hLimit(k) ) THEN
                0294            etaMin = MAX(       hLimit(k), hL(i,j,k)) - hL(i,j,k)
                0295            etaMax = MIN(hLimitNew(i,j,k), hR(i,j,k)) - hL(i,j,k)
                0296            kDonor = k
                0297            kRecvr = k+1
                0298           ELSE
                0299            etaMin = 0. _d 0
                0300            etaMax = MIN(hLimit(k), hR(i,j,k+1)) - hL(i,j,k+1)
                0301            kDonor = k+1
                0302            kRecvr = k
                0303           ENDIF
                0304 C     compute the area and volume to be moved
                0305           IF ( etaMax .GT. etaMin ) THEN
                0306            etam  = etaMax-etaMin
                0307            etap  = etaMax+etaMin
                0308            eta2  = 0.5*etam*etap
                0309            darea = g0(i,j,kDonor)*etam + g1(i,j,kDonor)*eta2
                0310 CML           dheff = g0(i,j,kDonor)*eta2
                0311 CML     &          +  g1(i,j,kDonor)*etam*(etap*etap-etaMax*etaMin)*third
                0312 CML     &          +  darea*hL(i,j,kDonor)
                0313            dheff = g0(i,j,kDonor)*eta2
                0314      &          +  g1(i,j,kDonor)*(etaMax**3-etaMin**3)*third
                0315      &          +  darea*hL(i,j,kDonor)
                0316           ENDIF
                0317 C     ... or shift entire category, if nearly all ice is to be shifted.
                0318 CMU          IF ( (darea .GT.AREAITD(i,j,kDonor,bi,bj)*oneMinusEps).OR.
                0319 CMU     &         (dheff .GT.HEFFITD(i,j,kDonor,bi,bj)*oneMinusEps) ) THEN
                0320           IF ( (darea .GT.AREAITD(i,j,kDonor,bi,bj)-SEAICE_eps).OR.
                0321      &         (dheff .GT.HEFFITD(i,j,kDonor,bi,bj)-SEAICE_eps) ) THEN
                0322            darea = AREAITD(i,j,kDonor,bi,bj)
                0323            dheff = HEFFITD(i,j,kDonor,bi,bj)
                0324           ENDIF
                0325 C     regularize: reset to zero, if there is too little ice to be shifted ...
                0326 CMU          IF ( (darea .LT. AREAITD(i,j,kDonor,bi,bj)*SEAICE_eps).OR.
                0327 CMU     &         (dheff .LT. HEFFITD(i,j,kDonor,bi,bj)*SEAICE_eps) ) THEN
                0328           IF ( (darea .LT. SEAICE_eps).OR.
                0329      &         (dheff .LT. SEAICE_eps) ) THEN
                0330            darea  = 0. _d 0
                0331            dheff  = 0. _d 0
                0332           ENDIF
                0333 C     snow scaled by area
                0334           IF ( AREAITD(i,j,kDonor,bi,bj) .GT. SEAICE_area_reg ) THEN
                0335 C     snow scaled by area (why not volume?), CICE also does it in this way
                0336            dhsnw = darea/AREAITD(i,j,kDonor,bi,bj)
                0337      &          * HSNOWITD(i,j,kDonor,bi,bj)
                0338 CMU          IF ( HEFFITD(i,j,kDonor,bi,bj) .GT. SEAICE_hice_reg ) THEN
dcd6ed0c75 Jean*0339 CMU           dhsnw = dheff/HEFFITD(i,j,kDonor,bi,bj)
ed2f6fecc4 Mart*0340 CMU     &         * HSNOWITD(i,j,kDonor,bi,bj)
                0341           ELSE
                0342            dhsnw = HSNOWITD(i,j,kDonor,bi,bj)
                0343           ENDIF
                0344 C     apply increments
                0345           HEFFITD(i,j,kRecvr,bi,bj) = HEFFITD(i,j,kRecvr,bi,bj) + dheff
                0346           HEFFITD(i,j,kDonor,bi,bj) = HEFFITD(i,j,kDonor,bi,bj) - dheff
                0347           AREAITD(i,j,kRecvr,bi,bj) = AREAITD(i,j,kRecvr,bi,bj) + darea
                0348           AREAITD(i,j,kDonor,bi,bj) = AREAITD(i,j,kDonor,bi,bj) - darea
                0349           HSNOWITD(i,j,kRecvr,bi,bj)=HSNOWITD(i,j,kRecvr,bi,bj) + dhsnw
                0350           HSNOWITD(i,j,kDonor,bi,bj)=HSNOWITD(i,j,kDonor,bi,bj) - dhsnw
                0351 C     end if doRemapping
                0352          ENDIF
                0353         ENDDO
                0354        ENDDO
                0355       ENDDO
                0356 
                0357       RETURN
                0358       END
                0359 
                0360 C---+-|--1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0361 
                0362 CBOP
                0363 C !ROUTINE: SEAICE_ITD_REMAP_LINEAR
                0364 
                0365 C !INTERFACE: ==========================================================
                0366       SUBROUTINE SEAICE_ITD_REMAP_LINEAR(
                0367      O     g0, g1,
                0368      U     hL, hR,
dcd6ed0c75 Jean*0369      I     hActual, area,
ed2f6fecc4 Mart*0370      I     SEAICE_area_reg, SEAICE_eps, doRemapping,
                0371      I     myTime, myIter, myThid )
                0372 
                0373 C !DESCRIPTION: \bv
                0374 C     *===========================================================*
                0375 C     | SUBROUTINE SEAICE_ITD_REMAP_LINEAR
dcd6ed0c75 Jean*0376 C     | o compute coefficients g0, g1 for piece-wise linear fit
ed2f6fecc4 Mart*0377 C     |    g(h) = g0 + g1*h
                0378 C     | o compute range boundaries hL, hR for this linear fit
                0379 C     |
                0380 C     | Martin Losch, May 2014, Martin.Losch@awi.de
                0381 C     *===========================================================*
                0382 C \ev
                0383 
                0384 C !USES: ===============================================================
                0385       IMPLICIT NONE
                0386 
                0387 #include "SIZE.h"
                0388 
                0389 C !INPUT PARAMETERS: ===================================================
                0390 C     === Routine arguments ===
                0391 C     myTime    :: current time
                0392 C     myIter    :: iteration number
                0393 C     myThid    :: Thread no. that called this routine.
                0394       _RL myTime
                0395       INTEGER myIter
                0396       INTEGER myThid
dcd6ed0c75 Jean*0397 C
ed2f6fecc4 Mart*0398 C     OUTPUT: coefficients for representing g(h)
                0399 C     g0 :: constant coefficient in g(h)
                0400 C     g1 :: linear  coefficient in g(h)
                0401 C     hL :: left end of range over which g(h) > 0
                0402 C     hL :: right end of range over which g(h) > 0
                0403       _RL g0 (1:sNx,1:sNy)
                0404       _RL g1 (1:sNx,1:sNy)
                0405       _RL hL (1:sNx,1:sNy)
                0406       _RL hR (1:sNx,1:sNy)
                0407 C     INPUT:
                0408 C     hActual :: ice thickness of current category
                0409 C     area    :: ice concentration of current category
                0410       _RL hActual (1:sNx,1:sNy)
                0411       _RL area    (1:sNx,1:sNy)
                0412 C     regularization constants
                0413       _RL SEAICE_area_reg
                0414       _RL SEAICE_eps
                0415 C     doRemapping :: mask where can be done, excludes points where
                0416 C                    new category limits are outside certain bounds
                0417       LOGICAL doRemapping (1:sNx,1:sNy)
                0418 
                0419 C !LOCAL VARIABLES: ====================================================
                0420 C     === Local variables ===
                0421 C     i,j       :: inner loop counters
                0422 C
                0423       INTEGER i, j
                0424 C     auxCoeff :: helper variable
                0425 C     recip_etaR :: reciprocal of range interval in eta space
                0426 C     etaNoR   :: ratio of distance to lower limit over etaR
                0427       _RL auxCoeff
                0428       _RL recip_etaR, etaNoR
                0429       _RL third, sixth
                0430       PARAMETER ( third = 0.333333333333333333333333333 _d 0 )
                0431       PARAMETER ( sixth = 0.666666666666666666666666666 _d 0 )
                0432 CEOP
dcd6ed0c75 Jean*0433 C
ed2f6fecc4 Mart*0434 C     initialisation of hL, hR is done outside this routine
                0435 C
                0436       DO j=1,sNy
                0437        DO i=1,sNx
                0438         g0(i,j) = 0. _d 0
                0439         g1(i,j) = 0. _d 0
                0440         IF ( doRemapping(i,j) .AND.
                0441      &       area(i,j) .GT. SEAICE_area_reg .AND.
                0442      &       hR(i,j) - hL(i,j) .GT. SEAICE_eps ) THEN
                0443 C     change hL and hR if hActual falls outside the central third of the range
                0444          IF ( hActual(i,j) .LT. (2. _d 0*hL(i,j) + hR(i,j))*third ) THEN
                0445           hR(i,j) = 3. _d 0 * hActual(i,j) - 2. _d 0 * hL(i,j)
                0446          ELSEIF ( hActual(i,j).GT.(hL(i,j)+2. _d 0*hR(i,j))*third ) THEN
                0447           hL(i,j) = 3. _d 0 * hActual(i,j) - 2. _d 0 * hR(i,j)
                0448          ENDIF
                0449 C     calculate new etaR = hR - hL;
                0450 C     catch the case of hR=hL, which can happen when hActual=hR or hL
                0451 C     before entering this routine; in this case g0=g1=0.
                0452          recip_etaR = 0. _d 0
                0453 CMU         IF ( hR(i,j) .GT. hL(i,j) ) ! crucial change; lets the model explode
                0454          IF ( hR(i,j) - hL(i,j) .GT. SEAICE_eps )
                0455      &        recip_etaR = 1. _d 0 / (hR(i,j) - hL(i,j))
                0456 C     some abbreviations to avoid computing the same thing multiple times
                0457          etaNoR     = (hActual(i,j) - hL(i,j))*recip_etaR
                0458          auxCoeff   = 6. _d 0 * area(i,j)*recip_etaR
                0459 C     equations (14) of Lipscomb (2001), JGR
                0460          g0(i,j) = auxCoeff*( sixth - etaNoR )
                0461          g1(i,j) = 2. _d 0 * auxCoeff*recip_etaR*( etaNoR - 0.5 _d 0 )
                0462         ELSE
                0463 C     not doRemapping
                0464 C     reset hL and hR
                0465          hL(i,j) = 0. _d 0
                0466          hR(i,j) = 0. _d 0
                0467         ENDIF
                0468        ENDDO
                0469       ENDDO
                0470 
                0471       RETURN
                0472       END
                0473 
                0474 C---+-|--1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0475 
e7af59f6fd Jean*0476 CBOP
ed2f6fecc4 Mart*0477 C !ROUTINE: SEAICE_ITD_REMAP_CHECK_BOUNDS
                0478 
                0479 C !INTERFACE: ==========================================================
                0480       SUBROUTINE SEAICE_ITD_REMAP_CHECK_BOUNDS(
                0481      I     AREAITD, hActual, hActualPre, hLimitNew, doRemapping,
                0482      I     bi, bj, myTime, myIter, myThid )
                0483 
                0484 C !DESCRIPTION: \bv
                0485 C     *===========================================================*
                0486 C     | SUBROUTINE SEAICE_ITD_REMAP_CHECK_BOUNDS
                0487 C     | o where doRemapping = .FALSE. print a warning
                0488 C     |
                0489 C     | Martin Losch, May 2014, Martin.Losch@awi.de
                0490 C     *===========================================================*
                0491 C \ev
                0492 
                0493 C !USES: ===============================================================
                0494       IMPLICIT NONE
                0495 
                0496 #include "SIZE.h"
                0497 #include "EEPARAMS.h"
                0498 #include "SEAICE_SIZE.h"
                0499 #include "SEAICE_PARAMS.h"
                0500 
                0501 C !INPUT PARAMETERS: ===================================================
                0502 C     === Routine arguments ===
                0503 C     bi, bj    :: outer loop counters
                0504 C     myTime    :: current time
                0505 C     myIter    :: iteration number
                0506 C     myThid    :: Thread no. that called this routine.
                0507       _RL myTime
                0508       INTEGER bi,bj
                0509       INTEGER myIter
                0510       INTEGER myThid
                0511 C     hActual :: ice thickness of current category
                0512       _RL hActual   (1:sNx,1:sNy,1:nITD)
                0513       _RL hActualPre(1:sNx,1:sNy,1:nITD)
                0514 C     hLimitNew :: new "advected" category boundaries after seaice_growth
                0515       _RL hLimitNew (1:sNx,1:sNy,0:nITD)
                0516 C     AREAITD :: ice concentration of current category
dcd6ed0c75 Jean*0517       _RL AREAITD   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nITD,nSx,nSy)
ed2f6fecc4 Mart*0518 C     doRemapping :: mask where can be done, excludes points where
                0519 C                    new category limits are outside certain bounds
                0520       LOGICAL doRemapping (1:sNx,1:sNy)
                0521 
                0522 C !LOCAL VARIABLES: ====================================================
                0523 C     === Local variables ===
                0524 C     i,j,k     :: inner loop counters
                0525 C
                0526       INTEGER i, j, k
                0527       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0528       CHARACTER*(39) tmpBuf
                0529 CEOP
dcd6ed0c75 Jean*0530 
ed2f6fecc4 Mart*0531        DO j=1,sNy
                0532         DO i=1,sNx
                0533          IF (.NOT.doRemapping(i,j) ) THEN
                0534           DO k=1,nITD-1
dcd6ed0c75 Jean*0535            WRITE(tmpBuf,'(A,2I5,A,I10)')
ed2f6fecc4 Mart*0536      &          ' at (', i, j, ') in timestep ', myIter
                0537            IF ( AREAITD(i,j,k,bi,bj).GT.SEAICE_area_reg .AND.
                0538      &          hActual(i,j,k) .GE. hLimitNew(i,j,k) ) THEN
                0539             WRITE(msgBuf,'(A,I3,A)')
                0540      &           'SEAICE_ITD_REMAP: hActual(k) >= hLimitNew(k) '//
                0541      &           'for category ', k, tmpBuf
                0542             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0543      &           SQUEEZE_RIGHT, myThid )
dcd6ed0c75 Jean*0544 CML            PRINT *, hActual(i,j,k),
ed2f6fecc4 Mart*0545 CML     &           hLimitNew(i,j,k), hLimit(k)
                0546            ENDIF
                0547            IF ( AREAITD(i,j,k+1,bi,bj).GT.SEAICE_area_reg .AND.
                0548      &          hActual(i,j,k+1) .LE. hLimitNew(i,j,k) ) THEN
dcd6ed0c75 Jean*0549             WRITE(msgBuf,'(A,I3,A)')
ed2f6fecc4 Mart*0550      &           'SEAICE_ITD_REMAP: hActual(k+1) <= hLimitNew(k) '//
                0551      &           'for category ', k, tmpBuf
                0552             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0553      &           SQUEEZE_RIGHT, myThid )
dcd6ed0c75 Jean*0554             PRINT '(8(1X,E10.4))',
ed2f6fecc4 Mart*0555      &           AREAITD(i,j,k+1,bi,bj), hActual(i,j,k+1),
                0556      &           hActualPre(i,j,k+1),
                0557      &           AREAITD(i,j,k,bi,bj), hActual(i,j,k),
                0558      &           hActualPre(i,j,k),
                0559      &           hLimitNew(i,j,k), hLimit(k)
                0560            ENDIF
                0561            IF ( hLimitNew(i,j,k) .GT. hLimit(k+1) ) THEN
dcd6ed0c75 Jean*0562             WRITE(msgBuf,'(A,I3,A)')
ed2f6fecc4 Mart*0563      &           'SEAICE_ITD_REMAP: hLimitNew(k) > hLimitNew(k+1) '//
                0564      &           'for category ', k, tmpBuf
                0565             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0566      &           SQUEEZE_RIGHT, myThid )
                0567            ENDIF
                0568            IF ( hLimitNew(i,j,k) .LT. hLimit(k-1) ) THEN
dcd6ed0c75 Jean*0569             WRITE(msgBuf,'(A,I3,A)')
ed2f6fecc4 Mart*0570      &           'SEAICE_ITD_REMAP: hLimitNew(k) < hLimitNew(k-1) '//
                0571      &           'for category ', k, tmpBuf
                0572             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0573      &           SQUEEZE_RIGHT, myThid )
                0574            ENDIF
                0575           ENDDO
                0576          ENDIF
                0577         ENDDO
                0578        ENDDO
                0579 
                0580 #endif /* SEAICE_ITD */
                0581 
                0582       RETURN
                0583       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!