Home

Context Navigation

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

Last change on this file since 4366 was 4366, checked in by raasch, 5 years ago
code vectorization for NEC Aurora: vectorized version of Temperton FFT, vectorization of Newtor iteration for calculating the Obukhov length
Property svn:keywords set to `Id`
File size: 78.2 KB

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