! source file: /Users/mtivig/REPOS/uvic_with_news/RUNS/RUN_UVic_2015_NEWS_BD_V2020/updates/MTivig_2020_01/source/ice/evp.F subroutine evp (is, ie, js, je) !======================================================================= ! calculate velocities using an "elastic-viscous-plastic" rheology ! see E. C. Hunke and J. K. Dukowicz. An elastic-viscous-plastic ! model for sea ice dynamics. J. Phys. Oceanogr., 1997. !======================================================================= implicit none integer i, j, k, ie, je, is, js include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "ice.h" include "evp.h" do j=1,jmt do i=1,imt xint(i,j) = 0.0 yint(i,j) = 0.0 enddo enddo call mass_prss call viscevp call stressprep do k=1,ndte call stressevp call stepu call embmbc (uice) call embmbc (vice) enddo return end subroutine viscevp ! $Id: dyn.F,v 1.3 1997/02/11 18:04:27 eclare Exp $ !.. Elastic-viscous-plastic sea ice dynamics model !.. Computes ice velocity !.. author Elizabeth C. Hunke !.. Fluid Dynamics Group, Los Alamos National Laboratory !.. See E. C. Hunke and J. K. Dukowicz. An elastic-viscous-plastic model !.. for sea ice dynamics. J. Phys. Oceanogr., 1997. !.. Copyright, 1997. The Regents of the University of California. !.. This software was produced under a U.S. Government contract !.. (W-7405-ENG-36) by Los Alamos National Laboratory, which is operated !.. by the University of California for the U.S. Department of Energy. !.. The U.S. Government is licensed to use, reproduce, and distribute this !.. software. Permission is granted to the public to copy and use this !.. software without charge, provided that this Notice and any statement !.. of authorship are reproduced on all copies. Neither the Government !.. nor the University makes any warranty, express or implied, or assumes !.. any liability or responsibility for the use of this software. !======================================================================= !.. Computes the rates of strain, and the bulk and shear viscosities !.. zeta and eta. The rates of strain (xi*) and viscosities are !.. calculated for each of the four triangles in each grid cell !.. (north, south, east, west). !======================================================================= implicit none integer n, j, i real prs, pf, cc, dd, xi11n, xi12n, xi22n, xi11e, xi12e, xi22e real xi11s, xi12s, xi22s, xi11w, xi12w, xi22w, deltan, deltae real deltas, deltaw, zetamax include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "grdvar.h" include "atm.h" include "ice.h" include "evp.h" do n=1,nseg do j=jsi(n),jei(n)+1 do i=2,imt prs = 0.5*pice(i,j) ! g/s^2 (P/2, varies with c*H) pf = prs/floor ! initializes zeta to large values ! initialize zeta ! g/s zetan(i,j) = pf zetae(i,j) = pf zetas(i,j) = pf zetaw(i,j) = pf ! rates of strain ! 1/s cc = (uice(i,j) + uice(i-1,j) - uice(i,j-1) & - uice(i-1,j-1))*dyt2r(j) dd = (vice(i,j) + vice(i,j-1) - vice(i-1,j) & - vice(i-1,j-1))*cstr(j)*dxt2r(i) xi11n = (uice(i,j) - uice(i-1,j))*csur(j)*dxur(i) xi12n = ((vice(i,j) - vice(i-1,j))*csur(j)*dxur(i) & + cc)*0.5 xi22n = (vice(i,j) + vice(i-1,j) - vice(i,j-1) & - vice(i-1,j-1))*dyt2r(j) xi11e = (uice(i,j) + uice(i,j-1) - uice(i-1,j) & - uice(i-1,j-1))*cstr(j)*dxt2r(i) xi12e = ((uice(i,j) - uice(i,j-1))*dyur(j) + dd)*0.5 xi22e = (vice(i,j) - vice(i,j-1))*dyur(j) xi11s = (uice(i,j-1) - uice(i-1,j-1))*csur(j-1)*dxur(i) xi12s = ((vice(i,j-1) - vice(i-1,j-1))* & csur(j-1)*dxur(i) + cc)*0.5 xi22s = xi22n xi11w = xi11e xi12w = ((uice(i-1,j) - uice(i-1,j-1))*dyur(j) & + dd)*0.5 xi22w = (vice(i-1,j) - vice(i-1,j-1))*dyur(j) ! Delta (in the denominator of zeta, eta) ! 1/s deltan = sqrt( (xi11n**2 + xi22n**2)*ecc2p & + 4.0*xi12n**2*ecc2 + xi11n*xi22n*ecc2m) deltae = sqrt( (xi11e**2 + xi22e**2)*ecc2p & + 4.0*xi12e**2*ecc2 + xi11e*xi22e*ecc2m) deltas = sqrt( (xi11s**2 + xi22s**2)*ecc2p & + 4.0*xi12s**2*ecc2 + xi11s*xi22s*ecc2m) deltaw = sqrt( (xi11w**2 + xi22w**2)*ecc2p & + 4.0*xi12w**2*ecc2 + xi11w*xi22w*ecc2m) deltan = amax1(1.e-20,deltan) deltae = amax1(1.e-20,deltae) deltas = amax1(1.e-20,deltas) deltaw = amax1(1.e-20,deltaw) ! bulk viscosity zeta, bounded by maximum, minimum values zetamax = 2.5e8*pice(i,j) !Hibler pg 819 ! g/s zetan(i,j) = prs/deltan zetan(i,j) = min(zetamax,zetan(i,j)) zetan(i,j) = max(zetamin,zetan(i,j)) zetae(i,j) = prs/deltae zetae(i,j) = min(zetamax,zetae(i,j)) zetae(i,j) = max(zetamin,zetae(i,j)) zetas(i,j) = prs/deltas zetas(i,j) = min(zetamax,zetas(i,j)) zetas(i,j) = max(zetamin,zetas(i,j)) zetaw(i,j) = prs/deltaw zetaw(i,j) = min(zetamax,zetaw(i,j)) zetaw(i,j) = max(zetamin,zetaw(i,j)) ! mask zeta (and therefore eta) zetan(i,j) = zetan(i,j)*tmsk(i,j) zetae(i,j) = zetae(i,j)*tmsk(i,j) zetas(i,j) = zetas(i,j)*tmsk(i,j) zetaw(i,j) = zetaw(i,j)*tmsk(i,j) ! shear viscosity eta ! g/s etan(i,j) = zetan(i,j)*ecc2 etae(i,j) = zetae(i,j)*ecc2 etas(i,j) = zetas(i,j)*ecc2 etaw(i,j) = zetaw(i,j)*ecc2 enddo enddo enddo return end !----------------------------------------------------------------------- subroutine stressprep !======================================================================= !.. Computes quantities used in the subroutines stress and stepu. !.. This subroutine is cryptic - I apologize - but it made the code !.. faster by almost a factor of 2 on a Cray-YMP. Here we compute !.. quantities needed in the stress tensor (sigma) and momentum (u) !.. equations, but which don't change during the subcycling. Many of !.. these variables are grouped in parties of four, one for each !.. triangle (north, east, south, west). !======================================================================= implicit none integer n, j, i, nc real costh, sinth, econst, ey, e2, en, ee, es, ew, zn, ze, zs real zw, c2n, c2e, c2s, c2w, c3n, c3e, c3s, c3w, d2n, d2e real d2s, d2w, a2n, a2e, a2s, a2w, prs include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "cembm.h" include "levind.h" include "grdvar.h" include "atm.h" include "ice.h" include "evp.h" include "csbc.h" costh = 0.9063 sinth = 0.4226 ! used to compute E = Econst*c*H Econst = 2.*eyc*rhoice*xyminevp*dtei**2 ! g/cm s^2 do n=1,nseg do j=jsi(n),jei(n)+1 do i=2,imt if (tmsk(i,j) .gt. floor) then ey = Econst*hice(i,j,2) ! E, g/s^2 ey = max(ey,floor) e2 = 0.5*ey ! E/2 edy(i,j) = e2*dytr(j) ! E/(2*dy) edx(i,j) = e2*cstr(j)*dxtr(i) ! E/(2*dx) eHN(i,j) = e2/(csu(j)*dxu(i)) ! E/(2*HTN) eHE(i,j) = e2/dyu(j) ! etc eHNm(i,j) = e2/(csu(j-1)*dxu(i)) eHEm(i,j) = e2/dyu(j) en = e2/etan(i,j) ! E ee = e2/etae(i,j) ! ------- es = e2/etas(i,j) ! 2*eta ew = e2/etaw(i,j) zn = e2/zetan(i,j) ! E ze = e2/zetae(i,j) ! -------- zs = e2/zetas(i,j) ! 2*zeta zw = e2/zetaw(i,j) c2n = dtei + en ! 1 E c2e = dtei + ee ! --- + ----- c2s = dtei + es ! dte 2*eta c2w = dtei + ew c3n = 0.5*(en - zn) ! E 1 1 c3e = 0.5*(ee - ze) ! - * ( --- - ---- ) c3s = 0.5*(es - zs) ! 4 eta zeta c3w = 0.5*(ew - zw) d2n = c2n - c3n ! 1 E E 1 1 d2e = c2e - c3e ! --- + ----- - - * (--- - ----) d2s = c2s - c3s ! dte 2*eta 4 eta zeta d2w = c2w - c3w h2n(i,j) = 1./c2n ! this rapidly gets out of hand h2e(i,j) = 1./c2e h2s(i,j) = 1./c2s h2w(i,j) = 1./c2w a2n = h2n(i,j)/(d2n - c3n) a2e = h2e(i,j)/(d2e - c3e) a2s = h2s(i,j)/(d2s - c3s) a2w = h2w(i,j)/(d2w - c3w) b2n(i,j) = a2n*d2n b2e(i,j) = a2e*d2e b2s(i,j) = a2s*d2s b2w(i,j) = a2w*d2w a2na(i,j) = a2n*c3n a2ea(i,j) = a2e*c3e a2sa(i,j) = a2s*c3s a2wa(i,j) = a2w*c3w prs = 0.5*pice(i,j) prssn(i,j) = prs*zn ! P*E prsse(i,j) = prs*ze ! -------- prsss(i,j) = prs*zs ! 4*zeta prssw(i,j) = prs*zw endif enddo enddo enddo ! for subroutine stepu do n=1,nseg do j=jsi(n),jei(n)+1 do i=2,imt HTN4(i,j) = 0.25/(csu(j)*dxu(i)) HTE4(i,j) = 0.25/dyu(j) dxt8(i,j) = 0.125/(cst(j)*dxt(i)) dyt8(i,j) = 0.125/dyt(j) enddo enddo enddo do n=1,nseg do j=jsi(n),jei(n) do i=2,imtm1 fmass(i,j) = fcor(i,j)*umass(i,j) ! Coriolis * mass sinth = sign(sinth,fmass(i,j)) ! for water stress waterx(i,j) = umsk(i,j)*(sbc(i,j,igu)*costh - & sbc(i,j,igv)*sinth) watery(i,j) = umsk(i,j)*(sbc(i,j,igv)*costh + & sbc(i,j,igu)*sinth) ! air stress and -mg*gradH_o term (tilt) strairx(i,j) = umsk(i,j)*(sbc(i,j,itaux) - & fmass(i,j)*sbc(i,j,igv)) strairy(i,j) = umsk(i,j)*(sbc(i,j,itauy) + & fmass(i,j)*sbc(i,j,igu)) enddo enddo enddo return end !----------------------------------------------------------------------- subroutine stressevp !======================================================================= !.. Calculates the internal stress components, sigma_ij, in the four !.. triangles of each cell. !======================================================================= implicit none integer n, j, i real dun, dus, due, duw, dvn, dvs, dve, dvw, cc, dd, xi11n, xi12n real xi22n, xi11e, xi12e, xi22e, xi11s, xi12s, xi22s, xi11w real xi12w, xi22w, c4n, c4e, c4s, c4w, c5n, c5e, c5s, c5w include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "levind.h" include "atm.h" include "ice.h" include "evp.h" do n=1,nseg do j=jsi(n), jei(n) + 1 do i=2,imt if (tmsk(i,j) .gt. floor) then dun = uice(i,j) - uice(i-1,j) dus = uice(i,j-1) - uice(i-1,j-1) due = uice(i,j) - uice(i,j-1) duw = uice(i-1,j) - uice(i-1,j-1) dvn = vice(i,j) - vice(i-1,j) dvs = vice(i,j-1) - vice(i-1,j-1) dve = vice(i,j) - vice(i,j-1) dvw = vice(i-1,j) - vice(i-1,j-1) cc = 0.5*edy(i,j)*(due + duw) dd = 0.5*edx(i,j)*(dvn + dvs) ! NOTE these are rates of strain * E xi11n = 2.0*dun*eHN(i,j) xi12n = dvn*eHN(i,j) + cc xi22n = edy(i,j)*(dve + dvw) xi11e = edx(i,j)*(dun + dus) xi12e = due*eHE(i,j) + dd xi22e = 2.0*dve*eHE(i,j) xi11s = 2.0*dus*eHNm(i,j) xi12s = dvs*eHNm(i,j) + cc xi22s = xi22n xi11w = xi11e xi12w = duw*eHEm(i,j) + dd xi22w = 2.0*dvw*eHEm(i,j) ! solve for the three components of sigma in each triangle c4n = dtei*sig11n(i,j) + xi11n - prssn(i,j) c4e = dtei*sig11e(i,j) + xi11e - prsse(i,j) c4s = dtei*sig11s(i,j) + xi11s - prsss(i,j) c4w = dtei*sig11w(i,j) + xi11w - prssw(i,j) c5n = dtei*sig22n(i,j) + xi22n - prssn(i,j) c5e = dtei*sig22e(i,j) + xi22e - prsse(i,j) c5s = dtei*sig22s(i,j) + xi22s - prsss(i,j) c5w = dtei*sig22w(i,j) + xi22w - prssw(i,j) sig11n(i,j) = a2na(i,j)*c5n + c4n*b2n(i,j) ! g/s^2 sig11e(i,j) = a2ea(i,j)*c5e + c4e*b2e(i,j) sig11s(i,j) = a2sa(i,j)*c5s + c4s*b2s(i,j) sig11w(i,j) = a2wa(i,j)*c5w + c4w*b2w(i,j) sig22n(i,j) = a2na(i,j)*c4n + c5n*b2n(i,j) sig22e(i,j) = a2ea(i,j)*c4e + c5e*b2e(i,j) sig22s(i,j) = a2sa(i,j)*c4s + c5s*b2s(i,j) sig22w(i,j) = a2wa(i,j)*c4w + c5w*b2w(i,j) sig12n(i,j) = h2n(i,j)*(xi12n + dtei*sig12n(i,j)) sig12e(i,j) = h2e(i,j)*(xi12e + dtei*sig12e(i,j)) sig12s(i,j) = h2s(i,j)*(xi12s + dtei*sig12s(i,j)) sig12w(i,j) = h2w(i,j)*(xi12w + dtei*sig12w(i,j)) endif enddo enddo enddo return end !----------------------------------------------------------------------- subroutine mass_prss !======================================================================= !.. Computes ice mass and pressure (strength) !======================================================================= implicit none integer n, j, i, nc real pstar, volm, area include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "cembm.h" include "atm.h" include "ice.h" include "evp.h" real tmass(imt,jmt) pstar = 2.75e5 ! total mass of ice and snow, centred in T-cell do n=1,nseg do j=jsi(n), jei(n) + 1 do i=1,imt if (tmsk(i,j) .ne. 0.0) then ! if you, like Hibler (1979), have only thick ice tmass(i,j) = rhoice*hice(i,j,2) ! g/cm^2 else tmass(i,j) = 0.0 endif enddo enddo enddo ! mass centred at velocity nodes (U-cells) do n=1,nseg do j=jsi(n),jei(n) do i=2,imtm1 umass(i,j) = 0.25*(tmass(i,j) + tmass(i+1,j) & + tmass(i,j+1) + tmass(i+1,j+1)) ! g/cm^2 enddo enddo enddo ! pressure P do n=1,nseg do j=jsi(n),jei(n)+1 do i=2,imt pice(i,j) = pstar*hice(i,j,2) & *exp(-20.0*(1.0 - aice(i,j,2))) ! g/s^2 enddo enddo enddo call embmbc (pice) return end !----------------------------------------------------------------------- subroutine stepu !======================================================================= !.. Calculation of the surface stresses !.. Integration of the momentum equation to find velocity (u,v) !======================================================================= implicit none integer n, j, i real costh, sinth, uorel, vorel, vrel, umassdtei, cca, ccb, ab2 real s11, s12, s21, s22, c1evp, c2evp include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "levind.h" include "atm.h" include "ice.h" include "evp.h" include "csbc.h" real s11ns(imt,jmt), s11ew(imt,jmt) real s22ew(imt,jmt), s22ns(imt,jmt) real s12ewi(imt,jmt), s12ns(imt,jmt) real s12nsj(imt,jmt), s12ew(imt,jmt) costh = 0.9063 sinth = 0.4226 ! some cryptic but useful arrays do n=1,nseg do j=jsi(n), jei(n) + 1 do i=2,imt s11ew(i,j) = dxt8(i,j)*(sig11e(i,j) + sig11w(i,j)) s22ns(i,j) = dyt8(i,j)*(sig22n(i,j) + sig22s(i,j)) s12ns(i,j) = dyt8(i,j)*(sig12n(i,j) + sig12s(i,j)) s12ew(i,j) = dxt8(i,j)*(sig12e(i,j) + sig12w(i,j)) enddo enddo enddo do n=1,nseg do j=jsi(n), jei(n) + 1 do i=2,imtm1 s22ew(i,j) = HTE4(i,j)*(sig22e(i,j) + sig22w(i+1,j)) s12ewi(i,j) = HTE4(i,j)*(sig12e(i,j) + sig12w(i+1,j)) enddo enddo enddo do n=1,nseg do j=jsi(n),jei(n) do i=2,imt s11ns(i,j) = HTN4(i,j)*(sig11s(i,j+1) + sig11n(i,j)) s12nsj(i,j) = HTN4(i,j)*(sig12s(i,j+1) + sig12n(i,j)) enddo enddo enddo ! integrate the momentum equation do n=1,nseg do j=jsi(n),jei(n) do i=2,imtm1 if (umsk(i,j) .gt. floor .and. umass(i,j) .gt. 0.01) then ! geostrophic ocean currents relative to ice velocity uorel = sbc(i,j,igu) - uice(i,j) vorel = sbc(i,j,igv) - vice(i,j) ! (magnitude of relative geostrophic ocean current)*rhow*drag vrel = 0.0055*1.03*sqrt(uorel**2 + vorel**2) ! cm/s ! vrel = dragw*sqrt(uorel**2 + vorel**2) ! cm/s ! alpha, beta are defined in Hunke and Dukowicz sec 3.2 umassdtei = umass(i,j)*dtei ! m/dte,alpha, beta, g/cm^2 s cca = umassdtei + vrel*costh ccb = fmass(i,j) + vrel*sign(sinth,fmass(i,j)) ab2 = cca**2 + ccb**2 ! more cryptic stuff s11 = - s11ns(i,j) + s11ns(i+1,j) + s11ew(i+1,j+1) & + s11ew(i+1,j) - s11ew(i,j+1) - s11ew(i,j) s12 = - s12ewi(i,j) + s12ewi(i,j+1) + s12ns(i+1,j+1) & + s12ns(i,j+1) - s12ns(i+1,j) - s12ns(i,j) s21 = - s12nsj(i,j) + s12nsj(i+1,j) + s12ew(i+1,j+1) & + s12ew(i+1,j) - s12ew(i,j+1) - s12ew(i,j) s22 = - s22ew(i,j) + s22ew(i,j+1) + s22ns(i+1,j+1) & + s22ns(i,j+1) - s22ns(i+1,j) - s22ns(i,j) xint(i,j) = s11 + s12 yint(i,j) = s21 + s22 ! finally, the velocity components (in cm/s) c1evp = xint(i,j) + strairx(i,j) + vrel*waterx(i,j) & + umassdtei*uice(i,j) c2evp = yint(i,j) + strairy(i,j) + vrel*watery(i,j) & + umassdtei*vice(i,j) uice(i,j) = (cca*c1evp + ccb*c2evp)/ab2 vice(i,j) = (cca*c2evp - ccb*c1evp)/ab2 else ! set velocity to zero on land and (nearly) open water uice(i,j) = 0. vice(i,j) = 0. endif enddo enddo enddo return end subroutine strain ! B-grid strain rate tensor for EVP model (lacks metric ! terms, see subroutine strain in vpadi.F for full calc) implicit none integer n, j, i include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "ice.h" include "evp.h" include "grdvar.h" real e11(imt,jmt),e22(imt,jmt),e12(imt,jmt) do n=1,nseg do j=jsi(n),jei(n) do i=2,imtm1 e11(i,j) = 0.5*cstr(j)*dxtr(i)*( uice(i,j) + uice(i,j-1) & - uice(i-1,j) - uice(i-1,j-1) ) e22(i,j) = 0.5*cstr(j)*dytr(j)*( & csu(j)*( vice(i,j) + vice(i-1,j) ) & - csu(j-1)*( vice(i,j-1) + vice(i-1,j-1) ) ) e12(i,j) = 0.25*cstr(j)*( dytr(j)*( & csu(j)*( uice(i,j) + uice(i-1,j) ) & - csu(j-1)*( uice(i,j-1) - uice(i-1,j-1) ) ) & + dxtr(i)*( vice(i,j) + vice(i,j-1) & - vice(i-1,j) - vice(i-1,j-1) ) ) eI(i,j)=e11(i,j)+e22(i,j) del(i,j) = (e11(i,j)**2 + e22(i,j)**2)*(1.0 + ecc2) & + 4.0*ecc2*e12(i,j)**2 + 2.0*e11(i,j)*e22(i,j)* & (1.0 - ecc2) del(i,j) = sqrt(del(i,j)) enddo enddo enddo call embmbc (del) call embmbc (eI) return end