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

Last change on this file since 1375 was 1375, checked in by raasch, 8 years ago

last commit documented

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