Home

Context Navigation

source: palm/trunk/SOURCE/large_scale_forcing_nudging_mod.f90 @ 3182

Last change on this file since 3182 was 3182, checked in by suehring, 6 years ago
New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed
Property svn:keywords set to `Id`
File size: 78.2 KB

Rev	Line
[2696]	1	!> @file large_scale_forcing_nudging_mod.f90
[2320]	2	!------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[2320]	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	!
[2718]	17	! Copyright 1997-2018 Leibniz Universitaet Hannover
[2320]	18	!------------------------------------------------------------------------------!
	19	!
	20	! Current revisions:
	21	! ------------------
[3182]	22	! * Adjustment to new Inifor version:
	23	! - No vertical interpolation/extrapolation of lateral boundary data required
	24	! any more (Inifor can treat grid stretching now
	25	! - Revise initialization in case of COSMO forcing
	26	! * Rename variables and subroutines for offline nesting
[2320]	27	!
	28	! Former revisions:
	29	! -----------------
	30	! $Id: large_scale_forcing_nudging_mod.f90 3182 2018-07-27 13:36:03Z suehring $
[3049]	31	! Error messages revised
	32	!
	33	! 3045 2018-05-28 07:55:41Z Giersch
[3045]	34	! Error messages revised
	35	!
	36	! 3026 2018-05-22 10:30:53Z schwenkel
[3026]	37	! Changed the name specific humidity to mixing ratio, since we are computing
	38	! mixing ratios.
	39	!
	40	! 2970 2018-04-13 15:09:23Z suehring
[2970]	41	! Bugfix in old large-scale forcing mode
	42	!
	43	! 2938 2018-03-27 15:52:42Z suehring
[2938]	44	! Further improvements for nesting in larger-scale model
	45	!
	46	! 2863 2018-03-08 11:36:25Z suehring
[2716]	47	! Corrected "Former revisions" section
	48	!
	49	! 2696 2017-12-14 17:12:51Z kanani
	50	! Change in file header (GPL part)
[2696]	51	! Forcing with larger-scale models implemented (MS)
	52	!
	53	! 2342 2017-08-08 11:00:43Z boeske
[2342]	54	! fixed check if surface forcing data is available until end of simulation
	55	!
	56	! 2320 2017-07-21 12:47:43Z suehring
	57	! initial revision
[2320]	58	!
	59	! Description:
	60	! ------------
	61	!> Calculates large scale forcings (geostrophic wind and subsidence velocity) as
	62	!> well as surfaces fluxes dependent on time given in an external file (LSF_DATA).
	63	!> Moreover, module contains nudging routines, where u, v, pt and q are nudged
	64	!> to given profiles on a relaxation timescale tnudge.
[2342]	65	!> Profiles are read in from NUDGING_DATA.
[2320]	66	!> Code is based on Neggers et al. (2012) and also in parts on DALES and UCLA-LES.
	67	!> @todo: Revise reading of ASCII-files
	68	!> @todo: Remove unused variables and control flags
	69	!> @todo: Revise large-scale facing of surface variables
	70	!> @todo: Revise control flags lsf_exception, lsf_surf, lsf_vert, etc.
	71	!--------------------------------------------------------------------------------!
	72	MODULE lsf_nudging_mod
	73
	74	USE arrays_3d, &
[3182]	75	ONLY: dzw, e, diss, heatflux_input_conversion, pt, pt_init, q, &
	76	q_init, s, tend, u, u_init, ug, v, v_init, vg, w, w_subs, &
[2696]	77	waterflux_input_conversion, zu, zw
[2320]	78
	79	USE control_parameters, &
[3182]	80	ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
	81	bc_lr, bc_ns, bc_pt_b, bc_q_b, constant_diffusion, &
[2970]	82	constant_heatflux, constant_waterflux, &
[3182]	83	data_output_pr, dt_3d, end_time, &
[2970]	84	humidity, initializing_actions, intermediate_timestep_count, &
	85	ibc_pt_b, ibc_q_b, &
[2320]	86	large_scale_forcing, large_scale_subsidence, lsf_surf, lsf_vert,&
[3182]	87	lsf_exception, message_string, nesting_offline, neutral, &
	88	nudging, passive_scalar, pt_surface, ocean, q_surface, &
	89	surface_heatflux, surface_pressure, surface_waterflux, &
	90	topography, use_subsidence_tendencies
[2320]	91
[2696]	92	USE grid_variables
[2320]	93
	94	USE indices, &
[2696]	95	ONLY: nbgp, ngp_sums_ls, nx, nxl, nxlg, nxlu, nxr, nxrg, ny, nys, &
	96	nysv, nysg, nyn, nyng, nzb, nz, nzt, wall_flags_0
[2320]	97
	98	USE kinds
	99
[3182]	100	USE netcdf_data_input_mod, &
	101	ONLY: nest_offl
	102
	103	USE pegrid
	104
[2320]	105	USE surface_mod, &
	106	ONLY: surf_def_h, surf_lsm_h, surf_usm_h
	107
	108	USE statistics, &
	109	ONLY: hom, statistic_regions, sums_ls_l, weight_substep
	110
	111	INTEGER(iwp) :: nlsf = 1000 !< maximum number of profiles in LSF_DATA (large scale forcing)
	112	INTEGER(iwp) :: ntnudge = 1000 !< maximum number of profiles in NUDGING_DATA (nudging)
	113
[2696]	114	REAL(wp) :: d_area_t
	115
[2320]	116	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ptnudge !< vertical profile of pot. temperature interpolated to vertical grid (nudging)
[3026]	117	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: qnudge !< vertical profile of water vapor mixing ratio interpolated to vertical grid (nudging)
[2320]	118	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tnudge !< vertical profile of nudging time scale interpolated to vertical grid (nudging)
	119	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: td_lsa_lpt !< temperature tendency due to large scale advection (large scale forcing)
[3026]	120	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: td_lsa_q !< water vapor mixing ratio tendency due to large scale advection (large scale forcing)
[2320]	121	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: td_sub_lpt !< temperature tendency due to subsidence/ascent (large scale forcing)
[3026]	122	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: td_sub_q !< water vapor mixing ratio tendency due to subsidence/ascent (large scale forcing)
[2320]	123	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ug_vert !< vertical profile of geostrophic wind component in x-direction interpolated to vertical grid (large scale forcing)
	124	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: unudge !< vertical profile of wind component in x-direction interpolated to vertical grid (nudging)
	125	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vnudge !< vertical profile of wind component in y-direction interpolated to vertical grid (nudging)
	126	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vg_vert !< vertical profile of geostrophic wind component in y-direction interpolated to vertical grid (large scale forcing)
	127	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: wnudge !< vertical profile of subsidence/ascent velocity interpolated to vertical grid (nudging) ???
	128	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: wsubs_vert !< vertical profile of wind component in z-direction interpolated to vertical grid (nudging) ???
	129
	130	REAL(wp), DIMENSION(:), ALLOCATABLE :: shf_surf !< time-dependent surface sensible heat flux (large scale forcing)
	131	REAL(wp), DIMENSION(:), ALLOCATABLE :: timenudge !< times at which vertical profiles are defined in NUDGING_DATA (nudging)
	132	REAL(wp), DIMENSION(:), ALLOCATABLE :: time_surf !< times at which surface values/fluxes are defined in LSF_DATA (large scale forcing)
	133	REAL(wp), DIMENSION(:), ALLOCATABLE :: time_vert !< times at which vertical profiles are defined in LSF_DATA (large scale forcing)
	134
	135	REAL(wp), DIMENSION(:), ALLOCATABLE :: tmp_tnudge !< current nudging time scale
	136
	137	REAL(wp), DIMENSION(:), ALLOCATABLE :: p_surf !< time-dependent surface pressure (large scale forcing)
	138	REAL(wp), DIMENSION(:), ALLOCATABLE :: pt_surf !< time-dependent surface temperature (large scale forcing)
	139	REAL(wp), DIMENSION(:), ALLOCATABLE :: qsws_surf !< time-dependent surface latent heat flux (large scale forcing)
[3026]	140	REAL(wp), DIMENSION(:), ALLOCATABLE :: q_surf !< time-dependent surface water vapor mixing ratio (large scale forcing)
[2320]	141
	142	SAVE
	143	PRIVATE
	144	!
	145	!-- Public subroutines
[3182]	146	PUBLIC calc_tnudge, ls_forcing_surf, ls_forcing_vert, ls_advec, lsf_init, &
[2320]	147	lsf_nudging_check_parameters, nudge_init, &
	148	lsf_nudging_check_data_output_pr, lsf_nudging_header, &
[3182]	149	lsf_nesting_offline, lsf_nesting_offline_mass_conservation, &
	150	nudge, nudge_ref
	151
[2320]	152	!
	153	!-- Public variables
	154	PUBLIC qsws_surf, shf_surf, td_lsa_lpt, td_lsa_q, td_sub_lpt, &
[3182]	155	td_sub_q, time_vert
[2320]	156
[2696]	157
[2320]	158	INTERFACE ls_advec
	159	MODULE PROCEDURE ls_advec
	160	MODULE PROCEDURE ls_advec_ij
	161	END INTERFACE ls_advec
	162
	163	INTERFACE nudge
	164	MODULE PROCEDURE nudge
	165	MODULE PROCEDURE nudge_ij
	166	END INTERFACE nudge
	167
	168	CONTAINS
	169
	170
	171	!------------------------------------------------------------------------------!
	172	! Description:
	173	! ------------
[3182]	174	!> In this subroutine a constant mass within the model domain is guaranteed.
	175	!> Larger-scale models may be based on a compressible equation system, which is
	176	!> not consistent with PALMs incompressible equation system. In order to avoid
	177	!> a decrease or increase of mass during the simulation, non-divergent flow
	178	!> through the lateral and top boundaries is compensated by the vertical wind
	179	!> component at the top boundary.
[2320]	180	!------------------------------------------------------------------------------!
[3182]	181	SUBROUTINE lsf_nesting_offline_mass_conservation
[2696]	182
	183	USE control_parameters, &
	184	ONLY: volume_flow
	185
	186	IMPLICIT NONE
	187
[3182]	188	INTEGER(iwp) :: i !< grid index in x-direction
	189	INTEGER(iwp) :: j !< grid index in y-direction
	190	INTEGER(iwp) :: k !< grid index in z-direction
[2696]	191
[3182]	192	REAL(wp) :: w_correct !< vertical velocity increment required to compensate non-divergent flow through the boundaries
	193	REAL(wp), DIMENSION(1:3) :: volume_flow_l !< local volume flow
[2696]	194
	195	volume_flow = 0.0_wp
	196	volume_flow_l = 0.0_wp
	197
	198	d_area_t = 1.0_wp / ( ( nx + 1 ) * dx * ( ny + 1 ) * dy )
	199
[3182]	200	IF ( bc_dirichlet_l ) THEN
[2696]	201	i = nxl
	202	DO j = nys, nyn
	203	DO k = nzb+1, nzt
	204	volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy &
	205	* MERGE( 1.0_wp, 0.0_wp, &
	206	BTEST( wall_flags_0(k,j,i), 1 ) )
	207	ENDDO
	208	ENDDO
	209	ENDIF
[3182]	210	IF ( bc_dirichlet_r ) THEN
[2696]	211	i = nxr+1
	212	DO j = nys, nyn
	213	DO k = nzb+1, nzt
	214	volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy &
	215	* MERGE( 1.0_wp, 0.0_wp, &
	216	BTEST( wall_flags_0(k,j,i), 1 ) )
	217	ENDDO
	218	ENDDO
	219	ENDIF
[3182]	220	IF ( bc_dirichlet_s ) THEN
[2696]	221	j = nys
	222	DO i = nxl, nxr
	223	DO k = nzb+1, nzt
	224	volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx &
	225	* MERGE( 1.0_wp, 0.0_wp, &
	226	BTEST( wall_flags_0(k,j,i), 2 ) )
	227	ENDDO
	228	ENDDO
	229	ENDIF
[3182]	230	IF ( bc_dirichlet_n ) THEN
[2696]	231	j = nyn+1
	232	DO i = nxl, nxr
	233	DO k = nzb+1, nzt
	234	volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx &
	235	* MERGE( 1.0_wp, 0.0_wp, &
	236	BTEST( wall_flags_0(k,j,i), 2 ) )
	237	ENDDO
	238	ENDDO
	239	ENDIF
	240	!
[2938]	241	!-- Top boundary
[2696]	242	k = nzt
	243	DO i = nxl, nxr
	244	DO j = nys, nyn
	245	volume_flow_l(3) = volume_flow_l(3) - w(k,j,i) * dx * dy
	246	ENDDO
	247	ENDDO
	248
	249	#if defined( __parallel )
	250	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	251	CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, 3, MPI_REAL, MPI_SUM, &
	252	comm2d, ierr )
	253	#else
	254	volume_flow = volume_flow_l
	255	#endif
	256
	257	w_correct = SUM( volume_flow ) * d_area_t
	258
	259	DO i = nxl, nxr
	260	DO j = nys, nyn
	261	DO k = nzt, nzt + 1
	262	w(k,j,i) = w(k,j,i) + w_correct
	263	ENDDO
	264	ENDDO
	265	ENDDO
	266
[3182]	267	END SUBROUTINE lsf_nesting_offline_mass_conservation
[2938]	268
[2696]	269
	270	!------------------------------------------------------------------------------!
	271	! Description:
	272	! ------------
[3182]	273	!> Set the lateral and top boundary conditions in case the PALM domain is
	274	!> nested offline in a mesoscale model.
[2696]	275	!------------------------------------------------------------------------------!
[3182]	276	SUBROUTINE lsf_nesting_offline
[2696]	277
	278	USE control_parameters, &
[3182]	279	ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, &
	280	bc_dirichlet_s, humidity, neutral, passive_scalar, rans_mode,&
	281	rans_tke_e, time_since_reference_point
[2696]	282
	283	IMPLICIT NONE
	284
	285	INTEGER(iwp) :: i !< running index x-direction
	286	INTEGER(iwp) :: j !< running index y-direction
	287	INTEGER(iwp) :: k !< running index z-direction
	288	INTEGER(iwp) :: t !< running index for time levels
	289
	290	REAL(wp) :: ddt_lsf !< inverse value of time resolution of forcing data
	291	REAL(wp) :: t_ref !< time past since last reference step
[3182]	292
[2696]	293	!
	294	!-- Calculate time interval of forcing data
[3182]	295	ddt_lsf = 1.0_wp / ( nest_offl%time(nest_offl%tind_p) - &
	296	nest_offl%time(nest_offl%tind) )
[2696]	297	!
	298	!-- Calculate reziproke time past since last reference step. Please note,
[3182]	299	!-- the time coordinate is still not updated, so that the actual time need
	300	!-- to be incremented by dt_3d. Moreover, note that the simulation time
	301	!-- passed since simulation start is time_since_reference_point, not
	302	!-- simulated_time!
	303	t_ref = time_since_reference_point + dt_3d - &
	304	nest_offl%time(nest_offl%tind)
	305
	306	IF ( bc_dirichlet_l ) THEN
[2696]	307
	308	DO j = nys, nyn
[3182]	309	DO k = nzb+1, nzt
	310	u(k,j,nxlg:nxl) = nest_offl%u_left(0,k,j) + ddt_lsf * t_ref * &
	311	( nest_offl%u_left(1,k,j) - nest_offl%u_left(0,k,j) ) * &
	312	MERGE( 1.0_wp, 0.0_wp, &
	313	BTEST( wall_flags_0(k,j,nxlg:nxl), 1 ) )
[2696]	314	ENDDO
	315	ENDDO
	316
	317	DO j = nys, nyn
[3182]	318	DO k = nzb+1, nzt-1
	319	w(k,j,nxlg:nxl-1) = nest_offl%w_left(0,k,j) + ddt_lsf * t_ref *&
	320	( nest_offl%w_left(1,k,j) - nest_offl%w_left(0,k,j) ) *&
	321	MERGE( 1.0_wp, 0.0_wp, &
	322	BTEST( wall_flags_0(k,j,nxlg:nxl-1), 3 ) )
[2696]	323	ENDDO
	324	ENDDO
	325
	326	DO j = nysv, nyn
[3182]	327	DO k = nzb+1, nzt
	328	v(k,j,nxlg:nxl-1) = nest_offl%v_left(0,k,j) + ddt_lsf * t_ref *&
	329	( nest_offl%v_left(1,k,j) - nest_offl%v_left(0,k,j) ) *&
	330	MERGE( 1.0_wp, 0.0_wp, &
	331	BTEST( wall_flags_0(k,j,nxlg:nxl-1), 2 ) )
[2696]	332	ENDDO
	333	ENDDO
	334
	335	IF ( .NOT. neutral ) THEN
	336	DO j = nys, nyn
[3182]	337	DO k = nzb+1, nzt
	338	pt(k,j,nxlg:nxl-1) = nest_offl%pt_left(0,k,j) + ddt_lsf * &
	339	t_ref * &
	340	( nest_offl%pt_left(1,k,j) - nest_offl%pt_left(0,k,j) )
[2696]	341
	342	ENDDO
	343	ENDDO
	344	ENDIF
	345
	346	IF ( humidity ) THEN
	347	DO j = nys, nyn
[3182]	348	DO k = nzb+1, nzt
	349	q(k,j,nxlg:nxl-1) = nest_offl%q_left(0,k,j) + ddt_lsf * &
	350	t_ref * &
	351	( nest_offl%q_left(1,k,j) - nest_offl%q_left(0,k,j) )
[2696]	352
	353	ENDDO
	354	ENDDO
	355	ENDIF
	356
	357	ENDIF
	358
[3182]	359	IF ( bc_dirichlet_r ) THEN
[2696]	360
	361	DO j = nys, nyn
[3182]	362	DO k = nzb+1, nzt
	363	u(k,j,nxr+1:nxrg) = nest_offl%u_right(0,k,j) + ddt_lsf * t_ref *&
	364	( nest_offl%u_right(1,k,j) - nest_offl%u_right(0,k,j) ) *&
	365	MERGE( 1.0_wp, 0.0_wp, &
	366	BTEST( wall_flags_0(k,j,nxr+1:nxrg), 1 ) )
[2696]	367
	368	ENDDO
	369	ENDDO
	370	DO j = nys, nyn
[3182]	371	DO k = nzb+1, nzt-1
	372	w(k,j,nxr+1:nxrg) = nest_offl%w_right(0,k,j) + ddt_lsf * t_ref *&
	373	( nest_offl%w_right(1,k,j) - nest_offl%w_right(0,k,j) ) *&
	374	MERGE( 1.0_wp, 0.0_wp, &
	375	BTEST( wall_flags_0(k,j,nxr+1:nxrg), 3 ) )
[2696]	376	ENDDO
	377	ENDDO
	378
	379	DO j = nysv, nyn
[3182]	380	DO k = nzb+1, nzt
	381	v(k,j,nxr+1:nxrg) = nest_offl%v_right(0,k,j) + ddt_lsf * t_ref *&
	382	( nest_offl%v_right(1,k,j) - nest_offl%v_right(0,k,j) ) *&
	383	MERGE( 1.0_wp, 0.0_wp, &
	384	BTEST( wall_flags_0(k,j,nxr+1:nxrg), 2 ) )
[2696]	385	ENDDO
	386	ENDDO
	387
	388	IF ( .NOT. neutral ) THEN
	389	DO j = nys, nyn
[3182]	390	DO k = nzb+1, nzt
	391	pt(k,j,nxr+1:nxrg) = nest_offl%pt_right(0,k,j) + ddt_lsf * &
	392	t_ref * &
	393	( nest_offl%pt_right(1,k,j) - nest_offl%pt_right(0,k,j) )
[2696]	394
	395	ENDDO
	396	ENDDO
	397	ENDIF
	398
	399	IF ( humidity ) THEN
	400	DO j = nys, nyn
[3182]	401	DO k = nzb+1, nzt
	402	q(k,j,nxr+1:nxrg) = nest_offl%q_right(0,k,j) + ddt_lsf * &
	403	t_ref * &
	404	( nest_offl%q_right(1,k,j) - nest_offl%q_right(0,k,j) )
[2696]	405
	406	ENDDO
	407	ENDDO
	408	ENDIF
	409
	410	ENDIF
	411
[3182]	412	IF ( bc_dirichlet_s ) THEN
[2696]	413
	414	DO i = nxl, nxr
[3182]	415	DO k = nzb+1, nzt
	416	v(k,nysg:nys,i) = nest_offl%v_south(0,k,i) + ddt_lsf * t_ref *&
	417	( nest_offl%v_south(1,k,i) - nest_offl%v_south(0,k,i) ) *&
	418	MERGE( 1.0_wp, 0.0_wp, &
	419	BTEST( wall_flags_0(k,nysg:nys,i), 2 ) )
[2696]	420	ENDDO
	421	ENDDO
	422
	423	DO i = nxl, nxr
[3182]	424	DO k = nzb+1, nzt-1
	425	w(k,nysg:nys-1,i) = nest_offl%w_south(0,k,i) + ddt_lsf * t_ref *&
	426	( nest_offl%w_south(1,k,i) - nest_offl%w_south(0,k,i) ) *&
	427	MERGE( 1.0_wp, 0.0_wp, &
[2696]	428	BTEST( wall_flags_0(k,nysg:nys-1,i), 3 ) )
	429	ENDDO
	430	ENDDO
	431
	432	DO i = nxlu, nxr
[3182]	433	DO k = nzb+1, nzt
	434	u(k,nysg:nys-1,i) = nest_offl%u_south(0,k,i) + ddt_lsf * t_ref *&
	435	( nest_offl%u_south(1,k,i) - nest_offl%u_south(0,k,i) ) *&
	436	MERGE( 1.0_wp, 0.0_wp, &
[2696]	437	BTEST( wall_flags_0(k,nysg:nys-1,i), 1 ) )
	438	ENDDO
	439	ENDDO
	440
	441	IF ( .NOT. neutral ) THEN
	442	DO i = nxl, nxr
[3182]	443	DO k = nzb+1, nzt
	444	pt(k,nysg:nys-1,i) = nest_offl%pt_south(0,k,i) + ddt_lsf * &
	445	t_ref * &
	446	( nest_offl%pt_south(1,k,i) - nest_offl%pt_south(0,k,i) )
[2696]	447
	448	ENDDO
	449	ENDDO
	450	ENDIF
	451
	452	IF ( humidity ) THEN
	453	DO i = nxl, nxr
[3182]	454	DO k = nzb+1, nzt
	455	q(k,nysg:nys-1,i) = nest_offl%q_south(0,k,i) + ddt_lsf * &
	456	t_ref * &
	457	( nest_offl%q_south(1,k,i) - nest_offl%q_south(0,k,i) )
[2696]	458
	459	ENDDO
	460	ENDDO
	461	ENDIF
	462
	463	ENDIF
	464
[3182]	465	IF ( bc_dirichlet_n ) THEN
[2696]	466
	467	DO i = nxl, nxr
[3182]	468	DO k = nzb+1, nzt
	469	v(k,nyn+1:nyng,i) = nest_offl%v_north(0,k,i) + ddt_lsf * t_ref *&
	470	( nest_offl%v_north(1,k,i) - nest_offl%v_north(0,k,i) ) *&
	471	MERGE( 1.0_wp, 0.0_wp, &
[2696]	472	BTEST( wall_flags_0(k,nyn+1:nyng,i), 2 ) )
	473	ENDDO
	474	ENDDO
	475	DO i = nxl, nxr
[3182]	476	DO k = nzb+1, nzt-1
	477	w(k,nyn+1:nyng,i) = nest_offl%w_north(0,k,i) + ddt_lsf * t_ref *&
	478	( nest_offl%w_north(1,k,i) - nest_offl%w_north(0,k,i) ) *&
	479	MERGE( 1.0_wp, 0.0_wp, &
[2696]	480	BTEST( wall_flags_0(k,nyn+1:nyng,i), 3 ) )
	481	ENDDO
	482	ENDDO
	483
	484	DO i = nxlu, nxr
[3182]	485	DO k = nzb+1, nzt
	486	u(k,nyn+1:nyng,i) = nest_offl%u_north(0,k,i) + ddt_lsf * t_ref *&
	487	( nest_offl%u_north(1,k,i) - nest_offl%u_north(0,k,i) ) *&
	488	MERGE( 1.0_wp, 0.0_wp, &
[2696]	489	BTEST( wall_flags_0(k,nyn+1:nyng,i), 1 ) )
	490
	491	ENDDO
	492	ENDDO
	493
	494	IF ( .NOT. neutral ) THEN
	495	DO i = nxl, nxr
[3182]	496	DO k = nzb+1, nzt
	497	pt(k,nyn+1:nyng,i) = nest_offl%pt_north(0,k,i) + ddt_lsf * &
	498	t_ref * &
	499	( nest_offl%pt_north(1,k,i) - nest_offl%pt_north(0,k,i) )
[2696]	500
	501	ENDDO
	502	ENDDO
	503	ENDIF
	504
	505	IF ( humidity ) THEN
	506	DO i = nxl, nxr
[3182]	507	DO k = nzb+1, nzt
	508	q(k,nyn+1:nyng,i) = nest_offl%q_north(0,k,i) + ddt_lsf * &
	509	t_ref * &
	510	( nest_offl%q_north(1,k,i) - nest_offl%q_north(0,k,i) )
[2696]	511
	512	ENDDO
	513	ENDDO
	514	ENDIF
	515
	516	ENDIF
	517	!
	518	!-- Top boundary.
[2938]	519	DO i = nxlu, nxr
	520	DO j = nys, nyn
[3182]	521	u(nzt+1,j,i) = nest_offl%u_top(0,j,i) + ddt_lsf * t_ref * &
	522	( nest_offl%u_top(1,j,i) - nest_offl%u_top(0,j,i) ) * &
	523	MERGE( 1.0_wp, 0.0_wp, &
[2938]	524	BTEST( wall_flags_0(nzt+1,j,i), 1 ) )
[2696]	525	ENDDO
[2938]	526	ENDDO
[2696]	527
[2938]	528	DO i = nxl, nxr
	529	DO j = nysv, nyn
[3182]	530	v(nzt+1,j,i) = nest_offl%v_top(0,j,i) + ddt_lsf * t_ref * &
	531	( nest_offl%v_top(1,j,i) - nest_offl%v_top(0,j,i) ) * &
[2938]	532	MERGE( 1.0_wp, 0.0_wp, &
	533	BTEST( wall_flags_0(nzt+1,j,i), 2 ) )
	534	ENDDO
	535	ENDDO
	536
	537	DO i = nxl, nxr
	538	DO j = nys, nyn
[3182]	539	w(nzt:nzt+1,j,i) = nest_offl%w_top(0,j,i) + ddt_lsf * t_ref * &
	540	( nest_offl%w_top(1,j,i) - nest_offl%w_top(0,j,i) ) * &
[2938]	541	MERGE( 1.0_wp, 0.0_wp, &
	542	BTEST( wall_flags_0(nzt:nzt+1,j,i), 3 ) )
	543	ENDDO
	544	ENDDO
	545
	546
	547	IF ( .NOT. neutral ) THEN
[2696]	548	DO i = nxl, nxr
[2938]	549	DO j = nys, nyn
[3182]	550	pt(nzt+1,j,i) = nest_offl%pt_top(0,j,i) + ddt_lsf * t_ref * &
	551	( nest_offl%pt_top(1,j,i) - nest_offl%pt_top(0,j,i) )
[2696]	552	ENDDO
	553	ENDDO
[2938]	554	ENDIF
[2696]	555
[2938]	556	IF ( humidity ) THEN
[2696]	557	DO i = nxl, nxr
	558	DO j = nys, nyn
[3182]	559	q(nzt+1,j,i) = nest_offl%q_top(0,j,i) + ddt_lsf * t_ref * &
	560	( nest_offl%q_top(1,j,i) - nest_offl%q_top(0,j,i) )
[2696]	561	ENDDO
	562	ENDDO
	563	ENDIF
	564	!
[2938]	565	!-- At the edges( left-south, left-north, right-south and right-north) set
	566	!-- data on ghost points.
[3182]	567	IF ( bc_dirichlet_l .AND. bc_dirichlet_s ) THEN
[2938]	568	DO i = 1, nbgp
	569	u(:,nys-i,nxlg:nxl) = u(:,nys,nxlg:nxl)
	570	w(:,nys-i,nxlg:nxl-1) = w(:,nys,nxlg:nxl-1)
	571	IF ( .NOT. neutral ) pt(:,nys-i,nxlg:nxl-1) = pt(:,nys,nxlg:nxl-1)
[3182]	572	IF ( humidity ) q(:,nys-i,nxlg:nxl-1) = q(:,nys,nxlg:nxl-1)
[2938]	573	ENDDO
	574	DO i = 1, nbgp+1
	575	v(:,nysv-i,nxlg:nxl-1) = v(:,nysv,nxlg:nxl-1)
	576	ENDDO
	577	ENDIF
[3182]	578	IF ( bc_dirichlet_l .AND. bc_dirichlet_n ) THEN
[2938]	579	DO i = 1, nbgp
	580	u(:,nyn+i,nxlg:nxl) = u(:,nyn,nxlg:nxl)
	581	v(:,nyn+i,nxlg:nxl-1) = v(:,nyn,nxlg:nxl-1)
	582	w(:,nyn+i,nxlg:nxl-1) = w(:,nyn,nxlg:nxl-1)
	583	IF ( .NOT. neutral ) pt(:,nyn+i,nxlg:nxl-1) = pt(:,nyn,nxlg:nxl-1)
[3182]	584	IF ( humidity ) q(:,nyn+i,nxlg:nxl-1) = q(:,nyn,nxlg:nxl-1)
[2938]	585	ENDDO
	586	ENDIF
[3182]	587	IF ( bc_dirichlet_r .AND. bc_dirichlet_s ) THEN
[2938]	588	DO i = 1, nbgp
	589	u(:,nys-i,nxr+1:nxrg) = u(:,nys,nxr+1:nxrg)
	590	w(:,nys-i,nxr+1:nxrg) = w(:,nys,nxr+1:nxrg)
	591	IF ( .NOT. neutral ) pt(:,nys-i,nxr+1:nxrg) = pt(:,nys,nxr+1:nxrg)
[3182]	592	IF ( humidity ) q(:,nys-i,nxr+1:nxrg) = q(:,nys,nxr+1:nxrg)
[2938]	593	ENDDO
	594	DO i = 1, nbgp+1
	595	v(:,nysv-i,nxr+1:nxrg) = v(:,nysv,nxr+1:nxrg)
	596	ENDDO
	597	ENDIF
[3182]	598	IF ( bc_dirichlet_r .AND. bc_dirichlet_n ) THEN
[2938]	599	DO i = 1, nbgp
	600	u(:,nyn+i,nxr+1:nxrg) = u(:,nyn,nxr+1:nxrg)
	601	v(:,nyn+i,nxr+1:nxrg) = v(:,nyn,nxr+1:nxrg)
	602	w(:,nyn+i,nxr+1:nxrg) = w(:,nyn,nxr+1:nxrg)
	603	IF ( .NOT. neutral ) pt(:,nyn+i,nxr+1:nxrg) = pt(:,nyn,nxr+1:nxrg)
[3182]	604	IF ( humidity ) q(:,nyn+i,nxr+1:nxrg) = q(:,nyn,nxr+1:nxrg)
[2938]	605	ENDDO
	606	ENDIF
	607	!
[3182]	608	!-- Moreover, set Neumann boundary condition for subgrid-scale TKE,
	609	!-- passive scalar, dissipation, and chemical species if required
	610	IF ( rans_mode .AND. rans_tke_e ) THEN
	611	IF ( bc_dirichlet_l ) diss(:,:,nxl-1) = diss(:,:,nxl)
	612	IF ( bc_dirichlet_r ) diss(:,:,nxr+1) = diss(:,:,nxr)
	613	IF ( bc_dirichlet_s ) diss(:,nys-1,:) = diss(:,nys,:)
	614	IF ( bc_dirichlet_n ) diss(:,nyn+1,:) = diss(:,nyn,:)
	615	ENDIF
[2696]	616	IF ( .NOT. constant_diffusion ) THEN
[3182]	617	IF ( bc_dirichlet_l ) e(:,:,nxl-1) = e(:,:,nxl)
	618	IF ( bc_dirichlet_r ) e(:,:,nxr+1) = e(:,:,nxr)
	619	IF ( bc_dirichlet_s ) e(:,nys-1,:) = e(:,nys,:)
	620	IF ( bc_dirichlet_n ) e(:,nyn+1,:) = e(:,nyn,:)
[2696]	621	e(nzt+1,:,:) = e(nzt,:,:)
	622	ENDIF
	623	IF ( passive_scalar ) THEN
[3182]	624	IF ( bc_dirichlet_l ) s(:,:,nxl-1) = s(:,:,nxl)
	625	IF ( bc_dirichlet_r ) s(:,:,nxr+1) = s(:,:,nxr)
	626	IF ( bc_dirichlet_s ) s(:,nys-1,:) = s(:,nys,:)
	627	IF ( bc_dirichlet_n ) s(:,nyn+1,:) = s(:,nyn,:)
[2696]	628	ENDIF
	629
	630
[3182]	631
[2696]	632	CALL exchange_horiz( u, nbgp )
	633	CALL exchange_horiz( v, nbgp )
	634	CALL exchange_horiz( w, nbgp )
	635	IF ( .NOT. neutral ) CALL exchange_horiz( pt, nbgp )
	636	IF ( humidity ) CALL exchange_horiz( q, nbgp )
	637
	638	!
	639	!-- Set surface pressure. Please note, time-dependent surface
	640	!-- pressure would require changes in anelastic approximation and
	641	!-- treatment of fluxes.
	642	!-- For the moment, comment this out!
[3182]	643	! surface_pressure = nest_offl%surface_pressure(nest_offl%tind) + &
[2696]	644	! ddt_lsf * t_ref * &
[3182]	645	! ( nest_offl%surface_pressure(nest_offl%tind_p) &
	646	! - nest_offl%surface_pressure(nest_offl%tind) )
[2696]	647
[3182]	648	END SUBROUTINE lsf_nesting_offline
[2696]	649
	650	!------------------------------------------------------------------------------!
	651	! Description:
	652	! ------------
	653	!> @todo Missing subroutine description.
	654	!------------------------------------------------------------------------------!
[2320]	655	SUBROUTINE lsf_nudging_check_parameters
	656
	657	IMPLICIT NONE
	658	!
	659	!-- Check nudging and large scale forcing from external file
	660	IF ( nudging .AND. ( .NOT. large_scale_forcing ) ) THEN
[3046]	661	message_string = 'Nudging requires large_scale_forcing = .T.. &'// &
	662	'Surface fluxes and geostrophic wind should be &'// &
[2320]	663	'prescribed in file LSF_DATA'
	664	CALL message( 'check_parameters', 'PA0374', 1, 2, 0, 6, 0 )
	665	ENDIF
	666
	667	IF ( large_scale_forcing .AND. ( bc_lr /= 'cyclic' .OR. &
[2970]	668	bc_ns /= 'cyclic' ) ) THEN
[3046]	669	message_string = 'Non-cyclic lateral boundaries do not allow for &'//&
[2320]	670	'the usage of large scale forcing from external file.'
	671	CALL message( 'check_parameters', 'PA0375', 1, 2, 0, 6, 0 )
	672	ENDIF
	673
	674	IF ( large_scale_forcing .AND. ( .NOT. humidity ) ) THEN
[3046]	675	message_string = 'The usage of large scale forcing from external &'//&
[2320]	676	'file LSF_DATA requires humidity = .T..'
	677	CALL message( 'check_parameters', 'PA0376', 1, 2, 0, 6, 0 )
	678	ENDIF
	679
	680	IF ( large_scale_forcing .AND. passive_scalar ) THEN
[3046]	681	message_string = 'The usage of large scale forcing from external &'// &
[2320]	682	'file LSF_DATA is not implemented for passive scalars'
	683	CALL message( 'check_parameters', 'PA0440', 1, 2, 0, 6, 0 )
	684	ENDIF
	685
	686	IF ( large_scale_forcing .AND. topography /= 'flat' &
	687	.AND. .NOT. lsf_exception ) THEN
[3046]	688	message_string = 'The usage of large scale forcing from external &'//&
[2320]	689	'file LSF_DATA is not implemented for non-flat topography'
	690	CALL message( 'check_parameters', 'PA0377', 1, 2, 0, 6, 0 )
	691	ENDIF
	692
	693	IF ( large_scale_forcing .AND. ocean ) THEN
[3046]	694	message_string = 'The usage of large scale forcing from external &'//&
[2320]	695	'file LSF_DATA is not implemented for ocean runs'
	696	CALL message( 'check_parameters', 'PA0378', 1, 2, 0, 6, 0 )
	697	ENDIF
	698
	699	END SUBROUTINE lsf_nudging_check_parameters
	700
	701	!------------------------------------------------------------------------------!
	702	! Description:
	703	! ------------
	704	!> Check data output of profiles for land surface model
	705	!------------------------------------------------------------------------------!
	706	SUBROUTINE lsf_nudging_check_data_output_pr( variable, var_count, unit, &
	707	dopr_unit )
	708
	709	USE profil_parameter
	710
	711	IMPLICIT NONE
	712
	713	CHARACTER (LEN=*) :: unit !<
	714	CHARACTER (LEN=*) :: variable !<
	715	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
	716
	717	INTEGER(iwp) :: user_pr_index !<
	718	INTEGER(iwp) :: var_count !<
	719
	720	SELECT CASE ( TRIM( variable ) )
	721
	722
	723	CASE ( 'td_lsa_lpt' )
	724	IF ( .NOT. large_scale_forcing ) THEN
	725	message_string = 'data_output_pr = ' // &
	726	TRIM( data_output_pr(var_count) ) // &
	727	' is not implemented for ' // &
	728	'large_scale_forcing = .FALSE.'
	729	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
	730	1, 2, 0, 6, 0 )
	731	ELSE
	732	dopr_index(var_count) = 81
	733	dopr_unit = 'K/s'
	734	unit = 'K/s'
	735	hom(:,2,81,:) = SPREAD( zu, 2, statistic_regions+1 )
	736	ENDIF
	737
	738	CASE ( 'td_lsa_q' )
	739	IF ( .NOT. large_scale_forcing ) THEN
	740	message_string = 'data_output_pr = ' // &
	741	TRIM( data_output_pr(var_count) ) // &
	742	' is not implemented for ' // &
	743	'large_scale_forcing = .FALSE.'
	744	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
	745	1, 2, 0, 6, 0 )
	746	ELSE
	747	dopr_index(var_count) = 82
	748	dopr_unit = 'kg/kgs'
	749	unit = 'kg/kgs'
	750	hom(:,2,82,:) = SPREAD( zu, 2, statistic_regions+1 )
	751	ENDIF
	752	CASE ( 'td_sub_lpt' )
	753	IF ( .NOT. large_scale_forcing ) THEN
	754	message_string = 'data_output_pr = ' // &
	755	TRIM( data_output_pr(var_count) ) // &
	756	' is not implemented for ' // &
	757	'large_scale_forcing = .FALSE.'
	758	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
	759	1, 2, 0, 6, 0 )
	760	ELSE
	761	dopr_index(var_count) = 83
	762	dopr_unit = 'K/s'
	763	unit = 'K/s'
	764	hom(:,2,83,:) = SPREAD( zu, 2, statistic_regions+1 )
	765	ENDIF
	766
	767	CASE ( 'td_sub_q' )
	768	IF ( .NOT. large_scale_forcing ) THEN
	769	message_string = 'data_output_pr = ' // &
	770	TRIM( data_output_pr(var_count) ) // &
	771	' is not implemented for ' // &
	772	'large_scale_forcing = .FALSE.'
	773	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0393', &
	774	1, 2, 0, 6, 0 )
	775	ELSE
	776	dopr_index(var_count) = 84
	777	dopr_unit = 'kg/kgs'
	778	unit = 'kg/kgs'
	779	hom(:,2,84,:) = SPREAD( zu, 2, statistic_regions+1 )
	780	ENDIF
	781
	782	CASE ( 'td_nud_lpt' )
	783	IF ( .NOT. nudging ) THEN
	784	message_string = 'data_output_pr = ' // &
	785	TRIM( data_output_pr(var_count) ) // &
	786	' is not implemented for ' // &
[3045]	787	'nudging = .FALSE.'
	788	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394', &
[2320]	789	1, 2, 0, 6, 0 )
	790	ELSE
	791	dopr_index(var_count) = 85
	792	dopr_unit = 'K/s'
	793	unit = 'K/s'
	794	hom(:,2,85,:) = SPREAD( zu, 2, statistic_regions+1 )
	795	ENDIF
	796
	797	CASE ( 'td_nud_q' )
	798	IF ( .NOT. nudging ) THEN
	799	message_string = 'data_output_pr = ' // &
	800	TRIM( data_output_pr(var_count) ) // &
	801	' is not implemented for ' // &
[3045]	802	'nudging = .FALSE.'
	803	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394', &
[2320]	804	1, 2, 0, 6, 0 )
	805	ELSE
	806	dopr_index(var_count) = 86
	807	dopr_unit = 'kg/kgs'
	808	unit = 'kg/kgs'
	809	hom(:,2,86,:) = SPREAD( zu, 2, statistic_regions+1 )
	810	ENDIF
	811
	812	CASE ( 'td_nud_u' )
	813	IF ( .NOT. nudging ) THEN
	814	message_string = 'data_output_pr = ' // &
	815	TRIM( data_output_pr(var_count) ) // &
	816	' is not implemented for ' // &
[3045]	817	'nudging = .FALSE.'
	818	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394', &
[2320]	819	1, 2, 0, 6, 0 )
	820	ELSE
	821	dopr_index(var_count) = 87
	822	dopr_unit = 'm/s2'
	823	unit = 'm/s2'
	824	hom(:,2,87,:) = SPREAD( zu, 2, statistic_regions+1 )
	825	ENDIF
	826
	827	CASE ( 'td_nud_v' )
	828	IF ( .NOT. nudging ) THEN
	829	message_string = 'data_output_pr = ' // &
	830	TRIM( data_output_pr(var_count) ) // &
	831	' is not implemented for ' // &
[3045]	832	'nudging = .FALSE.'
	833	CALL message( 'lsf_nudging_check_data_output_pr', 'PA0394', &
[2320]	834	1, 2, 0, 6, 0 )
	835	ELSE
	836	dopr_index(var_count) = 88
	837	dopr_unit = 'm/s2'
	838	unit = 'm/s2'
	839	hom(:,2,88,:) = SPREAD( zu, 2, statistic_regions+1 )
	840	ENDIF
	841
	842
	843	CASE DEFAULT
	844	unit = 'illegal'
	845
	846	END SELECT
	847
	848	END SUBROUTINE lsf_nudging_check_data_output_pr
	849
	850	!------------------------------------------------------------------------------!
	851	! Description:
	852	! ------------
	853	!> @todo Missing subroutine description.
	854	!------------------------------------------------------------------------------!
	855	SUBROUTINE lsf_nudging_header ( io )
	856
	857	IMPLICIT NONE
	858
	859	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
	860
	861	WRITE ( io, 1 )
	862	IF ( large_scale_forcing ) THEN
	863	WRITE ( io, 3 )
	864	WRITE ( io, 4 )
	865
	866	IF ( large_scale_subsidence ) THEN
	867	IF ( .NOT. use_subsidence_tendencies ) THEN
	868	WRITE ( io, 5 )
	869	ELSE
	870	WRITE ( io, 6 )
	871	ENDIF
	872	ENDIF
	873
	874	IF ( bc_pt_b == 'dirichlet' ) THEN
	875	WRITE ( io, 12 )
	876	ELSEIF ( bc_pt_b == 'neumann' ) THEN
	877	WRITE ( io, 13 )
	878	ENDIF
	879
	880	IF ( bc_q_b == 'dirichlet' ) THEN
	881	WRITE ( io, 14 )
	882	ELSEIF ( bc_q_b == 'neumann' ) THEN
	883	WRITE ( io, 15 )
	884	ENDIF
	885
	886	WRITE ( io, 7 )
	887	IF ( nudging ) THEN
	888	WRITE ( io, 10 )
	889	ENDIF
	890	ELSE
	891	WRITE ( io, 2 )
	892	WRITE ( io, 11 )
	893	ENDIF
	894	IF ( large_scale_subsidence ) THEN
	895	WRITE ( io, 8 )
	896	WRITE ( io, 9 )
	897	ENDIF
	898
	899
	900	1 FORMAT (//' Large scale forcing and nudging:'/ &
	901	' -------------------------------'/)
	902	2 FORMAT (' --> No large scale forcing from external is used (default) ')
	903	3 FORMAT (' --> Large scale forcing from external file LSF_DATA is used: ')
	904	4 FORMAT (' - large scale advection tendencies ')
	905	5 FORMAT (' - large scale subsidence velocity w_subs ')
	906	6 FORMAT (' - large scale subsidence tendencies ')
	907	7 FORMAT (' - and geostrophic wind components ug and vg')
	908	8 FORMAT (' --> Large-scale vertical motion is used in the ', &
	909	'prognostic equation(s) for')
	910	9 FORMAT (' the scalar(s) only')
	911	10 FORMAT (' --> Nudging is used')
	912	11 FORMAT (' --> No nudging is used (default) ')
	913	12 FORMAT (' - prescribed surface values for temperature')
	914	13 FORMAT (' - prescribed surface fluxes for temperature')
	915	14 FORMAT (' - prescribed surface values for humidity')
	916	15 FORMAT (' - prescribed surface fluxes for humidity')
	917
	918	END SUBROUTINE lsf_nudging_header
	919
	920	!------------------------------------------------------------------------------!
	921	! Description:
	922	! ------------
	923	!> @todo Missing subroutine description.
	924	!------------------------------------------------------------------------------!
	925	SUBROUTINE lsf_init
	926
[2938]	927	USE control_parameters, &
	928	ONLY: bc_lr_cyc, bc_ns_cyc
	929
[2696]	930	USE netcdf_data_input_mod, &
	931	ONLY: netcdf_data_input_lsf
	932
[2320]	933	IMPLICIT NONE
	934
	935	CHARACTER(100) :: chmess !<
	936	CHARACTER(1) :: hash !<
	937
	938	INTEGER(iwp) :: ierrn !<
	939	INTEGER(iwp) :: finput = 90 !<
	940	INTEGER(iwp) :: k !<
[2696]	941	INTEGER(iwp) :: nt !<
	942	INTEGER(iwp) :: t !< running index for time levels
[2320]	943
	944	REAL(wp) :: fac !<
	945	REAL(wp) :: highheight !<
	946	REAL(wp) :: highug_vert !<
	947	REAL(wp) :: highvg_vert !<
	948	REAL(wp) :: highwsubs_vert !<
	949	REAL(wp) :: lowheight !<
	950	REAL(wp) :: lowug_vert !<
	951	REAL(wp) :: lowvg_vert !<
	952	REAL(wp) :: lowwsubs_vert !<
	953	REAL(wp) :: high_td_lsa_lpt !<
	954	REAL(wp) :: low_td_lsa_lpt !<
	955	REAL(wp) :: high_td_lsa_q !<
	956	REAL(wp) :: low_td_lsa_q !<
	957	REAL(wp) :: high_td_sub_lpt !<
	958	REAL(wp) :: low_td_sub_lpt !<
	959	REAL(wp) :: high_td_sub_q !<
	960	REAL(wp) :: low_td_sub_q !<
	961	REAL(wp) :: r_dummy !<
	962
[3182]	963	IF ( nesting_offline ) THEN
[2696]	964	!
[2938]	965	!-- Allocate arrays for geostrophic wind components. Arrays will
	966	!-- incorporate 2 time levels in order to interpolate in between. Please
	967	!-- note, forcing using geostrophic wind components is only required in
	968	!-- case of cyclic boundary conditions.
	969	IF ( bc_lr_cyc .AND. bc_ns_cyc ) THEN
[3182]	970	ALLOCATE( nest_offl%ug(0:1,nzb:nzt+1) )
	971	ALLOCATE( nest_offl%vg(0:1,nzb:nzt+1) )
[2938]	972	ENDIF
	973	!
[2696]	974	!-- Allocate arrays for reading boundary values. Arrays will incorporate 2
	975	!-- time levels in order to interpolate in between.
[3182]	976	IF ( bc_dirichlet_l ) THEN
	977	ALLOCATE( nest_offl%u_left(0:1,nzb+1:nzt,nys:nyn) )
	978	ALLOCATE( nest_offl%v_left(0:1,nzb+1:nzt,nysv:nyn) )
	979	ALLOCATE( nest_offl%w_left(0:1,nzb+1:nzt-1,nys:nyn) )
	980	IF ( humidity ) ALLOCATE( nest_offl%q_left(0:1,nzb+1:nzt,nys:nyn) )
	981	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_left(0:1,nzb+1:nzt,nys:nyn) )
[2696]	982	ENDIF
[3182]	983	IF ( bc_dirichlet_r ) THEN
	984	ALLOCATE( nest_offl%u_right(0:1,nzb+1:nzt,nys:nyn) )
	985	ALLOCATE( nest_offl%v_right(0:1,nzb+1:nzt,nysv:nyn) )
	986	ALLOCATE( nest_offl%w_right(0:1,nzb+1:nzt-1,nys:nyn) )
	987	IF ( humidity ) ALLOCATE( nest_offl%q_right(0:1,nzb+1:nzt,nys:nyn) )
	988	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_right(0:1,nzb+1:nzt,nys:nyn) )
[2696]	989	ENDIF
[3182]	990	IF ( bc_dirichlet_n ) THEN
	991	ALLOCATE( nest_offl%u_north(0:1,nzb+1:nzt,nxlu:nxr) )
	992	ALLOCATE( nest_offl%v_north(0:1,nzb+1:nzt,nxl:nxr) )
	993	ALLOCATE( nest_offl%w_north(0:1,nzb+1:nzt-1,nxl:nxr) )
	994	IF ( humidity ) ALLOCATE( nest_offl%q_north(0:1,nzb+1:nzt,nxl:nxr) )
	995	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_north(0:1,nzb+1:nzt,nxl:nxr) )
[2696]	996	ENDIF
[3182]	997	IF ( bc_dirichlet_s ) THEN
	998	ALLOCATE( nest_offl%u_south(0:1,nzb+1:nzt,nxlu:nxr) )
	999	ALLOCATE( nest_offl%v_south(0:1,nzb+1:nzt,nxl:nxr) )
	1000	ALLOCATE( nest_offl%w_south(0:1,nzb+1:nzt-1,nxl:nxr) )
	1001	IF ( humidity ) ALLOCATE( nest_offl%q_south(0:1,nzb+1:nzt,nxl:nxr) )
	1002	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_south(0:1,nzb+1:nzt,nxl:nxr) )
[2696]	1003	ENDIF
	1004
[3182]	1005	ALLOCATE( nest_offl%u_top(0:1,nys:nyn,nxlu:nxr) )
	1006	ALLOCATE( nest_offl%v_top(0:1,nysv:nyn,nxl:nxr) )
	1007	ALLOCATE( nest_offl%w_top(0:1,nys:nyn,nxl:nxr) )
	1008	IF ( humidity ) ALLOCATE( nest_offl%q_top(0:1,nys:nyn,nxl:nxr) )
	1009	IF ( .NOT. neutral ) ALLOCATE( nest_offl%pt_top(0:1,nys:nyn,nxl:nxr) )
[2320]	1010
[2696]	1011	!
[3182]	1012	!-- Read COSMO data at lateral and top boundaries
[2696]	1013	CALL netcdf_data_input_lsf
	1014	!
[3182]	1015	!-- Write COSMO data at lateral and top boundaries
	1016	CALL lsf_nesting_offline
[2696]	1017	!
	1018	!-- After 3D data is initialized, ensure mass conservation
[3182]	1019	CALL lsf_nesting_offline_mass_conservation
[2696]	1020	!
	1021	!-- Initialize surface pressure. Please note, time-dependent surface
	1022	!-- pressure would require changes in anelastic approximation and
	1023	!-- treatment of fluxes.
	1024	!-- For the moment, comment this out!
[3182]	1025	! surface_pressure = nest_offl%surface_pressure(0)
[2320]	1026
[2696]	1027	ELSE
[2320]	1028
[2696]	1029	ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf), &
	1030	qsws_surf(0:nlsf), shf_surf(0:nlsf), &
	1031	td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf), &
	1032	td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf), &
	1033	time_vert(0:nlsf), time_surf(0:nlsf), &
	1034	ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf), &
	1035	wsubs_vert(nzb:nzt+1,0:nlsf) )
[2320]	1036
[2696]	1037	p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
	1038	shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
	1039	td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
	1040	time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
	1041	wsubs_vert = 0.0_wp
	1042
[2320]	1043	!
[2696]	1044	!-- Array for storing large scale forcing and nudging tendencies at each
	1045	!-- timestep for data output
	1046	ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
	1047	sums_ls_l = 0.0_wp
[2320]	1048
[2696]	1049	ngp_sums_ls = (nz+2)*6
[2320]	1050
[2696]	1051	OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
	1052	FORM='FORMATTED', IOSTAT=ierrn )
	1053
	1054	IF ( ierrn /= 0 ) THEN
	1055	message_string = 'file LSF_DATA does not exist'
	1056	CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
	1057	ENDIF
	1058
	1059	ierrn = 0
[2320]	1060	!
[2696]	1061	!-- First three lines of LSF_DATA contain header
	1062	READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
	1063	READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
	1064	READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
[2320]	1065
[2342]	1066	IF ( ierrn /= 0 ) THEN
[2696]	1067	message_string = 'errors in file LSF_DATA'
	1068	CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
	1069	ENDIF
[2320]	1070
[2696]	1071	!
	1072	!-- Surface values are read in
	1073	nt = 0
	1074	ierrn = 0
	1075
	1076	DO WHILE ( time_surf(nt) < end_time )
	1077	nt = nt + 1
	1078	READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt), &
	1079	qsws_surf(nt), pt_surf(nt), &
	1080	q_surf(nt), p_surf(nt)
	1081
	1082	IF ( ierrn /= 0 ) THEN
[3045]	1083	WRITE ( message_string, * ) 'No time dependent surface ' // &
[3046]	1084	'variables in & LSF_DATA for end of run found'
[2696]	1085
	1086	CALL message( 'ls_forcing', 'PA0363', 1, 2, 0, 6, 0 )
	1087	ENDIF
	1088	ENDDO
	1089
	1090	IF ( time_surf(1) > end_time ) THEN
	1091	WRITE ( message_string, * ) 'Time dependent surface variables in ' // &
[3046]	1092	'&LSF_DATA set in after end of ' , &
[2696]	1093	'simulation - lsf_surf is set to FALSE'
	1094	CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
	1095	lsf_surf = .FALSE.
[2320]	1096	ENDIF
	1097
[2696]	1098	!
	1099	!-- Go to the end of the list with surface variables
	1100	DO WHILE ( ierrn == 0 )
	1101	READ ( finput, *, IOSTAT = ierrn ) r_dummy
	1102	ENDDO
[2320]	1103
	1104	!
[2696]	1105	!-- Profiles of ug, vg and w_subs are read in (large scale forcing)
[2320]	1106
[2696]	1107	nt = 0
	1108	DO WHILE ( time_vert(nt) < end_time )
	1109	nt = nt + 1
	1110	hash = "#"
	1111	ierrn = 1 ! not zero
[2320]	1112	!
[2696]	1113	!-- Search for the next line consisting of "# time",
	1114	!-- from there onwards the profiles will be read
	1115	DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
	1116	READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
	1117	IF ( ierrn < 0 ) THEN
	1118	WRITE( message_string, * ) 'No time dependent vertical profiles',&
[3046]	1119	' in & LSF_DATA for end of run found'
[2696]	1120	CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
	1121	ENDIF
	1122	ENDDO
[2320]	1123
[2696]	1124	IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
	1125
	1126	READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert,&
	1127	lowwsubs_vert, low_td_lsa_lpt, &
	1128	low_td_lsa_q, low_td_sub_lpt, &
	1129	low_td_sub_q
	1130	IF ( ierrn /= 0 ) THEN
	1131	message_string = 'errors in file LSF_DATA'
	1132	CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
[2320]	1133	ENDIF
	1134
[2696]	1135	READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, &
	1136	highvg_vert, highwsubs_vert, &
	1137	high_td_lsa_lpt, high_td_lsa_q, &
	1138	high_td_sub_lpt, high_td_sub_q
	1139
	1140	IF ( ierrn /= 0 ) THEN
	1141	message_string = 'errors in file LSF_DATA'
	1142	CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
	1143	ENDIF
[2320]	1144
	1145
[2696]	1146	DO k = nzb, nzt+1
	1147	IF ( highheight < zu(k) ) THEN
	1148	lowheight = highheight
	1149	lowug_vert = highug_vert
	1150	lowvg_vert = highvg_vert
	1151	lowwsubs_vert = highwsubs_vert
	1152	low_td_lsa_lpt = high_td_lsa_lpt
	1153	low_td_lsa_q = high_td_lsa_q
	1154	low_td_sub_lpt = high_td_sub_lpt
	1155	low_td_sub_q = high_td_sub_q
[2320]	1156
[2696]	1157	ierrn = 0
	1158	READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, &
	1159	highvg_vert, highwsubs_vert,&
	1160	high_td_lsa_lpt, &
	1161	high_td_lsa_q, &
	1162	high_td_sub_lpt, high_td_sub_q
[2320]	1163
[2696]	1164	IF ( ierrn /= 0 ) THEN
[3045]	1165	WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm ', &
	1166	'is higher than the maximum height in LSF_DATA ', &
	1167	'which is ', lowheight, 'm. Interpolation on PALM ',&
[2696]	1168	'grid is not possible.'
	1169	CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
	1170	ENDIF
[2320]	1171
	1172	ENDIF
	1173
	1174	!
[2696]	1175	!-- Interpolation of prescribed profiles in space
	1176	fac = (highheight-zu(k))/(highheight - lowheight)
[2320]	1177
[2696]	1178	ug_vert(k,nt) = fac * lowug_vert &
	1179	+ ( 1.0_wp - fac ) * highug_vert
	1180	vg_vert(k,nt) = fac * lowvg_vert &
	1181	+ ( 1.0_wp - fac ) * highvg_vert
	1182	wsubs_vert(k,nt) = fac * lowwsubs_vert &
	1183	+ ( 1.0_wp - fac ) * highwsubs_vert
[2320]	1184
[2696]	1185	td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt &
	1186	+ ( 1.0_wp - fac ) * high_td_lsa_lpt
	1187	td_lsa_q(k,nt) = fac * low_td_lsa_q &
	1188	+ ( 1.0_wp - fac ) * high_td_lsa_q
	1189	td_sub_lpt(k,nt) = fac * low_td_sub_lpt &
	1190	+ ( 1.0_wp - fac ) * high_td_sub_lpt
	1191	td_sub_q(k,nt) = fac * low_td_sub_q &
	1192	+ ( 1.0_wp - fac ) * high_td_sub_q
[2320]	1193
[2696]	1194	ENDDO
[2320]	1195
[2696]	1196	ENDDO
[2320]	1197
	1198	!
[2696]	1199	!-- Large scale vertical velocity has to be zero at the surface
	1200	wsubs_vert(nzb,:) = 0.0_wp
	1201
	1202	IF ( time_vert(1) > end_time ) THEN
	1203	WRITE ( message_string, * ) 'Time dependent large scale profile ',&
	1204	'forcing from&LSF_DATA sets in after end of ' ,&
	1205	'simulation - lsf_vert is set to FALSE'
	1206	CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
	1207	lsf_vert = .FALSE.
	1208	ENDIF
[2320]	1209
[2696]	1210	CLOSE( finput )
[2320]	1211
[2696]	1212	ENDIF
[2320]	1213
	1214	END SUBROUTINE lsf_init
	1215
	1216	!------------------------------------------------------------------------------!
	1217	! Description:
	1218	! ------------
	1219	!> @todo Missing subroutine description.
	1220	!------------------------------------------------------------------------------!
	1221	SUBROUTINE ls_forcing_surf ( time )
	1222
	1223	IMPLICIT NONE
	1224
	1225	INTEGER(iwp) :: nt !<
	1226
	1227	REAL(wp) :: dum_surf_flux !<
	1228	REAL(wp) :: fac !<
	1229	REAL(wp), INTENT(in) :: time !<
	1230
	1231	!
	1232	!-- Interpolation in time of LSF_DATA at the surface
	1233	nt = 1
	1234	DO WHILE ( time > time_surf(nt) )
	1235	nt = nt + 1
	1236	ENDDO
	1237	IF ( time /= time_surf(nt) ) THEN
	1238	nt = nt - 1
	1239	ENDIF
	1240
	1241	fac = ( time -time_surf(nt) ) / ( time_surf(nt+1) - time_surf(nt) )
	1242
	1243	IF ( ibc_pt_b == 0 ) THEN
	1244	!
	1245	!-- In case of Dirichlet boundary condition shf must not
	1246	!-- be set - it is calculated via MOST in prandtl_fluxes
	1247	pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
	1248
	1249	ELSEIF ( ibc_pt_b == 1 ) THEN
	1250	!
	1251	!-- In case of Neumann boundary condition pt_surface is needed for
	1252	!-- calculation of reference density
	1253	dum_surf_flux = ( shf_surf(nt) + fac * &
	1254	( shf_surf(nt+1) - shf_surf(nt) ) &
	1255	) * heatflux_input_conversion(nzb)
	1256	!
	1257	!-- Save surface sensible heat flux on default, natural and urban surface
	1258	!-- type, if required
	1259	IF ( surf_def_h(0)%ns >= 1 ) surf_def_h(0)%shf(:) = dum_surf_flux
	1260	IF ( surf_lsm_h%ns >= 1 ) surf_lsm_h%shf(:) = dum_surf_flux
	1261	IF ( surf_usm_h%ns >= 1 ) surf_usm_h%shf(:) = dum_surf_flux
	1262
	1263	pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
	1264
	1265	ENDIF
	1266
	1267	IF ( ibc_q_b == 0 ) THEN
	1268	!
	1269	!-- In case of Dirichlet boundary condition qsws must not
	1270	!-- be set - it is calculated via MOST in prandtl_fluxes
	1271	q_surface = q_surf(nt) + fac * ( q_surf(nt+1) - q_surf(nt) )
	1272
	1273	ELSEIF ( ibc_q_b == 1 ) THEN
	1274	dum_surf_flux = ( qsws_surf(nt) + fac * &
	1275	( qsws_surf(nt+1) - qsws_surf(nt) ) &
	1276	) * waterflux_input_conversion(nzb)
	1277	!
	1278	!-- Save surface sensible heat flux on default, natural and urban surface
	1279	!-- type, if required
	1280	IF ( surf_def_h(0)%ns >= 1 ) surf_def_h(0)%qsws(:) = dum_surf_flux
	1281	IF ( surf_lsm_h%ns >= 1 ) surf_lsm_h%qsws(:) = dum_surf_flux
	1282	IF ( surf_usm_h%ns >= 1 ) surf_usm_h%qsws(:) = dum_surf_flux
	1283
	1284	ENDIF
[2970]	1285	!
	1286	!-- Surface heat- and waterflux will be written later onto surface elements
	1287	IF ( .NOT. neutral .AND. constant_heatflux .AND. &
	1288	TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
	1289	surface_heatflux = shf_surf(1)
	1290	ENDIF
	1291	IF ( humidity .AND. constant_waterflux .AND. &
	1292	TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
	1293	surface_waterflux = qsws_surf(1)
	1294	ENDIF
[2320]	1295
	1296	surface_pressure = p_surf(nt) + fac * ( p_surf(nt+1) - p_surf(nt) )
	1297
	1298	END SUBROUTINE ls_forcing_surf
	1299
	1300
	1301
	1302
	1303	!------------------------------------------------------------------------------!
	1304	! Description:
	1305	! ------------
	1306	!> @todo Missing subroutine description.
	1307	!------------------------------------------------------------------------------!
	1308	SUBROUTINE ls_forcing_vert ( time )
	1309
	1310
	1311	IMPLICIT NONE
	1312
	1313	INTEGER(iwp) :: nt !<
	1314
	1315	REAL(wp) :: fac !<
	1316	REAL(wp), INTENT(in) :: time !<
	1317
	1318	!
	1319	!-- Interpolation in time of LSF_DATA for ug, vg and w_subs
	1320	nt = 1
	1321	DO WHILE ( time > time_vert(nt) )
	1322	nt = nt + 1
	1323	ENDDO
	1324	IF ( time /= time_vert(nt) ) THEN
	1325	nt = nt - 1
	1326	ENDIF
	1327
	1328	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
	1329
	1330	ug = ug_vert(:,nt) + fac * ( ug_vert(:,nt+1) - ug_vert(:,nt) )
	1331	vg = vg_vert(:,nt) + fac * ( vg_vert(:,nt+1) - vg_vert(:,nt) )
	1332
	1333	IF ( large_scale_subsidence ) THEN
	1334	w_subs = wsubs_vert(:,nt) &
	1335	+ fac * ( wsubs_vert(:,nt+1) - wsubs_vert(:,nt) )
	1336	ENDIF
	1337
	1338	END SUBROUTINE ls_forcing_vert
	1339
	1340
	1341	!------------------------------------------------------------------------------!
	1342	! Description:
	1343	! ------------
	1344	!> Call for all grid points
	1345	!------------------------------------------------------------------------------!
	1346	SUBROUTINE ls_advec ( time, prog_var )
	1347
	1348
	1349	IMPLICIT NONE
	1350
	1351	CHARACTER (LEN=*) :: prog_var !<
	1352
	1353	REAL(wp), INTENT(in) :: time !<
	1354	REAL(wp) :: fac !<
	1355
	1356	INTEGER(iwp) :: i !<
	1357	INTEGER(iwp) :: j !<
	1358	INTEGER(iwp) :: k !<
	1359	INTEGER(iwp) :: nt !<
	1360
	1361	!
	1362	!-- Interpolation in time of LSF_DATA
	1363	nt = 1
	1364	DO WHILE ( time > time_vert(nt) )
	1365	nt = nt + 1
	1366	ENDDO
	1367	IF ( time /= time_vert(nt) ) THEN
	1368	nt = nt - 1
	1369	ENDIF
	1370
	1371	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
	1372
	1373	!
	1374	!-- Add horizontal large scale advection tendencies of pt and q
	1375	SELECT CASE ( prog_var )
	1376
	1377	CASE ( 'pt' )
	1378
	1379	DO i = nxl, nxr
	1380	DO j = nys, nyn
	1381	DO k = nzb+1, nzt
	1382	tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt) + fac * &
	1383	( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) ) *&
	1384	MERGE( 1.0_wp, 0.0_wp, &
	1385	BTEST( wall_flags_0(k,j,i), 0 ) )
	1386	ENDDO
	1387	ENDDO
	1388	ENDDO
	1389
	1390	CASE ( 'q' )
	1391
	1392	DO i = nxl, nxr
	1393	DO j = nys, nyn
	1394	DO k = nzb+1, nzt
	1395	tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt) + fac * &
	1396	( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) ) * &
	1397	MERGE( 1.0_wp, 0.0_wp, &
	1398	BTEST( wall_flags_0(k,j,i), 0 ) )
	1399	ENDDO
	1400	ENDDO
	1401	ENDDO
	1402
	1403	END SELECT
	1404
	1405	!
	1406	!-- Subsidence of pt and q with prescribed subsidence tendencies
	1407	IF ( large_scale_subsidence .AND. use_subsidence_tendencies ) THEN
	1408
	1409	SELECT CASE ( prog_var )
	1410
	1411	CASE ( 'pt' )
	1412
	1413	DO i = nxl, nxr
	1414	DO j = nys, nyn
	1415	DO k = nzb+1, nzt
	1416	tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac * &
	1417	( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )*&
	1418	MERGE( 1.0_wp, 0.0_wp, &
	1419	BTEST( wall_flags_0(k,j,i), 0 ) )
	1420	ENDDO
	1421	ENDDO
	1422	ENDDO
	1423
	1424	CASE ( 'q' )
	1425
	1426	DO i = nxl, nxr
	1427	DO j = nys, nyn
	1428	DO k = nzb+1, nzt
	1429	tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac * &
	1430	( td_sub_q(k,nt+1) - td_sub_q(k,nt) ) * &
	1431	MERGE( 1.0_wp, 0.0_wp, &
	1432	BTEST( wall_flags_0(k,j,i), 0 ) )
	1433	ENDDO
	1434	ENDDO
	1435	ENDDO
	1436
	1437	END SELECT
	1438
	1439	ENDIF
	1440
	1441	END SUBROUTINE ls_advec
	1442
	1443
	1444	!------------------------------------------------------------------------------!
	1445	! Description:
	1446	! ------------
	1447	!> Call for grid point i,j
	1448	!------------------------------------------------------------------------------!
	1449	SUBROUTINE ls_advec_ij ( i, j, time, prog_var )
	1450
	1451	IMPLICIT NONE
	1452
	1453	CHARACTER (LEN=*) :: prog_var !<
	1454
	1455	REAL(wp), INTENT(in) :: time !<
	1456	REAL(wp) :: fac !<
	1457
	1458	INTEGER(iwp) :: i !<
	1459	INTEGER(iwp) :: j !<
	1460	INTEGER(iwp) :: k !<
	1461	INTEGER(iwp) :: nt !<
	1462
	1463	!
	1464	!-- Interpolation in time of LSF_DATA
	1465	nt = 1
	1466	DO WHILE ( time > time_vert(nt) )
	1467	nt = nt + 1
	1468	ENDDO
	1469	IF ( time /= time_vert(nt) ) THEN
	1470	nt = nt - 1
	1471	ENDIF
	1472
	1473	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
	1474
	1475	!
	1476	!-- Add horizontal large scale advection tendencies of pt and q
	1477	SELECT CASE ( prog_var )
	1478
	1479	CASE ( 'pt' )
	1480
	1481	DO k = nzb+1, nzt
	1482	tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt) &
	1483	+ fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )*&
	1484	MERGE( 1.0_wp, 0.0_wp, &
	1485	BTEST( wall_flags_0(k,j,i), 0 ) )
	1486	ENDDO
	1487
	1488	CASE ( 'q' )
	1489
	1490	DO k = nzb+1, nzt
	1491	tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt) &
	1492	+ fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) ) * &
	1493	MERGE( 1.0_wp, 0.0_wp, &
	1494	BTEST( wall_flags_0(k,j,i), 0 ) )
	1495	ENDDO
	1496
	1497	END SELECT
	1498
	1499	!
	1500	!-- Subsidence of pt and q with prescribed profiles
	1501	IF ( large_scale_subsidence .AND. use_subsidence_tendencies ) THEN
	1502
	1503	SELECT CASE ( prog_var )
	1504
	1505	CASE ( 'pt' )
	1506
	1507	DO k = nzb+1, nzt
	1508	tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac * &
	1509	( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) ) * &
	1510	MERGE( 1.0_wp, 0.0_wp, &
	1511	BTEST( wall_flags_0(k,j,i), 0 ) )
	1512	ENDDO
	1513
	1514	CASE ( 'q' )
	1515
	1516	DO k = nzb+1, nzt
	1517	tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac * &
	1518	( td_sub_q(k,nt+1) - td_sub_q(k,nt) ) * &
	1519	MERGE( 1.0_wp, 0.0_wp, &
	1520	BTEST( wall_flags_0(k,j,i), 0 ) )
	1521	ENDDO
	1522
	1523	END SELECT
	1524
	1525	ENDIF
	1526
	1527	END SUBROUTINE ls_advec_ij
	1528
	1529
	1530	!------------------------------------------------------------------------------!
	1531	! Description:
	1532	! ------------
	1533	!> @todo Missing subroutine description.
	1534	!------------------------------------------------------------------------------!
	1535	SUBROUTINE nudge_init
	1536
	1537	IMPLICIT NONE
	1538
	1539
	1540	INTEGER(iwp) :: finput = 90 !<
	1541	INTEGER(iwp) :: ierrn !<
	1542	INTEGER(iwp) :: k !<
	1543	INTEGER(iwp) :: nt !<
	1544
	1545	CHARACTER(1) :: hash !<
	1546
	1547	REAL(wp) :: highheight !<
	1548	REAL(wp) :: highqnudge !<
	1549	REAL(wp) :: highptnudge !<
	1550	REAL(wp) :: highunudge !<
	1551	REAL(wp) :: highvnudge !<
	1552	REAL(wp) :: highwnudge !<
	1553	REAL(wp) :: hightnudge !<
	1554
	1555	REAL(wp) :: lowheight !<
	1556	REAL(wp) :: lowqnudge !<
	1557	REAL(wp) :: lowptnudge !<
	1558	REAL(wp) :: lowunudge !<
	1559	REAL(wp) :: lowvnudge !<
	1560	REAL(wp) :: lowwnudge !<
	1561	REAL(wp) :: lowtnudge !<
	1562
	1563	REAL(wp) :: fac !<
	1564
	1565	ALLOCATE( ptnudge(nzb:nzt+1,1:ntnudge), qnudge(nzb:nzt+1,1:ntnudge), &
	1566	tnudge(nzb:nzt+1,1:ntnudge), unudge(nzb:nzt+1,1:ntnudge), &
	1567	vnudge(nzb:nzt+1,1:ntnudge), wnudge(nzb:nzt+1,1:ntnudge) )
	1568
	1569	ALLOCATE( tmp_tnudge(nzb:nzt) )
	1570
	1571	ALLOCATE( timenudge(0:ntnudge) )
	1572
	1573	ptnudge = 0.0_wp; qnudge = 0.0_wp; tnudge = 0.0_wp; unudge = 0.0_wp
	1574	vnudge = 0.0_wp; wnudge = 0.0_wp; timenudge = 0.0_wp
	1575	!
	1576	!-- Initialize array tmp_nudge with a current nudging time scale of 6 hours
	1577	tmp_tnudge = 21600.0_wp
	1578
	1579	nt = 0
	1580	OPEN ( finput, FILE='NUDGING_DATA', STATUS='OLD', &
	1581	FORM='FORMATTED', IOSTAT=ierrn )
	1582
	1583	IF ( ierrn /= 0 ) THEN
	1584	message_string = 'file NUDGING_DATA does not exist'
	1585	CALL message( 'nudging', 'PA0365', 1, 2, 0, 6, 0 )
	1586	ENDIF
	1587
	1588	ierrn = 0
	1589
	1590	rloop:DO
	1591	nt = nt + 1
	1592	hash = "#"
	1593	ierrn = 1 ! not zero
	1594	!
	1595	!-- Search for the next line consisting of "# time",
	1596	!-- from there onwards the profiles will be read
	1597	DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
	1598
	1599	READ ( finput, *, IOSTAT=ierrn ) hash, timenudge(nt)
	1600	IF ( ierrn < 0 ) EXIT rloop
	1601
	1602	ENDDO
	1603
	1604	ierrn = 0
	1605	READ ( finput, *, IOSTAT=ierrn ) lowheight, lowtnudge, lowunudge, &
	1606	lowvnudge, lowwnudge , lowptnudge, &
	1607	lowqnudge
	1608
	1609	IF ( ierrn /= 0 ) THEN
	1610	message_string = 'errors in file NUDGING_DATA'
	1611	CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
	1612	ENDIF
	1613
	1614	ierrn = 0
	1615	READ ( finput, *, IOSTAT=ierrn ) highheight, hightnudge, highunudge, &
	1616	highvnudge, highwnudge , highptnudge, &
	1617	highqnudge
	1618
	1619	IF ( ierrn /= 0 ) THEN
	1620	message_string = 'errors in file NUDGING_DATA'
	1621	CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
	1622	ENDIF
	1623
	1624	DO k = nzb, nzt+1
	1625	DO WHILE ( highheight < zu(k) )
	1626	lowheight = highheight
	1627	lowtnudge = hightnudge
	1628	lowunudge = highunudge
	1629	lowvnudge = highvnudge
	1630	lowwnudge = highwnudge
	1631	lowptnudge = highptnudge
	1632	lowqnudge = highqnudge
	1633
	1634	ierrn = 0
[3045]	1635	READ ( finput, *, IOSTAT=ierrn ) highheight , hightnudge , &
	1636	highunudge , highvnudge , &
	1637	highwnudge , highptnudge, &
[2320]	1638	highqnudge
	1639	IF (ierrn /= 0 ) THEN
[3045]	1640	WRITE( message_string, * ) 'zu(',k,') = ', zu(k), 'm is ', &
	1641	'higher than the maximum height in NUDING_DATA which ',&
	1642	'is ', lowheight, 'm. Interpolation on PALM ', &
[2320]	1643	'grid is not possible.'
	1644	CALL message( 'nudging', 'PA0364', 1, 2, 0, 6, 0 )
	1645	ENDIF
	1646	ENDDO
	1647
	1648	!
	1649	!-- Interpolation of prescribed profiles in space
	1650
	1651	fac = ( highheight - zu(k) ) / ( highheight - lowheight )
	1652
	1653	tnudge(k,nt) = fac * lowtnudge + ( 1.0_wp - fac ) * hightnudge
	1654	unudge(k,nt) = fac * lowunudge + ( 1.0_wp - fac ) * highunudge
	1655	vnudge(k,nt) = fac * lowvnudge + ( 1.0_wp - fac ) * highvnudge
	1656	wnudge(k,nt) = fac * lowwnudge + ( 1.0_wp - fac ) * highwnudge
	1657	ptnudge(k,nt) = fac * lowptnudge + ( 1.0_wp - fac ) * highptnudge
	1658	qnudge(k,nt) = fac * lowqnudge + ( 1.0_wp - fac ) * highqnudge
	1659	ENDDO
	1660
	1661	ENDDO rloop
	1662
	1663	CLOSE ( finput )
	1664
	1665	!
	1666	!-- Overwrite initial profiles in case of nudging
	1667	IF ( nudging ) THEN
	1668	pt_init = ptnudge(:,1)
	1669	u_init = unudge(:,1)
	1670	v_init = vnudge(:,1)
	1671	IF ( humidity ) THEN ! is passive_scalar correct???
	1672	q_init = qnudge(:,1)
	1673	ENDIF
	1674
[3045]	1675	WRITE( message_string, * ) 'Initial profiles of u, v, pt and q ', &
	1676	'from NUDGING_DATA are used.'
[2320]	1677	CALL message( 'large_scale_forcing_nudging', 'PA0370', 0, 0, 0, 6, 0 )
	1678	ENDIF
	1679
	1680
	1681	END SUBROUTINE nudge_init
	1682
	1683	!------------------------------------------------------------------------------!
	1684	! Description:
	1685	! ------------
	1686	!> @todo Missing subroutine description.
	1687	!------------------------------------------------------------------------------!
	1688	SUBROUTINE calc_tnudge ( time )
	1689
	1690	IMPLICIT NONE
	1691
	1692
	1693	REAL(wp) :: dtm !<
	1694	REAL(wp) :: dtp !<
	1695	REAL(wp) :: time !<
	1696
	1697	INTEGER(iwp) :: k !<
	1698	INTEGER(iwp) :: nt !<
	1699
	1700	nt = 1
	1701	DO WHILE ( time > timenudge(nt) )
	1702	nt = nt+1
	1703	ENDDO
	1704	IF ( time /= timenudge(1) ) THEN
	1705	nt = nt-1
	1706	ENDIF
	1707
	1708	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
	1709	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
	1710
	1711	DO k = nzb, nzt
	1712	tmp_tnudge(k) = MAX( dt_3d, tnudge(k,nt) * dtp + tnudge(k,nt+1) * dtm )
	1713	ENDDO
	1714
	1715	END SUBROUTINE calc_tnudge
	1716
	1717	!------------------------------------------------------------------------------!
	1718	! Description:
	1719	! ------------
	1720	!> Call for all grid points
	1721	!------------------------------------------------------------------------------!
	1722	SUBROUTINE nudge ( time, prog_var )
	1723
	1724	IMPLICIT NONE
	1725
	1726	CHARACTER (LEN=*) :: prog_var !<
	1727
	1728	REAL(wp) :: tmp_tend !<
	1729	REAL(wp) :: dtm !<
	1730	REAL(wp) :: dtp !<
	1731	REAL(wp) :: time !<
	1732
	1733	INTEGER(iwp) :: i !<
	1734	INTEGER(iwp) :: j !<
	1735	INTEGER(iwp) :: k !<
	1736	INTEGER(iwp) :: nt !<
	1737
	1738
	1739	nt = 1
	1740	DO WHILE ( time > timenudge(nt) )
	1741	nt = nt+1
	1742	ENDDO
	1743	IF ( time /= timenudge(1) ) THEN
	1744	nt = nt-1
	1745	ENDIF
	1746
	1747	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
	1748	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
	1749
	1750	SELECT CASE ( prog_var )
	1751
	1752	CASE ( 'u' )
	1753
	1754	DO i = nxl, nxr
	1755	DO j = nys, nyn
	1756
	1757	DO k = nzb+1, nzt
	1758
	1759	tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp + &
	1760	unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
	1761
	1762	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
	1763	MERGE( 1.0_wp, 0.0_wp, &
	1764	BTEST( wall_flags_0(k,j,i), 1 ) )
	1765
	1766	sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend * &
	1767	weight_substep(intermediate_timestep_count)
	1768	ENDDO
	1769
	1770	sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
	1771
	1772	ENDDO
	1773	ENDDO
	1774
	1775	CASE ( 'v' )
	1776
	1777	DO i = nxl, nxr
	1778	DO j = nys, nyn
	1779
	1780	DO k = nzb+1, nzt
	1781
	1782	tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp + &
	1783	vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
	1784
	1785	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
	1786	MERGE( 1.0_wp, 0.0_wp, &
	1787	BTEST( wall_flags_0(k,j,i), 2 ) )
	1788
	1789	sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend * &
	1790	weight_substep(intermediate_timestep_count)
	1791	ENDDO
	1792
	1793	sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
	1794
	1795	ENDDO
	1796	ENDDO
	1797
	1798	CASE ( 'pt' )
	1799
	1800	DO i = nxl, nxr
	1801	DO j = nys, nyn
	1802
	1803	DO k = nzb+1, nzt
	1804
	1805	tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp + &
	1806	ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
	1807
	1808	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
	1809	MERGE( 1.0_wp, 0.0_wp, &
	1810	BTEST( wall_flags_0(k,j,i), 0 ) )
	1811
	1812	sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend * &
	1813	weight_substep(intermediate_timestep_count)
	1814	ENDDO
	1815
	1816	sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
	1817
	1818	ENDDO
	1819	ENDDO
	1820
	1821	CASE ( 'q' )
	1822
	1823	DO i = nxl, nxr
	1824	DO j = nys, nyn
	1825
	1826	DO k = nzb+1, nzt
	1827
	1828	tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp + &
	1829	qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
	1830
	1831	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
	1832	MERGE( 1.0_wp, 0.0_wp, &
	1833	BTEST( wall_flags_0(k,j,i), 0 ) )
	1834
	1835	sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend * &
	1836	weight_substep(intermediate_timestep_count)
	1837	ENDDO
	1838
	1839	sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
	1840
	1841	ENDDO
	1842	ENDDO
	1843
	1844	CASE DEFAULT
	1845	message_string = 'unknown prognostic variable "' // prog_var // '"'
	1846	CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
	1847
	1848	END SELECT
	1849
	1850	END SUBROUTINE nudge
	1851
	1852
	1853	!------------------------------------------------------------------------------!
	1854	! Description:
	1855	! ------------
	1856	!> Call for grid point i,j
	1857	!------------------------------------------------------------------------------!
	1858
	1859	SUBROUTINE nudge_ij( i, j, time, prog_var )
	1860
	1861	IMPLICIT NONE
	1862
	1863
	1864	CHARACTER (LEN=*) :: prog_var !<
	1865
	1866	REAL(wp) :: tmp_tend !<
	1867	REAL(wp) :: dtm !<
	1868	REAL(wp) :: dtp !<
	1869	REAL(wp) :: time !<
	1870
	1871	INTEGER(iwp) :: i !<
	1872	INTEGER(iwp) :: j !<
	1873	INTEGER(iwp) :: k !<
	1874	INTEGER(iwp) :: nt !<
	1875
	1876
	1877	nt = 1
	1878	DO WHILE ( time > timenudge(nt) )
	1879	nt = nt+1
	1880	ENDDO
	1881	IF ( time /= timenudge(1) ) THEN
	1882	nt = nt-1
	1883	ENDIF
	1884
	1885	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
	1886	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
	1887
	1888	SELECT CASE ( prog_var )
	1889
	1890	CASE ( 'u' )
	1891
	1892	DO k = nzb+1, nzt
	1893
	1894	tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp + &
	1895	unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
	1896
	1897	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
	1898	MERGE( 1.0_wp, 0.0_wp, &
	1899	BTEST( wall_flags_0(k,j,i), 1 ) )
	1900
	1901	sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend &
	1902	* weight_substep(intermediate_timestep_count)
	1903	ENDDO
	1904
	1905	sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
	1906
	1907	CASE ( 'v' )
	1908
	1909	DO k = nzb+1, nzt
	1910
	1911	tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp + &
	1912	vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
	1913
	1914	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
	1915	MERGE( 1.0_wp, 0.0_wp, &
	1916	BTEST( wall_flags_0(k,j,i), 2 ) )
	1917
	1918	sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend &
	1919	* weight_substep(intermediate_timestep_count)
	1920	ENDDO
	1921
	1922	sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
	1923
	1924	CASE ( 'pt' )
	1925
	1926	DO k = nzb+1, nzt
	1927
	1928	tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp + &
	1929	ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
	1930
	1931	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
	1932	MERGE( 1.0_wp, 0.0_wp, &
	1933	BTEST( wall_flags_0(k,j,i), 0 ) )
	1934
	1935	sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend &
	1936	* weight_substep(intermediate_timestep_count)
	1937	ENDDO
	1938
	1939	sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
	1940
	1941
	1942	CASE ( 'q' )
	1943
	1944	DO k = nzb+1, nzt
	1945
	1946	tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp + &
	1947	qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
	1948
	1949	tend(k,j,i) = tend(k,j,i) + tmp_tend * &
	1950	MERGE( 1.0_wp, 0.0_wp, &
	1951	BTEST( wall_flags_0(k,j,i), 0 ) )
	1952
	1953	sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend &
	1954	* weight_substep(intermediate_timestep_count)
	1955	ENDDO
	1956
	1957	sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
	1958
	1959	CASE DEFAULT
	1960	message_string = 'unknown prognostic variable "' // prog_var // '"'
	1961	CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
	1962
	1963	END SELECT
	1964
	1965
	1966	END SUBROUTINE nudge_ij
	1967
	1968
	1969	!------------------------------------------------------------------------------!
	1970	! Description:
	1971	! ------------
	1972	!> @todo Missing subroutine description.
	1973	!------------------------------------------------------------------------------!
	1974	SUBROUTINE nudge_ref ( time )
	1975
	1976	IMPLICIT NONE
	1977
	1978	INTEGER(iwp) :: nt !<
	1979
	1980	REAL(wp) :: fac !<
	1981	REAL(wp), INTENT(in) :: time !<
	1982
	1983	!
	1984	!-- Interpolation in time of NUDGING_DATA for pt_init and q_init. This is
	1985	!-- needed for correct upper boundary conditions for pt and q and in case that
	1986	! large scale subsidence as well as scalar Rayleigh-damping are used
	1987	nt = 1
	1988	DO WHILE ( time > time_vert(nt) )
	1989	nt = nt + 1
	1990	ENDDO
	1991	IF ( time /= time_vert(nt) ) THEN
	1992	nt = nt - 1
	1993	ENDIF
	1994
	1995	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
	1996
	1997	pt_init = ptnudge(:,nt) + fac * ( ptnudge(:,nt+1) - ptnudge(:,nt) )
	1998	q_init = qnudge(:,nt) + fac * ( qnudge(:,nt+1) - qnudge(:,nt) )
	1999	u_init = unudge(:,nt) + fac * ( unudge(:,nt+1) - unudge(:,nt) )
	2000	v_init = vnudge(:,nt) + fac * ( vnudge(:,nt+1) - vnudge(:,nt) )
	2001
	2002	END SUBROUTINE nudge_ref
	2003
[2696]	2004
[3026]	2005	END MODULE lsf_nudging_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |