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

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