Home

Context Navigation

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

Last change on this file since 4626 was 4626, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 67.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |