Home

Context Navigation

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

Last change on this file since 4498 was 4495, checked in by raasch, 5 years ago
restart data handling with MPI-IO added, first part
Property svn:keywords set to `Id`
File size: 60.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |