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

Last change on this file since 1318 was 1318, checked in by raasch, 8 years ago

former files/routines cpu_log and cpu_statistics combined to one module,
which also includes the former data module cpulog from the modules-file,
module interfaces removed

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