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