Home

Context Navigation

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

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