Home

Context Navigation

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

Last change on this file since 4810 was 4768, checked in by suehring, 4 years ago
Enable 3D data output also with 64-bit precision
Property svn:keywords set to `Id`
File size: 71.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |