Home

Context Navigation

source: palm/trunk/SOURCE/nesting_offl_mod.f90 @ 3987

Last change on this file since 3987 was 3987, checked in by kanani, 5 years ago
clean up location, debug and error messages
Property svn:keywords set to `Id`
File size: 58.2 KB

Rev	Line
[3347]	1	!> @file nesting_offl_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	!
[3655]	17	! Copyright 1997-2019 Leibniz Universitaet Hannover
[3347]	18	!------------------------------------------------------------------------------!
	19	!
	20	! Current revisions:
	21	! ------------------
	22	!
[3705]	23	!
[3347]	24	! Former revisions:
	25	! -----------------
[3413]	26	! $Id: nesting_offl_mod.f90 3987 2019-05-22 09:52:13Z kanani $
[3987]	27	! Introduce alternative switch for debug output during timestepping
	28	!
	29	! 3964 2019-05-09 09:48:32Z suehring
[3964]	30	! Ensure that veloctiy term in calculation of bulk Richardson number does not
	31	! become zero
	32	!
	33	! 3937 2019-04-29 15:09:07Z suehring
[3937]	34	! Set boundary conditon on upper-left and upper-south grid point for the u- and
	35	! v-component, respectively.
	36	!
	37	! 3891 2019-04-12 17:52:01Z suehring
[3891]	38	! Bugfix, do not overwrite lateral and top boundary data in case of restart
	39	! runs.
	40	!
	41	! 3885 2019-04-11 11:29:34Z kanani
[3885]	42	! Changes related to global restructuring of location messages and introduction
	43	! of additional debug messages
	44	!
	45	!
[3858]	46	! Do local data exchange for chemistry variables only when boundary data is
	47	! coming from dynamic file
	48	!
	49	! 3737 2019-02-12 16:57:06Z suehring
[3737]	50	! Introduce mesoscale nesting for chemical species
	51	!
	52	! 3705 2019-01-29 19:56:39Z suehring
[3705]	53	! Formatting adjustments
	54	!
	55	! 3704 2019-01-29 19:51:41Z suehring
[3579]	56	! Check implemented for offline nesting in child domain
	57	!
	58	! 3413 2018-10-24 10:28:44Z suehring
[3413]	59	! Keyword ID set
	60	!
	61	! 3404 2018-10-23 13:29:11Z suehring
[3404]	62	! Consider time-dependent geostrophic wind components in offline nesting
	63	!
	64	!
[3347]	65	! Initial Revision:
	66	! - separate offline nesting from large_scale_nudging_mod
	67	! - revise offline nesting, adjust for usage of synthetic turbulence generator
	68	! - adjust Rayleigh damping depending on the time-depending atmospheric
	69	! conditions
	70	!
	71	!
	72	! Description:
	73	! ------------
	74	!> Offline nesting in larger-scale models. Boundary conditions for the simulation
	75	!> are read from NetCDF file and are prescribed onto the respective arrays.
	76	!> Further, a mass-flux correction is performed to maintain the mass balance.
	77	!--------------------------------------------------------------------------------!
	78	MODULE nesting_offl_mod
	79
	80	USE arrays_3d, &
	81	ONLY: dzw, e, diss, pt, pt_init, q, q_init, s, u, u_init, ug, v, &
[3737]	82	v_init, vg, w, zu, zw
	83
[3876]	84	USE chem_modules, &
[3737]	85	ONLY: chem_species
[3347]	86
	87	USE control_parameters, &
[3737]	88	ONLY: air_chemistry, bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, &
[3885]	89	bc_dirichlet_s, dt_3d, dz, constant_diffusion, &
[3987]	90	debug_output_timestep, &
	91	humidity, &
	92	initializing_actions, &
[3737]	93	message_string, nesting_offline, neutral, passive_scalar, &
	94	rans_mode, rans_tke_e, time_since_reference_point, volume_flow
[3347]	95
	96	USE cpulog, &
	97	ONLY: cpu_log, log_point
	98
	99	USE grid_variables
	100
	101	USE indices, &
	102	ONLY: nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys, &
	103	nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
	104
	105	USE kinds
	106
	107	USE netcdf_data_input_mod, &
	108	ONLY: nest_offl
	109
	110	USE pegrid
	111
	112	REAL(wp) :: zi_ribulk = 0.0_wp !< boundary-layer depth according to bulk Richardson criterion, i.e. the height where Ri_bulk exceeds the critical
	113	!< bulk Richardson number of 0.25
	114
	115	SAVE
	116	PRIVATE
	117	!
	118	!-- Public subroutines
[3579]	119	PUBLIC nesting_offl_bc, nesting_offl_check_parameters, nesting_offl_header,&
	120	nesting_offl_init, nesting_offl_mass_conservation, nesting_offl_parin
[3347]	121	!
	122	!-- Public variables
	123	PUBLIC zi_ribulk
	124
	125	INTERFACE nesting_offl_bc
	126	MODULE PROCEDURE nesting_offl_bc
	127	END INTERFACE nesting_offl_bc
	128
[3579]	129	INTERFACE nesting_offl_check_parameters
	130	MODULE PROCEDURE nesting_offl_check_parameters
	131	END INTERFACE nesting_offl_check_parameters
	132
[3347]	133	INTERFACE nesting_offl_header
	134	MODULE PROCEDURE nesting_offl_header
	135	END INTERFACE nesting_offl_header
	136
	137	INTERFACE nesting_offl_init
	138	MODULE PROCEDURE nesting_offl_init
	139	END INTERFACE nesting_offl_init
	140
	141	INTERFACE nesting_offl_mass_conservation
	142	MODULE PROCEDURE nesting_offl_mass_conservation
	143	END INTERFACE nesting_offl_mass_conservation
	144
	145	INTERFACE nesting_offl_parin
	146	MODULE PROCEDURE nesting_offl_parin
	147	END INTERFACE nesting_offl_parin
	148
	149	CONTAINS
	150
	151
	152	!------------------------------------------------------------------------------!
	153	! Description:
	154	! ------------
	155	!> In this subroutine a constant mass within the model domain is guaranteed.
	156	!> Larger-scale models may be based on a compressible equation system, which is
	157	!> not consistent with PALMs incompressible equation system. In order to avoid
	158	!> a decrease or increase of mass during the simulation, non-divergent flow
	159	!> through the lateral and top boundaries is compensated by the vertical wind
	160	!> component at the top boundary.
	161	!------------------------------------------------------------------------------!
	162	SUBROUTINE nesting_offl_mass_conservation
	163
	164	IMPLICIT NONE
	165
	166	INTEGER(iwp) :: i !< grid index in x-direction
	167	INTEGER(iwp) :: j !< grid index in y-direction
	168	INTEGER(iwp) :: k !< grid index in z-direction
	169
	170	REAL(wp) :: d_area_t !< inverse of the total area of the horizontal model domain
	171	REAL(wp) :: w_correct !< vertical velocity increment required to compensate non-divergent flow through the boundaries
	172	REAL(wp), DIMENSION(1:3) :: volume_flow_l !< local volume flow
	173
[3987]	174
	175	IF ( debug_output_timestep ) CALL debug_message( 'nesting_offl_mass_conservation', 'start' )
	176
[3347]	177	CALL cpu_log( log_point(58), 'offline nesting', 'start' )
	178
	179	volume_flow = 0.0_wp
	180	volume_flow_l = 0.0_wp
	181
	182	d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
	183
	184	IF ( bc_dirichlet_l ) THEN
	185	i = nxl
	186	DO j = nys, nyn
	187	DO k = nzb+1, nzt
	188	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy &
	189	* MERGE( 1.0_wp, 0.0_wp, &
	190	BTEST( wall_flags_0(k,j,i), 1 ) )
	191	ENDDO
	192	ENDDO
	193	ENDIF
	194	IF ( bc_dirichlet_r ) THEN
	195	i = nxr+1
	196	DO j = nys, nyn
	197	DO k = nzb+1, nzt
	198	volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy &
	199	* MERGE( 1.0_wp, 0.0_wp, &
	200	BTEST( wall_flags_0(k,j,i), 1 ) )
	201	ENDDO
	202	ENDDO
	203	ENDIF
	204	IF ( bc_dirichlet_s ) THEN
	205	j = nys
	206	DO i = nxl, nxr
	207	DO k = nzb+1, nzt
	208	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx &
	209	* MERGE( 1.0_wp, 0.0_wp, &
	210	BTEST( wall_flags_0(k,j,i), 2 ) )
	211	ENDDO
	212	ENDDO
	213	ENDIF
	214	IF ( bc_dirichlet_n ) THEN
	215	j = nyn+1
	216	DO i = nxl, nxr
	217	DO k = nzb+1, nzt
	218	volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx &
	219	* MERGE( 1.0_wp, 0.0_wp, &
	220	BTEST( wall_flags_0(k,j,i), 2 ) )
	221	ENDDO
	222	ENDDO
	223	ENDIF
	224	!
	225	!-- Top boundary
	226	k = nzt
	227	DO i = nxl, nxr
	228	DO j = nys, nyn
	229	volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
	230	ENDDO
	231	ENDDO
	232
	233	#if defined( __parallel )
	234	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	235	CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM, &
	236	comm2d, ierr )
	237	#else
	238	volume_flow = volume_flow_l
	239	#endif
	240
	241	w_correct = SUM( volume_flow ) * d_area_t
	242
	243	DO i = nxl, nxr
	244	DO j = nys, nyn
	245	DO k = nzt, nzt + 1
	246	w(k,j,i) = w(k,j,i) + w_correct
	247	ENDDO
	248	ENDDO
	249	ENDDO
	250
	251	CALL cpu_log( log_point(58), 'offline nesting', 'stop' )
	252
[3987]	253	IF ( debug_output_timestep ) CALL debug_message( 'nesting_offl_mass_conservation', 'end' )
	254
[3347]	255	END SUBROUTINE nesting_offl_mass_conservation
	256
	257
	258	!------------------------------------------------------------------------------!
	259	! Description:
	260	! ------------
	261	!> Set the lateral and top boundary conditions in case the PALM domain is
	262	!> nested offline in a mesoscale model. Further, average boundary data and
	263	!> determine mean profiles, further used for correct damping in the sponge
	264	!> layer.
	265	!------------------------------------------------------------------------------!
	266	SUBROUTINE nesting_offl_bc
	267
	268	IMPLICIT NONE
	269
	270	INTEGER(iwp) :: i !< running index x-direction
	271	INTEGER(iwp) :: j !< running index y-direction
	272	INTEGER(iwp) :: k !< running index z-direction
[3737]	273	INTEGER(iwp) :: n !< running index for chemical species
[3347]	274
	275	REAL(wp) :: fac_dt !< interpolation factor
	276
	277	REAL(wp), DIMENSION(nzb:nzt+1) :: pt_ref !< reference profile for potential temperature
	278	REAL(wp), DIMENSION(nzb:nzt+1) :: pt_ref_l !< reference profile for potential temperature on subdomain
	279	REAL(wp), DIMENSION(nzb:nzt+1) :: q_ref !< reference profile for mixing ratio on subdomain
	280	REAL(wp), DIMENSION(nzb:nzt+1) :: q_ref_l !< reference profile for mixing ratio on subdomain
	281	REAL(wp), DIMENSION(nzb:nzt+1) :: u_ref !< reference profile for u-component on subdomain
	282	REAL(wp), DIMENSION(nzb:nzt+1) :: u_ref_l !< reference profile for u-component on subdomain
	283	REAL(wp), DIMENSION(nzb:nzt+1) :: v_ref !< reference profile for v-component on subdomain
	284	REAL(wp), DIMENSION(nzb:nzt+1) :: v_ref_l !< reference profile for v-component on subdomain
	285
[3885]	286
[3987]	287	IF ( debug_output_timestep ) CALL debug_message( 'nesting_offl_bc', 'start' )
	288
[3347]	289	CALL cpu_log( log_point(58), 'offline nesting', 'start' )
	290	!
	291	!-- Set mean profiles, derived from boundary data, to zero
	292	pt_ref = 0.0_wp
	293	q_ref = 0.0_wp
	294	u_ref = 0.0_wp
	295	v_ref = 0.0_wp
	296
	297	pt_ref_l = 0.0_wp
	298	q_ref_l = 0.0_wp
	299	u_ref_l = 0.0_wp
	300	v_ref_l = 0.0_wp
	301	!
	302	!-- Determine interpolation factor and limit it to 1. This is because
	303	!-- t+dt can slightly exceed time(tind_p) before boundary data is updated
	304	!-- again.
	305	fac_dt = ( time_since_reference_point - nest_offl%time(nest_offl%tind) &
	306	+ dt_3d ) / &
	307	( nest_offl%time(nest_offl%tind_p) - nest_offl%time(nest_offl%tind) )
	308	fac_dt = MIN( 1.0_wp, fac_dt )
	309	!
	310	!-- Set boundary conditions of u-, v-, w-component, as well as q, and pt.
	311	!-- Note, boundary values at the left boundary: i=-1 (v,w,pt,q) and
	312	!-- i=0 (u), at the right boundary: i=nxr+1 (all), at the south boundary:
	313	!-- j=-1 (u,w,pt,q) and j=0 (v), at the north boundary: j=nyn+1 (all).
	314	!-- Please note, at the left (for u) and south (for v) boundary, values
	315	!-- for u and v are set also at i/j=-1, since these values are used in
	316	!-- boundary_conditions() to restore prognostic values.
	317	!-- Further, sum up data to calculate mean profiles from boundary data,
	318	!-- used for Rayleigh damping.
	319	IF ( bc_dirichlet_l ) THEN
	320
	321	DO j = nys, nyn
	322	DO k = nzb+1, nzt
	323	u(k,j,0) = interpolate_in_time( nest_offl%u_left(0,k,j), &
	324	nest_offl%u_left(1,k,j), &
	325	fac_dt ) * &
	326	MERGE( 1.0_wp, 0.0_wp, &
	327	BTEST( wall_flags_0(k,j,0), 1 ) )
	328	u(k,j,-1) = u(k,j,0)
	329	ENDDO
	330	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,0)
	331	ENDDO
	332
	333	DO j = nys, nyn
	334	DO k = nzb+1, nzt-1
	335	w(k,j,-1) = interpolate_in_time( nest_offl%w_left(0,k,j), &
	336	nest_offl%w_left(1,k,j), &
	337	fac_dt ) * &
	338	MERGE( 1.0_wp, 0.0_wp, &
	339	BTEST( wall_flags_0(k,j,-1), 3 ) )
	340	ENDDO
	341	ENDDO
	342
	343	DO j = nysv, nyn
	344	DO k = nzb+1, nzt
	345	v(k,j,-1) = interpolate_in_time( nest_offl%v_left(0,k,j), &
	346	nest_offl%v_left(1,k,j), &
	347	fac_dt ) * &
	348	MERGE( 1.0_wp, 0.0_wp, &
	349	BTEST( wall_flags_0(k,j,-1), 2 ) )
	350	ENDDO
	351	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,-1)
	352	ENDDO
	353
	354	IF ( .NOT. neutral ) THEN
	355	DO j = nys, nyn
	356	DO k = nzb+1, nzt
	357	pt(k,j,-1) = interpolate_in_time( nest_offl%pt_left(0,k,j), &
	358	nest_offl%pt_left(1,k,j), &
	359	fac_dt )
	360
	361	ENDDO
	362	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,-1)
	363	ENDDO
	364	ENDIF
	365
	366	IF ( humidity ) THEN
	367	DO j = nys, nyn
	368	DO k = nzb+1, nzt
	369	q(k,j,-1) = interpolate_in_time( nest_offl%q_left(0,k,j), &
	370	nest_offl%q_left(1,k,j), &
	371	fac_dt )
	372
	373	ENDDO
	374	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,-1)
	375	ENDDO
	376	ENDIF
[3737]	377
	378	IF ( air_chemistry ) THEN
	379	DO n = 1, UBOUND( chem_species, 1 )
	380	IF ( nest_offl%chem_from_file_l(n) ) THEN
	381	DO j = nys, nyn
	382	DO k = nzb+1, nzt
	383	chem_species(n)%conc(k,j,-1) = interpolate_in_time( &
	384	nest_offl%chem_left(0,k,j,n),&
	385	nest_offl%chem_left(1,k,j,n),&
	386	fac_dt )
	387	ENDDO
	388	ENDDO
	389	ENDIF
	390	ENDDO
	391	ENDIF
[3347]	392
	393	ENDIF
	394
	395	IF ( bc_dirichlet_r ) THEN
	396
	397	DO j = nys, nyn
	398	DO k = nzb+1, nzt
	399	u(k,j,nxr+1) = interpolate_in_time( nest_offl%u_right(0,k,j), &
	400	nest_offl%u_right(1,k,j), &
	401	fac_dt ) * &
	402	MERGE( 1.0_wp, 0.0_wp, &
	403	BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
	404	ENDDO
	405	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,nxr+1)
	406	ENDDO
	407	DO j = nys, nyn
	408	DO k = nzb+1, nzt-1
	409	w(k,j,nxr+1) = interpolate_in_time( nest_offl%w_right(0,k,j), &
	410	nest_offl%w_right(1,k,j), &
	411	fac_dt ) * &
	412	MERGE( 1.0_wp, 0.0_wp, &
	413	BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
	414	ENDDO
	415	ENDDO
	416
	417	DO j = nysv, nyn
	418	DO k = nzb+1, nzt
	419	v(k,j,nxr+1) = interpolate_in_time( nest_offl%v_right(0,k,j), &
	420	nest_offl%v_right(1,k,j), &
	421	fac_dt ) * &
	422	MERGE( 1.0_wp, 0.0_wp, &
	423	BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
	424	ENDDO
	425	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,nxr+1)
	426	ENDDO
	427
	428	IF ( .NOT. neutral ) THEN
	429	DO j = nys, nyn
	430	DO k = nzb+1, nzt
	431	pt(k,j,nxr+1) = interpolate_in_time( &
	432	nest_offl%pt_right(0,k,j), &
	433	nest_offl%pt_right(1,k,j), &
	434	fac_dt )
	435	ENDDO
	436	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,nxr+1)
	437	ENDDO
	438	ENDIF
	439
	440	IF ( humidity ) THEN
	441	DO j = nys, nyn
	442	DO k = nzb+1, nzt
	443	q(k,j,nxr+1) = interpolate_in_time( &
	444	nest_offl%q_right(0,k,j), &
	445	nest_offl%q_right(1,k,j), &
	446	fac_dt )
	447
	448	ENDDO
	449	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,nxr+1)
	450	ENDDO
	451	ENDIF
[3737]	452
	453	IF ( air_chemistry ) THEN
	454	DO n = 1, UBOUND( chem_species, 1 )
	455	IF ( nest_offl%chem_from_file_r(n) ) THEN
	456	DO j = nys, nyn
	457	DO k = nzb+1, nzt
	458	chem_species(n)%conc(k,j,nxr+1) = interpolate_in_time(&
	459	nest_offl%chem_right(0,k,j,n),&
	460	nest_offl%chem_right(1,k,j,n),&
	461	fac_dt )
	462	ENDDO
	463	ENDDO
	464	ENDIF
	465	ENDDO
	466	ENDIF
[3347]	467
	468	ENDIF
	469
	470	IF ( bc_dirichlet_s ) THEN
	471
	472	DO i = nxl, nxr
	473	DO k = nzb+1, nzt
	474	v(k,0,i) = interpolate_in_time( nest_offl%v_south(0,k,i), &
	475	nest_offl%v_south(1,k,i), &
	476	fac_dt ) * &
	477	MERGE( 1.0_wp, 0.0_wp, &
	478	BTEST( wall_flags_0(k,0,i), 2 ) )
	479	v(k,-1,i) = v(k,0,i)
	480	ENDDO
	481	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,0,i)
	482	ENDDO
	483
	484	DO i = nxl, nxr
	485	DO k = nzb+1, nzt-1
	486	w(k,-1,i) = interpolate_in_time( nest_offl%w_south(0,k,i), &
	487	nest_offl%w_south(1,k,i), &
	488	fac_dt ) * &
	489	MERGE( 1.0_wp, 0.0_wp, &
	490	BTEST( wall_flags_0(k,-1,i), 3 ) )
	491	ENDDO
	492	ENDDO
	493
	494	DO i = nxlu, nxr
	495	DO k = nzb+1, nzt
	496	u(k,-1,i) = interpolate_in_time( nest_offl%u_south(0,k,i), &
	497	nest_offl%u_south(1,k,i), &
	498	fac_dt ) * &
	499	MERGE( 1.0_wp, 0.0_wp, &
	500	BTEST( wall_flags_0(k,-1,i), 1 ) )
	501	ENDDO
	502	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,-1,i)
	503	ENDDO
	504
	505	IF ( .NOT. neutral ) THEN
	506	DO i = nxl, nxr
	507	DO k = nzb+1, nzt
	508	pt(k,-1,i) = interpolate_in_time( &
	509	nest_offl%pt_south(0,k,i), &
	510	nest_offl%pt_south(1,k,i), &
	511	fac_dt )
	512
	513	ENDDO
	514	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,-1,i)
	515	ENDDO
	516	ENDIF
	517
	518	IF ( humidity ) THEN
	519	DO i = nxl, nxr
	520	DO k = nzb+1, nzt
	521	q(k,-1,i) = interpolate_in_time( &
	522	nest_offl%q_south(0,k,i), &
	523	nest_offl%q_south(1,k,i), &
	524	fac_dt )
	525
	526	ENDDO
	527	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,-1,i)
	528	ENDDO
	529	ENDIF
[3737]	530
	531	IF ( air_chemistry ) THEN
	532	DO n = 1, UBOUND( chem_species, 1 )
	533	IF ( nest_offl%chem_from_file_s(n) ) THEN
	534	DO i = nxl, nxr
	535	DO k = nzb+1, nzt
	536	chem_species(n)%conc(k,-1,i) = interpolate_in_time( &
	537	nest_offl%chem_south(0,k,i,n),&
	538	nest_offl%chem_south(1,k,i,n),&
	539	fac_dt )
	540	ENDDO
	541	ENDDO
	542	ENDIF
	543	ENDDO
	544	ENDIF
[3347]	545
	546	ENDIF
	547
	548	IF ( bc_dirichlet_n ) THEN
	549
	550	DO i = nxl, nxr
	551	DO k = nzb+1, nzt
	552	v(k,nyn+1,i) = interpolate_in_time( nest_offl%v_north(0,k,i), &
	553	nest_offl%v_north(1,k,i), &
	554	fac_dt ) * &
	555	MERGE( 1.0_wp, 0.0_wp, &
	556	BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
	557	ENDDO
	558	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,nyn+1,i)
	559	ENDDO
	560	DO i = nxl, nxr
	561	DO k = nzb+1, nzt-1
	562	w(k,nyn+1,i) = interpolate_in_time( nest_offl%w_north(0,k,i), &
	563	nest_offl%w_north(1,k,i), &
	564	fac_dt ) * &
	565	MERGE( 1.0_wp, 0.0_wp, &
	566	BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
	567	ENDDO
	568	ENDDO
	569
	570	DO i = nxlu, nxr
	571	DO k = nzb+1, nzt
	572	u(k,nyn+1,i) = interpolate_in_time( nest_offl%u_north(0,k,i), &
	573	nest_offl%u_north(1,k,i), &
	574	fac_dt ) * &
	575	MERGE( 1.0_wp, 0.0_wp, &
	576	BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
	577
	578	ENDDO
	579	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,nyn+1,i)
	580	ENDDO
	581
	582	IF ( .NOT. neutral ) THEN
	583	DO i = nxl, nxr
	584	DO k = nzb+1, nzt
	585	pt(k,nyn+1,i) = interpolate_in_time( &
	586	nest_offl%pt_north(0,k,i), &
	587	nest_offl%pt_north(1,k,i), &
	588	fac_dt )
	589
	590	ENDDO
	591	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,nyn+1,i)
	592	ENDDO
	593	ENDIF
	594
	595	IF ( humidity ) THEN
	596	DO i = nxl, nxr
	597	DO k = nzb+1, nzt
	598	q(k,nyn+1,i) = interpolate_in_time( &
	599	nest_offl%q_north(0,k,i), &
	600	nest_offl%q_north(1,k,i), &
	601	fac_dt )
	602
	603	ENDDO
	604	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,nyn+1,i)
	605	ENDDO
	606	ENDIF
[3737]	607
	608	IF ( air_chemistry ) THEN
	609	DO n = 1, UBOUND( chem_species, 1 )
	610	IF ( nest_offl%chem_from_file_n(n) ) THEN
	611	DO i = nxl, nxr
	612	DO k = nzb+1, nzt
	613	chem_species(n)%conc(k,nyn+1,i) = interpolate_in_time(&
	614	nest_offl%chem_north(0,k,i,n),&
	615	nest_offl%chem_north(1,k,i,n),&
	616	fac_dt )
	617	ENDDO
	618	ENDDO
	619	ENDIF
	620	ENDDO
	621	ENDIF
[3347]	622
	623	ENDIF
	624	!
	625	!-- Top boundary
	626	DO i = nxlu, nxr
	627	DO j = nys, nyn
	628	u(nzt+1,j,i) = interpolate_in_time( nest_offl%u_top(0,j,i), &
	629	nest_offl%u_top(1,j,i), &
	630	fac_dt ) * &
	631	MERGE( 1.0_wp, 0.0_wp, &
	632	BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
	633	u_ref_l(nzt+1) = u_ref_l(nzt+1) + u(nzt+1,j,i)
	634	ENDDO
	635	ENDDO
[3937]	636	!
	637	!-- For left boundary set boundary condition for u-component also at top
	638	!-- grid point.
	639	!-- Note, this has no effect on the numeric solution, only for data output.
	640	IF ( bc_dirichlet_l ) u(nzt+1,:,nxl) = u(nzt+1,:,nxlu)
[3347]	641
	642	DO i = nxl, nxr
	643	DO j = nysv, nyn
	644	v(nzt+1,j,i) = interpolate_in_time( nest_offl%v_top(0,j,i), &
	645	nest_offl%v_top(1,j,i), &
	646	fac_dt ) * &
	647	MERGE( 1.0_wp, 0.0_wp, &
	648	BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
	649	v_ref_l(nzt+1) = v_ref_l(nzt+1) + v(nzt+1,j,i)
	650	ENDDO
	651	ENDDO
[3937]	652	!
	653	!-- For south boundary set boundary condition for v-component also at top
	654	!-- grid point.
	655	!-- Note, this has no effect on the numeric solution, only for data output.
	656	IF ( bc_dirichlet_s ) v(nzt+1,nys,:) = v(nzt+1,nysv,:)
[3347]	657
	658	DO i = nxl, nxr
	659	DO j = nys, nyn
	660	w(nzt,j,i) = interpolate_in_time( nest_offl%w_top(0,j,i), &
	661	nest_offl%w_top(1,j,i), &
	662	fac_dt ) * &
	663	MERGE( 1.0_wp, 0.0_wp, &
	664	BTEST( wall_flags_0(nzt,j,i), 3 ) )
	665	w(nzt+1,j,i) = w(nzt,j,i)
	666	ENDDO
	667	ENDDO
	668
	669
	670	IF ( .NOT. neutral ) THEN
	671	DO i = nxl, nxr
	672	DO j = nys, nyn
	673	pt(nzt+1,j,i) = interpolate_in_time( nest_offl%pt_top(0,j,i), &
	674	nest_offl%pt_top(1,j,i), &
	675	fac_dt )
	676	pt_ref_l(nzt+1) = pt_ref_l(nzt+1) + pt(nzt+1,j,i)
	677	ENDDO
	678	ENDDO
	679	ENDIF
	680
	681	IF ( humidity ) THEN
	682	DO i = nxl, nxr
	683	DO j = nys, nyn
	684	q(nzt+1,j,i) = interpolate_in_time( nest_offl%q_top(0,j,i), &
	685	nest_offl%q_top(1,j,i), &
	686	fac_dt )
	687	q_ref_l(nzt+1) = q_ref_l(nzt+1) + q(nzt+1,j,i)
	688	ENDDO
	689	ENDDO
	690	ENDIF
[3737]	691
	692	IF ( air_chemistry ) THEN
	693	DO n = 1, UBOUND( chem_species, 1 )
	694	IF ( nest_offl%chem_from_file_t(n) ) THEN
	695	DO i = nxl, nxr
	696	DO j = nys, nyn
	697	chem_species(n)%conc(nzt+1,j,i) = interpolate_in_time( &
	698	nest_offl%chem_north(0,j,i,n), &
	699	nest_offl%chem_north(1,j,i,n), &
	700	fac_dt )
	701	ENDDO
	702	ENDDO
	703	ENDIF
	704	ENDDO
	705	ENDIF
[3347]	706	!
	707	!-- Moreover, set Neumann boundary condition for subgrid-scale TKE,
	708	!-- passive scalar, dissipation, and chemical species if required
	709	IF ( rans_mode .AND. rans_tke_e ) THEN
	710	IF ( bc_dirichlet_l ) diss(:,:,nxl-1) = diss(:,:,nxl)
	711	IF ( bc_dirichlet_r ) diss(:,:,nxr+1) = diss(:,:,nxr)
	712	IF ( bc_dirichlet_s ) diss(:,nys-1,:) = diss(:,nys,:)
	713	IF ( bc_dirichlet_n ) diss(:,nyn+1,:) = diss(:,nyn,:)
	714	ENDIF
	715	IF ( .NOT. constant_diffusion ) THEN
	716	IF ( bc_dirichlet_l ) e(:,:,nxl-1) = e(:,:,nxl)
	717	IF ( bc_dirichlet_r ) e(:,:,nxr+1) = e(:,:,nxr)
	718	IF ( bc_dirichlet_s ) e(:,nys-1,:) = e(:,nys,:)
	719	IF ( bc_dirichlet_n ) e(:,nyn+1,:) = e(:,nyn,:)
	720	e(nzt+1,:,:) = e(nzt,:,:)
	721	ENDIF
	722	IF ( passive_scalar ) THEN
	723	IF ( bc_dirichlet_l ) s(:,:,nxl-1) = s(:,:,nxl)
	724	IF ( bc_dirichlet_r ) s(:,:,nxr+1) = s(:,:,nxr)
	725	IF ( bc_dirichlet_s ) s(:,nys-1,:) = s(:,nys,:)
	726	IF ( bc_dirichlet_n ) s(:,nyn+1,:) = s(:,nyn,:)
	727	ENDIF
	728
	729	CALL exchange_horiz( u, nbgp )
	730	CALL exchange_horiz( v, nbgp )
	731	CALL exchange_horiz( w, nbgp )
	732	IF ( .NOT. neutral ) CALL exchange_horiz( pt, nbgp )
	733	IF ( humidity ) CALL exchange_horiz( q, nbgp )
[3737]	734	IF ( air_chemistry ) THEN
	735	DO n = 1, UBOUND( chem_species, 1 )
[3858]	736	!
	737	!-- Do local exchange only when necessary, i.e. when data is coming
	738	!-- from dynamic file.
	739	IF ( nest_offl%chem_from_file_t(n) ) &
	740	CALL exchange_horiz( chem_species(n)%conc, nbgp )
[3737]	741	ENDDO
	742	ENDIF
	743
[3347]	744	!
	745	!-- In case of Rayleigh damping, where the profiles u_init, v_init
	746	!-- q_init and pt_init are still used, update these profiles from the
	747	!-- averaged boundary data.
	748	!-- But first, average these data.
	749	#if defined( __parallel )
	750	CALL MPI_ALLREDUCE( u_ref_l, u_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM, &
	751	comm2d, ierr )
	752	CALL MPI_ALLREDUCE( v_ref_l, v_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM, &
	753	comm2d, ierr )
	754	IF ( humidity ) THEN
	755	CALL MPI_ALLREDUCE( q_ref_l, q_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM, &
	756	comm2d, ierr )
	757	ENDIF
	758	IF ( .NOT. neutral ) THEN
	759	CALL MPI_ALLREDUCE( pt_ref_l, pt_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,&
	760	comm2d, ierr )
	761	ENDIF
	762	#else
[3704]	763	u_ref = u_ref_l
	764	v_ref = v_ref_l
	765	IF ( humidity ) q_ref = q_ref_l
	766	IF ( .NOT. neutral ) pt_ref = pt_ref_l
[3347]	767	#endif
	768	!
[3704]	769	!-- Average data. Note, reference profiles up to nzt are derived from lateral
	770	!-- boundaries, at the model top it is derived from the top boundary. Thus,
	771	!-- number of input data is different from nzb:nzt compared to nzt+1.
	772	!-- Derived from lateral boundaries.
	773	u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx ), &
	774	KIND = wp )
	775	v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + nx + 1 ), &
	776	KIND = wp )
	777	IF ( humidity ) &
	778	q_ref(nzb:nzt) = q_ref(nzb:nzt) / REAL( 2.0_wp * &
	779	( ny + 1 + nx + 1 ), &
	780	KIND = wp )
	781	IF ( .NOT. neutral ) &
	782	pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp * &
	783	( ny + 1 + nx + 1 ), &
	784	KIND = wp )
[3347]	785	!
[3704]	786	!-- Derived from top boundary.
	787	u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx ), KIND = wp )
	788	v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny ) * ( nx + 1 ), KIND = wp )
	789	IF ( humidity ) &
	790	q_ref(nzt+1) = q_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), &
	791	KIND = wp )
	792	IF ( .NOT. neutral ) &
	793	pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), &
	794	KIND = wp )
[3347]	795	!
[3704]	796	!-- Write onto init profiles, which are used for damping
	797	u_init = u_ref
	798	v_init = v_ref
	799	IF ( humidity ) q_init = q_ref
	800	IF ( .NOT. neutral ) pt_init = pt_ref
[3347]	801	!
[3704]	802	!-- Set bottom boundary condition
	803	IF ( humidity ) q_init(nzb) = q_init(nzb+1)
	804	IF ( .NOT. neutral ) pt_init(nzb) = pt_init(nzb+1)
[3347]	805	!
[3704]	806	!-- Further, adjust Rayleigh damping height in case of time-changing conditions.
	807	!-- Therefore, calculate boundary-layer depth first.
	808	CALL calc_zi
	809	CALL adjust_sponge_layer
[3347]	810
	811	!
[3704]	812	!-- Update geostrophic wind components from dynamic input file.
	813	DO k = nzb+1, nzt
	814	ug(k) = interpolate_in_time( nest_offl%ug(0,k), nest_offl%ug(1,k), &
	815	fac_dt )
	816	vg(k) = interpolate_in_time( nest_offl%vg(0,k), nest_offl%vg(1,k), &
	817	fac_dt )
	818	ENDDO
	819	ug(nzt+1) = ug(nzt)
	820	vg(nzt+1) = vg(nzt)
[3347]	821
[3704]	822	CALL cpu_log( log_point(58), 'offline nesting', 'stop' )
[3347]	823
[3987]	824	IF ( debug_output_timestep ) CALL debug_message( 'nesting_offl_bc', 'end' )
	825
	826
[3347]	827	END SUBROUTINE nesting_offl_bc
	828
	829	!------------------------------------------------------------------------------!
	830	! Description:
	831	!------------------------------------------------------------------------------!
	832	!> Calculates the boundary-layer depth from the boundary data, according to
	833	!> bulk-Richardson criterion.
	834	!------------------------------------------------------------------------------!
	835	SUBROUTINE calc_zi
	836
	837	USE basic_constants_and_equations_mod, &
	838	ONLY: g
	839
	840	USE kinds
	841
	842	USE surface_mod, &
	843	ONLY: get_topography_top_index, get_topography_top_index_ji
	844
	845	IMPLICIT NONE
	846
	847	INTEGER(iwp) :: i !< loop index in x-direction
	848	INTEGER(iwp) :: j !< loop index in y-direction
	849	INTEGER(iwp) :: k !< loop index in z-direction
	850	INTEGER(iwp) :: k_surface !< topography top index in z-direction
	851
	852	REAL(wp) :: ri_bulk !< bulk Richardson number
	853	REAL(wp) :: ri_bulk_crit = 0.25_wp !< critical bulk Richardson number
	854	REAL(wp) :: topo_max !< maximum topography level in model domain
	855	REAL(wp) :: topo_max_l !< maximum topography level in subdomain
	856	REAL(wp) :: u_comp !< u-component
	857	REAL(wp) :: v_comp !< v-component
	858	REAL(wp) :: vpt_surface !< near-surface virtual potential temperature
	859	REAL(wp) :: zi_l !< mean boundary-layer depth on subdomain
	860	REAL(wp) :: zi_local !< local boundary-layer depth
	861
	862	REAL(wp), DIMENSION(nzb:nzt+1) :: vpt_col !< vertical profile of virtual potential temperature at (j,i)-grid point
	863
	864
	865	!
	866	!-- Calculate mean boundary-layer height from boundary data.
	867	!-- Start with the left and right boundaries.
	868	zi_l = 0.0_wp
	869	IF ( bc_dirichlet_l .OR. bc_dirichlet_r ) THEN
	870	!
	871	!-- Determine index along x. Please note, index indicates boundary
	872	!-- grid point for scalars.
	873	i = MERGE( -1, nxr + 1, bc_dirichlet_l )
	874
	875	DO j = nys, nyn
	876	!
	877	!-- Determine topography top index at current (j,i) index
	878	k_surface = get_topography_top_index_ji( j, i, 's' )
	879	!
	880	!-- Pre-compute surface virtual temperature. Therefore, use 2nd
	881	!-- prognostic level according to Heinze et al. (2017).
	882	IF ( humidity ) THEN
	883	vpt_surface = pt(k_surface+2,j,i) * &
	884	( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
	885	vpt_col = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
	886	ELSE
	887	vpt_surface = pt(k_surface+2,j,i)
	888	vpt_col = pt(:,j,i)
	889	ENDIF
	890	!
	891	!-- Calculate local boundary layer height from bulk Richardson number,
	892	!-- i.e. the height where the bulk Richardson number exceeds its
	893	!-- critical value of 0.25 (according to Heinze et al., 2017).
	894	!-- Note, no interpolation of u- and v-component is made, as both
	895	!-- are mainly mean inflow profiles with very small spatial variation.
[3964]	896	!-- Add a safety factor in case the velocity term becomes zero. This
	897	!-- may happen if overhanging 3D structures are directly located at
	898	!-- the boundary, where velocity inside the building is zero
	899	!-- (k_surface is the index of the lowest upward-facing surface).
[3347]	900	zi_local = 0.0_wp
	901	DO k = k_surface+1, nzt
	902	u_comp = MERGE( u(k,j,i+1), u(k,j,i), bc_dirichlet_l )
	903	v_comp = v(k,j,i)
	904	ri_bulk = zu(k) * g / vpt_surface * &
	905	( vpt_col(k) - vpt_surface ) / &
[3964]	906	( u_comp * u_comp + v_comp * v_comp + 1E-5_wp )
[3347]	907
	908	IF ( zi_local == 0.0_wp .AND. ri_bulk > ri_bulk_crit ) &
	909	zi_local = zu(k)
	910	ENDDO
	911	!
	912	!-- Assure that the minimum local boundary-layer depth is at least at
	913	!-- the second vertical grid level.
	914	zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )
	915
	916	ENDDO
	917
	918	ENDIF
	919	!
	920	!-- Do the same at the north and south boundaries.
	921	IF ( bc_dirichlet_s .OR. bc_dirichlet_n ) THEN
	922
	923	j = MERGE( -1, nyn + 1, bc_dirichlet_s )
	924
	925	DO i = nxl, nxr
	926	k_surface = get_topography_top_index_ji( j, i, 's' )
	927
	928	IF ( humidity ) THEN
	929	vpt_surface = pt(k_surface+2,j,i) * &
	930	( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
	931	vpt_col = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
	932	ELSE
	933	vpt_surface = pt(k_surface+2,j,i)
	934	vpt_col = pt(:,j,i)
	935	ENDIF
	936
	937	zi_local = 0.0_wp
	938	DO k = k_surface+1, nzt
	939	u_comp = u(k,j,i)
	940	v_comp = MERGE( v(k,j+1,i), v(k,j,i), bc_dirichlet_s )
	941	ri_bulk = zu(k) * g / vpt_surface * &
	942	( vpt_col(k) - vpt_surface ) / &
[3964]	943	( u_comp * u_comp + v_comp * v_comp + 1E-5_wp )
[3347]	944
	945	IF ( zi_local == 0.0_wp .AND. ri_bulk > 0.25_wp ) &
	946	zi_local = zu(k)
	947	ENDDO
	948	zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )
	949
	950	ENDDO
	951
	952	ENDIF
	953
	954	#if defined( __parallel )
	955	CALL MPI_ALLREDUCE( zi_l, zi_ribulk, 1, MPI_REAL, MPI_SUM, &
	956	comm2d, ierr )
	957	#else
	958	zi_ribulk = zi_l
	959	#endif
	960	zi_ribulk = zi_ribulk / REAL( 2 * nx + 2 * ny, KIND = wp )
	961	!
	962	!-- Finally, check if boundary layer depth is not below the any topography.
	963	!-- zi_ribulk will be used to adjust rayleigh damping height, i.e. the
	964	!-- lower level of the sponge layer, as well as to adjust the synthetic
	965	!-- turbulence generator accordingly. If Rayleigh damping would be applied
	966	!-- near buildings, etc., this would spoil the simulation results.
	967	topo_max_l = zw(MAXVAL( get_topography_top_index( 's' )))
	968
	969	#if defined( __parallel )
	970	CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX, &
	971	comm2d, ierr )
	972	#else
	973	topo_max = topo_max_l
	974	#endif
	975	zi_ribulk = MAX( zi_ribulk, topo_max )
[3937]	976
[3347]	977	END SUBROUTINE calc_zi
	978
	979
	980	!------------------------------------------------------------------------------!
	981	! Description:
	982	!------------------------------------------------------------------------------!
	983	!> Adjust the height where the rayleigh damping starts, i.e. the lower level
	984	!> of the sponge layer.
	985	!------------------------------------------------------------------------------!
	986	SUBROUTINE adjust_sponge_layer
	987
	988	USE arrays_3d, &
	989	ONLY: rdf, rdf_sc, zu
	990
	991	USE basic_constants_and_equations_mod, &
	992	ONLY: pi
	993
	994	USE control_parameters, &
	995	ONLY: rayleigh_damping_height, rayleigh_damping_factor
	996
	997	USE kinds
	998
	999	IMPLICIT NONE
	1000
	1001	INTEGER(iwp) :: k !< loop index in z-direction
	1002
	1003	REAL(wp) :: rdh !< updated Rayleigh damping height
	1004
	1005
	1006	IF ( rayleigh_damping_height > 0.0_wp .AND. &
	1007	rayleigh_damping_factor > 0.0_wp ) THEN
	1008	!
	1009	!-- Update Rayleigh-damping height and re-calculate height-depending
	1010	!-- damping coefficients.
	1011	!-- Assure that rayleigh damping starts well above the boundary layer.
	1012	rdh = MIN( MAX( zi_ribulk * 1.3_wp, 10.0_wp * dz(1) ), &
	1013	0.8_wp * zu(nzt), rayleigh_damping_height )
	1014	!
	1015	!-- Update Rayleigh damping factor
	1016	DO k = nzb+1, nzt
	1017	IF ( zu(k) >= rdh ) THEN
	1018	rdf(k) = rayleigh_damping_factor * &
	1019	( SIN( pi * 0.5_wp * ( zu(k) - rdh ) &
	1020	/ ( zu(nzt) - rdh ) ) &
	1021	)**2
	1022	ENDIF
	1023	ENDDO
	1024	rdf_sc = rdf
	1025
	1026	ENDIF
	1027
	1028	END SUBROUTINE adjust_sponge_layer
	1029
	1030	!------------------------------------------------------------------------------!
	1031	! Description:
	1032	! ------------
	1033	!> Performs consistency checks
	1034	!------------------------------------------------------------------------------!
	1035	SUBROUTINE nesting_offl_check_parameters
[3579]	1036
	1037	USE control_parameters, &
	1038	ONLY: child_domain, message_string, nesting_offline
[3347]	1039
	1040	IMPLICIT NONE
	1041	!
	1042	!-- Perform checks
[3579]	1043	IF ( nesting_offline .AND. child_domain ) THEN
	1044	message_string = 'Offline nesting is only applicable in root model.'
	1045	CALL message( 'stg_check_parameters', 'PA0622', 1, 2, 0, 6, 0 )
	1046	ENDIF
[3347]	1047
	1048
	1049	END SUBROUTINE nesting_offl_check_parameters
	1050
	1051	!------------------------------------------------------------------------------!
	1052	! Description:
	1053	! ------------
	1054	!> Reads the parameter list nesting_offl_parameters
	1055	!------------------------------------------------------------------------------!
	1056	SUBROUTINE nesting_offl_parin
	1057
	1058	IMPLICIT NONE
	1059
	1060	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
	1061
	1062
	1063	NAMELIST /nesting_offl_parameters/ nesting_offline
	1064
	1065	line = ' '
	1066
	1067	!
	1068	!-- Try to find stg package
	1069	REWIND ( 11 )
	1070	line = ' '
	1071	DO WHILE ( INDEX( line, '&nesting_offl_parameters' ) == 0 )
	1072	READ ( 11, '(A)', END=20 ) line
	1073	ENDDO
	1074	BACKSPACE ( 11 )
	1075
	1076	!
	1077	!-- Read namelist
	1078	READ ( 11, nesting_offl_parameters, ERR = 10, END = 20 )
	1079
	1080	GOTO 20
	1081
	1082	10 BACKSPACE( 11 )
	1083	READ( 11 , '(A)') line
	1084	CALL parin_fail_message( 'nesting_offl_parameters', line )
	1085
	1086	20 CONTINUE
	1087
	1088
	1089	END SUBROUTINE nesting_offl_parin
	1090
	1091	!------------------------------------------------------------------------------!
	1092	! Description:
	1093	! ------------
	1094	!> Writes information about offline nesting into HEADER file
	1095	!------------------------------------------------------------------------------!
	1096	SUBROUTINE nesting_offl_header ( io )
	1097
	1098	IMPLICIT NONE
	1099
	1100	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
	1101
	1102	WRITE ( io, 1 )
	1103	IF ( nesting_offline ) THEN
	1104	WRITE ( io, 3 )
	1105	ELSE
	1106	WRITE ( io, 2 )
	1107	ENDIF
	1108
	1109	1 FORMAT (//' Offline nesting in COSMO model:'/ &
	1110	' -------------------------------'/)
	1111	2 FORMAT (' --> No offlince nesting is used (default) ')
	1112	3 FORMAT (' --> Offlince nesting is used. Boundary data is read from dynamic input file ')
	1113
	1114	END SUBROUTINE nesting_offl_header
	1115
	1116	!------------------------------------------------------------------------------!
	1117	! Description:
	1118	! ------------
	1119	!> Allocate arrays used to read boundary data from NetCDF file and initialize
	1120	!> boundary data.
	1121	!------------------------------------------------------------------------------!
	1122	SUBROUTINE nesting_offl_init
	1123
	1124	USE netcdf_data_input_mod, &
	1125	ONLY: netcdf_data_input_offline_nesting
	1126
	1127	IMPLICIT NONE
[3737]	1128
	1129	INTEGER(iwp) :: n !< running index for chemical species
[3347]	1130
	1131
	1132	!-- Allocate arrays for geostrophic wind components. Arrays will
[3404]	1133	!-- incorporate 2 time levels in order to interpolate in between.
	1134	ALLOCATE( nest_offl%ug(0:1,1:nzt) )
	1135	ALLOCATE( nest_offl%vg(0:1,1:nzt) )
[3347]	1136	!
	1137	!-- Allocate arrays for reading boundary values. Arrays will incorporate 2
	1138	!-- time levels in order to interpolate in between.
	1139	IF ( bc_dirichlet_l ) THEN
	1140	ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn) )
	1141	ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
	1142	ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
	1143	IF ( humidity ) ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn) )
	1144	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
[3737]	1145	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_left(0:1,nzb+1:nzt,nys:nyn,&
	1146	1:UBOUND( chem_species, 1 )) )
[3347]	1147	ENDIF
	1148	IF ( bc_dirichlet_r ) THEN
	1149	ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn) )
	1150	ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
	1151	ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
	1152	IF ( humidity ) ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn) )
	1153	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
[3737]	1154	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_right(0:1,nzb+1:nzt,nys:nyn,&
	1155	1:UBOUND( chem_species, 1 )) )
[3347]	1156	ENDIF
	1157	IF ( bc_dirichlet_n ) THEN
	1158	ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
	1159	ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr) )
	1160	ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
	1161	IF ( humidity ) ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr) )
	1162	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
[3737]	1163	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_north(0:1,nzb+1:nzt,nxl:nxr,&
	1164	1:UBOUND( chem_species, 1 )) )
[3347]	1165	ENDIF
	1166	IF ( bc_dirichlet_s ) THEN
	1167	ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
	1168	ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr) )
	1169	ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr) )
	1170	IF ( humidity ) ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr) )
	1171	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
[3737]	1172	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_south(0:1,nzb+1:nzt,nxl:nxr,&
	1173	1:UBOUND( chem_species, 1 )) )
[3347]	1174	ENDIF
	1175
	1176	ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
	1177	ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
	1178	ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr) )
	1179	IF ( humidity ) ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr) )
	1180	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
[3737]	1181	IF ( air_chemistry ) ALLOCATE( nest_offl%chem_top(0:1,nys:nyn,nxl:nxr, &
	1182	1:UBOUND( chem_species, 1 )) )
[3347]	1183	!
[3737]	1184	!-- For chemical species, create the names of the variables. This is necessary
	1185	!-- to identify the respective variable and write it onto the correct array
	1186	!-- in the chem_species datatype.
	1187	IF ( air_chemistry ) THEN
	1188	ALLOCATE( nest_offl%chem_from_file_l(1:UBOUND( chem_species, 1 )) )
	1189	ALLOCATE( nest_offl%chem_from_file_n(1:UBOUND( chem_species, 1 )) )
	1190	ALLOCATE( nest_offl%chem_from_file_r(1:UBOUND( chem_species, 1 )) )
	1191	ALLOCATE( nest_offl%chem_from_file_s(1:UBOUND( chem_species, 1 )) )
	1192	ALLOCATE( nest_offl%chem_from_file_t(1:UBOUND( chem_species, 1 )) )
	1193
	1194	ALLOCATE( nest_offl%var_names_chem_l(1:UBOUND( chem_species, 1 )) )
	1195	ALLOCATE( nest_offl%var_names_chem_n(1:UBOUND( chem_species, 1 )) )
	1196	ALLOCATE( nest_offl%var_names_chem_r(1:UBOUND( chem_species, 1 )) )
	1197	ALLOCATE( nest_offl%var_names_chem_s(1:UBOUND( chem_species, 1 )) )
	1198	ALLOCATE( nest_offl%var_names_chem_t(1:UBOUND( chem_species, 1 )) )
	1199	!
	1200	!-- Initialize flags that indicate whether the variable is on file or
	1201	!-- not. Please note, this is only necessary for chemistry variables.
	1202	nest_offl%chem_from_file_l(:) = .FALSE.
	1203	nest_offl%chem_from_file_n(:) = .FALSE.
	1204	nest_offl%chem_from_file_r(:) = .FALSE.
	1205	nest_offl%chem_from_file_s(:) = .FALSE.
	1206	nest_offl%chem_from_file_t(:) = .FALSE.
	1207
	1208	DO n = 1, UBOUND( chem_species, 1 )
	1209	nest_offl%var_names_chem_l(n) = nest_offl%char_l // &
	1210	TRIM(chem_species(n)%name)
	1211	nest_offl%var_names_chem_n(n) = nest_offl%char_n // &
	1212	TRIM(chem_species(n)%name)
	1213	nest_offl%var_names_chem_r(n) = nest_offl%char_r // &
	1214	TRIM(chem_species(n)%name)
	1215	nest_offl%var_names_chem_s(n) = nest_offl%char_s // &
	1216	TRIM(chem_species(n)%name)
	1217	nest_offl%var_names_chem_t(n) = nest_offl%char_t // &
	1218	TRIM(chem_species(n)%name)
	1219	ENDDO
	1220	ENDIF
	1221	!
[3347]	1222	!-- Read COSMO data at lateral and top boundaries
	1223	CALL netcdf_data_input_offline_nesting
	1224	!
[3891]	1225	!-- Initialize boundary data. Please note, do not initialize boundaries in
	1226	!-- case of restart runs. This case the boundaries are already initialized
	1227	!-- and the boundary data from file would be on the wrong time level.
	1228	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
	1229	IF ( bc_dirichlet_l ) THEN
	1230	u(nzb+1:nzt,nys:nyn,0) = nest_offl%u_left(0,nzb+1:nzt,nys:nyn)
	1231	v(nzb+1:nzt,nysv:nyn,-1) = nest_offl%v_left(0,nzb+1:nzt,nysv:nyn)
	1232	w(nzb+1:nzt-1,nys:nyn,-1) = nest_offl%w_left(0,nzb+1:nzt-1,nys:nyn)
	1233	IF ( .NOT. neutral ) pt(nzb+1:nzt,nys:nyn,-1) = &
	1234	nest_offl%pt_left(0,nzb+1:nzt,nys:nyn)
	1235	IF ( humidity ) q(nzb+1:nzt,nys:nyn,-1) = &
	1236	nest_offl%q_left(0,nzb+1:nzt,nys:nyn)
	1237	IF ( air_chemistry ) THEN
	1238	DO n = 1, UBOUND( chem_species, 1 )
	1239	IF( nest_offl%chem_from_file_l(n) ) THEN
	1240	chem_species(n)%conc(nzb+1:nzt,nys:nyn,-1) = &
	1241	nest_offl%chem_left(0,nzb+1:nzt,nys:nyn,n)
	1242	ENDIF
	1243	ENDDO
	1244	ENDIF
[3737]	1245	ENDIF
[3891]	1246	IF ( bc_dirichlet_r ) THEN
	1247	u(nzb+1:nzt,nys:nyn,nxr+1) = nest_offl%u_right(0,nzb+1:nzt,nys:nyn)
	1248	v(nzb+1:nzt,nysv:nyn,nxr+1) = nest_offl%v_right(0,nzb+1:nzt,nysv:nyn)
	1249	w(nzb+1:nzt-1,nys:nyn,nxr+1) = nest_offl%w_right(0,nzb+1:nzt-1,nys:nyn)
	1250	IF ( .NOT. neutral ) pt(nzb+1:nzt,nys:nyn,nxr+1) = &
	1251	nest_offl%pt_right(0,nzb+1:nzt,nys:nyn)
	1252	IF ( humidity ) q(nzb+1:nzt,nys:nyn,nxr+1) = &
	1253	nest_offl%q_right(0,nzb+1:nzt,nys:nyn)
	1254	IF ( air_chemistry ) THEN
	1255	DO n = 1, UBOUND( chem_species, 1 )
	1256	IF( nest_offl%chem_from_file_r(n) ) THEN
	1257	chem_species(n)%conc(nzb+1:nzt,nys:nyn,nxr+1) = &
	1258	nest_offl%chem_right(0,nzb+1:nzt,nys:nyn,n)
	1259	ENDIF
	1260	ENDDO
	1261	ENDIF
[3737]	1262	ENDIF
[3891]	1263	IF ( bc_dirichlet_s ) THEN
	1264	u(nzb+1:nzt,-1,nxlu:nxr) = nest_offl%u_south(0,nzb+1:nzt,nxlu:nxr)
	1265	v(nzb+1:nzt,0,nxl:nxr) = nest_offl%v_south(0,nzb+1:nzt,nxl:nxr)
	1266	w(nzb+1:nzt-1,-1,nxl:nxr) = nest_offl%w_south(0,nzb+1:nzt-1,nxl:nxr)
	1267	IF ( .NOT. neutral ) pt(nzb+1:nzt,-1,nxl:nxr) = &
	1268	nest_offl%pt_south(0,nzb+1:nzt,nxl:nxr)
	1269	IF ( humidity ) q(nzb+1:nzt,-1,nxl:nxr) = &
	1270	nest_offl%q_south(0,nzb+1:nzt,nxl:nxr)
	1271	IF ( air_chemistry ) THEN
	1272	DO n = 1, UBOUND( chem_species, 1 )
	1273	IF( nest_offl%chem_from_file_s(n) ) THEN
	1274	chem_species(n)%conc(nzb+1:nzt,-1,nxl:nxr) = &
	1275	nest_offl%chem_south(0,nzb+1:nzt,nxl:nxr,n)
	1276	ENDIF
	1277	ENDDO
	1278	ENDIF
[3737]	1279	ENDIF
[3891]	1280	IF ( bc_dirichlet_n ) THEN
	1281	u(nzb+1:nzt,nyn+1,nxlu:nxr) = nest_offl%u_north(0,nzb+1:nzt,nxlu:nxr)
	1282	v(nzb+1:nzt,nyn+1,nxl:nxr) = nest_offl%v_north(0,nzb+1:nzt,nxl:nxr)
	1283	w(nzb+1:nzt-1,nyn+1,nxl:nxr) = nest_offl%w_north(0,nzb+1:nzt-1,nxl:nxr)
	1284	IF ( .NOT. neutral ) pt(nzb+1:nzt,nyn+1,nxl:nxr) = &
	1285	nest_offl%pt_north(0,nzb+1:nzt,nxl:nxr)
	1286	IF ( humidity ) q(nzb+1:nzt,nyn+1,nxl:nxr) = &
	1287	nest_offl%q_north(0,nzb+1:nzt,nxl:nxr)
	1288	IF ( air_chemistry ) THEN
	1289	DO n = 1, UBOUND( chem_species, 1 )
	1290	IF( nest_offl%chem_from_file_n(n) ) THEN
	1291	chem_species(n)%conc(nzb+1:nzt,nyn+1,nxl:nxr) = &
	1292	nest_offl%chem_north(0,nzb+1:nzt,nxl:nxr,n)
	1293	ENDIF
	1294	ENDDO
	1295	ENDIF
[3737]	1296	ENDIF
[3891]	1297	!
	1298	!-- Initialize geostrophic wind components. Actually this is already done in
	1299	!-- init_3d_model when initializing_action = 'inifor', however, in speical
	1300	!-- case of user-defined initialization this will be done here again, in
	1301	!-- order to have a consistent initialization.
	1302	ug(nzb+1:nzt) = nest_offl%ug(0,nzb+1:nzt)
	1303	vg(nzb+1:nzt) = nest_offl%vg(0,nzb+1:nzt)
	1304	!
	1305	!-- Set bottom and top boundary condition for geostrophic wind components
	1306	ug(nzt+1) = ug(nzt)
	1307	vg(nzt+1) = vg(nzt)
	1308	ug(nzb) = ug(nzb+1)
	1309	vg(nzb) = vg(nzb+1)
	1310	ENDIF
[3347]	1311	!
	1312	!-- After boundary data is initialized, mask topography at the
	1313	!-- boundaries for the velocity components.
	1314	u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
	1315	v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
	1316	w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
[3891]	1317	!
	1318	!-- Initial calculation of the boundary layer depth from the prescribed
	1319	!-- boundary data. This is requiered for initialize the synthetic turbulence
	1320	!-- generator correctly.
[3347]	1321	CALL calc_zi
	1322
	1323	!
[3891]	1324	!-- After boundary data is initialized, ensure mass conservation. Not
	1325	!-- necessary in restart runs.
	1326	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
	1327	CALL nesting_offl_mass_conservation
	1328	ENDIF
[3347]	1329
	1330	END SUBROUTINE nesting_offl_init
	1331
	1332	!------------------------------------------------------------------------------!
	1333	! Description:
	1334	!------------------------------------------------------------------------------!
	1335	!> Interpolation function, used to interpolate boundary data in time.
	1336	!------------------------------------------------------------------------------!
	1337	FUNCTION interpolate_in_time( var_t1, var_t2, fac )
	1338
	1339	USE kinds
	1340
	1341	IMPLICIT NONE
	1342
	1343	REAL(wp) :: interpolate_in_time !< time-interpolated boundary value
	1344	REAL(wp) :: var_t1 !< boundary value at t1
	1345	REAL(wp) :: var_t2 !< boundary value at t2
	1346	REAL(wp) :: fac !< interpolation factor
	1347
	1348	interpolate_in_time = ( 1.0_wp - fac ) * var_t1 + fac * var_t2
	1349
	1350	END FUNCTION interpolate_in_time
	1351
	1352
	1353
	1354	END MODULE nesting_offl_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |