m2d.F90 12.4 KB
Newer Older
gotm's avatar
gotm committed
1 2 3 4 5 6 7 8 9 10
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: m2d - depth integrated hydrodynamical model (2D)
!
! !INTERFACE:
   module m2d
!
! !DESCRIPTION:
hb's avatar
hb committed
11
!  This module contains declarations for all variables related to 2D
gotm's avatar
gotm committed
12
!  hydrodynamical calculations. Information about the calculation domain
hb's avatar
hb committed
13
!  is included from the {\tt domain} module.
gotm's avatar
gotm committed
14 15
!  The module contains public subroutines for initialisation, integration
!  and clean up of the 2D model component.
hb's avatar
hb committed
16 17 18
!  The actual calculation routines are called in {\tt integrate\_2d}
!  and are linked
!  in from the library {\tt lib2d.a}.
gotm's avatar
gotm committed
19 20
!
! !USES:
bjb's avatar
bjb committed
21
   use exceptions
kbk's avatar
kbk committed
22
   use time, only: julianday,secondsofday
gotm's avatar
gotm committed
23
   use parameters, only: avmmol
kbk's avatar
kbk committed
24
   use domain, only: imin,imax,jmin,jmax,az,au,av,H,HU,HV,min_depth
kb's avatar
kb committed
25 26
   use domain, only: ilg,ihg,jlg,jhg
   use domain, only: ill,ihl,jll,jhl
27
   use domain, only: openbdy,z0_method,z0_const,z0
bjb's avatar
bjb committed
28
   use domain, only: az,ax
kbk's avatar
kbk committed
29 30
   use halo_zones, only : update_2d_halo,wait_halo
   use halo_zones, only : U_TAG,V_TAG,H_TAG
kbk's avatar
kbk committed
31
   use variables_2d
gotm's avatar
gotm committed
32
   IMPLICIT NONE
bjb's avatar
bjb committed
33 34 35 36 37 38 39 40
! Temporary interface (should be read from module):
   interface
      subroutine get_2d_field(fn,varname,il,ih,jl,jh,f)
         character(len=*),intent(in)   :: fn,varname
         integer, intent(in)           :: il,ih,jl,jh
         REALTYPE, intent(out)         :: f(:,:)
      end subroutine get_2d_field
   end interface
gotm's avatar
gotm committed
41 42
!
! !PUBLIC DATA MEMBERS:
kbk's avatar
kbk committed
43
   logical                   :: have_boundaries
bjb's avatar
bjb committed
44 45 46 47 48 49 50
   REALTYPE                  :: dtm,Am=-_ONE_
!  method for specifying horizontal numerical diffusion coefficient
!     (0=const, 1=from named nc-file)
   integer                   :: An_method=0
   REALTYPE                  :: An_const=-_ONE_
   character(LEN = PATH_MAX) :: An_file

kbk's avatar
kbk committed
51
   integer                   :: MM=1,residual=-1
52
   integer                   :: sealevel_check=0
kbk's avatar
kbk committed
53 54 55
   logical                   :: bdy2d=.false.
   integer                   :: bdyfmt_2d,bdytype,bdyramp_2d=-1
   character(len=PATH_MAX)   :: bdyfile_2d
kbk's avatar
kbk committed
56
   REAL_4B                   :: bdy_data(1500)
kb's avatar
kb committed
57 58
   REAL_4B                   :: bdy_data_u(1500)
   REAL_4B                   :: bdy_data_v(1500)
kbk's avatar
kbk committed
59 60
   REAL_4B, allocatable      :: bdy_times(:)
   integer, parameter        :: comm_method=-1
gotm's avatar
gotm committed
61
!
kbk's avatar
kbk committed
62 63 64
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
gotm's avatar
gotm committed
65
! !LOCAL VARIABLES:
bjb's avatar
bjb committed
66 67
  integer                    :: num_neighbors
  REALTYPE                   :: An_sum
gotm's avatar
gotm committed
68 69 70 71 72 73 74 75 76
!
!EOP
!-----------------------------------------------------------------------

   contains

!-----------------------------------------------------------------------
!BOP
!
hb's avatar
hb committed
77
! !IROUTINE: init_2d - initialise 2D related stuff.
gotm's avatar
gotm committed
78 79 80 81 82 83
!
! !INTERFACE:
   subroutine init_2d(runtype,timestep,hotstart)
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
kbk's avatar
kbk committed
84 85 86
   integer, intent(in)                 :: runtype
   REALTYPE, intent(in)                :: timestep
   logical, intent(in)                 :: hotstart
gotm's avatar
gotm committed
87 88 89 90 91 92
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !DESCRIPTION:
hb's avatar
hb committed
93 94 95 96 97 98
!  Here, the {\tt m2d} namelist is read from {\tt getm.inp}, and the check
!  for the fulfilment of the CFL criterium for shallow water theory
!  {\tt cfl\_check} is called. A major part of this subroutine deals
!  then with the setting of local bathymetry values and initial surface
!  elevations in $u$- and $v$-points, also by calls to the subroutines 
!  {\tt uv\_depths} and {\tt depth\_update}.
gotm's avatar
gotm committed
99 100
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
101
   integer                   :: rc
kbk's avatar
kbk committed
102
   integer                   :: i,j
103
   integer                   :: vel_depth_method=0
kbk's avatar
kbk committed
104
   namelist /m2d/ &
bjb's avatar
bjb committed
105 106
          MM,vel_depth_method,Am,An_method,An_const,An_file,residual, &
          sealevel_check,bdy2d,bdyfmt_2d,bdyramp_2d,bdyfile_2d
gotm's avatar
gotm committed
107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'init_2d() # ',Ncall
#endif

   LEVEL1 'init_2d'

!  Read 2D-model specific things from the namelist.
   read(NAMLST,m2d)

   dtm = timestep

!  Allocates memory for the public data members - if not static
   call init_variables_2d(runtype)

126 127
#if defined(GETM_PARALLEL) || defined(NO_BAROTROPIC) 
!   STDERR 'Not calling cfl_check() - GETM_PARALLEL or NO_BAROTROPIC'
lars's avatar
lars committed
128
!   call cfl_check()
kbk's avatar
kbk committed
129
#else
130
   call cfl_check()
kbk's avatar
kbk committed
131
#endif
gotm's avatar
gotm committed
132 133

   if (Am .lt. _ZERO_) then
134 135
      LEVEL2 'Am < 0 --> horizontal momentum diffusion not included'
   end if
bjb's avatar
bjb committed
136 137 138 139

   select case (An_method)
      case(0)
         LEVEL2 'An_method=0 -> horizontal numerical diffusion not included'
kb's avatar
kb committed
140
         An = _ZERO_
bjb's avatar
bjb committed
141 142 143 144 145
      case(1)
         LEVEL2 'An_method=1 -> Using constant horizontal numerical diffusion'
         if (An_const .lt. _ZERO_) then
              call getm_error("init_2d()", &
                         "Constant horizontal numerical diffusion <0");
kb's avatar
kb committed
146
         else
bjb's avatar
bjb committed
147 148
            An  = An_const
            AnX = An_const
bjb's avatar
bjb committed
149 150 151
         end if
      case(2)
         LEVEL2 'An_method=2 -> Using space varying horizontal numerical diffusion'
kb's avatar
kb committed
152
         LEVEL2 '..  will read An from An_file ',trim(An_file)
bjb's avatar
bjb committed
153 154
         ! Initialize and then read field:
         An = _ZERO_
bjb's avatar
bjb committed
155
         call get_2d_field(trim(An_file),"An",ilg,ihg,jlg,jhg,An(ill:ihl,jll:jhl))
kb's avatar
kb committed
156 157
         call update_2d_halo(An,An,az,imin,jmin,imax,jmax,H_TAG)
         call wait_halo(H_TAG)
bjb's avatar
bjb committed
158 159
         ! Compute AnX (An in X-points) based on An and the X- and T- masks
         AnX = _ZERO_
Bjarne Buchmann's avatar
Bjarne Buchmann committed
160 161 162
         ! We loop over the X-points in the present domain.
         do j=jmin-1,jmax
            do i=imin-1,imax
bjb's avatar
bjb committed
163 164 165 166
               if (ax(i,j) .ge. 1) then
                  num_neighbors = 0
                  An_sum = _ZERO_
                  ! Each AnX should have up to 4 T-point neighbours.
Bjarne Buchmann's avatar
Bjarne Buchmann committed
167 168
                  if ( az(i,j) .ge. 1 ) then
                     An_sum        = An_sum + An(i,j)
bjb's avatar
bjb committed
169 170
                     num_neighbors = num_neighbors +1
                  end if
Bjarne Buchmann's avatar
Bjarne Buchmann committed
171 172
                  if ( az(i,j+1) .ge. 1 ) then
                     An_sum        = An_sum + An(i,j+1)
bjb's avatar
bjb committed
173 174
                     num_neighbors = num_neighbors +1
                  end if
Bjarne Buchmann's avatar
Bjarne Buchmann committed
175 176
                  if ( az(i+1,j) .ge. 1 ) then
                     An_sum        = An_sum + An(i+1,j)
bjb's avatar
bjb committed
177 178
                     num_neighbors = num_neighbors +1
                  end if
Bjarne Buchmann's avatar
Bjarne Buchmann committed
179 180
                  if ( az(i+1,j+1) .ge. 1 ) then
                     An_sum        = An_sum + An(i+1,j+1)
bjb's avatar
bjb committed
181 182 183 184 185 186 187 188 189 190 191
                     num_neighbors = num_neighbors +1
                  end if
                  ! Take average of actual neighbours:
                  if (num_neighbors .gt. 0) then
                     AnX(i,j) = An_sum/num_neighbors
                  end if
               end if
            end do
         end do
         call update_2d_halo(AnX,AnX,ax,imin,jmin,imax,jmax,H_TAG)
         call wait_halo(H_TAG)
192 193 194 195 196 197 198 199
         !
         ! If all An values are really zero, then we should not use An-smoothing at all...
         ! Note that smoothing may be on in other subdomains.
         if ( MAXVAL(ABS(An)) .eq. _ZERO_ ) then
            LEVEL2 '  All An is zero for this (sub)domain - switching to An_method=0'
            An_method=0
         end if
            
bjb's avatar
bjb committed
200 201 202 203
      case default
         call getm_error("init_2d()", &
                         "A non valid An method has been chosen");
   end select
gotm's avatar
gotm committed
204

205 206 207 208 209 210 211 212
   if (sealevel_check .eq. 0) then
      LEVEL2 'sealevel_check=0 --> NaN checks disabled'
   else if (sealevel_check .gt. 0) then
      LEVEL2 'sealevel_check>0 --> NaN values will result in error conditions'
   else
      LEVEL2 'sealevel_check<0 --> NaN values will result in warnings'
   end if

kb's avatar
kb committed
213
   if (.not. openbdy)  bdy2d=.false.
gotm's avatar
gotm committed
214 215
   LEVEL2 'Open boundary=',bdy2d
   if (bdy2d) then
216 217 218 219
      if (hotstart .and. bdyramp_2d .gt. 0) then
          LEVEL2 'WARNING: hotstart is .true. AND bdyramp_2d .gt. 0'
          LEVEL2 'WARNING: .. be sure you know what you are doing ..'
      end if
gotm's avatar
gotm committed
220 221 222 223
      LEVEL2 TRIM(bdyfile_2d)
      LEVEL2 'Format=',bdyfmt_2d
   end if

224
   call uv_depths(vel_depth_method)
gotm's avatar
gotm committed
225

kbk's avatar
kbk committed
226 227 228 229
   where ( -H+min_depth .gt. _ZERO_ )
      z = -H+min_depth
   end where
   zo=z
gotm's avatar
gotm committed
230

kbk's avatar
kbk committed
231 232 233
   where (-HU+min_depth .gt. _ZERO_ )
      zu = -HU+min_depth
   end where
gotm's avatar
gotm committed
234

kbk's avatar
kbk committed
235 236 237
   where (-HV+min_depth .gt. _ZERO_ )
      zv = -HV+min_depth
   end where
kbk's avatar
kbk committed
238 239 240 241 242 243 244

!  bottom roughness
   if (z0_method .eq. 0) then
      zub0 = z0_const
      zvb0 = z0_const
   end if
   if (z0_method .eq. 1) then
245 246
      do j=jmin-HALO,jmax+HALO
         do i=imin-HALO,imax+HALO-1
kbk's avatar
kbk committed
247
           if (au(i,j) .gt. 0) zub0(i,j) = 0.5*(z0(i,j)+z0(i+1,j))
248 249 250 251
         end do
      end do
      do j=jmin-HALO,jmax+HALO-1
         do i=imin-HALO,imax+HALO
kbk's avatar
kbk committed
252 253 254 255
           if (av(i,j) .gt. 0) zvb0(i,j) = 0.5*(z0(i,j)+z0(i,j+1))
         end do
      end do
   end if
kbk's avatar
kbk committed
256 257
   zub=zub0
   zvb=zvb0
gotm's avatar
gotm committed
258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275

   call depth_update()

#ifdef DEBUG
   write(debug,*) 'Leaving init_2d()'
   write(debug,*)
#endif
   return
   end subroutine init_2d
!EOC

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: integrate_2d - sequence of calls to do 2D model integration
!
! !INTERFACE:
   subroutine integrate_2d(runtype,loop,tausx,tausy,airp)
bjb's avatar
bjb committed
276
   use getm_timers, only: tic, toc, TIM_INTEGR2D
gotm's avatar
gotm committed
277 278 279
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
kbk's avatar
kbk committed
280 281 282 283
   integer, intent(in)                 :: runtype,loop
   REALTYPE, intent(in)                :: tausx(E2DFIELD)
   REALTYPE, intent(in)                :: tausy(E2DFIELD)
   REALTYPE, intent(in)                :: airp(E2DFIELD)
gotm's avatar
gotm committed
284 285 286 287 288 289
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !DESCRIPTION:
hb's avatar
hb committed
290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309
!  Here, all 2D related subroutines are called. The major calls and their
!  meaning are:
!
!  \vspace{0.5cm}
!
!  \begin{tabular}{ll}
!  {\tt call update\_2d\_bdy} & read in new lateral boundary conditions \\
!  {\tt call bottom\_friction} & update bottom friction\\
!  {\tt call uv\_advect} & calculate 2D advection terms\\
!  {\tt call uv\_diffusion} & calculate 2D  diffusion terms\\
!  {\tt call momentum} & iterate 2D momemtum equations\\
!  {\tt call sealevel} & update sea surface elevation\\
!  {\tt call depth\_update} & update water depths\\
!  {\tt call do\_residual} & calculate intermdediate values for residual currents
!  \end{tabular}
!
!  \vspace{0.5cm}
!
!  It should be noted that some of these calls may be excluded for certain
!  compiler options set in the {\tt Makefile} of the application.
gotm's avatar
gotm committed
310 311 312 313 314 315 316 317 318 319 320 321 322 323 324
!
! !LOCAL VARIABLES:
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'integrate_2d() # ',Ncall
#endif
   if (mod(loop-1,MM) .eq. 0) then        ! MacroMicro time step
#ifndef NO_BOTTFRIC
      call bottom_friction(runtype)
#endif
   end if
325
   UEx=_ZERO_ ; VEx=_ZERO_
gotm's avatar
gotm committed
326 327 328 329 330
#ifdef NO_ADVECT
   STDERR 'NO_ADVECT 2D'
#else
#ifndef UV_ADV_DIRECT
   call uv_advect()
bjb's avatar
bjb committed
331
   if (Am .gt. _ZERO_ .or. An_method .gt. 0) then
bjb's avatar
bjb committed
332
      call uv_diffusion(Am,An_method,An,AnX) ! Has to be called after uv_advect.
gotm's avatar
gotm committed
333
   end if
bjb's avatar
bjb committed
334
   call tic(TIM_INTEGR2D)
kbk's avatar
kbk committed
335 336
   call mirror_bdy_2d(UEx,U_TAG)
   call mirror_bdy_2d(VEx,V_TAG)
bjb's avatar
bjb committed
337
   call toc(TIM_INTEGR2D)
gotm's avatar
gotm committed
338 339 340 341
#endif
#endif
   call momentum(loop,tausx,tausy,airp)
   if (runtype .gt. 1) then
bjb's avatar
bjb committed
342
      call tic(TIM_INTEGR2D)
gotm's avatar
gotm committed
343 344
      Uint=Uint+U
      Vint=Vint+V
bjb's avatar
bjb committed
345
      call toc(TIM_INTEGR2D)
gotm's avatar
gotm committed
346
   end if
kbk's avatar
kbk committed
347
   if (have_boundaries) call update_2d_bdy(loop,bdyramp_2d)
348 349 350
   call sealevel()
   call depth_update()

gotm's avatar
gotm committed
351
   if(residual .gt. 0 .and. loop .ge. residual) then
bjb's avatar
bjb committed
352
      call tic(TIM_INTEGR2D)
gotm's avatar
gotm committed
353
      call do_residual(0)
bjb's avatar
bjb committed
354
      call toc(TIM_INTEGR2D)
gotm's avatar
gotm committed
355 356 357 358 359 360 361 362 363 364
   end if

#ifdef DEBUG
     write(debug,*) 'Leaving integrate_2d()'
     write(debug,*)
#endif
   return
   end subroutine integrate_2d
!EOC

bjb's avatar
bjb committed
365

gotm's avatar
gotm committed
366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: clean_2d - cleanup after 2D run.
!
! !INTERFACE:
   subroutine clean_2d()
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !DESCRIPTION:
hb's avatar
hb committed
382 383
!  This routine executes a final call to {\tt do\_residual} where the residual
!  current calculations are finished.
gotm's avatar
gotm committed
384 385 386 387 388 389 390 391 392 393 394 395
!
! !LOCAL VARIABLES:
!
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'clean_2d() # ',Ncall
#endif

kbk's avatar
kbk committed
396 397 398
   if(residual .gt. 0) then
      call do_residual(1)
   end if
gotm's avatar
gotm committed
399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414

#ifdef DEBUG
   write(debug,*) 'Leaving clean_2d()'
   write(debug,*)
#endif
   return
   end subroutine clean_2d
!EOC

!-----------------------------------------------------------------------

   end module m2d

!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding (BBH)         !
!-----------------------------------------------------------------------