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

Last change on this file since 1874 was 1874, checked in by maronga, 5 years ago

last commit documented

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