Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |