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

Last change on this file since 3436 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

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