Home

Context Navigation

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

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