source: palm/tags/release-3.8/SOURCE/advec_s_bc.f90 @ 4343

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

last commit documented

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