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

Last change on this file since 1011 was 1011, checked in by raasch, 9 years ago

last commit documented

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