Blame view

src/util/adv_center.F90 13.2 KB
33b83817   dumoda01   Depot du modele /...
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
412
413
414
!$Id: adv_center.F90,v 1.4 2006-11-06 13:36:46 hb Exp $
#include"cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: Advection schemes --- grid centers\label{sec:advectionMean}
!
! !INTERFACE:
   subroutine adv_center(N,dt,h,ho,ww,Bcup,Bcdw,Yup,Ydw,method,mode,Y)
!
! !DESCRIPTION:
!
! This subroutine solves a one-dimensional advection equation. There are two
! options, depending whether the advection should be conservative or not.
! Conservative advection has to be applied when settling of sediment or
! rising of phytoplankton is considered. In this case the advection is of
! the form
!  \begin{equation}
!   \label{Yadvection_cons}
!    \partder{Y}{t} = - \partder{F}{z}
!    \comma
!  \end{equation}
! where $F=wY$ is the flux caused by the advective velocity, $w$.
!
! Non-conservative advective transport has to be applied, when the water
! has a non-zero vertical velocity. In three-dimensional applications,
! this transport would be conservative, since vertical divergence would be
! compensated by horizontal convergence and vice versa. However, the
! key assumption of one-dimensional modelling is horizontal homogeneity, 
! such that we indeed have to apply a vertically non-conservative method,
! which is of the form
!  \begin{equation}
!   \label{Yadvection_noncons}
!    \partder{Y}{t} = - w\partder{Y}{z}
!                   = - \left(\partder{F}{z} - Y\partder{w}{z} \right).
!  \end{equation}
!
! The discretized form of \eq{Yadvection_cons} is
!  \begin{equation}
!   \label{advDiscretized_cons}
!   Y_i^{n+1} = Y_i^n
!   - \dfrac{\Delta t}{h_i}
!    \left( F^n_{i} - F^n_{i-1} \right)
!   \comma
!  \end{equation}
! where the integers $n$ and $i$ correspond to the present time and space
! level, respectively. 
!
! For the non-conservative form \eq{Yadvection_noncons}, 
! an extra term needs to be included:
!  \begin{equation}
!   \label{advDiscretized_noncons}
!   Y_i^{n+1} = Y_i^n
!   - \dfrac{\Delta t}{h_i}
!    \left( F^n_{i} - F^n_{i-1} -Y_i^n \left(w_k-w_{k-1}  \right)\right).
!  \end{equation}
!
!  Which advection method is applied is decided by the flag {\tt mode},
!  which gives conservative advection \eq{advDiscretized_cons}
!  for {\tt mode=1} and 
!  non-conservative advection \eq{advDiscretized_noncons} for {\tt mode=0}. 
!
! Fluxes are defined at the grid faces, the variable $Y_i$ is defined at the
!  grid centers. The fluxes are computed in an upstream-biased way,
!  \begin{equation}
!   \label{upstream}
!   F^n_{i} = \dfrac{1}{\Delta t}
!   \int_{z^\text{Face}_{i} - w \Delta t}^{z^\text{Face}_{i}} Y(z') dz'
!   \point
!  \end{equation}
!  For a third-order polynomial approximation of $Y$ (see \cite{Pietrzak98}),
!  these fluxes can be written the in so-called Lax-Wendroff form as
!  \begin{equation}
!   \label{fluxDiscretized}
!    \begin{array}{rcll}
!      F_{i} &=& w_{i} \left(Y_i +  \dfrac{1}{2} \Phi^+_{i}
!      \left(1-\magn{c_{i}} \right) \left( Y_{i+1} - Y_i \right) \right)
!      \quad & \text{for} \quad w_{i} > 0
!      \comma  \\[5mm]
!      F_{i} &=& w_{i} \left(Y_{i+1} +  \dfrac{1}{2} \Phi^-_{i}
!      \left(1-\magn{c_{i}} \right) \left( Y_i - Y_{i+1} \right) \right)
!      \quad & \text{for} \quad w_{i} < 0
!      \comma
!    \end{array}
!  \end{equation}
!  where $c_{i} = 2 w_{i} \Delta t / (h_i+h_{i+1})$ is the Courant number.
!  The factors appearing in \eq{fluxDiscretized} are defined as
!  \begin{equation}
!   \label{phiDiscretized}
!  \Phi^+_{i} =  \alpha_{i} +  \beta_{i}  r^+_{i}
!  \comma
!  \Phi^-_{i} =  \alpha_{i} +  \beta_{i}  r^-_{i}
!  \comma
!  \end{equation}
!  where
!  \begin{equation}
!   \label{alphaDiscretized}
!   \alpha_{i} = \dfrac{1}{2}
!    + \dfrac{1}{6} \left( 1- 2 \magn{c_{i}} \right) \comma
!   \beta_{i} = \dfrac{1}{2}
!    - \dfrac{1}{6} \left( 1- 2 \magn{c_{i}} \right)
!  \point
!  \end{equation}
!  The upstream and downstream slope parameters are
!  \begin{equation}
!   \label{slopeDiscretized}
!   r^+_{i} = \dfrac{Y_i - Y_{i-1}}{Y_{i+1}-Y_{i}}  \comma
!   r^-_{i} = \dfrac{Y_{i+2} - Y_{i+1}}{Y_{i+1}-Y_{i}}
!  \point
!  \end{equation}
!
!  To obtain monotonic and positive schemes also in the presence of strong
!  gradients, so-called slope limiters are aplied for the factors $\Phi^+_{i}$
!  and $\Phi^-_{i}$. The two most obvious cases are
!  the first-order upstream discretisation with $\Phi^+_{i}=\Phi^-_{i}=0$
!  and the Lax-Wendroff scheme with  $\Phi^+_{i}=\Phi^-_{i}=1$.
!  The subroutine {\tt adv\_center.F90} provides six different slope-limiters,
!  all discussed in detail by \cite{Pietrzak98}:
!
! \begin{itemize}
!  \item first-order upstream ({\tt method=UPSTREAM})
!  \item second-order upstream-biased polynomial scheme ({\tt method=P1},
!        not yet implemented)
!  \item third-order upstream-biased polynomial scheme ({\tt method=P2})
!  \item third-order scheme (TVD) with Superbee limiter ({\tt method=Superbee})
!  \item third-order scheme (TVD) with MUSCL limiter ({\tt method=MUSCL})
!  \item third-order scheme (TVD) with ULTIMATE QUICKEST limiter
!        ({\tt method=P2\_PDM})
! \end{itemize}
!
!
! If during a certain time step the maximum Courant number is larger
! than one, a split iteration will be carried out which guarantees that the
! split step Courant numbers are just smaller than 1.
!
! Several kinds of boundary conditions are implemented for the upper
! and lower boundaries. They are set by the integer values {\tt Bcup}
! and {\tt Bcdw}, that have to correspond to the parameters defined
! in the module {\tt util}, see \sect{sec:utils}. The
! following choices exist at the moment:
!
! For the value {\tt flux}, the boundary values {\tt Yup} and {\tt Ydw} are
! interpreted as specified fluxes at the uppermost and lowest interface.
! Fluxes into the boundary cells are counted positive by convention.
! For the value {\tt value}, {\tt Yup} and {\tt Ydw} specify the value
! of $Y$ at the interfaces, and the flux is computed by multiplying with
! the (known) speed  at the interface. For the value {\tt oneSided},
! {\tt Yup} and {\tt Ydw} are ignored and the flux is computed
! from a one-sided first-order upstream discretisation using the speed
! at the interface and the value of $Y$ at the center of the boundary cell.
! For the value {\tt zeroDivergence}, the fluxes into and out of the
! respective boundary cell are set equal.
! This corresponds to a zero-gradient formulation, or to zero
! flux divergence in the boundary cells.
!
! Be careful that your boundary conditions are mathematically well defined.
! For example, specifying an inflow into the boundary cell with the
! speed at the boundary being directed outward does not make sense.
!
!
! !USES:
   use util
   IMPLICIT NONE
!
! !INPUT PARAMETERS:

!  number of vertical layers
   integer,  intent(in)                :: N

!  time step (s)
   REALTYPE, intent(in)                :: dt

!  layer thickness (m)
   REALTYPE, intent(in)                :: h(0:N)

!  old layer thickness (m)
   REALTYPE, intent(in)                :: ho(0:N)

!  vertical advection speed
   REALTYPE, intent(in)                :: ww(0:N)

!  type of upper BC
   integer,  intent(in)                :: Bcup

!  type of lower BC
   integer,  intent(in)                :: Bcdw

!  value of upper BC
   REALTYPE, intent(in)                :: Yup

!  value of lower BC
   REALTYPE, intent(in)                :: Ydw

!  type of advection scheme
   integer,  intent(in)                :: method

!  advection mode (0: non-conservative, 1: conservative)
   integer,  intent(in)                :: mode
!
! !INPUT/OUTPUT PARAMETERS:
   REALTYPE                            :: Y(0:N)
!
! !DEFINED PARAMETERS:
   REALTYPE,     parameter             :: one6th=1.0d0/6.0d0
   integer,      parameter             :: itmax=100
!
! !REVISION HISTORY:
!  Original author(s): Lars Umlauf
!
!  $Log: adv_center.F90,v $
!  Revision 1.4  2006-11-06 13:36:46  hb
!  Option for conservative vertical advection added to adv_center
!
!  Revision 1.3  2006-03-20 09:06:38  kbk
!  removed explicit double precission dependency
!
!  Revision 1.2  2005/11/18 10:59:34  kbk
!  removed unused variables - some left in parameter lists
!
!  Revision 1.1  2005/06/27 10:54:33  kbk
!  new files needed
!
!
!
!EOP
!
! !LOCAL VARIABLES:
   integer                              :: i,k,it
   REALTYPE                             :: x,r,Phi,limit=_ZERO_
   REALTYPE                             :: Yu,Yc,Yd
   REALTYPE                             :: c,cmax
   REALTYPE                             :: cu(0:N)
!
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(*,*) 'adv_center # ',Ncall
#endif

!  initialize interface fluxes with zero
   cu   = _ZERO_

!  initialize maximum Courant number
   cmax = _ZERO_

!  compute maximum Courant number
   do k=1,N-1
      c=abs(ww(k))*dt/(0.5*(h(k)+h(k+1)))
      if (c.gt.cmax) cmax=c
   enddo

   it=min(itmax,int(cmax)+1)

!#ifdef DEBUG
   if (it .gt. 1) then
      STDERR 'In adv_center():'
      STDERR 'Maximum Courant number is ',cmax
      STDERR it,' iterations used for vertical advection'
   endif
!#endif

!  splitting loop
   do i=1,it

!     vertical loop
      do k=1,N-1

!        compute the slope ration
         if (ww(k) .gt. _ZERO_) then

!           compute Courant number
            c=ww(k)/float(it)*dt/(0.5*(h(k)+h(k+1)))

            if (k .gt. 1) then
               Yu=Y(k-1)                              ! upstream value
            else
               Yu=Y(k)
            end if
            Yc=Y(k  )                                 ! central value
            Yd=Y(k+1)                                 ! downstream value

!           compute slope ration
            if (abs(Yd-Yc) .gt. 1e-10) then
               r=(Yc-Yu)/(Yd-Yc)
            else
               r=(Yc-Yu)*1.e10
            end if

!        negative speed
         else

!           compute Courant number
            c=-ww(k)/float(it)*dt/(0.5*(h(k)+h(k+1)))

            if (k .lt. N-1) then
               Yu=Y(k+2)                              ! upstream value
            else
               Yu=Y(k+1)
            end if
            Yc=Y(k+1)                                 ! central value
            Yd=Y(k  )                                 ! downstream value


!           compute slope ratio
            if (abs(Yc-Yd) .gt. 1e-10) then
               r=(Yu-Yc)/(Yc-Yd)
            else
               r=(Yu-Yc)*1.e10
            end if

         end if

!        compute the flux-factor phi
         x    =  one6th*(1.-2.0*c)
         Phi  =  (0.5+x)+(0.5-x)*r

!        limit the flux according to different suggestions
         select case (method)
            case (UPSTREAM)
               limit=_ZERO_
            case (P1)
               FATAL "P1 advection method not yet implemented, choose other method"
               stop  "adv_center.F90"
            case ((P2),(P2_PDM))
               if (method.eq.P2) then
                  limit=Phi
               else
                  limit=max(_ZERO_,min(Phi,2./(1.-c),2.*r/(c+1.e-10)))
               end if
            case (Superbee)
               limit=max(_ZERO_, min(_ONE_, 2.0*r), min(r,2.*_ONE_) )
            case (MUSCL)
               limit=max(_ZERO_,min(2.*_ONE_,2.0*r,0.5*(1.0+r)))
            case default
               LEVEL3 method
               FATAL 'unkown advection method in adv_center()'
               stop
          end select

!        compute the limited flux
         cu(k)=ww(k)*(Yc+0.5*limit*(1-c)*(Yd-Yc))

      end do

!     do the upper boundary conditions
      select case (Bcup)
      case (flux)
         cu(N) = - Yup              ! flux into the domain is positive
      case (value)
         cu(N) =  ww(N)*Yup
      case (oneSided)
         if (ww(N).ge._ZERO_) then
            cu(N) =  ww(N)*Y(N)
         else
            cu(N) = _ZERO_
         end if
      case (zeroDivergence)
         cu(N) = cu(N-1)
      case default
         FATAL 'unkown upper boundary condition type in adv_center()'
         stop
      end select


!     do the lower boundary conditions
      select case (Bcdw)
      case (flux)
         cu(0) =   Ydw               ! flux into the domain is positive
      case (value)
         cu(0) =  ww(0)*Ydw
      case (oneSided)
         if(ww(0).le._ZERO_) then
            cu(0) =  ww(0)*Y(1)
         else
            cu(0) = _ZERO_
         end if
      case (zeroDivergence)
         cu(0) = cu(1)
      case default
         FATAL 'unkown lower boundary condition type in adv_center()'
         stop
      end select

!     do the vertical advection step which will be used for prescribed
!     vertical flow velocity and for settling of suspended matter.
 
      if (mode.eq.0) then ! non-conservative 
         do k=1,N
            Y(k)=Y(k)-1./float(it)*dt*((cu(k)-cu(k-1))/        &
                 h(k)-Y(k)*(ww(k)-ww(k-1))/h(k))
         enddo
      else                ! conservative  
         do k=1,N
            Y(k)=Y(k)-1./float(it)*dt*((cu(k)-cu(k-1))/h(k))
         enddo
      end if   

   end do ! end of the iteration loop


#ifdef DEBUG
   write(*,*) 'Leaving adv_center()'
   write(*,*)
#endif

   return
   end subroutine adv_center
!EOC

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