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

Last change on this file since 4631 was 4502, checked in by schwenkel, 5 years ago

Implementation of ice microphysics

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