advection_3d.F90 13.4 KB
Newer Older
gotm's avatar
gotm committed
1 2 3 4 5 6 7 8 9 10 11
!$id: advection_3d.F90,v 1.18 2001/09/19 13:53:08 bbh Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !MODULE:  3D advection
!
! !INTERFACE:
   module advection_3d
!
! !DESCRIPTION:
12
!
gotm's avatar
gotm committed
13 14 15 16
!  This module do advection of scalars.  The module follows the same
!  convention as the other modules in 'getm'. The module is initialised
!  by calling 'init\_advection\_3d()'. In the time-loop 'do\_advection\_3d' is
!  called. 'do\_advection\_3d' is a wrapper routine which - dependent on the
17 18 19
!  actual advection scheme chosen - makes calls to the appropriate 
!  subroutines, which may be done as one-step or multiple-step schemes.
!  The actual subroutines are coded in external FORTRAN files.
gotm's avatar
gotm committed
20 21 22
!  New advection schemes are easily implemented - at least from a program
!  point of view - since only this module needs to be changed.
!  Additional work arrays can easily be added following the stencil given
23 24 25 26 27 28 29 30 31
!  below. To add a new advection scheme three things must be done: 
!
!  \begin{enumerate}
!  \item define
!  a unique constant to identify the scheme (see e.g.\ {\tt UPSTREAM}
!  and {\tt TVD})
!  \item adopt the {\tt select case} in {\tt do\_advection\_3d} and 
!  \item  write the actual subroutine.
!  \end{enumerate}
gotm's avatar
gotm committed
32 33 34 35
!
! !USES:
   use domain, only: imin,imax,jmin,jmax
   use domain, only: iimin,iimax,jjmin,jjmax,kmax
36
   use halo_zones, only: update_3d_halo,wait_halo,D_TAG
gotm's avatar
gotm committed
37 38 39 40 41 42 43
   IMPLICIT NONE
!
   private
!
! !PUBLIC DATA MEMBERS:
   public init_advection_3d, do_advection_3d
#ifdef STATIC
44
   REALTYPE, public                    :: cu(I3DFIELD)
kbk's avatar
kbk committed
45 46
   REALTYPE, public                    :: hi(I3DFIELD)
   REALTYPE, public                    :: hio(I3DFIELD)
gotm's avatar
gotm committed
47
#else
48
   REALTYPE, public, dimension(:,:,:), allocatable       :: hi,hio,cu
gotm's avatar
gotm committed
49
#endif
50 51 52 53 54 55
   integer, public, parameter          :: UPSTREAM=1,UPSTREAM_SPLIT=2,P2=3
   integer, public, parameter          :: Superbee=4,MUSCL=5,P2_PDM=6,FCT=7
   REALTYPE, public, parameter         :: one6th=1./6.
   REALTYPE, public, parameter         :: ONE=_ONE_,TWO=2.*_ONE_
!
! !PRIVATE DATA MEMBERS:
gotm's avatar
gotm committed
56
!
kbk's avatar
kbk committed
57 58 59
! !REVISION HISTORY:
!  Original author(s): Karsten Bolding & Hans Burchard
!
gotm's avatar
gotm committed
60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
! !LOCAL VARIABLES:
   integer :: advection_method
!EOP
!-----------------------------------------------------------------------

   contains

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: init_advection_3d
!
! !INTERFACE:
   subroutine init_advection_3d(method)
!
! !DESCRIPTION:
76 77 78
!
! Here, memory for some variables is allocated, which are then initialised to
! zero.
gotm's avatar
gotm committed
79 80 81 82 83
!
! !USES
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
kbk's avatar
kbk committed
84
   integer, intent(in)                 :: method
gotm's avatar
gotm committed
85 86 87 88 89 90
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
91
   integer                   :: rc
gotm's avatar
gotm committed
92 93 94 95 96 97 98 99 100 101 102
!EOP
!-------------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'init_advection_3d() # ',Ncall
#endif

   LEVEL1 'init_advection_3d()'

kbk's avatar
kbk committed
103
#ifndef STATIC
gotm's avatar
gotm committed
104 105 106 107 108 109 110 111 112 113
   allocate(cu(I3DFIELD),stat=rc)    ! work array
   if (rc /= 0) stop 'init_advection_3d: Error allocating memory (cu)'

   allocate(hi(I3DFIELD),stat=rc)    ! work array
   if (rc /= 0) stop 'init_advection_3d: Error allocating memory (hi)'

   allocate(hio(I3DFIELD),stat=rc)    ! work array
   if (rc /= 0) stop 'init_advection_3d: Error allocating memory (hio)'
#endif

kbk's avatar
kbk committed
114
   cu  = _ZERO_ ; hi  = _ZERO_ ; hio = _ZERO_ 
gotm's avatar
gotm committed
115 116 117 118 119 120 121 122 123 124 125 126

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

!-----------------------------------------------------------------------
!BOP
!
127
! !IROUTINE:  do_advection_3d - 3D advection schemes \label{sec-do-advection-3d}
gotm's avatar
gotm committed
128 129 130
!
! !INTERFACE:
   subroutine do_advection_3d(dt,f,uu,vv,ww,hun,hvn,ho,hn,      &
kbk's avatar
kbk committed
131
                             delxu,delxv,delyu,delyv,area_inv,  &
132
                             az,au,av,hor_adv,ver_adv,adv_split,AH)
gotm's avatar
gotm committed
133 134 135
!
! !DESCRIPTION:
!
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
! Here, advection terms for all three-dimensional state variables are
! calculated by means of a finite-volume approach (an exception
! is the possibility to directly calculate the momentum advection
! by a one-step three-dimensional upstream scheme, 
! see {\tt uv\_advection\_3d}) and the advection step is carried out
! as a fractional advection time step. Those 3D variables may be defined on
! T-, U-, V- and W-points. The latter option is interesting for 
! turbulent quantities, but is not coded yet. Inside this advection
! routine, it does not matter, wehre the advected variable is located
! on the grid. All finite volume fluxes and geometric coefficients
! need to be calculated before {\tt advection\_3d} is called.
!
! Originally, this 3D advection routine has been written for tracer
! equations. There,
! after multiplying the layer-integrated and transformed to
! curvilinear coordinates tracer equation (\ref{C_Layer_IntCurvi}) 
! with $mn$, the advective 
! terms in this equation are discretised as follows. 
! 
! First advection term in (\ref{C_Layer_IntCurvi}):
! \begin{equation}\label{u_discr_advect}
! \left(mn\,\partial_{\cal X}\left(\frac{p_kc_k}{n}\right)\right)_{i,j}\approx 
! \frac{
! p_{i,j,k}\tilde c^u_{i,j,k}\Delta y^u_{i,j}-
! p_{i-1,j,k}\tilde c^u_{i-1,j,k}\Delta y^u_{i-1,j}
! }{\Delta x^c_{i,j}\Delta y^c_{i,j}}
! \end{equation}
! 
! Second advection term in (\ref{C_Layer_IntCurvi}):
! \begin{equation}\label{v_discr_advect}
! \left(mn\,\partial_{\cal Y}\left(\frac{q_kc_k}{m}\right)\right)_{i,j}\approx 
! \frac{
! q_{i,j,k}\tilde c^v_{i,j,k}\Delta y^v_{i,j}-
! q_{i,j-1,k}\tilde c^v_{i,j-1,k}\Delta y^v_{i,j-1}
! }{\Delta x^c_{i,j}\Delta y^c_{i,j}}
! \end{equation}
! 
! Vertical advective fluxes in (\ref{C_Layer_IntCurvi}):
! \begin{equation}\label{w_discr_advect}
! \left(\bar w_{k} \tilde c_{k}\right)_{i,j}\approx
! w_{i,j,k}\tilde c^w_{i,j,k}. 
! \end{equation}
! 
! The interfacial concentrations $\tilde c_{i,j,k}$ are calculated 
! according to upwind or higher order directional split
! schemes, which are discussed in detail below and in sections 
! \ref{sec-upstream-adv} - \ref{sec-fct-2dh-adv}. 
! 
! However, as said above, in the same way these routines may be applied 
! to quantities on
! U-, V-, and W-points, if the transports and geometric coefficients 
! are properly calculated.
! 
! There are various combinations of advection schemes possible. 
! The first selection is whether a one-step 3D first-order upstream method
! is cosen, or a fractional step method. 
!
! The next selection is (if a fractional step method is selected) 
! how to do the fractional steps (selection on {\tt adv\_split}). There
! are different options, 
!
! \begin{enumerate}
! \item directional split with subsequent full steps in $x$-, $y$- and 
! $z$-direction, 
! \item split with subsequent half steps in $x$-, and $y$-direction, a
! full step in $z$-direction, and half steps in $y$- and $x$-direction.
! \item directional split into a 2D horizontal step and a 1D vertical step.
! \end{enumerate}
!
! For the 1D directional-split schemes, first-order upstream, 
! ULTIMATE QUICKEST, and the Total Variation Diminishing (TVD) schemes
! Superbee, MUSCL, and P$_2$PDM are available. 

! For the 2D horizontal step, an upstream scheme and a Flux-Corrected
! Transport (FCT) scheme have been coded.
! 
! If the compiler option {\tt ITERATE\_VERT\_ADV} is chosen, the vertical
! advection is iterated as many times with reduced time step that 
! the CFL criterium for vertical advection is fulfilled, see the routine
! {\tt w\_split\_it\_adv}.
!
! With the compiler option {\tt SLICE\_MODEL}, the advection in
! $y$-direction is not executed.
!
!
gotm's avatar
gotm committed
221 222 223 224
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244
   REALTYPE, intent(in) :: uu(I3DFIELD)       ! layer-integrated x-transport
   REALTYPE, intent(in) :: vv(I3DFIELD)       ! layer-integrated y-transport
   REALTYPE, intent(in) :: ww(I3DFIELD)       ! grid-related vertical velocity
   REALTYPE, intent(in) :: ho(I3DFIELD)       ! old height of finite volume box
   REALTYPE, intent(in) :: hn(I3DFIELD)       ! new height of finite volume box
   REALTYPE, intent(in) :: hun(I3DFIELD)      ! height of x-interfaces
   REALTYPE, intent(in) :: hvn(I3DFIELD)      ! height of y-interfaces
   REALTYPE, intent(in) :: delxu(I2DFIELD)    ! dx centered on u-transport pt. 
   REALTYPE, intent(in) :: delxv(I2DFIELD)    ! length of y-interface
   REALTYPE, intent(in) :: delyu(I2DFIELD)    ! length of u-interface 
   REALTYPE, intent(in) :: delyv(I2DFIELD)    ! dy centered on v-transport pt.
   REALTYPE, intent(in) :: area_inv(I2DFIELD) ! inverse of horizontal box area
   REALTYPE, intent(in) :: dt                 ! advection time step
   REALTYPE, intent(in) :: AH                 ! constant horizontal diffusivity 
   integer, intent(in)  :: az(E2DFIELD)       ! mask for box centre (1: water)
   integer, intent(in)  :: au(E2DFIELD)       ! mask for u-transport (1: water)
   integer, intent(in)  :: av(E2DFIELD)       ! mask for v-transport (1: water)
   integer, intent(in)  :: hor_adv            ! selection for horizontal scheme
   integer, intent(in)  :: ver_adv            ! selection for vertical scheme
   integer, intent(in)  :: adv_split          ! selection for split mode
gotm's avatar
gotm committed
245 246
!
! !INPUT/OUTPUT PARAMETERS:
kbk's avatar
kbk committed
247
   REALTYPE, intent(inout)   :: f(I3DFIELD)
gotm's avatar
gotm committed
248 249 250 251
!
! !OUTPUT PARAMETERS:
!
! !LOCAL VARIABLES:
252
   REALTYPE, parameter       :: a1=0.5*ONE,a2=ONE
253
   integer         :: k
gotm's avatar
gotm committed
254 255 256 257 258 259 260 261 262 263 264 265
!
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'do_advection_3d() # ',Ncall
#endif

   select case (hor_adv)
      case (UPSTREAM)
kbk's avatar
kbk committed
266 267
         call upstream_adv(dt,f,uu,vv,ww,ho,hn, &
                           delxv,delyu,delxu,delyv,area_inv,az,AH)
268
      case ((UPSTREAM_SPLIT),(P2),(Superbee),(MUSCL),(P2_PDM),(FCT))
gotm's avatar
gotm committed
269
         hi=ho
270
         select case (adv_split)
gotm's avatar
gotm committed
271
            case (0)
kbk's avatar
kbk committed
272
               call u_split_adv(dt,f,uu,hun,delxu,delyu,area_inv,au,a2,&
gotm's avatar
gotm committed
273
                                hor_adv,az,AH)
274 275 276 277
               call update_3d_halo(f,f,az,& 
                                   iimin,jjmin,iimax,jjmax,kmax,D_TAG)
               call wait_halo(D_TAG)

278
#ifndef SLICE_MODEL
kbk's avatar
kbk committed
279
               call v_split_adv(dt,f,vv,hvn,delxv,delyv,area_inv,av,a2,&
gotm's avatar
gotm committed
280
                                hor_adv,az,AH)
kbk's avatar
kbk committed
281 282 283
               call update_3d_halo(f,f,az,& 
                                   iimin,jjmin,iimax,jjmax,kmax,D_TAG)
               call wait_halo(D_TAG)
284
#endif
kbk's avatar
kbk committed
285

gotm's avatar
gotm committed
286 287 288 289 290 291 292 293
               if (kmax.gt.1) then
#ifdef ITERATE_VERT_ADV
                  call w_split_it_adv(dt,f,ww,az,a2,ver_adv)
#else
                  call w_split_adv(dt,f,ww,az,a2,ver_adv)
#endif
               end if
            case (1)
kbk's avatar
kbk committed
294
               call u_split_adv(dt,f,uu,hun,delxu,delyu,area_inv,au,a1,&
gotm's avatar
gotm committed
295
                                hor_adv,az,AH)
296
               call update_3d_halo(f,f,az, &
297
                                   iimin,jjmin,iimax,jjmax,kmax,D_TAG)
298 299
               call wait_halo(D_TAG)

300
#ifndef SLICE_MODEL
kbk's avatar
kbk committed
301
               call v_split_adv(dt,f,vv,hvn,delxv,delyv,area_inv,av,a1,&
gotm's avatar
gotm committed
302
                                hor_adv,az,AH)
303 304 305
               call update_3d_halo(f,f,az, &
                                   iimin,jjmin,iimax,jjmax,kmax,D_TAG)
               call wait_halo(D_TAG)
306
#endif
307

gotm's avatar
gotm committed
308 309 310 311 312 313
               if (kmax.gt.1) then
#ifdef ITERATE_VERT_ADV
                  call w_split_it_adv(dt,f,ww,az,a2,ver_adv)
#else
                  call w_split_adv(dt,f,ww,az,a2,ver_adv)
#endif
kbk's avatar
kbk committed
314 315 316 317
                  call update_3d_halo(f,f,az, &
                                      iimin,jjmin,iimax,jjmax,kmax,D_TAG)
                  call wait_halo(D_TAG)

kbk's avatar
kbk committed
318
               end if
319
#ifndef SLICE_MODEL
kbk's avatar
kbk committed
320
               call v_split_adv(dt,f,vv,hvn,delxv,delyv,area_inv,av,a1,&
kbk's avatar
kbk committed
321
                                hor_adv,az,AH)
322 323 324
               call update_3d_halo(f,f,az, &
                                   iimin,jjmin,iimax,jjmax,kmax,D_TAG)
               call wait_halo(D_TAG)
325
#endif
326

kbk's avatar
kbk committed
327
               call u_split_adv(dt,f,uu,hun,delxu,delyu,area_inv,au,a1,&
kbk's avatar
kbk committed
328
                                hor_adv,az,AH)
kbk's avatar
kbk committed
329 330 331
               call update_3d_halo(f,f,az, &
                                   iimin,jjmin,iimax,jjmax,kmax,D_TAG)
               call wait_halo(D_TAG)
332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352

            case (2)
               select case (hor_adv)
                  case (UPSTREAM_SPLIT)
                     call upstream_2dh_adv(dt,f,uu,vv,ho,hn,hun,hvn, &
                               delxv,delyu,delxu,delyv,area_inv,az,AH)
                  case (FCT)
                     call fct_2dh_adv(dt,f,uu,vv,ho,hn,hun,hvn, &
                               delxv,delyu,delxu,delyv,area_inv,az,AH)
                  case default
                     FATAL 'For adv_split=2, hor_adv must be 2 (upstream) or 7 (fct)'
               end select
               if (kmax .gt. 1) then
#ifdef ITERATE_VERT_ADV
                  call w_split_it_adv(dt,f,ww,az,a2,ver_adv)
#else
                  call w_split_adv(dt,f,ww,az,a2,ver_adv)
#endif
               end if
            case default
               FATAL 'Not valid adv_split parameter'
gotm's avatar
gotm committed
353 354 355
         end select
      case default
         FATAL 'This is not so good - do_advection_3d()'
kbk's avatar
kbk committed
356
         stop
gotm's avatar
gotm committed
357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373
   end select

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

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

   end module advection_3d

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