Home

Context Navigation

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

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