Home

Context Navigation

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

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