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

Last change on this file since 1822 was 1818, checked in by maronga, 8 years ago

last commit documented / copyright update

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