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

Last change on this file since 4360 was 4360, checked in by suehring, 4 years ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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