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

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

last commit documented

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