Home

Context Navigation

surface_layer_fluxes_mod.f90 @ 4851

Last change on this file since 4851 was 4828, checked in by Giersch, 4 years ago
Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code
Property svn:keywords set to `Id`
File size: 83.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: