Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |