Home

Context Navigation

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

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