Home

Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |