Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |