source: palm/tags/release-3.1c/SOURCE/advec_s_bc.f90 @ 2516

Last change on this file since 2516 was 4, checked in by raasch, 17 years ago

Id keyword set as property for all *.f90 files

  • Property svn:keywords set to Id
File size: 42.6 KB
Line 
1 SUBROUTINE advec_s_bc( sk, sk_char )
2
3!-------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: advec_s_bc.f90 4 2007-02-13 11:33:16Z suehring $
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 cpulog
36    USE grid_variables
37    USE indices
38    USE interfaces
39    USE pegrid
40    USE statistics
41    USE control_parameters
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( __t3eb ) || defined( __t3eh ) || defined( __t3ej2 ) || defined( __t3ej5 )
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+2,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+1, 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 + 5 ) * ( 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-1
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-1,nxl-1:nxr+1),   a1(nzb+1:nzt-1,nxl-1:nxr+1),   &
163              a2(nzb+1:nzt-1,nxl-1:nxr+1),   a12(nzb+1:nzt-1,nxl-1:nxr+1),  &
164              a22(nzb+1:nzt-1,nxl-1:nxr+1),  immb(nzb+1:nzt-1,nxl-1:nxr+1), &
165              imme(nzb+1:nzt-1,nxl-1:nxr+1), impb(nzb+1:nzt-1,nxl-1:nxr+1), &
166              impe(nzb+1:nzt-1,nxl-1:nxr+1), ipmb(nzb+1:nzt-1,nxl-1:nxr+1), &
167              ipme(nzb+1:nzt-1,nxl-1:nxr+1), ippb(nzb+1:nzt-1,nxl-1:nxr+1), &
168              ippe(nzb+1:nzt-1,nxl-1:nxr+1), m1(nzb+1:nzt-1,nxl-2:nxr+2),   &
169              sw(nzb+1:nzt-1,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-1
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-1
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-1
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-1
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-1
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-1
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 + 5 ) * ( nyn - nys + 7 )
407    CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+5), 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-1
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-1
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-1,nys-1:nyn+1),   a1(nzb+1:nzt-1,nys-1:nyn+1),   &
460              a2(nzb+1:nzt-1,nys-1:nyn+1),   a12(nzb+1:nzt-1,nys-1:nyn+1),  &
461              a22(nzb+1:nzt-1,nys-1:nyn+1),  immb(nzb+1:nzt-1,nys-1:nyn+1), &
462              imme(nzb+1:nzt-1,nys-1:nyn+1), impb(nzb+1:nzt-1,nys-1:nyn+1), &
463              impe(nzb+1:nzt-1,nys-1:nyn+1), ipmb(nzb+1:nzt-1,nys-1:nyn+1), &
464              ipme(nzb+1:nzt-1,nys-1:nyn+1), ippb(nzb+1:nzt-1,nys-1:nyn+1), &
465              ippe(nzb+1:nzt-1,nys-1:nyn+1), m1(nzb+1:nzt-1,nys-2:nyn+2),   &
466              sw(nzb+1:nzt-1,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-1
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-1
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-1
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-1
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-1
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-1
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,j,i)   = sk_p(nzb+1,j,i)
723                sk_p(nzb-1,j,i) = sk_p(nzb+1,j,i)
724                sk_p(nzb-2,j,i) = sk_p(nzb+1,j,i)
725             ENDDO
726          ENDDO
727
728       ENDIF
729
730!
731!--    Temperature boundary condition at the top boundary
732       IF ( ibc_pt_t == 0 )  THEN
733!
734!--       Dirichlet
735          DO  i = nxl, nxr
736             DO  j = nys, nyn
737                sk_p(nzt+1,j,i)   = sk_p(nzt,j,i)
738                sk_p(nzt+2,j,i)   = sk_p(nzt,j,i)
739             ENDDO
740          ENDDO
741
742       ELSE
743!
744!--       Neumann: dzu(nzt+2) is not defined, thus instead dzu(nzt+1) is used
745          DO  i = nxl, nxr
746             DO  j = nys, nyn
747                sk_p(nzt,j,i)     = sk_p(nzt-1,j,i) + bc_pt_t_val * dzu(nzt)
748                sk_p(nzt+1,j,i)   = sk_p(nzt,j,i)   + bc_pt_t_val * dzu(nzt+1)
749                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_pt_t_val * dzu(nzt+1)
750             ENDDO
751          ENDDO
752
753       ENDIF
754
755    ELSEIF ( sk_char == 'q' )  THEN
756
757!
758!--    Temperature boundary condition at the bottom boundary
759       IF ( ibc_q_b == 0 )  THEN
760!
761!--       Dirichlet (fixed surface temperature)
762          DO  i = nxl, nxr
763             DO  j = nys, nyn
764                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
765                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
766             ENDDO
767          ENDDO
768
769       ELSE
770!
771!--       Neumann (i.e. here zero gradient)
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+1,j,i)
776                sk_p(nzb-2,j,i) = sk_p(nzb+1,j,i)
777             ENDDO
778          ENDDO
779
780       ENDIF
781
782!
783!--    Temperature boundary condition at the top boundary
784       IF ( ibc_q_t == 0 )  THEN
785!
786!--       Dirichlet
787          DO  i = nxl, nxr
788             DO  j = nys, nyn
789                sk_p(nzt+1,j,i)   = sk_p(nzt,j,i)
790                sk_p(nzt+2,j,i)   = sk_p(nzt,j,i)
791             ENDDO
792          ENDDO
793
794       ELSE
795!
796!--       Neumann: dzu(nzt+2) is not defined, thus instead dzu(nzt+1) is used
797          DO  i = nxl, nxr
798             DO  j = nys, nyn
799                sk_p(nzt,j,i)     = sk_p(nzt-1,j,i) + bc_q_t_val * dzu(nzt)
800                sk_p(nzt+1,j,i)   = sk_p(nzt,j,i)   + bc_q_t_val * dzu(nzt+1)
801                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_q_t_val * dzu(nzt+1)
802             ENDDO
803          ENDDO
804
805       ENDIF
806
807    ELSEIF ( sk_char == 'e' )  THEN
808
809!
810!--    TKE boundary condition at bottom and top boundary (generally Neumann)
811       sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
812       sk_p(nzb-1,j,i) = sk_p(nzb+1,j,i)
813       sk_p(nzb-2,j,i) = sk_p(nzb+1,j,i)
814       sk_p(nzt,j,i)   = sk_p(nzt-1,j,i)
815       sk_p(nzt+1,j,i) = sk_p(nzt-1,j,i)
816       sk_p(nzt+2,j,i) = sk_p(nzt-1,j,i)
817
818    ELSE
819
820       IF ( myid == 0 )  PRINT*,'+++ advec_s_bc:  no vertical boundary condi', &
821                                'tion for variable "', sk_char, '"'
822       CALL local_stop
823
824    ENDIF
825
826!
827!-- Determine the maxima of the first and second derivative in z-direction
828    fmax_l = 0.0
829    DO  i = nxl, nxr
830       DO  j = nys, nyn
831          DO  k = nzb, nzt
832             zaehler = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
833             nenner  = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) )
834             fmax_l(1) = MAX( fmax_l(1) , zaehler )
835             fmax_l(2) = MAX( fmax_l(2) , nenner  )
836          ENDDO
837       ENDDO
838    ENDDO
839#if defined( __parallel )
840    CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
841#else
842    fmax = fmax_l
843#endif
844
845    fmax = 0.04 * fmax
846
847!
848!-- Allocate temporary arrays
849    ALLOCATE( a0(nzb:nzt,nys:nyn),       a1(nzb:nzt,nys:nyn),       &
850              a2(nzb:nzt,nys:nyn),       a12(nzb:nzt,nys:nyn),      &
851              a22(nzb:nzt,nys:nyn),      immb(nzb+1:nzt-1,nys:nyn), &
852              imme(nzb+1:nzt-1,nys:nyn), impb(nzb+1:nzt-1,nys:nyn), &
853              impe(nzb+1:nzt-1,nys:nyn), ipmb(nzb+1:nzt-1,nys:nyn), &
854              ipme(nzb+1:nzt-1,nys:nyn), ippb(nzb+1:nzt-1,nys:nyn), &
855              ippe(nzb+1:nzt-1,nys:nyn), m1(nzb-1:nzt+1,nys:nyn),   &
856              sw(nzb:nzt,nys:nyn)                                   &
857            )
858    imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
859
860!
861!-- Outer loop of all i
862    DO  i = nxl, nxr
863
864!
865!--    Compute polynomial coefficients
866       DO  j = nys, nyn
867          DO  k = nzb, nzt
868             a12(k,j) = 0.5 * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
869             a22(k,j) = 0.5 * ( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) &
870                                              + sk_p(k-1,j,i) )
871             a0(k,j) = ( 9.0 * sk_p(k+2,j,i) - 116.0 * sk_p(k+1,j,i)    &
872                         + 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k-1,j,i) &
873                         + 9.0 * sk_p(k-2,j,i) ) * f1920
874             a1(k,j) = ( -5.0 * sk_p(k+2,j,i) + 34.0 * sk_p(k+1,j,i)  &
875                         - 34.0 * sk_p(k-1,j,i) + 5.0 * sk_p(k-2,j,i) &
876                       ) * f48
877             a2(k,j) = ( -3.0 * sk_p(k+2,j,i) + 36.0 * sk_p(k+1,j,i) &
878                         - 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k-1,j,i) &
879                         - 3.0 * sk_p(k-2,j,i) ) * f48
880          ENDDO
881       ENDDO
882
883!
884!--    Fluxes using the Bott scheme
885!--    *VOCL LOOP,UNROLL(2)
886       DO  j = nys, nyn
887          DO  k = nzb+1, nzt-1
888             cip  =  MAX( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
889             cim  = -MIN( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
890             cipf = 1.0 - 2.0 * cip
891             cimf = 1.0 - 2.0 * cim
892             ip   =   a0(k,j)   * f2  * ( 1.0 - cipf )             &
893                    + a1(k,j)   * f8  * ( 1.0 - cipf*cipf )        &
894                    + a2(k,j)   * f24 * ( 1.0 - cipf*cipf*cipf )
895             im   =   a0(k+1,j) * f2  * ( 1.0 - cimf )             &
896                    - a1(k+1,j) * f8  * ( 1.0 - cimf*cimf )        &
897                    + a2(k+1,j) * f24 * ( 1.0 - cimf*cimf*cimf )
898             ip   = MAX( ip, 0.0 )
899             im   = MAX( im, 0.0 )
900             ippb(k,j) = ip * MIN( 1.0, sk_p(k,j,i)   / (ip+im+1E-15) )
901             impb(k,j) = im * MIN( 1.0, sk_p(k+1,j,i) / (ip+im+1E-15) )
902
903             cip  =  MAX( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
904             cim  = -MIN( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
905             cipf = 1.0 - 2.0 * cip
906             cimf = 1.0 - 2.0 * cim
907             ip   =   a0(k-1,j) * f2  * ( 1.0 - cipf )             &
908                    + a1(k-1,j) * f8  * ( 1.0 - cipf*cipf )        &
909                    + a2(k-1,j) * f24 * ( 1.0 - cipf*cipf*cipf )
910             im   =   a0(k,j)   * f2  * ( 1.0 - cimf )             &
911                    - a1(k,j)   * f8  * ( 1.0 - cimf*cimf )        &
912                    + a2(k,j)   * f24 * ( 1.0 - cimf*cimf*cimf )
913             ip   = MAX( ip, 0.0 )
914             im   = MAX( im, 0.0 )
915             ipmb(k,j) = ip * MIN( 1.0, sk_p(k-1,j,i) / (ip+im+1E-15) )
916             immb(k,j) = im * MIN( 1.0, sk_p(k,j,i)   / (ip+im+1E-15) )
917          ENDDO
918       ENDDO
919
920!
921!--    Compute monitor function m1
922       DO  j = nys, nyn
923          DO  k = nzb-1, nzt+1
924             m1z = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
925             m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
926             IF ( m1n /= 0.0  .AND.  m1n >= m1z )  THEN
927                m1(k,j) = m1z / m1n
928                IF ( m1(k,j) /= 2.0  .AND.  m1n < fmax(2) )  m1(k,j) = 0.0
929             ELSEIF ( m1n < m1z )  THEN
930                m1(k,j) = -1.0
931             ELSE
932                m1(k,j) = 0.0
933             ENDIF
934          ENDDO
935       ENDDO
936
937!
938!--    Compute switch sw
939       sw = 0.0
940       DO  j = nys, nyn
941          DO  k = nzb, nzt
942             m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / &
943                  MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35 )
944             IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) )  m2 = 0.0
945
946             m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) / &
947                  MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35 )
948             IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) )  m3 = 0.0
949
950             t1 = 0.35
951             t2 = 0.35
952             IF ( m1(k,j) == -1.0 )  t2 = 0.12
953
954!--          *VOCL STMT,IF(10)
955             IF ( m1(k-1,j) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k+1,j) == 1.0 &
956                  .OR.  m2 > t2  .OR.  m3 > T2  .OR.                         &
957                  ( m1(k,j) > t1  .AND.  m1(k-1,j) /= -1.0  .AND.            &
958                    m1(k,j) /= -1.0  .AND.  m1(k+1,j) /= -1.0 )              &
959                )  sw(k,j) = 1.0
960          ENDDO
961       ENDDO
962
963!
964!--    Fluxes using exponential scheme
965       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
966       DO  j = nys, nyn
967          DO  k = nzb+1, nzt-1
968
969!--          *VOCL STMT,IF(10)
970             IF ( sw(k,j) == 1.0 )  THEN
971                snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i)
972                IF ( ABS( snenn ) < 1E-9 )  snenn = 1E-9
973                sterm = ( sk_p(k,j,i) - sk_p(k-1,j,i) ) / snenn
974                sterm = MIN( sterm, 0.9999 )
975                sterm = MAX( sterm, 0.0001 )
976
977                ix = INT( sterm * 1000 ) + 1
978
979                cip =  MAX( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
980
981                ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * (                    &
982                            aex(ix) * cip + bex(ix) / dex(ix) * (              &
983                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) &
984                                                                )              &
985                                                          )
986                IF ( sterm == 0.0001 )  ippe(k,j) = sk_p(k,j,i) * cip
987                IF ( sterm == 0.9999 )  ippe(k,j) = sk_p(k,j,i) * cip
988
989                snenn = sk_p(k-1,j,i) - sk_p(k+1,j,i)
990                IF ( ABS( snenn ) < 1E-9 )  snenn = 1E-9
991                sterm = ( sk_p(k,j,i) - sk_p(k+1,j,i) ) / snenn
992                sterm = MIN( sterm, 0.9999 )
993                sterm = MAX( sterm, 0.0001 )
994
995                ix = INT( sterm * 1000 ) + 1
996
997                cim = -MIN( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
998
999                imme(k,j) = sk_p(k+1,j,i) * cim + snenn * (                    &
1000                            aex(ix) * cim + bex(ix) / dex(ix) * (              &
1001                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) &
1002                                                                )              &
1003                                                          )
1004                IF ( sterm == 0.0001 )  imme(k,j) = sk_p(k,j,i) * cim
1005                IF ( sterm == 0.9999 )  imme(k,j) = sk_p(k,j,i) * cim
1006             ENDIF
1007
1008!--          *VOCL STMT,IF(10)
1009             IF ( sw(k+1,j) == 1.0 )  THEN
1010                snenn = sk_p(k,j,i) - sk_p(k+2,j,i)
1011                IF ( ABS( snenn ) .LT. 1E-9 )  snenn = 1E-9
1012                sterm = ( sk_p(k+1,j,i) - sk_p(k+2,j,i) ) / snenn
1013                sterm = MIN( sterm, 0.9999 )
1014                sterm = MAX( sterm, 0.0001 )
1015
1016                ix = INT( sterm * 1000 ) + 1
1017
1018                cim = -MIN( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
1019
1020                impe(k,j) = sk_p(k+2,j,i) * cim + snenn * (                    &
1021                            aex(ix) * cim + bex(ix) / dex(ix) * (              &
1022                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cim ) ) &
1023                                                                )              &
1024                                                          )
1025                IF ( sterm == 0.0001 )  impe(k,j) = sk_p(k+1,j,i) * cim
1026                IF ( sterm == 0.9999 )  impe(k,j) = sk_p(k+1,j,i) * cim
1027             ENDIF
1028
1029!--          *VOCL STMT,IF(10)
1030             IF ( sw(k-1,j) == 1.0 )  THEN
1031                snenn = sk_p(k,j,i) - sk_p(k-2,j,i)
1032                IF ( ABS( snenn ) < 1E-9 )  snenn = 1E-9
1033                sterm = ( sk_p(k-1,j,i) - sk_p(k-2,j,i) ) / snenn
1034                sterm = MIN( sterm, 0.9999 )
1035                sterm = MAX( sterm, 0.0001 )
1036
1037                ix = INT( sterm * 1000 ) + 1
1038
1039                cip = MAX( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
1040
1041                ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * (                    &
1042                            aex(ix) * cip + bex(ix) / dex(ix) * (              &
1043                            eex(ix) - EXP( dex(ix)*0.5 * ( 1.0 - 2.0 * cip ) ) &
1044                                                                )              &
1045                                                          )
1046                IF ( sterm == 0.0001 )  ipme(k,j) = sk_p(k-1,j,i) * cip
1047                IF ( sterm == 0.9999 )  ipme(k,j) = sk_p(k-1,j,i) * cip
1048             ENDIF
1049
1050          ENDDO
1051       ENDDO
1052       CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
1053
1054!
1055!--    Prognostic equation
1056       DO  j = nys, nyn
1057          DO  k = nzb+1, nzt-1
1058             fplus  = ( 1.0 - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
1059                    - ( 1.0 - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
1060             fminus = ( 1.0 - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) &
1061                    - ( 1.0 - sw(k,j)   ) * immb(k,j) - sw(k,j)   * imme(k,j)
1062             tendenz = fplus - fminus
1063!
1064!--           Removed in order to optimise speed
1065!             ffmax   = MAX( ABS( fplus ), ABS( fminus ), 1E-35 )
1066!             IF ( ( ABS( tendenz ) / ffmax ) < 1E-7 )  tendenz = 0.0
1067!
1068!--          Density correction because of possible remaining divergences
1069             d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k)
1070             sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendenz ) / &
1071                           ( 1.0 + d_new )
1072!
1073!--          Store heat flux for subsequent statistics output.
1074!--          array m1 is here used as temporary storage
1075             m1(k,j) = fplus / dt_3d * dzw(k)
1076          ENDDO
1077       ENDDO
1078
1079!
1080!--    Sum up heat flux in order to order to obtain horizontal averages
1081       IF ( sk_char == 'pt' )  THEN
1082          DO  sr = 0, statistic_regions
1083             DO  j = nys, nyn
1084                DO  k = nzb+1, nzt-1
1085                   sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) + &
1086                                          m1(k,j) * rmask(j,i,sr)
1087                ENDDO
1088             ENDDO
1089          ENDDO
1090       ENDIF
1091
1092    ENDDO   ! End of the advection in z-direction
1093    CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
1094    CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'stop' )
1095
1096!
1097!-- Deallocate temporary arrays
1098    DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, &
1099                ippb, ippe, m1, sw )
1100
1101!
1102!-- Store results as tendency and deallocate local array
1103    DO  i = nxl, nxr
1104       DO  j = nys, nyn
1105          DO  k = nzb+1, nzt-1
1106             tend(k,j,i) = tend(k,j,i) + ( sk_p(k,j,i) - sk(k,j,i) ) / dt_3d
1107          ENDDO
1108       ENDDO
1109    ENDDO
1110
1111    DEALLOCATE( sk_p )
1112
1113 END SUBROUTINE advec_s_bc
Note: See TracBrowser for help on using the repository browser.