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

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

last commit documented

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