source: palm/tags/release-3.10/SOURCE/advec_s_bc.f90 @ 3967

Last change on this file since 3967 was 1093, checked in by raasch, 11 years ago

last commit documented

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