Home

Context Navigation

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

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