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

Last change on this file since 1764 was 1764, checked in by raasch, 8 years ago

update of the nested domain system + some bugfixes

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