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

Last change on this file since 1925 was 1919, checked in by raasch, 8 years ago

last commit documented

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