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

Last change on this file since 4488 was 4488, checked in by raasch, 3 years ago

files re-formatted to follow the PALM coding standard

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