pkg/seaice/seaice_calc_ice_strength.F

0c32bd3cb0 Mart*0001 #include "SEAICE_OPTIONS.h"
772b2ed80e Gael*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
0c32bd3cb0 Mart*0005
                0006 CBOP
                0007 C !ROUTINE: SEAICE_CALC_ICE_STRENGTH
                0008 C !INTERFACE: ==========================================================
                0009       SUBROUTINE SEAICE_CALC_ICE_STRENGTH(
                0010      I     bi, bj, myTime, myIter, myThid )
                0011
                0012 C !DESCRIPTION: \bv
                0013 C     *===========================================================*
                0014 C     | SUBROUTINE SEAICE_CALC_ICE_STRENGTH
                0015 C     | o compute ice strengh PRESS0
73c2e960d7 Jean*0016 C     |   according to
0c32bd3cb0 Mart*0017 C     |   (a) Hibler (1979)
                0018 C     |   (b) Thorndyke et al (1975) and Hibler (1980)
                0019 C     |   (c) Bitz et al (2001) and Lipscomb et al (2007)
                0020 C     |
                0021 C     | Martin Losch, Apr. 2014, Martin.Losch@awi.de
                0022 C     *===========================================================*
                0023 C \ev
                0024
                0025 C !USES: ===============================================================
                0026       IMPLICIT NONE
                0027
                0028 #include "SIZE.h"
                0029 #include "EEPARAMS.h"
                0030 #include "PARAMS.h"
                0031 #include "GRID.h"
                0032 #include "SEAICE_SIZE.h"
                0033 #include "SEAICE_PARAMS.h"
3f0f10fc37 Mart*0034 #include "SEAICE_GRID.h"
0c32bd3cb0 Mart*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       INTEGER bi,bj
73c2e960d7 Jean*0044       _RL myTime
0c32bd3cb0 Mart*0045       INTEGER myIter
                0046       INTEGER myThid
73c2e960d7 Jean*0047 CEOP
353a8877c7 Mart*0048
45315406aa Mart*0049 #if ( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS )
353a8877c7 Mart*0050 C !LOCAL VARIABLES: ====================================================
                0051 C     === Local variables ===
                0052 C     i,j,k       :: inner loop counters
                0053 C     i/jMin/Max  :: loop boundaries
                0054 C
73c2e960d7 Jean*0055       INTEGER i, j
353a8877c7 Mart*0056       INTEGER iMin, iMax, jMin, jMax
                0057       _RL tmpscal1, tmpscal2
0c32bd3cb0 Mart*0058 #ifdef SEAICE_ITD
                0059 C     variables related to ridging schemes
                0060 C     ridgingModeNorm :: norm to ensure convervation (N in Lipscomb et al 2007)
                0061 C     partFunc   :: participation function (a_n in Lipscomb et al 2007)
                0062 C     ridgeRatio :: mean ridge thickness/ thickness of ridging ice
                0063 C     hrMin      :: min ridge thickness
                0064 C     hrMax      :: max ridge thickness   (SEAICEredistFunc = 0)
                0065 C     hrExp      :: ridge e-folding scale (SEAICEredistFunc = 1)
                0066 C     hActual    :: HEFFITD/AREAITD
73c2e960d7 Jean*0067       INTEGER k
0c32bd3cb0 Mart*0068       _RL ridgingModeNorm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0069       _RL partFunc        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,0:nITD)
                0070       _RL hrMin           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
                0071       _RL hrMax           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
                0072       _RL hrExp           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
                0073       _RL ridgeRatio      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
73c2e960d7 Jean*0074       _RL hActual         (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
0c32bd3cb0 Mart*0075 #endif /* SEAICE_ITD */
5c768f1941 Mart*0076 #ifdef SEAICE_CGRID
                0077 C     compute tensile strength
e4120ef5ff Jean*0078 c     _RL recip_tensilDepth
5c768f1941 Mart*0079 #endif /* SEAICE_CGRID */
0c32bd3cb0 Mart*0080 CEOP
                0081
                0082 C     loop boundaries
73c2e960d7 Jean*0083       iMin=1-OLx
                0084       iMax=sNx+OLx
                0085       jMin=1-OLy
                0086       jMax=sNy+OLy
0c32bd3cb0 Mart*0087
                0088 #ifdef SEAICE_ITD
f3136e1434 Mart*0089 C     compute the fraction of open water as early as possible, i.e.
                0090 C     before advection, but also before it is used in calculating the ice
                0091 C     strength according to Rothrock (1975), hidden in S/R seaice_repare_ridging
                0092       DO j=jMin,jMax
                0093        DO i=iMin,iMax
                0094         opnWtrFrac(i,j,bi,bj) = 1. _d 0 - AREA(i,j,bi,bj)
                0095        ENDDO
                0096       ENDDO
                0097
0c32bd3cb0 Mart*0098       IF ( useHibler79IceStrength ) THEN
                0099 #else
                0100       IF ( .TRUE. ) THEN
                0101 #endif /* SEAICE_ITD */
                0102        DO j=jMin,jMax
                0103         DO i=iMin,iMax
                0104 C--   now set up ice pressure and viscosities
                0105          IF ( (HEFF(i,j,bi,bj).LE.SEAICEpresH0).AND.
                0106      &        (SEAICEpresPow0.NE.1) ) THEN
                0107           tmpscal1=MAX(HEFF(i,j,bi,bj)/SEAICEpresH0,ZERO)
                0108           tmpscal2=SEAICEpresH0*(tmpscal1**SEAICEpresPow0)
                0109          ELSEIF ( (HEFF(i,j,bi,bj).GT.SEAICEpresH0).AND.
                0110      &         (SEAICEpresPow1.NE.1) ) THEN
                0111           tmpscal1=MAX(HEFF(i,j,bi,bj)/SEAICEpresH0,ZERO)
                0112           tmpscal2=SEAICEpresH0*(tmpscal1**SEAICEpresPow1)
                0113          ELSE
73c2e960d7 Jean*0114           tmpscal2=HEFF(i,j,bi,bj)
0c32bd3cb0 Mart*0115          ENDIF
3f0f10fc37 Mart*0116          PRESS0     (i,j,bi,bj) = SEAICE_strength*tmpscal2
ba6cfc5714 Mart*0117      &        *EXP(-SEAICE_cStar*(SEAICE_area_max-AREA(i,j,bi,bj)))
3f0f10fc37 Mart*0118          SEAICE_zMax(i,j,bi,bj) = SEAICE_zetaMaxFac*PRESS0(i,j,bi,bj)
                0119          SEAICE_zMin(i,j,bi,bj) = SEAICE_zetaMin
                0120          PRESS0     (i,j,bi,bj) = PRESS0(i,j,bi,bj)*HEFFM(i,j,bi,bj)
0c32bd3cb0 Mart*0121         ENDDO
73c2e960d7 Jean*0122        ENDDO
0c32bd3cb0 Mart*0123 #ifdef SEAICE_ITD
                0124       ELSE
                0125 C     not useHiber79IceStrength
                0126        DO j=jMin,jMax
                0127         DO i=iMin,iMax
                0128          PRESS0(i,j,bi,bj) = 0. _d 0
                0129         ENDDO
                0130        ENDDO
                0131        CALL SEAICE_PREPARE_RIDGING(
353a8877c7 Mart*0132      O      hActual,
0c32bd3cb0 Mart*0133      O      hrMin, hrMax, hrExp, ridgeRatio, ridgingModeNorm, partFunc,
                0134      I      iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
                0135        IF ( SEAICEredistFunc .EQ. 0 ) THEN
cd7aede93e Mart*0136         tmpscal1 = 1. _d 0 / 3. _d 0
0c32bd3cb0 Mart*0137         DO k = 1, nITD
                0138          DO j=jMin,jMax
                0139           DO i=iMin,iMax
cd7aede93e Mart*0140 C     replace (hrMax**3-hrMin**3)/(hrMax-hrMin) by identical
                0141 C     hrMax**2+hrMin**2 + hrMax*hrMin to avoid division by potentially
                0142 C     small number
73c2e960d7 Jean*0143            IF ( partFunc(i,j,k) .GT. 0. _d 0 )
0c32bd3cb0 Mart*0144      &          PRESS0(i,j,bi,bj) = PRESS0(i,j,bi,bj)
73c2e960d7 Jean*0145      &          + partFunc(i,j,k) * ( - hActual(i,j,k)**2
e4120ef5ff Jean*0146      &          + ( hrMax(i,j,k)**2 + hrMin(i,j,k)**2
cd7aede93e Mart*0147      &          + hrMax(i,j,k)*hrMin(i,j,k) )*tmpscal1
0c32bd3cb0 Mart*0148      &          / ridgeRatio(i,j,k) )
                0149           ENDDO
                0150          ENDDO
                0151         ENDDO
                0152        ELSEIF ( SEAICEredistFunc .EQ. 1 ) THEN
                0153         DO k = 1, nITD
                0154          DO j=jMin,jMax
                0155           DO i=iMin,iMax
                0156            PRESS0(i,j,bi,bj) = PRESS0(i,j,bi,bj)
73c2e960d7 Jean*0157      &          + partFunc(i,j,k) * ( - hActual(i,j,k)**2 +
0c32bd3cb0 Mart*0158      &          (           hrMin(i,j,k)*hrMin(i,j,k)
                0159      &          + 2. _d 0 * hrMin(i,j,k)*hrExp(i,j,k)
                0160      &          + 2. _d 0 * hrExp(i,j,k)*hrExp(i,j,k)
                0161      &          )/ridgeRatio(i,j,k) )
                0162           ENDDO
                0163          ENDDO
                0164         ENDDO
                0165        ENDIF
73c2e960d7 Jean*0166 C
1d045ed8ea Mart*0167        tmpscal1 = SEAICE_cf*0.5*gravity*(rhoConst-SEAICE_rhoIce)
                0168      &      * SEAICE_rhoIce/rhoConst
0c32bd3cb0 Mart*0169        DO j=jMin,jMax
                0170         DO i=iMin,iMax
8e32c48b8f Mart*0171          PRESS0(i,j,bi,bj)      = PRESS0(i,j,bi,bj)/ridgingModeNorm(i,j)
0c32bd3cb0 Mart*0172      &        *tmpscal1
3f0f10fc37 Mart*0173          SEAICE_zMax(i,j,bi,bj) = SEAICE_zetaMaxFac*PRESS0(i,j,bi,bj)
                0174          SEAICE_zMin(i,j,bi,bj) = SEAICE_zetaMin
                0175          PRESS0     (i,j,bi,bj) = PRESS0(i,j,bi,bj)*HEFFM(i,j,bi,bj)
0c32bd3cb0 Mart*0176         ENDDO
                0177        ENDDO
                0178 #endif /* SEAICE_ITD */
                0179       ENDIF
                0180
2f5e8addfd Mart*0181 CML#ifdef SEAICE_CGRID
                0182 CMLC     compute tensile strength factor k: tensileStrength = k*PRESS
e4120ef5ff Jean*0183 CMLC     can be done in initialisation phase as long as it depends only
2f5e8addfd Mart*0184 CMLC     on depth
                0185 CML      IF ( SEAICE_tensilFac .NE. 0. _d 0 ) THEN
                0186 CML       recip_tensilDepth = 0. _d 0
e4120ef5ff Jean*0187 CML       IF ( SEAICE_tensilDepth .GT. 0. _d 0 )
2f5e8addfd Mart*0188 CML     &      recip_tensilDepth = 1. _d 0 / SEAICE_tensilDepth
                0189 CML       DO j=jMin,jMax
                0190 CML        DO i=iMin,iMax
                0191 CML         tensileStrFac(i,j,bi,bj) = SEAICE_tensilFac
                0192 CML     &        *exp(-ABS(R_low(I,J,bi,bj))*recip_tensilDepth)
                0193 CML        ENDDO
                0194 CML       ENDDO
                0195 CML      ENDIF
                0196 CML#endif /* SEAICE_CGRID */
45315406aa Mart*0197 #endif /* SEAICE_CGRID or SEAICE_BGRID_DYNAMICS */
5c768f1941 Mart*0198
0c32bd3cb0 Mart*0199       RETURN
                0200       END