Home

Context Navigation

source: palm/trunk/SOURCE/biometeorology_pet_mod.f90 @ 3507

Last change on this file since 3507 was 3448, checked in by kanani, 6 years ago
Implementation of human thermal indices (from branch biomet_p2 at r3444)
Property svn:executable set to ``*
File size: 19.7 KB

Rev	Line
[3448]	1	!> @file biometeorology_PET.f90
	2	!------------------------------------------------------------------------------!
	3	! This file is part of PALM-4U.
	4	!
	5	! PALM-4U 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-4U 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 2019, Deutscher Wetterdienst (DWD) /
	18	! German Meteorological Service (DWD)
	19	!------------------------------------------------------------------------------!
	20	!
	21	! Current revisions: 003
	22	! -----------------
	23	!
	24	!
	25	! Former revisions: 001
	26	! -----------------
	27	! $Id$
	28	! Initial revision 001
	29	!
	30	!
	31	!
	32	! Authors:
	33	! --------
	34	! @author Peter HÃ¶ppe (original author)
	35	! @author Dominik Froehlich <dominik.froehlich@dwd.de> (PALM modification)
	36	!
	37	!------------------------------------------------------------------------------!
	38	MODULE biometeorology_pet_mod
	39
	40	!-- Load required variables from existing modules
	41	USE kinds !< to set precision of INTEGER and REAL according to PALM
	42
	43	IMPLICIT NONE
	44
	45	PRIVATE
	46
	47	!-- Internal constants:
	48	REAL(wp), PARAMETER :: po = 1013.25_wp !< Air pressure at sea level (hPa)
	49	REAL(wp), PARAMETER :: rob = 1.06_wp !<
	50	REAL(wp), PARAMETER :: cb = 3640._wp !<
	51	REAL(wp), PARAMETER :: food = 0._wp !< Heat gain by food (W)
	52	REAL(wp), PARAMETER :: emsk = 0.99_wp !< Longwave emission coef. of skin
	53	REAL(wp), PARAMETER :: emcl = 0.95_wp !< Longwave emission coef. of cloth
	54	REAL(wp), PARAMETER :: evap = 2.42_wp * 10._wp **6._wp !<
	55	REAL(wp), PARAMETER :: sigm = 5.67_wp * 10._wp **(-8._wp) !< Stefan-Boltzmann-Const.
	56	REAL(wp), PARAMETER :: cair = 1010._wp !<
	57	REAL(wp), PARAMETER :: cels_offs = 273.15_wp !< Kelvin Celsius Offset (K)
	58
	59	!-- Internal variables:
	60	REAL(wp) :: h !< Internal heat (W)
	61
	62	!-- MEMI configuration
	63	REAL(wp) :: age !< Persons age (a)
	64	REAL(wp) :: mbody !< Persons body mass (kg)
	65	REAL(wp) :: ht !< Persons height (m)
	66	REAL(wp) :: work !< Work load (W)
	67	REAL(wp) :: eta !< Work efficiency (dimensionless)
	68	REAL(wp) :: icl !< Clothing insulation index (clo)
	69	REAL(wp) :: fcl !< Surface area modification by clothing (factor)
	70	INTEGER(iwp) :: pos !< Posture: 1 = standing, 2 = sitting
	71	INTEGER(iwp) :: sex !< Sex: 1 = male, 2 = female
	72
	73	PUBLIC calculate_pet_static
	74
	75	!-- PALM interfaces:
	76	INTERFACE calculate_pet_static
	77	MODULE PROCEDURE calculate_pet_static
	78	END INTERFACE calculate_pet_static
	79
	80	CONTAINS
	81
	82
	83	!------------------------------------------------------------------------------!
	84	!
	85	! Description:
	86	! ------------
	87	!> Physiologically Equivalent Temperature (PET),
	88	! stationary (calculated based on MEMI),
	89	! Subroutine based on PETBER vers. 1.5.1996 by P. Hoeppe
	90	!------------------------------------------------------------------------------!
	91
	92	SUBROUTINE calculate_pet_static( ta, vpa, v, tmrt, p, tx )
	93
	94	IMPLICIT NONE
	95
	96	!-- Input arguments:
	97	REAL(wp), INTENT( IN ) :: ta !< Air temperature (Â°C)
	98	REAL(wp), INTENT( IN ) :: tmrt !< Mean radiant temperature (Â°C)
	99	REAL(wp), INTENT( IN ) :: v !< Wind speed (m/s)
	100	REAL(wp), INTENT( IN ) :: vpa !< Vapor pressure (hPa)
	101	REAL(wp), INTENT( IN ) :: p !< Air pressure (hPa)
	102
	103	!-- Output arguments:
	104	REAL(wp), INTENT ( OUT ) :: tx !< PET (Â°C)
	105
	106	!-- Internal variables:
	107	REAL(wp) :: acl
	108	REAL(wp) :: adu
	109	REAL(wp) :: aeff
	110	REAL(wp) :: ere
	111	REAL(wp) :: erel
	112	REAL(wp) :: esw
	113	REAL(wp) :: facl
	114	REAL(wp) :: feff
	115	REAL(wp) :: rdcl
	116	REAL(wp) :: rdsk
	117	REAL(wp) :: rtv
	118	REAL(wp) :: vpts
	119	REAL(wp) :: tsk
	120	REAL(wp) :: tcl
	121	REAL(wp) :: wetsk
	122
	123	!-- Configuration
	124	age = 35._wp
	125	mbody = 75._wp
	126	ht = 1.75_wp
	127	work = 80._wp
	128	eta = 0._wp
	129	icl = 0.9_wp
	130	fcl = 1.15_wp
	131
	132	!-- Call subfunctions
	133	CALL in_body ( ere, erel, p, rtv, ta, vpa )
	134
	135	CALL heat_exch ( acl, adu, aeff, ere, erel, esw, facl, feff, &
	136	p, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa, vpts, wetsk )
	137
	138	CALL pet_iteration ( acl, adu, aeff, esw, facl, feff, p, rdcl, &
	139	rdsk, rtv, ta, tcl, tsk, tx, vpts, wetsk )
	140
	141
	142	END SUBROUTINE calculate_pet_static
	143
	144
	145	!------------------------------------------------------------------------------!
	146	! Description:
	147	! ------------
	148	!> Calculate internal energy ballance
	149	!------------------------------------------------------------------------------!
	150	SUBROUTINE in_body ( ere, erel, p, rtv, ta, vpa )
	151
	152	!-- Input arguments:
	153	REAL(wp), INTENT( IN ) :: p !< air pressure (hPa)
	154	REAL(wp), INTENT( IN ) :: ta !< air temperature (Â°C)
	155	REAL(wp), INTENT( IN ) :: vpa !< vapor pressure (hPa)
	156
	157	!-- Output arguments:
	158	REAL(wp), INTENT( OUT ) :: ere !< energy ballance (W)
	159	REAL(wp), INTENT( OUT ) :: erel !< latent energy ballance (W)
	160	REAL(wp), INTENT( OUT ) :: rtv !<
	161
	162	!-- Internal variables:
	163	REAL(wp) :: eres
	164	REAL(wp) :: met
	165	REAL(wp) :: tex
	166	REAL(wp) :: vpex
	167
	168	met = 3.45_wp * mbody ** ( 3._wp / 4._wp ) * (1._wp + 0.004_wp * &
	169	( 30._wp - age) + 0.010_wp * ( ( ht * 100._wp / &
	170	( mbody ** ( 1._wp / 3._wp ) ) ) - 43.4_wp ) )
	171
	172	met = work + met
	173
	174	h = met * (1._wp - eta)
	175
	176	!-- Sensible respiration energy
	177	tex = 0.47_wp * ta + 21.0_wp
	178	rtv = 1.44_wp * 10._wp ** (-6._wp) * met
	179	eres = cair * (ta - tex) * rtv
	180
	181	!-- Latent respiration energy
	182	vpex = 6.11_wp * 10._wp ** ( 7.45_wp * tex / ( 235._wp + tex ) )
	183	erel = 0.623_wp * evap / p * ( vpa - vpex ) * rtv
	184
	185	!-- Sum of the results
	186	ere = eres + erel
	187
	188	END SUBROUTINE in_body
	189
	190
	191	!------------------------------------------------------------------------------!
	192	! Description:
	193	! ------------
	194	!> Calculate heat gain or loss
	195	!------------------------------------------------------------------------------!
	196	SUBROUTINE heat_exch ( acl, adu, aeff, ere, erel, esw, facl, feff, &
	197	p, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa, vpts, wetsk )
	198
	199
	200	!-- Input arguments:
	201	REAL(wp), INTENT( IN ) :: ere !< Energy ballance (W)
	202	REAL(wp), INTENT( IN ) :: erel !< Latent energy ballance (W)
	203	REAL(wp), INTENT( IN ) :: p !< Air pressure (hPa)
	204	REAL(wp), INTENT( IN ) :: ta !< Air temperature (Â°C)
	205	REAL(wp), INTENT( IN ) :: tmrt !< Mean radiant temperature (Â°C)
	206	REAL(wp), INTENT( IN ) :: v !< Wind speed (m/s)
	207	REAL(wp), INTENT( IN ) :: vpa !< Vapor pressure (hPa)
	208
	209	!-- Output arguments:
	210	REAL(wp), INTENT( OUT ) :: acl !< Clothing surface area (mÂ²)
	211	REAL(wp), INTENT( OUT ) :: adu !< Du-Bois area (mÂ²)
	212	REAL(wp), INTENT( OUT ) :: aeff !< Effective surface area (mÂ²)
	213	REAL(wp), INTENT( OUT ) :: esw !< Energy-loss through sweat evap. (W)
	214	REAL(wp), INTENT( OUT ) :: facl !< Surface area extension through clothing (factor)
	215	REAL(wp), INTENT( OUT ) :: feff !< Surface modification by posture (factor)
	216	REAL(wp), INTENT( OUT ) :: rdcl !< Diffusion resistence of clothing (factor)
	217	REAL(wp), INTENT( OUT ) :: rdsk !< Diffusion resistence of skin (factor)
	218	REAL(wp), INTENT( OUT ) :: tcl !< Clothing temperature (Â°C)
	219	REAL(wp), INTENT( OUT ) :: tsk !< Skin temperature (Â°C)
	220	REAL(wp), INTENT( OUT ) :: vpts !< Sat. vapor pressure over skin (hPa)
	221	REAL(wp), INTENT( OUT ) :: wetsk !< Fraction of wet skin (dimensionless)
	222
	223	!-- Internal variables
	224	REAL(wp) :: c(0:10)
	225	REAL(wp) :: cbare
	226	REAL(wp) :: cclo
	227	REAL(wp) :: csum
	228	REAL(wp) :: di
	229	REAL(wp) :: ed
	230	REAL(wp) :: enbal
	231	REAL(wp) :: enbal2
	232	REAL(wp) :: eswdif
	233	REAL(wp) :: eswphy
	234	REAL(wp) :: eswpot
	235	REAL(wp) :: fec
	236	REAL(wp) :: hc
	237	REAL(wp) :: he
	238	REAL(wp) :: htcl
	239	REAL(wp) :: r1
	240	REAL(wp) :: r2
	241	REAL(wp) :: rbare
	242	REAL(wp) :: rcl
	243	REAL(wp) :: rclo
	244	REAL(wp) :: rclo2
	245	REAL(wp) :: rsum
	246	REAL(wp) :: sw
	247	REAL(wp) :: swf
	248	REAL(wp) :: swm
	249	REAL(wp) :: tbody
	250	REAL(wp) :: tcore(1:7)
	251	REAL(wp) :: vb
	252	REAL(wp) :: vb1
	253	REAL(wp) :: vb2
	254	REAL(wp) :: wd
	255	REAL(wp) :: wr
	256	REAL(wp) :: ws
	257	REAL(wp) :: wsum
	258	REAL(wp) :: xx
	259	REAL(wp) :: y
	260	INTEGER(iwp) :: count1
	261	INTEGER(iwp) :: count3
	262	INTEGER(iwp) :: j
	263	INTEGER(iwp) :: i
	264	LOGICAL :: skipIncreaseCount
	265
	266	wetsk = 0._wp
	267	adu = 0.203_wp * mbody ** 0.425_wp * ht ** 0.725_wp
	268
	269	hc = 2.67_wp + ( 6.5_wp * v ** 0.67_wp)
	270	hc = hc * (p /po) ** 0.55_wp
	271	feff = 0.725_wp !< Posture: 0.725 for stading
	272	facl = (- 2.36_wp + 173.51_wp * icl - 100.76_wp * icl * icl + 19.28_wp &
	273	* (icl ** 3._wp)) / 100._wp
	274
	275	IF ( facl > 1._wp ) facl = 1._wp
	276	rcl = ( icl / 6.45_wp ) / facl
	277	IF ( icl >= 2._wp ) y = 1._wp
	278
	279	IF ( ( icl > 0.6_wp ) .AND. ( icl < 2._wp ) ) y = ( ht - 0.2_wp ) / ht
	280	IF ( ( icl <= 0.6_wp ) .AND. ( icl > 0.3_wp ) ) y = 0.5_wp
	281	IF ( ( icl <= 0.3_wp ) .AND. ( icl > 0._wp ) ) y = 0.1_wp
	282
	283	r2 = adu * (fcl - 1._wp + facl) / (2._wp * 3.14_wp * ht * y)
	284	r1 = facl * adu / (2._wp * 3.14_wp * ht * y)
	285
	286	di = r2 - r1
	287
	288	!-- Skin temperatur
	289	DO j = 1, 7
	290
	291	tsk = 34._wp
	292	count1 = 0_iwp
	293	tcl = ( ta + tmrt + tsk ) / 3._wp
	294	count3 = 1_iwp
	295	enbal2 = 0._wp
	296
	297	DO i = 1, 100 ! allow for 100 iterations max
	298	acl = adu * facl + adu * ( fcl - 1._wp )
	299	rclo2 = emcl * sigm * ( ( tcl + cels_offs )**4._wp - &
	300	( tmrt + cels_offs )** 4._wp ) * feff
	301	htcl = 6.28_wp * ht * y * di / ( rcl * LOG( r2 / r1 ) * acl )
	302	tsk = 1._wp / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl
	303
	304	!-- Radiation saldo
	305	aeff = adu * feff
	306	rbare = aeff * ( 1._wp - facl ) * emsk * sigm * &
	307	( ( tmrt + cels_offs ) 4._wp - ( tsk + cels_offs ) 4._wp )
	308	rclo = feff * acl * emcl * sigm * &
	309	( ( tmrt + cels_offs ) 4._wp - ( tcl + cels_offs ) 4._wp )
	310	rsum = rbare + rclo
	311
	312	!-- Convection
	313	cbare = hc * ( ta - tsk ) * adu * ( 1._wp - facl )
	314	cclo = hc * ( ta - tcl ) * acl
	315	csum = cbare + cclo
	316
	317	!-- Core temperature
	318	c(0) = h + ere
	319	c(1) = adu * rob * cb
	320	c(2) = 18._wp - 0.5_wp * tsk
	321	c(3) = 5.28_wp * adu * c(2)
	322	c(4) = 0.0208_wp * c(1)
	323	c(5) = 0.76075_wp * c(1)
	324	c(6) = c(3) - c(5) - tsk * c(4)
	325	c(7) = - c(0) * c(2) - tsk * c(3) + tsk * c(5)
	326	c(8) = c(6) * c(6) - 4._wp * c(4) * c(7)
	327	c(9) = 5.28_wp * adu - c(5) - c(4) * tsk
	328	c(10) = c(9) * c(9) - 4._wp * c(4) * &
	329	( c(5) * tsk - c(0) - 5.28_wp * adu * tsk )
	330
	331	IF ( tsk == 36._wp ) tsk = 36.01_wp
	332	tcore(7) = c(0) / ( 5.28_wp * adu + c(1) * 6.3_wp / 3600._wp ) + tsk
	333	tcore(3) = c(0) / ( 5.28_wp * adu + ( c(1) * 6.3_wp / 3600._wp ) / &
	334	( 1._wp + 0.5_wp * ( 34._wp - tsk ) ) ) + tsk
	335	IF ( c(10) >= 0._wp ) THEN
	336	tcore(6) = ( - c(9) - c(10)*0.5_wp ) / ( 2._wp c(4) )
	337	tcore(1) = ( - c(9) + c(10)*0.5_wp ) / ( 2._wp c(4) )
	338	END IF
	339
	340	IF ( c(8) >= 0._wp ) THEN
	341	tcore(2) = ( - c(6) + ABS( c(8) ) ** 0.5_wp ) / ( 2._wp * c(4) )
	342	tcore(5) = ( - c(6) - ABS( c(8) ) ** 0.5_wp ) / ( 2._wp * c(4) )
	343	tcore(4) = c(0) / ( 5.28_wp * adu + c(1) * 1._wp / 40._wp ) + tsk
	344	END IF
	345
	346	!-- Transpiration
	347	tbody = 0.1_wp * tsk + 0.9_wp * tcore(j)
	348	swm = 304.94_wp * ( tbody - 36.6_wp ) * adu / 3600000._wp
	349	vpts = 6.11_wp * 10._wp*( 7.45_wp tsk / ( 235._wp + tsk ) )
	350
	351	IF ( tbody <= 36.6_wp ) swm = 0._wp
	352
	353	sw = swm
	354	eswphy = - sw * evap
	355	he = 0.633_wp * hc / ( p * cair )
	356	fec = 1._wp / ( 1._wp + 0.92_wp * hc * rcl )
	357	eswpot = he * ( vpa - vpts ) * adu * evap * fec
	358	wetsk = eswphy / eswpot
	359
	360	IF ( wetsk > 1._wp ) wetsk = 1._wp
	361
	362	eswdif = eswphy - eswpot
	363
	364	IF ( eswdif <= 0._wp ) esw = eswpot
	365	IF ( eswdif > 0._wp ) esw = eswphy
	366	IF ( esw > 0._wp ) esw = 0._wp
	367
	368	!-- Diffusion
	369	rdsk = 0.79_wp * 10._wp ** 7._wp
	370	rdcl = 0._wp
	371	ed = evap / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( vpa - &
	372	vpts )
	373
	374	!-- Max vb
	375	vb1 = 34._wp - tsk
	376	vb2 = tcore(j) - 36.6_wp
	377
	378	IF ( vb2 < 0._wp ) vb2 = 0._wp
	379	IF ( vb1 < 0._wp ) vb1 = 0._wp
	380	vb = ( 6.3_wp + 75._wp * vb2 ) / ( 1._wp + 0.5_wp * vb1 )
	381
	382	!-- Energy ballence
	383	enbal = h + ed + ere + esw + csum + rsum + food
	384
	385	!-- Clothing temperature
	386	xx = 0.001_wp
	387	IF ( count1 == 0_iwp ) xx = 1._wp
	388	IF ( count1 == 1_iwp ) xx = 0.1_wp
	389	IF ( count1 == 2_iwp ) xx = 0.01_wp
	390	IF ( count1 == 3_iwp ) xx = 0.001_wp
	391
	392	IF ( enbal > 0._wp ) tcl = tcl + xx
	393	IF ( enbal < 0._wp ) tcl = tcl - xx
	394
	395	skipIncreaseCount = .FALSE.
	396	IF ( ( (enbal <= 0._wp ) .AND. (enbal2 > 0._wp ) ) .OR. &
	397	( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) ) THEN
	398	skipIncreaseCount = .TRUE.
	399	ELSE
	400	enbal2 = enbal
	401	count3 = count3 + 1_iwp
	402	END IF
	403
	404	IF ( ( count3 > 200_iwp ) .OR. skipIncreaseCount ) THEN
	405	IF ( count1 < 3_iwp ) THEN
	406	count1 = count1 + 1_iwp
	407	enbal2 = 0._wp
	408	ELSE
	409	EXIT
	410	END IF
	411	END IF
	412	END DO
	413
	414	IF ( count1 == 3_iwp ) THEN
	415	SELECT CASE ( j )
	416	CASE ( 2, 5)
	417	IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND. &
	418	( tsk <= 34.050_wp ) ) ) CYCLE
	419	CASE ( 6, 1 )
	420	IF ( c(10) < 0._wp ) CYCLE
	421	IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND. &
	422	( tsk > 33.850_wp ) ) ) CYCLE
	423	CASE ( 3 )
	424	IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND. &
	425	( tsk <= 34.000_wp ) ) ) CYCLE
	426	CASE ( 7 )
	427	IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND. &
	428	( tsk > 34.000_wp ) ) ) CYCLE
	429	CASE default
	430	END SELECT
	431	END IF
	432
	433	IF ( ( j /= 4_iwp ) .AND. ( vb >= 91._wp ) ) CYCLE
	434	IF ( ( j == 4_iwp ) .AND. ( vb < 89._wp ) ) CYCLE
	435	IF ( vb > 90._wp ) vb = 90._wp
	436
	437	!-- Loses by water
	438	ws = sw * 3600._wp * 1000._wp
	439	IF ( ws > 2000._wp ) ws = 2000._wp
	440	wd = ed / evap * 3600._wp * ( -1000._wp )
	441	wr = erel / evap * 3600._wp * ( -1000._wp )
	442
	443	wsum = ws + wr + wd
	444
	445	RETURN
	446	END DO
	447	END SUBROUTINE heat_exch
	448
	449	!------------------------------------------------------------------------------!
	450	! Description:
	451	! ------------
	452	!> Calculate PET
	453	!------------------------------------------------------------------------------!
	454	SUBROUTINE pet_iteration ( acl, adu, aeff, esw, facl, feff, p, rdcl, &
	455	rdsk, rtv, ta, tcl, tsk, tx, vpts, wetsk )
	456
	457	!-- Input arguments:
	458	REAL(wp), INTENT( IN ) :: acl !< clothing surface area (mÂ²)
	459	REAL(wp), INTENT( IN ) :: adu !< Du-Bois area (mÂ²)
	460	REAL(wp), INTENT( IN ) :: esw !< energy-loss through sweat evap. (W)
	461	REAL(wp), INTENT( IN ) :: facl !< surface area extension through clothing (factor)
	462	REAL(wp), INTENT( IN ) :: feff !< surface modification by posture (factor)
	463	REAL(wp), INTENT( IN ) :: p !< air pressure (hPa)
	464	REAL(wp), INTENT( IN ) :: rdcl !< diffusion resistence of clothing (factor)
	465	REAL(wp), INTENT( IN ) :: rdsk !< diffusion resistence of skin (factor)
	466	REAL(wp), INTENT( IN ) :: rtv !<
	467	REAL(wp), INTENT( IN ) :: ta !< air temperature (Â°C)
	468	REAL(wp), INTENT( IN ) :: tcl !< clothing temperature (Â°C)
	469	REAL(wp), INTENT( IN ) :: tsk !< skin temperature (Â°C)
	470	REAL(wp), INTENT( IN ) :: vpts !< sat. vapor pressure over skin (hPa)
	471	REAL(wp), INTENT( IN ) :: wetsk !< fraction of wet skin (dimensionless)
	472
	473	!-- Output arguments:
	474	REAL(wp), INTENT( OUT ) :: aeff !< effective surface area (mÂ²)
	475	REAL(wp), INTENT( OUT ) :: tx !< PET (Â°C)
	476
	477	!-- Internal variables
	478	REAL ( wp ) :: cbare
	479	REAL ( wp ) :: cclo
	480	REAL ( wp ) :: csum
	481	REAL ( wp ) :: ed
	482	REAL ( wp ) :: enbal
	483	REAL ( wp ) :: enbal2
	484	REAL ( wp ) :: ere
	485	REAL ( wp ) :: erel
	486	REAL ( wp ) :: eres
	487	REAL ( wp ) :: hc
	488	REAL ( wp ) :: rbare
	489	REAL ( wp ) :: rclo
	490	REAL ( wp ) :: rsum
	491	REAL ( wp ) :: tex
	492	REAL ( wp ) :: vpex
	493	REAL ( wp ) :: xx
	494
	495	INTEGER ( iwp ) :: count1
	496	INTEGER ( iwp ) :: i
	497
	498	tx = ta
	499	enbal2 = 0._wp
	500
	501	DO count1 = 0, 3
	502	DO i = 1, 125 ! 500 / 4
	503	hc = 2.67_wp + 6.5_wp * 0.1_wp ** 0.67_wp
	504	hc = hc * ( p / po ) ** 0.55_wp
	505
	506	!-- Radiation
	507	aeff = adu * feff
	508	rbare = aeff * ( 1._wp - facl ) * emsk * sigm * &
	509	( ( tx + cels_offs ) 4._wp - ( tsk + cels_offs ) 4._wp )
	510	rclo = feff * acl * emcl * sigm * &
	511	( ( tx + cels_offs ) 4._wp - ( tcl + cels_offs ) 4._wp )
	512	rsum = rbare + rclo
	513
	514	!-- Covection
	515	cbare = hc * ( tx - tsk ) * adu * ( 1._wp - facl )
	516	cclo = hc * ( tx - tcl ) * acl
	517	csum = cbare + cclo
	518
	519	!-- Diffusion
	520	ed = evap / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( 12._wp - &
	521	vpts )
	522
	523	!-- Respiration
	524	tex = 0.47_wp * tx + 21._wp
	525	eres = cair * ( tx - tex ) * rtv
	526	vpex = 6.11_wp * 10._wp ** ( 7.45_wp * tex / ( 235._wp + tex ) )
	527	erel = 0.623_wp * evap / p * ( 12._wp - vpex ) * rtv
	528	ere = eres + erel
	529
	530	!-- Energy ballance
	531	enbal = h + ed + ere + esw + csum + rsum
	532
	533	!-- Iteration concerning ta
	534	IF ( count1 == 0_iwp ) xx = 1._wp
	535	IF ( count1 == 1_iwp ) xx = 0.1_wp
	536	IF ( count1 == 2_iwp ) xx = 0.01_wp
	537	IF ( count1 == 3_iwp ) xx = 0.001_wp
	538	IF ( enbal > 0._wp ) tx = tx - xx
	539	IF ( enbal < 0._wp ) tx = tx + xx
	540	IF ( ( enbal <= 0._wp ) .AND. ( enbal2 > 0._wp ) ) EXIT
	541	IF ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) EXIT
	542
	543	enbal2 = enbal
	544	END DO
	545	END DO
	546	END SUBROUTINE pet_iteration
	547
	548	END MODULE biometeorology_pet_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |