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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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