Home

Context Navigation

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

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