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

Last change on this file since 1576 was 1576, checked in by raasch, 7 years ago

last commit documented

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