Home

Context Navigation

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

Last change on this file since 3691 was 3685, checked in by knoop, 6 years ago
Some interface calls moved to module_interface + cleanup
Property svn:keywords set to `Id`
File size: 75.6 KB

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