Home

Context Navigation

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

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