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

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

last commit documented

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