Home


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

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

Last change on this file since 547 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:

lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:

Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

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