Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |