Home

Context Navigation

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

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