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

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