Home

Context Navigation

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

Last change on this file since 3352 was 3347, checked in by suehring, 6 years ago
Offline nesting revised and separated from large_scale_forcing_mod; Time-dependent synthetic turbulence generator; bugfixes in USM/LSM radiation model and init_pegrid
File size: 45.4 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	!
	17	! Copyright 1997-2018 Leibniz Universitaet Hannover
	18	!------------------------------------------------------------------------------!
	19	!
	20	! Current revisions:
	21	! ------------------
	22	!
	23	!
	24	! Former revisions:
	25	! -----------------
	26	!
	27	! Initial Revision:
	28	! - separate offline nesting from large_scale_nudging_mod
	29	! - revise offline nesting, adjust for usage of synthetic turbulence generator
	30	! - adjust Rayleigh damping depending on the time-depending atmospheric
	31	! conditions
	32	!
	33	!
	34	! Description:
	35	! ------------
	36	!> Offline nesting in larger-scale models. Boundary conditions for the simulation
	37	!> are read from NetCDF file and are prescribed onto the respective arrays.
	38	!> Further, a mass-flux correction is performed to maintain the mass balance.
	39	!--------------------------------------------------------------------------------!
	40	MODULE nesting_offl_mod
	41
	42	USE arrays_3d, &
	43	ONLY: dzw, e, diss, pt, pt_init, q, q_init, s, u, u_init, ug, v, &
	44	v_init, vg, w, zu, zw
	45
	46	USE control_parameters, &
	47	ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
	48	dt_3d, dz, constant_diffusion, humidity, message_string, &
	49	nesting_offline, neutral, passive_scalar, rans_mode, rans_tke_e,&
	50	time_since_reference_point, volume_flow
	51
	52	USE cpulog, &
	53	ONLY: cpu_log, log_point
	54
	55	USE grid_variables
	56
	57	USE indices, &
	58	ONLY: nbgp, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys, &
	59	nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
	60
	61	USE kinds
	62
	63	USE netcdf_data_input_mod, &
	64	ONLY: nest_offl
	65
	66	USE pegrid
	67
	68	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
	69	!< bulk Richardson number of 0.25
	70
	71	SAVE
	72	PRIVATE
	73	!
	74	!-- Public subroutines
	75	PUBLIC nesting_offl_bc, nesting_offl_header, nesting_offl_init, &
	76	nesting_offl_mass_conservation, nesting_offl_parin
	77	!
	78	!-- Public variables
	79	PUBLIC zi_ribulk
	80
	81	INTERFACE nesting_offl_bc
	82	MODULE PROCEDURE nesting_offl_bc
	83	END INTERFACE nesting_offl_bc
	84
	85	INTERFACE nesting_offl_header
	86	MODULE PROCEDURE nesting_offl_header
	87	END INTERFACE nesting_offl_header
	88
	89	INTERFACE nesting_offl_init
	90	MODULE PROCEDURE nesting_offl_init
	91	END INTERFACE nesting_offl_init
	92
	93	INTERFACE nesting_offl_mass_conservation
	94	MODULE PROCEDURE nesting_offl_mass_conservation
	95	END INTERFACE nesting_offl_mass_conservation
	96
	97	INTERFACE nesting_offl_parin
	98	MODULE PROCEDURE nesting_offl_parin
	99	END INTERFACE nesting_offl_parin
	100
	101	CONTAINS
	102
	103
	104	!------------------------------------------------------------------------------!
	105	! Description:
	106	! ------------
	107	!> In this subroutine a constant mass within the model domain is guaranteed.
	108	!> Larger-scale models may be based on a compressible equation system, which is
	109	!> not consistent with PALMs incompressible equation system. In order to avoid
	110	!> a decrease or increase of mass during the simulation, non-divergent flow
	111	!> through the lateral and top boundaries is compensated by the vertical wind
	112	!> component at the top boundary.
	113	!------------------------------------------------------------------------------!
	114	SUBROUTINE nesting_offl_mass_conservation
	115
	116	IMPLICIT NONE
	117
	118	INTEGER(iwp) :: i !< grid index in x-direction
	119	INTEGER(iwp) :: j !< grid index in y-direction
	120	INTEGER(iwp) :: k !< grid index in z-direction
	121
	122	REAL(wp) :: d_area_t !< inverse of the total area of the horizontal model domain
	123	REAL(wp) :: w_correct !< vertical velocity increment required to compensate non-divergent flow through the boundaries
	124	REAL(wp), DIMENSION(1:3) :: volume_flow_l !< local volume flow
	125
	126
	127	CALL cpu_log( log_point(58), 'offline nesting', 'start' )
	128
	129	volume_flow = 0.0_wp
	130	volume_flow_l = 0.0_wp
	131
	132	d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
	133
	134	IF ( bc_dirichlet_l ) THEN
	135	i = nxl
	136	DO j = nys, nyn
	137	DO k = nzb+1, nzt
	138	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy &
	139	* MERGE( 1.0_wp, 0.0_wp, &
	140	BTEST( wall_flags_0(k,j,i), 1 ) )
	141	ENDDO
	142	ENDDO
	143	ENDIF
	144	IF ( bc_dirichlet_r ) THEN
	145	i = nxr+1
	146	DO j = nys, nyn
	147	DO k = nzb+1, nzt
	148	volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy &
	149	* MERGE( 1.0_wp, 0.0_wp, &
	150	BTEST( wall_flags_0(k,j,i), 1 ) )
	151	ENDDO
	152	ENDDO
	153	ENDIF
	154	IF ( bc_dirichlet_s ) THEN
	155	j = nys
	156	DO i = nxl, nxr
	157	DO k = nzb+1, nzt
	158	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx &
	159	* MERGE( 1.0_wp, 0.0_wp, &
	160	BTEST( wall_flags_0(k,j,i), 2 ) )
	161	ENDDO
	162	ENDDO
	163	ENDIF
	164	IF ( bc_dirichlet_n ) THEN
	165	j = nyn+1
	166	DO i = nxl, nxr
	167	DO k = nzb+1, nzt
	168	volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx &
	169	* MERGE( 1.0_wp, 0.0_wp, &
	170	BTEST( wall_flags_0(k,j,i), 2 ) )
	171	ENDDO
	172	ENDDO
	173	ENDIF
	174	!
	175	!-- Top boundary
	176	k = nzt
	177	DO i = nxl, nxr
	178	DO j = nys, nyn
	179	volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
	180	ENDDO
	181	ENDDO
	182
	183	#if defined( __parallel )
	184	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	185	CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM, &
	186	comm2d, ierr )
	187	#else
	188	volume_flow = volume_flow_l
	189	#endif
	190
	191	w_correct = SUM( volume_flow ) * d_area_t
	192
	193	DO i = nxl, nxr
	194	DO j = nys, nyn
	195	DO k = nzt, nzt + 1
	196	w(k,j,i) = w(k,j,i) + w_correct
	197	ENDDO
	198	ENDDO
	199	ENDDO
	200
	201	CALL cpu_log( log_point(58), 'offline nesting', 'stop' )
	202
	203	END SUBROUTINE nesting_offl_mass_conservation
	204
	205
	206	!------------------------------------------------------------------------------!
	207	! Description:
	208	! ------------
	209	!> Set the lateral and top boundary conditions in case the PALM domain is
	210	!> nested offline in a mesoscale model. Further, average boundary data and
	211	!> determine mean profiles, further used for correct damping in the sponge
	212	!> layer.
	213	!------------------------------------------------------------------------------!
	214	SUBROUTINE nesting_offl_bc
	215
	216	IMPLICIT NONE
	217
	218	INTEGER(iwp) :: i !< running index x-direction
	219	INTEGER(iwp) :: j !< running index y-direction
	220	INTEGER(iwp) :: k !< running index z-direction
	221
	222	REAL(wp) :: fac_dt !< interpolation factor
	223
	224	REAL(wp), DIMENSION(nzb:nzt+1) :: pt_ref !< reference profile for potential temperature
	225	REAL(wp), DIMENSION(nzb:nzt+1) :: pt_ref_l !< reference profile for potential temperature on subdomain
	226	REAL(wp), DIMENSION(nzb:nzt+1) :: q_ref !< reference profile for mixing ratio on subdomain
	227	REAL(wp), DIMENSION(nzb:nzt+1) :: q_ref_l !< reference profile for mixing ratio on subdomain
	228	REAL(wp), DIMENSION(nzb:nzt+1) :: u_ref !< reference profile for u-component on subdomain
	229	REAL(wp), DIMENSION(nzb:nzt+1) :: u_ref_l !< reference profile for u-component on subdomain
	230	REAL(wp), DIMENSION(nzb:nzt+1) :: v_ref !< reference profile for v-component on subdomain
	231	REAL(wp), DIMENSION(nzb:nzt+1) :: v_ref_l !< reference profile for v-component on subdomain
	232
	233	CALL cpu_log( log_point(58), 'offline nesting', 'start' )
	234	!
	235	!-- Set mean profiles, derived from boundary data, to zero
	236	pt_ref = 0.0_wp
	237	q_ref = 0.0_wp
	238	u_ref = 0.0_wp
	239	v_ref = 0.0_wp
	240
	241	pt_ref_l = 0.0_wp
	242	q_ref_l = 0.0_wp
	243	u_ref_l = 0.0_wp
	244	v_ref_l = 0.0_wp
	245	!
	246	!-- Determine interpolation factor and limit it to 1. This is because
	247	!-- t+dt can slightly exceed time(tind_p) before boundary data is updated
	248	!-- again.
	249	fac_dt = ( time_since_reference_point - nest_offl%time(nest_offl%tind) &
	250	+ dt_3d ) / &
	251	( nest_offl%time(nest_offl%tind_p) - nest_offl%time(nest_offl%tind) )
	252	fac_dt = MIN( 1.0_wp, fac_dt )
	253	!
	254	!-- Set boundary conditions of u-, v-, w-component, as well as q, and pt.
	255	!-- Note, boundary values at the left boundary: i=-1 (v,w,pt,q) and
	256	!-- i=0 (u), at the right boundary: i=nxr+1 (all), at the south boundary:
	257	!-- j=-1 (u,w,pt,q) and j=0 (v), at the north boundary: j=nyn+1 (all).
	258	!-- Please note, at the left (for u) and south (for v) boundary, values
	259	!-- for u and v are set also at i/j=-1, since these values are used in
	260	!-- boundary_conditions() to restore prognostic values.
	261	!-- Further, sum up data to calculate mean profiles from boundary data,
	262	!-- used for Rayleigh damping.
	263	IF ( bc_dirichlet_l ) THEN
	264
	265	DO j = nys, nyn
	266	DO k = nzb+1, nzt
	267	u(k,j,0) = interpolate_in_time( nest_offl%u_left(0,k,j), &
	268	nest_offl%u_left(1,k,j), &
	269	fac_dt ) * &
	270	MERGE( 1.0_wp, 0.0_wp, &
	271	BTEST( wall_flags_0(k,j,0), 1 ) )
	272	u(k,j,-1) = u(k,j,0)
	273	ENDDO
	274	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,0)
	275	ENDDO
	276
	277	DO j = nys, nyn
	278	DO k = nzb+1, nzt-1
	279	w(k,j,-1) = interpolate_in_time( nest_offl%w_left(0,k,j), &
	280	nest_offl%w_left(1,k,j), &
	281	fac_dt ) * &
	282	MERGE( 1.0_wp, 0.0_wp, &
	283	BTEST( wall_flags_0(k,j,-1), 3 ) )
	284	ENDDO
	285	ENDDO
	286
	287	DO j = nysv, nyn
	288	DO k = nzb+1, nzt
	289	v(k,j,-1) = interpolate_in_time( nest_offl%v_left(0,k,j), &
	290	nest_offl%v_left(1,k,j), &
	291	fac_dt ) * &
	292	MERGE( 1.0_wp, 0.0_wp, &
	293	BTEST( wall_flags_0(k,j,-1), 2 ) )
	294	ENDDO
	295	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,-1)
	296	ENDDO
	297
	298	IF ( .NOT. neutral ) THEN
	299	DO j = nys, nyn
	300	DO k = nzb+1, nzt
	301	pt(k,j,-1) = interpolate_in_time( nest_offl%pt_left(0,k,j), &
	302	nest_offl%pt_left(1,k,j), &
	303	fac_dt )
	304
	305	ENDDO
	306	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,-1)
	307	ENDDO
	308	ENDIF
	309
	310	IF ( humidity ) THEN
	311	DO j = nys, nyn
	312	DO k = nzb+1, nzt
	313	q(k,j,-1) = interpolate_in_time( nest_offl%q_left(0,k,j), &
	314	nest_offl%q_left(1,k,j), &
	315	fac_dt )
	316
	317	ENDDO
	318	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,-1)
	319	ENDDO
	320	ENDIF
	321
	322	ENDIF
	323
	324	IF ( bc_dirichlet_r ) THEN
	325
	326	DO j = nys, nyn
	327	DO k = nzb+1, nzt
	328	u(k,j,nxr+1) = interpolate_in_time( nest_offl%u_right(0,k,j), &
	329	nest_offl%u_right(1,k,j), &
	330	fac_dt ) * &
	331	MERGE( 1.0_wp, 0.0_wp, &
	332	BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
	333	ENDDO
	334	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,j,nxr+1)
	335	ENDDO
	336	DO j = nys, nyn
	337	DO k = nzb+1, nzt-1
	338	w(k,j,nxr+1) = interpolate_in_time( nest_offl%w_right(0,k,j), &
	339	nest_offl%w_right(1,k,j), &
	340	fac_dt ) * &
	341	MERGE( 1.0_wp, 0.0_wp, &
	342	BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
	343	ENDDO
	344	ENDDO
	345
	346	DO j = nysv, nyn
	347	DO k = nzb+1, nzt
	348	v(k,j,nxr+1) = interpolate_in_time( nest_offl%v_right(0,k,j), &
	349	nest_offl%v_right(1,k,j), &
	350	fac_dt ) * &
	351	MERGE( 1.0_wp, 0.0_wp, &
	352	BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
	353	ENDDO
	354	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,j,nxr+1)
	355	ENDDO
	356
	357	IF ( .NOT. neutral ) THEN
	358	DO j = nys, nyn
	359	DO k = nzb+1, nzt
	360	pt(k,j,nxr+1) = interpolate_in_time( &
	361	nest_offl%pt_right(0,k,j), &
	362	nest_offl%pt_right(1,k,j), &
	363	fac_dt )
	364	ENDDO
	365	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,j,nxr+1)
	366	ENDDO
	367	ENDIF
	368
	369	IF ( humidity ) THEN
	370	DO j = nys, nyn
	371	DO k = nzb+1, nzt
	372	q(k,j,nxr+1) = interpolate_in_time( &
	373	nest_offl%q_right(0,k,j), &
	374	nest_offl%q_right(1,k,j), &
	375	fac_dt )
	376
	377	ENDDO
	378	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,j,nxr+1)
	379	ENDDO
	380	ENDIF
	381
	382	ENDIF
	383
	384	IF ( bc_dirichlet_s ) THEN
	385
	386	DO i = nxl, nxr
	387	DO k = nzb+1, nzt
	388	v(k,0,i) = interpolate_in_time( nest_offl%v_south(0,k,i), &
	389	nest_offl%v_south(1,k,i), &
	390	fac_dt ) * &
	391	MERGE( 1.0_wp, 0.0_wp, &
	392	BTEST( wall_flags_0(k,0,i), 2 ) )
	393	v(k,-1,i) = v(k,0,i)
	394	ENDDO
	395	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,0,i)
	396	ENDDO
	397
	398	DO i = nxl, nxr
	399	DO k = nzb+1, nzt-1
	400	w(k,-1,i) = interpolate_in_time( nest_offl%w_south(0,k,i), &
	401	nest_offl%w_south(1,k,i), &
	402	fac_dt ) * &
	403	MERGE( 1.0_wp, 0.0_wp, &
	404	BTEST( wall_flags_0(k,-1,i), 3 ) )
	405	ENDDO
	406	ENDDO
	407
	408	DO i = nxlu, nxr
	409	DO k = nzb+1, nzt
	410	u(k,-1,i) = interpolate_in_time( nest_offl%u_south(0,k,i), &
	411	nest_offl%u_south(1,k,i), &
	412	fac_dt ) * &
	413	MERGE( 1.0_wp, 0.0_wp, &
	414	BTEST( wall_flags_0(k,-1,i), 1 ) )
	415	ENDDO
	416	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,-1,i)
	417	ENDDO
	418
	419	IF ( .NOT. neutral ) THEN
	420	DO i = nxl, nxr
	421	DO k = nzb+1, nzt
	422	pt(k,-1,i) = interpolate_in_time( &
	423	nest_offl%pt_south(0,k,i), &
	424	nest_offl%pt_south(1,k,i), &
	425	fac_dt )
	426
	427	ENDDO
	428	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,-1,i)
	429	ENDDO
	430	ENDIF
	431
	432	IF ( humidity ) THEN
	433	DO i = nxl, nxr
	434	DO k = nzb+1, nzt
	435	q(k,-1,i) = interpolate_in_time( &
	436	nest_offl%q_south(0,k,i), &
	437	nest_offl%q_south(1,k,i), &
	438	fac_dt )
	439
	440	ENDDO
	441	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,-1,i)
	442	ENDDO
	443	ENDIF
	444
	445	ENDIF
	446
	447	IF ( bc_dirichlet_n ) THEN
	448
	449	DO i = nxl, nxr
	450	DO k = nzb+1, nzt
	451	v(k,nyn+1,i) = interpolate_in_time( nest_offl%v_north(0,k,i), &
	452	nest_offl%v_north(1,k,i), &
	453	fac_dt ) * &
	454	MERGE( 1.0_wp, 0.0_wp, &
	455	BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
	456	ENDDO
	457	v_ref_l(nzb+1:nzt) = v_ref_l(nzb+1:nzt) + v(nzb+1:nzt,nyn+1,i)
	458	ENDDO
	459	DO i = nxl, nxr
	460	DO k = nzb+1, nzt-1
	461	w(k,nyn+1,i) = interpolate_in_time( nest_offl%w_north(0,k,i), &
	462	nest_offl%w_north(1,k,i), &
	463	fac_dt ) * &
	464	MERGE( 1.0_wp, 0.0_wp, &
	465	BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
	466	ENDDO
	467	ENDDO
	468
	469	DO i = nxlu, nxr
	470	DO k = nzb+1, nzt
	471	u(k,nyn+1,i) = interpolate_in_time( nest_offl%u_north(0,k,i), &
	472	nest_offl%u_north(1,k,i), &
	473	fac_dt ) * &
	474	MERGE( 1.0_wp, 0.0_wp, &
	475	BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
	476
	477	ENDDO
	478	u_ref_l(nzb+1:nzt) = u_ref_l(nzb+1:nzt) + u(nzb+1:nzt,nyn+1,i)
	479	ENDDO
	480
	481	IF ( .NOT. neutral ) THEN
	482	DO i = nxl, nxr
	483	DO k = nzb+1, nzt
	484	pt(k,nyn+1,i) = interpolate_in_time( &
	485	nest_offl%pt_north(0,k,i), &
	486	nest_offl%pt_north(1,k,i), &
	487	fac_dt )
	488
	489	ENDDO
	490	pt_ref_l(nzb+1:nzt) = pt_ref_l(nzb+1:nzt) + pt(nzb+1:nzt,nyn+1,i)
	491	ENDDO
	492	ENDIF
	493
	494	IF ( humidity ) THEN
	495	DO i = nxl, nxr
	496	DO k = nzb+1, nzt
	497	q(k,nyn+1,i) = interpolate_in_time( &
	498	nest_offl%q_north(0,k,i), &
	499	nest_offl%q_north(1,k,i), &
	500	fac_dt )
	501
	502	ENDDO
	503	q_ref_l(nzb+1:nzt) = q_ref_l(nzb+1:nzt) + q(nzb+1:nzt,nyn+1,i)
	504	ENDDO
	505	ENDIF
	506
	507	ENDIF
	508	!
	509	!-- Top boundary
	510	DO i = nxlu, nxr
	511	DO j = nys, nyn
	512	u(nzt+1,j,i) = interpolate_in_time( nest_offl%u_top(0,j,i), &
	513	nest_offl%u_top(1,j,i), &
	514	fac_dt ) * &
	515	MERGE( 1.0_wp, 0.0_wp, &
	516	BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
	517	u_ref_l(nzt+1) = u_ref_l(nzt+1) + u(nzt+1,j,i)
	518	ENDDO
	519	ENDDO
	520
	521	DO i = nxl, nxr
	522	DO j = nysv, nyn
	523	v(nzt+1,j,i) = interpolate_in_time( nest_offl%v_top(0,j,i), &
	524	nest_offl%v_top(1,j,i), &
	525	fac_dt ) * &
	526	MERGE( 1.0_wp, 0.0_wp, &
	527	BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
	528	v_ref_l(nzt+1) = v_ref_l(nzt+1) + v(nzt+1,j,i)
	529	ENDDO
	530	ENDDO
	531
	532	DO i = nxl, nxr
	533	DO j = nys, nyn
	534	w(nzt,j,i) = interpolate_in_time( nest_offl%w_top(0,j,i), &
	535	nest_offl%w_top(1,j,i), &
	536	fac_dt ) * &
	537	MERGE( 1.0_wp, 0.0_wp, &
	538	BTEST( wall_flags_0(nzt,j,i), 3 ) )
	539	w(nzt+1,j,i) = w(nzt,j,i)
	540	ENDDO
	541	ENDDO
	542
	543
	544	IF ( .NOT. neutral ) THEN
	545	DO i = nxl, nxr
	546	DO j = nys, nyn
	547	pt(nzt+1,j,i) = interpolate_in_time( nest_offl%pt_top(0,j,i), &
	548	nest_offl%pt_top(1,j,i), &
	549	fac_dt )
	550	pt_ref_l(nzt+1) = pt_ref_l(nzt+1) + pt(nzt+1,j,i)
	551	ENDDO
	552	ENDDO
	553	ENDIF
	554
	555	IF ( humidity ) THEN
	556	DO i = nxl, nxr
	557	DO j = nys, nyn
	558	q(nzt+1,j,i) = interpolate_in_time( nest_offl%q_top(0,j,i), &
	559	nest_offl%q_top(1,j,i), &
	560	fac_dt )
	561	q_ref_l(nzt+1) = q_ref_l(nzt+1) + q(nzt+1,j,i)
	562	ENDDO
	563	ENDDO
	564	ENDIF
	565	!
	566	!-- Moreover, set Neumann boundary condition for subgrid-scale TKE,
	567	!-- passive scalar, dissipation, and chemical species if required
	568	IF ( rans_mode .AND. rans_tke_e ) THEN
	569	IF ( bc_dirichlet_l ) diss(:,:,nxl-1) = diss(:,:,nxl)
	570	IF ( bc_dirichlet_r ) diss(:,:,nxr+1) = diss(:,:,nxr)
	571	IF ( bc_dirichlet_s ) diss(:,nys-1,:) = diss(:,nys,:)
	572	IF ( bc_dirichlet_n ) diss(:,nyn+1,:) = diss(:,nyn,:)
	573	ENDIF
	574	IF ( .NOT. constant_diffusion ) THEN
	575	IF ( bc_dirichlet_l ) e(:,:,nxl-1) = e(:,:,nxl)
	576	IF ( bc_dirichlet_r ) e(:,:,nxr+1) = e(:,:,nxr)
	577	IF ( bc_dirichlet_s ) e(:,nys-1,:) = e(:,nys,:)
	578	IF ( bc_dirichlet_n ) e(:,nyn+1,:) = e(:,nyn,:)
	579	e(nzt+1,:,:) = e(nzt,:,:)
	580	ENDIF
	581	IF ( passive_scalar ) THEN
	582	IF ( bc_dirichlet_l ) s(:,:,nxl-1) = s(:,:,nxl)
	583	IF ( bc_dirichlet_r ) s(:,:,nxr+1) = s(:,:,nxr)
	584	IF ( bc_dirichlet_s ) s(:,nys-1,:) = s(:,nys,:)
	585	IF ( bc_dirichlet_n ) s(:,nyn+1,:) = s(:,nyn,:)
	586	ENDIF
	587
	588	CALL exchange_horiz( u, nbgp )
	589	CALL exchange_horiz( v, nbgp )
	590	CALL exchange_horiz( w, nbgp )
	591	IF ( .NOT. neutral ) CALL exchange_horiz( pt, nbgp )
	592	IF ( humidity ) CALL exchange_horiz( q, nbgp )
	593	!
	594	!-- In case of Rayleigh damping, where the profiles u_init, v_init
	595	!-- q_init and pt_init are still used, update these profiles from the
	596	!-- averaged boundary data.
	597	!-- But first, average these data.
	598	#if defined( __parallel )
	599	CALL MPI_ALLREDUCE( u_ref_l, u_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM, &
	600	comm2d, ierr )
	601	CALL MPI_ALLREDUCE( v_ref_l, v_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM, &
	602	comm2d, ierr )
	603	IF ( humidity ) THEN
	604	CALL MPI_ALLREDUCE( q_ref_l, q_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM, &
	605	comm2d, ierr )
	606	ENDIF
	607	IF ( .NOT. neutral ) THEN
	608	CALL MPI_ALLREDUCE( pt_ref_l, pt_ref, nzt+1-nzb+1, MPI_REAL, MPI_SUM,&
	609	comm2d, ierr )
	610	ENDIF
	611	#else
	612	u_ref = u_ref_l
	613	v_ref = v_ref_l
	614	IF ( humidity ) q_ref = q_ref_l
	615	IF ( .NOT. neutral ) pt_ref = pt_ref_l
	616	#endif
	617	!
	618	!-- Average data. Note, reference profiles up to nzt are derived from lateral
	619	!-- boundaries, at the model top it is derived from the top boundary. Thus,
	620	!-- number of input data is different from nzb:nzt compared to nzt+1.
	621	!-- Derived from lateral boundaries.
	622	u_ref(nzb:nzt) = u_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx ), &
	623	KIND = wp )
	624	v_ref(nzb:nzt) = v_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + nx + 1 ), &
	625	KIND = wp )
	626	IF ( humidity ) &
	627	q_ref(nzb:nzt) = q_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), &
	628	KIND = wp )
	629	IF ( .NOT. neutral ) &
	630	pt_ref(nzb:nzt) = pt_ref(nzb:nzt) / REAL( 2.0_wp * ( ny + 1 + nx + 1 ), &
	631	KIND = wp )
	632	!
	633	!-- Derived from top boundary.
	634	u_ref(nzt+1) = u_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx ), KIND = wp )
	635	v_ref(nzt+1) = v_ref(nzt+1) / REAL( ( ny ) * ( nx + 1 ), KIND = wp )
	636	IF ( humidity ) &
	637	q_ref(nzt+1) = q_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp )
	638	IF ( .NOT. neutral ) &
	639	pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ), KIND = wp )
	640	!
	641	!-- Write onto init profiles, which are used for damping
	642	u_init = u_ref
	643	v_init = v_ref
	644	IF ( humidity ) q_init = q_ref
	645	IF ( .NOT. neutral ) pt_init = pt_ref
	646	!
	647	!-- Set bottom boundary condition
	648	IF ( humidity ) q_init(nzb) = q_init(nzb+1)
	649	IF ( .NOT. neutral ) pt_init(nzb) = pt_init(nzb+1)
	650	!
	651	!-- Further, adjust Rayleigh damping height in case of time-changing conditions.
	652	!-- Therefore, calculate boundary-layer depth first.
	653	CALL calc_zi
	654	CALL adjust_sponge_layer
	655
	656	!
	657	!-- Update geostrophic wind components from dynamic input file.
	658	!-- Note, at the moment this is not available. As a workaround, use reference
	659	!-- profiles derived from the boundary values.
	660	!-- This will be revised as soon as geostrophic wind components are provided by
	661	!-- dynamic input file.
	662	ug = u_init
	663	vg = v_init
	664
	665	CALL cpu_log( log_point(58), 'offline nesting', 'stop' )
	666
	667	END SUBROUTINE nesting_offl_bc
	668
	669	!------------------------------------------------------------------------------!
	670	! Description:
	671	!------------------------------------------------------------------------------!
	672	!> Calculates the boundary-layer depth from the boundary data, according to
	673	!> bulk-Richardson criterion.
	674	!------------------------------------------------------------------------------!
	675	SUBROUTINE calc_zi
	676
	677	USE basic_constants_and_equations_mod, &
	678	ONLY: g
	679
	680	USE kinds
	681
	682	USE surface_mod, &
	683	ONLY: get_topography_top_index, get_topography_top_index_ji
	684
	685	IMPLICIT NONE
	686
	687	INTEGER(iwp) :: i !< loop index in x-direction
	688	INTEGER(iwp) :: j !< loop index in y-direction
	689	INTEGER(iwp) :: k !< loop index in z-direction
	690	INTEGER(iwp) :: k_surface !< topography top index in z-direction
	691
	692	REAL(wp) :: ri_bulk !< bulk Richardson number
	693	REAL(wp) :: ri_bulk_crit = 0.25_wp !< critical bulk Richardson number
	694	REAL(wp) :: topo_max !< maximum topography level in model domain
	695	REAL(wp) :: topo_max_l !< maximum topography level in subdomain
	696	REAL(wp) :: u_comp !< u-component
	697	REAL(wp) :: v_comp !< v-component
	698	REAL(wp) :: vpt_surface !< near-surface virtual potential temperature
	699	REAL(wp) :: zi_l !< mean boundary-layer depth on subdomain
	700	REAL(wp) :: zi_local !< local boundary-layer depth
	701
	702	REAL(wp), DIMENSION(nzb:nzt+1) :: vpt_col !< vertical profile of virtual potential temperature at (j,i)-grid point
	703
	704
	705	!
	706	!-- Calculate mean boundary-layer height from boundary data.
	707	!-- Start with the left and right boundaries.
	708	zi_l = 0.0_wp
	709	IF ( bc_dirichlet_l .OR. bc_dirichlet_r ) THEN
	710	!
	711	!-- Determine index along x. Please note, index indicates boundary
	712	!-- grid point for scalars.
	713	i = MERGE( -1, nxr + 1, bc_dirichlet_l )
	714
	715	DO j = nys, nyn
	716	!
	717	!-- Determine topography top index at current (j,i) index
	718	k_surface = get_topography_top_index_ji( j, i, 's' )
	719	!
	720	!-- Pre-compute surface virtual temperature. Therefore, use 2nd
	721	!-- prognostic level according to Heinze et al. (2017).
	722	IF ( humidity ) THEN
	723	vpt_surface = pt(k_surface+2,j,i) * &
	724	( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
	725	vpt_col = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
	726	ELSE
	727	vpt_surface = pt(k_surface+2,j,i)
	728	vpt_col = pt(:,j,i)
	729	ENDIF
	730	!
	731	!-- Calculate local boundary layer height from bulk Richardson number,
	732	!-- i.e. the height where the bulk Richardson number exceeds its
	733	!-- critical value of 0.25 (according to Heinze et al., 2017).
	734	!-- Note, no interpolation of u- and v-component is made, as both
	735	!-- are mainly mean inflow profiles with very small spatial variation.
	736	zi_local = 0.0_wp
	737	DO k = k_surface+1, nzt
	738	u_comp = MERGE( u(k,j,i+1), u(k,j,i), bc_dirichlet_l )
	739	v_comp = v(k,j,i)
	740	ri_bulk = zu(k) * g / vpt_surface * &
	741	( vpt_col(k) - vpt_surface ) / &
	742	( u_comp * u_comp + v_comp * v_comp )
	743
	744	IF ( zi_local == 0.0_wp .AND. ri_bulk > ri_bulk_crit ) &
	745	zi_local = zu(k)
	746	ENDDO
	747	!
	748	!-- Assure that the minimum local boundary-layer depth is at least at
	749	!-- the second vertical grid level.
	750	zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )
	751
	752	ENDDO
	753
	754	ENDIF
	755	!
	756	!-- Do the same at the north and south boundaries.
	757	IF ( bc_dirichlet_s .OR. bc_dirichlet_n ) THEN
	758
	759	j = MERGE( -1, nyn + 1, bc_dirichlet_s )
	760
	761	DO i = nxl, nxr
	762	k_surface = get_topography_top_index_ji( j, i, 's' )
	763
	764	IF ( humidity ) THEN
	765	vpt_surface = pt(k_surface+2,j,i) * &
	766	( 1.0_wp + 0.61_wp * q(k_surface+2,j,i) )
	767	vpt_col = pt(:,j,i) * ( 1.0_wp + 0.61_wp * q(:,j,i) )
	768	ELSE
	769	vpt_surface = pt(k_surface+2,j,i)
	770	vpt_col = pt(:,j,i)
	771	ENDIF
	772
	773	zi_local = 0.0_wp
	774	DO k = k_surface+1, nzt
	775	u_comp = u(k,j,i)
	776	v_comp = MERGE( v(k,j+1,i), v(k,j,i), bc_dirichlet_s )
	777	ri_bulk = zu(k) * g / vpt_surface * &
	778	( vpt_col(k) - vpt_surface ) / &
	779	( u_comp * u_comp + v_comp * v_comp )
	780
	781	IF ( zi_local == 0.0_wp .AND. ri_bulk > 0.25_wp ) &
	782	zi_local = zu(k)
	783	ENDDO
	784	zi_l = zi_l + MAX( zi_local, zu(k_surface+2) )
	785
	786	ENDDO
	787
	788	ENDIF
	789
	790	#if defined( __parallel )
	791	CALL MPI_ALLREDUCE( zi_l, zi_ribulk, 1, MPI_REAL, MPI_SUM, &
	792	comm2d, ierr )
	793	#else
	794	zi_ribulk = zi_l
	795	#endif
	796	zi_ribulk = zi_ribulk / REAL( 2 * nx + 2 * ny, KIND = wp )
	797	!
	798	!-- Finally, check if boundary layer depth is not below the any topography.
	799	!-- zi_ribulk will be used to adjust rayleigh damping height, i.e. the
	800	!-- lower level of the sponge layer, as well as to adjust the synthetic
	801	!-- turbulence generator accordingly. If Rayleigh damping would be applied
	802	!-- near buildings, etc., this would spoil the simulation results.
	803	topo_max_l = zw(MAXVAL( get_topography_top_index( 's' )))
	804
	805	#if defined( __parallel )
	806	CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX, &
	807	comm2d, ierr )
	808	#else
	809	topo_max = topo_max_l
	810	#endif
	811
	812	zi_ribulk = MAX( zi_ribulk, topo_max )
	813
	814	END SUBROUTINE calc_zi
	815
	816
	817	!------------------------------------------------------------------------------!
	818	! Description:
	819	!------------------------------------------------------------------------------!
	820	!> Adjust the height where the rayleigh damping starts, i.e. the lower level
	821	!> of the sponge layer.
	822	!------------------------------------------------------------------------------!
	823	SUBROUTINE adjust_sponge_layer
	824
	825	USE arrays_3d, &
	826	ONLY: rdf, rdf_sc, zu
	827
	828	USE basic_constants_and_equations_mod, &
	829	ONLY: pi
	830
	831	USE control_parameters, &
	832	ONLY: rayleigh_damping_height, rayleigh_damping_factor
	833
	834	USE kinds
	835
	836	IMPLICIT NONE
	837
	838	INTEGER(iwp) :: k !< loop index in z-direction
	839
	840	REAL(wp) :: rdh !< updated Rayleigh damping height
	841
	842
	843	IF ( rayleigh_damping_height > 0.0_wp .AND. &
	844	rayleigh_damping_factor > 0.0_wp ) THEN
	845	!
	846	!-- Update Rayleigh-damping height and re-calculate height-depending
	847	!-- damping coefficients.
	848	!-- Assure that rayleigh damping starts well above the boundary layer.
	849	rdh = MIN( MAX( zi_ribulk * 1.3_wp, 10.0_wp * dz(1) ), &
	850	0.8_wp * zu(nzt), rayleigh_damping_height )
	851	!
	852	!-- Update Rayleigh damping factor
	853	DO k = nzb+1, nzt
	854	IF ( zu(k) >= rdh ) THEN
	855	rdf(k) = rayleigh_damping_factor * &
	856	( SIN( pi * 0.5_wp * ( zu(k) - rdh ) &
	857	/ ( zu(nzt) - rdh ) ) &
	858	)**2
	859	ENDIF
	860	ENDDO
	861	rdf_sc = rdf
	862
	863	ENDIF
	864
	865	END SUBROUTINE adjust_sponge_layer
	866
	867	!------------------------------------------------------------------------------!
	868	! Description:
	869	! ------------
	870	!> Performs consistency checks
	871	!------------------------------------------------------------------------------!
	872	SUBROUTINE nesting_offl_check_parameters
	873
	874	IMPLICIT NONE
	875	!
	876	!-- Perform checks
	877
	878
	879	END SUBROUTINE nesting_offl_check_parameters
	880
	881	!------------------------------------------------------------------------------!
	882	! Description:
	883	! ------------
	884	!> Reads the parameter list nesting_offl_parameters
	885	!------------------------------------------------------------------------------!
	886	SUBROUTINE nesting_offl_parin
	887
	888	IMPLICIT NONE
	889
	890	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
	891
	892
	893	NAMELIST /nesting_offl_parameters/ nesting_offline
	894
	895	line = ' '
	896
	897	!
	898	!-- Try to find stg package
	899	REWIND ( 11 )
	900	line = ' '
	901	DO WHILE ( INDEX( line, '&nesting_offl_parameters' ) == 0 )
	902	READ ( 11, '(A)', END=20 ) line
	903	ENDDO
	904	BACKSPACE ( 11 )
	905
	906	!
	907	!-- Read namelist
	908	READ ( 11, nesting_offl_parameters, ERR = 10, END = 20 )
	909
	910	GOTO 20
	911
	912	10 BACKSPACE( 11 )
	913	READ( 11 , '(A)') line
	914	CALL parin_fail_message( 'nesting_offl_parameters', line )
	915
	916	20 CONTINUE
	917
	918
	919	END SUBROUTINE nesting_offl_parin
	920
	921	!------------------------------------------------------------------------------!
	922	! Description:
	923	! ------------
	924	!> Writes information about offline nesting into HEADER file
	925	!------------------------------------------------------------------------------!
	926	SUBROUTINE nesting_offl_header ( io )
	927
	928	IMPLICIT NONE
	929
	930	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
	931
	932	WRITE ( io, 1 )
	933	IF ( nesting_offline ) THEN
	934	WRITE ( io, 3 )
	935	ELSE
	936	WRITE ( io, 2 )
	937	ENDIF
	938
	939	1 FORMAT (//' Offline nesting in COSMO model:'/ &
	940	' -------------------------------'/)
	941	2 FORMAT (' --> No offlince nesting is used (default) ')
	942	3 FORMAT (' --> Offlince nesting is used. Boundary data is read from dynamic input file ')
	943
	944	END SUBROUTINE nesting_offl_header
	945
	946	!------------------------------------------------------------------------------!
	947	! Description:
	948	! ------------
	949	!> Allocate arrays used to read boundary data from NetCDF file and initialize
	950	!> boundary data.
	951	!------------------------------------------------------------------------------!
	952	SUBROUTINE nesting_offl_init
	953
	954	USE netcdf_data_input_mod, &
	955	ONLY: netcdf_data_input_offline_nesting
	956
	957	IMPLICIT NONE
	958
	959
	960	!-- Allocate arrays for geostrophic wind components. Arrays will
	961	!-- incorporate 2 time levels in order to interpolate in between. Please
	962	!-- note, forcing using geostrophic wind components is only required in
	963	!-- case of cyclic boundary conditions.
	964	! IF ( bc_lr_cyc .AND. bc_ns_cyc ) THEN
	965	! ALLOCATE( nest_offl%ug(0:1,nzb:nzt+1) )
	966	! ALLOCATE( nest_offl%vg(0:1,nzb:nzt+1) )
	967	! ENDIF
	968	!
	969	!-- Allocate arrays for reading boundary values. Arrays will incorporate 2
	970	!-- time levels in order to interpolate in between.
	971	IF ( bc_dirichlet_l ) THEN
	972	ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn) )
	973	ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
	974	ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
	975	IF ( humidity ) ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn) )
	976	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
	977	ENDIF
	978	IF ( bc_dirichlet_r ) THEN
	979	ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn) )
	980	ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
	981	ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
	982	IF ( humidity ) ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn) )
	983	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
	984	ENDIF
	985	IF ( bc_dirichlet_n ) THEN
	986	ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
	987	ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr) )
	988	ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
	989	IF ( humidity ) ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr) )
	990	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
	991	ENDIF
	992	IF ( bc_dirichlet_s ) THEN
	993	ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
	994	ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr) )
	995	ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr) )
	996	IF ( humidity ) ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr) )
	997	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
	998	ENDIF
	999
	1000	ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
	1001	ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
	1002	ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr) )
	1003	IF ( humidity ) ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr) )
	1004	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
	1005
	1006	!
	1007	!-- Read COSMO data at lateral and top boundaries
	1008	CALL netcdf_data_input_offline_nesting
	1009	!
	1010	!-- Initialize boundary data
	1011	IF ( bc_dirichlet_l ) THEN
	1012	u(nzb+1:nzt,nys:nyn,0) = nest_offl%u_left(0,nzb+1:nzt,nys:nyn)
	1013	v(nzb+1:nzt,nysv:nyn,-1) = nest_offl%v_left(0,nzb+1:nzt,nysv:nyn)
	1014	w(nzb+1:nzt-1,nys:nyn,-1) = nest_offl%w_left(0,nzb+1:nzt-1,nys:nyn)
	1015	IF ( .NOT. neutral ) pt(nzb+1:nzt,nys:nyn,-1) = &
	1016	nest_offl%pt_left(0,nzb+1:nzt,nys:nyn)
	1017	IF ( humidity ) q(nzb+1:nzt,nys:nyn,-1) = &
	1018	nest_offl%q_left(0,nzb+1:nzt,nys:nyn)
	1019	ENDIF
	1020	IF ( bc_dirichlet_r ) THEN
	1021	u(nzb+1:nzt,nys:nyn,nxr+1) = nest_offl%u_right(0,nzb+1:nzt,nys:nyn)
	1022	v(nzb+1:nzt,nysv:nyn,nxr+1) = nest_offl%v_right(0,nzb+1:nzt,nysv:nyn)
	1023	w(nzb+1:nzt-1,nys:nyn,nxr+1) = nest_offl%w_right(0,nzb+1:nzt-1,nys:nyn)
	1024	IF ( .NOT. neutral ) pt(nzb+1:nzt,nys:nyn,nxr+1) = &
	1025	nest_offl%pt_right(0,nzb+1:nzt,nys:nyn)
	1026	IF ( humidity ) q(nzb+1:nzt,nys:nyn,nxr+1) = &
	1027	nest_offl%q_right(0,nzb+1:nzt,nys:nyn)
	1028	ENDIF
	1029	IF ( bc_dirichlet_s ) THEN
	1030	u(nzb+1:nzt,-1,nxlu:nxr) = nest_offl%u_south(0,nzb+1:nzt,nxlu:nxr)
	1031	v(nzb+1:nzt,0,nxl:nxr) = nest_offl%v_south(0,nzb+1:nzt,nxl:nxr)
	1032	w(nzb+1:nzt-1,-1,nxl:nxr) = nest_offl%w_south(0,nzb+1:nzt-1,nxl:nxr)
	1033	IF ( .NOT. neutral ) pt(nzb+1:nzt,-1,nxl:nxr) = &
	1034	nest_offl%pt_south(0,nzb+1:nzt,nxl:nxr)
	1035	IF ( humidity ) q(nzb+1:nzt,-1,nxl:nxr) = &
	1036	nest_offl%q_south(0,nzb+1:nzt,nxl:nxr)
	1037	ENDIF
	1038	IF ( bc_dirichlet_n ) THEN
	1039	u(nzb+1:nzt,nyn+1,nxlu:nxr) = nest_offl%u_north(0,nzb+1:nzt,nxlu:nxr)
	1040	v(nzb+1:nzt,nyn+1,nxl:nxr) = nest_offl%v_north(0,nzb+1:nzt,nxl:nxr)
	1041	w(nzb+1:nzt-1,nyn+1,nxl:nxr) = nest_offl%w_north(0,nzb+1:nzt-1,nxl:nxr)
	1042	IF ( .NOT. neutral ) pt(nzb+1:nzt,nyn+1,nxl:nxr) = &
	1043	nest_offl%pt_north(0,nzb+1:nzt,nxl:nxr)
	1044	IF ( humidity ) q(nzb+1:nzt,nyn+1,nxl:nxr) = &
	1045	nest_offl%q_north(0,nzb+1:nzt,nxl:nxr)
	1046	ENDIF
	1047	!
	1048	!-- Initial calculation of the boundary layer depth from the prescribed
	1049	!-- boundary data. This is requiered for initialize the synthetic turbulence
	1050	!-- generator correctly.
	1051	CALL calc_zi
	1052
	1053	!
	1054	!-- After boundary data is initialized, mask topography at the
	1055	!-- boundaries for the velocity components.
	1056	u = MERGE( u, 0.0_wp, BTEST( wall_flags_0, 1 ) )
	1057	v = MERGE( v, 0.0_wp, BTEST( wall_flags_0, 2 ) )
	1058	w = MERGE( w, 0.0_wp, BTEST( wall_flags_0, 3 ) )
	1059
	1060	CALL calc_zi
	1061
	1062	!
	1063	!-- After boundary data is initialized, ensure mass conservation
	1064	CALL nesting_offl_mass_conservation
	1065
	1066	END SUBROUTINE nesting_offl_init
	1067
	1068	!------------------------------------------------------------------------------!
	1069	! Description:
	1070	!------------------------------------------------------------------------------!
	1071	!> Interpolation function, used to interpolate boundary data in time.
	1072	!------------------------------------------------------------------------------!
	1073	FUNCTION interpolate_in_time( var_t1, var_t2, fac )
	1074
	1075	USE kinds
	1076
	1077	IMPLICIT NONE
	1078
	1079	REAL(wp) :: interpolate_in_time !< time-interpolated boundary value
	1080	REAL(wp) :: var_t1 !< boundary value at t1
	1081	REAL(wp) :: var_t2 !< boundary value at t2
	1082	REAL(wp) :: fac !< interpolation factor
	1083
	1084	interpolate_in_time = ( 1.0_wp - fac ) * var_t1 + fac * var_t2
	1085
	1086	END FUNCTION interpolate_in_time
	1087
	1088
	1089
	1090	END MODULE nesting_offl_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |