Home

Context Navigation

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

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