Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |