Home

Context Navigation

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

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