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

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