source: palm/trunk/SOURCE/advec_s_bc.f90 @ 1966

Last change on this file since 1966 was 1961, checked in by suehring, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 55.7 KB
Line 
1!> @file advec_s_bc.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: advec_s_bc.f90 1961 2016-07-12 16:37:58Z maronga $
26!
27! 1960 2016-07-12 16:34:24Z suehring
28! New CASE statement to treat scalars and humidity separately
29!
30! 1873 2016-04-18 14:50:06Z maronga
31! Module renamed (removed _mod)
32!
33!
34! 1850 2016-04-08 13:29:27Z maronga
35! Module renamed
36!
37!
38! 1815 2016-04-06 13:49:59Z raasch
39! comment change
40!
41! 1691 2015-10-26 16:17:44Z maronga
42! Formatting corrections
43!
44! 1682 2015-10-07 23:56:08Z knoop
45! Code annotations made doxygen readable
46!
47! 1517 2015-01-07 19:12:25Z hoffmann
48! interface added to advec_s_bc
49!
50! 1374 2014-04-25 12:55:07Z raasch
51! missing variables added to ONLY list
52!
53! 1361 2014-04-16 15:17:48Z hoffmann
54! nr and qr added
55!
56! 1353 2014-04-08 15:21:23Z heinze
57! REAL constants provided with KIND-attribute
58!
59! 1346 2014-03-27 13:18:20Z heinze
60! Bugfix: REAL constants provided with KIND-attribute especially in call of
61! intrinsic function like MAX, MIN, SIGN
62!
63! 1320 2014-03-20 08:40:49Z raasch
64! ONLY-attribute added to USE-statements,
65! kind-parameters added to all INTEGER and REAL declaration statements,
66! kinds are defined in new module kinds,
67! revision history before 2012 removed,
68! comment fields (!:) to be used for variable explanations added to
69! all variable declaration statements
70!
71! 1318 2014-03-17 13:35:16Z raasch
72! module interfaces removed
73!
74! 1092 2013-02-02 11:24:22Z raasch
75! unused variables removed
76!
77! 1036 2012-10-22 13:43:42Z raasch
78! code put under GPL (PALM 3.9)
79!
80! 1010 2012-09-20 07:59:54Z raasch
81! cpp switch __nopointer added for pointer free version
82!
83! Revision 1.1  1997/08/29 08:53:46  raasch
84! Initial revision
85!
86!
87! Description:
88! ------------
89!> Advection term for scalar quantities using the Bott-Chlond scheme.
90!> Computation in individual steps for each of the three dimensions.
91!> Limiting assumptions:
92!> So far the scheme has been assuming equidistant grid spacing. As this is not
93!> the case in the stretched portion of the z-direction, there dzw(k) is used as
94!> a substitute for a constant grid length. This certainly causes incorrect
95!> results; however, it is hoped that they are not too apparent for weakly
96!> stretched grids.
97!> NOTE: This is a provisional, non-optimised version!
98!------------------------------------------------------------------------------!
99MODULE advec_s_bc_mod
100 
101
102    PRIVATE
103    PUBLIC advec_s_bc
104
105    INTERFACE advec_s_bc
106       MODULE PROCEDURE advec_s_bc
107    END INTERFACE advec_s_bc
108
109 CONTAINS
110
111!------------------------------------------------------------------------------!
112! Description:
113! ------------
114!> @todo Missing subroutine description.
115!------------------------------------------------------------------------------!
116    SUBROUTINE advec_s_bc( sk, sk_char )
117
118       USE advection,                                                             &
119           ONLY:  aex, bex, dex, eex
120
121       USE arrays_3d,                                                             &
122           ONLY:  d, ddzw, dzu, dzw, tend, u, v, w
123
124       USE control_parameters,                                                    &
125           ONLY:  dt_3d, bc_pt_t_val, bc_q_t_val, bc_s_t_val, ibc_pt_b, ibc_pt_t, &
126                  ibc_q_t, ibc_s_t, message_string, pt_slope_offset,              &
127                  sloping_surface, u_gtrans, v_gtrans
128
129       USE cpulog,                                                                &
130           ONLY:  cpu_log, log_point_s
131
132       USE grid_variables,                                                        &
133           ONLY:  ddx, ddy
134
135       USE indices,                                                               &
136           ONLY:  nx, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
137
138       USE kinds
139
140       USE pegrid
141
142       USE statistics,                                                            &
143           ONLY:  rmask, statistic_regions, sums_wsts_bc_l
144
145
146       IMPLICIT NONE
147
148       CHARACTER (LEN=*) ::  sk_char !<
149
150       INTEGER(iwp) ::  i         !<
151       INTEGER(iwp) ::  ix        !<
152       INTEGER(iwp) ::  j         !<
153       INTEGER(iwp) ::  k         !<
154       INTEGER(iwp) ::  ngp       !<
155       INTEGER(iwp) ::  sr        !<
156       INTEGER(iwp) ::  type_xz_2 !<
157
158       REAL(wp) ::  cim    !<
159       REAL(wp) ::  cimf   !<
160       REAL(wp) ::  cip    !<
161       REAL(wp) ::  cipf   !<
162       REAL(wp) ::  d_new  !<
163       REAL(wp) ::  denomi !< denominator
164       REAL(wp) ::  fminus !<
165       REAL(wp) ::  fplus  !<
166       REAL(wp) ::  f2     !<
167       REAL(wp) ::  f4     !<
168       REAL(wp) ::  f8     !<
169       REAL(wp) ::  f12    !<
170       REAL(wp) ::  f24    !<
171       REAL(wp) ::  f48    !<
172       REAL(wp) ::  f1920  !<
173       REAL(wp) ::  im     !<
174       REAL(wp) ::  ip     !<
175       REAL(wp) ::  m2     !<
176       REAL(wp) ::  m3     !<
177       REAL(wp) ::  numera !< numerator
178       REAL(wp) ::  snenn  !<
179       REAL(wp) ::  sterm  !<
180       REAL(wp) ::  tendcy !<
181       REAL(wp) ::  t1     !<
182       REAL(wp) ::  t2     !<
183
184       REAL(wp) ::  fmax(2)   !<
185       REAL(wp) ::  fmax_l(2) !<
186       
187#if defined( __nopointer )
188       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  sk !<
189#else
190       REAL(wp), DIMENSION(:,:,:), POINTER ::  sk
191#endif
192
193       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a0   !<
194       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a1   !<
195       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a12  !<
196       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a2   !<
197       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  a22  !<
198       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  immb !<
199       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  imme !<
200       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  impb !<
201       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  impe !<
202       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ipmb !<
203       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ipme !<
204       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ippb !<
205       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ippe !<
206       
207       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  sk_p !<
208
209#if defined( __nec )
210       REAL(sp) ::  m1n, m1z  !< important for optimisation of division
211       REAL(sp), DIMENSION(:,:), ALLOCATABLE :: m1, sw !<
212#else
213       REAL(wp) ::  m1n, m1z 
214       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: m1, sw
215#endif
216
217
218!
219!--    Array sk_p requires 2 extra elements for each dimension
220       ALLOCATE( sk_p(nzb-2:nzt+3,nys-3:nyn+3,nxl-3:nxr+3) )
221       sk_p = 0.0_wp
222
223!
224!--    Assign reciprocal values in order to avoid divisions later
225       f2    = 0.5_wp
226       f4    = 0.25_wp
227       f8    = 0.125_wp
228       f12   = 0.8333333333333333E-01_wp
229       f24   = 0.4166666666666666E-01_wp
230       f48   = 0.2083333333333333E-01_wp
231       f1920 = 0.5208333333333333E-03_wp
232
233!
234!--    Advection in x-direction:
235
236!
237!--    Save the quantity to be advected in a local array
238!--    add an enlarged boundary in x-direction
239       DO  i = nxl-1, nxr+1
240          DO  j = nys, nyn
241             DO  k = nzb, nzt+1
242                sk_p(k,j,i) = sk(k,j,i)
243             ENDDO
244          ENDDO
245       ENDDO
246#if defined( __parallel )
247       ngp = 2 * ( nzt - nzb + 6 ) * ( nyn - nys + 7 )
248       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'start' )
249!
250!--    Send left boundary, receive right boundary
251       CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxl+1), ngp, MPI_REAL, pleft,  0,      &
252                          sk_p(nzb-2,nys-3,nxr+2), ngp, MPI_REAL, pright, 0,      &
253                          comm2d, status, ierr )
254!
255!--    Send right boundary, receive left boundary
256       CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxr-2), ngp, MPI_REAL, pright, 1,      &
257                          sk_p(nzb-2,nys-3,nxl-3), ngp, MPI_REAL, pleft,  1,      &
258                          comm2d, status, ierr )
259       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' )
260#else
261
262!
263!--    Cyclic boundary conditions
264       sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxr-2)
265       sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxr-1)
266       sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxl+1)
267       sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxl+2)
268#endif
269
270!
271!--    In case of a sloping surface, the additional gridpoints in x-direction
272!--    of the temperature field at the left and right boundary of the total
273!--    domain must be adjusted by the temperature difference between this distance
274       IF ( sloping_surface  .AND.  sk_char == 'pt' )  THEN
275          IF ( nxl ==  0 )  THEN
276             sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxl-3) - pt_slope_offset
277             sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxl-2) - pt_slope_offset
278          ENDIF
279          IF ( nxr == nx )  THEN
280             sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxr+2) + pt_slope_offset
281             sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxr+3) + pt_slope_offset
282          ENDIF
283       ENDIF
284
285!
286!--    Initialise control density
287       d = 0.0_wp
288
289!
290!--    Determine maxima of the first and second derivative in x-direction
291       fmax_l = 0.0_wp
292       DO  i = nxl, nxr
293          DO  j = nys, nyn
294             DO  k = nzb+1, nzt
295                numera = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) )
296                denomi  = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
297                fmax_l(1) = MAX( fmax_l(1) , numera )
298                fmax_l(2) = MAX( fmax_l(2) , denomi  )
299             ENDDO
300          ENDDO
301       ENDDO
302#if defined( __parallel )
303       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
304       CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
305#else
306       fmax = fmax_l
307#endif
308
309       fmax = 0.04_wp * fmax
310
311!
312!--    Allocate temporary arrays
313       ALLOCATE( a0(nzb+1:nzt,nxl-1:nxr+1),   a1(nzb+1:nzt,nxl-1:nxr+1),          &
314                 a2(nzb+1:nzt,nxl-1:nxr+1),   a12(nzb+1:nzt,nxl-1:nxr+1),         &
315                 a22(nzb+1:nzt,nxl-1:nxr+1),  immb(nzb+1:nzt,nxl-1:nxr+1),        &
316                 imme(nzb+1:nzt,nxl-1:nxr+1), impb(nzb+1:nzt,nxl-1:nxr+1),        &
317                 impe(nzb+1:nzt,nxl-1:nxr+1), ipmb(nzb+1:nzt,nxl-1:nxr+1),        &
318                 ipme(nzb+1:nzt,nxl-1:nxr+1), ippb(nzb+1:nzt,nxl-1:nxr+1),        &
319                 ippe(nzb+1:nzt,nxl-1:nxr+1), m1(nzb+1:nzt,nxl-2:nxr+2),          &
320                 sw(nzb+1:nzt,nxl-1:nxr+1)                                        &
321               )
322       imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
323
324!
325!--    Initialise point of time measuring of the exponential portion (this would
326!--    not work if done locally within the loop)
327       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'start' )
328       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
329
330!
331!--    Outer loop of all j
332       DO  j = nys, nyn
333
334!
335!--       Compute polynomial coefficients
336          DO  i = nxl-1, nxr+1
337             DO  k = nzb+1, nzt
338                a12(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
339                a22(k,i) = 0.5_wp * ( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i)        &
340                                                    + sk_p(k,j,i-1) )
341                a0(k,i) = ( 9.0_wp * sk_p(k,j,i+2)    - 116.0_wp * sk_p(k,j,i+1)  &
342                            + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j,i-1)  &
343                            + 9.0_wp * sk_p(k,j,i-2) ) * f1920
344                a1(k,i) = ( -5.0_wp * sk_p(k,j,i+2)   + 34.0_wp * sk_p(k,j,i+1)   &
345                            - 34.0_wp * sk_p(k,j,i-1) + 5.0_wp * sk_p(k,j,i-2)    &
346                          ) * f48
347                a2(k,i) = ( -3.0_wp * sk_p(k,j,i+2) + 36.0_wp * sk_p(k,j,i+1)     &
348                            - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j,i-1)     &
349                            - 3.0_wp * sk_p(k,j,i-2) ) * f48
350             ENDDO
351          ENDDO
352
353!
354!--       Fluxes using the Bott scheme
355!--       *VOCL LOOP,UNROLL(2)
356          DO  i = nxl, nxr
357             DO  k = nzb+1, nzt
358                cip  =  MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
359                cim  = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
360                cipf = 1.0_wp - 2.0_wp * cip
361                cimf = 1.0_wp - 2.0_wp * cim
362                ip   =   a0(k,i)   * f2  * ( 1.0_wp - cipf )                      &
363                       + a1(k,i)   * f8  * ( 1.0_wp - cipf*cipf )                 &
364                       + a2(k,i)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
365                im   =   a0(k,i+1) * f2  * ( 1.0_wp - cimf )                      &
366                       - a1(k,i+1) * f8  * ( 1.0_wp - cimf*cimf )                 &
367                       + a2(k,i+1) * f24 * ( 1.0_wp - cimf*cimf*cimf )
368                ip   = MAX( ip, 0.0_wp )
369                im   = MAX( im, 0.0_wp )
370                ippb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
371                impb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i+1) / (ip+im+1E-15_wp) )
372
373                cip  =  MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
374                cim  = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
375                cipf = 1.0_wp - 2.0_wp * cip
376                cimf = 1.0_wp - 2.0_wp * cim
377                ip   =   a0(k,i-1) * f2  * ( 1.0_wp - cipf )                      &
378                       + a1(k,i-1) * f8  * ( 1.0_wp - cipf*cipf )                 &
379                       + a2(k,i-1) * f24 * ( 1.0_wp - cipf*cipf*cipf )
380                im   =   a0(k,i)   * f2  * ( 1.0_wp - cimf )                      &
381                       - a1(k,i)   * f8  * ( 1.0_wp - cimf*cimf )                 &
382                       + a2(k,i)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
383                ip   = MAX( ip, 0.0_wp )
384                im   = MAX( im, 0.0_wp )
385                ipmb(k,i) = ip * MIN( 1.0_wp, sk_p(k,j,i-1) / (ip+im+1E-15_wp) )
386                immb(k,i) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
387             ENDDO
388          ENDDO
389
390!
391!--       Compute monitor function m1
392          DO  i = nxl-2, nxr+2
393             DO  k = nzb+1, nzt
394                m1z = ABS( sk_p(k,j,i+1) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j,i-1) )
395                m1n = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
396                IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
397                   m1(k,i) = m1z / m1n
398                   IF ( m1(k,i) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,i) = 0.0_wp
399                ELSEIF ( m1n < m1z )  THEN
400                   m1(k,i) = -1.0_wp
401                ELSE
402                   m1(k,i) = 0.0_wp
403                ENDIF
404             ENDDO
405          ENDDO
406
407!
408!--       Compute switch sw
409          sw = 0.0_wp
410          DO  i = nxl-1, nxr+1
411             DO  k = nzb+1, nzt
412                m2 = 2.0_wp * ABS( a1(k,i) - a12(k,i) ) /                         &
413                     MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35_wp )
414                IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) )  m2 = 0.0_wp
415
416                m3 = 2.0_wp * ABS( a2(k,i) - a22(k,i) ) /                         &
417                     MAX( ABS( a2(k,i) + a22(k,i) ), 1E-35_wp )
418                IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) )  m3 = 0.0_wp
419
420                t1 = 0.35_wp
421                t2 = 0.35_wp
422                IF ( m1(k,i) == -1.0_wp )  t2 = 0.12_wp
423
424!--             *VOCL STMT,IF(10)
425                IF ( m1(k,i-1) == 1.0_wp .OR. m1(k,i) == 1.0_wp                   &
426                     .OR. m1(k,i+1) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
427                     ( m1(k,i) > t1  .AND.  m1(k,i-1) /= -1.0_wp  .AND.           &
428                       m1(k,i) /= -1.0_wp  .AND.  m1(k,i+1) /= -1.0_wp )          &
429                   )  sw(k,i) = 1.0_wp
430             ENDDO
431          ENDDO
432
433!
434!--       Fluxes using the exponential scheme
435          CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
436          DO  i = nxl, nxr
437             DO  k = nzb+1, nzt
438
439!--             *VOCL STMT,IF(10)
440                IF ( sw(k,i) == 1.0_wp )  THEN
441                   snenn = sk_p(k,j,i+1) - sk_p(k,j,i-1)
442                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
443                   sterm = ( sk_p(k,j,i) - sk_p(k,j,i-1) ) / snenn
444                   sterm = MIN( sterm, 0.9999_wp )
445                   sterm = MAX( sterm, 0.0001_wp )
446
447                   ix = INT( sterm * 1000 ) + 1
448
449                   cip =  MAX( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
450
451                   ippe(k,i) = sk_p(k,j,i-1) * cip + snenn * (                    &
452                               aex(ix) * cip + bex(ix) / dex(ix) * (              &
453                               eex(ix) -                                          &
454                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
455                                                                   )              &
456                                                             )
457                   IF ( sterm == 0.0001_wp )  ippe(k,i) = sk_p(k,j,i) * cip
458                   IF ( sterm == 0.9999_wp )  ippe(k,i) = sk_p(k,j,i) * cip
459
460                   snenn = sk_p(k,j,i-1) - sk_p(k,j,i+1)
461                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
462                   sterm = ( sk_p(k,j,i) - sk_p(k,j,i+1) ) / snenn
463                   sterm = MIN( sterm, 0.9999_wp )
464                   sterm = MAX( sterm, 0.0001_wp )
465
466                   ix = INT( sterm * 1000 ) + 1
467
468                   cim = -MIN( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
469
470                   imme(k,i) = sk_p(k,j,i+1) * cim + snenn * (                    &
471                               aex(ix) * cim + bex(ix) / dex(ix) * (              &
472                               eex(ix) -                                          &
473                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
474                                                                   )              &
475                                                             )
476                   IF ( sterm == 0.0001_wp )  imme(k,i) = sk_p(k,j,i) * cim
477                   IF ( sterm == 0.9999_wp )  imme(k,i) = sk_p(k,j,i) * cim
478                ENDIF
479
480!--             *VOCL STMT,IF(10)
481                IF ( sw(k,i+1) == 1.0_wp )  THEN
482                   snenn = sk_p(k,j,i) - sk_p(k,j,i+2)
483                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
484                   sterm = ( sk_p(k,j,i+1) - sk_p(k,j,i+2) ) / snenn
485                   sterm = MIN( sterm, 0.9999_wp )
486                   sterm = MAX( sterm, 0.0001_wp )
487
488                   ix = INT( sterm * 1000 ) + 1
489
490                   cim = -MIN( 0.0_wp, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
491
492                   impe(k,i) = sk_p(k,j,i+2) * cim + snenn * (                    &
493                               aex(ix) * cim + bex(ix) / dex(ix) * (              &
494                               eex(ix) -                                          &
495                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
496                                                                   )              &
497                                                             )
498                   IF ( sterm == 0.0001_wp )  impe(k,i) = sk_p(k,j,i+1) * cim
499                   IF ( sterm == 0.9999_wp )  impe(k,i) = sk_p(k,j,i+1) * cim
500                ENDIF
501
502!--             *VOCL STMT,IF(10)
503                IF ( sw(k,i-1) == 1.0_wp )  THEN
504                   snenn = sk_p(k,j,i) - sk_p(k,j,i-2)
505                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
506                   sterm = ( sk_p(k,j,i-1) - sk_p(k,j,i-2) ) / snenn
507                   sterm = MIN( sterm, 0.9999_wp )
508                   sterm = MAX( sterm, 0.0001_wp )
509
510                   ix = INT( sterm * 1000 ) + 1
511
512                   cip = MAX( 0.0_wp, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
513
514                   ipme(k,i) = sk_p(k,j,i-2) * cip + snenn * (                    &
515                               aex(ix) * cip + bex(ix) / dex(ix) * (              &
516                               eex(ix) -                                          &
517                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
518                                                                   )              &
519                                                             )
520                   IF ( sterm == 0.0001_wp )  ipme(k,i) = sk_p(k,j,i-1) * cip
521                   IF ( sterm == 0.9999_wp )  ipme(k,i) = sk_p(k,j,i-1) * cip
522                ENDIF
523
524             ENDDO
525          ENDDO
526          CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
527
528!
529!--       Prognostic equation
530          DO  i = nxl, nxr
531             DO  k = nzb+1, nzt
532                fplus  = ( 1.0_wp - sw(k,i)   ) * ippb(k,i) + sw(k,i)   * ippe(k,i)  &
533                       - ( 1.0_wp - sw(k,i+1) ) * impb(k,i) - sw(k,i+1) * impe(k,i)
534                fminus = ( 1.0_wp - sw(k,i-1) ) * ipmb(k,i) + sw(k,i-1) * ipme(k,i)  &
535                       - ( 1.0_wp - sw(k,i)   ) * immb(k,i) - sw(k,i)   * imme(k,i)
536                tendcy = fplus - fminus
537!
538!--             Removed in order to optimize speed
539!                ffmax   = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp )
540!                IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp )  tendcy = 0.0
541!
542!--             Density correction because of possible remaining divergences
543                d_new = d(k,j,i) - ( u(k,j,i+1) - u(k,j,i) ) * dt_3d * ddx
544                sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /    &
545                              ( 1.0_wp + d_new )
546                d(k,j,i)  = d_new
547             ENDDO
548          ENDDO
549
550       ENDDO   ! End of the advection in x-direction
551
552!
553!--    Deallocate temporary arrays
554       DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme,      &
555                   ippb, ippe, m1, sw )
556
557
558!
559!--    Enlarge boundary of local array cyclically in y-direction
560#if defined( __parallel )
561       ngp = ( nzt - nzb + 6 ) * ( nyn - nys + 7 )
562       CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+6), ngp, MPI_REAL,             &
563                             type_xz_2, ierr )
564       CALL MPI_TYPE_COMMIT( type_xz_2, ierr )
565!
566!--    Send front boundary, receive rear boundary
567       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' )
568       CALL MPI_SENDRECV( sk_p(nzb-2,nys,nxl-3),   1, type_xz_2, psouth, 0,       &
569                          sk_p(nzb-2,nyn+1,nxl-3), 1, type_xz_2, pnorth, 0,       &
570                          comm2d, status, ierr )
571!
572!--    Send rear boundary, receive front boundary
573       CALL MPI_SENDRECV( sk_p(nzb-2,nyn-2,nxl-3), 1, type_xz_2, pnorth, 1,       &
574                          sk_p(nzb-2,nys-3,nxl-3), 1, type_xz_2, psouth, 1,       &
575                          comm2d, status, ierr )
576       CALL MPI_TYPE_FREE( type_xz_2, ierr )
577       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' )
578#else
579       DO  i = nxl, nxr
580          DO  k = nzb+1, nzt
581             sk_p(k,nys-1,i) = sk_p(k,nyn,i)
582             sk_p(k,nys-2,i) = sk_p(k,nyn-1,i)
583             sk_p(k,nys-3,i) = sk_p(k,nyn-2,i)
584             sk_p(k,nyn+1,i) = sk_p(k,nys,i)
585             sk_p(k,nyn+2,i) = sk_p(k,nys+1,i)
586             sk_p(k,nyn+3,i) = sk_p(k,nys+2,i)
587          ENDDO
588       ENDDO
589#endif
590
591!
592!--    Determine the maxima of the first and second derivative in y-direction
593       fmax_l = 0.0_wp
594       DO  i = nxl, nxr
595          DO  j = nys, nyn
596             DO  k = nzb+1, nzt
597                numera = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) )
598                denomi  = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
599                fmax_l(1) = MAX( fmax_l(1) , numera )
600                fmax_l(2) = MAX( fmax_l(2) , denomi  )
601             ENDDO
602          ENDDO
603       ENDDO
604#if defined( __parallel )
605       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
606       CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
607#else
608       fmax = fmax_l
609#endif
610
611       fmax = 0.04_wp * fmax
612
613!
614!--    Allocate temporary arrays
615       ALLOCATE( a0(nzb+1:nzt,nys-1:nyn+1),   a1(nzb+1:nzt,nys-1:nyn+1),          &
616                 a2(nzb+1:nzt,nys-1:nyn+1),   a12(nzb+1:nzt,nys-1:nyn+1),         &
617                 a22(nzb+1:nzt,nys-1:nyn+1),  immb(nzb+1:nzt,nys-1:nyn+1),        &
618                 imme(nzb+1:nzt,nys-1:nyn+1), impb(nzb+1:nzt,nys-1:nyn+1),        &
619                 impe(nzb+1:nzt,nys-1:nyn+1), ipmb(nzb+1:nzt,nys-1:nyn+1),        &
620                 ipme(nzb+1:nzt,nys-1:nyn+1), ippb(nzb+1:nzt,nys-1:nyn+1),        &
621                 ippe(nzb+1:nzt,nys-1:nyn+1), m1(nzb+1:nzt,nys-2:nyn+2),          &
622                 sw(nzb+1:nzt,nys-1:nyn+1)                                        &
623               )
624       imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
625
626!
627!--    Outer loop of all i
628       DO  i = nxl, nxr
629
630!
631!--       Compute polynomial coefficients
632          DO  j = nys-1, nyn+1
633             DO  k = nzb+1, nzt
634                a12(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
635                a22(k,j) = 0.5_wp * ( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i)        &
636                                                    + sk_p(k,j-1,i) )
637                a0(k,j) = ( 9.0_wp * sk_p(k,j+2,i)    - 116.0_wp * sk_p(k,j+1,i)  &
638                            + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k,j-1,i)  &
639                            + 9.0_wp * sk_p(k,j-2,i) ) * f1920
640                a1(k,j) = ( -5.0_wp   * sk_p(k,j+2,i) + 34.0_wp * sk_p(k,j+1,i)   &
641                            - 34.0_wp * sk_p(k,j-1,i) + 5.0_wp  * sk_p(k,j-2,i)   &
642                          ) * f48
643                a2(k,j) = ( -3.0_wp * sk_p(k,j+2,i) + 36.0_wp * sk_p(k,j+1,i)     &
644                            - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k,j-1,i)     &
645                            - 3.0_wp * sk_p(k,j-2,i) ) * f48
646             ENDDO
647          ENDDO
648
649!
650!--       Fluxes using the Bott scheme
651!--       *VOCL LOOP,UNROLL(2)
652          DO  j = nys, nyn
653             DO  k = nzb+1, nzt
654                cip  =  MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
655                cim  = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
656                cipf = 1.0_wp - 2.0_wp * cip
657                cimf = 1.0_wp - 2.0_wp * cim
658                ip   =   a0(k,j)   * f2  * ( 1.0_wp - cipf )                      &
659                       + a1(k,j)   * f8  * ( 1.0_wp - cipf*cipf )                 &
660                       + a2(k,j)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
661                im   =   a0(k,j+1) * f2  * ( 1.0_wp - cimf )                      &
662                       - a1(k,j+1) * f8  * ( 1.0_wp - cimf*cimf )                 &
663                       + a2(k,j+1) * f24 * ( 1.0_wp - cimf*cimf*cimf )
664                ip   = MAX( ip, 0.0_wp )
665                im   = MAX( im, 0.0_wp )
666                ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
667                impb(k,j) = im * MIN( 1.0_wp, sk_p(k,j+1,i) / (ip+im+1E-15_wp) )
668
669                cip  =  MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
670                cim  = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
671                cipf = 1.0_wp - 2.0_wp * cip
672                cimf = 1.0_wp - 2.0_wp * cim
673                ip   =   a0(k,j-1) * f2  * ( 1.0_wp - cipf )                      &
674                       + a1(k,j-1) * f8  * ( 1.0_wp - cipf*cipf )                 &
675                       + a2(k,j-1) * f24 * ( 1.0_wp - cipf*cipf*cipf )
676                im   =   a0(k,j)   * f2  * ( 1.0_wp - cimf )                      &
677                       - a1(k,j)   * f8  * ( 1.0_wp - cimf*cimf )                 &
678                       + a2(k,j)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
679                ip   = MAX( ip, 0.0_wp )
680                im   = MAX( im, 0.0_wp )
681                ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j-1,i) / (ip+im+1E-15_wp) )
682                immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
683             ENDDO
684          ENDDO
685
686!
687!--       Compute monitor function m1
688          DO  j = nys-2, nyn+2
689             DO  k = nzb+1, nzt
690                m1z = ABS( sk_p(k,j+1,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k,j-1,i) )
691                m1n = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
692                IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
693                   m1(k,j) = m1z / m1n
694                   IF ( m1(k,j) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0_wp
695                ELSEIF ( m1n < m1z )  THEN
696                   m1(k,j) = -1.0_wp
697                ELSE
698                   m1(k,j) = 0.0_wp
699                ENDIF
700             ENDDO
701          ENDDO
702
703!
704!--       Compute switch sw
705          sw = 0.0_wp
706          DO  j = nys-1, nyn+1
707             DO  k = nzb+1, nzt
708                m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) /                            &
709                     MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp )
710                IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0_wp
711
712                m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) /                            &
713                     MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp )
714                IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0_wp
715
716                t1 = 0.35_wp
717                t2 = 0.35_wp
718                IF ( m1(k,j) == -1.0_wp )  t2 = 0.12_wp
719
720!--             *VOCL STMT,IF(10)
721                IF ( m1(k,j-1) == 1.0_wp .OR. m1(k,j) == 1.0_wp                   &
722                     .OR. m1(k,j+1) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
723                     ( m1(k,j) > t1  .AND.  m1(k,j-1) /= -1.0_wp  .AND.           &
724                       m1(k,j) /= -1.0_wp  .AND.  m1(k,j+1) /= -1.0_wp )          &
725                   )  sw(k,j) = 1.0_wp
726             ENDDO
727          ENDDO
728
729!
730!--       Fluxes using exponential scheme
731          CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
732          DO  j = nys, nyn
733             DO  k = nzb+1, nzt
734
735!--             *VOCL STMT,IF(10)
736                IF ( sw(k,j) == 1.0_wp )  THEN
737                   snenn = sk_p(k,j+1,i) - sk_p(k,j-1,i)
738                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
739                   sterm = ( sk_p(k,j,i) - sk_p(k,j-1,i) ) / snenn
740                   sterm = MIN( sterm, 0.9999_wp )
741                   sterm = MAX( sterm, 0.0001_wp )
742
743                   ix = INT( sterm * 1000 ) + 1
744
745                   cip =  MAX( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
746
747                   ippe(k,j) = sk_p(k,j-1,i) * cip + snenn * (                    &
748                               aex(ix) * cip + bex(ix) / dex(ix) * (              &
749                               eex(ix) -                                          &
750                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
751                                                                   )              &
752                                                             )
753                   IF ( sterm == 0.0001_wp )  ippe(k,j) = sk_p(k,j,i) * cip
754                   IF ( sterm == 0.9999_wp )  ippe(k,j) = sk_p(k,j,i) * cip
755
756                   snenn = sk_p(k,j-1,i) - sk_p(k,j+1,i)
757                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
758                   sterm = ( sk_p(k,j,i) - sk_p(k,j+1,i) ) / snenn
759                   sterm = MIN( sterm, 0.9999_wp )
760                   sterm = MAX( sterm, 0.0001_wp )
761
762                   ix = INT( sterm * 1000 ) + 1
763
764                   cim = -MIN( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
765
766                   imme(k,j) = sk_p(k,j+1,i) * cim + snenn * (                    &
767                               aex(ix) * cim + bex(ix) / dex(ix) * (              &
768                               eex(ix) -                                          &
769                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
770                                                                   )              &
771                                                             )
772                   IF ( sterm == 0.0001_wp )  imme(k,j) = sk_p(k,j,i) * cim
773                   IF ( sterm == 0.9999_wp )  imme(k,j) = sk_p(k,j,i) * cim
774                ENDIF
775
776!--             *VOCL STMT,IF(10)
777                IF ( sw(k,j+1) == 1.0_wp )  THEN
778                   snenn = sk_p(k,j,i) - sk_p(k,j+2,i)
779                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
780                   sterm = ( sk_p(k,j+1,i) - sk_p(k,j+2,i) ) / snenn
781                   sterm = MIN( sterm, 0.9999_wp )
782                   sterm = MAX( sterm, 0.0001_wp )
783
784                   ix = INT( sterm * 1000 ) + 1
785
786                   cim = -MIN( 0.0_wp, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
787
788                   impe(k,j) = sk_p(k,j+2,i) * cim + snenn * (                    &
789                               aex(ix) * cim + bex(ix) / dex(ix) * (              &
790                               eex(ix) -                                          &
791                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
792                                                                   )              &
793                                                             )
794                   IF ( sterm == 0.0001_wp )  impe(k,j) = sk_p(k,j+1,i) * cim
795                   IF ( sterm == 0.9999_wp )  impe(k,j) = sk_p(k,j+1,i) * cim
796                ENDIF
797
798!--             *VOCL STMT,IF(10)
799                IF ( sw(k,j-1) == 1.0_wp )  THEN
800                   snenn = sk_p(k,j,i) - sk_p(k,j-2,i)
801                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
802                   sterm = ( sk_p(k,j-1,i) - sk_p(k,j-2,i) ) / snenn
803                   sterm = MIN( sterm, 0.9999_wp )
804                   sterm = MAX( sterm, 0.0001_wp )
805
806                   ix = INT( sterm * 1000 ) + 1
807
808                   cip = MAX( 0.0_wp, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
809
810                   ipme(k,j) = sk_p(k,j-2,i) * cip + snenn * (                    &
811                               aex(ix) * cip + bex(ix) / dex(ix) * (              &
812                               eex(ix) -                                          &
813                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
814                                                                   )              &
815                                                             )
816                   IF ( sterm == 0.0001_wp )  ipme(k,j) = sk_p(k,j-1,i) * cip
817                   IF ( sterm == 0.9999_wp )  ipme(k,j) = sk_p(k,j-1,i) * cip
818                ENDIF
819
820             ENDDO
821          ENDDO
822          CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
823
824!
825!--       Prognostic equation
826          DO  j = nys, nyn
827             DO  k = nzb+1, nzt
828                fplus  = ( 1.0_wp - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
829                       - ( 1.0_wp - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j)
830                fminus = ( 1.0_wp - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) &
831                       - ( 1.0_wp - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
832                tendcy = fplus - fminus
833!
834!--             Removed in order to optimise speed
835!                ffmax   = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp )
836!                 IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp )  tendcy = 0.0
837!
838!--             Density correction because of possible remaining divergences
839                d_new = d(k,j,i) - ( v(k,j+1,i) - v(k,j,i) ) * dt_3d * ddy
840                sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /  &
841                              ( 1.0_wp + d_new )
842                d(k,j,i)  = d_new
843             ENDDO
844          ENDDO
845
846       ENDDO   ! End of the advection in y-direction
847       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' )
848       CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'stop' )
849
850!
851!--    Deallocate temporary arrays
852       DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme,      &
853                   ippb, ippe, m1, sw )
854
855
856!
857!--    Initialise for the computation of heat fluxes (see below; required in
858!--    UP flow_statistics)
859       IF ( sk_char == 'pt' )  sums_wsts_bc_l = 0.0_wp
860
861!
862!--    Add top and bottom boundaries according to the relevant boundary conditions
863       IF ( sk_char == 'pt' )  THEN
864
865!
866!--       Temperature boundary condition at the bottom boundary
867          IF ( ibc_pt_b == 0 )  THEN
868!
869!--       Dirichlet (fixed surface temperature)
870             DO  i = nxl, nxr
871                DO  j = nys, nyn
872                   sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
873                   sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
874                ENDDO
875             ENDDO
876
877          ELSE
878!
879!--          Neumann (i.e. here zero gradient)
880             DO  i = nxl, nxr
881                DO  j = nys, nyn
882                   sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
883                   sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
884                   sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
885                ENDDO
886             ENDDO
887
888          ENDIF
889
890!
891!--       Temperature boundary condition at the top boundary
892          IF ( ibc_pt_t == 0  .OR.  ibc_pt_t == 1 )  THEN
893!
894!--          Dirichlet or Neumann (zero gradient)
895             DO  i = nxl, nxr
896                DO  j = nys, nyn
897                   sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
898                   sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
899                ENDDO
900             ENDDO
901
902          ELSEIF ( ibc_pt_t == 2 )  THEN
903!
904!--          Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
905             DO  i = nxl, nxr
906                DO  j = nys, nyn
907                   sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_pt_t_val * dzu(nzt+1)
908                   sk_p(nzt+3,j,i)   = sk_p(nzt+2,j,i) + bc_pt_t_val * dzu(nzt+1)
909                ENDDO
910             ENDDO
911
912          ENDIF
913
914       ELSEIF ( sk_char == 'sa' )  THEN
915
916!
917!--       Salinity boundary condition at the bottom boundary.
918!--       So far, always Neumann (i.e. here zero gradient) is used
919          DO  i = nxl, nxr
920             DO  j = nys, nyn
921                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
922                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
923                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
924             ENDDO
925          ENDDO
926
927!
928!--       Salinity boundary condition at the top boundary.
929!--       Dirichlet or Neumann (zero gradient)
930          DO  i = nxl, nxr
931             DO  j = nys, nyn
932                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
933                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
934             ENDDO
935          ENDDO
936
937       ELSEIF ( sk_char == 'q' )  THEN
938
939!
940!--       Specific humidity boundary condition at the bottom boundary.
941!--       Dirichlet (fixed surface humidity) or Neumann (i.e. zero gradient)
942          DO  i = nxl, nxr
943             DO  j = nys, nyn
944                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
945                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
946                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
947             ENDDO
948          ENDDO
949
950!
951!--       Specific humidity boundary condition at the top boundary
952          IF ( ibc_q_t == 0 )  THEN
953!
954!--          Dirichlet
955             DO  i = nxl, nxr
956                DO  j = nys, nyn
957                   sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
958                   sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
959                ENDDO
960             ENDDO
961
962          ELSE
963!
964!--          Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
965             DO  i = nxl, nxr
966                DO  j = nys, nyn
967                   sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_q_t_val * dzu(nzt+1)
968                   sk_p(nzt+3,j,i)   = sk_p(nzt+2,j,i) + bc_q_t_val * dzu(nzt+1)
969                ENDDO
970             ENDDO
971
972          ENDIF
973
974       ELSEIF ( sk_char == 's' )  THEN
975!
976!--       Specific scalar boundary condition at the bottom boundary.
977!--       Dirichlet (fixed surface humidity) or Neumann (i.e. zero gradient)
978          DO  i = nxl, nxr
979             DO  j = nys, nyn
980                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
981                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
982                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
983             ENDDO
984          ENDDO
985
986!
987!--       Specific scalar boundary condition at the top boundary
988          IF ( ibc_s_t == 0 )  THEN
989!
990!--          Dirichlet
991             DO  i = nxl, nxr
992                DO  j = nys, nyn
993                   sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
994                   sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
995                ENDDO
996             ENDDO
997
998          ELSE
999!
1000!--          Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
1001             DO  i = nxl, nxr
1002                DO  j = nys, nyn
1003                   sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_s_t_val * dzu(nzt+1)
1004                   sk_p(nzt+3,j,i)   = sk_p(nzt+2,j,i) + bc_s_t_val * dzu(nzt+1)
1005                ENDDO
1006             ENDDO
1007
1008          ENDIF
1009
1010       ELSEIF ( sk_char == 'qr' )  THEN
1011
1012!
1013!--       Rain water content boundary condition at the bottom boundary:
1014!--       Dirichlet (fixed surface rain water content).
1015          DO  i = nxl, nxr
1016             DO  j = nys, nyn
1017                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
1018                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
1019                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
1020             ENDDO
1021          ENDDO
1022
1023!
1024!--       Rain water content boundary condition at the top boundary: Dirichlet
1025          DO  i = nxl, nxr
1026             DO  j = nys, nyn
1027                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
1028                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
1029             ENDDO
1030          ENDDO
1031
1032       ELSEIF ( sk_char == 'nr' )  THEN
1033
1034!
1035!--       Rain drop concentration boundary condition at the bottom boundary:
1036!--       Dirichlet (fixed surface rain drop concentration).
1037          DO  i = nxl, nxr
1038             DO  j = nys, nyn
1039                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
1040                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
1041                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
1042             ENDDO
1043          ENDDO
1044
1045!
1046!--       Rain drop concentration boundary condition at the top boundary: Dirichlet
1047          DO  i = nxl, nxr
1048             DO  j = nys, nyn
1049                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
1050                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
1051             ENDDO
1052          ENDDO
1053
1054       ELSEIF ( sk_char == 'e' )  THEN
1055
1056!
1057!--       TKE boundary condition at bottom and top boundary (generally Neumann)
1058          DO  i = nxl, nxr
1059             DO  j = nys, nyn
1060                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
1061                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
1062                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
1063                sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
1064                sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
1065             ENDDO
1066          ENDDO
1067
1068       ELSE
1069
1070          WRITE( message_string, * ) 'no vertical boundary condi',                &
1071                                   'tion for variable "', sk_char, '"'
1072          CALL message( 'advec_s_bc', 'PA0158', 1, 2, 0, 6, 0 )
1073         
1074       ENDIF
1075
1076!
1077!--    Determine the maxima of the first and second derivative in z-direction
1078       fmax_l = 0.0_wp
1079       DO  i = nxl, nxr
1080          DO  j = nys, nyn
1081             DO  k = nzb, nzt+1
1082                numera = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) )
1083                denomi  = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) )
1084                fmax_l(1) = MAX( fmax_l(1) , numera )
1085                fmax_l(2) = MAX( fmax_l(2) , denomi  )
1086             ENDDO
1087          ENDDO
1088       ENDDO
1089#if defined( __parallel )
1090       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1091       CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
1092#else
1093       fmax = fmax_l
1094#endif
1095
1096       fmax = 0.04_wp * fmax
1097
1098!
1099!--    Allocate temporary arrays
1100       ALLOCATE( a0(nzb:nzt+1,nys:nyn),   a1(nzb:nzt+1,nys:nyn),                  &
1101                 a2(nzb:nzt+1,nys:nyn),   a12(nzb:nzt+1,nys:nyn),                 &
1102                 a22(nzb:nzt+1,nys:nyn),  immb(nzb+1:nzt,nys:nyn),                &
1103                 imme(nzb+1:nzt,nys:nyn), impb(nzb+1:nzt,nys:nyn),                &
1104                 impe(nzb+1:nzt,nys:nyn), ipmb(nzb+1:nzt,nys:nyn),                &
1105                 ipme(nzb+1:nzt,nys:nyn), ippb(nzb+1:nzt,nys:nyn),                &
1106                 ippe(nzb+1:nzt,nys:nyn), m1(nzb-1:nzt+2,nys:nyn),                &
1107                 sw(nzb:nzt+1,nys:nyn)                                            &
1108               )
1109       imme = 0.0_wp; impe = 0.0_wp; ipme = 0.0_wp; ippe = 0.0_wp
1110
1111!
1112!--    Outer loop of all i
1113       DO  i = nxl, nxr
1114
1115!
1116!--       Compute polynomial coefficients
1117          DO  j = nys, nyn
1118             DO  k = nzb, nzt+1
1119                a12(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
1120                a22(k,j) = 0.5_wp * ( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i)        &
1121                                                    + sk_p(k-1,j,i) )
1122                a0(k,j) = ( 9.0_wp * sk_p(k+2,j,i)    - 116.0_wp * sk_p(k+1,j,i)  &
1123                            + 2134.0_wp * sk_p(k,j,i) - 116.0_wp * sk_p(k-1,j,i)  &
1124                            + 9.0_wp * sk_p(k-2,j,i) ) * f1920
1125                a1(k,j) = ( -5.0_wp   * sk_p(k+2,j,i) + 34.0_wp * sk_p(k+1,j,i)   &
1126                            - 34.0_wp * sk_p(k-1,j,i) + 5.0_wp  * sk_p(k-2,j,i)   &
1127                          ) * f48
1128                a2(k,j) = ( -3.0_wp * sk_p(k+2,j,i) + 36.0_wp * sk_p(k+1,j,i)     &
1129                            - 66.0_wp * sk_p(k,j,i) + 36.0_wp * sk_p(k-1,j,i)     &
1130                            - 3.0_wp * sk_p(k-2,j,i) ) * f48
1131             ENDDO
1132          ENDDO
1133
1134!
1135!--       Fluxes using the Bott scheme
1136!--       *VOCL LOOP,UNROLL(2)
1137          DO  j = nys, nyn
1138             DO  k = nzb+1, nzt
1139                cip  =  MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
1140                cim  = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
1141                cipf = 1.0_wp - 2.0_wp * cip
1142                cimf = 1.0_wp - 2.0_wp * cim
1143                ip   =   a0(k,j)   * f2  * ( 1.0_wp - cipf )                      &
1144                       + a1(k,j)   * f8  * ( 1.0_wp - cipf*cipf )                 &
1145                       + a2(k,j)   * f24 * ( 1.0_wp - cipf*cipf*cipf )
1146                im   =   a0(k+1,j) * f2  * ( 1.0_wp - cimf )                      &
1147                       - a1(k+1,j) * f8  * ( 1.0_wp - cimf*cimf )                 &
1148                       + a2(k+1,j) * f24 * ( 1.0_wp - cimf*cimf*cimf )
1149                ip   = MAX( ip, 0.0_wp )
1150                im   = MAX( im, 0.0_wp )
1151                ippb(k,j) = ip * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
1152                impb(k,j) = im * MIN( 1.0_wp, sk_p(k+1,j,i) / (ip+im+1E-15_wp) )
1153
1154                cip  =  MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
1155                cim  = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
1156                cipf = 1.0_wp - 2.0_wp * cip
1157                cimf = 1.0_wp - 2.0_wp * cim
1158                ip   =   a0(k-1,j) * f2  * ( 1.0_wp - cipf )                      &
1159                       + a1(k-1,j) * f8  * ( 1.0_wp - cipf*cipf )                 &
1160                       + a2(k-1,j) * f24 * ( 1.0_wp - cipf*cipf*cipf )
1161                im   =   a0(k,j)   * f2  * ( 1.0_wp - cimf )                      &
1162                       - a1(k,j)   * f8  * ( 1.0_wp - cimf*cimf )                 &
1163                       + a2(k,j)   * f24 * ( 1.0_wp - cimf*cimf*cimf )
1164                ip   = MAX( ip, 0.0_wp )
1165                im   = MAX( im, 0.0_wp )
1166                ipmb(k,j) = ip * MIN( 1.0_wp, sk_p(k-1,j,i) / (ip+im+1E-15_wp) )
1167                immb(k,j) = im * MIN( 1.0_wp, sk_p(k,j,i)   / (ip+im+1E-15_wp) )
1168             ENDDO
1169          ENDDO
1170
1171!
1172!--       Compute monitor function m1
1173          DO  j = nys, nyn
1174             DO  k = nzb-1, nzt+2
1175                m1z = ABS( sk_p(k+1,j,i) - 2.0_wp * sk_p(k,j,i) + sk_p(k-1,j,i) )
1176                m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
1177                IF ( m1n /= 0.0_wp  .AND.  m1n >= m1z )  THEN
1178                   m1(k,j) = m1z / m1n
1179                   IF ( m1(k,j) /= 2.0_wp  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0_wp
1180                ELSEIF ( m1n < m1z )  THEN
1181                   m1(k,j) = -1.0_wp
1182                ELSE
1183                   m1(k,j) = 0.0_wp
1184                ENDIF
1185             ENDDO
1186          ENDDO
1187
1188!
1189!--       Compute switch sw
1190          sw = 0.0_wp
1191          DO  j = nys, nyn
1192             DO  k = nzb, nzt+1
1193                m2 = 2.0_wp * ABS( a1(k,j) - a12(k,j) ) /                         &
1194                     MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35_wp )
1195                IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0_wp
1196
1197                m3 = 2.0_wp * ABS( a2(k,j) - a22(k,j) ) /                         &
1198                     MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35_wp )
1199                IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0_wp
1200
1201                t1 = 0.35_wp
1202                t2 = 0.35_wp
1203                IF ( m1(k,j) == -1.0_wp )  t2 = 0.12_wp
1204
1205!--             *VOCL STMT,IF(10)
1206                IF ( m1(k-1,j) == 1.0_wp .OR. m1(k,j) == 1.0_wp                   &
1207                     .OR. m1(k+1,j) == 1.0_wp .OR.  m2 > t2  .OR.  m3 > t2  .OR.  &
1208                     ( m1(k,j) > t1  .AND.  m1(k-1,j) /= -1.0_wp  .AND.           &
1209                       m1(k,j) /= -1.0_wp  .AND.  m1(k+1,j) /= -1.0_wp )          &
1210                   )  sw(k,j) = 1.0_wp
1211             ENDDO
1212          ENDDO
1213
1214!
1215!--       Fluxes using exponential scheme
1216          CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
1217          DO  j = nys, nyn
1218             DO  k = nzb+1, nzt
1219
1220!--             *VOCL STMT,IF(10)
1221                IF ( sw(k,j) == 1.0_wp )  THEN
1222                   snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i)
1223                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
1224                   sterm = ( sk_p(k,j,i) - sk_p(k-1,j,i) ) / snenn
1225                   sterm = MIN( sterm, 0.9999_wp )
1226                   sterm = MAX( sterm, 0.0001_wp )
1227
1228                   ix = INT( sterm * 1000 ) + 1
1229
1230                   cip =  MAX( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
1231
1232                   ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * (                    &
1233                               aex(ix) * cip + bex(ix) / dex(ix) * (              &
1234                               eex(ix) -                                          &
1235                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
1236                                                                   )              &
1237                                                             )
1238                   IF ( sterm == 0.0001_wp )  ippe(k,j) = sk_p(k,j,i) * cip
1239                   IF ( sterm == 0.9999_wp )  ippe(k,j) = sk_p(k,j,i) * cip
1240
1241                   snenn = sk_p(k-1,j,i) - sk_p(k+1,j,i)
1242                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
1243                   sterm = ( sk_p(k,j,i) - sk_p(k+1,j,i) ) / snenn
1244                   sterm = MIN( sterm, 0.9999_wp )
1245                   sterm = MAX( sterm, 0.0001_wp )
1246
1247                   ix = INT( sterm * 1000 ) + 1
1248
1249                   cim = -MIN( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
1250
1251                   imme(k,j) = sk_p(k+1,j,i) * cim + snenn * (                    &
1252                               aex(ix) * cim + bex(ix) / dex(ix) * (              &
1253                               eex(ix) -                                          &
1254                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) ) &
1255                                                                   )              &
1256                                                             )
1257                   IF ( sterm == 0.0001_wp )  imme(k,j) = sk_p(k,j,i) * cim
1258                   IF ( sterm == 0.9999_wp )  imme(k,j) = sk_p(k,j,i) * cim
1259                ENDIF
1260
1261!--             *VOCL STMT,IF(10)
1262                IF ( sw(k+1,j) == 1.0_wp )  THEN
1263                   snenn = sk_p(k,j,i) - sk_p(k+2,j,i)
1264                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
1265                   sterm = ( sk_p(k+1,j,i) - sk_p(k+2,j,i) ) / snenn
1266                   sterm = MIN( sterm, 0.9999_wp )
1267                   sterm = MAX( sterm, 0.0001_wp )
1268
1269                   ix = INT( sterm * 1000 ) + 1
1270
1271                   cim = -MIN( 0.0_wp, w(k,j,i) * dt_3d * ddzw(k) )
1272
1273                   impe(k,j) = sk_p(k+2,j,i) * cim + snenn * (                    &
1274                               aex(ix) * cim + bex(ix) / dex(ix) * (              &
1275                               eex(ix) -                                          &
1276                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cim ) )  &
1277                                                                   )              &
1278                                                             )
1279                   IF ( sterm == 0.0001_wp )  impe(k,j) = sk_p(k+1,j,i) * cim
1280                   IF ( sterm == 0.9999_wp )  impe(k,j) = sk_p(k+1,j,i) * cim
1281                ENDIF
1282
1283!--             *VOCL STMT,IF(10)
1284                IF ( sw(k-1,j) == 1.0_wp )  THEN
1285                   snenn = sk_p(k,j,i) - sk_p(k-2,j,i)
1286                   IF ( ABS( snenn ) < 1E-9_wp )  snenn = 1E-9_wp
1287                   sterm = ( sk_p(k-1,j,i) - sk_p(k-2,j,i) ) / snenn
1288                   sterm = MIN( sterm, 0.9999_wp )
1289                   sterm = MAX( sterm, 0.0001_wp )
1290
1291                   ix = INT( sterm * 1000 ) + 1
1292
1293                   cip = MAX( 0.0_wp, w(k-1,j,i) * dt_3d * ddzw(k) )
1294
1295                   ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * (                    &
1296                               aex(ix) * cip + bex(ix) / dex(ix) * (              &
1297                               eex(ix) -                                          &
1298                               EXP( dex(ix)*0.5_wp * ( 1.0_wp - 2.0_wp * cip ) )  &
1299                                                                   )              &
1300                                                             )
1301                   IF ( sterm == 0.0001_wp )  ipme(k,j) = sk_p(k-1,j,i) * cip
1302                   IF ( sterm == 0.9999_wp )  ipme(k,j) = sk_p(k-1,j,i) * cip
1303                ENDIF
1304
1305             ENDDO
1306          ENDDO
1307          CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
1308
1309!
1310!--       Prognostic equation
1311          DO  j = nys, nyn
1312             DO  k = nzb+1, nzt
1313                fplus  = ( 1.0_wp - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
1314                       - ( 1.0_wp - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
1315                fminus = ( 1.0_wp - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) &
1316                       - ( 1.0_wp - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
1317                tendcy = fplus - fminus
1318!
1319!--              Removed in order to optimise speed
1320!                ffmax   = MAX( ABS( fplus ), ABS( fminus ), 1E-35_wp )
1321!                IF ( ( ABS( tendcy ) / ffmax ) < 1E-7_wp )  tendcy = 0.0
1322!
1323!--             Density correction because of possible remaining divergences
1324                d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k)
1325                sk_p(k,j,i) = ( ( 1.0_wp + d(k,j,i) ) * sk_p(k,j,i) - tendcy ) /  &
1326                              ( 1.0_wp + d_new )
1327!
1328!--             Store heat flux for subsequent statistics output.
1329!--             array m1 is here used as temporary storage
1330                m1(k,j) = fplus / dt_3d * dzw(k)
1331             ENDDO
1332          ENDDO
1333
1334!
1335!--       Sum up heat flux in order to order to obtain horizontal averages
1336          IF ( sk_char == 'pt' )  THEN
1337             DO  sr = 0, statistic_regions
1338                DO  j = nys, nyn
1339                   DO  k = nzb+1, nzt
1340                      sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) +               &
1341                                             m1(k,j) * rmask(j,i,sr)
1342                   ENDDO
1343                ENDDO
1344             ENDDO
1345          ENDIF
1346
1347       ENDDO   ! End of the advection in z-direction
1348       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
1349       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'stop' )
1350
1351!
1352!--    Deallocate temporary arrays
1353       DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme,      &
1354                   ippb, ippe, m1, sw )
1355
1356!
1357!--    Store results as tendency and deallocate local array
1358       DO  i = nxl, nxr
1359          DO  j = nys, nyn
1360             DO  k = nzb+1, nzt
1361                tend(k,j,i) = tend(k,j,i) + ( sk_p(k,j,i) - sk(k,j,i) ) / dt_3d
1362             ENDDO
1363          ENDDO
1364       ENDDO
1365
1366       DEALLOCATE( sk_p )
1367
1368    END SUBROUTINE advec_s_bc
1369
1370 END MODULE advec_s_bc_mod
Note: See TracBrowser for help on using the repository browser.