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

Last change on this file since 1823 was 1823, checked in by hoffmann, 8 years ago

last commit documented

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