Home


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

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

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

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