!!****f* HYCOM_2.1.03/EVP_MPI/evp_stress !! !! NAME !! stress !! !! SYNOPSIS !! !! !! !! DESCRIPTION !! Computes strain rates and internal stress components. !! Computes the rates of strain and internal stress components for !! each of the four corners on each T-grid cell !! !! !! AUTHOR !! Elizabeth C. Hunke, Fluid Dynamics Group, Los Alamos National Laboratory !! !! HISTORY !! Oct 2006 - Modified to use HYCOM tiling logic - Knut Liseter !! 12.03.2007 - Added OMP directives !! 28.07.2009 - Corrected a few bugs and cleaned the code - Dany Dumont !! !! INPUT !! ksub -- Subcycling time step !! SIDE EFFECTS !! Sets up the following variables for evp_stepu !! stressm_[1-4] !! stressp_[1-4] !! stress12_[1-4] !! Keeps the following variables for diagnostics and mechanical redist !! Divu !! Delta !! shear !! Because a grid difference is involved, the effective margin is !! reduced by one !! !! SOURCE !! subroutine evp_stress_miz_2(ksub) use mod_evp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! A few options here: !! (1) with WAVES & with/without ICE_DYN_DIAG #if defined (WAVES) !use mod_common_ice, only : dfloe,dfloe_miz,MIZ_MASK use mod_common_wavesice, only : dfloe,dfloe_miz #if defined(ICE_DYN_DIAG) use mod_common_ice, only : MIZ_MASK,strainI,strainII #else use mod_common_ice, only : MIZ_MASK #endif !! (2) without WAVES & with ICE_DYN_DIAG #elif defined(ICE_DYN_DIAG) /*no WAVES, yes ICE_DYN_DIAG*/ use mod_common_ice, only : MIZ_MASK,strainI,strainII !! (3) without WAVES & without ICE_DYN_DIAG #else /*no WAVES, no ICE_DYN_DIAG*/ use mod_common_ice, only : MIZ_MASK #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! implicit none integer, intent(in) :: ksub ! subcycling step integer :: i, j real :: & divune, divunw, divuse, divusw ! divergence &, tensionne, tensionnw, tensionse, tensionsw ! tension &, shearne, shearnw, shearse, shearsw ! shearing &, msq_tension,msq_shear,divu_mean ! mean square tension, shear; mean of divu &, Deltane, Deltanw, Deltase, Deltasw ! Delta &, c0ne, c0nw, c0se, c0sw ! useful combinations &, c1ne, c1nw, c1se, c1sw &, detne, detnw, detse, detsw &, K2ne, K2nw, K2se, K2sw &, K3ne, K3nw, K3se, K3sw &, vdne, vdnw, vdse, vdsw &, etane, etanw, etase, etasw &, zetane, zetanw, zetase, zetasw &, prssne, prssnw, prssse, prsssw %, Eyoung, fac1, fac2, fac, pi, gamp, sep, hfloe,wec,wecci &, maxstressm, maxstressp,maxstress12 &, minstressm, minstressp,minstress12 &, minvpD, maxvpD, cice, minvp integer :: & ij ! loop index, combination of i and j loops logical, parameter :: diag=.true. c --- From the collisional rheology model of Shen et al. (1987) real, parameter :: e=0.1 ! restitution coefficient real, parameter :: diam=30. ! diameter of floes real, parameter :: C0=1.1 ! max compactness of hexagonal packing real, parameter :: K1=(1-e**2)*0.5 ! factor fo ert. velocity assumption real, parameter :: rhoi=930. ! density of sea ice real, parameter :: cs=0.89 ! cutoff concentration to turn on MIZ real, parameter :: hmin=1.0 ! minimum ice thickness for VP at c=1 real, parameter :: hmax=2.0 ! maximum MIZ thickness at any concentration c --- MIZ viscosity ++ real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & mizzeta &, mizeta &, mizprss &, mizvd !include 'common_blocks.h'!![TW 2013.03.10]: want itest/jtest pi = 4.d0*atan(1.d0) maxstressm = -1e14 maxstressp = -1e14 maxstress12 = -1e14 maxvpD = -1e14 minstressm = 1e14 minstressp = 1e14 minstress12 = 1e14 minvpD = 1e14 !$OMP PARALLEL DO PRIVATE(j,i, !$OMP& divune, divunw, divusw, divuse, !$OMP& tensionne, tensionnw,tensionsw,tensionse, !$OMP& shearne, shearnw, shearsw, shearse, !$OMP& Deltane, Deltanw, Deltasw, Deltase, !$OMP& c0ne, c0nw, c0sw, c0se, !$OMP& c1ne, c1nw, c1sw, c1se, !$OMP& detne, detnw, detsw, detse, !$OMP& K2ne, K2nw, K2sw, K2se, !$OMP& K3ne, K3nw, K3sw, K3se, !$OMP& vdne, vdnw, vdse, vdsw, !$OMP& etane, etanw, etase, etasw, !$OMP& zetane, zetanw, zetase, zetasw, !$OMP& prssne, prssnw, prssse, prsssw, !$OMP& fac , dfloe, Eyoung, fac1, fac2,cice) !$OMP& SCHEDULE(STATIC,jblk) MIZ_MASK = 0.0!!test variable stored in mod_common_ice.F msq_tension = 0.0 msq_shear = 0.0 divu_mean = 0.0 do j=1-imargin,jj+imargin do i=1-imargin,ii+imargin if (icetmask(i,j)) then !----------------------------------------------------------------- ! strain rates ! NOTE these are actually strain rates * area (m^2/s) !----------------------------------------------------------------- ! divergence = e_11 + e_22 divune = cyp(i,j)*uvel(i ,j ) - dyt(i,j)*uvel(i-1,j ) & + cxp(i,j)*vvel(i ,j ) - dxt(i,j)*vvel(i ,j-1) divunw = cym(i,j)*uvel(i-1,j ) + dyt(i,j)*uvel(i ,j ) & + cxp(i,j)*vvel(i-1,j ) - dxt(i,j)*vvel(i-1,j-1) divusw = cym(i,j)*uvel(i-1,j-1) + dyt(i,j)*uvel(i ,j-1) & + cxm(i,j)*vvel(i-1,j-1) + dxt(i,j)*vvel(i-1,j ) divuse = cyp(i,j)*uvel(i ,j-1) - dyt(i,j)*uvel(i-1,j-1) & + cxm(i,j)*vvel(i ,j-1) + dxt(i,j)*vvel(i ,j ) ! tension strain rate = e_11 - e_22 tensionne = -cym(i,j)*uvel(i ,j ) - dyt(i,j)*uvel(i-1,j ) & + cxm(i,j)*vvel(i ,j ) + dxt(i,j)*vvel(i ,j-1) tensionnw = -cyp(i,j)*uvel(i-1,j ) + dyt(i,j)*uvel(i ,j ) & + cxm(i,j)*vvel(i-1,j ) + dxt(i,j)*vvel(i-1,j-1) tensionsw = -cyp(i,j)*uvel(i-1,j-1) + dyt(i,j)*uvel(i ,j-1) & + cxp(i,j)*vvel(i-1,j-1) - dxt(i,j)*vvel(i-1,j ) tensionse = -cym(i,j)*uvel(i ,j-1) - dyt(i,j)*uvel(i-1,j-1) & + cxp(i,j)*vvel(i ,j-1) - dxt(i,j)*vvel(i ,j ) ! shearing strain rate = 2*e_12 shearne = -cym(i,j)*vvel(i ,j ) - dyt(i,j)*vvel(i-1,j ) & - cxm(i,j)*uvel(i ,j ) - dxt(i,j)*uvel(i ,j-1) shearnw = -cyp(i,j)*vvel(i-1,j ) + dyt(i,j)*vvel(i ,j ) & - cxm(i,j)*uvel(i-1,j ) - dxt(i,j)*uvel(i-1,j-1) shearsw = -cyp(i,j)*vvel(i-1,j-1) + dyt(i,j)*vvel(i ,j-1) & - cxp(i,j)*uvel(i-1,j-1) + dxt(i,j)*uvel(i-1,j ) shearse = -cym(i,j)*vvel(i ,j-1) - dyt(i,j)*vvel(i-1,j-1) & - cxp(i,j)*uvel(i ,j-1) + dxt(i,j)*uvel(i ,j ) ! Define and store strain rate invariants for diagnostics #if defined(ICE_DYN_DIAG) !!TW change: 20140307 !!averaging of corners done weirdly !! if (ksub.eq.ndte) then !! strainI(i,j) = 0.5*0.25*tarear(i,j)*(divune + divunw !! & + divuse + divusw) !! strainII(i,j) = 0.5*0.25*tarear(i,j)*sqrt( !! & (tensionne + tensionnw + tensionse + tensionsw)**2 !! & + ( shearne + shearnw + shearse + shearsw)**2) !! endif if (ksub.eq.ndte) then divu_mean = .25*(divune + divunw & + divuse + divusw) msq_tension = .25*(tensionne**2+tensionnw**2 & +tensionse**2+tensionsw**2); msq_shear = .25*(shearne**2+shearnw**2 & +shearse**2+shearsw**2); !!principal strains are strainI +/- strainII !!NB tension=e11-e22, shear=2*e12 strainI(i,j) = 0.5*tarear(i,j)*divu_mean strainII(i,j) = 0.5*tarear(i,j)*sqrt(msq_tension+msq_shear) endif #endif ! Ice concentration cice = max(aice(i,j),0.01) ! always > 0.01 hfloe = max(vice(i,j)/cice,0.10) ! always > 0.10 !! if (mnproc==1) then !! print*,'TWt2: c,h',cice,hfloe !! end if #if defined (WAVES) ! VERSION 3 : with explicit inclusion of waves ! We use the floe size to alter the major to minor axis ratio (shear strength) ! NB max value is wec=2 ! - this is the EVP value (ecc=4) ! NB wec -> sqrt(ecc) in EVP ! wecci = 1/wec/wec -> 1/ecc = ecci in EVP ! in evp_init.F, evp_stress.F ! dfloe_miz=200.0 (defined in mod_common_ice.F) !! Formula is completely arbitrary !! wec = max(1.8+10*exp(-dfloe(i,j)/dfloe_miz),2.0)!![TW/JB change 20140311] 200m is too big wec = max(1.8+10*exp(-dfloe(i,j)/40),2.0) wecci = 1./(wec*wec) #else ! Version without WAVES !! Formula is completely arbitrary wec = max(1.8+10*exp(-(cice/cs+hfloe/hmax)/4.0),2.0) wecci = 1./(wec*wec) #endif ! Viscous-plastic regime computed with the EVP scheme ! Delta (in the denominator of zeta, eta) Deltane = sqrt(divune**2 + wecci*(tensionne**2 + shearne**2)) Deltanw = sqrt(divunw**2 + wecci*(tensionnw**2 + shearnw**2)) Deltase = sqrt(divuse**2 + wecci*(tensionse**2 + shearse**2)) Deltasw = sqrt(divusw**2 + wecci*(tensionsw**2 + shearsw**2)) if (evp_damping) then ! enforce damping criterion c0ne=min(prss(i,j)/max(Deltane,4.*tinyarea(i,j)),rcon_evp(i,j)) c0nw=min(prss(i,j)/max(Deltanw,4.*tinyarea(i,j)),rcon_evp(i,j)) c0sw=min(prss(i,j)/max(Deltasw,4.*tinyarea(i,j)),rcon_evp(i,j)) c0se=min(prss(i,j)/max(Deltase,4.*tinyarea(i,j)),rcon_evp(i,j)) prs_sig(i,j) = prss(i,j)*Deltane/max(Deltane,4.*tinyarea(i,j)) ! ne else ! original version c0ne = prss(i,j)/max(Deltane,tinyarea(i,j)) c0nw = prss(i,j)/max(Deltanw,tinyarea(i,j)) c0sw = prss(i,j)/max(Deltasw,tinyarea(i,j)) c0se = prss(i,j)/max(Deltase,tinyarea(i,j)) prs_sig(i,j) = c0ne*Deltane ! northeast endif c1ne = c0ne*dte2T c1nw = c0nw*dte2T c1sw = c0sw*dte2T c1se = c0se*dte2T !----------------------------------------------------------------- ! the stresses ! kg/s^2 ! (1) northeast, (2) northwest, (3) southwest, (4) southeast !----------------------------------------------------------------- #if defined (EVP_CHG) !!bug found by Sylvain to prevent features aligned with grid stressp_1(i,j) = (stressp_1(i,j) + c1ne*(divune - Deltane))*denom stressp_2(i,j) = (stressp_2(i,j) + c1nw*(divunw - Deltanw))*denom stressp_3(i,j) = (stressp_3(i,j) + c1sw*(divusw - Deltasw))*denom stressp_4(i,j) = (stressp_4(i,j) + c1se*(divuse - Deltase))*denom stressm_1(i,j) = (stressm_1(i,j) + c1ne*wecci*tensionne)*denom stressm_2(i,j) = (stressm_2(i,j) + c1nw*wecci*tensionnw)*denom stressm_3(i,j) = (stressm_3(i,j) + c1sw*wecci*tensionsw)*denom stressm_4(i,j) = (stressm_4(i,j) + c1se*wecci*tensionse)*denom stress12_1(i,j) = (stress12_1(i,j) + c1ne*wecci*shearne*.5)*denom stress12_2(i,j) = (stress12_2(i,j) + c1nw*wecci*shearnw*.5)*denom stress12_3(i,j) = (stress12_3(i,j) + c1sw*wecci*shearsw*.5)*denom stress12_4(i,j) = (stress12_4(i,j) + c1se*wecci*shearse*.5)*denom #else stressp_1(i,j) = (stressp_1(i,j) + c1ne*(divune - Deltane))*denom1 stressp_2(i,j) = (stressp_2(i,j) + c1nw*(divunw - Deltanw))*denom1 stressp_3(i,j) = (stressp_3(i,j) + c1sw*(divusw - Deltasw))*denom1 stressp_4(i,j) = (stressp_4(i,j) + c1se*(divuse - Deltase))*denom1 stressm_1(i,j) = (stressm_1(i,j) + c1ne*tensionne)*denom2 stressm_2(i,j) = (stressm_2(i,j) + c1nw*tensionnw)*denom2 stressm_3(i,j) = (stressm_3(i,j) + c1sw*tensionsw)*denom2 stressm_4(i,j) = (stressm_4(i,j) + c1se*tensionse)*denom2 stress12_1(i,j) = (stress12_1(i,j) + c1ne*shearne*.5)*denom2 stress12_2(i,j) = (stress12_2(i,j) + c1nw*shearnw*.5)*denom2 stress12_3(i,j) = (stress12_3(i,j) + c1sw*shearsw*.5)*denom2 stress12_4(i,j) = (stress12_4(i,j) + c1se*shearse*.5)*denom2 #endif endif!icetmask(i,j).eq.(.true.) enddo enddo !$OMP END PARALLEL DO ! if (diag) then ! call xcmaxr(maxstressm) ! call xcminr(minstressm) ! if (mnproc==1) print *,'min/max stressm:',minstressm,maxstressm ! call xcmaxr(maxstressp) ! call xcminr(minstressp) ! if (mnproc==1) print *,'min/max stressp:',minstressp,maxstressp ! call xcmaxr(maxstress12) ! call xcminr(minstress12) ! if (mnproc==1) print *,'min/max stress12:',minstress12,maxstress12 ! call xcmaxr(maxvpD) ! call xcminr(minvpD) ! if (mnproc==1) print *,'min/max vpD:',minvpD,maxvpD ! end if ! !call xcstop('(evp_stress_miz)') end subroutine evp_stress_miz_2 !!******