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

Last change on this file since 1917 was 1917, checked in by witha, 8 years ago

Corrected documented changes in the file headers

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