! source file: /Users/csomes/Research/Models/UVic_ESCM/2.9/sim/nonredfield_dom/clean_up_code/seaice_wind_bug/afs_high/updates/isopyc.F subroutine isopi (error, am, ah) !======================================================================= ! Initialization for isopycnal mixing scheme ! -Disopycmix gives either the full or small angle isopycnal ! mixing tensor ! -Disopycmix -Dgent_mcwilliams gives the isopycnal mixing tensor ! plus the advective velocity parameterization of Gent_McWilliams ! input: ! error = logical to signal problems ! am = background mixing coeff for momentum ! ah = horizontal mixing coeff for tracers ! slmx = max slope of isopycnals ! ahisop = isopycnal tracer diffusivity(cm**2/sec) ! athkdf = isopycnal thickness diffusivity (cm**2/sec) ! output: ! The above input can be reset via namelist !======================================================================= implicit none character(120) :: fname, new_file_name integer i, k, j, ip, kr, jq, io, i_delta, j_delta, k_delta, jrow integer iou, n, ib(10), ic(10) logical exists, inqvardef, error real slmx, s_dm, ah, am, ft, delta1, delta2 include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "accel.h" include "coord.h" include "grdvar.h" include "iounit.h" include "isopyc.h" include "levind.h" include "mw.h" include "scalar.h" include "switch.h" include "vmixc.h" real tmpijk(imtm2,jmtm2,km) namelist /isopyc/ slmx, ahisop, athkdf, del_dm, s_dm write (stdout,'(/,20x,a,/,20x,a,/)') & 'I S O P Y C M I X I N I T I A L I Z A T I O N' &,'Isopycnal mixing tensor + GM thickness diffusion' !----------------------------------------------------------------------- ! USER INPUT ! initialize variables (all mixing units are cm**2/sec.) !----------------------------------------------------------------------- ! maximum isopycnal slope slmx = 1.0/100.0 ! isopycnal diffusion coefficient ahisop = 1.e7 ! isopycnal thickness diffusion coefficient athkdf = 1.0e7 ! transition for scaling diffusion coefficient del_dm = 4.0/1000.0 ! half width scaling for diffusion coefficient s_dm = 1.0/1000.0 !----------------------------------------------------------------------- ! provide for namelist over-ride of above settings + documentation !----------------------------------------------------------------------- call getunit (io, 'control.in' &, 'formatted sequential rewind') read (io,isopyc,end=100) 100 continue write (stdout,isopyc) call relunit (io) ! set reciprocal of maximum isopycnal slope and scaling half width slmxr = c1/slmx s_dmr = c1/s_dm !----------------------------------------------------------------------- ! check for problems !----------------------------------------------------------------------- write (stdout,*) & '==> Note: small angle approximation to isopycnal tensor is' &,' being used. ' write (stdout,*) & '==> Note: Re-scaling of Gerdes et al is being used to ' &,' reduce mixing coeff within areas of steeply ' &,' sloping isopycnals ' !----------------------------------------------------------------------- ! print out tracer diffusion coefficients in isopycnal case !----------------------------------------------------------------------- write(stdout,9102) ahisop, athkdf, ah, am 9102 format(/' ahisop = ',e12.6,' along isopyncal tracer mixing ' &, '(cm**2/sec) '/' athkdf = ',e12.6,' isopycnal thickness ' &,'diffusion (cm**2/sec) '/' ah = ',e12.6,' cm**2/sec for ' &,' am = ',e12.6,' cm**2/sec for horizontal mixing of momentum'/) ft = c1/(4.0*ahisop*dtts) delta_iso = dxt(1)*cst(jmt/2)*dzt(1)*ft i_delta = 1 j_delta = 1 k_delta = 1 do jrow=2,jmt-1 do i=2,imt-1 do k=1,km delta1 = dxt(i)*cst(jrow)*dzt(k)*ft delta2 = dyt(jrow)*dzt(k)*ft if (delta_iso .ge. delta1 .or. delta_iso .ge. delta2) then i_delta = i j_delta = jrow k_delta = k delta_iso = min(delta1,delta2) endif enddo enddo enddo if (delta_iso .lt. p5) then s_minus = (c1 - sqrt(c1 - 4.0*delta_iso**2))/(c2*delta_iso) s_plus = c1/s_minus else s_minus = c0 s_plus = c0 endif write(stdout,'(a)') &'------------------------------------------------' write(stdout,'(a,e14.7)') &'The isopycnal diffusion grid factor delta_iso =',delta_iso write(stdout,'(a)') &'was determined at the grid point' write(stdout,'(a,i4,a,e14.7)') &'dxt(',i_delta,') = ',dxt(i_delta) write(stdout,'(a,i4,a,e14.7)') &'dyt(',j_delta,') = ',dyt(j_delta) write(stdout,'(a,i4,a,e14.7)') &'dzt(',k_delta,') = ',dzt(k_delta) write(stdout,'(a,i4,a,e14.7)') &'cst(',j_delta,') = ',cst(j_delta) write(stdout,'(a,e14.7)') &'dtts = ',dtts write(stdout,'(a,e14.7)') &'ahisop = ',ahisop write(stdout,'(/a,e14.7/)') & 'Maximum allowable isopycnal slope is specified as slmx = ',slmx write(stdout,'(a)') &'------------------------------------------------' !----------------------------------------------------------------------- ! store the square root of the tracer timestep acceleration values ! into variable "dtxsqr" for use in isopycnal mixing !----------------------------------------------------------------------- do k=1,km dtxsqr(k) = sqrt(dtxcel(k)) enddo write (stdout,'(a)') ' Acceleration with depth "dtxsqr(k)"=' write (stdout,'(5(1x,e12.6))') (dtxsqr(k),k=1,km) !----------------------------------------------------------------------- ! read ocean diffusion factor !----------------------------------------------------------------------- fisop(:,:,:) = 1. fname = new_file_name ("O_diffac.nc") inquire (file=trim(fname), exist=exists) if (exists) then ib(:) = 1 ic(:) = 1 ic(1) = imtm2 ic(2) = jmtm2 ic(3) = km call openfile (fname, iou) exists = inqvardef('O_diffac', iou) if (exists) then write(*,*) "reading ocean diffusion factor from: ",fname call getvara ('O_diffac', iou, imtm2*jmtm2*km, ib, ic, tmpijk &, c1, c0) fisop(2:imtm1,2:jmtm1,1:km) = tmpijk(1:imtm2,1:jmtm2,1:km) endif endif !----------------------------------------------------------------------- ! initialize arrays !----------------------------------------------------------------------- do j=1,jmw do k=1,km do i=1,imt alphai(i,k,j) = c0 betai(i,k,j) = c0 enddo enddo enddo do j=jsmw,jemw do n=1,2 do k=1,km do i=1,imt ddxt(i,k,j,n) = c0 enddo enddo enddo enddo do j=1,jemw do n=1,2 do k=1,km do i=1,imt ddyt(i,k,j,n) = c0 enddo enddo enddo enddo do j=1,jmw do n=1,2 do k=0,km do i=1,imt ddzt(i,k,j,n) = c0 enddo enddo enddo enddo do j=jsmw,jemw do k=1,km do i=1,imt do kr=0,1 do ip=0,1 Ai_ez(i,k,j,ip,kr) = c0 enddo enddo do kr=0,1 do ip=0,1 Ai_bx(i,k,j,ip,kr) = c0 enddo enddo do kr=0,1 do jq=0,1 Ai_by(i,k,j,jq,kr) = c0 enddo enddo K33(i,k,j) = c0 K11(i,k,j) = c0 adv_vetiso(i,k,j) = c0 enddo enddo enddo do j=1,jemw do k=1,km do i=1,imt do kr=0,1 do jq=0,1 Ai_nz(i,k,j,jq,kr) = c0 enddo enddo K22(i,k,j) = c0 adv_vntiso(i,k,j) = c0 enddo enddo enddo do j=jsmw,jemw do k=0,km do i=1,imt adv_vbtiso(i,k,j) = c0 adv_fbiso(i,k,j) = c0 enddo enddo enddo return end subroutine elements (joff, js, je, is, ie) !======================================================================= ! Estimate alpha, beta, and normal gradients on faces of T cells !======================================================================= implicit none integer i, k, j, ip, kr, jq, js, je, is, ie, n, kp1, jrow, joff real dens, tq, sq, drodt, drods, tprime, sprime include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "accel.h" include "grdvar.h" include "levind.h" include "mw.h" include "state.h" include "isopyc.h" include "dens.h" !----------------------------------------------------------------------- ! alpha and beta at centers of T cells !----------------------------------------------------------------------- do j=js,je do k=1,km do i=is,ie tprime = t(i,k,j,1,taum1)-to(k) sprime = t(i,k,j,2,taum1)-so(k) alphai(i,k,j) = drodt (tprime, sprime, k) betai(i,k,j) = drods (tprime, sprime, k) enddo enddo call setbcx (alphai(1,1,j), imt, km) call setbcx (betai(1,1,j), imt, km) enddo !----------------------------------------------------------------------- ! gradients at bottom face of T cells !----------------------------------------------------------------------- do j=js,je do n=1,2 do k=1,km kp1 = min(k+1,km) do i=is,ie ddzt(i,k,j,n) = tmask(i,kp1,j)*dzwr(k)* & (t(i,k,j,n,taum1) - t(i,kp1,j,n,taum1)) enddo enddo do i=is,ie ddzt(i,0,j,n) = c0 enddo call setbcx (ddzt(1,0,j,n), imt, km+1) enddo enddo !----------------------------------------------------------------------- ! gradients at eastern face of T cells !----------------------------------------------------------------------- do j=max(js-1,2), je-1 jrow = j + joff do n=1,2 do k=1,km do i=is,ie ddxt(i,k,j,n) = tmask(i,k,j)*tmask(i+1,k,j)*cstr(jrow)* & dxur(i)*(t(i+1,k,j,n,taum1) - t(i,k,j,n,taum1)) enddo enddo call setbcx (ddxt(1,1,j,n), imt, km) enddo enddo !----------------------------------------------------------------------- ! gradients at northern face of T cells !----------------------------------------------------------------------- do j=max(js-1,1), je-1 jrow = j + joff do n=1,2 do k=1,km kp1 = min(k+1,km) do i=is,ie ddyt(i,k,j,n) = tmask(i,k,j)*tmask(i,k,j+1)*dyur(jrow)* & (t(i,k,j+1,n,taum1) - t(i,k,j,n,taum1)) enddo enddo call setbcx (ddyt(1,1,j,n), imt, km) enddo enddo return end subroutine isopyc (joff, js, je, is, ie) !======================================================================= ! Compute isopycnal diffusion coefficients Ai_ez, Ai_nz, Ai_bx, ! Ai_by, and the Gent/McWilliams advection velocities. ! input: ! joff = offset relating row "j" in the MW to latitude "jrow" ! js = starting row within the MW for calculations ! je = ending row within the MW for calculations ! is = starting index longitude within the MW ! ie = ending index longitude within the MW ! output: ! Ai_ez = diffusion coefficient centered on east face of T cells ! Ai_nz = diffusion coefficient centered on north face of T cells ! Ai_bx = diffusion coefficient centered on bottom face of T cells ! Ai_by = diffusion coefficient centered on bottom face of T cells ! adv_vetiso = isopycnal advective vel on east face of "T" cell ! adv_vntiso = isopycnal advective vel on north face of "T" cell ! (Note: this includes the cosine factor as in "adv_vnt") ! adv_vbtiso = isopycnal advective vel on bottom face of "T" cell !======================================================================= implicit none integer istrt, is, iend, ie, joff, js, je include "size.h" include "param.h" include "pconst.h" include "stdunits.h" !----------------------------------------------------------------------- ! set local constants !----------------------------------------------------------------------- istrt = max(2,is) iend = min(imt-1,ie) !----------------------------------------------------------------------- ! estimate alpha, beta, and gradients on sides of T cells !----------------------------------------------------------------------- call elements (joff, js, je, is, ie) !----------------------------------------------------------------------- ! compute Ai_ez centered on eastern face of T cells ! 1 <= MW < last: calculate for rows jsmw..jemw ! last MW : calculate for rows jsmw..<=jemw !----------------------------------------------------------------------- call ai_east (joff, max(js-1,2), je-1, is, ie) !----------------------------------------------------------------------- ! compute Ai_nz centered on the northern face of T cells ! first MW : calculate for rows 1..jemw ! 1 < MW < last: calculate for rows jsmw..jemw ! last MW : calculate for rows jsmw..<=jemw !----------------------------------------------------------------------- call ai_north (joff, max(js-1,1), je-1, is, ie) !----------------------------------------------------------------------- ! evaluate Ai_bx & Ai_by centered on bottom face of T cells ! 1 <= MW < last: calculate for rows jsmw..jemw ! last MW : calculate for rows jsmw..<=jemw !----------------------------------------------------------------------- call ai_bottom (joff, max(js-1,2), je-1, is, ie) !----------------------------------------------------------------------- ! compute isopycnal advective velocities for tracers ! first MW : calculate (viso) for rows 1..jemw ! (uiso,wiso) for rows 2..jemw ! 1 < MW < last: calculate (uiso,viso,wiso) for rows jsmw..jemw ! last MW : calculate (uiso,viso,wiso) for rows jsmw..<=jemw !----------------------------------------------------------------------- call isopyc_adv (joff, max(js-1,1), je-1, is, ie) return end subroutine ai_east (joff, js, je, is, ie) !======================================================================= ! compute "Ai_ez" & "K11" at the center of eastern face of T cells. !======================================================================= implicit none integer i, k, j, ip, kr, jq, js, je, jrow, joff, ie, is real sc, dzt4r, ai0, sumz, sxe, sumy, facty, ck ! JG - ck include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "accel.h" include "coord.h" include "grdvar.h" include "hmixc.h" include "isopyc.h" include "mw.h" !----------------------------------------------------------------------- ! compute "Ai_ez" on east face of T cell. Note re-scaling factor ! which reduces mixing coefficient "Ai" where slope "sxe" ! exceeds the critical slope "sc" for the small slope approx. For ! the full tensor, it is re-scaled if outside the stable range. !----------------------------------------------------------------------- !JG -- start addisop(1:imtm1,js:je,1:km) = 0; ! Pacific do i=30,83 do j=44,47 do k=4,km addisop(i,k,j) = 1.25e8*j - 5.5e9 c print*,"j=",j,", addisop=",addisop(i,k,j) enddo enddo do j=48,53 do k=4,km addisop(i,k,j) = 5.0e8 enddo enddo do j=54,57 do k=4,km addisop(i,k,j) = -1.25e8*j + 7.125e9 c print*,"j=",j,", addisop=",addisop(i,k,j) enddo enddo enddo ! Indic do i=5,29 do j=44,47 do k=4,km addisop(i,k,j) = 1.25e8*j - 5.5e9 c print*,"j=",j,", addisop=",addisop(i,k,j) enddo enddo do j=48,53 do k=4,km addisop(i,k,j) = 5.0e8 enddo enddo do j=54,57 do k=4,km addisop(i,k,j) = -1.25e8*j + 7.125e9 c print*,"j=",j,", addisop=",addisop(i,k,j) enddo enddo enddo ! Atlantic do i=1,5 do j=44,47 do k=4,km addisop(i,k,j) = 5.e7*j - 2.2e9 c print*,"j=",j,", addisop=",addisop(i,k,j) enddo enddo do j=48,53 do k=4,km addisop(i,k,j) = 2.0e8 enddo enddo do j=54,57 do k=4,km addisop(i,k,j) = -5.e7*j + 2.85e9 c print*,"j=",j,", addisop=",addisop(i,k,j) enddo enddo enddo do i=84,imtm1 do j=44,47 do k=4,km addisop(i,k,j) = 5.e7*j - 2.2e9 c print*,"j=",j,", addisop=",addisop(i,k,j) enddo enddo do j=48,53 do k=4,km addisop(i,k,j) = 2.0e8 enddo enddo do j=54,57 do k=4,km addisop(i,k,j) = -5.e7*j + 2.85e9 c print*,"j=",j,", addisop=",addisop(i,k,j) enddo enddo enddo !JG -- end do j=js,je jrow = j + joff do k=1,km sc = c1/(slmxr*dtxsqr(k)) dzt4r = p5*dzt2r(k) do i=2,imtm1 ! JG - start ck = addisop(i,k,jrow) Ai0 =(ahisop+ck)*p5*(fisop(i,jrow,k)+fisop(i+1,jrow,k)) ! JG - end sumz = c0 do kr=0,1 do ip=0,1 sxe = abs(drodxe(i,k,j,ip)/(drodze(i,k,j,ip,kr)+epsln)) if (sxe .gt. sc) then Ai_ez(i,k,j,ip,kr) = Ai0*tmask(i,k,j)*tmask(i+1,k,j) & *(sc/(sxe + epsln))**2 else Ai_ez(i,k,j,ip,kr) = Ai0*tmask(i,k,j)*tmask(i+1,k,j) endif sumz = sumz + dzw(k-1+kr)*Ai_ez(i,k,j,ip,kr) enddo enddo K11(i,k,j) = dzt4r*sumz enddo enddo call setbcx (Ai_ez(1,1,j,0,0), imt, km) call setbcx (Ai_ez(1,1,j,1,0), imt, km) call setbcx (Ai_ez(1,1,j,0,1), imt, km) call setbcx (Ai_ez(1,1,j,1,1), imt, km) call setbcx (K11(1,1,j), imt, km) enddo return end subroutine ai_north (joff, js, je, is, ie) !======================================================================= ! compute "Ai_nz" & "K22" at the center of northern face of T cells. !======================================================================= implicit none integer i, k, j, ip, kr, jq, js, je, jrow, joff, ie, is real sc, dzt4r, ai0, sumz, syn, sumx, ck include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "accel.h" include "coord.h" include "grdvar.h" include "hmixc.h" include "isopyc.h" include "mw.h" !----------------------------------------------------------------------- ! compute "Ai_nz" on north face of T cell. Note re-scaling factor ! which reduces mixing coefficient "Ai" where slope "syn" ! exceeds the critical slope "sc" for the small slope approx. For ! the full tensor, it is re-scaled if outside the stable range. !----------------------------------------------------------------------- do j=js,je jrow = j + joff do k=1,km sc = c1/(slmxr*dtxsqr(k)) dzt4r = p5*dzt2r(k) do i=2,imtm1 ! print*, "ADDISO", ck Ai0 = ahisop*p5*(fisop(i,jrow,k) + fisop(i,jrow+1,k)) sumz = c0 do kr=0,1 do jq=0,1 syn =abs(drodyn(i,k,j,jq)/(drodzn(i,k,j,jq,kr)+epsln)) if (syn .gt. sc) then Ai_nz(i,k,j,jq,kr) = Ai0*tmask(i,k,j)*tmask(i,k,j+1) & *(sc/(syn + epsln))**2 else Ai_nz(i,k,j,jq,kr) = Ai0*tmask(i,k,j)*tmask(i,k,j+1) endif sumz = sumz + dzw(k-1+kr)*Ai_nz(i,k,j,jq,kr) enddo enddo K22(i,k,j) = dzt4r*sumz enddo enddo call setbcx (Ai_nz(1,1,j,0,0), imt, km) call setbcx (Ai_nz(1,1,j,1,0), imt, km) call setbcx (Ai_nz(1,1,j,0,1), imt, km) call setbcx (Ai_nz(1,1,j,1,1), imt, km) call setbcx (K22(1,1,j), imt, km) enddo return end subroutine ai_bottom (joff, js, je, is, ie) !======================================================================= ! compute Ai_bx, Ai_by, and K33 at the center of the bottom face of ! T cells. !======================================================================= implicit none integer i, k, j, ip, kr, jq, js, je, jrow, joff, ie, is real sc, kp1, ai0, sumx, sxb, sumy, facty, syb, ck include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "accel.h" include "grdvar.h" include "hmixc.h" include "isopyc.h" include "mw.h" !----------------------------------------------------------------------- ! compute "Ai_bx", "Ai_by", & K33 on bottom face of T cell. Note ! re-scaling factor to reduce mixing coefficient "Ai" where slopes ! "sxb" and "syb" exceeds the critical slope "sc" for the small slope approx. For ! the full tensor, it is re-scaled if outside the stable range. !----------------------------------------------------------------------- do j=js,je jrow = j + joff do k=1,kmm1 sc = c1/(slmxr*dtxsqr(k)) kp1 = min(k+1,km) do i=2,imtm1 Ai0 = ahisop*p5*(fisop(i,jrow,k+1) + fisop(i,jrow,k)) ! eastward slopes at the base of T cells sumx = c0 do ip=0,1 do kr=0,1 sxb = abs(drodxb(i,k,j,ip,kr)/(drodzb(i,k,j,kr)+epsln)) if (sxb .gt. sc) then Ai_bx(i,k,j,ip,kr) = Ai0*tmask(i,k+1,j) & *(sc/(sxb + epsln))**2 else Ai_bx(i,k,j,ip,kr) = Ai0*tmask(i,k+1,j) endif sumx = sumx + dxu(i-1+ip)*Ai_bx(i,k,j,ip,kr)*sxb**2 enddo enddo ! northward slopes at the base of T cells sumy = c0 do jq=0,1 facty = csu(jrow-1+jq)*dyu(jrow-1+jq) do kr=0,1 syb = abs(drodyb(i,k,j,jq,kr)/(drodzb(i,k,j,kr)+epsln)) if (syb .gt. sc) then Ai_by(i,k,j,jq,kr) = Ai0*tmask(i,k+1,j) & *(sc/(syb + epsln))**2 else Ai_by(i,k,j,jq,kr) = Ai0*tmask(i,k+1,j) endif sumy = sumy + facty*Ai_by(i,k,j,jq,kr)*syb**2 enddo enddo K33(i,k,j) = dxt4r(i)*sumx + dyt4r(jrow)*cstr(jrow)*sumy enddo enddo call setbcx (Ai_bx(1,1,j,1,0), imt, km) call setbcx (Ai_bx(1,1,j,0,0), imt, km) call setbcx (Ai_bx(1,1,j,1,1), imt, km) call setbcx (Ai_bx(1,1,j,0,1), imt, km) call setbcx (Ai_by(1,1,j,1,0), imt, km) call setbcx (Ai_by(1,1,j,0,0), imt, km) call setbcx (Ai_by(1,1,j,1,1), imt, km) call setbcx (Ai_by(1,1,j,0,1), imt, km) call setbcx (K33(1,1,j), imt, km) enddo return end subroutine isoflux (joff, js, je, is, ie, n) !======================================================================= ! isopycnal diffusive tracer fluxes are computed. !======================================================================= implicit none integer i, k, j, ip, kr, jq, js, je, jrow, joff, km1kr, kpkr integer n, ie, is real dzt4r, sumz, flux_x, ai0, sumy, facty, csu_dzt4r, flux_y real sumx, ck include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "coord.h" include "grdvar.h" include "hmixc.h" include "isopyc.h" include "levind.h" include "mw.h" include "vmixc.h" !----------------------------------------------------------------------- ! construct total isopycnal tracer flux at east face of "T" cells !----------------------------------------------------------------------- do j=js,je jrow = j + joff do k=1,km dzt4r = p5*dzt2r(k) do i=2,imtm1 sumz = c0 do kr=0,1 km1kr = max(k-1+kr,1) kpkr = min(k+kr,km) do ip=0,1 sumz = sumz - Ai_ez(i,k,j,ip,kr) & *(t(i+ip,km1kr,j,n,taum1)-t(i+ip,kpkr,j,n,taum1)) & *drodxe(i,k,j,ip) & /(drodze(i,k,j,ip,kr) + epsln) enddo enddo flux_x = dzt4r*sumz diff_fe(i,k,j) = diff_fe(i,k,j) + K11(i,k,j)*cstdxur(i,j)* & (t(i+1,k,j,n,taum1) - t(i,k,j,n,taum1)) & + flux_x enddo enddo call setbcx (diff_fe(1,1,j), imt, km) enddo !----------------------------------------------------------------------- ! construct total isopycnal tracer flux at north face of "T" cells !----------------------------------------------------------------------- do j=js-1,je jrow = j + joff do k=1,km csu_dzt4r = csu(jrow)*p5*dzt2r(k) do i=2,imtm1 sumz = c0 do kr=0,1 km1kr = max(k-1+kr,1) kpkr = min(k+kr,km) do jq=0,1 sumz = sumz - Ai_nz(i,k,j,jq,kr) & *(t(i,km1kr,j+jq,n,taum1)-t(i,kpkr,j+jq,n,taum1)) & *drodyn(i,k,j,jq) & /(drodzn(i,k,j,jq,kr) + epsln) enddo enddo flux_y = csu_dzt4r*sumz diff_fn(i,k,j) = diff_fn(i,k,j)+K22(i,k,j)*csu_dyur(jrow)* & (t(i,k,j+1,n,taum1)-t(i,k,j,n,taum1)) & + flux_y enddo enddo call setbcx (diff_fn(1,1,j), imt, km) enddo !----------------------------------------------------------------------- ! compute the vertical tracer flux "diff_fbiso" containing the K31 ! and K32 components which are to be solved explicitly. The K33 ! component will be treated implicitly. Note that there are some ! cancellations of dxu(i-1+ip) and dyu(jrow-1+jq) !----------------------------------------------------------------------- do j=js,je jrow = j + joff do k=1,kmm1 do i=2,imtm1 sumx = c0 do ip=0,1 do kr=0,1 sumx = sumx & - Ai_bx(i,k,j,ip,kr)*cstr(jrow)* & (t(i+ip,k+kr,j,n,taum1) - t(i-1+ip,k+kr,j,n,taum1)) & *drodxb(i,k,j,ip,kr) & /(drodzb(i,k,j,kr) + epsln) enddo enddo sumy = c0 do jq=0,1 do kr=0,1 sumy = sumy & - Ai_by(i,k,j,jq,kr)*csu(jrow-1+jq)* & (t(i,k+kr,j+jq,n,taum1)-t(i,k+kr,j-1+jq,n,taum1)) & *drodyb(i,k,j,jq,kr) & /(drodzb(i,k,j,kr) + epsln) enddo enddo diff_fbiso(i,k,j) = dxt4r(i)*sumx & + dyt4r(jrow)*cstr(jrow)*sumy enddo enddo do i=2,imtm1 diff_fbiso(i,0,j) = c0 diff_fbiso(i,km,j) = c0 enddo call setbcx (diff_fbiso(1,0,j), imt, km+1) enddo !----------------------------------------------------------------------- ! compute advective tracer flux at the center of the bottom face of ! the "T" cells !----------------------------------------------------------------------- do j=js,je do k=1,kmm1 do i=2,imt-1 adv_fbiso(i,k,j) = adv_vbtiso(i,k,j)* & (t(i,k,j,n,taum1) + t(i,k+1,j,n,taum1)) enddo enddo enddo ! now consider the top and bottom boundaries do j=js,je do i=2,imt-1 adv_fbiso(i,0,j) = c0 adv_fbiso(i,km,j) = c0 enddo enddo return end subroutine isopyc_adv (joff, js, je, is, ie) !======================================================================= ! compute isopycnal transport velocities. !======================================================================= implicit none integer i, k, j, ip, kr, jq, js, je, jrow, joff, km1, kp1 integer jstrt, ie, is real sc, ath0, at, bt, stn, ab, bb, sbn, absstn, abssbn real ath_t, ath_b, ste, sbe, absste, abssbe include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "accel.h" include "coord.h" include "grdvar.h" include "isopyc.h" include "levind.h" include "mw.h" include "switch.h" include "timeavgs.h" real top_bc(km), bot_bc(km) do k=1,km top_bc(k) = c1 bot_bc(k) = c1 enddo top_bc(1) = c0 bot_bc(km) = c0 !----------------------------------------------------------------------- ! compute the meridional component of the isopycnal mixing velocity ! at the center of the northern face of the "T" cells. !----------------------------------------------------------------------- do j=js,je jrow = j + joff do k=1,km sc = c1/(slmxr*dtxsqr(k)) km1 = max(k-1,1) kp1 = min(k+1,km) do i=1,imt Ath0 = athkdf*p5*(fisop(i,jrow,k) + fisop(i,jrow+1,k)) at = (alphai(i,k,j) + alphai(i,k,j+1) + alphai(i,km1,j) & + alphai(i,km1,j+1)) bt = (betai(i,k,j) + betai(i,k,j+1) + betai(i,km1,j) & + betai(i,km1,j+1)) stn = -(at*(ddyt(i,k,j,1) + ddyt(i,km1,j,1)) & + bt*(ddyt(i,k,j,2) + ddyt(i,km1,j,2))) / & (at*(ddzt(i,km1,j,1) + ddzt(i,km1,j+1,1)) & + bt*(ddzt(i,km1,j,2) + ddzt(i,km1,j+1,2))+epsln) ab = (alphai(i,k,j) + alphai(i,k,j+1) + alphai(i,kp1,j) & + alphai(i,kp1,j+1)) bb = (betai(i,k,j) + betai(i,k,j+1) + betai(i,kp1,j) & + betai(i,kp1,j+1)) sbn = -(ab*(ddyt(i,k,j,1) + ddyt(i,kp1,j,1)) & + bb*(ddyt(i,k,j,2) + ddyt(i,kp1,j,2))) / & (ab*(ddzt(i,k,j,1) + ddzt(i,k,j+1,1)) & + bb*(ddzt(i,k,j,2) + ddzt(i,k,j+1,2))+epsln) absstn = abs(stn) abssbn = abs(sbn) if (absstn .gt. sc) then ath_t = Ath0*tmask(i,k,j)*tmask(i,k,j+1) & *(sc/(absstn + epsln))**2 else ath_t = Ath0*tmask(i,k,j)*tmask(i,k,j+1) endif if (abssbn .gt. sc) then ath_b = Ath0*tmask(i,kp1,j)*tmask(i,kp1,j+1) & *(sc/(abssbn + epsln))**2 else ath_b = Ath0*tmask(i,kp1,j)*tmask(i,kp1,j+1) endif adv_vntiso(i,k,j) = -(ath_t*stn*top_bc(k) - & ath_b*sbn*bot_bc(k))*dztr(k)*csu(jrow) enddo enddo enddo !----------------------------------------------------------------------- ! compute the zonal component of the isopycnal mixing velocity ! at the center of the eastern face of the "T" grid cell. !----------------------------------------------------------------------- jstrt = max(js,jsmw) do j=jstrt,je jrow = j + joff do k=1,km sc = c1/(slmxr*dtxsqr(k)) km1 = max(k-1,1) kp1 = min(k+1,km) do i=1,imtm1 Ath0 = athkdf*p5*(fisop(i,jrow,k) + fisop(i+1,jrow,k)) at = (alphai(i,k,j) + alphai(i+1,k,j) + alphai(i,km1,j) & + alphai(i+1,km1,j)) bt = (betai(i,k,j) + betai(i+1,k,j) + betai(i,km1,j) & + betai(i+1,km1,j)) ste =-(at*(ddxt(i,k,j,1) + ddxt(i,km1,j,1)) & + bt*(ddxt(i,k,j,2) + ddxt(i,km1,j,2))) & / (at*(ddzt(i,km1,j,1) + ddzt(i+1,km1,j,1)) & + bt*(ddzt(i,km1,j,2) + ddzt(i+1,km1,j,2))+epsln) ab = (alphai(i,k,j) + alphai(i+1,k,j) + alphai(i,kp1,j) & + alphai(i+1,kp1,j)) bb = (betai(i,k,j) + betai(i+1,k,j) + betai(i,kp1,j) & + betai(i+1,kp1,j)) sbe =-(ab*(ddxt(i,k,j,1) + ddxt(i,kp1,j,1)) & + bb*(ddxt(i,k,j,2) + ddxt(i,kp1,j,2))) / & (ab*(ddzt(i,k,j,1) + ddzt(i+1,k,j,1)) & + bb*(ddzt(i,k,j,2) + ddzt(i+1,k,j,2))+epsln) absste = abs(ste) abssbe = abs(sbe) if (absste .gt. sc) then ath_t = Ath0*tmask(i,k,j)*tmask(i+1,k,j) & *(sc/(absste + epsln))**2 else ath_t = Ath0*tmask(i,k,j)*tmask(i+1,k,j) endif if (abssbe .gt. sc) then ath_b = Ath0*tmask(i,kp1,j)*tmask(i+1,kp1,j) & *(sc/(abssbe + epsln))**2 else ath_b = Ath0*tmask(i,kp1,j)*tmask(i+1,kp1,j) endif adv_vetiso(i,k,j) = -(ath_t*ste*top_bc(k) - & ath_b*sbe*bot_bc(k))*dztr(k) enddo enddo enddo ! set the boundary conditions do j=jstrt,je call setbcx (adv_vetiso(1,1,j), imt, km) enddo !----------------------------------------------------------------------- ! compute the vertical component of the isopycnal mixing velocity ! at the center of the bottom face of the "T" cells, using the ! continuity equation for the isopycnal mixing velocities !----------------------------------------------------------------------- do j=jstrt,je do i=1,imt adv_vbtiso(i,0,j) = c0 enddo enddo do j=jstrt,je jrow = j + joff do k=1,kmm1 do i=2,imt adv_vbtiso(i,k,j) = dzt(k)*cstr(jrow)*( & (adv_vetiso(i,k,j) - adv_vetiso(i-1,k,j))*dxtr(i) + & (adv_vntiso(i,k,j) - adv_vntiso(i,k,j-1))*dytr(jrow)) enddo enddo enddo do j=jstrt,je do k=1,kmm1 do i=2,imt adv_vbtiso(i,k,j) = adv_vbtiso(i,k,j) + adv_vbtiso(i,k-1,j) enddo enddo enddo do j=jstrt,je jrow = j + joff do i=2,imt adv_vbtiso(i,kmt(i,jrow),j) = c0 enddo enddo ! set the boundary conditions do j=jstrt,je call setbcx (adv_vbtiso(1,0,j), imt, km+1) enddo !----------------------------------------------------------------------- ! accumulate time average gm velocities !----------------------------------------------------------------------- if (timavgperts .and. .not. euler2) then do j=jstrt,je jrow = j + joff do k=1,kmm1 do i=2,imt ta_vetiso(i,k,jrow) = ta_vetiso(i,k,jrow) + & adv_vetiso(i,k,j) ta_vntiso(i,k,jrow) = ta_vntiso(i,k,jrow) + & adv_vntiso(i,k,j) ta_vbtiso(i,k,jrow) = ta_vbtiso(i,k,jrow) + & adv_vbtiso(i,k,j) enddo enddo enddo endif return end