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

Last change on this file since 1361 was 1361, checked in by hoffmann, 10 years ago

improved version of two-moment cloud physics

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