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

Last change on this file since 1330 was 1330, checked in by suehring, 11 years ago

In case of SGS-particle velocity advection of TKE is also allowed with 5th-order scheme.

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