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

Last change on this file since 2292 was 2292, checked in by schwenkel, 4 years ago

implementation of new bulk microphysics scheme

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