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

Last change on this file since 241 was 226, checked in by raasch, 15 years ago

preparations for the next release

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