subroutine mxkrtm(m,n) use mod_xc ! HYCOM communication interface implicit none c include 'common_blocks.h' c integer m,n c c --- hycom version 1.0 (adapted from micom version 2.8) c real, save, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & sdot c integer i,j,k,l real delp,q,thk ccc integer kmax ccc real totem,tosal,tndcyt,tndcys,work(3) c !$OMP PARALLEL DO PRIVATE(j,k,l,i) !$OMP& SCHEDULE(STATIC,jblk) do 1 j=1-margin,jj+margin do 1 k=1,kk do 1 l=1,isp(j) do 1 i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) p(i,j,k+1)=p(i,j,k)+dp(i,j,k,n) 1 continue !$OMP END PARALLEL DO c 103 format (i9,2i5,a/(33x,i3,2f8.3,f8.3,0p,f8.2,f8.1)) cdiag if (itest.gt.0 .and. jtest.gt.0) write (lp,103) nstep,itest,jtest, cdiag. ' entering mxlayr: temp saln dens thkns dpth', cdiag. (k,temp(itest,jtest,k,n),saln(itest,jtest,k,n), cdiag. th3d(itest,jtest,k,n)+thbase,dp(itest,jtest,k,n)*qonem, cdiag. p(itest,jtest,k+1)*qonem,k=1,kk) c if (thermo .or. sstflg.gt.0 .or. srelax) then c c --- ----------------------------------- c --- mixed layer entrainment/detrainment c --- ----------------------------------- c !$OMP PARALLEL DO PRIVATE(j) !$OMP& SHARED(m,n,sdot) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin call mxkrtmaj(m,n, sdot, j) enddo !$OMP END PARALLEL DO c else !.not.thermo ... c !$OMP PARALLEL DO PRIVATE(j,l,i) !$OMP& SCHEDULE(STATIC,jblk) do 87 j=1-margin,jj+margin do 87 l=1,isp(j) do 87 i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) surflx(i,j)=0. salflx(i,j)=0. sdot(i,j)=dp(i,j,1,n) 87 continue !$OMP END PARALLEL DO c end if !thermo.or.sstflg.gt.0.or.srelax:else c cdiag if (itest.gt.0.and.jtest.gt.0.and.turgen(itest,jtest).lt.0.) cdiag. write (lp,'(i9,2i5,a,f8.2)') nstep,itest,jtest, cdiag. ' monin-obukhov length (m):',sdot(itest,jtest)*qonem c c --- store 'old' t/s column integral in totem/tosal (diagnostic use only) ccc totem=0. ccc tosal=0. ccc do k=1,kk ccc if (max(dp(itest,jtest,1,n)+sdot(itest,jtest),thkmin*onem).gt. ccc . p(itest,jtest,k) .or. max(th3d(itest,jtest,1,m),th3d(itest, ccc . jtest,1,n)) +sigjmp.ge.th3d(i,j,k,n)) then ccc kmax=k ccc totem=totem+temp(itest,jtest,k,n)*dp(itest,jtest,k,n) ccc tosal=tosal+saln(itest,jtest,k,n)*dp(itest,jtest,k,n) ccc end if ccc end do c !$OMP PARALLEL DO PRIVATE(j) !$OMP& SHARED(m,n,sdot) !$OMP& SCHEDULE(STATIC,jblk) do j=1-margin,jj+margin call mxkrtmbj(m,n, sdot, j) enddo !$OMP END PARALLEL DO c c --- compare 'old' with 'new' t/s column integral (diagnostic use only) c ccc tndcyt=-totem ccc tndcys=-tosal ccc do k=kmax,1,-1 ccc tndcyt=tndcyt+temp(itest,jtest,k,n)*dp(itest,jtest,k,n) ccc tndcys=tndcys+saln(itest,jtest,k,n)*dp(itest,jtest,k,n) ccc end do ccc write (lp,'(i9,2i5,i3,3x,a,1p,3e10.2/25x,a,3e10.2)') nstep,itest, ccc . jtest,kmax,'total saln,srf.flux,tndcy:',tosal/g,salflx(itest, ccc . jtest)*delt1,tndcys/g,'total temp,srf.flux,tndcy:',totem/g, ccc . surflx(itest,jtest)*delt1,tndcyt*spcifh/g c c --- store 'old' interface pressures in -pu,pv- c !$OMP PARALLEL DO PRIVATE(j,k,l,i) !$OMP& SCHEDULE(STATIC,jblk) do 882 j=1-margin,jj+margin do 882 k=2,kk+1 c do 881 l=1,isu(j) do 881 i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) 881 pu(i,j,k)=min(depthu(i,j),.5*(p(i,j,k)+p(i-1,j,k))) c do 882 l=1,isv(j) do 882 i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) pv(i,j,k)=min(depthv(i,j),.5*(p(i,j,k)+p(i,j-1,k))) 882 continue !$OMP END PARALLEL DO c c --- store 'new' layer thicknesses in -dpu,dpv- c !$OMP PARALLEL DO PRIVATE(j,k,l,i) !$OMP& SCHEDULE(STATIC,jblk) do j=1,jj ! no margin because p's halo is updated in dpudpv do k=1,kk do l=1,isp(j) do i=max(1,ifp(j,l)),min(ii,ilp(j,l)) p(i,j,k+1)=p(i,j,k)+dp(i,j,k,n) enddo enddo enddo enddo !$OMP END PARALLEL DO c call dpudpv(dpu(1-nbdy,1-nbdy,1,n), & dpv(1-nbdy,1-nbdy,1,n), & p,depthu,depthv, margin) ! p's halo updated by dpudpv c c --- redistribute momentum in the vertical. c --- homogenize (u,v) over depth range defined in -util1,util2- c c --- thk>0 activates momentum diffusion across mixed-layer interface thk=vertmx*onem*delt1 c !$OMP PARALLEL DO PRIVATE(j,l,i,k,delp,q) !$OMP& SCHEDULE(STATIC,jblk) do 97 j=1-margin,jj+margin c do 83 l=1,isu(j) c do 822 i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) util1(i,j)=max(dpu(i,j,1,n),pu(i,j,2)+thk) uflux(i,j)=0. 822 util3(i,j)=0. c do 82 k=1,kk do 82 i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) delp=max(0.,min(util1(i,j),pu(i,j,k+1)) . -min(util1(i,j),pu(i,j,k ))) uflux(i,j)=uflux(i,j)+u(i,j,k,n)*delp 82 util3(i,j)=util3(i,j) +delp c do 83 i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) u(i,j,1,n)=uflux(i,j)/util3(i,j) 83 continue c do 84 l=1,isv(j) c do 844 i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) util2(i,j)=max(dpv(i,j,1,n),pv(i,j,2)+thk) vflux(i,j)=0. 844 util4(i,j)=0. c do 80 k=1,kk do 80 i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) delp=max(0.,min(util2(i,j),pv(i,j,k+1)) . -min(util2(i,j),pv(i,j,k ))) vflux(i,j)=vflux(i,j)+v(i,j,k,n)*delp 80 util4(i,j)=util4(i,j) +delp c do 84 i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) v(i,j,1,n)=vflux(i,j)/util4(i,j) 84 continue c do 97 k=2,kk c do 96 l=1,isu(j) do 96 i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l)) pu(i,j,k)=pu(i,j,k-1)+dpu(i,j,k-1,n) q=max(0.,min(1.,(util1(i,j)-pu(i,j,k))/(dpu(i,j,k,n)+epsil))) 96 u(i,j,k,n)=u(i,j,1,n)*q+u(i,j,k,n)*(1.-q) c do 97 l=1,isv(j) do 97 i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l)) pv(i,j,k)=pv(i,j,k-1)+dpv(i,j,k-1,n) q=max(0.,min(1.,(util2(i,j)-pv(i,j,k))/(dpv(i,j,k,n)+epsil))) v(i,j,k,n)=v(i,j,1,n)*q+v(i,j,k,n)*(1.-q) 97 continue !$OMP END PARALLEL DO c cdiag if (itest.gt.0 .and. jtest.gt.0) write (lp,103) nstep,itest,jtest, cdiag. ' exiting mxlayr: temp saln dens thkns dpth', cdiag. (k,temp(itest,jtest,k,n),saln(itest,jtest,k,n), cdiag. th3d(itest,jtest,k,n)+thbase,dp(itest,jtest,k,n)*qonem, cdiag. p(itest,jtest,k+1)*qonem,k=1,kk) return end subroutine mxkrtmaj(m,n, sdot, j) use mod_xc ! HYCOM communication interface implicit none c include 'common_blocks.h' c integer m,n, j real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & sdot c c --- hycom version 1.0 (adapted from micom version 2.8) c integer i,k,ka,l c real thknss,ustar3,dpth,ekminv,obuinv,buoyfl,dsgdt,tmn,smn, . ex,alf1,alf2,cp1,cp3,ape,cc4,spe,pnew,alfadt,betads,thet c real ea1, ea2, em1, em2, em3, em4, em5 data ea1, ea2, em1, em2, em3, em4, em5 . /0.60,0.30,0.45,2.60,1.90,2.30,0.60/ ! Gaspar coefficients c include 'stmt_fns.h' c locsig=.true. c c --- ----------------------------------- c --- mixed layer entrainment/detrainment c --- ----------------------------------- c do 85 l=1,isp(j) c do 86 i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) c c --- determine turb.kin.energy generation due to wind stirring c --- ustar computed in subr. -thermf- c --- buoyancy flux (m**2/sec**3), all fluxes into the ocean c --- note: surface density increases (column is destabilized) if buoyfl < 0 thknss=dp(i,j,1,n) ustar3=ustar(i,j)**3 tmn=.5*(temp(i,j,1,m)+temp(i,j,1,n)) smn=.5*(saln(i,j,1,m)+saln(i,j,1,n)) dsgdt=dsigdt(tmn,smn) buoyfl=-g*thref*(dsigds(tmn,smn)*salflx(i,j)*thref+ & dsgdt *surflx(i,j)*thref/spcifh) c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c --- option 1 : k r a u s - t u r n e r mixed-layer t.k.e. closure c ccc em=0.8*exp(-p(i,j,2)/(50.*onem)) ! hadley centre choice (orig.: 1.25) ccc en=0.15 ! hadley centre choice (orig.: 0.4) ccc thermg=-.5*g*((en+1.)*buoyfl+(en-1.)*abs(buoyfl))*qthref ccc turgen(i,j)=delt1*(2.*em*g*ustar3*qthref+thknss*thermg)*qthref**2 c c --- find monin-obukhov length in case of receding mixed layer (turgen < 0). c --- the monin-obukhov length is found by stipulating turgen = 0. c --- store temporarily in 'sdot'. c ccc if (turgen(i,j).lt.0.) then ccc sdot(i,j)=-2.*em*g*ustar3/min(-epsil,thref*thermg) ccc else ccc sdot(i,j)=thknss ccc end if c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c --- option 2 : g a s p a r mixed-layer t.k.e. closure c dpth=thknss*qonem ekminv=1./hekman(i,j) obuinv=buoyfl/max(epsil,ustar3) ex=exp(min(50.,dpth*obuinv)) alf1=ea1+ea2*max(1.,2.5*dpth*ekminv)*ex alf2=ea1+ea2*ex cp1=((1.-em5)*(alf1/alf2)+.5*em4)*athird cp3=max(0.,(em4*(em2+em3)-(alf1/alf2)*(em2+em3-em3*em5))*athird) ape=cp3*ustar3-cp1*dpth*buoyfl c if(ape.lt.0.) then ! detrainment turgen(i,j)=(g*delt1*qthref**3)*ape sdot(i,j)=max(thkmin*onem,min(thknss,g*cp3/ .(thref*cp1*max(epsil,obuinv)))) c else ! entrainment cc4=2.*em4/(em1*em1) * alf1*alf1 spe=(em2+em3)*ustar3-0.5*dpth*buoyfl turgen(i,j)=(g*delt1*qthref**3)*(sqrt((.5*ape-cp1*spe)**2 . +2.*cc4*ape*spe)-(.5*ape+cp1*spe))/(cc4-cp1) sdot(i,j)=thknss end if c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c --- util1,util2 are used to evaluate pot.energy changes during entrainment util1(i,j)=th3d(i,j,1,n)*thknss 86 util2(i,j)=th3d(i,j,1,n)*thknss**2 c c --- find pnew in case of mixed layer deepening (turgen > 0). store in 'sdot'. c --- entrain as many layers as needed to deplete -turgen-. c do 85 k=2,kk ka=k-1 do 85 i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) if (k.eq.2) then thstar(i,j,ka,1)=th3d(i,j,ka,n) endif if (locsig) then alfadt=0.5* & (dsiglocdt(temp(i,j,ka,n),saln(i,j,ka,n),p(i,j,k))+ & dsiglocdt(temp(i,j,k ,n),saln(i,j,k ,n),p(i,j,k)))* & (temp(i,j,ka,n)-temp(i,j,k,n)) betads=0.5* & (dsiglocds(temp(i,j,ka,n),saln(i,j,ka,n),p(i,j,k))+ & dsiglocds(temp(i,j,k ,n),saln(i,j,k ,n),p(i,j,k)))* & (saln(i,j,ka,n)-saln(i,j,k,n)) thstar(i,j,k,1)=thstar(i,j,ka,1)-alfadt-betads thet=thstar(i,j,k,1) else thet=th3d(i,j,k,n) endif pnew=(2.*turgen(i,j)+thet*p(i,j,k)**2-util2(i,j))/ . max(epsil,thet*p(i,j,k) -util1(i,j)) c --- stop iterating for 'pnew' as soon as pnew < k-th interface pressure if (pnew.lt.p(i,j,k)) pnew=sdot(i,j) c --- substitute 'pnew' for monin-obukhov length if mixed layer is deepening if (turgen(i,j).ge.0.) sdot(i,j)=pnew c util1(i,j)=util1(i,j)+thet*dp(i,j,k,n) util2(i,j)=util2(i,j)+thet*(p(i,j,k+1)**2-p(i,j,k)**2) 85 continue return end subroutine mxkrtmbj(m,n, sdot, j) use mod_xc ! HYCOM communication interface implicit none c include 'common_blocks.h' c integer m,n, j real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: & sdot c c --- hycom version 1.0 (adapted from micom version 2.8) c integer i,k,ktr,l,num c real tdp(idm),sdp(idm) real pnew,thknss,t1,s1,tmxl,smxl, . dpn,sn,tn,dtemp,dsaln,tnew,snew,z,s_up,a,e,b,f,d,c1msig, . cc0,cc3,cc1,cc2,x c real ccubq,ccubr,ccubqr,ccubs1,ccubs2,ccubrl,ccubim,root,root1, . root2,root3 c include 'stmt_fns.h' c c --- cubic eqn. solver used in mixed-layer detrainment ccubq(s)=athird*(cc1/cc3-athird*(cc2/cc3)**2) ccubr(s)=athird*(.5*(cc1*cc2)/(cc3*cc3)-1.5*cc0/cc3) . -(athird*cc2/cc3)**3 ccubqr(s)=sqrt(abs(ccubq(s)**3+ccubr(s)**2)) ccubs1(s)=sign(abs(ccubr(s)+ccubqr(s))**athird,ccubr(s)+ccubqr(s)) ccubs2(s)=sign(abs(ccubr(s)-ccubqr(s))**athird,ccubr(s)-ccubqr(s)) root(s)=ccubs1(s)+ccubs2(s)-athird*cc2/cc3 ccubrl(s)=sqrt(max(0.,-ccubq(s))) . *cos(athird*atan2(ccubqr(s),ccubr(s))) ccubim(s)=sqrt(max(0.,-ccubq(s))) . *sin(athird*atan2(ccubqr(s),ccubr(s))) root1(s)=2.*ccubrl(s)-athird*cc2/cc3 root2(s)=-ccubrl(s)+sqrt(3.)*ccubim(s)-athird*cc2/cc3 root3(s)=-ccubrl(s)-sqrt(3.)*ccubim(s)-athird*cc2/cc3 c do 26 l=1,isp(j) c do 42 i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) c --- store (pnew - pold) in 'sdot'. c --- don't allow mixed layer to get too deep or too shallow. sdot(i,j)=min(p(i,j,kk+1),max(thkmin*onem,sdot(i,j)))- . dp(i,j,1,n) klist(i,j)=2 tdp(i)=0. 42 sdp(i)=0. c do 43 k=2,kk do 43 i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) pnew=dp(i,j,1,n)+sdot(i,j) c --- 'tdp,sdp' will be needed for temp./salin. mixing during entrainment tdp(i)=tdp(i)+temp(i,j,k,n)*(min(pnew,p(i,j,k+1)) . -min(pnew,p(i,j,k ))) sdp(i)=sdp(i)+saln(i,j,k,n)*(min(pnew,p(i,j,k+1)) . -min(pnew,p(i,j,k ))) c c --- if sdot > 0, remove water from layers about to be entrained. dpold(i,j,k)=dp(i,j,k,n) ! diapyc.flux dp(i,j,k,n)=max(p(i,j,k+1),pnew)-max(p(i,j,k),pnew) diaflx(i,j,k)=diaflx(i,j,k)+(dp(i,j,k,n)-dpold(i,j,k)) ! diapyc.flux if (pnew.ge.p(i,j,k+1)) then do ktr= 1,ntracr tracer(i,j,k,n,ktr)=0. enddo endif c c --- if sdot < 0, mixed layer water will be detrained into isopycnic layer c --- defined in -klist-. to prevent odd/even time step decoupling of mixed- c --- layer depth, determine -klist- from layer one -th3d- at 2 consecutive c --- time levels c if (max(th3d(i,j,1,m),th3d(i,j,1,n))+sigjmp.ge.th3d(i,j,k,n)) . klist(i,j)=k+1 c c --- set t/s in massless layers. step 1: copy salinity from layer(s) above c saln(i,j,k,n)=(saln(i,j,k,n)*dp(i,j,k,n)+saln(i,j,k-1,n)*epsil)/ . ( dp(i,j,k,n)+ epsil) 43 continue c c --- set t/s in massless layers. step 2: copy salinity from layer(s) below c do 44 k=kk-1,2,-1 do 44 i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) 44 saln(i,j,k,n)=(saln(i,j,k,n)*dp(i,j,k,n)+saln(i,j,k+1,n)*epsil)/ . ( dp(i,j,k,n)+ epsil) c c --- set t/s in massless layers. step 3: increase salinity where water c --- is too fresh to fit into layer k c do 45 k=2,kk do 45 i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) if (saln(i,j,k,n).lt.salmin(k)) then saln(i,j,k,n)=salmin(k) temp(i,j,k,n)=tofsig(th3d(i,j,k,n)+thbase,saln(i,j,k,n)) end if 45 continue c c --- redistribute temp. and salin. during both de- and entrainment c do 26 i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l)) thknss=dp(i,j,1,n) pnew=thknss+sdot(i,j) t1=temp(i,j,1,n) s1=saln(i,j,1,n) c tmxl=t1+surflx(i,j)*delt1*g/(spcifh*thknss) smxl=s1+salflx(i,j)*delt1*g/ thknss c cdiag if (i.eq.itest.and.j.eq.jtest) write (lp,'(i9,2i5,a,3f7.3,f8.2)') cdiag. nstep,i,j,' t,s,sig,dp after diab.forcing',tmxl,smxl, cdiag. sig(tmxl,smxl),thknss*qonem c if (sdot(i,j).ge.0.) then c c --- (mixed layer d e e p e n s) c cdiag if (i.eq.itest.and.j.eq.jtest) write (lp,'(i9,2i5,a,f9.3,a)') cdiag. nstep,i,j,' entrain',sdot(i,j)*qonem,' m of water' c tmxl=(tmxl*thknss+tdp(i))/pnew smxl=(smxl*thknss+sdp(i))/pnew dp(i,j,1,n)=pnew diaflx(i,j,1)=diaflx(i,j,1)+sdot(i,j) ! diapyc.flux c else if (sdot(i,j).lt.-onecm.and.surflx(i,j).ge.0.) then ! sdot < 0 c c --- (mixed layer r e c e d e s) c k=klist(i,j) if (k.gt.kk) go to 27 c cdiag if (i.eq.itest.and.j.eq.jtest) cdiag. write (lp,'(i9,2i5,a,i2,a,3p,2f7.3)') nstep,i,j, cdiag. ' sig\*(1),sig\*(',k,') =',th3d(i,j,1,n)+thbase, cdiag. th3d(i,j,k,n)+thbase c dpn=max(dp(i,j,k,n),0.) sn=saln(i,j,k,n) tn=temp(i,j,k,n) c cdiag if (i.eq.itest.and.j.eq.jtest) cdiag. write (lp,'(i9,2i5,i3,a,2f9.4,f8.2)') nstep,i,j,k, cdiag. ' t,s,dp before detrainment',tn,sn,dpn*qonem c c --- distribute last time step's heating and freshwater flux over depth range c --- 'pnew' (monin-obukhov length). split fossil mixed layer (depth= -sdot= c --- thknss-pnew) into lower part ('lo') of depth z cooled and detrained into c --- layer k, and an upper part ('up') heated to match temperature rise in c --- mixed layer. transfer as much salinity as possible from sublayer 'up' to c --- sublayer 'lo' without creating new maxima/minima in water column. c dtemp=delt1*g*surflx(i,j)/(spcifh*pnew) dsaln=delt1*g*salflx(i,j)/ pnew c tnew= t1+dtemp snew=max(s1+dsaln,0.0) !must be non-negative c if (s1.le.sn .and. t1.gt.tn) then c c --- scenario 1: transfer t/s so as to achieve t_lo = t_k, s_lo = s_k c z=-sdot(i,j)*min(1.,dtemp/max(epsil,tnew-tn))*qonem s_up=s1+(s1-sn)*dtemp/max(epsil*dtemp,t1-tn) c --- is scenario 1 feasible? if (s_up.ge.min(snew,s1)) go to 24 end if ! s_1 < s_n c c --- scenario 2: (t_lo,s_lo) differ from (tn,sn). main problem now is in c --- maintaining density in layer k during detrainment. This requires solving c --- 3rd deg. polynomial cc3*z**3 + cc2*z**2 + cc1*z + cc0 = 0 for z. c s_up=min(s1,snew) c --- new (t,s) in layer k will be t=(a*z+b)/(z+d), s=(e*z+f)/(z+d). a=tnew e=s_up b=(tn*dpn+ dtemp*sdot(i,j))*qonem f=(sn*dpn+(s_up-s1)*sdot(i,j))*qonem d=dpn*qonem c c1msig=c1-(th3d(i,j,k,n)+thbase) cc0=d*d*(d*c1msig+b*c2+f*c3)+b*(d*f*c5+b*(d*c4+b*c6+f*c7)) cc3= ( c1msig+a*c2+e*c3)+a*( e*c5+a*( c4+a*c6+e*c7)) cc1=d*(3. *d*c1msig+(2.*b +a*d)*c2+(2. *f+d*e)*c3)+b*((2.*a*d . +b )*c4+3.*a*b*c6+(2.*a*f+b*e)*c7)+(a*d*f+b*(d*e+ f))*c5 cc2= (3. *d*c1msig+(2.*a*d+b )*c2+(2.*d*e+ f)*c3)+a*((2.*b . +a*d)*c4+3.*a*b*c6+(2.*b*e+a*f)*c7)+(b *e+a*( f+d*e))*c5 c --- bound cc3 away from zero cc3=sign(max(1.e-6,abs(cc3)),cc3) c x=0.0 ! dummy argument that is never used if (ccubq(x)**3+ccubr(x)**2.gt.0.) then c --- one real root num=1 z=root(x) else c --- three real roots num=3 z=root1(x) end if c cdiag if (i.eq.itest.and.j.eq.jtest) then cdiag work(1)=z cdiag if (num.eq.3) then cdiag work(2)=root2(x) cdiag work(3)=root3(x) cdiag end if cdiag write (lp,100) nstep,i,j,' t,s,dp( 1)=',tnew,snew, cdiag. thknss*qonem,'sdot,z=',sdot(i,j)*qonem,z,'t,s,dp(',k,')=',tn, cdiag. sn,dpn*qonem,'real root(s):',(work(nu),nu=1,num) cdiag end if 100 format (i9,2i5,a,2f7.3,f8.2,3x,a,2f8.2/20x,a,i2,a,2f7.3,f8.2, . 3x,a,1p3e11.4) c c --- does root fall into appropriate range? if (z.le.0.005) go to 27 c c --- ready to detrain lowest 'z' meters from mixed layer c temp(i,j,k,n)=(a*z+b)/(z+d) saln(i,j,k,n)=(e*z+f)/(z+d) c 24 sdot(i,j)=max(sdot(i,j),-z*onem) dp(i,j,1,n)=thknss+sdot(i,j) dp(i,j,k,n) =dpn -sdot(i,j) smxl=(snew*pnew+s_up*(dp(i,j,1,n)-pnew))/dp(i,j,1,n) tmxl=tnew diaflx(i,j,1)=diaflx(i,j,1)+sdot(i,j) ! diapyc.flux diaflx(i,j,k)=diaflx(i,j,k)-sdot(i,j) ! diapyc.flux c c --- inject 'ventilation' tracer into layer k do ktr= 1,ntracr tracer(i,j,k,n,ktr)=(tracer(i,j,k,n,ktr)*dpn-sdot(i,j)) & /(dpn-sdot(i,j)) enddo c cdiag if (i.eq.itest.and.j.eq.jtest) cdiag. write (lp,'(i9,2i5,i3,a,2f9.4,f8.2)') nstep,i,j,k, cdiag. ' t,s,dp after detrainment',temp(i,j,k,n),saln(i,j,k,n), cdiag. dp(i,j,k,n)*qonem c end if ! sdot > or < 0 c 27 temp(i,j,1,n)=tmxl saln(i,j,1,n)=smxl th3d(i,j,1,n)=sig(tmxl,smxl)-thbase do ktr= 1,ntracr tracer(i,j,1,n,ktr)=1.0 enddo c dpmixl(i,j,n)=dp( i,j,1,n) dpbl( i,j) =dp( i,j,1,n) tmix( i,j) =temp(i,j,1,n) smix( i,j) =saln(i,j,1,n) thmix( i,j) =th3d(i,j,1,n) c cdiag if (i.eq.itest.and.j.eq.jtest) write cdiag. (lp,'(i9,2i5,i3,a,2f9.4,f8.2)') nstep,i,j,1, cdiag. ' final mixed-layer t,s,dp ',tmxl,smxl,dp(i,j,1,n)*qonem c 26 continue return end c c c> Revision history: c> c> June 1995 - removed restriction 'klist(i,j) .le. kk' c> June 1995 - added code for setting t/s in massless layers below mix.layer c> Oct. 1995 - removed bug created while changing klist (June 1995 revision): c> 'if (k.gt.kk) go to 26' now reads 'if (k.gt.kk) go to 27' c> May 1997 - changed -sdot- into local array c> Mar. 1998 - added -th3d- c> Nov. 1998 - fixed bug in computing tnew,snew in situations where z < 0.005 c> Dec. 1998 - replaced dsaln by (s_up-s1) in definition of 'f' c> Feb. 1999 - limited 'tofsig' call in loop 45 to cases where saln < salmin c> Aug. 2000 - adapted from micom 2.8 to run within hycom 1.0 c> May 2002 - buoyfl (into the ocean), calculated here