pkg/seaice/seaice_calc_viscosities.F

3ed8349f04 Mart*0001 #include "SEAICE_OPTIONS.h"
772b2ed80e Gael*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
3ed8349f04 Mart*0005
ad2b15ed4c Mart*0006 CBOP
3ed8349f04 Mart*0007 CStartOfInterface
b4949dd6db Jean*0008       SUBROUTINE SEAICE_CALC_VISCOSITIES(
ec0d7df165 Mart*0009      I     e11, e22, e12, zMin, zMax, HEFFM, press0, tnsFac,
0fd6b2de1a Mart*0010      O     eta, etaZ, zeta, zetaZ, press, deltaC,
b4949dd6db Jean*0011      I     iStep, myTime, myIter, myThid )
e826e89d88 Jean*0012 C     *==========================================================*
3ed8349f04 Mart*0013 C     | SUBROUTINE  SEAICE_CALC_VISCOSITIES                      |
                0014 C     | o compute shear and bulk viscositites eta, zeta and the  |
037351a1a6 Mart*0015 C     |   corrected ice strength P                               |
3ed8349f04 Mart*0016 C     |   (see Zhang and Hibler,   JGR, 102, 8691-8702, 1997)    |
e826e89d88 Jean*0017 C     *==========================================================*
c512e371cc drin*0018 C     | started by Martin Losch, Mar 2006                        |
e826e89d88 Jean*0019 C     *==========================================================*
3ed8349f04 Mart*0020       IMPLICIT NONE
                0021
                0022 C     === Global variables ===
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
                0026 #include "GRID.h"
c8d9ddcff2 Patr*0027 #include "SEAICE_SIZE.h"
3ed8349f04 Mart*0028 #include "SEAICE_PARAMS.h"
                0029
                0030 C     === Routine arguments ===
b4949dd6db Jean*0031 C     iStep  :: Sub-time-step number
                0032 C     myTime :: Simulation time
                0033 C     myIter :: Simulation timestep number
                0034 C     myThid :: My Thread Id. number
                0035       INTEGER iStep
                0036       _RL     myTime
                0037       INTEGER myIter
3ed8349f04 Mart*0038       INTEGER myThid
037351a1a6 Mart*0039 C     strain rate tensor
                0040       _RL e11   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0041       _RL e22   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0042       _RL e12   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
b4949dd6db Jean*0043 C
e826e89d88 Jean*0044       _RL zMin  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0045       _RL zMax  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
ec0d7df165 Mart*0046       _RL HEFFM (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
b4949dd6db Jean*0047 C
e826e89d88 Jean*0048       _RL press0(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0049       _RL press (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
36fa289698 Mart*0050 C
                0051       _RL deltaC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
2f5e8addfd Mart*0052 C     factor k to compute tensile strength from k*press0
                0053       _RL tnsFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
037351a1a6 Mart*0054 C     bulk viscosity
e826e89d88 Jean*0055       _RL  eta  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
ad2b15ed4c Mart*0056       _RL  etaZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
037351a1a6 Mart*0057 C     shear viscosity
e826e89d88 Jean*0058       _RL zeta  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0fd6b2de1a Mart*0059       _RL zetaZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
3ed8349f04 Mart*0060 CEndOfInterface
                0061
45315406aa Mart*0062 #ifdef SEAICE_CGRID
3ed8349f04 Mart*0063 C     === Local variables ===
                0064 C     i,j,bi,bj - Loop counters
                0065 C     e11, e12, e22 - components of strain rate tensor
314f7b043c Mart*0066 C     recip_e2      - inverse of square of ratio of yield curve main axes
                0067 C     ep            - e11+e22 (abbreviation)
                0068 C     em            - e11-e22 (abbreviation)
3ed8349f04 Mart*0069       INTEGER i, j, bi, bj
c512e371cc drin*0070       _RL recip_e2, deltaCsq, deltaMinSq, tmp, ep, em, smallNbr
                0071       _RL recip_efr4, recip_efr2, oneThird, oneNinth, smallNbrSq
                0072       _RL e12Csq     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0073       _RL deltaCreg  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0074       _RL recip_shear(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0075       _RL shearDefSq, shearDef
                0076 #ifdef SEAICE_ALLOW_TEARDROP
                0077       _RL k,ksq,x,a,ap1,etamax_TD,zetamax_TD,cyc,xma, etapr
                0078 #endif /* SEAICE_ALLOW_TEARDROP */
                0079 #if (defined SEAICE_ALLOW_MCE || defined SEAICE_ALLOW_TEM)
                0080       _RL etaMaxFac
                0081 #endif /* SEAICE_ALLOW_MCE or SEAICE_ALLOW_TEM */
0fd6b2de1a Mart*0082       _RL sumNorm, maskZ
6a337ce369 Mart*0083 #ifdef SEAICE_ZETA_SMOOTHREG
                0084       _RL argTmp
                0085 #endif /* SEAICE_ZETA_SMOOTHREG */
e826e89d88 Jean*0086 CEOP
3ed8349f04 Mart*0087
ec0d7df165 Mart*0088 C--   basic constants
c512e371cc drin*0089       oneThird = 1. _d 0 / 3. _d 0
                0090       oneNinth = 1. _d 0 / 9. _d 0
                0091       smallNbr = 1. _d -20
                0092       smallNbrSq = smallNbr * smallNbr
                0093
314f7b043c Mart*0094       recip_e2=0. _d 0
                0095       IF ( SEAICE_eccen .NE. 0. _d 0 ) recip_e2=ONE/(SEAICE_eccen**2)
c512e371cc drin*0096       recip_efr2=0. _d 0
                0097       recip_efr4=0. _d 0
                0098       IF( SEAICE_eccfr .NE. 0. _d 0) THEN
                0099        recip_efr2=ONE/(SEAICE_eccfr**2)
                0100        recip_efr4=(SEAICE_eccen**2) / (SEAICE_eccfr**4)
                0101       ENDIF
3daf25222c Mart*0102       deltaMinSq = SEAICE_deltaMin**2
3ed8349f04 Mart*0103 C
                0104       DO bj=myByLo(myThid),myByHi(myThid)
                0105        DO bi=myBxLo(myThid),myBxHi(myThid)
8377b8ee87 Mart*0106 #ifdef ALLOW_AUTODIFF
e7b5767f0b Mart*0107         DO j=1-OLy,sNy+OLy
                0108          DO i=1-OLx,sNx+OLx
c512e371cc drin*0109           deltaCreg  (i,j) = SEAICE_deltaMin
                0110           e12Csq     (i,j) = 0. _d 0
                0111           recip_shear(i,j) = 0. _d 0
e7b5767f0b Mart*0112          ENDDO
                0113         ENDDO
8377b8ee87 Mart*0114 #endif /* ALLOW_AUTODIFF */
03cd49feda Mart*0115 C     need to do this beforehand for easier vectorization after
                0116 C     TAFization
dc66f404d9 Mart*0117         IF ( SEAICEetaZmethod .EQ. 0 ) THEN
e826e89d88 Jean*0118          DO j=1-OLy+1,sNy+OLy-1
                0119           DO i=1-OLx+1,sNx+OLx-1
ad2b15ed4c Mart*0120            tmp = 0.25 *
c512e371cc drin*0121      &          ( e12(i,j  ,bi,bj) + e12(i+1,j  ,bi,bj)
                0122      &          + e12(i,j+1,bi,bj) + e12(i+1,j+1,bi,bj) )
f3224168cd Mart*0123            e12Csq(i,j) = tmp*tmp
ad2b15ed4c Mart*0124           ENDDO
03cd49feda Mart*0125          ENDDO
dc66f404d9 Mart*0126         ELSEIF ( SEAICEetaZmethod .EQ. 3 ) THEN
e826e89d88 Jean*0127          DO j=1-OLy+1,sNy+OLy-1
                0128           DO i=1-OLx+1,sNx+OLx-1
dc66f404d9 Mart*0129 C     area weighted average of the squares of e12 is more accurate
                0130 C     (and energy conserving)
c512e371cc drin*0131            e12Csq(i,j) = 0.25 _d 0 * recip_rA(i,j,bi,bj) *
                0132      &          ( rAz(i  ,j  ,bi,bj)*e12(i  ,j  ,bi,bj)**2
                0133      &          + rAz(i+1,j  ,bi,bj)*e12(i+1,j  ,bi,bj)**2
                0134      &          + rAz(i  ,j+1,bi,bj)*e12(i  ,j+1,bi,bj)**2
                0135      &          + rAz(i+1,j+1,bi,bj)*e12(i+1,j+1,bi,bj)**2 )
ad2b15ed4c Mart*0136           ENDDO
                0137          ENDDO
                0138         ENDIF
e826e89d88 Jean*0139         DO j=1-OLy+1,sNy+OLy-1
                0140          DO i=1-OLx+1,sNx+OLx-1
314f7b043c Mart*0141           em = e11(i,j,bi,bj)-e22(i,j,bi,bj)
c512e371cc drin*0142           shearDefSq = em*em + 4. _d 0 * e12Csq(i,j)
                0143 #ifdef SEAICE_DELTA_SMOOTHREG
                0144           recip_shear(i,j) = 1. _d 0
                0145      &                / SQRT( shearDefSq + smallNbrSq )
                0146 C         recip_shear(i,j) = 1. _d 0 / ( shearDef + smallNbr )
                0147 #else
8377b8ee87 Mart*0148 # ifdef ALLOW_AUTODIFF
c512e371cc drin*0149 C     avoid sqrt of 0
                0150           shearDef = 0. _d 0
                0151           IF ( shearDefSq .GT. 0. _d 0 ) shearDef = SQRT(shearDefSq)
                0152 # else
                0153           shearDef = SQRT(shearDefSq)
8377b8ee87 Mart*0154 # endif /* ALLOW_AUTODIFF */
c512e371cc drin*0155           recip_shear(i,j) = 1. _d 0 / MAX( ShearDef, smallNbr )
                0156 #endif /* SEAICE_DELTA_SMOOTHREG */
                0157          ENDDO
                0158         ENDDO
                0159
                0160 C     select between different yield curves; default: elliptical yield curve
                0161         IF ( .FALSE. ) THEN
                0162 C     do nothing (hack to be able add code with different cpp-flag)
                0163 #ifdef SEAICE_ALLOW_TEARDROP
                0164         ELSEIF ( SEAICEuseTD ) THEN
                0165
                0166          DO j=1-OLy+1,sNy+OLy-1
                0167           DO i=1-OLx+1,sNx+OLx-1
                0168            ep = e11(i,j,bi,bj)+e22(i,j,bi,bj)
                0169            em = e11(i,j,bi,bj)-e22(i,j,bi,bj)
                0170 C            shearDef = SQRT(em*em + 4. _d 0 * e12Csq(i,j))
                0171
                0172            k   = ep * recip_shear(i,j)
                0173            ksq = k * k
                0174            a   = tnsFac(i,j,bi,bj)
                0175
                0176 C     cyc: coordinate of the maximum point of the yield curve,
                0177 C     where eI changes sign
                0178            cyc=( TWO - tnsFac(i,j,bi,bj) ) * oneThird
                0179            ap1 = a + ONE
                0180
                0181 C     handle floating point errors when k is negative and very large
                0182            IF ( k .LT. -1. _d 3 ) THEN
                0183             x = (3. _d 0 * ap1)
                0184      &        * ( - 1. _d 0 + 0.125 _d 0 *(3. _d 0 * ap1 ) / ksq )
                0185            ELSE
                0186             x = TWO * ksq + TWO *k*SQRT(ksq + 3. _d 0 * ap1 )
                0187            ENDIF
                0188            x = (x - (6. _d 0 * ap1) ) * oneNinth + a
                0189 C     original formulation of Zhang 2005:
                0190 C          x = ( TWO * ksq + TWO *k*sqrt(ksq + 3. _d 0 * ap1 )
                0191 C      &          - (6. _d 0 * ap1) ) * oneNinth + a
                0192
                0193 C     Teardrop yield curve:
                0194 C     capping x at a is probably better that setting it to a for x>1,
                0195 C     but still not ideal (because non-differentiable)
                0196 C     possible alternative: log(exp(rho*x)+exp(rho*a))/rho, with rho=O(100)
                0197            x = MIN( x, a )
                0198 C          x = MAX( x,-1 _d 0 ) ! to be sure
                0199
                0200 #ifdef SEAICE_DELTA_SMOOTHREG
                0201            zeta(i,j,bi,bj) = (x + cyc) * press0(i,j,bi,bj)
                0202      &          / SIGN( SQRT(ep*ep+smallNbrSq), ep)
                0203 #else
                0204            zeta(i,j,bi,bj) = (x + cyc) * press0(i,j,bi,bj)
                0205      &          / SIGN( MAX( ABS(ep),smallNbr ), ep)
                0206 #endif /* SEAICE_DELTA_SMOOTHREG */
                0207 C      original formulation of Zhang 2005
                0208 C          zeta(i,j,bi,bj) = (x + 0.5 _d 0) * press0(i,j,bi,bj)
                0209 C      &        / SIGN( ABS(ep) + SEAICE_deltaMin, ep)
                0210 C      &        / SIGN( MAX( ABS(ep),deltaMinSq ), 0.5 _d 0 * ep)
                0211
                0212            eta(i,j,bi,bj) = -(x-a) * SQRT(ONE+x) * press0(i,j,bi,bj)
                0213      &                     * recip_shear(i,j)
                0214
                0215 C      define the max values for eta and zeta
8e32c48b8f Mart*0216            etamax_TD = zMax(i,j,bi,bj) * MIN(eta(i,j,bi,bj)
c512e371cc drin*0217      &               /(ABS(zeta(i,j,bi,bj)) + smallNbr), ONE)
8e32c48b8f Mart*0218            zetamax_TD =zMax(i,j,bi,bj) * MIN(zeta(i,j,bi,bj)
c512e371cc drin*0219      &               /(ABS(eta(i,j,bi,bj)) + smallNbr), ONE )
                0220
                0221 C      apply the max values
                0222            zeta(i,j,bi,bj) = MIN(zeta(i,j,bi,bj), zetamax_TD)
                0223            eta (i,j,bi,bj) = MIN( eta(i,j,bi,bj),  etamax_TD)
                0224 C      original formulation of Zhang 2005
8e32c48b8f Mart*0225 C          zeta(i,j,bi,bj) = MIN(zeta(i,j,bi,bj), zMax(i,j,bi,bj))
                0226 C          eta (i,j,bi,bj) = MIN( eta(i,j,bi,bj), zMax(i,j,bi,bj))
c512e371cc drin*0227
                0228 C      compute the replacement pressure
                0229            press(i,j,bi,bj) = TWO * cyc * (
                0230      &          press0(i,j,bi,bj)*( ONE - SEAICEpressReplFac )
                0231      &         + SEAICEpressReplFac * ep * zeta(i,j,bi,bj)
                0232      &         / SIGN(MAX(ABS(x+cyc), smallNbr), x+cyc) )
                0233
                0234 C      reduce eta according to replacement pressure
                0235            etapr = -(x-a) * SQRT(ONE + x) * press(i,j,bi,bj)/(ONE-a)
                0236      &                    * recip_shear(i,j)
                0237            etapr = MIN( etapr, etamax_TD )
                0238            eta(i,j,bi,bj) =  SEAICEpressReplFac * etapr
                0239      &          + eta(i,j,bi,bj) * ( ONE - SEAICEpressReplFac )
                0240
                0241           ENDDO ! x loop
                0242          ENDDO ! y loop
                0243
                0244         ELSEIF ( SEAICEusePL ) THEN
                0245
                0246          DO j=1-OLy+1,sNy+OLy-1
                0247           DO i=1-OLx+1,sNx+OLx-1
                0248
                0249            ep = e11(i,j,bi,bj)+e22(i,j,bi,bj)
                0250            k = ep * recip_shear(i,j)
                0251            a = tnsFac(i,j,bi,bj)
                0252            x = 0.5 _d 0 * (k - 1. _d 0 + a)
                0253            x = MAX( x, -1. _d 0 )
                0254            x = MIN( x, a )
                0255            xma = x + 0.5 _d 0 *(1. _d 0 - a )
                0256 C     compute zeta
                0257 #ifdef SEAICE_DELTA_SMOOTHREG
                0258            zeta(i,j,bi,bj) = xma * press0(i,j,bi,bj)
                0259      &                    / SIGN( SQRT(ep*ep + smallNbrSq) , ep)
                0260 #else
                0261            zeta(i,j,bi,bj) = xma * press0(i,j,bi,bj)
                0262      &                    / SIGN( MAX(ABS(ep),smallNbr) , ep)
                0263 #endif /* SEAICE_DELTA_SMOOTHREG */
                0264 C      original formulation of Zhang 2005
                0265 C          zeta(i,j,bi,bj) = (x + 0.5 _d 0) * press0(i,j,bi,bj)
                0266 C      &        / SIGN( ABS(ep) + SEAICE_deltaMin, ep)
                0267 C      &        / SIGN( MAX( ABS(ep),deltaMinSq ), 0.5 _d 0 * ep)
                0268
                0269 C     compute eta
                0270            eta(i,j,bi,bj) = -(x-a) * (ONE + x) * press0(i,j,bi,bj)
                0271      &                    * recip_shear(i,j)
                0272
                0273 C     maximum of eta and zeta
8e32c48b8f Mart*0274            etamax_TD = zMax(i,j,bi,bj) * MIN(eta(i,j,bi,bj)
c512e371cc drin*0275      &               /(ABS(zeta(i,j,bi,bj)) + smallNbr),ONE )
8e32c48b8f Mart*0276            zetamax_TD =zMax(i,j,bi,bj) * MIN(zeta(i,j,bi,bj)
c512e371cc drin*0277      &               /(ABS(eta(i,j,bi,bj))  + smallNbr),ONE )
                0278
                0279 C     apply the maxiums on zeta and eta
                0280            zeta(i,j,bi,bj) = MIN(zeta(i,j,bi,bj),zetamax_TD)
                0281            eta (i,j,bi,bj) = MIN( eta(i,j,bi,bj), etamax_TD)
                0282
                0283 C     replacement pressure
                0284            press(i,j,bi,bj) =
                0285      &          ( press0(i,j,bi,bj)*( ONE - SEAICEpressReplFac )
                0286      &          + zeta(i,j,bi,bj) * ep * SEAICEpressReplFac
                0287      &          / SIGN(ABS(xma) + smallNbr, xma ) )
                0288      &          * ( ONE - a )
                0289
                0290 C     change eta according to replacement pressure
                0291            etapr = -(x-a) * (ONE + x) * press(i,j,bi,bj)/(ONE-a)
                0292      &                    * recip_shear(i,j)
                0293            etapr = MIN( etapr, etamax_TD )
                0294            eta(i,j,bi,bj) =  SEAICEpressReplFac * etapr
                0295      &          + eta(i,j,bi,bj) * ( ONE - SEAICEpressReplFac )
                0296           ENDDO !i loop
                0297          ENDDO !j loop
                0298 #endif /* SEAICE_ALLOW_TEARDROP */
                0299 #ifdef SEAICE_ALLOW_MCS
                0300         ELSEIF ( SEAICEuseMCS ) THEN
                0301 C     Full Mohr-Coulomb following IP et al. 1991
                0302          DO j=1-OLy+1,sNy+OLy-1
                0303           DO i=1-OLx+1,sNx+OLx-1
                0304
                0305 C     compute eI = ep
                0306            ep = e11(i,j,bi,bj)+e22(i,j,bi,bj)
                0307
                0308 C     compute zeta
                0309 #ifdef SEAICE_DELTA_SMOOTHREG
                0310            zeta(i,j,bi,bj) = press0(i,j,bi,bj)*(ONE+tnsFac(i,j,bi,bj))
                0311      &                      / (TWO * SQRT(ep*ep+deltaMinSq) )
                0312 #else
                0313            zeta(i,j,bi,bj) = press0(i,j,bi,bj)*(ONE+tnsFac(i,j,bi,bj))
                0314      &                      / (TWO * MAX(ABS(ep),SEAICE_deltaMin) )
                0315 #endif /* SEAICE_DELTA_SMOOTHREG */
                0316
                0317 C     replacement pressure
                0318            press(i,j,bi,bj) = ( ONE - tnsFac(i,j,bi,bj) )
                0319      &      * ( press0(i,j,bi,bj) * ( ONE - SEAICEpressReplFac )
                0320      &      + SEAICEpressReplFac * TWO * zeta(i,j,bi,bj) * ABS(ep)
                0321      &      / ( ONE + tnsFac(i,j,bi,bj) ) )
                0322
8e32c48b8f Mart*0323 C     compute eta  (eMax=zMax)
c512e371cc drin*0324            eta(i,j,bi,bj) = SEAICEmcMU * (0.5 _d 0 * press(i,j,bi,bj)
                0325      &      - zeta(i,j,bi,bj)*ep+press0(i,j,bi,bj)*tnsFac(i,j,bi,bj))
                0326      &       * recip_shear(i,j)
                0327
                0328 C     maximum for eta (high compressive stresses)
8e32c48b8f Mart*0329            eta(i,j,bi,bj) = MIN(eta(i,j,bi,bj) , zMax(i,j,bi,bj))
c512e371cc drin*0330
                0331           ENDDO
                0332          ENDDO
                0333 #endif /* SEAICE_ALLOW_MCS */
                0334         ELSE
                0335 C     For all elliptic yield curves, the computation of deltaC and zeta
                0336 C     is identical so we do it first.
                0337          DO j=1-OLy+1,sNy+OLy-1
                0338           DO i=1-OLx+1,sNx+OLx-1
                0339            ep = e11(i,j,bi,bj)+e22(i,j,bi,bj)
                0340            em = e11(i,j,bi,bj)-e22(i,j,bi,bj)
                0341            shearDefSq = em*em + 4. _d 0*e12Csq(i,j)
                0342            deltaCsq   = ep*ep + recip_efr4*shearDefSq
314f7b043c Mart*0343 CML The old formulation does not ensure that deltaC**2 is always positive,
                0344 CML but in case you need it to reproduce old results, here it is:
c512e371cc drin*0345 CML          deltaCsq =
314f7b043c Mart*0346 CML     &         (e11(i,j,bi,bj)**2+e22(i,j,bi,bj)**2)*(ONE+recip_e2)
f3224168cd Mart*0347 CML     &         + 4. _d 0*recip_e2*e12Csq(i,j)
314f7b043c Mart*0348 CML     &         + 2. _d 0*e11(i,j,bi,bj)*e22(i,j,bi,bj)*(ONE-recip_e2)
8377b8ee87 Mart*0349 #ifdef ALLOW_AUTODIFF
e29ea25404 Mart*0350 C     avoid sqrt of 0
c512e371cc drin*0351            deltaC(i,j,bi,bj) = 0. _d 0
                0352            IF ( deltaCsq .GT. 0. _d 0 )
                0353      &          deltaC(i,j,bi,bj) = SQRT(deltaCsq)
e29ea25404 Mart*0354 #else
c512e371cc drin*0355            deltaC(i,j,bi,bj) = SQRT(deltaCsq)
8377b8ee87 Mart*0356 #endif /* ALLOW_AUTODIFF */
f3224168cd Mart*0357 #ifdef SEAICE_DELTA_SMOOTHREG
                0358 C     smooth regularization (without max-function) of delta for
                0359 C     better differentiability
c512e371cc drin*0360            deltaCreg(i,j) = SQRT(deltaCsq + deltaMinSq)
                0361 CML          deltaCreg(i,j) = deltaC(i,j,bi,bj) + SEAICE_deltaMin
f3224168cd Mart*0362 #else
c512e371cc drin*0363            deltaCreg(i,j) = MAX(deltaC(i,j,bi,bj),SEAICE_deltaMin)
f3224168cd Mart*0364 #endif /* SEAICE_DELTA_SMOOTHREG */
510b9c42ef Mart*0365 #ifdef SEAICE_ZETA_SMOOTHREG
e826e89d88 Jean*0366 C     regularize zeta to zmax with a smooth tanh-function instead
510b9c42ef Mart*0367 C     of a min(zeta,zmax). This improves convergence of iterative
                0368 C     solvers (Lemieux and Tremblay 2009, JGR). No effect on EVP
c512e371cc drin*0369            argTmp = exp(-1. _d 0/(deltaCreg(i,j)*SEAICE_zetaMaxFac))
8e32c48b8f Mart*0370            zeta (i,j,bi,bj) = zMax(i,j,bi,bj)
c512e371cc drin*0371      &          *(1. _d 0 - argTmp)/(1. _d 0 + argTmp)
                0372      &          *(1. _d 0 + tnsFac(i,j,bi,bj) )
510b9c42ef Mart*0373 #else
c512e371cc drin*0374            zeta (i,j,bi,bj) = HALF*( press0(i,j,bi,bj)
                0375      &          * ( 1. _d 0 + tnsFac(i,j,bi,bj) )
                0376      &          )/deltaCreg(i,j)
52b1b58d3c Mart*0377 C     put min and max viscosities in
8e32c48b8f Mart*0378            zeta (i,j,bi,bj) = MIN(zMax(i,j,bi,bj),zeta(i,j,bi,bj))
510b9c42ef Mart*0379 #endif /*  SEAICE_ZETA_SMOOTHREG */
8e32c48b8f Mart*0380            zeta (i,j,bi,bj) = MAX(zMin(i,j,bi,bj),zeta(i,j,bi,bj))
ec0d7df165 Mart*0381 C     set viscosities to zero at HEFFM flow pts
c512e371cc drin*0382            zeta (i,j,bi,bj) = zeta(i,j,bi,bj)*HEFFM(i,j,bi,bj)
52b1b58d3c Mart*0383 C     "replacement pressure"
c512e371cc drin*0384            press(i,j,bi,bj) =
                0385      &          ( press0(i,j,bi,bj)*( 1. _d 0 - SEAICEpressReplFac )
                0386      &          + TWO*zeta(i,j,bi,bj)*deltaC(i,j,bi,bj)
                0387      &          * SEAICEpressReplFac/( 1. _d 0 + tnsFac(i,j,bi,bj) )
                0388      &          ) * ( 1. _d 0 - tnsFac(i,j,bi,bj) )
                0389 CML          press(i,j,bi,bj) = press0(i,j,bi,bj) *
2f5e8addfd Mart*0390 CML     &         ( 1. _d 0 + SEAICEpressReplFac
c512e371cc drin*0391 CML     &                  * ( deltaC(i,j,bi,bj)/deltaCreg(i,j) - 1. _d 0 )
                0392 CML     &         ) * ( 1. _d 0 - tnsFac(i,j,bi,bj) )
dadd13178c Mart*0393           ENDDO
                0394          ENDDO
c512e371cc drin*0395 C     The elliptical yield curves differ in the way eta is computed.
                0396 #if (defined SEAICE_ALLOW_MCE) || (defined SEAICE_ALLOW_TEM)
                0397          IF ( SEAICEuseMCE .OR. SEAICEuseTEM ) THEN
                0398 C     MC yield curve with non-normal flow rule (elliptical plastic potential)
                0399 C     or Truncated Ellipse method with Coulombic limbs
                0400           DO j=1-OLy+1,sNy+OLy-1
                0401            DO i=1-OLx+1,sNx+OLx-1
                0402             ep = e11(i,j,bi,bj)+e22(i,j,bi,bj)
                0403 C     compute eta:
                0404 C     In principle, we need to recompute zeta here, because we need the
                0405 C     unlimited version
                0406 C          zetaLoc = HALF*( press0(i,j,bi,bj)
                0407 C    &          * ( 1. _d 0 + tnsFac(i,j,bi,bj) )
                0408 C    &          )/deltaCreg(i,j)
                0409 C           eta(i,j,bi,bj) = SEAICEmcMU * (0.5 * press0(i,j,bi,bj)
                0410 C    &           * ( 1. _d 0 - tnsFac(i,j,bi,bj) ) - zetaLoc * ep
                0411 C    &           + press0(i,j,bi,bj)*tnsFac(i,j,bi,bj)) / shearDef
                0412 C     but in this formulation, only ep/deltaCreg(i,j) remains
                0413             eta(i,j,bi,bj) = SEAICEmcMU * HALF * press0(i,j,bi,bj)
                0414      &           * ( 1. _d 0 + tnsFac(i,j,bi,bj) )
                0415      &           * ( 1. _d 0 - ep/deltaCreg(i,j) ) * recip_shear(i,j)
8e32c48b8f Mart*0416 C     compute etaMaxFac = zMax/zetaLoc
c512e371cc drin*0417             etaMaxFac = SEAICE_zetaMaxFac * 2. _d 0 * deltaCreg(i,j)
                0418      &           / ( 1. _d 0 + tnsFac(i,j,bi,bj) )
                0419 C     apply maximum from Mohr-Coulomb limbs
                0420             eta(i,j,bi,bj) = eta(i,j,bi,bj) * MIN( 1. _d 0, etaMaxFac )
                0421 C     apply maximum from Elliptical formulation (to be sure)
                0422             eta(i,j,bi,bj) =
8e32c48b8f Mart*0423      &           MIN( eta(i,j,bi,bj), zMax(i,j,bi,bj)*recip_efr2 )
c512e371cc drin*0424            ENDDO
                0425           ENDDO
                0426 C     close the coulombic limbs with the ellipse for compression
                0427           IF ( SEAICEuseTEM) THEN
                0428            DO j=1-OLy+1,sNy+OLy-1
                0429             DO i=1-OLx+1,sNx+OLx-1
                0430              eta(i,j,bi,bj) = MIN( eta(i,j,bi,bj),
                0431      &                             zeta(i,j,bi,bj) * recip_efr2 )
                0432             ENDDO
                0433            ENDDO
                0434           ENDIF
                0435          ELSE
                0436 #endif /* SEAICE_ALLOW_MCE or SEAICE_ALLOW_TEM */
                0437 C     default elliptical yield curve
                0438           DO j=1-OLy+1,sNy+OLy-1
                0439            DO i=1-OLx+1,sNx+OLx-1
                0440             eta(i,j,bi,bj) = zeta(i,j,bi,bj) * recip_efr2
                0441            ENDDO
                0442           ENDDO
                0443 #if (defined SEAICE_ALLOW_MCE) || (defined SEAICE_ALLOW_TEM)
                0444 C     end of elliptical yield curve options
                0445          ENDIF
                0446 #endif /* SEAICE_ALLOW_MCE or SEAICE_ALLOW_TEM */
                0447 C     end of yield curve options
dadd13178c Mart*0448         ENDIF
c512e371cc drin*0449
dc66f404d9 Mart*0450 C     compute eta at Z-points by simple averaging
e826e89d88 Jean*0451         DO j=1-OLy+1,sNy+OLy-1
                0452          DO i=1-OLx+1,sNx+OLx-1
c512e371cc drin*0453           sumNorm  = HEFFM(i,j,  bi,bj)+HEFFM(i-1,j,  bi,bj)
                0454      &             + HEFFM(i,j-1,bi,bj)+HEFFM(i-1,j-1,bi,bj)
dc66f404d9 Mart*0455           IF ( sumNorm.GT.0. _d 0 ) sumNorm = 1. _d 0 / sumNorm
c512e371cc drin*0456           etaZ(i,j,bi,bj) = sumNorm *
                0457      &         ( eta (i,j  ,bi,bj)  + eta (i-1,j  ,bi,bj)
                0458      &         + eta (i,j-1,bi,bj)  + eta (i-1,j-1,bi,bj) )
                0459           zetaZ(i,j,bi,bj) = sumNorm *
                0460      &         ( zeta(i,j  ,bi,bj)  + zeta(i-1,j  ,bi,bj)
                0461      &         + zeta(i,j-1,bi,bj)  + zeta(i-1,j-1,bi,bj) )
ad2b15ed4c Mart*0462          ENDDO
dc66f404d9 Mart*0463         ENDDO
52b1b58d3c Mart*0464 C     free-slip means no lateral stress, which is best achieved by masking
ad2b15ed4c Mart*0465 C     eta on vorticity(=Z)-points; from now on we only need to worry
                0466 C     about the no-slip boundary conditions
                0467         IF (.NOT.SEAICE_no_slip) THEN
c512e371cc drin*0468          DO j=1-OLy+1,sNy+OLy-1
                0469           DO i=1-OLx+1,sNx+OLx-1
                0470            maskZ = HEFFM(i,j,  bi,bj)*HEFFM(i-1,j,  bi,bj)
                0471      &          *  HEFFM(i,j-1,bi,bj)*HEFFM(i-1,j-1,bi,bj)
                0472            etaZ (i,j,bi,bj) =  etaZ(i,j,bi,bj) * maskZ
                0473            zetaZ(i,j,bi,bj) = zetaZ(i,j,bi,bj) * maskZ
ad2b15ed4c Mart*0474           ENDDO
                0475          ENDDO
                0476         ENDIF
3ed8349f04 Mart*0477        ENDDO
                0478       ENDDO
                0479
45315406aa Mart*0480 #endif /* SEAICE_CGRID */
3ed8349f04 Mart*0481       RETURN
                0482       END