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

Last change on this file since 1971 was 1971, checked in by suehring, 8 years ago

last commit documented

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