Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |