Home

Context Navigation

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

Last change on this file since 3676 was 3668, checked in by maronga, 6 years ago
removed most_methods circular and lookup. added improved version of palm_csd
Property svn:keywords set to `Id`
File size: 75.4 KB

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