Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

Last change on this file since 1329 was 1321, checked in by raasch, 10 years ago

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