Home

Context Navigation

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

Last change on this file since 4586 was 4562, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 78.1 KB

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