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

Last change on this file since 1727 was 1692, checked in by maronga, 9 years ago

last commit documented

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