Home

Context Navigation

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

Last change on this file since 4366 was 4366, checked in by raasch, 4 years ago
code vectorization for NEC Aurora: vectorized version of Temperton FFT, vectorization of Newtor iteration for calculating the Obukhov length
Property svn:keywords set to `Id`
File size: 78.2 KB

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