source: palm/tags/release-4.0/SOURCE/check_parameters.f90 @ 2704

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

last commit documented

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