Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |