Commit ea351e1f authored by kbk's avatar kbk
Browse files

removed old turbulence files

parent 4c2479d2
#$Id: Makefile,v 1.1 2002-05-02 14:01:56 gotm Exp $
#
# Makefile to build the turbulence library - libturb.a.
#
include ../Rules.make
LIB = $(LIBDIR)/libturb${buildtype}.a
MODSRC = turb.F90
LIBSRC = dissipation.F90 pran_kol.F90 production.F90 tke_calc.F90 eddyviscosity.F90
SRC = $(MODSRC) $(LIBSRC)
MOD = \
${LIB}(turb.o)
OBJ = \
${LIB}(dissipation.o) \
${LIB}(pran_kol.o) \
${LIB}(production.o) \
${LIB}(tke_calc.o) \
${LIB}(eddyviscosity.o)
all: modules objects
modules: $(MOD)
objects: $(OBJ)
doc:
$(PROTEX) $(SRC) > $(DOCDIR)/turb.tex
clean:
$(RM) $(LIB) $(MODDIR)/turb.{m,mod}
realclean: clean
$(RM) *.o
distclean: realclean
#-----------------------------------------------------------------------
# Copyright (C) 2001 - Hans Burchard and Karsten Bolding (BBH) !
#-----------------------------------------------------------------------
C
C comturb.h created Mon Jun 8 12:18:54 MET DST 1998
C
C Created using /home/kbk/scr/newdoc by kbk
C
C
C
C $Id: comturb.h,v 1.1 2002-05-02 14:01:56 gotm Exp $
C
C
C
C
C $Log: comturb.h,v $
C Revision 1.1 2002-05-02 14:01:56 gotm
C Initial revision
C
C Revision 1.1.1.1 2001/04/17 08:43:08 bbh
C initial import into CVS
C
C
C
double precision cde,sige ! Parameters for k-epsilon model
common /comturb/
& cde,sige
double precision cmue,tkemin,epsmin,sigk,ce1,ce2,ce3,z0s
parameter(cmue=0.5625d0)
parameter(tkemin=1.0d-7)
parameter(epsmin=5.0d-10)
parameter(sigk=1.0)
parameter(ce1=1.44)
parameter(ce2=1.92)
parameter(ce3=-0.9)
parameter(z0s=0.1)
!$Id: dissipation.F90,v 1.1 2002-05-02 14:01:56 gotm Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: dissipation() - calculates turbulent dissipation.
!
! !INTERFACE:
subroutine dissipation
!
! !DESCRIPTION:
!
! !USES:
use parameters, only: kappa
use domain, only: iimin,iimax,jjmin,jjmax,kmax,az
use m2d, only: zub,zvb
use m3d, only: kmin,hn,ho,num
use m3d, only: dt,taub,eps,P,B,tke,tkeo
use turb, only: ce1,ce2,ce3,z0s,cde,sige,epsmin
IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !REVISION HISTORY:
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: dissipation.F90,v $
! Revision 1.1 2002-05-02 14:01:56 gotm
! Initial revision
!
! Revision 1.1.1.1 2001/04/17 08:43:08 bbh
! initial import into CVS
!
!
! !LOCAL VARIABLES:
integer :: i,j,k
REALTYPE :: prod,buoy,diss,dtl,cee3
REALTYPE :: pplus(1:kmax-1),pminus(1:kmax-1),dif(1:kmax),aux(1:kmax)
REALTYPE :: a1(0:kmax),a2(0:kmax),a3(0:kmax),a4(0:kmax),res(0:kmax)
!
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'dissipation() # ',Ncall
#endif
do i=iimin,iimax
do j=jjmin,jjmax
dtl=dt*az(i,j) ! Applying land mask to time step
do k=kmin(i,j)+1,kmax
dif(k)=0.5*(num(i,j,k-1)+num(i,j,k))/sige
end do
dif(kmin(i,j))=taub(i,j)**2 &
/(0.5*(eps(i,j,kmin(i,j))+eps(i,j,kmin(i,j)-1)))/sige
do k=kmin(i,j),kmax
aux(k)=dtl*dif(k)/hn(i,j,k)
end do
do k=kmin(i,j),kmax-1
prod=ce1*eps(i,j,k)/tkeo(i,j,k)*P(i,j,k)
if (B(i,j,k).gt.0) then
cee3=1.
else
cee3=ce3
end if
buoy=cee3*eps(i,j,k)/tkeo(i,j,k)*B(i,j,k)
diss=ce2*eps(i,j,k)/tkeo(i,j,k)*eps(i,j,k)
if (prod+buoy.gt.0) then
pplus(k)=prod+buoy
pminus(k)=diss
else
pplus(k)=prod
pminus(k)=diss-buoy
end if
end do
do k=kmin(i,j),kmax-1
a3(k)=-2.*aux(k+1)/(hn(i,j,k+1)+hn(i,j,k))
a1(k)=-2.*aux(k )/(hn(i,j,k+1)+hn(i,j,k))
a2(k)=1-a1(k)-a3(k)+dtl*pminus(k)/eps(i,j,k)
a4(k)=(ho(i,j,k+1)+ho(i,j,k))/(hn(i,j,k+1)+hn(i,j,k))*eps(i,j,k) &
+dtl*pplus(k)
end do
a3(kmin(i,j)-1)=0
a2(kmin(i,j)-1)=1.
a4(kmin(i,j)-1)=cde*(tke(i,j,kmin(i,j)-1)**1.5)/kappa &
/(0.25*(zub(i-1,j)+zub(i,j)+zvb(i,j-1)+zvb(i,j)))
a2(kmax)=1.
a1(kmax)=0
a4(kmax)=cde*(tke(i,j,kmax)**1.5)/kappa/z0s
call getm_tridiagonal(kmax,kmin(i,j)-1,kmax,a1,a2,a3,a4,res)
do k=kmin(i,j)-1,kmax
eps(i,j,k)=max(epsmin,res(k))
end do
end do
end do
#ifdef DEBUG
write(debug,*) 'Leaving dissipation()'
write(debug,*)
#endif
return
end subroutine dissipation
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
!$Id: eddyviscosity.F90,v 1.1 2002-05-02 14:01:57 gotm Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: eddyviscosity() - calculates the eddy-viscosity.
!
! !INTERFACE:
subroutine eddyviscosity
!
! !DESCRIPTION:
!
! !USES:
use domain, only: iimin,iimax,jjmin,jjmax,H,min_depth,crit_depth
use meteo, only: tausx,tausy
use m3d, only: kmax,kmin,hn,ssen,num,nuh,taub
use m3d, only: kappa
IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! $Log: eddyviscosity.F90,v $
! Revision 1.1 2002-05-02 14:01:57 gotm
! Initial revision
!
! Revision 1.1.1.1 2001/04/17 08:43:08 bbh
! initial import into CVS
!
!
! !LOCAL VARIABLES:
integer :: i,j,k
REALTYPE :: zz,DD,dry,taus
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'eddyviscosity() # ',Ncall
#endif
do j=jjmin,jjmax
do i=iimin,iimax
zz=0
DD=ssen(i,j)+H(i,j)
dry=max(_ZERO_,min(_ONE_,(DD-min_depth)/(crit_depth-min_depth)))
taus=sqrt(tausx(i,j)**2+tausy(i,j)**2)
do k=kmin(i,j),kmax-1
zz=zz+hn(i,j,k)
num(i,j,k)=kappa*zz*(1-zz/DD)*sqrt(max(taus,taub(i,j))) &
+(1-dry)*1.e-2
end do
end do
end do
#ifdef DEBUG
write(debug,*) 'Leaving eddyviscosity()'
write(debug,*)
#endif
return
end subroutine eddyviscosity
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
!$Id: pran_kol.F90,v 1.1 2002-05-02 14:01:57 gotm Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: pran_kol() - using the cont. eq. to get the sealevel.
!
! !INTERFACE:
subroutine pran_kol
!
! !DESCRIPTION:
!
! !USES:
use domain, only: iimin,iimax,jjmin,jjmax,kmax,min_depth,crit_depth
use m2d, only: D
use m3d, only: kmin,num,tke,eps
use turb, only: cde,cmue
IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! $Log: pran_kol.F90,v $
! Revision 1.1 2002-05-02 14:01:57 gotm
! Initial revision
!
! Revision 1.1.1.1 2001/04/17 08:43:08 bbh
! initial import into CVS
!
!
! !LOCAL VARIABLES:
REALTYPE :: dry
integer :: i,j,k
!
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'pran_kol() # ',Ncall
#endif
!
do i=iimin,iimax
do j=jjmin,jjmax
dry=max(_ZERO_,min(_ONE_,(D(i,j)-min_depth)/(crit_depth-min_depth)))
do k=kmin(i,j),kmax-1
num(i,j,k)=cmue*cde*tke(i,j,k)**2/eps(i,j,k)+(1-dry)*1.e-2
end do
end do
end do
!
#ifdef DEBUG
write(debug,*) 'Leaving pran_kol()'
write(debug,*)
#endif
return
end subroutine pran_kol
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
!$Id: production.F90,v 1.1 2002-05-02 14:01:57 gotm Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: production() -
!
! !INTERFACE:
subroutine production
!
! !DESCRIPTION:
!
! !USES:
use domain, only: iimin,iimax,jjmin,jjmax,kmax
use m3d, only: kmin,kumin,uu,hun,kvmin,vv,hvn,num,P,B
IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! $Log: production.F90,v $
! Revision 1.1 2002-05-02 14:01:57 gotm
! Initial revision
!
! Revision 1.1.1.1 2001/04/17 08:43:08 bbh
! initial import into CVS
!
!
! !LOCAL VARIABLES:
integer :: i,j,k,rc
REALTYPE :: dzu(iimin-1:iimax,jjmin :jjmax,1:kmax-1)
REALTYPE :: dzv(iimin :iimax,jjmin-1:jjmax,1:kmax-1)
!
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'production() # ',Ncall
#endif
do j=jjmin,jjmax
do i=iimin,iimax-1
do k=1,kumin(i,j)-1
dzu(i,j,k)=0
end do
do k=kumin(i,j),kmax-1
dzu(i,j,k)=uu(i,j,k+1)/hun(i,j,k+1)-uu(i,j,k)/hun(i,j,k)
dzu(i,j,k)=dzu(i,j,k)/(0.5*(hun(i,j,k+1)+hun(i,j,k)))
end do
end do
do k=kumin(iimin-1,j),kmax-1
dzu(iimin-1,j,k)=dzu(iimin ,j,k)
end do
do k=kumin(iimax,j),kmax-1
dzu(iimax ,j,k)=dzu(iimax-1,j,k)
end do
end do
do i=iimin,iimax
do j=jjmin,jjmax-1
do k=1,kvmin(i,j)-1
dzv(i,j,k)=0
end do
do k=kvmin(i,j),kmax-1
dzv(i,j,k)=vv(i,j,k+1)/hvn(i,j,k+1)-vv(i,j,k)/hvn(i,j,k)
dzv(i,j,k)=dzv(i,j,k)/(0.5*(hvn(i,j,k+1)+hvn(i,j,k)))
end do
end do
do k=kvmin(i,jjmin-1),kmax-1
dzv(i,jjmin-1,k)=dzv(i,jjmin ,k)
end do
do k=kvmin(i,jjmax),kmax-1
dzv(i,jjmax ,k)=dzv(i,jjmax-1,k)
end do
end do
do i=iimin+1,iimax-1
do j=jjmin+1,jjmax-1
do k=1,kmin(i,j)-1
P(i,j,k)=0
B(i,j,k)=0
end do
do k=kmin(i,j),kmax-1
P(i,j,k)=0.5*( &
0.5*(num(i-1,j ,k)+num(i ,j ,k))*dzu(i-1,j ,k)**2+ &
0.5*(num(i ,j ,k)+num(i+1,j ,k))*dzu(i ,j ,k)**2+ &
0.5*(num(i ,j-1,k)+num(i ,j ,k))*dzv(i ,j-1,k)**2+ &
0.5*(num(i ,j ,k)+num(i ,j+1,k))*dzv(i ,j ,k)**2)
B(i,j,k)=0
end do
end do
end do
i=iimin
do j=jjmin+1,jjmax-1
do k=kmin(i,j)+1,kmax-1
P(i,j,k)=0.5*( &
num(i ,j ,k) *dzu(i-1,j ,k)**2+ &
0.5*(num(i ,j ,k)+num(i+1,j ,k))*dzu(i ,j ,k)**2+ &
0.5*(num(i ,j-1,k)+num(i ,j ,k))*dzv(i ,j-1,k)**2+ &
0.5*(num(i ,j ,k)+num(i ,j+1,k))*dzv(i ,j ,k)**2)
end do
end do
i=iimax
do j=jjmin+1,jjmax-1
do k=kmin(i,j)+1,kmax-1
P(i,j,k)=0.5*( &
0.5*(num(i-1,j ,k)+num(i ,j ,k))*dzu(i-1,j ,k)**2+ &
num(i ,j ,k) *dzu(i ,j ,k)**2+ &
0.5*(num(i ,j-1,k)+num(i ,j ,k))*dzv(i ,j-1,k)**2+ &
0.5*(num(i ,j ,k)+num(i ,j+1,k))*dzv(i ,j ,k)**2)
end do
end do
j=jjmin
do i=iimin+1,iimax-1
do k=kmin(i,j)+1,kmax-1
P(i,j,k)=0.5*( &
0.5*(num(i-1,j ,k)+num(i ,j ,k))*dzu(i-1,j ,k)**2+ &
0.5*(num(i ,j ,k)+num(i+1,j ,k))*dzu(i ,j ,k)**2+ &
num(i ,j ,k) *dzv(i ,j-1,k)**2+ &
0.5*(num(i ,j ,k)+num(i ,j+1,k))*dzv(i ,j ,k)**2)
end do
end do
j=jjmax
do i=iimin+1,iimax-1
do k=kmin(i,j)+1,kmax-1
P(i,j,k)=0.5*( &
0.5*(num(i-1,j ,k)+num(i ,j ,k))*dzu(i-1,j ,k)**2+ &
0.5*(num(i ,j ,k)+num(i+1,j ,k))*dzu(i ,j ,k)**2+ &
0.5*(num(i ,j-1,k)+num(i ,j ,k))*dzv(i ,j-1,k)**2+ &
num(i ,j ,k) *dzv(i ,j ,k)**2)
end do
end do
i=iimin
j=jjmin
do k=kmin(i,j)+1,kmax-1
P(i,j,k)=0.5*( &
num(i ,j ,k) *dzu(i-1,j ,k)**2+ &
0.5*(num(i ,j ,k)+num(i+1,j ,k))*dzu(i ,j ,k)**2+ &
num(i ,j ,k) *dzv(i ,j-1,k)**2+ &
0.5*(num(i ,j ,k)+num(i ,j+1,k))*dzv(i ,j ,k)**2)
end do
i=iimax
j=jjmin
do k=kmin(i,j)+1,kmax-1
P(i,j,k)=0.5*( &
0.5*(num(i-1,j ,k)+num(i ,j ,k))*dzu(i-1,j ,k)**2+ &
num(i ,j ,k) *dzu(i ,j ,k)**2+ &
num(i ,j ,k) *dzv(i ,j-1,k)**2+ &
0.5*(num(i ,j ,k)+num(i ,j+1,k))*dzv(i ,j ,k)**2)
end do
i=iimin
j=jjmax
do k=kmin(i,j)+1,kmax-1
P(i,j,k)=0.5*( &
num(i ,j ,k) *dzu(i-1,j ,k)**2+ &
0.5*(num(i ,j ,k)+num(i+1,j ,k))*dzu(i ,j ,k)**2+ &
0.5*(num(i ,j-1,k)+num(i ,j ,k))*dzv(i ,j-1,k)**2+ &
num(i ,j ,k) *dzv(i ,j ,k)**2)
end do
i=iimax
j=jjmax
do k=kmin(i,j)+1,kmax-1
P(i,j,k)=0.5*( &
0.5*(num(i-1,j ,k)+num(i ,j ,k))*dzu(i-1,j ,k)**2+ &
num(i ,j ,k) *dzu(i ,j ,k)**2+ &
0.5*(num(i ,j-1,k)+num(i ,j ,k))*dzv(i ,j-1,k)**2+ &
num(i ,j ,k) *dzv(i ,j ,k)**2)
end do
#ifdef DEBUG
write(debug,*) 'Leaving production()'
write(debug,*)
#endif
return
end subroutine production
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
!$Id: tke_calc.F90,v 1.1 2002-05-02 14:01:57 gotm Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: tke_calc() -
!
! !INTERFACE:
subroutine tke_calc
!
! !DESCRIPTION:
!
! !USES:
use domain, only: iimin,iimax,jjmin,jjmax,kmax,az
use m3d, only: dt,kmin,hn,ho,num,taub,P,B,eps,tke,tkeo
use turb, only: sigk,cmue,cde,tkemin
IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! $Log: tke_calc.F90,v $
! Revision 1.1 2002-05-02 14:01:57 gotm
! Initial revision
!
! Revision 1.1.1.1 2001/04/17 08:43:08 bbh
! initial import into CVS
!
!
! !LOCAL VARIABLES:
integer :: i,j,k
REALTYPE :: prod,buoy,diss,dtl