Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |