Home

Context Navigation

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

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