C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_prepare_ridging.F,v 1.5 2014/10/20 03:20:58 gforget Exp $ C $Name: $ #include "SEAICE_OPTIONS.h" #ifdef ALLOW_AUTODIFF # include "AUTODIFF_OPTIONS.h" #endif CBOP C !ROUTINE: SEAICE_PREPARE_RIDGING C !INTERFACE: ========================================================== SUBROUTINE SEAICE_PREPARE_RIDGING( #ifdef SEAICE_ITD O hActual, O hrMin, hrMax, hrExp, ridgeRatio, ridgingModeNorm, partFunc, #endif /* SEAICE_ITD */ I iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid ) C !DESCRIPTION: \bv C *===========================================================* C | SUBROUTINE SEAICE_PREPARE_RIDGING C | o compute ridging parameters according to Thorndyke et al C | (1975), Hibler (1980), Bitz et al (2001) or C | Lipscomb et al (2007) C | o this routine is called from s/r seaice_do_ridging and C | from s/r seaice_calc_ice_strength C | C | Martin Losch, Apr. 2014, Martin.Losch@awi.de C *===========================================================* C \ev C !USES: =============================================================== IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "SEAICE_SIZE.h" #include "SEAICE_PARAMS.h" #include "SEAICE.h" #ifdef ALLOW_AUTODIFF_TAMC # include "tamc.h" #endif 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. C i/jMin/Max:: loop boundaries _RL myTime INTEGER bi,bj INTEGER myIter INTEGER myThid INTEGER iMin, iMax, jMin, jMax #ifdef SEAICE_ITD C ridgingModeNorm :: norm to ensure convervation (N in Lipscomb et al 2007) C partFunc :: participation function (a_n in Lipscomb et al 2007) C ridgeRatio :: mean ridge thickness/ thickness of ridging ice C hrMin :: min ridge thickness C hrMax :: max ridge thickness (SEAICEredistFunc = 0) C hrExp :: ridge e-folding scale (SEAICEredistFunc = 1) C hActual :: HEFFITD/AREAITD, regularized _RL ridgingModeNorm (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL partFunc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,0:nITD) _RL hrMin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD) _RL hrMax (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD) _RL hrExp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD) _RL ridgeRatio (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD) _RL hActual (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD) CEndOfInterface C !LOCAL VARIABLES: ==================================================== C === Local variables === C i,j,k :: inner loop counters C INTEGER i, j INTEGER k C variables related to ridging schemes C gSum :: cumulative distribution function G _RL gSum (1-OLx:sNx+OLx,1-OLy:sNy+OLy,-1:nITD) _RL recip_gStar, recip_aStar, tmp C Regularization values squared _RL area_reg_sq, hice_reg_sq CEOP C---+-|--1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C regularization constants area_reg_sq = SEAICE_area_reg * SEAICE_area_reg hice_reg_sq = SEAICE_hice_reg * SEAICE_hice_reg DO k=1,nITD DO j=jMin,jMax DO i=iMin,iMax hActual(i,j,k) = 0. _d 0 CML IF ( AREAITD(i,j,k,bi,bj) .GT. SEAICE_area_reg ) THEN CML hActual(i,j,k) = HEFFITD(i,j,k,bi,bj)/AREAITD(i,j,k,bi,bj) CML ENDIF IF ( HEFFITD(i,j,k,bi,bj) .GT. 0. _d 0 ) THEN C regularize as in seaice_growth: compute hActual with regularized C AREA and regularize from below with a minimum thickness tmp = HEFFITD(i,j,k,bi,bj) & /SQRT( AREAITD(i,j,k,bi,bj)**2 + area_reg_sq ) hActual(i,j,k) = SQRT(tmp * tmp + hice_reg_sq) ENDIF ENDDO ENDDO ENDDO C---+-|--1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C compute the cumulative thickness distribution function gSum DO j=jMin,jMax DO i=iMin,iMax gSum(i,j,-1) = 0. _d 0 gSum(i,j,0) = 0. _d 0 IF ( opnWtrFrac(i,j,bi,bj) .GT. SEAICE_area_floor ) & gSum(i,j,0) = opnWtrFrac(i,j,bi,bj) ENDDO ENDDO DO k = 1, nITD DO j=jMin,jMax DO i=iMin,iMax gSum(i,j,k) = gSum(i,j,k-1) IF ( AREAITD(i,j,k,bi,bj) .GT. SEAICE_area_floor ) & gSum(i,j,k) = gSum(i,j,k) + AREAITD(i,j,k,bi,bj) ENDDO ENDDO ENDDO C normalize DO k = 0, nITD DO j=jMin,jMax DO i=iMin,iMax IF ( gSum(i,j,nITD).NE.0. _d 0 ) & gSum(i,j,k) = gSum(i,j,k) / gSum(i,j,nITD) ENDDO ENDDO ENDDO C Compute the participation function C area lost from category n due to ridging/closing C partFunc(n) = -------------------------------------------------- C total area lost due to ridging/closing IF ( SEAICEpartFunc .EQ. 0 ) THEN C Thorndike et al. (1975) discretize b(h) = (2/Gstar) * (1 - G(h)/Gstar) C The expressions for the partition function partFunc are found by C integrating b(h)g(h) between the category boundaries. recip_gStar = 1. _d 0 / SEAICEgStar DO k = 0, nITD DO j=jMin,jMax DO i=iMin,iMax partFunc(i,j,k) = 0. _d 0 IF ( gSum(i,j,k) .LT. SEAICEgStar ) THEN partFunc(i,j,k) = & (gSum(i,j,k)-gSum(i,j,k-1)) * recip_gStar & *( 2. _d 0 - (gSum(i,j,k-1)+gSum(i,j,k))*recip_gStar) ELSEIF ( gSum(i,j,k-1) .LT. SEAICEgStar & .AND. gSum(i,j,k) .GE. SEAICEgStar ) THEN partFunc(i,j,k) = & (SEAICEgStar-gSum(i,j,k-1)) * recip_gStar & *( 2. _d 0 - (gSum(i,j,k-1)+SEAICEgStar)*recip_gStar) ENDIF ENDDO ENDDO ENDDO ELSEIF ( SEAICEpartFunc .EQ. 1 ) THEN C Lipscomb et al. (2007) discretize b(h) = exp(-G(h)/astar) into C partFunc(n) = [exp(-G(n-1)/astar - exp(-G(n)/astar] / [1-exp(-1/astar)]. C The expression is found by integrating b(h)g(h) between the category C boundaries. recip_astar = 1. _d 0 / SEAICEaStar tmp = 1. _d 0 / ( 1. _d 0 - EXP( -recip_astar ) ) C abuse gSum as a work array k = -1 DO j=jMin,jMax DO i=iMin,iMax gSum(i,j,k) = EXP(-gSum(i,j,k)*recip_astar) * tmp ENDDO ENDDO DO k = 0, nITD DO j=jMin,jMax DO i=iMin,iMax gSum(i,j,k) = EXP(-gSum(i,j,k)*recip_astar) * tmp partFunc(i,j,k) = gSum(i,j,k-1) - gSum(i,j,k) ENDDO ENDDO ENDDO ELSE STOP 'Ooops: SEAICEpartFunc > 1 not implemented' ENDIF C Compute variables of ITD of ridged ice C ridgeRatio :: mean ridge thickness/ thickness of ridging ice C hrMin :: min ridge thickness C hrMax :: max ridge thickness (SEAICEredistFunc = 0) C hrExp :: ridge e-folding scale (SEAICEredistFunc = 1) DO k = 1, nITD DO j=jMin,jMax DO i=iMin,iMax hrMin(i,j,k) = 0. _d 0 hrMax(i,j,k) = 0. _d 0 hrExp(i,j,k) = 0. _d 0 C avoid divisions by zero ridgeRatio(i,j,k) = 1. _d 0 ENDDO ENDDO ENDDO IF ( SEAICEredistFunc .EQ. 0 ) THEN C Assume ridged ice is uniformly distributed between hrmin and hrmax. C (Hibler, 1980) DO k = 1, nITD DO j=jMin,jMax DO i=iMin,iMax IF ( hActual(i,j,k) .GT. 0. _d 0 ) THEN C This is the original Hibler (1980) scheme: hrMin(i,j,k) = 2. _d 0 * hActual(i,j,k) hrMax(i,j,k) = 2. _d 0 * SQRT(hActual(i,j,k)*SEAICEhStar) C CICE does this in addition, so that thick ridging ice is not required C to raft: hrMin(i,j,k) = MIN(hrMin(i,j,k),hActual(i,j,k)+SEAICEmaxRaft) hrMax(i,j,k) = MAX(hrMax(i,j,k),hrMin(i,j,k)+SEAICE_hice_reg) C ridgeRatio(i,j,k) = & 0.5 _d 0 * (hrMax(i,j,k)+hrMin(i,j,k))/hActual(i,j,k) ENDIF ENDDO ENDDO ENDDO ELSEIF ( SEAICEredistFunc .EQ. 1 ) THEN C Follow Lipscomb et al. (2007) and model ridge ITD as an exponentially C decaying function DO k = 1, nITD DO j=jMin,jMax DO i=iMin,iMax IF ( hActual(i,j,k) .GT. 0. _d 0 ) THEN C regularization is only required in this case but already done above CML tmp = MAX(hActual(i,j,k), SEAICE_hice_reg) tmp = hActual(i,j,k) hrMin(i,j,k) = MIN(2.D0 * tmp, tmp+SEAICEmaxRaft) hrExp(i,j,k) = SEAICEmuRidging*SQRT(tmp) C arent we missing a factor 0.5 here? ridgeRatio(i,j,k)=(hrMin(i,j,k)+hrExp(i,j,k))/tmp ENDIF ENDDO ENDDO ENDDO ELSE STOP 'Ooops: SEAICEredistFunc > 1 not implemented' ENDIF C Compute the norm of the ridging mode N (in Lipscomp et al 2007) C or omega (in Bitz et al 2001): C rigdingModeNorm = net ice area removed / total area participating. C For instance, if a unit area of ice with thickness = 1 participates in C ridging to form a ridge with a = 1/3 and thickness = 3, then C rigdingModeNorm = 1 - 1/3 = 2/3. DO j=jMin,jMax DO i=iMin,iMax ridgingModeNorm(i,j) = partFunc(i,j,0) ENDDO ENDDO DO k = 1, nITD DO j=jMin,jMax DO i=iMin,iMax partFunc(i,j,k) = partFunc(i,j,k) * heffM(i,j,bi,bj) ridgingModeNorm(i,j) = ridgingModeNorm(i,j) & + partFunc(i,j,k)*( 1. _d 0 - 1. _d 0/ridgeRatio(i,j,k) ) ENDDO ENDDO ENDDO C avoid division by zero DO j=jMin,jMax DO i=iMin,iMax IF ( ridgingModeNorm(i,j) .LE. 0. _d 0 ) & ridgingModeNorm(i,j) = 1. _d 0 ENDDO ENDDO #endif /* SEAICE_ITD */ RETURN END