Home

Context Navigation

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

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