Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |