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

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

Initial repository layout and content

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