Home

Context Navigation

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

Last change on this file since 4680 was 4627, checked in by raasch, 4 years ago
bugfix for r4646
Property svn:keywords set to `Id`
File size: 68.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |