Home

Context Navigation

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

Last change on this file since 3790 was 3787, checked in by raasch, 6 years ago
unused variables removed
Property svn:keywords set to `Id`
File size: 74.9 KB

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