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

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

last commit documented / copyright update

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