Home

Context Navigation

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

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