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

Last change on this file since 4178 was 3761, checked in by raasch, 2 years ago

unused variables removed, OpenACC directives re-formatted, statements added to avoid compiler warnings

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