Home

Context Navigation

source: palm/trunk/SOURCE/dynamics_mod.f90 @ 4844

Last change on this file since 4844 was 4843, checked in by raasch, 4 years ago
local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates
Property svn:keywords set to `Id`
File size: 72.2 KB

Line
1	!> @file dynamics_mod.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: dynamics_mod.f90 4843 2021-01-15 15:22:11Z raasch $
26	! local namelist parameter added to switch off the module although the respective module namelist
27	! appears in the namelist file
28	!
29	! 4842 2021-01-14 10:42:28Z raasch
30	! reading of namelist file and actions in case of namelist errors revised so that statement labels
31	! and goto statements are not required any more
32	!
33	! 4828 2021-01-05 11:21:41Z Giersch
34	! Enable 3D data output also with 64-bit precision
35	!
36	! 4760 2020-10-26 13:26:47Z schwenkel
37	! Implement relative humidity as diagnostic output quantity
38	!
39	! 4731 2020-10-07 13:25:11Z schwenkel
40	! Move exchange_horiz from time_integration to modules
41	!
42	! 4627 2020-07-26 10:14:44Z raasch
43	! bugfix for r4626
44	!
45	! 4626 2020-07-26 09:49:48Z raasch
46	! file re-formatted to follow the PALM coding standard
47	!
48	! 4517 2020-05-03 14:29:30Z raasch
49	! added restart with MPI-IO for reading local arrays
50	!
51	! 4505 2020-04-20 15:37:15Z schwenkel
52	! Add flag for saturation check
53	!
54	! 4495 2020-04-13 20:11:20Z resler
55	! restart data handling with MPI-IO added
56	!
57	! 4360 2020-01-07 11:25:50Z suehring
58	! Bugfix for last commit.
59	!
60	! 4359 2019-12-30 13:36:50Z suehring
61	! Refine post-initialization check for realistically inital values of mixing ratio. Give an error
62	! message for faulty initial values, but only a warning in a restart run.
63	!
64	! 4347 2019-12-18 13:18:33Z suehring
65	! Implement post-initialization check for realistically inital values of mixing ratio
66	!
67	! 4281 2019-10-29 15:15:39Z schwenkel
68	! Moved boundary conditions in dynamics module
69	!
70	! 4097 2019-07-15 11:59:11Z suehring
71	! Avoid overlong lines - limit is 132 characters per line
72	!
73	! 4047 2019-06-21 18:58:09Z knoop
74	! Initial introduction of the dynamics module with only dynamics_swap_timelevel implemented
75	!
76	!
77	! Description:
78	! ------------
79	!> This module contains the dynamics of PALM.
80	!--------------------------------------------------------------------------------------------------!
81	MODULE dynamics_mod
82
83
84	USE arrays_3d, &
85	ONLY: c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l, &
86	diss, &
87	diss_p, &
88	dzu, &
89	e, &
90	e_p, &
91	exner, &
92	hyp, &
93	pt, pt_1, pt_2, pt_init, pt_p, &
94	q, q_1, q_2, q_p, &
95	s, s_1, s_2, s_p, &
96	u, u_1, u_2, u_init, u_p, u_m_l, u_m_n, u_m_r, u_m_s, &
97	v, v_1, v_2, v_p, v_init, v_m_l, v_m_n, v_m_r, v_m_s, &
98	w, w_1, w_2, w_p, w_m_l, w_m_n, w_m_r, w_m_s, &
99	zu
100
101	USE basic_constants_and_equations_mod, &
102	ONLY: magnus, &
103	rd_d_rv
104
105	USE control_parameters, &
106	ONLY: bc_dirichlet_l, &
107	bc_dirichlet_s, &
108	bc_radiation_l, &
109	bc_radiation_n, &
110	bc_radiation_r, &
111	bc_radiation_s, &
112	bc_pt_t_val, &
113	bc_q_t_val, &
114	bc_s_t_val, &
115	check_realistic_q, &
116	child_domain, &
117	coupling_mode, &
118	constant_diffusion, &
119	dt_3d, &
120	humidity, &
121	ibc_pt_b, &
122	ibc_pt_t, &
123	ibc_q_b, &
124	ibc_q_t, &
125	ibc_s_b, &
126	ibc_s_t, &
127	ibc_uv_b, &
128	ibc_uv_t, &
129	initializing_actions, &
130	intermediate_timestep_count, &
131	length, &
132	message_string, &
133	nesting_offline, &
134	neutral, &
135	nudging, &
136	passive_scalar, &
137	restart_string, &
138	rans_mode, &
139	rans_tke_e, &
140	tsc, &
141	use_cmax
142
143	USE exchange_horiz_mod, &
144	ONLY: exchange_horiz
145
146
147	USE grid_variables, &
148	ONLY: ddx, &
149	ddy, &
150	dx, &
151	dy
152
153	USE indices, &
154	ONLY: nbgp, &
155	nx, &
156	nxl, &
157	nxlg, &
158	nxr, &
159	nxrg, &
160	ny, &
161	nys, &
162	nysg, &
163	nyn, &
164	nyng, &
165	nzb, &
166	nzt
167
168	USE kinds
169
170	USE pegrid
171
172	USE pmc_interface, &
173	ONLY : nesting_mode
174
175	! USE restart_data_mpi_io_mod, &
176	! ONLY:
177
178	USE surface_mod, &
179	ONLY : bc_h
180
181
182	IMPLICIT NONE
183
184	LOGICAL :: dynamics_module_enabled = .FALSE. !<
185
186	SAVE
187
188	PRIVATE
189
190	!
191	!-- Public functions
192	PUBLIC &
193	dynamics_parin, &
194	dynamics_check_parameters, &
195	dynamics_check_data_output_ts, &
196	dynamics_check_data_output_pr, &
197	dynamics_check_data_output, &
198	dynamics_init_masks, &
199	dynamics_define_netcdf_grid, &
200	dynamics_init_arrays, &
201	dynamics_init, &
202	dynamics_init_checks, &
203	dynamics_header, &
204	dynamics_actions, &
205	dynamics_non_advective_processes, &
206	dynamics_exchange_horiz, &
207	dynamics_prognostic_equations, &
208	dynamics_boundary_conditions, &
209	dynamics_swap_timelevel, &
210	dynamics_3d_data_averaging, &
211	dynamics_data_output_2d, &
212	dynamics_data_output_3d, &
213	dynamics_statistics, &
214	dynamics_rrd_global, &
215	dynamics_rrd_local, &
216	dynamics_wrd_global, &
217	dynamics_wrd_local, &
218	dynamics_last_actions
219
220	!
221	!-- Public parameters, constants and initial values
222	PUBLIC &
223	dynamics_module_enabled
224
225	INTERFACE dynamics_parin
226	MODULE PROCEDURE dynamics_parin
227	END INTERFACE dynamics_parin
228
229	INTERFACE dynamics_check_parameters
230	MODULE PROCEDURE dynamics_check_parameters
231	END INTERFACE dynamics_check_parameters
232
233	INTERFACE dynamics_check_data_output_ts
234	MODULE PROCEDURE dynamics_check_data_output_ts
235	END INTERFACE dynamics_check_data_output_ts
236
237	INTERFACE dynamics_check_data_output_pr
238	MODULE PROCEDURE dynamics_check_data_output_pr
239	END INTERFACE dynamics_check_data_output_pr
240
241	INTERFACE dynamics_check_data_output
242	MODULE PROCEDURE dynamics_check_data_output
243	END INTERFACE dynamics_check_data_output
244
245	INTERFACE dynamics_init_masks
246	MODULE PROCEDURE dynamics_init_masks
247	END INTERFACE dynamics_init_masks
248
249	INTERFACE dynamics_define_netcdf_grid
250	MODULE PROCEDURE dynamics_define_netcdf_grid
251	END INTERFACE dynamics_define_netcdf_grid
252
253	INTERFACE dynamics_init_arrays
254	MODULE PROCEDURE dynamics_init_arrays
255	END INTERFACE dynamics_init_arrays
256
257	INTERFACE dynamics_init
258	MODULE PROCEDURE dynamics_init
259	END INTERFACE dynamics_init
260
261	INTERFACE dynamics_init_checks
262	MODULE PROCEDURE dynamics_init_checks
263	END INTERFACE dynamics_init_checks
264
265	INTERFACE dynamics_header
266	MODULE PROCEDURE dynamics_header
267	END INTERFACE dynamics_header
268
269	INTERFACE dynamics_actions
270	MODULE PROCEDURE dynamics_actions
271	MODULE PROCEDURE dynamics_actions_ij
272	END INTERFACE dynamics_actions
273
274	INTERFACE dynamics_non_advective_processes
275	MODULE PROCEDURE dynamics_non_advective_processes
276	MODULE PROCEDURE dynamics_non_advective_processes_ij
277	END INTERFACE dynamics_non_advective_processes
278
279	INTERFACE dynamics_exchange_horiz
280	MODULE PROCEDURE dynamics_exchange_horiz
281	END INTERFACE dynamics_exchange_horiz
282
283	INTERFACE dynamics_prognostic_equations
284	MODULE PROCEDURE dynamics_prognostic_equations
285	MODULE PROCEDURE dynamics_prognostic_equations_ij
286	END INTERFACE dynamics_prognostic_equations
287
288	INTERFACE dynamics_boundary_conditions
289	MODULE PROCEDURE dynamics_boundary_conditions
290	END INTERFACE dynamics_boundary_conditions
291
292	INTERFACE dynamics_swap_timelevel
293	MODULE PROCEDURE dynamics_swap_timelevel
294	END INTERFACE dynamics_swap_timelevel
295
296	INTERFACE dynamics_3d_data_averaging
297	MODULE PROCEDURE dynamics_3d_data_averaging
298	END INTERFACE dynamics_3d_data_averaging
299
300	INTERFACE dynamics_data_output_2d
301	MODULE PROCEDURE dynamics_data_output_2d
302	END INTERFACE dynamics_data_output_2d
303
304	INTERFACE dynamics_data_output_3d
305	MODULE PROCEDURE dynamics_data_output_3d
306	END INTERFACE dynamics_data_output_3d
307
308	INTERFACE dynamics_statistics
309	MODULE PROCEDURE dynamics_statistics
310	END INTERFACE dynamics_statistics
311
312	INTERFACE dynamics_rrd_global
313	MODULE PROCEDURE dynamics_rrd_global_ftn
314	MODULE PROCEDURE dynamics_rrd_global_mpi
315	END INTERFACE dynamics_rrd_global
316
317	INTERFACE dynamics_rrd_local
318	MODULE PROCEDURE dynamics_rrd_local_ftn
319	MODULE PROCEDURE dynamics_rrd_local_mpi
320	END INTERFACE dynamics_rrd_local
321
322	INTERFACE dynamics_wrd_global
323	MODULE PROCEDURE dynamics_wrd_global
324	END INTERFACE dynamics_wrd_global
325
326	INTERFACE dynamics_wrd_local
327	MODULE PROCEDURE dynamics_wrd_local
328	END INTERFACE dynamics_wrd_local
329
330	INTERFACE dynamics_last_actions
331	MODULE PROCEDURE dynamics_last_actions
332	END INTERFACE dynamics_last_actions
333
334
335	CONTAINS
336
337
338	!--------------------------------------------------------------------------------------------------!
339	! Description:
340	! ------------
341	!> Read module-specific namelist
342	!--------------------------------------------------------------------------------------------------!
343	SUBROUTINE dynamics_parin
344
345	CHARACTER(LEN=100) :: line !< dummy string that contains the current line of the parameter
346	!< file
347	INTEGER(iwp) :: io_status !< status after reading the namelist file
348
349	LOGICAL :: switch_off_module = .FALSE. !< local namelist parameter to switch off the module
350	!< although the respective module namelist appears in
351	!< the namelist file
352
353	NAMELIST /dynamics_parameters/ switch_off_module
354
355
356	!
357	!-- For the time beeing (unless the dynamics module is further developed), set default module
358	!-- switch to true.
359	dynamics_module_enabled = .TRUE.
360
361	!
362	!-- Move to the beginning of the namelist file and try to find and read the namelist.
363	REWIND( 11 )
364	READ( 11, dynamics_parameters, IOSTAT=io_status )
365
366	!
367	!-- Action depending on the READ status
368	IF ( io_status == 0 ) THEN
369	!
370	!-- dynamics_parameters namelist was found and read correctly.
371	IF ( .NOT. switch_off_module ) dynamics_module_enabled = .TRUE.
372
373	ELSEIF ( io_status > 0 ) THEN
374	!
375	!-- dynamics_parameters namelist was found, but contained errors. Print an error message
376	!-- including the line that caused the problem.
377	BACKSPACE( 11 )
378	READ( 11 , '(A)') line
379	CALL parin_fail_message( 'dynamics_parameters', line )
380
381	ENDIF
382
383	END SUBROUTINE dynamics_parin
384
385
386	!--------------------------------------------------------------------------------------------------!
387	! Description:
388	! ------------
389	!> Check control parameters and deduce further quantities.
390	!--------------------------------------------------------------------------------------------------!
391	SUBROUTINE dynamics_check_parameters
392
393
394	END SUBROUTINE dynamics_check_parameters
395
396
397	!--------------------------------------------------------------------------------------------------!
398	! Description:
399	! ------------
400	!> Set module-specific timeseries units and labels
401	!--------------------------------------------------------------------------------------------------!
402	SUBROUTINE dynamics_check_data_output_ts( dots_max, dots_num, dots_label, dots_unit )
403
404	INTEGER(iwp), INTENT(IN) :: dots_max
405
406	CHARACTER (LEN=*), DIMENSION(dots_max), INTENT(INOUT) :: dots_label
407	CHARACTER (LEN=*), DIMENSION(dots_max), INTENT(INOUT) :: dots_unit
408
409	INTEGER(iwp), INTENT(INOUT) :: dots_num
410
411	!
412	!-- Next line is to avoid compiler warning about unused variables. Please remove.
413	IF ( dots_num == 0 .OR. dots_label(1)(1:1) == ' ' .OR. dots_unit(1)(1:1) == ' ' ) CONTINUE
414
415
416	END SUBROUTINE dynamics_check_data_output_ts
417
418
419	!--------------------------------------------------------------------------------------------------!
420	! Description:
421	! ------------
422	!> Set the unit of module-specific profile output quantities. For those variables not recognized,
423	!> the parameter unit is set to "illegal", which tells the calling routine that the output variable
424	!> is not defined and leads to a program abort.
425	!--------------------------------------------------------------------------------------------------!
426	SUBROUTINE dynamics_check_data_output_pr( variable, var_count, unit, dopr_unit )
427
428
429	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
430	CHARACTER (LEN=*) :: unit !<
431	CHARACTER (LEN=*) :: variable !<
432
433	INTEGER(iwp) :: var_count !<
434
435	!
436	!-- Next line is to avoid compiler warning about unused variables. Please remove.
437	IF ( unit(1:1) == ' ' .OR. dopr_unit(1:1) == ' ' .OR. var_count == 0 ) CONTINUE
438
439	SELECT CASE ( TRIM( variable ) )
440
441	! CASE ( 'var_name' )
442
443	CASE DEFAULT
444	unit = 'illegal'
445
446	END SELECT
447
448
449	END SUBROUTINE dynamics_check_data_output_pr
450
451
452	!--------------------------------------------------------------------------------------------------!
453	! Description:
454	! ------------
455	!> Set the unit of module-specific output quantities. For those variables not recognized,
456	!> the parameter unit is set to "illegal", which tells the calling routine that the output variable
457	!< is not defined and leads to a program abort.
458	!--------------------------------------------------------------------------------------------------!
459	SUBROUTINE dynamics_check_data_output( variable, unit )
460
461
462	CHARACTER (LEN=*) :: unit !<
463	CHARACTER (LEN=*) :: variable !<
464
465	SELECT CASE ( TRIM( variable ) )
466
467	! CASE ( 'u2' )
468
469	CASE DEFAULT
470	unit = 'illegal'
471
472	END SELECT
473
474
475	END SUBROUTINE dynamics_check_data_output
476
477
478	!--------------------------------------------------------------------------------------------------!
479	!
480	! Description:
481	! ------------
482	!> Initialize module-specific masked output
483	!--------------------------------------------------------------------------------------------------!
484	SUBROUTINE dynamics_init_masks( variable, unit )
485
486
487	CHARACTER (LEN=*) :: unit !<
488	CHARACTER (LEN=*) :: variable !<
489
490
491	SELECT CASE ( TRIM( variable ) )
492
493	! CASE ( 'u2' )
494
495	CASE DEFAULT
496	unit = 'illegal'
497
498	END SELECT
499
500
501	END SUBROUTINE dynamics_init_masks
502
503
504	!--------------------------------------------------------------------------------------------------!
505	! Description:
506	! ------------
507	!> Initialize module-specific arrays
508	!--------------------------------------------------------------------------------------------------!
509	SUBROUTINE dynamics_init_arrays
510
511
512	END SUBROUTINE dynamics_init_arrays
513
514
515	!--------------------------------------------------------------------------------------------------!
516	! Description:
517	! ------------
518	!> Execution of module-specific initializing actions
519	!--------------------------------------------------------------------------------------------------!
520	SUBROUTINE dynamics_init
521
522
523	END SUBROUTINE dynamics_init
524
525
526	!--------------------------------------------------------------------------------------------------!
527	! Description:
528	! ------------
529	!> Perform module-specific post-initialization checks
530	!--------------------------------------------------------------------------------------------------!
531	SUBROUTINE dynamics_init_checks
532
533	INTEGER(iwp) :: i !< loop index in x-direction
534	INTEGER(iwp) :: j !< loop index in y-direction
535	INTEGER(iwp) :: k !< loop index in z-direction
536
537	LOGICAL :: realistic_q = .TRUE. !< flag indicating realistic mixing ratios
538
539	REAL(wp) :: e_s !< saturation water vapor pressure
540	REAL(wp) :: q_s !< saturation mixing ratio
541	REAL(wp) :: t_l !< actual temperature
542	REAL(wp) :: rh_check = 9999999.9_wp !< relative humidity
543	REAL(wp) :: rh_min = 9999999.9_wp !< max relative humidity
544	REAL(wp) :: height = 9999999.9_wp !< height of supersaturated regions
545	REAL(wp) :: min_height = 9999999.9_wp !< height of supersaturated regions
546
547	!
548	!-- Check for realistic initial mixing ratio. This must be in a realistic phyiscial range and must
549	!-- not exceed the saturation mixing ratio by more than 2 percent. Please note, the check is
550	!-- performed for each grid point (not just for a vertical profile), in order to cover also
551	!-- three-dimensional initialization. Note, this check gives an error only for the initial run not
552	!-- for a restart run. In case there are no cloud physics considered, the mixing ratio can exceed
553	!-- the saturation moisture. This case a warning is given.
554	IF ( humidity .AND. .NOT. neutral .AND. check_realistic_q ) THEN
555	DO i = nxl, nxr
556	DO j = nys, nyn
557	DO k = nzb+1, nzt
558	!
559	!-- Calculate actual temperature, water vapor saturation pressure, and based on this
560	!-- the saturation mixing ratio.
561	t_l = exner(k) * pt(k,j,i)
562	e_s = magnus( t_l )
563	q_s = rd_d_rv * e_s / ( hyp(k) - e_s )
564
565	IF ( q(k,j,i) > 1.02_wp * q_s ) THEN
566	realistic_q = .FALSE.
567	rh_check = q(k,j,i) / q_s * 100.0_wp
568	height = zu(k)
569	ENDIF
570	ENDDO
571	ENDDO
572	ENDDO
573	!
574	!-- Since the check is performed locally, merge the logical flag from all mpi ranks,
575	!-- in order to do not print the error message multiple times.
576	#if defined( __parallel )
577	CALL MPI_ALLREDUCE( MPI_IN_PLACE, realistic_q, 1, MPI_LOGICAL, MPI_LAND, comm2d, ierr)
578	CALL MPI_ALLREDUCE( rh_check, rh_min, 1, MPI_REAL, MPI_MIN, comm2d, ierr )
579	CALL MPI_ALLREDUCE( height, min_height, 1, MPI_REAL, MPI_MIN, comm2d, ierr )
580	#endif
581
582	IF ( .NOT. realistic_q .AND. &
583	TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
584	WRITE( message_string, * ) 'The initial mixing ratio exceeds the saturation mixing' // &
585	'ratio, with rh =', rh_min, '% in a height of', min_height, 'm for the first time'
586	CALL message( 'dynamic_init_checks', 'PA0697', 2, 2, 0, 6, 0 )
587	ELSEIF ( .NOT. realistic_q .AND. &
588	TRIM( initializing_actions ) == 'read_restart_data' ) THEN
589	WRITE( message_string, * ) 'The initial mixing ratio exceeds the saturation mixing' // &
590	'ratio, with rh =', rh_min, '% in a height of', min_height, 'm for the first time'
591	CALL message( 'dynamic_init_checks', 'PA0697', 0, 1, 0, 6, 0 )
592	ENDIF
593	ENDIF
594
595	END SUBROUTINE dynamics_init_checks
596
597
598	!--------------------------------------------------------------------------------------------------!
599	! Description:
600	! ------------
601	!> Set the grids on which module-specific output quantities are defined. Allowed values for
602	!> grid_x are "x" and "xu", for grid_y "y" and "yv", and for grid_z "zu" and "zw".
603	!--------------------------------------------------------------------------------------------------!
604	SUBROUTINE dynamics_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
605
606
607	CHARACTER (LEN=*) :: grid_x !<
608	CHARACTER (LEN=*) :: grid_y !<
609	CHARACTER (LEN=*) :: grid_z !<
610	CHARACTER (LEN=*) :: variable !<
611
612	LOGICAL :: found !<
613
614
615	SELECT CASE ( TRIM( variable ) )
616
617	! CASE ( 'u2' )
618
619	CASE DEFAULT
620	found = .FALSE.
621	grid_x = 'none'
622	grid_y = 'none'
623	grid_z = 'none'
624
625	END SELECT
626
627
628	END SUBROUTINE dynamics_define_netcdf_grid
629
630
631	!--------------------------------------------------------------------------------------------------!
632	! Description:
633	! ------------
634	!> Print a header with module-specific information.
635	!--------------------------------------------------------------------------------------------------!
636	SUBROUTINE dynamics_header( io )
637
638
639	INTEGER(iwp) :: io !<
640
641	!
642	!-- If no module-specific variables are read from the namelist-file, no information will be printed.
643	IF ( .NOT. dynamics_module_enabled ) THEN
644	WRITE ( io, 100 )
645	RETURN
646	ENDIF
647
648	!
649	!-- Printing the information.
650	WRITE ( io, 110 )
651
652	!
653	!-- Format-descriptors
654	100 FORMAT (//' *** dynamic module disabled'/)
655	110 FORMAT (//1X,78('#') &
656	//' User-defined variables and actions:'/ &
657	' -----------------------------------'//)
658
659	END SUBROUTINE dynamics_header
660
661
662	!--------------------------------------------------------------------------------------------------!
663	! Description:
664	! ------------
665	!> Execute module-specific actions for all grid points
666	!--------------------------------------------------------------------------------------------------!
667	SUBROUTINE dynamics_actions( location )
668
669
670	CHARACTER (LEN=*) :: location !<
671
672	! INTEGER(iwp) :: i !<
673	! INTEGER(iwp) :: j !<
674	! INTEGER(iwp) :: k !<
675
676	!
677	!-- Here the user-defined actions follow
678	!-- No calls for single grid points are allowed at locations before and after the timestep, since
679	!-- these calls are not within an i,j-loop
680	SELECT CASE ( location )
681
682	CASE ( 'before_timestep' )
683
684
685	CASE ( 'before_prognostic_equations' )
686
687
688	CASE ( 'after_integration' )
689
690
691	CASE ( 'after_timestep' )
692
693
694	CASE ( 'u-tendency' )
695
696
697	CASE ( 'v-tendency' )
698
699
700	CASE ( 'w-tendency' )
701
702
703	CASE ( 'pt-tendency' )
704
705
706	CASE ( 'sa-tendency' )
707
708
709	CASE ( 'e-tendency' )
710
711
712	CASE ( 'q-tendency' )
713
714
715	CASE ( 's-tendency' )
716
717
718	CASE DEFAULT
719	CONTINUE
720
721	END SELECT
722
723	END SUBROUTINE dynamics_actions
724
725
726	!--------------------------------------------------------------------------------------------------!
727	! Description:
728	! ------------
729	!> Execute module-specific actions for grid point i,j
730	!--------------------------------------------------------------------------------------------------!
731	SUBROUTINE dynamics_actions_ij( i, j, location )
732
733
734	CHARACTER (LEN=*) :: location
735
736	INTEGER(iwp) :: i
737	INTEGER(iwp) :: j
738
739	!
740	!-- Here the user-defined actions follow
741	SELECT CASE ( location )
742
743	CASE ( 'u-tendency' )
744
745	!
746	!-- Next line is to avoid compiler warning about unused variables. Please remove.
747	IF ( i + j < 0 ) CONTINUE
748
749	CASE ( 'v-tendency' )
750
751
752	CASE ( 'w-tendency' )
753
754
755	CASE ( 'pt-tendency' )
756
757
758	CASE ( 'sa-tendency' )
759
760
761	CASE ( 'e-tendency' )
762
763
764	CASE ( 'q-tendency' )
765
766
767	CASE ( 's-tendency' )
768
769
770	CASE DEFAULT
771	CONTINUE
772
773	END SELECT
774
775	END SUBROUTINE dynamics_actions_ij
776
777
778	!--------------------------------------------------------------------------------------------------!
779	! Description:
780	! ------------
781	!> Compute module-specific non-advective processes for all grid points
782	!--------------------------------------------------------------------------------------------------!
783	SUBROUTINE dynamics_non_advective_processes
784
785
786
787	END SUBROUTINE dynamics_non_advective_processes
788
789
790	!--------------------------------------------------------------------------------------------------!
791	! Description:
792	! ------------
793	!> Compute module-specific non-advective processes for grid points i,j
794	!--------------------------------------------------------------------------------------------------!
795	SUBROUTINE dynamics_non_advective_processes_ij( i, j )
796
797
798	INTEGER(iwp) :: i !<
799	INTEGER(iwp) :: j !<
800
801	!
802	!-- Next line is just to avoid compiler warnings about unused variables. You may remove it.
803	IF ( i + j < 0 ) CONTINUE
804
805
806	END SUBROUTINE dynamics_non_advective_processes_ij
807
808
809	!--------------------------------------------------------------------------------------------------!
810	! Description:
811	! ------------
812	!> Perform module-specific horizontal boundary exchange
813	!--------------------------------------------------------------------------------------------------!
814	SUBROUTINE dynamics_exchange_horiz( location )
815
816	CHARACTER (LEN=*), INTENT(IN) :: location !< call location string
817
818	SELECT CASE ( location )
819
820	CASE ( 'before_prognostic_equation' )
821
822	CASE ( 'after_prognostic_equation' )
823
824	CALL exchange_horiz( u_p, nbgp )
825	CALL exchange_horiz( v_p, nbgp )
826	CALL exchange_horiz( w_p, nbgp )
827	CALL exchange_horiz( pt_p, nbgp )
828	IF ( .NOT. constant_diffusion ) CALL exchange_horiz( e_p, nbgp )
829	IF ( rans_tke_e ) CALL exchange_horiz( diss_p, nbgp )
830	IF ( humidity ) CALL exchange_horiz( q_p, nbgp )
831	IF ( passive_scalar ) CALL exchange_horiz( s_p, nbgp )
832
833	CASE ( 'after_anterpolation' )
834
835	CALL exchange_horiz( u, nbgp )
836	CALL exchange_horiz( v, nbgp )
837	CALL exchange_horiz( w, nbgp )
838	IF ( .NOT. neutral ) CALL exchange_horiz( pt, nbgp )
839	IF ( humidity ) CALL exchange_horiz( q, nbgp )
840	IF ( passive_scalar ) CALL exchange_horiz( s, nbgp )
841	IF ( .NOT. constant_diffusion ) CALL exchange_horiz( e, nbgp )
842	IF ( .NOT. constant_diffusion .AND. rans_mode .AND. rans_tke_e ) THEN
843	CALL exchange_horiz( diss, nbgp )
844	ENDIF
845
846	END SELECT
847
848	END SUBROUTINE dynamics_exchange_horiz
849
850
851	!--------------------------------------------------------------------------------------------------!
852	! Description:
853	! ------------
854	!> Compute module-specific prognostic equations for all grid points
855	!--------------------------------------------------------------------------------------------------!
856	SUBROUTINE dynamics_prognostic_equations
857
858
859
860	END SUBROUTINE dynamics_prognostic_equations
861
862
863	!--------------------------------------------------------------------------------------------------!
864	! Description:
865	! ------------
866	!> Compute module-specific prognostic equations for grid point i,j
867	!--------------------------------------------------------------------------------------------------!
868	SUBROUTINE dynamics_prognostic_equations_ij( i, j, i_omp_start, tn )
869
870
871	INTEGER(iwp), INTENT(IN) :: i !< grid index in x-direction
872	INTEGER(iwp), INTENT(IN) :: i_omp_start !< first loop index of i-loop in prognostic_equations
873	INTEGER(iwp), INTENT(IN) :: j !< grid index in y-direction
874	INTEGER(iwp), INTENT(IN) :: tn !< task number of openmp task
875
876	!
877	!-- Next line is just to avoid compiler warnings about unused variables. You may remove it.
878	IF ( i + j + i_omp_start + tn < 0 ) CONTINUE
879
880	END SUBROUTINE dynamics_prognostic_equations_ij
881
882
883	!--------------------------------------------------------------------------------------------------!
884	! Description:
885	! ------------
886	!> Compute boundary conditions of dynamics model
887	!--------------------------------------------------------------------------------------------------!
888	SUBROUTINE dynamics_boundary_conditions
889
890	IMPLICIT NONE
891
892	INTEGER(iwp) :: i !< grid index x direction
893	INTEGER(iwp) :: j !< grid index y direction
894	INTEGER(iwp) :: k !< grid index z direction
895	INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls
896	INTEGER(iwp) :: m !< running index surface elements
897
898	REAL(wp) :: c_max !< maximum phase velocity allowed by CFL criterion, used for outflow boundary condition
899	REAL(wp) :: denom !< horizontal gradient of velocity component normal to the outflow boundary
900
901	!
902	!-- Bottom boundary
903	IF ( ibc_uv_b == 1 ) THEN
904	u_p(nzb,:,:) = u_p(nzb+1,:,:)
905	v_p(nzb,:,:) = v_p(nzb+1,:,:)
906	ENDIF
907	!
908	!-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case
909	!-- of downward-facing surfaces.
910	DO l = 0, 1
911	!$OMP PARALLEL DO PRIVATE( i, j, k )
912	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
913	!$ACC PRESENT(bc_h, w_p)
914	DO m = 1, bc_h(l)%ns
915	i = bc_h(l)%i(m)
916	j = bc_h(l)%j(m)
917	k = bc_h(l)%k(m)
918	w_p(k+bc_h(l)%koff,j,i) = 0.0_wp
919	ENDDO
920	ENDDO
921
922	!
923	!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
924	IF ( ibc_uv_t == 0 ) THEN
925	!$ACC KERNELS PRESENT(u_p, v_p, u_init, v_init)
926	u_p(nzt+1,:,:) = u_init(nzt+1)
927	v_p(nzt+1,:,:) = v_init(nzt+1)
928	!$ACC END KERNELS
929	ELSEIF ( ibc_uv_t == 1 ) THEN
930	u_p(nzt+1,:,:) = u_p(nzt,:,:)
931	v_p(nzt+1,:,:) = v_p(nzt,:,:)
932	ENDIF
933
934	!
935	!-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
936	IF ( .NOT. child_domain .AND. .NOT. nesting_offline .AND. &
937	TRIM(coupling_mode) /= 'vnested_fine' ) THEN
938	!$ACC KERNELS PRESENT(w_p)
939	w_p(nzt:nzt+1,:,:) = 0.0_wp !< nzt is not a prognostic level (but cf. pres)
940	!$ACC END KERNELS
941	ENDIF
942
943	!
944	!-- Temperature at bottom and top boundary.
945	!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by the sea surface temperature
946	!-- of the coupled ocean model.
947	!-- Dirichlet
948	IF ( .NOT. neutral ) THEN
949	IF ( ibc_pt_b == 0 ) THEN
950	DO l = 0, 1
951	!$OMP PARALLEL DO PRIVATE( i, j, k )
952	DO m = 1, bc_h(l)%ns
953	i = bc_h(l)%i(m)
954	j = bc_h(l)%j(m)
955	k = bc_h(l)%k(m)
956	pt_p(k+bc_h(l)%koff,j,i) = pt(k+bc_h(l)%koff,j,i)
957	ENDDO
958	ENDDO
959	!
960	!-- Neumann, zero-gradient
961	ELSEIF ( ibc_pt_b == 1 ) THEN
962	DO l = 0, 1
963	!$OMP PARALLEL DO PRIVATE( i, j, k )
964	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
965	!$ACC PRESENT(bc_h, pt_p)
966	DO m = 1, bc_h(l)%ns
967	i = bc_h(l)%i(m)
968	j = bc_h(l)%j(m)
969	k = bc_h(l)%k(m)
970	pt_p(k+bc_h(l)%koff,j,i) = pt_p(k,j,i)
971	ENDDO
972	ENDDO
973	ENDIF
974
975	!
976	!-- Temperature at top boundary
977	IF ( ibc_pt_t == 0 ) THEN
978	pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
979	!
980	!-- In case of nudging adjust top boundary to pt which is
981	!-- read in from NUDGING-DATA
982	IF ( nudging ) THEN
983	pt_p(nzt+1,:,:) = pt_init(nzt+1)
984	ENDIF
985	ELSEIF ( ibc_pt_t == 1 ) THEN
986	pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
987	ELSEIF ( ibc_pt_t == 2 ) THEN
988	!$ACC KERNELS PRESENT(pt_p, dzu)
989	pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
990	!$ACC END KERNELS
991	ENDIF
992	ENDIF
993	!
994	!-- Boundary conditions for total water content, bottom and top boundary (see also temperature)
995	IF ( humidity ) THEN
996	!
997	!-- Surface conditions for constant_humidity_flux
998	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype the k coordinate
999	!-- belongs to the atmospheric grid point, therefore, set q_p at k-1
1000	IF ( ibc_q_b == 0 ) THEN
1001
1002	DO l = 0, 1
1003	!$OMP PARALLEL DO PRIVATE( i, j, k )
1004	DO m = 1, bc_h(l)%ns
1005	i = bc_h(l)%i(m)
1006	j = bc_h(l)%j(m)
1007	k = bc_h(l)%k(m)
1008	q_p(k+bc_h(l)%koff,j,i) = q(k+bc_h(l)%koff,j,i)
1009	ENDDO
1010	ENDDO
1011
1012	ELSE
1013
1014	DO l = 0, 1
1015	!$OMP PARALLEL DO PRIVATE( i, j, k )
1016	DO m = 1, bc_h(l)%ns
1017	i = bc_h(l)%i(m)
1018	j = bc_h(l)%j(m)
1019	k = bc_h(l)%k(m)
1020	q_p(k+bc_h(l)%koff,j,i) = q_p(k,j,i)
1021	ENDDO
1022	ENDDO
1023	ENDIF
1024	!
1025	!-- Top boundary
1026	IF ( ibc_q_t == 0 ) THEN
1027	q_p(nzt+1,:,:) = q(nzt+1,:,:)
1028	ELSEIF ( ibc_q_t == 1 ) THEN
1029	q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
1030	ENDIF
1031	ENDIF
1032	!
1033	!-- Boundary conditions for scalar, bottom and top boundary (see also temperature)
1034	IF ( passive_scalar ) THEN
1035	!
1036	!-- Surface conditions for constant_humidity_flux
1037	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype the k coordinate
1038	!-- belongs to the atmospheric grid point, therefore, set s_p at k-1
1039	IF ( ibc_s_b == 0 ) THEN
1040
1041	DO l = 0, 1
1042	!$OMP PARALLEL DO PRIVATE( i, j, k )
1043	DO m = 1, bc_h(l)%ns
1044	i = bc_h(l)%i(m)
1045	j = bc_h(l)%j(m)
1046	k = bc_h(l)%k(m)
1047	s_p(k+bc_h(l)%koff,j,i) = s(k+bc_h(l)%koff,j,i)
1048	ENDDO
1049	ENDDO
1050
1051	ELSE
1052
1053	DO l = 0, 1
1054	!$OMP PARALLEL DO PRIVATE( i, j, k )
1055	DO m = 1, bc_h(l)%ns
1056	i = bc_h(l)%i(m)
1057	j = bc_h(l)%j(m)
1058	k = bc_h(l)%k(m)
1059	s_p(k+bc_h(l)%koff,j,i) = s_p(k,j,i)
1060	ENDDO
1061	ENDDO
1062	ENDIF
1063	!
1064	!-- Top boundary condition
1065	IF ( ibc_s_t == 0 ) THEN
1066	s_p(nzt+1,:,:) = s(nzt+1,:,:)
1067	ELSEIF ( ibc_s_t == 1 ) THEN
1068	s_p(nzt+1,:,:) = s_p(nzt,:,:)
1069	ELSEIF ( ibc_s_t == 2 ) THEN
1070	s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
1071	ENDIF
1072
1073	ENDIF
1074	!
1075	!-- In case of inflow or nest boundary at the south boundary the boundary for v is at nys and in
1076	!-- case of inflow or nest boundary at the left boundary the boundary for u is at nxl. Since in
1077	!-- prognostic_equations (cache optimized version) these levels are handled as a prognostic level,
1078	!-- boundary values have to be restored here.
1079	IF ( bc_dirichlet_s ) THEN
1080	v_p(:,nys,:) = v_p(:,nys-1,:)
1081	ELSEIF ( bc_dirichlet_l ) THEN
1082	u_p(:,:,nxl) = u_p(:,:,nxl-1)
1083	ENDIF
1084
1085	!
1086	!-- The same restoration for u at i=nxl and v at j=nys as above must be made in case of nest
1087	!-- boundaries. This must not be done in case of vertical nesting mode as in that case the lateral
1088	!-- boundaries are actually cyclic.
1089	!-- Lateral oundary conditions for TKE and dissipation are set in tcm_boundary_conds.
1090	IF ( nesting_mode /= 'vertical' .OR. nesting_offline ) THEN
1091	IF ( bc_dirichlet_s ) THEN
1092	v_p(:,nys,:) = v_p(:,nys-1,:)
1093	ENDIF
1094	IF ( bc_dirichlet_l ) THEN
1095	u_p(:,:,nxl) = u_p(:,:,nxl-1)
1096	ENDIF
1097	ENDIF
1098
1099	!
1100	!-- Lateral boundary conditions for scalar quantities at the outflow.
1101	!-- Lateral oundary conditions for TKE and dissipation are set in tcm_boundary_conds.
1102	IF ( bc_radiation_s ) THEN
1103	pt_p(:,nys-1,:) = pt_p(:,nys,:)
1104	IF ( humidity ) THEN
1105	q_p(:,nys-1,:) = q_p(:,nys,:)
1106	ENDIF
1107	IF ( passive_scalar ) s_p(:,nys-1,:) = s_p(:,nys,:)
1108	ELSEIF ( bc_radiation_n ) THEN
1109	pt_p(:,nyn+1,:) = pt_p(:,nyn,:)
1110	IF ( humidity ) THEN
1111	q_p(:,nyn+1,:) = q_p(:,nyn,:)
1112	ENDIF
1113	IF ( passive_scalar ) s_p(:,nyn+1,:) = s_p(:,nyn,:)
1114	ELSEIF ( bc_radiation_l ) THEN
1115	pt_p(:,:,nxl-1) = pt_p(:,:,nxl)
1116	IF ( humidity ) THEN
1117	q_p(:,:,nxl-1) = q_p(:,:,nxl)
1118	ENDIF
1119	IF ( passive_scalar ) s_p(:,:,nxl-1) = s_p(:,:,nxl)
1120	ELSEIF ( bc_radiation_r ) THEN
1121	pt_p(:,:,nxr+1) = pt_p(:,:,nxr)
1122	IF ( humidity ) THEN
1123	q_p(:,:,nxr+1) = q_p(:,:,nxr)
1124	ENDIF
1125	IF ( passive_scalar ) s_p(:,:,nxr+1) = s_p(:,:,nxr)
1126	ENDIF
1127
1128	!
1129	!-- Radiation boundary conditions for the velocities at the respective outflow.
1130	!-- The phase velocity is either assumed to the maximum phase velocity that ensures numerical
1131	!-- stability (CFL-condition) or calculated after Orlanski(1976) and averaged along the outflow
1132	!-- boundary.
1133	IF ( bc_radiation_s ) THEN
1134
1135	IF ( use_cmax ) THEN
1136	u_p(:,-1,:) = u(:,0,:)
1137	v_p(:,0,:) = v(:,1,:)
1138	w_p(:,-1,:) = w(:,0,:)
1139	ELSEIF ( .NOT. use_cmax ) THEN
1140
1141	c_max = dy / dt_3d
1142
1143	c_u_m_l = 0.0_wp
1144	c_v_m_l = 0.0_wp
1145	c_w_m_l = 0.0_wp
1146
1147	c_u_m = 0.0_wp
1148	c_v_m = 0.0_wp
1149	c_w_m = 0.0_wp
1150
1151	!
1152	!-- Calculate the phase speeds for u, v, and w, first local and then average along the outflow
1153	!-- boundary.
1154	DO k = nzb+1, nzt+1
1155	DO i = nxl, nxr
1156
1157	denom = u_m_s(k,0,i) - u_m_s(k,1,i)
1158
1159	IF ( denom /= 0.0_wp ) THEN
1160	c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) )
1161	IF ( c_u(k,i) < 0.0_wp ) THEN
1162	c_u(k,i) = 0.0_wp
1163	ELSEIF ( c_u(k,i) > c_max ) THEN
1164	c_u(k,i) = c_max
1165	ENDIF
1166	ELSE
1167	c_u(k,i) = c_max
1168	ENDIF
1169
1170	denom = v_m_s(k,1,i) - v_m_s(k,2,i)
1171
1172	IF ( denom /= 0.0_wp ) THEN
1173	c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) )
1174	IF ( c_v(k,i) < 0.0_wp ) THEN
1175	c_v(k,i) = 0.0_wp
1176	ELSEIF ( c_v(k,i) > c_max ) THEN
1177	c_v(k,i) = c_max
1178	ENDIF
1179	ELSE
1180	c_v(k,i) = c_max
1181	ENDIF
1182
1183	denom = w_m_s(k,0,i) - w_m_s(k,1,i)
1184
1185	IF ( denom /= 0.0_wp ) THEN
1186	c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
1187	IF ( c_w(k,i) < 0.0_wp ) THEN
1188	c_w(k,i) = 0.0_wp
1189	ELSEIF ( c_w(k,i) > c_max ) THEN
1190	c_w(k,i) = c_max
1191	ENDIF
1192	ELSE
1193	c_w(k,i) = c_max
1194	ENDIF
1195
1196	c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
1197	c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
1198	c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
1199
1200	ENDDO
1201	ENDDO
1202
1203	#if defined( __parallel )
1204	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
1205	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, &
1206	ierr )
1207	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
1208	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, &
1209	ierr )
1210	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
1211	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, &
1212	ierr )
1213	#else
1214	c_u_m = c_u_m_l
1215	c_v_m = c_v_m_l
1216	c_w_m = c_w_m_l
1217	#endif
1218
1219	c_u_m = c_u_m / (nx+1)
1220	c_v_m = c_v_m / (nx+1)
1221	c_w_m = c_w_m / (nx+1)
1222
1223	!
1224	!-- Save old timelevels for the next timestep
1225	IF ( intermediate_timestep_count == 1 ) THEN
1226	u_m_s(:,:,:) = u(:,0:1,:)
1227	v_m_s(:,:,:) = v(:,1:2,:)
1228	w_m_s(:,:,:) = w(:,0:1,:)
1229	ENDIF
1230
1231	!
1232	!-- Calculate the new velocities
1233	DO k = nzb+1, nzt+1
1234	DO i = nxlg, nxrg
1235	u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) * &
1236	( u(k,-1,i) - u(k,0,i) ) * ddy
1237
1238	v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v_m(k) * &
1239	( v(k,0,i) - v(k,1,i) ) * ddy
1240
1241	w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) * &
1242	( w(k,-1,i) - w(k,0,i) ) * ddy
1243	ENDDO
1244	ENDDO
1245
1246	!
1247	!-- Bottom boundary at the outflow
1248	IF ( ibc_uv_b == 0 ) THEN
1249	u_p(nzb,-1,:) = 0.0_wp
1250	v_p(nzb,0,:) = 0.0_wp
1251	ELSE
1252	u_p(nzb,-1,:) = u_p(nzb+1,-1,:)
1253	v_p(nzb,0,:) = v_p(nzb+1,0,:)
1254	ENDIF
1255	w_p(nzb,-1,:) = 0.0_wp
1256
1257	!
1258	!-- Top boundary at the outflow
1259	IF ( ibc_uv_t == 0 ) THEN
1260	u_p(nzt+1,-1,:) = u_init(nzt+1)
1261	v_p(nzt+1,0,:) = v_init(nzt+1)
1262	ELSE
1263	u_p(nzt+1,-1,:) = u_p(nzt,-1,:)
1264	v_p(nzt+1,0,:) = v_p(nzt,0,:)
1265	ENDIF
1266	w_p(nzt:nzt+1,-1,:) = 0.0_wp
1267
1268	ENDIF
1269
1270	ENDIF
1271
1272	IF ( bc_radiation_n ) THEN
1273
1274	IF ( use_cmax ) THEN
1275	u_p(:,ny+1,:) = u(:,ny,:)
1276	v_p(:,ny+1,:) = v(:,ny,:)
1277	w_p(:,ny+1,:) = w(:,ny,:)
1278	ELSEIF ( .NOT. use_cmax ) THEN
1279
1280	c_max = dy / dt_3d
1281
1282	c_u_m_l = 0.0_wp
1283	c_v_m_l = 0.0_wp
1284	c_w_m_l = 0.0_wp
1285
1286	c_u_m = 0.0_wp
1287	c_v_m = 0.0_wp
1288	c_w_m = 0.0_wp
1289
1290	!
1291	!-- Calculate the phase speeds for u, v, and w, first local and then average along the outflow
1292	!-- boundary.
1293	DO k = nzb+1, nzt+1
1294	DO i = nxl, nxr
1295
1296	denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
1297
1298	IF ( denom /= 0.0_wp ) THEN
1299	c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
1300	IF ( c_u(k,i) < 0.0_wp ) THEN
1301	c_u(k,i) = 0.0_wp
1302	ELSEIF ( c_u(k,i) > c_max ) THEN
1303	c_u(k,i) = c_max
1304	ENDIF
1305	ELSE
1306	c_u(k,i) = c_max
1307	ENDIF
1308
1309	denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
1310
1311	IF ( denom /= 0.0_wp ) THEN
1312	c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
1313	IF ( c_v(k,i) < 0.0_wp ) THEN
1314	c_v(k,i) = 0.0_wp
1315	ELSEIF ( c_v(k,i) > c_max ) THEN
1316	c_v(k,i) = c_max
1317	ENDIF
1318	ELSE
1319	c_v(k,i) = c_max
1320	ENDIF
1321
1322	denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
1323
1324	IF ( denom /= 0.0_wp ) THEN
1325	c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
1326	IF ( c_w(k,i) < 0.0_wp ) THEN
1327	c_w(k,i) = 0.0_wp
1328	ELSEIF ( c_w(k,i) > c_max ) THEN
1329	c_w(k,i) = c_max
1330	ENDIF
1331	ELSE
1332	c_w(k,i) = c_max
1333	ENDIF
1334
1335	c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
1336	c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
1337	c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
1338
1339	ENDDO
1340	ENDDO
1341
1342	#if defined( __parallel )
1343	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
1344	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, &
1345	ierr )
1346	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
1347	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, &
1348	ierr )
1349	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
1350	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx, &
1351	ierr )
1352	#else
1353	c_u_m = c_u_m_l
1354	c_v_m = c_v_m_l
1355	c_w_m = c_w_m_l
1356	#endif
1357
1358	c_u_m = c_u_m / (nx+1)
1359	c_v_m = c_v_m / (nx+1)
1360	c_w_m = c_w_m / (nx+1)
1361
1362	!
1363	!-- Save old timelevels for the next timestep
1364	IF ( intermediate_timestep_count == 1 ) THEN
1365	u_m_n(:,:,:) = u(:,ny-1:ny,:)
1366	v_m_n(:,:,:) = v(:,ny-1:ny,:)
1367	w_m_n(:,:,:) = w(:,ny-1:ny,:)
1368	ENDIF
1369
1370	!
1371	!-- Calculate the new velocities
1372	DO k = nzb+1, nzt+1
1373	DO i = nxlg, nxrg
1374	u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) * &
1375	( u(k,ny+1,i) - u(k,ny,i) ) * ddy
1376
1377	v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v_m(k) * &
1378	( v(k,ny+1,i) - v(k,ny,i) ) * ddy
1379
1380	w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) * &
1381	( w(k,ny+1,i) - w(k,ny,i) ) * ddy
1382	ENDDO
1383	ENDDO
1384
1385	!
1386	!-- Bottom boundary at the outflow
1387	IF ( ibc_uv_b == 0 ) THEN
1388	u_p(nzb,ny+1,:) = 0.0_wp
1389	v_p(nzb,ny+1,:) = 0.0_wp
1390	ELSE
1391	u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:)
1392	v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:)
1393	ENDIF
1394	w_p(nzb,ny+1,:) = 0.0_wp
1395
1396	!
1397	!-- Top boundary at the outflow
1398	IF ( ibc_uv_t == 0 ) THEN
1399	u_p(nzt+1,ny+1,:) = u_init(nzt+1)
1400	v_p(nzt+1,ny+1,:) = v_init(nzt+1)
1401	ELSE
1402	u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
1403	v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
1404	ENDIF
1405	w_p(nzt:nzt+1,ny+1,:) = 0.0_wp
1406
1407	ENDIF
1408
1409	ENDIF
1410
1411	IF ( bc_radiation_l ) THEN
1412
1413	IF ( use_cmax ) THEN
1414	u_p(:,:,0) = u(:,:,1)
1415	v_p(:,:,-1) = v(:,:,0)
1416	w_p(:,:,-1) = w(:,:,0)
1417	ELSEIF ( .NOT. use_cmax ) THEN
1418
1419	c_max = dx / dt_3d
1420
1421	c_u_m_l = 0.0_wp
1422	c_v_m_l = 0.0_wp
1423	c_w_m_l = 0.0_wp
1424
1425	c_u_m = 0.0_wp
1426	c_v_m = 0.0_wp
1427	c_w_m = 0.0_wp
1428
1429	!
1430	!-- Calculate the phase speeds for u, v, and w, first local and then average along the outflow
1431	!-- boundary.
1432	DO k = nzb+1, nzt+1
1433	DO j = nys, nyn
1434
1435	denom = u_m_l(k,j,1) - u_m_l(k,j,2)
1436
1437	IF ( denom /= 0.0_wp ) THEN
1438	c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
1439	IF ( c_u(k,j) < 0.0_wp ) THEN
1440	c_u(k,j) = 0.0_wp
1441	ELSEIF ( c_u(k,j) > c_max ) THEN
1442	c_u(k,j) = c_max
1443	ENDIF
1444	ELSE
1445	c_u(k,j) = c_max
1446	ENDIF
1447
1448	denom = v_m_l(k,j,0) - v_m_l(k,j,1)
1449
1450	IF ( denom /= 0.0_wp ) THEN
1451	c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
1452	IF ( c_v(k,j) < 0.0_wp ) THEN
1453	c_v(k,j) = 0.0_wp
1454	ELSEIF ( c_v(k,j) > c_max ) THEN
1455	c_v(k,j) = c_max
1456	ENDIF
1457	ELSE
1458	c_v(k,j) = c_max
1459	ENDIF
1460
1461	denom = w_m_l(k,j,0) - w_m_l(k,j,1)
1462
1463	IF ( denom /= 0.0_wp ) THEN
1464	c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
1465	IF ( c_w(k,j) < 0.0_wp ) THEN
1466	c_w(k,j) = 0.0_wp
1467	ELSEIF ( c_w(k,j) > c_max ) THEN
1468	c_w(k,j) = c_max
1469	ENDIF
1470	ELSE
1471	c_w(k,j) = c_max
1472	ENDIF
1473
1474	c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
1475	c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
1476	c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
1477
1478	ENDDO
1479	ENDDO
1480
1481	#if defined( __parallel )
1482	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
1483	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, &
1484	ierr )
1485	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
1486	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, &
1487	ierr )
1488	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
1489	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, &
1490	ierr )
1491	#else
1492	c_u_m = c_u_m_l
1493	c_v_m = c_v_m_l
1494	c_w_m = c_w_m_l
1495	#endif
1496
1497	c_u_m = c_u_m / (ny+1)
1498	c_v_m = c_v_m / (ny+1)
1499	c_w_m = c_w_m / (ny+1)
1500
1501	!
1502	!-- Save old timelevels for the next timestep
1503	IF ( intermediate_timestep_count == 1 ) THEN
1504	u_m_l(:,:,:) = u(:,:,1:2)
1505	v_m_l(:,:,:) = v(:,:,0:1)
1506	w_m_l(:,:,:) = w(:,:,0:1)
1507	ENDIF
1508
1509	!
1510	!-- Calculate the new velocities
1511	DO k = nzb+1, nzt+1
1512	DO j = nysg, nyng
1513	u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) * &
1514	( u(k,j,0) - u(k,j,1) ) * ddx
1515
1516	v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) * &
1517	( v(k,j,-1) - v(k,j,0) ) * ddx
1518
1519	w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) * &
1520	( w(k,j,-1) - w(k,j,0) ) * ddx
1521	ENDDO
1522	ENDDO
1523
1524	!
1525	!-- Bottom boundary at the outflow
1526	IF ( ibc_uv_b == 0 ) THEN
1527	u_p(nzb,:,0) = 0.0_wp
1528	v_p(nzb,:,-1) = 0.0_wp
1529	ELSE
1530	u_p(nzb,:,0) = u_p(nzb+1,:,0)
1531	v_p(nzb,:,-1) = v_p(nzb+1,:,-1)
1532	ENDIF
1533	w_p(nzb,:,-1) = 0.0_wp
1534
1535	!
1536	!-- Top boundary at the outflow
1537	IF ( ibc_uv_t == 0 ) THEN
1538	u_p(nzt+1,:,0) = u_init(nzt+1)
1539	v_p(nzt+1,:,-1) = v_init(nzt+1)
1540	ELSE
1541	u_p(nzt+1,:,0) = u_p(nzt,:,0)
1542	v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
1543	ENDIF
1544	w_p(nzt:nzt+1,:,-1) = 0.0_wp
1545
1546	ENDIF
1547
1548	ENDIF
1549
1550	IF ( bc_radiation_r ) THEN
1551
1552	IF ( use_cmax ) THEN
1553	u_p(:,:,nx+1) = u(:,:,nx)
1554	v_p(:,:,nx+1) = v(:,:,nx)
1555	w_p(:,:,nx+1) = w(:,:,nx)
1556	ELSEIF ( .NOT. use_cmax ) THEN
1557
1558	c_max = dx / dt_3d
1559
1560	c_u_m_l = 0.0_wp
1561	c_v_m_l = 0.0_wp
1562	c_w_m_l = 0.0_wp
1563
1564	c_u_m = 0.0_wp
1565	c_v_m = 0.0_wp
1566	c_w_m = 0.0_wp
1567
1568	!
1569	!-- Calculate the phase speeds for u, v, and w, first local and then average along the outflow
1570	!-- boundary.
1571	DO k = nzb+1, nzt+1
1572	DO j = nys, nyn
1573
1574	denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
1575
1576	IF ( denom /= 0.0_wp ) THEN
1577	c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
1578	IF ( c_u(k,j) < 0.0_wp ) THEN
1579	c_u(k,j) = 0.0_wp
1580	ELSEIF ( c_u(k,j) > c_max ) THEN
1581	c_u(k,j) = c_max
1582	ENDIF
1583	ELSE
1584	c_u(k,j) = c_max
1585	ENDIF
1586
1587	denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
1588
1589	IF ( denom /= 0.0_wp ) THEN
1590	c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
1591	IF ( c_v(k,j) < 0.0_wp ) THEN
1592	c_v(k,j) = 0.0_wp
1593	ELSEIF ( c_v(k,j) > c_max ) THEN
1594	c_v(k,j) = c_max
1595	ENDIF
1596	ELSE
1597	c_v(k,j) = c_max
1598	ENDIF
1599
1600	denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
1601
1602	IF ( denom /= 0.0_wp ) THEN
1603	c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
1604	IF ( c_w(k,j) < 0.0_wp ) THEN
1605	c_w(k,j) = 0.0_wp
1606	ELSEIF ( c_w(k,j) > c_max ) THEN
1607	c_w(k,j) = c_max
1608	ENDIF
1609	ELSE
1610	c_w(k,j) = c_max
1611	ENDIF
1612
1613	c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
1614	c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
1615	c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
1616
1617	ENDDO
1618	ENDDO
1619
1620	#if defined( __parallel )
1621	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
1622	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, &
1623	ierr )
1624	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
1625	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, &
1626	ierr )
1627	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
1628	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy, &
1629	ierr )
1630	#else
1631	c_u_m = c_u_m_l
1632	c_v_m = c_v_m_l
1633	c_w_m = c_w_m_l
1634	#endif
1635
1636	c_u_m = c_u_m / (ny+1)
1637	c_v_m = c_v_m / (ny+1)
1638	c_w_m = c_w_m / (ny+1)
1639
1640	!
1641	!-- Save old timelevels for the next timestep
1642	IF ( intermediate_timestep_count == 1 ) THEN
1643	u_m_r(:,:,:) = u(:,:,nx-1:nx)
1644	v_m_r(:,:,:) = v(:,:,nx-1:nx)
1645	w_m_r(:,:,:) = w(:,:,nx-1:nx)
1646	ENDIF
1647
1648	!
1649	!-- Calculate the new velocities
1650	DO k = nzb+1, nzt+1
1651	DO j = nysg, nyng
1652	u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) * &
1653	( u(k,j,nx+1) - u(k,j,nx) ) * ddx
1654
1655	v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) * &
1656	( v(k,j,nx+1) - v(k,j,nx) ) * ddx
1657
1658	w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) * &
1659	( w(k,j,nx+1) - w(k,j,nx) ) * ddx
1660	ENDDO
1661	ENDDO
1662
1663	!
1664	!-- Bottom boundary at the outflow
1665	IF ( ibc_uv_b == 0 ) THEN
1666	u_p(nzb,:,nx+1) = 0.0_wp
1667	v_p(nzb,:,nx+1) = 0.0_wp
1668	ELSE
1669	u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1)
1670	v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1)
1671	ENDIF
1672	w_p(nzb,:,nx+1) = 0.0_wp
1673
1674	!
1675	!-- Top boundary at the outflow
1676	IF ( ibc_uv_t == 0 ) THEN
1677	u_p(nzt+1,:,nx+1) = u_init(nzt+1)
1678	v_p(nzt+1,:,nx+1) = v_init(nzt+1)
1679	ELSE
1680	u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
1681	v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
1682	ENDIF
1683	w_p(nzt:nzt+1,:,nx+1) = 0.0_wp
1684
1685	ENDIF
1686
1687	ENDIF
1688
1689	END SUBROUTINE dynamics_boundary_conditions
1690	!--------------------------------------------------------------------------------------------------!
1691	! Description:
1692	! ------------
1693	!> Swap timelevels of module-specific array pointers
1694	!--------------------------------------------------------------------------------------------------!
1695	SUBROUTINE dynamics_swap_timelevel ( mod_count )
1696
1697
1698	INTEGER, INTENT(IN) :: mod_count
1699
1700
1701	SELECT CASE ( mod_count )
1702
1703	CASE ( 0 )
1704
1705	u => u_1; u_p => u_2
1706	v => v_1; v_p => v_2
1707	w => w_1; w_p => w_2
1708	IF ( .NOT. neutral ) THEN
1709	pt => pt_1; pt_p => pt_2
1710	ENDIF
1711	IF ( humidity ) THEN
1712	q => q_1; q_p => q_2
1713	ENDIF
1714	IF ( passive_scalar ) THEN
1715	s => s_1; s_p => s_2
1716	ENDIF
1717
1718	CASE ( 1 )
1719
1720	u => u_2; u_p => u_1
1721	v => v_2; v_p => v_1
1722	w => w_2; w_p => w_1
1723	IF ( .NOT. neutral ) THEN
1724	pt => pt_2; pt_p => pt_1
1725	ENDIF
1726	IF ( humidity ) THEN
1727	q => q_2; q_p => q_1
1728	ENDIF
1729	IF ( passive_scalar ) THEN
1730	s => s_2; s_p => s_1
1731	ENDIF
1732
1733	END SELECT
1734
1735	END SUBROUTINE dynamics_swap_timelevel
1736
1737
1738	!--------------------------------------------------------------------------------------------------!
1739	! Description:
1740	! ------------
1741	!> Sum up and time-average module-specific output quantities as well as allocate the array necessary
1742	!> for storing the average.
1743	!--------------------------------------------------------------------------------------------------!
1744	SUBROUTINE dynamics_3d_data_averaging( mode, variable )
1745
1746
1747	CHARACTER (LEN=*) :: mode !<
1748	CHARACTER (LEN=*) :: variable !<
1749
1750
1751	IF ( mode == 'allocate' ) THEN
1752
1753	SELECT CASE ( TRIM( variable ) )
1754
1755	! CASE ( 'u2' )
1756
1757	CASE DEFAULT
1758	CONTINUE
1759
1760	END SELECT
1761
1762	ELSEIF ( mode == 'sum' ) THEN
1763
1764	SELECT CASE ( TRIM( variable ) )
1765
1766	! CASE ( 'u2' )
1767
1768	CASE DEFAULT
1769	CONTINUE
1770
1771	END SELECT
1772
1773	ELSEIF ( mode == 'average' ) THEN
1774
1775	SELECT CASE ( TRIM( variable ) )
1776
1777	! CASE ( 'u2' )
1778
1779	END SELECT
1780
1781	ENDIF
1782
1783
1784	END SUBROUTINE dynamics_3d_data_averaging
1785
1786
1787	!--------------------------------------------------------------------------------------------------!
1788	! Description:
1789	! ------------
1790	!> Resorts the module-specific output quantity with indices (k,j,i) to a temporary array with
1791	!> indices (i,j,k) and sets the grid on which it is defined.
1792	!> Allowed values for grid are "zu" and "zw".
1793	!--------------------------------------------------------------------------------------------------!
1794	SUBROUTINE dynamics_data_output_2d( av, variable, found, grid, mode, local_pf, two_d, nzb_do, &
1795	nzt_do, fill_value )
1796
1797
1798	CHARACTER (LEN=*) :: grid !<
1799	CHARACTER (LEN=*), INTENT(IN) :: mode !< either 'xy', 'xz' or 'yz'
1800	CHARACTER (LEN=*) :: variable !<
1801
1802	INTEGER(iwp) :: av !< flag to control data output of instantaneous or time-averaged data
1803	! INTEGER(iwp) :: i !< grid index along x-direction
1804	! INTEGER(iwp) :: j !< grid index along y-direction
1805	! INTEGER(iwp) :: k !< grid index along z-direction
1806	! INTEGER(iwp) :: m !< running index surface elements
1807	INTEGER(iwp) :: nzb_do !< lower limit of the domain (usually nzb)
1808	INTEGER(iwp) :: nzt_do !< upper limit of the domain (usually nzt+1)
1809
1810	LOGICAL :: found !<
1811	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
1812
1813	REAL(wp), INTENT(IN) :: fill_value
1814
1815	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
1816
1817	!
1818	!-- Next line is just to avoid compiler warnings about unused variables. You may remove it.
1819	IF ( two_d .AND. av + LEN( mode ) + local_pf(nxl,nys,nzb_do) + fill_value < 0.0 ) CONTINUE
1820
1821	found = .TRUE.
1822
1823	SELECT CASE ( TRIM( variable ) )
1824
1825	! CASE ( 'u2_xy', 'u2_xz', 'u2_yz' )
1826
1827	CASE DEFAULT
1828	found = .FALSE.
1829	grid = 'none'
1830
1831	END SELECT
1832
1833
1834	END SUBROUTINE dynamics_data_output_2d
1835
1836
1837	!--------------------------------------------------------------------------------------------------!
1838	! Description:
1839	! ------------
1840	!> Resorts the module-specific output quantity with indices (k,j,i) to a temporary array with
1841	!> indices (i,j,k).
1842	!--------------------------------------------------------------------------------------------------!
1843	SUBROUTINE dynamics_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
1844
1845
1846	CHARACTER (LEN=*) :: variable !<
1847
1848	INTEGER(iwp) :: av !<
1849	! INTEGER(iwp) :: i !<
1850	! INTEGER(iwp) :: j !<
1851	! INTEGER(iwp) :: k !<
1852	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
1853	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
1854
1855	LOGICAL :: found !<
1856
1857	REAL(wp), INTENT(IN) :: fill_value !< value for the _FillValue attribute
1858
1859	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
1860
1861	!
1862	!-- Next line is to avoid compiler warning about unused variables. Please remove.
1863	IF ( av + local_pf(nxl,nys,nzb_do) + fill_value < 0.0 ) CONTINUE
1864
1865
1866	found = .TRUE.
1867
1868	SELECT CASE ( TRIM( variable ) )
1869
1870	! CASE ( 'u2' )
1871
1872	CASE DEFAULT
1873	found = .FALSE.
1874
1875	END SELECT
1876
1877
1878	END SUBROUTINE dynamics_data_output_3d
1879
1880
1881	!--------------------------------------------------------------------------------------------------!
1882	! Description:
1883	! ------------
1884	!> Calculation of module-specific statistics, i.e. horizontally averaged profiles and time series.
1885	!> This is called for every statistic region sr, but at least for the region "total domain" (sr=0).
1886	!--------------------------------------------------------------------------------------------------!
1887	SUBROUTINE dynamics_statistics( mode, sr, tn )
1888
1889
1890	CHARACTER (LEN=*) :: mode !<
1891	! INTEGER(iwp) :: i !<
1892	! INTEGER(iwp) :: j !<
1893	! INTEGER(iwp) :: k !<
1894	INTEGER(iwp) :: sr !<
1895	INTEGER(iwp) :: tn !<
1896
1897	!
1898	!-- Next line is to avoid compiler warning about unused variables. Please remove.
1899	IF ( sr == 0 .OR. tn == 0 ) CONTINUE
1900
1901	IF ( mode == 'profiles' ) THEN
1902
1903	ELSEIF ( mode == 'time_series' ) THEN
1904
1905	ENDIF
1906
1907	END SUBROUTINE dynamics_statistics
1908
1909
1910	!--------------------------------------------------------------------------------------------------!
1911	! Description:
1912	! ------------
1913	!> Read module-specific global restart data (Fortran binary format).
1914	!--------------------------------------------------------------------------------------------------!
1915	SUBROUTINE dynamics_rrd_global_ftn( found )
1916
1917	LOGICAL, INTENT(OUT) :: found
1918
1919
1920	found = .TRUE.
1921
1922
1923	SELECT CASE ( restart_string(1:length) )
1924
1925	CASE ( 'global_paramter' )
1926	! READ ( 13 ) global_parameter
1927
1928	CASE DEFAULT
1929
1930	found = .FALSE.
1931
1932	END SELECT
1933
1934
1935	END SUBROUTINE dynamics_rrd_global_ftn
1936
1937
1938	!--------------------------------------------------------------------------------------------------!
1939	! Description:
1940	! ------------
1941	!> Read module-specific global restart data (MPI-IO).
1942	!--------------------------------------------------------------------------------------------------!
1943	SUBROUTINE dynamics_rrd_global_mpi
1944
1945	! CALL rrd_mpi_io( 'global_parameter', global_parameter )
1946	CONTINUE
1947
1948	END SUBROUTINE dynamics_rrd_global_mpi
1949
1950
1951	!--------------------------------------------------------------------------------------------------!
1952	! Description:
1953	! ------------
1954	!> Read module-specific local restart data arrays (Fortran binary format).
1955	!> Subdomain index limits on file are given by nxl_on_file, etc.
1956	!> Indices nxlc, etc. indicate the range of gridpoints to be mapped from the subdomain on file (f)
1957	!> to the subdomain of the current PE (c). They have been calculated in routine rrd_local.
1958	!--------------------------------------------------------------------------------------------------!
1959	SUBROUTINE dynamics_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, &
1960	nync, nyn_on_file, nysf, nysc, nys_on_file, tmp_2d, tmp_3d, &
1961	found )
1962
1963
1964	INTEGER(iwp) :: k !<
1965	INTEGER(iwp) :: nxlc !<
1966	INTEGER(iwp) :: nxlf !<
1967	INTEGER(iwp) :: nxl_on_file !<
1968	INTEGER(iwp) :: nxrc !<
1969	INTEGER(iwp) :: nxrf !<
1970	INTEGER(iwp) :: nxr_on_file !<
1971	INTEGER(iwp) :: nync !<
1972	INTEGER(iwp) :: nynf !<
1973	INTEGER(iwp) :: nyn_on_file !<
1974	INTEGER(iwp) :: nysc !<
1975	INTEGER(iwp) :: nysf !<
1976	INTEGER(iwp) :: nys_on_file !<
1977
1978	LOGICAL, INTENT(OUT) :: found
1979
1980	REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d !<
1981	REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !<
1982
1983	!
1984	!-- Next line is to avoid compiler warning about unused variables. Please remove.
1985	IF ( k + nxlc + nxlf + nxrc + nxrf + nync + nynf + nysc + nysf + &
1986	tmp_2d(nys_on_file,nxl_on_file) + &
1987	tmp_3d(nzb,nys_on_file,nxl_on_file) < 0.0 ) CONTINUE
1988	!
1989	!-- Here the reading of user-defined restart data follows:
1990	!-- Sample for user-defined output
1991
1992	found = .TRUE.
1993
1994	SELECT CASE ( restart_string(1:length) )
1995
1996	! CASE ( 'u2_av' )
1997
1998	CASE DEFAULT
1999
2000	found = .FALSE.
2001
2002	END SELECT
2003
2004	END SUBROUTINE dynamics_rrd_local_ftn
2005
2006
2007	!--------------------------------------------------------------------------------------------------!
2008	! Description:
2009	! ------------
2010	!> Read module-specific local restart data arrays (MPI-IO).
2011	!--------------------------------------------------------------------------------------------------!
2012	SUBROUTINE dynamics_rrd_local_mpi
2013
2014	IMPLICIT NONE
2015
2016	! LOGICAL :: array_found !<
2017
2018
2019	! CALL rd_mpi_io_check_array( 'u2_av' , found = array_found )
2020	! IF ( array_found ) THEN
2021	! IF ( .NOT. ALLOCATED( u2_av ) ) ALLOCATE( u2_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2022	! CALL rrd_mpi_io( 'u2_av', u2_av )
2023	! ENDIF
2024
2025	CONTINUE
2026
2027	END SUBROUTINE dynamics_rrd_local_mpi
2028
2029
2030
2031	!--------------------------------------------------------------------------------------------------!
2032	! Description:
2033	! ------------
2034	!> Writes global module-specific restart data into binary file(s) for restart runs.
2035	!--------------------------------------------------------------------------------------------------!
2036	SUBROUTINE dynamics_wrd_global
2037
2038
2039	END SUBROUTINE dynamics_wrd_global
2040
2041
2042	!--------------------------------------------------------------------------------------------------!
2043	! Description:
2044	! ------------
2045	!> Writes processor specific and module-specific restart data into binary file(s) for restart runs.
2046	!--------------------------------------------------------------------------------------------------!
2047	SUBROUTINE dynamics_wrd_local
2048
2049
2050	END SUBROUTINE dynamics_wrd_local
2051
2052
2053	!--------------------------------------------------------------------------------------------------!
2054	! Description:
2055	! ------------
2056	!> Execute module-specific actions at the very end of the program.
2057	!--------------------------------------------------------------------------------------------------!
2058	SUBROUTINE dynamics_last_actions
2059
2060
2061	END SUBROUTINE dynamics_last_actions
2062
2063	END MODULE dynamics_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |