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

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

output of location messages to stdout

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