source: palm/trunk/SOURCE/data_output_2d.f90 @ 399

Last change on this file since 399 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: 48.8 KB
Line 
1 SUBROUTINE data_output_2d( mode, av )
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: data_output_2d.f90 392 2009-09-24 10:39:14Z weinreis $
11!
12! 367 2009-08-25 08:35:52Z maronga
13! simulated_time in NetCDF output replaced by time_since_reference_point.
14! Output of NetCDF messages with aid of message handling routine.
15! Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0)
16! Output of messages replaced by message handling routine.
17! Output of user defined 2D (XY) arrays at z=nzb+1 is now possible
18! Bugfix: to_be_resorted => s_av for time-averaged scalars
19! Calculation of shf* and qsws* added.
20!
21! 215 2008-11-18 09:54:31Z raasch
22! Bugfix: no output of particle concentration and radius unless particles
23! have been started
24!
25! 96 2007-06-04 08:07:41Z raasch
26! Output of density and salinity
27!
28! 75 2007-03-22 09:54:05Z raasch
29! Output of precipitation amount/rate and roughness length,
30! 2nd+3rd argument removed from exchange horiz
31!
32! RCS Log replace by Id keyword, revision history cleaned up
33!
34! Revision 1.5  2006/08/22 13:50:29  raasch
35! xz and yz cross sections now up to nzt+1
36!
37! Revision 1.2  2006/02/23 10:19:22  raasch
38! Output of time-averaged data, output of averages along x, y, or z,
39! output of user-defined quantities,
40! section data are copied from local_pf to local_2d before they are output,
41! output of particle concentration and mean radius,
42! Former subroutine plot_2d renamed data_output_2d, pl2d.. renamed do2d..,
43! anz renamed ngp, ebene renamed section, pl2d_.._anz renamed do2d_.._n
44!
45! Revision 1.1  1997/08/11 06:24:09  raasch
46! Initial revision
47!
48!
49! Description:
50! ------------
51! Data output of horizontal cross-sections in NetCDF format or binary format
52! compatible to old graphic software iso2d.
53! Attention: The position of the sectional planes is still not always computed
54! ---------  correctly. (zu is used always)!
55!------------------------------------------------------------------------------!
56
57    USE arrays_3d
58    USE averaging
59    USE cloud_parameters
60    USE control_parameters
61    USE cpulog
62    USE grid_variables
63    USE indices
64    USE interfaces
65    USE netcdf_control
66    USE particle_attributes
67    USE pegrid
68
69    IMPLICIT NONE
70
71    CHARACTER (LEN=2)  ::  do2d_mode, mode
72    CHARACTER (LEN=4)  ::  grid
73    CHARACTER (LEN=25) ::  section_chr
74    CHARACTER (LEN=50) ::  rtext
75    INTEGER ::  av, ngp, file_id, i, if, is, j, k, l, layer_xy, n, psi, s, &
76                sender, &
77                ind(4)
78    LOGICAL ::  found, resorted, two_d
79    REAL    ::  mean_r, s_r3, s_r4
80    REAL, DIMENSION(:), ALLOCATABLE ::      level_z
81    REAL, DIMENSION(:,:), ALLOCATABLE ::    local_2d, local_2d_l
82    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  local_pf
83#if defined( __parallel )
84    REAL, DIMENSION(:,:),   ALLOCATABLE ::  total_2d
85#endif
86    REAL, DIMENSION(:,:,:), POINTER ::  to_be_resorted
87
88    NAMELIST /LOCAL/  rtext
89
90    CALL cpu_log (log_point(3),'data_output_2d','start')
91
92!
93!-- Immediate return, if no output is requested (no respective sections
94!-- found in parameter data_output)
95    IF ( mode == 'xy'  .AND.  .NOT. data_output_xy(av) )  RETURN
96    IF ( mode == 'xz'  .AND.  .NOT. data_output_xz(av) )  RETURN
97    IF ( mode == 'yz'  .AND.  .NOT. data_output_yz(av) )  RETURN
98
99    two_d = .FALSE.    ! local variable to distinguish between output of pure 2D
100                       ! arrays and cross-sections of 3D arrays.
101
102!
103!-- Depending on the orientation of the cross-section, the respective output
104!-- files have to be opened.
105    SELECT CASE ( mode )
106
107       CASE ( 'xy' )
108
109          s = 1
110          ALLOCATE( level_z(0:nzt+1), local_2d(nxl-1:nxr+1,nys-1:nyn+1) )
111
112#if defined( __netcdf )
113          IF ( myid == 0  .AND.  netcdf_output )  CALL check_open( 101+av*10 )
114#endif
115
116          IF ( data_output_2d_on_each_pe )  THEN
117             CALL check_open( 21 )
118          ELSE
119             IF ( myid == 0 )  THEN
120                IF ( iso2d_output )  CALL check_open( 21 )
121#if defined( __parallel )
122                ALLOCATE( total_2d(-1:nx+1,-1:ny+1) )
123#endif
124             ENDIF
125          ENDIF
126
127       CASE ( 'xz' )
128
129          s = 2
130          ALLOCATE( local_2d(nxl-1:nxr+1,nzb:nzt+1) )
131
132#if defined( __netcdf )
133          IF ( myid == 0  .AND.  netcdf_output )  CALL check_open( 102+av*10 )
134#endif
135
136          IF ( data_output_2d_on_each_pe )  THEN
137             CALL check_open( 22 )
138          ELSE
139             IF ( myid == 0 )  THEN
140                IF ( iso2d_output )  CALL check_open( 22 )
141#if defined( __parallel )
142                ALLOCATE( total_2d(-1:nx+1,nzb:nzt+1) )
143#endif
144             ENDIF
145          ENDIF
146
147       CASE ( 'yz' )
148
149          s = 3
150          ALLOCATE( local_2d(nys-1:nyn+1,nzb:nzt+1) )
151
152#if defined( __netcdf )
153          IF ( myid == 0  .AND.  netcdf_output )  CALL check_open( 103+av*10 )
154#endif
155
156          IF ( data_output_2d_on_each_pe )  THEN
157             CALL check_open( 23 )
158          ELSE
159             IF ( myid == 0 )  THEN
160                IF ( iso2d_output )  CALL check_open( 23 )
161#if defined( __parallel )
162                ALLOCATE( total_2d(-1:ny+1,nzb:nzt+1) )
163#endif
164             ENDIF
165          ENDIF
166
167       CASE DEFAULT
168
169          message_string = 'unknown cross-section: ' // TRIM( mode )
170          CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 )
171
172    END SELECT
173
174!
175!-- Allocate a temporary array for resorting (kji -> ijk).
176    ALLOCATE( local_pf(nxl-1:nxr+1,nys-1:nyn+1,nzb:nzt+1) )
177
178!
179!-- Loop of all variables to be written.
180!-- Output dimensions chosen
181    if = 1
182    l = MAX( 2, LEN_TRIM( do2d(av,if) ) )
183    do2d_mode = do2d(av,if)(l-1:l)
184
185    DO  WHILE ( do2d(av,if)(1:1) /= ' ' )
186
187       IF ( do2d_mode == mode )  THEN
188!
189!--       Store the array chosen on the temporary array.
190          resorted = .FALSE.
191          SELECT CASE ( TRIM( do2d(av,if) ) )
192
193             CASE ( 'e_xy', 'e_xz', 'e_yz' )
194                IF ( av == 0 )  THEN
195                   to_be_resorted => e
196                ELSE
197                   to_be_resorted => e_av
198                ENDIF
199                IF ( mode == 'xy' )  level_z = zu
200
201             CASE ( 'lwp*_xy' )        ! 2d-array
202                IF ( av == 0 )  THEN
203                   DO  i = nxl-1, nxr+1
204                      DO  j = nys-1, nyn+1
205                         local_pf(i,j,nzb+1) = SUM( ql(nzb:nzt,j,i) * &
206                                                    dzw(1:nzt+1) )
207                      ENDDO
208                   ENDDO
209                ELSE
210                   DO  i = nxl-1, nxr+1
211                      DO  j = nys-1, nyn+1
212                         local_pf(i,j,nzb+1) = lwp_av(j,i)
213                      ENDDO
214                   ENDDO
215                ENDIF
216                resorted = .TRUE.
217                two_d = .TRUE.
218                level_z(nzb+1) = zu(nzb+1)
219
220             CASE ( 'p_xy', 'p_xz', 'p_yz' )
221                IF ( av == 0 )  THEN
222                   to_be_resorted => p
223                ELSE
224                   to_be_resorted => p_av
225                ENDIF
226                IF ( mode == 'xy' )  level_z = zu
227
228             CASE ( 'pc_xy', 'pc_xz', 'pc_yz' )  ! particle concentration
229                IF ( av == 0 )  THEN
230                   IF ( simulated_time >= particle_advection_start )  THEN
231                      tend = prt_count
232                      CALL exchange_horiz( tend )
233                   ELSE
234                      tend = 0.0
235                   ENDIF
236                   DO  i = nxl-1, nxr+1
237                      DO  j = nys-1, nyn+1
238                         DO  k = nzb, nzt+1
239                            local_pf(i,j,k) = tend(k,j,i)
240                         ENDDO
241                      ENDDO
242                   ENDDO
243                   resorted = .TRUE.
244                ELSE
245                   CALL exchange_horiz( pc_av )
246                   to_be_resorted => pc_av
247                ENDIF
248
249             CASE ( 'pr_xy', 'pr_xz', 'pr_yz' )  ! mean particle radius
250                IF ( av == 0 )  THEN
251                   IF ( simulated_time >= particle_advection_start )  THEN
252                      DO  i = nxl, nxr
253                         DO  j = nys, nyn
254                            DO  k = nzb, nzt+1
255                               psi = prt_start_index(k,j,i)
256                               s_r3 = 0.0
257                               s_r4 = 0.0
258                               DO  n = psi, psi+prt_count(k,j,i)-1
259                                  s_r3 = s_r3 + particles(n)%radius**3
260                                  s_r4 = s_r4 + particles(n)%radius**4
261                               ENDDO
262                               IF ( s_r3 /= 0.0 )  THEN
263                                  mean_r = s_r4 / s_r3
264                               ELSE
265                                  mean_r = 0.0
266                               ENDIF
267                               tend(k,j,i) = mean_r
268                            ENDDO
269                         ENDDO
270                      ENDDO
271                      CALL exchange_horiz( tend )
272                   ELSE
273                      tend = 0.0
274                   ENDIF
275                   DO  i = nxl-1, nxr+1
276                      DO  j = nys-1, nyn+1
277                         DO  k = nzb, nzt+1
278                            local_pf(i,j,k) = tend(k,j,i)
279                         ENDDO
280                      ENDDO
281                   ENDDO
282                   resorted = .TRUE.
283                ELSE
284                   CALL exchange_horiz( pr_av )
285                   to_be_resorted => pr_av
286                ENDIF
287
288             CASE ( 'pra*_xy' )        ! 2d-array / integral quantity => no av
289                CALL exchange_horiz_2d( precipitation_amount )
290                DO  i = nxl-1, nxr+1
291                   DO  j = nys-1, nyn+1
292                      local_pf(i,j,nzb+1) =  precipitation_amount(j,i)
293                   ENDDO
294                ENDDO
295                precipitation_amount = 0.0   ! reset for next integ. interval
296                resorted = .TRUE.
297                two_d = .TRUE.
298                level_z(nzb+1) = zu(nzb+1)
299
300             CASE ( 'prr*_xy' )        ! 2d-array
301                IF ( av == 0 )  THEN
302                   CALL exchange_horiz_2d( precipitation_rate )
303                   DO  i = nxl-1, nxr+1
304                      DO  j = nys-1, nyn+1 
305                         local_pf(i,j,nzb+1) =  precipitation_rate(j,i)
306                      ENDDO
307                   ENDDO
308                ELSE
309                   CALL exchange_horiz_2d( precipitation_rate_av )
310                   DO  i = nxl-1, nxr+1
311                      DO  j = nys-1, nyn+1 
312                         local_pf(i,j,nzb+1) =  precipitation_rate_av(j,i)
313                      ENDDO
314                   ENDDO
315                ENDIF
316                resorted = .TRUE.
317                two_d = .TRUE.
318                level_z(nzb+1) = zu(nzb+1)
319
320             CASE ( 'pt_xy', 'pt_xz', 'pt_yz' )
321                IF ( av == 0 )  THEN
322                   IF ( .NOT. cloud_physics ) THEN
323                      to_be_resorted => pt
324                   ELSE
325                      DO  i = nxl-1, nxr+1
326                         DO  j = nys-1, nyn+1
327                            DO  k = nzb, nzt+1
328                               local_pf(i,j,k) = pt(k,j,i) + l_d_cp *    &
329                                                             pt_d_t(k) * &
330                                                             ql(k,j,i)
331                            ENDDO
332                         ENDDO
333                      ENDDO
334                      resorted = .TRUE.
335                   ENDIF
336                ELSE
337                   to_be_resorted => pt_av
338                ENDIF
339                IF ( mode == 'xy' )  level_z = zu
340
341             CASE ( 'q_xy', 'q_xz', 'q_yz' )
342                IF ( av == 0 )  THEN
343                   to_be_resorted => q
344                ELSE
345                   to_be_resorted => q_av
346                ENDIF
347                IF ( mode == 'xy' )  level_z = zu
348
349             CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
350                IF ( av == 0 )  THEN
351                   to_be_resorted => ql
352                ELSE
353                   to_be_resorted => ql_av
354                ENDIF
355                IF ( mode == 'xy' )  level_z = zu
356
357             CASE ( 'ql_c_xy', 'ql_c_xz', 'ql_c_yz' )
358                IF ( av == 0 )  THEN
359                   to_be_resorted => ql_c
360                ELSE
361                   to_be_resorted => ql_c_av
362                ENDIF
363                IF ( mode == 'xy' )  level_z = zu
364
365             CASE ( 'ql_v_xy', 'ql_v_xz', 'ql_v_yz' )
366                IF ( av == 0 )  THEN
367                   to_be_resorted => ql_v
368                ELSE
369                   to_be_resorted => ql_v_av
370                ENDIF
371                IF ( mode == 'xy' )  level_z = zu
372
373             CASE ( 'ql_vp_xy', 'ql_vp_xz', 'ql_vp_yz' )
374                IF ( av == 0 )  THEN
375                   to_be_resorted => ql_vp
376                ELSE
377                   to_be_resorted => ql_vp_av
378                ENDIF
379                IF ( mode == 'xy' )  level_z = zu
380
381             CASE ( 'qsws*_xy' )        ! 2d-array
382                IF ( av == 0 ) THEN
383                   DO  i = nxl-1, nxr+1
384                      DO  j = nys-1, nyn+1 
385                         local_pf(i,j,nzb+1) =  qsws(j,i)
386                      ENDDO
387                   ENDDO
388                ELSE
389                   DO  i = nxl-1, nxr+1
390                      DO  j = nys-1, nyn+1 
391                         local_pf(i,j,nzb+1) =  qsws_av(j,i)
392                      ENDDO
393                   ENDDO
394                ENDIF
395                resorted = .TRUE.
396                two_d = .TRUE.
397                level_z(nzb+1) = zu(nzb+1)
398
399             CASE ( 'qv_xy', 'qv_xz', 'qv_yz' )
400                IF ( av == 0 )  THEN
401                   DO  i = nxl-1, nxr+1
402                      DO  j = nys-1, nyn+1
403                         DO  k = nzb, nzt+1
404                            local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
405                         ENDDO
406                      ENDDO
407                   ENDDO
408                   resorted = .TRUE.
409                ELSE
410                   to_be_resorted => qv_av
411                ENDIF
412                IF ( mode == 'xy' )  level_z = zu
413
414             CASE ( 'rho_xy', 'rho_xz', 'rho_yz' )
415                IF ( av == 0 )  THEN
416                   to_be_resorted => rho
417                ELSE
418                   to_be_resorted => rho_av
419                ENDIF
420
421             CASE ( 's_xy', 's_xz', 's_yz' )
422                IF ( av == 0 )  THEN
423                   to_be_resorted => q
424                ELSE
425                   to_be_resorted => s_av
426                ENDIF
427
428             CASE ( 'sa_xy', 'sa_xz', 'sa_yz' )
429                IF ( av == 0 )  THEN
430                   to_be_resorted => sa
431                ELSE
432                   to_be_resorted => sa_av
433                ENDIF
434
435             CASE ( 'shf*_xy' )        ! 2d-array
436                IF ( av == 0 ) THEN
437                   DO  i = nxl-1, nxr+1
438                      DO  j = nys-1, nyn+1 
439                         local_pf(i,j,nzb+1) =  shf(j,i)
440                      ENDDO
441                   ENDDO
442                ELSE
443                   DO  i = nxl-1, nxr+1
444                      DO  j = nys-1, nyn+1 
445                         local_pf(i,j,nzb+1) =  shf_av(j,i)
446                      ENDDO
447                   ENDDO
448                ENDIF
449                resorted = .TRUE.
450                two_d = .TRUE.
451                level_z(nzb+1) = zu(nzb+1)
452
453             CASE ( 't*_xy' )        ! 2d-array
454                IF ( av == 0 )  THEN
455                   DO  i = nxl-1, nxr+1
456                      DO  j = nys-1, nyn+1
457                         local_pf(i,j,nzb+1) = ts(j,i)
458                      ENDDO
459                   ENDDO
460                ELSE
461                   DO  i = nxl-1, nxr+1
462                      DO  j = nys-1, nyn+1
463                         local_pf(i,j,nzb+1) = ts_av(j,i)
464                      ENDDO
465                   ENDDO
466                ENDIF
467                resorted = .TRUE.
468                two_d = .TRUE.
469                level_z(nzb+1) = zu(nzb+1)
470
471             CASE ( 'u_xy', 'u_xz', 'u_yz' )
472                IF ( av == 0 )  THEN
473                   to_be_resorted => u
474                ELSE
475                   to_be_resorted => u_av
476                ENDIF
477                IF ( mode == 'xy' )  level_z = zu
478!
479!--             Substitute the values generated by "mirror" boundary condition
480!--             at the bottom boundary by the real surface values.
481                IF ( do2d(av,if) == 'u_xz'  .OR.  do2d(av,if) == 'u_yz' )  THEN
482                   IF ( ibc_uv_b == 0 )  local_pf(:,:,nzb) = 0.0
483                ENDIF
484
485             CASE ( 'u*_xy' )        ! 2d-array
486                IF ( av == 0 )  THEN
487                   DO  i = nxl-1, nxr+1
488                      DO  j = nys-1, nyn+1
489                         local_pf(i,j,nzb+1) = us(j,i)
490                      ENDDO
491                   ENDDO
492                ELSE
493                   DO  i = nxl-1, nxr+1
494                      DO  j = nys-1, nyn+1
495                         local_pf(i,j,nzb+1) = us_av(j,i)
496                      ENDDO
497                   ENDDO
498                ENDIF
499                resorted = .TRUE.
500                two_d = .TRUE.
501                level_z(nzb+1) = zu(nzb+1)
502
503             CASE ( 'v_xy', 'v_xz', 'v_yz' )
504                IF ( av == 0 )  THEN
505                   to_be_resorted => v
506                ELSE
507                   to_be_resorted => v_av
508                ENDIF
509                IF ( mode == 'xy' )  level_z = zu
510!
511!--             Substitute the values generated by "mirror" boundary condition
512!--             at the bottom boundary by the real surface values.
513                IF ( do2d(av,if) == 'v_xz'  .OR.  do2d(av,if) == 'v_yz' )  THEN
514                   IF ( ibc_uv_b == 0 )  local_pf(:,:,nzb) = 0.0
515                ENDIF
516
517             CASE ( 'vpt_xy', 'vpt_xz', 'vpt_yz' )
518                IF ( av == 0 )  THEN
519                   to_be_resorted => vpt
520                ELSE
521                   to_be_resorted => vpt_av
522                ENDIF
523                IF ( mode == 'xy' )  level_z = zu
524
525             CASE ( 'w_xy', 'w_xz', 'w_yz' )
526                IF ( av == 0 )  THEN
527                   to_be_resorted => w
528                ELSE
529                   to_be_resorted => w_av
530                ENDIF
531                IF ( mode == 'xy' )  level_z = zw
532
533             CASE ( 'z0*_xy' )        ! 2d-array
534                IF ( av == 0 ) THEN
535                   DO  i = nxl-1, nxr+1
536                      DO  j = nys-1, nyn+1 
537                         local_pf(i,j,nzb+1) =  z0(j,i)
538                      ENDDO
539                   ENDDO
540                ELSE
541                   DO  i = nxl-1, nxr+1
542                      DO  j = nys-1, nyn+1 
543                         local_pf(i,j,nzb+1) =  z0_av(j,i)
544                      ENDDO
545                   ENDDO
546                ENDIF
547                resorted = .TRUE.
548                two_d = .TRUE.
549                level_z(nzb+1) = zu(nzb+1)
550
551             CASE DEFAULT
552!
553!--             User defined quantity
554                CALL user_data_output_2d( av, do2d(av,if), found, grid, &
555                                          local_pf, two_d )
556                resorted = .TRUE.
557
558                IF ( grid == 'zu' )  THEN
559                   IF ( mode == 'xy' )  level_z = zu
560                ELSEIF ( grid == 'zw' )  THEN
561                   IF ( mode == 'xy' )  level_z = zw
562                ELSEIF ( grid == 'zu1' ) THEN
563                   IF ( mode == 'xy' )  level_z(nzb+1) = zu(nzb+1)
564                ENDIF
565
566                IF ( .NOT. found )  THEN
567                   message_string = 'no output provided for: ' //    &
568                                    TRIM( do2d(av,if) )
569                   CALL message( 'data_output_2d', 'PA0181', 0, 0, 0, 6, 0 )
570                ENDIF
571
572          END SELECT
573
574!
575!--       Resort the array to be output, if not done above
576          IF ( .NOT. resorted )  THEN
577             DO  i = nxl-1, nxr+1
578                DO  j = nys-1, nyn+1
579                   DO  k = nzb, nzt+1
580                      local_pf(i,j,k) = to_be_resorted(k,j,i)
581                   ENDDO
582                ENDDO
583             ENDDO
584          ENDIF
585
586!
587!--       Output of the individual cross-sections, depending on the cross-
588!--       section mode chosen.
589          is = 1
590   loop1: DO  WHILE ( section(is,s) /= -9999  .OR.  two_d )
591
592             SELECT CASE ( mode )
593
594                CASE ( 'xy' )
595!
596!--                Determine the cross section index
597                   IF ( two_d )  THEN
598                      layer_xy = nzb+1
599                   ELSE
600                      layer_xy = section(is,s)
601                   ENDIF
602
603!
604!--                Update the NetCDF xy cross section time axis
605                   IF ( myid == 0 )  THEN
606                      IF ( simulated_time /= do2d_xy_last_time(av) )  THEN
607                         do2d_xy_time_count(av) = do2d_xy_time_count(av) + 1
608                         do2d_xy_last_time(av)  = simulated_time
609                         IF ( .NOT. data_output_2d_on_each_pe  .AND. &
610                              netcdf_output )  THEN
611#if defined( __netcdf )
612                            nc_stat = NF90_PUT_VAR( id_set_xy(av),             &
613                                                    id_var_time_xy(av),        &
614                                             (/ time_since_reference_point /), &
615                                         start = (/ do2d_xy_time_count(av) /), &
616                                                    count = (/ 1 /) )
617                         CALL handle_netcdf_error( 'data_output_2d', 53 )
618#endif
619                         ENDIF
620                      ENDIF
621                   ENDIF
622!
623!--                If required, carry out averaging along z
624                   IF ( section(is,s) == -1  .AND.  .NOT. two_d )  THEN
625
626                      local_2d = 0.0
627!
628!--                   Carry out the averaging (all data are on the PE)
629                      DO  k = nzb, nzt+1
630                         DO  j = nys-1, nyn+1
631                            DO  i = nxl-1, nxr+1
632                               local_2d(i,j) = local_2d(i,j) + local_pf(i,j,k)
633                            ENDDO
634                         ENDDO
635                      ENDDO
636
637                      local_2d = local_2d / ( nzt -nzb + 2.0 )
638
639                   ELSE
640!
641!--                   Just store the respective section on the local array
642                      local_2d = local_pf(:,:,layer_xy)
643
644                   ENDIF
645
646#if defined( __parallel )
647                   IF ( data_output_2d_on_each_pe )  THEN
648!
649!--                   Output of partial arrays on each PE
650#if defined( __netcdf )
651                      IF ( netcdf_output  .AND.  myid == 0 )  THEN
652                         WRITE ( 21 )  simulated_time, do2d_xy_time_count(av), &
653                                       av
654                      ENDIF
655#endif
656                      WRITE ( 21 )  nxl-1, nxr+1, nys-1, nyn+1
657                      WRITE ( 21 )  local_2d
658
659                   ELSE
660!
661!--                   PE0 receives partial arrays from all processors and then
662!--                   outputs them. Here a barrier has to be set, because
663!--                   otherwise "-MPI- FATAL: Remote protocol queue full" may
664!--                   occur.
665                      CALL MPI_BARRIER( comm2d, ierr )
666
667                      ngp = ( nxr-nxl+3 ) * ( nyn-nys+3 )
668                      IF ( myid == 0 )  THEN
669!
670!--                      Local array can be relocated directly.
671                         total_2d(nxl-1:nxr+1,nys-1:nyn+1) = local_2d
672!
673!--                      Receive data from all other PEs.
674                         DO  n = 1, numprocs-1
675!
676!--                         Receive index limits first, then array.
677!--                         Index limits are received in arbitrary order from
678!--                         the PEs.
679                            CALL MPI_RECV( ind(1), 4, MPI_INTEGER,            &
680                                           MPI_ANY_SOURCE, 0, comm2d, status, &
681                                           ierr )
682                            sender = status(MPI_SOURCE)
683                            DEALLOCATE( local_2d )
684                            ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
685                            CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp,      &
686                                           MPI_REAL, sender, 1, comm2d,       &
687                                           status, ierr )
688                            total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
689                         ENDDO
690!
691!--                      Output of the total cross-section.
692                         IF ( iso2d_output ) WRITE (21)  total_2d(0:nx+1,0:ny+1)
693!
694!--                      Relocate the local array for the next loop increment
695                         DEALLOCATE( local_2d )
696                         ALLOCATE( local_2d(nxl-1:nxr+1,nys-1:nyn+1) )
697
698#if defined( __netcdf )
699                         IF ( netcdf_output )  THEN
700                            IF ( two_d ) THEN
701                               nc_stat = NF90_PUT_VAR( id_set_xy(av),          &
702                                                       id_var_do2d(av,if),     &
703                                                      total_2d(0:nx+1,0:ny+1), &
704                                start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
705                                                count = (/ nx+2, ny+2, 1, 1 /) )
706                            ELSE
707                               nc_stat = NF90_PUT_VAR( id_set_xy(av),          &
708                                                       id_var_do2d(av,if),     &
709                                                      total_2d(0:nx+1,0:ny+1), &
710                               start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
711                                                count = (/ nx+2, ny+2, 1, 1 /) )
712                            ENDIF
713                            CALL handle_netcdf_error( 'data_output_2d', 54 )
714                         ENDIF
715#endif
716
717                      ELSE
718!
719!--                      First send the local index limits to PE0
720                         ind(1) = nxl-1; ind(2) = nxr+1
721                         ind(3) = nys-1; ind(4) = nyn+1
722                         CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, comm2d, &
723                                        ierr )
724!
725!--                      Send data to PE0
726                         CALL MPI_SEND( local_2d(nxl-1,nys-1), ngp, MPI_REAL, &
727                                        0, 1, comm2d, ierr )
728                      ENDIF
729!
730!--                   A barrier has to be set, because otherwise some PEs may
731!--                   proceed too fast so that PE0 may receive wrong data on
732!--                   tag 0
733                      CALL MPI_BARRIER( comm2d, ierr )
734                   ENDIF
735#else
736                   IF ( iso2d_output )  THEN
737                      WRITE (21)  local_2d(nxl:nxr+1,nys:nyn+1)
738                   ENDIF
739#if defined( __netcdf )
740                   IF ( netcdf_output )  THEN
741                      IF ( two_d ) THEN
742                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
743                                                 id_var_do2d(av,if),           &
744                                                local_2d(nxl:nxr+1,nys:nyn+1), &
745                                start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
746                                              count = (/ nx+2, ny+2, 1, 1 /) )
747                      ELSE
748                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
749                                                 id_var_do2d(av,if),           &
750                                                local_2d(nxl:nxr+1,nys:nyn+1), &
751                               start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
752                                              count = (/ nx+2, ny+2, 1, 1 /) )
753                      ENDIF
754                      CALL handle_netcdf_error( 'data_output_2d', 55 )
755                   ENDIF
756#endif
757#endif
758                   do2d_xy_n = do2d_xy_n + 1
759!
760!--                Write LOCAL parameter set for ISO2D.
761                   IF ( myid == 0  .AND.  iso2d_output )  THEN
762                      IF ( section(is,s) /= -1 )  THEN
763                         WRITE ( section_chr, '(''z = '',F7.2,'' m  (GP '',I3, &
764                                               &'')'')'                        &
765                               )  level_z(layer_xy), layer_xy
766                      ELSE
767                         section_chr = 'averaged along z'
768                      ENDIF
769                      IF ( av == 0 )  THEN
770                         rtext = TRIM( do2d(av,if) ) // '  t = ' //    &
771                                 TRIM( simulated_time_chr ) // '  ' // &
772                                 TRIM( section_chr )
773                      ELSE
774                         rtext = TRIM( do2d(av,if) ) // '  averaged t = ' // &
775                                 TRIM( simulated_time_chr ) // '  ' //       &
776                                 TRIM( section_chr )
777                      ENDIF
778                      WRITE (27,LOCAL)
779                   ENDIF
780!
781!--                For 2D-arrays (e.g. u*) only one cross-section is available.
782!--                Hence exit loop of output levels.
783                   IF ( two_d )  THEN
784                      two_d = .FALSE.
785                      EXIT loop1
786                   ENDIF
787
788                CASE ( 'xz' )
789!
790!--                Update the NetCDF xz cross section time axis
791                   IF ( myid == 0 )  THEN
792                      IF ( simulated_time /= do2d_xz_last_time(av) )  THEN
793                         do2d_xz_time_count(av) = do2d_xz_time_count(av) + 1
794                         do2d_xz_last_time(av)  = simulated_time
795                         IF ( .NOT. data_output_2d_on_each_pe  .AND. &
796                              netcdf_output )  THEN
797#if defined( __netcdf )
798                            nc_stat = NF90_PUT_VAR( id_set_xz(av),             &
799                                                    id_var_time_xz(av),        &
800                                             (/ time_since_reference_point /), &
801                                         start = (/ do2d_xz_time_count(av) /), &
802                                                    count = (/ 1 /) )
803                             CALL handle_netcdf_error( 'data_output_2d', 56 )
804#endif
805                         ENDIF
806                      ENDIF
807                   ENDIF
808!
809!--                If required, carry out averaging along y
810                   IF ( section(is,s) == -1 )  THEN
811
812                      ALLOCATE( local_2d_l(nxl-1:nxr+1,nzb:nzt+1) )
813                      local_2d_l = 0.0
814                      ngp = ( nxr-nxl+3 ) * ( nzt-nzb+2 )
815!
816!--                   First local averaging on the PE
817                      DO  k = nzb, nzt+1
818                         DO  j = nys, nyn
819                            DO  i = nxl-1, nxr+1
820                               local_2d_l(i,k) = local_2d_l(i,k) + &
821                                                 local_pf(i,j,k)
822                            ENDDO
823                         ENDDO
824                      ENDDO
825#if defined( __parallel )
826!
827!--                   Now do the averaging over all PEs along y
828                      CALL MPI_ALLREDUCE( local_2d_l(nxl-1,nzb),              &
829                                          local_2d(nxl-1,nzb), ngp, MPI_REAL, &
830                                          MPI_SUM, comm1dy, ierr )
831#else
832                      local_2d = local_2d_l
833#endif
834                      local_2d = local_2d / ( ny + 1.0 )
835
836                      DEALLOCATE( local_2d_l )
837
838                   ELSE
839!
840!--                   Just store the respective section on the local array
841!--                   (but only if it is available on this PE!)
842                      IF ( section(is,s) >= nys  .AND.  section(is,s) <= nyn ) &
843                      THEN
844                         local_2d = local_pf(:,section(is,s),nzb:nzt+1)
845                      ENDIF
846
847                   ENDIF
848
849#if defined( __parallel )
850                   IF ( data_output_2d_on_each_pe )  THEN
851!
852!--                   Output of partial arrays on each PE. If the cross section
853!--                   does not reside on the PE, output special index values.
854#if defined( __netcdf )
855                      IF ( netcdf_output  .AND.  myid == 0 )  THEN
856                         WRITE ( 22 )  simulated_time, do2d_xz_time_count(av), &
857                                       av
858                      ENDIF
859#endif
860                      IF ( ( section(is,s)>=nys .AND. section(is,s)<=nyn ) .OR.&
861                           ( section(is,s) == -1  .AND.  nys-1 == -1 ) )       &
862                      THEN
863                         WRITE (22)  nxl-1, nxr+1, nzb, nzt+1
864                         WRITE (22)  local_2d
865                      ELSE
866                         WRITE (22)  -1, -1, -1, -1
867                      ENDIF
868
869                   ELSE
870!
871!--                   PE0 receives partial arrays from all processors of the
872!--                   respective cross section and outputs them. Here a
873!--                   barrier has to be set, because otherwise
874!--                   "-MPI- FATAL: Remote protocol queue full" may occur.
875                      CALL MPI_BARRIER( comm2d, ierr )
876
877                      ngp = ( nxr-nxl+3 ) * ( nzt-nzb+2 )
878                      IF ( myid == 0 )  THEN
879!
880!--                      Local array can be relocated directly.
881                         IF ( ( section(is,s)>=nys .AND. section(is,s)<=nyn )  &
882                            .OR. ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
883                         THEN
884                            total_2d(nxl-1:nxr+1,nzb:nzt+1) = local_2d
885                         ENDIF
886!
887!--                      Receive data from all other PEs.
888                         DO  n = 1, numprocs-1
889!
890!--                         Receive index limits first, then array.
891!--                         Index limits are received in arbitrary order from
892!--                         the PEs.
893                            CALL MPI_RECV( ind(1), 4, MPI_INTEGER,            &
894                                           MPI_ANY_SOURCE, 0, comm2d, status, &
895                                           ierr )
896!
897!--                         Not all PEs have data for XZ-cross-section.
898                            IF ( ind(1) /= -9999 )  THEN
899                               sender = status(MPI_SOURCE)
900                               DEALLOCATE( local_2d )
901                               ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
902                               CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
903                                              MPI_REAL, sender, 1, comm2d,  &
904                                              status, ierr )
905                               total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
906                            ENDIF
907                         ENDDO
908!
909!--                      Output of the total cross-section.
910                         IF ( iso2d_output )  THEN
911                            WRITE (22)  total_2d(0:nx+1,nzb:nzt+1)
912                         ENDIF
913!
914!--                      Relocate the local array for the next loop increment
915                         DEALLOCATE( local_2d )
916                         ALLOCATE( local_2d(nxl-1:nxr+1,nzb:nzt+1) )
917
918#if defined( __netcdf )
919                         IF ( netcdf_output )  THEN
920                            nc_stat = NF90_PUT_VAR( id_set_xz(av),             &
921                                                    id_var_do2d(av,if),        &
922                                                    total_2d(0:nx+1,nzb:nzt+1),&
923                               start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
924                                                count = (/ nx+2, 1, nz+2, 1 /) )
925                            CALL handle_netcdf_error( 'data_output_2d', 57 )
926                         ENDIF
927#endif
928
929                      ELSE
930!
931!--                      If the cross section resides on the PE, send the
932!--                      local index limits, otherwise send -9999 to PE0.
933                         IF ( ( section(is,s)>=nys .AND. section(is,s)<=nyn )  &
934                            .OR. ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
935                         THEN
936                            ind(1) = nxl-1; ind(2) = nxr+1
937                            ind(3) = nzb;   ind(4) = nzt+1
938                         ELSE
939                            ind(1) = -9999; ind(2) = -9999
940                            ind(3) = -9999; ind(4) = -9999
941                         ENDIF
942                         CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, comm2d, &
943                                        ierr )
944!
945!--                      If applicable, send data to PE0.
946                         IF ( ind(1) /= -9999 )  THEN
947                            CALL MPI_SEND( local_2d(nxl-1,nzb), ngp, MPI_REAL, &
948                                           0, 1, comm2d, ierr )
949                         ENDIF
950                      ENDIF
951!
952!--                   A barrier has to be set, because otherwise some PEs may
953!--                   proceed too fast so that PE0 may receive wrong data on
954!--                   tag 0
955                      CALL MPI_BARRIER( comm2d, ierr )
956                   ENDIF
957#else
958                   IF ( iso2d_output )  THEN
959                      WRITE (22)  local_2d(nxl:nxr+1,nzb:nzt+1)
960                   ENDIF
961#if defined( __netcdf )
962                   IF ( netcdf_output )  THEN
963                      nc_stat = NF90_PUT_VAR( id_set_xz(av),                   &
964                                              id_var_do2d(av,if),              &
965                                              local_2d(nxl:nxr+1,nzb:nzt+1),   &
966                               start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
967                                              count = (/ nx+2, 1, nz+2, 1 /) )
968                      CALL handle_netcdf_error( 'data_output_2d', 58 )
969                   ENDIF
970#endif
971#endif
972                   do2d_xz_n = do2d_xz_n + 1
973!
974!--                Write LOCAL-parameter set for ISO2D.
975                   IF ( myid == 0  .AND.  iso2d_output )  THEN
976                      IF ( section(is,s) /= -1 )  THEN
977                         WRITE ( section_chr, '(''y = '',F8.2,'' m  (GP '',I3, &
978                                               &'')'')'                        &
979                               )  section(is,s)*dy, section(is,s)
980                      ELSE
981                         section_chr = 'averaged along y'
982                      ENDIF
983                      IF ( av == 0 )  THEN
984                         rtext = TRIM( do2d(av,if) ) // '  t = ' //    &
985                                 TRIM( simulated_time_chr ) // '  ' // &
986                                 TRIM( section_chr )
987                      ELSE
988                         rtext = TRIM( do2d(av,if) ) // '  averaged t = ' // &
989                                 TRIM( simulated_time_chr ) // '  ' //       &
990                                 TRIM( section_chr )
991                      ENDIF
992                      WRITE (28,LOCAL)
993                   ENDIF
994
995                CASE ( 'yz' )
996!
997!--                Update the NetCDF xy cross section time axis
998                   IF ( myid == 0 )  THEN
999                      IF ( simulated_time /= do2d_yz_last_time(av) )  THEN
1000                         do2d_yz_time_count(av) = do2d_yz_time_count(av) + 1
1001                         do2d_yz_last_time(av)  = simulated_time
1002                         IF ( .NOT. data_output_2d_on_each_pe  .AND. &
1003                              netcdf_output )  THEN
1004#if defined( __netcdf )
1005                            nc_stat = NF90_PUT_VAR( id_set_yz(av),             &
1006                                                    id_var_time_yz(av),        &
1007                                             (/ time_since_reference_point /), &
1008                                         start = (/ do2d_yz_time_count(av) /), &
1009                                                    count = (/ 1 /) )
1010                            CALL handle_netcdf_error( 'data_output_2d', 59 )
1011#endif
1012                         ENDIF
1013                      ENDIF
1014                   ENDIF
1015!
1016!--                If required, carry out averaging along x
1017                   IF ( section(is,s) == -1 )  THEN
1018
1019                      ALLOCATE( local_2d_l(nys-1:nyn+1,nzb:nzt+1) )
1020                      local_2d_l = 0.0
1021                      ngp = ( nyn-nys+3 ) * ( nzt-nzb+2 )
1022!
1023!--                   First local averaging on the PE
1024                      DO  k = nzb, nzt+1
1025                         DO  j = nys-1, nyn+1
1026                            DO  i = nxl, nxr
1027                               local_2d_l(j,k) = local_2d_l(j,k) + &
1028                                                 local_pf(i,j,k)
1029                            ENDDO
1030                         ENDDO
1031                      ENDDO
1032#if defined( __parallel )
1033!
1034!--                   Now do the averaging over all PEs along x
1035                      CALL MPI_ALLREDUCE( local_2d_l(nys-1,nzb),              &
1036                                          local_2d(nys-1,nzb), ngp, MPI_REAL, &
1037                                          MPI_SUM, comm1dx, ierr )
1038#else
1039                      local_2d = local_2d_l
1040#endif
1041                      local_2d = local_2d / ( nx + 1.0 )
1042
1043                      DEALLOCATE( local_2d_l )
1044
1045                   ELSE
1046!
1047!--                   Just store the respective section on the local array
1048!--                   (but only if it is available on this PE!)
1049                      IF ( section(is,s) >= nxl  .AND.  section(is,s) <= nxr ) &
1050                      THEN
1051                         local_2d = local_pf(section(is,s),:,nzb:nzt+1)
1052                      ENDIF
1053
1054                   ENDIF
1055
1056#if defined( __parallel )
1057                   IF ( data_output_2d_on_each_pe )  THEN
1058!
1059!--                   Output of partial arrays on each PE. If the cross section
1060!--                   does not reside on the PE, output special index values.
1061#if defined( __netcdf )
1062                      IF ( netcdf_output  .AND.  myid == 0 )  THEN
1063                         WRITE ( 23 )  simulated_time, do2d_yz_time_count(av), &
1064                                       av
1065                      ENDIF
1066#endif
1067                      IF ( ( section(is,s)>=nxl .AND. section(is,s)<=nxr ) .OR.&
1068                           ( section(is,s) ==  -1  .AND.  nxl-1 == -1 ) )      &
1069                      THEN
1070                         WRITE (23)  nys-1, nyn+1, nzb, nzt+1
1071                         WRITE (23)  local_2d
1072                      ELSE
1073                         WRITE (23)  -1, -1, -1, -1
1074                      ENDIF
1075
1076                   ELSE
1077!
1078!--                   PE0 receives partial arrays from all processors of the
1079!--                   respective cross section and outputs them. Here a
1080!--                   barrier has to be set, because otherwise
1081!--                   "-MPI- FATAL: Remote protocol queue full" may occur.
1082                      CALL MPI_BARRIER( comm2d, ierr )
1083
1084                      ngp = ( nyn-nys+3 ) * ( nzt-nzb+2 )
1085                      IF ( myid == 0 )  THEN
1086!
1087!--                      Local array can be relocated directly.
1088                         IF ( ( section(is,s)>=nxl .AND. section(is,s)<=nxr )  &
1089                           .OR. ( section(is,s) ==  -1  .AND.  nxl-1 == -1 ) ) &
1090                         THEN
1091                            total_2d(nys-1:nyn+1,nzb:nzt+1) = local_2d
1092                         ENDIF
1093!
1094!--                      Receive data from all other PEs.
1095                         DO  n = 1, numprocs-1
1096!
1097!--                         Receive index limits first, then array.
1098!--                         Index limits are received in arbitrary order from
1099!--                         the PEs.
1100                            CALL MPI_RECV( ind(1), 4, MPI_INTEGER,            &
1101                                           MPI_ANY_SOURCE, 0, comm2d, status, &
1102                                           ierr )
1103!
1104!--                         Not all PEs have data for YZ-cross-section.
1105                            IF ( ind(1) /= -9999 )  THEN
1106                               sender = status(MPI_SOURCE)
1107                               DEALLOCATE( local_2d )
1108                               ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
1109                               CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1110                                              MPI_REAL, sender, 1, comm2d,  &
1111                                              status, ierr )
1112                               total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
1113                            ENDIF
1114                         ENDDO
1115!
1116!--                      Output of the total cross-section.
1117                         IF ( iso2d_output )  THEN
1118                            WRITE (23)  total_2d(0:ny+1,nzb:nzt+1)
1119                         ENDIF
1120!
1121!--                      Relocate the local array for the next loop increment
1122                         DEALLOCATE( local_2d )
1123                         ALLOCATE( local_2d(nys-1:nyn+1,nzb:nzt+1) )
1124
1125#if defined( __netcdf )
1126                         IF ( netcdf_output )  THEN
1127                            nc_stat = NF90_PUT_VAR( id_set_yz(av),             &
1128                                                    id_var_do2d(av,if),        &
1129                                                    total_2d(0:ny+1,nzb:nzt+1),&
1130                               start = (/ is, 1, 1, do2d_yz_time_count(av) /), &
1131                                                count = (/ 1, ny+2, nz+2, 1 /) )
1132                            CALL handle_netcdf_error( 'data_output_2d', 60 )
1133                         ENDIF
1134#endif
1135
1136                      ELSE
1137!
1138!--                      If the cross section resides on the PE, send the
1139!--                      local index limits, otherwise send -9999 to PE0.
1140                         IF ( ( section(is,s)>=nxl .AND. section(is,s)<=nxr )  &
1141                           .OR. ( section(is,s) ==  -1  .AND.  nxl-1 == -1 ) ) &
1142                         THEN
1143                            ind(1) = nys-1; ind(2) = nyn+1
1144                            ind(3) = nzb;   ind(4) = nzt+1
1145                         ELSE
1146                            ind(1) = -9999; ind(2) = -9999
1147                            ind(3) = -9999; ind(4) = -9999
1148                         ENDIF
1149                         CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, comm2d, &
1150                                        ierr )
1151!
1152!--                      If applicable, send data to PE0.
1153                         IF ( ind(1) /= -9999 )  THEN
1154                            CALL MPI_SEND( local_2d(nys-1,nzb), ngp, MPI_REAL, &
1155                                           0, 1, comm2d, ierr )
1156                         ENDIF
1157                      ENDIF
1158!
1159!--                   A barrier has to be set, because otherwise some PEs may
1160!--                   proceed too fast so that PE0 may receive wrong data on
1161!--                   tag 0
1162                      CALL MPI_BARRIER( comm2d, ierr )
1163                   ENDIF
1164#else
1165                   IF ( iso2d_output )  THEN
1166                      WRITE (23)  local_2d(nys:nyn+1,nzb:nzt+1)
1167                   ENDIF
1168#if defined( __netcdf )
1169                   IF ( netcdf_output )  THEN
1170                      nc_stat = NF90_PUT_VAR( id_set_yz(av),                   &
1171                                              id_var_do2d(av,if),              &
1172                                              local_2d(nys:nyn+1,nzb:nzt+1),   &
1173                               start = (/ is, 1, 1, do2d_xz_time_count(av) /), &
1174                                              count = (/ 1, ny+2, nz+2, 1 /) )
1175                      CALL handle_netcdf_error( 'data_output_2d', 61 )
1176                   ENDIF
1177#endif
1178#endif
1179                   do2d_yz_n = do2d_yz_n + 1
1180!
1181!--                Write LOCAL-parameter set for ISO2D.
1182                   IF ( myid == 0  .AND.  iso2d_output )  THEN
1183                      IF ( section(is,s) /= -1 )  THEN
1184                         WRITE ( section_chr, '(''x = '',F8.2,'' m  (GP '',I3, &
1185                                               &'')'')'                        &
1186                               )  section(is,s)*dx, section(is,s)
1187                      ELSE
1188                         section_chr = 'averaged along x'
1189                      ENDIF
1190                      IF ( av == 0 )  THEN
1191                         rtext = TRIM( do2d(av,if) ) // '  t = ' //    &
1192                                 TRIM( simulated_time_chr ) // '  ' // &
1193                                 TRIM( section_chr )
1194                      ELSE
1195                         rtext = TRIM( do2d(av,if) ) // '  averaged t = ' // &
1196                                 TRIM( simulated_time_chr ) // '  ' //       &
1197                                 TRIM( section_chr )
1198                      ENDIF
1199                      WRITE (29,LOCAL)
1200                   ENDIF
1201
1202             END SELECT
1203
1204             is = is + 1
1205          ENDDO loop1
1206
1207       ENDIF
1208
1209       if = if + 1
1210       l = MAX( 2, LEN_TRIM( do2d(av,if) ) )
1211       do2d_mode = do2d(av,if)(l-1:l)
1212
1213    ENDDO
1214
1215!
1216!-- Deallocate temporary arrays.
1217    IF ( ALLOCATED( level_z ) )  DEALLOCATE( level_z )
1218    DEALLOCATE( local_pf, local_2d )
1219#if defined( __parallel )
1220    IF ( .NOT.  data_output_2d_on_each_pe  .AND.  myid == 0 )  THEN
1221       DEALLOCATE( total_2d )
1222    ENDIF
1223#endif
1224
1225!
1226!-- Close plot output file.
1227    file_id = 20 + s
1228
1229    IF ( data_output_2d_on_each_pe )  THEN
1230       CALL close_file( file_id )
1231    ELSE
1232       IF ( myid == 0 )  CALL close_file( file_id )
1233    ENDIF
1234
1235
1236    CALL cpu_log (log_point(3),'data_output_2d','stop','nobarrier')
1237
1238 END SUBROUTINE data_output_2d
Note: See TracBrowser for help on using the repository browser.