Home

Context Navigation

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

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