source: palm/trunk/SOURCE/biometeorology_ipt_mod.f90 @ 3518

Last change on this file since 3518 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: 21.0 KB
Line 
1!> @file biometeorology_ipt_mod.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 1997-2019, Deutscher Wetterdienst (DWD) /
18! German Meteorological Service (DWD) and
19! Leibniz Universitaet Hannover
20!--------------------------------------------------------------------------------!
21!
22! Current revisions: 001
23! -----------------
24!
25!
26! Former revisions: 001
27! -----------------
28! $Id$
29! Initial revision 001
30!
31!
32!
33! Authors:
34! --------
35! @author Dominik Froehlich <dominik.froehlich@mailbox.org>
36!
37!
38! Description:
39! ------------
40!> Module for the calculation of the dynamic thermal index iPT
41!>
42!> @todo
43!>
44!> @note nothing now
45!>
46!> @bug  no known bugs by now
47!------------------------------------------------------------------------------!
48 MODULE biometeorology_ipt_mod
49!
50!-- Load required variables from existing modules
51    USE kinds  !< to set precision of INTEGER and REAL arrays according to PALM
52
53    USE biometeorology_pt_mod,                                                 &
54       ONLY: saturation_vapor_pressure, deltapmv, dpmv_adj, dpmv_cold,         &
55             calc_sultr, pt_regression, ireq_neutral, iso_ridder
56
57    IMPLICIT NONE
58
59    PRIVATE  !-- No private interfaces --!
60
61    PUBLIC ipt_cycle, ipt_init
62
63!-- PALM interfaces:
64    INTERFACE ipt_cycle
65       MODULE PROCEDURE ipt_cycle
66    END INTERFACE ipt_cycle
67
68    INTERFACE ipt_init
69       MODULE PROCEDURE ipt_init
70    END INTERFACE ipt_init
71
72
73 CONTAINS
74
75
76
77!------------------------------------------------------------------------------!
78! Description:
79! ------------
80!> The SUBROUTINE surface area calculates the surface area of the individual
81!> according to its height (m), weight (kg), and age (y)
82!------------------------------------------------------------------------------!
83 SUBROUTINE surface_area ( height_cm, weight, age, surf )
84
85    IMPLICIT NONE
86
87    REAL(wp)    , INTENT(in)  :: weight
88    REAL(wp)    , INTENT(in)  :: height_cm
89    INTEGER(iwp), INTENT(in)  :: age
90    REAL(wp)    , INTENT(out) :: surf
91    REAL(wp)                  :: height
92
93    height = height_cm * 100._wp
94
95!-- According to Gehan-George, for children
96    IF ( age < 19_iwp ) THEN
97       IF ( age < 5_iwp ) THEN
98          surf = 0.02667_wp * height**0.42246_wp * weight**0.51456_wp
99          RETURN
100       END IF
101       surf = 0.03050_wp * height**0.35129_wp * weight**0.54375_wp
102       RETURN
103    END IF
104
105!-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if
106!   height and weight be known. In: Arch. Int. Med.. 17, 1916, S. 863?871.
107    surf = 0.007184_wp * height**0.725_wp * weight**0.425_wp
108    RETURN
109
110 END SUBROUTINE surface_area
111
112!------------------------------------------------------------------------------!
113! Description:
114! ------------
115!> The SUBROUTINE persdat calculates
116!>  - the total internal energy production = metabolic + workload,
117!>  - the total internal energy production for a standardized surface (actlev)
118!>  - the DuBois - area (a_surf [m2])
119!> from
120!>  - the persons age (years),
121!>  - weight (kg),
122!>  - height (m),
123!>  - sex (1 = male, 2 = female),
124!>  - work load (W)
125!> for a sample human.
126!------------------------------------------------------------------------------!
127 SUBROUTINE persdat ( age, weight, height, sex, work, a_surf, actlev )
128!
129    IMPLICIT NONE
130
131    REAL(wp), INTENT(in) ::  age
132    REAL(wp), INTENT(in) ::  weight
133    REAL(wp), INTENT(in) ::  height
134    REAL(wp), INTENT(in) ::  work
135    INTEGER(iwp), INTENT(in) ::  sex
136    REAL(wp), INTENT(out) ::  actlev
137    REAL(wp) ::  a_surf
138    REAL(wp) ::  energieumsatz
139    REAL(wp) ::  s
140    REAL(wp) ::  factor
141    REAL(wp) ::  heightundumsatz
142
143    CALL surface_area ( height, weight, INT( age ), a_surf )
144    s = height * 100._wp / ( weight ** ( 1._wp / 3._wp ) )
145    factor = 1._wp + .004_wp  * ( 30._wp - age )
146    heightundumsatz = 0.
147    IF ( sex == 1_iwp ) THEN
148       heightundumsatz = 3.45_wp * weight ** ( 3._wp / 4._wp ) * ( factor +    &
149                     .01_wp  * ( s - 43.4_wp ) )
150    ELSE IF ( sex == 2_iwp ) THEN
151       heightundumsatz = 3.19_wp * weight ** ( 3._wp / 4._wp ) * ( factor +    &
152                    .018_wp * ( s - 42.1_wp ) )
153    END IF
154
155    energieumsatz = work + heightundumsatz
156    actlev = energieumsatz / a_surf
157
158 END SUBROUTINE persdat
159
160
161!------------------------------------------------------------------------------!
162! Description:
163! ------------
164!> SUBROUTINE ipt_init
165!> initializes the instationary perceived temperature
166!------------------------------------------------------------------------------!
167
168 SUBROUTINE ipt_init (age, weight, height, sex, work, actlev, clo,             &
169     ta, vp, vau1m, tmrt, pb, dt, storage, t_clothing,                         &
170     ipt )
171
172    IMPLICIT NONE
173
174!-- Input parameters
175    REAL(wp), INTENT(in) ::  age        !< Persons age          (years)
176    REAL(wp), INTENT(in) ::  weight     !< Persons weight       (kg)
177    REAL(wp), INTENT(in) ::  height     !< Persons height       (m)
178    REAL(wp), INTENT(in) ::  work       !< Current workload     (W)
179    REAL(wp), INTENT(in) ::  ta         !< Air Temperature      (°C)
180    REAL(wp), INTENT(in) ::  vp         !< Vapor pressure       (hPa)
181    REAL(wp), INTENT(in) ::  vau1m      !< Wind speed in approx. 1.1m (m/s)
182    REAL(wp), INTENT(in) ::  tmrt       !< Mean radiant temperature   (°C)
183    REAL(wp), INTENT(in) ::  pb         !< Air pressure         (hPa)
184    REAL(wp), INTENT(in) ::  dt         !< Timestep             (s)
185    INTEGER(iwp), INTENT(in)  :: sex    !< Persons sex (1 = male, 2 = female)
186
187!-- Output parameters
188    REAL(wp), INTENT(out) ::  actlev
189    REAL(wp), INTENT(out) ::  clo
190    REAL(wp), INTENT(out) ::  storage
191    REAL(wp), INTENT(out) ::  t_clothing
192    REAL(wp), INTENT(out) ::  ipt
193
194!-- Internal variables
195    REAL(wp), PARAMETER :: eps = 0.0005_wp
196    REAL(wp), PARAMETER :: eta = 0._wp
197    REAL(wp) ::  uclo
198    REAL(wp) ::  oclo
199    REAL(wp) ::  d_pmv
200    REAL(wp) ::  svp_tt
201    REAL(wp) ::  sult_lim
202    REAL(wp) ::  dgtcm
203    REAL(wp) ::  dgtcstd
204    REAL(wp) ::  clon
205    REAL(wp) ::  ireq_minimal
206    REAL(wp) ::  clo_fanger
207    REAL(wp) ::  pmv_o
208    REAL(wp) ::  pmv_u
209    REAL(wp) ::  pmva
210    REAL(wp) ::  ptc
211    REAL(wp) ::  d_std
212    REAL(wp) ::  pmvs
213    REAL(wp) ::  top
214    REAL(wp) ::  a_surf
215    REAL(wp) ::  acti
216    INTEGER(iwp) ::  nzaehl
217    INTEGER(iwp) ::  nerr_kalt
218    INTEGER(iwp) ::  nerr
219
220    LOGICAL ::  sultrieness
221
222    storage = 0._wp
223    CALL persdat ( age, weight, height, sex, work, a_surf, actlev )
224
225!-- Initialise
226    t_clothing = -999.0_wp
227    ipt        = -999.0_wp
228    nerr       = 0_wp
229    nzaehl     = 0_wp
230    sultrieness    = .FALSE.
231!-- Tresholds: clothing insulation (account for model inaccuracies)
232!-- Summer clothing
233    uclo     = 0.44453_wp
234!-- Winter clothing
235    oclo     = 1.76267_wp
236
237!-- Decision: firstly calculate for winter or summer clothing
238    IF ( ta <= 10._wp ) THEN
239
240!--    First guess: winter clothing insulation: cold stress
241       clo = oclo
242       t_clothing = -999.0_wp  ! force initial run
243       CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,        &
244          t_clothing, storage, dt, pmva )
245       pmv_o = pmva
246
247       IF ( pmva > 0._wp ) THEN
248
249!--       Case summer clothing insulation: heat load ?           
250          clo = uclo
251          t_clothing = -999.0_wp  ! force initial run
252          CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,     &
253             t_clothing, storage, dt, pmva )
254          pmv_u = pmva
255          IF ( pmva <= 0._wp ) THEN
256!--          Case: comfort achievable by varying clothing insulation
257!            between winter and summer set values
258             CALL iso_ridder ( ta, tmrt, vp, vau1m, pb, actlev, eta , uclo,    &
259                            pmv_u, oclo, pmv_o, eps, pmva, top, nzaehl, clo )
260             IF ( nzaehl < 0_iwp ) THEN
261                nerr = -1_iwp
262                RETURN
263             ENDIF
264          ELSE IF ( pmva > 0.06_wp ) THEN
265             clo = 0.5_wp
266             t_clothing = -999.0_wp
267             CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,  &
268                t_clothing, storage, dt, pmva )
269          ENDIF
270       ELSE IF ( pmva < -0.11_wp ) THEN
271          clo = 1.75_wp
272          t_clothing = -999.0_wp
273          CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,     &
274             t_clothing, storage, dt, pmva )
275       ENDIF
276
277    ELSE
278
279!--    First guess: summer clothing insulation: heat load
280       clo = uclo
281       t_clothing = -999.0_wp
282       CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,        &
283          t_clothing, storage, dt, pmva )
284       pmv_u = pmva
285
286       IF ( pmva < 0._wp ) THEN
287
288!--       Case winter clothing insulation: cold stress ?
289          clo = oclo
290          t_clothing = -999.0_wp
291          CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,     &
292             t_clothing, storage, dt, pmva )
293          pmv_o = pmva
294
295          IF ( pmva >= 0._wp ) THEN
296
297!--          Case: comfort achievable by varying clothing insulation
298!            between winter and summer set values
299             CALL iso_ridder ( ta, tmrt, vp, vau1m, pb, actlev, eta, uclo,     &
300                               pmv_u, oclo, pmv_o, eps, pmva, top, nzaehl, clo )
301             IF ( nzaehl < 0_wp ) THEN
302                nerr = -1_iwp
303                RETURN
304             ENDIF
305          ELSE IF ( pmva < -0.11_wp ) THEN
306             clo = 1.75_wp
307             t_clothing = -999.0_wp
308             CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,  &
309                t_clothing, storage, dt, pmva )
310          ENDIF
311       ELSE IF ( pmva > 0.06_wp ) THEN
312          clo = 0.5_wp
313          t_clothing = -999.0_wp
314          CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,     &
315             t_clothing, storage, dt, pmva )
316       ENDIF
317
318    ENDIF
319
320!-- Determine perceived temperature by regression equation + adjustments
321    pmvs = pmva
322    CALL pt_regression ( pmva, clo, ipt )
323    ptc = ipt
324    IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
325!--    Adjust for cold conditions according to Gagge 1986
326       CALL dpmv_cold ( pmva, ta, vau1m, tmrt, nerr_kalt, d_pmv )
327       IF ( nerr_kalt > 0_iwp ) nerr = -5_iwp
328       pmvs = pmva - d_pmv
329       IF ( pmvs > -0.11_wp ) THEN
330          d_pmv  = 0._wp
331          pmvs   = -0.11_wp
332       ENDIF
333       CALL pt_regression ( pmvs, clo, ipt )
334    ENDIF
335    clo_fanger = clo
336    clon = clo
337    IF ( clo > 0.5_wp .AND. ipt <= 8.73_wp ) THEN
338!--    Required clothing insulation (ireq) is exclusively defined for
339!      operative temperatures (top) less 10 (C) for a
340!      reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
341       clon = ireq_neutral ( ipt, ireq_minimal, nerr )
342       clo = clon
343    ENDIF
344    CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
345    sultrieness    = .FALSE.
346    d_std      = -99._wp
347    IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
348!--    Adjust for warm/humid conditions according to Gagge 1986
349       CALL saturation_vapor_pressure ( ta, svp_tt )
350       d_pmv  = deltapmv ( pmva, ta, vp, svp_tt, tmrt, vau1m, nerr )
351       pmvs   = pmva + d_pmv
352       CALL pt_regression ( pmvs, clo, ipt )
353       IF ( sult_lim < 99._wp ) THEN
354          IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE.
355       ENDIF
356    ENDIF
357
358 
359 END SUBROUTINE ipt_init
360 
361!------------------------------------------------------------------------------!
362! Description:
363! ------------
364!> SUBROUTINE ipt_cycle
365!> Calculates one timestep for the instationary version of perceived
366!> temperature (iPT, °C) for
367!>  - standard measured/predicted meteorological values and TMRT
368!>    as input;
369!>  - regressions for determination of PT;
370!>  - adjustment to Gagge's PMV* (2-node-model, 1986) as base of PT
371!>    under warm/humid conditions (Icl= 0.50 clo) and under cold
372!>    conditions (Icl= 1.75 clo)
373!>
374!------------------------------------------------------------------------------!
375 SUBROUTINE ipt_cycle( ta, vp, vau1m, tmrt, pb, dt, storage, t_clothing, clo,  &
376     actlev, work, ipt )
377
378    IMPLICIT NONE
379
380!-- Type of input of the argument list
381    REAL(wp), INTENT ( IN )  ::  ta      !< Air temperature             (°C)
382    REAL(wp), INTENT ( IN )  ::  vp      !< Vapor pressure              (hPa)
383    REAL(wp), INTENT ( IN )  ::  tmrt    !< Mean radiant temperature    (°C)
384    REAL(wp), INTENT ( IN )  ::  vau1m   !< Wind speed                  (m/s)
385    REAL(wp), INTENT ( IN )  ::  pb      !< Air pressure                (hPa)
386    REAL(wp), INTENT ( IN )  ::  dt      !< Timestep                    (s)
387    REAL(wp), INTENT ( IN )  ::  clo     !< Clothing index              (no dim)
388    REAL(wp), INTENT ( IN )  ::  actlev  !< Internal heat production    (W)
389    REAL(wp), INTENT ( IN )  ::  work    !< Mechanical work load        (W)
390
391!-- In and output parameters
392    REAL(wp), INTENT (INOUT) ::  storage     !< Heat storage            (W)
393    REAL(wp), INTENT (INOUT) ::  t_clothing  !< Clothig temperature     (°C)
394
395!-- Type of output of the argument list
396    REAL(wp), INTENT ( OUT ) ::  ipt  !< Instationary perceived temperature (°C)
397
398!-- Type of program variables
399    REAL(wp), PARAMETER ::  eps = 0.0005_wp
400    REAL(wp) ::  d_pmv
401    REAL(wp) ::  pmv_o
402    REAL(wp) ::  pmv_s
403    REAL(wp) ::  storage_rate
404    REAL(wp) ::  ener_in
405    REAL(wp) ::  svp_tt
406    REAL(wp) ::  sult_lim
407    REAL(wp) ::  dgtcm
408    REAL(wp) ::  dgtcstd
409    REAL(wp) ::  clon
410    REAL(wp) ::  ireq_minimal
411    REAL(wp) ::  clo_fanger
412    REAL(wp) ::  pmv_u
413    REAL(wp) ::  half_hour_frac
414    REAL(wp) ::  storage_delta
415    REAL(wp) ::  storage_delta_max
416    REAL(wp) ::  pmva
417    REAL(wp) ::  ptc
418    REAL(wp) ::  d_std
419    REAL(wp) ::  pmvs
420    REAL(wp) ::  top
421    INTEGER(iwp) ::  nzaehl
422    INTEGER(iwp) ::  nerr_kalt
423    INTEGER(iwp) ::  nerr
424
425    LOGICAL ::  sultrieness
426
427!-- Initialise
428    ipt = -999.0_wp
429
430    nerr     = 0_iwp
431    nzaehl   = 0_iwp
432    sultrieness  = .FALSE.
433
434!-- Determine pmv_adjusted for current conditions
435    CALL fanger_s_acti ( ta, tmrt, vp, vau1m, pb, clo, actlev, work,           &
436       t_clothing, storage, dt, pmva )
437    pmv_s = pmva  !< PMV considering storage
438
439!-- Determine perceived temperature by regression equation + adjustments
440    CALL pt_regression ( pmva, clo, ipt )
441
442    IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
443!--    Adjust for cold conditions according to Gagge 1986
444       CALL dpmv_cold ( pmva, ta, vau1m, tmrt, nerr_kalt, d_pmv )
445       IF ( nerr_kalt > 0_iwp ) nerr = -5_iwp
446       pmvs = pmva - d_pmv
447       IF ( pmvs > -0.11_wp ) THEN
448          d_pmv  = 0._wp
449          pmvs   = -0.11_wp
450       ENDIF
451       CALL pt_regression ( pmvs, clo, ipt )
452    ENDIF
453
454!-- Consider sultriness if appropriate
455    ptc = ipt
456    CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
457    sultrieness    = .FALSE.
458    d_std      = -99._wp
459    IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
460!--    Adjust for warm/humid conditions according to Gagge 1986
461       CALL saturation_vapor_pressure ( ta, svp_tt )
462       d_pmv  = deltapmv ( pmva, ta, vp, svp_tt, tmrt, vau1m, nerr )
463       pmvs   = pmva + d_pmv
464       CALL pt_regression ( pmvs, clo, ipt )
465       IF ( sult_lim < 99._wp ) THEN
466          IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE.
467       ENDIF
468    ENDIF
469
470 END SUBROUTINE ipt_cycle
471
472!------------------------------------------------------------------------------!
473! Description:
474! ------------
475!> SUBROUTINE fanger_s calculates the
476!> actual Predicted Mean Vote (dimensionless) according
477!> to Fanger corresponding to meteorological (tt,tmrt,pa,vau,pb)
478!> and individual variables (clo, actlev, eta) considering a storage
479!> and clothing temperature for a given timestep.
480!------------------------------------------------------------------------------!
481 SUBROUTINE fanger_s_acti( tt, tmrt, pa, in_vau, pb, in_clo, actlev,           &
482    activity, t_cloth, s, dt, pmva )
483
484    IMPLICIT NONE
485
486!--  Input argument types
487    REAL(wp), INTENT ( IN )  ::  tt       !< Air temperature          (°C)
488    REAL(wp), INTENT ( IN )  ::  tmrt     !< Mean radiant temperature (°C)
489    REAL(wp), INTENT ( IN )  ::  pa       !< Vapour pressure          (hPa)
490    REAL(wp), INTENT ( IN )  ::  pb       !< Air pressure             (hPa)
491    REAL(wp), INTENT ( IN )  ::  in_vau   !< Wind speed               (m/s)
492    REAL(wp), INTENT ( IN )  ::  actlev   !< Metabolic + work energy  (W/m²)
493    REAL(wp), INTENT ( IN )  ::  dt       !< Timestep                 (s)
494    REAL(wp), INTENT ( IN )  ::  activity !< Work load                (W/m²)
495    REAL(wp), INTENT ( IN )  ::  in_clo   !< Clothing index (clo)     (no dim)
496
497!-- Output argument types
498    REAL(wp), INTENT ( OUT ) ::  pmva  !< actual Predicted Mean Vote (no dim)
499
500    REAL(wp), INTENT (INOUT) ::  s  !< storage var. of energy balance (W/m2)
501    REAL(wp), INTENT (INOUT) ::  t_cloth  !< clothing temperature (°C)
502
503!-- Internal variables
504    REAL(wp) ::  f_cl
505    REAL(wp) ::  heat_convection
506    REAL(wp) ::  t_skin_aver
507    REAL(wp) ::  bc
508    REAL(wp) ::  cc
509    REAL(wp) ::  dc
510    REAL(wp) ::  ec
511    REAL(wp) ::  gc
512    REAL(wp) ::  hr
513    REAL(wp) ::  clo
514    REAL(wp) ::  vau
515    REAL(wp) ::  z1
516    REAL(wp) ::  z2
517    REAL(wp) ::  z3
518    REAL(wp) ::  z4
519    REAL(wp) ::  z5
520    REAL(wp) ::  z6
521    REAL(wp) ::  en
522    REAL(wp) ::  d_s
523    REAL(wp) ::  t_clothing   !< Clothing index value   (dimensionless)
524    REAL(wp) ::  adjustrate
525    REAL(wp) ::  adjustrate_cloth
526    REAL(wp), PARAMETER  ::  time_equil = 7200._wp
527    INTEGER(iwp) ::  i
528    INTEGER(iwp) ::  niter  !< Running index
529
530
531!-- Clo must be > 0. to avoid div. by 0!
532    clo = in_clo
533    IF ( clo < 001._wp ) clo = .001_wp
534
535!-- Increase in surface due to clothing
536    f_cl = 1._wp + .15_wp * clo
537
538!-- Case of free convection (vau < 0.1 m/s ) not considered
539    vau = in_vau
540    IF ( vau < .1_wp ) THEN
541       vau = .1_wp
542    ENDIF
543
544!-- Heat_convection = forced convection
545    heat_convection = 12.1_wp * SQRT ( vau * pb / 1013.25_wp )
546
547!-- Average skin temperature
548    t_skin_aver = 35.7_wp - .0275_wp * activity
549
550!-- Calculation of constants for evaluation below
551    bc = .155_wp * clo * 3.96_wp * 10._wp**( -8._wp ) * f_cl
552    cc = f_cl * heat_convection
553    ec = .155_wp * clo
554    dc = ( 1._wp + ec * cc ) / bc
555    gc = ( t_skin_aver + bc * ( tmrt + 273.2_wp )**4._wp + ec * cc * tt ) / bc
556
557!-- Calculation of clothing surface temperature (t_clothing) based on
558!   newton-approximation with air temperature as initial guess
559    niter = dt
560    adjustrate = 1._wp - EXP ( -1._wp * ( 10._wp / time_equil ) * dt )
561    IF ( adjustrate >= 1._wp ) adjustrate = 1._wp
562    adjustrate_cloth = adjustrate * 30._wp
563    t_clothing = t_cloth
564
565    IF ( t_cloth <= -998.0_wp ) THEN  ! If initial run
566       niter = 3_wp
567       adjustrate = 1._wp
568       adjustrate_cloth = 1._wp
569       t_clothing = tt
570    ENDIF
571
572    DO i = 1, niter
573       t_clothing = t_clothing - adjustrate_cloth * ( ( t_clothing +           &
574          273.2_wp )**4._wp + t_clothing *                                     &
575          dc - gc ) / ( 4._wp * ( t_clothing + 273.2_wp )**3._wp + dc )
576    ENDDO
577
578!-- Empiric factor for the adaption of the heat ballance equation
579!   to the psycho-physical scale (Equ. 40 in FANGER)
580    z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp )
581
582!-- Water vapour diffution through the skin
583    z2 = .31_wp * ( 57.3_wp - .07_wp * activity-pa )
584
585!-- Sweat evaporation from the skin surface
586    z3 = .42_wp * ( activity - 58._wp )
587
588!-- Loss of latent heat through respiration
589    z4 = .0017_wp * actlev * ( 58.7_wp - pa ) + .0014_wp * actlev *            &
590      ( 34._wp - tt )
591
592!-- Loss of radiational heat
593    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + 273.2_wp )**4 - ( tmrt +         &
594       273.2_wp )**4 )
595
596!-- Heat loss through forced convection
597    z6 = cc * ( t_clothing - tt )
598
599!-- Write together as energy ballance
600    en = activity - z2 - z3 - z4 - z5 - z6
601
602!-- Manage storage
603    d_s = adjustrate * en + ( 1._wp - adjustrate ) * s
604
605!-- Predicted Mean Vote
606    pmva = z1 * d_s
607
608!-- Update storage
609    s = d_s
610    t_cloth = t_clothing
611
612 END SUBROUTINE fanger_s_acti
613
614
615 END MODULE biometeorology_ipt_mod
Note: See TracBrowser for help on using the repository browser.