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

Last change on this file since 4782 was 4742, checked in by schwenkel, 4 years ago

Implement snow and graupel (bulk microphysics)

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