Home

Context Navigation

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

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