Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |