Home

Context Navigation

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

Last change on this file since 4574 was 4562, checked in by raasch, 5 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 78.1 KB

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