source: palm/trunk/SOURCE/advec_s_bc_mod.f90 @ 1851

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

last commit documented

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