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

Last change on this file since 2716 was 2716, checked in by kanani, 6 years ago

Correction of "Former revisions" section

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