Home

Context Navigation

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

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