Home

Context Navigation

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

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