Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |