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

Last change on this file since 1374 was 1374, checked in by raasch, 7 years ago

bugfixes:
missing variables added to ONLY lists in USE statements (advec_s_bc, advec_s_pw, advec_s_up, advec_ws, buoyancy, diffusion_e, diffusion_s, fft_xy, flow_statistics, palm, prognostic_equations)
missing module kinds added (cuda_fft_interfaces)
dpk renamed dp (fft_xy)
missing dependency for check_open added (Makefile)
variables removed from acc-present-list (diffusion_e, diffusion_w, diffusivities, production_e, wall_fluxes)
syntax errors removed from openacc-branch (flow_statistics)
USE-statement for nopointer-case added (swap_timelevel)

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