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

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