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

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

revised renaming of modules

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