Home

Context Navigation

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

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

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