Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-28 05:09:44 UTC

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