updategrid.F90 9.38 KB
!$Id: updategrid.F90,v 1.19 2007-01-06 11:49:16 kbk Exp$
#include"cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: The vertical grid \label{sec:updategrid}
!
! !INTERFACE:
subroutine updategrid(nlev,dt,zeta)
!
! !DESCRIPTION:
!  This subroutine calculates for each time step new layer thicknesses
!  in order to fit them to the changing water depth.
!  Three different grids can be specified:
!  \begin{enumerate}
!  \item Equidistant grid with possible zooming towards surface and bottom.
!  The number of layers, {\tt nlev}, and the zooming factors,
!  {\tt ddu}=$d_u$ and  {\tt ddl}=$d_l$,
!  are specified in {\tt gotmmean.nml}.
!  Zooming is applied according to the formula
!  \begin{equation}\label{formula_Antoine}
!    h_k = D\frac{\mbox{tanh}\left( (d_l+d_u)\frac{k}{M}-d_l\right)
!    +\mbox{tanh}(d_l)}{\mbox{tanh}(d_l)+\mbox{tanh}(d_u)}-1
!   \point
!  \end{equation}
!
!  From this formula, the following grids are constructed:
!  \begin{itemize}
!    \item $d_l=d_u=0$ results in equidistant discretisations.
!    \item $d_l>0, d_u=0$ results in zooming near the bottom.
!    \item $d_l=0, d_u>0$ results in zooming near the surface.
!    \item $d_l>0, d_u>0$ results in double zooming nea both,
!          the surface and the bottom.
!  \end{itemize}
!
!  \item Sigma-layers. The fraction that every layer occupies is
!  read-in from file, see {\tt gotmmean.nml}.
!  \item Cartesian layers. The height of every layer is read in from file,
!  see {\tt gotmmean.nml}.
!  This method is not recommended when a varying sea surface is considered.
!  \end{enumerate}
!
!  Furthermore, vertical velocity profiles are calculated here, if
!  {\tt w\_adv\_method} is 1 or 2, which has to be chosen in the
!  {\tt w\_advspec} namelist in {\tt obs.nml}. The profiles of vertical
!  velocity are determined by two values,
!  the height of maximum absolute value of vertical velocity, {\tt w\_height},
!  and the vertical velocity at this height, {\tt w\_adv}. From {\tt w\_height},
!  the vertical velocity is linearly decreasing towards the surface and
!  the bottom, where is value is zero.

!
! !USES:
use meanflow,     only: depth0,depth
use meanflow,     only: ga,z,h,ho,ddu,ddl,grid_method
use meanflow,     only: NN,SS,w_grid,grid_file,w
IMPLICIT NONE
!
! !INPUT PARAMETERS:
integer, intent(in)                 :: nlev
REALTYPE, intent(in)                :: dt,zeta
!
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!  $Log: updategrid.F90,v$
!  Revision 1.19  2007-01-06 11:49:16  kbk
!  namelist file extension changed .inp --> .nml
!
!  Revision 1.18  2006-11-27 15:26:37  kbk
!  initialise grid depending on grid_ready
!
!  Revision 1.17  2006-11-24 15:13:41  kbk
!  de-allocate memory and close open files
!
!  Revision 1.16  2005-11-18 10:59:35  kbk
!  removed unused variables - some left in parameter lists
!
!  Revision 1.15  2005/11/15 11:39:32  lars
!  documentation finish for print
!
!  Revision 1.14  2005/08/25 19:41:33  hb
!  small deviations between depth and depth0 tolerated now
!
!  Revision 1.13  2005/08/15 20:23:40  hb
!  Vertical advection profiles triangle-shaped also for
!  temporally constant vertical velocity
!
!  Revision 1.12  2005/06/27 13:44:07  kbk
!  modified + removed traling blanks
!
!  Revision 1.11  2004/08/18 11:46:19  lars
!  updated documentation
!
!  Revision 1.10  2003/07/23 10:52:52  hb
!  proper initialisation of gridinit + cleaning
!
!  Revision 1.9  2003/03/28 09:20:35  kbk
!
!  Revision 1.8  2003/03/28 08:56:56  kbk
!  removed tabs
!
!  Revision 1.7  2003/03/10 13:43:42  lars
!  double definitions removed - to conform with DEC compiler
!
!  Revision 1.6  2003/03/10 08:50:08  gotm
!  Improved documentation and cleaned up code
!
!  Revision 1.5  2002/02/08 08:33:44  gotm
!
!  Revision 1.4  2001/11/27 19:51:49  gotm
!  Cleaned
!  Revision 1.3  2001/11/27 15:38:06  gotm
!  Possible to read coordinate distribution from file
!
!  Revision 1.1.1.1  2001/02/12 15:55:57  gotm
!  initial import into CVS
!
!EOP
!
! !LOCAL VARIABLES:
integer                   :: i,rc,j,nlayers
REALTYPE                  :: zi(0:nlev)
integer, parameter        :: grid_unit = 101
!-----------------------------------------------------------------------
!BOC
if (.not. grid_ready) then ! Build up dimensionless grid (0<=ga<=1)
select case (grid_method)
case(0) !Equidistant grid with possible zooming to surface and bottom
LEVEL2 "sigma coordinates (zooming possible)"
if (ddu .le. 0 .and. ddl .le. 0) then
do i=1,nlev
ga(i)=ga(i-1)+1/float(nlev)
end do
else
do i=1,nlev ! This zooming routine is from Antoine Garapon, ICCH, DK
ga(i)=tanh((ddl+ddu)*i/nlev-ddl)+tanh(ddl)
ga(i)=ga(i)/(tanh(ddl)+tanh(ddu))
end do
end if
depth = depth0 + Zeta
do i=1,nlev
h(i)=(ga(i)-ga(i-1))*depth
end do
case(1) !Sigma, the fraction each layer occupies is specified.
LEVEL2 "external specified sigma coordinates"
open (grid_unit,FILE =grid_file,status='unknown',ERR=100)
if (nlayers /= nlev) then
FATAL "number of layers spefified in file <> # of model layers"
stop 'updategrid'
end if
depth = _ZERO_
j = 0
do i=nlev,1,-1 !The first layer to be read is at the surface
depth = depth + ga(i)
j=j+1
end do
if (j /= nlayers) then
FATAL "number of layers read from file <> # of model layers"
stop 'updategrid'
end if
close (grid_unit)
if (abs(depth-1.).gt.1.e-8) then
FATAL "sum of all layers in grid_file should be 1."
stop 'updategrid'
end if
case(2) !Cartesian, the layer thickness is read from file
LEVEL2 "external specified cartesian coordinates"
open (grid_unit,FILE =grid_file,ERR=100)
! Observations is called after meanflow is initialised, and we don#t have
! zeta_method
!        if (zeta_method /= 0) then
!          stop "You are using Cartesian coordinates with varying surface elevation"
!        end if
if(nlayers /= nlev) then
FATAL "nlev must be equal to the number of layers in: ", &
trim(grid_file)
stop 'updategrid'
end if
depth = _ZERO_
j=0
do i=nlev,1,-1 !The first layer read is the surface
depth = depth + h(i)
j=j+1
end do
if (j /= nlayers) then
FATAL "number of layers read from file <> # of model layers"
stop 'updategrid'
end if
close (grid_unit)
if (abs(depth-depth0).gt.1.e-5) then
FATAL "sum of all layers should be equal to the total depth",depth0,depth
stop 'updategrid'
end if
ga(0)=-1.
if (ddu.le.0.and.ddl.le.0) then
do i=1,nlev
ga(i)=ga(i-1)+1/float(nlev)
end do
else
do i=1,nlev ! This zooming is from Antoine Garapon, ICCH, DK
ga(i)=tanh((ddl+ddu)*i/nlev-ddl)+tanh(ddl)
ga(i)=ga(i)/(tanh(ddl)+tanh(ddu))-1.
end do
end if
depth = depth0 + Zeta
do i=1,nlev
h(i)  = (ga(i)-ga(i-1)) * depth
end do
case default
stop "updategrid: No valid grid_method specified"
end select

grid_ready=.true.  !  Grid is now initialised !
end if

depth = depth0 + zeta

select case(grid_method)
case (0)
do i=1,nlev
ho(i) = h(i)
h(i)  = (ga(i)-ga(i-1)) * depth
end do
case (1)
ho = h
h = ga *depth
case (2)
ho=h
case default
stop "updategrid: No valid grid_method specified"
end select

z(1)=-depth0+0.5*h(1)
do i=2,nlev
z(i)=z(i-1)+0.5*(h(i-1)+h(i))
end do

!  Vertical velocity calculation:

case(0)
case(1,2)
! linearly varying advection velocity with peak at "w_height"
zi(0)=-depth0
do i=1,nlev
zi(i)=zi(i-1)+h(i)
end do
if (w_height.lt.0.01*(zi(nlev)-zi(0))) w_height=0.5*(zi(0)-zi(nlev))
do i=1,nlev-1
if (zi(i).gt.w_height) then
else
end if
end do
w(0)    =_ZERO_
w(nlev) =_ZERO_
case default
end select

return

100 FATAL 'Unable to open ',trim(grid_file),' for reading'
stop 'updategrid'
101 FATAL 'Error reading grid file ',trim(grid_file)
stop 'updategrid'

end subroutine updategrid
!EOC

!-----------------------------------------------------------------------
! Copyright by the GOTM-team under the GNU Public License - www.gnu.org
!-----------------------------------------------------------------------