Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |