Home

Context Navigation

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

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