#include"cppdefs.h" !------------------------------------------------------------------------- !BOP ! ! !ROUTINE: Algebraic length-scale with two master scales \label{sec:potentialml} ! ! !INTERFACE: subroutine potentialml(nlev,z0b,z0s,h,depth,NN) ! !DESCRIPTION: ! Computes the length scale by defining two master ! length scales $l_u$ and $l_d$ ! \begin{equation} ! \begin{array}{l} ! \int_{z_0}^{z_0+l_u(z_0)} (b(z_0)-b(z)) dz =k(z_0) \comma \\[4mm] ! \int_{z_0-l_d(z_0)}^{z_0} (b(z)-b(z_0)) dz =k(z_0) ! \end{array} ! \end{equation} ! ! From $l_u$ and $l_d$ two length--scales are defined: $l_k$, ! a characteristic mixing length, ! and $l_\epsilon$, a characteristic dissipation length. ! They are computed according to ! \begin{equation} ! \begin{array}{l} ! l_k(z_0)= \text{Min} ( l_d(z_0),l_u(z_0)) \comma \\[4mm] ! l_{\epsilon}(z_0)=\left( l_d(z_0)l_u(z_0)\right)^\frac{1}{2} ! \point ! \end{array} ! \end{equation} ! ! $l_k$ is used in {\tt kolpran()} to compute eddy viscosity/difussivity. ! $l_{\epsilon}$ is used to compute the dissipation rate, $\epsilon$ ! according to ! \begin{equation} ! \epsilon=C_{\epsilon} k^{3/2} l_{\epsilon}^{-1} ! \comma ! C_{\epsilon}=0.7 ! \point ! \end{equation} ! ! !USES: use turbulence, only: L,eps,tke,k_min,eps_min use turbulence, only: cde,galp,kappa,length_lim IMPLICIT NONE ! ! !INPUT PARAMETERS: ! number of vertical layers integer, intent(in) :: nlev ! bottom and surface roughness (m) REALTYPE, intent(in) :: z0b,z0s ! layer thickness (m) REALTYPE, intent(in) :: h(0:nlev) ! local depth (m) REALTYPE, intent(in) :: depth ! buoyancy frequency (1/s^2) REALTYPE, intent(in) :: NN(0:nlev) ! ! !REVISION HISTORY: ! Original author(s): Manuel Ruiz Villarreal, Hans Burchard ! ! $Log: potentialml.F90,v $ ! Revision 20.0 2013/12/14 00:13:48 fms ! Merged revision 1.1.2.1 onto trunk ! ! Revision 1.1.2.1 2012/05/15 16:00:53 smg ! initial cvs ci for these modules to mom5. ! AUTHOR:Griffies ! REVIEWERS: ! TEST STATUS: ! CHANGES PUBLIC INTERFACES? ! CHANGES ANSWERS? ! ! Revision 1.1.2.1.390.1 2012/04/23 20:30:29 smg ! updated to the gotm-2012.03.09 CVS tag. ! AUTHOR:Martin Schmidt ! REVIEWERS:Griffies ! TEST STATUS: ! CHANGES PUBLIC INTERFACES? ! CHANGES ANSWERS? ! ! Revision 1.6 2005-11-15 11:35:02 lars ! documentation finish for print ! ! Revision 1.5 2005/06/27 13:44:07 kbk ! modified + removed traling blanks ! ! Revision 1.4 2003/03/28 09:20:35 kbk ! added new copyright to files ! ! Revision 1.3 2003/03/10 09:02:05 gotm ! Added new Generic Turbulence Model + ! improved documentation and cleaned up code ! ! Revision 1.2 2002/02/08 08:59:59 gotm ! Revision 1.1.1.1 2001/02/12 15:55:58 gotm ! initial import into CVS ! !EOP !------------------------------------------------------------------------- ! ! !LOCAL VARIABLES: integer :: i,j REALTYPE :: ds(0:nlev),db(0:nlev) REALTYPE :: lu(0:nlev),ld(0:nlev) REALTYPE :: lk(0:nlev),leps(0:nlev) REALTYPE :: Lcrit,buoydiff,integral,ceps REALTYPE, parameter :: NNmin=1.e-8 ! !------------------------------------------------------------------------- !BOC db(0)=0. ds(nlev)=0. do i=1,nlev-1 db(i)=db(i-1)+h(i) ! distance of intercace i from bottom ds(i)=depth-db(i) ! distance of intercace i from surface end do ! ! Calculation of lu and ld by solving the integral equation following ! Gaspar (1990). Some other approximations of the integral equation ! are possible. ! ! Computation of lupward ! do i=1,nlev-1 lu(i)=0. integral=0. buoydiff=0. do j=i+1,nlev buoydiff=buoydiff+NN(j-1)*0.5*(h(j)+h(j-1)) integral=integral+buoydiff*h(j) if (integral.ge.tke(i)) then if(j.ne.nlev) then if(j.ne.i+1) then lu(i)=lu(i)-(integral-tke(i))/buoydiff else ! To avoid lu(i) from becoming too large if NN(i) is too small if(NN(i).gt.NNmin) then lu(i)=sqrt(2.)*sqrt(tke(i))/sqrt(NN(i)) else lu(i)=h(i) end if end if goto 600 end if end if lu(i)=lu(i)+h(j) end do 600 continue ! Implicitely done in the do loop: if (lu(i).gt.ds(i)) lu(i)=ds(i) ! lu limited by distance to surface end do ! Computation of ldownward do i=nlev-1,1,-1 ld(i)=0. integral=0. buoydiff=0. do j=i-1,1,-1 buoydiff=buoydiff+NN(j)*0.5*(h(j+1)+h(j)) integral=integral-buoydiff*h(j) if (integral.ge.tke(i)) then if(j.ne.0) then if(j.ne.i-1) then ld(i)=ld(i)-(integral-tke(i))/buoydiff else ! To avoid ld(i) from becoming too large if NN(i) is too small if(NN(i).gt.NNmin) then ld(i)=sqrt(2.)*sqrt(tke(i))/sqrt(NN(i)) else ld(i)=h(i) end if end if goto 610 end if end if ld(i)=ld(i)+h(j) end do 610 continue ! if (ld(i).gt.db(i)) ld(i)=db(i) !ld limited by distance to bottom end do ! Calculation of lk and leps, mixing and dissipation lengths do i=nlev-1,1,-1 ! Suggested by Gaspar: lk(i) = min(lu(i),ld(i)) lk(i)=sqrt(lu(i)*ld(i)) leps(i) = sqrt(lu(i)*ld(i)) end do ! We set L=lk because it is the one we use to calculate num and nuh ceps=0.7 do i=1,nlev-1 L(i)=lk(i) end do ! do the boundaries assuming linear log-law length-scale L(0)=kappa*z0b L(nlev)=kappa*z0s do i=0,nlev ! clip the length-scale at the Galperin et al. (1988) value ! under stable stratifcitation if ((NN(i).gt.0).and.(length_lim)) then Lcrit=sqrt(2*galp*galp*tke(i)/NN(i)) if (L(i).gt.Lcrit) L(i)=Lcrit end if ! compute the dissipation rate eps(i)=cde*sqrt(tke(i)*tke(i)*tke(i))/L(i) ! substitute minimum value if (eps(i).lt.eps_min) then eps(i) = eps_min L(i) = cde*sqrt(tke(i)*tke(i)*tke(i))/eps_min endif enddo return end subroutine potentialml !EOC !----------------------------------------------------------------------- ! Copyright by the GOTM-team under the GNU Public License - www.gnu.org !-----------------------------------------------------------------------