m2d.F90 14.3 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
!  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
97
!  elevations in $u$- and $v$-points, also by calls to the subroutines
hb's avatar
hb committed
98
!  {\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
#if defined(GETM_PARALLEL) || defined(NO_BAROTROPIC)
127
!   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
         !
         ! 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
199

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

   call depth_update()

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

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
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: postinit_2d - re-initialise some 2D after hotstart read.
!
! !INTERFACE:
   subroutine postinit_2d(runtype)
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   integer, intent(in)                 :: runtype
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !DESCRIPTION:
!  This routine provides possibility to reset/initialize 2D variables to 
!  ensure that velocities are correctly set on land cells after read 
!  of a hotstart file.
!
! !LOCAL VARIABLES:
   integer                   :: i,j, ischange
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'postinit_2d() # ',Ncall
#endif

   LEVEL1 'postinit_2d'

   ischange = 0

! The first two loops are pure diagnostics, logging where changes will actually take place
! (and if there is something to do at all, to be able to skip the second part)
   do j=jmin-HALO,jmax+HALO
      do i=imin-HALO,imax+HALO
         if ( au(i,j).eq.0 .and. U(i,j).ne._ZERO_ ) then
            LEVEL3 'postinit_2d: Reset to mask(au), U=0 for i,j=',i,j
            ischange = 1
         end if
      end do
   end do
   do j=jmin-HALO,jmax+HALO
      do i=imin-HALO,imax+HALO
         if ( av(i,j).eq.0 .and. V(i,j).ne._ZERO_ ) then
            LEVEL3 'postinit_2d: Reset to mask(av), V=0 for i,j=',i,j
            ischange = 1
         end if
      end do
   end do

! The actual reset is below here - independent of the above diagnostics (except for the if)
   if (ischange.ne.0) then
      where (au .eq. 0)
         U     = _ZERO_
         Uint  = _ZERO_
         Uinto = _ZERO_
      end where
      where (av .eq. 0)
         V     = _ZERO_
         Vint  = _ZERO_
         Vinto = _ZERO_
      end where
      where (az .eq. 0)
         z  = _ZERO_
         zo = _ZERO_
      end where
   end if

   return
   end subroutine postinit_2d
!EOC

gotm's avatar
gotm committed
346 347 348 349 350 351 352
!-----------------------------------------------------------------------
!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
353
   use getm_timers, only: tic, toc, TIM_INTEGR2D
gotm's avatar
gotm committed
354 355 356
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
kbk's avatar
kbk committed
357 358 359 360
   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
361 362 363 364 365 366
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !DESCRIPTION:
hb's avatar
hb committed
367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386
!  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
387 388 389 390 391 392 393 394 395 396 397 398 399 400 401
!
! !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
402
   UEx=_ZERO_ ; VEx=_ZERO_
gotm's avatar
gotm committed
403 404 405 406 407
#ifdef NO_ADVECT
   STDERR 'NO_ADVECT 2D'
#else
#ifndef UV_ADV_DIRECT
   call uv_advect()
bjb's avatar
bjb committed
408
   if (Am .gt. _ZERO_ .or. An_method .gt. 0) then
bjb's avatar
bjb committed
409
      call uv_diffusion(Am,An_method,An,AnX) ! Has to be called after uv_advect.
gotm's avatar
gotm committed
410
   end if
bjb's avatar
bjb committed
411
   call tic(TIM_INTEGR2D)
kbk's avatar
kbk committed
412 413
   call mirror_bdy_2d(UEx,U_TAG)
   call mirror_bdy_2d(VEx,V_TAG)
bjb's avatar
bjb committed
414
   call toc(TIM_INTEGR2D)
gotm's avatar
gotm committed
415 416 417 418
#endif
#endif
   call momentum(loop,tausx,tausy,airp)
   if (runtype .gt. 1) then
bjb's avatar
bjb committed
419
      call tic(TIM_INTEGR2D)
gotm's avatar
gotm committed
420 421
      Uint=Uint+U
      Vint=Vint+V
bjb's avatar
bjb committed
422
      call toc(TIM_INTEGR2D)
gotm's avatar
gotm committed
423
   end if
kbk's avatar
kbk committed
424
   if (have_boundaries) call update_2d_bdy(loop,bdyramp_2d)
425 426 427
   call sealevel()
   call depth_update()

gotm's avatar
gotm committed
428
   if(residual .gt. 0 .and. loop .ge. residual) then
bjb's avatar
bjb committed
429
      call tic(TIM_INTEGR2D)
gotm's avatar
gotm committed
430
      call do_residual(0)
bjb's avatar
bjb committed
431
      call toc(TIM_INTEGR2D)
gotm's avatar
gotm committed
432 433 434 435 436 437 438 439 440 441
   end if

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

bjb's avatar
bjb committed
442

gotm's avatar
gotm committed
443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458
!-----------------------------------------------------------------------
!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
459 460
!  This routine executes a final call to {\tt do\_residual} where the residual
!  current calculations are finished.
gotm's avatar
gotm committed
461 462 463 464 465 466 467 468 469 470 471 472
!
! !LOCAL VARIABLES:
!
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'clean_2d() # ',Ncall
#endif

kbk's avatar
kbk committed
473 474 475
   if(residual .gt. 0) then
      call do_residual(1)
   end if
gotm's avatar
gotm committed
476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491

#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)         !
!-----------------------------------------------------------------------