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

Last change on this file since 224 was 216, checked in by raasch, 16 years ago

reading mechanism for restart files completely revised

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