Home

Context Navigation

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

Last change on this file since 3614 was 3597, checked in by maronga, 6 years ago
revised calculation of near surface air potential temperature
Property svn:keywords set to `Id`
File size: 89.8 KB

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