Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

Last change on this file since 3456 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
Line 
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
548END MODULE biometeorology_pet_mod
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |