Home

Context Navigation

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

Last change on this file since 4842 was 4842, checked in by raasch, 4 years ago
reading of namelist file and actions in case of namelist errors revised so that statement labels and goto statements are not required any more, deprecated namelists removed
Property svn:keywords set to `Id`
File size: 71.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |