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

Last change on this file since 2001 was 2001, checked in by knoop, 8 years ago

last commit documented

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