Home

Context Navigation

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

Last change on this file since 4464 was 4370, checked in by raasch, 5 years ago
bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added
Property svn:keywords set to `Id`
File size: 78.8 KB

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