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

Last change on this file since 1623 was 1518, checked in by hoffmann, 10 years ago

last commit documented

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