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

Last change on this file since 1833 was 1833, checked in by raasch, 5 years ago

spectrum renamed spactra_par and further modularized, POINTER-attributes added in coupler-routines to avoid gfortran error messages

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