C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_evp.F,v 1.73 2017/04/28 20:04:57 mlosch Exp $ C $Name: $ #include "SEAICE_OPTIONS.h" #ifdef ALLOW_OBCS # include "OBCS_OPTIONS.h" #else # define OBCS_UVICE_OLD #endif #ifdef ALLOW_AUTODIFF # include "AUTODIFF_OPTIONS.h" #endif CBOP C !ROUTINE: SEAICE_EVP C !INTERFACE: SUBROUTINE SEAICE_EVP( myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE SEAICE_EVP C | o Ice dynamics using an EVP solver following C | E. C. Hunke and J. K. Dukowicz. An C | Elastic-Viscous-Plastic Model for Sea Ice Dynamics, C | J. Phys. Oceanogr., 27, 1849-1867 (1997). C *==========================================================* C | written by Martin Losch, March 2006 C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "DYNVARS.h" #include "GRID.h" #include "SEAICE_SIZE.h" #include "SEAICE_PARAMS.h" #include "SEAICE.h" #ifdef ALLOW_AUTODIFF_TAMC # include "tamc.h" #endif C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myTime :: Simulation time C myIter :: Simulation timestep number C myThid :: my Thread Id. number _RL myTime INTEGER myIter INTEGER myThid #if ( defined (SEAICE_CGRID) defined (SEAICE_ALLOW_EVP) defined (SEAICE_ALLOW_DYNAMICS) ) C === Local variables === C i,j,bi,bj :: Loop counters C kSrf :: vertical index of surface layer C nEVPstep :: number of timesteps within the EVP solver C iEVPstep :: Loop counter C SIN/COSWAT :: sine/cosine of turning angle C (recip_)ecc2 :: (one over) eccentricity squared C recip_evpAlpha :: 1/SEAICE_evpAlpha C recip_deltaT :: 1/SEAICE_deltaTdyn C evpStarFac :: 1 if SEAICEuseEVPstar = .true., 0 otherwise C betaFac :: SEAICE_evpBeta/SEAICE_deltaTdyn=1/SEAICE_deltaTevp C betaFacP1 :: betaFac + evpStarFac/SEAICE_deltaTdyn C e11,e12,e22 :: components of strain rate tensor C seaice_div :: divergence strain rates at C-points times P C /divided by Delta minus 1 C seaice_tension :: tension strain rates at C-points times P C /divided by Delta C seaice_shear :: shear strain rates, defined at Z-points times P C /divided by Delta C sig11, sig22 :: sum and difference of diagonal terms of stress tensor C modification for adaptive alpha and beta C (see Kimmritz, Danilov, Losch 2015 for gamma << alpha beta) C EVPcFac :: SEAICE_deltaTdyn*SEAICEaEVPcStar*(SEAICEaEVPcoeff*PI)**2 C with C SEAICEaEVPcStar:: multiple of stabilty factor: alpha*beta = cstar*gamma C SEAICEaEVPcoeff:: largest stabilized frequency according to C gamma = zeta * (cfac/cellarea)*deltaT/m C with (cfac/cellarea) <= pi**2/cellarea C evpAlphaC/Z :: alpha field on C points and on Z points C := sqrt(cstar gamma) C evpBetaU/V :: beta field on u and on v points C := sqrt(cstar gamma) C evpAlphaMin :: lower limit of alpha and beta, regularisation C to prevent singularities of system matrix, C e.g. when ice concentration is too low. C betaFacP1U/V :: = betaFacP1 in standard case, C with varying beta in the adaptive case C on u and on v point C betaFacU/V :: analog betaFacP1U/V INTEGER i, j, bi, bj INTEGER kSrf INTEGER nEVPstep, iEVPstep #ifdef ALLOW_AUTODIFF_TAMC INTEGER ikeyloc #else INTEGER nEVPstepMax #endif _RL COSWAT _RS SINWAT _RL ecc2, recip_ecc2, recip_evpAlpha, recip_deltaT _RL betaFacP1, betaFac, evpStarFac, evpRevFac, recip_evpRevFac _RL seaice_div (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL seaice_tension(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL seaice_shear (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL sig11 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL sig22 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) C fractional area at velocity points _RL areaW (1:sNx,1:sNy,nSx,nSy) _RL areaS (1:sNx,1:sNy,nSx,nSy) C auxilliary variables _RL ep (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL em (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL e12Csq (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL pressC (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL zetaC (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL deltaZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy) C _RL zetaZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy) #ifdef SEAICE_ALLOW_MOM_ADVECTION C tendency due to advection of momentum _RL gUmom (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL gVmom (1-OLx:sNx+OLx,1-OLy:sNy+OLy) #endif /* SEAICE_ALLOW_MOM_ADVECTION */ _RL deltaCreg, deltaSq, deltaMinSq, tmp #ifdef SEAICE_ALLOW_TEM _RL etaDenC, zetaMaxC, etaDenZ, zetaMaxZ #endif /* SEAICE_ALLOW_TEM */ #ifdef SEAICE_ALLOW_CLIPZETA _RL zMaxZ, zMinZ, fac #endif /* SEAICE_ALLOW_CLIPZETA */ _RL denom1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL denom2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL sumNorm, denomU, denomV _RL locMaskU, locMaskV _RL EVPcFac _RL evpAlphaC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL evpAlphaZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL evpBetaU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL evpBetaV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL betaFacP1U, betaFacP1V _RL betaFacU, betaFacV LOGICAL useAdaptiveEVP #ifdef ALLOW_SEAICE_EVP_RESIDUAL _RL resTile(nSx,nSy) _RL resLoc _RL uIcePm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vIcePm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL sig11pm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL sig22pm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL sig12pm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) LOGICAL printResidual C CHARACTER*(10) suff C !FUNCTIONS: LOGICAL DIFFERENT_MULTIPLE EXTERNAL #endif /* ALLOW_SEAICE_EVP_RESIDUAL */ CEOP C set tuning parameters for adaptive EVP useAdaptiveEVP = .FALSE. IF ( SEAICEaEvpCoeff .NE. UNSET_RL ) useAdaptiveEVP = .TRUE. EVPcFac = 0. _d 0 IF ( useAdaptiveEVP ) & EVPcFac = SEAICE_deltaTdyn*SEAICEaEVPcStar & * (SEAICEaEvpCoeff * PI)**2 #ifdef ALLOW_SEAICE_EVP_RESIDUAL printResidual = debugLevel.GE.debLevA & .AND. DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) #endif /* ALLOW_SEAICE_EVP_RESIDUAL */ C surface level kSrf = 1 C-- introduce turning angles SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad) COSWAT=COS(SEAICE_waterTurnAngle*deg2rad) C abbreviation eccentricity squared ecc2=SEAICE_eccen**2 recip_ecc2 = 0. _d 0 IF ( ecc2 .NE. 0. _d 0 ) recip_ecc2=ONE/ecc2 deltaMinSq = SEAICE_deltaMin**2 C determine number of internal time steps nEVPstep = INT(SEAICE_deltaTdyn/SEAICE_deltaTevp) IF ( SEAICEnEVPstarSteps.NE.UNSET_I ) nEVPstep=SEAICEnEVPstarSteps C SEAICE_evpAlpha = 2. * SEAICE_evpTauRelax/SEAICE_deltaTevp DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx denom1(I,J,bi,bj) = 1. _d 0 / ( SEAICE_evpAlpha + 1. _d 0 ) denom2(I,J,bi,bj) = 1. _d 0 / ( SEAICE_evpAlpha + ecc2 ) ENDDO ENDDO ENDDO ENDDO recip_deltaT = 1. _d 0 / SEAICE_deltaTdyn recip_evpAlpha = 0. _d 0 IF ( SEAICE_evpAlpha .GT. 0. _d 0 ) & recip_evpAlpha = 1. _d 0 / SEAICE_evpAlpha evpStarFac = 0. _d 0 evpRevFac = 0. _d 0 recip_evpRevFac = 1. _d 0 IF ( SEAICEuseEVPstar ) evpStarFac = 1. _d 0 IF ( SEAICEuseEVPrev ) THEN C the Bouillon et al. (2013) discretization in time has more C explicit terms evpRevFac = 1. _d 0 evpStarFac = 1. _d 0 recip_evpRevFac = recip_ecc2 DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx denom1(I,J,bi,bj) = 1. _d 0 / SEAICE_evpAlpha denom2(I,J,bi,bj) = denom1(i,j,bi,bj) ENDDO ENDDO ENDDO ENDDO ENDIF C betaFac = SEAICE_evpBeta*recip_deltaT betaFacU = betaFac betaFacV = betaFac C betaFacP1 = betaFac + evpStarFac*recip_deltaT betaFacP1U = betaFacP1 betaFacP1V = betaFacP1 #ifndef ALLOW_AUTODIFF_TAMC nEVPstepMax = nEVPstep #endif DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx C use u/vIce as work arrays: copy previous time step to u/vIceNm1 uIceNm1(I,J,bi,bj) = uIce(I,J,bi,bj) vIceNm1(I,J,bi,bj) = vIce(I,J,bi,bj) C initialise strain rates e11 (I,J,bi,bj) = 0. _d 0 e22 (I,J,bi,bj) = 0. _d 0 e12 (I,J,bi,bj) = 0. _d 0 C initialise adative-EVP-specific fields evpAlphaC(I,J,bi,bj) = SEAICE_evpAlpha evpAlphaZ(I,J,bi,bj) = SEAICE_evpAlpha evpBetaU (I,J,bi,bj) = SEAICE_evpBeta evpBetaV (I,J,bi,bj) = SEAICE_evpBeta ENDDO ENDDO C initialise fractional areas at velocity points IF ( SEAICEscaleSurfStress ) THEN DO J=1,sNy DO I=1,sNx areaW(I,J,bi,bj) = & 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I-1,J,bi,bj)) areaS(I,J,bi,bj) = & 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I,J-1,bi,bj)) ENDDO ENDDO ELSE DO J=1,sNy DO I=1,sNx areaW(I,J,bi,bj) = 1. _d 0 areaS(I,J,bi,bj) = 1. _d 0 ENDDO ENDDO ENDIF ENDDO ENDDO #ifdef SEAICE_ALLOW_CLIPZETA C damping constraint (Hunke, J.Comp.Phys.,2001) IF ( SEAICE_evpDampC .GT. 0. _d 0 ) THEN CML fac = HALF * SEAICE_evpDampC * SEAICE_evpTauRelax CML & /SEAICE_deltaTevp**2 fac = 0.25 _d 0 * SEAICE_evpDampC * SEAICE_evpAlpha & /SEAICE_deltaTevp DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx zMax (I,J,bi,bj) = _rA(I,J,bi,bj) * fac ENDDO ENDDO ENDDO ENDDO ENDIF #endif /* SEAICE_ALLOW_CLIPZETA */ C C start of the main time loop DO iEVPstep = 1, nEVPstepMax IF (iEVPstep.LE.nEVPstep) THEN C #ifdef ALLOW_AUTODIFF_TAMC ikeyloc = iEVPstep + & (ikey_dynamics-1)*nEVPstepMax CADJ STORE uice = comlev1_evp, key = ikeyloc, byte = isbyte CADJ STORE vice = comlev1_evp, key = ikeyloc, byte = isbyte CADJ STORE seaice_sigma1 = comlev1_evp, key = ikeyloc, byte = isbyte CADJ STORE seaice_sigma2 = comlev1_evp, key = ikeyloc, byte = isbyte CADJ STORE seaice_sigma12 = comlev1_evp, key = ikeyloc, byte = isbyte CADJ STORE evpAlphaC = comlev1_evp, key = ikeyloc, byte = isbyte CADJ STORE evpBetaU = comlev1_evp, key = ikeyloc, byte = isbyte CADJ STORE evpBetaV = comlev1_evp, key = ikeyloc, byte = isbyte #endif /* ALLOW_AUTODIFF_TAMC */ C C first calculate strain rates and bulk moduli/viscosities C CALL SEAICE_CALC_STRAINRATES( I uIce, vIce, O e11, e22, e12, I iEVPstep, myTime, myIter, myThid ) #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE e11 = comlev1_evp,key = ikeyloc, byte = isbyte CADJ STORE e12 = comlev1_evp,key = ikeyloc, byte = isbyte CADJ STORE e22 = comlev1_evp,key = ikeyloc, byte = isbyte #endif /* ALLOW_AUTODIFF_TAMC */ 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 = ikeyloc - 1 iicekey = (act1 + 1) + act2*max1 & + act3*max1*max2 & + act4*max1*max2*max3 #endif /* ALLOW_AUTODIFF_TAMC */ #ifdef ALLOW_SEAICE_EVP_RESIDUAL C save previous (p-1) iteration IF ( printResidual ) THEN DO j=1,sNy DO i=1,sNx sig11Pm1(I,J,bi,bj) = seaice_sigma1(I,J,bi,bj) sig22Pm1(I,J,bi,bj) = seaice_sigma2(I,J,bi,bj) sig12Pm1(I,J,bi,bj) = seaice_sigma12(I,J,bi,bj) uIcePm1 (I,J,bi,bj) = uIce(i,j,bi,bj) vIcePm1 (I,J,bi,bj) = vIce(i,j,bi,bj) ENDDO ENDDO ENDIF #endif /* ALLOW_SEAICE_EVP_RESIDUAL */ DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx seaice_div (I,J) = 0. _d 0 seaice_tension(I,J) = 0. _d 0 seaice_shear (I,J) = 0. _d 0 pressC (I,J) = 0. _d 0 e12Csq (I,J) = 0. _d 0 zetaC (I,J) = 0. _d 0 deltaZ (I,J) = 0. _d 0 zetaZ (I,J,bi,bj) = 0. _d 0 deltaC (I,J,bi,bj) = 0. _d 0 ENDDO ENDDO DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx ep(i,j) = e11(i,j,bi,bj) + e22(i,j,bi,bj) em(i,j) = e11(i,j,bi,bj) - e22(i,j,bi,bj) ENDDO ENDDO C need to do this beforehand for easier vectorization after C TAFization C average strain rates to C points 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) according to Bouillon et al. 2013, eq (11) 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=0,sNy+1 DO i=0,sNx+1 deltaSq = ep(I,J)**2 + recip_ecc2 * em(I,J)**2 & + recip_ecc2 * 4. _d 0 * e12Csq(I,J) #ifdef ALLOW_AUTODIFF_TAMC C avoid sqrt of 0 deltaC(I,J,bi,bj) = 0. _d 0 IF ( deltaSq .GT. 0. _d 0 ) & deltaC(I,J,bi,bj) = SQRT(deltaSq) #else deltaC(I,J,bi,bj) = SQRT(deltaSq) #endif /* ALLOW_AUTODIFF_TAMC */ #ifdef SEAICE_DELTA_SMOOTHREG C smooth regularization (without max-function) of delta for C better differentiability CML deltaCreg = SQRT(deltaSq + deltaMinSq) deltaCreg = deltaC(I,J,bi,bj) + SEAICE_deltaMin #else deltaCreg = MAX(deltaC(I,J,bi,bj),SEAICE_deltaMin) #endif /* SEAICE_DELTA_SMOOTHREG */ zetaC(I,J) = HALF*( press0(I,J,bi,bj) & * ( 1. _d 0 + tensileStrFac(I,J,bi,bj) ) & )/deltaCreg ENDDO ENDDO IF ( useAdaptiveEVP ) THEN DO j=0,sNy+1 DO i=0,sNx+1 CML I do not like these hidden regularisations, why do we need to CML divide by mass? evpAlphaC(I,J,bi,bj) = SQRT(zetaC(I,J) & * EVPcFac / MAX(seaiceMassC(I,J,bi,bj), 1.D-04) & * recip_rA(I,J,bi,bj) ) * maskC(I,J,kSrf,bi,bj) evpAlphaC(I,J,bi,bj) = & MAX(evpAlphaC(I,J,bi,bj),SEAICEaEVPalphaMin) ENDDO ENDDO ENDIF C compute zetaZ and deltaZ by simple averaging (following C Bouillon et al., 2013) DO J=1,sNy+1 DO I=1,sNx+1 sumNorm = maskC(I,J, kSrf,bi,bj)+maskC(I-1,J, kSrf,bi,bj) & + maskC(I,J-1,kSrf,bi,bj)+maskC(I-1,J-1,kSrf,bi,bj) IF ( sumNorm.GT.0. _d 0 ) sumNorm = 1. _d 0 / sumNorm zetaZ(I,J,bi,bj) = sumNorm * & ( zetaC(I, J) + zetaC(I-1,J-1) & + zetaC(I-1,J) + zetaC(I, J-1) ) deltaZ(I,J) = sumNorm * & ( deltaC(I, J,bi,bj) + deltaC(I-1,J-1,bi,bj) & + deltaC(I-1,J,bi,bj) + deltaC(I, J-1,bi,bj) ) ENDDO ENDDO #ifdef SEAICE_ALLOW_CLIPZETA #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE zetac = comlev1_bibj_evp, CADJ & key = iiceloc, byte = isbyte CADJ STORE zetaz(:,:,bi,bj) = comlev1_bibj_evp, CADJ & key = iicekey, byte = isbyte #endif /* ALLOW_AUTODIFF_TAMC */ C regularize zeta if necessary DO j=0,sNy+1 DO i=0,sNx+1 zetaC(I,J) = MAX(zMin(I,J,bi,bj),MIN(zMax(I,J,bi,bj) & ,zetaC(I,J))) CML zetaC(I,J) = zetaC(I,J)*hEffM(I,J,bi,bj) C C zMin, zMax are defined at C-points, make sure that they are not C masked by boundaries/land points zMaxZ = MAX( & MAX(zMax(I, J,bi,bj),zMax(I, J-1,bi,bj)), & MAX(zMax(I-1,J,bi,bj),zMax(I-1,J-1,bi,bj)) ) zMinZ = MAX( & MAX(zMin(I, J,bi,bj),zMin(I, J-1,bi,bj)), & MAX(zMin(I-1,J,bi,bj),zMin(I-1,J-1,bi,bj)) ) zetaZ(I,J,bi,bj) = MAX(zMinZ,MIN(zMaxZ,zetaZ(I,J,bi,bj))) ENDDO ENDDO #endif /* SEAICE_ALLOW_CLIPZETA */ C recompute pressure DO j=0,sNy+1 DO i=0,sNx+1 pressC(I,J) = & ( press0(I,J,bi,bj) * ( 1. _d 0 - SEAICEpressReplFac ) & + TWO*zetaC(I,J)*deltaC(I,J,bi,bj)*SEAICEpressReplFac & /(1. _d 0 + tensileStrFac(I,J,bi,bj)) & ) * (1. _d 0 - tensileStrFac(I,J,bi,bj)) ENDDO ENDDO #ifdef ALLOW_DIAGNOSTICS IF ( useDiagnostics ) THEN C save eta, zeta, and pressure for diagnostics DO j=1,sNy DO i=1,sNx press(I,J,bi,bj) = pressC(I,J) zeta (I,J,bi,bj) = zetaC(I,J) eta (I,J,bi,bj) = zetaC(I,J)*recip_ecc2 ENDDO ENDDO ENDIF #endif /* ALLOW_DIAGNOSTICS */ C Calculate the RHS of the stress equations. Do this now in order to C avoid multiple divisions by Delta C P * ( D_d/Delta - 1 ) = 2*zeta*D_d - P = 2*zeta*D_d - 2*zeta*Delta C P * ( D_t/Delta ) = 2*zeta*D_t C P * ( D_s/Delta ) = 2*zeta*D_s #ifdef SEAICE_ALLOW_TEM IF ( SEAICEuseTEM ) THEN DO j=0,sNy DO i=0,sNx etaDenC = em(I,J)**2 + 4. _d 0 * e12Csq(I,J) etaDenC = SQRT(MAX(deltaMinSq,etaDenC)) zetaMaxC = ecc2*zetaC(I,J) & *(deltaC(I,J,bi,bj)-ep(I,J))/etaDenC #ifdef ALLOW_DIAGNOSTICS C save new eta for diagnostics eta(I,J,bi,bj) = MIN(zetaC(I,J),zetaMaxC)*recip_ecc2 #endif /* ALLOW_DIAGNOSTICS */ seaice_div (I,J) = & ( 2. _d 0 *zetaC(I,J)*ep(I,J) - pressC(I,J) & ) * hEffM(I,J,bi,bj) seaice_tension(I,J) = 2. _d 0*MIN(zetaC(I,J),zetaMaxC) & * em(I,J) * hEffM(I,J,bi,bj) ENDDO ENDDO DO j=1,sNy+1 DO i=1,sNx+1 sumNorm = 0.25 _d 0 CML sumNorm = 1.0 _d 0 CML & / MAX(1.D0,maskC(I, J, kSrf,bi,bj) CML & + maskC(I-1,J, kSrf,bi,bj) CML & + maskC(I, J-1,kSrf,bi,bj) CML & + maskC(I-1,J-1,kSrf,bi,bj) ) C Averaging the squares gives more accurate viscous-plastic behavior C than squaring the averages etaDenZ = & sumNorm * recip_rAz(I,J,bi,bj) * & ( _rA(I ,J ,bi,bj) * em(I, J )**2 & + _rA(I-1,J-1,bi,bj) * em(I-1,J-1)**2 & + _rA(I-1,J ,bi,bj) * em(I-1,J )**2 & + _rA(I ,J-1,bi,bj) * em(I, J-1)**2 ) & + 4. _d 0*e12(I,J,bi,bj)**2 etaDenZ = SQRT(MAX(deltaMinSq,etaDenZ)) zetaMaxZ = ecc2*zetaZ(I,J,bi,bj) * ( deltaZ(I,J) & - sumNorm * ( ep(I,J ) + ep(I-1,J ) & + ep(I,J-1) + ep(I-1,J-1) ) & )/etaDenZ seaice_shear (I,J) = & 2. _d 0*MIN(zetaZ(I,J,bi,bj),zetaMaxZ) & * 2. _d 0*e12(I,J,bi,bj) ENDDO ENDDO ELSE #else IF (.TRUE. ) THEN #endif /* SEAICE_ALLOW_TEM */ DO J=0,sNy DO I=0,sNx seaice_div (I,J) = & ( 2. _d 0 *zetaC(I,J)*ep(I,J) - pressC(I,J) & ) * hEffM(I,J,bi,bj) seaice_tension(I,J) = 2. _d 0*zetaC(I,J) & * em(I,J) * hEffM(I,J,bi,bj) ENDDO ENDDO DO J=1,sNy+1 DO I=1,sNx+1 seaice_shear (I,J) = & 2. _d 0*zetaZ(I,J,bi,bj)*e12(I,J,bi,bj) ENDDO ENDDO ENDIF C C first step stress equations C #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE denom1(:,:,bi,bj) = comlev1_bibj_evp, CADJ & key=iicekey, byte=isbyte CADJ STORE denom2(:,:,bi,bj) = comlev1_bibj_evp, CADJ & key=iicekey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ IF ( useAdaptiveEVP ) THEN DO j=0,sNy DO i=0,sNx denom1(I,J,bi,bj) = 1. _d 0 / evpAlphaC(I,J,bi,bj) denom2(I,J,bi,bj) = denom1(I,J,bi,bj) ENDDO ENDDO ENDIF DO j=0,sNy DO i=0,sNx C sigma1 and sigma2 are computed on C points seaice_sigma1 (I,J,bi,bj) = ( seaice_sigma1 (I,J,bi,bj) & * ( evpAlphaC(I,J,bi,bj) - evpRevFac ) & + seaice_div(I,J) & ) * denom1(I,J,bi,bj) & *hEffM(I,J,bi,bj) seaice_sigma2 (I,J,bi,bj) = ( seaice_sigma2 (I,J,bi,bj) & * ( evpAlphaC(I,J,bi,bj) - evpRevFac ) & + seaice_tension(I,J)*recip_evpRevFac & ) * denom2(I,J,bi,bj) & *hEffM(I,J,bi,bj) #ifdef SEAICE_EVP_ELIMINATE_UNDERFLOWS C Code to avoid very small numbers that can degrade performance. C Many compilers can handle this more efficiently with the help of C a flag (copied from CICE after correspondence with Elizabeth Hunke) seaice_sigma1(I,J,bi,bj) = SIGN(MAX( & ABS( seaice_sigma1(I,J,bi,bj) ), SEAICE_EPS ), & seaice_sigma1(I,J,bi,bj) ) seaice_sigma2(I,J,bi,bj) = SIGN(MAX( & ABS( seaice_sigma2(I,J,bi,bj) ), SEAICE_EPS ), & seaice_sigma2(I,J,bi,bj) ) #endif /* SEAICE_EVP_ELIMINATE_UNDERFLOWS */ C recover sigma11 and sigma22 sig11(I,J) = 0.5 _d 0 * & ( seaice_sigma1(I,J,bi,bj)+seaice_sigma2(I,J,bi,bj) ) sig22(I,J) = 0.5 _d 0 * & ( seaice_sigma1(I,J,bi,bj)-seaice_sigma2(I,J,bi,bj) ) ENDDO ENDDO #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE denom2 (:,:,bi,bj) = comlev1_bibj_evp, CADJ & key=iicekey, byte=isbyte CADJ STORE evpAlphaZ(:,:,bi,bj) = comlev1_bibj_evp, CADJ & key=iicekey, byte=isbyte #endif /* ALLOW_AUTODIFF_TAMC */ C sigma12 is computed on Z points IF ( useAdaptiveEVP ) THEN DO j=1,sNy+1 DO i=1,sNx+1 evpAlphaZ(I,J,bi,bj) = 0.25 _d 0 * & ( evpAlphaC(I, J,bi,bj)+evpAlphaC(I-1,J-1,bi,bj) & + evpAlphaC(I-1,J,bi,bj)+evpAlphaC(I, J-1,bi,bj) ) denom2(I,J,bi,bj) = 1. _d 0 / evpAlphaZ(I,J,bi,bj) ENDDO ENDDO ENDIF DO j=1,sNy+1 DO i=1,sNx+1 seaice_sigma12(I,J,bi,bj) = ( seaice_sigma12(I,J,bi,bj) & * ( evpAlphaZ(I,J,bi,bj) - evpRevFac ) & + seaice_shear(I,J)*recip_evpRevFac & ) * denom2(I,J,bi,bj) #ifdef SEAICE_EVP_ELIMINATE_UNDERFLOWS seaice_sigma12(I,J,bi,bj) = SIGN(MAX( & ABS( seaice_sigma12(I,J,bi,bj) ), SEAICE_EPS ), & seaice_sigma12(I,J,bi,bj) ) #endif /* SEAICE_EVP_ELIMINATE_UNDERFLOWS */ ENDDO ENDDO C C compute divergence of stress tensor C DO J=1,sNy DO I=1,sNx stressDivergenceX(I,J,bi,bj) = & ( sig11(I ,J ) * _dyF(I ,J,bi,bj) & - sig11(I-1,J ) * _dyF(I-1,J,bi,bj) & + seaice_sigma12(I,J+1,bi,bj) * _dxV(I,J+1,bi,bj) & - seaice_sigma12(I,J ,bi,bj) * _dxV(I,J ,bi,bj) & ) * recip_rAw(I,J,bi,bj) stressDivergenceY(I,J,bi,bj) = & ( sig22(I,J ) * _dxF(I,J ,bi,bj) & - sig22(I,J-1) * _dxF(I,J-1,bi,bj) & + seaice_sigma12(I+1,J,bi,bj) * _dyU(I+1,J,bi,bj) & - seaice_sigma12(I ,J,bi,bj) * _dyU(I ,J,bi,bj) & ) * recip_rAs(I,J,bi,bj) ENDDO ENDDO ENDDO ENDDO #ifdef ALLOW_SEAICE_EVP_RESIDUAL IF ( printResidual ) THEN DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) resTile(bi,bj) = 0. _d 0 DO j=1,sNy DO i=1,sNx sig11Pm1(i,j,bi,bj) = & seaice_sigma1(i,j,bi,bj)-sig11pm1(i,j,bi,bj) sig22Pm1(i,j,bi,bj) = & seaice_sigma2(i,j,bi,bj)-sig22pm1(i,j,bi,bj) sig12Pm1(i,j,bi,bj) = & seaice_sigma12(i,j,bi,bj)-sig12Pm1(i,j,bi,bj) sig11Pm1(i,j,bi,bj) = & evpAlphaC(I,J,bi,bj) * sig11Pm1(i,j,bi,bj) sig22Pm1(i,j,bi,bj) = & evpAlphaC(I,J,bi,bj) * sig22Pm1(i,j,bi,bj) sig12Pm1(i,j,bi,bj) = & evpAlphaZ(I,J,bi,bj) * sig12Pm1(i,j,bi,bj) ENDDO ENDDO IF ( .NOT. SEAICEscaleSurfStress ) THEN C multiply with mask (concentration) to only count ice contributions DO j=1,sNy DO i=1,sNx resTile(bi,bj) = resTile(bi,bj) + AREA(i,j,bi,bj) * & ( sig11Pm1(i,j,bi,bj)*sig11Pm1(i,j,bi,bj) & + sig22Pm1(i,j,bi,bj)*sig22Pm1(i,j,bi,bj) & + sig12Pm1(i,j,bi,bj)*sig12Pm1(i,j,bi,bj) ) ENDDO ENDDO ELSE C in this case the scaling with AREA is already done DO j=1,sNy DO i=1,sNx resTile(bi,bj) = resTile(bi,bj) & + sig11Pm1(i,j,bi,bj)*sig11Pm1(i,j,bi,bj) & + sig22Pm1(i,j,bi,bj)*sig22Pm1(i,j,bi,bj) & + sig12Pm1(i,j,bi,bj)*sig12Pm1(i,j,bi,bj) ENDDO ENDDO ENDIF ENDDO ENDDO resloc = 0. _d 0 CALL GLOBAL_SUM_TILE_RL( resTile, resloc, myThid ) resloc = SQRT(resloc) WRITE(standardMessageUnit,'(A,1X,I6,1PE16.8)') & ' SEAICE_EVP: iEVPstep, residual sigma = ', iEVPstep, resLoc ENDIF #endif /* ALLOW_SEAICE_EVP_RESIDUAL */ #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE stressDivergenceX = comlev1_evp, CADJ & key = ikeyloc, byte = isbyte CADJ STORE stressDivergenceY = comlev1_evp, CADJ & key = ikeyloc, byte = isbyte # ifdef SEAICE_DYN_STABLE_ADJOINT Cgf zero out adjoint fields to stabilize pkg/seaice dyna. adjoint CALL ZERO_ADJ( 1, stressDivergenceX, myThid) CALL ZERO_ADJ( 1, stressDivergenceY, myThid) # endif /* SEAICE_DYN_STABLE_ADJOINT */ #endif /* ALLOW_AUTODIFF_TAMC */ C C set up rhs for stepping the velocity field C CALL SEAICE_OCEANDRAG_COEFFS( I uIce, vIce, O DWATN, I iEVPstep, myTime, myIter, myThid ) #ifdef SEAICE_ALLOW_BOTTOMDRAG IF (SEAICEbasalDragK2.GT.0. _d 0) CALL SEAICE_BOTTOMDRAG_COEFFS( I uIce, vIce, #ifdef SEAICE_ITD I HEFFITD, AREAITD, AREA, #else I HEFF, AREA, #endif O CbotC, I iEVPstep, myTime, myIter, myThid ) #endif /* SEAICE_ALLOW_BOTTOMDRAG */ DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1,sNy DO I=1,sNx C over open water, all terms that contain sea ice mass drop out and C the balance is determined by the atmosphere-ice and ice-ocean stress; C the staggering of uIce and vIce can cause stripes in the open ocean C solution when the turning angle is non-zero (SINWAT.NE.0); C we mask this term here in order to avoid the stripes and because C over open ocean, u/vIce do not advect anything, so that the associated C error is small and most likely only confined to the ice edge but may C propagate into the ice covered regions. locMaskU = seaiceMassU(I,J,bi,bj) locMaskV = seaiceMassV(I,J,bi,bj) IF ( locMaskU .NE. 0. _d 0 ) locMaskU = 1. _d 0 IF ( locMaskV .NE. 0. _d 0 ) locMaskV = 1. _d 0 C to recover old results replace the above lines with the line below C locMaskU = 1. _d 0 C locMaskV = 1. _d 0 C set up anti symmetric drag force and add in ice ocean stress C ( remember to average to correct velocity points ) FORCEX(I,J,bi,bj)=FORCEX0(I,J,bi,bj)+ & ( 0.5 _d 0 * ( DWATN(I,J,bi,bj)+DWATN(I-1,J,bi,bj) ) * & COSWAT * uVel(I,J,kSrf,bi,bj) & - SIGN(SINWAT, _fCori(I,J,bi,bj))* 0.5 _d 0 * & ( DWATN(I ,J,bi,bj) * 0.5 _d 0 * & (vVel(I ,J ,kSrf,bi,bj)-vIce(I ,J ,bi,bj) & +vVel(I ,J+1,kSrf,bi,bj)-vIce(I ,J+1,bi,bj)) & + DWATN(I-1,J,bi,bj) * 0.5 _d 0 * & (vVel(I-1,J ,kSrf,bi,bj)-vIce(I-1,J ,bi,bj) & +vVel(I-1,J+1,kSrf,bi,bj)-vIce(I-1,J+1,bi,bj)) & )*locMaskU ) * areaW(I,J,bi,bj) FORCEY(I,J,bi,bj)=FORCEY0(I,J,bi,bj)+ & ( 0.5 _d 0 * ( DWATN(I,J,bi,bj)+DWATN(I,J-1,bi,bj) ) * & COSWAT * vVel(I,J,kSrf,bi,bj) & + SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 * & ( DWATN(I,J ,bi,bj) * 0.5 _d 0 * & (uVel(I ,J ,kSrf,bi,bj)-uIce(I ,J ,bi,bj) & +uVel(I+1,J ,kSrf,bi,bj)-uIce(I+1,J ,bi,bj)) & + DWATN(I,J-1,bi,bj) * 0.5 _d 0 * & (uVel(I ,J-1,kSrf,bi,bj)-uIce(I ,J-1,bi,bj) & +uVel(I+1,J-1,kSrf,bi,bj)-uIce(I+1,J-1,bi,bj)) & )*locMaskV ) * areaS(I,J,bi,bj) C coriols terms FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj) + 0.5 _d 0*( & seaiceMassC(I ,J,bi,bj) * _fCori(I ,J,bi,bj) & * 0.5 _d 0*( vIce(I ,J,bi,bj)+vIce(I ,J+1,bi,bj) ) & + seaiceMassC(I-1,J,bi,bj) * _fCori(I-1,J,bi,bj) & * 0.5 _d 0*( vIce(I-1,J,bi,bj)+vIce(I-1,J+1,bi,bj) ) & ) FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj) - 0.5 _d 0*( & seaiceMassC(I,J ,bi,bj) * _fCori(I,J ,bi,bj) & * 0.5 _d 0*( uIce(I,J ,bi,bj)+uIce(I+1, J,bi,bj) ) & + seaiceMassC(I,J-1,bi,bj) * _fCori(I,J-1,bi,bj) & * 0.5 _d 0*( uIce(I,J-1,bi,bj)+uIce(I+1,J-1,bi,bj) ) & ) ENDDO ENDDO #ifdef SEAICE_ALLOW_MOM_ADVECTION IF ( SEAICEmomAdvection ) THEN DO J=1-Oly,sNy+Oly DO I=1-Olx,sNx+Olx gUmom(I,J) = 0. _d 0 gVmom(I,J) = 0. _d 0 ENDDO ENDDO CALL SEAICE_MOM_ADVECTION( I bi,bj,1,sNx,1,sNy, I uIce, vIce, O gUmom, gVmom, I myTime, myIter, myThid ) DO J=1,sNy DO I=1,sNx FORCEX(I,J,bi,bj) = FORCEX(I,J,bi,bj) + gUmom(I,J) FORCEY(I,J,bi,bj) = FORCEY(I,J,bi,bj) + gVmom(I,J) ENDDO ENDDO ENDIF #endif /* SEAICE_ALLOW_MOM_ADVECTION */ C C step momentum equations with ice-ocean stress term treated implicitly C IF ( useAdaptiveEVP ) THEN DO J=1,sNy DO I=1,sNx C compute and adjust parameters that are constant otherwise evpBetaU(I,J,bi,bj) = 0.5 _d 0*(evpAlphaC(I-1,J,bi,bj) & +evpAlphaC(I, J,bi,bj)) evpBetaV(I,J,bi,bj) = 0.5 _d 0*(evpAlphaC(I,J-1,bi,bj) & +evpAlphaC(I,J, bi,bj)) ENDDO ENDDO ENDIF DO J=1,sNy DO I=1,sNx betaFacU = evpBetaU(I,J,bi,bj)*recip_deltaT betaFacV = evpBetaV(I,J,bi,bj)*recip_deltaT tmp = evpStarFac*recip_deltaT betaFacP1V = betaFacV + tmp betaFacP1U = betaFacU + tmp denomU = seaiceMassU(I,J,bi,bj)*betaFacP1U & + 0.5 _d 0*( DWATN(I,J,bi,bj) + DWATN(I-1,J,bi,bj) ) & * COSWAT * areaW(I,J,bi,bj) denomV = seaiceMassV(I,J,bi,bj)*betaFacP1V & + 0.5 _d 0*( DWATN(I,J,bi,bj) + DWATN(I,J-1,bi,bj) ) & * COSWAT * areaS(I,J,bi,bj) #ifdef SEAICE_ALLOW_BOTTOMDRAG denomU = denomU + areaW(I,J,bi,bj) & * 0.5 _d 0*( CbotC(I,J,bi,bj) + CbotC(I-1,J,bi,bj) ) denomV = denomV + areaS(I,J,bi,bj) & * 0.5 _d 0*( CbotC(I,J,bi,bj) + CbotC(I,J-1,bi,bj) ) #endif /* SEAICE_ALLOW_BOTTOMDRAG */ IF ( denomU .EQ. 0. _d 0 ) denomU = 1. _d 0 IF ( denomV .EQ. 0. _d 0 ) denomV = 1. _d 0 uIce(I,J,bi,bj) = seaiceMaskU(I,J,bi,bj) * & ( seaiceMassU(I,J,bi,bj)*betaFacU & * uIce(I,J,bi,bj) & + seaiceMassU(I,J,bi,bj)*recip_deltaT*evpStarFac & * uIceNm1(I,J,bi,bj) & + FORCEX(I,J,bi,bj) & + stressDivergenceX(I,J,bi,bj) ) / denomU vIce(I,J,bi,bj) = seaiceMaskV(I,J,bi,bj) * & ( seaiceMassV(I,J,bi,bj)*betaFacV & * vIce(I,J,bi,bj) & + seaiceMassV(I,J,bi,bj)*recip_deltaT*evpStarFac & * vIceNm1(I,J,bi,bj) & + FORCEY(I,J,bi,bj) & + stressDivergenceY(I,J,bi,bj) ) / denomV C-- to change results of lab_sea.hb87 test exp. (only preserve 2 digits for cg2d) c uIce(i,j,bi,bj) = uIceNm1(i,j,bi,bj) c & +( uIce(i,j,bi,bj) - uIceNm1(i,j,bi,bj) ) c vIce(i,j,bi,bj) = vIceNm1(i,j,bi,bj) c & +( vIce(i,j,bi,bj) - vIceNm1(i,j,bi,bj) ) ENDDO ENDDO #ifndef OBCS_UVICE_OLD DO j=1,sNy DO i=1,sNx locMaskU = maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj) locMaskV = maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) uIce(i,j,bi,bj) = uIce(i,j,bi,bj)*locMaskU & + uIceNm1(i,j,bi,bj)*(ONE-locMaskU) vIce(i,j,bi,bj) = vIce(i,j,bi,bj)*locMaskV & + vIceNm1(i,j,bi,bj)*(ONE-locMaskV) ENDDO ENDDO #endif /* OBCS_UVICE_OLD */ ENDDO ENDDO #ifdef ALLOW_AUTODIFF_TAMC CADJ STORE uIce = comlev1_evp, key = ikeyloc, byte = isbyte CADJ STORE vIce = comlev1_evp, key = ikeyloc, byte = isbyte #endif /* ALLOW_AUTODIFF_TAMC */ CALL EXCH_UV_XY_RL(uIce,vIce,.TRUE.,myThid) #ifdef ALLOW_SEAICE_EVP_RESIDUAL IF ( printResidual ) THEN DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) resTile(bi,bj) = 0. _d 0 DO J=1,sNy DO I=1,sNx uIcePm1(I,J,bi,bj) = seaiceMaskU(I,J,bi,bj) * & ( uIce(I,J,bi,bj)-uIcePm1(i,j,bi,bj) ) & * evpBetaU(I,J,bi,bj) vIcePm1(I,J,bi,bj) = seaiceMaskV(I,J,bi,bj) * & ( vIce(I,J,bi,bj)-vIcePm1(i,j,bi,bj) ) & * evpBetaV(I,J,bi,bj) ENDDO ENDDO IF ( .NOT. SEAICEscaleSurfStress ) THEN C multiply with mask (concentration) to only count ice contributions DO j=1,sNy DO i=1,sNx resTile(bi,bj) = resTile(bi,bj) + AREA(i,j,bi,bj) * & ( uIcePm1(I,J,bi,bj)*uIcePm1(I,J,bi,bj) & + vIcePm1(I,J,bi,bj)*vIcePm1(I,J,bi,bj) ) ENDDO ENDDO ELSE DO j=1,sNy DO i=1,sNx resTile(bi,bj) = resTile(bi,bj) & + uIcePm1(I,J,bi,bj)*uIcePm1(I,J,bi,bj) & + vIcePm1(I,J,bi,bj)*vIcePm1(I,J,bi,bj) ENDDO ENDDO ENDIF ENDDO ENDDO resloc = 0. _d 0 CALL GLOBAL_SUM_TILE_RL( resTile, resloc, myThid ) resloc = SQRT(resloc) WRITE(standardMessageUnit,'(A,1X,I6,1PE16.8)') & ' SEAICE_EVP: iEVPstep, residual U = ', iEVPstep, resLoc ENDIF CML WRITE(suff,'(I10.10)') myIter*100000+iEvpStep CML CALL WRITE_FLD_XY_RL( 'DELTA.',suff,deltaC, CML & myIter*100000+iEvpStep,myThid) CML CALL WRITE_FLD_XY_RL( 'RSIG1.',suff,sig11pm1, CML & myIter*100000+iEvpStep,myThid) CML CALL WRITE_FLD_XY_RL( 'RSIG2.',suff,sig22pm1, CML & myIter*100000+iEvpStep,myThid) CML CALL WRITE_FLD_XY_RL( 'RSIG12.',suff,sig12pm1, CML & myIter*100000+iEvpStep,myThid) CML CALL WRITE_FLD_XY_RL( 'RUICE.',suff,uIcePm1, CML & myIter*100000+iEvpStep,myThid) CML CALL WRITE_FLD_XY_RL( 'RVICE.',suff,vIcePm1, CML & myIter*100000+iEvpStep,myThid) #endif /* ALLOW_SEAICE_EVP_RESIDUAL */ ENDIF C end of the main time loop ENDDO #endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID and SEAICE_ALLOW_EVP */ RETURN END