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

Last change on this file since 354 was 247, checked in by heinze, 15 years ago

Output of messages replaced by message handling routin

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