C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_calc_viscosities.F,v 1.36 2016/05/17 15:26:46 mlosch Exp $ C $Name: $ #include "SEAICE_OPTIONS.h" #ifdef ALLOW_AUTODIFF # include "AUTODIFF_OPTIONS.h" #endif CBOP CStartOfInterface SUBROUTINE SEAICE_CALC_VISCOSITIES( I e11, e22, e12, zMin, zMax, hEffM, press0, tnsFac, O eta, etaZ, zeta, zetaZ, press, deltaC, I iStep, myTime, myIter, myThid ) C *==========================================================* C | SUBROUTINE SEAICE_CALC_VISCOSITIES | C | o compute shear and bulk viscositites eta, zeta and the | C | corrected ice strength P | C | (see Zhang and Hibler, JGR, 102, 8691-8702, 1997) | C *==========================================================* C | written by Martin Losch, Mar 2006 | C *==========================================================* IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "SEAICE_SIZE.h" #include "SEAICE_PARAMS.h" C === Routine arguments === C iStep :: Sub-time-step number C myTime :: Simulation time C myIter :: Simulation timestep number C myThid :: My Thread Id. number INTEGER iStep _RL myTime INTEGER myIter INTEGER myThid C strain rate tensor _RL e11 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL e22 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL e12 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C _RL zMin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL zMax (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL hEffM (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C _RL press0(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL press (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C _RL deltaC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C factor k to compute tensile strength from k*press0 _RL tnsFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C bulk viscosity _RL eta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C shear viscosity _RL zeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL zetaZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) CEndOfInterface #if ( defined (SEAICE_CGRID) defined (SEAICE_ALLOW_DYNAMICS) ) C === Local variables === C i,j,bi,bj - Loop counters C e11, e12, e22 - components of strain rate tensor C recip_e2 - inverse of square of ratio of yield curve main axes C ep - e11+e22 (abbreviation) C em - e11-e22 (abbreviation) INTEGER i, j, bi, bj _RL recip_e2, deltaCreg, deltaCsq, deltaMinSq, tmp, ep, em _RL e12Csq(1-OLx:sNx+OLx,1-OLy:sNy+OLy) #ifdef SEAICE_ALLOW_TEM _RL etaMax, etaDen #endif /* SEAICE_ALLOW_TEM */ INTEGER k _RL sumNorm, maskZ #ifdef SEAICE_ZETA_SMOOTHREG _RL argTmp #endif /* SEAICE_ZETA_SMOOTHREG */ CEOP C-- FIRST SET UP BASIC CONSTANTS k=1 recip_e2=0. _d 0 IF ( SEAICE_eccen .NE. 0. _d 0 ) recip_e2=ONE/(SEAICE_eccen**2) deltaMinSq = SEAICE_deltaMin**2 C DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) #ifdef ALLOW_AUTODIFF_TAMC DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx e12Csq(i,j) = 0. _d 0 ENDDO ENDDO #endif /* ALLOW_AUTODIFF_TAMC */ C need to do this beforehand for easier vectorization after C TAFization IF ( SEAICEetaZmethod .EQ. 0 ) THEN DO j=1-OLy+1,sNy+OLy-1 DO i=1-OLx+1,sNx+OLx-1 tmp = 0.25 * & ( e12(I,J ,bi,bj) + e12(I+1,J ,bi,bj) & + e12(I,J+1,bi,bj) + e12(I+1,J+1,bi,bj) ) e12Csq(i,j) = tmp*tmp ENDDO ENDDO ELSEIF ( SEAICEetaZmethod .EQ. 3 ) THEN DO j=1-OLy+1,sNy+OLy-1 DO i=1-OLx+1,sNx+OLx-1 C area weighted average of the squares of e12 is more accurate C (and energy conserving) e12Csq(i,j) = 0.25 _d 0 * recip_rA(I,J,bi,bj) * & ( rAz(I ,J ,bi,bj)*e12(I ,J ,bi,bj)**2 & + rAz(I+1,J ,bi,bj)*e12(I+1,J ,bi,bj)**2 & + rAz(I ,J+1,bi,bj)*e12(I ,J+1,bi,bj)**2 & + rAz(I+1,J+1,bi,bj)*e12(I+1,J+1,bi,bj)**2 ) ENDDO ENDDO ENDIF DO j=1-OLy+1,sNy+OLy-1 DO i=1-OLx+1,sNx+OLx-1 ep = e11(i,j,bi,bj)+e22(i,j,bi,bj) em = e11(i,j,bi,bj)-e22(i,j,bi,bj) deltaCsq = ep*ep+recip_e2*(em*em+4. _d 0*e12Csq(i,j)) CML The old formulation does not ensure that deltaC**2 is always positive, CML but in case you need it to reproduce old results, here it is: CML deltaCsq = CML & (e11(i,j,bi,bj)**2+e22(i,j,bi,bj)**2)*(ONE+recip_e2) CML & + 4. _d 0*recip_e2*e12Csq(i,j) CML & + 2. _d 0*e11(i,j,bi,bj)*e22(i,j,bi,bj)*(ONE-recip_e2) #ifdef ALLOW_AUTODIFF_TAMC C avoid sqrt of 0 deltaC(I,J,bi,bj) = 0. _d 0 IF ( deltaCsq .GT. 0. _d 0 ) & deltaC(I,J,bi,bj) = SQRT(deltaCsq) #else deltaC(I,J,bi,bj) = SQRT(deltaCsq) #endif /* ALLOW_AUTODIFF_TAMC */ #ifdef SEAICE_DELTA_SMOOTHREG C smooth regularization (without max-function) of delta for C better differentiability deltaCreg = SQRT(deltaCsq + deltaMinSq) CML deltaCreg = deltaC(I,J,bi,bj) + SEAICE_deltaMin #else deltaCreg = MAX(deltaC(I,J,bi,bj),SEAICE_deltaMin) #endif /* SEAICE_DELTA_SMOOTHREG */ #ifdef SEAICE_ZETA_SMOOTHREG C regularize zeta to zmax with a smooth tanh-function instead C of a min(zeta,zmax). This improves convergence of iterative C solvers (Lemieux and Tremblay 2009, JGR). No effect on EVP argTmp = exp(-1. _d 0/(deltaCreg*SEAICE_zetaMaxFac)) zeta (I,J,bi,bj) = ZMAX(I,J,bi,bj) & *(1. _d 0 - argTmp)/(1. _d 0 + argTmp) #else zeta (I,J,bi,bj) = HALF*( press0(I,J,bi,bj) & * ( 1. _d 0 + tnsFac(I,J,bi,bj) ) & )/deltaCreg C put min and max viscosities in zeta (I,J,bi,bj) = MIN(ZMAX(I,J,bi,bj),zeta(I,J,bi,bj)) #endif /* SEAICE_ZETA_SMOOTHREG */ zeta (I,J,bi,bj) = MAX(ZMIN(I,J,bi,bj),zeta(I,J,bi,bj)) C set viscosities to zero at hEffM flow pts zeta (I,J,bi,bj) = zeta(I,J,bi,bj)*HEFFM(I,J,bi,bj) eta (I,J,bi,bj) = recip_e2*zeta(I,J,bi,bj) C "replacement pressure" press(I,J,bi,bj) = & ( press0(I,J,bi,bj)*( 1. _d 0 - SEAICEpressReplFac ) & + TWO*zeta(I,J,bi,bj)*deltaC(I,J,bi,bj) & * SEAICEpressReplFac/( 1. _d 0 + tnsFac(I,J,bi,bj) ) & ) * ( 1. _d 0 - tnsFac(I,J,bi,bj) ) CML press(I,J,bi,bj) = press0(I,J,bi,bj) * CML & ( 1. _d 0 + SEAICEpressReplFac CML & * ( deltaC(I,J,bi,bj)/deltaCreg - 1. _d 0 ) CML & ) * ( 1. _d 0 - tnsFac(I,J,bi,bj) ) ENDDO ENDDO #ifdef SEAICE_ALLOW_TEM IF ( SEAICEuseTEM ) THEN DO j=1-OLy+1,sNy+OLy-1 DO i=1-OLx+1,sNx+OLx-1 ep = e11(i,j,bi,bj)+e22(i,j,bi,bj) em = e11(i,j,bi,bj)-e22(i,j,bi,bj) etaDen = em*em + 4. _d 0 * e12Csq(i,j) etaDen = SQRT(MAX(deltaMinSq,etaDen)) etaMax = ( 0.5 _d 0*press(I,J,bi,bj)-zeta(I,J,bi,bj)*ep & )/etaDen eta(I,J,bi,bj) = MIN(eta(I,J,bi,bj),etaMax) ENDDO ENDDO ENDIF #endif /* SEAICE_ALLOW_TEM */ C compute eta at Z-points by simple averaging DO j=1-OLy+1,sNy+OLy-1 DO i=1-OLx+1,sNx+OLx-1 sumNorm = maskC(I,J, k,bi,bj)+maskC(I-1,J, k,bi,bj) & + maskC(I,J-1,k,bi,bj)+maskC(I-1,J-1,k,bi,bj) IF ( sumNorm.GT.0. _d 0 ) sumNorm = 1. _d 0 / sumNorm etaZ(I,J,bi,bj) = sumNorm * & ( eta (I,J ,bi,bj) + eta (I-1,J ,bi,bj) & + eta (I,J-1,bi,bj) + eta (I-1,J-1,bi,bj) ) zetaZ(I,J,bi,bj) = sumNorm * & ( zeta(I,J ,bi,bj) + zeta(I-1,J ,bi,bj) & + zeta(I,J-1,bi,bj) + zeta(I-1,J-1,bi,bj) ) ENDDO ENDDO C free-slip means no lateral stress, which is best achieved by masking C eta on vorticity(=Z)-points; from now on we only need to worry C about the no-slip boundary conditions IF (.NOT.SEAICE_no_slip) THEN DO J=1-OLy+1,sNy+OLy-1 DO I=1-OLx+1,sNx+OLx-1 maskZ = maskC(I,J, k,bi,bj)*maskC(I-1,J, k,bi,bj) & * maskC(I,J-1,k,bi,bj)*maskC(I-1,J-1,k,bi,bj) etaZ (I,J,bi,bj) = etaZ(I,J,bi,bj) * maskZ zetaZ(I,J,bi,bj) = zetaZ(I,J,bi,bj) * maskZ ENDDO ENDDO ENDIF ENDDO ENDDO #endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID */ RETURN END