Home

Context Navigation

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

Last change on this file since 4686 was 4627, checked in by raasch, 4 years ago
bugfix for r4646
Property svn:keywords set to `Id`
File size: 68.0 KB

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