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

Last change on this file since 4853 was 4853, checked in by raasch, 3 years ago

array sizes for output profiles increased from 300 to 400

  • Property svn:keywords set to Id
  • Property svn:mergeinfo set to (toggle deleted branches)
    /palm/branches/chemistry/SOURCE/check_parameters.f902047-3190,​3218-3297
    /palm/branches/rans/SOURCE/check_parameters.f902078-3128
    /palm/branches/resler/SOURCE/check_parameters.f902023-4492
    /palm/branches/salsa/SOURCE/check_parameters.f902503-3581
    /palm/branches/forwind/SOURCE/check_parameters.f901564-1913
    /palm/branches/mosaik_M2/check_parameters.f902360-3471
    /palm/branches/palm4u/SOURCE/check_parameters.f902540-2692
File size: 143.0 KB
Line 
1!> @file check_parameters.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2021 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: check_parameters.f90 4853 2021-01-25 07:56:09Z raasch $
26! number of allowed profiles increased from 300 to 400
27!
28! 4844 2021-01-18 09:22:26Z raasch
29! size check for profile output added
30!
31! 4828 2021-01-05 11:21:41Z Giersch
32! allow prescribing initial uv_profiles AND geostrophic wind
33! by changing error message PA0347 to an informative message
34!
35! 4685 2020-09-18 11:49:50Z gronemeier
36! filter data_output_pr and data_output list for duplicate entries
37!
38! 4633 2020-08-05 14:21:14Z suehring
39! todo added
40!
41! 4565 2020-06-15 08:30:38Z oliver.maas
42! added check for pt_surface_heating_rate
43!
44! 4564 2020-06-12 14:03:36Z raasch
45! Vertical nesting method of Huq et al. (2019) removed
46!
47! 4562 2020-06-12 08:38:47Z raasch
48! bugfix: revised error message for exceeding allow number of time series
49!
50! 4559 2020-06-11 08:51:48Z raasch
51! file re-formatted to follow the PALM coding standard
52!
53! 4536 2020-05-17 17:24:13Z raasch
54! unneccessary query for restart data format removed
55!
56! 4534 2020-05-14 18:35:22Z raasch
57! adjustments for I/O on reduced number of cores using shared memory MPI
58!
59! 4514 2020-04-30 16:29:59Z suehring
60! Enable output of qsurf and ssurf
61!
62! 4513 2020-04-30 13:45:47Z raasch
63! unused modules removed
64!
65! 4511 2020-04-30 12:20:40Z raasch
66! call of chem_boundary_conds removed (respective settings are now done in the chemistry module)
67!
68! 4495 2020-04-13 20:11:20Z raasch
69! check new restart_data_format parameters
70!
71! 4493 2020-04-10 09:49:43Z pavelkrc
72! Kolmogorov length scale eta added to profile output
73!
74! 4444 2020-03-05 15:59:50Z raasch
75! bugfix: cpp-directives for serial mode added
76!
77! 4392 2020-01-31 16:14:57Z pavelkrc
78! Some error numbers revised to prevent double usage
79!
80! 11:55:33Z oliver.maas
81! Checks for closed channel flow implemented
82!
83! 11:55:33Z oliver.maas
84! Move 2-m potential temperature output to diagnostic_output_quantities
85!
86! 11:55:33Z oliver.maas
87! removed message PA0421, concerning old parameter recycling_yshift
88!
89! 11:55:33Z oliver.maas
90! adjust message to the modified parameter recycling_yshift
91!
92! 11:55:33Z oliver.maas
93! Check if a cross section is specified if any output cross-section quantity
94! is given
95!
96! 11:55:33Z oliver.maas
97! Overwrite rotation_angle from namelist by value from static driver
98!
99! 11:55:33Z oliver.maas
100! removed conversion from recycle_absolute_quantities to raq, added check and
101! error message for correct input of recycling_method_for_thermodynamic_quantities
102!
103! 11:55:33Z oliver.maas
104! Corrected "Former revisions" section
105!
106! 11:55:33Z oliver.maas
107! bugfix error message: replaced PA184 by PA0184
108!
109! 11:55:33Z oliver.maas
110! added conversion from recycle_absolute_quantities to raq for recycling of
111! absolute quantities and added error message PA184 for not implemented quantities
112!
113! 4142 2019-08-05 12:38:31Z suehring
114! Consider spinup in number of output timesteps for averaged 2D output (merge
115! from branch resler).
116!
117! 4069 2019-07-01 14:05:51Z Giersch
118! Masked output running index mid has been introduced as a local variable to
119! avoid runtime error (Loop variable has been modified) in time_integration
120!
121! 4048 2019-06-21 21:00:21Z knoop
122! Moved tcm_check_data_output to module_interface
123!
124! 4039 2019-06-18 10:32:41Z suehring
125! Modularize diagnostic output
126!
127! 4017 2019-06-06 12:16:46Z schwenkel
128! output of turbulence intensity added
129!
130! 3933 2019-04-25 12:33:20Z kanani
131! Alphabetical resorting in CASE, condense settings for theta_2m* into one IF clause
132!
133! 3885 2019-04-11 11:29:34Z kanani
134! Changes related to global restructuring of location messages and introduction
135! of additional debug messages
136!
137! 3766 2019-02-26 16:23:41Z raasch
138! trim added to avoid truncation compiler warnings
139!
140! 3761 2019-02-25 15:31:42Z raasch
141! unused variables removed
142!
143! 3735 2019-02-12 09:52:40Z dom_dwd_user
144! Passing variable j (averaged output?) to
145! module_interface.f90:chem_check_data_output.
146!
147! 3705 2019-01-29 19:56:39Z suehring
148! bugfix: renamed thetav_t to vtheta_t
149!
150! 3702 2019-01-28 13:19:30Z gronemeier
151! most_method removed
152!
153! 3655 2019-01-07 16:51:22Z knoop
154! Formatting
155!
156! Revision 1.1  1997/08/26 06:29:23  raasch
157! Initial revision
158!
159!
160! Description:
161! ------------
162!> Check control parameters and deduce further quantities.
163!
164!> @todo Increase character length of unit and corresponding characters to LEN>=8 in order to allow
165!>       units like degree_C (05.08.2020)
166!--------------------------------------------------------------------------------------------------!
167 SUBROUTINE check_parameters
168
169
170    USE arrays_3d
171
172    USE basic_constants_and_equations_mod
173
174    USE bulk_cloud_model_mod,                                                                      &
175        ONLY:  bulk_cloud_model
176
177    USE control_parameters
178
179    USE grid_variables
180
181    USE kinds
182
183    USE indices
184
185    USE model_1d_mod,                                                                              &
186        ONLY:  damp_level_1d, damp_level_ind_1d
187
188    USE module_interface,                                                                          &
189        ONLY:  module_interface_check_data_output,                                                 &
190               module_interface_check_data_output_pr,                                              &
191               module_interface_check_data_output_ts,                                              &
192               module_interface_check_parameters
193
194    USE netcdf_data_input_mod,                                                                     &
195        ONLY:  init_model, input_pids_static, netcdf_data_input_check_dynamic,                     &
196               netcdf_data_input_check_static
197
198    USE netcdf_interface,                                                                          &
199        ONLY:  do2d_unit, do3d_unit, dopr_unit, dots_label, dots_max, dots_num, dots_unit,         &
200               heatflux_output_unit, momentumflux_output_unit, netcdf_data_format,                 &
201               netcdf_data_format_string, waterflux_output_unit
202
203    USE particle_attributes,                                                                       &
204        ONLY:  particle_advection, use_sgs_for_particles
205
206    USE pegrid
207
208    USE pmc_interface,                                                                             &
209        ONLY:  cpl_id, nested_run
210
211    USE profil_parameter
212
213    USE statistics
214
215    USE subsidence_mod
216
217    USE transpose_indices
218
219    IMPLICIT NONE
220
221    CHARACTER (LEN=varnamelength)  ::  var           !< variable name
222    CHARACTER (LEN=7)   ::  unit                     !< unit of variable
223    CHARACTER (LEN=8)   ::  date                     !< current date string
224    CHARACTER (LEN=10)  ::  time                     !< current time string
225    CHARACTER (LEN=20)  ::  ensemble_string          !< string containing number of ensemble member
226    CHARACTER (LEN=15)  ::  nest_string              !< string containing id of nested domain
227    CHARACTER (LEN=40)  ::  coupling_string          !< string containing type of coupling
228    CHARACTER (LEN=100) ::  action                   !< flag string
229
230    CHARACTER (LEN=varnamelength), DIMENSION(:), ALLOCATABLE ::  data_output_filtered  !< filtered list of output variables
231
232    INTEGER(iwp) ::  i                               !< loop index
233    INTEGER(iwp) ::  ilen                            !< string length
234    INTEGER(iwp) ::  j                               !< loop index
235    INTEGER(iwp) ::  k                               !< loop index
236    INTEGER(iwp) ::  kk                              !< loop index
237    INTEGER(iwp) ::  mid                             !< masked output running index
238    INTEGER(iwp) ::  netcdf_data_format_save         !< initial value of netcdf_data_format
239    INTEGER(iwp) ::  position                        !< index position of string
240
241    LOGICAL     ::  found                            !< flag, true if output variable is already marked for averaging
242
243    REAL(wp)    ::  gradient                         !< local gradient
244#if defined( __parallel )
245    REAL(wp)    ::  dt_spinup_max                    !< maximum spinup timestep in nested domains
246    REAL(wp)    ::  remote = 0.0_wp                  !< MPI id of remote processor
247    REAL(wp)    ::  spinup_time_max                  !< maximum spinup time in nested domains
248    REAL(wp)    ::  time_to_be_simulated_from_reference_point  !< time to be simulated from reference point
249#endif
250
251
252    CALL location_message( 'checking parameters', 'start' )
253!
254!-- At first, check static and dynamic input for consistency.
255    CALL netcdf_data_input_check_dynamic
256    CALL netcdf_data_input_check_static
257!
258!-- Check for overlap combinations, which are not realized yet
259    IF ( transpose_compute_overlap  .AND. numprocs == 1 )  THEN
260          message_string = 'transpose-compute-overlap not implemented for single PE runs'
261          CALL message( 'check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
262    ENDIF
263
264!
265!-- Check the coupling mode
266    IF ( coupling_mode /= 'uncoupled'            .AND.                                             &
267         coupling_mode /= 'precursor_atmos'      .AND.                                             &
268         coupling_mode /= 'precursor_ocean'      .AND.                                             &
269         coupling_mode /= 'atmosphere_to_ocean'  .AND.                                             &
270         coupling_mode /= 'ocean_to_atmosphere' )  THEN
271       message_string = 'illegal coupling mode: ' // TRIM( coupling_mode )
272       CALL message( 'check_parameters', 'PA0002', 1, 2, 0, 6, 0 )
273    ENDIF
274
275!
276!-- Check if humidity is set to .TRUE. in case of the atmospheric run (for coupled runs)
277    IF ( coupling_mode == 'atmosphere_to_ocean' .AND. .NOT. humidity) THEN
278       message_string = ' Humidity has to be set to .T. in the _p3d file ' //                      &
279                        'for coupled runs between ocean and atmosphere.'
280       CALL message( 'check_parameters', 'PA0476', 1, 2, 0, 6, 0 )
281    ENDIF
282
283!
284!-- Check and set the restart data format variables
285    IF ( TRIM( restart_data_format ) /= 'fortran_binary'  .AND.                                    &
286         TRIM( restart_data_format ) /= 'mpi'             .AND.                                    &
287         TRIM( restart_data_format ) /= 'mpi_shared_memory' )  THEN
288       message_string = 'illegal restart data format "' // TRIM( restart_data_format ) // '"'
289       CALL message( 'check_parameters', 'PA0724', 1, 2, 0, 6, 0 )
290    ENDIF
291
292    IF ( TRIM( restart_data_format_input ) == 'undefined' )  THEN
293       restart_data_format_input = restart_data_format
294    ENDIF
295    IF ( TRIM( restart_data_format_output ) == 'undefined' )  THEN
296       restart_data_format_output = restart_data_format
297    ENDIF
298
299    IF ( TRIM( restart_data_format_input ) /= 'fortran_binary'  .AND.                              &
300         TRIM( restart_data_format_input ) /= 'mpi'             .AND.                              &
301         TRIM( restart_data_format_input ) /= 'mpi_shared_memory' )  THEN
302       message_string = 'illegal restart input data format "' //                                   &
303                        TRIM( restart_data_format_input ) // '"'
304       CALL message( 'check_parameters', 'PA0725', 1, 2, 0, 6, 0 )
305    ENDIF
306    IF ( TRIM( restart_data_format_output ) /= 'fortran_binary'  .AND.                             &
307         TRIM( restart_data_format_output ) /= 'mpi'             .AND.                             &
308         TRIM( restart_data_format_output ) /= 'mpi_shared_memory' )  THEN
309       message_string = 'illegal restart output data format "' //                                  &
310                        TRIM( restart_data_format_output ) // '"'
311       CALL message( 'check_parameters', 'PA0726', 1, 2, 0, 6, 0 )
312    ENDIF
313
314!
315!-- Set flag for including total domain boundaries in the restart data (MPI-IO) in case of
316!-- non-cyclic boundary conditions
317    include_total_domain_boundaries = .NOT. ( bc_lr_cyc  .AND.  bc_ns_cyc )
318
319!
320!-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny
321    IF ( coupling_mode /= 'uncoupled'        .AND.                                                 &
322         coupling_mode /= 'precursor_atmos'  .AND.                                                 &
323         coupling_mode /= 'precursor_ocean' )  THEN
324
325       IF ( dt_coupling == 9999999.9_wp )  THEN
326          message_string = 'dt_coupling is not set but required for coupling mode "' //            &
327                           TRIM( coupling_mode ) // '"'
328          CALL message( 'check_parameters', 'PA0003', 1, 2, 0, 6, 0 )
329       ENDIF
330
331#if defined( __parallel )
332
333
334       IF ( myid == 0 ) THEN
335          CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, ierr )
336          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, status, ierr )
337       ENDIF
338       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
339
340       IF ( dt_coupling /= remote )  THEN
341          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), '": dt_coupling = ',&
342                 dt_coupling, '& is not equal to ', 'dt_coupling_remote = ', remote
343          CALL message( 'check_parameters', 'PA0004', 1, 2, 0, 6, 0 )
344       ENDIF
345       IF ( dt_coupling <= 0.0_wp )  THEN
346
347          IF ( myid == 0  ) THEN
348             CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
349             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, status, ierr )
350          ENDIF
351          CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
352
353          dt_coupling = MAX( dt_max, remote )
354          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                     &
355                 '": dt_coupling <= 0.0 & is not allowed and is reset to ', 'MAX(dt_max(A,O)) = ', &
356                 dt_coupling
357          CALL message( 'check_parameters', 'PA0005', 0, 1, 0, 6, 0 )
358       ENDIF
359
360       IF ( myid == 0 ) THEN
361          CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, ierr )
362          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, status, ierr )
363       ENDIF
364       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
365
366       IF ( restart_time /= remote )  THEN
367          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                     &
368                 '": restart_time = ', restart_time, '& is not equal to ',                         &
369                 'restart_time_remote = ', remote
370          CALL message( 'check_parameters', 'PA0006', 1, 2, 0, 6, 0 )
371       ENDIF
372
373       IF ( myid == 0 ) THEN
374          CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, ierr )
375          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, status, ierr )
376       ENDIF
377       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
378
379       IF ( dt_restart /= remote )  THEN
380          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), '": dt_restart = ', &
381                 dt_restart, '& is not equal to ', 'dt_restart_remote = ', remote
382          CALL message( 'check_parameters', 'PA0007', 1, 2, 0, 6, 0 )
383       ENDIF
384
385       time_to_be_simulated_from_reference_point = end_time-coupling_start_time
386
387       IF ( myid == 0 ) THEN
388          CALL MPI_SEND( time_to_be_simulated_from_reference_point, 1, MPI_REAL, target_id, 14,    &
389                         comm_inter, ierr )
390          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, status, ierr )
391       ENDIF
392       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
393
394       IF ( time_to_be_simulated_from_reference_point /= remote )  THEN
395          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                     &
396                 '": time_to_be_simulated_from_reference_point = ',                                &
397                 time_to_be_simulated_from_reference_point, '& is not equal ',                     &
398                 'to time_to_be_simulated_from_reference_point_remote = ',                         &
399                 remote
400          CALL message( 'check_parameters', 'PA0008', 1, 2, 0, 6, 0 )
401       ENDIF
402
403       IF ( myid == 0 ) THEN
404          CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
405          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, status, ierr )
406       ENDIF
407       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
408
409
410       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
411
412          IF ( dx < remote ) THEN
413             WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                  &
414                    '": dx in Atmosphere is not equal to or not larger than dx in ocean'
415             CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
416          ENDIF
417
418          IF ( (nx_a+1)*dx /= (nx_o+1)*remote )  THEN
419             WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                  &
420                    '": Domain size in x-direction is not equal in ocean and atmosphere'
421             CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
422          ENDIF
423
424       ENDIF
425
426       IF ( myid == 0) THEN
427          CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )
428          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, status, ierr )
429       ENDIF
430       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
431
432       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
433
434          IF ( dy < remote )  THEN
435             WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                  &
436                    '": dy in Atmosphere is not equal to or not larger than dy in ocean'
437             CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
438          ENDIF
439
440          IF ( (ny_a+1)*dy /= (ny_o+1)*remote )  THEN
441             WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                  &
442                    '": Domain size in y-direction is not equal in ocean and atmosphere'
443             CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )
444          ENDIF
445
446          IF ( MOD( nx_o+1, nx_a+1 ) /= 0 )  THEN
447             WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                  &
448                    '": nx+1 in ocean is not divisible by nx+1 in', ' atmosphere without remainder'
449             CALL message( 'check_parameters', 'PA0339', 1, 2, 0, 6, 0 )
450          ENDIF
451
452          IF ( MOD(ny_o+1,ny_a+1) /= 0 )  THEN
453             WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ),                  &
454                    '": ny+1 in ocean is not divisible by ny+1 in', ' atmosphere without remainder'
455
456             CALL message( 'check_parameters', 'PA0340', 1, 2, 0, 6, 0 )
457          ENDIF
458
459       ENDIF
460#else
461       WRITE( message_string, * ) 'coupling requires PALM to be compiled with',                    &
462              ' cpp-option "-D__parallel"'
463       CALL message( 'check_parameters', 'PA0141', 1, 2, 0, 6, 0 )
464#endif
465    ENDIF
466
467#if defined( __parallel )
468!
469!-- Exchange via intercommunicator
470    IF ( coupling_mode == 'atmosphere_to_ocean'  .AND.  myid == 0 )  THEN
471       CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter, ierr )
472    ELSEIF ( coupling_mode == 'ocean_to_atmosphere'  .AND.  myid == 0)  THEN
473       CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19, comm_inter, status, ierr )
474    ENDIF
475    CALL MPI_BCAST( humidity_remote, 1, MPI_LOGICAL, 0, comm2d, ierr)
476
477#endif
478
479!
480!-- User settings for restart times requires that "restart" has been given as file activation
481!-- string. Otherwise, binary output would not be saved by palmrun.
482    IF ( ( restart_time /= 9999999.9_wp  .OR.  dt_restart /= 9999999.9_wp )                        &
483         .AND.  .NOT. write_binary )  THEN
484       WRITE( message_string, * ) 'manual restart settings requires file ',                        &
485                                  'activation string "restart"'
486       CALL message( 'check_parameters', 'PA0001', 1, 2, 0, 6, 0 )
487    ENDIF
488
489
490!
491!-- Generate the file header which is used as a header for most of PALM's output files
492    CALL DATE_AND_TIME( date, time, run_zone )
493    run_date = date(1:4) // '-' // date(5:6) // '-' // date(7:8)
494    run_time = time(1:2) // ':' // time(3:4) // ':' // time(5:6)
495    IF ( coupling_mode == 'uncoupled' )  THEN
496       coupling_string = ''
497    ELSEIF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
498       coupling_string = ' coupled (atmosphere)'
499    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
500       coupling_string = ' coupled (ocean)'
501    ENDIF
502    IF ( ensemble_member_nr /= 0 )  THEN
503       WRITE( ensemble_string, '(2X,A,I2.2)' )  'en-no: ', ensemble_member_nr
504    ELSE
505       ensemble_string = ''
506    ENDIF
507    IF ( nested_run )  THEN
508       WRITE( nest_string, '(2X,A,I2.2)' )  'nest-id: ', cpl_id
509    ELSE
510       nest_string = ''
511    ENDIF
512
513    WRITE ( run_description_header, '(A,2X,A,2X,A,A,A,I2.2,A,A,A,2X,A,A,2X,A,1X,A)' )              &
514          TRIM( version ), TRIM( revision ), 'run: ', TRIM( run_identifier ), '.', runnr,          &
515          TRIM( coupling_string ), TRIM( nest_string ), TRIM( ensemble_string), 'host: ',          &
516          TRIM( host ), run_date, run_time
517
518!
519!-- Check the general loop optimization method
520    SELECT CASE ( TRIM( loop_optimization ) )
521
522       CASE ( 'cache', 'vector' )
523          CONTINUE
524
525       CASE DEFAULT
526          message_string = 'illegal value given for loop_optimization: "' //                       &
527                           TRIM( loop_optimization ) // '"'
528          CALL message( 'check_parameters', 'PA0013', 1, 2, 0, 6, 0 )
529
530    END SELECT
531
532!
533!-- Check topography setting (check for illegal parameter combinations)
534    IF ( topography /= 'flat' )  THEN
535       action = ' '
536       IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' )  THEN
537          WRITE( action, '(A,A)' )  'scalar_advec = ', scalar_advec
538       ENDIF
539       IF ( momentum_advec /= 'pw-scheme'  .AND.  momentum_advec /= 'ws-scheme' )  THEN
540          WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
541       ENDIF
542       IF ( psolver == 'sor' )  THEN
543          WRITE( action, '(A,A)' )  'psolver = ', psolver
544       ENDIF
545       IF ( sloping_surface )  THEN
546          WRITE( action, '(A)' )  'sloping surface = .TRUE.'
547       ENDIF
548       IF ( galilei_transformation )  THEN
549          WRITE( action, '(A)' )  'galilei_transformation = .TRUE.'
550       ENDIF
551       IF ( cloud_droplets )  THEN
552          WRITE( action, '(A)' )  'cloud_droplets = .TRUE.'
553       ENDIF
554       IF ( .NOT. constant_flux_layer  .AND.  topography /= 'closed_channel' )  THEN
555          WRITE( action, '(A)' )  'constant_flux_layer = .FALSE.'
556       ENDIF
557       IF ( action /= ' ' )  THEN
558          message_string = 'The specified topography does not allow ' // TRIM( action )
559          CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
560       ENDIF
561!
562!--    Check illegal/untested parameter combinations for closed channel
563       If ( topography == 'closed_channel' )  THEN
564          symmetry_flag = 1
565          message_string = 'Bottom and top boundary are treated equal'
566          CALL message( 'check_parameters', 'PA0699', 0, 0, 0, 6, 0 )
567
568          IF ( dz(1) /= dz(COUNT( dz /= -1.0_wp ))  .OR.  dz_stretch_level /= -9999999.9_wp)  THEN
569             WRITE( message_string, * )  'dz should be equal close to the ' //                     &
570                                         'boundaries due to symmetrical problem'
571             CALL message( 'check_parameters', 'PA0700', 1, 2, 0, 6, 0 )
572          ENDIF
573
574          IF ( constant_flux_layer )  THEN
575             WRITE( message_string, * )  'A constant flux layer is not ' //                        &
576                                         'allowed if a closed channel shall be used'
577             CALL message( 'check_parameters', 'PA0701', 1, 2, 0, 6, 0 )
578          ENDIF
579
580          IF ( ocean_mode )  THEN
581             WRITE( message_string, * )  'The ocean mode is not allowed if ' //                    &
582                                         'a closed channel shall be used'
583             CALL message( 'check_parameters', 'PA0702', 1, 2, 0, 6, 0 )
584          ENDIF
585
586          IF ( momentum_advec /= 'ws-scheme'  .OR.                                                 &
587               scalar_advec /= 'ws-scheme' )  THEN
588             WRITE( message_string, * )  'A closed channel require the upwind scheme of Wicker' // &
589                                         ' and Skamarock as the advection scheme'
590             CALL message( 'check_parameters', 'PA0703', 1, 2, 0, 6, 0 )
591          ENDIF
592       ENDIF
593    ENDIF
594
595!
596!-- Check approximation
597    IF ( TRIM( approximation ) /= 'boussinesq'  .AND.  TRIM( approximation ) /= 'anelastic' )  THEN
598       message_string = 'unknown approximation: approximation = "' // TRIM( approximation ) // '"'
599       CALL message( 'check_parameters', 'PA0446', 1, 2, 0, 6, 0 )
600    ENDIF
601
602!
603!-- Check approximation requirements
604    IF ( TRIM( approximation ) == 'anelastic'  .AND.  TRIM( momentum_advec ) /= 'ws-scheme' )  THEN
605       message_string = 'Anelastic approximation requires momentum_advec = "ws-scheme"'
606       CALL message( 'check_parameters', 'PA0447', 1, 2, 0, 6, 0 )
607    ENDIF
608    IF ( TRIM( approximation ) == 'anelastic'  .AND.  TRIM( psolver ) == 'multigrid' )  THEN
609       message_string = 'Anelastic approximation currently only supports psolver = "poisfft", ' // &
610                        'psolver = "sor" and psolver = "multigrid_noopt"'
611       CALL message( 'check_parameters', 'PA0448', 1, 2, 0, 6, 0 )
612    ENDIF
613    IF ( TRIM( approximation ) == 'anelastic'  .AND.  conserve_volume_flow )  THEN
614       message_string = 'Anelastic approximation is not allowed with ' //                          &
615                        'conserve_volume_flow = .TRUE.'
616       CALL message( 'check_parameters', 'PA0449', 1, 2, 0, 6, 0 )
617    ENDIF
618
619!
620!-- Check flux input mode
621    IF ( TRIM( flux_input_mode ) /= 'dynamic'  .AND.  TRIM( flux_input_mode ) /= 'kinematic'       &
622         .AND.  TRIM( flux_input_mode ) /= 'approximation-specific' )  THEN
623       message_string = 'unknown flux input mode: flux_input_mode = "' //                          &
624                        TRIM( flux_input_mode ) // '"'
625       CALL message( 'check_parameters', 'PA0450', 1, 2, 0, 6, 0 )
626    ENDIF
627!
628!-- Set flux input mode according to approximation if applicable
629    IF ( TRIM( flux_input_mode ) == 'approximation-specific' )  THEN
630       IF ( TRIM( approximation ) == 'anelastic' )  THEN
631          flux_input_mode = 'dynamic'
632       ELSEIF ( TRIM( approximation ) == 'boussinesq' )  THEN
633          flux_input_mode = 'kinematic'
634       ENDIF
635    ENDIF
636
637!
638!-- Check flux output mode
639    IF ( TRIM( flux_output_mode ) /= 'dynamic'  .AND.  TRIM( flux_output_mode ) /= 'kinematic'     &
640         .AND.  TRIM( flux_output_mode ) /= 'approximation-specific' )  THEN
641       message_string = 'unknown flux output mode: flux_output_mode = "' //                        &
642                        TRIM( flux_output_mode ) // '"'
643       CALL message( 'check_parameters', 'PA0451', 1, 2, 0, 6, 0 )
644    ENDIF
645!
646!-- Set flux output mode according to approximation if applicable
647    IF ( TRIM( flux_output_mode ) == 'approximation-specific' )  THEN
648       IF ( TRIM( approximation ) == 'anelastic' )  THEN
649          flux_output_mode = 'dynamic'
650       ELSEIF ( TRIM( approximation ) == 'boussinesq' )  THEN
651          flux_output_mode = 'kinematic'
652       ENDIF
653    ENDIF
654
655
656!
657!-- When the land- or urban-surface model is used, the flux output must be dynamic.
658    IF ( land_surface  .OR.  urban_surface )  THEN
659       flux_output_mode = 'dynamic'
660    ENDIF
661
662!
663!-- Set the flux output units according to flux_output_mode
664    IF ( TRIM( flux_output_mode ) == 'kinematic' )  THEN
665        heatflux_output_unit              = 'K m/s'
666        waterflux_output_unit             = 'kg/kg m/s'
667        momentumflux_output_unit          = 'm2/s2'
668    ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' )  THEN
669        heatflux_output_unit              = 'W/m2'
670        waterflux_output_unit             = 'W/m2'
671        momentumflux_output_unit          = 'N/m2'
672    ENDIF
673
674!
675!-- set time series output units for fluxes
676    dots_unit(14:16) = TRIM( heatflux_output_unit )
677    dots_unit(21)    = TRIM( waterflux_output_unit )
678    dots_unit(19:20) = TRIM( momentumflux_output_unit )
679
680!
681!-- Add other module specific timeseries
682    CALL module_interface_check_data_output_ts( dots_max, dots_num, dots_label, dots_unit )
683
684!
685!-- Check if maximum number of allowed timeseries is exceeded
686    IF ( dots_num > dots_max )  THEN
687       WRITE( message_string, * ) 'number of time series quantities exceeds',                      &
688                                  ' its maximum of dots_max = ', dots_max,                         &
689                                  '&Please increase dots_max in netcdf_interface_mod.f90.'
690       CALL message( 'check_parameters', 'PA0194', 1, 2, 0, 6, 0 )
691    ENDIF
692
693!
694!-- Check whether there are any illegal values
695!-- Pressure solver:
696    IF ( psolver /= 'poisfft'  .AND.  psolver /= 'sor'  .AND.  psolver /= 'multigrid'  .AND.       &
697         psolver /= 'multigrid_noopt' )  THEN
698       message_string = 'unknown solver for perturbation pressure: psolver' //                     &
699                        ' = "' // TRIM( psolver ) // '"'
700       CALL message( 'check_parameters', 'PA0016', 1, 2, 0, 6, 0 )
701    ENDIF
702
703    IF ( psolver(1:9) == 'multigrid' )  THEN
704       IF ( cycle_mg == 'w' )  THEN
705          gamma_mg = 2
706       ELSEIF ( cycle_mg == 'v' )  THEN
707          gamma_mg = 1
708       ELSE
709          message_string = 'unknown multigrid cycle: cycle_mg = "' //  TRIM( cycle_mg ) // '"'
710          CALL message( 'check_parameters', 'PA0020', 1, 2, 0, 6, 0 )
711       ENDIF
712    ENDIF
713
714    IF ( fft_method /= 'singleton-algorithm'  .AND.  fft_method /= 'temperton-algorithm'  .AND.    &
715         fft_method /= 'fftw'                 .AND.  fft_method /= 'system-specific' )  THEN
716       message_string = 'unknown fft-algorithm: fft_method = "' // TRIM( fft_method ) // '"'
717       CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
718    ENDIF
719
720    IF( momentum_advec == 'ws-scheme' .AND.  .NOT. call_psolver_at_all_substeps  ) THEN
721        message_string = 'psolver must be called at each RK3 substep when "'//                     &
722                         TRIM(momentum_advec) // ' "is used for momentum_advec'
723        CALL message( 'check_parameters', 'PA0344', 1, 2, 0, 6, 0 )
724    END IF
725!
726!-- Advection schemes:
727    IF ( momentum_advec /= 'pw-scheme'  .AND.  momentum_advec /= 'ws-scheme'  .AND.                &
728         momentum_advec /= 'up-scheme' )  THEN
729       message_string = 'unknown advection scheme: momentum_advec = "' //                          &
730                        TRIM( momentum_advec ) // '"'
731       CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
732    ENDIF
733    IF ( ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' )                        &
734         .AND. ( timestep_scheme == 'euler' .OR.  timestep_scheme == 'runge-kutta-2' ) )  THEN
735       message_string = 'momentum_advec or scalar_advec = "' // TRIM( momentum_advec ) //          &
736                        '" is not allowed with timestep_scheme = "' //                             &
737                        TRIM( timestep_scheme ) // '"'
738       CALL message( 'check_parameters', 'PA0023', 1, 2, 0, 6, 0 )
739    ENDIF
740    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme'  .AND.                    &
741         scalar_advec /= 'bc-scheme'  .AND.  scalar_advec /= 'up-scheme' )  THEN
742       message_string = 'unknown advection scheme: scalar_advec = "' // TRIM( scalar_advec ) // '"'
743       CALL message( 'check_parameters', 'PA0024', 1, 2, 0, 6, 0 )
744    ENDIF
745    IF ( scalar_advec == 'bc-scheme'  .AND.  loop_optimization == 'cache' )  THEN
746       message_string = 'advection_scheme scalar_advec = "' // TRIM( scalar_advec ) //             &
747                        '" not implemented for loop_optimization = "' //                           &
748                        TRIM( loop_optimization ) // '"'
749       CALL message( 'check_parameters', 'PA0026', 1, 2, 0, 6, 0 )
750    ENDIF
751
752    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets  .AND.  .NOT. use_upstream_for_tke     &
753         .AND.  scalar_advec /= 'ws-scheme' )  THEN
754       use_upstream_for_tke = .TRUE.
755       message_string = 'use_upstream_for_tke is set to .TRUE. because ' //                        &
756                        'use_sgs_for_particles = .TRUE. and scalar_advec /= ws-scheme'
757       CALL message( 'check_parameters', 'PA0025', 0, 1, 0, 6, 0 )
758    ENDIF
759
760!
761!-- Set LOGICAL switches to enhance performance
762    IF ( momentum_advec == 'ws-scheme' )  ws_scheme_mom = .TRUE.
763    IF ( scalar_advec   == 'ws-scheme' )  ws_scheme_sca = .TRUE.
764
765
766!
767!-- Timestep schemes:
768    SELECT CASE ( TRIM( timestep_scheme ) )
769
770       CASE ( 'euler' )
771          intermediate_timestep_count_max = 1
772
773       CASE ( 'runge-kutta-2' )
774          intermediate_timestep_count_max = 2
775
776       CASE ( 'runge-kutta-3' )
777          intermediate_timestep_count_max = 3
778
779       CASE DEFAULT
780          message_string = 'unknown timestep scheme: timestep_scheme = "' //                       &
781                           TRIM( timestep_scheme ) // '"'
782          CALL message( 'check_parameters', 'PA0027', 1, 2, 0, 6, 0 )
783
784    END SELECT
785
786    IF ( ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' )                     &
787         .AND.  timestep_scheme(1:5) == 'runge' ) THEN
788       message_string = 'momentum advection scheme "' // TRIM( momentum_advec ) //                 &
789                        '" & does not work with timestep_scheme "' // TRIM( timestep_scheme )      &
790                        // '"'
791       CALL message( 'check_parameters', 'PA0029', 1, 2, 0, 6, 0 )
792    ENDIF
793!
794!-- Check for proper settings for microphysics
795    IF ( bulk_cloud_model  .AND.  cloud_droplets )  THEN
796       message_string = 'bulk_cloud_model = .TRUE. is not allowed with cloud_droplets = .TRUE.'
797       CALL message( 'check_parameters', 'PA0442', 1, 2, 0, 6, 0 )
798    ENDIF
799
800!
801!-- Initializing actions must have been set by the user
802    IF ( TRIM( initializing_actions ) == '' )  THEN
803       message_string = 'no value specified for initializing_actions'
804       CALL message( 'check_parameters', 'PA0149', 1, 2, 0, 6, 0 )
805    ENDIF
806
807    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.                                &
808         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
809!
810!--    No restart run: several initialising actions are possible
811       action = initializing_actions
812       DO  WHILE ( TRIM( action ) /= '' )
813          position = INDEX( action, ' ' )
814          SELECT CASE ( action(1:position-1) )
815
816             CASE ( 'set_constant_profiles', 'set_1d-model_profiles', 'by_user',                   &
817                    'initialize_vortex', 'initialize_ptanom', 'initialize_bubble', 'inifor' )
818                action = action(position+1:)
819
820             CASE DEFAULT
821                message_string = 'initializing_action = "' //                                      &
822                                 TRIM( action ) // '" unknown or not allowed'
823                CALL message( 'check_parameters', 'PA0030', 1, 2, 0, 6, 0 )
824
825          END SELECT
826       ENDDO
827    ENDIF
828
829    IF ( TRIM( initializing_actions ) == 'initialize_vortex'  .AND.  conserve_volume_flow ) THEN
830         message_string = 'initializing_actions = "initialize_vortex"' //                          &
831                          ' is not allowed with conserve_volume_flow = .T.'
832       CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
833    ENDIF
834
835
836    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.                        &
837         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
838       message_string = 'initializing_actions = "set_constant_profiles"' //                        &
839                        ' and "set_1d-model_profiles" are not allowed simultaneously'
840       CALL message( 'check_parameters', 'PA0031', 1, 2, 0, 6, 0 )
841    ENDIF
842
843    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.                        &
844         INDEX( initializing_actions, 'by_user' ) /= 0 )  THEN
845       message_string = 'initializing_actions = "set_constant_profiles"' //                        &
846                        ' and "by_user" are not allowed simultaneously'
847       CALL message( 'check_parameters', 'PA0032', 1, 2, 0, 6, 0 )
848    ENDIF
849
850    IF ( INDEX( initializing_actions, 'by_user' ) /= 0  .AND.                                      &
851         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
852       message_string = 'initializing_actions = "by_user" and ' //                                 &
853                        '"set_1d-model_profiles" are not allowed simultaneously'
854       CALL message( 'check_parameters', 'PA0033', 1, 2, 0, 6, 0 )
855    ENDIF
856!
857!-- In case of spinup and nested run, spinup end time must be identical in order to have
858!-- synchronously running simulations.
859    IF ( nested_run )  THEN
860#if defined( __parallel )
861       CALL MPI_ALLREDUCE( spinup_time, spinup_time_max, 1, MPI_REAL, MPI_MAX, MPI_COMM_WORLD,     &
862                           ierr )
863       CALL MPI_ALLREDUCE( dt_spinup,   dt_spinup_max,   1, MPI_REAL, MPI_MAX, MPI_COMM_WORLD,     &
864                           ierr )
865
866       IF ( spinup_time /= spinup_time_max  .OR.  dt_spinup /= dt_spinup_max )  THEN
867          message_string = 'In case of nesting, spinup_time and ' //                               &
868                           'dt_spinup must be identical in all parent and child domains.'
869          CALL message( 'check_parameters', 'PA0489', 3, 2, 0, 6, 0 )
870       ENDIF
871#endif
872    ENDIF
873
874    IF ( bulk_cloud_model  .AND.  .NOT. humidity )  THEN
875       WRITE( message_string, * ) 'bulk_cloud_model = ', bulk_cloud_model,                         &
876              ' is not allowed with humidity = ', humidity
877       CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 )
878    ENDIF
879
880    IF ( humidity  .AND.  sloping_surface )  THEN
881       message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' //                       &
882                        'are not allowed simultaneously'
883       CALL message( 'check_parameters', 'PA0036', 1, 2, 0, 6, 0 )
884    ENDIF
885
886!-- Check the module settings
887    CALL module_interface_check_parameters
888
889!
890!-- In case of no restart run, check initialising parameters and calculate further quantities
891    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
892
893!
894!--    Initial profiles for 1D and 3D model, respectively (u,v further below)
895       pt_init = pt_surface
896       IF ( humidity       )  q_init  = q_surface
897       IF ( passive_scalar )  s_init  = s_surface
898
899!--
900!--    If required, compute initial profile of the geostrophic wind (component ug)
901       i = 1
902       gradient = 0.0_wp
903
904       IF ( .NOT. ocean_mode )  THEN
905
906          ug_vertical_gradient_level_ind(1) = 0
907          ug(0) = ug_surface
908          DO  k = 1, nzt+1
909             IF ( i < 11 )  THEN
910                IF ( ug_vertical_gradient_level(i) < zu(k)  .AND.                                  &
911                     ug_vertical_gradient_level(i) >= 0.0_wp )  THEN
912                   gradient = ug_vertical_gradient(i) / 100.0_wp
913                   ug_vertical_gradient_level_ind(i) = k - 1
914                   i = i + 1
915                ENDIF
916             ENDIF
917             IF ( gradient /= 0.0_wp )  THEN
918                IF ( k /= 1 )  THEN
919                   ug(k) = ug(k-1) + dzu(k) * gradient
920                ELSE
921                   ug(k) = ug_surface + dzu(k) * gradient
922                ENDIF
923             ELSE
924                ug(k) = ug(k-1)
925             ENDIF
926          ENDDO
927
928       ELSE
929
930          ug_vertical_gradient_level_ind(1) = nzt+1
931          ug(nzt+1) = ug_surface
932          DO  k = nzt, nzb, -1
933             IF ( i < 11 )  THEN
934                IF ( ug_vertical_gradient_level(i) > zu(k)  .AND.                                  &
935                     ug_vertical_gradient_level(i) <= 0.0_wp )  THEN
936                   gradient = ug_vertical_gradient(i) / 100.0_wp
937                   ug_vertical_gradient_level_ind(i) = k + 1
938                   i = i + 1
939                ENDIF
940             ENDIF
941             IF ( gradient /= 0.0_wp )  THEN
942                IF ( k /= nzt )  THEN
943                   ug(k) = ug(k+1) - dzu(k+1) * gradient
944                ELSE
945                   ug(k)   = ug_surface - 0.5_wp * dzu(k+1) * gradient
946                   ug(k+1) = ug_surface + 0.5_wp * dzu(k+1) * gradient
947                ENDIF
948             ELSE
949                ug(k) = ug(k+1)
950             ENDIF
951          ENDDO
952
953       ENDIF
954
955!
956!--    In case of no given gradients for ug, choose a zero gradient
957       IF ( ug_vertical_gradient_level(1) == -9999999.9_wp )  THEN
958          ug_vertical_gradient_level(1) = 0.0_wp
959       ENDIF
960
961!
962!--
963!--    If required, compute initial profile of the geostrophic wind (component vg)
964       i = 1
965       gradient = 0.0_wp
966
967       IF ( .NOT. ocean_mode )  THEN
968
969          vg_vertical_gradient_level_ind(1) = 0
970          vg(0) = vg_surface
971          DO  k = 1, nzt+1
972             IF ( i < 11 )  THEN
973                IF ( vg_vertical_gradient_level(i) < zu(k)  .AND.                                  &
974                     vg_vertical_gradient_level(i) >= 0.0_wp )  THEN
975                   gradient = vg_vertical_gradient(i) / 100.0_wp
976                   vg_vertical_gradient_level_ind(i) = k - 1
977                   i = i + 1
978                ENDIF
979             ENDIF
980             IF ( gradient /= 0.0_wp )  THEN
981                IF ( k /= 1 )  THEN
982                   vg(k) = vg(k-1) + dzu(k) * gradient
983                ELSE
984                   vg(k) = vg_surface + dzu(k) * gradient
985                ENDIF
986             ELSE
987                vg(k) = vg(k-1)
988             ENDIF
989          ENDDO
990
991       ELSE
992
993          vg_vertical_gradient_level_ind(1) = nzt+1
994          vg(nzt+1) = vg_surface
995          DO  k = nzt, nzb, -1
996             IF ( i < 11 )  THEN
997                IF ( vg_vertical_gradient_level(i) > zu(k)  .AND.                                  &
998                     vg_vertical_gradient_level(i) <= 0.0_wp )  THEN
999                   gradient = vg_vertical_gradient(i) / 100.0_wp
1000                   vg_vertical_gradient_level_ind(i) = k + 1
1001                   i = i + 1
1002                ENDIF
1003             ENDIF
1004             IF ( gradient /= 0.0_wp )  THEN
1005                IF ( k /= nzt )  THEN
1006                   vg(k) = vg(k+1) - dzu(k+1) * gradient
1007                ELSE
1008                   vg(k)   = vg_surface - 0.5_wp * dzu(k+1) * gradient
1009                   vg(k+1) = vg_surface + 0.5_wp * dzu(k+1) * gradient
1010                ENDIF
1011             ELSE
1012                vg(k) = vg(k+1)
1013             ENDIF
1014          ENDDO
1015
1016       ENDIF
1017
1018!
1019!--    In case of no given gradients for vg, choose a zero gradient
1020       IF ( vg_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1021          vg_vertical_gradient_level(1) = 0.0_wp
1022       ENDIF
1023
1024!
1025!--    Let the initial wind profiles be the calculated ug/vg profiles or interpolate them from wind
1026!--    profile data (if given)
1027       IF ( u_profile(1) == 9999999.9_wp  .AND.  v_profile(1) == 9999999.9_wp )  THEN
1028
1029          u_init = ug
1030          v_init = vg
1031
1032       ELSEIF ( u_profile(1) == 0.0_wp  .AND.  v_profile(1) == 0.0_wp )  THEN
1033
1034          IF ( uv_heights(1) /= 0.0_wp )  THEN
1035             message_string = 'uv_heights(1) must be 0.0'
1036             CALL message( 'check_parameters', 'PA0345', 1, 2, 0, 6, 0 )
1037          ENDIF
1038
1039          IF ( omega /= 0.0_wp )  THEN
1040             message_string = 'Coriolis force is switched on (omega /= 0.0) and' //                &
1041                              ' initial profiles u_profile and v_profile are prescribed'
1042             CALL message( 'check_parameters', 'PA0347', 0, 0, 0, 6, 0 )
1043          ENDIF
1044
1045          use_prescribed_profile_data = .TRUE.
1046
1047          kk = 1
1048          u_init(0) = 0.0_wp
1049          v_init(0) = 0.0_wp
1050
1051          DO  k = 1, nz+1
1052
1053             IF ( kk < 200 )  THEN
1054                DO  WHILE ( uv_heights(kk+1) <= zu(k) )
1055                   kk = kk + 1
1056                   IF ( kk == 200 )  EXIT
1057                ENDDO
1058             ENDIF
1059
1060             IF ( kk < 200  .AND.  uv_heights(kk+1) /= 9999999.9_wp )  THEN
1061                u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /                           &
1062                                            ( uv_heights(kk+1) - uv_heights(kk) ) *                &
1063                                            ( u_profile(kk+1) - u_profile(kk) )
1064                v_init(k) = v_profile(kk) + ( zu(k) - uv_heights(kk) ) /                           &
1065                                            ( uv_heights(kk+1) - uv_heights(kk) ) *                &
1066                                            ( v_profile(kk+1) - v_profile(kk) )
1067             ELSE
1068                u_init(k) = u_profile(kk)
1069                v_init(k) = v_profile(kk)
1070             ENDIF
1071
1072          ENDDO
1073
1074       ELSE
1075
1076          message_string = 'u_profile(1) and v_profile(1) must be 0.0'
1077          CALL message( 'check_parameters', 'PA0346', 1, 2, 0, 6, 0 )
1078
1079       ENDIF
1080
1081!
1082!--    Compute initial temperature profile using the given temperature gradients
1083       IF (  .NOT.  neutral )  THEN
1084          CALL init_vertical_profiles( pt_vertical_gradient_level_ind, pt_vertical_gradient_level, &
1085                                       pt_vertical_gradient, pt_init, pt_surface, bc_pt_t_val )
1086       ENDIF
1087!
1088!--    Compute initial humidity profile using the given humidity gradients
1089       IF ( humidity )  THEN
1090          CALL init_vertical_profiles( q_vertical_gradient_level_ind, q_vertical_gradient_level,   &
1091                                       q_vertical_gradient, q_init, q_surface, bc_q_t_val )
1092       ENDIF
1093!
1094!--    Compute initial scalar profile using the given scalar gradients
1095       IF ( passive_scalar )  THEN
1096          CALL init_vertical_profiles( s_vertical_gradient_level_ind, s_vertical_gradient_level,   &
1097                                       s_vertical_gradient, s_init, s_surface, bc_s_t_val )
1098       ENDIF
1099!
1100!--    TODO
1101!--    Compute initial chemistry profile using the given chemical species gradients
1102!--    Russo: Is done in chem_init --> kanani: Revise
1103
1104    ENDIF
1105
1106!
1107!-- Check if the control parameter use_subsidence_tendencies is used correctly
1108    IF ( use_subsidence_tendencies  .AND.  .NOT.  large_scale_subsidence )  THEN
1109       message_string = 'The usage of use_subsidence_tendencies ' //                               &
1110                        'requires large_scale_subsidence = .T..'
1111       CALL message( 'check_parameters', 'PA0396', 1, 2, 0, 6, 0 )
1112    ELSEIF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_forcing )  THEN
1113       message_string = 'The usage of use_subsidence_tendencies ' //           &
1114                        'requires large_scale_forcing = .T..'
1115       CALL message( 'check_parameters', 'PA0397', 1, 2, 0, 6, 0 )
1116    ENDIF
1117
1118!
1119!-- Initialize large scale subsidence if required
1120    If ( large_scale_subsidence )  THEN
1121       IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp  .AND. .NOT. large_scale_forcing )    &
1122       THEN
1123          CALL init_w_subsidence
1124       ENDIF
1125!
1126!--    In case large_scale_forcing is used, profiles for subsidence velocity are read in from file
1127!--    LSF_DATA
1128
1129       IF ( subs_vertical_gradient_level(1) == -9999999.9_wp  .AND. .NOT. large_scale_forcing )    &
1130       THEN
1131          message_string = 'There is no default large scale vertical velocity profile set. ' //    &
1132                           'Specify the subsidence velocity profile via subs_vertical_gradient' // &
1133                           ' and subs_vertical_gradient_level.'
1134          CALL message( 'check_parameters', 'PA0380', 1, 2, 0, 6, 0 )
1135       ENDIF
1136    ELSE
1137        IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp )  THEN
1138           message_string = 'Enable usage of large scale subsidence by ' //                        &
1139                            'setting large_scale_subsidence = .T..'
1140          CALL message( 'check_parameters', 'PA0381', 1, 2, 0, 6, 0 )
1141        ENDIF
1142    ENDIF
1143
1144!
1145!-- Overwrite parameters from namelist if necessary and compute Coriolis parameter.
1146!-- @todo - move initialization of f and fs to coriolis_mod.
1147    IF ( input_pids_static )  THEN
1148       latitude       = init_model%latitude
1149       longitude      = init_model%longitude
1150       rotation_angle = init_model%rotation_angle
1151    ENDIF
1152
1153    f  = 2.0_wp * omega * SIN( latitude / 180.0_wp * pi )
1154    fs = 2.0_wp * omega * COS( latitude / 180.0_wp * pi )
1155
1156!
1157!-- Check and set buoyancy related parameters and switches
1158    IF ( reference_state == 'horizontal_average' )  THEN
1159       CONTINUE
1160    ELSEIF ( reference_state == 'initial_profile' )  THEN
1161       use_initial_profile_as_reference = .TRUE.
1162    ELSEIF ( reference_state == 'single_value' )  THEN
1163       use_single_reference_value = .TRUE.
1164       IF ( pt_reference == 9999999.9_wp )  pt_reference = pt_surface
1165       vpt_reference = pt_reference * ( 1.0_wp + 0.61_wp * q_surface )
1166    ELSE
1167       message_string = 'illegal value for reference_state: "' // TRIM( reference_state ) // '"'
1168       CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 )
1169    ENDIF
1170
1171!
1172!-- In case of a given slope, compute the relevant quantities
1173    IF ( alpha_surface /= 0.0_wp )  THEN
1174       IF ( ABS( alpha_surface ) > 90.0_wp )  THEN
1175          WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface, ' ) must be < 90.0'
1176          CALL message( 'check_parameters', 'PA0043', 1, 2, 0, 6, 0 )
1177       ENDIF
1178       sloping_surface = .TRUE.
1179       cos_alpha_surface = COS( alpha_surface / 180.0_wp * pi )
1180       sin_alpha_surface = SIN( alpha_surface / 180.0_wp * pi )
1181    ENDIF
1182
1183!
1184!-- Check time step and cfl_factor
1185    IF ( dt /= -1.0_wp )  THEN
1186       IF ( dt <= 0.0_wp )  THEN
1187          WRITE( message_string, * ) 'dt = ', dt , ' <= 0.0'
1188          CALL message( 'check_parameters', 'PA0044', 1, 2, 0, 6, 0 )
1189       ENDIF
1190       dt_3d = dt
1191       dt_fixed = .TRUE.
1192    ENDIF
1193
1194    IF ( cfl_factor <= 0.0_wp  .OR.  cfl_factor > 1.0_wp )  THEN
1195       IF ( cfl_factor == -1.0_wp )  THEN
1196          IF ( timestep_scheme == 'runge-kutta-2' )  THEN
1197             cfl_factor = 0.8_wp
1198          ELSEIF ( timestep_scheme == 'runge-kutta-3' )  THEN
1199             cfl_factor = 0.9_wp
1200          ELSE
1201             cfl_factor = 0.9_wp
1202          ENDIF
1203       ELSE
1204          WRITE( message_string, * ) 'cfl_factor = ', cfl_factor,                                  &
1205                 ' out of range &0.0 < cfl_factor <= 1.0 is required'
1206          CALL message( 'check_parameters', 'PA0045', 1, 2, 0, 6, 0 )
1207       ENDIF
1208    ENDIF
1209
1210!
1211!-- Store simulated time at begin
1212    simulated_time_at_begin = simulated_time
1213
1214!
1215!-- Store reference time for coupled runs and change the coupling flag, if ...
1216    IF ( simulated_time == 0.0_wp )  THEN
1217       IF ( coupling_start_time == 0.0_wp )  THEN
1218          time_since_reference_point = 0.0_wp
1219       ELSEIF ( time_since_reference_point < 0.0_wp )  THEN
1220          run_coupled = .FALSE.
1221       ENDIF
1222    ENDIF
1223
1224!
1225!-- Set wind speed in the Galilei-transformed system
1226    IF ( galilei_transformation )  THEN
1227       IF ( use_ug_for_galilei_tr                    .AND.                                         &
1228            ug_vertical_gradient_level(1) == 0.0_wp  .AND.                                         &
1229            ug_vertical_gradient(1)       == 0.0_wp  .AND.                                         &
1230            vg_vertical_gradient_level(1) == 0.0_wp  .AND.                                         &
1231            vg_vertical_gradient(1)       == 0.0_wp )  THEN
1232          u_gtrans = ug_surface * 0.6_wp
1233          v_gtrans = vg_surface * 0.6_wp
1234       ELSEIF ( use_ug_for_galilei_tr  .AND.  ( ug_vertical_gradient_level(1) /= 0.0_wp .OR.       &
1235                                                ug_vertical_gradient(1) /= 0.0_wp ) )  THEN
1236          message_string = 'baroclinity (ug) not allowed simultaneously' //                        &
1237                           ' with galilei transformation'
1238          CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 )
1239       ELSEIF ( use_ug_for_galilei_tr  .AND.  ( vg_vertical_gradient_level(1) /= 0.0_wp  .OR.      &
1240                                                vg_vertical_gradient(1) /= 0.0_wp ) )  THEN
1241          message_string = 'baroclinity (vg) not allowed simultaneously' //                        &
1242                           ' with galilei transformation'
1243          CALL message( 'check_parameters', 'PA0047', 1, 2, 0, 6, 0 )
1244       ELSE
1245          message_string = 'variable translation speed used for Galilei-transformation, which ' // &
1246                           'may cause & instabilities in stably stratified regions'
1247          CALL message( 'check_parameters', 'PA0048', 0, 1, 0, 6, 0 )
1248       ENDIF
1249    ENDIF
1250
1251!
1252!-- In case of using a prandtl-layer, calculated (or prescribed) surface fluxes have to be used in
1253!-- the diffusion-terms
1254    IF ( constant_flux_layer )  use_surface_fluxes = .TRUE.
1255
1256!
1257!-- Check boundary conditions and set internal variables:
1258!-- Attention: the lateral boundary conditions have been already checked in parin
1259!
1260!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-Willimas or
1261!-- Wicker - Skamarock advection scheme. Several schemes and tools do not work with non-cyclic
1262!-- boundary conditions.
1263    IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
1264       IF ( psolver(1:9) /= 'multigrid' )  THEN
1265          message_string = 'non-cyclic lateral boundaries do not allow ' // 'psolver = "' //       &
1266                           TRIM( psolver ) // '"'
1267          CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
1268       ENDIF
1269       IF ( momentum_advec /= 'pw-scheme'  .AND.  momentum_advec /= 'ws-scheme' )  THEN
1270
1271          message_string = 'non-cyclic lateral boundaries do not allow momentum_advec = "' //      &
1272                           TRIM( momentum_advec ) // '"'
1273          CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
1274       ENDIF
1275       IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' )  THEN
1276          message_string = 'non-cyclic lateral boundaries do not allow scalar_advec = "' //        &
1277                           TRIM( scalar_advec ) // '"'
1278          CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 )
1279       ENDIF
1280       IF ( galilei_transformation )  THEN
1281          message_string = 'non-cyclic lateral boundaries do not allow galilei_transformation = .T.'
1282          CALL message( 'check_parameters', 'PA0054', 1, 2, 0, 6, 0 )
1283       ENDIF
1284    ENDIF
1285
1286!
1287!-- Bottom boundary condition for the turbulent Kinetic energy
1288    IF ( bc_e_b == 'neumann' )  THEN
1289       ibc_e_b = 1
1290    ELSEIF ( bc_e_b == '(u*)**2+neumann' )  THEN
1291       ibc_e_b = 2
1292       IF ( .NOT. constant_flux_layer )  THEN
1293          bc_e_b = 'neumann'
1294          ibc_e_b = 1
1295          message_string = 'boundary condition bc_e_b changed to "' // TRIM( bc_e_b ) // '"'
1296          CALL message( 'check_parameters', 'PA0057', 0, 1, 0, 6, 0 )
1297       ENDIF
1298    ELSE
1299       message_string = 'unknown boundary condition: bc_e_b = "' // TRIM( bc_e_b ) // '"'
1300       CALL message( 'check_parameters', 'PA0058', 1, 2, 0, 6, 0 )
1301    ENDIF
1302
1303!
1304!-- Boundary conditions for perturbation pressure
1305    IF ( bc_p_b == 'dirichlet' )  THEN
1306       ibc_p_b = 0
1307    ELSEIF ( bc_p_b == 'neumann' )  THEN
1308       ibc_p_b = 1
1309    ELSE
1310       message_string = 'unknown boundary condition: bc_p_b = "' // TRIM( bc_p_b ) // '"'
1311       CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 )
1312    ENDIF
1313
1314    IF ( bc_p_t == 'dirichlet' )  THEN
1315       ibc_p_t = 0
1316!-- TO_DO: later set bc_p_t to neumann before, in case of nested domain
1317    ELSEIF ( bc_p_t == 'neumann' .OR. bc_p_t == 'nested' )  THEN
1318       ibc_p_t = 1
1319    ELSE
1320       message_string = 'unknown boundary condition: bc_p_t = "' // TRIM( bc_p_t ) // '"'
1321       CALL message( 'check_parameters', 'PA0061', 1, 2, 0, 6, 0 )
1322    ENDIF
1323
1324!
1325!-- Boundary conditions for potential temperature
1326    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
1327       ibc_pt_b = 2
1328    ELSE
1329       IF ( bc_pt_b == 'dirichlet' )  THEN
1330          ibc_pt_b = 0
1331       ELSEIF ( bc_pt_b == 'neumann' )  THEN
1332          ibc_pt_b = 1
1333       ELSE
1334          message_string = 'unknown boundary condition: bc_pt_b = "' // TRIM( bc_pt_b ) // '"'
1335          CALL message( 'check_parameters', 'PA0062', 1, 2, 0, 6, 0 )
1336       ENDIF
1337    ENDIF
1338
1339    IF ( bc_pt_t == 'dirichlet' )  THEN
1340       ibc_pt_t = 0
1341    ELSEIF ( bc_pt_t == 'neumann' )  THEN
1342       ibc_pt_t = 1
1343    ELSEIF ( bc_pt_t == 'initial_gradient' )  THEN
1344       ibc_pt_t = 2
1345    ELSEIF ( bc_pt_t == 'nested'  .OR.  bc_pt_t == 'nesting_offline' )  THEN
1346       ibc_pt_t = 3
1347    ELSE
1348       message_string = 'unknown boundary condition: bc_pt_t = "' // TRIM( bc_pt_t ) // '"'
1349       CALL message( 'check_parameters', 'PA0063', 1, 2, 0, 6, 0 )
1350    ENDIF
1351
1352    IF ( ANY( wall_heatflux /= 0.0_wp )  .AND.  surface_heatflux == 9999999.9_wp )  THEN
1353       message_string = 'wall_heatflux additionally requires setting of surface_heatflux'
1354       CALL message( 'check_parameters', 'PA0443', 1, 2, 0, 6, 0 )
1355    ENDIF
1356
1357!
1358!   This IF clause needs revision, got too complex!!
1359    IF ( surface_heatflux == 9999999.9_wp  )  THEN
1360       constant_heatflux = .FALSE.
1361       IF ( large_scale_forcing  .OR.  land_surface  .OR.  urban_surface )  THEN
1362          IF ( ibc_pt_b == 0 )  THEN
1363             constant_heatflux = .FALSE.
1364          ELSEIF ( ibc_pt_b == 1 )  THEN
1365             constant_heatflux = .TRUE.
1366             surface_heatflux = 0.0_wp
1367          ENDIF
1368       ENDIF
1369    ELSE
1370       constant_heatflux = .TRUE.
1371    ENDIF
1372
1373    IF ( top_heatflux     == 9999999.9_wp )  constant_top_heatflux = .FALSE.
1374
1375    IF ( neutral )  THEN
1376
1377       IF ( surface_heatflux /= 0.0_wp  .AND.  surface_heatflux /= 9999999.9_wp )  THEN
1378          message_string = 'heatflux must not be set for pure neutral flow'
1379          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
1380       ENDIF
1381
1382       IF ( top_heatflux /= 0.0_wp  .AND.  top_heatflux /= 9999999.9_wp )  THEN
1383          message_string = 'heatflux must not be set for pure neutral flow'
1384          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
1385       ENDIF
1386
1387    ENDIF
1388
1389    IF ( top_momentumflux_u /= 9999999.9_wp  .AND.  top_momentumflux_v /= 9999999.9_wp )  THEN
1390       constant_top_momentumflux = .TRUE.
1391    ELSEIF ( .NOT. ( top_momentumflux_u == 9999999.9_wp  .AND.                                     &
1392           top_momentumflux_v == 9999999.9_wp ) )  THEN
1393       message_string = 'both, top_momentumflux_u AND top_momentumflux_v must be set'
1394       CALL message( 'check_parameters', 'PA0064', 1, 2, 0, 6, 0 )
1395    ENDIF
1396
1397!
1398!-- A given surface temperature implies Dirichlet boundary condition for temperature. In this case
1399!-- specification of a constant heat flux is forbidden.
1400    IF ( ibc_pt_b == 0  .AND.  constant_heatflux  .AND.  surface_heatflux /= 0.0_wp )  THEN
1401       message_string = 'boundary_condition: bc_pt_b = "' // TRIM( bc_pt_b ) //                    &
1402                        '& is not allowed with constant_heatflux = .TRUE.'
1403       CALL message( 'check_parameters', 'PA0065', 1, 2, 0, 6, 0 )
1404    ENDIF
1405    IF ( constant_heatflux  .AND.  pt_surface_initial_change /= 0.0_wp )  THEN
1406       WRITE ( message_string, * )  'constant_heatflux = .TRUE. is not allo',                      &
1407               'wed with pt_surface_initial_change (/=0) = ', pt_surface_initial_change
1408       CALL message( 'check_parameters', 'PA0066', 1, 2, 0, 6, 0 )
1409    ENDIF
1410    IF ( constant_heatflux  .AND.  pt_surface_heating_rate /= 0.0_wp )  THEN
1411       WRITE ( message_string, * )  'constant_heatflux = .TRUE. is not allo',                      &
1412               'wed with pt_surface_heating_rate (/=0) = ', pt_surface_heating_rate
1413       CALL message( 'check_parameters', 'PA0497', 1, 2, 0, 6, 0 )
1414    ENDIF
1415
1416!
1417!-- A given temperature at the top implies Dirichlet boundary condition for temperature. In this
1418!-- case specification of a constant heat flux is forbidden.
1419    IF ( ibc_pt_t == 0  .AND.  constant_top_heatflux  .AND.  top_heatflux /= 0.0_wp )  THEN
1420       message_string = 'boundary_condition: bc_pt_t = "' // TRIM( bc_pt_t ) //                    &
1421                        '" is not allowed with constant_top_heatflux = .TRUE.'
1422       CALL message( 'check_parameters', 'PA0067', 1, 2, 0, 6, 0 )
1423    ENDIF
1424
1425!
1426!-- Set boundary conditions for total water content
1427    IF ( humidity )  THEN
1428
1429       IF ( ANY( wall_humidityflux /= 0.0_wp )  .AND.  surface_waterflux == 9999999.9_wp )  THEN
1430          message_string = 'wall_humidityflux additionally requires setting of surface_waterflux'
1431          CALL message( 'check_parameters', 'PA0444', 1, 2, 0, 6, 0 )
1432       ENDIF
1433
1434       CALL set_bc_scalars( 'q', bc_q_b, bc_q_t, ibc_q_b, ibc_q_t, 'PA0071', 'PA0072' )
1435
1436       IF ( surface_waterflux == 9999999.9_wp  )  THEN
1437          constant_waterflux = .FALSE.
1438          IF ( large_scale_forcing .OR. land_surface )  THEN
1439             IF ( ibc_q_b == 0 )  THEN
1440                constant_waterflux = .FALSE.
1441             ELSEIF ( ibc_q_b == 1 )  THEN
1442                constant_waterflux = .TRUE.
1443             ENDIF
1444          ENDIF
1445       ELSE
1446          constant_waterflux = .TRUE.
1447       ENDIF
1448
1449       CALL check_bc_scalars( 'q', bc_q_b, ibc_q_b, 'PA0073', 'PA0074', constant_waterflux,        &
1450                              q_surface_initial_change )
1451
1452    ENDIF
1453
1454    IF ( passive_scalar )  THEN
1455
1456       IF ( ANY( wall_scalarflux /= 0.0_wp )  .AND.  surface_scalarflux == 9999999.9_wp )  THEN
1457          message_string = 'wall_scalarflux additionally requires setting of surface_scalarflux'
1458          CALL message( 'check_parameters', 'PA0445', 1, 2, 0, 6, 0 )
1459       ENDIF
1460
1461       IF ( surface_scalarflux == 9999999.9_wp )  constant_scalarflux = .FALSE.
1462
1463       CALL set_bc_scalars( 's', bc_s_b, bc_s_t, ibc_s_b, ibc_s_t, 'PA0071', 'PA0072' )
1464
1465       CALL check_bc_scalars( 's', bc_s_b, ibc_s_b, 'PA0073', 'PA0074', constant_scalarflux,       &
1466                              s_surface_initial_change )
1467
1468       IF ( top_scalarflux == 9999999.9_wp )  constant_top_scalarflux = .FALSE.
1469!
1470!--    A fixed scalar concentration at the top implies Dirichlet boundary condition for scalar.
1471!--    Hence, in this case specification of a constant scalar flux is forbidden.
1472       IF ( ( ibc_s_t == 0 .OR. ibc_s_t == 2 )  .AND.  constant_top_scalarflux  .AND.              &
1473              top_scalarflux /= 0.0_wp )  THEN
1474          message_string = 'boundary condition: bc_s_t = "' // TRIM( bc_s_t ) //                   &
1475                           '" is not allowed with top_scalarflux /= 0.0'
1476          CALL message( 'check_parameters', 'PA0441', 1, 2, 0, 6, 0 )
1477       ENDIF
1478    ENDIF
1479
1480!
1481!-- Boundary conditions for horizontal components of wind speed
1482    IF ( bc_uv_b == 'dirichlet' )  THEN
1483       ibc_uv_b = 0
1484    ELSEIF ( bc_uv_b == 'neumann' )  THEN
1485       ibc_uv_b = 1
1486       IF ( constant_flux_layer )  THEN
1487          message_string = 'boundary condition: bc_uv_b = "' // TRIM( bc_uv_b ) //                 &
1488                           '" is not allowed with constant_flux_layer = .TRUE.'
1489          CALL message( 'check_parameters', 'PA0075', 1, 2, 0, 6, 0 )
1490       ENDIF
1491    ELSE
1492       message_string = 'unknown boundary condition: bc_uv_b = "' // TRIM( bc_uv_b ) // '"'
1493       CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )
1494    ENDIF
1495!
1496!-- In case of coupled simulations u and v at the ground in atmosphere will be assigned with the u
1497!-- and v values of the ocean surface
1498    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
1499       ibc_uv_b = 2
1500    ENDIF
1501
1502    IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
1503       bc_uv_t = 'neumann'
1504       ibc_uv_t = 1
1505    ELSE
1506       IF ( bc_uv_t == 'dirichlet' .OR. bc_uv_t == 'dirichlet_0' )  THEN
1507          ibc_uv_t = 0
1508          IF ( bc_uv_t == 'dirichlet_0' )  THEN
1509!
1510!--          Velocities for the initial u,v-profiles are set zero at the top in case of dirichlet_0
1511!--          conditions
1512             u_init(nzt+1)    = 0.0_wp
1513             v_init(nzt+1)    = 0.0_wp
1514          ENDIF
1515       ELSEIF ( bc_uv_t == 'neumann' )  THEN
1516          ibc_uv_t = 1
1517       ELSEIF ( bc_uv_t == 'nested'  .OR.  bc_uv_t == 'nesting_offline' )  THEN
1518          ibc_uv_t = 3
1519       ELSE
1520          message_string = 'unknown boundary condition: bc_uv_t = "' // TRIM( bc_uv_t ) // '"'
1521          CALL message( 'check_parameters', 'PA0077', 1, 2, 0, 6, 0 )
1522       ENDIF
1523    ENDIF
1524
1525!
1526!-- Compute and check, respectively, the Rayleigh Damping parameter
1527    IF ( rayleigh_damping_factor == -1.0_wp )  THEN
1528       rayleigh_damping_factor = 0.0_wp
1529    ELSE
1530       IF ( rayleigh_damping_factor < 0.0_wp  .OR.  rayleigh_damping_factor > 1.0_wp )  THEN
1531          WRITE( message_string, * )  'rayleigh_damping_factor = ', rayleigh_damping_factor,       &
1532                 ' out of range [0.0,1.0]'
1533          CALL message( 'check_parameters', 'PA0078', 1, 2, 0, 6, 0 )
1534       ENDIF
1535    ENDIF
1536
1537    IF ( rayleigh_damping_height == -1.0_wp )  THEN
1538       IF ( .NOT. ocean_mode )  THEN
1539          rayleigh_damping_height = 0.66666666666_wp * zu(nzt)
1540       ELSE
1541          rayleigh_damping_height = 0.66666666666_wp * zu(nzb)
1542       ENDIF
1543    ELSE
1544       IF ( .NOT. ocean_mode )  THEN
1545          IF ( rayleigh_damping_height < 0.0_wp  .OR.  rayleigh_damping_height > zu(nzt) )  THEN
1546             WRITE( message_string, * )  'rayleigh_damping_height = ',  rayleigh_damping_height,   &
1547                    ' out of range [0.0,', zu(nzt), ']'
1548             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
1549          ENDIF
1550       ELSE
1551          IF ( rayleigh_damping_height > 0.0_wp  .OR.  rayleigh_damping_height < zu(nzb) )  THEN
1552             WRITE( message_string, * )  'rayleigh_damping_height = ', rayleigh_damping_height,    &
1553                    ' out of range [0.0,', zu(nzb), ']'
1554             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
1555          ENDIF
1556       ENDIF
1557    ENDIF
1558
1559!
1560!-- Check number of chosen statistic regions
1561    IF ( statistic_regions < 0 )  THEN
1562       WRITE ( message_string, * ) 'number of statistic_regions = ', statistic_regions+1,          &
1563               ' is not allowed'
1564       CALL message( 'check_parameters', 'PA0082', 1, 2, 0, 6, 0 )
1565    ENDIF
1566    IF ( normalizing_region > statistic_regions  .OR.  normalizing_region < 0)  THEN
1567       WRITE ( message_string, * ) 'normalizing_region = ', normalizing_region,                    &
1568               ' must be >= 0 and <= ',statistic_regions, ' (value of statistic_regions)'
1569       CALL message( 'check_parameters', 'PA0083', 1, 2, 0, 6, 0 )
1570    ENDIF
1571
1572!
1573!-- Set the default intervals for data output, if necessary
1574!-- NOTE: dt_dosp has already been set in spectra_parin
1575    IF ( dt_data_output /= 9999999.9_wp )  THEN
1576       IF ( dt_dopr           == 9999999.9_wp )  dt_dopr           = dt_data_output
1577       IF ( dt_dopts          == 9999999.9_wp )  dt_dopts          = dt_data_output
1578       IF ( dt_do2d_xy        == 9999999.9_wp )  dt_do2d_xy        = dt_data_output
1579       IF ( dt_do2d_xz        == 9999999.9_wp )  dt_do2d_xz        = dt_data_output
1580       IF ( dt_do2d_yz        == 9999999.9_wp )  dt_do2d_yz        = dt_data_output
1581       IF ( dt_do3d           == 9999999.9_wp )  dt_do3d           = dt_data_output
1582       IF ( dt_data_output_av == 9999999.9_wp )  dt_data_output_av = dt_data_output
1583       DO  mid = 1, max_masks
1584          IF ( dt_domask(mid) == 9999999.9_wp )  dt_domask(mid)    = dt_data_output
1585       ENDDO
1586    ENDIF
1587
1588!
1589!-- Set the default skip time intervals for data output, if necessary
1590    IF ( skip_time_dopr    == 9999999.9_wp )  skip_time_dopr    = skip_time_data_output
1591    IF ( skip_time_do2d_xy == 9999999.9_wp )  skip_time_do2d_xy = skip_time_data_output
1592    IF ( skip_time_do2d_xz == 9999999.9_wp )  skip_time_do2d_xz = skip_time_data_output
1593    IF ( skip_time_do2d_yz == 9999999.9_wp )  skip_time_do2d_yz = skip_time_data_output
1594    IF ( skip_time_do3d    == 9999999.9_wp )  skip_time_do3d    = skip_time_data_output
1595    IF ( skip_time_data_output_av == 9999999.9_wp )                                                &
1596                                       skip_time_data_output_av = skip_time_data_output
1597    DO  mid = 1, max_masks
1598       IF ( skip_time_domask(mid) == 9999999.9_wp )                                                &
1599                                       skip_time_domask(mid)    = skip_time_data_output
1600    ENDDO
1601
1602!
1603!-- Check the average intervals (first for 3d-data, then for profiles)
1604    IF ( averaging_interval > dt_data_output_av )  THEN
1605       WRITE( message_string, * )  'averaging_interval = ', averaging_interval,                    &
1606              ' must be <= dt_data_output_av = ', dt_data_output_av
1607       CALL message( 'check_parameters', 'PA0085', 1, 2, 0, 6, 0 )
1608    ENDIF
1609
1610    IF ( averaging_interval_pr == 9999999.9_wp )  THEN
1611       averaging_interval_pr = averaging_interval
1612    ENDIF
1613
1614    IF ( averaging_interval_pr > dt_dopr )  THEN
1615       WRITE( message_string, * )  'averaging_interval_pr = ', averaging_interval_pr,              &
1616              ' must be <= dt_dopr = ', dt_dopr
1617       CALL message( 'check_parameters', 'PA0086', 1, 2, 0, 6, 0 )
1618    ENDIF
1619
1620!
1621!-- Set the default interval for profiles entering the temporal average
1622    IF ( dt_averaging_input_pr == 9999999.9_wp )  THEN
1623       dt_averaging_input_pr = dt_averaging_input
1624    ENDIF
1625
1626!
1627!-- Set the default interval for the output of timeseries to a reasonable value (tries to minimize
1628!-- the number of calls of flow_statistics)
1629    IF ( dt_dots == 9999999.9_wp )  THEN
1630       IF ( averaging_interval_pr == 0.0_wp )  THEN
1631          dt_dots = MIN( dt_run_control, dt_dopr )
1632       ELSE
1633          dt_dots = MIN( dt_run_control, dt_averaging_input_pr )
1634       ENDIF
1635    ENDIF
1636
1637!
1638!-- Check the sample rate for averaging (first for 3d-data, then for profiles)
1639    IF ( dt_averaging_input > averaging_interval )  THEN
1640       WRITE( message_string, * )  'dt_averaging_input = ', dt_averaging_input,                    &
1641              ' must be <= averaging_interval = ', averaging_interval
1642       CALL message( 'check_parameters', 'PA0088', 1, 2, 0, 6, 0 )
1643    ENDIF
1644
1645    IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
1646       WRITE( message_string, * )  'dt_averaging_input_pr = ', dt_averaging_input_pr,              &
1647              ' must be <= averaging_interval_pr = ', averaging_interval_pr
1648       CALL message( 'check_parameters', 'PA0089', 1, 2, 0, 6, 0 )
1649    ENDIF
1650
1651!
1652!-- Check data-output variable list for any duplicates and remove them
1653    ALLOCATE( data_output_filtered(1:UBOUND( data_output_pr, DIM=1 )) )
1654    CALL filter_duplicate_strings( varnamelength, data_output_pr, data_output_filtered )
1655    data_output_pr = data_output_filtered
1656    DEALLOCATE( data_output_filtered )
1657!
1658!-- Determine the number of output profiles and check whether they are permissible
1659    DO  WHILE ( data_output_pr(dopr_n+1) /= '          ' )
1660
1661       dopr_n = dopr_n + 1
1662       i = dopr_n
1663
1664!
1665!--    Size check. Some arrays in the netcdf interface module are internally limited to 300 elements
1666!--    and would need to be increased for output of more than 300 profiles.
1667       IF ( dopr_n > 400 )  THEN
1668          message_string = 'output of more than 300 profiles requested'
1669          CALL message( 'check_parameters', 'PA0017', 1, 2, 0, 6, 0 )
1670       ENDIF
1671
1672!
1673!--    Determine internal profile number (for hom, homs) and store height levels
1674       SELECT CASE ( TRIM( data_output_pr(i) ) )
1675
1676          CASE ( 'u', '#u' )
1677             dopr_index(i) = 1
1678             dopr_unit(i)  = 'm/s'
1679             hom(:,2,1,:)  = SPREAD( zu, 2, statistic_regions+1 )
1680             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1681                dopr_initial_index(i) = 5
1682                hom(:,2,5,:)          = SPREAD( zu, 2, statistic_regions+1 )
1683                data_output_pr(i)     = data_output_pr(i)(2:)
1684             ENDIF
1685
1686          CASE ( 'v', '#v' )
1687             dopr_index(i) = 2
1688             dopr_unit(i)  = 'm/s'
1689             hom(:,2,2,:)  = SPREAD( zu, 2, statistic_regions+1 )
1690             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1691                dopr_initial_index(i) = 6
1692                hom(:,2,6,:)          = SPREAD( zu, 2, statistic_regions+1 )
1693                data_output_pr(i)     = data_output_pr(i)(2:)
1694             ENDIF
1695
1696          CASE ( 'w' )
1697             dopr_index(i) = 3
1698             dopr_unit(i)  = 'm/s'
1699             hom(:,2,3,:)  = SPREAD( zw, 2, statistic_regions+1 )
1700
1701          CASE ( 'theta', '#theta' )
1702             IF ( .NOT. bulk_cloud_model )  THEN
1703                dopr_index(i) = 4
1704                dopr_unit(i)  = 'K'
1705                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
1706                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1707                   dopr_initial_index(i) = 7
1708                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
1709                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
1710                   data_output_pr(i)     = data_output_pr(i)(2:)
1711                ENDIF
1712             ELSE
1713                dopr_index(i) = 43
1714                dopr_unit(i)  = 'K'
1715                hom(:,2,43,:)  = SPREAD( zu, 2, statistic_regions+1 )
1716                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1717                   dopr_initial_index(i) = 28
1718                   hom(:,2,28,:)         = SPREAD( zu, 2, statistic_regions+1 )
1719                   hom(nzb,2,28,:)       = 0.0_wp    ! because zu(nzb) is negative
1720                   data_output_pr(i)     = data_output_pr(i)(2:)
1721                ENDIF
1722             ENDIF
1723
1724          CASE ( 'e', '#e' )
1725             dopr_index(i)  = 8
1726             dopr_unit(i)   = 'm2/s2'
1727             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
1728             hom(nzb,2,8,:) = 0.0_wp
1729             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1730                dopr_initial_index(i) = 8
1731                hom(:,2,8,:)          = SPREAD( zu, 2, statistic_regions+1 )
1732                data_output_pr(i)     = data_output_pr(i)(2:)
1733             ENDIF
1734
1735          CASE ( 'km', '#km' )
1736             dopr_index(i)  = 9
1737             dopr_unit(i)   = 'm2/s'
1738             hom(:,2,9,:)   = SPREAD( zu, 2, statistic_regions+1 )
1739             hom(nzb,2,9,:) = 0.0_wp
1740             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1741                dopr_initial_index(i) = 23
1742                hom(:,2,23,:)         = hom(:,2,9,:)
1743                data_output_pr(i)     = data_output_pr(i)(2:)
1744             ENDIF
1745
1746          CASE ( 'kh', '#kh' )
1747             dopr_index(i)   = 10
1748             dopr_unit(i)    = 'm2/s'
1749             hom(:,2,10,:)   = SPREAD( zu, 2, statistic_regions+1 )
1750             hom(nzb,2,10,:) = 0.0_wp
1751             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1752                dopr_initial_index(i) = 24
1753                hom(:,2,24,:)         = hom(:,2,10,:)
1754                data_output_pr(i)     = data_output_pr(i)(2:)
1755             ENDIF
1756
1757          CASE ( 'l', '#l' )
1758             dopr_index(i)   = 11
1759             dopr_unit(i)    = 'm'
1760             hom(:,2,11,:)   = SPREAD( zu, 2, statistic_regions+1 )
1761             hom(nzb,2,11,:) = 0.0_wp
1762             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1763                dopr_initial_index(i) = 25
1764                hom(:,2,25,:)         = hom(:,2,11,:)
1765                data_output_pr(i)     = data_output_pr(i)(2:)
1766             ENDIF
1767
1768          CASE ( 'w"u"' )
1769             dopr_index(i) = 12
1770             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
1771             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
1772             IF ( constant_flux_layer )  hom(nzb,2,12,:) = zu(1)
1773
1774          CASE ( 'w*u*' )
1775             dopr_index(i) = 13
1776             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
1777             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
1778
1779          CASE ( 'w"v"' )
1780             dopr_index(i) = 14
1781             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
1782             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
1783             IF ( constant_flux_layer )  hom(nzb,2,14,:) = zu(1)
1784
1785          CASE ( 'w*v*' )
1786             dopr_index(i) = 15
1787             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
1788             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
1789
1790          CASE ( 'w"theta"' )
1791             dopr_index(i) = 16
1792             dopr_unit(i)  = TRIM ( heatflux_output_unit )
1793             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
1794
1795          CASE ( 'w*theta*' )
1796             dopr_index(i) = 17
1797             dopr_unit(i)  = TRIM ( heatflux_output_unit )
1798             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
1799
1800          CASE ( 'wtheta' )
1801             dopr_index(i) = 18
1802             dopr_unit(i)  = TRIM ( heatflux_output_unit )
1803             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
1804
1805          CASE ( 'wu' )
1806             dopr_index(i) = 19
1807             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
1808             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
1809             IF ( constant_flux_layer )  hom(nzb,2,19,:) = zu(1)
1810
1811          CASE ( 'wv' )
1812             dopr_index(i) = 20
1813             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
1814             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
1815             IF ( constant_flux_layer )  hom(nzb,2,20,:) = zu(1)
1816
1817          CASE ( 'w*theta*BC' )
1818             dopr_index(i) = 21
1819             dopr_unit(i)  = TRIM ( heatflux_output_unit )
1820             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
1821
1822          CASE ( 'wthetaBC' )
1823             dopr_index(i) = 22
1824             dopr_unit(i)  = TRIM ( heatflux_output_unit )
1825             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
1826
1827          CASE ( 'u*2' )
1828             dopr_index(i) = 30
1829             dopr_unit(i)  = 'm2/s2'
1830             hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
1831
1832          CASE ( 'v*2' )
1833             dopr_index(i) = 31
1834             dopr_unit(i)  = 'm2/s2'
1835             hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
1836
1837          CASE ( 'w*2' )
1838             dopr_index(i) = 32
1839             dopr_unit(i)  = 'm2/s2'
1840             hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
1841
1842          CASE ( 'theta*2' )
1843             dopr_index(i) = 33
1844             dopr_unit(i)  = 'K2'
1845             hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
1846
1847          CASE ( 'e*' )
1848             dopr_index(i) = 34
1849             dopr_unit(i)  = 'm2/s2'
1850             hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
1851
1852          CASE ( 'w*2theta*' )
1853             dopr_index(i) = 35
1854             dopr_unit(i)  = 'K m2/s2'
1855             hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
1856
1857          CASE ( 'w*theta*2' )
1858             dopr_index(i) = 36
1859             dopr_unit(i)  = 'K2 m/s'
1860             hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
1861
1862          CASE ( 'w*e*' )
1863             dopr_index(i) = 37
1864             dopr_unit(i)  = 'm3/s3'
1865             hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
1866
1867          CASE ( 'w*3' )
1868             dopr_index(i) = 38
1869             dopr_unit(i)  = 'm3/s3'
1870             hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
1871
1872          CASE ( 'Sw' )
1873             dopr_index(i) = 39
1874             dopr_unit(i)  = 'none'
1875             hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
1876
1877          CASE ( 'p' )
1878             dopr_index(i) = 40
1879             dopr_unit(i)  = 'Pa'
1880             hom(:,2,40,:) = SPREAD( zu, 2, statistic_regions+1 )
1881
1882          CASE ( 'q', '#q' )
1883             IF ( .NOT. humidity )  THEN
1884                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
1885                                 ' is not implemented for humidity = .FALSE.'
1886                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
1887             ELSE
1888                dopr_index(i) = 41
1889                dopr_unit(i)  = 'kg/kg'
1890                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
1891                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1892                   dopr_initial_index(i) = 26
1893                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
1894                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
1895                   data_output_pr(i)     = data_output_pr(i)(2:)
1896                ENDIF
1897             ENDIF
1898
1899          CASE ( 's', '#s' )
1900             IF ( .NOT. passive_scalar )  THEN
1901                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
1902                                 ' is not implemented for passive_scalar = .FALSE.'
1903                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
1904             ELSE
1905                dopr_index(i) = 115
1906                dopr_unit(i)  = 'kg/m3'
1907                hom(:,2,115,:) = SPREAD( zu, 2, statistic_regions+1 )
1908                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1909                   dopr_initial_index(i) = 121
1910                   hom(:,2,121,:)        = SPREAD( zu, 2, statistic_regions+1 )
1911                   hom(nzb,2,121,:)      = 0.0_wp    ! because zu(nzb) is negative
1912                   data_output_pr(i)     = data_output_pr(i)(2:)
1913                ENDIF
1914             ENDIF
1915
1916          CASE ( 'qv', '#qv' )
1917             IF ( .NOT. bulk_cloud_model ) THEN
1918                dopr_index(i) = 41
1919                dopr_unit(i)  = 'kg/kg'
1920                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
1921                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1922                   dopr_initial_index(i) = 26
1923                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
1924                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
1925                   data_output_pr(i)     = data_output_pr(i)(2:)
1926                ENDIF
1927             ELSE
1928                dopr_index(i) = 42
1929                dopr_unit(i)  = 'kg/kg'
1930                hom(:,2,42,:) = SPREAD( zu, 2, statistic_regions+1 )
1931                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1932                   dopr_initial_index(i) = 27
1933                   hom(:,2,27,:)         = SPREAD( zu, 2, statistic_regions+1 )
1934                   hom(nzb,2,27,:)       = 0.0_wp   ! because zu(nzb) is negative
1935                   data_output_pr(i)     = data_output_pr(i)(2:)
1936                ENDIF
1937             ENDIF
1938
1939          CASE ( 'thetal', '#thetal' )
1940             IF ( .NOT. bulk_cloud_model ) THEN
1941                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
1942                                 ' is not implemented for bulk_cloud_model = .FALSE.'
1943                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
1944             ELSE
1945                dopr_index(i) = 4
1946                dopr_unit(i)  = 'K'
1947                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
1948                IF ( data_output_pr(i)(1:1) == '#' )  THEN
1949                   dopr_initial_index(i) = 7
1950                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
1951                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
1952                   data_output_pr(i)     = data_output_pr(i)(2:)
1953                ENDIF
1954             ENDIF
1955
1956          CASE ( 'thetav', '#thetav' )
1957             dopr_index(i) = 44
1958             dopr_unit(i)  = 'K'
1959             hom(:,2,44,:) = SPREAD( zu, 2, statistic_regions+1 )
1960             IF ( data_output_pr(i)(1:1) == '#' )  THEN
1961                dopr_initial_index(i) = 29
1962                hom(:,2,29,:)         = SPREAD( zu, 2, statistic_regions+1 )
1963                hom(nzb,2,29,:)       = 0.0_wp    ! because zu(nzb) is negative
1964                data_output_pr(i)     = data_output_pr(i)(2:)
1965             ENDIF
1966
1967          CASE ( 'w"thetav"' )
1968             dopr_index(i) = 45
1969             dopr_unit(i)  = TRIM ( heatflux_output_unit )
1970             hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
1971
1972          CASE ( 'w*thetav*' )
1973             dopr_index(i) = 46
1974             dopr_unit(i)  = TRIM ( heatflux_output_unit )
1975             hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
1976
1977          CASE ( 'wthetav' )
1978             dopr_index(i) = 47
1979             dopr_unit(i)  = TRIM ( heatflux_output_unit )
1980             hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
1981
1982          CASE ( 'w"q"' )
1983             IF ( .NOT. humidity )  THEN
1984                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
1985                                 ' is not implemented for humidity = .FALSE.'
1986                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
1987             ELSE
1988                dopr_index(i) = 48
1989                dopr_unit(i)  = TRIM ( waterflux_output_unit )
1990                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
1991             ENDIF
1992
1993          CASE ( 'w*q*' )
1994             IF ( .NOT. humidity )  THEN
1995                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
1996                                 ' is not implemented for humidity = .FALSE.'
1997                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
1998             ELSE
1999                dopr_index(i) = 49
2000                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2001                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2002             ENDIF
2003
2004          CASE ( 'wq' )
2005             IF ( .NOT. humidity )  THEN
2006                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
2007                                 ' is not implemented for humidity = .FALSE.'
2008                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2009             ELSE
2010                dopr_index(i) = 50
2011                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2012                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2013             ENDIF
2014
2015          CASE ( 'w"s"' )
2016             IF ( .NOT. passive_scalar )  THEN
2017                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
2018                                 ' is not implemented for passive_scalar = .FALSE.'
2019                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2020             ELSE
2021                dopr_index(i) = 117
2022                dopr_unit(i)  = 'kg/m3 m/s'
2023                hom(:,2,117,:) = SPREAD( zw, 2, statistic_regions+1 )
2024             ENDIF
2025
2026          CASE ( 'w*s*' )
2027             IF ( .NOT. passive_scalar )  THEN
2028                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
2029                                 ' is not implemented for passive_scalar = .FALSE.'
2030                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2031             ELSE
2032                dopr_index(i) = 114
2033                dopr_unit(i)  = 'kg/m3 m/s'
2034                hom(:,2,114,:) = SPREAD( zw, 2, statistic_regions+1 )
2035             ENDIF
2036
2037          CASE ( 'ws' )
2038             IF ( .NOT. passive_scalar )  THEN
2039                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
2040                                 ' is not implemented for passive_scalar = .FALSE.'
2041                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2042             ELSE
2043                dopr_index(i) = 118
2044                dopr_unit(i)  = 'kg/m3 m/s'
2045                hom(:,2,118,:) = SPREAD( zw, 2, statistic_regions+1 )
2046             ENDIF
2047
2048          CASE ( 'w"qv"' )
2049             IF ( humidity  .AND.  .NOT. bulk_cloud_model )  THEN
2050                dopr_index(i) = 48
2051                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2052                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2053             ELSEIF ( humidity  .AND.  bulk_cloud_model )  THEN
2054                dopr_index(i) = 51
2055                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2056                hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
2057             ELSE
2058                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
2059                                 ' is not implemented for bulk_cloud_model = .FALSE. ' //          &
2060                                 'and humidity = .FALSE.'
2061                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2062             ENDIF
2063
2064          CASE ( 'w*qv*' )
2065             IF ( humidity  .AND.  .NOT. bulk_cloud_model )  THEN
2066                dopr_index(i) = 49
2067                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2068                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2069             ELSEIF( humidity  .AND.  bulk_cloud_model )  THEN
2070                dopr_index(i) = 52
2071                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2072                hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
2073             ELSE
2074                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
2075                                 ' is not implemented for bulk_cloud_model = .FALSE. ' //          &
2076                                 'and humidity = .FALSE.'
2077                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2078             ENDIF
2079
2080          CASE ( 'wqv' )
2081             IF ( humidity  .AND.  .NOT. bulk_cloud_model )  THEN
2082                dopr_index(i) = 50
2083                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2084                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2085             ELSEIF ( humidity  .AND.  bulk_cloud_model )  THEN
2086                dopr_index(i) = 53
2087                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2088                hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
2089             ELSE
2090                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
2091                                 ' is not implemented for bulk_cloud_model = .FALSE. ' //          &
2092                                 'and humidity = .FALSE.'
2093                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2094             ENDIF
2095
2096          CASE ( 'ql' )
2097             IF ( .NOT. bulk_cloud_model  .AND.  .NOT. cloud_droplets )  THEN
2098                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
2099                                 ' is not implemented for bulk_cloud_model = .FALSE. ' //          &
2100                                 'and cloud_droplets = .FALSE.'
2101                CALL message( 'check_parameters', 'PA0096', 1, 2, 0, 6, 0 )
2102             ELSE
2103                dopr_index(i) = 54
2104                dopr_unit(i)  = 'kg/kg'
2105                hom(:,2,54,:)  = SPREAD( zu, 2, statistic_regions+1 )
2106             ENDIF
2107
2108          CASE ( 'w*u*u*:dz' )
2109             dopr_index(i) = 55
2110             dopr_unit(i)  = 'm2/s3'
2111             hom(:,2,55,:) = SPREAD( zu, 2, statistic_regions+1 )
2112
2113          CASE ( 'w*p*:dz' )
2114             dopr_index(i) = 56
2115             dopr_unit(i)  = 'm2/s3'
2116             hom(:,2,56,:) = SPREAD( zw, 2, statistic_regions+1 )
2117
2118          CASE ( 'w"e:dz' )
2119             dopr_index(i) = 57
2120             dopr_unit(i)  = 'm2/s3'
2121             hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
2122
2123          CASE ( 'u"theta"' )
2124             dopr_index(i) = 58
2125             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2126             hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
2127
2128          CASE ( 'u*theta*' )
2129             dopr_index(i) = 59
2130             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2131             hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
2132
2133          CASE ( 'utheta_t' )
2134             dopr_index(i) = 60
2135             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2136             hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
2137
2138          CASE ( 'v"theta"' )
2139             dopr_index(i) = 61
2140             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2141             hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
2142
2143          CASE ( 'v*theta*' )
2144             dopr_index(i) = 62
2145             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2146             hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
2147
2148          CASE ( 'vtheta_t' )
2149             dopr_index(i) = 63
2150             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2151             hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
2152
2153          CASE ( 'w*p*' )
2154             dopr_index(i) = 68
2155             dopr_unit(i)  = 'm3/s3'
2156             hom(:,2,68,:) = SPREAD( zu, 2, statistic_regions+1 )
2157
2158          CASE ( 'w"e' )
2159             dopr_index(i) = 69
2160             dopr_unit(i)  = 'm3/s3'
2161             hom(:,2,69,:) = SPREAD( zu, 2, statistic_regions+1 )
2162
2163          CASE ( 'q*2' )
2164             IF (  .NOT.  humidity )  THEN
2165                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
2166                                 ' is not implemented for humidity = .FALSE.'
2167                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2168             ELSE
2169                dopr_index(i) = 70
2170                dopr_unit(i)  = 'kg2/kg2'
2171                hom(:,2,70,:) = SPREAD( zu, 2, statistic_regions+1 )
2172             ENDIF
2173
2174          CASE ( 'hyp' )
2175             dopr_index(i) = 72
2176             dopr_unit(i)  = 'hPa'
2177             hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
2178
2179          CASE ( 'rho' )
2180             dopr_index(i)  = 119
2181             dopr_unit(i)   = 'kg/m3'
2182             hom(:,2,119,:) = SPREAD( zu, 2, statistic_regions+1 )
2183
2184          CASE ( 'rho_zw' )
2185             dopr_index(i)  = 120
2186             dopr_unit(i)   = 'kg/m3'
2187             hom(:,2,120,:) = SPREAD( zw, 2, statistic_regions+1 )
2188
2189          CASE ( 'ug' )
2190             dopr_index(i) = 78
2191             dopr_unit(i)  = 'm/s'
2192             hom(:,2,78,:) = SPREAD( zu, 2, statistic_regions+1 )
2193
2194          CASE ( 'vg' )
2195             dopr_index(i) = 79
2196             dopr_unit(i)  = 'm/s'
2197             hom(:,2,79,:) = SPREAD( zu, 2, statistic_regions+1 )
2198
2199          CASE ( 'w_subs' )
2200             IF (  .NOT.  large_scale_subsidence )  THEN
2201                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
2202                                 ' is not implemented for large_scale_subsidence = .FALSE.'
2203                CALL message( 'check_parameters', 'PA0382', 1, 2, 0, 6, 0 )
2204             ELSE
2205                dopr_index(i) = 80
2206                dopr_unit(i)  = 'm/s'
2207                hom(:,2,80,:) = SPREAD( zu, 2, statistic_regions+1 )
2208             ENDIF
2209
2210          CASE ( 's*2' )
2211             IF (  .NOT.  passive_scalar )  THEN
2212                message_string = 'data_output_pr = ' // TRIM( data_output_pr(i) ) //               &
2213                                 ' is not implemented for passive_scalar = .FALSE.'
2214                CALL message( 'check_parameters', 'PA0185', 1, 2, 0, 6, 0 )
2215             ELSE
2216                dopr_index(i) = 116
2217                dopr_unit(i)  = 'kg2/m6'
2218                hom(:,2,116,:) = SPREAD( zu, 2, statistic_regions+1 )
2219             ENDIF
2220
2221           CASE ( 'eta' )
2222              dopr_index(i) = 121
2223              dopr_unit(i)  = 'mm'
2224              hom(:,2,121,:) = SPREAD( zu, 2, statistic_regions+1 )
2225
2226              kolmogorov_length_scale = .TRUE.
2227
2228          CASE DEFAULT
2229             unit = 'illegal'
2230!
2231!--          Check for other modules
2232             CALL module_interface_check_data_output_pr( data_output_pr(i), i, unit, dopr_unit(i) )
2233
2234!
2235!--          No valid quantity found
2236             IF ( unit == 'illegal' )  THEN
2237                IF ( data_output_pr_user(1) /= ' ' )  THEN
2238                   message_string = 'illegal value for data_output_pr or ' //                      &
2239                                    'data_output_pr_user = "' // TRIM( data_output_pr(i) ) // '"'
2240                   CALL message( 'check_parameters', 'PA0097', 1, 2, 0, 6, 0 )
2241                ELSE
2242                   message_string = 'illegal value for data_output_pr = "' //                      &
2243                                    TRIM( data_output_pr(i) ) // '"'
2244                   CALL message( 'check_parameters', 'PA0098', 1, 2, 0, 6, 0 )
2245                ENDIF
2246             ENDIF
2247
2248       END SELECT
2249
2250    ENDDO
2251
2252
2253!
2254!-- Append user-defined data output variables to the standard data output
2255    IF ( data_output_user(1) /= ' ' )  THEN
2256       i = 1
2257       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 500 )
2258          i = i + 1
2259       ENDDO
2260       j = 1
2261       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 500 )
2262          IF ( i > 500 )  THEN
2263             message_string = 'number of output quantitities given by data' //                     &
2264                              '_output and data_output_user exceeds the limit of 500'
2265             CALL message( 'check_parameters', 'PA0102', 1, 2, 0, 6, 0 )
2266          ENDIF
2267          data_output(i) = data_output_user(j)
2268          i = i + 1
2269          j = j + 1
2270       ENDDO
2271    ENDIF
2272!
2273!-- Check data-output variable list for any duplicates and remove them
2274    ALLOCATE( data_output_filtered(1:UBOUND( data_output, DIM=1 )) )
2275    CALL filter_duplicate_strings( varnamelength, data_output, data_output_filtered )
2276    data_output = data_output_filtered
2277    DEALLOCATE( data_output_filtered )
2278!
2279!-- Check and set steering parameters for 2d/3d data output and averaging
2280    i   = 1
2281    DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 500 )
2282!
2283!--    Check for data averaging
2284       ilen = LEN_TRIM( data_output(i) )
2285       j = 0                                                 ! no data averaging
2286       IF ( ilen > 3 )  THEN
2287          IF ( data_output(i)(ilen-2:ilen) == '_av' )  THEN
2288             j = 1                                           ! data averaging
2289             data_output(i) = data_output(i)(1:ilen-3)
2290          ENDIF
2291       ENDIF
2292!
2293!--    Check for cross section or volume data
2294       ilen = LEN_TRIM( data_output(i) )
2295       k = 0                                                   ! 3d data
2296       var = data_output(i)(1:ilen)
2297       IF ( ilen > 3 )  THEN
2298          IF ( data_output(i)(ilen-2:ilen) == '_xy'  .OR.  data_output(i)(ilen-2:ilen) == '_xz'    &
2299               .OR.  data_output(i)(ilen-2:ilen) == '_yz' )  THEN
2300             k = 1                                             ! 2d data
2301             var = data_output(i)(1:ilen-3)
2302          ENDIF
2303       ENDIF
2304
2305!
2306!--    Check for allowed value and set units
2307       SELECT CASE ( TRIM( var ) )
2308
2309          CASE ( 'e' )
2310             IF ( constant_diffusion )  THEN
2311                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
2312                                 'constant_diffusion = .FALSE.'
2313                CALL message( 'check_parameters', 'PA0103', 1, 2, 0, 6, 0 )
2314             ENDIF
2315             unit = 'm2/s2'
2316
2317          CASE ( 'thetal' )
2318             IF ( .NOT. bulk_cloud_model )  THEN
2319                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
2320                                 'bulk_cloud_model = .TRUE.'
2321                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2322             ENDIF
2323             unit = 'K'
2324
2325          CASE ( 'pc', 'pr' )
2326             IF ( .NOT. particle_advection )  THEN
2327                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
2328                                 'a "particle_parameters"-NAMELIST in the parameter file (PARIN)'
2329                CALL message( 'check_parameters', 'PA0104', 1, 2, 0, 6, 0 )
2330             ENDIF
2331             IF ( TRIM( var ) == 'pc' )  unit = 'number'
2332             IF ( TRIM( var ) == 'pr' )  unit = 'm'
2333
2334          CASE ( 'q', 'thetav' )
2335             IF ( .NOT. humidity )  THEN
2336                message_string = 'output of "' // TRIM( var ) // '" requires humidity = .TRUE.'
2337                CALL message( 'check_parameters', 'PA0105', 1, 2, 0, 6, 0 )
2338             ENDIF
2339             IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
2340             IF ( TRIM( var ) == 'thetav' )  unit = 'K'
2341
2342          CASE ( 'ql' )
2343             IF ( .NOT.  ( bulk_cloud_model  .OR.  cloud_droplets ) )  THEN
2344                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
2345                                 'bulk_cloud_model = .TRUE. or cloud_droplets = .TRUE.'
2346                CALL message( 'check_parameters', 'PA0106', 1, 2, 0, 6, 0 )
2347             ENDIF
2348             unit = 'kg/kg'
2349
2350          CASE ( 'ql_c', 'ql_v', 'ql_vp' )
2351             IF ( .NOT. cloud_droplets )  THEN
2352                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
2353                                 'cloud_droplets = .TRUE.'
2354                CALL message( 'check_parameters', 'PA0107', 1, 2, 0, 6, 0 )
2355             ENDIF
2356             IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
2357             IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
2358             IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
2359
2360          CASE ( 'qv' )
2361             IF ( .NOT. bulk_cloud_model )  THEN
2362                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
2363                                 'bulk_cloud_model = .TRUE.'
2364                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2365             ENDIF
2366             unit = 'kg/kg'
2367
2368          CASE ( 's' )
2369             IF ( .NOT. passive_scalar )  THEN
2370                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
2371                                 'passive_scalar = .TRUE.'
2372                CALL message( 'check_parameters', 'PA0110', 1, 2, 0, 6, 0 )
2373             ENDIF
2374             unit = 'kg/m3'
2375
2376          CASE ( 'p', 'theta', 'u', 'v', 'w' )
2377             IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
2378             IF ( TRIM( var ) == 'theta' )  unit = 'K'
2379             IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
2380             IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
2381             IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
2382             CONTINUE
2383
2384          CASE ( 'ghf*', 'lwp*', 'ol*', 'qsurf*', 'qsws*', 'r_a*', 'shf*', 'ssurf*', 'ssws*', 't*',&
2385                 'tsurf*', 'us*', 'z0*', 'z0h*', 'z0q*' )
2386             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
2387                message_string = 'illegal value for data_output: "' // TRIM( var ) //              &
2388                                 '" & only 2d-horizontal cross sections are allowed for this value'
2389                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
2390             ENDIF
2391
2392             IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. bulk_cloud_model )  THEN
2393                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
2394                                 'bulk_cloud_model = .TRUE.'
2395                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
2396             ENDIF
2397             IF ( TRIM( var ) == 'qsws*'  .AND.  .NOT.  humidity )  THEN
2398                message_string = 'output of "' // TRIM( var ) // '" requires humidity = .TRUE.'
2399                CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
2400             ENDIF
2401
2402             IF ( TRIM( var ) == 'ghf*'  .AND.  .NOT.  land_surface )  THEN
2403                message_string = 'output of "' // TRIM( var ) // '" requires land_surface = .TRUE.'
2404                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
2405             ENDIF
2406
2407             IF ( (  TRIM( var ) == 'r_a*' .OR. TRIM( var ) == 'ghf*' )  .AND.  .NOT. land_surface &
2408                   .AND.  .NOT. urban_surface )  THEN
2409                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
2410                                 'land_surface = .TRUE. or ' // 'urban_surface = .TRUE.'
2411                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
2412             ENDIF
2413
2414             IF ( TRIM( var ) == 'ssws*'  .AND.  .NOT. passive_scalar )  THEN
2415                message_string = 'output of "' // TRIM( var ) // '" requires ' //                  &
2416                                 'passive_scalar = .TRUE.'
2417                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
2418             ENDIF
2419
2420             IF ( TRIM( var ) == 'ghf*'   )  unit = 'W/m2'
2421             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/m2'
2422             IF ( TRIM( var ) == 'ol*'    )  unit = 'm'
2423             IF ( TRIM( var ) == 'qsurf*' )  unit = 'kg/kg'
2424             IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
2425             IF ( TRIM( var ) == 'r_a*'   )  unit = 's/m'
2426             IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
2427             IF ( TRIM( var ) == 'ssurf*' )  unit = 'kg/kg'
2428             IF ( TRIM( var ) == 'ssws*'  )  unit = 'kg/m2*s'
2429             IF ( TRIM( var ) == 't*'     )  unit = 'K'
2430             IF ( TRIM( var ) == 'tsurf*' )  unit = 'K'
2431             IF ( TRIM( var ) == 'us*'    )  unit = 'm/s'
2432             IF ( TRIM( var ) == 'z0*'    )  unit = 'm'
2433             IF ( TRIM( var ) == 'z0h*'   )  unit = 'm'
2434!
2435!--          Output of surface latent and sensible heat flux will be in W/m2 in case of natural- and
2436!--          urban-type surfaces, even if flux_output_mode is set to kinematic units.
2437             IF ( land_surface  .OR.  urban_surface )  THEN
2438                IF ( TRIM( var ) == 'shf*'   )  unit = 'W/m2'
2439                IF ( TRIM( var ) == 'qsws*'  )  unit = 'W/m2'
2440             ENDIF
2441
2442          CASE DEFAULT
2443!
2444!--          Check for other modules
2445             CALL module_interface_check_data_output( var, unit, i, j, ilen, k )
2446
2447             IF ( unit == 'illegal' )  THEN
2448                IF ( data_output_user(1) /= ' ' )  THEN
2449                   message_string = 'illegal value for data_output or ' //                         &
2450                                    'data_output_user = "' // TRIM( data_output(i) ) // '"'
2451                   CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
2452                ELSE
2453                   message_string = 'illegal value for data_output = "' //                         &
2454                                    TRIM( data_output(i) ) // '"'
2455                   CALL message( 'check_parameters', 'PA0115', 1, 2, 0, 6, 0 )
2456                ENDIF
2457             ENDIF
2458
2459       END SELECT
2460!
2461!--    Set the internal steering parameters appropriately
2462       IF ( k == 0 )  THEN
2463          do3d_no(j)              = do3d_no(j) + 1
2464          do3d(j,do3d_no(j))      = data_output(i)
2465          do3d_unit(j,do3d_no(j)) = unit
2466       ELSE
2467          do2d_no(j)              = do2d_no(j) + 1
2468          do2d(j,do2d_no(j))      = data_output(i)
2469          do2d_unit(j,do2d_no(j)) = unit
2470          IF ( data_output(i)(ilen-2:ilen) == '_xy' )  THEN
2471             data_output_xy(j) = .TRUE.
2472          ENDIF
2473          IF ( data_output(i)(ilen-2:ilen) == '_xz' )  THEN
2474             data_output_xz(j) = .TRUE.
2475          ENDIF
2476          IF ( data_output(i)(ilen-2:ilen) == '_yz' )  THEN
2477             data_output_yz(j) = .TRUE.
2478          ENDIF
2479       ENDIF
2480
2481       IF ( j == 1 )  THEN
2482!
2483!--       Check, if variable is already subject to averaging
2484          found = .FALSE.
2485          DO  k = 1, doav_n
2486             IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
2487          ENDDO
2488
2489          IF ( .NOT. found )  THEN
2490             doav_n = doav_n + 1
2491             doav(doav_n) = var
2492          ENDIF
2493       ENDIF
2494
2495       i = i + 1
2496    ENDDO
2497
2498!
2499!-- Averaged 2d or 3d output requires that an averaging interval has been set
2500    IF ( doav_n > 0  .AND.  averaging_interval == 0.0_wp )  THEN
2501       WRITE( message_string, * )  'output of averaged quantity "', TRIM( doav(1) ),               &
2502                                   '_av" requires to set a ', 'non-zero averaging interval'
2503       CALL message( 'check_parameters', 'PA0323', 1, 2, 0, 6, 0 )
2504    ENDIF
2505
2506!
2507!-- Check sectional planes and store them in one shared array
2508    IF ( ANY( section_xy > nz + 1 ) )  THEN
2509       WRITE( message_string, * )  'section_xy must be <= nz + 1 = ', nz + 1
2510       CALL message( 'check_parameters', 'PA0319', 1, 2, 0, 6, 0 )
2511    ENDIF
2512    IF ( ANY( section_xz > ny + 1 ) )  THEN
2513       WRITE( message_string, * )  'section_xz must be <= ny + 1 = ', ny + 1
2514       CALL message( 'check_parameters', 'PA0320', 1, 2, 0, 6, 0 )
2515    ENDIF
2516    IF ( ANY( section_yz > nx + 1 ) )  THEN
2517       WRITE( message_string, * )  'section_yz must be <= nx + 1 = ', nx + 1
2518       CALL message( 'check_parameters', 'PA0321', 1, 2, 0, 6, 0 )
2519    ENDIF
2520    section(:,1) = section_xy
2521    section(:,2) = section_xz
2522    section(:,3) = section_yz
2523
2524    IF ( ANY( data_output_xy )  .AND.  .NOT. ANY( section(:,1) /= -9999 ) )  THEN
2525       WRITE( message_string, * )  'section_xy not defined for requested xy-cross section ' //     &
2526                                   'output.&At least one cross section must be given.'
2527       CALL message( 'check_parameters', 'PA0681', 1, 2, 0, 6, 0 )
2528    ENDIF
2529    IF ( ANY( data_output_xz )  .AND.  .NOT. ANY( section(:,2) /= -9999 ) )  THEN
2530       WRITE( message_string, * )  'section_xz not defined for requested xz-cross section ' //     &
2531                                   'output.&At least one cross section must be given.'
2532       CALL message( 'check_parameters', 'PA0681', 1, 2, 0, 6, 0 )
2533    ENDIF
2534    IF ( ANY( data_output_yz )  .AND.  .NOT. ANY( section(:,3) /= -9999 ) )  THEN
2535       WRITE( message_string, * )  'section_yz not defined for requested yz-cross section ' //     &
2536                                   'output.&At least one cross section must be given.'
2537       CALL message( 'check_parameters', 'PA0681', 1, 2, 0, 6, 0 )
2538    ENDIF
2539!
2540!-- Upper plot limit for 3D arrays
2541    IF ( nz_do3d == -9999 )  nz_do3d = nzt + 1
2542
2543!
2544!-- Set output format string (used in header)
2545    SELECT CASE ( netcdf_data_format )
2546       CASE ( 1 )
2547          netcdf_data_format_string = 'netCDF classic'
2548       CASE ( 2 )
2549          netcdf_data_format_string = 'netCDF 64bit offset'
2550       CASE ( 3 )
2551          netcdf_data_format_string = 'netCDF4/HDF5'
2552       CASE ( 4 )
2553          netcdf_data_format_string = 'netCDF4/HDF5 classic'
2554       CASE ( 5 )
2555          netcdf_data_format_string = 'parallel netCDF4/HDF5'
2556       CASE ( 6 )
2557          netcdf_data_format_string = 'parallel netCDF4/HDF5 classic'
2558
2559    END SELECT
2560
2561!
2562!-- Check mask conditions
2563    DO mid = 1, max_masks
2564       IF ( data_output_masks(mid,1) /= ' '  .OR.  data_output_masks_user(mid,1) /= ' ' )  THEN
2565          masks = masks + 1
2566       ENDIF
2567    ENDDO
2568
2569    IF ( masks < 0  .OR.  masks > max_masks )  THEN
2570       WRITE( message_string, * )  'illegal value: masks must be >= 0 and ', '<= ', max_masks,     &
2571              ' (=max_masks)'
2572       CALL message( 'check_parameters', 'PA0325', 1, 2, 0, 6, 0 )
2573    ENDIF
2574    IF ( masks > 0 )  THEN
2575       mask_scale(1) = mask_scale_x
2576       mask_scale(2) = mask_scale_y
2577       mask_scale(3) = mask_scale_z
2578       IF ( ANY( mask_scale <= 0.0_wp ) )  THEN
2579          WRITE( message_string, * )  'illegal value: mask_scale_x, mask_scale_y and mask_scale_z',&
2580                 'must be > 0.0'
2581          CALL message( 'check_parameters', 'PA0326', 1, 2, 0, 6, 0 )
2582       ENDIF
2583!
2584!--    Generate masks for masked data output.
2585!--    Parallel netcdf output is not tested so far for masked data, hence netcdf_data_format is
2586!--    switched back to non-parallel output.
2587       netcdf_data_format_save = netcdf_data_format
2588       IF ( netcdf_data_format > 4 )  THEN
2589          IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
2590          IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
2591          message_string = 'netCDF file formats '// '5 (parallel netCDF 4) and 6 (parallel ' //    &
2592                           'netCDF 4 Classic model) & are currently not supported (not yet ' //    &
2593                           'tested) for masked data. &Using respective non-parallel' //            &
2594                           ' output for masked data.'
2595          CALL message( 'check_parameters', 'PA0383', 0, 0, 0, 6, 0 )
2596       ENDIF
2597       CALL init_masks
2598       netcdf_data_format = netcdf_data_format_save
2599    ENDIF
2600
2601!
2602!-- Check the NetCDF data format
2603    IF ( netcdf_data_format > 2 )  THEN
2604#if defined( __netcdf4 )
2605       CONTINUE
2606#else
2607       message_string = 'netCDF: netCDF4 format requested but no ' //                              &
2608                        'cpp-directive __netcdf4 given & switch back to 64-bit offset format'
2609       CALL message( 'check_parameters', 'PA0171', 0, 1, 0, 6, 0 )
2610       netcdf_data_format = 2
2611#endif
2612    ENDIF
2613    IF ( netcdf_data_format > 4 )  THEN
2614#if defined( __netcdf4 ) && defined( __netcdf4_parallel )
2615       CONTINUE
2616#else
2617       message_string = 'netCDF: netCDF4 parallel output requested but no ' //                     &
2618                        'cpp-directive __netcdf4_parallel given & switch ' //                      &
2619                        'back to netCDF4 non-parallel output'
2620       CALL message( 'check_parameters', 'PA0099', 0, 1, 0, 6, 0 )
2621       netcdf_data_format = netcdf_data_format - 2
2622#endif
2623    ENDIF
2624
2625!
2626!-- Calculate fixed number of output time levels for parallel netcdf output.
2627!-- The time dimension has to be defined as limited for parallel output, because otherwise the I/O
2628!-- performance drops significantly.
2629    IF ( netcdf_data_format > 4 )  THEN
2630
2631!
2632!--    Check if any of the follwoing data output interval is 0.0s, which is not allowed for parallel
2633!--    output.
2634       CALL check_dt_do( dt_do3d,           'dt_do3d'           )
2635       CALL check_dt_do( dt_do2d_xy,        'dt_do2d_xy'        )
2636       CALL check_dt_do( dt_do2d_xz,        'dt_do2d_xz'        )
2637       CALL check_dt_do( dt_do2d_yz,        'dt_do2d_yz'        )
2638       CALL check_dt_do( dt_data_output_av, 'dt_data_output_av' )
2639
2640!--    Set needed time levels (ntdim) to saved time levels + to be saved time levels.
2641       ntdim_3d(0) = do3d_time_count(0) + CEILING(                                                 &
2642                     ( end_time - MAX(                                                             &
2643                         MERGE( skip_time_do3d, skip_time_do3d + spinup_time,                      &
2644                                data_output_during_spinup ),                                       &
2645                         simulated_time_at_begin )                                                 &
2646                     ) / dt_do3d )
2647       IF ( do3d_at_begin ) ntdim_3d(0) = ntdim_3d(0) + 1
2648
2649       ntdim_3d(1) = do3d_time_count(1) + CEILING(                                                 &
2650                     ( end_time - MAX(                                                             &
2651                         MERGE( skip_time_data_output_av, skip_time_data_output_av + spinup_time,  &
2652                                data_output_during_spinup ),                                       &
2653                         simulated_time_at_begin )                                                 &
2654                     ) / dt_data_output_av )
2655
2656       ntdim_2d_xy(0) = do2d_xy_time_count(0) + CEILING(                                           &
2657                        ( end_time - MAX(                                                          &
2658                            MERGE( skip_time_do2d_xy, skip_time_do2d_xy + spinup_time,             &
2659                                   data_output_during_spinup ),                                    &
2660                            simulated_time_at_begin )                                              &
2661                        ) / dt_do2d_xy )
2662
2663       ntdim_2d_xz(0) = do2d_xz_time_count(0) + CEILING(                                           &
2664                        ( end_time - MAX(                                                          &
2665                            MERGE( skip_time_do2d_xz, skip_time_do2d_xz + spinup_time,             &
2666                                   data_output_during_spinup ),                                    &
2667                            simulated_time_at_begin )                                              &
2668                        ) / dt_do2d_xz )
2669
2670       ntdim_2d_yz(0) = do2d_yz_time_count(0) + CEILING(                                           &
2671                        ( end_time - MAX(                                                          &
2672                            MERGE( skip_time_do2d_yz, skip_time_do2d_yz + spinup_time,             &
2673                                   data_output_during_spinup ),                                    &
2674                            simulated_time_at_begin )                                              &
2675                        ) / dt_do2d_yz )
2676
2677       IF ( do2d_at_begin )  THEN
2678          ntdim_2d_xy(0) = ntdim_2d_xy(0) + 1
2679          ntdim_2d_xz(0) = ntdim_2d_xz(0) + 1
2680          ntdim_2d_yz(0) = ntdim_2d_yz(0) + 1
2681       ENDIF
2682!
2683!--    Please note, for averaged 2D data skip_time_data_output_av is the relavant output control
2684!--    parameter.
2685       ntdim_2d_xy(1) = do2d_xy_time_count(1) + CEILING(                                           &
2686                        ( end_time - MAX( MERGE( skip_time_data_output_av,                         &
2687                                                 skip_time_data_output_av + spinup_time,           &
2688                                                 data_output_during_spinup ),                      &
2689                                          simulated_time_at_begin )                                &
2690                        ) / dt_data_output_av )
2691
2692       ntdim_2d_xz(1) = do2d_xz_time_count(1) + CEILING(                                           &
2693                        ( end_time - MAX( MERGE( skip_time_data_output_av,                         &
2694                                                 skip_time_data_output_av + spinup_time,           &
2695                                                 data_output_during_spinup ),                      &
2696                                          simulated_time_at_begin )                                &
2697                        ) / dt_data_output_av )
2698
2699       ntdim_2d_yz(1) = do2d_yz_time_count(1) + CEILING(                                           &
2700                        ( end_time - MAX( MERGE( skip_time_data_output_av,                         &
2701                                                 skip_time_data_output_av + spinup_time,           &
2702                                                 data_output_during_spinup ),                      &
2703                                          simulated_time_at_begin )                                &
2704                        ) / dt_data_output_av )
2705
2706    ENDIF
2707
2708!
2709!-- Check, whether a constant diffusion coefficient shall be used
2710    IF ( km_constant /= -1.0_wp )  THEN
2711       IF ( km_constant < 0.0_wp )  THEN
2712          WRITE( message_string, * )  'km_constant = ', km_constant, ' < 0.0'
2713          CALL message( 'check_parameters', 'PA0121', 1, 2, 0, 6, 0 )
2714       ELSE
2715          IF ( prandtl_number < 0.0_wp )  THEN
2716             WRITE( message_string, * )  'prandtl_number = ', prandtl_number, ' < 0.0'
2717             CALL message( 'check_parameters', 'PA0122', 1, 2, 0, 6, 0 )
2718          ENDIF
2719          constant_diffusion = .TRUE.
2720
2721          IF ( constant_flux_layer )  THEN
2722             message_string = 'constant_flux_layer is not allowed with fixed value of km'
2723             CALL message( 'check_parameters', 'PA0123', 1, 2, 0, 6, 0 )
2724          ENDIF
2725       ENDIF
2726    ENDIF
2727
2728!
2729!-- In case of non-cyclic lateral boundaries and a damping layer for the potential temperature,
2730!-- check the width of the damping layer
2731    IF ( bc_lr /= 'cyclic' ) THEN
2732       IF ( pt_damping_width < 0.0_wp  .OR.  pt_damping_width > REAL( (nx+1) * dx ) )  THEN
2733          message_string = 'pt_damping_width out of range'
2734          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
2735       ENDIF
2736    ENDIF
2737
2738    IF ( bc_ns /= 'cyclic' )  THEN
2739       IF ( pt_damping_width < 0.0_wp  .OR.  pt_damping_width > REAL( (ny+1) * dy ) )  THEN
2740          message_string = 'pt_damping_width out of range'
2741          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
2742       ENDIF
2743    ENDIF
2744
2745!
2746!-- Check value range for zeta = z/L
2747    IF ( zeta_min >= zeta_max )  THEN
2748       WRITE( message_string, * )  'zeta_min = ', zeta_min, ' must be less ', 'than zeta_max = ',  &
2749              zeta_max
2750       CALL message( 'check_parameters', 'PA0125', 1, 2, 0, 6, 0 )
2751    ENDIF
2752
2753!
2754!-- Check random generator
2755    IF ( (random_generator /= 'system-specific'      .AND.                                         &
2756          random_generator /= 'random-parallel'   )  .AND.                                         &
2757          random_generator /= 'numerical-recipes' )  THEN
2758       message_string = 'unknown random generator: random_generator = "' //                        &
2759                        TRIM( random_generator ) // '"'
2760       CALL message( 'check_parameters', 'PA0135', 1, 2, 0, 6, 0 )
2761    ENDIF
2762
2763!
2764!-- Determine upper and lower hight level indices for random perturbations
2765    IF ( disturbance_level_b == -9999999.9_wp )  THEN
2766       IF ( ocean_mode )  THEN
2767          disturbance_level_b     = zu((nzt*2)/3)
2768          disturbance_level_ind_b = ( nzt * 2 ) / 3
2769       ELSE
2770          disturbance_level_b     = zu(nzb+3)
2771          disturbance_level_ind_b = nzb + 3
2772       ENDIF
2773    ELSEIF ( disturbance_level_b < zu(3) )  THEN
2774       WRITE( message_string, * )  'disturbance_level_b = ', disturbance_level_b, ' must be >= ',  &
2775              zu(3), '(zu(3))'
2776       CALL message( 'check_parameters', 'PA0126', 1, 2, 0, 6, 0 )
2777    ELSEIF ( disturbance_level_b > zu(nzt-2) )  THEN
2778       WRITE( message_string, * )  'disturbance_level_b = ', disturbance_level_b, ' must be <= ',  &
2779              zu(nzt-2), '(zu(nzt-2))'
2780       CALL message( 'check_parameters', 'PA0127', 1, 2, 0, 6, 0 )
2781    ELSE
2782       DO  k = 3, nzt-2
2783          IF ( disturbance_level_b <= zu(k) )  THEN
2784             disturbance_level_ind_b = k
2785             EXIT
2786          ENDIF
2787       ENDDO
2788    ENDIF
2789
2790    IF ( disturbance_level_t == -9999999.9_wp )  THEN
2791       IF ( ocean_mode )  THEN
2792          disturbance_level_t     = zu(nzt-3)
2793          disturbance_level_ind_t = nzt - 3
2794       ELSE
2795          disturbance_level_t     = zu(nzt/3)
2796          disturbance_level_ind_t = nzt / 3
2797       ENDIF
2798    ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
2799       WRITE( message_string, * )  'disturbance_level_t = ', disturbance_level_t, ' must be <= ',  &
2800              zu(nzt-2), '(zu(nzt-2))'
2801       CALL message( 'check_parameters', 'PA0128', 1, 2, 0, 6, 0 )
2802    ELSEIF ( disturbance_level_t < disturbance_level_b )  THEN
2803       WRITE( message_string, * )  'disturbance_level_t = ', disturbance_level_t,                  &
2804             ' must be >= disturbance_level_b = ', disturbance_level_b
2805       CALL message( 'check_parameters', 'PA0129', 1, 2, 0, 6, 0 )
2806    ELSE
2807       DO  k = 3, nzt-2
2808          IF ( disturbance_level_t <= zu(k) )  THEN
2809             disturbance_level_ind_t = k
2810             EXIT
2811          ENDIF
2812       ENDDO
2813    ENDIF
2814
2815!
2816!-- Check again whether the levels determined this way are ok.
2817!-- Error may occur at automatic determination and too few grid points in z-direction.
2818    IF ( disturbance_level_ind_t < disturbance_level_ind_b )  THEN
2819       WRITE( message_string, * )  'disturbance_level_ind_t = ', disturbance_level_ind_t,          &
2820              ' must be >= ', 'disturbance_level_ind_b = ', disturbance_level_ind_b
2821       CALL message( 'check_parameters', 'PA0130', 1, 2, 0, 6, 0 )
2822    ENDIF
2823
2824!
2825!-- Determine the horizontal index range for random perturbations.
2826!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed near the inflow and
2827!-- the perturbation area is further limited to ...(1) after the initial phase of the flow.
2828
2829    IF ( bc_lr /= 'cyclic' )  THEN
2830       IF ( inflow_disturbance_begin == -1 )  THEN
2831          inflow_disturbance_begin = MIN( 10, nx/2 )
2832       ENDIF
2833       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > nx )  THEN
2834          message_string = 'inflow_disturbance_begin out of range'
2835          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
2836       ENDIF
2837       IF ( inflow_disturbance_end == -1 )  THEN
2838          inflow_disturbance_end = MIN( 100, 3*nx/4 )
2839       ENDIF
2840       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > nx )  THEN
2841          message_string = 'inflow_disturbance_end out of range'
2842          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
2843       ENDIF
2844    ELSEIF ( bc_ns /= 'cyclic' )  THEN
2845       IF ( inflow_disturbance_begin == -1 )  THEN
2846          inflow_disturbance_begin = MIN( 10, ny/2 )
2847       ENDIF
2848       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > ny )  THEN
2849          message_string = 'inflow_disturbance_begin out of range'
2850          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
2851       ENDIF
2852       IF ( inflow_disturbance_end == -1 )  THEN
2853          inflow_disturbance_end = MIN( 100, 3*ny/4 )
2854       ENDIF
2855       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > ny )  THEN
2856          message_string = 'inflow_disturbance_end out of range'
2857          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
2858       ENDIF
2859    ENDIF
2860
2861    IF ( random_generator == 'random-parallel' )  THEN
2862       dist_nxl = nxl;  dist_nxr = nxr
2863       dist_nys = nys;  dist_nyn = nyn
2864       IF ( bc_lr == 'radiation/dirichlet' )  THEN
2865          dist_nxr    = MIN( nx - inflow_disturbance_begin, nxr )
2866          dist_nxl(1) = MAX( nx - inflow_disturbance_end, nxl )
2867       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
2868          dist_nxl    = MAX( inflow_disturbance_begin, nxl )
2869          dist_nxr(1) = MIN( inflow_disturbance_end, nxr )
2870       ELSEIF ( bc_lr == 'nested'  .OR.  bc_lr == 'nesting_offline' )  THEN
2871          dist_nxl    = MAX( inflow_disturbance_begin, nxl )
2872          dist_nxr    = MIN( nx - inflow_disturbance_begin, nxr )
2873       ENDIF
2874       IF ( bc_ns == 'dirichlet/radiation' )  THEN
2875          dist_nyn    = MIN( ny - inflow_disturbance_begin, nyn )
2876          dist_nys(1) = MAX( ny - inflow_disturbance_end, nys )
2877       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
2878          dist_nys    = MAX( inflow_disturbance_begin, nys )
2879          dist_nyn(1) = MIN( inflow_disturbance_end, nyn )
2880       ELSEIF ( bc_ns == 'nested'  .OR.  bc_ns == 'nesting_offline' )  THEN
2881          dist_nys    = MAX( inflow_disturbance_begin, nys )
2882          dist_nyn    = MIN( ny - inflow_disturbance_begin, nyn )
2883       ENDIF
2884    ELSE
2885       dist_nxl = 0;  dist_nxr = nx
2886       dist_nys = 0;  dist_nyn = ny
2887       IF ( bc_lr == 'radiation/dirichlet' )  THEN
2888          dist_nxr    = nx - inflow_disturbance_begin
2889          dist_nxl(1) = nx - inflow_disturbance_end
2890       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
2891          dist_nxl    = inflow_disturbance_begin
2892          dist_nxr(1) = inflow_disturbance_end
2893       ELSEIF ( bc_lr == 'nested'  .OR.  bc_lr == 'nesting_offline' )  THEN
2894          dist_nxr    = nx - inflow_disturbance_begin
2895          dist_nxl    = inflow_disturbance_begin
2896       ENDIF
2897       IF ( bc_ns == 'dirichlet/radiation' )  THEN
2898          dist_nyn    = ny - inflow_disturbance_begin
2899          dist_nys(1) = ny - inflow_disturbance_end
2900       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
2901          dist_nys    = inflow_disturbance_begin
2902          dist_nyn(1) = inflow_disturbance_end
2903       ELSEIF ( bc_ns == 'nested'  .OR.  bc_ns == 'nesting_offline' )  THEN
2904          dist_nyn    = ny - inflow_disturbance_begin
2905          dist_nys    = inflow_disturbance_begin
2906       ENDIF
2907    ENDIF
2908
2909!
2910!-- A turbulent inflow requires Dirichlet conditions at the respective inflow boundary (so far, a
2911!-- turbulent inflow is realized from the left side only).
2912    IF ( turbulent_inflow  .AND.  bc_lr /= 'dirichlet/radiation' )  THEN
2913       message_string = 'turbulent_inflow = .T. requires a Dirichlet ' //                          &
2914                        'condition at the inflow boundary'
2915       CALL message( 'check_parameters', 'PA0133', 1, 2, 0, 6, 0 )
2916    ENDIF
2917
2918!
2919!-- Turbulent inflow requires that 3d arrays have been cyclically filled with data from prerun in
2920!-- the first main run
2921    IF ( turbulent_inflow  .AND.  initializing_actions /= 'cyclic_fill'  .AND.                     &
2922         initializing_actions /= 'read_restart_data' )  THEN
2923       message_string = 'turbulent_inflow = .T. requires ' //                                      &
2924                        'initializing_actions = ''cyclic_fill'' or ' //                            &
2925                        'initializing_actions = ''read_restart_data'' '
2926       CALL message( 'check_parameters', 'PA0055', 1, 2, 0, 6, 0 )
2927    ENDIF
2928
2929!
2930!-- In case of turbulent inflow
2931    IF ( turbulent_inflow )  THEN
2932
2933!
2934!--    Calculate the index of the recycling plane
2935       IF ( recycling_width <= dx  .OR.  recycling_width >= nx * dx )  THEN
2936          WRITE( message_string, * )  'illegal value for recycling_width: ', recycling_width
2937          CALL message( 'check_parameters', 'PA0134', 1, 2, 0, 6, 0 )
2938       ENDIF
2939!
2940!--    Calculate the index
2941       recycling_plane = recycling_width / dx
2942!
2943!--   Check for correct input of recycling method for thermodynamic quantities
2944       IF ( TRIM( recycling_method_for_thermodynamic_quantities ) /= 'turbulent_fluctuation'  .AND.&
2945            TRIM( recycling_method_for_thermodynamic_quantities ) /= 'absolute_value' )  THEN
2946          WRITE( message_string, * )  'unknown recycling method for thermodynamic quantities: ',   &
2947               TRIM( recycling_method_for_thermodynamic_quantities )
2948          CALL message( 'check_parameters', 'PA0184', 1, 2, 0, 6, 0 )
2949       ENDIF
2950
2951    ENDIF
2952
2953
2954    IF ( turbulent_outflow )  THEN
2955!
2956!--    Turbulent outflow requires Dirichlet conditions at the respective inflow boundary (so far, a
2957!--    turbulent outflow is realized at the right side only).
2958       IF ( bc_lr /= 'dirichlet/radiation' )  THEN
2959          message_string = 'turbulent_outflow = .T. requires bc_lr = "dirichlet/radiation"'
2960          CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )
2961       ENDIF
2962!
2963!--    The ouflow-source plane must lay inside the model domain
2964       IF ( outflow_source_plane < dx  .OR.  outflow_source_plane > nx * dx )  THEN
2965          WRITE( message_string, * )  'illegal value for outflow_source'// '_plane: ',             &
2966                 outflow_source_plane
2967          CALL message( 'check_parameters', 'PA0145', 1, 2, 0, 6, 0 )
2968       ENDIF
2969    ENDIF
2970
2971!
2972!-- Determine damping level index for 1D model
2973    IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
2974       IF ( damp_level_1d == -1.0_wp )  THEN
2975          damp_level_1d     = zu(nzt+1)
2976          damp_level_ind_1d = nzt + 1
2977       ELSEIF ( damp_level_1d < 0.0_wp  .OR.  damp_level_1d > zu(nzt+1) )  THEN
2978          WRITE( message_string, * )  'damp_level_1d = ', damp_level_1d, ' must be >= 0.0 and <= ',&
2979                 zu(nzt+1), '(zu(nzt+1))'
2980          CALL message( 'check_parameters', 'PA0136', 1, 2, 0, 6, 0 )
2981       ELSE
2982          DO  k = 1, nzt+1
2983             IF ( damp_level_1d <= zu(k) )  THEN
2984                damp_level_ind_1d = k
2985                EXIT
2986             ENDIF
2987          ENDDO
2988       ENDIF
2989    ENDIF
2990
2991!
2992!-- Check some other 1d-model parameters
2993    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                                       &
2994         TRIM( mixing_length_1d ) /= 'blackadar' )  THEN
2995       message_string = 'mixing_length_1d = "' // TRIM( mixing_length_1d ) // '" is unknown'
2996       CALL message( 'check_parameters', 'PA0137', 1, 2, 0, 6, 0 )
2997    ENDIF
2998    IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model'  .AND.                                         &
2999         TRIM( dissipation_1d ) /= 'detering'        .AND.                                         &
3000         TRIM( dissipation_1d ) /= 'prognostic' )  THEN
3001       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) // '" is unknown'
3002       CALL message( 'check_parameters', 'PA0138', 1, 2, 0, 6, 0 )
3003    ENDIF
3004    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                                       &
3005         TRIM( dissipation_1d ) == 'as_in_3d_model' )  THEN
3006       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) //                          &
3007                        '" requires mixing_length_1d = "as_in_3d_model"'
3008       CALL message( 'check_parameters', 'PA0485', 1, 2, 0, 6, 0 )
3009    ENDIF
3010
3011!
3012!-- Set time for the next user defined restart (time_restart is the internal parameter for steering
3013!-- restart events)
3014    IF ( restart_time /= 9999999.9_wp )  THEN
3015       IF ( restart_time > time_since_reference_point )  THEN
3016          time_restart = restart_time
3017       ENDIF
3018    ELSE
3019!
3020!--    In case of a restart run, set internal parameter to default (no restart) if the
3021!--    NAMELIST-parameter restart_time is omitted
3022       time_restart = 9999999.9_wp
3023    ENDIF
3024
3025!
3026!-- Check pressure gradient conditions
3027    IF ( dp_external  .AND.  conserve_volume_flow )  THEN
3028       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo',                     &
3029              'w are .TRUE. but one of them must be .FALSE.'
3030       CALL message( 'check_parameters', 'PA0150', 1, 2, 0, 6, 0 )
3031    ENDIF
3032    IF ( dp_external )  THEN
3033       IF ( dp_level_b < zu(nzb)  .OR.  dp_level_b > zu(nzt) )  THEN
3034          WRITE( message_string, * )  'dp_level_b = ', dp_level_b, ' is out ',                     &
3035                 ' of range [zu(nzb), zu(nzt)]'
3036          CALL message( 'check_parameters', 'PA0151', 1, 2, 0, 6, 0 )
3037       ENDIF
3038       IF ( .NOT. ANY( dpdxy /= 0.0_wp ) )  THEN
3039          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is ze',                     &
3040                 'ro, i.e. the external pressure gradient will not be applied'
3041          CALL message( 'check_parameters', 'PA0152', 0, 1, 0, 6, 0 )
3042       ENDIF
3043    ENDIF
3044    IF ( ANY( dpdxy /= 0.0_wp )  .AND.  .NOT.  dp_external )  THEN
3045       WRITE( message_string, * )  'dpdxy is nonzero but dp_external is ',                         &
3046              '.FALSE., i.e. the external pressure gradient & will not be applied'
3047       CALL message( 'check_parameters', 'PA0153', 0, 1, 0, 6, 0 )
3048    ENDIF
3049    IF ( conserve_volume_flow )  THEN
3050       IF ( TRIM( conserve_volume_flow_mode ) == 'default' )  THEN
3051
3052          conserve_volume_flow_mode = 'initial_profiles'
3053
3054       ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles'  .AND.                     &
3055                TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' )  THEN
3056          WRITE( message_string, * )  'unknown conserve_volume_flow_mode: ',                       &
3057                 conserve_volume_flow_mode
3058          CALL message( 'check_parameters', 'PA0154', 1, 2, 0, 6, 0 )
3059       ENDIF
3060       IF ( ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic')  .AND.                                   &
3061            TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
3062          WRITE( message_string, * )  'non-cyclic boundary conditions ',                           &
3063                 'require  conserve_volume_flow_mode = ''initial_profiles'''
3064          CALL message( 'check_parameters', 'PA0155', 1, 2, 0, 6, 0 )
3065       ENDIF
3066    ENDIF
3067    IF ( ( u_bulk /= 0.0_wp .OR. v_bulk /= 0.0_wp )  .AND.                                         &
3068         ( .NOT. conserve_volume_flow .OR. TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) )&
3069    THEN
3070       WRITE( message_string, * )  'nonzero bulk velocity requires ',                              &
3071              'conserve_volume_flow = .T. and ', 'conserve_volume_flow_mode = ''bulk_velocity'''
3072       CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 )
3073    ENDIF
3074
3075!
3076!-- Prevent empty time records in volume, cross-section and masked data in case of non-parallel
3077!-- netcdf-output in restart runs
3078    IF ( netcdf_data_format < 5 )  THEN
3079       IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
3080          do3d_time_count    = 0
3081          do2d_xy_time_count = 0
3082          do2d_xz_time_count = 0
3083          do2d_yz_time_count = 0
3084          domask_time_count  = 0
3085       ENDIF
3086    ENDIF
3087
3088
3089!
3090!-- Check roughness length, which has to be smaller than dz/2
3091    IF ( ( constant_flux_layer .OR. INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  &
3092         .AND.  roughness_length >= 0.5 * dz(1) )  THEN
3093       message_string = 'roughness_length must be smaller than dz/2'
3094       CALL message( 'check_parameters', 'PA0424', 1, 2, 0, 6, 0 )
3095    ENDIF
3096
3097!
3098!-- Check if topography is read from file in case of complex terrain simulations
3099    IF ( complex_terrain  .AND.  TRIM( topography ) /= 'read_from_file' )  THEN
3100       message_string = 'complex_terrain requires topography = ''read_from_file'''
3101       CALL message( 'check_parameters', 'PA0295', 1, 2, 0, 6, 0 )
3102    ENDIF
3103
3104!
3105!-- Check if vertical grid stretching is switched off in case of complex terrain simulations
3106    IF ( complex_terrain  .AND.  ANY( dz_stretch_level_start /= -9999999.9_wp ) )  THEN
3107       message_string = 'Vertical grid stretching is not allowed for complex_terrain = .TRUE.'
3108       CALL message( 'check_parameters', 'PA0473', 1, 2, 0, 6, 0 )
3109    ENDIF
3110
3111    CALL location_message( 'checking parameters', 'finished' )
3112
3113 CONTAINS
3114
3115!--------------------------------------------------------------------------------------------------!
3116! Description:
3117! ------------
3118!> Check the length of data output intervals. In case of parallel NetCDF output the time levels of
3119!> the output files need to be fixed. Therefore setting the output interval to 0.0s (usually used to
3120!> output each timestep) is not possible as long as a non-fixed timestep is used.
3121!--------------------------------------------------------------------------------------------------!
3122
3123    SUBROUTINE check_dt_do( dt_do, dt_do_name )
3124
3125       IMPLICIT NONE
3126
3127       CHARACTER (LEN=*), INTENT (IN) :: dt_do_name !< parin variable name
3128
3129       REAL(wp), INTENT (INOUT)       :: dt_do      !< data output interval
3130
3131       IF ( dt_do == 0.0_wp )  THEN
3132          IF ( dt_fixed )  THEN
3133             WRITE( message_string, '(A,F9.4,A)' )  'Output at every timestep is wanted (' //      &
3134                    dt_do_name // ' = 0.0).&'//                                                    &
3135                    'The output interval is set to the fixed timestep dt '// '= ', dt, 's.'
3136             CALL message( 'check_parameters', 'PA0060', 0, 0, 0, 6, 0 )
3137             dt_do = dt
3138          ELSE
3139             message_string = dt_do_name // ' = 0.0 while using a ' //                             &
3140                              'variable timestep and parallel netCDF4 is not allowed.'
3141             CALL message( 'check_parameters', 'PA0081', 1, 2, 0, 6, 0 )
3142          ENDIF
3143       ENDIF
3144
3145    END SUBROUTINE check_dt_do
3146
3147
3148
3149!--------------------------------------------------------------------------------------------------!
3150! Description:
3151! ------------
3152!> Set the bottom and top boundary conditions for humidity and scalars.
3153!--------------------------------------------------------------------------------------------------!
3154
3155    SUBROUTINE set_bc_scalars( sq, bc_b, bc_t, ibc_b, ibc_t, err_nr_b, err_nr_t )
3156
3157
3158       IMPLICIT NONE
3159
3160       CHARACTER (LEN=*)   ::  bc_b       !< bottom boundary condition
3161       CHARACTER (LEN=*)   ::  bc_t       !< top boundary condition
3162       CHARACTER (LEN=*)   ::  err_nr_b   !< error number if bottom bc is unknown
3163       CHARACTER (LEN=*)   ::  err_nr_t   !< error number if top bc is unknown
3164       CHARACTER (LEN=1)   ::  sq         !< name of scalar quantity
3165
3166
3167       INTEGER(iwp)        ::  ibc_b      !< index for bottom boundary condition
3168       INTEGER(iwp)        ::  ibc_t      !< index for top boundary condition
3169
3170!
3171!--    Set Integer flags and check for possilbe errorneous settings for bottom boundary condition
3172       IF ( bc_b == 'dirichlet' )  THEN
3173          ibc_b = 0
3174       ELSEIF ( bc_b == 'neumann' )  THEN
3175          ibc_b = 1
3176       ELSE
3177          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // '_b ="' //           &
3178                           TRIM( bc_b ) // '"'
3179          CALL message( 'check_parameters', err_nr_b, 1, 2, 0, 6, 0 )
3180       ENDIF
3181!
3182!--    Set Integer flags and check for possilbe errorneous settings for top boundary condition
3183       IF ( bc_t == 'dirichlet' )  THEN
3184          ibc_t = 0
3185       ELSEIF ( bc_t == 'neumann' )  THEN
3186          ibc_t = 1
3187       ELSEIF ( bc_t == 'initial_gradient' )  THEN
3188          ibc_t = 2
3189       ELSEIF ( bc_t == 'nested'  .OR.  bc_t == 'nesting_offline' )  THEN
3190          ibc_t = 3
3191       ELSE
3192          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) // '_t ="' //           &
3193                           TRIM( bc_t ) // '"'
3194          CALL message( 'check_parameters', err_nr_t, 1, 2, 0, 6, 0 )
3195       ENDIF
3196
3197
3198    END SUBROUTINE set_bc_scalars
3199
3200
3201
3202!--------------------------------------------------------------------------------------------------!
3203! Description:
3204! ------------
3205!> Check for consistent settings of bottom boundary conditions for humidity and scalars.
3206!--------------------------------------------------------------------------------------------------!
3207
3208    SUBROUTINE check_bc_scalars( sq, bc_b, ibc_b, err_nr_1, err_nr_2, constant_flux,               &
3209                                 surface_initial_change )
3210
3211
3212       IMPLICIT NONE
3213
3214       CHARACTER (LEN=*)   ::  bc_b                     !< bottom boundary condition
3215       CHARACTER (LEN=*)   ::  err_nr_1                 !< error number of first error
3216       CHARACTER (LEN=*)   ::  err_nr_2                 !< error number of second error
3217       CHARACTER (LEN=1)   ::  sq                       !< name of scalar quantity
3218
3219
3220       INTEGER(iwp)        ::  ibc_b                    !< index of bottom boundary condition
3221
3222       LOGICAL             ::  constant_flux            !< flag for constant-flux layer
3223
3224       REAL(wp)            ::  surface_initial_change   !< value of initial change at the surface
3225
3226!
3227!--    A given surface value implies Dirichlet boundary condition for the respective quantity. In
3228!--    this case specification of a constant flux is forbidden. However, an exception is made for
3229!--    large-scale forcing as well as land-surface model.
3230       IF ( .NOT. land_surface  .AND.  .NOT. large_scale_forcing )  THEN
3231          IF ( ibc_b == 0  .AND.  constant_flux )  THEN
3232             message_string = 'boundary condition: bc_' // TRIM( sq ) // '_b = "' //               &
3233                              TRIM( bc_b ) // '" is not allowed with prescribed surface flux'
3234             CALL message( 'check_parameters', err_nr_1, 1, 2, 0, 6, 0 )
3235          ENDIF
3236       ENDIF
3237       IF ( constant_flux  .AND.  surface_initial_change /= 0.0_wp )  THEN
3238          WRITE( message_string, * )  'a prescribed surface flux is not allo', 'wed with ', sq,    &
3239                '_surface_initial_change (/=0) = ', surface_initial_change
3240          CALL message( 'check_parameters', err_nr_2, 1, 2, 0, 6, 0 )
3241       ENDIF
3242
3243
3244    END SUBROUTINE check_bc_scalars
3245
3246
3247
3248!--------------------------------------------------------------------------------------------------!
3249! Description:
3250! ------------
3251!> Filter a list of strings for duplicates. The output list is of same length but does not contain
3252!> any duplicates.
3253!--------------------------------------------------------------------------------------------------!
3254
3255    SUBROUTINE filter_duplicate_strings( string_length, string_list, string_list_filtered )
3256
3257       IMPLICIT NONE
3258
3259       INTEGER, INTENT(IN) :: string_length  !< length of a single string in list
3260
3261       CHARACTER (LEN=*),             DIMENSION(:), INTENT(IN)               ::  string_list           !< non-filtered list
3262       CHARACTER (LEN=string_length), DIMENSION(:), INTENT(OUT), ALLOCATABLE ::  string_list_filtered  !< filtered list
3263
3264       CHARACTER (LEN=string_length) ::  string_to_check  !< string to be check for duplicates
3265
3266       INTEGER(iwp) ::  i         !< loop index
3267       INTEGER(iwp) ::  j         !< loop index
3268       INTEGER(iwp) ::  k         !< loop index
3269       INTEGER(iwp) ::  nstrings  !< number of strings in list
3270
3271       LOGICAL ::  is_duplicate  !< true if string has duplicates in list
3272
3273
3274       nstrings = UBOUND( string_list, DIM=1 )
3275       ALLOCATE( string_list_filtered(1:nstrings) )
3276       string_list_filtered(:) = ' '
3277
3278       k = 0
3279       DO  i = 1, nstrings
3280          string_to_check = string_list(i)
3281
3282          IF ( string_to_check == ' ' )  EXIT
3283
3284          is_duplicate = .FALSE.
3285          DO  j = 1, nstrings
3286             IF ( string_list_filtered(j) == ' ' )  THEN
3287                EXIT
3288             ELSEIF ( string_to_check == string_list_filtered(j) )  THEN
3289                is_duplicate = .TRUE.
3290                EXIT
3291             ENDIF
3292          ENDDO
3293
3294          IF ( .NOT. is_duplicate )  THEN
3295             k = k + 1
3296             string_list_filtered(k) = string_to_check
3297          ENDIF
3298
3299       ENDDO
3300
3301    END SUBROUTINE filter_duplicate_strings
3302
3303
3304 END SUBROUTINE check_parameters
Note: See TracBrowser for help on using the repository browser.