Home

Context Navigation

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

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