Home

Context Navigation

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

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