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

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

updating comments and rc-file

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