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

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