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

Last change on this file since 1553 was 1553, checked in by maronga, 9 years ago

further adjustments in land surface model

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