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

Last change on this file since 622 was 622, checked in by raasch, 11 years ago

New:
---

Optional barriers included in order to speed up collective operations
MPI_ALLTOALL and MPI_ALLREDUCE. This feature is controlled with new initial
parameter collective_wait. Default is .FALSE, but .TRUE. on SGI-type
systems. (advec_particles, advec_s_bc, buoyancy, check_for_restart,
cpu_statistics, data_output_2d, data_output_ptseries, flow_statistics,
global_min_max, inflow_turbulence, init_3d_model, init_particles, init_pegrid,
init_slope, parin, pres, poismg, set_particle_attributes, timestep,
read_var_list, user_statistics, write_compressed, write_var_list)

Adjustments for Kyushu Univ. (lcrte, ibmku). Concerning hybrid
(MPI/openMP) runs, the number of openMP threads per MPI tasks can now
be given as an argument to mrun-option -O. (mbuild, mrun, subjob)

Changed:


Initialization of the module command changed for SGI-ICE/lcsgi (mbuild, subjob)

Errors:


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