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"
                0034 #include "SEAICE.h"
                0035
                0036 C !INPUT PARAMETERS: ===================================================
                0037 C     === Routine arguments ===
                0038 C     bi, bj    :: outer loop counters
                0039 C     myTime    :: current time
                0040 C     myIter    :: iteration number
                0041 C     myThid    :: Thread no. that called this routine.
                0042       INTEGER bi,bj
73c2e960d7 Jean*0043       _RL myTime
0c32bd3cb0 Mart*0044       INTEGER myIter
                0045       INTEGER myThid
73c2e960d7 Jean*0046 CEOP
353a8877c7 Mart*0047
45315406aa Mart*0048 #if ( defined SEAICE_CGRID || defined SEAICE_BGRID_DYNAMICS )
353a8877c7 Mart*0049 C !LOCAL VARIABLES: ====================================================
                0050 C     === Local variables ===
                0051 C     i,j,k       :: inner loop counters
                0052 C     i/jMin/Max  :: loop boundaries
                0053 C
73c2e960d7 Jean*0054       INTEGER i, j
353a8877c7 Mart*0055       INTEGER iMin, iMax, jMin, jMax
                0056       _RL tmpscal1, tmpscal2
0c32bd3cb0 Mart*0057 #ifdef SEAICE_ITD
                0058 C     variables related to ridging schemes
                0059 C     ridgingModeNorm :: norm to ensure convervation (N in Lipscomb et al 2007)
                0060 C     partFunc   :: participation function (a_n in Lipscomb et al 2007)
                0061 C     ridgeRatio :: mean ridge thickness/ thickness of ridging ice
                0062 C     hrMin      :: min ridge thickness
                0063 C     hrMax      :: max ridge thickness   (SEAICEredistFunc = 0)
                0064 C     hrExp      :: ridge e-folding scale (SEAICEredistFunc = 1)
                0065 C     hActual    :: HEFFITD/AREAITD
73c2e960d7 Jean*0066       INTEGER k
0c32bd3cb0 Mart*0067       _RL ridgingModeNorm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0068       _RL partFunc        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,0:nITD)
                0069       _RL hrMin           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
                0070       _RL hrMax           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
                0071       _RL hrExp           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
                0072       _RL ridgeRatio      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
73c2e960d7 Jean*0073       _RL hActual         (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:nITD)
0c32bd3cb0 Mart*0074 #endif /* SEAICE_ITD */
5c768f1941 Mart*0075 #ifdef SEAICE_CGRID
                0076 C     compute tensile strength
e4120ef5ff Jean*0077 c     _RL recip_tensilDepth
5c768f1941 Mart*0078 #endif /* SEAICE_CGRID */
0c32bd3cb0 Mart*0079 CEOP
                0080
                0081 C     loop boundaries
73c2e960d7 Jean*0082       iMin=1-OLx
                0083       iMax=sNx+OLx
                0084       jMin=1-OLy
                0085       jMax=sNy+OLy
0c32bd3cb0 Mart*0086
                0087 #ifdef SEAICE_ITD
f3136e1434 Mart*0088 C     compute the fraction of open water as early as possible, i.e.
                0089 C     before advection, but also before it is used in calculating the ice
                0090 C     strength according to Rothrock (1975), hidden in S/R seaice_repare_ridging
                0091       DO j=jMin,jMax
                0092        DO i=iMin,iMax
                0093         opnWtrFrac(i,j,bi,bj) = 1. _d 0 - AREA(i,j,bi,bj)
                0094        ENDDO
                0095       ENDDO
                0096
0c32bd3cb0 Mart*0097       IF ( useHibler79IceStrength ) THEN
                0098 #else
                0099       IF ( .TRUE. ) THEN
                0100 #endif /* SEAICE_ITD */
                0101        DO j=jMin,jMax
                0102         DO i=iMin,iMax
                0103 C--   now set up ice pressure and viscosities
                0104          IF ( (HEFF(i,j,bi,bj).LE.SEAICEpresH0).AND.
                0105      &        (SEAICEpresPow0.NE.1) ) THEN
                0106           tmpscal1=MAX(HEFF(i,j,bi,bj)/SEAICEpresH0,ZERO)
                0107           tmpscal2=SEAICEpresH0*(tmpscal1**SEAICEpresPow0)
                0108          ELSEIF ( (HEFF(i,j,bi,bj).GT.SEAICEpresH0).AND.
                0109      &         (SEAICEpresPow1.NE.1) ) THEN
                0110           tmpscal1=MAX(HEFF(i,j,bi,bj)/SEAICEpresH0,ZERO)
                0111           tmpscal2=SEAICEpresH0*(tmpscal1**SEAICEpresPow1)
                0112          ELSE
73c2e960d7 Jean*0113           tmpscal2=HEFF(i,j,bi,bj)
0c32bd3cb0 Mart*0114          ENDIF
8e32c48b8f Mart*0115          PRESS0     (I,J,bi,bj) = SEAICE_strength*tmpscal2
ba6cfc5714 Mart*0116      &        *EXP(-SEAICE_cStar*(SEAICE_area_max-AREA(i,j,bi,bj)))
8e32c48b8f Mart*0117          SEAICE_zMax(I,J,bi,bj) = SEAICE_zetaMaxFac*PRESS0(I,J,bi,bj)
                0118          SEAICE_zMin(I,J,bi,bj) = SEAICE_zetaMin
                0119          PRESS0     (I,J,bi,bj) = PRESS0(I,J,bi,bj)*HEFFM(I,J,bi,bj)
0c32bd3cb0 Mart*0120         ENDDO
73c2e960d7 Jean*0121        ENDDO
0c32bd3cb0 Mart*0122 #ifdef SEAICE_ITD
                0123       ELSE
                0124 C     not useHiber79IceStrength
                0125        DO j=jMin,jMax
                0126         DO i=iMin,iMax
                0127          PRESS0(i,j,bi,bj) = 0. _d 0
                0128         ENDDO
                0129        ENDDO
                0130        CALL SEAICE_PREPARE_RIDGING(
353a8877c7 Mart*0131      O      hActual,
0c32bd3cb0 Mart*0132      O      hrMin, hrMax, hrExp, ridgeRatio, ridgingModeNorm, partFunc,
                0133      I      iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
                0134        IF ( SEAICEredistFunc .EQ. 0 ) THEN
cd7aede93e Mart*0135         tmpscal1 = 1. _d 0 / 3. _d 0
0c32bd3cb0 Mart*0136         DO k = 1, nITD
                0137          DO j=jMin,jMax
                0138           DO i=iMin,iMax
cd7aede93e Mart*0139 C     replace (hrMax**3-hrMin**3)/(hrMax-hrMin) by identical
                0140 C     hrMax**2+hrMin**2 + hrMax*hrMin to avoid division by potentially
                0141 C     small number
73c2e960d7 Jean*0142            IF ( partFunc(i,j,k) .GT. 0. _d 0 )
0c32bd3cb0 Mart*0143      &          PRESS0(i,j,bi,bj) = PRESS0(i,j,bi,bj)
73c2e960d7 Jean*0144      &          + partFunc(i,j,k) * ( - hActual(i,j,k)**2
e4120ef5ff Jean*0145      &          + ( hrMax(i,j,k)**2 + hrMin(i,j,k)**2
cd7aede93e Mart*0146      &          + hrMax(i,j,k)*hrMin(i,j,k) )*tmpscal1
0c32bd3cb0 Mart*0147      &          / ridgeRatio(i,j,k) )
                0148           ENDDO
                0149          ENDDO
                0150         ENDDO
                0151        ELSEIF ( SEAICEredistFunc .EQ. 1 ) THEN
                0152         DO k = 1, nITD
                0153          DO j=jMin,jMax
                0154           DO i=iMin,iMax
                0155            PRESS0(i,j,bi,bj) = PRESS0(i,j,bi,bj)
73c2e960d7 Jean*0156      &          + partFunc(i,j,k) * ( - hActual(i,j,k)**2 +
0c32bd3cb0 Mart*0157      &          (           hrMin(i,j,k)*hrMin(i,j,k)
                0158      &          + 2. _d 0 * hrMin(i,j,k)*hrExp(i,j,k)
                0159      &          + 2. _d 0 * hrExp(i,j,k)*hrExp(i,j,k)
                0160      &          )/ridgeRatio(i,j,k) )
                0161           ENDDO
                0162          ENDDO
                0163         ENDDO
                0164        ENDIF
73c2e960d7 Jean*0165 C
1d045ed8ea Mart*0166        tmpscal1 = SEAICE_cf*0.5*gravity*(rhoConst-SEAICE_rhoIce)
                0167      &      * SEAICE_rhoIce/rhoConst
0c32bd3cb0 Mart*0168        DO j=jMin,jMax
                0169         DO i=iMin,iMax
8e32c48b8f Mart*0170          PRESS0(i,j,bi,bj)      = PRESS0(i,j,bi,bj)/ridgingModeNorm(i,j)
0c32bd3cb0 Mart*0171      &        *tmpscal1
8e32c48b8f Mart*0172          SEAICE_zMax(I,J,bi,bj) = SEAICE_zetaMaxFac*PRESS0(I,J,bi,bj)
                0173          SEAICE_zMin(I,J,bi,bj) = SEAICE_zetaMin
                0174          PRESS0     (I,J,bi,bj) = PRESS0(I,J,bi,bj)*HEFFM(I,J,bi,bj)
0c32bd3cb0 Mart*0175         ENDDO
                0176        ENDDO
                0177 #endif /* SEAICE_ITD */
                0178       ENDIF
                0179
2f5e8addfd Mart*0180 CML#ifdef SEAICE_CGRID
                0181 CMLC     compute tensile strength factor k: tensileStrength = k*PRESS
e4120ef5ff Jean*0182 CMLC     can be done in initialisation phase as long as it depends only
2f5e8addfd Mart*0183 CMLC     on depth
                0184 CML      IF ( SEAICE_tensilFac .NE. 0. _d 0 ) THEN
                0185 CML       recip_tensilDepth = 0. _d 0
e4120ef5ff Jean*0186 CML       IF ( SEAICE_tensilDepth .GT. 0. _d 0 )
2f5e8addfd Mart*0187 CML     &      recip_tensilDepth = 1. _d 0 / SEAICE_tensilDepth
                0188 CML       DO j=jMin,jMax
                0189 CML        DO i=iMin,iMax
                0190 CML         tensileStrFac(i,j,bi,bj) = SEAICE_tensilFac
                0191 CML     &        *exp(-ABS(R_low(I,J,bi,bj))*recip_tensilDepth)
                0192 CML        ENDDO
                0193 CML       ENDDO
                0194 CML      ENDIF
                0195 CML#endif /* SEAICE_CGRID */
45315406aa Mart*0196 #endif /* SEAICE_CGRID or SEAICE_BGRID_DYNAMICS */
5c768f1941 Mart*0197
0c32bd3cb0 Mart*0198       RETURN
                0199       END