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