Home

Context Navigation

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

Last change on this file since 2004 was 2001, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 44.5 KB

Line
1	!> @file surface_layer_fluxes_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2016 Leibniz Universitaet Hannover
18	!
19	!------------------------------------------------------------------------------!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: surface_layer_fluxes_mod.f90 2001 2016-08-20 18:41:22Z suehring $
27	!
28	! 2000 2016-08-20 18:09:15Z knoop
29	! Forced header and separation lines into 80 columns
30	!
31	! 1992 2016-08-12 15:14:59Z suehring
32	! Minor bug, declaration of look-up index as INTEGER
33	!
34	! 1960 2016-07-12 16:34:24Z suehring
35	! Treat humidity and passive scalar separately
36	!
37	! 1929 2016-06-09 16:25:25Z suehring
38	! Bugfix: avoid segmentation fault in case one grid point is horizontally
39	! completely surrounded by topography
40	!
41	! 1920 2016-05-30 10:50:15Z suehring
42	! Avoid segmentation fault (see change in 1915) by different initialization of
43	! us instead of adding a very small number in the denominator
44	!
45	! 1915 2016-05-27 11:05:02Z suehring
46	! Bugfix: avoid segmentation fault in case of most_method = 'circular' at first
47	! timestep
48	!
49	! 1850 2016-04-08 13:29:27Z maronga
50	! Module renamed
51	!
52	!
53	! 1822 2016-04-07 07:49:42Z hoffmann
54	! icloud_scheme replaced by microphysics_*
55	!
56	! 1788 2016-03-10 11:01:04Z maronga
57	! Added parameter z0q which replaces z0h in the similarity functions for
58	! humidity.
59	! Syntax layout improved.
60	!
61	! 1757 2016-02-22 15:49:32Z maronga
62	! Minor fixes.
63	!
64	! 1749 2016-02-09 12:19:56Z raasch
65	! further OpenACC adjustments
66	!
67	! 1747 2016-02-08 12:25:53Z raasch
68	! adjustments for OpenACC usage
69	!
70	! 1709 2015-11-04 14:47:01Z maronga
71	! Bugfix: division by zero could occur when calculating rib at low wind speeds
72	! Bugfix: calculation of uv_total for neutral = .T., initial value for ol for
73	! neutral = .T.
74	!
75	! 1705 2015-11-02 14:28:56Z maronga
76	! Typo removed
77	!
78	! 1697 2015-10-28 17:14:10Z raasch
79	! FORTRAN and OpenMP errors removed
80	!
81	! 1696 2015-10-27 10:03:34Z maronga
82	! Modularized and completely re-written version of prandtl_fluxes.f90. In the
83	! course of the re-writing two additional methods have been implemented. See
84	! updated description.
85	!
86	! 1551 2015-03-03 14:18:16Z maronga
87	! Removed land surface model part. The surface fluxes are now always calculated
88	! within prandtl_fluxes, based on the given surface temperature/humidity (which
89	! is either provided by the land surface model, by large scale forcing data, or
90	! directly prescribed by the user.
91	!
92	! 1496 2014-12-02 17:25:50Z maronga
93	! Adapted for land surface model
94	!
95	! 1494 2014-11-21 17:14:03Z maronga
96	! Bugfixes: qs is now calculated before calculation of Rif. calculation of
97	! buoyancy flux in Rif corrected (added missing humidity term), allow use of
98	! topography for coupled runs (not tested)
99	!
100	! 1361 2014-04-16 15:17:48Z hoffmann
101	! Bugfix: calculation of turbulent fluxes of rain water content (qrsws) and rain
102	! drop concentration (nrsws) added
103	!
104	! 1340 2014-03-25 19:45:13Z kanani
105	! REAL constants defined as wp-kind
106	!
107	! 1320 2014-03-20 08:40:49Z raasch
108	! ONLY-attribute added to USE-statements,
109	! kind-parameters added to all INTEGER and REAL declaration statements,
110	! kinds are defined in new module kinds,
111	! old module precision_kind is removed,
112	! revision history before 2012 removed,
113	! comment fields (!:) to be used for variable explanations added to
114	! all variable declaration statements
115	!
116	! 1276 2014-01-15 13:40:41Z heinze
117	! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
118	!
119	! 1257 2013-11-08 15:18:40Z raasch
120	! openACC "kernels do" replaced by "kernels loop", "loop independent" added
121	!
122	! 1036 2012-10-22 13:43:42Z raasch
123	! code put under GPL (PALM 3.9)
124	!
125	! 1015 2012-09-27 09:23:24Z raasch
126	! OpenACC statements added
127	!
128	! 978 2012-08-09 08:28:32Z fricke
129	! roughness length for scalar quantities z0h added
130	!
131	! Revision 1.1 1998/01/23 10:06:06 raasch
132	! Initial revision
133	!
134	!
135	! Description:
136	! ------------
137	!> Diagnostic computation of vertical fluxes in the constant flux layer from the
138	!> values of the variables at grid point k=1. Three different methods are
139	!> available:
140	!> 1) the "old" version (most_method = 'circular') which is fast, but inaccurate
141	!> 2) a Newton iteration method (most_method = 'newton'), which is accurate, but
142	!> slower
143	!> 3) a method using a lookup table which is fast and accurate. Note, however,
144	!> that this method cannot be used in case of roughness heterogeneity
145	!>
146	!> @todo (re)move large_scale_forcing actions
147	!> @todo check/optimize OpenMP and OpenACC directives
148	!------------------------------------------------------------------------------!
149	MODULE surface_layer_fluxes_mod
150
151	USE arrays_3d, &
152	ONLY: e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, &
153	s, shf, ss, ssws, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, &
154	z0h, z0q
155
156	USE cloud_parameters, &
157	ONLY: l_d_cp, pt_d_t
158
159	USE constants, &
160	ONLY: pi
161
162	USE cpulog
163
164	USE control_parameters, &
165	ONLY: cloud_physics, constant_heatflux, constant_scalarflux, &
166	constant_waterflux, coupling_mode, g, humidity, ibc_e_b, &
167	ibc_pt_b, initializing_actions, kappa, &
168	intermediate_timestep_count, &
169	intermediate_timestep_count_max, large_scale_forcing, lsf_surf, &
170	message_string, microphysics_seifert, most_method, neutral, &
171	passive_scalar, pt_surface, q_surface, run_coupled, &
172	surface_pressure, simulated_time, terminate_run, zeta_max, &
173	zeta_min
174
175	USE indices, &
176	ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_s_inner, &
177	nzb_u_inner, nzb_v_inner
178
179	USE kinds
180
181	USE pegrid
182
183	USE land_surface_model_mod, &
184	ONLY: land_surface, skip_time_do_lsm
185
186	IMPLICIT NONE
187
188	INTEGER(iwp) :: i !< loop index x direction
189	INTEGER(iwp) :: j !< loop index y direction
190	INTEGER(iwp) :: k !< loop index z direction
191	INTEGER(iwp) :: l_bnd = 7500 !< Lookup table index of the last time step
192
193	INTEGER(iwp), PARAMETER :: num_steps = 15000 !< number of steps in the lookup table
194
195	LOGICAL :: coupled_run !< Flag for coupled atmosphere-ocean runs
196
197	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt1, & !< Potential temperature at first grid level (required for cloud_physics = .T.)
198	qv1, & !< Specific humidity at first grid level (required for cloud_physics = .T.)
199	uv_total !< Total velocity at first grid level
200
201	REAL(wp), DIMENSION(0:num_steps-1) :: rib_tab, & !< Lookup table bulk Richardson number
202	ol_tab !< Lookup table values of L
203
204	REAL(wp) :: e_s, & !< Saturation water vapor pressure
205	ol_max = 1.0E6_wp, & !< Maximum Obukhov length
206	rib_max, & !< Maximum Richardson number in lookup table
207	rib_min, & !< Minimum Richardson number in lookup table
208	z_mo !< Height of the constant flux layer where MOST is assumed
209
210
211	SAVE
212
213	PRIVATE
214
215	PUBLIC init_surface_layer_fluxes, pt1, qv1, surface_layer_fluxes, uv_total
216
217	INTERFACE init_surface_layer_fluxes
218	MODULE PROCEDURE init_surface_layer_fluxes
219	END INTERFACE init_surface_layer_fluxes
220
221	INTERFACE surface_layer_fluxes
222	MODULE PROCEDURE surface_layer_fluxes
223	END INTERFACE surface_layer_fluxes
224
225
226	CONTAINS
227
228
229	!------------------------------------------------------------------------------!
230	! Description:
231	! ------------
232	!> Main routine to compute the surface fluxes
233	!------------------------------------------------------------------------------!
234	SUBROUTINE surface_layer_fluxes
235
236	IMPLICIT NONE
237
238	!
239	!-- In case cloud physics is used, it is required to derive potential
240	!-- temperature and specific humidity at first grid level from the fields pt
241	!-- and q
242	IF ( cloud_physics ) THEN
243	CALL calc_pt_q
244	ENDIF
245
246	!
247	!-- First, calculate the new Obukhov length, then new friction velocity,
248	!-- followed by the new scaling parameters (th, q, etc.), and the new
249	!-- surface fluxes if required. The old routine ("circular") requires a
250	!-- different order of calls as the scaling parameters from the previous time
251	!-- steps are used to calculate the Obukhov length
252
253	!
254	!-- Depending on setting of most_method use the "old" routine
255	IF ( most_method == 'circular' ) THEN
256
257	CALL calc_scaling_parameters
258
259	CALL calc_uv_total
260
261	IF ( .NOT. neutral ) THEN
262	CALL calc_ol
263	ENDIF
264
265	CALL calc_us
266
267	CALL calc_surface_fluxes
268
269	!
270	!-- Use either Newton iteration or a lookup table for the bulk Richardson
271	!-- number to calculate the Obukhov length
272	ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' ) THEN
273
274	CALL calc_uv_total
275
276	IF ( .NOT. neutral ) THEN
277	CALL calc_ol
278	ENDIF
279
280	CALL calc_us
281
282	CALL calc_scaling_parameters
283
284	CALL calc_surface_fluxes
285
286	ENDIF
287
288	END SUBROUTINE surface_layer_fluxes
289
290
291	!------------------------------------------------------------------------------!
292	! Description:
293	! ------------
294	!> Initializing actions for the surface layer routine. Basically, this involves
295	!> the preparation of a lookup table for the the bulk Richardson number vs
296	!> Obukhov length L when using the lookup table method.
297	!------------------------------------------------------------------------------!
298	SUBROUTINE init_surface_layer_fluxes
299
300	IMPLICIT NONE
301
302	INTEGER(iwp) :: l, & !< Index for loop to create lookup table
303	num_steps_n !< Number of non-stretched zeta steps
304
305	LOGICAL :: terminate_run_l = .FALSE. !< Flag to terminate run (global)
306
307	REAL(wp), PARAMETER :: zeta_stretch = -10.0_wp !< Start of stretching in the free convection limit
308
309	REAL(wp), DIMENSION(:), ALLOCATABLE :: zeta_tmp
310
311
312	REAL(wp) :: zeta_step, & !< Increment of zeta
313	regr = 1.01_wp, & !< Stretching factor of zeta_step in the free convection limit
314	regr_old = 1.0E9_wp, & !< Stretching factor of last iteration step
315	z0h_min = 0.0_wp, & !< Minimum value of z0h to create table
316	z0_min = 0.0_wp !< Minimum value of z0 to create table
317	!
318	!-- When cloud physics is used, arrays for storing potential temperature and
319	!-- specific humidity at first grid level are required
320	IF ( cloud_physics ) THEN
321	ALLOCATE ( pt1(nysg:nyng,nxlg:nxrg) )
322	ALLOCATE ( qv1(nysg:nyng,nxlg:nxrg) )
323	ENDIF
324
325	!
326	!-- Allocate field for storing the horizontal velocity
327	ALLOCATE ( uv_total(nysg:nyng,nxlg:nxrg) )
328
329
330	!
331	!-- In case of runs with neutral statification, set Obukhov length to a
332	!-- large value
333	IF ( neutral ) ol = 1.0E10_wp
334
335	IF ( most_method == 'lookup' ) THEN
336
337	!
338	!-- Check for roughness heterogeneity. In that case terminate run and
339	!-- inform user
340	IF ( MINVAL( z0h ) /= MAXVAL( z0h ) .OR. &
341	MINVAL( z0 ) /= MAXVAL( z0 ) ) THEN
342	terminate_run_l = .TRUE.
343	ENDIF
344
345	#if defined( __parallel )
346	!
347	!-- Make a logical OR for all processes. Force termiation of model if result
348	!-- is TRUE
349	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
350	CALL MPI_ALLREDUCE( terminate_run_l, terminate_run, 1, MPI_LOGICAL, &
351	MPI_LOR, comm2d, ierr )
352	#else
353	terminate_run = terminate_run_l
354	#endif
355
356	IF ( terminate_run ) THEN
357	message_string = 'most_method = "lookup" cannot be used in ' // &
358	'combination with a prescribed roughness ' // &
359	'heterogeneity'
360	CALL message( 'surface_layer_fluxes', 'PA0417', 1, 2, 0, 6, 0 )
361	ENDIF
362
363	ALLOCATE( zeta_tmp(0:num_steps-1) )
364
365	!
366	!-- Use the lowest possible value for z_mo
367	k = MINVAL(nzb_s_inner)
368	z_mo = zu(k+1) - zw(k)
369
370	!
371	!-- Calculate z/L range from zeta_stretch to zeta_max using 90% of the
372	!-- available steps (num_steps). The calculation is done with negative
373	!-- values of zeta in order to simplify the stretching in the free
374	!-- convection limit for the remaining 10% of steps.
375	zeta_tmp(0) = - zeta_max
376	num_steps_n = ( num_steps * 9 / 10 ) - 1
377	zeta_step = (zeta_max - zeta_stretch) / REAL(num_steps_n)
378
379	DO l = 1, num_steps_n
380	zeta_tmp(l) = zeta_tmp(l-1) + zeta_step
381	ENDDO
382
383	!
384	!-- Calculate stretching factor for the free convection range
385	DO WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
386	regr_old = regr
387	regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr ) &
388	)**( 10.0_wp / REAL(num_steps) )
389	ENDDO
390
391	!
392	!-- Calculate z/L range from zeta_min to zeta_stretch
393	DO l = num_steps_n+1, num_steps-1
394	zeta_tmp(l) = zeta_tmp(l-1) + zeta_step
395	zeta_step = zeta_step * regr
396	ENDDO
397
398	!
399	!-- Save roughness lengths to temporary variables
400	z0h_min = z0h(nys,nxl)
401	z0_min = z0(nys,nxl)
402
403	!
404	!-- Calculate lookup table for the Richardson number versus Obukhov length
405	!-- The Richardson number (rib) is defined depending on the choice of
406	!-- boundary conditions for temperature
407	IF ( ibc_pt_b == 1 ) THEN
408	DO l = 0, num_steps-1
409	ol_tab(l) = - z_mo / zeta_tmp(num_steps-1-l)
410	rib_tab(l) = z_mo / ol_tab(l) / ( LOG( z_mo / z0_min ) &
411	- psi_m( z_mo / ol_tab(l) ) &
412	+ psi_m( z0_min / ol_tab(l) ) &
413	)**3
414	ENDDO
415	ELSE
416	DO l = 0, num_steps-1
417	ol_tab(l) = - z_mo / zeta_tmp(num_steps-1-l)
418	rib_tab(l) = z_mo / ol_tab(l) * ( LOG( z_mo / z0h_min ) &
419	- psi_h( z_mo / ol_tab(l) ) &
420	+ psi_h( z0h_min / ol_tab(l) ) &
421	) &
422	/ ( LOG( z_mo / z0_min ) &
423	- psi_m( z_mo / ol_tab(l) ) &
424	+ psi_m( z0_min / ol_tab(l) ) &
425	)**2
426	ENDDO
427	ENDIF
428
429	!
430	!-- Determine minimum values of rib in the lookup table. Set upper limit
431	!-- to critical Richardson number (0.25)
432	rib_min = MINVAL(rib_tab)
433	rib_max = 0.25 !MAXVAL(rib_tab)
434
435	DEALLOCATE( zeta_tmp )
436	ENDIF
437
438	END SUBROUTINE init_surface_layer_fluxes
439
440
441	!------------------------------------------------------------------------------!
442	! Description:
443	! ------------
444	!> Compute the absolute value of the horizontal velocity (relative to the
445	!> surface). This is required by all methods
446	!------------------------------------------------------------------------------!
447	SUBROUTINE calc_uv_total
448
449	IMPLICIT NONE
450
451
452	!$OMP PARALLEL DO PRIVATE( k )
453	!$acc kernels loop present( nzb_s_inner, u, uv_total, v ) private( j, k )
454	DO i = nxl, nxr
455	DO j = nys, nyn
456
457	k = nzb_s_inner(j,i)
458	uv_total(j,i) = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) &
459	- u(k,j,i) - u(k,j,i+1) ) )**2 + &
460	( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) &
461	- v(k,j,i) - v(k,j+1,i) ) )**2 )
462
463	!
464	!-- For too small values of the local wind, MOST does not work. A
465	!-- threshold value is thus set if required
466	! uv_total(j,i) = MAX(0.01_wp,uv_total(j,i))
467
468	ENDDO
469	ENDDO
470
471	!
472	!-- Values of uv_total need to be exchanged at the ghost boundaries
473	!$acc update host( uv_total )
474	CALL exchange_horiz_2d( uv_total )
475	!$acc update device( uv_total )
476
477	END SUBROUTINE calc_uv_total
478
479
480	!------------------------------------------------------------------------------!
481	! Description:
482	! ------------
483	!> Calculate the Obukhov length (L) and Richardson flux number (z/L)
484	!------------------------------------------------------------------------------!
485	SUBROUTINE calc_ol
486
487	IMPLICIT NONE
488
489	INTEGER(iwp) :: iter, & !< Newton iteration step
490	l !< look index
491
492	REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) :: rib !< Bulk Richardson number
493
494	REAL(wp) :: f, & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
495	f_d_ol, & !< Derivative of f
496	ol_l, & !< Lower bound of L for Newton iteration
497	ol_m, & !< Previous value of L for Newton iteration
498	ol_old, & !< Previous time step value of L
499	ol_u !< Upper bound of L for Newton iteration
500
501	IF ( TRIM( most_method ) /= 'circular' ) THEN
502
503	!$acc data present( nzb_s_inner, pt, q, qsws, rib, shf, uv_total, vpt, zu, zw )
504
505	!$OMP PARALLEL DO PRIVATE( k, z_mo )
506	!$acc kernels loop private( j, k, z_mo )
507	DO i = nxl, nxr
508	DO j = nys, nyn
509
510	k = nzb_s_inner(j,i)
511	z_mo = zu(k+1) - zw(k)
512
513	!
514	!-- Evaluate bulk Richardson number (calculation depends on
515	!-- definition based on setting of boundary conditions
516	IF ( ibc_pt_b /= 1 ) THEN
517	IF ( humidity ) THEN
518	rib(j,i) = g * z_mo * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
519	/ ( uv_total(j,i)*2 vpt(k+1,j,i) + 1.0E-20_wp )
520	ELSE
521	rib(j,i) = g * z_mo * ( pt(k+1,j,i) - pt(k,j,i) ) &
522	/ ( uv_total(j,i)*2 pt(k+1,j,i) + 1.0E-20_wp )
523	ENDIF
524	ELSE
525	!
526	!-- When using Neumann boundary conditions, the buoyancy flux
527	!-- is required but cannot be calculated at the surface, as pt
528	!-- and q are not known at the surface. Hence the values at
529	!-- first grid level are used to estimate the buoyancy flux
530	IF ( humidity ) THEN
531	rib(j,i) = - g * z_mo * ( ( 1.0_wp + 0.61_wp &
532	* q(k+1,j,i) ) * shf(j,i) + 0.61_wp &
533	* pt(k+1,j,i) * qsws(j,i) ) &
534	/ ( uv_total(j,i)*3 vpt(k+1,j,i) * kappa**2&
535	+ 1.0E-20_wp)
536	ELSE
537	rib(j,i) = - g * z_mo * shf(j,i) &
538	/ ( uv_total(j,i)*3 pt(k+1,j,i) * kappa**2 &
539	+ 1.0E-20_wp )
540	ENDIF
541	ENDIF
542
543	ENDDO
544	ENDDO
545	!$acc end data
546
547	ENDIF
548
549	!
550	!-- Calculate the Obukhov length either using a Newton iteration
551	!-- method, via a lookup table, or using the old circular way
552	IF ( TRIM( most_method ) == 'newton' ) THEN
553
554	!$OMP PARALLEL DO PRIVATE( k, z_mo )
555	!# WARNING: does not work on GPU so far because of DO-loop with
556	!# undetermined iterations
557	!!!!!!$acc kernels loop
558	DO i = nxl, nxr
559	DO j = nys, nyn
560
561	k = nzb_s_inner(j,i)
562	z_mo = zu(k+1) - zw(k)
563
564	!
565	!-- Store current value in case the Newton iteration fails
566	ol_old = ol(j,i)
567
568	!
569	!-- Ensure that the bulk Richardson number and the Obukhov
570	!-- lengtH have the same sign
571	IF ( rib(j,i) * ol(j,i) < 0.0_wp .OR. &
572	ABS( ol(j,i) ) == ol_max ) THEN
573	IF ( rib(j,i) > 0.0_wp ) ol(j,i) = 0.01_wp
574	IF ( rib(j,i) < 0.0_wp ) ol(j,i) = -0.01_wp
575	ENDIF
576	!
577	!-- Iteration to find Obukhov length
578	iter = 0
579	DO
580	iter = iter + 1
581	!
582	!-- In case of divergence, use the value of the previous time step
583	IF ( iter > 1000 ) THEN
584	ol(j,i) = ol_old
585	EXIT
586	ENDIF
587
588	ol_m = ol(j,i)
589	ol_l = ol_m - 0.001_wp * ol_m
590	ol_u = ol_m + 0.001_wp * ol_m
591
592
593	IF ( ibc_pt_b /= 1 ) THEN
594	!
595	!-- Calculate f = Ri - [...]/[...]^2 = 0
596	f = rib(j,i) - ( z_mo / ol_m ) * ( LOG( z_mo / z0h(j,i) )&
597	- psi_h( z_mo / ol_m ) &
598	+ psi_h( z0h(j,i) / ol_m ) &
599	) &
600	/ ( LOG( z_mo / z0(j,i) ) &
601	- psi_m( z_mo / ol_m ) &
602	+ psi_m( z0(j,i) / ol_m ) &
603	)**2
604
605	!
606	!-- Calculate df/dL
607	f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / z0h(j,i) ) &
608	- psi_h( z_mo / ol_u ) &
609	+ psi_h( z0h(j,i) / ol_u ) &
610	) &
611	/ ( LOG( z_mo / z0(j,i) ) &
612	- psi_m( z_mo / ol_u ) &
613	+ psi_m( z0(j,i) / ol_u ) &
614	)**2 &
615	+ ( z_mo / ol_l ) * ( LOG( z_mo / z0h(j,i) ) &
616	- psi_h( z_mo / ol_l ) &
617	+ psi_h( z0h(j,i) / ol_l ) &
618	) &
619	/ ( LOG( z_mo / z0(j,i) ) &
620	- psi_m( z_mo / ol_l ) &
621	+ psi_m( z0(j,i) / ol_l ) &
622	)**2 &
623	) / ( ol_u - ol_l )
624	ELSE
625	!
626	!-- Calculate f = Ri - 1 /[...]^3 = 0
627	f = rib(j,i) - ( z_mo / ol_m ) / ( LOG( z_mo / z0(j,i) )&
628	- psi_m( z_mo / ol_m ) &
629	+ psi_m( z0(j,i) / ol_m ) &
630	)**3
631
632	!
633	!-- Calculate df/dL
634	f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / z0(j,i) ) &
635	- psi_m( z_mo / ol_u ) &
636	+ psi_m( z0(j,i) / ol_u ) &
637	)**3 &
638	+ ( z_mo / ol_l ) / ( LOG( z_mo / z0(j,i) ) &
639	- psi_m( z_mo / ol_l ) &
640	+ psi_m( z0(j,i) / ol_l ) &
641	)**3 &
642	) / ( ol_u - ol_l )
643	ENDIF
644	!
645	!-- Calculate new L
646	ol(j,i) = ol_m - f / f_d_ol
647
648	!
649	!-- Ensure that the bulk Richardson number and the Obukhov
650	!-- length have the same sign and ensure convergence.
651	IF ( ol(j,i) * ol_m < 0.0_wp ) ol(j,i) = ol_m * 0.5_wp
652
653	!
654	!-- If unrealistic value occurs, set L to the maximum
655	!-- value that is allowed
656	IF ( ABS( ol(j,i) ) > ol_max ) THEN
657	ol(j,i) = ol_max
658	EXIT
659	ENDIF
660	!
661	!-- Check for convergence
662	IF ( ABS( ( ol(j,i) - ol_m ) / ol(j,i) ) < 1.0E-4_wp ) THEN
663	EXIT
664	ELSE
665	CYCLE
666	ENDIF
667
668	ENDDO
669
670	ENDDO
671	ENDDO
672
673	ELSEIF ( TRIM( most_method ) == 'lookup' ) THEN
674
675	!$OMP PARALLEL DO PRIVATE( k, z_mo )
676	!# WARNING: does not work on GPU so far because of DO WHILE construct
677	!!!!!!$acc kernels loop
678	DO i = nxl, nxr
679	DO j = nys, nyn
680
681	!
682	!-- If the bulk Richardson number is outside the range of the lookup
683	!-- table, set it to the exceeding threshold value
684	IF ( rib(j,i) < rib_min ) rib(j,i) = rib_min
685	IF ( rib(j,i) > rib_max ) rib(j,i) = rib_max
686
687	!
688	!-- Find the correct index bounds for linear interpolation. As the
689	!-- Richardson number will not differ very much from time step to
690	!-- time step , use the index from the last step and search in the
691	!-- correct direction
692	l = l_bnd
693	IF ( rib_tab(l) - rib(j,i) > 0.0_wp ) THEN
694	DO WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp .AND. l > 0 )
695	l = l-1
696	ENDDO
697	ELSE
698	DO WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp &
699	.AND. l < num_steps-1 )
700	l = l+1
701	ENDDO
702	ENDIF
703	l_bnd = l
704
705	!
706	!-- Linear interpolation to find the correct value of z/L
707	ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) ) &
708	/ ( rib_tab(l) - rib_tab(l-1) ) &
709	* ( rib(j,i) - rib_tab(l-1) ) )
710
711	ENDDO
712	ENDDO
713
714	ELSEIF ( TRIM( most_method ) == 'circular' ) THEN
715
716	!$OMP PARALLEL DO PRIVATE( k, z_mo )
717	!$acc kernels loop present( nzb_s_inner, ol, pt, pt1, q, ql, qs, qv1, ts, us, vpt, zu, zw ) private( j, k, z_mo )
718	DO i = nxl, nxr
719	DO j = nys, nyn
720
721	k = nzb_s_inner(j,i)
722	z_mo = zu(k+1) - zw(k)
723
724	IF ( .NOT. humidity ) THEN
725	ol(j,i) = ( pt(k+1,j,i) * us(j,i)*2 ) / ( kappa g &
726	* ts(j,i) + 1E-30_wp )
727	ELSEIF ( cloud_physics ) THEN
728
729	ol(j,i) = ( vpt(k+1,j,i) * us(j,i)*2 ) / ( kappa g &
730	* ( ts(j,i) + 0.61_wp * pt1(j,i) * qs(j,i) &
731	+ 0.61_wp * qv1(j,i) * ts(j,i) - ts(j,i) &
732	* ql(k+1,j,i) ) + 1E-30_wp )
733	ELSE
734	ol(j,i) = ( vpt(k+1,j,i) * us(j,i)*2 ) / ( kappa g &
735	* ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i) &
736	+ 0.61_wp * q(k+1,j,i) * ts(j,i) ) + 1E-30_wp )
737	ENDIF
738	!
739	!-- Limit the value range of the Obukhov length.
740	!-- This is necessary for very small velocities (u,v --> 0), because
741	!-- the absolute value of ol can then become very small, which in
742	!-- consequence would result in very large shear stresses and very
743	!-- small momentum fluxes (both are generally unrealistic).
744	IF ( ( z_mo / ( ol(j,i) + 1E-30_wp ) ) < zeta_min ) &
745	ol(j,i) = z_mo / zeta_min
746	IF ( ( z_mo / ( ol(j,i) + 1E-30_wp ) ) > zeta_max ) &
747	ol(j,i) = z_mo / zeta_max
748	ENDDO
749	ENDDO
750
751	ENDIF
752
753	!
754	!-- Values of ol at ghost point locations are needed for the evaluation
755	!-- of usws and vsws.
756	!$acc update host( ol )
757	CALL exchange_horiz_2d( ol )
758	!$acc update device( ol )
759
760	END SUBROUTINE calc_ol
761
762	!
763	!-- Calculate friction velocity u*
764	SUBROUTINE calc_us
765
766	IMPLICIT NONE
767
768	!$OMP PARALLEL DO PRIVATE( k, z_mo )
769	!$acc kernels loop present( nzb_s_inner, ol, us, uv_total, zu, zw, z0 ) private( j, k, z_mo )
770	DO i = nxlg, nxrg
771	DO j = nysg, nyng
772
773	k = nzb_s_inner(j,i)+1
774	z_mo = zu(k+1) - zw(k)
775
776	!
777	!-- Compute u* at the scalars' grid points
778	us(j,i) = kappa * uv_total(j,i) / ( LOG( z_mo / z0(j,i) ) &
779	- psi_m( z_mo / ol(j,i) ) &
780	+ psi_m( z0(j,i) / ol(j,i) ) )
781	ENDDO
782	ENDDO
783
784	END SUBROUTINE calc_us
785
786	!
787	!-- Calculate potential temperature and specific humidity at first grid level
788	SUBROUTINE calc_pt_q
789
790	IMPLICIT NONE
791
792	!$acc kernels loop present( nzb_s_inner, pt, pt1, pt_d_t, q, ql, qv1 ) private( j, k )
793	DO i = nxlg, nxrg
794	DO j = nysg, nyng
795	k = nzb_s_inner(j,i)+1
796	pt1(j,i) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
797	qv1(j,i) = q(k,j,i) - ql(k,j,i)
798	ENDDO
799	ENDDO
800
801	END SUBROUTINE calc_pt_q
802
803	!
804	!-- Calculate the other MOST scaling parameters theta, q, (qr, nr)
805	SUBROUTINE calc_scaling_parameters
806
807	IMPLICIT NONE
808
809	!
810	!-- Data information for accelerators
811	!$acc data present( e, nrsws, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt ) &
812	!$acc present( q, qs, qsws, qrsws, shf, ts, u, us, usws, v ) &
813	!$acc present( vpt, vsws, zu, zw, z0, z0h )
814	!
815	!-- Compute theta*
816	IF ( constant_heatflux ) THEN
817
818	!
819	!-- For a given heat flux in the surface layer:
820	!$OMP PARALLEL DO
821	!$acc kernels loop private( j )
822	DO i = nxlg, nxrg
823	DO j = nysg, nyng
824	ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30_wp )
825	!
826	!-- ts must be limited, because otherwise overflow may occur in case
827	!-- of us=0 when computing ol further below
828	IF ( ts(j,i) < -1.05E5_wp ) ts(j,i) = -1.0E5_wp
829	IF ( ts(j,i) > 1.0E5_wp ) ts(j,i) = 1.0E5_wp
830	ENDDO
831	ENDDO
832
833	ELSE
834	!
835	!-- For a given surface temperature:
836	IF ( large_scale_forcing .AND. lsf_surf ) THEN
837	!$OMP PARALLEL DO
838	!$acc kernels loop private( j, k )
839	DO i = nxlg, nxrg
840	DO j = nysg, nyng
841	k = nzb_s_inner(j,i)
842	pt(k,j,i) = pt_surface
843	ENDDO
844	ENDDO
845	ENDIF
846
847	!$OMP PARALLEL DO PRIVATE( k, z_mo )
848	!$acc kernels loop present( nzb_s_inner, ol, pt, pt1, ts, zu, zw, z0h ) private( j, k, z_mo )
849	DO i = nxlg, nxrg
850	DO j = nysg, nyng
851
852	k = nzb_s_inner(j,i)
853	z_mo = zu(k+1) - zw(k)
854
855	IF ( cloud_physics ) THEN
856	ts(j,i) = kappa * ( pt1(j,i) - pt(k,j,i) ) &
857	/ ( LOG( z_mo / z0h(j,i) ) &
858	- psi_h( z_mo / ol(j,i) ) &
859	+ psi_h( z0h(j,i) / ol(j,i) ) )
860	ELSE
861	ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) &
862	/ ( LOG( z_mo / z0h(j,i) ) &
863	- psi_h( z_mo / ol(j,i) ) &
864	+ psi_h( z0h(j,i) / ol(j,i) ) )
865	ENDIF
866
867	ENDDO
868	ENDDO
869	ENDIF
870
871	!
872	!-- If required compute q*
873	IF ( humidity ) THEN
874	IF ( constant_waterflux ) THEN
875	!
876	!-- For a given water flux in the surface layer
877	!$OMP PARALLEL DO
878	!$acc kernels loop private( j )
879	DO i = nxlg, nxrg
880	DO j = nysg, nyng
881	qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp )
882	ENDDO
883	ENDDO
884
885	ELSE
886	coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. &
887	run_coupled )
888
889	IF ( large_scale_forcing .AND. lsf_surf ) THEN
890	!$OMP PARALLEL DO
891	!$acc kernels loop private( j, k )
892	DO i = nxlg, nxrg
893	DO j = nysg, nyng
894	k = nzb_s_inner(j,i)
895	q(k,j,i) = q_surface
896	ENDDO
897	ENDDO
898	ENDIF
899
900	!$OMP PARALLEL DO PRIVATE( e_s, k, z_mo )
901	!$acc kernels loop independent present( nzb_s_inner, ol, pt, q, qs, qv1, zu, zw, z0q ) private( e_s, j, k, z_mo )
902	DO i = nxlg, nxrg
903	!$acc loop independent
904	DO j = nysg, nyng
905
906	k = nzb_s_inner(j,i)
907	z_mo = zu(k+1) - zw(k)
908
909	!
910	!-- Assume saturation for atmosphere coupled to ocean (but not
911	!-- in case of precursor runs)
912	IF ( coupled_run ) THEN
913	e_s = 6.1_wp * &
914	EXP( 0.07_wp * ( MIN(pt(k,j,i),pt(k+1,j,i)) &
915	- 273.15_wp ) )
916	q(k,j,i) = 0.622_wp * e_s / ( surface_pressure - e_s )
917	ENDIF
918
919	IF ( cloud_physics ) THEN
920	qs(j,i) = kappa * ( qv1(j,i) - q(k,j,i) ) &
921	/ ( LOG( z_mo / z0q(j,i) ) &
922	- psi_h( z_mo / ol(j,i) ) &
923	+ psi_h( z0q(j,i) / ol(j,i) ) )
924
925	ELSE
926	qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) &
927	/ ( LOG( z_mo / z0q(j,i) ) &
928	- psi_h( z_mo / ol(j,i) ) &
929	+ psi_h( z0q(j,i) / ol(j,i) ) )
930	ENDIF
931
932	ENDDO
933	ENDDO
934	ENDIF
935	ENDIF
936
937	!
938	!-- If required compute s*
939	IF ( passive_scalar ) THEN
940	IF ( constant_scalarflux ) THEN
941	!
942	!-- For a given water flux in the surface layer
943	!$OMP PARALLEL DO
944	!$acc kernels loop private( j )
945	DO i = nxlg, nxrg
946	DO j = nysg, nyng
947	ss(j,i) = -ssws(j,i) / ( us(j,i) + 1E-30_wp )
948	ENDDO
949	ENDDO
950	ENDIF
951	ENDIF
952
953
954	!
955	!-- If required compute qr* and nr*
956	IF ( cloud_physics .AND. microphysics_seifert ) &
957	THEN
958
959	!$OMP PARALLEL DO PRIVATE( k, z_mo )
960	!$acc kernels loop independent present( nr, nrs, nzb_s_inner, ol, qr, qrs, zu, zw, z0q ) private( j, k, z_mo )
961	DO i = nxlg, nxrg
962	!$acc loop independent
963	DO j = nysg, nyng
964
965	k = nzb_s_inner(j,i)
966	z_mo = zu(k+1) - zw(k)
967
968	qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) &
969	/ ( LOG( z_mo / z0q(j,i) ) &
970	- psi_h( z_mo / ol(j,i) ) &
971	+ psi_h( z0q(j,i) / ol(j,i) ) )
972
973	nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) &
974	/ ( LOG( z_mo / z0q(j,i) ) &
975	- psi_h( z_mo / ol(j,i) ) &
976	+ psi_h( z0q(j,i) / ol(j,i) ) )
977	ENDDO
978	ENDDO
979
980	ENDIF
981	!$acc end data
982
983	END SUBROUTINE calc_scaling_parameters
984
985
986
987	!
988	!-- Calculate surface fluxes usws, vsws, shf, qsws, (qrsws, nrsws)
989	SUBROUTINE calc_surface_fluxes
990
991	IMPLICIT NONE
992
993	REAL(wp) :: ol_mid !< Grid-interpolated L
994
995	!
996	!-- Compute u'w' for the total model domain.
997	!-- First compute the corresponding component of u* and square it.
998	!$OMP PARALLEL DO PRIVATE( k, ol_mid, z_mo )
999	!$acc kernels loop present( nzb_u_inner, ol, u, us, usws, zu, zw, z0 ) private( j, k, z_mo )
1000	DO i = nxl, nxr
1001	DO j = nys, nyn
1002
1003	k = nzb_u_inner(j,i)
1004	z_mo = zu(k+1) - zw(k)
1005	!
1006	!-- Compute bulk Obukhov length for this point
1007	ol_mid = 0.5_wp * ( ol(j,i-1) + ol(j,i) )
1008
1009	IF ( ol_mid == 0.0_wp ) THEN
1010	ol_mid = MIN(ol(j,i-1), ol(j,i))
1011	ENDIF
1012
1013	usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) ) &
1014	/ ( LOG( z_mo / z0(j,i) ) &
1015	- psi_m( z_mo / ol_mid ) &
1016	+ psi_m( z0(j,i) / ol_mid ) )
1017
1018	usws(j,i) = -usws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
1019	ENDDO
1020	ENDDO
1021
1022	!
1023	!-- Compute v'w' for the total model domain.
1024	!-- First compute the corresponding component of u* and square it.
1025	!$OMP PARALLEL DO PRIVATE( k, ol_mid, z_mo )
1026	!$acc kernels loop present( nzb_v_inner, ol, v, us, vsws, zu, zw, z0 ) private( j, k, ol_mid, z_mo )
1027	DO i = nxl, nxr
1028	DO j = nys, nyn
1029
1030	k = nzb_v_inner(j,i)
1031	z_mo = zu(k+1) - zw(k)
1032	!
1033	!-- Compute bulk Obukhov length for this point
1034	ol_mid = 0.5_wp * ( ol(j-1,i) + ol(j,i) )
1035
1036	IF ( ol_mid == 0.0_wp ) THEN
1037	ol_mid = MIN(ol(j-1,i), ol(j-1,i))
1038	ENDIF
1039
1040	vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) &
1041	/ ( LOG( z_mo / z0(j,i) ) &
1042	- psi_m( z_mo / ol_mid ) &
1043	+ psi_m( z0(j,i) / ol_mid ) )
1044
1045	vsws(j,i) = -vsws(j,i) * 0.5_wp * ( us(j,i-1) + us(j,i) )
1046
1047	ENDDO
1048	ENDDO
1049
1050	!
1051	!-- Exchange the boundaries for the momentum fluxes (is this still required?)
1052	!$acc update host( usws, vsws )
1053	CALL exchange_horiz_2d( usws )
1054	CALL exchange_horiz_2d( vsws )
1055	!$acc update device( usws, vsws )
1056
1057	!
1058	!-- Compute the vertical kinematic heat flux
1059	IF ( .NOT. constant_heatflux .AND. ( simulated_time <= &
1060	skip_time_do_lsm .OR. .NOT. land_surface ) ) THEN
1061	!$OMP PARALLEL DO
1062	!$acc kernels loop independent present( shf, ts, us )
1063	DO i = nxlg, nxrg
1064	!$acc loop independent
1065	DO j = nysg, nyng
1066	shf(j,i) = -ts(j,i) * us(j,i)
1067	ENDDO
1068	ENDDO
1069
1070	ENDIF
1071
1072	!
1073	!-- Compute the vertical water flux
1074	IF ( .NOT. constant_waterflux .AND. humidity .AND. &
1075	( simulated_time <= skip_time_do_lsm &
1076	.OR. .NOT. land_surface ) ) THEN
1077	!$OMP PARALLEL DO
1078	!$acc kernels loop independent present( qs, qsws, us )
1079	DO i = nxlg, nxrg
1080	!$acc loop independent
1081	DO j = nysg, nyng
1082	qsws(j,i) = -qs(j,i) * us(j,i)
1083	ENDDO
1084	ENDDO
1085
1086	ENDIF
1087
1088	!
1089	!-- Compute the vertical scalar flux
1090	IF ( .NOT. constant_scalarflux .AND. passive_scalar .AND. &
1091	( simulated_time <= skip_time_do_lsm &
1092	.OR. .NOT. land_surface ) ) THEN
1093	!$OMP PARALLEL DO
1094	!$acc kernels loop independent present( qs, qsws, us )
1095	DO i = nxlg, nxrg
1096	!$acc loop independent
1097	DO j = nysg, nyng
1098	ssws(j,i) = -ss(j,i) * us(j,i)
1099	ENDDO
1100	ENDDO
1101
1102	ENDIF
1103
1104	!
1105	!-- Compute (turbulent) fluxes of rain water content and rain drop conc.
1106	IF ( cloud_physics .AND. microphysics_seifert ) THEN
1107	!$OMP PARALLEL DO
1108	!$acc kernels loop independent present( nrs, nrsws, qrs, qrsws, us )
1109	DO i = nxlg, nxrg
1110	!$acc loop independent
1111	DO j = nysg, nyng
1112	qrsws(j,i) = -qrs(j,i) * us(j,i)
1113	nrsws(j,i) = -nrs(j,i) * us(j,i)
1114	ENDDO
1115	ENDDO
1116	ENDIF
1117
1118	!
1119	!-- Bottom boundary condition for the TKE
1120	IF ( ibc_e_b == 2 ) THEN
1121	!$OMP PARALLEL DO
1122	!$acc kernels loop independent present( e, nzb_s_inner, us )
1123	DO i = nxlg, nxrg
1124	!$acc loop independent
1125	DO j = nysg, nyng
1126	k = nzb_s_inner(j,i)
1127	e(k+1,j,i) = ( us(j,i) / 0.1_wp )**2
1128	!
1129	!-- As a test: cm = 0.4
1130	! e(k+1,j,i) = ( us(j,i) / 0.4_wp )**2
1131	e(k,j,i) = e(k+1,j,i)
1132	ENDDO
1133	ENDDO
1134	ENDIF
1135
1136	END SUBROUTINE calc_surface_fluxes
1137
1138
1139	!
1140	!-- Integrated stability function for momentum
1141	FUNCTION psi_m( zeta )
1142
1143	USE kinds
1144
1145	IMPLICIT NONE
1146
1147	REAL(wp) :: psi_m !< Integrated similarity function result
1148	REAL(wp) :: zeta !< Stability parameter z/L
1149	REAL(wp) :: x !< dummy variable
1150
1151	REAL(wp), PARAMETER :: a = 1.0_wp !< constant
1152	REAL(wp), PARAMETER :: b = 0.66666666666_wp !< constant
1153	REAL(wp), PARAMETER :: c = 5.0_wp !< constant
1154	REAL(wp), PARAMETER :: d = 0.35_wp !< constant
1155	REAL(wp), PARAMETER :: c_d_d = c / d !< constant
1156	REAL(wp), PARAMETER :: bc_d_d = b * c / d !< constant
1157
1158
1159	IF ( zeta < 0.0_wp ) THEN
1160	x = SQRT( SQRT( 1.0_wp - 16.0_wp * zeta ) )
1161	psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2 &
1162	* ( 1.0_wp + x*2 ) 0.125_wp )
1163	ELSE
1164
1165	psi_m = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - a * zeta &
1166	- bc_d_d
1167	!
1168	!-- Old version for stable conditions (only valid for z/L < 0.5)
1169	!-- psi_m = - 5.0_wp * zeta
1170
1171	ENDIF
1172
1173	END FUNCTION psi_m
1174
1175
1176	!
1177	!-- Integrated stability function for heat and moisture
1178	FUNCTION psi_h( zeta )
1179
1180	USE kinds
1181
1182	IMPLICIT NONE
1183
1184	REAL(wp) :: psi_h !< Integrated similarity function result
1185	REAL(wp) :: zeta !< Stability parameter z/L
1186	REAL(wp) :: x !< dummy variable
1187
1188	REAL(wp), PARAMETER :: a = 1.0_wp !< constant
1189	REAL(wp), PARAMETER :: b = 0.66666666666_wp !< constant
1190	REAL(wp), PARAMETER :: c = 5.0_wp !< constant
1191	REAL(wp), PARAMETER :: d = 0.35_wp !< constant
1192	REAL(wp), PARAMETER :: c_d_d = c / d !< constant
1193	REAL(wp), PARAMETER :: bc_d_d = b * c / d !< constant
1194
1195
1196	IF ( zeta < 0.0_wp ) THEN
1197	x = SQRT( 1.0_wp - 16.0_wp * zeta )
1198	psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
1199	ELSE
1200	psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp &
1201	+ 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d &
1202	+ 1.0_wp
1203	!
1204	!-- Old version for stable conditions (only valid for z/L < 0.5)
1205	!-- psi_h = - 5.0_wp * zeta
1206	ENDIF
1207
1208	END FUNCTION psi_h
1209
1210	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 |