source: palm/trunk/SOURCE/advec_s_bc_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 8 years ago

added _mod string to several filenames to meet the naming convection for modules

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