Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

Last change on this file since 92 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:

General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:

Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

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