Home

Context Navigation

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

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