Home

Context Navigation

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

Last change on this file since 2303 was 2299, checked in by maronga, 7 years ago
improvements for spinup mechanism
Property svn:keywords set to `Id`
File size: 80.3 KB

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