Home

Context Navigation

source: palm/trunk/SOURCE/surface_layer_fluxes_mod.f90 @ 4759

Last change on this file since 4759 was 4750, checked in by suehring, 4 years ago
Indoor-model new: Namelist parameter added to switch-off/on the indoor model during wall/soil spinup; bugfixes concerning indoor model: bugfix in window-wall treatment during spinup - in the urban-surface model the window fraction is set to zero during spinup, so it is done here also; bugfix in wall treatment - inner wall temperature was too low due to wrong weighting of wall/green/window fractions; Revision of 10-cm temperature at vertical walls - assume grid-cell temperature rather than employ MOST; call hourly-based indoor model only once per hour during spinup, not every timestep; add missing dependency in Makefile; urban-surface model: bugfix in openmp directive
Property svn:keywords set to `Id`
File size: 83.1 KB

Rev	Line
[1850]	1	!> @file surface_layer_fluxes_mod.f90
[4562]	2	!--------------------------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1691]	4	!
[4562]	5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
	6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
	7	! (at your option) any later version.
[1691]	8	!
[4562]	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
	10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	11	! Public License for more details.
[1691]	12	!
[4562]	13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
	14	! <http://www.gnu.org/licenses/>.
[1691]	15	!
[4360]	16	! Copyright 1997-2020 Leibniz Universitaet Hannover
[4562]	17	!--------------------------------------------------------------------------------------------------!
[1691]	18	!
[2696]	19	!
[1691]	20	! Current revisions:
[4562]	21	! -----------------
[4593]	22	!
[4691]	23	!
[1692]	24	! Former revisions:
	25	! -----------------
	26	! $Id: surface_layer_fluxes_mod.f90 4750 2020-10-16 14:27:48Z suehring $
[4750]	27	! Revision of 10-cm temperature at vertical walls - assume grid-cell temperature rather than employ
	28	! MOST
	29	!
	30	! 4717 2020-09-30 22:27:40Z pavelkrc
[4717]	31	! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
	32	! directives (J. Resler)
	33	!
	34	! 4691 2020-09-22 14:38:38Z suehring
[4691]	35	! Bugfix for commit 4593 in vector branch of calc_ol
	36	!
	37	! 4671 2020-09-09 20:27:58Z pavelkrc
[4671]	38	! Implementation of downward facing USM and LSM surfaces
	39	!
	40	! 4594 2020-07-09 15:01:00Z suehring
[4594]	41	! Include k index in OMP PRIVATE statements
	42	!
	43	! 4593 2020-07-09 12:48:18Z suehring
[4593]	44	! - Pre-calculate ln(z/z0) at each timestep in order to reduce the number of log-calculations
	45	! - Bugfix - add missing density to fluxes of passive-scalars, chemistry and cloud-phyiscal
	46	! quantities at upward-facing surfaces
	47	! - Move if-statement out of inner loop
	48	! - Remove unnecessary index referencing
	49	!
	50	! 4562 2020-06-12 08:38:47Z raasch
[4562]	51	! File re-formatted to follow the PALM coding standard
	52	!
	53	! 4519 2020-05-05 17:33:30Z suehring
[4519]	54	! Add missing computation of passive scalar scaling parameter
[4562]	55	!
[4519]	56	! 4370 2020-01-10 14:00:44Z raasch
[4562]	57	! Bugfix: openacc porting for vector version of OL calculation added
	58	!
[4370]	59	! 4366 2020-01-09 08:12:43Z raasch
[4562]	60	! Vector version for calculation of Obukhov length via Newton iteration added
	61	!
[4366]	62	! 4360 2020-01-07 11:25:50Z suehring
[4562]	63	! Calculation of diagnostic-only 2-m potential temperature moved to diagnostic_output_quantities.
	64	!
[4331]	65	! 4298 2019-11-21 15:59:16Z suehring
[4562]	66	! Calculation of 2-m temperature adjusted to the case the 2-m level is above the first grid point.
	67	!
[4298]	68	! 4258 2019-10-07 13:29:08Z suehring
[4258]	69	! Initialization of Obukhov lenght also at vertical surfaces (if allocated).
[4562]	70	!
[4258]	71	! 4237 2019-09-25 11:33:42Z knoop
[4237]	72	! Added missing OpenMP directives
[4562]	73	!
[4237]	74	! 4186 2019-08-23 16:06:14Z suehring
[4562]	75	! - To enable limitation of Obukhov length, move it before exit-cycle construct.
	76	! Further, change the limit to 10E-5 in order to get rid-off unrealistic peaks in the heat fluxes
	77	! during nighttime
[4186]	78	! - Unused variable removed
[4562]	79	!
[4186]	80	! 4182 2019-08-22 15:20:23Z scharf
[4182]	81	! Corrected "Former revisions" section
[4562]	82	!
[4182]	83	! 3987 2019-05-22 09:52:13Z kanani
[3987]	84	! Introduce alternative switch for debug output during timestepping
[4562]	85	!
[3987]	86	! 3885 2019-04-11 11:29:34Z kanani
[4562]	87	! Changes related to global restructuring of location messages and introduction of additional debug
	88	! messages
	89	!
[3885]	90	! 3881 2019-04-10 09:31:22Z suehring
[3881]	91	! Assure that Obukhov length does not become zero
[4562]	92	!
[3881]	93	! 3834 2019-03-28 15:40:15Z forkel
[4562]	94	! Added USE chem_gasphase_mod
	95	!
[3833]	96	! 3787 2019-03-07 08:43:54Z raasch
[4562]	97	! Unused variables removed
	98	!
[3787]	99	! 3745 2019-02-15 18:57:56Z suehring
[4562]	100	! Bugfix, missing calculation of 10cm temperature at vertical building walls, required for indoor
	101	! model
	102	!
[3745]	103	! 3744 2019-02-15 18:38:58Z suehring
[3685]	104	! Some interface calls moved to module_interface + cleanup
[4562]	105	!
[3685]	106	! 3668 2019-01-14 12:49:24Z maronga
[3668]	107	! Removed methods "circular" and "lookup"
[4562]	108	!
[3668]	109	! 3655 2019-01-07 16:51:22Z knoop
[3634]	110	! OpenACC port for SPEC
[1692]	111	!
[4182]	112	! Revision 1.1 1998/01/23 10:06:06 raasch
	113	! Initial revision
	114	!
	115	!
[1691]	116	! Description:
	117	! ------------
[4562]	118	!> Diagnostic computation of vertical fluxes in the constant flux layer from the values of the
	119	!> variables at grid point k=1 based on Newton iteration.
[1691]	120	!>
[4562]	121	!> @todo (Re)move large_scale_forcing actions
	122	!> @todo Check/optimize OpenMP directives
	123	!> @todo Simplify if conditions (which flux need to be computed in which case)
	124	!--------------------------------------------------------------------------------------------------!
[1691]	125	MODULE surface_layer_fluxes_mod
	126
[4562]	127	USE arrays_3d, &
	128	ONLY: d_exner, &
	129	drho_air_zw, &
	130	e, &
	131	kh, &
	132	nc, &
	133	nr, &
	134	pt, &
	135	q, &
	136	ql, &
	137	qc, &
	138	qr, &
	139	s, &
	140	u, &
	141	v, &
	142	vpt, &
	143	w, &
	144	zu, &
	145	zw, &
	146	rho_air_zw
[1691]	147
[3274]	148
[4562]	149	USE basic_constants_and_equations_mod, &
	150	ONLY: g, &
	151	kappa, &
	152	lv_d_cp, &
	153	pi, &
	154	rd_d_rv
	155
	156	USE chem_gasphase_mod, &
[3833]	157	ONLY: nvar
	158
[4562]	159	USE chem_modules, &
[3834]	160	ONLY: constant_csflux
[2696]	161
[1691]	162	USE cpulog
	163
[4562]	164	USE control_parameters, &
	165	ONLY: air_chemistry, &
	166	cloud_droplets, &
	167	constant_heatflux, &
	168	constant_scalarflux, &
	169	constant_waterflux, &
	170	coupling_mode, &
	171	debug_output_timestep, &
	172	humidity, &
	173	ibc_e_b, &
	174	ibc_pt_b, &
	175	indoor_model, &
	176	land_surface, &
	177	large_scale_forcing, &
	178	loop_optimization, &
	179	lsf_surf, &
	180	message_string, &
	181	neutral, &
	182	passive_scalar, &
	183	pt_surface, &
	184	q_surface, &
	185	run_coupled, &
	186	surface_pressure, &
	187	simulated_time, &
	188	time_since_reference_point, &
	189	urban_surface, &
	190	use_free_convection_scaling, &
	191	zeta_max, &
	192	zeta_min
[1691]	193
[4562]	194	USE grid_variables, &
	195	ONLY: dx, &
	196	dy
[2232]	197
[4562]	198	USE indices, &
[4298]	199	ONLY: nzt
[1691]	200
	201	USE kinds
	202
[4562]	203	USE bulk_cloud_model_mod, &
	204	ONLY: bulk_cloud_model, &
	205	microphysics_morrison, &
	206	microphysics_seifert
[3274]	207
[1691]	208	USE pegrid
	209
[4562]	210	USE land_surface_model_mod, &
	211	ONLY: aero_resist_kray, &
	212	skip_time_do_lsm
[2011]	213
[4562]	214	USE surface_mod, &
	215	ONLY : surf_def_h, &
	216	surf_def_v, &
	217	surf_lsm_h, &
	218	surf_lsm_v, &
	219	surf_type, &
	220	surf_usm_h, &
	221	surf_usm_v
[1691]	222
[4562]	223
[1691]	224	IMPLICIT NONE
	225
[4562]	226	INTEGER(iwp) :: i !< loop index x direction
	227	INTEGER(iwp) :: j !< loop index y direction
	228	INTEGER(iwp) :: k !< loop index z direction
	229	INTEGER(iwp) :: l !< loop index for surf type
[1691]	230
[4562]	231	LOGICAL :: coupled_run !< Flag for coupled atmosphere-ocean runs
	232	LOGICAL :: downward = .FALSE. !< Flag indicating downward-facing horizontal surface
	233	LOGICAL :: mom_uv = .FALSE. !< Flag indicating calculation of usvs and vsus at vertical surfaces
	234	LOGICAL :: mom_w = .FALSE. !< Flag indicating calculation of wsus and wsvs at vertical surfaces
	235	LOGICAL :: mom_tke = .FALSE. !< Flag indicating calculation of momentum fluxes at vertical surfaces used for TKE production
	236	LOGICAL :: surf_vertical !< Flag indicating vertical surfaces
[1691]	237
[4562]	238	REAL(wp) :: e_s !< Saturation water vapor pressure
	239	REAL(wp) :: ol_max = 1.0E6_wp !< Maximum Obukhov length
	240	REAL(wp) :: z_mo !< Height of the constant flux layer where MOST is assumed
[1691]	241
[4562]	242	TYPE(surf_type), POINTER :: surf !< surf-type array, used to generalize subroutines
[1691]	243
[2232]	244
[1691]	245	SAVE
	246
	247	PRIVATE
	248
[4562]	249	PUBLIC init_surface_layer_fluxes, &
	250	phi_m, &
	251	psi_h, &
	252	psi_m, &
[4331]	253	surface_layer_fluxes
[1691]	254
	255	INTERFACE init_surface_layer_fluxes
	256	MODULE PROCEDURE init_surface_layer_fluxes
	257	END INTERFACE init_surface_layer_fluxes
	258
[3130]	259	INTERFACE phi_m
	260	MODULE PROCEDURE phi_m
	261	END INTERFACE phi_m
	262
[4331]	263	INTERFACE psi_h
	264	MODULE PROCEDURE psi_h
	265	END INTERFACE psi_h
	266
	267	INTERFACE psi_m
	268	MODULE PROCEDURE psi_m
	269	END INTERFACE psi_m
	270
[1691]	271	INTERFACE surface_layer_fluxes
	272	MODULE PROCEDURE surface_layer_fluxes
	273	END INTERFACE surface_layer_fluxes
	274
	275
	276	CONTAINS
	277
	278
[4562]	279	!--------------------------------------------------------------------------------------------------!
[1691]	280	! Description:
	281	! ------------
[4562]	282	!> Main routine to compute the surface fluxes.
	283	!--------------------------------------------------------------------------------------------------!
	284	SUBROUTINE surface_layer_fluxes
[1691]	285
[4562]	286	IMPLICIT NONE
[1691]	287
[3885]	288
[4562]	289	IF ( debug_output_timestep ) CALL debug_message( 'surface_layer_fluxes', 'start' )
[3885]	290
[4562]	291	surf_vertical = .FALSE. !< flag indicating vertically orientated surface elements
	292	downward = .FALSE. !< flag indicating downward-facing surface elements
[1691]	293	!
[4593]	294	!-- First, precalculate ln(z/z0) for all surfaces. This is done each timestep, in order to account
	295	!-- for time-dependent roughness or user-modifications.
	296	DO l = 0, 1
	297	IF ( surf_def_h(l)%ns >= 1 ) THEN
	298	surf => surf_def_h(l)
	299	CALL calc_ln
	300	ENDIF
[4671]	301	IF ( surf_lsm_h(l)%ns >= 1 ) THEN
	302	surf => surf_lsm_h(l)
	303	CALL calc_ln
	304	ENDIF
	305	IF ( surf_usm_h(l)%ns >= 1 ) THEN
	306	surf => surf_usm_h(l)
	307	CALL calc_ln
	308	ENDIF
[4593]	309	ENDDO
	310
	311	DO l = 0, 3
	312	IF ( surf_def_v(l)%ns >= 1 ) THEN
	313	surf => surf_def_v(l)
	314	CALL calc_ln
	315	ENDIF
	316	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	317	surf => surf_lsm_v(l)
	318	CALL calc_ln
	319	ENDIF
	320	IF ( surf_usm_v(l)%ns >= 1 ) THEN
	321	surf => surf_usm_v(l)
	322	CALL calc_ln
	323	ENDIF
	324	ENDDO
	325	!
[4562]	326	!-- Derive potential temperature and specific humidity at first grid level from the fields pt and q
	327	DO l = 0, 1
[4671]	328	!-- First call for horizontal default-type surfaces (l=0 - upward facing, l=1 - downward facing)
[4562]	329	IF ( surf_def_h(l)%ns >= 1 ) THEN
	330	surf => surf_def_h(l)
	331	CALL calc_pt_q
	332	IF ( .NOT. neutral ) THEN
	333	CALL calc_pt_surface
	334	IF ( humidity ) THEN
	335	CALL calc_q_surface
	336	CALL calc_vpt_surface
[3146]	337	ENDIF
[2696]	338	ENDIF
[4562]	339	ENDIF
[2232]	340	!
[4671]	341	!-- Call for natural-type horizontal surfaces
	342	IF ( surf_lsm_h(l)%ns >= 1 ) THEN
	343	surf => surf_lsm_h(l)
	344	CALL calc_pt_q
	345	ENDIF
[4562]	346	!
[4671]	347	!-- Call for urban-type horizontal surfaces
	348	IF ( surf_usm_h(l)%ns >= 1 ) THEN
	349	surf => surf_usm_h(l)
	350	CALL calc_pt_q
	351	ENDIF
	352	ENDDO
[4562]	353	!
	354	!-- Call for natural-type vertical surfaces
	355	DO l = 0, 3
	356	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	357	surf => surf_lsm_v(l)
[2696]	358	CALL calc_pt_q
	359	ENDIF
	360
[4562]	361	!-- Call for urban-type vertical surfaces
	362	IF ( surf_usm_v(l)%ns >= 1 ) THEN
	363	surf => surf_usm_v(l)
[2696]	364	CALL calc_pt_q
	365	ENDIF
[4562]	366	ENDDO
[2696]	367
	368	!
[4593]	369	!-- First, calculate the new Obukhov length from precalculated values of log(z/z0) and wind speeds.
	370	!-- As a second step, then calculate new friction velocity, followed by the new scaling
[4562]	371	!-- parameters (th, q, etc.), and the new surface fluxes, if required. Note, each routine is called
	372	!-- for different surface types. First call for default-type horizontal surfaces, for natural- and
	373	!-- urban-type surfaces. Note, here only upward-facing horizontal surfaces are treated.
[4593]	374	!-- Note, calculation of log(z/z0) is redone each timestep, in order to account for time-dependent
	375	!-- values.
	376	!-- Start with default-type upward-facing horizontal surfaces
[4562]	377	IF ( surf_def_h(0)%ns >= 1 ) THEN
	378	surf => surf_def_h(0)
	379	CALL calc_uvw_abs
	380	IF ( .NOT. neutral ) CALL calc_ol
	381	CALL calc_us
	382	CALL calc_scaling_parameters
	383	CALL calc_surface_fluxes
	384	ENDIF
	385	!
	386	!-- Natural-type horizontal surfaces
[4671]	387	IF ( surf_lsm_h(0)%ns >= 1 ) THEN
	388	surf => surf_lsm_h(0)
[4562]	389	CALL calc_uvw_abs
	390	IF ( .NOT. neutral ) CALL calc_ol
	391	CALL calc_us
	392	CALL calc_scaling_parameters
	393	CALL calc_surface_fluxes
	394	ENDIF
	395	!
	396	!-- Urban-type horizontal surfaces
[4671]	397	IF ( surf_usm_h(0)%ns >= 1 ) THEN
	398	surf => surf_usm_h(0)
[4562]	399	CALL calc_uvw_abs
	400	IF ( .NOT. neutral ) CALL calc_ol
	401	CALL calc_us
	402	CALL calc_scaling_parameters
	403	CALL calc_surface_fluxes
	404	!
	405	!-- Calculate 10cm temperature, required in indoor model
	406	IF ( indoor_model ) CALL calc_pt_near_surface ( '10cm' )
	407	ENDIF
[3668]	408
[2232]	409	!
[4671]	410	!-- Treat downward-facing horizontal surfaces.
[4562]	411	!-- Stratification is not considered in this case, hence, no further distinction between different
	412	!-- most_method is required.
[4671]	413	downward = .TRUE.
	414	!-- Default type.
[4562]	415	IF ( surf_def_h(1)%ns >= 1 ) THEN
	416	surf => surf_def_h(1)
	417	CALL calc_uvw_abs
	418	CALL calc_us
	419	CALL calc_surface_fluxes
	420	ENDIF
[4671]	421	!-- Natural surface type.
	422	IF ( surf_lsm_h(1)%ns >= 1 ) THEN
	423	surf => surf_lsm_h(1)
	424	CALL calc_uvw_abs
	425	CALL calc_us
	426	CALL calc_surface_fluxes
	427	ENDIF
	428	!-- Urban surface type.
	429	IF ( surf_usm_h(1)%ns >= 1 ) THEN
	430	surf => surf_usm_h(1)
	431	CALL calc_uvw_abs
	432	CALL calc_us
	433	CALL calc_surface_fluxes
	434	ENDIF
	435	downward = .FALSE.
[1691]	436	!
[4562]	437	!-- Calculate surfaces fluxes at vertical surfaces for momentum and subgrid-scale TKE. No stability
	438	!-- is considered. Therefore, scaling parameters and Obukhov length do not need to be calculated and
	439	!-- no distinction in 'circular', 'Newton' or 'lookup' is necessary so far. Note, this will change
	440	!-- if stability is once considered.
	441	surf_vertical = .TRUE.
[2232]	442	!
[4562]	443	!-- Calculate horizontal momentum fluxes at north- and south-facing surfaces(usvs).
	444	!-- For default-type surfaces
	445	mom_uv = .TRUE.
	446	DO l = 0, 1
	447	IF ( surf_def_v(l)%ns >= 1 ) THEN
	448	surf => surf_def_v(l)
[3744]	449	!
[4562]	450	!-- Compute surface-parallel velocity
	451	CALL calc_uvw_abs_v_ugrid
[2232]	452	!
[4562]	453	!-- Compute respective friction velocity on staggered grid
[1691]	454	CALL calc_us
[4562]	455	!
	456	!-- Compute respective surface fluxes for momentum and TKE
[1691]	457	CALL calc_surface_fluxes
	458	ENDIF
[4562]	459	ENDDO
[2232]	460	!
[4562]	461	!-- For natural-type surfaces. Please note, even though stability is not considered for the
	462	!-- calculation of momentum fluxes at vertical surfaces, scaling parameters and Obukhov length are
	463	!-- calculated nevertheless in this case. This is due to the requirement of ts in parameterization
	464	!-- of heat flux in land-surface model in case that aero_resist_kray is not true.
	465	IF ( .NOT. aero_resist_kray ) THEN
	466	DO l = 0, 1
	467	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	468	surf => surf_lsm_v(l)
[2232]	469	!
[4562]	470	!-- Compute surface-parallel velocity
	471	CALL calc_uvw_abs_v_ugrid
	472	!
	473	!-- Compute Obukhov length
	474	IF ( .NOT. neutral ) CALL calc_ol
	475	!
	476	!-- Compute respective friction velocity on staggered grid
	477	CALL calc_us
	478	!
	479	!-- Compute scaling parameters
	480	CALL calc_scaling_parameters
	481	!
	482	!-- Compute respective surface fluxes for momentum and TKE
	483	CALL calc_surface_fluxes
	484	ENDIF
	485	ENDDO
	486	!
	487	!-- No ts is required, so scaling parameters and Obukhov length do not need to be computed.
	488	ELSE
[2232]	489	DO l = 0, 1
[4562]	490	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	491	surf => surf_lsm_v(l)
[2232]	492	!
	493	!-- Compute surface-parallel velocity
	494	CALL calc_uvw_abs_v_ugrid
	495	!
	496	!-- Compute respective friction velocity on staggered grid
	497	CALL calc_us
	498	!
	499	!-- Compute respective surface fluxes for momentum and TKE
	500	CALL calc_surface_fluxes
	501	ENDIF
	502	ENDDO
[4562]	503	ENDIF
[2232]	504	!
[4562]	505	!-- For urban-type surfaces
	506	DO l = 0, 1
	507	IF ( surf_usm_v(l)%ns >= 1 ) THEN
	508	surf => surf_usm_v(l)
[2232]	509	!
[4562]	510	!-- Compute surface-parallel velocity
	511	CALL calc_uvw_abs_v_ugrid
[2232]	512	!
[4562]	513	!-- Compute respective friction velocity on staggered grid
	514	CALL calc_us
[2232]	515	!
[4562]	516	!-- Compute respective surface fluxes for momentum and TKE
	517	CALL calc_surface_fluxes
[2232]	518	!
[4562]	519	!-- Calculate 10cm temperature, required in indoor model
	520	IF ( indoor_model ) CALL calc_pt_near_surface ( '10cm' )
	521	ENDIF
	522	ENDDO
[2232]	523	!
[4562]	524	!-- Calculate horizontal momentum fluxes at east- and west-facing surfaces (vsus).
	525	!-- For default-type surfaces
	526	DO l = 2, 3
	527	IF ( surf_def_v(l)%ns >= 1 ) THEN
	528	surf => surf_def_v(l)
[2232]	529	!
[4562]	530	!-- Compute surface-parallel velocity
	531	CALL calc_uvw_abs_v_vgrid
[2232]	532	!
[4562]	533	!-- Compute respective friction velocity on staggered grid
	534	CALL calc_us
[2232]	535	!
[4562]	536	!-- Compute respective surface fluxes for momentum and TKE
	537	CALL calc_surface_fluxes
[2232]	538	ENDIF
[4562]	539	ENDDO
[2232]	540	!
[4562]	541	!-- For natural-type surfaces. Please note, even though stability is not considered for the
	542	!-- calculation of momentum fluxes at vertical surfaces, scaling parameters and Obukov length are
	543	!-- calculated nevertheless in this case. This is due to the requirement of ts in parameterization
	544	!-- of heat flux in land-surface model in case that aero_resist_kray is not true.
	545	IF ( .NOT. aero_resist_kray ) THEN
	546	DO l = 2, 3
	547	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	548	surf => surf_lsm_v(l)
[2232]	549	!
	550	!-- Compute surface-parallel velocity
[4562]	551	CALL calc_uvw_abs_v_vgrid
[2232]	552	!
[4562]	553	!-- Compute Obukhov length
	554	IF ( .NOT. neutral ) CALL calc_ol
	555	!
[2232]	556	!-- Compute respective friction velocity on staggered grid
	557	CALL calc_us
	558	!
[4562]	559	!-- Compute scaling parameters
	560	CALL calc_scaling_parameters
	561	!
[2232]	562	!-- Compute respective surface fluxes for momentum and TKE
	563	CALL calc_surface_fluxes
	564	ENDIF
	565	ENDDO
[4562]	566	ELSE
[2232]	567	DO l = 2, 3
[4562]	568	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	569	surf => surf_lsm_v(l)
[2232]	570	!
	571	!-- Compute surface-parallel velocity
	572	CALL calc_uvw_abs_v_vgrid
	573	!
	574	!-- Compute respective friction velocity on staggered grid
	575	CALL calc_us
	576	!
	577	!-- Compute respective surface fluxes for momentum and TKE
	578	CALL calc_surface_fluxes
	579	ENDIF
	580	ENDDO
[4562]	581	ENDIF
[2232]	582	!
[4562]	583	!-- For urban-type surfaces
	584	DO l = 2, 3
	585	IF ( surf_usm_v(l)%ns >= 1 ) THEN
	586	surf => surf_usm_v(l)
[2232]	587	!
[4562]	588	!-- Compute surface-parallel velocity
	589	CALL calc_uvw_abs_v_vgrid
[2232]	590	!
[4562]	591	!-- Compute respective friction velocity on staggered grid
	592	CALL calc_us
[2232]	593	!
[4562]	594	!-- Compute respective surface fluxes for momentum and TKE
	595	CALL calc_surface_fluxes
[2232]	596	!
[4562]	597	!-- Calculate 10cm temperature, required in indoor model
	598	IF ( indoor_model ) CALL calc_pt_near_surface ( '10cm' )
	599	ENDIF
	600	ENDDO
	601	mom_uv = .FALSE.
[2232]	602	!
[4562]	603	!-- Calculate horizontal momentum fluxes of w (wsus and wsvs) at vertial surfaces.
	604	mom_w = .TRUE.
[2232]	605	!
[4562]	606	!-- Default-type surfaces
	607	DO l = 0, 3
	608	IF ( surf_def_v(l)%ns >= 1 ) THEN
	609	surf => surf_def_v(l)
	610	CALL calc_uvw_abs_v_wgrid
	611	CALL calc_us
	612	CALL calc_surface_fluxes
	613	ENDIF
	614	ENDDO
[2232]	615	!
[4562]	616	!-- Natural-type surfaces
	617	DO l = 0, 3
	618	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	619	surf => surf_lsm_v(l)
	620	CALL calc_uvw_abs_v_wgrid
	621	CALL calc_us
	622	CALL calc_surface_fluxes
[2232]	623	ENDIF
[4562]	624	ENDDO
[2232]	625	!
[4562]	626	!-- Urban-type surfaces
	627	DO l = 0, 3
	628	IF ( surf_usm_v(l)%ns >= 1 ) THEN
	629	surf => surf_usm_v(l)
	630	CALL calc_uvw_abs_v_wgrid
	631	CALL calc_us
	632	CALL calc_surface_fluxes
	633	ENDIF
	634	ENDDO
	635	mom_w = .FALSE.
[2232]	636	!
[4562]	637	!-- Calculate momentum fluxes usvs, vsus, wsus and wsvs at vertical surfaces for TKE production.
	638	!-- Note, here, momentum fluxes are defined at grid center and are not staggered as before.
	639	mom_tke = .TRUE.
[2232]	640	!
[4562]	641	!-- Default-type surfaces
	642	DO l = 0, 3
	643	IF ( surf_def_v(l)%ns >= 1 ) THEN
	644	surf => surf_def_v(l)
	645	CALL calc_uvw_abs_v_sgrid
	646	CALL calc_us
	647	CALL calc_surface_fluxes
	648	ENDIF
	649	ENDDO
[2232]	650	!
[4562]	651	!-- Natural-type surfaces
	652	DO l = 0, 3
	653	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	654	surf => surf_lsm_v(l)
	655	CALL calc_uvw_abs_v_sgrid
	656	CALL calc_us
	657	CALL calc_surface_fluxes
	658	ENDIF
	659	ENDDO
[3744]	660	!
[4562]	661	!-- Urban-type surfaces
	662	DO l = 0, 3
	663	IF ( surf_usm_v(l)%ns >= 1 ) THEN
	664	surf => surf_usm_v(l)
	665	CALL calc_uvw_abs_v_sgrid
	666	CALL calc_us
	667	CALL calc_surface_fluxes
	668	ENDIF
	669	ENDDO
	670	mom_tke = .FALSE.
[1691]	671
[4562]	672	IF ( debug_output_timestep ) CALL debug_message( 'surface_layer_fluxes', 'end' )
[3885]	673
[4562]	674	END SUBROUTINE surface_layer_fluxes
[1691]	675
	676
[4562]	677	!--------------------------------------------------------------------------------------------------!
[1691]	678	! Description:
	679	! ------------
[4258]	680	!> Initializing actions for the surface layer routine.
[4562]	681	!--------------------------------------------------------------------------------------------------!
	682	SUBROUTINE init_surface_layer_fluxes
[1691]	683
[4562]	684	IMPLICIT NONE
[1691]	685
[4562]	686	INTEGER(iwp) :: l !< running index for vertical surface orientation
[1691]	687
[4562]	688	CALL location_message( 'initializing surface layer', 'start' )
[1709]	689
	690	!
[4562]	691	!-- In case of runs with neutral statification, set Obukhov length to a large value
	692	IF ( neutral ) THEN
[4671]	693	DO l = 0, 1
	694	IF ( surf_def_h(l)%ns >= 1 .AND. &
	695	ALLOCATED( surf_def_h(l)%ol ) ) surf_def_h(l)%ol = 1.0E10_wp
	696	IF ( surf_lsm_h(l)%ns >= 1 .AND. &
	697	ALLOCATED( surf_lsm_h(l)%ol ) ) surf_lsm_h(l)%ol = 1.0E10_wp
	698	IF ( surf_usm_h(l)%ns >= 1 .AND. &
	699	ALLOCATED( surf_usm_h(l)%ol ) ) surf_usm_h(l)%ol = 1.0E10_wp
	700	ENDDO
[4562]	701	DO l = 0, 3
	702	IF ( surf_def_v(l)%ns >= 1 .AND. &
	703	ALLOCATED( surf_def_v(l)%ol ) ) surf_def_v(l)%ol = 1.0E10_wp
	704	IF ( surf_lsm_v(l)%ns >= 1 .AND. &
	705	ALLOCATED( surf_lsm_v(l)%ol ) ) surf_lsm_v(l)%ol = 1.0E10_wp
	706	IF ( surf_usm_v(l)%ns >= 1 .AND. &
	707	ALLOCATED( surf_usm_v(l)%ol ) ) surf_usm_v(l)%ol = 1.0E10_wp
	708	ENDDO
[3685]	709
[4562]	710	ENDIF
[1691]	711
[4562]	712	CALL location_message( 'initializing surface layer', 'finished' )
[1691]	713
[4562]	714	END SUBROUTINE init_surface_layer_fluxes
	715
	716
	717	!--------------------------------------------------------------------------------------------------!
[1691]	718	! Description:
	719	! ------------
[4593]	720	!> Compute ln(z/z0).
	721	!--------------------------------------------------------------------------------------------------!
	722	SUBROUTINE calc_ln
	723
	724	INTEGER(iwp) :: m !< running index surface elements
	725
	726	!
	727	!-- Note, ln(z/z0h) and ln(z/z0q) is also calculated even if neural simulations are applied.
	728	!-- This is because the scalar coefficients are also used for other scalars such as passive scalars,
	729	!-- chemistry and aerosols.
[4717]	730	!$OMP PARALLEL DO PRIVATE( z_mo )
[4593]	731	!$ACC PARALLEL LOOP PRIVATE(z_mo) &
	732	!$ACC PRESENT(surf)
	733	DO m = 1, surf%ns
	734	z_mo = surf%z_mo(m)
	735	surf%ln_z_z0(m) = LOG( z_mo / surf%z0(m) )
	736	surf%ln_z_z0h(m) = LOG( z_mo / surf%z0h(m) )
	737	surf%ln_z_z0q(m) = LOG( z_mo / surf%z0q(m) )
	738	ENDDO
	739
	740	END SUBROUTINE calc_ln
	741
	742	!--------------------------------------------------------------------------------------------------!
	743	! Description:
	744	! ------------
[4562]	745	!> Compute the absolute value of the horizontal velocity (relative to the surface) for horizontal
	746	!> surface elements. This is required by all methods.
	747	!--------------------------------------------------------------------------------------------------!
	748	SUBROUTINE calc_uvw_abs
[1691]	749
[4562]	750	IMPLICIT NONE
[1691]	751
[4562]	752	INTEGER(iwp) :: i !< running index x direction
	753	INTEGER(iwp) :: ibit !< flag to mask computation of relative velocity in case of downward-facing surfaces
	754	INTEGER(iwp) :: j !< running index y direction
	755	INTEGER(iwp) :: k !< running index z direction
	756	INTEGER(iwp) :: m !< running index surface elements
	757
	758	REAL(wp) :: w_lfc !< local free convection velocity scale
[2232]	759	!
[4562]	760	!-- ibit is 1 for upward-facing surfaces, zero for downward-facing surfaces.
	761	ibit = MERGE( 1, 0, .NOT. downward )
[1691]	762
[4593]	763	IF ( use_free_convection_scaling ) THEN
	764	!$OMP PARALLEL DO PRIVATE(i, j, k, w_lfc)
	765	!$ACC PARALLEL LOOP PRIVATE(i, j, k, w_lfc) &
	766	!$ACC PRESENT(surf, u, v)
	767	DO m = 1, surf%ns
	768	i = surf%i(m)
	769	j = surf%j(m)
	770	k = surf%k(m)
[1691]	771
	772	!
[4593]	773	!-- Calculate free convection velocity scale w_lfc is use_free_convection_scaling = .T.. This
	774	!-- will maintain a horizontal velocity even for very weak wind convective conditions. SIGN is
	775	!-- used to set w_lfc to zero under stable conditions.
[4562]	776	w_lfc = ABS(g / surf%pt1(m) * surf%z_mo(m) * surf%shf(m))
	777	w_lfc = ( 0.5_wp * ( w_lfc + SIGN(w_lfc,surf%shf(m)) ) )**(0.33333_wp)
[3157]	778	!
[4593]	779	!-- Compute the absolute value of the horizontal velocity. (relative to the surface in case the
	780	!-- lower surface is the ocean). Please note, in new surface modelling concept the index values
	781	!-- changed, i.e. the reference grid point is not the surface-grid point itself but the first
	782	!-- grid point outside of the topography. Note, in case of coupled ocean-atmosphere simulations
	783	!-- relative velocity with respect to the ocean surface is used, hence, (k-1,j,i) values are used
	784	!-- to calculate the absolute velocity. However, this does not apply for downward-facing walls.
	785	!-- To mask this, use ibit, which checks for upward/downward-facing surfaces.
	786	surf%uvw_abs(m) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) - ( u(k-1,j,i) + u(k-1,j,i+1) &
	787	) * ibit ) &
	788	)**2 &
	789	+ ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) - ( v(k-1,j,i) + v(k-1,j+1,i) &
	790	) * ibit ) &
[4562]	791	)2 + w_lfc2 )
[4593]	792	ENDDO
	793	ELSE
	794	!$OMP PARALLEL DO PRIVATE(i, j, k)
	795	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	796	!$ACC PRESENT(surf, u, v)
	797	DO m = 1, surf%ns
	798	i = surf%i(m)
	799	j = surf%j(m)
	800	k = surf%k(m)
	801	!
	802	!-- Compute the absolute value of the horizontal velocity. (relative to the surface in case the
	803	!-- lower surface is the ocean). Please note, in new surface modelling concept the index values
	804	!-- changed, i.e. the reference grid point is not the surface-grid point itself but the first
	805	!-- grid point outside of the topography. Note, in case of coupled ocean-atmosphere simulations
	806	!-- relative velocity with respect to the ocean surface is used, hence, (k-1,j,i) values are used
	807	!-- to calculate the absolute velocity. However, this does not apply for downward-facing walls.
	808	!-- To mask this, use ibit, which checks for upward/downward-facing surfaces.
	809	surf%uvw_abs(m) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) - ( u(k-1,j,i) + u(k-1,j,i+1) &
	810	) * ibit ) &
	811	)**2 &
	812	+ ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) - ( v(k-1,j,i) + v(k-1,j+1,i) &
	813	) * ibit ) &
	814	)**2 )
	815	ENDDO
	816	ENDIF
[2232]	817
[4562]	818	END SUBROUTINE calc_uvw_abs
[3148]	819
[1691]	820
[4562]	821	!--------------------------------------------------------------------------------------------------!
[2232]	822	! Description:
	823	! ------------
[4562]	824	!> Compute the absolute value of the horizontal velocity (relative to the surface) for horizontal
	825	!> surface elements. This is required by all methods.
	826	!--------------------------------------------------------------------------------------------------!
	827	SUBROUTINE calc_uvw_abs_v_ugrid
[2232]	828
[4562]	829	IMPLICIT NONE
[2232]	830
[4562]	831	INTEGER(iwp) :: i !< running index x direction
	832	INTEGER(iwp) :: j !< running index y direction
	833	INTEGER(iwp) :: k !< running index z direction
	834	INTEGER(iwp) :: m !< running index surface elements
[2232]	835
[4562]	836	REAL(wp) :: u_i !< u-component on xu-grid
	837	REAL(wp) :: w_i !< w-component on xu-grid
[2232]	838
	839
[4562]	840	DO m = 1, surf%ns
	841	i = surf%i(m)
	842	j = surf%j(m)
	843	k = surf%k(m)
[1691]	844	!
[4562]	845	!-- Compute the absolute value of the surface parallel velocity on u-grid.
	846	u_i = u(k,j,i)
	847	w_i = 0.25_wp * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) )
[1691]	848
[4562]	849	surf%uvw_abs(m) = SQRT( u_i2 + w_i2 )
	850	ENDDO
[1709]	851
[4562]	852	END SUBROUTINE calc_uvw_abs_v_ugrid
[1709]	853
[4562]	854	!--------------------------------------------------------------------------------------------------!
[1709]	855	! Description:
	856	! ------------
[4562]	857	!> Compute the absolute value of the horizontal velocity (relative to the surface) for horizontal
	858	!> surface elements. This is required by all methods.
	859	!--------------------------------------------------------------------------------------------------!
	860	SUBROUTINE calc_uvw_abs_v_vgrid
[2232]	861
[4562]	862	IMPLICIT NONE
[2232]	863
[4562]	864	INTEGER(iwp) :: i !< running index x direction
	865	INTEGER(iwp) :: j !< running index y direction
	866	INTEGER(iwp) :: k !< running index z direction
	867	INTEGER(iwp) :: m !< running index surface elements
[2232]	868
[4562]	869	REAL(wp) :: v_i !< v-component on yv-grid
	870	REAL(wp) :: w_i !< w-component on yv-grid
[2232]	871
	872
[4562]	873	DO m = 1, surf%ns
	874	i = surf%i(m)
	875	j = surf%j(m)
	876	k = surf%k(m)
[2232]	877
[4562]	878	v_i = u(k,j,i)
	879	w_i = 0.25_wp * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) )
[2232]	880
[4562]	881	surf%uvw_abs(m) = SQRT( v_i2 + w_i2 )
	882	ENDDO
[2232]	883
[4562]	884	END SUBROUTINE calc_uvw_abs_v_vgrid
[2232]	885
[4562]	886	!--------------------------------------------------------------------------------------------------!
[2232]	887	! Description:
	888	! ------------
[4562]	889	!> Compute the absolute value of the horizontal velocity (relative to the surface) for horizontal
	890	!> surface elements. This is required by all methods.
	891	!--------------------------------------------------------------------------------------------------!
	892	SUBROUTINE calc_uvw_abs_v_wgrid
[2232]	893
[4562]	894	IMPLICIT NONE
[2232]	895
[4562]	896	INTEGER(iwp) :: i !< running index x direction
	897	INTEGER(iwp) :: j !< running index y direction
	898	INTEGER(iwp) :: k !< running index z direction
	899	INTEGER(iwp) :: m !< running index surface elements
[2232]	900
[4562]	901	REAL(wp) :: u_i !< u-component on x-zw-grid
	902	REAL(wp) :: v_i !< v-component on y-zw-grid
	903	REAL(wp) :: w_i !< w-component on zw-grid
[2232]	904	!
[4562]	905	!-- North- (l=0) and south-facing (l=1) surfaces
	906	IF ( l == 0 .OR. l == 1 ) THEN
	907	DO m = 1, surf%ns
	908	i = surf%i(m)
	909	j = surf%j(m)
	910	k = surf%k(m)
[2232]	911
[4562]	912	u_i = 0.25_wp * ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
	913	v_i = 0.0_wp
	914	w_i = w(k,j,i)
[2232]	915
[4562]	916	surf%uvw_abs(m) = SQRT( u_i2 + v_i2 + w_i**2 )
	917	ENDDO
[2232]	918	!
[4562]	919	!-- East- (l=2) and west-facing (l=3) surfaces
	920	ELSE
	921	DO m = 1, surf%ns
	922	i = surf%i(m)
	923	j = surf%j(m)
	924	k = surf%k(m)
[2232]	925
[4562]	926	u_i = 0.0_wp
	927	v_i = 0.25_wp * ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
	928	w_i = w(k,j,i)
[2232]	929
[4562]	930	surf%uvw_abs(m) = SQRT( u_i2 + v_i2 + w_i**2 )
	931	ENDDO
	932	ENDIF
[2232]	933
[4562]	934	END SUBROUTINE calc_uvw_abs_v_wgrid
[2232]	935
[4562]	936	!--------------------------------------------------------------------------------------------------!
[2232]	937	! Description:
	938	! ------------
[4562]	939	!> Compute the absolute value of the horizontal velocity (relative to the surface) for horizontal
	940	!> surface elements. This is required by all methods.
	941	!--------------------------------------------------------------------------------------------------!
	942	SUBROUTINE calc_uvw_abs_v_sgrid
[2232]	943
[4562]	944	IMPLICIT NONE
[2232]	945
[4562]	946	INTEGER(iwp) :: i !< running index x direction
	947	INTEGER(iwp) :: j !< running index y direction
	948	INTEGER(iwp) :: k !< running index z direction
	949	INTEGER(iwp) :: m !< running index surface elements
[2232]	950
[4562]	951	REAL(wp) :: u_i !< u-component on scalar grid
	952	REAL(wp) :: v_i !< v-component on scalar grid
	953	REAL(wp) :: w_i !< w-component on scalar grid
[2232]	954
	955	!
[4562]	956	!-- North- (l=0) and south-facing (l=1) walls
	957	IF ( l == 0 .OR. l == 1 ) THEN
	958	DO m = 1, surf%ns
	959	i = surf%i(m)
	960	j = surf%j(m)
	961	k = surf%k(m)
[2232]	962
[4562]	963	u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
	964	v_i = 0.0_wp
	965	w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
[2232]	966
[4562]	967	surf%uvw_abs(m) = SQRT( u_i2 + v_i2 + w_i**2 )
	968	ENDDO
[2232]	969	!
[4562]	970	!-- East- (l=2) and west-facing (l=3) walls
	971	ELSE
	972	DO m = 1, surf%ns
	973	i = surf%i(m)
	974	j = surf%j(m)
	975	k = surf%k(m)
[2232]	976
[4562]	977	u_i = 0.0_wp
	978	v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
	979	w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
[2232]	980
[4562]	981	surf%uvw_abs(m) = SQRT( u_i2 + v_i2 + w_i**2 )
	982	ENDDO
	983	ENDIF
[2232]	984
[4562]	985	END SUBROUTINE calc_uvw_abs_v_sgrid
[2232]	986
	987
[4562]	988	!--------------------------------------------------------------------------------------------------!
[2232]	989	! Description:
	990	! ------------
[1709]	991	!> Calculate the Obukhov length (L) and Richardson flux number (z/L)
[4562]	992	!--------------------------------------------------------------------------------------------------!
	993	SUBROUTINE calc_ol
[1709]	994
[4562]	995	IMPLICIT NONE
[1709]	996
[4562]	997	INTEGER(iwp) :: iter !< Newton iteration step
	998	INTEGER(iwp) :: m !< loop variable over all horizontal wall elements
[1709]	999
[4562]	1000	LOGICAL, DIMENSION(surf%ns) :: convergence_reached !< convergence switch for vectorization
	1001	!$ACC DECLARE CREATE( convergence_reached )
[4366]	1002
[4562]	1003	REAL(wp) :: f !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
	1004	REAL(wp) :: f_d_ol !< Derivative of f
	1005	REAL(wp) :: ol_l !< Lower bound of L for Newton iteration
	1006	REAL(wp) :: ol_m !< Previous value of L for Newton iteration
	1007	REAL(wp) :: ol_old !< Previous time step value of L
	1008	REAL(wp) :: ol_u !< Upper bound of L for Newton iteration
[1709]	1009
[4562]	1010	REAL(wp), DIMENSION(surf%ns) :: ol_old_vec !< temporary array required for vectorization
	1011	REAL(wp), DIMENSION(surf%ns) :: z_mo_vec !< temporary array required for vectorization
	1012	!$ACC DECLARE CREATE( ol_old_vec, z_mo_vec )
[4366]	1013
[2232]	1014	!
[4562]	1015	!-- Evaluate bulk Richardson number (calculation depends on definition based on setting of boundary
	1016	!-- conditions)
	1017	IF ( ibc_pt_b /= 1 ) THEN
	1018	IF ( humidity ) THEN
	1019	!$OMP PARALLEL DO PRIVATE( z_mo )
	1020	DO m = 1, surf%ns
	1021	z_mo = surf%z_mo(m)
	1022	surf%rib(m) = g * z_mo * ( surf%vpt1(m) - surf%vpt_surface(m) ) / &
	1023	( surf%uvw_abs(m)*2 surf%vpt1(m) + 1.0E-20_wp )
	1024	ENDDO
[3668]	1025	ELSE
[4562]	1026	!$OMP PARALLEL DO PRIVATE( z_mo )
	1027	DO m = 1, surf%ns
	1028	z_mo = surf%z_mo(m)
	1029	surf%rib(m) = g * z_mo * ( surf%pt1(m) - surf%pt_surface(m) ) / &
	1030	( surf%uvw_abs(m)*2 surf%pt1(m) + 1.0E-20_wp )
	1031	ENDDO
	1032	ENDIF
	1033	ELSE
	1034	IF ( humidity ) THEN
	1035	!$OMP PARALLEL DO PRIVATE( k, z_mo )
	1036	DO m = 1, surf%ns
	1037	k = surf%k(m)
	1038	z_mo = surf%z_mo(m)
	1039	surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp * surf%qv1(m) ) * &
	1040	surf%shf(m) + 0.61_wp * surf%pt1(m) * surf%qsws(m) ) * &
	1041	drho_air_zw(k-1) / ( surf%uvw_abs(m)*3 surf%vpt1(m) * kappa**2 &
[1709]	1042	+ 1.0E-20_wp )
[4562]	1043	ENDDO
	1044	ELSE
	1045	!$OMP PARALLEL DO PRIVATE( k, z_mo )
	1046	!$ACC PARALLEL LOOP PRIVATE(k, z_mo) &
	1047	!$ACC PRESENT(surf, drho_air_zw)
	1048	DO m = 1, surf%ns
	1049	k = surf%k(m)
	1050	z_mo = surf%z_mo(m)
	1051	surf%rib(m) = - g * z_mo * surf%shf(m) * drho_air_zw(k-1) / &
	1052	( surf%uvw_abs(m)*3 surf%pt1(m) * kappa**2 + 1.0E-20_wp )
	1053	ENDDO
[1691]	1054	ENDIF
[4562]	1055	ENDIF
[1691]	1056
[4562]	1057	IF ( loop_optimization == 'cache' ) THEN
[1691]	1058	!
[4562]	1059	!-- Calculate the Obukhov length using Newton iteration
[4717]	1060	!$OMP PARALLEL DO PRIVATE(i, j, z_mo, ol_old, iter, ol_m, ol_l, ol_u, f, f_d_ol)
[4562]	1061	!$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
	1062	!$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
	1063	!$ACC PRESENT(surf)
	1064	DO m = 1, surf%ns
	1065	i = surf%i(m)
	1066	j = surf%j(m)
[1691]	1067
[4562]	1068	z_mo = surf%z_mo(m)
[1691]	1069
	1070	!
[4562]	1071	!-- Store current value in case the Newton iteration fails
	1072	ol_old = surf%ol(m)
[1691]	1073
	1074	!
[4562]	1075	!-- Ensure that the bulk Richardson number and the Obukhov length have the same sign
	1076	IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. ABS( surf%ol(m) ) == ol_max ) THEN
	1077	IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp
	1078	IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
	1079	ENDIF
[1691]	1080	!
[4562]	1081	!-- Iteration to find Obukhov length
	1082	iter = 0
	1083	DO
	1084	iter = iter + 1
[4366]	1085	!
[4562]	1086	!-- In case of divergence, use the value of the previous time step
	1087	IF ( iter > 1000 ) THEN
	1088	surf%ol(m) = ol_old
	1089	EXIT
	1090	ENDIF
[4366]	1091
[4562]	1092	ol_m = surf%ol(m)
	1093	ol_l = ol_m - 0.001_wp * ol_m
	1094	ol_u = ol_m + 0.001_wp * ol_m
[4366]	1095
	1096
[4562]	1097	IF ( ibc_pt_b /= 1 ) THEN
[4366]	1098	!
[4562]	1099	!-- Calculate f = Ri - [...]/[...]^2 = 0
[4593]	1100	f = surf%rib(m) - ( z_mo / ol_m ) * ( surf%ln_z_z0h(m) &
[4562]	1101	- psi_h( z_mo / ol_m ) &
	1102	+ psi_h( surf%z0h(m) / ol_m ) ) / &
[4593]	1103	( surf%ln_z_z0(m) - psi_m( z_mo / ol_m ) &
	1104	+ psi_m( surf%z0(m) / ol_m ) )**2
[4366]	1105
	1106	!
[4562]	1107	!-- Calculate df/dL
[4593]	1108	f_d_ol = ( - ( z_mo / ol_u ) * ( surf%ln_z_z0h(m) &
[4562]	1109	- psi_h( z_mo / ol_u ) &
	1110	+ psi_h( surf%z0h(m) / ol_u ) ) / &
[4593]	1111	( surf%ln_z_z0(m) - psi_m( z_mo / ol_u ) &
	1112	+ psi_m( surf%z0(m) / ol_u ) )**2 &
	1113	+ ( z_mo / ol_l ) * ( surf%ln_z_z0h(m) - psi_h( z_mo / ol_l ) &
	1114	+ psi_h( surf%z0h(m) / ol_l ) ) /&
	1115	( surf%ln_z_z0(m) - psi_m( z_mo / ol_l ) &
	1116	+ psi_m( surf%z0(m) / ol_l ) )**2 ) / ( ol_u - ol_l )
[4562]	1117	ELSE
[4366]	1118	!
[4562]	1119	!-- Calculate f = Ri - 1 /[...]^3 = 0
	1120	f = surf%rib(m) - ( z_mo / ol_m ) / &
[4593]	1121	( surf%ln_z_z0(m) - psi_m( z_mo / ol_m ) + psi_m( surf%z0(m) / ol_m ) )**3
[4366]	1122
	1123	!
[4562]	1124	!-- Calculate df/dL
[4593]	1125	f_d_ol = ( - ( z_mo / ol_u ) / ( surf%ln_z_z0(m) &
[4562]	1126	- psi_m( z_mo / ol_u ) &
	1127	+ psi_m( surf%z0(m) / ol_u ) &
	1128	)**3 &
[4593]	1129	+ ( z_mo / ol_l ) / ( surf%ln_z_z0(m) &
[4562]	1130	- psi_m( z_mo / ol_l ) &
	1131	+ psi_m( surf%z0(m) / ol_l ) &
	1132	)**3 &
	1133	) / ( ol_u - ol_l )
	1134	ENDIF
[4366]	1135	!
[4562]	1136	!-- Calculate new L
	1137	surf%ol(m) = ol_m - f / f_d_ol
[4366]	1138
	1139	!
[4562]	1140	!-- Ensure that the bulk Richardson number and the Obukhov length have the same sign and
	1141	!-- ensure convergence.
	1142	IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp
[4366]	1143	!
[4562]	1144	!-- If unrealistic value occurs, set L to the maximum value that is allowed
	1145	IF ( ABS( surf%ol(m) ) > ol_max ) THEN
	1146	surf%ol(m) = ol_max
	1147	EXIT
	1148	ENDIF
[4366]	1149	!
[4562]	1150	!-- Assure that Obukhov length does not become zero. If the limit is reached, exit the loop.
	1151	IF ( ABS( surf%ol(m) ) < 1E-5_wp ) THEN
	1152	surf%ol(m) = SIGN( 1E-5_wp, surf%ol(m) )
	1153	EXIT
	1154	ENDIF
[4366]	1155	!
[4562]	1156	!-- Check for convergence
	1157	IF ( ABS( ( surf%ol(m) - ol_m ) / surf%ol(m) ) < 1.0E-4_wp ) EXIT
[4366]	1158	ENDDO
[4562]	1159	ENDDO
[4366]	1160
	1161	!
[4562]	1162	!-- Vector Version
	1163	ELSE
[4366]	1164	!
[4562]	1165	!-- Calculate the Obukhov length using Newton iteration
	1166	!-- First set arrays required for vectorization
	1167	!$ACC PARALLEL LOOP &
	1168	!$ACC PRESENT(surf)
	1169	DO m = 1, surf%ns
	1170	z_mo_vec(m) = surf%z_mo(m)
[4366]	1171	!
[4562]	1172	!-- Store current value in case the Newton iteration fails
	1173	ol_old_vec(m) = surf%ol(m)
[4366]	1174	!
[4562]	1175	!-- Ensure that the bulk Richardson number and the Obukhov length have the same sign
	1176	IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. ABS( surf%ol(m) ) == ol_max ) THEN
	1177	IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp
	1178	IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
	1179	ENDIF
[4370]	1180	!
[4562]	1181	!-- Initialize convergence flag
	1182	convergence_reached(m) = .FALSE.
	1183	ENDDO
[4366]	1184
	1185	!
[4562]	1186	!-- Iteration to find Obukhov length
	1187	iter = 0
	1188	DO
	1189	iter = iter + 1
[1691]	1190	!
[4562]	1191	!-- In case of divergence, use the value(s) of the previous time step
	1192	IF ( iter > 1000 ) THEN
	1193	!$ACC PARALLEL LOOP &
[4370]	1194	!$ACC PRESENT(surf)
[4366]	1195	DO m = 1, surf%ns
[4562]	1196	IF ( .NOT. convergence_reached(m) ) surf%ol(m) = ol_old_vec(m)
	1197	ENDDO
	1198	EXIT
	1199	ENDIF
[1691]	1200
[4562]	1201	!$ACC PARALLEL LOOP PRIVATE(ol_m, ol_l, ol_u, f, f_d_ol) &
	1202	!$ACC PRESENT(surf)
	1203	DO m = 1, surf%ns
	1204	IF ( convergence_reached(m) ) CYCLE
[4366]	1205
[4562]	1206	ol_m = surf%ol(m)
	1207	ol_l = ol_m - 0.001_wp * ol_m
	1208	ol_u = ol_m + 0.001_wp * ol_m
[4366]	1209
	1210
[4562]	1211	IF ( ibc_pt_b /= 1 ) THEN
[1691]	1212	!
[4562]	1213	!-- Calculate f = Ri - [...]/[...]^2 = 0
[4690]	1214	f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) * ( surf%ln_z_z0h(m) &
[4562]	1215	- psi_h( z_mo_vec(m) / ol_m ) &
	1216	+ psi_h( surf%z0h(m) / ol_m ) &
	1217	) / &
[4593]	1218	( surf%ln_z_z0(m) &
[4562]	1219	- psi_m( z_mo_vec(m) / ol_m ) &
	1220	+ psi_m( surf%z0(m) / ol_m ) &
	1221	)**2
[1691]	1222
	1223	!
[4562]	1224	!-- Calculate df/dL
[4593]	1225	f_d_ol = ( - ( z_mo_vec(m) / ol_u ) * ( surf%ln_z_z0h(m) &
[4562]	1226	- psi_h( z_mo_vec(m) / ol_u ) &
	1227	+ psi_h( surf%z0h(m) / ol_u ) &
	1228	) / &
[4593]	1229	( surf%ln_z_z0(m) &
[4562]	1230	- psi_m( z_mo_vec(m) / ol_u ) &
	1231	+ psi_m( surf%z0(m) / ol_u ) &
	1232	)**2 &
[4593]	1233	+ ( z_mo_vec(m) / ol_l ) * ( surf%ln_z_z0h(m) &
[4562]	1234	- psi_h( z_mo_vec(m) / ol_l ) &
	1235	+ psi_h( surf%z0h(m) / ol_l ) &
	1236	) / &
[4593]	1237	( surf%ln_z_z0(m) &
[4562]	1238	- psi_m( z_mo_vec(m) / ol_l ) &
	1239	+ psi_m( surf%z0(m) / ol_l ) &
	1240	)**2 &
	1241	) / ( ol_u - ol_l )
	1242	ELSE
[1691]	1243	!
[4562]	1244	!-- Calculate f = Ri - 1 /[...]^3 = 0
[4593]	1245	f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) / ( surf%ln_z_z0(m) &
[4562]	1246	- psi_m( z_mo_vec(m) / ol_m ) &
	1247	+ psi_m( surf%z0(m) / ol_m ) &
	1248	)**3
[1691]	1249
	1250	!
[4562]	1251	!-- Calculate df/dL
[4593]	1252	f_d_ol = ( - ( z_mo_vec(m) / ol_u ) / ( surf%ln_z_z0(m) &
[4562]	1253	- psi_m( z_mo_vec(m) / ol_u ) &
	1254	+ psi_m( surf%z0(m) / ol_u ) &
	1255	)**3 &
[4593]	1256	+ ( z_mo_vec(m) / ol_l ) / ( surf%ln_z_z0(m) &
[4562]	1257	- psi_m( z_mo_vec(m) / ol_l ) &
	1258	+ psi_m( surf%z0(m) / ol_l ) &
	1259	)**3 &
	1260	) / ( ol_u - ol_l )
	1261	ENDIF
[1691]	1262	!
[4562]	1263	!-- Calculate new L
	1264	surf%ol(m) = ol_m - f / f_d_ol
[1691]	1265
	1266	!
[4562]	1267	!-- Ensure that the bulk Richardson number and the Obukhov length have the same sign and
	1268	!-- ensure convergence.
	1269	IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp
[4366]	1270
[1691]	1271	!
[4562]	1272	!-- Check for convergence
	1273	!-- This check does not modify surf%ol, therefore this is done first
	1274	IF ( ABS( ( surf%ol(m) - ol_m ) / surf%ol(m) ) < 1.0E-4_wp ) THEN
	1275	convergence_reached(m) = .TRUE.
	1276	ENDIF
[1691]	1277	!
[4562]	1278	!-- If unrealistic value occurs, set L to the maximum allowed value
	1279	IF ( ABS( surf%ol(m) ) > ol_max ) THEN
	1280	surf%ol(m) = ol_max
	1281	convergence_reached(m) = .TRUE.
	1282	ENDIF
	1283	ENDDO
[4186]	1284	!
[4562]	1285	!-- Assure that Obukhov length does not become zero
	1286	!$ACC PARALLEL LOOP &
	1287	!$ACC PRESENT(surf)
	1288	DO m = 1, surf%ns
	1289	IF ( convergence_reached(m) ) CYCLE
	1290	IF ( ABS( surf%ol(m) ) < 1E-5_wp ) THEN
	1291	surf%ol(m) = SIGN( 10E-6_wp, surf%ol(m) )
	1292	convergence_reached(m) = .TRUE.
	1293	ENDIF
	1294	ENDDO
[1691]	1295
[4562]	1296	IF ( ALL( convergence_reached ) ) EXIT
[1691]	1297
[4562]	1298	ENDDO ! End of iteration loop
[4366]	1299
[4562]	1300	ENDIF ! End of vector branch
[4366]	1301
[4562]	1302	END SUBROUTINE calc_ol
[1691]	1303
	1304
[4562]	1305	!--------------------------------------------------------------------------------------------------!
	1306	! Description:
	1307	! ------------
	1308	!> Calculate friction velocity u*.
	1309	!--------------------------------------------------------------------------------------------------!
	1310	SUBROUTINE calc_us
[1691]	1311
[4562]	1312	IMPLICIT NONE
[1691]	1313
[4562]	1314	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
[1691]	1315
	1316	!
[4562]	1317	!-- Compute u* at horizontal surfaces at the scalars' grid points
	1318	IF ( .NOT. surf_vertical ) THEN
[2232]	1319	!
[4562]	1320	!-- Compute u* at upward-facing surfaces
	1321	IF ( .NOT. downward ) THEN
	1322	!$OMP PARALLEL DO PRIVATE( z_mo )
	1323	!$ACC PARALLEL LOOP PRIVATE(z_mo) &
	1324	!$ACC PRESENT(surf)
	1325	DO m = 1, surf%ns
	1326	z_mo = surf%z_mo(m)
[2232]	1327	!
[4562]	1328	!-- Compute u* at the scalars' grid points
[4593]	1329	surf%us(m) = kappa * surf%uvw_abs(m) / ( surf%ln_z_z0(m) &
[4562]	1330	- psi_m( z_mo / surf%ol(m) ) &
	1331	+ psi_m( surf%z0(m) / surf%ol(m) ) )
	1332	ENDDO
[2232]	1333	!
[4562]	1334	!-- Compute u* at downward-facing surfaces. This case, do not consider any stability.
[2232]	1335	ELSE
[4594]	1336	!$OMP PARALLEL DO
[4593]	1337	!$ACC PARALLEL LOOP &
[3634]	1338	!$ACC PRESENT(surf)
[2232]	1339	DO m = 1, surf%ns
[4562]	1340	!
	1341	!-- Compute u* at the scalars' grid points
[4593]	1342	surf%us(m) = kappa * surf%uvw_abs(m) / surf%ln_z_z0(m)
[1691]	1343	ENDDO
[2232]	1344	ENDIF
[4562]	1345	!
	1346	!-- Compute u* at vertical surfaces at the u/v/v grid, respectively.
	1347	!-- No stability is considered in this case.
	1348	ELSE
[4593]	1349	!$OMP PARALLEL DO
	1350	!$ACC PARALLEL LOOP &
[4562]	1351	!$ACC PRESENT(surf)
	1352	DO m = 1, surf%ns
[4593]	1353	surf%us(m) = kappa * surf%uvw_abs(m) / surf%ln_z_z0(m)
[4562]	1354	ENDDO
	1355	ENDIF
[1691]	1356
[4562]	1357	END SUBROUTINE calc_us
[1691]	1358
[4562]	1359	!--------------------------------------------------------------------------------------------------!
	1360	! Description:
	1361	! ------------
	1362	!> Calculate potential temperature, specific humidity, and virtual potential temperature at first
	1363	!> grid level.
	1364	!--------------------------------------------------------------------------------------------------!
	1365	SUBROUTINE calc_pt_q
[1691]	1366
[4562]	1367	IMPLICIT NONE
[1691]	1368
[4562]	1369	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
[2232]	1370
[4562]	1371	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1372	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	1373	!$ACC PRESENT(surf, pt)
	1374	DO m = 1, surf%ns
	1375	i = surf%i(m)
	1376	j = surf%j(m)
	1377	k = surf%k(m)
[2232]	1378
[3634]	1379	#ifndef _OPENACC
[4562]	1380	IF ( bulk_cloud_model ) THEN
	1381	surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
	1382	surf%qv1(m) = q(k,j,i) - ql(k,j,i)
	1383	ELSEIF( cloud_droplets ) THEN
	1384	surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
	1385	surf%qv1(m) = q(k,j,i)
	1386	ELSE
[3634]	1387	#endif
[4562]	1388	surf%pt1(m) = pt(k,j,i)
[3634]	1389	#ifndef _OPENACC
[4562]	1390	IF ( humidity ) THEN
	1391	surf%qv1(m) = q(k,j,i)
	1392	ELSE
[3634]	1393	#endif
[4562]	1394	surf%qv1(m) = 0.0_wp
[3634]	1395	#ifndef _OPENACC
[4562]	1396	ENDIF
	1397	ENDIF
[2232]	1398
[4562]	1399	IF ( humidity ) THEN
	1400	surf%vpt1(m) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
	1401	ENDIF
[3634]	1402	#endif
[4562]	1403	ENDDO
[1691]	1404
[4562]	1405	END SUBROUTINE calc_pt_q
[1691]	1406
[2696]	1407
[4562]	1408	!--------------------------------------------------------------------------------------------------!
	1409	! Description:
	1410	! ------------
	1411	!> Set potential temperature at surface grid level( only for upward-facing surfs ).
	1412	!--------------------------------------------------------------------------------------------------!
	1413	SUBROUTINE calc_pt_surface
[2696]	1414
[4562]	1415	IMPLICIT NONE
[2696]	1416
[4562]	1417	INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
	1418	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
[2696]	1419
[4562]	1420	k_off = surf%koff
	1421	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1422	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	1423	!$ACC PRESENT(surf, pt)
	1424	DO m = 1, surf%ns
	1425	i = surf%i(m)
	1426	j = surf%j(m)
	1427	k = surf%k(m)
	1428	surf%pt_surface(m) = pt(k+k_off,j,i)
	1429	ENDDO
[2696]	1430
[4562]	1431	END SUBROUTINE calc_pt_surface
[2696]	1432
	1433	!
[3152]	1434	!-- Set mixing ratio at surface grid level. ( Only for upward-facing surfs. )
[4562]	1435	SUBROUTINE calc_q_surface
[3152]	1436
[4562]	1437	IMPLICIT NONE
[3152]	1438
[4562]	1439	INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
	1440	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
[3152]	1441
[4562]	1442	k_off = surf%koff
	1443	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1444	DO m = 1, surf%ns
	1445	i = surf%i(m)
	1446	j = surf%j(m)
	1447	k = surf%k(m)
	1448	surf%q_surface(m) = q(k+k_off,j,i)
	1449	ENDDO
[3152]	1450
[4562]	1451	END SUBROUTINE calc_q_surface
[3152]	1452
[4562]	1453	!--------------------------------------------------------------------------------------------------!
	1454	! Description:
	1455	! ------------
	1456	!> Set virtual potential temperature at surface grid level ( only for upward-facing surfs ).
	1457	!--------------------------------------------------------------------------------------------------!
	1458	SUBROUTINE calc_vpt_surface
[3152]	1459
[4562]	1460	IMPLICIT NONE
[3146]	1461
[4562]	1462	INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
	1463	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
[3146]	1464
[4562]	1465	k_off = surf%koff
	1466	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1467	DO m = 1, surf%ns
	1468	i = surf%i(m)
	1469	j = surf%j(m)
	1470	k = surf%k(m)
	1471	surf%vpt_surface(m) = vpt(k+k_off,j,i)
[3146]	1472
[4562]	1473	ENDDO
[3146]	1474
[4562]	1475	END SUBROUTINE calc_vpt_surface
[3146]	1476
[4562]	1477	!--------------------------------------------------------------------------------------------------!
	1478	! Description:
	1479	! ------------
	1480	!> Calculate the other MOST scaling parameters theta, q, (qc, qr, nc, nr)
	1481	!--------------------------------------------------------------------------------------------------!
	1482	SUBROUTINE calc_scaling_parameters
[3146]	1483
[4562]	1484	IMPLICIT NONE
[1691]	1485
	1486
[4562]	1487	INTEGER(iwp) :: lsp !< running index for chemical species
	1488	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
	1489	!
	1490	!-- Compute theta* at horizontal surfaces
	1491	IF ( constant_heatflux .AND. .NOT. surf_vertical ) THEN
	1492	!
	1493	!-- For a given heat flux in the surface layer:
[2232]	1494
[4562]	1495	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1496	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	1497	!$ACC PRESENT(surf, drho_air_zw)
	1498	DO m = 1, surf%ns
	1499	i = surf%i(m)
	1500	j = surf%j(m)
	1501	k = surf%k(m)
	1502	surf%ts(m) = -surf%shf(m) * drho_air_zw(k-1) / ( surf%us(m) + 1E-30_wp )
[1691]	1503	!
[4562]	1504	!-- ts must be limited, because otherwise overflow may occur in case of us=0 when computing
	1505	!-- ol further below
	1506	IF ( surf%ts(m) < -1.05E5_wp ) surf%ts(m) = -1.0E5_wp
	1507	IF ( surf%ts(m) > 1.0E5_wp ) surf%ts(m) = 1.0E5_wp
	1508	ENDDO
[2232]	1509
[4562]	1510	ELSEIF ( .NOT. surf_vertical ) THEN
	1511	!
	1512	!-- For a given surface temperature:
	1513	IF ( large_scale_forcing .AND. lsf_surf ) THEN
	1514
[2232]	1515	!$OMP PARALLEL DO PRIVATE( i, j, k )
[4562]	1516	DO m = 1, surf%ns
	1517	i = surf%i(m)
[2232]	1518	j = surf%j(m)
	1519	k = surf%k(m)
[4562]	1520	pt(k-1,j,i) = pt_surface
	1521	ENDDO
	1522	ENDIF
[2232]	1523
[4562]	1524	!$OMP PARALLEL DO PRIVATE( z_mo )
	1525	DO m = 1, surf%ns
	1526	z_mo = surf%z_mo(m)
	1527	surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) ) &
[4593]	1528	/ ( surf%ln_z_z0h(m) - psi_h( z_mo / surf%ol(m) ) &
	1529	+ psi_h( surf%z0h(m) / surf%ol(m) ) )
[4562]	1530	ENDDO
[2232]	1531
[4562]	1532	ENDIF
[1691]	1533	!
[4562]	1534	!-- Compute theta* at vertical surfaces. This is only required in case of land-surface model, in
	1535	!-- order to compute aerodynamical resistance.
	1536	IF ( surf_vertical ) THEN
	1537	!$OMP PARALLEL DO PRIVATE( i, j )
	1538	DO m = 1, surf%ns
	1539	i = surf%i(m)
	1540	j = surf%j(m)
	1541	surf%ts(m) = -surf%shf(m) / ( surf%us(m) + 1E-30_wp )
	1542	!
	1543	!-- ts must be limited, because otherwise overflow may occur in case of us=0 when computing ol
	1544	!-- further below
	1545	IF ( surf%ts(m) < -1.05E5_wp ) surf%ts(m) = -1.0E5_wp
	1546	IF ( surf%ts(m) > 1.0E5_wp ) surf%ts(m) = 1.0E5_wp
	1547	ENDDO
	1548	ENDIF
[2232]	1549
[4562]	1550	!
	1551	!-- If required compute q* at horizontal surfaces
	1552	IF ( humidity ) THEN
	1553	IF ( constant_waterflux .AND. .NOT. surf_vertical ) THEN
	1554	!
	1555	!-- For a given water flux in the surface layer
	1556	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1557	DO m = 1, surf%ns
	1558	i = surf%i(m)
	1559	j = surf%j(m)
	1560	k = surf%k(m)
	1561	surf%qs(m) = -surf%qsws(m) * drho_air_zw(k-1) / ( surf%us(m) + 1E-30_wp )
[1691]	1562	ENDDO
	1563
[4562]	1564	ELSEIF ( .NOT. surf_vertical ) THEN
	1565	coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. run_coupled )
	1566
[1788]	1567	IF ( large_scale_forcing .AND. lsf_surf ) THEN
[4562]	1568	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1569	DO m = 1, surf%ns
	1570	i = surf%i(m)
	1571	j = surf%j(m)
	1572	k = surf%k(m)
	1573	q(k-1,j,i) = q_surface
[2232]	1574
[4562]	1575	ENDDO
	1576	ENDIF
	1577
	1578	!
	1579	!-- Assume saturation for atmosphere coupled to ocean (but not in case of precursor runs)
	1580	IF ( coupled_run ) THEN
	1581	!$OMP PARALLEL DO PRIVATE( i, j, k, e_s )
	1582	DO m = 1, surf%ns
	1583	i = surf%i(m)
[2232]	1584	j = surf%j(m)
	1585	k = surf%k(m)
[4562]	1586	e_s = 6.1_wp * EXP( 0.07_wp * ( MIN( pt(k-1,j,i), pt(k,j,i) ) - 273.15_wp ) )
	1587	q(k-1,j,i) = rd_d_rv * e_s / ( surface_pressure - e_s )
	1588	ENDDO
	1589	ENDIF
[2232]	1590
[4562]	1591	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
[4717]	1592	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
[4562]	1593	DO m = 1, surf%ns
	1594	i = surf%i(m)
	1595	j = surf%j(m)
	1596	k = surf%k(m)
	1597	z_mo = surf%z_mo(m)
	1598	surf%qs(m) = kappa * ( surf%qv1(m) - surf%q_surface(m) ) &
[4593]	1599	/ ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) ) &
	1600	+ psi_h( surf%z0q(m) / surf%ol(m) ) )
[1691]	1601	ENDDO
[4562]	1602	ELSE
	1603	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	1604	DO m = 1, surf%ns
	1605	i = surf%i(m)
	1606	j = surf%j(m)
	1607	k = surf%k(m)
	1608	z_mo = surf%z_mo(m)
	1609	surf%qs(m) = kappa * ( q(k,j,i) - q(k-1,j,i) ) &
[4593]	1610	/ ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) ) &
	1611	+ psi_h( surf%z0q(m) / surf%ol(m) ) )
[4562]	1612	ENDDO
[1691]	1613	ENDIF
[2232]	1614	ENDIF
[4562]	1615	!
	1616	!-- Compute q* at vertical surfaces
[2232]	1617	IF ( surf_vertical ) THEN
[2281]	1618	!$OMP PARALLEL DO PRIVATE( i, j )
[4562]	1619	DO m = 1, surf%ns
[2232]	1620
[4562]	1621	i = surf%i(m)
	1622	j = surf%j(m)
	1623	surf%qs(m) = -surf%qsws(m) / ( surf%us(m) + 1E-30_wp )
[2232]	1624
[1691]	1625	ENDDO
	1626	ENDIF
[4562]	1627	ENDIF
[1691]	1628
	1629	!
[4562]	1630	!-- If required compute s*
	1631	IF ( passive_scalar ) THEN
[1691]	1632	!
[4562]	1633	!-- At horizontal surfaces
	1634	IF ( constant_scalarflux .AND. .NOT. surf_vertical ) THEN
[1691]	1635	!
[4562]	1636	!-- For a given scalar flux in the surface layer
	1637	!$OMP PARALLEL DO PRIVATE( i, j )
	1638	DO m = 1, surf%ns
	1639	i = surf%i(m)
	1640	j = surf%j(m)
	1641	surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
	1642	ENDDO
	1643	ELSEIF ( .NOT. surf_vertical ) THEN
[1691]	1644
[4562]	1645	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	1646	DO m = 1, surf%ns
	1647	i = surf%i(m)
	1648	j = surf%j(m)
	1649	k = surf%k(m)
	1650	z_mo = surf%z_mo(m)
[1691]	1651
[4562]	1652	surf%ss(m) = kappa * ( s(k,j,i) - s(k-1,j,i) ) &
[4593]	1653	/ ( surf%ln_z_z0h(m) - psi_h( z_mo / surf%ol(m) ) &
	1654	+ psi_h( surf%z0h(m) / surf%ol(m) ) )
[4562]	1655	ENDDO
[1691]	1656	ENDIF
[1960]	1657	!
[4562]	1658	!-- At vertical surfaces
	1659	IF ( surf_vertical ) THEN
	1660	!$OMP PARALLEL DO PRIVATE( i, j )
	1661	DO m = 1, surf%ns
	1662	i = surf%i(m)
	1663	j = surf%j(m)
	1664	surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
	1665	ENDDO
	1666	ENDIF
	1667	ENDIF
	1668
[1960]	1669	!
[4562]	1670	!-- If required compute cs* (chemical species)
	1671	IF ( air_chemistry ) THEN
[2232]	1672	!
[4562]	1673	!-- At horizontal surfaces
	1674	DO lsp = 1, nvar
	1675	IF ( constant_csflux(lsp) .AND. .NOT. surf_vertical ) THEN
	1676	!-- For a given chemical species' flux in the surface layer
[2281]	1677	!$OMP PARALLEL DO PRIVATE( i, j )
[4562]	1678	DO m = 1, surf%ns
	1679	i = surf%i(m)
	1680	j = surf%j(m)
	1681	surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
[1960]	1682	ENDDO
	1683	ENDIF
[4562]	1684	ENDDO
[2232]	1685	!
[4562]	1686	!-- At vertical surfaces
	1687	IF ( surf_vertical ) THEN
	1688	DO lsp = 1, nvar
[2281]	1689	!$OMP PARALLEL DO PRIVATE( i, j )
[4562]	1690	DO m = 1, surf%ns
	1691	i = surf%i(m)
	1692	j = surf%j(m)
	1693	surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
[2232]	1694	ENDDO
[2696]	1695	ENDDO
	1696	ENDIF
[4562]	1697	ENDIF
[2696]	1698
	1699	!
[4562]	1700	!-- If required compute qc* and nc*
	1701	IF ( bulk_cloud_model .AND. microphysics_morrison .AND. .NOT. surf_vertical ) THEN
	1702	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	1703	DO m = 1, surf%ns
	1704	i = surf%i(m)
	1705	j = surf%j(m)
	1706	k = surf%k(m)
[1691]	1707
[4562]	1708	z_mo = surf%z_mo(m)
[2292]	1709
[4562]	1710	surf%qcs(m) = kappa * ( qc(k,j,i) - qc(k-1,j,i) ) &
[4593]	1711	/ ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) ) &
	1712	+ psi_h( surf%z0q(m) / surf%ol(m) ) )
[2292]	1713
[4562]	1714	surf%ncs(m) = kappa * ( nc(k,j,i) - nc(k-1,j,i) ) &
[4593]	1715	/ ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) ) &
	1716	+ psi_h( surf%z0q(m) / surf%ol(m) ) )
[4562]	1717	ENDDO
[2292]	1718
[4562]	1719	ENDIF
[2292]	1720
[1691]	1721	!
[4562]	1722	!-- If required compute qr* and nr*
	1723	IF ( bulk_cloud_model .AND. microphysics_seifert .AND. .NOT. surf_vertical ) THEN
	1724	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	1725	DO m = 1, surf%ns
	1726	i = surf%i(m)
	1727	j = surf%j(m)
	1728	k = surf%k(m)
[1691]	1729
[4562]	1730	z_mo = surf%z_mo(m)
[1691]	1731
[4562]	1732	surf%qrs(m) = kappa * ( qr(k,j,i) - qr(k-1,j,i) ) &
[4593]	1733	/ ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) ) &
	1734	+ psi_h( surf%z0q(m) / surf%ol(m) ) )
[1691]	1735
[4562]	1736	surf%nrs(m) = kappa * ( nr(k,j,i) - nr(k-1,j,i) ) &
[4593]	1737	/ ( surf%ln_z_z0q(m) - psi_h( z_mo / surf%ol(m) ) &
	1738	+ psi_h( surf%z0q(m) / surf%ol(m) ) )
[4562]	1739	ENDDO
[1691]	1740
[4562]	1741	ENDIF
[1691]	1742
[4562]	1743	END SUBROUTINE calc_scaling_parameters
[1691]	1744
	1745
	1746
[4562]	1747	!--------------------------------------------------------------------------------------------------!
	1748	! Description:
	1749	! ------------
	1750	!> Calculate surface fluxes usws, vsws, shf, qsws, (qcsws, qrsws, ncsws, nrsws)
	1751	!--------------------------------------------------------------------------------------------------!
	1752	SUBROUTINE calc_surface_fluxes
[1691]	1753
[4562]	1754	IMPLICIT NONE
[1691]	1755
[4562]	1756	INTEGER(iwp) :: lsp !< running index for chemical species
	1757	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
[1691]	1758
[4562]	1759	REAL(wp) :: dum !< dummy to precalculate logarithm
	1760	REAL(wp) :: flag_u !< flag indicating u-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
	1761	REAL(wp) :: flag_v !< flag indicating v-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
[1691]	1762
[4562]	1763	REAL(wp), DIMENSION(:), ALLOCATABLE :: u_i !< u-component interpolated onto scalar grid point, required for momentum fluxes
	1764	!< at vertical surfaces
	1765	REAL(wp), DIMENSION(:), ALLOCATABLE :: v_i !< v-component interpolated onto scalar grid point, required for momentum fluxes
	1766	!< at vertical surfaces
	1767	REAL(wp), DIMENSION(:), ALLOCATABLE :: w_i !< w-component interpolated onto scalar grid point, required for momentum fluxes
	1768	!< at vertical surfaces
[1691]	1769
	1770	!
[4562]	1771	!-- Calcuate surface fluxes at horizontal walls
	1772	IF ( .NOT. surf_vertical ) THEN
[2232]	1773	!
[4562]	1774	!-- Compute u'w' for the total model domain at upward-facing surfaces. First compute the
	1775	!-- corresponding component of u* and square it.
	1776	IF ( .NOT. downward ) THEN
	1777	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	1778	!$ACC PARALLEL LOOP PRIVATE(i, j, k, z_mo) &
	1779	!$ACC PRESENT(surf, u, rho_air_zw)
	1780	DO m = 1, surf%ns
	1781	i = surf%i(m)
	1782	j = surf%j(m)
	1783	k = surf%k(m)
[1691]	1784
[4562]	1785	z_mo = surf%z_mo(m)
[1691]	1786
[4562]	1787	surf%usws(m) = kappa * ( u(k,j,i) - u(k-1,j,i) ) &
[4593]	1788	/ ( surf%ln_z_z0(m) - psi_m( z_mo / surf%ol(m) ) &
	1789	+ psi_m( surf%z0(m) / surf%ol(m) ) )
[2232]	1790	!
[4562]	1791	!-- Please note, the computation of usws is not fully accurate. Actually a further
	1792	!-- interpolation of us onto the u-grid, where usws is defined, is required. However, this
	1793	!-- is not done as this would require several data transfers between 2D-grid and the
	1794	!-- surf-type. The impact of the missing interpolation is negligible as several tests have
	1795	!-- shown. Same also for ol.
	1796	surf%usws(m) = -surf%usws(m) * surf%us(m) * rho_air_zw(k-1)
	1797	ENDDO
[1691]	1798	!
[4562]	1799	!-- At downward-facing surfaces
	1800	ELSE
[4593]	1801	!$OMP PARALLEL DO PRIVATE( i, j, k )
[4562]	1802	DO m = 1, surf%ns
	1803	i = surf%i(m)
	1804	j = surf%j(m)
	1805	k = surf%k(m)
[1691]	1806
[4593]	1807	surf%usws(m) = kappa * u(k,j,i) / surf%ln_z_z0(m)
[4562]	1808	surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k)
	1809	ENDDO
	1810	ENDIF
[1691]	1811
[2232]	1812	!
[4562]	1813	!-- Compute v'w' for the total model domain. First compute the corresponding component of u* and
	1814	!-- square it.
	1815	!-- Upward-facing surfaces
	1816	IF ( .NOT. downward ) THEN
	1817	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	1818	!$ACC PARALLEL LOOP PRIVATE(i, j, k, z_mo) &
	1819	!$ACC PRESENT(surf, v, rho_air_zw)
	1820	DO m = 1, surf%ns
	1821	i = surf%i(m)
	1822	j = surf%j(m)
	1823	k = surf%k(m)
[1691]	1824
[4562]	1825	z_mo = surf%z_mo(m)
[1691]	1826
[4562]	1827	surf%vsws(m) = kappa * ( v(k,j,i) - v(k-1,j,i) ) &
[4593]	1828	/ ( surf%ln_z_z0(m) - psi_m( z_mo / surf%ol(m) ) &
	1829	+ psi_m( surf%z0(m) / surf%ol(m) ) )
[1691]	1830	!
[4562]	1831	!-- Please note, the computation of vsws is not fully accurate. Actually a further
	1832	!-- interpolation of us onto the v-grid, where vsws is defined, is required. However, this
	1833	!-- is not done as this would require several data transfers between 2D-grid and the
	1834	!-- surf-type. The impact of the missing interpolation is negligible as several tests have
	1835	!-- shown. Same also for ol.
	1836	surf%vsws(m) = -surf%vsws(m) * surf%us(m) * rho_air_zw(k-1)
	1837	ENDDO
[2232]	1838	!
[4562]	1839	!-- Downward-facing surfaces
	1840	ELSE
[4593]	1841	!$OMP PARALLEL DO PRIVATE( i, j, k )
[4562]	1842	DO m = 1, surf%ns
	1843	i = surf%i(m)
	1844	j = surf%j(m)
	1845	k = surf%k(m)
[1691]	1846
[4593]	1847	surf%vsws(m) = kappa * v(k,j,i) / surf%ln_z_z0(m)
[4562]	1848	surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k)
	1849	ENDDO
	1850	ENDIF
[1691]	1851	!
[4593]	1852	!-- Compute the vertical kinematic heat flux. Note, only upward-facing surfaces are considered,
	1853	!-- at downward-facing surfaces the flux is not parametrized with a scaling parameter.
[4562]	1854	IF ( .NOT. constant_heatflux .AND. ( ( time_since_reference_point <= skip_time_do_lsm &
	1855	.AND. simulated_time > 0.0_wp ) .OR. .NOT. land_surface ) .AND. &
	1856	.NOT. urban_surface .AND. .NOT. downward ) THEN
[4593]	1857	!$OMP PARALLEL DO PRIVATE( k )
[4562]	1858	DO m = 1, surf%ns
	1859	k = surf%k(m)
	1860	surf%shf(m) = -surf%ts(m) * surf%us(m) * rho_air_zw(k-1)
	1861	ENDDO
	1862	ENDIF
[2232]	1863	!
[4562]	1864	!-- Compute the vertical water flux
	1865	IF ( .NOT. constant_waterflux .AND. humidity .AND. &
	1866	( ( time_since_reference_point <= skip_time_do_lsm .AND. simulated_time > 0.0_wp ) &
	1867	.OR. .NOT. land_surface ) .AND. .NOT. urban_surface .AND. .NOT. downward ) &
	1868	THEN
[4593]	1869	!$OMP PARALLEL DO PRIVATE( k )
[4562]	1870	DO m = 1, surf%ns
[4594]	1871	k = surf%k(m)
[4562]	1872	surf%qsws(m) = -surf%qs(m) * surf%us(m) * rho_air_zw(k-1)
	1873	ENDDO
	1874	ENDIF
	1875	!
	1876	!-- Compute the vertical scalar flux
	1877	IF ( .NOT. constant_scalarflux .AND. passive_scalar .AND. .NOT. downward ) THEN
[4594]	1878	!$OMP PARALLEL DO PRIVATE( k )
[4562]	1879	DO m = 1, surf%ns
[4594]	1880	k = surf%k(m)
[4593]	1881	surf%ssws(m) = -surf%ss(m) * surf%us(m) * rho_air_zw(k-1)
[4562]	1882	ENDDO
	1883	ENDIF
	1884	!
	1885	!-- Compute the vertical chemical species' flux
	1886	DO lsp = 1, nvar
	1887	IF ( .NOT. constant_csflux(lsp) .AND. air_chemistry .AND. .NOT. downward ) THEN
[4594]	1888	!$OMP PARALLEL DO PRIVATE( k )
[4562]	1889	DO m = 1, surf%ns
[4594]	1890	k = surf%k(m)
[4593]	1891	surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m) * rho_air_zw(k-1)
[2232]	1892	ENDDO
	1893	ENDIF
[4562]	1894	ENDDO
	1895
[2232]	1896	!
[4562]	1897	!-- Compute (turbulent) fluxes of cloud water content and cloud drop conc.
	1898	IF ( bulk_cloud_model .AND. microphysics_morrison .AND. .NOT. downward) THEN
[4594]	1899	!$OMP PARALLEL DO PRIVATE( k )
[4562]	1900	DO m = 1, surf%ns
[4594]	1901	k = surf%k(m)
[4593]	1902	surf%qcsws(m) = -surf%qcs(m) * surf%us(m) * rho_air_zw(k-1)
	1903	surf%ncsws(m) = -surf%ncs(m) * surf%us(m) * rho_air_zw(k-1)
[4562]	1904	ENDDO
	1905	ENDIF
	1906	!
	1907	!-- Compute (turbulent) fluxes of rain water content and rain drop conc.
	1908	IF ( bulk_cloud_model .AND. microphysics_seifert .AND. .NOT. downward) THEN
[4594]	1909	!$OMP PARALLEL DO PRIVATE( k )
[4562]	1910	DO m = 1, surf%ns
[4594]	1911	k = surf%k(m)
[4593]	1912	surf%qrsws(m) = -surf%qrs(m) * surf%us(m) * rho_air_zw(k-1)
	1913	surf%nrsws(m) = -surf%nrs(m) * surf%us(m) * rho_air_zw(k-1)
[2696]	1914	ENDDO
[4562]	1915	ENDIF
[2696]	1916
	1917	!
[4562]	1918	!-- Bottom boundary condition for the TKE.
	1919	IF ( ibc_e_b == 2 ) THEN
	1920	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1921	DO m = 1, surf%ns
	1922	i = surf%i(m)
	1923	j = surf%j(m)
	1924	k = surf%k(m)
[2292]	1925
[4562]	1926	e(k,j,i) = ( surf%us(m) / 0.1_wp )**2
	1927	!
	1928	!-- As a test: cm = 0.4
	1929	! e(k,j,i) = ( us(j,i) / 0.4_wp )**2
	1930	e(k-1,j,i) = e(k,j,i)
[2292]	1931
[4562]	1932	ENDDO
	1933	ENDIF
[2292]	1934	!
[4562]	1935	!-- Calcuate surface fluxes at vertical surfaces. No stability is considered.
[4593]	1936	!-- Further, no density needs to be considered here.
[4562]	1937	ELSE
	1938	!
	1939	!-- Compute usvs l={0,1} and vsus l={2,3}
	1940	IF ( mom_uv ) THEN
	1941	!
	1942	!-- Generalize computation by introducing flags. At north- and south-facing surfaces
	1943	!-- u-component is used, at east- and west-facing surfaces v-component is used.
	1944	flag_u = MERGE( 1.0_wp, 0.0_wp, l == 0 .OR. l == 1 )
	1945	flag_v = MERGE( 1.0_wp, 0.0_wp, l == 2 .OR. l == 3 )
[4593]	1946	!$OMP PARALLEL DO PRIVATE( i, j, k )
[4562]	1947	DO m = 1, surf%ns
	1948	i = surf%i(m)
	1949	j = surf%j(m)
	1950	k = surf%k(m)
[2232]	1951
[4562]	1952	surf%mom_flux_uv(m) = kappa * ( flag_u * u(k,j,i) + flag_v * v(k,j,i) ) / &
[4593]	1953	surf%ln_z_z0(m)
[1691]	1954
[4562]	1955	surf%mom_flux_uv(m) = - surf%mom_flux_uv(m) * surf%us(m)
	1956	ENDDO
	1957	ENDIF
[1960]	1958	!
[4562]	1959	!-- Compute wsus l={0,1} and wsvs l={2,3}
	1960	IF ( mom_w ) THEN
[4717]	1961	!$OMP PARALLEL DO PRIVATE( i, j, k )
[4562]	1962	DO m = 1, surf%ns
	1963	i = surf%i(m)
	1964	j = surf%j(m)
	1965	k = surf%k(m)
[2232]	1966
[4593]	1967	surf%mom_flux_w(m) = kappa * w(k,j,i) / surf%ln_z_z0(m)
[2232]	1968
[4562]	1969	surf%mom_flux_w(m) = - surf%mom_flux_w(m) * surf%us(m)
	1970	ENDDO
	1971	ENDIF
[2232]	1972	!
[4562]	1973	!-- Compute momentum fluxes used for subgrid-scale TKE production at vertical surfaces. In
	1974	!-- constrast to the calculated momentum fluxes at vertical surfaces before, which are defined on
	1975	!-- the u/v/w-grid, respectively), the TKE fluxes are defined at the scalar grid.
	1976	!--
	1977	IF ( mom_tke ) THEN
[2232]	1978	!
[4562]	1979	!-- Precalculate velocity components at scalar grid point.
	1980	ALLOCATE( u_i(1:surf%ns) )
	1981	ALLOCATE( v_i(1:surf%ns) )
	1982	ALLOCATE( w_i(1:surf%ns) )
	1983
	1984	IF ( l == 0 .OR. l == 1 ) THEN
	1985	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1986	DO m = 1, surf%ns
	1987	i = surf%i(m)
[2232]	1988	j = surf%j(m)
	1989	k = surf%k(m)
[1691]	1990
[4562]	1991	u_i(m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
	1992	v_i(m) = 0.0_wp
	1993	w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
[2232]	1994	ENDDO
[4562]	1995	ELSE
	1996	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1997	DO m = 1, surf%ns
	1998	i = surf%i(m)
[2232]	1999	j = surf%j(m)
	2000	k = surf%k(m)
	2001
[4562]	2002	u_i(m) = 0.0_wp
	2003	v_i(m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
	2004	w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
[1691]	2005	ENDDO
[2232]	2006	ENDIF
[1691]	2007
[4593]	2008	!$OMP PARALLEL DO PRIVATE( i, j, dum )
[4562]	2009	DO m = 1, surf%ns
	2010	i = surf%i(m)
	2011	j = surf%j(m)
[2232]	2012
[4593]	2013	dum = kappa / surf%ln_z_z0(m)
[1691]	2014	!
[4562]	2015	!-- usvs (l=0,1) and vsus (l=2,3)
	2016	surf%mom_flux_tke(0,m) = dum * ( u_i(m) + v_i(m) )
[1691]	2017	!
[4562]	2018	!-- wsvs (l=0,1) and wsus (l=2,3)
	2019	surf%mom_flux_tke(1,m) = dum * w_i(m)
[2232]	2020
[4562]	2021	surf%mom_flux_tke(0:1,m) = - surf%mom_flux_tke(0:1,m) * surf%us(m)
	2022	ENDDO
[2232]	2023	!
[4562]	2024	!-- Deallocate temporary arrays
	2025	DEALLOCATE( u_i )
	2026	DEALLOCATE( v_i )
	2027	DEALLOCATE( w_i )
[1691]	2028	ENDIF
[4562]	2029	ENDIF
[1691]	2030
[4562]	2031	END SUBROUTINE calc_surface_fluxes
[1691]	2032
[4562]	2033
	2034	!--------------------------------------------------------------------------------------------------!
[3597]	2035	! Description:
	2036	! ------------
[4562]	2037	!> Calculates temperature near surface (10 cm) for indoor model or 2 m temperature for output.
	2038	!--------------------------------------------------------------------------------------------------!
	2039	SUBROUTINE calc_pt_near_surface ( z_char )
[1691]	2040
[4562]	2041	IMPLICIT NONE
[3597]	2042
[4562]	2043	CHARACTER(LEN = *), INTENT(IN) :: z_char !< string identifier to identify z level
[3597]	2044
[4562]	2045	INTEGER(iwp) :: m !< running index for surface elements
[3597]	2046
[4562]	2047	SELECT CASE ( z_char)
[3597]	2048
[4562]	2049	CASE ( '10cm' )
[4750]	2050	!
	2051	!-- For horizontal upward-facing surfaces 10-cm temperature can be calculated using MOST.
	2052	IF ( .NOT. downward .AND. .NOT. surf_vertical ) THEN
	2053	DO m = 1, surf%ns
	2054	surf%pt_10cm(m) = surf%pt_surface(m) + surf%ts(m) / kappa &
[4562]	2055	* ( LOG( 0.1_wp / surf%z0h(m) ) - psi_h( 0.1_wp / surf%ol(m) ) &
	2056	+ psi_h( surf%z0h(m) / surf%ol(m) ) )
[4750]	2057	ENDDO
	2058	!
	2059	!-- At vertical surfaces 10-cm temperature cannot be calculated via MOST as the Obukhov length
	2060	!-- and temperature scaling parameter are not calculated. Hence, set 10-cm temperature to
	2061	!-- the grid-cell temperature.
	2062	ELSE
	2063	DO m = 1, surf%ns
	2064	surf%pt_10cm(m) = pt(surf%k(m)+surf%koff,surf%j(m)+surf%joff,surf%i(m)+surf%ioff)
	2065	ENDDO
	2066	ENDIF
[3597]	2067
[4562]	2068	END SELECT
[1691]	2069
[4562]	2070	END SUBROUTINE calc_pt_near_surface
[1691]	2071
	2072
[4562]	2073	!--------------------------------------------------------------------------------------------------!
	2074	! Description:
	2075	! ------------
	2076	!> Integrated stability function for momentum.
	2077	!--------------------------------------------------------------------------------------------------!
	2078	FUNCTION psi_m( zeta )
	2079	!$ACC ROUTINE SEQ
[1691]	2080
[4562]	2081	USE kinds
[1691]	2082
[4562]	2083	IMPLICIT NONE
[1691]	2084
[4562]	2085	REAL(wp) :: psi_m !< Integrated similarity function result
	2086	REAL(wp) :: zeta !< Stability parameter z/L
	2087	REAL(wp) :: x !< dummy variable
[1691]	2088
[4562]	2089	REAL(wp), PARAMETER :: a = 1.0_wp !< constant
	2090	REAL(wp), PARAMETER :: b = 0.66666666666_wp !< constant
	2091	REAL(wp), PARAMETER :: c = 5.0_wp !< constant
	2092	REAL(wp), PARAMETER :: d = 0.35_wp !< constant
	2093	REAL(wp), PARAMETER :: c_d_d = c / d !< constant
	2094	REAL(wp), PARAMETER :: bc_d_d = b * c / d !< constant
[1691]	2095
	2096
[4562]	2097	IF ( zeta < 0.0_wp ) THEN
	2098	x = SQRT( SQRT( 1.0_wp - 16.0_wp * zeta ) )
	2099	psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2 &
	2100	* ( 1.0_wp + x*2 ) 0.125_wp )
	2101	ELSE
[1691]	2102
[4562]	2103	psi_m = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - a * zeta - bc_d_d
[1691]	2104	!
[4562]	2105	!-- Old version for stable conditions (only valid for z/L < 0.5) psi_m = - 5.0_wp * zeta
[1691]	2106
[4562]	2107	ENDIF
[1691]	2108
[4562]	2109	END FUNCTION psi_m
[1691]	2110
	2111
[4562]	2112	!--------------------------------------------------------------------------------------------------!
	2113	! Description:
	2114	!------------
	2115	!> Integrated stability function for heat and moisture.
	2116	!--------------------------------------------------------------------------------------------------!
	2117	FUNCTION psi_h( zeta )
	2118	!$ACC ROUTINE SEQ
[1691]	2119
[4562]	2120	USE kinds
	2121
	2122	IMPLICIT NONE
	2123
	2124	REAL(wp) :: psi_h !< Integrated similarity function result
	2125	REAL(wp) :: zeta !< Stability parameter z/L
	2126	REAL(wp) :: x !< dummy variable
	2127
	2128	REAL(wp), PARAMETER :: a = 1.0_wp !< constant
	2129	REAL(wp), PARAMETER :: b = 0.66666666666_wp !< constant
	2130	REAL(wp), PARAMETER :: c = 5.0_wp !< constant
	2131	REAL(wp), PARAMETER :: d = 0.35_wp !< constant
	2132	REAL(wp), PARAMETER :: c_d_d = c / d !< constant
	2133	REAL(wp), PARAMETER :: bc_d_d = b * c / d !< constant
	2134
	2135
	2136	IF ( zeta < 0.0_wp ) THEN
	2137	x = SQRT( 1.0_wp - 16.0_wp * zeta )
	2138	psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
	2139	ELSE
	2140	psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp &
	2141	+ 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d + 1.0_wp
[1691]	2142	!
[4562]	2143	!-- Old version for stable conditions (only valid for z/L < 0.5)
	2144	!-- psi_h = - 5.0_wp * zeta
	2145	ENDIF
[1691]	2146
[4562]	2147	END FUNCTION psi_h
[1691]	2148
[3130]	2149
[4562]	2150	!--------------------------------------------------------------------------------------------------!
[3130]	2151	! Description:
	2152	! ------------
	2153	!> Calculates stability function for momentum
	2154	!>
	2155	!> @author Hauke Wurps
[4562]	2156	!--------------------------------------------------------------------------------------------------!
	2157	FUNCTION phi_m( zeta )
	2158	!$ACC ROUTINE SEQ
[3130]	2159
[4562]	2160	IMPLICIT NONE
	2161
	2162	REAL(wp) :: phi_m !< Value of the function
	2163	REAL(wp) :: zeta !< Stability parameter z/L
	2164
	2165	REAL(wp), PARAMETER :: a = 16.0_wp !< constant
	2166	REAL(wp), PARAMETER :: c = 5.0_wp !< constant
	2167
	2168	IF ( zeta < 0.0_wp ) THEN
	2169	phi_m = 1.0_wp / SQRT( SQRT( 1.0_wp - a * zeta ) )
	2170	ELSE
	2171	phi_m = 1.0_wp + c * zeta
	2172	ENDIF
	2173
	2174	END FUNCTION phi_m
	2175
[1697]	2176	END MODULE surface_layer_fluxes_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |