Home

Context Navigation

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

Last change on this file since 3661 was 3655, checked in by knoop, 6 years ago
Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update
Property svn:keywords set to `Id`
File size: 91.0 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	!
[3147]	24	!
[1692]	25	! Former revisions:
	26	! -----------------
	27	! $Id: surface_layer_fluxes_mod.f90 3655 2019-01-07 16:51:22Z suehring $
[3634]	28	! OpenACC port for SPEC
	29	!
	30	! 3597 2018-12-04 08:40:18Z maronga
[3597]	31	! Added routine for calculating near surface air potential temperature (moved
	32	! from urban_surface_mod)
	33	!
	34	! 3547 2018-11-21 13:21:24Z suehring
[3547]	35	! variables documented
	36	!
	37	! 3361 2018-10-16 20:39:37Z knoop
[3274]	38	! Modularization of all bulk cloud physics code components
	39	!
	40	! 3271 2018-09-24 08:20:34Z suehring
[3271]	41	! Comment revised
	42	!
	43	! 3157 2018-07-19 21:08:49Z maronga
[3157]	44	! Added local free convection velocity scale w_lfc in calculation of uvw_abs.
	45	! This can be switche on/off by the user via the flag namelist parameter
	46	! use_free_convection_scaling
	47	!
	48	! 3152 2018-07-19 13:26:52Z suehring
[3152]	49	! q_surface is now part of surface structure
	50	!
	51	! 3149 2018-07-19 05:48:09Z maronga
[3149]	52	! Correct accidental last commit
	53	!
	54	! 3148 2018-07-19 05:45:25Z maronga
[3147]	55	! Major bugfix in calculation of Obukhov length
	56	!
	57	! 3146 2018-07-18 22:36:19Z maronga
[3130]	58	! move phi_m from turbulence_closure_mod
	59	!
	60	! 3045 2018-05-28 07:55:41Z Giersch
[3045]	61	! Error message revised
	62	!
	63	! 2766 2018-01-22 17:17:47Z kanani
[2766]	64	! Removed preprocessor directive __chem
	65	!
	66	! 2718 2018-01-02 08:49:38Z maronga
[2716]	67	! Corrected "Former revisions" section
	68	!
	69	! 2696 2017-12-14 17:12:51Z kanani
	70	! - Change in file header (GPL part)
[2696]	71	! - Implementation of chemistry module (FK)
	72	! - Added calculation of pt1 and qv1 for all surface types. Added calculation of
	73	! pt_surface for default-type surfaces (BM)
	74	! - Add flag to disable computation of qsws in case of urban surface (MS)
	75	!
	76	! 2547 2017-10-16 12:41:56Z schwenkel
[2547]	77	! extended by cloud_droplets option
	78	!
	79	! 2321 2017-07-24 15:57:07Z schwenkel
[2321]	80	! Bugfix: Correct index in lookup table for Obukhov length
	81	!
	82	! 2299 2017-06-29 10:14:38Z suehring
[2299]	83	! Adjusted for allow separate spinups of LSM and atmosphere code
	84	!
	85	! 2292 2017-06-20 09:51:42Z schwenkel
[2292]	86	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
	87	! includes two more prognostic equations for cloud drop concentration (nc)
	88	! and cloud water content (qc).
	89	!
	90	! 2281 2017-06-13 11:34:50Z suehring
[2281]	91	! Clean-up unnecessary index access to surface type
	92	!
	93	! 2233 2017-05-30 18:08:54Z suehring
[1692]	94	!
[2233]	95	! 2232 2017-05-30 17:47:52Z suehring
	96	! Adjustments to new surface concept
	97	! OpenMP bugfix
	98	!
[2119]	99	! 2118 2017-01-17 16:38:49Z raasch
	100	! OpenACC directives and related code removed
	101	!
[2092]	102	! 2091 2016-12-21 16:38:18Z suehring
	103	! Bugfix in calculation of vsws ( incorrect linear interpolation of us )
	104	!
[2077]	105	! 2076 2016-12-02 13:54:20Z raasch
	106	! further openmp bugfix for lookup method
	107	!
[2074]	108	! 2073 2016-11-30 14:34:05Z raasch
	109	! openmp bugfix for lookup method
	110	!
[2038]	111	! 2037 2016-10-26 11:15:40Z knoop
	112	! Anelastic approximation implemented
	113	!
[2012]	114	! 2011 2016-09-19 17:29:57Z kanani
	115	! Flag urban_surface is now defined in module control_parameters.
	116	!
[2008]	117	! 2007 2016-08-24 15:47:17Z kanani
	118	! Account for urban surface model in computation of vertical kinematic heatflux
	119	!
[2001]	120	! 2000 2016-08-20 18:09:15Z knoop
	121	! Forced header and separation lines into 80 columns
	122	!
[1993]	123	! 1992 2016-08-12 15:14:59Z suehring
	124	! Minor bug, declaration of look-up index as INTEGER
	125	!
[1961]	126	! 1960 2016-07-12 16:34:24Z suehring
	127	! Treat humidity and passive scalar separately
	128	!
[1930]	129	! 1929 2016-06-09 16:25:25Z suehring
	130	! Bugfix: avoid segmentation fault in case one grid point is horizontally
	131	! completely surrounded by topography
	132	!
[1921]	133	! 1920 2016-05-30 10:50:15Z suehring
	134	! Avoid segmentation fault (see change in 1915) by different initialization of
	135	! us instead of adding a very small number in the denominator
	136	!
[1916]	137	! 1915 2016-05-27 11:05:02Z suehring
	138	! Bugfix: avoid segmentation fault in case of most_method = 'circular' at first
	139	! timestep
	140	!
[1851]	141	! 1850 2016-04-08 13:29:27Z maronga
	142	! Module renamed
	143	!
	144	!
[1823]	145	! 1822 2016-04-07 07:49:42Z hoffmann
	146	! icloud_scheme replaced by microphysics_*
	147	!
[1789]	148	! 1788 2016-03-10 11:01:04Z maronga
	149	! Added parameter z0q which replaces z0h in the similarity functions for
	150	! humidity.
	151	! Syntax layout improved.
	152	!
[1758]	153	! 1757 2016-02-22 15:49:32Z maronga
	154	! Minor fixes.
	155	!
[1750]	156	! 1749 2016-02-09 12:19:56Z raasch
	157	! further OpenACC adjustments
	158	!
[1748]	159	! 1747 2016-02-08 12:25:53Z raasch
	160	! adjustments for OpenACC usage
	161	!
[1710]	162	! 1709 2015-11-04 14:47:01Z maronga
	163	! Bugfix: division by zero could occur when calculating rib at low wind speeds
	164	! Bugfix: calculation of uv_total for neutral = .T., initial value for ol for
	165	! neutral = .T.
	166	!
[1706]	167	! 1705 2015-11-02 14:28:56Z maronga
	168	! Typo removed
	169	!
[1698]	170	! 1697 2015-10-28 17:14:10Z raasch
	171	! FORTRAN and OpenMP errors removed
	172	!
[1697]	173	! 1696 2015-10-27 10:03:34Z maronga
[1691]	174	! Modularized and completely re-written version of prandtl_fluxes.f90. In the
	175	! course of the re-writing two additional methods have been implemented. See
	176	! updated description.
	177	!
	178	! 1551 2015-03-03 14:18:16Z maronga
	179	! Removed land surface model part. The surface fluxes are now always calculated
	180	! within prandtl_fluxes, based on the given surface temperature/humidity (which
	181	! is either provided by the land surface model, by large scale forcing data, or
	182	! directly prescribed by the user.
	183	!
	184	! 1496 2014-12-02 17:25:50Z maronga
	185	! Adapted for land surface model
	186	!
	187	! 1494 2014-11-21 17:14:03Z maronga
	188	! Bugfixes: qs is now calculated before calculation of Rif. calculation of
	189	! buoyancy flux in Rif corrected (added missing humidity term), allow use of
	190	! topography for coupled runs (not tested)
	191	!
	192	! 1361 2014-04-16 15:17:48Z hoffmann
	193	! Bugfix: calculation of turbulent fluxes of rain water content (qrsws) and rain
	194	! drop concentration (nrsws) added
	195	!
	196	! 1340 2014-03-25 19:45:13Z kanani
	197	! REAL constants defined as wp-kind
	198	!
	199	! 1320 2014-03-20 08:40:49Z raasch
	200	! ONLY-attribute added to USE-statements,
	201	! kind-parameters added to all INTEGER and REAL declaration statements,
	202	! kinds are defined in new module kinds,
	203	! old module precision_kind is removed,
	204	! revision history before 2012 removed,
	205	! comment fields (!:) to be used for variable explanations added to
	206	! all variable declaration statements
	207	!
	208	! 1276 2014-01-15 13:40:41Z heinze
	209	! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
	210	!
	211	! 1257 2013-11-08 15:18:40Z raasch
	212	! openACC "kernels do" replaced by "kernels loop", "loop independent" added
	213	!
	214	! 1036 2012-10-22 13:43:42Z raasch
	215	! code put under GPL (PALM 3.9)
	216	!
	217	! 1015 2012-09-27 09:23:24Z raasch
	218	! OpenACC statements added
	219	!
	220	! 978 2012-08-09 08:28:32Z fricke
	221	! roughness length for scalar quantities z0h added
	222	!
	223	! Revision 1.1 1998/01/23 10:06:06 raasch
	224	! Initial revision
	225	!
	226	!
	227	! Description:
	228	! ------------
	229	!> Diagnostic computation of vertical fluxes in the constant flux layer from the
	230	!> values of the variables at grid point k=1. Three different methods are
	231	!> available:
	232	!> 1) the "old" version (most_method = 'circular') which is fast, but inaccurate
	233	!> 2) a Newton iteration method (most_method = 'newton'), which is accurate, but
	234	!> slower
	235	!> 3) a method using a lookup table which is fast and accurate. Note, however,
	236	!> that this method cannot be used in case of roughness heterogeneity
	237	!>
	238	!> @todo (re)move large_scale_forcing actions
[2118]	239	!> @todo check/optimize OpenMP directives
[2696]	240	!> @todo simplify if conditions (which flux need to be computed in which case)
[1691]	241	!------------------------------------------------------------------------------!
	242	MODULE surface_layer_fluxes_mod
	243
	244	USE arrays_3d, &
[2292]	245	ONLY: e, kh, nc, nr, pt, q, ql, qc, qr, s, u, v, vpt, w, zu, zw, &
[3274]	246	drho_air_zw, rho_air_zw, d_exner
[1691]	247
[3274]	248	USE basic_constants_and_equations_mod, &
[3361]	249	ONLY: g, kappa, lv_d_cp, pi, rd_d_rv
[3274]	250
[2696]	251	USE chem_modules, &
	252	ONLY: constant_csflux, nvar
	253
[1691]	254	USE cpulog
	255
	256	USE control_parameters, &
[3274]	257	ONLY: air_chemistry, cloud_droplets, &
	258	constant_heatflux, constant_scalarflux, &
[3597]	259	constant_waterflux, coupling_mode, do_output_at_2m, humidity, &
	260	ibc_e_b, ibc_pt_b, indoor_model, initializing_actions, &
[2232]	261	intermediate_timestep_count, intermediate_timestep_count_max, &
[3274]	262	land_surface, large_scale_forcing, lsf_surf, message_string, &
[2292]	263	most_method, neutral, passive_scalar, pt_surface, q_surface, &
	264	run_coupled, surface_pressure, simulated_time, terminate_run, &
[3157]	265	time_since_reference_point, urban_surface, &
	266	use_free_convection_scaling, zeta_max, zeta_min
[1691]	267
[2232]	268	USE grid_variables, &
	269	ONLY: dx, dy
	270
[1691]	271	USE indices, &
[2232]	272	ONLY: nxl, nxr, nys, nyn, nzb
[1691]	273
	274	USE kinds
	275
[3274]	276	USE bulk_cloud_model_mod, &
	277	ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
	278
[1691]	279	USE pegrid
	280
	281	USE land_surface_model_mod, &
[2232]	282	ONLY: aero_resist_kray, skip_time_do_lsm
[2011]	283
[2232]	284	USE surface_mod, &
	285	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_type, &
	286	surf_usm_h, surf_usm_v
[2007]	287
[1691]	288
	289	IMPLICIT NONE
	290
[1992]	291	INTEGER(iwp) :: i !< loop index x direction
	292	INTEGER(iwp) :: j !< loop index y direction
	293	INTEGER(iwp) :: k !< loop index z direction
[2232]	294	INTEGER(iwp) :: l !< loop index for surf type
	295	INTEGER(iwp) :: li_bnd = 7500 !< Lookup table index of the last time step
[1691]	296
[2232]	297	INTEGER(iwp), PARAMETER :: num_steps = 15000 !< number of steps in the lookup table
[1691]	298
[2232]	299	LOGICAL :: coupled_run !< Flag for coupled atmosphere-ocean runs
	300	LOGICAL :: downward = .FALSE.!< Flag indicating downward-facing horizontal surface
	301	LOGICAL :: mom_uv = .FALSE. !< Flag indicating calculation of usvs and vsus at vertical surfaces
	302	LOGICAL :: mom_w = .FALSE. !< Flag indicating calculation of wsus and wsvs at vertical surfaces
	303	LOGICAL :: mom_tke = .FALSE. !< Flag indicating calculation of momentum fluxes at vertical surfaces used for TKE production
	304	LOGICAL :: surf_vertical !< Flag indicating vertical surfaces
[1691]	305
	306	REAL(wp), DIMENSION(0:num_steps-1) :: rib_tab, & !< Lookup table bulk Richardson number
	307	ol_tab !< Lookup table values of L
	308
	309	REAL(wp) :: e_s, & !< Saturation water vapor pressure
	310	ol_max = 1.0E6_wp, & !< Maximum Obukhov length
	311	rib_max, & !< Maximum Richardson number in lookup table
	312	rib_min, & !< Minimum Richardson number in lookup table
	313	z_mo !< Height of the constant flux layer where MOST is assumed
	314
[2232]	315	TYPE(surf_type), POINTER :: surf !< surf-type array, used to generalize subroutines
[1691]	316
[2232]	317
[1691]	318	SAVE
	319
	320	PRIVATE
	321
[3130]	322	PUBLIC init_surface_layer_fluxes, phi_m, surface_layer_fluxes
[1691]	323
	324	INTERFACE init_surface_layer_fluxes
	325	MODULE PROCEDURE init_surface_layer_fluxes
	326	END INTERFACE init_surface_layer_fluxes
	327
[3130]	328	INTERFACE phi_m
	329	MODULE PROCEDURE phi_m
	330	END INTERFACE phi_m
	331
[1691]	332	INTERFACE surface_layer_fluxes
	333	MODULE PROCEDURE surface_layer_fluxes
	334	END INTERFACE surface_layer_fluxes
	335
	336
	337	CONTAINS
	338
	339
	340	!------------------------------------------------------------------------------!
	341	! Description:
	342	! ------------
	343	!> Main routine to compute the surface fluxes
	344	!------------------------------------------------------------------------------!
	345	SUBROUTINE surface_layer_fluxes
	346
	347	IMPLICIT NONE
	348
[3547]	349	surf_vertical = .FALSE. !< flag indicating vertically orientated surface elements
	350	downward = .FALSE. !< flag indicating downward-facing surface elements
[1691]	351	!
[2696]	352	!-- Derive potential temperature and specific humidity at first grid level
	353	!-- from the fields pt and q
[2232]	354	!
[2696]	355	!-- First call for horizontal default-type surfaces (l=0 - upward facing,
	356	!-- l=1 - downward facing)
	357	DO l = 0, 1
	358	IF ( surf_def_h(l)%ns >= 1 ) THEN
	359	surf => surf_def_h(l)
	360	CALL calc_pt_q
[3146]	361	IF ( .NOT. neutral ) THEN
	362	CALL calc_pt_surface
	363	IF ( humidity ) THEN
[3152]	364	CALL calc_q_surface
[3146]	365	CALL calc_vpt_surface
	366	ENDIF
	367	ENDIF
[2696]	368	ENDIF
	369	ENDDO
[2232]	370	!
[2696]	371	!-- Call for natural-type horizontal surfaces
	372	IF ( surf_lsm_h%ns >= 1 ) THEN
	373	surf => surf_lsm_h
	374	CALL calc_pt_q
	375	ENDIF
	376
	377	!
	378	!-- Call for urban-type horizontal surfaces
	379	IF ( surf_usm_h%ns >= 1 ) THEN
	380	surf => surf_usm_h
	381	CALL calc_pt_q
	382	ENDIF
	383
	384	!
	385	!-- Call for natural-type vertical surfaces
	386	DO l = 0, 3
	387	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	388	surf => surf_lsm_v(l)
[2232]	389	CALL calc_pt_q
	390	ENDIF
[2696]	391
[3146]	392	!-- Call for urban-type vertical surfaces
[2696]	393	IF ( surf_usm_v(l)%ns >= 1 ) THEN
	394	surf => surf_usm_v(l)
[2232]	395	CALL calc_pt_q
	396	ENDIF
[2696]	397	ENDDO
[1691]	398
	399	!
	400	!-- First, calculate the new Obukhov length, then new friction velocity,
	401	!-- followed by the new scaling parameters (th, q, etc.), and the new
	402	!-- surface fluxes if required. The old routine ("circular") requires a
	403	!-- different order of calls as the scaling parameters from the previous time
	404	!-- steps are used to calculate the Obukhov length
	405
	406	!
	407	!-- Depending on setting of most_method use the "old" routine
[2232]	408	!-- Note, each routine is called for different surface types.
	409	!-- First call for default-type horizontal surfaces, for natural- and
	410	!-- urban-type surfaces. Note, at this place only upward-facing horizontal
	411	!-- surfaces are treted.
[1691]	412	IF ( most_method == 'circular' ) THEN
[2232]	413	!
	414	!-- Default-type upward-facing horizontal surfaces
	415	IF ( surf_def_h(0)%ns >= 1 ) THEN
	416	surf => surf_def_h(0)
	417	CALL calc_scaling_parameters
	418	CALL calc_uvw_abs
	419	IF ( .NOT. neutral ) CALL calc_ol
	420	CALL calc_us
	421	CALL calc_surface_fluxes
[3597]	422
	423	IF ( do_output_at_2m ) THEN
	424	CALL calc_pt_near_surface ( '2m' )
	425	ENDIF
[1691]	426	ENDIF
	427	!
[2232]	428	!-- Natural-type horizontal surfaces
	429	IF ( surf_lsm_h%ns >= 1 ) THEN
	430	surf => surf_lsm_h
	431	CALL calc_scaling_parameters
	432	CALL calc_uvw_abs
	433	IF ( .NOT. neutral ) CALL calc_ol
	434	CALL calc_us
	435	CALL calc_surface_fluxes
[3597]	436	IF ( do_output_at_2m ) THEN
	437	CALL calc_pt_near_surface ( '2m' )
	438	ENDIF
[2232]	439	ENDIF
	440	!
	441	!-- Urban-type horizontal surfaces
	442	IF ( surf_usm_h%ns >= 1 ) THEN
	443	surf => surf_usm_h
	444	CALL calc_scaling_parameters
	445	CALL calc_uvw_abs
	446	IF ( .NOT. neutral ) CALL calc_ol
	447	CALL calc_us
	448	CALL calc_surface_fluxes
[3597]	449	IF ( do_output_at_2m ) THEN
	450	CALL calc_pt_near_surface ( '2m' )
	451	ENDIF
	452	IF ( indoor_model ) THEN
	453	CALL calc_pt_near_surface ( '10cm' )
	454	ENDIF
[2232]	455	ENDIF
	456	!
[1691]	457	!-- Use either Newton iteration or a lookup table for the bulk Richardson
	458	!-- number to calculate the Obukhov length
[1788]	459	ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' ) THEN
[2232]	460	!
	461	!-- Default-type upward-facing horizontal surfaces
	462	IF ( surf_def_h(0)%ns >= 1 ) THEN
	463	surf => surf_def_h(0)
	464	CALL calc_uvw_abs
	465	IF ( .NOT. neutral ) CALL calc_ol
	466	CALL calc_us
	467	CALL calc_scaling_parameters
	468	CALL calc_surface_fluxes
[3597]	469	IF ( do_output_at_2m ) THEN
	470	CALL calc_pt_near_surface ( '2m' )
	471	ENDIF
[1691]	472	ENDIF
[2232]	473	!
	474	!-- Natural-type horizontal surfaces
	475	IF ( surf_lsm_h%ns >= 1 ) THEN
	476	surf => surf_lsm_h
	477	CALL calc_uvw_abs
	478	IF ( .NOT. neutral ) CALL calc_ol
	479	CALL calc_us
	480	CALL calc_scaling_parameters
	481	CALL calc_surface_fluxes
[3597]	482	IF ( do_output_at_2m ) THEN
	483	CALL calc_pt_near_surface ( '2m' )
	484	ENDIF
[2232]	485	ENDIF
	486	!
	487	!-- Urban-type horizontal surfaces
	488	IF ( surf_usm_h%ns >= 1 ) THEN
	489	surf => surf_usm_h
	490	CALL calc_uvw_abs
	491	IF ( .NOT. neutral ) CALL calc_ol
	492	CALL calc_us
	493	CALL calc_scaling_parameters
	494	CALL calc_surface_fluxes
[3597]	495	IF ( do_output_at_2m ) THEN
	496	CALL calc_pt_near_surface ( '2m' )
	497	ENDIF
	498	IF ( indoor_model ) THEN
	499	CALL calc_pt_near_surface ( '10cm' )
	500	ENDIF
[2232]	501	ENDIF
[1691]	502
[2232]	503	ENDIF
	504	!
	505	!-- Treat downward-facing horizontal surfaces. Note, so far, these are
	506	!-- always default type. Stratification is not considered
	507	!-- in this case, hence, no further distinction between different
	508	!-- most_method is required.
	509	IF ( surf_def_h(1)%ns >= 1 ) THEN
	510	downward = .TRUE.
	511	surf => surf_def_h(1)
	512	CALL calc_uvw_abs
[1691]	513	CALL calc_us
	514	CALL calc_surface_fluxes
[2232]	515	downward = .FALSE.
[1691]	516	ENDIF
[2232]	517	!
	518	!-- Calculate surfaces fluxes at vertical surfaces for momentum
	519	!-- and subgrid-scale TKE.
	520	!-- No stability is considered. Therefore, scaling parameters and Obukhov-
	521	!-- length do not need to be calculated and no distinction in 'circular',
	522	!-- 'Newton' or 'lookup' is necessary so far.
	523	!-- Note, this will change if stability is once considered.
	524	surf_vertical = .TRUE.
	525	!
	526	!-- Calculate horizontal momentum fluxes at north- and south-facing
	527	!-- surfaces(usvs).
	528	!-- For default-type surfaces
	529	mom_uv = .TRUE.
	530	DO l = 0, 1
	531	IF ( surf_def_v(l)%ns >= 1 ) THEN
	532	surf => surf_def_v(l)
	533	!
	534	!-- Compute surface-parallel velocity
	535	CALL calc_uvw_abs_v_ugrid
	536	!
	537	!-- Compute respective friction velocity on staggered grid
	538	CALL calc_us
	539	!
	540	!-- Compute respective surface fluxes for momentum and TKE
	541	CALL calc_surface_fluxes
	542	ENDIF
	543	ENDDO
	544	!
	545	!-- For natural-type surfaces. Please note, even though stability is not
	546	!-- considered for the calculation of momentum fluxes at vertical surfaces,
	547	!-- scaling parameters and Obukov length are calculated nevertheless in this
	548	!-- case. This is due to the requirement of ts in parameterization of heat
	549	!-- flux in land-surface model in case of aero_resist_kray is not true.
	550	IF ( .NOT. aero_resist_kray ) THEN
	551	IF ( most_method == 'circular' ) THEN
	552	DO l = 0, 1
	553	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	554	surf => surf_lsm_v(l)
	555	!
	556	!-- Compute scaling parameters
	557	CALL calc_scaling_parameters
	558	!
	559	!-- Compute surface-parallel velocity
	560	CALL calc_uvw_abs_v_ugrid
	561	!
	562	!-- Compute Obukhov length
	563	IF ( .NOT. neutral ) CALL calc_ol
	564	!
	565	!-- Compute respective friction velocity on staggered grid
	566	CALL calc_us
	567	!
	568	!-- Compute respective surface fluxes for momentum and TKE
	569	CALL calc_surface_fluxes
	570	ENDIF
	571	ENDDO
	572	ELSE
	573	DO l = 0, 1
	574	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	575	surf => surf_lsm_v(l)
	576	!
	577	!-- Compute surface-parallel velocity
	578	CALL calc_uvw_abs_v_ugrid
	579	!
	580	!-- Compute Obukhov length
	581	IF ( .NOT. neutral ) CALL calc_ol
	582	!
	583	!-- Compute respective friction velocity on staggered grid
	584	CALL calc_us
	585	!
	586	!-- Compute scaling parameters
	587	CALL calc_scaling_parameters
	588	!
	589	!-- Compute respective surface fluxes for momentum and TKE
	590	CALL calc_surface_fluxes
	591	ENDIF
	592	ENDDO
	593	ENDIF
	594	!
	595	!-- No ts is required, so scaling parameters and Obukhov length do not need
	596	!-- to be computed.
	597	ELSE
	598	DO l = 0, 1
	599	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	600	surf => surf_lsm_v(l)
	601	!
	602	!-- Compute surface-parallel velocity
	603	CALL calc_uvw_abs_v_ugrid
	604	!
	605	!-- Compute respective friction velocity on staggered grid
	606	CALL calc_us
	607	!
	608	!-- Compute respective surface fluxes for momentum and TKE
	609	CALL calc_surface_fluxes
	610	ENDIF
	611	ENDDO
	612	ENDIF
	613	!
	614	!-- For urban-type surfaces
	615	DO l = 0, 1
	616	IF ( surf_usm_v(l)%ns >= 1 ) THEN
	617	surf => surf_usm_v(l)
	618	!
	619	!-- Compute surface-parallel velocity
	620	CALL calc_uvw_abs_v_ugrid
	621	!
	622	!-- Compute respective friction velocity on staggered grid
	623	CALL calc_us
	624	!
	625	!-- Compute respective surface fluxes for momentum and TKE
	626	CALL calc_surface_fluxes
	627	ENDIF
	628	ENDDO
	629	!
	630	!-- Calculate horizontal momentum fluxes at east- and west-facing
	631	!-- surfaces (vsus).
	632	!-- For default-type surfaces
	633	DO l = 2, 3
	634	IF ( surf_def_v(l)%ns >= 1 ) THEN
	635	surf => surf_def_v(l)
	636	!
	637	!-- Compute surface-parallel velocity
	638	CALL calc_uvw_abs_v_vgrid
	639	!
	640	!-- Compute respective friction velocity on staggered grid
	641	CALL calc_us
	642	!
	643	!-- Compute respective surface fluxes for momentum and TKE
	644	CALL calc_surface_fluxes
	645	ENDIF
	646	ENDDO
	647	!
	648	!-- For natural-type surfaces. Please note, even though stability is not
	649	!-- considered for the calculation of momentum fluxes at vertical surfaces,
	650	!-- scaling parameters and Obukov length are calculated nevertheless in this
	651	!-- case. This is due to the requirement of ts in parameterization of heat
	652	!-- flux in land-surface model in case of aero_resist_kray is not true.
	653	IF ( .NOT. aero_resist_kray ) THEN
	654	IF ( most_method == 'circular' ) THEN
	655	DO l = 2, 3
	656	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	657	surf => surf_lsm_v(l)
	658	!
	659	!-- Compute scaling parameters
	660	CALL calc_scaling_parameters
	661	!
	662	!-- Compute surface-parallel velocity
	663	CALL calc_uvw_abs_v_vgrid
	664	!
	665	!-- Compute Obukhov length
	666	IF ( .NOT. neutral ) CALL calc_ol
	667	!
	668	!-- Compute respective friction velocity on staggered grid
	669	CALL calc_us
	670	!
	671	!-- Compute respective surface fluxes for momentum and TKE
	672	CALL calc_surface_fluxes
	673	ENDIF
	674	ENDDO
	675	ELSE
	676	DO l = 2, 3
	677	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	678	surf => surf_lsm_v(l)
	679	!
	680	!-- Compute surface-parallel velocity
	681	CALL calc_uvw_abs_v_vgrid
	682	!
	683	!-- Compute Obukhov length
	684	IF ( .NOT. neutral ) CALL calc_ol
	685	!
	686	!-- Compute respective friction velocity on staggered grid
	687	CALL calc_us
	688	!
	689	!-- Compute scaling parameters
	690	CALL calc_scaling_parameters
	691	!
	692	!-- Compute respective surface fluxes for momentum and TKE
	693	CALL calc_surface_fluxes
	694	ENDIF
	695	ENDDO
	696	ENDIF
	697	ELSE
	698	DO l = 2, 3
	699	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	700	surf => surf_lsm_v(l)
	701	!
	702	!-- Compute surface-parallel velocity
	703	CALL calc_uvw_abs_v_vgrid
	704	!
	705	!-- Compute respective friction velocity on staggered grid
	706	CALL calc_us
	707	!
	708	!-- Compute respective surface fluxes for momentum and TKE
	709	CALL calc_surface_fluxes
	710	ENDIF
	711	ENDDO
	712	ENDIF
	713	!
	714	!-- For urban-type surfaces
	715	DO l = 2, 3
	716	IF ( surf_usm_v(l)%ns >= 1 ) THEN
	717	surf => surf_usm_v(l)
	718	!
	719	!-- Compute surface-parallel velocity
	720	CALL calc_uvw_abs_v_vgrid
	721	!
	722	!-- Compute respective friction velocity on staggered grid
	723	CALL calc_us
	724	!
	725	!-- Compute respective surface fluxes for momentum and TKE
	726	CALL calc_surface_fluxes
	727	ENDIF
	728	ENDDO
	729	mom_uv = .FALSE.
	730	!
	731	!-- Calculate horizontal momentum fluxes of w (wsus and wsvs) at vertial
	732	!-- surfaces.
	733	mom_w = .TRUE.
	734	!
	735	!-- Default-type surfaces
	736	DO l = 0, 3
	737	IF ( surf_def_v(l)%ns >= 1 ) THEN
	738	surf => surf_def_v(l)
	739	CALL calc_uvw_abs_v_wgrid
	740	CALL calc_us
	741	CALL calc_surface_fluxes
	742	ENDIF
	743	ENDDO
	744	!
	745	!-- Natural-type surfaces
	746	DO l = 0, 3
	747	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	748	surf => surf_lsm_v(l)
	749	CALL calc_uvw_abs_v_wgrid
	750	CALL calc_us
	751	CALL calc_surface_fluxes
	752	ENDIF
	753	ENDDO
	754	!
	755	!-- Urban-type surfaces
	756	DO l = 0, 3
	757	IF ( surf_usm_v(l)%ns >= 1 ) THEN
	758	surf => surf_usm_v(l)
	759	CALL calc_uvw_abs_v_wgrid
	760	CALL calc_us
	761	CALL calc_surface_fluxes
	762	ENDIF
	763	ENDDO
	764	mom_w = .FALSE.
	765	!
	766	!-- Calculate momentum fluxes usvs, vsus, wsus and wsvs at vertical
	767	!-- surfaces for TKE production. Note, here, momentum fluxes are defined
	768	!-- at grid center and are not staggered as before.
	769	mom_tke = .TRUE.
	770	!
	771	!-- Default-type surfaces
	772	DO l = 0, 3
	773	IF ( surf_def_v(l)%ns >= 1 ) THEN
	774	surf => surf_def_v(l)
	775	CALL calc_uvw_abs_v_sgrid
	776	CALL calc_us
	777	CALL calc_surface_fluxes
	778	ENDIF
	779	ENDDO
	780	!
	781	!-- Natural-type surfaces
	782	DO l = 0, 3
	783	IF ( surf_lsm_v(l)%ns >= 1 ) THEN
	784	surf => surf_lsm_v(l)
	785	CALL calc_uvw_abs_v_sgrid
	786	CALL calc_us
	787	CALL calc_surface_fluxes
	788	ENDIF
	789	ENDDO
	790	!
	791	!-- Urban-type surfaces
	792	DO l = 0, 3
	793	IF ( surf_usm_v(l)%ns >= 1 ) THEN
	794	surf => surf_usm_v(l)
	795	CALL calc_uvw_abs_v_sgrid
	796	CALL calc_us
	797	CALL calc_surface_fluxes
	798	ENDIF
	799	ENDDO
	800	mom_tke = .FALSE.
	801
[1691]	802
	803	END SUBROUTINE surface_layer_fluxes
	804
	805
	806	!------------------------------------------------------------------------------!
	807	! Description:
	808	! ------------
	809	!> Initializing actions for the surface layer routine. Basically, this involves
	810	!> the preparation of a lookup table for the the bulk Richardson number vs
	811	!> Obukhov length L when using the lookup table method.
	812	!------------------------------------------------------------------------------!
	813	SUBROUTINE init_surface_layer_fluxes
	814
	815	IMPLICIT NONE
	816
[2232]	817	INTEGER(iwp) :: li, & !< Index for loop to create lookup table
[1691]	818	num_steps_n !< Number of non-stretched zeta steps
	819
	820	LOGICAL :: terminate_run_l = .FALSE. !< Flag to terminate run (global)
	821
	822	REAL(wp), PARAMETER :: zeta_stretch = -10.0_wp !< Start of stretching in the free convection limit
	823
	824	REAL(wp), DIMENSION(:), ALLOCATABLE :: zeta_tmp
	825
	826
	827	REAL(wp) :: zeta_step, & !< Increment of zeta
	828	regr = 1.01_wp, & !< Stretching factor of zeta_step in the free convection limit
	829	regr_old = 1.0E9_wp, & !< Stretching factor of last iteration step
	830	z0h_min = 0.0_wp, & !< Minimum value of z0h to create table
	831	z0_min = 0.0_wp !< Minimum value of z0 to create table
	832
	833
[1709]	834
[2232]	835
[1709]	836	!
	837	!-- In case of runs with neutral statification, set Obukhov length to a
	838	!-- large value
[2232]	839	IF ( neutral ) THEN
	840	IF ( surf_def_h(0)%ns >= 1 ) surf_def_h(0)%ol = 1.0E10_wp
	841	IF ( surf_lsm_h%ns >= 1 ) surf_lsm_h%ol = 1.0E10_wp
	842	IF ( surf_usm_h%ns >= 1 ) surf_usm_h%ol = 1.0E10_wp
	843	ENDIF
[1709]	844
[1691]	845	IF ( most_method == 'lookup' ) THEN
	846
	847	!
	848	!-- Check for roughness heterogeneity. In that case terminate run and
[2232]	849	!-- inform user. Check for both, natural and non-natural walls.
	850	IF ( surf_def_h(0)%ns >= 1 ) THEN
	851	IF ( MINVAL( surf_def_h(0)%z0h ) /= MAXVAL( surf_def_h(0)%z0h ) .OR. &
	852	MINVAL( surf_def_h(0)%z0 ) /= MAXVAL( surf_def_h(0)%z0 ) ) THEN
	853	terminate_run_l = .TRUE.
	854	ENDIF
[1691]	855	ENDIF
[2232]	856	IF ( surf_lsm_h%ns >= 1 ) THEN
	857	IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h ) .OR. &
	858	MINVAL( surf_lsm_h%z0 ) /= MAXVAL( surf_lsm_h%z0 ) ) THEN
	859	terminate_run_l = .TRUE.
	860	ENDIF
	861	ENDIF
	862	IF ( surf_usm_h%ns >= 1 ) THEN
	863	IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_usm_h%z0h ) .OR. &
	864	MINVAL( surf_usm_h%z0 ) /= MAXVAL( surf_usm_h%z0 ) ) THEN
	865	terminate_run_l = .TRUE.
	866	ENDIF
	867	ENDIF
	868	!
	869	!-- Check roughness homogeneity between differt surface types.
	870	IF ( surf_lsm_h%ns >= 1 .AND. surf_def_h(0)%ns >= 1 ) THEN
	871	IF ( MINVAL( surf_lsm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h ) .OR. &
	872	MINVAL( surf_lsm_h%z0 ) /= MAXVAL( surf_def_h(0)%z0 ) ) THEN
	873	terminate_run_l = .TRUE.
	874	ENDIF
	875	ENDIF
	876	IF ( surf_usm_h%ns >= 1 .AND. surf_def_h(0)%ns >= 1 ) THEN
	877	IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_def_h(0)%z0h ) .OR. &
	878	MINVAL( surf_usm_h%z0 ) /= MAXVAL( surf_def_h(0)%z0 ) ) THEN
	879	terminate_run_l = .TRUE.
	880	ENDIF
	881	ENDIF
	882	IF ( surf_usm_h%ns >= 1 .AND. surf_lsm_h%ns >= 1 ) THEN
	883	IF ( MINVAL( surf_usm_h%z0h ) /= MAXVAL( surf_lsm_h%z0h ) .OR. &
	884	MINVAL( surf_usm_h%z0 ) /= MAXVAL( surf_lsm_h%z0 ) ) THEN
	885	terminate_run_l = .TRUE.
	886	ENDIF
	887	ENDIF
[1691]	888
[2232]	889
[1691]	890	#if defined( __parallel )
	891	!
	892	!-- Make a logical OR for all processes. Force termiation of model if result
	893	!-- is TRUE
	894	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	895	CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL, &
	896	MPI_LOR, comm2d, ierr )
	897	#else
	898	terminate_run = terminate_run_l
	899	#endif
	900
	901	IF ( terminate_run ) THEN
	902	message_string = 'most_method = "lookup" cannot be used in ' // &
	903	'combination with a prescribed roughness ' // &
	904	'heterogeneity'
[3045]	905	CALL message( 'surface_layer_fluxes', 'PA0116', 1, 2, 0, 6, 0 )
[1691]	906	ENDIF
	907
	908	ALLOCATE( zeta_tmp(0:num_steps-1) )
	909
	910	!
	911	!-- Use the lowest possible value for z_mo
[2232]	912	k = nzb
[1691]	913	z_mo = zu(k+1) - zw(k)
	914
	915	!
	916	!-- Calculate z/L range from zeta_stretch to zeta_max using 90% of the
	917	!-- available steps (num_steps). The calculation is done with negative
	918	!-- values of zeta in order to simplify the stretching in the free
	919	!-- convection limit for the remaining 10% of steps.
	920	zeta_tmp(0) = - zeta_max
	921	num_steps_n = ( num_steps * 9 / 10 ) - 1
	922	zeta_step = (zeta_max - zeta_stretch) / REAL(num_steps_n)
	923
[2232]	924	DO li = 1, num_steps_n
	925	zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
[1691]	926	ENDDO
	927
	928	!
	929	!-- Calculate stretching factor for the free convection range
[1788]	930	DO WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
[1691]	931	regr_old = regr
	932	regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr ) &
	933	)**( 10.0_wp / REAL(num_steps) )
	934	ENDDO
	935
	936	!
	937	!-- Calculate z/L range from zeta_min to zeta_stretch
[2232]	938	DO li = num_steps_n+1, num_steps-1
	939	zeta_tmp(li) = zeta_tmp(li-1) + zeta_step
[1691]	940	zeta_step = zeta_step * regr
	941	ENDDO
	942
	943	!
[1757]	944	!-- Save roughness lengths to temporary variables
[2232]	945	IF ( surf_def_h(0)%ns >= 1 ) THEN
	946	z0h_min = surf_def_h(0)%z0h(1)
	947	z0_min = surf_def_h(0)%z0(1)
	948	ELSEIF ( surf_lsm_h%ns >= 1 ) THEN
	949	z0h_min = surf_lsm_h%z0h(1)
	950	z0_min = surf_lsm_h%z0(1)
	951	ELSEIF ( surf_usm_h%ns >= 1 ) THEN
	952	z0h_min = surf_usm_h%z0h(1)
	953	z0_min = surf_usm_h%z0(1)
	954	ENDIF
[1691]	955	!
	956	!-- Calculate lookup table for the Richardson number versus Obukhov length
	957	!-- The Richardson number (rib) is defined depending on the choice of
	958	!-- boundary conditions for temperature
	959	IF ( ibc_pt_b == 1 ) THEN
[2232]	960	DO li = 0, num_steps-1
	961	ol_tab(li) = - z_mo / zeta_tmp(num_steps-1-li)
	962	rib_tab(li) = z_mo / ol_tab(li) / ( LOG( z_mo / z0_min ) &
	963	- psi_m( z_mo / ol_tab(li) ) &
	964	+ psi_m( z0_min / ol_tab(li) ) &
[1691]	965	)**3
	966	ENDDO
	967	ELSE
[2232]	968	DO li = 0, num_steps-1
	969	ol_tab(li) = - z_mo / zeta_tmp(num_steps-1-li)
	970	rib_tab(li) = z_mo / ol_tab(li) * ( LOG( z_mo / z0h_min ) &
	971	- psi_h( z_mo / ol_tab(li) ) &
	972	+ psi_h( z0h_min / ol_tab(li) ) &
[1691]	973	) &
	974	/ ( LOG( z_mo / z0_min ) &
[2232]	975	- psi_m( z_mo / ol_tab(li) ) &
	976	+ psi_m( z0_min / ol_tab(li) ) &
[1691]	977	)**2
	978	ENDDO
	979	ENDIF
	980
	981	!
	982	!-- Determine minimum values of rib in the lookup table. Set upper limit
	983	!-- to critical Richardson number (0.25)
	984	rib_min = MINVAL(rib_tab)
	985	rib_max = 0.25 !MAXVAL(rib_tab)
	986
	987	DEALLOCATE( zeta_tmp )
	988	ENDIF
	989
	990	END SUBROUTINE init_surface_layer_fluxes
	991
	992
	993	!------------------------------------------------------------------------------!
	994	! Description:
	995	! ------------
[1709]	996	!> Compute the absolute value of the horizontal velocity (relative to the
[2232]	997	!> surface) for horizontal surface elements. This is required by all methods.
[1691]	998	!------------------------------------------------------------------------------!
[2232]	999	SUBROUTINE calc_uvw_abs
[3157]	1000
[1691]	1001	IMPLICIT NONE
	1002
[2232]	1003	INTEGER(iwp) :: i !< running index x direction
	1004	INTEGER(iwp) :: ibit !< flag to mask computation of relative velocity in case of downward-facing surfaces
	1005	INTEGER(iwp) :: j !< running index y direction
	1006	INTEGER(iwp) :: k !< running index z direction
	1007	INTEGER(iwp) :: m !< running index surface elements
[1691]	1008
[3157]	1009	REAL(wp) :: w_lfc !< local free convection velocity scale
[2232]	1010	!
	1011	!-- ibit is 1 for upward-facing surfaces, zero for downward-facing surfaces.
	1012	ibit = MERGE( 1, 0, .NOT. downward )
[1691]	1013
[3634]	1014	!$ACC PARALLEL LOOP PRIVATE(i, j, k, w_lfc) &
	1015	!$ACC PRESENT(surf, u, v)
[2232]	1016	DO m = 1, surf%ns
[1691]	1017
[2232]	1018	i = surf%i(m)
	1019	j = surf%j(m)
	1020	k = surf%k(m)
[3157]	1021
[1691]	1022	!
[3157]	1023	!-- Calculate free convection velocity scale w_lfc is
	1024	!-- use_free_convection_scaling = .T.. This will maintain a horizontal
	1025	!-- velocity even for very weak wind convective conditions. SIGN is used
	1026	!-- to set w_lfc to zero under stable conditions.
	1027	IF ( use_free_convection_scaling ) THEN
	1028	w_lfc = ABS(g / surf%pt1(m) * surf%z_mo(m) * surf%shf(m))
	1029	w_lfc = ( 0.5_wp * ( w_lfc + SIGN(w_lfc,surf%shf(m)) ) )**(0.33333_wp)
	1030	ELSE
	1031	w_lfc = 0.0_wp
	1032	ENDIF
	1033
	1034	!
[2232]	1035	!-- Compute the absolute value of the horizontal velocity.
	1036	!-- (relative to the surface in case the lower surface is the ocean).
	1037	!-- Please note, in new surface modelling concept the index values changed,
	1038	!-- i.e. the reference grid point is not the surface-grid point itself but
	1039	!-- the first grid point outside of the topography.
	1040	!-- Note, in case of coupled ocean-atmosphere simulations relative velocity
	1041	!-- with respect to the ocean surface is used, hence, (k-1,j,i) values
	1042	!-- are used to calculate the absolute velocity.
	1043	!-- However, this do not apply for downward-facing walls. To mask this,
	1044	!-- use ibit, which checks for upward/downward-facing surfaces.
	1045	surf%uvw_abs(m) = SQRT( &
	1046	( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) &
	1047	- ( u(k-1,j,i) + u(k-1,j,i+1) &
	1048	) * ibit &
	1049	) &
	1050	)**2 + &
	1051	( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) &
	1052	- ( v(k-1,j,i) + v(k-1,j+1,i) &
	1053	) * ibit &
	1054	) &
[3157]	1055	)2 + w_lfc2 &
[2232]	1056	)
	1057
[3148]	1058
	1059
[1691]	1060	ENDDO
	1061
[2232]	1062	END SUBROUTINE calc_uvw_abs
	1063
	1064
	1065	!------------------------------------------------------------------------------!
	1066	! Description:
	1067	! ------------
	1068	!> Compute the absolute value of the horizontal velocity (relative to the
	1069	!> surface) for horizontal surface elements. This is required by all methods.
	1070	!------------------------------------------------------------------------------!
	1071	SUBROUTINE calc_uvw_abs_v_ugrid
	1072
	1073	IMPLICIT NONE
	1074
[3547]	1075	INTEGER(iwp) :: i !< running index x direction
	1076	INTEGER(iwp) :: j !< running index y direction
	1077	INTEGER(iwp) :: k !< running index z direction
	1078	INTEGER(iwp) :: m !< running index surface elements
[2232]	1079
[3547]	1080	REAL(wp) :: u_i !< u-component on xu-grid
	1081	REAL(wp) :: w_i !< w-component on xu-grid
[2232]	1082
	1083
	1084	DO m = 1, surf%ns
	1085	i = surf%i(m)
	1086	j = surf%j(m)
	1087	k = surf%k(m)
[1691]	1088	!
[2232]	1089	!-- Compute the absolute value of the surface parallel velocity on u-grid.
	1090	u_i = u(k,j,i)
	1091	w_i = 0.25_wp * ( w(k-1,j,i-1) + w(k-1,j,i) + &
	1092	w(k,j,i-1) + w(k,j,i) )
[1691]	1093
[2232]	1094	surf%uvw_abs(m) = SQRT( u_i2 + w_i2 )
[1709]	1095
[2232]	1096	ENDDO
[1709]	1097
[2232]	1098	END SUBROUTINE calc_uvw_abs_v_ugrid
	1099
[1709]	1100	!------------------------------------------------------------------------------!
	1101	! Description:
	1102	! ------------
[2232]	1103	!> Compute the absolute value of the horizontal velocity (relative to the
	1104	!> surface) for horizontal surface elements. This is required by all methods.
	1105	!------------------------------------------------------------------------------!
	1106	SUBROUTINE calc_uvw_abs_v_vgrid
	1107
	1108	IMPLICIT NONE
	1109
[3547]	1110	INTEGER(iwp) :: i !< running index x direction
	1111	INTEGER(iwp) :: j !< running index y direction
	1112	INTEGER(iwp) :: k !< running index z direction
	1113	INTEGER(iwp) :: m !< running index surface elements
[2232]	1114
[3547]	1115	REAL(wp) :: v_i !< v-component on yv-grid
	1116	REAL(wp) :: w_i !< w-component on yv-grid
[2232]	1117
	1118
	1119	DO m = 1, surf%ns
	1120	i = surf%i(m)
	1121	j = surf%j(m)
	1122	k = surf%k(m)
	1123
	1124	v_i = u(k,j,i)
	1125	w_i = 0.25_wp * ( w(k-1,j-1,i) + w(k-1,j,i) + &
	1126	w(k,j-1,i) + w(k,j,i) )
	1127
	1128	surf%uvw_abs(m) = SQRT( v_i2 + w_i2 )
	1129
	1130	ENDDO
	1131
	1132	END SUBROUTINE calc_uvw_abs_v_vgrid
	1133
	1134	!------------------------------------------------------------------------------!
	1135	! Description:
	1136	! ------------
	1137	!> Compute the absolute value of the horizontal velocity (relative to the
	1138	!> surface) for horizontal surface elements. This is required by all methods.
	1139	!------------------------------------------------------------------------------!
	1140	SUBROUTINE calc_uvw_abs_v_wgrid
	1141
	1142	IMPLICIT NONE
	1143
[3547]	1144	INTEGER(iwp) :: i !< running index x direction
	1145	INTEGER(iwp) :: j !< running index y direction
	1146	INTEGER(iwp) :: k !< running index z direction
	1147	INTEGER(iwp) :: m !< running index surface elements
[2232]	1148
[3547]	1149	REAL(wp) :: u_i !< u-component on x-zw-grid
	1150	REAL(wp) :: v_i !< v-component on y-zw-grid
	1151	REAL(wp) :: w_i !< w-component on zw-grid
[2232]	1152	!
	1153	!-- North- (l=0) and south-facing (l=1) surfaces
	1154	IF ( l == 0 .OR. l == 1 ) THEN
	1155	DO m = 1, surf%ns
	1156	i = surf%i(m)
	1157	j = surf%j(m)
	1158	k = surf%k(m)
	1159
	1160	u_i = 0.25_wp * ( u(k+1,j,i+1) + u(k+1,j,i) + &
	1161	u(k,j,i+1) + u(k,j,i) )
	1162	v_i = 0.0_wp
	1163	w_i = w(k,j,i)
	1164
	1165	surf%uvw_abs(m) = SQRT( u_i2 + v_i2 + w_i**2 )
	1166	ENDDO
	1167	!
	1168	!-- East- (l=2) and west-facing (l=3) surfaces
	1169	ELSE
	1170	DO m = 1, surf%ns
	1171	i = surf%i(m)
	1172	j = surf%j(m)
	1173	k = surf%k(m)
	1174
	1175	u_i = 0.0_wp
	1176	v_i = 0.25_wp * ( v(k+1,j+1,i) + v(k+1,j,i) + &
	1177	v(k,j+1,i) + v(k,j,i) )
	1178	w_i = w(k,j,i)
	1179
	1180	surf%uvw_abs(m) = SQRT( u_i2 + v_i2 + w_i**2 )
	1181	ENDDO
	1182	ENDIF
	1183
	1184	END SUBROUTINE calc_uvw_abs_v_wgrid
	1185
	1186	!------------------------------------------------------------------------------!
	1187	! Description:
	1188	! ------------
	1189	!> Compute the absolute value of the horizontal velocity (relative to the
	1190	!> surface) for horizontal surface elements. This is required by all methods.
	1191	!------------------------------------------------------------------------------!
	1192	SUBROUTINE calc_uvw_abs_v_sgrid
	1193
	1194	IMPLICIT NONE
	1195
[3547]	1196	INTEGER(iwp) :: i !< running index x direction
	1197	INTEGER(iwp) :: j !< running index y direction
	1198	INTEGER(iwp) :: k !< running index z direction
	1199	INTEGER(iwp) :: m !< running index surface elements
[2232]	1200
[3547]	1201	REAL(wp) :: u_i !< u-component on scalar grid
	1202	REAL(wp) :: v_i !< v-component on scalar grid
	1203	REAL(wp) :: w_i !< w-component on scalar grid
[2232]	1204
	1205	!
	1206	!-- North- (l=0) and south-facing (l=1) walls
	1207	IF ( l == 0 .OR. l == 1 ) THEN
	1208	DO m = 1, surf%ns
	1209	i = surf%i(m)
	1210	j = surf%j(m)
	1211	k = surf%k(m)
	1212
	1213	u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
	1214	v_i = 0.0_wp
	1215	w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
	1216
	1217	surf%uvw_abs(m) = SQRT( u_i2 + v_i2 + w_i**2 )
	1218	ENDDO
	1219	!
	1220	!-- East- (l=2) and west-facing (l=3) walls
	1221	ELSE
	1222	DO m = 1, surf%ns
	1223	i = surf%i(m)
	1224	j = surf%j(m)
	1225	k = surf%k(m)
	1226
	1227	u_i = 0.0_wp
	1228	v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
	1229	w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
	1230
	1231	surf%uvw_abs(m) = SQRT( u_i2 + v_i2 + w_i**2 )
	1232	ENDDO
	1233	ENDIF
	1234
	1235	END SUBROUTINE calc_uvw_abs_v_sgrid
	1236
	1237
	1238	!------------------------------------------------------------------------------!
	1239	! Description:
	1240	! ------------
[1709]	1241	!> Calculate the Obukhov length (L) and Richardson flux number (z/L)
	1242	!------------------------------------------------------------------------------!
	1243	SUBROUTINE calc_ol
	1244
	1245	IMPLICIT NONE
	1246
[2232]	1247	INTEGER(iwp) :: iter !< Newton iteration step
	1248	INTEGER(iwp) :: li !< look index
	1249	INTEGER(iwp) :: m !< loop variable over all horizontal wall elements
[1709]	1250
	1251	REAL(wp) :: f, & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
	1252	f_d_ol, & !< Derivative of f
	1253	ol_l, & !< Lower bound of L for Newton iteration
	1254	ol_m, & !< Previous value of L for Newton iteration
	1255	ol_old, & !< Previous time step value of L
	1256	ol_u !< Upper bound of L for Newton iteration
	1257
[1691]	1258	IF ( TRIM( most_method ) /= 'circular' ) THEN
[2232]	1259	!
	1260	!-- Evaluate bulk Richardson number (calculation depends on
	1261	!-- definition based on setting of boundary conditions
	1262	IF ( ibc_pt_b /= 1 ) THEN
	1263	IF ( humidity ) THEN
[3146]	1264	!$OMP PARALLEL DO PRIVATE( z_mo )
[2232]	1265	DO m = 1, surf%ns
[1691]	1266
[2232]	1267	z_mo = surf%z_mo(m)
	1268
[3146]	1269	surf%rib(m) = g * z_mo &
	1270	* ( surf%vpt1(m) - surf%vpt_surface(m) ) &
	1271	/ ( surf%uvw_abs(m)*2 surf%vpt1(m) &
	1272	+ 1.0E-20_wp )
[2232]	1273	ENDDO
	1274	ELSE
[3146]	1275	!$OMP PARALLEL DO PRIVATE( z_mo )
[2232]	1276	DO m = 1, surf%ns
	1277
	1278	z_mo = surf%z_mo(m)
	1279
[3146]	1280	surf%rib(m) = g * z_mo &
	1281	* ( surf%pt1(m) - surf%pt_surface(m) ) &
	1282	/ ( surf%uvw_abs(m)*2 surf%pt1(m) &
	1283	+ 1.0E-20_wp )
[2232]	1284	ENDDO
	1285	ENDIF
	1286	ELSE
	1287	IF ( humidity ) THEN
[3146]	1288	!$OMP PARALLEL DO PRIVATE( k, z_mo )
[2232]	1289	DO m = 1, surf%ns
	1290
	1291	k = surf%k(m)
	1292
	1293	z_mo = surf%z_mo(m)
	1294
	1295	surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp &
[3146]	1296	* surf%qv1(m) ) * surf%shf(m) + 0.61_wp &
	1297	* surf%pt1(m) * surf%qsws(m) ) * &
[2232]	1298	drho_air_zw(k-1) / &
[3146]	1299	( surf%uvw_abs(m)*3 surf%vpt1(m) * kappa**2 &
[1709]	1300	+ 1.0E-20_wp )
[2232]	1301	ENDDO
	1302	ELSE
[3146]	1303	!$OMP PARALLEL DO PRIVATE( k, z_mo )
[3634]	1304	!$ACC PARALLEL LOOP PRIVATE(k, z_mo) &
	1305	!$ACC PRESENT(surf, drho_air_zw)
[2232]	1306	DO m = 1, surf%ns
[1691]	1307
[2232]	1308	k = surf%k(m)
	1309
	1310	z_mo = surf%z_mo(m)
	1311
	1312	surf%rib(m) = - g * z_mo * surf%shf(m) * &
	1313	drho_air_zw(k-1) / &
[3146]	1314	( surf%uvw_abs(m)*3 surf%pt1(m) * kappa**2 &
[2232]	1315	+ 1.0E-20_wp )
	1316	ENDDO
	1317	ENDIF
	1318	ENDIF
	1319
[1691]	1320	ENDIF
	1321
[2232]	1322
[1691]	1323	!
	1324	!-- Calculate the Obukhov length either using a Newton iteration
	1325	!-- method, via a lookup table, or using the old circular way
	1326	IF ( TRIM( most_method ) == 'newton' ) THEN
	1327
[3634]	1328	!$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
	1329	!$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
	1330	!$ACC PRESENT(surf)
[2232]	1331	DO m = 1, surf%ns
[1691]	1332
[2232]	1333	i = surf%i(m)
	1334	j = surf%j(m)
[1691]	1335
[2232]	1336	z_mo = surf%z_mo(m)
	1337
[1691]	1338	!
[2232]	1339	!-- Store current value in case the Newton iteration fails
	1340	ol_old = surf%ol(m)
[1691]	1341
	1342	!
[2232]	1343	!-- Ensure that the bulk Richardson number and the Obukhov
	1344	!-- length have the same sign
	1345	IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. &
	1346	ABS( surf%ol(m) ) == ol_max ) THEN
	1347	IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp
	1348	IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
	1349	ENDIF
[1691]	1350	!
[2232]	1351	!-- Iteration to find Obukhov length
	1352	iter = 0
	1353	DO
	1354	iter = iter + 1
[1691]	1355	!
[2232]	1356	!-- In case of divergence, use the value of the previous time step
	1357	IF ( iter > 1000 ) THEN
	1358	surf%ol(m) = ol_old
	1359	EXIT
	1360	ENDIF
[1691]	1361
[2232]	1362	ol_m = surf%ol(m)
	1363	ol_l = ol_m - 0.001_wp * ol_m
	1364	ol_u = ol_m + 0.001_wp * ol_m
[1691]	1365
	1366
[2232]	1367	IF ( ibc_pt_b /= 1 ) THEN
[1691]	1368	!
[2232]	1369	!-- Calculate f = Ri - [...]/[...]^2 = 0
	1370	f = surf%rib(m) - ( z_mo / ol_m ) * ( &
	1371	LOG( z_mo / surf%z0h(m) ) &
	1372	- psi_h( z_mo / ol_m ) &
	1373	+ psi_h( surf%z0h(m) / &
	1374	ol_m ) &
	1375	) &
	1376	/ ( LOG( z_mo / surf%z0(m) ) &
	1377	- psi_m( z_mo / ol_m ) &
	1378	+ psi_m( surf%z0(m) / &
	1379	ol_m ) &
	1380	)**2
[1691]	1381
	1382	!
[2232]	1383	!-- Calculate df/dL
	1384	f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / &
	1385	surf%z0h(m) ) &
	1386	- psi_h( z_mo / ol_u ) &
	1387	+ psi_h( surf%z0h(m) / ol_u ) &
	1388	) &
	1389	/ ( LOG( z_mo / surf%z0(m) ) &
	1390	- psi_m( z_mo / ol_u ) &
	1391	+ psi_m( surf%z0(m) / ol_u ) &
	1392	)**2 &
	1393	+ ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) ) &
	1394	- psi_h( z_mo / ol_l ) &
	1395	+ psi_h( surf%z0h(m) / ol_l ) &
	1396	) &
	1397	/ ( LOG( z_mo / surf%z0(m) ) &
	1398	- psi_m( z_mo / ol_l ) &
	1399	+ psi_m( surf%z0(m) / ol_l ) &
	1400	)**2 &
	1401	) / ( ol_u - ol_l )
	1402	ELSE
[1691]	1403	!
[2232]	1404	!-- Calculate f = Ri - 1 /[...]^3 = 0
	1405	f = surf%rib(m) - ( z_mo / ol_m ) / &
	1406	( LOG( z_mo / surf%z0(m) ) &
	1407	- psi_m( z_mo / ol_m ) &
	1408	+ psi_m( surf%z0(m) / ol_m ) &
	1409	)**3
[1691]	1410
	1411	!
[2232]	1412	!-- Calculate df/dL
	1413	f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / &
	1414	surf%z0(m) ) &
	1415	- psi_m( z_mo / ol_u ) &
	1416	+ psi_m( surf%z0(m) / ol_u ) &
	1417	)**3 &
	1418	+ ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) ) &
	1419	- psi_m( z_mo / ol_l ) &
	1420	+ psi_m( surf%z0(m) / ol_l ) &
	1421	)**3 &
	1422	) / ( ol_u - ol_l )
	1423	ENDIF
[1691]	1424	!
[2232]	1425	!-- Calculate new L
	1426	surf%ol(m) = ol_m - f / f_d_ol
[1691]	1427
	1428	!
[2232]	1429	!-- Ensure that the bulk Richardson number and the Obukhov
	1430	!-- length have the same sign and ensure convergence.
	1431	IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp
[1691]	1432	!
[2232]	1433	!-- If unrealistic value occurs, set L to the maximum
	1434	!-- value that is allowed
	1435	IF ( ABS( surf%ol(m) ) > ol_max ) THEN
	1436	surf%ol(m) = ol_max
	1437	EXIT
	1438	ENDIF
[1691]	1439	!
[2232]	1440	!-- Check for convergence
	1441	IF ( ABS( ( surf%ol(m) - ol_m ) / &
	1442	surf%ol(m) ) < 1.0E-4_wp ) THEN
	1443	EXIT
	1444	ELSE
	1445	CYCLE
	1446	ENDIF
[1691]	1447
	1448	ENDDO
	1449	ENDDO
	1450
	1451	ELSEIF ( TRIM( most_method ) == 'lookup' ) THEN
	1452
[2281]	1453	!$OMP PARALLEL DO PRIVATE( i, j, z_mo, li ) FIRSTPRIVATE( li_bnd ) LASTPRIVATE( li_bnd )
[2232]	1454	DO m = 1, surf%ns
[1691]	1455
[2232]	1456	i = surf%i(m)
	1457	j = surf%j(m)
[1691]	1458	!
[2232]	1459	!-- If the bulk Richardson number is outside the range of the lookup
	1460	!-- table, set it to the exceeding threshold value
	1461	IF ( surf%rib(m) < rib_min ) surf%rib(m) = rib_min
	1462	IF ( surf%rib(m) > rib_max ) surf%rib(m) = rib_max
[1691]	1463
	1464	!
[2232]	1465	!-- Find the correct index bounds for linear interpolation. As the
	1466	!-- Richardson number will not differ very much from time step to
	1467	!-- time step , use the index from the last step and search in the
	1468	!-- correct direction
	1469	li = li_bnd
[2321]	1470	IF ( rib_tab(li) - surf%rib(m) > 0.0_wp ) THEN
	1471	DO WHILE ( rib_tab(li-1) - surf%rib(m) > 0.0_wp .AND. li > 0 )
[2232]	1472	li = li-1
	1473	ENDDO
	1474	ELSE
	1475	DO WHILE ( rib_tab(li) - surf%rib(m) < 0.0_wp &
	1476	.AND. li < num_steps-1 )
	1477	li = li+1
	1478	ENDDO
	1479	ENDIF
	1480	li_bnd = li
[1691]	1481
	1482	!
[2232]	1483	!-- Linear interpolation to find the correct value of z/L
	1484	surf%ol(m) = ( ol_tab(li-1) + ( ol_tab(li) - ol_tab(li-1) ) &
	1485	/ ( rib_tab(li) - rib_tab(li-1) ) &
	1486	* ( surf%rib(m) - rib_tab(li-1) ) )
[1691]	1487
	1488	ENDDO
	1489
	1490	ELSEIF ( TRIM( most_method ) == 'circular' ) THEN
	1491
[2232]	1492	IF ( .NOT. humidity ) THEN
[3146]	1493	!$OMP PARALLEL DO PRIVATE( z_mo )
[2232]	1494	DO m = 1, surf%ns
[1691]	1495
[2232]	1496	z_mo = surf%z_mo(m)
[1691]	1497
[3146]	1498	surf%ol(m) = ( surf%pt1(m) * surf%us(m)**2 ) / &
[2232]	1499	( kappa * g * &
	1500	surf%ts(m) + 1E-30_wp )
[1691]	1501	!
	1502	!-- Limit the value range of the Obukhov length.
[2232]	1503	!-- This is necessary for very small velocities (u,v --> 1), because
[1691]	1504	!-- the absolute value of ol can then become very small, which in
	1505	!-- consequence would result in very large shear stresses and very
	1506	!-- small momentum fluxes (both are generally unrealistic).
[2232]	1507	IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min ) &
	1508	surf%ol(m) = z_mo / zeta_min
	1509	IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max ) &
	1510	surf%ol(m) = z_mo / zeta_max
	1511
[1691]	1512	ENDDO
[3274]	1513	ELSEIF ( bulk_cloud_model .OR. cloud_droplets ) THEN
[2232]	1514	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	1515	DO m = 1, surf%ns
[1691]	1516
[2232]	1517	i = surf%i(m)
	1518	j = surf%j(m)
	1519	k = surf%k(m)
[1691]	1520
[2232]	1521	z_mo = surf%z_mo(m)
	1522
	1523
[3146]	1524	surf%ol(m) = ( surf%vpt1(m) * surf%us(m)**2 ) / &
[2232]	1525	( kappa * g * ( surf%ts(m) + &
	1526	0.61_wp * surf%pt1(m) * surf%us(m) &
	1527	+ 0.61_wp * surf%qv1(m) * surf%ts(m) - &
	1528	surf%ts(m) * ql(k,j,i) ) + 1E-30_wp )
[1691]	1529	!
[2232]	1530	!-- Limit the value range of the Obukhov length.
	1531	!-- This is necessary for very small velocities (u,v --> 1), because
	1532	!-- the absolute value of ol can then become very small, which in
	1533	!-- consequence would result in very large shear stresses and very
	1534	!-- small momentum fluxes (both are generally unrealistic).
	1535	IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min ) &
	1536	surf%ol(m) = z_mo / zeta_min
	1537	IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max ) &
	1538	surf%ol(m) = z_mo / zeta_max
[1691]	1539
[2232]	1540	ENDDO
	1541	ELSE
	1542
[3146]	1543	!$OMP PARALLEL DO PRIVATE( z_mo )
[2232]	1544	DO m = 1, surf%ns
	1545
	1546	z_mo = surf%z_mo(m)
	1547
[3146]	1548	surf%ol(m) = ( surf%vpt1(m) * surf%us(m)**2 ) / &
	1549	( kappa * g * ( surf%ts(m) + 0.61_wp * surf%pt1(m) * &
	1550	surf%qs(m) + 0.61_wp * surf%qv1(m) * &
[2232]	1551	surf%ts(m) ) + 1E-30_wp )
	1552
	1553	!
	1554	!-- Limit the value range of the Obukhov length.
	1555	!-- This is necessary for very small velocities (u,v --> 1), because
	1556	!-- the absolute value of ol can then become very small, which in
	1557	!-- consequence would result in very large shear stresses and very
	1558	!-- small momentum fluxes (both are generally unrealistic).
	1559	IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) < zeta_min ) &
	1560	surf%ol(m) = z_mo / zeta_min
	1561	IF ( ( z_mo / ( surf%ol(m) + 1E-30_wp ) ) > zeta_max ) &
	1562	surf%ol(m) = z_mo / zeta_max
	1563
	1564	ENDDO
	1565
	1566	ENDIF
	1567
	1568	ENDIF
	1569
[1691]	1570	END SUBROUTINE calc_ol
	1571
	1572	!
	1573	!-- Calculate friction velocity u*
	1574	SUBROUTINE calc_us
	1575
	1576	IMPLICIT NONE
	1577
[2232]	1578	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
[1691]	1579
[2232]	1580	!
	1581	!-- Compute u* at horizontal surfaces at the scalars' grid points
	1582	IF ( .NOT. surf_vertical ) THEN
	1583	!
	1584	!-- Compute u* at upward-facing surfaces
	1585	IF ( .NOT. downward ) THEN
[2281]	1586	!$OMP PARALLEL DO PRIVATE( z_mo )
[3634]	1587	!$ACC PARALLEL LOOP PRIVATE(z_mo) &
	1588	!$ACC PRESENT(surf)
[2232]	1589	DO m = 1, surf%ns
[1691]	1590
[2232]	1591	z_mo = surf%z_mo(m)
[1691]	1592	!
[2232]	1593	!-- Compute u* at the scalars' grid points
	1594	surf%us(m) = kappa * surf%uvw_abs(m) / &
	1595	( LOG( z_mo / surf%z0(m) ) &
	1596	- psi_m( z_mo / surf%ol(m) ) &
	1597	+ psi_m( surf%z0(m) / surf%ol(m) ) )
	1598
	1599	ENDDO
	1600	!
	1601	!-- Compute u* at downward-facing surfaces. This case, do not consider
	1602	!-- any stability.
	1603	ELSE
[2281]	1604	!$OMP PARALLEL DO PRIVATE( z_mo )
[3634]	1605	!$ACC PARALLEL LOOP PRIVATE(z_mo) &
	1606	!$ACC PRESENT(surf)
[2232]	1607	DO m = 1, surf%ns
	1608
	1609	z_mo = surf%z_mo(m)
	1610	!
	1611	!-- Compute u* at the scalars' grid points
	1612	surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
	1613
	1614	ENDDO
	1615	ENDIF
	1616	!
	1617	!-- Compute u* at vertical surfaces at the u/v/v grid, respectively.
	1618	!-- No stability is considered in this case.
	1619	ELSE
	1620	!$OMP PARALLEL DO PRIVATE( z_mo )
[3634]	1621	!$ACC PARALLEL LOOP PRIVATE(z_mo) &
	1622	!$ACC PRESENT(surf)
[2232]	1623	DO m = 1, surf%ns
	1624	z_mo = surf%z_mo(m)
	1625
	1626	surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) )
[1691]	1627	ENDDO
[2232]	1628	ENDIF
[1691]	1629
	1630	END SUBROUTINE calc_us
	1631
	1632	!
[3146]	1633	!-- Calculate potential temperature, specific humidity, and virtual potential
	1634	!-- temperature at first grid level
[1691]	1635	SUBROUTINE calc_pt_q
	1636
	1637	IMPLICIT NONE
	1638
[2232]	1639	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
	1640
	1641	!$OMP PARALLEL DO PRIVATE( i, j, k )
[3634]	1642	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	1643	!$ACC PRESENT(surf, pt)
[2232]	1644	DO m = 1, surf%ns
	1645
	1646	i = surf%i(m)
	1647	j = surf%j(m)
	1648	k = surf%k(m)
	1649
[3634]	1650	#ifndef _OPENACC
[3274]	1651	IF ( bulk_cloud_model ) THEN
	1652	surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
[2547]	1653	surf%qv1(m) = q(k,j,i) - ql(k,j,i)
	1654	ELSEIF( cloud_droplets ) THEN
[3274]	1655	surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)
[2547]	1656	surf%qv1(m) = q(k,j,i)
[2696]	1657	ELSE
[3634]	1658	#endif
[2696]	1659	surf%pt1(m) = pt(k,j,i)
[3634]	1660	#ifndef _OPENACC
[2696]	1661	IF ( humidity ) THEN
	1662	surf%qv1(m) = q(k,j,i)
	1663	ELSE
[3634]	1664	#endif
[2696]	1665	surf%qv1(m) = 0.0_wp
[3634]	1666	#ifndef _OPENACC
[2696]	1667	ENDIF
[2547]	1668	ENDIF
[2232]	1669
[3146]	1670	IF ( humidity ) THEN
	1671	surf%vpt1(m) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
	1672	ENDIF
[3634]	1673	#endif
[3146]	1674
[1691]	1675	ENDDO
	1676
	1677	END SUBROUTINE calc_pt_q
	1678
[2696]	1679
[1691]	1680	!
[3152]	1681	!-- Set potential temperature at surface grid level.
[2696]	1682	!-- ( only for upward-facing surfs )
	1683	SUBROUTINE calc_pt_surface
	1684
	1685	IMPLICIT NONE
	1686
[3146]	1687	INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
[2696]	1688	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
	1689
[3146]	1690	k_off = surf%koff
[2696]	1691	!$OMP PARALLEL DO PRIVATE( i, j, k )
[3634]	1692	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	1693	!$ACC PRESENT(surf, pt)
[2696]	1694	DO m = 1, surf%ns
	1695
	1696	i = surf%i(m)
	1697	j = surf%j(m)
	1698	k = surf%k(m)
	1699
[3146]	1700	surf%pt_surface(m) = pt(k+k_off,j,i)
[2696]	1701
	1702	ENDDO
	1703
	1704	END SUBROUTINE calc_pt_surface
	1705
	1706	!
[3152]	1707	!-- Set mixing ratio at surface grid level. ( Only for upward-facing surfs. )
	1708	SUBROUTINE calc_q_surface
	1709
	1710	IMPLICIT NONE
	1711
	1712	INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
	1713	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
	1714
	1715	k_off = surf%koff
	1716	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1717	DO m = 1, surf%ns
	1718
	1719	i = surf%i(m)
	1720	j = surf%j(m)
	1721	k = surf%k(m)
	1722
	1723	surf%q_surface(m) = q(k+k_off,j,i)
	1724
	1725	ENDDO
	1726
	1727	END SUBROUTINE calc_q_surface
	1728
	1729	!
	1730	!-- Set virtual potential temperature at surface grid level.
[3146]	1731	!-- ( only for upward-facing surfs )
	1732	SUBROUTINE calc_vpt_surface
	1733
	1734	IMPLICIT NONE
	1735
	1736	INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
	1737	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
	1738
	1739	k_off = surf%koff
	1740	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1741	DO m = 1, surf%ns
	1742
	1743	i = surf%i(m)
	1744	j = surf%j(m)
	1745	k = surf%k(m)
	1746
	1747	surf%vpt_surface(m) = vpt(k+k_off,j,i)
	1748
	1749	ENDDO
	1750
	1751	END SUBROUTINE calc_vpt_surface
	1752
	1753	!
[2292]	1754	!-- Calculate the other MOST scaling parameters theta, q, (qc, qr, nc, nr)
[1691]	1755	SUBROUTINE calc_scaling_parameters
	1756
	1757	IMPLICIT NONE
	1758
[2232]	1759
	1760	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
[2696]	1761	INTEGER(iwp) :: lsp !< running index for chemical species
[1691]	1762	!
[2232]	1763	!-- Compute theta* at horizontal surfaces
	1764	IF ( constant_heatflux .AND. .NOT. surf_vertical ) THEN
[1691]	1765	!
	1766	!-- For a given heat flux in the surface layer:
[2232]	1767
	1768	!$OMP PARALLEL DO PRIVATE( i, j, k )
[3634]	1769	!$ACC PARALLEL LOOP PRIVATE(i, j, k) &
	1770	!$ACC PRESENT(surf, drho_air_zw)
[2232]	1771	DO m = 1, surf%ns
	1772
	1773	i = surf%i(m)
	1774	j = surf%j(m)
	1775	k = surf%k(m)
	1776
	1777	surf%ts(m) = -surf%shf(m) * drho_air_zw(k-1) / &
	1778	( surf%us(m) + 1E-30_wp )
	1779
[1691]	1780	!
[2232]	1781	!-- ts must be limited, because otherwise overflow may occur in case
	1782	!-- of us=0 when computing ol further below
	1783	IF ( surf%ts(m) < -1.05E5_wp ) surf%ts(m) = -1.0E5_wp
	1784	IF ( surf%ts(m) > 1.0E5_wp ) surf%ts(m) = 1.0E5_wp
	1785
[1691]	1786	ENDDO
	1787
[2232]	1788	ELSEIF ( .NOT. surf_vertical ) THEN
[1691]	1789	!
	1790	!-- For a given surface temperature:
[1788]	1791	IF ( large_scale_forcing .AND. lsf_surf ) THEN
[2232]	1792
	1793	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1794	DO m = 1, surf%ns
	1795	i = surf%i(m)
	1796	j = surf%j(m)
	1797	k = surf%k(m)
	1798
	1799	pt(k-1,j,i) = pt_surface
[1691]	1800	ENDDO
	1801	ENDIF
	1802
[2696]	1803	!$OMP PARALLEL DO PRIVATE( z_mo )
	1804	DO m = 1, surf%ns
[1691]	1805
[2696]	1806	z_mo = surf%z_mo(m)
[1691]	1807
[2696]	1808	surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) ) &
	1809	/ ( LOG( z_mo / surf%z0h(m) ) &
	1810	- psi_h( z_mo / surf%ol(m) ) &
	1811	+ psi_h( surf%z0h(m) / surf%ol(m) ) )
[1691]	1812
[2696]	1813	ENDDO
[2232]	1814
	1815	ENDIF
	1816	!
	1817	!-- Compute theta* at vertical surfaces. This is only required in case of
	1818	!-- land-surface model, in order to compute aerodynamical resistance.
	1819	IF ( surf_vertical ) THEN
[2281]	1820	!$OMP PARALLEL DO PRIVATE( i, j )
[2232]	1821	DO m = 1, surf%ns
	1822
	1823	i = surf%i(m)
	1824	j = surf%j(m)
	1825	surf%ts(m) = -surf%shf(m) / ( surf%us(m) + 1E-30_wp )
	1826	!
	1827	!-- ts must be limited, because otherwise overflow may occur in case
	1828	!-- of us=0 when computing ol further below
	1829	IF ( surf%ts(m) < -1.05E5_wp ) surf%ts(m) = -1.0E5_wp
	1830	IF ( surf%ts(m) > 1.0E5_wp ) surf%ts(m) = 1.0E5_wp
	1831
[1691]	1832	ENDDO
	1833	ENDIF
	1834
	1835	!
[2232]	1836	!-- If required compute q* at horizontal surfaces
[1960]	1837	IF ( humidity ) THEN
[2232]	1838	IF ( constant_waterflux .AND. .NOT. surf_vertical ) THEN
[1691]	1839	!
[1788]	1840	!-- For a given water flux in the surface layer
[2232]	1841	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1842	DO m = 1, surf%ns
	1843
	1844	i = surf%i(m)
	1845	j = surf%j(m)
	1846	k = surf%k(m)
	1847	surf%qs(m) = -surf%qsws(m) * drho_air_zw(k-1) / &
	1848	( surf%us(m) + 1E-30_wp )
	1849
[1691]	1850	ENDDO
	1851
[2232]	1852	ELSEIF ( .NOT. surf_vertical ) THEN
[1788]	1853	coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. &
[1691]	1854	run_coupled )
	1855
[1788]	1856	IF ( large_scale_forcing .AND. lsf_surf ) THEN
[2232]	1857	!$OMP PARALLEL DO PRIVATE( i, j, k )
	1858	DO m = 1, surf%ns
	1859
	1860	i = surf%i(m)
	1861	j = surf%j(m)
	1862	k = surf%k(m)
	1863	q(k-1,j,i) = q_surface
	1864
[1691]	1865	ENDDO
	1866	ENDIF
	1867
	1868	!
[2232]	1869	!-- Assume saturation for atmosphere coupled to ocean (but not
	1870	!-- in case of precursor runs)
	1871	IF ( coupled_run ) THEN
	1872	!$OMP PARALLEL DO PRIVATE( i, j, k, e_s )
	1873	DO m = 1, surf%ns
	1874	i = surf%i(m)
	1875	j = surf%j(m)
	1876	k = surf%k(m)
	1877	e_s = 6.1_wp * &
	1878	EXP( 0.07_wp * ( MIN(pt(k-1,j,i),pt(k,j,i)) &
[1691]	1879	- 273.15_wp ) )
[3361]	1880	q(k-1,j,i) = rd_d_rv * e_s / ( surface_pressure - e_s )
[2232]	1881	ENDDO
	1882	ENDIF
[1691]	1883
[3274]	1884	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
[2232]	1885	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	1886	DO m = 1, surf%ns
[1691]	1887
[2232]	1888	i = surf%i(m)
	1889	j = surf%j(m)
	1890	k = surf%k(m)
	1891
	1892	z_mo = surf%z_mo(m)
[1691]	1893
[3152]	1894	surf%qs(m) = kappa * ( surf%qv1(m) - surf%q_surface(m) ) &
[2232]	1895	/ ( LOG( z_mo / surf%z0q(m) ) &
	1896	- psi_h( z_mo / surf%ol(m) ) &
	1897	+ psi_h( surf%z0q(m) / &
	1898	surf%ol(m) ) )
[1691]	1899	ENDDO
[2232]	1900	ELSE
	1901	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	1902	DO m = 1, surf%ns
	1903
	1904	i = surf%i(m)
	1905	j = surf%j(m)
	1906	k = surf%k(m)
	1907
	1908	z_mo = surf%z_mo(m)
	1909
	1910	surf%qs(m) = kappa * ( q(k,j,i) - q(k-1,j,i) ) &
	1911	/ ( LOG( z_mo / surf%z0q(m) ) &
	1912	- psi_h( z_mo / surf%ol(m) ) &
	1913	+ psi_h( surf%z0q(m) / &
	1914	surf%ol(m) ) )
	1915	ENDDO
	1916	ENDIF
	1917	ENDIF
	1918	!
	1919	!-- Compute q* at vertical surfaces
	1920	IF ( surf_vertical ) THEN
[2281]	1921	!$OMP PARALLEL DO PRIVATE( i, j )
[2232]	1922	DO m = 1, surf%ns
	1923
	1924	i = surf%i(m)
	1925	j = surf%j(m)
	1926	surf%qs(m) = -surf%qsws(m) / ( surf%us(m) + 1E-30_wp )
	1927
[1691]	1928	ENDDO
	1929	ENDIF
	1930	ENDIF
[1960]	1931
	1932	!
	1933	!-- If required compute s*
	1934	IF ( passive_scalar ) THEN
	1935	!
[2232]	1936	!-- At horizontal surfaces
	1937	IF ( constant_scalarflux .AND. .NOT. surf_vertical ) THEN
	1938	!
	1939	!-- For a given scalar flux in the surface layer
[2281]	1940	!$OMP PARALLEL DO PRIVATE( i, j )
[2232]	1941	DO m = 1, surf%ns
	1942	i = surf%i(m)
	1943	j = surf%j(m)
	1944	surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
[1960]	1945	ENDDO
	1946	ENDIF
[2232]	1947	!
	1948	!-- At vertical surfaces
	1949	IF ( surf_vertical ) THEN
[2281]	1950	!$OMP PARALLEL DO PRIVATE( i, j )
[2232]	1951	DO m = 1, surf%ns
	1952	i = surf%i(m)
	1953	j = surf%j(m)
	1954	surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )
	1955	ENDDO
	1956	ENDIF
[1960]	1957	ENDIF
[1691]	1958
[2292]	1959	!
[2696]	1960	!-- If required compute cs* (chemical species)
	1961	IF ( air_chemistry ) THEN
	1962	!
	1963	!-- At horizontal surfaces
	1964	DO lsp = 1, nvar
	1965	IF ( constant_csflux(lsp) .AND. .NOT. surf_vertical ) THEN
	1966	!-- For a given chemical species' flux in the surface layer
	1967	!$OMP PARALLEL DO PRIVATE( i, j )
	1968	DO m = 1, surf%ns
	1969	i = surf%i(m)
	1970	j = surf%j(m)
	1971	surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
	1972	ENDDO
	1973	ENDIF
	1974	ENDDO
	1975	!
	1976	!-- At vertical surfaces
	1977	IF ( surf_vertical ) THEN
	1978	DO lsp = 1, nvar
	1979	!$OMP PARALLEL DO PRIVATE( i, j )
	1980	DO m = 1, surf%ns
	1981	i = surf%i(m)
	1982	j = surf%j(m)
	1983	surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
	1984	ENDDO
	1985	ENDDO
	1986	ENDIF
	1987	ENDIF
	1988
	1989	!
[2292]	1990	!-- If required compute qc* and nc*
[3274]	1991	IF ( bulk_cloud_model .AND. microphysics_morrison .AND. &
[2292]	1992	.NOT. surf_vertical ) THEN
	1993	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	1994	DO m = 1, surf%ns
[1691]	1995
[2292]	1996	i = surf%i(m)
	1997	j = surf%j(m)
	1998	k = surf%k(m)
	1999
	2000	z_mo = surf%z_mo(m)
	2001
	2002	surf%qcs(m) = kappa * ( qc(k,j,i) - qc(k-1,j,i) ) &
	2003	/ ( LOG( z_mo / surf%z0q(m) ) &
	2004	- psi_h( z_mo / surf%ol(m) ) &
	2005	+ psi_h( surf%z0q(m) / surf%ol(m) ) )
	2006
	2007	surf%ncs(m) = kappa * ( nc(k,j,i) - nc(k-1,j,i) ) &
	2008	/ ( LOG( z_mo / surf%z0q(m) ) &
	2009	- psi_h( z_mo / surf%ol(m) ) &
	2010	+ psi_h( surf%z0q(m) / surf%ol(m) ) )
	2011	ENDDO
	2012
	2013	ENDIF
	2014
[1691]	2015	!
	2016	!-- If required compute qr* and nr*
[3274]	2017	IF ( bulk_cloud_model .AND. microphysics_seifert .AND. &
[2232]	2018	.NOT. surf_vertical ) THEN
	2019	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	2020	DO m = 1, surf%ns
[1691]	2021
[2232]	2022	i = surf%i(m)
	2023	j = surf%j(m)
	2024	k = surf%k(m)
[1691]	2025
[2232]	2026	z_mo = surf%z_mo(m)
[1691]	2027
[2232]	2028	surf%qrs(m) = kappa * ( qr(k,j,i) - qr(k-1,j,i) ) &
	2029	/ ( LOG( z_mo / surf%z0q(m) ) &
	2030	- psi_h( z_mo / surf%ol(m) ) &
	2031	+ psi_h( surf%z0q(m) / surf%ol(m) ) )
[1691]	2032
[2232]	2033	surf%nrs(m) = kappa * ( nr(k,j,i) - nr(k-1,j,i) ) &
	2034	/ ( LOG( z_mo / surf%z0q(m) ) &
	2035	- psi_h( z_mo / surf%ol(m) ) &
	2036	+ psi_h( surf%z0q(m) / surf%ol(m) ) )
[1691]	2037	ENDDO
	2038
	2039	ENDIF
	2040
	2041	END SUBROUTINE calc_scaling_parameters
	2042
	2043
	2044
	2045	!
[2292]	2046	!-- Calculate surface fluxes usws, vsws, shf, qsws, (qcsws, qrsws, ncsws, nrsws)
[1691]	2047	SUBROUTINE calc_surface_fluxes
	2048
	2049	IMPLICIT NONE
	2050
[2696]	2051	INTEGER(iwp) :: m !< loop variable over all horizontal surf elements
	2052	INTEGER(iwp) :: lsp !< running index for chemical species
[1691]	2053
[2232]	2054	REAL(wp) :: dum !< dummy to precalculate logarithm
	2055	REAL(wp) :: flag_u !< flag indicating u-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
	2056	REAL(wp) :: flag_v !< flag indicating v-grid, used for calculation of horizontal momentum fluxes at vertical surfaces
	2057	REAL(wp), DIMENSION(:), ALLOCATABLE :: u_i !< u-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
	2058	REAL(wp), DIMENSION(:), ALLOCATABLE :: v_i !< v-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
	2059	REAL(wp), DIMENSION(:), ALLOCATABLE :: w_i !< w-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces
[1691]	2060
	2061	!
[2232]	2062	!-- Calcuate surface fluxes at horizontal walls
	2063	IF ( .NOT. surf_vertical ) THEN
	2064	!
	2065	!-- Compute u'w' for the total model domain at upward-facing surfaces.
	2066	!-- First compute the corresponding component of u* and square it.
	2067	IF ( .NOT. downward ) THEN
	2068	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
[3634]	2069	!$ACC PARALLEL LOOP PRIVATE(i, j, k, z_mo) &
	2070	!$ACC PRESENT(surf, u, rho_air_zw)
[2232]	2071	DO m = 1, surf%ns
	2072
	2073	i = surf%i(m)
	2074	j = surf%j(m)
	2075	k = surf%k(m)
[1691]	2076
[2232]	2077	z_mo = surf%z_mo(m)
[1691]	2078
[2232]	2079	surf%usws(m) = kappa * ( u(k,j,i) - u(k-1,j,i) ) &
	2080	/ ( LOG( z_mo / surf%z0(m) ) &
	2081	- psi_m( z_mo / surf%ol(m) ) &
	2082	+ psi_m( surf%z0(m) / surf%ol(m) ) )
	2083	!
	2084	!-- Please note, the computation of usws is not fully accurate. Actually
	2085	!-- a further interpolation of us onto the u-grid, where usws is defined,
	2086	!-- is required. However, this is not done as this would require several
	2087	!-- data transfers between 2D-grid and the surf-type.
	2088	!-- The impact of the missing interpolation is negligible as several
	2089	!-- tests had shown.
	2090	!-- Same also for ol.
	2091	surf%usws(m) = -surf%usws(m) * surf%us(m) * rho_air_zw(k-1)
[1691]	2092
[2232]	2093	ENDDO
[1691]	2094	!
[2232]	2095	!-- At downward-facing surfaces
	2096	ELSE
	2097	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	2098	DO m = 1, surf%ns
	2099
	2100	i = surf%i(m)
	2101	j = surf%j(m)
	2102	k = surf%k(m)
[1691]	2103
[2232]	2104	z_mo = surf%z_mo(m)
	2105
	2106	surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) )
	2107	surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k)
[1691]	2108
[2232]	2109	ENDDO
	2110	ENDIF
[1691]	2111
[2232]	2112	!
	2113	!-- Compute v'w' for the total model domain.
	2114	!-- First compute the corresponding component of u* and square it.
	2115	!-- Upward-facing surfaces
	2116	IF ( .NOT. downward ) THEN
	2117	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
[3634]	2118	!$ACC PARALLEL LOOP PRIVATE(i, j, k, z_mo) &
	2119	!$ACC PRESENT(surf, v, rho_air_zw)
[2232]	2120	DO m = 1, surf%ns
	2121	i = surf%i(m)
	2122	j = surf%j(m)
	2123	k = surf%k(m)
[1691]	2124
[2232]	2125	z_mo = surf%z_mo(m)
[1691]	2126
[2232]	2127	surf%vsws(m) = kappa * ( v(k,j,i) - v(k-1,j,i) ) &
	2128	/ ( LOG( z_mo / surf%z0(m) ) &
	2129	- psi_m( z_mo / surf%ol(m) ) &
	2130	+ psi_m( surf%z0(m) / surf%ol(m) ) )
[1691]	2131	!
[2232]	2132	!-- Please note, the computation of vsws is not fully accurate. Actually
	2133	!-- a further interpolation of us onto the v-grid, where vsws is defined,
	2134	!-- is required. However, this is not done as this would require several
	2135	!-- data transfers between 2D-grid and the surf-type.
	2136	!-- The impact of the missing interpolation is negligible as several
	2137	!-- tests had shown.
	2138	!-- Same also for ol.
	2139	surf%vsws(m) = -surf%vsws(m) * surf%us(m) * rho_air_zw(k-1)
	2140	ENDDO
	2141	!
	2142	!-- Downward-facing surfaces
	2143	ELSE
	2144	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	2145	DO m = 1, surf%ns
	2146	i = surf%i(m)
	2147	j = surf%j(m)
	2148	k = surf%k(m)
[1691]	2149
[2232]	2150	z_mo = surf%z_mo(m)
	2151
	2152	surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) )
	2153	surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k)
	2154	ENDDO
	2155	ENDIF
[1691]	2156	!
[2232]	2157	!-- Compute the vertical kinematic heat flux
[2299]	2158	IF ( .NOT. constant_heatflux .AND. ( ( time_since_reference_point&
	2159	<= skip_time_do_lsm .AND. simulated_time > 0.0_wp ) .OR. &
[2696]	2160	.NOT. land_surface ) .AND. .NOT. urban_surface .AND. &
[2299]	2161	.NOT. downward ) THEN
[2232]	2162	!$OMP PARALLEL DO PRIVATE( i, j, k )
	2163	DO m = 1, surf%ns
	2164	i = surf%i(m)
	2165	j = surf%j(m)
	2166	k = surf%k(m)
	2167	surf%shf(m) = -surf%ts(m) * surf%us(m) * rho_air_zw(k-1)
[1691]	2168	ENDDO
[2232]	2169	ENDIF
	2170	!
	2171	!-- Compute the vertical water flux
	2172	IF ( .NOT. constant_waterflux .AND. humidity .AND. &
[2299]	2173	( ( time_since_reference_point <= skip_time_do_lsm .AND. &
[2696]	2174	simulated_time > 0.0_wp ) .OR. .NOT. land_surface ) .AND. &
	2175	.NOT. urban_surface .AND. .NOT. downward ) THEN
[2232]	2176	!$OMP PARALLEL DO PRIVATE( i, j, k )
	2177	DO m = 1, surf%ns
	2178	i = surf%i(m)
	2179	j = surf%j(m)
	2180	k = surf%k(m)
	2181	surf%qsws(m) = -surf%qs(m) * surf%us(m) * rho_air_zw(k-1)
	2182	ENDDO
	2183	ENDIF
	2184	!
	2185	!-- Compute the vertical scalar flux
	2186	IF ( .NOT. constant_scalarflux .AND. passive_scalar .AND. &
	2187	.NOT. downward ) THEN
[2281]	2188	!$OMP PARALLEL DO PRIVATE( i, j )
[2232]	2189	DO m = 1, surf%ns
[1691]	2190
[2232]	2191	i = surf%i(m)
	2192	j = surf%j(m)
	2193	surf%ssws(m) = -surf%ss(m) * surf%us(m)
[1691]	2194
[2232]	2195	ENDDO
[2292]	2196	ENDIF
[1691]	2197	!
[2696]	2198	!-- Compute the vertical chemical species' flux
	2199	DO lsp = 1, nvar
	2200	IF ( .NOT. constant_csflux(lsp) .AND. air_chemistry .AND. &
	2201	.NOT. downward ) THEN
	2202	!$OMP PARALLEL DO PRIVATE( i, j )
	2203	DO m = 1, surf%ns
	2204	i = surf%i(m)
	2205	j = surf%j(m)
	2206	surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m)
	2207	ENDDO
	2208	ENDIF
	2209	ENDDO
	2210
	2211	!
[2292]	2212	!-- Compute (turbulent) fluxes of cloud water content and cloud drop conc.
[3274]	2213	IF ( bulk_cloud_model .AND. microphysics_morrison .AND. &
[2292]	2214	.NOT. downward) THEN
	2215	!$OMP PARALLEL DO PRIVATE( i, j )
	2216	DO m = 1, surf%ns
	2217
	2218	i = surf%i(m)
	2219	j = surf%j(m)
	2220
	2221	surf%qcsws(m) = -surf%qcs(m) * surf%us(m)
	2222	surf%ncsws(m) = -surf%ncs(m) * surf%us(m)
	2223	ENDDO
	2224	ENDIF
	2225	!
[2232]	2226	!-- Compute (turbulent) fluxes of rain water content and rain drop conc.
[3274]	2227	IF ( bulk_cloud_model .AND. microphysics_seifert .AND. &
[2232]	2228	.NOT. downward) THEN
[2281]	2229	!$OMP PARALLEL DO PRIVATE( i, j )
[2232]	2230	DO m = 1, surf%ns
	2231
	2232	i = surf%i(m)
	2233	j = surf%j(m)
	2234
	2235	surf%qrsws(m) = -surf%qrs(m) * surf%us(m)
	2236	surf%nrsws(m) = -surf%nrs(m) * surf%us(m)
[1691]	2237	ENDDO
[2232]	2238	ENDIF
[1691]	2239
[1960]	2240	!
[2232]	2241	!-- Bottom boundary condition for the TKE.
	2242	IF ( ibc_e_b == 2 ) THEN
	2243	!$OMP PARALLEL DO PRIVATE( i, j, k )
	2244	DO m = 1, surf%ns
	2245
	2246	i = surf%i(m)
	2247	j = surf%j(m)
	2248	k = surf%k(m)
	2249
	2250	e(k,j,i) = ( surf%us(m) / 0.1_wp )**2
	2251	!
	2252	!-- As a test: cm = 0.4
	2253	! e(k,j,i) = ( us(j,i) / 0.4_wp )**2
	2254	e(k-1,j,i) = e(k,j,i)
	2255
[1960]	2256	ENDDO
[2232]	2257	ENDIF
	2258	!
	2259	!-- Calcuate surface fluxes at vertical surfaces. No stability is considered.
	2260	ELSE
	2261	!
	2262	!-- Compute usvs l={0,1} and vsus l={2,3}
	2263	IF ( mom_uv ) THEN
	2264	!
	2265	!-- Generalize computation by introducing flags. At north- and south-
	2266	!-- facing surfaces u-component is used, at east- and west-facing
	2267	!-- surfaces v-component is used.
	2268	flag_u = MERGE( 1.0_wp, 0.0_wp, l == 0 .OR. l == 1 )
	2269	flag_v = MERGE( 1.0_wp, 0.0_wp, l == 2 .OR. l == 3 )
	2270	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	2271	DO m = 1, surf%ns
	2272	i = surf%i(m)
	2273	j = surf%j(m)
	2274	k = surf%k(m)
[1691]	2275
[2232]	2276	z_mo = surf%z_mo(m)
[1960]	2277
[2232]	2278	surf%mom_flux_uv(m) = kappa * &
	2279	( flag_u * u(k,j,i) + flag_v * v(k,j,i) ) / &
	2280	LOG( z_mo / surf%z0(m) )
	2281
	2282	surf%mom_flux_uv(m) = &
	2283	- surf%mom_flux_uv(m) * surf%us(m)
	2284	ENDDO
	2285	ENDIF
[1691]	2286	!
[2232]	2287	!-- Compute wsus l={0,1} and wsvs l={2,3}
	2288	IF ( mom_w ) THEN
	2289	!$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
	2290	DO m = 1, surf%ns
	2291	i = surf%i(m)
	2292	j = surf%j(m)
	2293	k = surf%k(m)
	2294
	2295	z_mo = surf%z_mo(m)
	2296
	2297	surf%mom_flux_w(m) = kappa * w(k,j,i) / LOG( z_mo / surf%z0(m) )
	2298
	2299	surf%mom_flux_w(m) = &
	2300	- surf%mom_flux_w(m) * surf%us(m)
[1691]	2301	ENDDO
[2232]	2302	ENDIF
	2303	!
	2304	!-- Compute momentum fluxes used for subgrid-scale TKE production at
	2305	!-- vertical surfaces. In constrast to the calculated momentum fluxes at
	2306	!-- vertical surfaces before, which are defined on the u/v/w-grid,
	2307	!-- respectively), the TKE fluxes are defined at the scalar grid.
	2308	!--
	2309	IF ( mom_tke ) THEN
	2310	!
	2311	!-- Precalculate velocity components at scalar grid point.
	2312	ALLOCATE( u_i(1:surf%ns) )
	2313	ALLOCATE( v_i(1:surf%ns) )
	2314	ALLOCATE( w_i(1:surf%ns) )
[1691]	2315
[2232]	2316	IF ( l == 0 .OR. l == 1 ) THEN
	2317	!$OMP PARALLEL DO PRIVATE( i, j, k )
	2318	DO m = 1, surf%ns
	2319	i = surf%i(m)
	2320	j = surf%j(m)
	2321	k = surf%k(m)
	2322
	2323	u_i(m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
	2324	v_i(m) = 0.0_wp
	2325	w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
	2326	ENDDO
	2327	ELSE
	2328	!$OMP PARALLEL DO PRIVATE( i, j, k )
	2329	DO m = 1, surf%ns
	2330	i = surf%i(m)
	2331	j = surf%j(m)
	2332	k = surf%k(m)
	2333
	2334	u_i(m) = 0.0_wp
	2335	v_i(m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
	2336	w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
	2337	ENDDO
	2338	ENDIF
	2339
[2281]	2340	!$OMP PARALLEL DO PRIVATE( i, j, dum, z_mo )
[2232]	2341	DO m = 1, surf%ns
	2342	i = surf%i(m)
	2343	j = surf%j(m)
	2344
	2345	z_mo = surf%z_mo(m)
	2346
	2347	dum = kappa / LOG( z_mo / surf%z0(m) )
[1691]	2348	!
[2232]	2349	!-- usvs (l=0,1) and vsus (l=2,3)
	2350	surf%mom_flux_tke(0,m) = dum * ( u_i(m) + v_i(m) )
[1691]	2351	!
[2232]	2352	!-- wsvs (l=0,1) and wsus (l=2,3)
	2353	surf%mom_flux_tke(1,m) = dum * w_i(m)
	2354
	2355	surf%mom_flux_tke(0:1,m) = &
	2356	- surf%mom_flux_tke(0:1,m) * surf%us(m)
[1691]	2357	ENDDO
[2232]	2358	!
	2359	!-- Deallocate temporary arrays
	2360	DEALLOCATE( u_i )
	2361	DEALLOCATE( v_i )
	2362	DEALLOCATE( w_i )
	2363	ENDIF
[1691]	2364	ENDIF
	2365
	2366	END SUBROUTINE calc_surface_fluxes
	2367
[3597]	2368
	2369	!------------------------------------------------------------------------------!
	2370	! Description:
	2371	! ------------
	2372	!> Calculates temperature near surface (10 cm) for indoor model or 2 m
	2373	!> temperature for output
	2374	!------------------------------------------------------------------------------!
	2375	SUBROUTINE calc_pt_near_surface ( z_char )
[1691]	2376
[3597]	2377	IMPLICIT NONE
	2378
	2379	CHARACTER (LEN = *), INTENT(IN) :: z_char !< string identifier to identify z level
	2380	INTEGER(iwp) :: i, j, k, l, m !< running indices
	2381
	2382
	2383	SELECT CASE ( z_char)
	2384
	2385
	2386	CASE ( '10cm' )
	2387
	2388	DO m = 1, surf%ns
	2389
	2390	i = surf%i(m)
	2391	j = surf%j(m)
	2392	k = surf%k(m)
	2393
	2394	surf%pt_10cm(m) = surf%pt_surface(m) + surf%ts(m) / kappa &
	2395	* ( log( 0.1_wp / surf%z0h(m) ) &
	2396	- psi_h( 0.1_wp / surf%ol(m) ) &
	2397	+ psi_h( surf%z0h(m) / surf%ol(m) ) )
	2398
	2399	ENDDO
	2400
	2401
	2402	CASE ( '2m' )
	2403
	2404	DO m = 1, surf%ns
	2405
	2406	i = surf%i(m)
	2407	j = surf%j(m)
	2408	k = surf%k(m)
	2409
	2410	surf%pt_2m(m) = surf%pt_surface(m) + surf%ts(m) / kappa &
	2411	* ( log( 2.0_wp / surf%z0h(m) ) &
	2412	- psi_h( 2.0_wp / surf%ol(m) ) &
	2413	+ psi_h( surf%z0h(m) / surf%ol(m) ) )
	2414
	2415	ENDDO
	2416
	2417
	2418	END SELECT
	2419
	2420	END SUBROUTINE calc_pt_near_surface
	2421
	2422
[1691]	2423	!
	2424	!-- Integrated stability function for momentum
	2425	FUNCTION psi_m( zeta )
[3634]	2426	!$ACC ROUTINE SEQ
[1691]	2427
	2428	USE kinds
	2429
	2430	IMPLICIT NONE
	2431
	2432	REAL(wp) :: psi_m !< Integrated similarity function result
	2433	REAL(wp) :: zeta !< Stability parameter z/L
	2434	REAL(wp) :: x !< dummy variable
	2435
	2436	REAL(wp), PARAMETER :: a = 1.0_wp !< constant
	2437	REAL(wp), PARAMETER :: b = 0.66666666666_wp !< constant
	2438	REAL(wp), PARAMETER :: c = 5.0_wp !< constant
	2439	REAL(wp), PARAMETER :: d = 0.35_wp !< constant
	2440	REAL(wp), PARAMETER :: c_d_d = c / d !< constant
	2441	REAL(wp), PARAMETER :: bc_d_d = b * c / d !< constant
	2442
	2443
	2444	IF ( zeta < 0.0_wp ) THEN
[1788]	2445	x = SQRT( SQRT( 1.0_wp - 16.0_wp * zeta ) )
[1691]	2446	psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2 &
	2447	* ( 1.0_wp + x*2 ) 0.125_wp )
	2448	ELSE
	2449
	2450	psi_m = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - a * zeta &
	2451	- bc_d_d
	2452	!
	2453	!-- Old version for stable conditions (only valid for z/L < 0.5)
	2454	!-- psi_m = - 5.0_wp * zeta
	2455
	2456	ENDIF
	2457
	2458	END FUNCTION psi_m
	2459
	2460
	2461	!
	2462	!-- Integrated stability function for heat and moisture
	2463	FUNCTION psi_h( zeta )
[3634]	2464	!$ACC ROUTINE SEQ
[1691]	2465
	2466	USE kinds
	2467
	2468	IMPLICIT NONE
	2469
	2470	REAL(wp) :: psi_h !< Integrated similarity function result
	2471	REAL(wp) :: zeta !< Stability parameter z/L
	2472	REAL(wp) :: x !< dummy variable
	2473
	2474	REAL(wp), PARAMETER :: a = 1.0_wp !< constant
	2475	REAL(wp), PARAMETER :: b = 0.66666666666_wp !< constant
	2476	REAL(wp), PARAMETER :: c = 5.0_wp !< constant
	2477	REAL(wp), PARAMETER :: d = 0.35_wp !< constant
	2478	REAL(wp), PARAMETER :: c_d_d = c / d !< constant
	2479	REAL(wp), PARAMETER :: bc_d_d = b * c / d !< constant
	2480
	2481
	2482	IF ( zeta < 0.0_wp ) THEN
[1788]	2483	x = SQRT( 1.0_wp - 16.0_wp * zeta )
[1691]	2484	psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
	2485	ELSE
	2486	psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp &
	2487	+ 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d &
	2488	+ 1.0_wp
	2489	!
	2490	!-- Old version for stable conditions (only valid for z/L < 0.5)
	2491	!-- psi_h = - 5.0_wp * zeta
	2492	ENDIF
	2493
	2494	END FUNCTION psi_h
	2495
[3130]	2496
	2497	!------------------------------------------------------------------------------!
	2498	! Description:
	2499	! ------------
	2500	!> Calculates stability function for momentum
	2501	!>
	2502	!> @author Hauke Wurps
	2503	!------------------------------------------------------------------------------!
	2504	FUNCTION phi_m( zeta )
[3634]	2505	!$ACC ROUTINE SEQ
[3130]	2506
	2507	IMPLICIT NONE
	2508
	2509	REAL(wp) :: phi_m !< Value of the function
	2510	REAL(wp) :: zeta !< Stability parameter z/L
	2511
	2512	REAL(wp), PARAMETER :: a = 16.0_wp !< constant
	2513	REAL(wp), PARAMETER :: c = 5.0_wp !< constant
	2514
	2515	IF ( zeta < 0.0_wp ) THEN
	2516	phi_m = 1.0_wp / SQRT( SQRT( 1.0_wp - a * zeta ) )
	2517	ELSE
	2518	phi_m = 1.0_wp + c * zeta
	2519	ENDIF
	2520
	2521	END FUNCTION phi_m
	2522
[1697]	2523	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 |