Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |