C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_do_ridging.F,v 1.8 2014/10/20 03:20:57 gforget Exp $
C $Name: $
#include "SEAICE_OPTIONS.h"
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif
CBOP
C !ROUTINE: SEAICE_DO_RIDGING
C !INTERFACE: ==========================================================
SUBROUTINE SEAICE_DO_RIDGING(
I bi, bj, myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *===========================================================*
C | SUBROUTINE SEAICE_DO_RIDGING
C | o compute mechanical redistribution of thin (undeformed) into
C | thick (deformed, i.e. ridged) ice categories
C | according to Thorndyke et al (1975) and Hibler (1980)
C | or Bitz et al (2001) and Lipscomb et al (2007)
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.
_RL myTime
INTEGER bi,bj
INTEGER myIter
INTEGER myThid
CEndOfInterface
C !LOCAL VARIABLES: ====================================================
C === Local variables ===
C i,j,k :: inner loop counters
C openWater :: open water area fraction
C
INTEGER i, j
INTEGER iMin, iMax, jMin, jMax
#ifdef SEAICE_ITD
INTEGER k, l, n
_RL ridgingModeNorm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL partFunc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,0:nITD)
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)
_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)
C computed and returned by S/R seaice_prepare_ridging, but not needed here
_RL hActual (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
C
_RL openWater (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL netArea (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C variables related to ridging schemes
_RL openingRate (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL closingRate (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL grossClosing (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C amount of ice that participates in ridging (according to partFunc)
_RL ridgingArea (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL ridgingHeff (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL ridgingHsnw (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C fractions of deformed/ridged ice
_RL areaFraction (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL volFraction (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C absolute area/concentration of deformed/ridged ice
_RL ridgedArea (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
LOGICAL doRidging (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
LOGICAL doRidgeAgain, areaTooLarge
C
_RL recip_deltaT, convergence, divergence, shear, divAdv
_RL tmp, tmpFac, hL, hR, expL, expR
_RL areaPR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
_RL heffPR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
_RL hsnwPR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
C
CHARACTER*(MAX_LEN_MBUF) msgBuf
#endif /* SEAICE_ITD */
CEOP
C---+-|--1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
iMin=1
iMax=sNx
jMin=1
jMax=sNy
#ifndef SEAICE_ITD
C Hiblers "ridging function" for single category ice
DO j=jMin,jMax
DO i=iMin,iMax
AREA(I,J,bi,bj) = MIN(AREA(I,J,bi,bj),SEAICE_area_max)
ENDDO
ENDDO
#else
C calculate area of open water
DO j=jMin,jMax
DO i=iMin,iMax
openWater(i,j) = ONE
ENDDO
ENDDO
DO k=1,nITD
DO j=jMin,jMax
DO i=iMin,iMax
openWater(i,j) = openWater(i,j) - AREAITD(i,j,k,bi,bj)
ENDDO
ENDDO
ENDDO
IF ( SEAICEsimpleRidging ) THEN
C-- Hibler-type "ridging", i.e. cut back excessive ice area fraction ---
C in case ice concentration exceeds 100% assume that
C convergence of floe field has eliminated all open water
C and eventual rafting occured in thinnest category:
DO j=jMin,jMax
DO i=iMin,iMax
IF (openWater(i,j) .lt. 0.0)
& AREAITD(i,j,1,bi,bj) = openWater(i,j)+AREAITD(i,j,1,bi,bj)
ENDDO
ENDDO
C---+-|--1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
ELSE
C initialisation
DO j=jMin,jMax
DO i=iMin,iMax
ridgingArea(i,j) = 0. _d 0
ridgingHeff(i,j) = 0. _d 0
ridgingHsnw(i,j) = 0. _d 0
areaFraction(i,j) = 0. _d 0
volFraction(i,j) = 0. _d 0
fw2ObyRidge(i,j,bi,bj) = 0. _d 0
ridgedArea(i,j) = 0. _d 0
ENDDO
ENDDO
CALL SEAICE_PREPARE_RIDGING(
O hActual,
O hrMin, hrMax, hrExp, ridgeRatio, ridgingModeNorm, partFunc,
I iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
C Compute the first strain rate invariant epsilonI (divergence)
C energy dissipation by convergence = -min (divergence, 0)
C energy dissipation by shearing = (1/2) * (Delta - abs(divergence))
DO j=jMin,jMax
DO i=iMin,iMax
divergence = e11(i,j,bi,bj) + e22(i,j,bi,bj)
shear = 0.5 _d 0 * ( deltaC(i,j,bi,bj) - ABS(divergence) )
convergence = - MIN(divergence, 0.D0)
closingRate(i,j) = SEAICEshearParm*shear + convergence
ENDDO
ENDDO
C we need a new estimate of the total AREA (including the open water
C fraction, but for computational reason it is not included here)
DO j=jMin,jMax
DO i=iMin,iMax
netArea(i,j) = 0. _d 0
ENDDO
ENDDO
DO k=1,nITD
DO j=jMin,jMax
DO i=iMin,iMax
netArea(i,j) = netArea(i,j) + AREAITD(i,j,k,bi,bj)
ENDDO
ENDDO
ENDDO
recip_DeltaT = 1. _d 0/SEAICE_deltaTtherm
DO j=jMin,jMax
DO i=iMin,iMax
C divergence rate due to advection; this term need not be zero due
C to numerical effects
C (this is copied from CICE but I am not sure about that)
divAdv = (1. _d 0-netArea(i,j)-opnWtrFrac(i,j,bi,bj))
& *recip_deltaT
IF (divAdv .LT. 0. _d 0)
& closingRate(i,j) = MAX(closingRate(i,j), -divAdv)
C finally compute a non-negative opening rate that will lead to
C a net area of 1
openingRate(i,j) = closingRate(i,j) + divAdv
ENDDO
ENDDO
C
C start of the ridging loop
C
doRidgeAgain = .TRUE.
n = 1
DO WHILE (doRidgeAgain)
C save pre-ridging ice concentration and ridged ice volume
DO k=1,nITD
DO j=jMin,jMax
DO i=iMin,iMax
areaPR(i,j,k) = AREAITD(i,j,k,bi,bj)
heffPR(i,j,k) = HEFFITD(i,j,k,bi,bj)
hsnwPR(i,j,k) = HSNOWITD(i,j,k,bi,bj)
C ridgeFrac(i,j,k) = 0. _d 0
C IF (HEFFITD(i,j,k,bi,bj) .GT. 0. _d 0 )
C & ridgeFrac(i,j,k) = ridgedHeff/HEFF
ENDDO
ENDDO
ENDDO
C
DO j=jMin,jMax
DO i=iMin,iMax
C Based on the ITD of ridging and ridged ice, convert the net
C closing rate to a gross closing rate times deltaT.
C NOTE: 0 < ridgingModeNorm <= 1
grossClosing(i,j) = closingRate(i,j)*SEAICE_deltaTtherm
& /ridgingModeNorm(i,j)
C reduce rates in case more than 100% of open water would be removed
IF ( partFunc(i,j,0) .GT. 0. _d 0 ) THEN
tmp = partFunc(i,j,0)*grossClosing(i,j)
IF ( tmp .GT. opnWtrFrac(i,j,bi,bj) ) THEN
tmpFac = opnWtrFrac(i,j,bi,bj)/tmp
grossClosing(i,j) = grossClosing(i,j) * tmpFac
openingRate(i,j) = openingRate(i,j) * tmpFac
ENDIF
ENDIF
ENDDO
ENDDO
DO k=1,nITD
DO j=jMin,jMax
DO i=iMin,iMax
C reduce rates in case more than 100% of any ice categroy would be removed
IF ( areaPR(i,j,k) .GT. SEAICE_area_reg
& .AND. partFunc(i,j,k) .GT. 0. _d 0 ) THEN
tmp = partFunc(i,j,k)*grossClosing(i,j)
IF ( tmp .GT. AREAITD(i,j,k,bi,bj) ) THEN
tmpFac = AREAITD(i,j,k,bi,bj)/tmp
grossClosing(i,j) = grossClosing(i,j) * tmpFac
openingRate(i,j) = openingRate(i,j) * tmpFac
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
C
C start redistribution
C
DO j=jMin,jMax
DO i=iMin,iMax
C open water first
opnWtrFrac(i,j,bi,bj) = opnWtrFrac(i,j,bi,bj)
& - partFunc(i,j,0)*grossClosing(i,j)
& + openingRate(i,j)*SEAICE_deltaTtherm
C need to catch openWater << 0 properly
C negative open water it not allowed
opnWtrFrac(i,j,bi,bj) = MAX( 0. _d 0, opnWtrFrac(i,j,bi,bj) )
ENDDO
ENDDO
C
DO k=1,nITD
C need to catch partFunc*grossClosing > AREAITD (or areaPR)
DO j=jMin,jMax
DO i=iMin,iMax
doRidging(i,j) = areaPR(i,j,k) .GT. SEAICE_area_reg
& .AND. partFunc(i,j,k) .GT. 0. _d 0
& .AND. grossClosing(i,j) .GT. 0. _d 0
& .AND. HEFFM(i,j,bi,bj) .GT. 0. _d 0
C this would be safety catch only
C & .AND. netArea(i,j) .GT. 1. _d 0
IF ( doRidging(i,j) ) THEN
CML ridgingArea(i,j) = MIN(partFunc(i,j,k)*grossClosing(i,j),
CML & areaPR(i,j,k))
ridgingArea(i,j) = partFunc(i,j,k)*grossClosing(i,j)
IF ( ridgingArea(i,j) .GT. areaPR(i,j,k) ) THEN
ridgingArea(i,j) = areaPR(i,j,k)
ENDIF
areaFraction(i,j) = ridgingArea(i,j)/areaPR(i,j,k)
ridgedArea(i,j) = ridgingArea(i,j)/ridgeRatio(i,j,k)
C compute ice volume (HEFF) and snow volume to be removed from this
C ridging category
ridgingHEFF(i,j) = heffPR(i,j,k) * areaFraction(i,j)
ridgingHsnw(i,j) = hsnwPR(i,j,k) * areaFraction(i,j)
C part of the snow mass is pushed into the ocean during ridging;
C this freshwater flux will be added to the net feshwater flux into
C the ocean in seaice_growth
fw2ObyRidge(i,j,bi,bj) = fw2ObyRidge(i,j,bi,bj)
& + SEAICE_rhoSnow*ridgingHsnw(i,j)
& *(1. _d 0 - SEAICEsnowFracRidge)
C reduce the snow volume that is left for redistribution
ridgingHsnw(i,j) = ridgingHsnw(i,j) * SEAICEsnowFracRidge
C remove ice concentration, volume (HEFF), and snow volume from
C this ridging category
AREAITD(i,j,k,bi,bj) = AREAITD(i,j,k,bi,bj) - ridgingArea(i,j)
HEFFITD(i,j,k,bi,bj) = HEFFITD(i,j,k,bi,bj) - ridgingHeff(i,j)
HSNOWITD(i,j,k,bi,bj)=HSNOWITD(i,j,k,bi,bj) - ridgingHsnw(i,j)
ENDIF
ENDDO
ENDDO
C inner loop over categories: distribute what has been removed from the
C kth category to all categories according to area/volFraction
DO l=1,nITD
C initialising these is essential, because here the ridging-mask doRidging
C to area/volFraction, and applied via these fields
DO j=jMin,jMax
DO i=iMin,iMax
areaFraction(i,j) = 0. _d 0
volFraction (i,j) = 0. _d 0
ENDDO
ENDDO
IF ( SEAICEredistFunc .EQ. 0 ) THEN
C Assume ridged ice is uniformly distributed between hrmin and hrmax
C (Hibler, 1980), see also s/r seaice_prepare_ridging.
DO j=jMin,jMax
DO i=iMin,iMax
IF ( doRidging(i,j) ) THEN
IF ( hrMin(i,j,k) .GE. hLimit(l) .OR.
& hrMax(i,j,k) .LE. hLimit(l-1) ) THEN
CML hL = 0. _d 0
CML hR = 0. _d 0
areaFraction(i,j) = 0. _d 0
volFraction (i,j) = 0. _d 0
ELSE
hL = MAX(hrMin(i,j,k), hLimit(l-1))
hR = MIN(hrMax(i,j,k), hLimit(l))
areaFraction(i,j) = ( hR - hL )
& / ( hrMax(i,j,k) - hrMin(i,j,k) )
CML volFraction (i,j) = ( hR*hR - hL*hL )
CML & / ( hrMax(i,j,k)**2 - hrMin(i,j,k)**2 )
volFraction (i,j) = areaFraction(i,j)*( hR + hL )
& / ( hrMax(i,j,k) + hrMin(i,j,k) )
ENDIF
ENDIF
ENDDO
ENDDO
ELSEIF ( SEAICEredistFunc .EQ. 1 ) THEN
C Follow Lipscomb et al. (2007) and model ridge ITD as an exponentially
C decaying function, see also s/r seaice_prepare_ridging.
IF ( l.LT.nITD ) THEN
DO j=jMin,jMax
DO i=iMin,iMax
IF ( doRidging(i,j)
& .AND. hrMin(i,j,k) .LT. hLimit(l)
& .AND. hrExp(i,j,k) .NE. 0. _d 0 ) THEN
hL = MAX( hrMin(i,j,k), hLimit(l-1) )
hR = hLimit(l)
expL = EXP(-( hL - hrMin(i,j,k) )/hrExp(i,j,k) )
expR = EXP(-( hR - hrMin(i,j,k) )/hrExp(i,j,k) )
areaFraction(i,j) = expL - expR
volFraction (i,j) =
& ( ( hL + hrExp(i,j,k) ) * expL
& - ( hR + hrExp(i,j,k) ) * expR )
& / ( hrMin(i,j,k) + hrExp(i,j,k) )
ENDIF
ENDDO
ENDDO
ELSE
DO j=jMin,jMax
DO i=iMin,iMax
IF ( doRidging(i,j) .AND. hrExp(i,j,k) .NE. 0. _d 0 ) THEN
hL = MAX( hrMin(i,j,k), hLimit(l-1) )
expL = EXP(-( hL - hrMin(i,j,k) )/hrExp(i,j,k) )
areaFraction(i,j) = expL
volFraction (i,j) = ( hL + hrExp(i,j,k) ) * expL
& / ( hrMin(i,j,k) + hrExp(i,j,k) )
ENDIF
ENDDO
ENDDO
ENDIF
ENDIF
C after computing the fraction ridged ice for this category, apply it
DO j=jMin,jMax
DO i=iMin,iMax
AREAITD(i,j,l,bi,bj) = AREAITD(i,j,l,bi,bj)
& +areaFraction(i,j)*ridgedArea(i,j)
HEFFITD(i,j,l,bi,bj) = HEFFITD(i,j,l,bi,bj)
& +volFraction(i,j)*ridgingHeff(i,j)
HSNOWITD(i,j,l,bi,bj) = HSNOWITD(i,j,l,bi,bj)
& +volFraction(i,j)*ridgingHsnw(i,j)*SEAICEsnowFracRidge
ENDDO
ENDDO
C category l-loop
ENDDO
C category k-loop
ENDDO
C determine if the ridging process needs to be repeated
C we need a new estimate of the total AREA
DO j=jMin,jMax
DO i=iMin,iMax
netArea(i,j) = 0. _d 0
ENDDO
ENDDO
DO k=1,nITD
DO j=jMin,jMax
DO i=iMin,iMax
netArea(i,j) = netArea(i,j) + AREAITD(i,j,k,bi,bj)
ENDDO
ENDDO
ENDDO
doRidgeAgain = .FALSE.
DO j=jMin,jMax
DO i=iMin,iMax
tmp = netArea(i,j)+opnWtrFrac(i,j,bi,bj)
areaTooLarge = tmp - 1. _d 0 .GT. 1. _d -11
IF ( HEFFM(i,j,bi,bj) .GT. 0. _d 0 .AND. areaTooLarge ) THEN
doRidging(i,j) = .TRUE.
doRidgeAgain = .TRUE.
divAdv = (1. _d 0-tmp)*recip_deltaT
closingRate(i,j) = MAX( 0. _d 0, -divAdv)
openingRate(i,j) = MAX( 0. _d 0, divAdv)
ELSE
C set to zero avoid going through this grid point again
closingRate(i,j) = 0. _d 0
openingRate(i,j) = 0. _d 0
doRidging(i,j) = .FALSE.
ENDIF
ENDDO
ENDDO
IF ( doRidgeAgain .AND. n.GE.SEAICEridgingIterMax ) THEN
C some debugging information
WRITE(msgBuf,'(A)') 'SEAICE_DO_RIDGING: *** WARNING ***'
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT, myThid )
WRITE(msgBuf,'(A,I2,A)') 'SEAICE_DO_RIDGING: '//
& 'did not converge in SEAICEridgingIterMax = ',
& SEAICEridgingIterMax, ' iterations.'
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT, myThid )
CML CALL PRINT_ERROR( msgBuf, myThid )
ENDIF
doRidgeAgain = doRidgeAgain .AND. n.LT.SEAICEridgingIterMax
IF ( doRidgeAgain .AND. debugLevel .GE. debLevA ) THEN
C some debugging information
WRITE(msgBuf,'(A,I2,A,I10)')
& 'SEAICE_DO_RIDGING: Repeat ridging after iteration ',
& n, ' in timestep ', myIter
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT, myThid )
ENDIF
IF ( doRidgeAgain ) CALL SEAICE_PREPARE_RIDGING(
O hActual,
O hrMin, hrMax, hrExp, ridgeRatio, ridgingModeNorm, partFunc,
I iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
n = n + 1
C ridging iteration
ENDDO
C .not. SEAICEsimpleRidging
ENDIF
#endif /* SEAICE_ITD */
C after ridging is complete, the critical variables need to be
C regularized and checked for consistency. This is done in a separate
C routine
CML CALL SEAICE_REGULARIZE( bi, bj, myTime, myIter, myThid )
RETURN
END