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

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

Separate balance equations for humidity and passive_scalar

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