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

Last change on this file since 1517 was 1517, checked in by hoffmann, 7 years ago

bug-fix: Bott-Chlond scheme

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