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

Last change on this file since 3746 was 3655, checked in by knoop, 6 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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