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

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