pkg/seaice/seaice_prepare_ridging.F

0c32bd3cb0 Mart*0001 #include "SEAICE_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: SEAICE_PREPARE_RIDGING
                0005 C !INTERFACE: ==========================================================
                0006       SUBROUTINE SEAICE_PREPARE_RIDGING(
                0007 #ifdef SEAICE_ITD
353a8877c7 Mart*0008      O     hActual,
0c32bd3cb0 Mart*0009      O     hrMin, hrMax, hrExp, ridgeRatio, ridgingModeNorm, partFunc,
                0010 #endif /* SEAICE_ITD */
                0011      I     iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
                0012
                0013 C !DESCRIPTION: \bv
                0014 C     *===========================================================*
                0015 C     | SUBROUTINE SEAICE_PREPARE_RIDGING
4e4ad91a39 Jean*0016 C     | o compute ridging parameters according to Thorndyke et al
0c32bd3cb0 Mart*0017 C     |   (1975), Hibler (1980), Bitz et al (2001) or
                0018 C     |   Lipscomb et al (2007)
4e4ad91a39 Jean*0019 C     | o this routine is called from s/r seaice_do_ridging and
0c32bd3cb0 Mart*0020 C     |   from s/r seaice_calc_ice_strength
4e4ad91a39 Jean*0021 C     |
0c32bd3cb0 Mart*0022 C     | Martin Losch, Apr. 2014, Martin.Losch@awi.de
                0023 C     *===========================================================*
                0024 C \ev
                0025
                0026 C !USES: ===============================================================
                0027       IMPLICIT NONE
                0028
                0029 #include "SIZE.h"
                0030 #include "EEPARAMS.h"
                0031 #include "PARAMS.h"
                0032 #include "GRID.h"
                0033 #include "SEAICE_SIZE.h"
                0034 #include "SEAICE_PARAMS.h"
3f0f10fc37 Mart*0035 #include "SEAICE_GRID.h"
0c32bd3cb0 Mart*0036 #include "SEAICE.h"
                0037
                0038 C !INPUT PARAMETERS: ===================================================
                0039 C     === Routine arguments ===
                0040 C     bi, bj    :: outer loop counters
                0041 C     myTime    :: current time
                0042 C     myIter    :: iteration number
                0043 C     myThid    :: Thread no. that called this routine.
                0044 C     i/jMin/Max:: loop boundaries
                0045       _RL myTime
                0046       INTEGER bi,bj
                0047       INTEGER myIter
                0048       INTEGER myThid
                0049       INTEGER iMin, iMax, jMin, jMax
                0050 #ifdef SEAICE_ITD
                0051 C     ridgingModeNorm :: norm to ensure convervation (N in Lipscomb et al 2007)
                0052 C     partFunc   :: participation function (a_n in Lipscomb et al 2007)
                0053 C     ridgeRatio :: mean ridge thickness/ thickness of ridging ice
                0054 C     hrMin      :: min ridge thickness
                0055 C     hrMax      :: max ridge thickness   (SEAICEredistFunc = 0)
                0056 C     hrExp      :: ridge e-folding scale (SEAICEredistFunc = 1)
353a8877c7 Mart*0057 C     hActual    :: HEFFITD/AREAITD, regularized
0c32bd3cb0 Mart*0058       _RL ridgingModeNorm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0059       _RL partFunc        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,0:nITD)
                0060       _RL hrMin           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
                0061       _RL hrMax           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
                0062       _RL hrExp           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
                0063       _RL ridgeRatio      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
4e4ad91a39 Jean*0064       _RL hActual         (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
0c32bd3cb0 Mart*0065 CEndOfInterface
                0066
                0067 C !LOCAL VARIABLES: ====================================================
                0068 C     === Local variables ===
                0069 C     i,j,k       :: inner loop counters
                0070 C
                0071       INTEGER i, j
                0072       INTEGER k
                0073 C     variables related to ridging schemes
                0074 C     gSum        :: cumulative distribution function G
                0075       _RL gSum            (1-OLx:sNx+OLx,1-OLy:sNy+OLy,-1:nITD)
                0076       _RL recip_gStar, recip_aStar, tmp
353a8877c7 Mart*0077 C     Regularization values squared
                0078       _RL area_reg_sq, hice_reg_sq
0c32bd3cb0 Mart*0079 CEOP
                0080
                0081 C---+-|--1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0082
353a8877c7 Mart*0083 C     regularization constants
                0084       area_reg_sq = SEAICE_area_reg * SEAICE_area_reg
                0085       hice_reg_sq = SEAICE_hice_reg * SEAICE_hice_reg
0c32bd3cb0 Mart*0086       DO k=1,nITD
                0087        DO j=jMin,jMax
                0088         DO i=iMin,iMax
                0089          hActual(i,j,k) = 0. _d 0
353a8877c7 Mart*0090 CML         IF ( AREAITD(i,j,k,bi,bj) .GT. SEAICE_area_reg ) THEN
                0091 CML          hActual(i,j,k) = HEFFITD(i,j,k,bi,bj)/AREAITD(i,j,k,bi,bj)
                0092 CML         ENDIF
                0093          IF ( HEFFITD(i,j,k,bi,bj) .GT. 0. _d 0 ) THEN
4e4ad91a39 Jean*0094 C     regularize as in seaice_growth: compute hActual with regularized
353a8877c7 Mart*0095 C     AREA and regularize from below with a minimum thickness
                0096           tmp = HEFFITD(i,j,k,bi,bj)
e4863bf4ee Mart*0097      &         /SQRT( AREAITD(i,j,k,bi,bj)**2 + area_reg_sq )
353a8877c7 Mart*0098           hActual(i,j,k) = SQRT(tmp * tmp + hice_reg_sq)
0c32bd3cb0 Mart*0099          ENDIF
                0100         ENDDO
                0101        ENDDO
                0102       ENDDO
4e4ad91a39 Jean*0103
0c32bd3cb0 Mart*0104 C---+-|--1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0105
                0106 C     compute the cumulative thickness distribution function gSum
                0107       DO j=jMin,jMax
                0108        DO i=iMin,iMax
                0109         gSum(i,j,-1) = 0. _d 0
                0110         gSum(i,j,0)  = 0. _d 0
4e4ad91a39 Jean*0111         IF ( opnWtrFrac(i,j,bi,bj) .GT. SEAICE_area_floor )
8d0c6ad10c Mart*0112      &       gSum(i,j,0) = opnWtrFrac(i,j,bi,bj)
0c32bd3cb0 Mart*0113        ENDDO
                0114       ENDDO
                0115       DO k = 1, nITD
                0116        DO j=jMin,jMax
                0117         DO i=iMin,iMax
                0118          gSum(i,j,k) = gSum(i,j,k-1)
                0119          IF ( AREAITD(i,j,k,bi,bj) .GT. SEAICE_area_floor )
                0120      &        gSum(i,j,k) = gSum(i,j,k) + AREAITD(i,j,k,bi,bj)
                0121         ENDDO
                0122        ENDDO
                0123       ENDDO
                0124 C     normalize
                0125       DO k = 0, nITD
                0126        DO j=jMin,jMax
                0127         DO i=iMin,iMax
4e4ad91a39 Jean*0128          IF ( gSum(i,j,nITD).NE.0. _d 0 )
0c32bd3cb0 Mart*0129      &        gSum(i,j,k) = gSum(i,j,k) / gSum(i,j,nITD)
                0130         ENDDO
                0131        ENDDO
                0132       ENDDO
                0133
                0134 C     Compute the participation function
                0135 C                    area lost from category n due to ridging/closing
                0136 C     partFunc(n) = --------------------------------------------------
                0137 C                       total area lost due to ridging/closing
                0138
                0139       IF ( SEAICEpartFunc .EQ. 0 ) THEN
                0140 C     Thorndike et al. (1975) discretize b(h) = (2/Gstar) * (1 - G(h)/Gstar)
                0141 C     The expressions for the partition function partFunc are found by
                0142 C     integrating b(h)g(h) between the category boundaries.
                0143        recip_gStar = 1. _d 0 / SEAICEgStar
                0144        DO k = 0, nITD
                0145         DO j=jMin,jMax
                0146          DO i=iMin,iMax
                0147           partFunc(i,j,k) = 0. _d 0
                0148           IF ( gSum(i,j,k) .LT. SEAICEgStar ) THEN
4e4ad91a39 Jean*0149            partFunc(i,j,k) =
0c32bd3cb0 Mart*0150      &          (gSum(i,j,k)-gSum(i,j,k-1)) * recip_gStar
                0151      &          *( 2. _d 0 - (gSum(i,j,k-1)+gSum(i,j,k))*recip_gStar)
4e4ad91a39 Jean*0152           ELSEIF (  gSum(i,j,k-1) .LT. SEAICEgStar
0c32bd3cb0 Mart*0153      &          .AND. gSum(i,j,k) .GE. SEAICEgStar ) THEN
4e4ad91a39 Jean*0154            partFunc(i,j,k) =
0c32bd3cb0 Mart*0155      &          (SEAICEgStar-gSum(i,j,k-1)) * recip_gStar
                0156      &          *( 2. _d 0 - (gSum(i,j,k-1)+SEAICEgStar)*recip_gStar)
                0157           ENDIF
                0158          ENDDO
                0159         ENDDO
                0160        ENDDO
                0161       ELSEIF  ( SEAICEpartFunc .EQ. 1 ) THEN
                0162 C     Lipscomb et al. (2007) discretize b(h) = exp(-G(h)/astar) into
4e4ad91a39 Jean*0163 C     partFunc(n) = [exp(-G(n-1)/astar - exp(-G(n)/astar] / [1-exp(-1/astar)].
                0164 C     The expression is found by integrating b(h)g(h) between the category
0c32bd3cb0 Mart*0165 C     boundaries.
                0166        recip_astar = 1. _d 0 / SEAICEaStar
                0167        tmp = 1. _d 0 / ( 1. _d 0 - EXP( -recip_astar ) )
                0168 C     abuse gSum as a work array
                0169        k = -1
                0170        DO j=jMin,jMax
                0171         DO i=iMin,iMax
                0172          gSum(i,j,k)     = EXP(-gSum(i,j,k)*recip_astar) * tmp
                0173         ENDDO
                0174        ENDDO
                0175        DO k = 0, nITD
                0176         DO j=jMin,jMax
                0177          DO i=iMin,iMax
                0178           gSum(i,j,k)     = EXP(-gSum(i,j,k)*recip_astar) * tmp
                0179           partFunc(i,j,k) = gSum(i,j,k-1) - gSum(i,j,k)
                0180          ENDDO
                0181         ENDDO
                0182        ENDDO
                0183       ELSE
                0184        STOP 'Ooops: SEAICEpartFunc > 1 not implemented'
                0185       ENDIF
                0186
                0187 C     Compute variables of ITD of ridged ice
                0188 C     ridgeRatio :: mean ridge thickness/ thickness of ridging ice
                0189 C     hrMin      :: min ridge thickness
                0190 C     hrMax      :: max ridge thickness   (SEAICEredistFunc = 0)
                0191 C     hrExp      :: ridge e-folding scale (SEAICEredistFunc = 1)
                0192       DO k = 1, nITD
                0193        DO j=jMin,jMax
                0194         DO i=iMin,iMax
                0195          hrMin(i,j,k)      = 0. _d 0
                0196          hrMax(i,j,k)      = 0. _d 0
                0197          hrExp(i,j,k)      = 0. _d 0
                0198 C     avoid divisions by zero
                0199          ridgeRatio(i,j,k) = 1. _d 0
                0200         ENDDO
                0201        ENDDO
                0202       ENDDO
                0203       IF ( SEAICEredistFunc .EQ. 0 ) THEN
                0204 C     Assume ridged ice is uniformly distributed between hrmin and hrmax.
                0205 C     (Hibler, 1980)
                0206        DO k = 1, nITD
                0207         DO j=jMin,jMax
                0208          DO i=iMin,iMax
                0209           IF ( hActual(i,j,k) .GT. 0. _d 0 ) THEN
                0210 C     This is the original Hibler (1980) scheme:
                0211            hrMin(i,j,k) = 2. _d 0 * hActual(i,j,k)
                0212            hrMax(i,j,k) = 2. _d 0 * SQRT(hActual(i,j,k)*SEAICEhStar)
                0213 C     CICE does this in addition, so that thick ridging ice is not required
4e4ad91a39 Jean*0214 C     to raft:
0c32bd3cb0 Mart*0215            hrMin(i,j,k) = MIN(hrMin(i,j,k),hActual(i,j,k)+SEAICEmaxRaft)
                0216            hrMax(i,j,k) = MAX(hrMax(i,j,k),hrMin(i,j,k)+SEAICE_hice_reg)
                0217 C
                0218            ridgeRatio(i,j,k) =
                0219      &          0.5 _d 0 * (hrMax(i,j,k)+hrMin(i,j,k))/hActual(i,j,k)
                0220           ENDIF
                0221          ENDDO
                0222         ENDDO
                0223        ENDDO
                0224       ELSEIF ( SEAICEredistFunc .EQ. 1 ) THEN
                0225 C     Follow Lipscomb et al. (2007) and model ridge ITD as an exponentially
                0226 C     decaying function
                0227        DO k = 1, nITD
                0228         DO j=jMin,jMax
                0229          DO i=iMin,iMax
                0230           IF ( hActual(i,j,k) .GT. 0. _d 0 ) THEN
353a8877c7 Mart*0231 C     regularization is only required in this case but already done above
                0232 CML           tmp = MAX(hActual(i,j,k), SEAICE_hice_reg)
                0233            tmp = hActual(i,j,k)
                0234            hrMin(i,j,k) = MIN(2.D0 * tmp, tmp+SEAICEmaxRaft)
                0235            hrExp(i,j,k) = SEAICEmuRidging*SQRT(tmp)
0c32bd3cb0 Mart*0236 C     arent we missing a factor 0.5 here?
353a8877c7 Mart*0237            ridgeRatio(i,j,k)=(hrMin(i,j,k)+hrExp(i,j,k))/tmp
0c32bd3cb0 Mart*0238           ENDIF
                0239          ENDDO
                0240         ENDDO
                0241        ENDDO
                0242       ELSE
4e4ad91a39 Jean*0243        STOP 'Ooops: SEAICEredistFunc > 1 not implemented'
0c32bd3cb0 Mart*0244       ENDIF
                0245
4e4ad91a39 Jean*0246 C     Compute the norm of the ridging mode N (in Lipscomp et al 2007)
                0247 C     or omega (in Bitz et al 2001):
0c32bd3cb0 Mart*0248 C     rigdingModeNorm = net ice area removed / total area participating.
                0249 C     For instance, if a unit area of ice with thickness = 1 participates in
                0250 C     ridging to form a ridge with a = 1/3 and thickness = 3, then
                0251 C     rigdingModeNorm = 1 - 1/3 = 2/3.
                0252       DO j=jMin,jMax
                0253        DO i=iMin,iMax
                0254         ridgingModeNorm(i,j) = partFunc(i,j,0)
                0255        ENDDO
                0256       ENDDO
                0257       DO k = 1, nITD
                0258        DO j=jMin,jMax
                0259         DO i=iMin,iMax
8d0c6ad10c Mart*0260          partFunc(i,j,k) = partFunc(i,j,k) * heffM(i,j,bi,bj)
0c32bd3cb0 Mart*0261          ridgingModeNorm(i,j) = ridgingModeNorm(i,j)
                0262      &        + partFunc(i,j,k)*( 1. _d 0 - 1. _d 0/ridgeRatio(i,j,k) )
                0263         ENDDO
                0264        ENDDO
                0265       ENDDO
8d0c6ad10c Mart*0266 C     avoid division by zero
                0267       DO j=jMin,jMax
                0268        DO i=iMin,iMax
                0269         IF ( ridgingModeNorm(i,j) .LE. 0. _d 0 )
                0270      &       ridgingModeNorm(i,j) = 1. _d 0
                0271        ENDDO
                0272       ENDDO
0c32bd3cb0 Mart*0273
                0274 #endif /* SEAICE_ITD */
                0275
                0276       RETURN
                0277       END