Home

Context Navigation

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

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