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

Last change on this file since 1754 was 1692, checked in by maronga, 9 years ago

last commit documented

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