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

Last change on this file since 1826 was 1826, checked in by maronga, 8 years ago

further modularization of radiation model and plant canopy model

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