Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |