Home

Context Navigation

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

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