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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

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