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

Last change on this file since 4444 was 4429, checked in by raasch, 4 years ago

serial (non-MPI) test case added, several bugfixes for the serial mode

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