Home


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

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

Last change on this file since 1347 was 1347, checked in by heinze, 10 years ago

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