Home

Context Navigation

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

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