Home

Context Navigation

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

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

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