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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

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