C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_reg_ridge.F,v 1.5 2014/10/20 03:20:58 gforget Exp $
C $Name: $
#include "SEAICE_OPTIONS.h"
#ifdef ALLOW_EXF
# include "EXF_OPTIONS.h"
#endif
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif
CBOP
C !ROUTINE: SEAICE_REG_RIDGE
C !INTERFACE:
SUBROUTINE SEAICE_REG_RIDGE( myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *=================================================================*
C | SUBROUTINE seaice_reg_ridge
C | o this routine has two purposes:
C | (1) clean up after advection (undershoots etc.);
C | after advection, the sea ice variables may have unphysical
C | values, e.g. < 0 or very thin ice, that are regularized
C | here.
C | (2) driver for ice ridging;
C | concentration as a special case may be > 1 in convergent
C | motion and a ridging algorithm redistributes the ice to
C | limit the concentration to 1.
C | o called after S/R seaice_advdiff
C *=================================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"
#include "SEAICE_TRACER.h"
#ifdef ALLOW_EXF
# include "EXF_FIELDS.h"
#endif
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif /* ALLOW_AUTODIFF_TAMC */
C !INPUT/OUTPUT PARAMETERS:
C === Routine arguments ===
C myTime :: Simulation time
C myIter :: Simulation timestep number
C myThid :: Thread no. that called this routine.
_RL myTime
INTEGER myIter, myThid
#ifdef ALLOW_SEAICE
C !LOCAL VARIABLES:
C === Local variables ===
C i,j,bi,bj :: Loop counters
INTEGER i, j, bi, bj
C number of surface interface layer
C IT :: ice thickness category index (ITD and SEAICE_multDim code)
INTEGER IT
C reciprocal of time step
_RL recip_deltaTtherm
C temporary variables available for the various computations
_RL tmpscal1, tmpscal2
#ifdef SEAICE_ITD
_RL tmpscal1itd(1:sNx,1:sNy), tmpscal2itd(1:sNx,1:sNy)
_RL tmpscal3itd(1:sNx,1:sNy)
C reciprocal number of ice classes nITD
_RL recip_nitd
#endif /* SEAICE_ITD */
#ifdef ALLOW_DIAGNOSTICS
C Helper variables for diagnostics
_RL DIAGarrayA (1:sNx,1:sNy)
#endif /* ALLOW_DIAGNOSTICS */
CEOP
C
C === Routine body ===
C
recip_deltaTtherm = ONE / SEAICE_deltaTtherm
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef ALLOW_AUTODIFF_TAMC
act1 = bi - myBxLo(myThid)
max1 = myBxHi(myThid) - myBxLo(myThid) + 1
act2 = bj - myByLo(myThid)
max2 = myByHi(myThid) - myByLo(myThid) + 1
act3 = myThid - 1
max3 = nTx*nTy
act4 = ikey_dynamics - 1
iicekey = (act1 + 1) + act2*max1
& + act3*max1*max2
& + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */
DO J=1-OLy,sNy+OLy
DO I=1-OLx,sNx+OLx
d_HEFFbyNEG(I,J,bi,bj) = 0.0 _d 0
d_HSNWbyNEG(I,J,bi,bj) = 0.0 _d 0
#ifdef EXF_SEAICE_FRACTION
d_AREAbyRLX(I,J,bi,bj) = 0.0 _d 0
d_HEFFbyRLX(I,J,bi,bj) = 0.0 _d 0
#endif /* EXF_SEAICE_FRACTION */
#ifdef SEAICE_VARIABLE_SALINITY
saltFluxAdjust(I,J,bi,bj) = 0.0 _d 0
#endif /* SEAICE_VARIABLE_SALINITY */
ENDDO
ENDDO
C =====================================================================
C ========== PART 1: treat pathological cases (post advdiff) ==========
C =====================================================================
#if (defined ALLOW_AUTODIFF_TAMC defined SEAICE_MODIFY_GROWTH_ADJ)
Cgf no dependency through pathological cases treatment
IF ( SEAICEadjMODE.EQ.0 ) THEN
#endif
#ifdef EXF_SEAICE_FRACTION
CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
C-- (0) relax sea ice concentration towards observation
IF ( SEAICE_tauAreaObsRelax .GT. zeroRL ) THEN
DO J=1,sNy
DO I=1,sNx
IF ( exf_iceFraction(I,J,bi,bj).GT.AREA(I,J,bi,bj) ) THEN
d_AREAbyRLX(i,j,bi,bj) =
& SEAICE_deltaTtherm/SEAICE_tauAreaObsRelax
& * (exf_iceFraction(I,J,bi,bj) - AREA(I,J,bi,bj))
ENDIF
IF ( exf_iceFraction(I,J,bi,bj).GT.zeroRS .AND.
& AREA(I,J,bi,bj).EQ.0. _d 0) THEN
C d_HEFFbyRLX(i,j,bi,bj) = 1. _d 1 * siEps * d_AREAbyRLX(i,j,bi,bj)
d_HEFFbyRLX(i,j,bi,bj) = 1. _d 1 * siEps
ENDIF
#ifdef SEAICE_ITD
AREAITD(I,J,1,bi,bj) = AREAITD(I,J,1,bi,bj)
& + d_AREAbyRLX(i,j,bi,bj)
HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj)
& + d_HEFFbyRLX(i,j,bi,bj)
#endif /* SEAICE_ITD */
AREA(I,J,bi,bj) = AREA(I,J,bi,bj) + d_AREAbyRLX(i,j,bi,bj)
HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + d_HEFFbyRLX(i,j,bi,bj)
ENDDO
ENDDO
ENDIF
#endif /* EXF_SEAICE_FRACTION */
C-- (1) treat the case of negative values:
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
tmpscal1=0. _d 0
tmpscal2=0. _d 0
tmpscal1=MAX(-HEFFITD(I,J,IT,bi,bj),0. _d 0)
HEFFITD(I,J,IT,bi,bj)=HEFFITD(I,J,IT,bi,bj)+tmpscal1
d_HEFFbyNEG(I,J,bi,bj)=d_HEFFbyNEG(I,J,bi,bj)+tmpscal1
tmpscal2=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)
HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal2
d_HSNWbyNEG(I,J,bi,bj)=d_HSNWbyNEG(I,J,bi,bj)+tmpscal2
AREAITD(I,J,IT,bi,bj)=MAX(AREAITD(I,J,IT,bi,bj),0. _d 0)
C AREA, HEFF, and HSNOW will be updated at end of PART 1
C by calling SEAICE_ITD_SUM
ENDDO
ENDDO
ENDDO
C update mean thicknesses HEFF and HSNOW and total ice
C concentration AREA to match single category values
CALL SEAICE_ITD_SUM ( bi, bj, myTime, myIter, myThid )
#else /* ndef SEAICE_ITD */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
DO J=1,sNy
DO I=1,sNx
d_HEFFbyNEG(I,J,bi,bj)=MAX(-HEFF(I,J,bi,bj),0. _d 0)
HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J,bi,bj)
d_HSNWbyNEG(I,J,bi,bj)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)
HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J,bi,bj)
AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)
ENDDO
ENDDO
#endif /* SEAICE_ITD */
C-- (2) treat the case of very thin ice:
#ifdef SEAICE_ITD
C Here we risk that even though HEFF may be larger than siEps (=1e-5)
C HEFFITD can have classes with very small (< siEps) non-zero ice volume.
C We avoid applying the correction to each class because that leads to
C funny structures in the net heat and freshwater flux into the ocean.
C Let us keep our fingers crossed, that the model will be benign!
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
IF (HEFF(I,J,bi,bj).LE.siEps) THEN
HEFFITD(I,J,IT,bi,bj) = 0. _d 0
HSNOWITD(I,J,IT,bi,bj) = 0. _d 0
ENDIF
ENDDO
ENDDO
ENDDO
#endif /* SEAICE_ITD */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
DO J=1,sNy
DO I=1,sNx
tmpscal1=0. _d 0
tmpscal2=0. _d 0
IF (HEFF(I,J,bi,bj).LE.siEps) THEN
tmpscal1=-HEFF(I,J,bi,bj)
tmpscal2=-HSNOW(I,J,bi,bj)
DO IT=1,SEAICE_multDim
TICES(I,J,IT,bi,bj)=celsius2K
ENDDO
ENDIF
HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+tmpscal1
HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal2
d_HEFFbyNEG(I,J,bi,bj)=d_HEFFbyNEG(I,J,bi,bj)+tmpscal1
d_HSNWbyNEG(I,J,bi,bj)=d_HSNWbyNEG(I,J,bi,bj)+tmpscal2
ENDDO
ENDDO
C-- (3) treat the case of area but no ice/snow:
#ifdef SEAICE_ITD
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
IF ( (HEFFITD(I,J,IT,bi,bj) .EQ.0. _d 0).AND.
& (HSNOWITD(I,J,IT,bi,bj).EQ.0. _d 0))
& AREAITD(I,J,IT,bi,bj)=0. _d 0
ENDDO
ENDDO
ENDDO
#else /* ndef SEAICE_ITD */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
DO J=1,sNy
DO I=1,sNx
IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
& (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
ENDDO
ENDDO
#endif /* SEAICE_ITD */
C-- (4) treat the case of very small area:
#ifndef DISABLE_AREA_FLOOR
#ifdef SEAICE_ITD
recip_nitd = 1. _d 0 / float(SEAICE_multDim)
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
IF ((HEFFITD(I,J,IT,bi,bj).GT.0).OR.
& (HSNOWITD(I,J,IT,bi,bj).GT.0)) THEN
C SEAICE_area_floor*SEAICE_multDim cannot be allowed to exceed 1
C hence use SEAICE_area_floor devided by SEAICE_multDim
C (or install a warning in e.g. seaice_readparms.F)
AREAITD(I,J,IT,bi,bj)=
& MAX(AREAITD(I,J,IT,bi,bj),SEAICE_area_floor*recip_nitd)
ENDIF
ENDDO
ENDDO
ENDDO
#else /* ndef SEAICE_ITD */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
DO J=1,sNy
DO I=1,sNx
IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN
AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),SEAICE_area_floor)
ENDIF
ENDDO
ENDDO
#endif /* SEAICE_ITD */
#endif /* DISABLE_AREA_FLOOR */
C (5) treat sea ice salinity pathological cases
#ifdef SEAICE_VARIABLE_SALINITY
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
DO J=1,sNy
DO I=1,sNx
IF ( (HSALT(I,J,bi,bj) .LT. 0.0).OR.
& (HEFF(I,J,bi,bj) .EQ. 0.0) ) THEN
saltFluxAdjust(I,J,bi,bj) = - HEFFM(I,J,bi,bj) *
& HSALT(I,J,bi,bj) * recip_deltaTtherm
HSALT(I,J,bi,bj) = 0.0 _d 0
ENDIF
ENDDO
ENDDO
#endif /* SEAICE_VARIABLE_SALINITY */
C =====================================================================
C ========== PART 2: ridging algorithm ===============================
C =====================================================================
C treat case of excessive ice cover, e.g., due to ridging:
#ifdef SEAICE_ITD
C catch up with item (2) that involves category sums AREA and HEFF
DO J=1,sNy
DO I=1,sNx
tmpscal1itd(i,j) = 0. _d 0
tmpscal2itd(i,j) = 0. _d 0
tmpscal3itd(i,j) = 0. _d 0
ENDDO
ENDDO
DO IT=1,SEAICE_multDim
DO J=1,sNy
DO I=1,sNx
C TICES was changed above (item 2), now update TICE as ice volume
C weighted average of TICES
tmpscal1itd(i,j)=tmpscal1itd(i,j)
& + TICES(I,J,IT,bi,bj) * HEFFITD(I,J,IT,bi,bj)
tmpscal2itd(i,j)=tmpscal2itd(i,j) + HEFFITD(I,J,IT,bi,bj)
C also compute total of AREAITD for diagnostics and SItrArea
tmpscal3itd(i,j)=tmpscal3itd(i,j) + AREAITD(I,J,IT,bi,bj)
ENDDO
ENDDO
ENDDO
DO J=1,sNy
DO I=1,sNx
C save pre-ridging ice concentration for diagnostics:
C these lines are executed before "ridging" is applied to AREA
C hence we execute them here before SEAICE_ITD_REDIST is called
C although this means that AREA has not been completely regularized
#ifdef ALLOW_DIAGNOSTICS
DIAGarrayA(I,J) = tmpscal3itd(i,j)
#endif
#ifdef ALLOW_SITRACER
SItrAREA(I,J,bi,bj,1)=tmpscal3itd(i,j)
#endif
ENDDO
ENDDO
C ridge ice according to Lipscomb et al. (2007), Bitz et al. (2001)
C Thorndyke et al. (1975), Hibler (1980)
CALL SEAICE_DO_RIDGING( bi, bj, myTime, myIter, myThid )
C check that all ice thickness categories meet their limits
C (includes Hibler-type ridging)
CALL SEAICE_ITD_REDIST( bi, bj, myTime, myIter, myThid )
C update mean thicknesses HEFF and HSNOW and total ice
C concentration AREA to match single category values
CALL SEAICE_ITD_SUM ( bi, bj, myTime, myIter, myThid )
#else /* ifndef SEAICE_ITD */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
DO J=1,sNy
DO I=1,sNx
C save pre-ridging ice concentration for diagnostics
#ifdef ALLOW_DIAGNOSTICS
DIAGarrayA(I,J) = AREA(I,J,bi,bj)
#endif /* ALLOW_DIAGNOSTICS */
#ifdef ALLOW_SITRACER
SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)
#endif /* ALLOW_SITRACER */
C this is the simple Hibler (1979)-type ridging (capping of
C concentrations > 1) for the non-ITD sea ice model
AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),SEAICE_area_max)
ENDDO
ENDDO
#endif /* SEAICE_ITD */
#if (defined ALLOW_AUTODIFF_TAMC defined SEAICE_MODIFY_GROWTH_ADJ)
C end SEAICEadjMODE.EQ.0 statement:
ENDIF
#endif
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIareaPR',0,1,3,bi,bj,myThid)
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
C close bi,bj loops
ENDDO
ENDDO
#endif /* ALLOW_SEAICE */
RETURN
END