Commit f5fd2b58 authored by Knut's avatar Knut
Browse files

split ss_nn (not ready yet)

parent 400a3e5c
......@@ -14,7 +14,7 @@ TEXSRC = m3d.F90 variables_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 ip_stelling_vankester.F90 \
bdy_3d.F90 rivers.F90 spm.F90 \
start_macro.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 slow_bottom_friction.F90 slow_terms.F90 stop_macro.F90 \
ss_nn.F90 stresses_3d.F90 gotm.F90 tke_eps_advect_3d.F90 numerical_mixing.F90 physical_mixing.F90 structure_friction_3d.F90
shear_frequency.F90 buoyancy_frequency.F90 stresses_3d.F90 gotm.F90 tke_eps_advect_3d.F90 numerical_mixing.F90 physical_mixing.F90 structure_friction_3d.F90
MOD = \
$(LIB)(variables_3d.o) \
......@@ -71,6 +71,7 @@ $(LIB)(ip_song_wright.o) \
$(LIB)(ip_shchepetkin_mcwilliams.o) \
$(LIB)(ip_stelling_vankester.o) \
$(LIB)(ip_chu_fan.o) \
$(LIB)(buoyancy_frequency.o) \
$(LIB)(numerical_mixing.o) \
$(LIB)(physical_mixing.o)
endif
......@@ -80,7 +81,7 @@ $(LIB)(slow_bottom_friction.o) \
$(LIB)(slow_terms.o) \
$(LIB)(adv_split_w.o) \
$(LIB)(stresses_3d.o) \
$(LIB)(ss_nn.o) \
$(LIB)(shear_frequency.o) \
$(LIB)(gotm.o) \
$(LIB)(check_h.o) \
$(LIB)(tke_eps_advect_3d.o) \
......
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: buoyancy_frequency - calculates the buoyancy frequency
! \label{sec-buoyancy-frequency}
!
! !INTERFACE:
subroutine buoyancy_frequency()
!
! !DESCRIPTION:
!
! Here, the buoyancy frequency squared, $N^2=\partial_z b$, with buoyancy $b$ from
! (\ref{bdef}) is calculated.
! Two alternative methods are coded.
! The straight-forward method which is explained first,
! has the disadvantage of generating numerical instabilities.
!
! The straight-forward discretisation of
! $N^2$ is given by
!
! \begin{equation}\label{Nstraight}
! \begin{array}{l}
! \displaystyle
! (N^2)_{i,j,k}\approx
! \frac{b_{i,j,k+1}-b_{i,j,k}}{\frac12(h^t_{i,j,k+1}+h^t_{i,j,k})}.
! \end{array}
! \end{equation}
!
! In some cases, together with the straight-forward discretisation
! of the shear squared, (\ref{ShearSquaredOld}), this
! did not produce stable numerical results. The reason for this might be that
! the velocities involved in the calculation for the shear squared do depend
! on the buoyancies in the two
! neighbouring T-points such that the straight-forward
! method (\ref{Nstraight}) leads to an inconsistency.
! However, other experiments with the energy-conserving discretisation
! of the shear stress squared, (\ref{ShearSquaredNew})
! and the straight-forward discretisation of
! $N^2$, (\ref{Nstraight}), produced numerically stable results.
!
! Most stable results have been obtained with a weighted average
! for the $N^2$ calculation:
!
! \begin{equation}\label{Naveraged}
! \begin{array}{l}
! \displaystyle
! (N^2)_{i,j,k}\approx
! \frac16 \Bigg(
! 2\frac{b_{i,j,k+1}-b_{i,j,k}}{\frac12(h^t_{i,j,k+1}+h^t_{i,j,k})}
! \\ \\ \qquad\qquad\qquad
! \displaystyle
! +
! \frac{b_{i+1,j,k+1}-b_{i+1,j,k}}{\frac12(h^t_{i+1,j,k+1}+h^t_{i+1,j,k})}
! +
! \frac{b_{i-1,j,k+1}-b_{i-1,j,k}}{\frac12(h^t_{i-1,j,k+1}+h^t_{i-1,j,k})}
! \\ \\ \qquad\qquad\qquad
! \displaystyle
! +
! \frac{b_{i,j+1,k+1}-b_{i,j+1,k}}{\frac12(h^t_{i,j+1,k+1}+h^t_{i,j+1,k})}
! +
! \frac{b_{i,j-1,k+1}-b_{i,j-1,k}}{\frac12(h^t_{i,j-1,k+1}+h^t_{i,j-1,k})}
! \Bigg).
! \end{array}
! \end{equation}
!
! These stability issues need to be further investigated in the future.
!
! !USES:
use domain, only: imin,imax,jmin,jmax,kmax,au,av,az
use variables_3d, only: hn
use parameters, only: g,rho_0
use variables_3d, only: NN,buoy,T,S
#ifndef _OLD_BVF_
use variables_3d, only: alpha,beta
#endif
use getm_timers, only: tic, toc, TIM_NN
!$ use omp_lib
IMPLICIT NONE
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
integer :: i,j,k,nb
REALTYPE :: dz,NNc,ttt
REALTYPE :: NNe,NNw,NNn,NNs
REALTYPE, parameter :: small_bvf = 1.d-10
#ifdef _SMOOTH_BVF_VERT_
REALTYPE :: below,center,above
#endif
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'buoyancy_frequency() # ',Ncall
#endif
call tic(TIM_NN)
! BJB-TODO: There are a lot of single- (or lower-) precision constants
! in this routine. Enough that it significantly influences the result
! if they are correctly converted to full (double) precision.
! A conversion should be done, but I will not do it as part of
! the OMP-implementation as it muddles the comparisons of old and
! new results. BJB 2009-09-22
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,nb,dz,NNc,NNe,NNw,NNn,NNs)
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-1,jmax+1
do i=imin-1,imax+1
if (az(i,j) .ge. 1 ) then
NN(i,j,kmax) = small_bvf
do k=kmax-1,1,-1
dz=_HALF_*(hn(i,j,k+1)+hn(i,j,k))
#ifdef _OLD_BVF_
NNc =(buoy(i,j,k+1)-buoy(i,j,k))/dz
#else
NNc = -g / rho_0
NNc = NNc/dz *(alpha(i,j,k) *(T(i,j,k+1)-T(i,j,k)) &
+ beta(i,j,k) *(S(i,j,k+1)-S(i,j,k)))
#endif
if (abs(NNc) .lt. small_bvf ) then
NNc = sign(_ONE_,NNc) * small_bvf
end if
NN(i,j,k)= NNc
end do
#ifdef _SMOOTH_BVF_VERT_
if ( kmax .ge. 4 ) then
below = NN(i,j,1)
center = NN(i,j,2)
above = NN(i,j,3)
do k= 2,kmax-2
center = _HALF_ * center + _QUART_ * (below+above)
below = NN(i,j,k)
NN(i,j,k) = center
center = NN(i,j,k+1)
above = NN(i,j,k+2)
end do
end if
#endif
end if
end do
end do
!$OMP END DO
#ifdef SMOOTH_BVF_HORI
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin,jmax
do i=imin,imax
if (az(i,j) .ge. 1 ) then
do k=kmax-1,1,-1
NNc = NN(i,j,k)
if (az(i+1,j) .ge. 1) then
NNe= NN(i+1,j,k)
else
NNe=NNc
end if
if (az(i-1,j) .ge. 1) then
NNw= NN(i-1,j,k)
else
NNw=NNc
end if
if (az(i,j+1) .ge. 1) then
NNn= NN(i,j+1,k)
else
NNn=NNc
end if
if (az(i,j-1) .ge. 1) then
NNs= NN(i,j-1,k)
else
NNs=NNc
end if
NN(i,j,k) = (_ONE_/3)*NNc+(_ONE_/6)*(NNe+NNw+NNn+NNs)
end do
end if
end do
end do
!$OMP END DO
#endif
#ifdef SLICE_MODEL
do k=1,kmax-1
do i = imin,imax
if (az(i,2) .ge. 1 ) then
NN(i,3,k)=NN(i,2,k)
end if
end do
end do
#endif
!$OMP END PARALLEL
call toc(TIM_NN)
#ifdef DEBUG
write(debug,*) 'Leaving buoyancy_frequency()'
write(debug,*)
#endif
return
end subroutine buoyancy_frequency
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
......@@ -529,6 +529,8 @@
ufirst=.true.
end if
! HERE WE HAVE TO UPDATE SS !!!
#ifndef MUDFLAT
if (kmax .gt. 1) then
if (vert_cord .eq. _ADAPTIVE_COORDS_) call ss_nn()
......@@ -585,6 +587,7 @@
#ifndef PECS
call do_eqstate()
#endif
! HERE WE HAVE TO UPDATE NN !!!
end if
#endif
......
......@@ -2,21 +2,20 @@
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: ss_nn - calculates shear and buoyancy frequency
! \label{sec-ss-nn}
! !ROUTINE: shear_frequency - calculates the shear frequency
! \label{sec-shear-frequency}
!
! !INTERFACE:
subroutine ss_nn()
subroutine shear_frequency()
!
! !DESCRIPTION:
!
! Here, the shear frequency squared,
! $M^2=\left(\partial_z u\right)^2+\left(\partial_z v\right)^2$,
! and the buoyancy frequency squared, $N^2=\partial_z b$, with buoyancy $b$ from
! (\ref{bdef}) are calculated.
! For both calculations, two alternative methods are coded.
! The two straight-forward methods which are explained first, do
! both have the disadvantage of generating numerical instabilities.
! is calculated.
! Two alternative methods are coded.
! The straight-forward method which is explained first,
! has the disadvantage of generating numerical instabilities.
! The straight-forward way for calculating $M^2$ is as follows:
!
! \begin{equation}\label{ShearSquaredOld}
......@@ -76,19 +75,8 @@
! \end{array}
! \end{equation}
!
! The straight-forward discretisation of
! $N^2$ is given by
!
! \begin{equation}\label{Nstraight}
! \begin{array}{l}
! \displaystyle
! (N^2)_{i,j,k}\approx
! \frac{b_{i,j,k+1}-b_{i,j,k}}{\frac12(h^t_{i,j,k+1}+h^t_{i,j,k})}.
! \end{array}
! \end{equation}
!
! In some cases, together with the straight-forward discretisation
! of the shear squared, (\ref{ShearSquaredOld}), this
! of the buoyancy frequency squared, (\ref{Nstraight}), this
! did not produce stable numerical results. The reason for this might be that
! the velocities involved in the calculation for the shear squared do depend
! on the buoyancies in the two
......@@ -99,44 +87,12 @@
! and the straight-forward discretisation of
! $N^2$, (\ref{Nstraight}), produced numerically stable results.
!
! Most stable results have been obtained with a weighted average
! for the $N^2$ calculation:
!
! \begin{equation}\label{Naveraged}
! \begin{array}{l}
! \displaystyle
! (N^2)_{i,j,k}\approx
! \frac16 \Bigg(
! 2\frac{b_{i,j,k+1}-b_{i,j,k}}{\frac12(h^t_{i,j,k+1}+h^t_{i,j,k})}
! \\ \\ \qquad\qquad\qquad
! \displaystyle
! +
! \frac{b_{i+1,j,k+1}-b_{i+1,j,k}}{\frac12(h^t_{i+1,j,k+1}+h^t_{i+1,j,k})}
! +
! \frac{b_{i-1,j,k+1}-b_{i-1,j,k}}{\frac12(h^t_{i-1,j,k+1}+h^t_{i-1,j,k})}
! \\ \\ \qquad\qquad\qquad
! \displaystyle
! +
! \frac{b_{i,j+1,k+1}-b_{i,j+1,k}}{\frac12(h^t_{i,j+1,k+1}+h^t_{i,j+1,k})}
! +
! \frac{b_{i,j-1,k+1}-b_{i,j-1,k}}{\frac12(h^t_{i,j-1,k+1}+h^t_{i,j-1,k})}
! \Bigg).
! \end{array}
! \end{equation}
!
! These stability issues need to be further investigated in the future.
!
! !USES:
use domain, only: imin,imax,jmin,jmax,kmax,au,av,az
use variables_3d, only: kmin,kumin,hn,uu,hun,kvmin,vv,hvn,SS,num
use parameters, only: g,rho_0
#ifndef NO_BAROCLINIC
use variables_3d, only: NN,buoy,T,S
#ifndef _OLD_BVF_
use variables_3d, only: alpha,beta
#endif
#endif
use getm_timers, only: tic, toc, TIM_SSNN
use variables_3d, only: hn,uu,hun,vv,hvn,SS,num
use getm_timers, only: tic, toc, TIM_SS
!$ use omp_lib
IMPLICIT NONE
!
......@@ -144,22 +100,16 @@
! Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
integer :: i,j,k,nb
REALTYPE :: dz,NNc,ttt
REALTYPE :: NNe,NNw,NNn,NNs
REALTYPE, parameter :: small_bvf = 1.d-10
#ifdef _SMOOTH_BVF_VERT_
REALTYPE :: below,center,above
#endif
integer :: i,j,k
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'ss_nn() # ',Ncall
write(debug,*) 'shear_frequency() # ',Ncall
#endif
call tic(TIM_SSNN)
call tic(TIM_SS)
! BJB-TODO: There are a lot of single- (or lower-) precision constants
! in this routine. Enough that it significantly influences the result
......@@ -168,7 +118,7 @@
! the OMP-implementation as it muddles the comparisons of old and
! new results. BJB 2009-09-22
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,nb,dz,NNc,NNe,NNw,NNn,NNs)
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k)
! Prandtl frequency
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin,jmax
......@@ -211,89 +161,11 @@
end do
!$OMP END DO
#ifndef NO_BAROCLINIC
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-1,jmax+1
do i=imin-1,imax+1
if (az(i,j) .ge. 1 ) then
NN(i,j,kmax) = small_bvf
do k=kmax-1,1,-1
dz=_HALF_*(hn(i,j,k+1)+hn(i,j,k))
#ifdef _OLD_BVF_
NNc =(buoy(i,j,k+1)-buoy(i,j,k))/dz
#else
NNc = -g / rho_0
NNc = NNc/dz *(alpha(i,j,k) *(T(i,j,k+1)-T(i,j,k)) &
+ beta(i,j,k) *(S(i,j,k+1)-S(i,j,k)))
#endif
if (abs(NNc) .lt. small_bvf ) then
NNc = sign(_ONE_,NNc) * small_bvf
end if
NN(i,j,k)= NNc
end do
#ifdef _SMOOTH_BVF_VERT_
if ( kmax .ge. 4 ) then
below = NN(i,j,1)
center = NN(i,j,2)
above = NN(i,j,3)
do k= 2,kmax-2
center = _HALF_ * center + _QUART_ * (below+above)
below = NN(i,j,k)
NN(i,j,k) = center
center = NN(i,j,k+1)
above = NN(i,j,k+2)
end do
end if
#endif
end if
end do
end do
!$OMP END DO
#ifdef SMOOTH_BVF_HORI
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin,jmax
do i=imin,imax
if (az(i,j) .ge. 1 ) then
do k=kmax-1,1,-1
NNc = NN(i,j,k)
if (az(i+1,j) .ge. 1) then
NNe= NN(i+1,j,k)
else
NNe=NNc
end if
if (az(i-1,j) .ge. 1) then
NNw= NN(i-1,j,k)
else
NNw=NNc
end if
if (az(i,j+1) .ge. 1) then
NNn= NN(i,j+1,k)
else
NNn=NNc
end if
if (az(i,j-1) .ge. 1) then
NNs= NN(i,j-1,k)
else
NNs=NNc
end if
NN(i,j,k) = (_ONE_/3)*NNc+(_ONE_/6)*(NNe+NNw+NNn+NNs)
end do
end if
end do
end do
!$OMP END DO
#endif
#endif
#ifdef SLICE_MODEL
do k=1,kmax-1
do i = imin,imax
if (az(i,2) .ge. 1 ) then
SS(i,3,k)=SS(i,2,k)
#ifndef NO_BAROCLINIC
NN(i,3,k)=NN(i,2,k)
#endif
end if
end do
end do
......@@ -301,15 +173,14 @@ end do
!$OMP END PARALLEL
call toc(TIM_SSNN)
call toc(TIM_SS)
#ifdef DEBUG
write(debug,*) 'Leaving ss_nn()'
write(debug,*) 'Leaving shear_frequency()'
write(debug,*)
#endif
return
end subroutine ss_nn
end subroutine shear_frequency
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
......@@ -47,7 +47,8 @@
integer, parameter :: TIM_UUMOMENTUMH = 37 ! 3d uu_momentum_3d - halo part only
integer, parameter :: TIM_WWMOMENTUM = 38 ! 3d ww_momentum_3d
integer, parameter :: TIM_WWMOMENTUMH = 39 ! 3d ww_momentum_3d - halo part only
integer, parameter :: TIM_SSNN = 40 ! 3d ss_nn
integer, parameter :: TIM_SS = 40 ! 3d shear_frequency
integer, parameter :: TIM_NN = 41 ! 3d buoyancy_frequency
integer, parameter :: TIM_BOTTFRICT3D = 42 ! 3d bottom_friction_3d
integer, parameter :: TIM_STRESSES3D = 44 ! 3d stresses_3d
integer, parameter :: TIM_STRESSES3DH = 45 ! 3d stresses_3d - halo part only
......@@ -197,7 +198,8 @@
timernames(TIM_VVMOMENTUM) = 'vv_momentum_3d'
timernames(TIM_UUMOMENTUM) = 'uu_momentum_3d'
timernames(TIM_WWMOMENTUM) = 'ww_momentum_3d'
timernames(TIM_SSNN) = 'ss_nn'
timernames(TIM_SS) = 'shear_frequency'
timernames(TIM_NN) = 'buoyancy_frequency'
timernames(TIM_BOTTFRICT3D) = 'bottom_friction_3d'
timernames(TIM_STRESSES3D) = 'stresses_3d'
timernames(TIM_SLOWTERMS) = 'slow_terms'
......
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