Commit fc7d61f2 authored by hb's avatar hb
Browse files

Turbulence advection enabled

parent 7da6990f
#$Id: Makefile,v 1.18 2010-02-23 08:23:34 kb Exp $
#$Id: Makefile,v 1.19 2010-03-02 19:39:32 hb Exp $
#
# Makefile to build the 3D specific library - libm3d.a
#
......@@ -11,9 +11,9 @@ LIB = $(LIBDIR)/lib3d${buildtype}.a
MODSRC = m3d.F90 variables_3d.F90 advection_3d.F90 eqstate.F90 \
temperature.F90 salinity.F90 spm.F90 bdy_3d.F90 rivers.F90
LIBSRC = start_macro.F90 hcc_check.F90 bdy_3d.F90 coordinates.F90 sigma_coordinates.F90 general_coordinates.F90 hybrid_coordinates.F90 check_h.F90 bottom_friction_3d.F90 internal_pressure.F90 uu_momentum_3d.F90 vv_momentum_3d.F90 ww_momentum_3d.F90 structure_friction_3d.F90 uv_advect_3d.F90 uv_diffusion_3d.F90 stresses_3d.F90 ss_nn.F90 gotm.F90 slow_bottom_friction.F90 slow_advection.F90 slow_diffusion.F90 slow_terms.F90 stop_macro.F90
LIBSRC = start_macro.F90 hcc_check.F90 bdy_3d.F90 coordinates.F90 sigma_coordinates.F90 general_coordinates.F90 hybrid_coordinates.F90 check_h.F90 bottom_friction_3d.F90 internal_pressure.F90 uu_momentum_3d.F90 vv_momentum_3d.F90 ww_momentum_3d.F90 structure_friction_3d.F90 uv_advect_3d.F90 uv_diffusion_3d.F90 tke_eps_advect_3d.F90 stresses_3d.F90 ss_nn.F90 gotm.F90 slow_bottom_friction.F90 slow_advection.F90 slow_diffusion.F90 slow_terms.F90 stop_macro.F90
TEXSRC = m3d.F90 variables_3d.F90 coordinates.F90 sigma_coordinates.F90 general_coordinates.F90 hybrid_coordinates.F90 adaptive_coordinates.F90 check_h.F90 hcc_check.F90 uu_momentum_3d.F90 vv_momentum_3d.F90 ww_momentum_3d.F90 structure_friction_3d.F90 uv_advect_3d.F90 uv_diffusion_3d.F90 bottom_friction_3d.F90 internal_pressure.F90 ip_blumberg_mellor.F90 ip_blumberg_mellor_lin.F90 ip_z_interpol.F90 ip_song_wright.F90 ip_chu_fan.F90 ip_shchepetkin_mcwilliams.F90 advection_3d.F90 upstream_adv.F90 upstream_2dh_adv.F90 u_split_adv.F90 v_split_adv.F90 w_split_adv.F90 w_split_it_adv.F90 fct_2dh_adv.F90 temperature.F90 salinity.F90 spm.F90 eqstate.F90 ss_nn.F90 stresses_3d.F90 gotm.F90 rivers.F90 bdy_3d.F90 slow_bottom_friction.F90 slow_advection.F90 slow_diffusion.F90 slow_terms.F90 start_macro.F90 stop_macro.F90
TEXSRC = m3d.F90 variables_3d.F90 coordinates.F90 sigma_coordinates.F90 general_coordinates.F90 hybrid_coordinates.F90 adaptive_coordinates.F90 check_h.F90 hcc_check.F90 uu_momentum_3d.F90 vv_momentum_3d.F90 ww_momentum_3d.F90 structure_friction_3d.F90 uv_advect_3d.F90 uv_diffusion_3d.F90 tke_eps_advect_3d.F90 bottom_friction_3d.F90 internal_pressure.F90 ip_blumberg_mellor.F90 ip_blumberg_mellor_lin.F90 ip_z_interpol.F90 ip_song_wright.F90 ip_chu_fan.F90 ip_shchepetkin_mcwilliams.F90 advection_3d.F90 upstream_adv.F90 upstream_2dh_adv.F90 u_split_adv.F90 v_split_adv.F90 w_split_adv.F90 w_split_it_adv.F90 fct_2dh_adv.F90 temperature.F90 salinity.F90 spm.F90 eqstate.F90 ss_nn.F90 stresses_3d.F90 gotm.F90 rivers.F90 bdy_3d.F90 slow_bottom_friction.F90 slow_advection.F90 slow_diffusion.F90 slow_terms.F90 start_macro.F90 stop_macro.F90
SRC = $(MODSRC) $(LIBSRC)
......@@ -88,7 +88,9 @@ ifeq ($(turbulence),gotm)
OBJ += \
$(LIB)(stresses_3d.o) \
$(LIB)(ss_nn.o) \
$(LIB)(gotm.o)
$(LIB)(gotm.o) \
$(LIB)(tke_eps_advect_3d.o)
endif
......
!$Id: m3d.F90,v 1.47 2010-02-23 08:23:35 kb Exp $
!$Id: m3d.F90,v 1.48 2010-03-02 19:39:32 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -492,6 +492,9 @@
if (vert_cord .ne. 5) call ss_nn()
#endif
call gotm()
#ifdef TURB_ADV
call tke_eps_advect_3d(vel_hor_adv,vel_ver_adv,vel_adv_split)
#endif
#endif
end if
#ifndef NO_BAROCLINIC
......
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: tke_eps_advect_3d - 3D turbulence advection
!
! !INTERFACE:
subroutine tke_eps_advect_3d(hor_adv,ver_adv,adv_split)
!
! !DESCRIPTION:
!
! This routine carries out advection of the prognostic turbulence quantities
! {\tt tke} (turbuent kinetic energy, $k$) and {\tt eps} (lenght scale related
! turbulence quantity, e.g.\ dissipation rate of $k$, $\varepsilon$, or
! turbulent frequency, $\omega=\varepsilon/k$. Here, the TVD advection
! schemes are used which are also used for the momentum advection.
!
! !
! !USES:
use domain, only: imin,imax,jmin,jmax,kmax,az,au,av
#if defined(SPHERICAL) || defined(CURVILINEAR)
use domain, only: dxu,dxv,dyu,dyv,arcd1
#else
use domain, only: dx,dy,ard1
#endif
use variables_3d, only: dt,kumin,kvmin,uu,vv,ww,hun,hvn,ho,hn,uuEx,vvEx
#ifdef UV_TVD
use variables_3d, only: uadv,vadv,wadv,huadv,hvadv,hoadv,hnadv
#endif
use variables_3d, only: tke,eps
use advection_3d, only: do_advection_3d
use halo_zones, only: update_3d_halo,wait_halo,U_TAG,V_TAG
IMPLICIT NONE
!
! !INPUT PARAMETERS:
integer, intent(in) :: hor_adv,ver_adv,adv_split
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
integer :: i,j,k
#ifdef UV_TVD
REALTYPE :: dxuadv(I2DFIELD)
REALTYPE :: dxvadv(I2DFIELD)
REALTYPE :: dyuadv(I2DFIELD)
REALTYPE :: dyvadv(I2DFIELD)
REALTYPE :: area_inv(I2DFIELD)
REALTYPE :: AH=_ZERO_
REALTYPE :: dxdyi
#endif
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'tke_eps_advect_3d() # ',Ncall
#endif
#ifndef UV_TVD
STDERR 'Do not use tke_eps_advect_3d() without the compiler option UV_TVD'
stop 'tke_eps_advect_3d()'
#else
#if defined(SPHERICAL) || defined(CURVILINEAR)
#else
dxdyi = _ONE_/(dx*dy)
#endif
! Here begins dimensional split advection for tke and eps
do k=1,kmax-1
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
uadv(i,j,k)=_HALF_*(uu(i,j,k)+uu(i,j,k+1))
vadv(i,j,k)=_HALF_*(vv(i,j,k)+vv(i,j,k+1))
wadv(i,j,k)=_HALF_*(ww(i,j,k)+ww(i,j,k+1))
huadv(i,j,k)=_HALF_*(hun(i,j,k)+hun(i,j,k+1))
hvadv(i,j,k)=_HALF_*(hvn(i,j,k)+hvn(i,j,k+1))
hoadv(i,j,k)=_HALF_*(ho(i,j,k)+ho(i,j,k+1))
hnadv(i,j,k)=_HALF_*(hn(i,j,k)+hn(i,j,k+1))
end do
end do
end do
k=kmax ! Only needed to allow for advection of k=kmax
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
uadv(i,j,k)=_HALF_*(uu(i,j,k))
vadv(i,j,k)=_HALF_*(vv(i,j,k))
wadv(i,j,k)=_HALF_*(ww(i,j,k))
huadv(i,j,k)=_HALF_*(hun(i,j,k))
hvadv(i,j,k)=_HALF_*(hvn(i,j,k))
hoadv(i,j,k)=_HALF_*(ho(i,j,k))
hnadv(i,j,k)=_HALF_*(hn(i,j,k))
end do
end do
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
#if defined(SPHERICAL) || defined(CURVILINEAR)
dxuadv(i,j)=dxu(i,j)
dxvadv(i,j)=dxv(i,j)
dyuadv(i,j)=dyu(i,j)
dyvadv(i,j)=dyv(i,j)
area_inv(i,j)=arcd1(i,j)
#else
dxuadv(i,j)=dx
dxvadv(i,j)=dx
dyuadv(i,j)=dy
dyvadv(i,j)=dy
area_inv(i,j)=dxdyi
#endif
end do
end do
call do_advection_3d(dt,tke,uadv,vadv,wadv,huadv,hvadv,hoadv,hnadv, &
dxuadv,dxvadv,dyuadv,dyvadv,area_inv, &
az,au,av,hor_adv,ver_adv,adv_split,AH)
call do_advection_3d(dt,eps,uadv,vadv,wadv,huadv,hvadv,hoadv,hnadv, &
dxuadv,dxvadv,dyuadv,dyvadv,area_inv, &
az,au,av,hor_adv,ver_adv,adv_split,AH)
#endif
#ifdef DEBUG
write(debug,*) 'Leaving tke_eps_advect_3d()'
write(debug,*)
#endif
return
end subroutine tke_eps_advect_3d
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2010 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment