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: grid_ready
   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
   use observations, only: zeta_method,w_adv_method
   use observations, only: w_adv,w_height,w_adv_discr
   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
!  added new copyright to files
!
!  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
!  Manuel added support for reading grid distribution from file
!
!  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)
         read (grid_unit,*) nlayers
         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
            read(grid_unit,*,ERR=101,END=101) ga(i)
            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
         read (grid_unit,*) nlayers
         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
            read(grid_unit,*,ERR=101) h(i)
            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
     case(3) ! Adaptive grid
          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:

   select case(w_adv_method)
      case(0)
         ! no vertical advection
      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
               w(i)=(zi(nlev)-zi(i))/(zi(nlev)-w_height)*w_adv
            else
               w(i)=(zi(0)-zi(i))/(zi(0)-w_height)*w_adv
            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
!-----------------------------------------------------------------------