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

Last change on this file since 2008 was 2008, checked in by kanani, 5 years ago

last commit documented

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