Home

Context Navigation

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

Last change on this file since 2675 was 2547, checked in by schwenkel, 7 years ago
extended by cloud_droplets option
Property svn:keywords set to `Id`
File size: 80.8 KB

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