seagrass.F90 11.6 KB
Newer Older
dumoda01's avatar
dumoda01 committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411
!$Id: seagrass.F90,v 1.9 2007-01-06 11:49:15 kbk Exp $
#include"cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: seagrass --- sea grass dynamics \label{sec:seagrass}
!
! !INTERFACE:
   module seagrass
!
! !DESCRIPTION:
! In this module, seagrass canopies are treated as Lagrangian tracers,
! which either advect passively with the horizontal current speed or
! rest at their excursion limits and thus exert friction on the mean flow,
! see \cite{VerduinBackhaus2000}.
! Turbulence generation due to seagrass friction is possible, see
! namelist file {\tt seagrass.nml}. The extra production term
! in the balance of TKE, \eq{tkeA}, is included as described in
! \sect{sec:production}.

!
! !USES:
!
!  default: all is private.
   private
!
! !PUBLIC MEMBER FUNCTIONS:
   public init_seagrass, do_seagrass, save_seagrass, end_seagrass
   logical, public                     :: seagrass_calc=.false.
!
! !REVISION HISTORY:!
!  Original author(s): Hans Burchard & Karsten Bolding
!  $Log: seagrass.F90,v $
!  Revision 1.9  2007-01-06 11:49:15  kbk
!  namelist file extension changed .inp --> .nml
!
!  Revision 1.8  2006-12-03 13:54:22  hb
!  No extra production above seagrass
!
!  Revision 1.7  2006-11-21 15:21:56  kbk
!  seagrass working again
!
!  Revision 1.6  2005-12-02 21:10:25  hb
!  Documentation updated
!
!  Revision 1.5  2005/06/27 13:44:07  kbk
!  modified + removed traling blanks
!
!  Revision 1.4  2003/03/28 09:20:34  kbk
!  added new copyright to files
!
!  Revision 1.3  2003/03/28 08:28:36  kbk
!  removed tabs
!
!  Revision 1.2  2003/03/10 09:13:09  gotm
!  Improved documentation
!
!  Revision 1.1.1.1  2001/02/12 15:55:57  gotm
!  initial import into CVS
!
   REALTYPE, dimension(:), allocatable :: xx,yy
   REALTYPE, dimension(:), allocatable :: exc,vfric,grassz,xxP
!  from a namelist
   character(len=PATH_MAX)   :: grassfile='seagrass.dat'
   REALTYPE                  :: XP_rat
   integer                   :: grassind
   integer                   :: grassn
   integer                   :: out_unit
   integer                   :: maxn

!EOP
!-----------------------------------------------------------------------

   contains

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Initialise the sea grass module
!
! !INTERFACE:
   subroutine init_seagrass(namlst,fname,unit,nlev,h)
!
! !DESCRIPTION:
! Here, the seagrass namelist {\tt seagrass.nml} is read
! and memory is allocated
! for some relevant vectors. Afterwards, excursion limits and friction
! coefficients are read from a file. The uppermost grid related index
! for the seagrass canopy is then calculated.
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   integer,          intent(in)   :: namlst
   character(len=*), intent(in)   :: fname
   integer,          intent(in)   :: unit
   integer,          intent(in)   :: nlev
   REALTYPE,         intent(in)   :: h(0:nlev)
!
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
! !LOCAL VARIABLES:
   integer                   :: i,rc
   REALTYPE                  :: z
   namelist /canopy/  seagrass_calc,grassfile,XP_rat
!EOP
!-----------------------------------------------------------------------
!BOC
   LEVEL1 'init_seagrass'

   maxn=nlev

!  Open and read the namelist
   open(namlst,file=fname,action='read',status='old',err=98)
   read(namlst,nml=canopy,err=99)
   close(namlst)

   if (seagrass_calc) then
      out_unit=unit

      open(unit,status='unknown',file=grassfile)

      read(unit,*) grassn

      allocate(xx(0:nlev),stat=rc)
      if (rc /= 0) STOP 'init_seagrass: Error allocating (xx)'

      allocate(yy(0:nlev),stat=rc)
      if (rc /= 0) STOP 'init_seagrass: Error allocating (yy)'

      allocate(xxP(0:nlev),stat=rc)
      if (rc /= 0) STOP 'init_seagrass: Error allocating (xxP)'

      allocate(exc(0:grassn),stat=rc)
      if (rc /= 0) STOP 'init_seagrass: Error allocating (exc)'

      allocate(vfric(0:grassn),stat=rc)
      if (rc /= 0) STOP 'init_seagrass: Error allocating (vfric)'

      allocate(grassz(0:grassn),stat=rc)
      if (rc /= 0) STOP 'init_seagrass: Error allocating (grassz)'

      xx = _ZERO_
      yy = _ZERO_
      xxP= _ZERO_

      do i=1,grassn
         read(unit,*) grassz(i),exc(i),vfric(i)
      end do

      z=0.5*h(1)
      do i=2,nlev
         z=z+0.5*(h(i-1)+h(i))
         if (grassz(grassn).gt.z) grassind=i
      end do


      close(unit)

      LEVEL2 'seagrass initialised ...'

   end if
   return

98 LEVEL2 'I could not open seagrass.nml'
   LEVEL2 'Ill continue but set seagrass_calc to false.'
   LEVEL2 'If thats not what you want you have to supply seagrass.nml'
   LEVEL2 'See the Seagrass example on www.gotm.net for a working seagrass.nml'
   seagrass_calc = .false.
   return
99 FATAL 'I could not read seagrass.nml'
   stop 'init_seagrass'
   end subroutine init_seagrass
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Update the sea grass model
!
! !INTERFACE:
   subroutine do_seagrass(nlev,dt)
!
! !DESCRIPTION:
!
!  Here the time depending seagrass equation suggested by
!  \cite{VerduinBackhaus2000} is calculated. In order to
!  explain the basic principle, an idealised example is examined here
!  with a simplified momentum equation,
!
!  \begin{equation}
!  \partial_t u - \partial_z(\nu_t \partial_z u) = -g\partial_x\zeta-C_fu|u|
!  \comma
!  \end{equation}
!  and the Lagrangian tracer equation for seagrass,
!  \begin{equation}
!  \partial_t X =
!  \left\{
!  \begin{array}{ll}
!  u & \mbox{ for } |X|< X_{\max} \mbox{ or } X \cdot u <0,\\
!  0 & \mbox{ else}
!  \comma
!  \end{array}
!  \right.
!  \end{equation}
!  where $X$ is the Langrangian horizontal excursion of the seagrass.
!  The seagrass friction coefficient, $C_f$, is only non--zero at heights
!  where seagrass tracers are at their excursion limits:
!
!  \begin{equation}
!  C_f =
!  \left\{
!  \begin{array}{ll}
!  C_f^{\max} & \mbox{ for } |X|=X_{\max} \comma \\
!  0 & \mbox{ else} \point
!  \end{array}
!  \right.
!  \end{equation}
!
!  The maximum excursion limits $X_{\max}$ and the friction coefficients
!  $C_f^{\max}$ are read from a file.
!
!  The production of turbulence is calculated here as the sum of shear
!  production and friction loss at the seagrass leaves,
!  \begin{equation}
!   \label{sgProduction}
!    X_P = \alpha_{sg} C_f |u|^3
!   \comma
!  \end{equation}
!  which is added to the usual shear--production term as illustrated in
!  \eq{computeP}. The efficiency coefficient of turbulence production
!  by sea--grass friction, $\alpha_{sg}$, is denoted as {\tt xP\_rat}
!  in the code. It has to be read--in from the {\tt canopy} namelist.
!  For details and example calculations, see \cite{BurchardBolding2000}.
!
! !USES:
   use meanflow, only:     u,v,h,drag,xP
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   integer,  intent(in)     :: nlev
   REALTYPE, intent(in)     :: dt
!
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
! !LOCAL VARIABLES:
   integer                   :: i
   REALTYPE                  :: dist
   REALTYPE                  :: grassfric(0:nlev)
   REALTYPE                  :: excur(0:nlev)
   REALTYPE                  :: z(0:nlev)
!EOP
!-----------------------------------------------------------------------
   if (seagrass_calc) then

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

!     Interpolate excursion limits and friction to actual grid.

      call gridinterpol(grassn,1,grassz,exc,nlev,z,excur)
      call gridinterpol(grassn,1,grassz,vfric,nlev,z,grassfric)

      do i=1,grassind
         xx(i)=xx(i)+dt*u(i)                ! Motion of seagrass elements with
         yy(i)=yy(i)+dt*v(i)                ! mean flow.
         dist=sqrt(xx(i)*xx(i)+yy(i)*yy(i))
         if (dist .gt. excur(i)) then       ! Excursion limit reached
            xx(i)= excur(i)/dist * xx(i)
            yy(i)= excur(i)/dist * yy(i)

!           Increased drag by seagrass friction
            drag(i)=drag(i)+grassfric(i)

!           Extra turbulence production by seagrass friction
            xxP(i)=xP_rat*grassfric(i)*(sqrt(u(i)**2+v(i)**2))**3
         else
            xxP(i)=_ZERO_
         end if
      end do

!     Interpolate onto turbulence grid points
      do i=1,nlev-1
         xP(i)=0.5*(xxP(i)+xxP(i+1))
      end do

   end if
   return
   end subroutine do_seagrass
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Storing the results
!
! !INTERFACE:
   subroutine save_seagrass
!
! !DESCRIPTION:
!  Here, storing of the seagrass profiles to an ascii or a
!  netCDF file is managed.
!
! !USES:
   use meanflow, only:     h
   use output, only: out_fmt,ascii_unit,ts
#ifdef NETCDF_FMT
   use ncdfout, only: ncid
   use ncdfout, only: lon_dim,lat_dim,z_dim,time_dim,dims
   use ncdfout, only: define_mode,new_nc_variable,set_attributes,store_data
#endif
   IMPLICIT NONE

#ifdef NETCDF_FMT
#include "netcdf.inc"
#endif

!
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
! !LOCAL VARIABLES:
   logical, save             :: first=.true.
   integer, save             :: x_excur_id,y_excur_id,n
   integer                   :: i,iret
   REALTYPE                  :: zz
   REALTYPE                  :: miss_val
!EOP
!-----------------------------------------------------------------------
!BOC

   select case (out_fmt)
      case (ASCII)
         write(ascii_unit,*) 'Seagrass ',grassind
         write(ascii_unit,110) 'height','x-excur','y-excur'
         zz = _ZERO_
         do i=1,grassind
            zz=zz+0.5*h(i)
            write(ascii_unit,111) zz,xx(i),yy(i)
            zz=zz+0.5*h(i)
         end do
      case (NETCDF)
#ifdef NETCDF_FMT
         if(first) then
            dims(1) = lon_dim
            dims(2) = lat_dim
            dims(3) = z_dim
            dims(4) = time_dim
            miss_val = -999.0
            iret = define_mode(ncid,.true.)
            iret = new_nc_variable(ncid,'x-excur',NF_REAL,4,dims,x_excur_id)
            iret = set_attributes(ncid,x_excur_id,units='m',    &
                   long_name='seagrass excursion(x)',missing_value=miss_val)
            iret = new_nc_variable(ncid,'y-excur',NF_REAL,4,dims,y_excur_id)
            iret = set_attributes(ncid,y_excur_id,units='m',    &
                   long_name='seagrass excursion(y)',missing_value=miss_val)
            iret = define_mode(ncid,.false.)
            n = ubound(xx,1)
            first = .false.
         end if
         xx(grassind+1:maxn)=miss_val
         yy(grassind+1:maxn)=miss_val
         iret = store_data(ncid,x_excur_id,XYZT_SHAPE,n,array=xx)
         iret = store_data(ncid,y_excur_id,XYZT_SHAPE,n,array=yy)
#endif
      case default
         FATAL 'A non valid output format has been chosen'
         stop 'save_seagrass'
   end select
110 format(3(1x,A12))
111 format(3(1x,E12.3E2))
   return
   end subroutine save_seagrass
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Finish the sea grass calculations
!
! !INTERFACE:
   subroutine end_seagrass
!
! !DESCRIPTION:
!  Nothing done yet --- supplied for completeness.
!
! !USES:
   IMPLICIT NONE
!
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
!EOP
!-----------------------------------------------------------------------
!BOC

   return
   end subroutine end_seagrass
!EOC
!-----------------------------------------------------------------------

   end module seagrass

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