Home

Context Navigation

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

Last change on this file since 2701 was 2696, checked in by kanani, 7 years ago
Merge of branch palm4u into trunk
Property svn:keywords set to `Id`
File size: 83.4 KB

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