Home

Context Navigation

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

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