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

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