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

Last change on this file since 1586 was 1586, checked in by maronga, 10 years ago

last commit documented

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