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

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

last commit documented

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