source: palm/trunk/SOURCE/check_parameters.f90 @ 1276

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

Usage of Dirichlet bottom boundary condition for scalars in conjunction with large scale forcing enabled

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