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

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

small adjustments for RRTMG

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