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

Last change on this file since 1353 was 1353, checked in by heinze, 7 years ago

REAL constants provided with KIND-attribute

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