source: palm/tags/release-3.10/SOURCE/check_parameters.f90 @ 1296

Last change on this file since 1296 was 1277, checked in by heinze, 10 years ago

Last commit documented

  • Property svn:keywords set to Id
File size: 143.9 KB
Line 
1 SUBROUTINE check_parameters
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: check_parameters.f90 1277 2014-01-15 13:43:45Z raasch $
27!
28! 1276 2014-01-15 13:40:41Z heinze
29! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
30!
31! 1241 2013-10-30 11:36:58Z heinze
32! output for profiles of ug and vg added
33! set surface_heatflux and surface_waterflux also in dependence on
34! large_scale_forcing
35! checks for nudging and large scale forcing from external file
36!
37! 1236 2013-09-27 07:21:13Z raasch
38! check number of spectra levels
39!
40! 1216 2013-08-26 09:31:42Z raasch
41! check for transpose_compute_overlap (temporary)
42!
43! 1214 2013-08-21 12:29:17Z kanani
44! additional check for simultaneous use of vertical grid stretching
45! and particle advection
46!
47! 1212 2013-08-15 08:46:27Z raasch
48! checks for poisfft_hybrid removed
49!
50! 1210 2013-08-14 10:58:20Z raasch
51! check for fftw
52!
53! 1179 2013-06-14 05:57:58Z raasch
54! checks and settings of buoyancy parameters and switches revised,
55! initial profile for rho added to hom (id=77)
56!
57! 1174 2013-05-31 10:28:08Z gryschka
58! Bugfix in computing initial profiles for ug, vg, lad, q in case of Atmosphere
59!
60! 1159 2013-05-21 11:58:22Z fricke
61! bc_lr/ns_dirneu/neudir removed
62!
63! 1115 2013-03-26 18:16:16Z hoffmann
64! unused variables removed
65! drizzle can be used without precipitation
66!
67! 1111 2013-03-08 23:54:10Z raasch
68! ibc_p_b = 2 removed
69!
70! 1103 2013-02-20 02:15:53Z raasch
71! Bugfix: turbulent inflow must not require cyclic fill in restart runs
72!
73! 1092 2013-02-02 11:24:22Z raasch
74! unused variables removed
75!
76! 1069 2012-11-28 16:18:43Z maronga
77! allow usage of topography in combination with cloud physics
78!
79! 1065 2012-11-22 17:42:36Z hoffmann
80! Bugfix: It is not allowed to use cloud_scheme = seifert_beheng without
81!         precipitation in order to save computational resources.
82!
83! 1060 2012-11-21 07:19:51Z raasch
84! additional check for parameter turbulent_inflow
85!
86! 1053 2012-11-13 17:11:03Z hoffmann
87! necessary changes for the new two-moment cloud physics scheme added:
88! - check cloud physics scheme (Kessler or Seifert and Beheng)
89! - plant_canopy is not allowed
90! - currently, only cache loop_optimization is allowed
91! - initial profiles of nr, qr
92! - boundary condition of nr, qr
93! - check output quantities (qr, nr, prr)
94!
95! 1036 2012-10-22 13:43:42Z raasch
96! code put under GPL (PALM 3.9)
97!
98! 1031/1034 2012-10-22 11:32:49Z raasch
99! check of netcdf4 parallel file support
100!
101! 1019 2012-09-28 06:46:45Z raasch
102! non-optimized version of prognostic_equations not allowed any more
103!
104! 1015 2012-09-27 09:23:24Z raasch
105! acc allowed for loop optimization,
106! checks for adjustment of mixing length to the Prandtl mixing length removed
107!
108! 1003 2012-09-14 14:35:53Z raasch
109! checks for cases with unequal subdomain sizes removed
110!
111! 1001 2012-09-13 14:08:46Z raasch
112! all actions concerning leapfrog- and upstream-spline-scheme removed
113!
114! 996 2012-09-07 10:41:47Z raasch
115! little reformatting
116
117! 978 2012-08-09 08:28:32Z fricke
118! setting of bc_lr/ns_dirneu/neudir
119! outflow damping layer removed
120! check for z0h*
121! check for pt_damping_width
122!
123! 964 2012-07-26 09:14:24Z raasch
124! check of old profil-parameters removed
125!
126! 940 2012-07-09 14:31:00Z raasch
127! checks for parameter neutral
128!
129! 924 2012-06-06 07:44:41Z maronga
130! Bugfix: preprocessor directives caused error during compilation
131!
132! 892 2012-05-02 13:51:44Z maronga
133! Bugfix for parameter file check ( excluding __netcdf4 )
134!
135! 866 2012-03-28 06:44:41Z raasch
136! use only 60% of the geostrophic wind as translation speed in case of Galilean
137! transformation and use_ug_for_galilei_tr = .T. in order to mimimize the
138! timestep
139!
140! 861 2012-03-26 14:18:34Z suehring
141! Check for topography and ws-scheme removed.
142! Check for loop_optimization = 'vector' and ws-scheme removed.
143!
144! 845 2012-03-07 10:23:05Z maronga
145! Bugfix: exclude __netcdf4 directive part from namelist file check compilation
146!
147! 828 2012-02-21 12:00:36Z raasch
148! check of collision_kernel extended
149!
150! 825 2012-02-19 03:03:44Z raasch
151! check for collision_kernel and curvature_solution_effects
152!
153! 809 2012-01-30 13:32:58Z maronga
154! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
155!
156! 807 2012-01-25 11:53:51Z maronga
157! New cpp directive "__check" implemented which is used by check_namelist_files
158!
159! 774 2011-10-27 13:34:16Z letzel
160! bugfix for prescribed u,v-profiles
161!
162! 767 2011-10-14 06:39:12Z raasch
163! Calculating u,v-profiles from given profiles by linear interpolation.
164! bugfix: dirichlet_0 conditions for ug/vg moved from init_3d_model to here
165!
166! 707 2011-03-29 11:39:40Z raasch
167! setting of bc_lr/ns_dirrad/raddir
168!
169! 689 2011-02-20 19:31:12z gryschka
170! Bugfix for some logical expressions
171! (syntax was not compatible with all compilers)
172!
173! 680 2011-02-04 23:16:06Z gryschka
174! init_vortex is not allowed with volume_flow_control
175!
176! 673 2011-01-18 16:19:48Z suehring
177! Declaration of ws_scheme_sca and ws_scheme_mom added (moved from advec_ws).
178!
179! 667 2010-12-23 12:06:00Z suehring/gryschka
180! Exchange of parameters between ocean and atmosphere via PE0
181! Check for illegal combination of ws-scheme and timestep scheme.
182! Check for topography and ws-scheme.
183! Check for not cyclic boundary conditions in combination with ws-scheme and
184! loop_optimization = 'vector'.
185! Check for call_psolver_at_all_substeps and ws-scheme for momentum_advec.
186! Different processor/grid topology in atmosphere and ocean is now allowed!
187! Bugfixes in checking for conserve_volume_flow_mode
188! 600 2010-11-24 16:10:51Z raasch
189! change due to new default value of surface_waterflux
190! 580 2010-10-05 13:59:11Z heinze
191! renaming of ws_vertical_gradient_level to subs_vertical_gradient_level
192!
193! 567 2010-10-01 10:46:30Z helmke
194! calculating masks changed
195!
196! 564 2010-09-30 13:18:59Z helmke
197! palm message identifiers of masked output changed, 20 replaced by max_masks
198!
199! 553 2010-09-01 14:09:06Z weinreis
200! masks is calculated and removed from inipar
201!
202! 531 2010-04-21 06:47:21Z heinze
203! Bugfix: unit of hyp changed to dbar
204!
205! 524 2010-03-30 02:04:51Z raasch
206! Bugfix: "/" in netcdf profile variable names replaced by ":"
207!
208! 493 2010-03-01 08:30:24Z raasch
209! netcdf_data_format is checked
210!
211! 411 2009-12-11 14:15:58Z heinze
212! Enabled passive scalar/humidity wall fluxes for non-flat topography
213! Initialization of large scale vertical motion (subsidence/ascent)
214!
215! 410 2009-12-04 17:05:40Z letzel
216! masked data output
217!
218! 388 2009-09-23 09:40:33Z raasch
219! Check profiles fpr prho and hyp.
220! Bugfix: output of averaged 2d/3d quantities requires that an avaraging
221! interval has been set, respective error message is included
222! bc_lr_cyc and bc_ns_cyc are set,
223! initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill'
224! Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
225! Coupling with independent precursor runs.
226! Check particle_color, particle_dvrpsize, color_interval, dvrpsize_interval
227! Bugfix: pressure included for profile output
228! Check pressure gradient conditions
229! topography_grid_convention moved from user_check_parameters
230! 'single_street_canyon'
231! Added shf* and qsws* to the list of available output data
232!
233! 222 2009-01-12 16:04:16Z letzel
234! +user_check_parameters
235! Output of messages replaced by message handling routine.
236! Implementation of an MPI-1 coupling: replaced myid with target_id,
237! deleted __mpi2 directives
238! Check that PALM is called with mrun -K parallel for coupling
239!
240! 197 2008-09-16 15:29:03Z raasch
241! Bug fix: Construction of vertical profiles when 10 gradients have been
242! specified in the parameter list (ug, vg, pt, q, sa, lad)
243!   
244! Strict grid matching along z is not needed for mg-solver.
245! Leaf area density (LAD) explicitly set to its surface value at k=0
246! Case of reading data for recycling included in initializing_actions,
247! check of turbulent_inflow and calculation of recycling_plane.
248! q*2 profile added
249!
250! 138 2007-11-28 10:03:58Z letzel
251! Plant canopy added
252! Allow new case bc_uv_t = 'dirichlet_0' for channel flow.
253! Multigrid solver allows topography, checking of dt_sort_particles
254! Bugfix: initializing u_init and v_init in case of ocean runs
255!
256! 109 2007-08-28 15:26:47Z letzel
257! Check coupling_mode and set default (obligatory) values (like boundary
258! conditions for temperature and fluxes) in case of coupled runs.
259! +profiles for w*p* and w"e
260! Bugfix: Error message concerning output of particle concentration (pc)
261! modified
262! More checks and more default values for coupled runs
263! allow data_output_pr= q, wq, w"q", w*q* for humidity = .T. (instead of
264! cloud_physics = .T.)
265! Rayleigh damping for ocean fixed.
266! Check and, if necessary, set default value for dt_coupling
267!
268! 97 2007-06-21 08:23:15Z raasch
269! Initial salinity profile is calculated, salinity boundary conditions are
270! checked,
271! z_max_do1d is checked only in case of ocean = .f.,
272! +initial temperature and geostrophic velocity profiles for the ocean version,
273! use_pt_reference renamed use_reference
274!
275! 89 2007-05-25 12:08:31Z raasch
276! Check for user-defined profiles
277!
278! 75 2007-03-22 09:54:05Z raasch
279! "by_user" allowed as initializing action, -data_output_ts,
280! leapfrog with non-flat topography not allowed any more, loop_optimization
281! and pt_reference are checked, moisture renamed humidity,
282! output of precipitation amount/rate and roughnes length + check
283! possible negative humidities are avoided in initial profile,
284! dirichlet/neumann changed to dirichlet/radiation, etc.,
285! revision added to run_description_header
286!
287! 20 2007-02-26 00:12:32Z raasch
288! Temperature and humidity gradients at top are now calculated for nzt+1,
289! top_heatflux and respective boundary condition bc_pt_t is checked
290!
291! RCS Log replace by Id keyword, revision history cleaned up
292!
293! Revision 1.61  2006/08/04 14:20:25  raasch
294! do2d_unit and do3d_unit now defined as 2d-arrays, check of
295! use_upstream_for_tke, default value for dt_dopts,
296! generation of file header moved from routines palm and header to here
297!
298! Revision 1.1  1997/08/26 06:29:23  raasch
299! Initial revision
300!
301!
302! Description:
303! ------------
304! Check control parameters and deduce further quantities.
305!------------------------------------------------------------------------------!
306
307    USE arrays_3d
308    USE cloud_parameters
309    USE constants
310    USE control_parameters
311    USE dvrp_variables
312    USE grid_variables
313    USE indices
314    USE model_1d
315    USE netcdf_control
316    USE particle_attributes
317    USE pegrid
318    USE profil_parameter
319    USE spectrum
320    USE statistics
321    USE subsidence_mod
322    USE statistics
323    USE transpose_indices
324
325    IMPLICIT NONE
326
327    CHARACTER (LEN=1)   ::  sq
328    CHARACTER (LEN=6)   ::  var
329    CHARACTER (LEN=7)   ::  unit
330    CHARACTER (LEN=8)   ::  date
331    CHARACTER (LEN=10)  ::  time
332    CHARACTER (LEN=40)  ::  coupling_string
333    CHARACTER (LEN=100) ::  action
334
335    INTEGER ::  i, ilen, iremote = 0, j, k, kk, position, prec
336    LOGICAL ::  found, ldum
337    REAL    ::  gradient, remote = 0.0, simulation_time_since_reference
338
339!
340!-- Check for overlap combinations, which are not realized yet
341    IF ( transpose_compute_overlap )  THEN
342       IF ( numprocs == 1 )  STOP '+++ transpose-compute-overlap not implemented for single PE runs'
343#if defined( __openacc )
344       STOP '+++ transpose-compute-overlap not implemented for GPU usage'
345#endif
346    ENDIF
347
348!
349!-- Warning, if host is not set
350    IF ( host(1:1) == ' ' )  THEN
351       message_string = '"host" is not set. Please check that environment ' // &
352                        'variable "localhost" & is set before running PALM'
353       CALL message( 'check_parameters', 'PA0001', 0, 0, 0, 6, 0 )
354    ENDIF
355
356!
357!-- Check the coupling mode
358    IF ( coupling_mode /= 'uncoupled'            .AND.  &
359         coupling_mode /= 'atmosphere_to_ocean'  .AND.  &
360         coupling_mode /= 'ocean_to_atmosphere' )  THEN
361       message_string = 'illegal coupling mode: ' // TRIM( coupling_mode )
362       CALL message( 'check_parameters', 'PA0002', 1, 2, 0, 6, 0 )
363    ENDIF
364
365!
366!-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny
367    IF ( coupling_mode /= 'uncoupled')  THEN
368
369       IF ( dt_coupling == 9999999.9 )  THEN
370          message_string = 'dt_coupling is not set but required for coup' // &
371                           'ling mode "' //  TRIM( coupling_mode ) // '"'
372          CALL message( 'check_parameters', 'PA0003', 1, 2, 0, 6, 0 )
373       ENDIF
374
375#if defined( __parallel )
376
377!
378!--    NOTE: coupled runs have not been implemented in the check_namelist_files
379!--    program.
380!--    check_namelist_files will need the following information of the other
381!--    model (atmosphere/ocean).
382!       dt_coupling = remote
383!       dt_max = remote
384!       restart_time = remote
385!       dt_restart= remote
386!       simulation_time_since_reference = remote
387!       dx = remote
388
389
390#if ! defined( __check )
391       IF ( myid == 0 ) THEN
392          CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, &
393                         ierr )
394          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, &
395                         status, ierr )
396       ENDIF
397       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
398#endif     
399       IF ( dt_coupling /= remote )  THEN
400          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
401                 '": dt_coupling = ', dt_coupling, '& is not equal to ',       &
402                 'dt_coupling_remote = ', remote
403          CALL message( 'check_parameters', 'PA0004', 1, 2, 0, 6, 0 )
404       ENDIF
405       IF ( dt_coupling <= 0.0 )  THEN
406#if ! defined( __check )
407          IF ( myid == 0  ) THEN
408             CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
409             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, &
410                            status, ierr )
411          ENDIF   
412          CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
413#endif         
414          dt_coupling = MAX( dt_max, remote )
415          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
416                 '": dt_coupling <= 0.0 & is not allowed and is reset to ',    &
417                 'MAX(dt_max(A,O)) = ', dt_coupling
418          CALL message( 'check_parameters', 'PA0005', 0, 1, 0, 6, 0 )
419       ENDIF
420#if ! defined( __check )
421       IF ( myid == 0 ) THEN
422          CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &
423                         ierr )
424          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, &
425                         status, ierr )
426       ENDIF
427       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
428#endif     
429       IF ( restart_time /= remote )  THEN
430          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
431                 '": restart_time = ', restart_time, '& is not equal to ',     &
432                 'restart_time_remote = ', remote
433          CALL message( 'check_parameters', 'PA0006', 1, 2, 0, 6, 0 )
434       ENDIF
435#if ! defined( __check )
436       IF ( myid == 0 ) THEN
437          CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, &
438                         ierr )
439          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, &
440                         status, ierr )
441       ENDIF   
442       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
443#endif     
444       IF ( dt_restart /= remote )  THEN
445          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
446                 '": dt_restart = ', dt_restart, '& is not equal to ',         &
447                 'dt_restart_remote = ', remote
448          CALL message( 'check_parameters', 'PA0007', 1, 2, 0, 6, 0 )
449       ENDIF
450
451       simulation_time_since_reference = end_time - coupling_start_time
452#if ! defined( __check )
453       IF  ( myid == 0 ) THEN
454          CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, &
455                         14, comm_inter, ierr )
456          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, &
457                         status, ierr )   
458       ENDIF
459       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
460#endif     
461       IF ( simulation_time_since_reference /= remote )  THEN
462          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
463                 '": simulation_time_since_reference = ',                      &
464                 simulation_time_since_reference, '& is not equal to ',        &
465                 'simulation_time_since_reference_remote = ', remote
466          CALL message( 'check_parameters', 'PA0008', 1, 2, 0, 6, 0 )
467       ENDIF
468
469#if ! defined( __check )
470       IF ( myid == 0 ) THEN
471          CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
472          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, &
473                                                             status, ierr )
474       ENDIF
475       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
476
477#endif
478       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
479
480          IF ( dx < remote ) THEN
481             WRITE( message_string, * ) 'coupling mode "', &
482                   TRIM( coupling_mode ),                  &
483           '": dx in Atmosphere is not equal to or not larger then dx in ocean'
484             CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
485          ENDIF
486
487          IF ( (nx_a+1)*dx /= (nx_o+1)*remote )  THEN
488             WRITE( message_string, * ) 'coupling mode "', &
489                    TRIM( coupling_mode ), &
490             '": Domain size in x-direction is not equal in ocean and atmosphere'
491             CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
492          ENDIF
493
494       ENDIF
495
496#if ! defined( __check )
497       IF ( myid == 0) THEN
498          CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )
499          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, &
500                         status, ierr )
501       ENDIF
502       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
503#endif
504       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
505
506          IF ( dy < remote )  THEN
507             WRITE( message_string, * ) 'coupling mode "', &
508                    TRIM( coupling_mode ), &
509                 '": dy in Atmosphere is not equal to or not larger then dy in ocean'
510             CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
511          ENDIF
512
513          IF ( (ny_a+1)*dy /= (ny_o+1)*remote )  THEN
514             WRITE( message_string, * ) 'coupling mode "', &
515                   TRIM( coupling_mode ), &
516             '": Domain size in y-direction is not equal in ocean and atmosphere'
517             CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )
518          ENDIF
519
520          IF ( MOD(nx_o+1,nx_a+1) /= 0 )  THEN
521             WRITE( message_string, * ) 'coupling mode "', &
522                   TRIM( coupling_mode ), &
523             '": nx+1 in ocean is not divisible without remainder with nx+1 in', & 
524             ' atmosphere'
525             CALL message( 'check_parameters', 'PA0339', 1, 2, 0, 6, 0 )
526          ENDIF
527
528          IF ( MOD(ny_o+1,ny_a+1) /= 0 )  THEN
529             WRITE( message_string, * ) 'coupling mode "', &
530                   TRIM( coupling_mode ), &
531             '": ny+1 in ocean is not divisible without remainder with ny+1 in', & 
532             ' atmosphere'
533             CALL message( 'check_parameters', 'PA0340', 1, 2, 0, 6, 0 )
534          ENDIF
535
536       ENDIF
537#else
538       WRITE( message_string, * ) 'coupling requires PALM to be called with', &
539            ' ''mrun -K parallel'''
540       CALL message( 'check_parameters', 'PA0141', 1, 2, 0, 6, 0 )
541#endif
542    ENDIF
543
544#if defined( __parallel ) && ! defined ( __check )
545!
546!-- Exchange via intercommunicator
547    IF ( coupling_mode == 'atmosphere_to_ocean' .AND. myid == 0 )  THEN
548       CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter, &
549                      ierr )
550    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' .AND. myid == 0)  THEN
551       CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19, &
552                      comm_inter, status, ierr )
553    ENDIF
554    CALL MPI_BCAST( humidity_remote, 1, MPI_LOGICAL, 0, comm2d, ierr)
555   
556#endif
557
558
559!
560!-- Generate the file header which is used as a header for most of PALM's
561!-- output files
562    CALL DATE_AND_TIME( date, time )
563    run_date = date(7:8)//'-'//date(5:6)//'-'//date(3:4)
564    run_time = time(1:2)//':'//time(3:4)//':'//time(5:6)
565    IF ( coupling_mode == 'uncoupled' )  THEN
566       coupling_string = ''
567    ELSEIF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
568       coupling_string = ' coupled (atmosphere)'
569    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
570       coupling_string = ' coupled (ocean)'
571    ENDIF       
572
573    WRITE ( run_description_header,                                        &
574                             '(A,2X,A,2X,A,A,A,I2.2,A,2X,A,A,2X,A,1X,A)' ) &
575              TRIM( version ), TRIM( revision ), 'run: ',                  &
576              TRIM( run_identifier ), '.', runnr, TRIM( coupling_string ), &
577              'host: ', TRIM( host ), run_date, run_time
578
579!
580!-- Check the general loop optimization method
581    IF ( loop_optimization == 'default' )  THEN
582       IF ( host(1:3) == 'nec' )  THEN
583          loop_optimization = 'vector'
584       ELSE
585          loop_optimization = 'cache'
586       ENDIF
587    ENDIF
588
589    SELECT CASE ( TRIM( loop_optimization ) )
590
591       CASE ( 'acc', 'cache', 'vector' )
592          CONTINUE
593
594       CASE DEFAULT
595          message_string = 'illegal value given for loop_optimization: "' // &
596                           TRIM( loop_optimization ) // '"'
597          CALL message( 'check_parameters', 'PA0013', 1, 2, 0, 6, 0 )
598
599    END SELECT
600
601!
602!-- Check if vertical grid stretching is used together with particles
603    IF ( dz_stretch_level < 100000.0 .AND. particle_advection )  THEN
604       message_string = 'Vertical grid stretching is not allowed together ' // &
605                        'with particle advection.'
606       CALL message( 'check_parameters', 'PA0017', 1, 2, 0, 6, 0 )
607    ENDIF
608
609!
610!-- Check topography setting (check for illegal parameter combinations)
611    IF ( topography /= 'flat' )  THEN
612       action = ' '
613       IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme')  THEN
614          WRITE( action, '(A,A)' )  'scalar_advec = ', scalar_advec
615       ENDIF
616       IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' ) &
617       THEN
618          WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
619       ENDIF
620       IF ( psolver == 'sor' )  THEN
621          WRITE( action, '(A,A)' )  'psolver = ', psolver
622       ENDIF
623       IF ( sloping_surface )  THEN
624          WRITE( action, '(A)' )  'sloping surface = .TRUE.'
625       ENDIF
626       IF ( galilei_transformation )  THEN
627          WRITE( action, '(A)' )  'galilei_transformation = .TRUE.'
628       ENDIF
629       IF ( cloud_physics )  THEN
630          WRITE( action, '(A)' )  'cloud_physics = .TRUE.'
631       ENDIF
632       IF ( cloud_droplets )  THEN
633          WRITE( action, '(A)' )  'cloud_droplets = .TRUE.'
634       ENDIF
635       IF ( .NOT. prandtl_layer )  THEN
636          WRITE( action, '(A)' )  'prandtl_layer = .FALSE.'
637       ENDIF
638       IF ( action /= ' ' )  THEN
639          message_string = 'a non-flat topography does not allow ' // &
640                           TRIM( action )
641          CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
642       ENDIF
643!
644!--    In case of non-flat topography, check whether the convention how to
645!--    define the topography grid has been set correctly, or whether the default
646!--    is applicable. If this is not possible, abort.
647       IF ( TRIM( topography_grid_convention ) == ' ' )  THEN
648          IF ( TRIM( topography ) /= 'single_building' .AND.  &
649               TRIM( topography ) /= 'single_street_canyon' .AND.  &
650               TRIM( topography ) /= 'read_from_file' )  THEN
651!--          The default value is not applicable here, because it is only valid
652!--          for the two standard cases 'single_building' and 'read_from_file'
653!--          defined in init_grid.
654             WRITE( message_string, * )  &
655                  'The value for "topography_grid_convention" ',  &
656                  'is not set. Its default value is & only valid for ',  &
657                  '"topography" = ''single_building'', ',  &
658                  '''single_street_canyon'' & or ''read_from_file''.',  &
659                  ' & Choose ''cell_edge'' or ''cell_center''.'
660             CALL message( 'user_check_parameters', 'PA0239', 1, 2, 0, 6, 0 )
661          ELSE
662!--          The default value is applicable here.
663!--          Set convention according to topography.
664             IF ( TRIM( topography ) == 'single_building' .OR.  &
665                  TRIM( topography ) == 'single_street_canyon' )  THEN
666                topography_grid_convention = 'cell_edge'
667             ELSEIF ( TRIM( topography ) == 'read_from_file' )  THEN
668                topography_grid_convention = 'cell_center'
669             ENDIF
670          ENDIF
671       ELSEIF ( TRIM( topography_grid_convention ) /= 'cell_edge' .AND.  &
672                TRIM( topography_grid_convention ) /= 'cell_center' )  THEN
673          WRITE( message_string, * )  &
674               'The value for "topography_grid_convention" is ', &
675               'not recognized. & Choose ''cell_edge'' or ''cell_center''.'
676          CALL message( 'user_check_parameters', 'PA0240', 1, 2, 0, 6, 0 )
677       ENDIF
678
679    ENDIF
680
681!
682!-- Check ocean setting
683    IF ( ocean )  THEN
684
685       action = ' '
686       IF ( action /= ' ' )  THEN
687          message_string = 'ocean = .T. does not allow ' // TRIM( action )
688          CALL message( 'check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
689       ENDIF
690
691    ELSEIF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.  &
692             TRIM( coupling_char ) == '_O' )  THEN
693
694!
695!--    Check whether an (uncoupled) atmospheric run has been declared as an
696!--    ocean run (this setting is done via mrun-option -y)
697
698       message_string = 'ocean = .F. does not allow coupling_char = "' // &
699                        TRIM( coupling_char ) // '" set by mrun-option "-y"'
700       CALL message( 'check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
701
702    ENDIF
703!
704!-- Check cloud scheme
705    IF ( cloud_scheme == 'seifert_beheng' )  THEN
706       icloud_scheme = 0
707    ELSEIF ( cloud_scheme == 'kessler' )  THEN
708       icloud_scheme = 1
709    ELSE
710       message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
711                        TRIM( cloud_scheme ) // '"'
712       CALL message( 'check_parameters', 'PA0357', 1, 2, 0, 6, 0 )
713    ENDIF
714!
715!-- Check whether there are any illegal values
716!-- Pressure solver:
717    IF ( psolver /= 'poisfft'  .AND. &
718         psolver /= 'sor'  .AND.  psolver /= 'multigrid' )  THEN
719       message_string = 'unknown solver for perturbation pressure: psolver' // &
720                        ' = "' // TRIM( psolver ) // '"'
721       CALL message( 'check_parameters', 'PA0016', 1, 2, 0, 6, 0 )
722    ENDIF
723
724    IF ( psolver == 'multigrid' )  THEN
725       IF ( cycle_mg == 'w' )  THEN
726          gamma_mg = 2
727       ELSEIF ( cycle_mg == 'v' )  THEN
728          gamma_mg = 1
729       ELSE
730          message_string = 'unknown multigrid cycle: cycle_mg = "' // &
731                           TRIM( cycle_mg ) // '"'
732          CALL message( 'check_parameters', 'PA0020', 1, 2, 0, 6, 0 )
733       ENDIF
734    ENDIF
735
736    IF ( fft_method /= 'singleton-algorithm'  .AND.  &
737         fft_method /= 'temperton-algorithm'  .AND.  &
738         fft_method /= 'fftw'                 .AND.  &
739         fft_method /= 'system-specific' )  THEN
740       message_string = 'unknown fft-algorithm: fft_method = "' // &
741                        TRIM( fft_method ) // '"'
742       CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
743    ENDIF
744   
745    IF( momentum_advec == 'ws-scheme' .AND. & 
746        .NOT. call_psolver_at_all_substeps  ) THEN
747        message_string = 'psolver must be called at each RK3 substep when "'//&
748                      TRIM(momentum_advec) // ' "is used for momentum_advec'
749        CALL message( 'check_parameters', 'PA0344', 1, 2, 0, 6, 0 )
750    END IF
751!
752!-- Advection schemes:
753    IF ( momentum_advec /= 'pw-scheme'  .AND.  momentum_advec /= 'ws-scheme' ) &
754    THEN
755       message_string = 'unknown advection scheme: momentum_advec = "' // &
756                        TRIM( momentum_advec ) // '"'
757       CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
758    ENDIF
759    IF ( ( momentum_advec == 'ws-scheme' .OR.  scalar_advec == 'ws-scheme' )   &
760           .AND. ( timestep_scheme == 'euler' .OR.                             &
761                   timestep_scheme == 'runge-kutta-2' ) )                      &
762    THEN
763       message_string = 'momentum_advec or scalar_advec = "' &
764         // TRIM( momentum_advec ) // '" is not allowed with timestep_scheme = "' // &
765         TRIM( timestep_scheme ) // '"'
766       CALL message( 'check_parameters', 'PA0023', 1, 2, 0, 6, 0 )
767    ENDIF
768    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' .AND. &
769         scalar_advec /= 'bc-scheme' )                                         &
770    THEN
771       message_string = 'unknown advection scheme: scalar_advec = "' // &
772                        TRIM( scalar_advec ) // '"'
773       CALL message( 'check_parameters', 'PA0024', 1, 2, 0, 6, 0 )
774    ENDIF
775    IF ( scalar_advec == 'bc-scheme'  .AND.  loop_optimization == 'cache' ) &
776    THEN
777       message_string = 'advection_scheme scalar_advec = "' &
778         // TRIM( scalar_advec ) // '" not implemented for & loop_optimization = "' // &
779         TRIM( loop_optimization ) // '"'
780       CALL message( 'check_parameters', 'PA0026', 1, 2, 0, 6, 0 )
781    ENDIF
782
783    IF ( use_sgs_for_particles  .AND.  .NOT. use_upstream_for_tke )  THEN
784       use_upstream_for_tke = .TRUE.
785       message_string = 'use_upstream_for_tke set .TRUE. because ' // &
786                        'use_sgs_for_particles = .TRUE.'
787       CALL message( 'check_parameters', 'PA0025', 0, 1, 0, 6, 0 )
788    ENDIF
789
790    IF ( use_sgs_for_particles  .AND.  curvature_solution_effects )  THEN
791       message_string = 'use_sgs_for_particles = .TRUE. not allowed with ' // &
792                        'curvature_solution_effects = .TRUE.'
793       CALL message( 'check_parameters', 'PA0349', 1, 2, 0, 6, 0 )
794    ENDIF
795
796!
797!-- Set LOGICAL switches to enhance performance
798    IF ( momentum_advec == 'ws-scheme' )    ws_scheme_mom = .TRUE.
799    IF ( scalar_advec   == 'ws-scheme'   )  ws_scheme_sca = .TRUE.
800
801!
802!-- Timestep schemes:
803    SELECT CASE ( TRIM( timestep_scheme ) )
804
805       CASE ( 'euler' )
806          intermediate_timestep_count_max = 1
807
808       CASE ( 'runge-kutta-2' )
809          intermediate_timestep_count_max = 2
810
811       CASE ( 'runge-kutta-3' )
812          intermediate_timestep_count_max = 3
813
814       CASE DEFAULT
815          message_string = 'unknown timestep scheme: timestep_scheme = "' // &
816                           TRIM( timestep_scheme ) // '"'
817          CALL message( 'check_parameters', 'PA0027', 1, 2, 0, 6, 0 )
818
819    END SELECT
820
821    IF ( (momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme') &
822         .AND. timestep_scheme(1:5) == 'runge' ) THEN
823       message_string = 'momentum advection scheme "' // &
824                        TRIM( momentum_advec ) // '" & does not work with ' // &
825                        'timestep_scheme "' // TRIM( timestep_scheme ) // '"'
826       CALL message( 'check_parameters', 'PA0029', 1, 2, 0, 6, 0 )
827    ENDIF
828
829!
830!-- Collision kernels:
831    SELECT CASE ( TRIM( collision_kernel ) )
832
833       CASE ( 'hall', 'hall_fast' )
834          hall_kernel = .TRUE.
835
836       CASE ( 'palm' )
837          palm_kernel = .TRUE.
838
839       CASE ( 'wang', 'wang_fast' )
840          wang_kernel = .TRUE.
841
842       CASE ( 'none' )
843
844
845       CASE DEFAULT
846          message_string = 'unknown collision kernel: collision_kernel = "' // &
847                           TRIM( collision_kernel ) // '"'
848          CALL message( 'check_parameters', 'PA0350', 1, 2, 0, 6, 0 )
849
850    END SELECT
851    IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
852
853    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.  &
854         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
855!
856!--    No restart run: several initialising actions are possible
857       action = initializing_actions
858       DO WHILE ( TRIM( action ) /= '' )
859          position = INDEX( action, ' ' )
860          SELECT CASE ( action(1:position-1) )
861
862             CASE ( 'set_constant_profiles', 'set_1d-model_profiles', &
863                    'by_user', 'initialize_vortex',     'initialize_ptanom' )
864                action = action(position+1:)
865
866             CASE DEFAULT
867                message_string = 'initializing_action = "' // &
868                                 TRIM( action ) // '" unkown or not allowed'
869                CALL message( 'check_parameters', 'PA0030', 1, 2, 0, 6, 0 )
870
871          END SELECT
872       ENDDO
873    ENDIF
874
875    IF ( TRIM( initializing_actions ) == 'initialize_vortex' .AND. &
876         conserve_volume_flow ) THEN
877         message_string = 'initializing_actions = "initialize_vortex"' // &
878                        ' ist not allowed with conserve_volume_flow = .T.'
879       CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
880    ENDIF       
881
882
883    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND. &
884         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
885       message_string = 'initializing_actions = "set_constant_profiles"' // &
886                        ' and "set_1d-model_profiles" are not allowed ' //  &
887                        'simultaneously'
888       CALL message( 'check_parameters', 'PA0031', 1, 2, 0, 6, 0 )
889    ENDIF
890
891    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND. &
892         INDEX( initializing_actions, 'by_user' ) /= 0 )  THEN
893       message_string = 'initializing_actions = "set_constant_profiles"' // &
894                        ' and "by_user" are not allowed simultaneously'
895       CALL message( 'check_parameters', 'PA0032', 1, 2, 0, 6, 0 )
896    ENDIF
897
898    IF ( INDEX( initializing_actions, 'by_user' ) /= 0  .AND. &
899         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
900       message_string = 'initializing_actions = "by_user" and ' // &
901                        '"set_1d-model_profiles" are not allowed simultaneously'
902       CALL message( 'check_parameters', 'PA0033', 1, 2, 0, 6, 0 )
903    ENDIF
904
905    IF ( cloud_physics  .AND.  .NOT. humidity )  THEN
906       WRITE( message_string, * ) 'cloud_physics = ', cloud_physics, ' is ', &
907              'not allowed with humidity = ', humidity
908       CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 )
909    ENDIF
910
911    IF ( precipitation  .AND.  .NOT.  cloud_physics )  THEN
912       WRITE( message_string, * ) 'precipitation = ', precipitation, ' is ', &
913              'not allowed with cloud_physics = ', cloud_physics
914       CALL message( 'check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
915    ENDIF
916
917    IF ( humidity  .AND.  sloping_surface )  THEN
918       message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' // &
919                        'are not allowed simultaneously'
920       CALL message( 'check_parameters', 'PA0036', 1, 2, 0, 6, 0 )
921    ENDIF
922
923    IF ( passive_scalar  .AND.  humidity )  THEN
924       message_string = 'humidity = .TRUE. and passive_scalar = .TRUE. ' // &
925                        'is not allowed simultaneously'
926       CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )
927    ENDIF
928
929    IF ( plant_canopy .AND. ( drag_coefficient == 0.0 ) ) THEN
930       message_string = 'plant_canopy = .TRUE. requires a non-zero drag ' // &
931                        'coefficient & given value is drag_coefficient = 0.0'
932       CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
933    ENDIF
934
935    IF ( plant_canopy  .AND.  cloud_physics  .AND.  icloud_scheme == 0 ) THEN
936       message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' //  &
937                        ' seifert_beheng'
938       CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
939    ENDIF
940
941    IF ( loop_optimization /= 'cache' .AND.  cloud_physics  .AND.            &
942         icloud_scheme == 0 ) THEN
943       message_string = 'cloud_scheme = seifert_beheng requires ' // &
944                        'loop_optimization = cache'
945       CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
946    ENDIF 
947
948!    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
949!         .NOT. precipitation  .AND.  .NOT. drizzle ) THEN
950!       message_string = 'cloud_scheme = seifert_beheng requires ' // &
951!                        'precipitation = .TRUE. or drizzle = .TRUE.'
952!       CALL message( 'check_parameters', 'PA0363', 1, 2, 0, 6, 0 )
953!    ENDIF
954
955!
956!-- In case of no model continuation run, check initialising parameters and
957!-- deduce further quantities
958    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
959
960!
961!--    Initial profiles for 1D and 3D model, respectively (u,v further below)
962       pt_init = pt_surface
963       IF ( humidity )  THEN
964          q_init  = q_surface
965       ENDIF
966       IF ( ocean )           sa_init = sa_surface
967       IF ( passive_scalar )  q_init  = s_surface
968       IF ( plant_canopy )    lad = 0.0
969
970!
971!--
972!--    If required, compute initial profile of the geostrophic wind
973!--    (component ug)
974       i = 1
975       gradient = 0.0
976
977       IF ( .NOT. ocean )  THEN
978
979          ug_vertical_gradient_level_ind(1) = 0
980          ug(0) = ug_surface
981          DO  k = 1, nzt+1
982             IF ( i < 11 ) THEN
983                IF ( ug_vertical_gradient_level(i) < zu(k)  .AND. &
984                     ug_vertical_gradient_level(i) >= 0.0 )  THEN
985                   gradient = ug_vertical_gradient(i) / 100.0
986                   ug_vertical_gradient_level_ind(i) = k - 1
987                   i = i + 1
988                ENDIF
989             ENDIF       
990             IF ( gradient /= 0.0 )  THEN
991                IF ( k /= 1 )  THEN
992                   ug(k) = ug(k-1) + dzu(k) * gradient
993                ELSE
994                   ug(k) = ug_surface + dzu(k) * gradient
995                ENDIF
996             ELSE
997                ug(k) = ug(k-1)
998             ENDIF
999          ENDDO
1000
1001       ELSE
1002
1003          ug_vertical_gradient_level_ind(1) = nzt+1
1004          ug(nzt+1) = ug_surface
1005          DO  k = nzt, nzb, -1
1006             IF ( i < 11 ) THEN
1007                IF ( ug_vertical_gradient_level(i) > zu(k)  .AND. &
1008                     ug_vertical_gradient_level(i) <= 0.0 )  THEN
1009                   gradient = ug_vertical_gradient(i) / 100.0
1010                   ug_vertical_gradient_level_ind(i) = k + 1
1011                   i = i + 1
1012                ENDIF
1013             ENDIF
1014             IF ( gradient /= 0.0 )  THEN
1015                IF ( k /= nzt )  THEN
1016                   ug(k) = ug(k+1) - dzu(k+1) * gradient
1017                ELSE
1018                   ug(k)   = ug_surface - 0.5 * dzu(k+1) * gradient
1019                   ug(k+1) = ug_surface + 0.5 * dzu(k+1) * gradient
1020                ENDIF
1021             ELSE
1022                ug(k) = ug(k+1)
1023             ENDIF
1024          ENDDO
1025
1026       ENDIF
1027
1028!
1029!--    In case of no given gradients for ug, choose a zero gradient
1030       IF ( ug_vertical_gradient_level(1) == -9999999.9 )  THEN
1031          ug_vertical_gradient_level(1) = 0.0
1032       ENDIF 
1033
1034!
1035!--
1036!--    If required, compute initial profile of the geostrophic wind
1037!--    (component vg)
1038       i = 1
1039       gradient = 0.0
1040
1041       IF ( .NOT. ocean )  THEN
1042
1043          vg_vertical_gradient_level_ind(1) = 0
1044          vg(0) = vg_surface
1045          DO  k = 1, nzt+1
1046             IF ( i < 11 ) THEN
1047                IF ( vg_vertical_gradient_level(i) < zu(k)  .AND. &
1048                     vg_vertical_gradient_level(i) >= 0.0 )  THEN
1049                   gradient = vg_vertical_gradient(i) / 100.0
1050                   vg_vertical_gradient_level_ind(i) = k - 1
1051                   i = i + 1
1052                ENDIF
1053             ENDIF
1054             IF ( gradient /= 0.0 )  THEN
1055                IF ( k /= 1 )  THEN
1056                   vg(k) = vg(k-1) + dzu(k) * gradient
1057                ELSE
1058                   vg(k) = vg_surface + dzu(k) * gradient
1059                ENDIF
1060             ELSE
1061                vg(k) = vg(k-1)
1062             ENDIF
1063          ENDDO
1064
1065       ELSE
1066
1067          vg_vertical_gradient_level_ind(1) = nzt+1
1068          vg(nzt+1) = vg_surface
1069          DO  k = nzt, nzb, -1
1070             IF ( i < 11 ) THEN
1071                IF ( vg_vertical_gradient_level(i) > zu(k)  .AND. &
1072                     vg_vertical_gradient_level(i) <= 0.0 )  THEN
1073                   gradient = vg_vertical_gradient(i) / 100.0
1074                   vg_vertical_gradient_level_ind(i) = k + 1
1075                   i = i + 1
1076                ENDIF
1077             ENDIF
1078             IF ( gradient /= 0.0 )  THEN
1079                IF ( k /= nzt )  THEN
1080                   vg(k) = vg(k+1) - dzu(k+1) * gradient
1081                ELSE
1082                   vg(k)   = vg_surface - 0.5 * dzu(k+1) * gradient
1083                   vg(k+1) = vg_surface + 0.5 * dzu(k+1) * gradient
1084                ENDIF
1085             ELSE
1086                vg(k) = vg(k+1)
1087             ENDIF
1088          ENDDO
1089
1090       ENDIF
1091
1092!
1093!--    In case of no given gradients for vg, choose a zero gradient
1094       IF ( vg_vertical_gradient_level(1) == -9999999.9 )  THEN
1095          vg_vertical_gradient_level(1) = 0.0
1096       ENDIF
1097
1098!
1099!--    Let the initial wind profiles be the calculated ug/vg profiles or
1100!--    interpolate them from wind profile data (if given)
1101       IF ( u_profile(1) == 9999999.9  .AND.  v_profile(1) == 9999999.9 )  THEN
1102
1103          u_init = ug
1104          v_init = vg
1105
1106       ELSEIF ( u_profile(1) == 0.0  .AND.  v_profile(1) == 0.0 )  THEN
1107
1108          IF ( uv_heights(1) /= 0.0 )  THEN
1109             message_string = 'uv_heights(1) must be 0.0'
1110             CALL message( 'check_parameters', 'PA0345', 1, 2, 0, 6, 0 )
1111          ENDIF
1112
1113          use_prescribed_profile_data = .TRUE.
1114
1115          kk = 1
1116          u_init(0) = 0.0
1117          v_init(0) = 0.0
1118
1119          DO  k = 1, nz+1
1120
1121             IF ( kk < 100 )  THEN
1122                DO WHILE ( uv_heights(kk+1) <= zu(k) )
1123                   kk = kk + 1
1124                   IF ( kk == 100 )  EXIT
1125                ENDDO
1126             ENDIF
1127
1128             IF ( kk < 100 .AND. uv_heights(kk+1) /= 9999999.9 )  THEN
1129                u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
1130                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
1131                                       ( u_profile(kk+1) - u_profile(kk) )
1132                v_init(k) = v_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
1133                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
1134                                       ( v_profile(kk+1) - v_profile(kk) )
1135             ELSE
1136                u_init(k) = u_profile(kk)
1137                v_init(k) = v_profile(kk)
1138             ENDIF
1139
1140          ENDDO
1141
1142       ELSE
1143
1144          message_string = 'u_profile(1) and v_profile(1) must be 0.0'
1145          CALL message( 'check_parameters', 'PA0346', 1, 2, 0, 6, 0 )
1146
1147       ENDIF
1148
1149!
1150!--    Compute initial temperature profile using the given temperature gradients
1151       IF ( .NOT. neutral )  THEN
1152
1153          i = 1
1154          gradient = 0.0
1155
1156          IF ( .NOT. ocean )  THEN
1157
1158             pt_vertical_gradient_level_ind(1) = 0
1159             DO  k = 1, nzt+1
1160                IF ( i < 11 ) THEN
1161                   IF ( pt_vertical_gradient_level(i) < zu(k)  .AND. &
1162                        pt_vertical_gradient_level(i) >= 0.0 )  THEN
1163                      gradient = pt_vertical_gradient(i) / 100.0
1164                      pt_vertical_gradient_level_ind(i) = k - 1
1165                      i = i + 1
1166                   ENDIF
1167                ENDIF
1168                IF ( gradient /= 0.0 )  THEN
1169                   IF ( k /= 1 )  THEN
1170                      pt_init(k) = pt_init(k-1) + dzu(k) * gradient
1171                   ELSE
1172                      pt_init(k) = pt_surface   + dzu(k) * gradient
1173                   ENDIF
1174                ELSE
1175                   pt_init(k) = pt_init(k-1)
1176                ENDIF
1177             ENDDO
1178
1179          ELSE
1180
1181             pt_vertical_gradient_level_ind(1) = nzt+1
1182             DO  k = nzt, 0, -1
1183                IF ( i < 11 ) THEN
1184                   IF ( pt_vertical_gradient_level(i) > zu(k)  .AND. &
1185                        pt_vertical_gradient_level(i) <= 0.0 )  THEN
1186                      gradient = pt_vertical_gradient(i) / 100.0
1187                      pt_vertical_gradient_level_ind(i) = k + 1
1188                      i = i + 1
1189                   ENDIF
1190                ENDIF
1191                IF ( gradient /= 0.0 )  THEN
1192                   IF ( k /= nzt )  THEN
1193                      pt_init(k) = pt_init(k+1) - dzu(k+1) * gradient
1194                   ELSE
1195                      pt_init(k)   = pt_surface - 0.5 * dzu(k+1) * gradient
1196                      pt_init(k+1) = pt_surface + 0.5 * dzu(k+1) * gradient
1197                   ENDIF
1198                ELSE
1199                   pt_init(k) = pt_init(k+1)
1200                ENDIF
1201             ENDDO
1202
1203          ENDIF
1204
1205       ENDIF
1206
1207!
1208!--    In case of no given temperature gradients, choose gradient of neutral
1209!--    stratification
1210       IF ( pt_vertical_gradient_level(1) == -9999999.9 )  THEN
1211          pt_vertical_gradient_level(1) = 0.0
1212       ENDIF
1213
1214!
1215!--    Store temperature gradient at the top boundary for possible Neumann
1216!--    boundary condition
1217       bc_pt_t_val = ( pt_init(nzt+1) - pt_init(nzt) ) / dzu(nzt+1)
1218
1219!
1220!--    If required, compute initial humidity or scalar profile using the given
1221!--    humidity/scalar gradient. In case of scalar transport, initially store
1222!--    values of the scalar parameters on humidity parameters
1223       IF ( passive_scalar )  THEN
1224          bc_q_b                    = bc_s_b
1225          bc_q_t                    = bc_s_t
1226          q_surface                 = s_surface
1227          q_surface_initial_change  = s_surface_initial_change
1228          q_vertical_gradient       = s_vertical_gradient
1229          q_vertical_gradient_level = s_vertical_gradient_level
1230          surface_waterflux         = surface_scalarflux
1231          wall_humidityflux         = wall_scalarflux
1232       ENDIF
1233
1234       IF ( humidity  .OR.  passive_scalar )  THEN
1235
1236          i = 1
1237          gradient = 0.0
1238          q_vertical_gradient_level_ind(1) = 0
1239          DO  k = 1, nzt+1
1240             IF ( i < 11 ) THEN
1241                IF ( q_vertical_gradient_level(i) < zu(k)  .AND. &
1242                     q_vertical_gradient_level(i) >= 0.0 )  THEN
1243                   gradient = q_vertical_gradient(i) / 100.0
1244                   q_vertical_gradient_level_ind(i) = k - 1
1245                   i = i + 1
1246                ENDIF
1247             ENDIF
1248             IF ( gradient /= 0.0 )  THEN
1249                IF ( k /= 1 )  THEN
1250                   q_init(k) = q_init(k-1) + dzu(k) * gradient
1251                ELSE
1252                   q_init(k) = q_init(k-1) + dzu(k) * gradient
1253                ENDIF
1254             ELSE
1255                q_init(k) = q_init(k-1)
1256             ENDIF
1257!
1258!--          Avoid negative humidities
1259             IF ( q_init(k) < 0.0 )  THEN
1260                q_init(k) = 0.0
1261             ENDIF
1262          ENDDO
1263
1264!
1265!--       In case of no given humidity gradients, choose zero gradient
1266!--       conditions
1267          IF ( q_vertical_gradient_level(1) == -1.0 )  THEN
1268             q_vertical_gradient_level(1) = 0.0
1269          ENDIF
1270!
1271!--       Store humidity, rain water content and rain drop concentration
1272!--       gradient at the top boundary for possile Neumann boundary condition
1273          bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
1274       ENDIF
1275
1276!
1277!--    If required, compute initial salinity profile using the given salinity
1278!--    gradients
1279       IF ( ocean )  THEN
1280
1281          i = 1
1282          gradient = 0.0
1283
1284          sa_vertical_gradient_level_ind(1) = nzt+1
1285          DO  k = nzt, 0, -1
1286             IF ( i < 11 ) THEN
1287                IF ( sa_vertical_gradient_level(i) > zu(k)  .AND. &
1288                     sa_vertical_gradient_level(i) <= 0.0 )  THEN
1289                   gradient = sa_vertical_gradient(i) / 100.0
1290                   sa_vertical_gradient_level_ind(i) = k + 1
1291                   i = i + 1
1292                ENDIF
1293             ENDIF
1294             IF ( gradient /= 0.0 )  THEN
1295                IF ( k /= nzt )  THEN
1296                   sa_init(k) = sa_init(k+1) - dzu(k+1) * gradient
1297                ELSE
1298                   sa_init(k)   = sa_surface - 0.5 * dzu(k+1) * gradient
1299                   sa_init(k+1) = sa_surface + 0.5 * dzu(k+1) * gradient
1300                ENDIF
1301             ELSE
1302                sa_init(k) = sa_init(k+1)
1303             ENDIF
1304          ENDDO
1305
1306       ENDIF
1307
1308!
1309!--    If required compute the profile of leaf area density used in the plant
1310!--    canopy model
1311       IF ( plant_canopy ) THEN
1312       
1313          i = 1
1314          gradient = 0.0
1315
1316          IF ( .NOT. ocean ) THEN
1317
1318             lad(0) = lad_surface
1319 
1320             lad_vertical_gradient_level_ind(1) = 0
1321             DO k = 1, pch_index
1322                IF ( i < 11 ) THEN
1323                   IF ( lad_vertical_gradient_level(i) < zu(k) .AND.  &
1324                        lad_vertical_gradient_level(i) >= 0.0 ) THEN
1325                      gradient = lad_vertical_gradient(i)
1326                      lad_vertical_gradient_level_ind(i) = k - 1
1327                      i = i + 1
1328                   ENDIF
1329                ENDIF
1330                IF ( gradient /= 0.0 ) THEN
1331                   IF ( k /= 1 ) THEN
1332                      lad(k) = lad(k-1) + dzu(k) * gradient
1333                   ELSE
1334                      lad(k) = lad_surface + dzu(k) *gradient
1335                   ENDIF
1336                ELSE
1337                   lad(k) = lad(k-1)
1338                ENDIF
1339             ENDDO
1340
1341          ENDIF
1342
1343!
1344!--       In case of no given leaf area density gradients, choose a vanishing
1345!--       gradient
1346          IF ( lad_vertical_gradient_level(1) == -9999999.9 ) THEN
1347             lad_vertical_gradient_level(1) = 0.0
1348          ENDIF
1349
1350       ENDIF
1351         
1352    ENDIF
1353
1354!
1355!-- Initialize large scale subsidence if required
1356    IF ( subs_vertical_gradient_level(1) /= -9999999.9 )  THEN
1357       large_scale_subsidence = .TRUE.
1358       CALL init_w_subsidence
1359    ENDIF
1360
1361!
1362!-- Compute Coriolis parameter
1363    f  = 2.0 * omega * SIN( phi / 180.0 * pi )
1364    fs = 2.0 * omega * COS( phi / 180.0 * pi )
1365
1366!
1367!-- Check and set buoyancy related parameters and switches
1368    IF ( reference_state == 'horizontal_average' )  THEN
1369       CONTINUE
1370    ELSEIF ( reference_state == 'initial_profile' )  THEN
1371       use_initial_profile_as_reference = .TRUE.
1372    ELSEIF ( reference_state == 'single_value' )  THEN
1373       use_single_reference_value = .TRUE.
1374       IF ( pt_reference == 9999999.9 )  pt_reference = pt_surface
1375       vpt_reference = pt_reference * ( 1.0 + 0.61 * q_surface )
1376    ELSE
1377       message_string = 'illegal value for reference_state: "' // &
1378                        TRIM( reference_state ) // '"'
1379       CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 )
1380    ENDIF
1381
1382!
1383!-- Ocean runs always use reference values in the buoyancy term
1384    IF ( ocean )  THEN
1385       reference_state = 'single_value'
1386       use_single_reference_value = .TRUE.
1387    ENDIF
1388
1389!
1390!-- Sign of buoyancy/stability terms
1391    IF ( ocean )  atmos_ocean_sign = -1.0
1392
1393!
1394!-- Ocean version must use flux boundary conditions at the top
1395    IF ( ocean .AND. .NOT. use_top_fluxes )  THEN
1396       message_string = 'use_top_fluxes must be .TRUE. in ocean version'
1397       CALL message( 'check_parameters', 'PA0042', 1, 2, 0, 6, 0 )
1398    ENDIF
1399
1400!
1401!-- In case of a given slope, compute the relevant quantities
1402    IF ( alpha_surface /= 0.0 )  THEN
1403       IF ( ABS( alpha_surface ) > 90.0 )  THEN
1404          WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface, &
1405                                     ' ) must be < 90.0'
1406          CALL message( 'check_parameters', 'PA0043', 1, 2, 0, 6, 0 )
1407       ENDIF
1408       sloping_surface = .TRUE.
1409       cos_alpha_surface = COS( alpha_surface / 180.0 * pi )
1410       sin_alpha_surface = SIN( alpha_surface / 180.0 * pi )
1411    ENDIF
1412
1413!
1414!-- Check time step and cfl_factor
1415    IF ( dt /= -1.0 )  THEN
1416       IF ( dt <= 0.0  .AND.  dt /= -1.0 )  THEN
1417          WRITE( message_string, * ) 'dt = ', dt , ' <= 0.0'
1418          CALL message( 'check_parameters', 'PA0044', 1, 2, 0, 6, 0 )
1419       ENDIF
1420       dt_3d = dt
1421       dt_fixed = .TRUE.
1422    ENDIF
1423
1424    IF ( cfl_factor <= 0.0  .OR.  cfl_factor > 1.0 )  THEN
1425       IF ( cfl_factor == -1.0 )  THEN
1426          IF ( timestep_scheme == 'runge-kutta-2' )  THEN
1427             cfl_factor = 0.8
1428          ELSEIF ( timestep_scheme == 'runge-kutta-3' )  THEN
1429             cfl_factor = 0.9
1430          ELSE
1431             cfl_factor = 0.9
1432          ENDIF
1433       ELSE
1434          WRITE( message_string, * ) 'cfl_factor = ', cfl_factor, &
1435                 ' out of range & 0.0 < cfl_factor <= 1.0 is required'
1436          CALL message( 'check_parameters', 'PA0045', 1, 2, 0, 6, 0 )
1437       ENDIF
1438    ENDIF
1439
1440!
1441!-- Store simulated time at begin
1442    simulated_time_at_begin = simulated_time
1443
1444!
1445!-- Store reference time for coupled runs and change the coupling flag,
1446!-- if ...
1447    IF ( simulated_time == 0.0 )  THEN
1448       IF ( coupling_start_time == 0.0 )  THEN
1449          time_since_reference_point = 0.0
1450       ELSEIF ( time_since_reference_point < 0.0 )  THEN
1451          run_coupled = .FALSE.
1452       ENDIF
1453    ENDIF
1454
1455!
1456!-- Set wind speed in the Galilei-transformed system
1457    IF ( galilei_transformation )  THEN
1458       IF ( use_ug_for_galilei_tr .AND.                  &
1459            ug_vertical_gradient_level(1) == 0.0  .AND.  &
1460            ug_vertical_gradient(1) == 0.0  .AND.        & 
1461            vg_vertical_gradient_level(1) == 0.0  .AND.  &
1462            vg_vertical_gradient(1) == 0.0 )  THEN
1463          u_gtrans = ug_surface * 0.6
1464          v_gtrans = vg_surface * 0.6
1465       ELSEIF ( use_ug_for_galilei_tr  .AND.                  &
1466                ( ug_vertical_gradient_level(1) /= 0.0  .OR.  &
1467                ug_vertical_gradient(1) /= 0.0 ) )  THEN
1468          message_string = 'baroclinicity (ug) not allowed simultaneously' // &
1469                           ' with galilei transformation'
1470          CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 )
1471       ELSEIF ( use_ug_for_galilei_tr  .AND.                  &
1472                ( vg_vertical_gradient_level(1) /= 0.0  .OR.  &
1473                vg_vertical_gradient(1) /= 0.0 ) )  THEN
1474          message_string = 'baroclinicity (vg) not allowed simultaneously' // &
1475                           ' with galilei transformation'
1476          CALL message( 'check_parameters', 'PA0047', 1, 2, 0, 6, 0 )
1477       ELSE
1478          message_string = 'variable translation speed used for galilei-' // &
1479             'transformation, which may cause & instabilities in stably ' // &
1480             'stratified regions'
1481          CALL message( 'check_parameters', 'PA0048', 0, 1, 0, 6, 0 )
1482       ENDIF
1483    ENDIF
1484
1485!
1486!-- In case of using a prandtl-layer, calculated (or prescribed) surface
1487!-- fluxes have to be used in the diffusion-terms
1488    IF ( prandtl_layer )  use_surface_fluxes = .TRUE.
1489
1490!
1491!-- Check boundary conditions and set internal variables:
1492!-- Lateral boundary conditions
1493    IF ( bc_lr /= 'cyclic'  .AND.  bc_lr /= 'dirichlet/radiation'  .AND. &
1494         bc_lr /= 'radiation/dirichlet' )  THEN
1495       message_string = 'unknown boundary condition: bc_lr = "' // &
1496                        TRIM( bc_lr ) // '"'
1497       CALL message( 'check_parameters', 'PA0049', 1, 2, 0, 6, 0 )
1498    ENDIF
1499    IF ( bc_ns /= 'cyclic'  .AND.  bc_ns /= 'dirichlet/radiation'  .AND. &
1500         bc_ns /= 'radiation/dirichlet' )  THEN
1501       message_string = 'unknown boundary condition: bc_ns = "' // &
1502                        TRIM( bc_ns ) // '"'
1503       CALL message( 'check_parameters', 'PA0050', 1, 2, 0, 6, 0 )
1504    ENDIF
1505
1506!
1507!-- Internal variables used for speed optimization in if clauses
1508    IF ( bc_lr /= 'cyclic' )               bc_lr_cyc    = .FALSE.
1509    IF ( bc_lr == 'dirichlet/radiation' )  bc_lr_dirrad = .TRUE.
1510    IF ( bc_lr == 'radiation/dirichlet' )  bc_lr_raddir = .TRUE.
1511    IF ( bc_ns /= 'cyclic' )               bc_ns_cyc    = .FALSE.
1512    IF ( bc_ns == 'dirichlet/radiation' )  bc_ns_dirrad = .TRUE.
1513    IF ( bc_ns == 'radiation/dirichlet' )  bc_ns_raddir = .TRUE.
1514
1515!
1516!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-
1517!-- Willimas or Wicker - Skamarock advection scheme. Several schemes
1518!-- and tools do not work with non-cyclic boundary conditions.
1519    IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
1520       IF ( psolver /= 'multigrid' )  THEN
1521          message_string = 'non-cyclic lateral boundaries do not allow ' // &
1522                           'psolver = "' // TRIM( psolver ) // '"'
1523          CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
1524       ENDIF
1525       IF ( momentum_advec /= 'pw-scheme' .AND. &
1526            momentum_advec /= 'ws-scheme')  THEN
1527          message_string = 'non-cyclic lateral boundaries do not allow ' // &
1528                           'momentum_advec = "' // TRIM( momentum_advec ) // '"'
1529          CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
1530       ENDIF
1531       IF ( scalar_advec /= 'pw-scheme' .AND. &
1532            scalar_advec /= 'ws-scheme' )  THEN
1533          message_string = 'non-cyclic lateral boundaries do not allow ' // &
1534                           'scalar_advec = "' // TRIM( scalar_advec ) // '"'
1535          CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 )
1536       ENDIF
1537       IF ( galilei_transformation )  THEN
1538          message_string = 'non-cyclic lateral boundaries do not allow ' // &
1539                           'galilei_transformation = .T.'
1540          CALL message( 'check_parameters', 'PA0054', 1, 2, 0, 6, 0 )
1541       ENDIF
1542    ENDIF
1543
1544!
1545!-- Bottom boundary condition for the turbulent Kinetic energy
1546    IF ( bc_e_b == 'neumann' )  THEN
1547       ibc_e_b = 1
1548    ELSEIF ( bc_e_b == '(u*)**2+neumann' )  THEN
1549       ibc_e_b = 2
1550       IF ( prandtl_layer )  THEN
1551          message_string = 'adjust mixing length = FALSE and bc_e_b = "' // &
1552                           TRIM( bc_e_b ) // '"'
1553          CALL message( 'check_parameters', 'PA0056', 0, 1, 0, 6, 0 )
1554       ENDIF
1555       IF ( .NOT. prandtl_layer )  THEN
1556          bc_e_b = 'neumann'
1557          ibc_e_b = 1
1558          message_string = 'boundary condition bc_e_b changed to "' // &
1559                           TRIM( bc_e_b ) // '"'
1560          CALL message( 'check_parameters', 'PA0057', 0, 1, 0, 6, 0 )
1561       ENDIF
1562    ELSE
1563       message_string = 'unknown boundary condition: bc_e_b = "' // &
1564                        TRIM( bc_e_b ) // '"'
1565       CALL message( 'check_parameters', 'PA0058', 1, 2, 0, 6, 0 )
1566    ENDIF
1567
1568!
1569!-- Boundary conditions for perturbation pressure
1570    IF ( bc_p_b == 'dirichlet' )  THEN
1571       ibc_p_b = 0
1572    ELSEIF ( bc_p_b == 'neumann' )  THEN
1573       ibc_p_b = 1
1574    ELSE
1575       message_string = 'unknown boundary condition: bc_p_b = "' // &
1576                        TRIM( bc_p_b ) // '"'
1577       CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 )
1578    ENDIF
1579
1580    IF ( bc_p_t == 'dirichlet' )  THEN
1581       ibc_p_t = 0
1582    ELSEIF ( bc_p_t == 'neumann' )  THEN
1583       ibc_p_t = 1
1584    ELSE
1585       message_string = 'unknown boundary condition: bc_p_t = "' // &
1586                        TRIM( bc_p_t ) // '"'
1587       CALL message( 'check_parameters', 'PA0061', 1, 2, 0, 6, 0 )
1588    ENDIF
1589
1590!
1591!-- Boundary conditions for potential temperature
1592    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
1593       ibc_pt_b = 2
1594    ELSE
1595       IF ( bc_pt_b == 'dirichlet' )  THEN
1596          ibc_pt_b = 0
1597       ELSEIF ( bc_pt_b == 'neumann' )  THEN
1598          ibc_pt_b = 1
1599       ELSE
1600          message_string = 'unknown boundary condition: bc_pt_b = "' // &
1601                           TRIM( bc_pt_b ) // '"'
1602          CALL message( 'check_parameters', 'PA0062', 1, 2, 0, 6, 0 )
1603       ENDIF
1604    ENDIF
1605
1606    IF ( bc_pt_t == 'dirichlet' )  THEN
1607       ibc_pt_t = 0
1608    ELSEIF ( bc_pt_t == 'neumann' )  THEN
1609       ibc_pt_t = 1
1610    ELSEIF ( bc_pt_t == 'initial_gradient' )  THEN
1611       ibc_pt_t = 2
1612    ELSE
1613       message_string = 'unknown boundary condition: bc_pt_t = "' // &
1614                        TRIM( bc_pt_t ) // '"'
1615       CALL message( 'check_parameters', 'PA0063', 1, 2, 0, 6, 0 )
1616    ENDIF
1617
1618    IF ( surface_heatflux == 9999999.9  )  THEN
1619       constant_heatflux = .FALSE.
1620       IF ( large_scale_forcing )  THEN
1621          IF ( ibc_pt_b == 0 )  THEN
1622             constant_heatflux = .FALSE.
1623          ELSEIF ( ibc_pt_b == 1 )  THEN
1624             constant_heatflux = .TRUE.
1625             IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1626                surface_heatflux = shf_surf(1)
1627             ENDIF
1628          ENDIF
1629       ENDIF
1630    ELSE
1631        constant_heatflux = .TRUE.
1632        IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
1633               large_scale_forcing ) THEN
1634           surface_heatflux = shf_surf(1)
1635        ENDIF
1636    ENDIF
1637
1638    IF ( top_heatflux     == 9999999.9 )  constant_top_heatflux = .FALSE.
1639
1640    IF ( neutral )  THEN
1641
1642       IF ( surface_heatflux /= 0.0  .AND.  surface_heatflux /= 9999999.9 ) &
1643       THEN
1644          message_string = 'heatflux must not be set for pure neutral flow'
1645          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
1646       ENDIF
1647
1648       IF ( top_heatflux /= 0.0  .AND.  top_heatflux /= 9999999.9 ) &
1649       THEN
1650          message_string = 'heatflux must not be set for pure neutral flow'
1651          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
1652       ENDIF
1653
1654    ENDIF
1655
1656    IF ( top_momentumflux_u /= 9999999.9  .AND.  &
1657         top_momentumflux_v /= 9999999.9 )  THEN
1658       constant_top_momentumflux = .TRUE.
1659    ELSEIF (  .NOT. ( top_momentumflux_u == 9999999.9  .AND.  &
1660           top_momentumflux_v == 9999999.9 ) )  THEN
1661       message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' // &
1662                        'must be set'
1663       CALL message( 'check_parameters', 'PA0064', 1, 2, 0, 6, 0 )
1664    ENDIF
1665
1666!
1667!-- A given surface temperature implies Dirichlet boundary condition for
1668!-- temperature. In this case specification of a constant heat flux is
1669!-- forbidden.
1670    IF ( ibc_pt_b == 0  .AND.   constant_heatflux  .AND. &
1671         surface_heatflux /= 0.0 )  THEN
1672       message_string = 'boundary_condition: bc_pt_b = "' // TRIM( bc_pt_b ) //&
1673                        '& is not allowed with constant_heatflux = .TRUE.'
1674       CALL message( 'check_parameters', 'PA0065', 1, 2, 0, 6, 0 )
1675    ENDIF
1676    IF ( constant_heatflux  .AND.  pt_surface_initial_change /= 0.0 )  THEN
1677       WRITE ( message_string, * )  'constant_heatflux = .TRUE. is not allo', &
1678               'wed with pt_surface_initial_change (/=0) = ', &
1679               pt_surface_initial_change
1680       CALL message( 'check_parameters', 'PA0066', 1, 2, 0, 6, 0 )
1681    ENDIF
1682
1683!
1684!-- A given temperature at the top implies Dirichlet boundary condition for
1685!-- temperature. In this case specification of a constant heat flux is
1686!-- forbidden.
1687    IF ( ibc_pt_t == 0  .AND.   constant_top_heatflux  .AND. &
1688         top_heatflux /= 0.0 )  THEN
1689       message_string = 'boundary_condition: bc_pt_t = "' // TRIM( bc_pt_t ) //&
1690                        '" is not allowed with constant_top_heatflux = .TRUE.'
1691       CALL message( 'check_parameters', 'PA0067', 1, 2, 0, 6, 0 )
1692    ENDIF
1693
1694!
1695!-- Boundary conditions for salinity
1696    IF ( ocean )  THEN
1697       IF ( bc_sa_t == 'dirichlet' )  THEN
1698          ibc_sa_t = 0
1699       ELSEIF ( bc_sa_t == 'neumann' )  THEN
1700          ibc_sa_t = 1
1701       ELSE
1702          message_string = 'unknown boundary condition: bc_sa_t = "' // &
1703                           TRIM( bc_sa_t ) // '"'
1704          CALL message( 'check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
1705       ENDIF
1706
1707       IF ( top_salinityflux == 9999999.9 )  constant_top_salinityflux = .FALSE.
1708       IF ( ibc_sa_t == 1  .AND.   top_salinityflux == 9999999.9 )  THEN
1709          message_string = 'boundary condition: bc_sa_t = "' // &
1710                           TRIM( bc_sa_t ) // '" requires to set ' // &
1711                           'top_salinityflux'
1712          CALL message( 'check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
1713       ENDIF
1714
1715!
1716!--    A fixed salinity at the top implies Dirichlet boundary condition for
1717!--    salinity. In this case specification of a constant salinity flux is
1718!--    forbidden.
1719       IF ( ibc_sa_t == 0  .AND.   constant_top_salinityflux  .AND. &
1720            top_salinityflux /= 0.0 )  THEN
1721          message_string = 'boundary condition: bc_sa_t = "' // &
1722                           TRIM( bc_sa_t ) // '" is not allowed with ' // &
1723                           'constant_top_salinityflux = .TRUE.'
1724          CALL message( 'check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
1725       ENDIF
1726
1727    ENDIF
1728
1729!
1730!-- In case of humidity or passive scalar, set boundary conditions for total
1731!-- water content / scalar
1732    IF ( humidity  .OR.  passive_scalar ) THEN
1733       IF ( humidity )  THEN
1734          sq = 'q'
1735       ELSE
1736          sq = 's'
1737       ENDIF
1738       IF ( bc_q_b == 'dirichlet' )  THEN
1739          ibc_q_b = 0
1740       ELSEIF ( bc_q_b == 'neumann' )  THEN
1741          ibc_q_b = 1
1742       ELSE
1743          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // &
1744                           '_b ="' // TRIM( bc_q_b ) // '"'
1745          CALL message( 'check_parameters', 'PA0071', 1, 2, 0, 6, 0 )
1746       ENDIF
1747       IF ( bc_q_t == 'dirichlet' )  THEN
1748          ibc_q_t = 0
1749       ELSEIF ( bc_q_t == 'neumann' )  THEN
1750          ibc_q_t = 1
1751       ELSE
1752          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // &
1753                           '_t ="' // TRIM( bc_q_t ) // '"'
1754          CALL message( 'check_parameters', 'PA0072', 1, 2, 0, 6, 0 )
1755       ENDIF
1756
1757       IF ( surface_waterflux == 9999999.9  )  THEN
1758          constant_waterflux = .FALSE.
1759          IF ( large_scale_forcing )  THEN
1760             IF ( ibc_q_b == 0 )  THEN
1761                constant_waterflux = .FALSE.
1762             ELSEIF ( ibc_q_b == 1 )  THEN
1763                constant_waterflux = .TRUE.
1764                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1765                   surface_waterflux = qsws_surf(1)
1766                ENDIF
1767             ENDIF
1768          ENDIF
1769       ELSE
1770          constant_waterflux = .TRUE.
1771          IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
1772                 large_scale_forcing ) THEN
1773             surface_waterflux = qsws_surf(1)
1774          ENDIF
1775       ENDIF
1776
1777!
1778!--    A given surface humidity implies Dirichlet boundary condition for
1779!--    humidity. In this case specification of a constant water flux is
1780!--    forbidden.
1781       IF ( ibc_q_b == 0  .AND.  constant_waterflux )  THEN
1782          message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b ' // &
1783                           '= "' // TRIM( bc_q_b ) // '" is not allowed wi' // &
1784                           'th prescribed surface flux'
1785          CALL message( 'check_parameters', 'PA0073', 1, 2, 0, 6, 0 )
1786       ENDIF
1787       IF ( constant_waterflux  .AND.  q_surface_initial_change /= 0.0 )  THEN
1788          WRITE( message_string, * )  'a prescribed surface flux is not allo', &
1789                 'wed with ', sq, '_surface_initial_change (/=0) = ', &
1790                 q_surface_initial_change
1791          CALL message( 'check_parameters', 'PA0074', 1, 2, 0, 6, 0 )
1792       ENDIF
1793
1794    ENDIF
1795!
1796!-- Boundary conditions for horizontal components of wind speed
1797    IF ( bc_uv_b == 'dirichlet' )  THEN
1798       ibc_uv_b = 0
1799    ELSEIF ( bc_uv_b == 'neumann' )  THEN
1800       ibc_uv_b = 1
1801       IF ( prandtl_layer )  THEN
1802          message_string = 'boundary condition: bc_uv_b = "' // &
1803               TRIM( bc_uv_b ) // '" is not allowed with prandtl_layer = .TRUE.'
1804          CALL message( 'check_parameters', 'PA0075', 1, 2, 0, 6, 0 )
1805       ENDIF
1806    ELSE
1807       message_string = 'unknown boundary condition: bc_uv_b = "' // &
1808                        TRIM( bc_uv_b ) // '"'
1809       CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )
1810    ENDIF
1811!
1812!-- In case of coupled simulations u and v at the ground in atmosphere will be
1813!-- assigned with the u and v values of the ocean surface
1814    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
1815       ibc_uv_b = 2
1816    ENDIF
1817
1818    IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
1819       bc_uv_t = 'neumann'
1820       ibc_uv_t = 1
1821    ELSE
1822       IF ( bc_uv_t == 'dirichlet' .OR. bc_uv_t == 'dirichlet_0' )  THEN
1823          ibc_uv_t = 0
1824          IF ( bc_uv_t == 'dirichlet_0' )  THEN
1825!
1826!--          Velocities for the initial u,v-profiles are set zero at the top
1827!--          in case of dirichlet_0 conditions
1828             u_init(nzt+1)    = 0.0
1829             v_init(nzt+1)    = 0.0
1830          ENDIF
1831       ELSEIF ( bc_uv_t == 'neumann' )  THEN
1832          ibc_uv_t = 1
1833       ELSE
1834          message_string = 'unknown boundary condition: bc_uv_t = "' // &
1835                           TRIM( bc_uv_t ) // '"'
1836          CALL message( 'check_parameters', 'PA0077', 1, 2, 0, 6, 0 )
1837       ENDIF
1838    ENDIF
1839
1840!
1841!-- Compute and check, respectively, the Rayleigh Damping parameter
1842    IF ( rayleigh_damping_factor == -1.0 )  THEN
1843       rayleigh_damping_factor = 0.0
1844    ELSE
1845       IF ( rayleigh_damping_factor < 0.0 .OR. rayleigh_damping_factor > 1.0 ) &
1846       THEN
1847          WRITE( message_string, * )  'rayleigh_damping_factor = ', &
1848                              rayleigh_damping_factor, ' out of range [0.0,1.0]'
1849          CALL message( 'check_parameters', 'PA0078', 1, 2, 0, 6, 0 )
1850       ENDIF
1851    ENDIF
1852
1853    IF ( rayleigh_damping_height == -1.0 )  THEN
1854       IF ( .NOT. ocean )  THEN
1855          rayleigh_damping_height = 0.66666666666 * zu(nzt)
1856       ELSE
1857          rayleigh_damping_height = 0.66666666666 * zu(nzb)
1858       ENDIF
1859    ELSE
1860       IF ( .NOT. ocean )  THEN
1861          IF ( rayleigh_damping_height < 0.0  .OR. &
1862               rayleigh_damping_height > zu(nzt) )  THEN
1863             WRITE( message_string, * )  'rayleigh_damping_height = ', &
1864                   rayleigh_damping_height, ' out of range [0.0,', zu(nzt), ']'
1865             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
1866          ENDIF
1867       ELSE
1868          IF ( rayleigh_damping_height > 0.0  .OR. &
1869               rayleigh_damping_height < zu(nzb) )  THEN
1870             WRITE( message_string, * )  'rayleigh_damping_height = ', &
1871                   rayleigh_damping_height, ' out of range [0.0,', zu(nzb), ']'
1872             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
1873          ENDIF
1874       ENDIF
1875    ENDIF
1876
1877!
1878!-- Check number of chosen statistic regions. More than 10 regions are not
1879!-- allowed, because so far no more than 10 corresponding output files can
1880!-- be opened (cf. check_open)
1881    IF ( statistic_regions > 9  .OR.  statistic_regions < 0 )  THEN
1882       WRITE ( message_string, * ) 'number of statistic_regions = ', &
1883                   statistic_regions+1, ' but only 10 regions are allowed'
1884       CALL message( 'check_parameters', 'PA0082', 1, 2, 0, 6, 0 )
1885    ENDIF
1886    IF ( normalizing_region > statistic_regions  .OR. &
1887         normalizing_region < 0)  THEN
1888       WRITE ( message_string, * ) 'normalizing_region = ', &
1889                normalizing_region, ' must be >= 0 and <= ',statistic_regions, &
1890                ' (value of statistic_regions)'
1891       CALL message( 'check_parameters', 'PA0083', 1, 2, 0, 6, 0 )
1892    ENDIF
1893
1894!
1895!-- Check the interval for sorting particles.
1896!-- Using particles as cloud droplets requires sorting after each timestep.
1897    IF ( dt_sort_particles /= 0.0  .AND.  cloud_droplets )  THEN
1898       dt_sort_particles = 0.0
1899       message_string = 'dt_sort_particles is reset to 0.0 because of cloud' //&
1900                        '_droplets = .TRUE.'
1901       CALL message( 'check_parameters', 'PA0084', 0, 1, 0, 6, 0 )
1902    ENDIF
1903
1904!
1905!-- Set the default intervals for data output, if necessary
1906!-- NOTE: dt_dosp has already been set in package_parin
1907    IF ( dt_data_output /= 9999999.9 )  THEN
1908       IF ( dt_dopr           == 9999999.9 )  dt_dopr           = dt_data_output
1909       IF ( dt_dopts          == 9999999.9 )  dt_dopts          = dt_data_output
1910       IF ( dt_do2d_xy        == 9999999.9 )  dt_do2d_xy        = dt_data_output
1911       IF ( dt_do2d_xz        == 9999999.9 )  dt_do2d_xz        = dt_data_output
1912       IF ( dt_do2d_yz        == 9999999.9 )  dt_do2d_yz        = dt_data_output
1913       IF ( dt_do3d           == 9999999.9 )  dt_do3d           = dt_data_output
1914       IF ( dt_data_output_av == 9999999.9 )  dt_data_output_av = dt_data_output
1915       DO  mid = 1, max_masks
1916          IF ( dt_domask(mid) == 9999999.9 )  dt_domask(mid)    = dt_data_output
1917       ENDDO
1918    ENDIF
1919
1920!
1921!-- Set the default skip time intervals for data output, if necessary
1922    IF ( skip_time_dopr    == 9999999.9 ) &
1923                                       skip_time_dopr    = skip_time_data_output
1924    IF ( skip_time_dosp    == 9999999.9 ) &
1925                                       skip_time_dosp    = skip_time_data_output
1926    IF ( skip_time_do2d_xy == 9999999.9 ) &
1927                                       skip_time_do2d_xy = skip_time_data_output
1928    IF ( skip_time_do2d_xz == 9999999.9 ) &
1929                                       skip_time_do2d_xz = skip_time_data_output
1930    IF ( skip_time_do2d_yz == 9999999.9 ) &
1931                                       skip_time_do2d_yz = skip_time_data_output
1932    IF ( skip_time_do3d    == 9999999.9 ) &
1933                                       skip_time_do3d    = skip_time_data_output
1934    IF ( skip_time_data_output_av == 9999999.9 ) &
1935                                skip_time_data_output_av = skip_time_data_output
1936    DO  mid = 1, max_masks
1937       IF ( skip_time_domask(mid) == 9999999.9 ) &
1938                                skip_time_domask(mid)    = skip_time_data_output
1939    ENDDO
1940
1941!
1942!-- Check the average intervals (first for 3d-data, then for profiles and
1943!-- spectra)
1944    IF ( averaging_interval > dt_data_output_av )  THEN
1945       WRITE( message_string, * )  'averaging_interval = ', &
1946             averaging_interval, ' must be <= dt_data_output = ', dt_data_output
1947       CALL message( 'check_parameters', 'PA0085', 1, 2, 0, 6, 0 )
1948    ENDIF
1949
1950    IF ( averaging_interval_pr == 9999999.9 )  THEN
1951       averaging_interval_pr = averaging_interval
1952    ENDIF
1953
1954    IF ( averaging_interval_pr > dt_dopr )  THEN
1955       WRITE( message_string, * )  'averaging_interval_pr = ', &
1956             averaging_interval_pr, ' must be <= dt_dopr = ', dt_dopr
1957       CALL message( 'check_parameters', 'PA0086', 1, 2, 0, 6, 0 )
1958    ENDIF
1959
1960    IF ( averaging_interval_sp == 9999999.9 )  THEN
1961       averaging_interval_sp = averaging_interval
1962    ENDIF
1963
1964    IF ( averaging_interval_sp > dt_dosp )  THEN
1965       WRITE( message_string, * )  'averaging_interval_sp = ', &
1966             averaging_interval_sp, ' must be <= dt_dosp = ', dt_dosp
1967       CALL message( 'check_parameters', 'PA0087', 1, 2, 0, 6, 0 )
1968    ENDIF
1969
1970!
1971!-- Set the default interval for profiles entering the temporal average
1972    IF ( dt_averaging_input_pr == 9999999.9 )  THEN
1973       dt_averaging_input_pr = dt_averaging_input
1974    ENDIF
1975
1976!
1977!-- Set the default interval for the output of timeseries to a reasonable
1978!-- value (tries to minimize the number of calls of flow_statistics)
1979    IF ( dt_dots == 9999999.9 )  THEN
1980       IF ( averaging_interval_pr == 0.0 )  THEN
1981          dt_dots = MIN( dt_run_control, dt_dopr )
1982       ELSE
1983          dt_dots = MIN( dt_run_control, dt_averaging_input_pr )
1984       ENDIF
1985    ENDIF
1986
1987!
1988!-- Check the sample rate for averaging (first for 3d-data, then for profiles)
1989    IF ( dt_averaging_input > averaging_interval )  THEN
1990       WRITE( message_string, * )  'dt_averaging_input = ', &
1991                dt_averaging_input, ' must be <= averaging_interval = ', &
1992                averaging_interval
1993       CALL message( 'check_parameters', 'PA0088', 1, 2, 0, 6, 0 )
1994    ENDIF
1995
1996    IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
1997       WRITE( message_string, * )  'dt_averaging_input_pr = ', &
1998                dt_averaging_input_pr, ' must be <= averaging_interval_pr = ', &
1999                averaging_interval_pr
2000       CALL message( 'check_parameters', 'PA0089', 1, 2, 0, 6, 0 )
2001    ENDIF
2002
2003!
2004!-- Set the default value for the integration interval of precipitation amount
2005    IF ( precipitation )  THEN
2006       IF ( precipitation_amount_interval == 9999999.9 )  THEN
2007          precipitation_amount_interval = dt_do2d_xy
2008       ELSE
2009          IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
2010             WRITE( message_string, * )  'precipitation_amount_interval = ', &
2011                 precipitation_amount_interval, ' must not be larger than ', &
2012                 'dt_do2d_xy = ', dt_do2d_xy
2013             CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
2014          ENDIF
2015       ENDIF
2016    ENDIF
2017
2018!
2019!-- Determine the number of output profiles and check whether they are
2020!-- permissible
2021    DO  WHILE ( data_output_pr(dopr_n+1) /= '          ' )
2022
2023       dopr_n = dopr_n + 1
2024       i = dopr_n
2025
2026!
2027!--    Determine internal profile number (for hom, homs)
2028!--    and store height levels
2029       SELECT CASE ( TRIM( data_output_pr(i) ) )
2030
2031          CASE ( 'u', '#u' )
2032             dopr_index(i) = 1
2033             dopr_unit(i)  = 'm/s'
2034             hom(:,2,1,:)  = SPREAD( zu, 2, statistic_regions+1 )
2035             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2036                dopr_initial_index(i) = 5
2037                hom(:,2,5,:)          = SPREAD( zu, 2, statistic_regions+1 )
2038                data_output_pr(i)     = data_output_pr(i)(2:)
2039             ENDIF
2040
2041          CASE ( 'v', '#v' )
2042             dopr_index(i) = 2
2043             dopr_unit(i)  = 'm/s'
2044             hom(:,2,2,:)  = SPREAD( zu, 2, statistic_regions+1 )
2045             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2046                dopr_initial_index(i) = 6
2047                hom(:,2,6,:)          = SPREAD( zu, 2, statistic_regions+1 )
2048                data_output_pr(i)     = data_output_pr(i)(2:)
2049             ENDIF
2050
2051          CASE ( 'w' )
2052             dopr_index(i) = 3
2053             dopr_unit(i)  = 'm/s'
2054             hom(:,2,3,:)  = SPREAD( zw, 2, statistic_regions+1 )
2055
2056          CASE ( 'pt', '#pt' )
2057             IF ( .NOT. cloud_physics ) THEN
2058                dopr_index(i) = 4
2059                dopr_unit(i)  = 'K'
2060                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2061                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2062                   dopr_initial_index(i) = 7
2063                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2064                   hom(nzb,2,7,:)        = 0.0    ! because zu(nzb) is negative
2065                   data_output_pr(i)     = data_output_pr(i)(2:)
2066                ENDIF
2067             ELSE
2068                dopr_index(i) = 43
2069                dopr_unit(i)  = 'K'
2070                hom(:,2,43,:)  = SPREAD( zu, 2, statistic_regions+1 )
2071                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2072                   dopr_initial_index(i) = 28
2073                   hom(:,2,28,:)         = SPREAD( zu, 2, statistic_regions+1 )
2074                   hom(nzb,2,28,:)       = 0.0    ! because zu(nzb) is negative
2075                   data_output_pr(i)     = data_output_pr(i)(2:)
2076                ENDIF
2077             ENDIF
2078
2079          CASE ( 'e' )
2080             dopr_index(i)  = 8
2081             dopr_unit(i)   = 'm2/s2'
2082             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
2083             hom(nzb,2,8,:) = 0.0
2084
2085          CASE ( 'km', '#km' )
2086             dopr_index(i)  = 9
2087             dopr_unit(i)   = 'm2/s'
2088             hom(:,2,9,:)   = SPREAD( zu, 2, statistic_regions+1 )
2089             hom(nzb,2,9,:) = 0.0
2090             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2091                dopr_initial_index(i) = 23
2092                hom(:,2,23,:)         = hom(:,2,9,:)
2093                data_output_pr(i)     = data_output_pr(i)(2:)
2094             ENDIF
2095
2096          CASE ( 'kh', '#kh' )
2097             dopr_index(i)   = 10
2098             dopr_unit(i)    = 'm2/s'
2099             hom(:,2,10,:)   = SPREAD( zu, 2, statistic_regions+1 )
2100             hom(nzb,2,10,:) = 0.0
2101             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2102                dopr_initial_index(i) = 24
2103                hom(:,2,24,:)         = hom(:,2,10,:)
2104                data_output_pr(i)     = data_output_pr(i)(2:)
2105             ENDIF
2106
2107          CASE ( 'l', '#l' )
2108             dopr_index(i)   = 11
2109             dopr_unit(i)    = 'm'
2110             hom(:,2,11,:)   = SPREAD( zu, 2, statistic_regions+1 )
2111             hom(nzb,2,11,:) = 0.0
2112             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2113                dopr_initial_index(i) = 25
2114                hom(:,2,25,:)         = hom(:,2,11,:)
2115                data_output_pr(i)     = data_output_pr(i)(2:)
2116             ENDIF
2117
2118          CASE ( 'w"u"' )
2119             dopr_index(i) = 12
2120             dopr_unit(i)  = 'm2/s2'
2121             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
2122             IF ( prandtl_layer )  hom(nzb,2,12,:) = zu(1)
2123
2124          CASE ( 'w*u*' )
2125             dopr_index(i) = 13
2126             dopr_unit(i)  = 'm2/s2'
2127             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
2128
2129          CASE ( 'w"v"' )
2130             dopr_index(i) = 14
2131             dopr_unit(i)  = 'm2/s2'
2132             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
2133             IF ( prandtl_layer )  hom(nzb,2,14,:) = zu(1)
2134
2135          CASE ( 'w*v*' )
2136             dopr_index(i) = 15
2137             dopr_unit(i)  = 'm2/s2'
2138             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
2139
2140          CASE ( 'w"pt"' )
2141             dopr_index(i) = 16
2142             dopr_unit(i)  = 'K m/s'
2143             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
2144
2145          CASE ( 'w*pt*' )
2146             dopr_index(i) = 17
2147             dopr_unit(i)  = 'K m/s'
2148             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
2149
2150          CASE ( 'wpt' )
2151             dopr_index(i) = 18
2152             dopr_unit(i)  = 'K m/s'
2153             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
2154
2155          CASE ( 'wu' )
2156             dopr_index(i) = 19
2157             dopr_unit(i)  = 'm2/s2'
2158             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
2159             IF ( prandtl_layer )  hom(nzb,2,19,:) = zu(1)
2160
2161          CASE ( 'wv' )
2162             dopr_index(i) = 20
2163             dopr_unit(i)  = 'm2/s2'
2164             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
2165             IF ( prandtl_layer )  hom(nzb,2,20,:) = zu(1)
2166
2167          CASE ( 'w*pt*BC' )
2168             dopr_index(i) = 21
2169             dopr_unit(i)  = 'K m/s'
2170             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
2171
2172          CASE ( 'wptBC' )
2173             dopr_index(i) = 22
2174             dopr_unit(i)  = 'K m/s'
2175             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
2176
2177          CASE ( 'sa', '#sa' )
2178             IF ( .NOT. ocean )  THEN
2179                message_string = 'data_output_pr = ' // &
2180                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2181                                 'lemented for ocean = .FALSE.'
2182                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2183             ELSE
2184                dopr_index(i) = 23
2185                dopr_unit(i)  = 'psu'
2186                hom(:,2,23,:) = SPREAD( zu, 2, statistic_regions+1 )
2187                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2188                   dopr_initial_index(i) = 26
2189                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2190                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
2191                   data_output_pr(i)     = data_output_pr(i)(2:)
2192                ENDIF
2193             ENDIF
2194
2195          CASE ( 'u*2' )
2196             dopr_index(i) = 30
2197             dopr_unit(i)  = 'm2/s2'
2198             hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
2199
2200          CASE ( 'v*2' )
2201             dopr_index(i) = 31
2202             dopr_unit(i)  = 'm2/s2'
2203             hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
2204
2205          CASE ( 'w*2' )
2206             dopr_index(i) = 32
2207             dopr_unit(i)  = 'm2/s2'
2208             hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
2209
2210          CASE ( 'pt*2' )
2211             dopr_index(i) = 33
2212             dopr_unit(i)  = 'K2'
2213             hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
2214
2215          CASE ( 'e*' )
2216             dopr_index(i) = 34
2217             dopr_unit(i)  = 'm2/s2'
2218             hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
2219
2220          CASE ( 'w*2pt*' )
2221             dopr_index(i) = 35
2222             dopr_unit(i)  = 'K m2/s2'
2223             hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
2224
2225          CASE ( 'w*pt*2' )
2226             dopr_index(i) = 36
2227             dopr_unit(i)  = 'K2 m/s'
2228             hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
2229
2230          CASE ( 'w*e*' )
2231             dopr_index(i) = 37
2232             dopr_unit(i)  = 'm3/s3'
2233             hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
2234
2235          CASE ( 'w*3' )
2236             dopr_index(i) = 38
2237             dopr_unit(i)  = 'm3/s3'
2238             hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
2239
2240          CASE ( 'Sw' )
2241             dopr_index(i) = 39
2242             dopr_unit(i)  = 'none'
2243             hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
2244
2245          CASE ( 'p' )
2246             dopr_index(i) = 40
2247             dopr_unit(i)  = 'Pa'
2248             hom(:,2,40,:) = SPREAD( zu, 2, statistic_regions+1 )
2249
2250          CASE ( 'q', '#q' )
2251             IF ( .NOT. humidity )  THEN
2252                message_string = 'data_output_pr = ' // &
2253                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2254                                 'lemented for humidity = .FALSE.'
2255                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2256             ELSE
2257                dopr_index(i) = 41
2258                dopr_unit(i)  = 'kg/kg'
2259                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2260                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2261                   dopr_initial_index(i) = 26
2262                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2263                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
2264                   data_output_pr(i)     = data_output_pr(i)(2:)
2265                ENDIF
2266             ENDIF
2267
2268          CASE ( 's', '#s' )
2269             IF ( .NOT. passive_scalar )  THEN
2270                message_string = 'data_output_pr = ' // &
2271                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2272                                 'lemented for passive_scalar = .FALSE.'
2273                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2274             ELSE
2275                dopr_index(i) = 41
2276                dopr_unit(i)  = 'kg/m3'
2277                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2278                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2279                   dopr_initial_index(i) = 26
2280                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2281                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
2282                   data_output_pr(i)     = data_output_pr(i)(2:)
2283                ENDIF
2284             ENDIF
2285
2286          CASE ( 'qv', '#qv' )
2287             IF ( .NOT. cloud_physics ) THEN
2288                dopr_index(i) = 41
2289                dopr_unit(i)  = 'kg/kg'
2290                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2291                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2292                   dopr_initial_index(i) = 26
2293                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2294                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
2295                   data_output_pr(i)     = data_output_pr(i)(2:)
2296                ENDIF
2297             ELSE
2298                dopr_index(i) = 42
2299                dopr_unit(i)  = 'kg/kg'
2300                hom(:,2,42,:) = SPREAD( zu, 2, statistic_regions+1 )
2301                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2302                   dopr_initial_index(i) = 27
2303                   hom(:,2,27,:)         = SPREAD( zu, 2, statistic_regions+1 )
2304                   hom(nzb,2,27,:)       = 0.0    ! weil zu(nzb) negativ ist
2305                   data_output_pr(i)     = data_output_pr(i)(2:)
2306                ENDIF
2307             ENDIF
2308
2309          CASE ( 'lpt', '#lpt' )
2310             IF ( .NOT. cloud_physics ) THEN
2311                message_string = 'data_output_pr = ' // &
2312                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2313                                 'lemented for cloud_physics = .FALSE.'
2314                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2315             ELSE
2316                dopr_index(i) = 4
2317                dopr_unit(i)  = 'K'
2318                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2319                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2320                   dopr_initial_index(i) = 7
2321                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2322                   hom(nzb,2,7,:)        = 0.0    ! weil zu(nzb) negativ ist
2323                   data_output_pr(i)     = data_output_pr(i)(2:)
2324                ENDIF
2325             ENDIF
2326
2327          CASE ( 'vpt', '#vpt' )
2328             dopr_index(i) = 44
2329             dopr_unit(i)  = 'K'
2330             hom(:,2,44,:) = SPREAD( zu, 2, statistic_regions+1 )
2331             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2332                dopr_initial_index(i) = 29
2333                hom(:,2,29,:)         = SPREAD( zu, 2, statistic_regions+1 )
2334                hom(nzb,2,29,:)       = 0.0    ! weil zu(nzb) negativ ist
2335                data_output_pr(i)     = data_output_pr(i)(2:)
2336             ENDIF
2337
2338          CASE ( 'w"vpt"' )
2339             dopr_index(i) = 45
2340             dopr_unit(i)  = 'K m/s'
2341             hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
2342
2343          CASE ( 'w*vpt*' )
2344             dopr_index(i) = 46
2345             dopr_unit(i)  = 'K m/s'
2346             hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
2347
2348          CASE ( 'wvpt' )
2349             dopr_index(i) = 47
2350             dopr_unit(i)  = 'K m/s'
2351             hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
2352
2353          CASE ( 'w"q"' )
2354             IF ( .NOT. humidity )  THEN
2355                message_string = 'data_output_pr = ' // &
2356                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2357                                 'lemented for humidity = .FALSE.'
2358                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2359             ELSE
2360                dopr_index(i) = 48
2361                dopr_unit(i)  = 'kg/kg m/s'
2362                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2363             ENDIF
2364
2365          CASE ( 'w*q*' )
2366             IF ( .NOT. humidity )  THEN
2367                message_string = 'data_output_pr = ' // &
2368                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2369                                 'lemented for humidity = .FALSE.'
2370                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2371             ELSE
2372                dopr_index(i) = 49
2373                dopr_unit(i)  = 'kg/kg m/s'
2374                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2375             ENDIF
2376
2377          CASE ( 'wq' )
2378             IF ( .NOT. humidity )  THEN
2379                message_string = 'data_output_pr = ' // &
2380                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2381                                 'lemented for humidity = .FALSE.'
2382                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2383             ELSE
2384                dopr_index(i) = 50
2385                dopr_unit(i)  = 'kg/kg m/s'
2386                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2387             ENDIF
2388
2389          CASE ( 'w"s"' )
2390             IF ( .NOT. passive_scalar ) THEN
2391                message_string = 'data_output_pr = ' // &
2392                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2393                                 'lemented for passive_scalar = .FALSE.'
2394                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2395             ELSE
2396                dopr_index(i) = 48
2397                dopr_unit(i)  = 'kg/m3 m/s'
2398                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2399             ENDIF
2400
2401          CASE ( 'w*s*' )
2402             IF ( .NOT. passive_scalar ) THEN
2403                message_string = 'data_output_pr = ' // &
2404                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2405                                 'lemented for passive_scalar = .FALSE.'
2406                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2407             ELSE
2408                dopr_index(i) = 49
2409                dopr_unit(i)  = 'kg/m3 m/s'
2410                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2411             ENDIF
2412
2413          CASE ( 'ws' )
2414             IF ( .NOT. passive_scalar ) THEN
2415                message_string = 'data_output_pr = ' // &
2416                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2417                                 'lemented for passive_scalar = .FALSE.'
2418                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2419             ELSE
2420                dopr_index(i) = 50
2421                dopr_unit(i)  = 'kg/m3 m/s'
2422                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2423             ENDIF
2424
2425          CASE ( 'w"qv"' )
2426             IF ( humidity  .AND.  .NOT. cloud_physics ) &
2427             THEN
2428                dopr_index(i) = 48
2429                dopr_unit(i)  = 'kg/kg m/s'
2430                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2431             ELSEIF( humidity .AND. cloud_physics ) THEN
2432                dopr_index(i) = 51
2433                dopr_unit(i)  = 'kg/kg m/s'
2434                hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
2435             ELSE
2436                message_string = 'data_output_pr = ' // &
2437                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2438                                 'lemented for cloud_physics = .FALSE. an&' // &
2439                                 'd humidity = .FALSE.'
2440                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2441             ENDIF
2442
2443          CASE ( 'w*qv*' )
2444             IF ( humidity  .AND.  .NOT. cloud_physics ) &
2445             THEN
2446                dopr_index(i) = 49
2447                dopr_unit(i)  = 'kg/kg m/s'
2448                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2449             ELSEIF( humidity .AND. cloud_physics ) THEN
2450                dopr_index(i) = 52
2451                dopr_unit(i)  = 'kg/kg m/s'
2452                hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
2453             ELSE
2454                message_string = 'data_output_pr = ' // &
2455                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2456                                 'lemented for cloud_physics = .FALSE. an&' // &
2457                                 'd humidity = .FALSE.'
2458                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2459             ENDIF
2460
2461          CASE ( 'wqv' )
2462             IF ( humidity  .AND.  .NOT. cloud_physics ) &
2463             THEN
2464                dopr_index(i) = 50
2465                dopr_unit(i)  = 'kg/kg m/s'
2466                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2467             ELSEIF( humidity .AND. cloud_physics ) THEN
2468                dopr_index(i) = 53
2469                dopr_unit(i)  = 'kg/kg m/s'
2470                hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
2471             ELSE
2472                message_string = 'data_output_pr = ' // &
2473                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2474                                 'lemented for cloud_physics = .FALSE. an&' // &
2475                                 'd humidity = .FALSE.'
2476                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2477             ENDIF
2478
2479          CASE ( 'ql' )
2480             IF ( .NOT. cloud_physics  .AND.  .NOT. cloud_droplets )  THEN
2481                message_string = 'data_output_pr = ' // &
2482                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2483                                 'lemented for cloud_physics = .FALSE. or'  // &
2484                                 '&cloud_droplets = .FALSE.'
2485                CALL message( 'check_parameters', 'PA0096', 1, 2, 0, 6, 0 )
2486             ELSE
2487                dopr_index(i) = 54
2488                dopr_unit(i)  = 'kg/kg'
2489                hom(:,2,54,:)  = SPREAD( zu, 2, statistic_regions+1 )
2490             ENDIF
2491
2492          CASE ( 'w*u*u*:dz' )
2493             dopr_index(i) = 55
2494             dopr_unit(i)  = 'm2/s3'
2495             hom(:,2,55,:) = SPREAD( zu, 2, statistic_regions+1 )
2496
2497          CASE ( 'w*p*:dz' )
2498             dopr_index(i) = 56
2499             dopr_unit(i)  = 'm2/s3'
2500             hom(:,2,56,:) = SPREAD( zw, 2, statistic_regions+1 )
2501
2502          CASE ( 'w"e:dz' )
2503             dopr_index(i) = 57
2504             dopr_unit(i)  = 'm2/s3'
2505             hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
2506
2507
2508          CASE ( 'u"pt"' )
2509             dopr_index(i) = 58
2510             dopr_unit(i)  = 'K m/s'
2511             hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
2512
2513          CASE ( 'u*pt*' )
2514             dopr_index(i) = 59
2515             dopr_unit(i)  = 'K m/s'
2516             hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
2517
2518          CASE ( 'upt_t' )
2519             dopr_index(i) = 60
2520             dopr_unit(i)  = 'K m/s'
2521             hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
2522
2523          CASE ( 'v"pt"' )
2524             dopr_index(i) = 61
2525             dopr_unit(i)  = 'K m/s'
2526             hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
2527             
2528          CASE ( 'v*pt*' )
2529             dopr_index(i) = 62
2530             dopr_unit(i)  = 'K m/s'
2531             hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
2532
2533          CASE ( 'vpt_t' )
2534             dopr_index(i) = 63
2535             dopr_unit(i)  = 'K m/s'
2536             hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
2537
2538          CASE ( 'rho' )
2539             IF ( .NOT. ocean ) THEN
2540                message_string = 'data_output_pr = ' // &
2541                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2542                                 'lemented for ocean = .FALSE.'
2543                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2544             ELSE
2545                dopr_index(i) = 64
2546                dopr_unit(i)  = 'kg/m3'
2547                hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
2548                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2549                   dopr_initial_index(i) = 77
2550                   hom(:,2,77,:)         = SPREAD( zu, 2, statistic_regions+1 )
2551                   hom(nzb,2,77,:)       = 0.0    ! because zu(nzb) is negative
2552                   data_output_pr(i)     = data_output_pr(i)(2:)
2553                ENDIF
2554             ENDIF
2555
2556          CASE ( 'w"sa"' )
2557             IF ( .NOT. ocean ) THEN
2558                message_string = 'data_output_pr = ' // &
2559                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2560                                 'lemented for ocean = .FALSE.'
2561                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2562             ELSE
2563                dopr_index(i) = 65
2564                dopr_unit(i)  = 'psu m/s'
2565                hom(:,2,65,:) = SPREAD( zw, 2, statistic_regions+1 )
2566             ENDIF
2567
2568          CASE ( 'w*sa*' )
2569             IF ( .NOT. ocean ) THEN
2570                message_string = 'data_output_pr = ' // &
2571                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2572                                 'lemented for ocean = .FALSE.'
2573                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2574             ELSE
2575                dopr_index(i) = 66
2576                dopr_unit(i)  = 'psu m/s'
2577                hom(:,2,66,:) = SPREAD( zw, 2, statistic_regions+1 )
2578             ENDIF
2579
2580          CASE ( 'wsa' )
2581             IF ( .NOT. ocean ) THEN
2582                message_string = 'data_output_pr = ' // &
2583                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2584                                 'lemented for ocean = .FALSE.'
2585                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2586             ELSE
2587                dopr_index(i) = 67
2588                dopr_unit(i)  = 'psu m/s'
2589                hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
2590             ENDIF
2591
2592          CASE ( 'w*p*' )
2593             dopr_index(i) = 68
2594             dopr_unit(i)  = 'm3/s3'
2595             hom(:,2,68,:) = SPREAD( zu, 2, statistic_regions+1 )
2596
2597          CASE ( 'w"e' )
2598             dopr_index(i) = 69
2599             dopr_unit(i)  = 'm3/s3'
2600             hom(:,2,69,:) = SPREAD( zu, 2, statistic_regions+1 )
2601
2602          CASE ( 'q*2' )
2603             IF ( .NOT. humidity )  THEN
2604                message_string = 'data_output_pr = ' // &
2605                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2606                                 'lemented for humidity = .FALSE.'
2607                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2608             ELSE
2609                dopr_index(i) = 70
2610                dopr_unit(i)  = 'kg2/kg2'
2611                hom(:,2,70,:) = SPREAD( zu, 2, statistic_regions+1 )
2612             ENDIF
2613
2614          CASE ( 'prho' )
2615             IF ( .NOT. ocean ) THEN
2616                message_string = 'data_output_pr = ' // &
2617                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2618                                 'lemented for ocean = .FALSE.'
2619                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2620             ELSE
2621                dopr_index(i) = 71
2622                dopr_unit(i)  = 'kg/m3'
2623                hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 )
2624             ENDIF
2625
2626          CASE ( 'hyp' )
2627             dopr_index(i) = 72
2628             dopr_unit(i)  = 'dbar'
2629             hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
2630
2631          CASE ( 'nr' )
2632             IF ( .NOT. cloud_physics )  THEN
2633                message_string = 'data_output_pr = ' // &
2634                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2635                                 'lemented for cloud_physics = .FALSE.'
2636                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2637             ELSEIF ( icloud_scheme /= 0 )  THEN
2638                message_string = 'data_output_pr = ' // &
2639                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2640                                 'lemented for cloud_scheme /= seifert_beheng'
2641                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2642             ELSEIF ( .NOT. precipitation )  THEN
2643                message_string = 'data_output_pr = ' // &
2644                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2645                                 'lemented for precipitation = .FALSE.'
2646                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
2647             ELSE
2648                dopr_index(i) = 73
2649                dopr_unit(i)  = '1/m3'
2650                hom(:,2,73,:)  = SPREAD( zu, 2, statistic_regions+1 )
2651             ENDIF
2652
2653          CASE ( 'qr' )
2654             IF ( .NOT. cloud_physics )  THEN
2655                message_string = 'data_output_pr = ' // &
2656                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2657                                 'lemented for cloud_physics = .FALSE.'
2658                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2659             ELSEIF ( icloud_scheme /= 0 )  THEN
2660                message_string = 'data_output_pr = ' // &
2661                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2662                                 'lemented for cloud_scheme /= seifert_beheng'
2663                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2664             ELSEIF ( .NOT. precipitation )  THEN
2665                message_string = 'data_output_pr = ' // &
2666                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2667                                 'lemented for precipitation = .FALSE.'
2668                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
2669             ELSE
2670                dopr_index(i) = 74
2671                dopr_unit(i)  = 'kg/kg'
2672                hom(:,2,74,:)  = SPREAD( zu, 2, statistic_regions+1 )
2673             ENDIF
2674
2675          CASE ( 'qc' )
2676             IF ( .NOT. cloud_physics )  THEN
2677                message_string = 'data_output_pr = ' // &
2678                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2679                                 'lemented for cloud_physics = .FALSE.'
2680                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2681             ELSEIF ( icloud_scheme /= 0 )  THEN
2682                message_string = 'data_output_pr = ' // &
2683                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2684                                 'lemented for cloud_scheme /= seifert_beheng'
2685                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2686             ELSE
2687                dopr_index(i) = 75
2688                dopr_unit(i)  = 'kg/kg'
2689                hom(:,2,75,:)  = SPREAD( zu, 2, statistic_regions+1 )
2690             ENDIF
2691
2692          CASE ( 'prr' )
2693             IF ( .NOT. cloud_physics )  THEN
2694                message_string = 'data_output_pr = ' // &
2695                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2696                                 'lemented for cloud_physics = .FALSE.'
2697                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2698             ELSEIF ( icloud_scheme /= 0 )  THEN
2699                message_string = 'data_output_pr = ' // &
2700                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2701                                 'lemented for cloud_scheme /= seifert_beheng'
2702                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
2703             ELSEIF ( .NOT. precipitation )  THEN
2704                message_string = 'data_output_pr = ' // &
2705                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2706                                 'lemented for precipitation = .FALSE.'
2707                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
2708
2709             ELSE
2710                dopr_index(i) = 76
2711                dopr_unit(i)  = 'kg/kg m/s'
2712                hom(:,2,76,:)  = SPREAD( zu, 2, statistic_regions+1 )
2713             ENDIF
2714
2715          CASE ( 'ug' )
2716             dopr_index(i) = 78
2717             dopr_unit(i)  = 'm/s'
2718             hom(:,2,78,:) = SPREAD( zu, 2, statistic_regions+1 )
2719
2720          CASE ( 'vg' )
2721             dopr_index(i) = 79
2722             dopr_unit(i)  = 'm/s'
2723             hom(:,2,79,:) = SPREAD( zu, 2, statistic_regions+1 )
2724
2725          CASE DEFAULT
2726
2727             CALL user_check_data_output_pr( data_output_pr(i), i, unit )
2728
2729             IF ( unit == 'illegal' )  THEN
2730                IF ( data_output_pr_user(1) /= ' ' )  THEN
2731                   message_string = 'illegal value for data_output_pr or ' // &
2732                                    'data_output_pr_user = "' // &
2733                                    TRIM( data_output_pr(i) ) // '"'
2734                   CALL message( 'check_parameters', 'PA0097', 1, 2, 0, 6, 0 )
2735                ELSE
2736                   message_string = 'illegal value for data_output_pr = "' // &
2737                                    TRIM( data_output_pr(i) ) // '"'
2738                   CALL message( 'check_parameters', 'PA0098', 1, 2, 0, 6, 0 )
2739                ENDIF
2740             ENDIF
2741
2742       END SELECT
2743
2744    ENDDO
2745
2746
2747!
2748!-- Append user-defined data output variables to the standard data output
2749    IF ( data_output_user(1) /= ' ' )  THEN
2750       i = 1
2751       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
2752          i = i + 1
2753       ENDDO
2754       j = 1
2755       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 100 )
2756          IF ( i > 100 )  THEN
2757             message_string = 'number of output quantitities given by data' // &
2758                '_output and data_output_user exceeds the limit of 100'
2759             CALL message( 'check_parameters', 'PA0102', 1, 2, 0, 6, 0 )
2760          ENDIF
2761          data_output(i) = data_output_user(j)
2762          i = i + 1
2763          j = j + 1
2764       ENDDO
2765    ENDIF
2766
2767!
2768!-- Check and set steering parameters for 2d/3d data output and averaging
2769    i   = 1
2770    DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 100 )
2771!
2772!--    Check for data averaging
2773       ilen = LEN_TRIM( data_output(i) )
2774       j = 0                                                 ! no data averaging
2775       IF ( ilen > 3 )  THEN
2776          IF ( data_output(i)(ilen-2:ilen) == '_av' )  THEN
2777             j = 1                                           ! data averaging
2778             data_output(i) = data_output(i)(1:ilen-3)
2779          ENDIF
2780       ENDIF
2781!
2782!--    Check for cross section or volume data
2783       ilen = LEN_TRIM( data_output(i) )
2784       k = 0                                                   ! 3d data
2785       var = data_output(i)(1:ilen)
2786       IF ( ilen > 3 )  THEN
2787          IF ( data_output(i)(ilen-2:ilen) == '_xy'  .OR. &
2788               data_output(i)(ilen-2:ilen) == '_xz'  .OR. &
2789               data_output(i)(ilen-2:ilen) == '_yz' )  THEN
2790             k = 1                                             ! 2d data
2791             var = data_output(i)(1:ilen-3)
2792          ENDIF
2793       ENDIF
2794!
2795!--    Check for allowed value and set units
2796       SELECT CASE ( TRIM( var ) )
2797
2798          CASE ( 'e' )
2799             IF ( constant_diffusion )  THEN
2800                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2801                                 'res constant_diffusion = .FALSE.'
2802                CALL message( 'check_parameters', 'PA0103', 1, 2, 0, 6, 0 )
2803             ENDIF
2804             unit = 'm2/s2'
2805
2806          CASE ( 'lpt' )
2807             IF ( .NOT. cloud_physics )  THEN
2808                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2809                         'res cloud_physics = .TRUE.'
2810                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2811             ENDIF
2812             unit = 'K'
2813
2814          CASE ( 'nr' )
2815             IF ( .NOT. cloud_physics )  THEN
2816                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2817                         'res cloud_physics = .TRUE.'
2818                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2819             ELSEIF ( icloud_scheme /= 0 )  THEN
2820                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2821                         'res cloud_scheme = seifert_beheng'
2822                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
2823             ENDIF
2824             unit = '1/m3'
2825
2826          CASE ( 'pc', 'pr' )
2827             IF ( .NOT. particle_advection )  THEN
2828                message_string = 'output of "' // TRIM( var ) // '" requir' // &
2829                   'es a "particles_par"-NAMELIST in the parameter file (PARIN)'
2830                CALL message( 'check_parameters', 'PA0104', 1, 2, 0, 6, 0 )
2831             ENDIF
2832             IF ( TRIM( var ) == 'pc' )  unit = 'number'
2833             IF ( TRIM( var ) == 'pr' )  unit = 'm'
2834
2835          CASE ( 'prr' )
2836             IF ( .NOT. cloud_physics )  THEN
2837                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2838                         'res cloud_physics = .TRUE.'
2839                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2840             ELSEIF ( icloud_scheme /= 0 )  THEN
2841                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2842                         'res cloud_scheme = seifert_beheng'
2843                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
2844             ELSEIF ( .NOT. precipitation )  THEN
2845                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2846                                 'res precipitation = .TRUE.'
2847                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
2848             ENDIF
2849             unit = 'kg/kg m/s'
2850
2851          CASE ( 'q', 'vpt' )
2852             IF ( .NOT. humidity )  THEN
2853                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2854                                 'res humidity = .TRUE.'
2855                CALL message( 'check_parameters', 'PA0105', 1, 2, 0, 6, 0 )
2856             ENDIF
2857             IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
2858             IF ( TRIM( var ) == 'vpt' )  unit = 'K'
2859
2860          CASE ( 'qc' )
2861             IF ( .NOT. cloud_physics )  THEN
2862                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2863                         'res cloud_physics = .TRUE.'
2864                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2865             ELSEIF ( icloud_scheme /= 0 ) THEN
2866                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2867                         'res cloud_scheme = seifert_beheng'
2868                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
2869             ENDIF
2870             unit = 'kg/kg'
2871
2872          CASE ( 'ql' )
2873             IF ( .NOT. ( cloud_physics  .OR.  cloud_droplets ) )  THEN
2874                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2875                         'res cloud_physics = .TRUE. or cloud_droplets = .TRUE.'
2876                CALL message( 'check_parameters', 'PA0106', 1, 2, 0, 6, 0 )
2877             ENDIF
2878             unit = 'kg/kg'
2879
2880          CASE ( 'ql_c', 'ql_v', 'ql_vp' )
2881             IF ( .NOT. cloud_droplets )  THEN
2882                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2883                                 'res cloud_droplets = .TRUE.'
2884                CALL message( 'check_parameters', 'PA0107', 1, 2, 0, 6, 0 )
2885             ENDIF
2886             IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
2887             IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
2888             IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
2889
2890          CASE ( 'qr' )
2891             IF ( .NOT. cloud_physics )  THEN
2892                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2893                         'res cloud_physics = .TRUE.'
2894                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2895             ELSEIF ( icloud_scheme /= 0 ) THEN
2896                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2897                         'res cloud_scheme = seifert_beheng'
2898                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
2899             ELSEIF ( .NOT. precipitation )  THEN
2900                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2901                                 'res precipitation = .TRUE.'
2902                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
2903             ENDIF
2904             unit = 'kg/kg'
2905
2906          CASE ( 'qv' )
2907             IF ( .NOT. cloud_physics )  THEN
2908                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2909                                 'res cloud_physics = .TRUE.'
2910                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2911             ENDIF
2912             unit = 'kg/kg'
2913
2914          CASE ( 'rho' )
2915             IF ( .NOT. ocean )  THEN
2916                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2917                                 'res ocean = .TRUE.'
2918                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
2919             ENDIF
2920             unit = 'kg/m3'
2921
2922          CASE ( 's' )
2923             IF ( .NOT. passive_scalar )  THEN
2924                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2925                                 'res passive_scalar = .TRUE.'
2926                CALL message( 'check_parameters', 'PA0110', 1, 2, 0, 6, 0 )
2927             ENDIF
2928             unit = 'conc'
2929
2930          CASE ( 'sa' )
2931             IF ( .NOT. ocean )  THEN
2932                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2933                                 'res ocean = .TRUE.'
2934                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
2935             ENDIF
2936             unit = 'psu'
2937
2938          CASE ( 'u*', 't*', 'lwp*', 'pra*', 'prr*', 'qsws*', 'shf*', 'z0*', 'z0h*' )
2939             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
2940                message_string = 'illegal value for data_output: "' // &
2941                                 TRIM( var ) // '" & only 2d-horizontal ' // &
2942                                 'cross sections are allowed for this value'
2943                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
2944             ENDIF
2945             IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. cloud_physics )  THEN
2946                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2947                                 'res cloud_physics = .TRUE.'
2948                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2949             ENDIF
2950             IF ( TRIM( var ) == 'pra*'  .AND.  .NOT. precipitation )  THEN
2951                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2952                                 'res precipitation = .TRUE.'
2953                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
2954             ENDIF
2955             IF ( TRIM( var ) == 'pra*'  .AND.  j == 1 )  THEN
2956                message_string = 'temporal averaging of precipitation ' // &
2957                          'amount "' // TRIM( var ) // '" is not possible'
2958                CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
2959             ENDIF
2960             IF ( TRIM( var ) == 'prr*'  .AND.  .NOT. precipitation )  THEN
2961                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2962                                 'res precipitation = .TRUE.'
2963                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
2964             ENDIF
2965             IF ( TRIM( var ) == 'qsws*'  .AND.  .NOT. humidity )  THEN
2966                message_string = 'output of "' // TRIM( var ) // '" requi' // &
2967                                 'res humidity = .TRUE.'
2968                CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
2969             ENDIF
2970
2971             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/kg*m'
2972             IF ( TRIM( var ) == 'pra*'   )  unit = 'mm'
2973             IF ( TRIM( var ) == 'prr*'   )  unit = 'mm/s'
2974             IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
2975             IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
2976             IF ( TRIM( var ) == 't*'     )  unit = 'K'
2977             IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
2978             IF ( TRIM( var ) == 'z0*'    )  unit = 'm'
2979             IF ( TRIM( var ) == 'z0h*'   )  unit = 'm'
2980
2981
2982          CASE ( 'p', 'pt', 'u', 'v', 'w' )
2983             IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
2984             IF ( TRIM( var ) == 'pt' )  unit = 'K'
2985             IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
2986             IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
2987             IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
2988             CONTINUE
2989
2990          CASE DEFAULT
2991             CALL user_check_data_output( var, unit )
2992
2993             IF ( unit == 'illegal' )  THEN
2994                IF ( data_output_user(1) /= ' ' )  THEN
2995                   message_string = 'illegal value for data_output or ' // &
2996                         'data_output_user = "' // TRIM( data_output(i) ) // '"'
2997                   CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
2998                ELSE
2999                   message_string = 'illegal value for data_output =' // &
3000                                    TRIM( data_output(i) ) // '"'
3001                   CALL message( 'check_parameters', 'PA0115', 1, 2, 0, 6, 0 )
3002                ENDIF
3003             ENDIF
3004
3005       END SELECT
3006!
3007!--    Set the internal steering parameters appropriately
3008       IF ( k == 0 )  THEN
3009          do3d_no(j)              = do3d_no(j) + 1
3010          do3d(j,do3d_no(j))      = data_output(i)
3011          do3d_unit(j,do3d_no(j)) = unit
3012       ELSE
3013          do2d_no(j)              = do2d_no(j) + 1
3014          do2d(j,do2d_no(j))      = data_output(i)
3015          do2d_unit(j,do2d_no(j)) = unit
3016          IF ( data_output(i)(ilen-2:ilen) == '_xy' )  THEN
3017             data_output_xy(j) = .TRUE.
3018          ENDIF
3019          IF ( data_output(i)(ilen-2:ilen) == '_xz' )  THEN
3020             data_output_xz(j) = .TRUE.
3021          ENDIF
3022          IF ( data_output(i)(ilen-2:ilen) == '_yz' )  THEN
3023             data_output_yz(j) = .TRUE.
3024          ENDIF
3025       ENDIF
3026
3027       IF ( j == 1 )  THEN
3028!
3029!--       Check, if variable is already subject to averaging
3030          found = .FALSE.
3031          DO  k = 1, doav_n
3032             IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
3033          ENDDO
3034
3035          IF ( .NOT. found )  THEN
3036             doav_n = doav_n + 1
3037             doav(doav_n) = var
3038          ENDIF
3039       ENDIF
3040
3041       i = i + 1
3042    ENDDO
3043
3044!
3045!-- Averaged 2d or 3d output requires that an averaging interval has been set
3046    IF ( doav_n > 0  .AND.  averaging_interval == 0.0 )  THEN
3047       WRITE( message_string, * )  'output of averaged quantity "',            &
3048                                   TRIM( doav(1) ), '_av" requires to set a ', &
3049                                   'non-zero & averaging interval'
3050       CALL message( 'check_parameters', 'PA0323', 1, 2, 0, 6, 0 )
3051    ENDIF
3052
3053!
3054!-- Check sectional planes and store them in one shared array
3055    IF ( ANY( section_xy > nz + 1 ) )  THEN
3056       WRITE( message_string, * )  'section_xy must be <= nz + 1 = ', nz + 1
3057       CALL message( 'check_parameters', 'PA0319', 1, 2, 0, 6, 0 )
3058    ENDIF
3059    IF ( ANY( section_xz > ny + 1 ) )  THEN
3060       WRITE( message_string, * )  'section_xz must be <= ny + 1 = ', ny + 1
3061       CALL message( 'check_parameters', 'PA0320', 1, 2, 0, 6, 0 )
3062    ENDIF
3063    IF ( ANY( section_yz > nx + 1 ) )  THEN
3064       WRITE( message_string, * )  'section_yz must be <= nx + 1 = ', nx + 1
3065       CALL message( 'check_parameters', 'PA0321', 1, 2, 0, 6, 0 )
3066    ENDIF
3067    section(:,1) = section_xy
3068    section(:,2) = section_xz
3069    section(:,3) = section_yz
3070
3071!
3072!-- Upper plot limit for 2D vertical sections
3073    IF ( z_max_do2d == -1.0 )  z_max_do2d = zu(nzt)
3074    IF ( z_max_do2d < zu(nzb+1)  .OR.  z_max_do2d > zu(nzt) )  THEN
3075       WRITE( message_string, * )  'z_max_do2d = ', z_max_do2d, &
3076                    ' must be >= ', zu(nzb+1), '(zu(nzb+1)) and <= ', zu(nzt), &
3077                    ' (zu(nzt))'
3078       CALL message( 'check_parameters', 'PA0116', 1, 2, 0, 6, 0 )
3079    ENDIF
3080
3081!
3082!-- Upper plot limit for 3D arrays
3083    IF ( nz_do3d == -9999 )  nz_do3d = nzt + 1
3084
3085!
3086!-- Determine and check accuracy for compressed 3D plot output
3087    IF ( do3d_compress )  THEN
3088!
3089!--    Compression only permissible on T3E machines
3090       IF ( host(1:3) /= 't3e' )  THEN
3091          message_string = 'do3d_compress = .TRUE. not allowed on host "' // &
3092                           TRIM( host ) // '"'
3093          CALL message( 'check_parameters', 'PA0117', 1, 2, 0, 6, 0 )
3094       ENDIF
3095
3096       i = 1
3097       DO  WHILE ( do3d_comp_prec(i) /= ' ' )
3098
3099          ilen = LEN_TRIM( do3d_comp_prec(i) )
3100          IF ( LLT( do3d_comp_prec(i)(ilen:ilen), '0' ) .OR. &
3101               LGT( do3d_comp_prec(i)(ilen:ilen), '9' ) )  THEN
3102             WRITE( message_string, * )  'illegal precision: do3d_comp_prec', &
3103                                   '(', i, ') = "', TRIM(do3d_comp_prec(i)),'"'
3104             CALL message( 'check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
3105          ENDIF
3106
3107          prec = IACHAR( do3d_comp_prec(i)(ilen:ilen) ) - IACHAR( '0' )
3108          var = do3d_comp_prec(i)(1:ilen-1)
3109
3110          SELECT CASE ( var )
3111
3112             CASE ( 'u' )
3113                j = 1
3114             CASE ( 'v' )
3115                j = 2
3116             CASE ( 'w' )
3117                j = 3
3118             CASE ( 'p' )
3119                j = 4
3120             CASE ( 'pt' )
3121                j = 5
3122
3123             CASE DEFAULT
3124                WRITE( message_string, * )  'unknown variable "', &
3125                     TRIM( do3d_comp_prec(i) ), '" given for do3d_comp_prec(', &
3126                     i, ')'
3127                CALL message( 'check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
3128
3129          END SELECT
3130
3131          plot_3d_precision(j)%precision = prec
3132          i = i + 1
3133
3134       ENDDO
3135    ENDIF
3136
3137!
3138!-- Check the data output format(s)
3139    IF ( data_output_format(1) == ' ' )  THEN
3140!
3141!--    Default value
3142       netcdf_output = .TRUE.
3143    ELSE
3144       i = 1
3145       DO  WHILE ( data_output_format(i) /= ' ' )
3146
3147          SELECT CASE ( data_output_format(i) )
3148
3149             CASE ( 'netcdf' )
3150                netcdf_output = .TRUE.
3151             CASE ( 'iso2d' )
3152                iso2d_output  = .TRUE.
3153             CASE ( 'avs' )
3154                avs_output    = .TRUE.
3155
3156             CASE DEFAULT
3157                message_string = 'unknown value for data_output_format "' // &
3158                                 TRIM( data_output_format(i) ) // '"'
3159                CALL message( 'check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
3160
3161          END SELECT
3162
3163          i = i + 1
3164          IF ( i > 10 )  EXIT
3165
3166       ENDDO
3167    ENDIF
3168
3169!
3170!-- Set output format string (used in header)
3171    IF ( netcdf_output )  THEN
3172
3173       SELECT CASE ( netcdf_data_format )
3174          CASE ( 1 )
3175             output_format_netcdf = 'netCDF classic'
3176          CASE ( 2 )
3177             output_format_netcdf = 'netCDF 64bit offset'
3178          CASE ( 3 )
3179             output_format_netcdf = 'netCDF4/HDF5'
3180          CASE ( 4 )
3181             output_format_netcdf = 'netCDF4/HDF5 classic'
3182          CASE ( 5 )
3183             output_format_netcdf = 'parallel netCDF4/HDF5'
3184          CASE ( 6 )
3185             output_format_netcdf = 'parallel netCDF4/HDF5 classic'
3186
3187       END SELECT
3188
3189    ENDIF
3190
3191#if defined( __spectra )
3192!
3193!-- Check the number of spectra level to be output
3194    i = 1
3195    DO WHILE ( comp_spectra_level(i) /= 999999  .AND.  i <= 100 )
3196       i = i + 1
3197    ENDDO
3198    i = i - 1
3199    IF ( i == 0 )  THEN
3200       WRITE( message_string, * )  'no spectra levels given'
3201       CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
3202    ENDIF
3203#endif
3204
3205!
3206!-- Check mask conditions
3207    DO mid = 1, max_masks
3208       IF ( data_output_masks(mid,1) /= ' ' .OR.   &
3209            data_output_masks_user(mid,1) /= ' ' ) THEN
3210          masks = masks + 1
3211       ENDIF
3212    ENDDO
3213   
3214    IF ( masks < 0 .OR. masks > max_masks )  THEN
3215       WRITE( message_string, * )  'illegal value: masks must be >= 0 and ', &
3216            '<= ', max_masks, ' (=max_masks)'
3217       CALL message( 'check_parameters', 'PA0325', 1, 2, 0, 6, 0 )
3218    ENDIF
3219    IF ( masks > 0 )  THEN
3220       mask_scale(1) = mask_scale_x
3221       mask_scale(2) = mask_scale_y
3222       mask_scale(3) = mask_scale_z
3223       IF ( ANY( mask_scale <= 0.0 ) )  THEN
3224          WRITE( message_string, * )  &
3225               'illegal value: mask_scale_x, mask_scale_y and mask_scale_z', &
3226               'must be > 0.0'
3227          CALL message( 'check_parameters', 'PA0326', 1, 2, 0, 6, 0 )
3228       ENDIF
3229!
3230!--    Generate masks for masked data output
3231       CALL init_masks
3232    ENDIF
3233
3234!
3235!-- Check the NetCDF data format
3236#if ! defined ( __check )
3237    IF ( netcdf_data_format > 2 )  THEN
3238#if defined( __netcdf4 )
3239       CONTINUE
3240#else
3241       message_string = 'netCDF: netCDF4 format requested but no ' // &
3242                        'cpp-directive __netcdf4 given & switch '  // &
3243                        'back to 64-bit offset format'
3244       CALL message( 'check_parameters', 'PA0171', 0, 1, 0, 6, 0 )
3245       netcdf_data_format = 2
3246#endif
3247    ENDIF
3248    IF ( netcdf_data_format > 4 )  THEN
3249#if defined( __netcdf4 ) && defined( __netcdf4_parallel )
3250       CONTINUE
3251#else
3252       message_string = 'netCDF: netCDF4 parallel output requested but no ' // &
3253                        'cpp-directive __netcdf4_parallel given & switch '  // &
3254                        'back to netCDF4 non-parallel output'
3255       CALL message( 'check_parameters', 'PA0099', 0, 1, 0, 6, 0 )
3256       netcdf_data_format = netcdf_data_format - 2
3257#endif
3258    ENDIF
3259#endif
3260
3261#if ! defined( __check )
3262!
3263!-- Check netcdf precison
3264    ldum = .FALSE.
3265    CALL define_netcdf_header( 'ch', ldum, 0 )
3266#endif
3267!
3268!-- Check, whether a constant diffusion coefficient shall be used
3269    IF ( km_constant /= -1.0 )  THEN
3270       IF ( km_constant < 0.0 )  THEN
3271          WRITE( message_string, * )  'km_constant = ', km_constant, ' < 0.0'
3272          CALL message( 'check_parameters', 'PA0121', 1, 2, 0, 6, 0 )
3273       ELSE
3274          IF ( prandtl_number < 0.0 )  THEN
3275             WRITE( message_string, * )  'prandtl_number = ', prandtl_number, &
3276                                         ' < 0.0'
3277             CALL message( 'check_parameters', 'PA0122', 1, 2, 0, 6, 0 )
3278          ENDIF
3279          constant_diffusion = .TRUE.
3280
3281          IF ( prandtl_layer )  THEN
3282             message_string = 'prandtl_layer is not allowed with fixed ' // &
3283                              'value of km'
3284             CALL message( 'check_parameters', 'PA0123', 1, 2, 0, 6, 0 )
3285          ENDIF
3286       ENDIF
3287    ENDIF
3288
3289!
3290!-- In case of non-cyclic lateral boundaries and a damping layer for the
3291!-- potential temperature, check the width of the damping layer
3292    IF ( bc_lr /= 'cyclic' ) THEN
3293       IF ( pt_damping_width < 0.0  .OR.  pt_damping_width > REAL( nx * dx ) )  THEN
3294          message_string = 'pt_damping_width out of range'
3295          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3296       ENDIF
3297    ENDIF
3298
3299    IF ( bc_ns /= 'cyclic' )  THEN
3300       IF ( pt_damping_width < 0.0  .OR.  pt_damping_width > REAL( ny * dy ) )  THEN
3301          message_string = 'pt_damping_width out of range'
3302          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3303       ENDIF
3304    ENDIF
3305
3306!
3307!-- Check value range for rif
3308    IF ( rif_min >= rif_max )  THEN
3309       WRITE( message_string, * )  'rif_min = ', rif_min, ' must be less ', &
3310                                   'than rif_max = ', rif_max
3311       CALL message( 'check_parameters', 'PA0125', 1, 2, 0, 6, 0 )
3312    ENDIF
3313
3314!
3315!-- Determine upper and lower hight level indices for random perturbations
3316    IF ( disturbance_level_b == -9999999.9 )  THEN
3317       IF ( ocean ) THEN
3318          disturbance_level_b     = zu((nzt*2)/3)
3319          disturbance_level_ind_b = ( nzt * 2 ) / 3
3320       ELSE
3321          disturbance_level_b     = zu(nzb+3)
3322          disturbance_level_ind_b = nzb + 3
3323       ENDIF
3324    ELSEIF ( disturbance_level_b < zu(3) )  THEN
3325       WRITE( message_string, * )  'disturbance_level_b = ', &
3326                           disturbance_level_b, ' must be >= ', zu(3), '(zu(3))'
3327       CALL message( 'check_parameters', 'PA0126', 1, 2, 0, 6, 0 )
3328    ELSEIF ( disturbance_level_b > zu(nzt-2) )  THEN
3329       WRITE( message_string, * )  'disturbance_level_b = ', &
3330                   disturbance_level_b, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3331       CALL message( 'check_parameters', 'PA0127', 1, 2, 0, 6, 0 )
3332    ELSE
3333       DO  k = 3, nzt-2
3334          IF ( disturbance_level_b <= zu(k) )  THEN
3335             disturbance_level_ind_b = k
3336             EXIT
3337          ENDIF
3338       ENDDO
3339    ENDIF
3340
3341    IF ( disturbance_level_t == -9999999.9 )  THEN
3342       IF ( ocean )  THEN
3343          disturbance_level_t     = zu(nzt-3)
3344          disturbance_level_ind_t = nzt - 3
3345       ELSE
3346          disturbance_level_t     = zu(nzt/3)
3347          disturbance_level_ind_t = nzt / 3
3348       ENDIF
3349    ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
3350       WRITE( message_string, * )  'disturbance_level_t = ', &
3351                   disturbance_level_t, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3352       CALL message( 'check_parameters', 'PA0128', 1, 2, 0, 6, 0 )
3353    ELSEIF ( disturbance_level_t < disturbance_level_b )  THEN
3354       WRITE( message_string, * )  'disturbance_level_t = ', &
3355                   disturbance_level_t, ' must be >= disturbance_level_b = ', &
3356                   disturbance_level_b
3357       CALL message( 'check_parameters', 'PA0129', 1, 2, 0, 6, 0 )
3358    ELSE
3359       DO  k = 3, nzt-2
3360          IF ( disturbance_level_t <= zu(k) )  THEN
3361             disturbance_level_ind_t = k
3362             EXIT
3363          ENDIF
3364       ENDDO
3365    ENDIF
3366
3367!
3368!-- Check again whether the levels determined this way are ok.
3369!-- Error may occur at automatic determination and too few grid points in
3370!-- z-direction.
3371    IF ( disturbance_level_ind_t < disturbance_level_ind_b )  THEN
3372       WRITE( message_string, * )  'disturbance_level_ind_t = ', &
3373                disturbance_level_ind_t, ' must be >= disturbance_level_b = ', &
3374                disturbance_level_b
3375       CALL message( 'check_parameters', 'PA0130', 1, 2, 0, 6, 0 )
3376    ENDIF
3377
3378!
3379!-- Determine the horizontal index range for random perturbations.
3380!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed
3381!-- near the inflow and the perturbation area is further limited to ...(1)
3382!-- after the initial phase of the flow.
3383    dist_nxl = 0;  dist_nxr = nx
3384    dist_nys = 0;  dist_nyn = ny
3385    IF ( bc_lr /= 'cyclic' )  THEN
3386       IF ( inflow_disturbance_begin == -1 )  THEN
3387          inflow_disturbance_begin = MIN( 10, nx/2 )
3388       ENDIF
3389       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > nx )&
3390       THEN
3391          message_string = 'inflow_disturbance_begin out of range'
3392          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
3393       ENDIF
3394       IF ( inflow_disturbance_end == -1 )  THEN
3395          inflow_disturbance_end = MIN( 100, 3*nx/4 )
3396       ENDIF
3397       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > nx )    &
3398       THEN
3399          message_string = 'inflow_disturbance_end out of range'
3400          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
3401       ENDIF
3402    ELSEIF ( bc_ns /= 'cyclic' )  THEN
3403       IF ( inflow_disturbance_begin == -1 )  THEN
3404          inflow_disturbance_begin = MIN( 10, ny/2 )
3405       ENDIF
3406       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > ny )&
3407       THEN
3408          message_string = 'inflow_disturbance_begin out of range'
3409          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
3410       ENDIF
3411       IF ( inflow_disturbance_end == -1 )  THEN
3412          inflow_disturbance_end = MIN( 100, 3*ny/4 )
3413       ENDIF
3414       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > ny )    &
3415       THEN
3416          message_string = 'inflow_disturbance_end out of range'
3417          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
3418       ENDIF
3419    ENDIF
3420
3421    IF ( bc_lr == 'radiation/dirichlet' )  THEN
3422       dist_nxr    = nx - inflow_disturbance_begin
3423       dist_nxl(1) = nx - inflow_disturbance_end
3424    ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
3425       dist_nxl    = inflow_disturbance_begin
3426       dist_nxr(1) = inflow_disturbance_end
3427    ENDIF
3428    IF ( bc_ns == 'dirichlet/radiation' )  THEN
3429       dist_nyn    = ny - inflow_disturbance_begin
3430       dist_nys(1) = ny - inflow_disturbance_end
3431    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
3432       dist_nys    = inflow_disturbance_begin
3433       dist_nyn(1) = inflow_disturbance_end
3434    ENDIF
3435
3436!
3437!-- A turbulent inflow requires Dirichlet conditions at the respective inflow
3438!-- boundary (so far, a turbulent inflow is realized from the left side only)
3439    IF ( turbulent_inflow  .AND.  bc_lr /= 'dirichlet/radiation' )  THEN
3440       message_string = 'turbulent_inflow = .T. requires a Dirichlet ' // &
3441                        'condition at the inflow boundary'
3442       CALL message( 'check_parameters', 'PA0133', 1, 2, 0, 6, 0 )
3443    ENDIF
3444
3445!
3446!-- Turbulent inflow requires that 3d arrays have been cyclically filled with
3447!-- data from prerun in the first main run
3448    IF ( turbulent_inflow  .AND.  initializing_actions /= 'cyclic_fill'  .AND. &
3449         initializing_actions /= 'read_restart_data' )  THEN
3450       message_string = 'turbulent_inflow = .T. requires ' // &
3451                        'initializing_actions = ''cyclic_fill'' '
3452       CALL message( 'check_parameters', 'PA0055', 1, 2, 0, 6, 0 )
3453    ENDIF
3454
3455!
3456!-- In case of turbulent inflow calculate the index of the recycling plane
3457    IF ( turbulent_inflow )  THEN
3458       IF ( recycling_width == 9999999.9 )  THEN
3459!
3460!--       Set the default value for the width of the recycling domain
3461          recycling_width = 0.1 * nx * dx
3462       ELSE
3463          IF ( recycling_width < dx  .OR.  recycling_width > nx * dx )  THEN
3464             WRITE( message_string, * )  'illegal value for recycling_width:', &
3465                                         ' ', recycling_width
3466             CALL message( 'check_parameters', 'PA0134', 1, 2, 0, 6, 0 )
3467          ENDIF
3468       ENDIF
3469!
3470!--    Calculate the index
3471       recycling_plane = recycling_width / dx
3472    ENDIF
3473
3474!
3475!-- Check random generator
3476    IF ( random_generator /= 'system-specific'  .AND. &
3477         random_generator /= 'numerical-recipes' )  THEN
3478       message_string = 'unknown random generator: random_generator = "' // &
3479                        TRIM( random_generator ) // '"'
3480       CALL message( 'check_parameters', 'PA0135', 1, 2, 0, 6, 0 )
3481    ENDIF
3482
3483!
3484!-- Determine damping level index for 1D model
3485    IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
3486       IF ( damp_level_1d == -1.0 )  THEN
3487          damp_level_1d     = zu(nzt+1)
3488          damp_level_ind_1d = nzt + 1
3489       ELSEIF ( damp_level_1d < 0.0  .OR.  damp_level_1d > zu(nzt+1) )  THEN
3490          WRITE( message_string, * )  'damp_level_1d = ', damp_level_1d, &
3491                 ' must be > 0.0 and < ', zu(nzt+1), '(zu(nzt+1))'
3492          CALL message( 'check_parameters', 'PA0136', 1, 2, 0, 6, 0 )
3493       ELSE
3494          DO  k = 1, nzt+1
3495             IF ( damp_level_1d <= zu(k) )  THEN
3496                damp_level_ind_1d = k
3497                EXIT
3498             ENDIF
3499          ENDDO
3500       ENDIF
3501    ENDIF
3502
3503!
3504!-- Check some other 1d-model parameters
3505    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND. &
3506         TRIM( mixing_length_1d ) /= 'blackadar' )  THEN
3507       message_string = 'mixing_length_1d = "' // TRIM( mixing_length_1d ) // &
3508                        '" is unknown'
3509       CALL message( 'check_parameters', 'PA0137', 1, 2, 0, 6, 0 )
3510    ENDIF
3511    IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model'  .AND. &
3512         TRIM( dissipation_1d ) /= 'detering' )  THEN
3513       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) // &
3514                        '" is unknown'
3515       CALL message( 'check_parameters', 'PA0138', 1, 2, 0, 6, 0 )
3516    ENDIF
3517
3518!
3519!-- Set time for the next user defined restart (time_restart is the
3520!-- internal parameter for steering restart events)
3521    IF ( restart_time /= 9999999.9 )  THEN
3522       IF ( restart_time > time_since_reference_point )  THEN
3523          time_restart = restart_time
3524       ENDIF
3525    ELSE
3526!
3527!--    In case of a restart run, set internal parameter to default (no restart)
3528!--    if the NAMELIST-parameter restart_time is omitted
3529       time_restart = 9999999.9
3530    ENDIF
3531
3532!
3533!-- Set default value of the time needed to terminate a model run
3534    IF ( termination_time_needed == -1.0 )  THEN
3535       IF ( host(1:3) == 'ibm' )  THEN
3536          termination_time_needed = 300.0
3537       ELSE
3538          termination_time_needed = 35.0
3539       ENDIF
3540    ENDIF
3541
3542!
3543!-- Check the time needed to terminate a model run
3544    IF ( host(1:3) == 't3e' )  THEN
3545!
3546!--    Time needed must be at least 30 seconds on all CRAY machines, because
3547!--    MPP_TREMAIN gives the remaining CPU time only in steps of 30 seconds
3548       IF ( termination_time_needed <= 30.0 )  THEN
3549          WRITE( message_string, * )  'termination_time_needed = ', &
3550                 termination_time_needed, ' must be > 30.0 on host "', &
3551                 TRIM( host ), '"'
3552          CALL message( 'check_parameters', 'PA0139', 1, 2, 0, 6, 0 )
3553       ENDIF
3554    ELSEIF ( host(1:3) == 'ibm' )  THEN
3555!
3556!--    On IBM-regatta machines the time should be at least 300 seconds,
3557!--    because the job time consumed before executing palm (for compiling,
3558!--    copying of files, etc.) has to be regarded
3559       IF ( termination_time_needed < 300.0 )  THEN
3560          WRITE( message_string, * )  'termination_time_needed = ', &
3561                 termination_time_needed, ' should be >= 300.0 on host "', &
3562                 TRIM( host ), '"'
3563          CALL message( 'check_parameters', 'PA0140', 1, 2, 0, 6, 0 )
3564       ENDIF
3565    ENDIF
3566
3567!
3568!-- Check pressure gradient conditions
3569    IF ( dp_external .AND. conserve_volume_flow )  THEN
3570       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo', &
3571            'w are .TRUE. but one of them must be .FALSE.'
3572       CALL message( 'check_parameters', 'PA0150', 1, 2, 0, 6, 0 )
3573    ENDIF
3574    IF ( dp_external )  THEN
3575       IF ( dp_level_b < zu(nzb) .OR. dp_level_b > zu(nzt) )  THEN
3576          WRITE( message_string, * )  'dp_level_b = ', dp_level_b, ' is out ', &
3577               ' of range'
3578          CALL message( 'check_parameters', 'PA0151', 1, 2, 0, 6, 0 )
3579       ENDIF
3580       IF ( .NOT. ANY( dpdxy /= 0.0 ) )  THEN
3581          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is ze', &
3582               'ro, i.e. the external pressure gradient & will not be applied'
3583          CALL message( 'check_parameters', 'PA0152', 0, 1, 0, 6, 0 )
3584       ENDIF
3585    ENDIF
3586    IF ( ANY( dpdxy /= 0.0 ) .AND. .NOT. dp_external )  THEN
3587       WRITE( message_string, * )  'dpdxy is nonzero but dp_external is ', &
3588            '.FALSE., i.e. the external pressure gradient & will not be applied'
3589       CALL message( 'check_parameters', 'PA0153', 0, 1, 0, 6, 0 )
3590    ENDIF
3591    IF ( conserve_volume_flow )  THEN
3592       IF ( TRIM( conserve_volume_flow_mode ) == 'default' )  THEN
3593
3594          conserve_volume_flow_mode = 'initial_profiles'
3595
3596       ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles' .AND.  &
3597            TRIM( conserve_volume_flow_mode ) /= 'inflow_profile' .AND.  &
3598            TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' )  THEN
3599          WRITE( message_string, * )  'unknown conserve_volume_flow_mode: ', &
3600               conserve_volume_flow_mode
3601          CALL message( 'check_parameters', 'PA0154', 1, 2, 0, 6, 0 )
3602       ENDIF
3603       IF ( (bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic')  .AND. &
3604          TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
3605          WRITE( message_string, * )  'non-cyclic boundary conditions ', &
3606               'require  conserve_volume_flow_mode = ''initial_profiles'''
3607          CALL message( 'check_parameters', 'PA0155', 1, 2, 0, 6, 0 )
3608       ENDIF
3609       IF ( bc_lr == 'cyclic'  .AND.  bc_ns == 'cyclic'  .AND.  &
3610            TRIM( conserve_volume_flow_mode ) == 'inflow_profile' )  THEN
3611          WRITE( message_string, * )  'cyclic boundary conditions ', &
3612               'require conserve_volume_flow_mode = ''initial_profiles''', &
3613               ' or ''bulk_velocity'''
3614          CALL message( 'check_parameters', 'PA0156', 1, 2, 0, 6, 0 )
3615       ENDIF
3616    ENDIF
3617    IF ( ( u_bulk /= 0.0 .OR. v_bulk /= 0.0 ) .AND.  &
3618         ( .NOT. conserve_volume_flow .OR.  &
3619         TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) )  THEN
3620       WRITE( message_string, * )  'nonzero bulk velocity requires ', &
3621            'conserve_volume_flow = .T. and ', &
3622            'conserve_volume_flow_mode = ''bulk_velocity'''
3623       CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 )
3624    ENDIF
3625
3626!
3627!-- Check particle attributes
3628    IF ( particle_color /= 'none' )  THEN
3629       IF ( particle_color /= 'absuv'  .AND.  particle_color /= 'pt*'  .AND.  &
3630            particle_color /= 'z' )  THEN
3631          message_string = 'illegal value for parameter particle_color: ' // &
3632                           TRIM( particle_color)
3633          CALL message( 'check_parameters', 'PA0313', 1, 2, 0, 6, 0 )
3634       ELSE
3635          IF ( color_interval(2) <= color_interval(1) )  THEN
3636             message_string = 'color_interval(2) <= color_interval(1)'
3637             CALL message( 'check_parameters', 'PA0315', 1, 2, 0, 6, 0 )
3638          ENDIF
3639       ENDIF
3640    ENDIF
3641
3642    IF ( particle_dvrpsize /= 'none' )  THEN
3643       IF ( particle_dvrpsize /= 'absw' )  THEN
3644          message_string = 'illegal value for parameter particle_dvrpsize:' // &
3645                           ' ' // TRIM( particle_color)
3646          CALL message( 'check_parameters', 'PA0314', 1, 2, 0, 6, 0 )
3647       ELSE
3648          IF ( dvrpsize_interval(2) <= dvrpsize_interval(1) )  THEN
3649             message_string = 'dvrpsize_interval(2) <= dvrpsize_interval(1)'
3650             CALL message( 'check_parameters', 'PA0316', 1, 2, 0, 6, 0 )
3651          ENDIF
3652       ENDIF
3653    ENDIF
3654
3655!
3656!-- Check nudging and large scale forcing from external file
3657    IF ( nudging .AND. ( .NOT. large_scale_forcing ) )  THEN
3658       message_string = 'Nudging requires large_scale_forcing = .T.. &'// &
3659                        'Surface fluxes and geostrophic wind should be &'// &
3660                        'prescribed in file LSF_DATA'
3661       CALL message( 'check_parameters', 'PA0374', 1, 2, 0, 6, 0 )
3662    ENDIF
3663
3664    IF ( large_scale_forcing .AND. ( bc_lr /= 'cyclic'  .OR. &
3665                                    bc_ns /= 'cyclic' ) )  THEN
3666       message_string = 'Non-cyclic lateral boundaries do not allow for &' // &
3667                        'the usage of large scale forcing from external file.'
3668       CALL message( 'check_parameters', 'PA0375', 1, 2, 0, 6, 0 )
3669     ENDIF
3670
3671    IF ( large_scale_forcing .AND. ( .NOT. humidity ) )  THEN
3672       message_string = 'The usage of large scale forcing from external &'//& 
3673                        'file LSF_DATA requires humidity = .T..'
3674       CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
3675     ENDIF
3676
3677    IF ( large_scale_forcing .AND. topography /= 'flat' )  THEN
3678       message_string = 'The usage of large scale forcing from external &'//& 
3679                        'file LSF_DATA is not implemented for non-flat topography'
3680       CALL message( 'check_parameters', 'PA0377', 1, 2, 0, 6, 0 )
3681    ENDIF
3682
3683    IF ( large_scale_forcing .AND.  ocean  )  THEN
3684       message_string = 'The usage of large scale forcing from external &'//& 
3685                        'file LSF_DATA is not implemented for ocean runs'
3686       CALL message( 'check_parameters', 'PA0378', 1, 2, 0, 6, 0 )
3687    ENDIF
3688!
3689!-- Check &userpar parameters
3690    CALL user_check_parameters
3691
3692
3693
3694 END SUBROUTINE check_parameters
Note: See TracBrowser for help on using the repository browser.