Home

Context Navigation

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

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