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

Last change on this file since 1827 was 1827, checked in by maronga, 8 years ago

last commit documented

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