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

Last change on this file since 2101 was 2101, checked in by suehring, 5 years ago

last commit documented

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