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

Last change on this file since 2416 was 2300, checked in by raasch, 7 years ago

NEC related code partly removed, host variable partly removed, host specific code completely removed, default values for host, loop_optimization and termination time_needed changed

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