Home

Context Navigation

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

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