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

Last change on this file since 63 was 63, checked in by raasch, 14 years ago

preliminary changes concerning update of BC-scheme

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