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

Last change on this file since 1379 was 1375, checked in by raasch, 11 years ago

last commit documented

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