Commit 84abb257 authored by hb's avatar hb
Browse files

Inclusion of friction due to structures in the water column

parent 0c66dbf1
#$Id: Makefile,v 1.15 2007-11-09 10:43:43 kb Exp $
#$Id: Makefile,v 1.16 2008-03-26 13:25:51 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 adaptive_coordinates.F90 check_h.F90 bottom_friction_3d.F90 internal_pressure.F90 uu_momentum_3d.F90 vv_momentum_3d.F90 ww_momentum_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 adaptive_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
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 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 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)
......@@ -50,7 +50,12 @@ $(LIB)(adaptive_coordinates.o) \
$(LIB)(bottom_friction_3d.o) \
$(LIB)(uu_momentum_3d.o) \
$(LIB)(vv_momentum_3d.o) \
$(LIB)(ww_momentum_3d.o) \
$(LIB)(ww_momentum_3d.o)
ifeq ($(GETM_STRUCTURE_FRICTION),true)
OBJ += \
$(LIB)(structure_friction_3d.o)
endif
OBJ += \
$(LIB)(uv_advect_3d.o) \
$(LIB)(uv_diffusion_3d.o)
ifneq ($(GETM_NO_BAROCLINIC),true)
......
......@@ -7,6 +7,11 @@
allocate(ww(I3DFIELD),stat=rc) ! 3D field for w-velocity
if (rc /= 0) stop 'init_3d: Error allocating memory (ww)'
#ifdef STRUCTURE_FRICTION
allocate(sf(I3DFIELD),stat=rc) ! 3D field for velocity in T-points
if (rc /= 0) stop 'init_3d: Error allocating memory (sf)'
#endif
allocate(ho(I3DFIELD),stat=rc) ! 3D field for old box height (z-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (ho)'
......
......@@ -10,6 +10,9 @@
integer, dimension(:,:), allocatable:: kmin_pmz,kumin_pmz,kvmin_pmz
REALTYPE, dimension(:,:,:), allocatable :: uu,vv,ww
#ifdef STRUCTURE_FRICTION
REALTYPE, dimension(:,:,:), allocatable :: sf
#endif
REALTYPE, dimension(:,:,:), allocatable :: ho,hn
REALTYPE, dimension(:,:,:), allocatable :: huo,hun
REALTYPE, dimension(:,:,:), allocatable :: hvo,hvn
......
!$Id: gotm.F90,v 1.15 2007-06-07 10:25:19 kbk Exp $
!$Id: gotm.F90,v 1.16 2008-03-26 13:25:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -31,6 +31,15 @@
! vertical vectors {\tt tke1d}, {\tt eps1d}, {\tt num1d} and {\tt nuh1d}
! are transformed back to 3D fields.
!
! In case that the compiler option {\tt STRUCTURE\_FRICTION} is switched on,
! the additional turbulence production by structures in the water column is calculated
! by calculating the total production as
! \begin{equation}
! P_{tot} = P +C \left(u^2+v^2\right)^{3/2},
! \end{equation}
! with the shear production $P$, and the structure friction coefficient $C$. The
! latter is calculated in the routine {\tt structure\_friction\_3d.F90}.
!
! There are furthermore a number of compiler options provided, e.g.\
! for an older GOTM version, for barotropic calcuations,
! and for simple parabolic viscosity profiles circumventing the GOTM
......@@ -45,6 +54,9 @@
use variables_3d, only: NN,nuh
#endif
use variables_3d, only: avmback,avhback
#ifdef STRUCTURE_FRICTION
use variables_3d, only: uu,vv,hun,hvn,sf
#endif
use turbulence, only: do_turbulence,cde
use turbulence, only: tke1d => tke, eps1d => eps, L1d => L
use turbulence, only: num1d => num, nuh1d => nuh
......@@ -65,6 +77,10 @@
REALTYPE :: h(0:kmax),dry,zz
REALTYPE :: NN1d(0:kmax),SS1d(0:kmax)
REALTYPE :: xP(0:kmax)
#ifdef STRUCTURE_FRICTION
REALTYPE :: cds=0.01
#endif
!
!EOP
!-----------------------------------------------------------------------
......@@ -81,6 +97,15 @@
if (az(i,j) .eq. 1 ) then
#ifdef STRUCTURE_FRICTION
do k=1,kmax
xP(k)= 0.25*( &
(uu(i ,j ,k)/hun(i ,j ,k))**2*(sf(i ,j ,k)+sf(i+1,j ,k)) &
+(uu(i-1,j ,k)/hun(i-1,j ,k))**2*(sf(i-1,j ,k)+sf(i ,j ,k)) &
+(vv(i ,j ,k)/hvn(i ,j ,k))**2*(sf(i ,j ,k)+sf(i ,j+1,k)) &
+(vv(i ,j-1,k)/hvn(i ,j-1,k))**2*(sf(i ,j-1,k)+sf(i ,j ,k)))
end do
#endif
u_taus = sqrt(taus(i,j))
u_taub = sqrt(taub(i,j))
......
!$Id: m3d.F90,v 1.40 2008-02-07 08:27:44 kb Exp $
!$Id: m3d.F90,v 1.41 2008-03-26 13:25:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -400,6 +400,9 @@
hvo=hvn
ip_fac=_ONE_
if (ip_ramp .gt. 0) ip_fac=min( _ONE_ , n*_ONE_/ip_ramp)
#ifdef STRUCTURE_FRICTION
call structure_friction_3d
#endif
if (ufirst) then
call uu_momentum_3d(n,bdy3d)
call vv_momentum_3d(n,bdy3d)
......
!$Id: slow_terms.F90,v 1.8 2007-06-07 10:25:19 kbk Exp $
!$Id: slow_terms.F90,v 1.9 2008-03-26 13:25:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -28,6 +28,10 @@
#ifndef NO_BAROCLINIC
use variables_3d, only: idpdx,idpdy
#endif
#ifdef STRUCTURE_FRICTION
use variables_3d, only: sf
#endif
IMPLICIT NONE
!
! !INPUT PARAMETERS:
......@@ -147,6 +151,11 @@
Slru(i,j)=-Uint(i,j)/(0.5*(ssuo(i,j)+ssun(i,j)) &
+HU(i,j))*ru(i,j) &
+uu(i,j,k)/(0.5*(huo(i,j,k)+hun(i,j,k)))*rru(i,j)
#endif
#ifdef STRUCTURE_FRICTION
do k=1,kmax
Slru(i,j)=Slru(i,j)+uu(i,j,k)*0.5*(sf(i,j,k)+sf(i-1,j,k))
end do
#endif
else
Slru(i,j)= _ZERO_
......@@ -167,6 +176,11 @@
Slrv(i,j)=-Vint(i,j)/(0.5*(ssvo(i,j)+ssvn(i,j)) &
+HV(i,j))*rv(i,j) &
+vv(i,j,k)/(0.5*(hvo(i,j,k)+hvn(i,j,k)))*rrv(i,j)
#endif
#ifdef STRUCTURE_FRICTION
do k=1,kmax
Slrv(i,j)=Slrv(i,j)+vv(i,j,k)*0.5*(sf(i,j,k)+sf(i,j-1,k))
end do
#endif
else
Slrv(i,j)=_ZERO_
......
......@@ -23,6 +23,9 @@
REALTYPE :: uu(I3DFIELD)
REALTYPE :: vv(I3DFIELD)
REALTYPE :: ww(I3DFIELD)
#ifdef STRUCTURE_FRICTION
REALTYPE :: sf(I3DFIELD)
#endif
REALTYPE :: ho(I3DFIELD)
REALTYPE :: hn(I3DFIELD)
REALTYPE :: huo(I3DFIELD)
......
!$Id: structure_friction_3d.F90,v 1.1 2008-03-26 13:25:51 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: structure_friction_3d -
! \label{sec-structure_friction-3d}
!
! !INTERFACE:
subroutine structure_friction_3d()
!
! !DESCRIPTION:
! Here, the quadratic friction term resulting from a structure in the water column is
! calculated. This term will be added as additional forcing to the
! three-dimensional momentum equations, where it is treated numerically
! implicitly. Therefore here, only the following terms is calculated:
! \begin{equation}
! \mbox{\tt sf} = C(z) \sqrt{u(z)^2+v(z)^2},
! \end{equation}
! with the friction coefficient $C$ bearing the physical unit [1/m].
!
! !USES:
use domain, only: imin,imax,jmin,jmax,kmax
use variables_3d, only: uu,vv,sf,huo,hvo
! #define CALC_HALO_WW
#ifndef CALC_HALO_WW
use domain, only: az
use halo_zones, only: update_3d_halo,wait_halo,z_TAG
#endif
IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
REALTYPE :: dtm1
integer :: i,j,k
#ifdef STRUCTURE_FRICTION
REALTYPE :: cds(I2DFIELD)
#endif
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'structure_friction_3d() # ',Ncall
#endif
! cds is the friction coefficient which a structure imposes to a
! flow field. It has the unit [1/m] and should be a 3d field. For the
! first, a dependency on the vertical coordinate is not assumed.
! cds should be later read from a netcdf input file.
cds=0
do j=jmin,jmax
do i=imin,imax
if ((i.eq.90).and.(j.eq.6)) then
cds(i,j)=0.01
end if
end do
end do
sf=_ZERO_
do k=1,kmax
#ifdef CALC_HALO_WW
do j=jmin-1,jmax+1
do i=imin-1,imax+1
#else
do j=jmin,jmax
do i=imin,imax
#endif
sf(i,j,k)=cds(i,j)*sqrt(0.5*( &
(uu(i ,j ,k)/(huo(i ,j ,k)))**2 &
+(uu(i-1,j ,k)/(huo(i-1,j ,k)))**2 &
+(vv(i ,j ,k)/(hvo(i ,j ,k)))**2 &
+(vv(i ,j-1,k)/(hvo(i ,j-1,k)))**2 &
))
end do
end do
end do
#ifndef CALC_HALO_WW
call update_3d_halo(sf,sf,az,imin,jmin,imax,jmax,kmax,z_TAG)
call wait_halo(z_TAG)
#endif
#ifdef DEBUG
write(debug,*) 'Leaving structure_friction_3d()'
write(debug,*)
#endif
return
end subroutine structure_friction_3d
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
!$Id: uu_momentum_3d.F90,v 1.15 2007-06-07 10:25:19 kbk Exp $
!$Id: uu_momentum_3d.F90,v 1.16 2008-03-26 13:25:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -29,6 +29,13 @@
! Afterwards, the matrix is set up for each water column, and it is solved
! by means of a tri-diagonal matrix solver.
!
! In case that the compiler option {\tt STRUCTURE\_FRICTION} is switched on,
! the frictional effect of structures in the water column is calculated
! by adding the quadratic frictional term $C u \sqrt{u^2+v^2}$ (with a minus sign on
! the right hand side) numerically implicitly to the $u$-equation,
! with the friction coefficient $C$. The explicit part of this term, $C\sqrt{u^2+v^2}$,
! is calculated in the routine {\tt structure\_friction\_3d.F90}.
!
! Finally, the new velocity profile is shifted such that its vertical
! integral is identical to the time integral of the vertically integrated
! transport.
......@@ -55,6 +62,9 @@
use variables_3d, only: dt,cnpar,kumin,uu,vv,huo,hun,hvo,uuEx,ww,hvn
use variables_3d, only: num,nuh,sseo,ssun,rru
use variables_3d, only: ssuo
#ifdef STRUCTURE_FRICTION
use variables_3d, only: sf
#endif
#ifndef NO_BAROCLINIC
use variables_3d, only: idpdx
#endif
......@@ -176,6 +186,12 @@
+dt*ex(k) &
-dt*0.5*(huo(i,j,k)+hun(i,j,k))*g*zx
#ifdef STRUCTURE_FRICTION
do k=kumin(i,j),kmax
a2(k)=a2(k)+dt*0.5*(sf(i,j,k)+sf(i+1,j,k))
end do
#endif
call getm_tridiagonal(kmax,kumin(i,j),kmax,a1,a2,a3,a4,Res)
! Transport correction: the integral of the new velocities has to
......
!$Id: vv_momentum_3d.F90,v 1.18 2007-06-07 10:25:19 kbk Exp $
!$Id: vv_momentum_3d.F90,v 1.19 2008-03-26 13:25:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -28,6 +28,15 @@
! section \ref{Section_dry}).
! Afterwards, the matrix is set up for each water column, and it is solved
! by means of a tri-diagonal matrix solver.
!
! In case that the compiler option {\tt STRUCTURE\_FRICTION} is switched on,
! the frictional effect of structures in the water column is calculated
! by adding the quadratic frictional term $C v \sqrt{u^2+v^2}$ (with a minus
! sign on
! the right hand side) numerically implicitly to the $v$-equation,
! with the friction coefficient $C$. The explicit part of this term,
! $C\sqrt{u^2+v^2}$,
! is calculated in the routine {\tt structure\_friction\_3d.F90}.
!
! Finally, the new velocity profile is shifted such that its vertical
! integral is identical to the time integral of the vertically integrated
......@@ -61,6 +70,9 @@
#ifdef XZ_PLUME_TEST
use variables_3d, only: buoy
#endif
#ifdef STRUCTURE_FRICTION
use variables_3d, only: sf
#endif
#ifndef NO_BAROCLINIC
use variables_3d, only: idpdy
#endif
......@@ -191,6 +203,12 @@
+dt*ex(k) &
-dt*0.5*(hvo(i,j,k)+hvn(i,j,k))*g*zy
#ifdef STRUCTURE_FRICTION
do k=kvmin(i,j),kmax
a2(k)=a2(k)+dt*0.5*(sf(i,j,k)+sf(i,j+1,k))
end do
#endif
call getm_tridiagonal(kmax,kvmin(i,j),kmax,a1,a2,a3,a4,Res)
! Transport correction: the integral of the new velocities has to
......
#$Id: Rules.make,v 1.15 2007-10-24 11:23:33 kbk Exp $
#$Id: Rules.make,v 1.16 2008-03-26 13:25:52 hb Exp $
#
# This file contains rules which are shared between multiple Makefiles.
# This file is quite complicated - all compilation options are set in this
......@@ -41,6 +41,11 @@ ifeq ($(GETM_SPM),true)
DEFINES += -DSPM
endif
# Structure friction
ifeq ($(GETM_STRUCTURE_FRICTION),true)
DEFINES += -DSTRUCTURE_FRICTION
endif
# Bio-geochemical component
ifeq ($(GETM_BIO),true)
DEFINES += -DGETM_BIO
......
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