Home

Context Navigation

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

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