Changeset 3753 for palm


Ignore:
Timestamp:
Feb 19, 2019 2:48:54 PM (5 years ago)
Author:
dom_dwd_user
Message:

biometeorology_mod.f90:
(C) Added automatic setting of mrt_nlevels in case it was not part of radiation_parameters namelist (or set to 0 accidentially).
(C) Minor speed improvoemnts in perceived temperature calculations.
(C) Perceived temperature regression arrays now declared as PARAMETERs.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/biometeorology_mod.f90

    r3750 r3753  
    2727! -----------------
    2828! $Id$
     29! - Added automatic setting of mrt_nlevels in case it was not part of
     30! radiation_parameters namelist (or set to 0 accidentially).
     31! - Minor speed improvoemnts in perceived temperature calculations.
     32! - Perceived temperature regression arrays now declared as PARAMETERs.
     33!
     34! 3750 2019-02-19 07:29:39Z dom_dwd_user
    2935! - Added addittional safety meassures to bio_calculate_thermal_index_maps.
    3036! - Replaced several REAL (un-)equality comparisons.
     
    836842!--    Further checks if thermal comfort output is desired.
    837843    IF ( thermal_comfort  .AND.  unit == 'degree_C' )  THEN
    838 
    839 !
    840 !--    Break if required modules "radiation" and "humidity" are not avalable.
     844!
     845!--    Break if required modules "radiation" is not avalable.
    841846       IF ( .NOT.  radiation )  THEN
    842847          message_string = 'output of "' // TRIM( var ) // '" require'         &
     
    845850          unit = 'illegal'
    846851       ENDIF
    847        IF ( .NOT.  humidity )  THEN
    848           message_string = 'The estimation of thermal comfort requires '    // &
    849                            'air humidity information, but humidity module ' // &
    850                            'is disabled!'
    851           CALL message( 'check_parameters', 'PA0561', 1, 2, 0, 6, 0 )
    852           unit = 'illegal'
    853        ENDIF
    854        IF ( mrt_nlevels == 0 )  THEN
    855           message_string = 'output of "' // TRIM( var ) // '" require'         &
    856                            // 's mrt_nlevels > 0'
    857           CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
    858           unit = 'illegal'
     852!
     853!--    All "thermal_comfort" outputs except from 'bio_mrt' will also need
     854!--    humidity input. Check also for that.
     855       IF ( TRIM( var ) /= 'bio_mrt' )  THEN
     856          IF ( .NOT.  humidity )  THEN
     857             message_string = 'The estimation of thermal comfort '    //       &
     858                              'requires air humidity information, but ' //     &
     859                              'humidity module is disabled!'
     860             CALL message( 'check_parameters', 'PA0561', 1, 2, 0, 6, 0 )
     861             unit = 'illegal'
     862          ENDIF
    859863       ENDIF
    860864
     
    12141218    height = 0.0_wp
    12151219
    1216     bio_cell_level = INT ( 1.099_wp / dz(1) )
     1220    bio_cell_level = INT( 1.099_wp / dz(1) )
    12171221    bio_output_height = bio_output_height + bio_cell_level * dz(1)
    1218 
     1222!
     1223!-- Set radiation level if not done by user
     1224    IF ( mrt_nlevels == 0 )  THEN
     1225       mrt_nlevels = bio_cell_level + 1_iwp
     1226    ENDIF
    12191227!
    12201228!-- Init UVEM and load lookup tables
     
    15741582!-- Determine cell level closest to 1.1m above ground
    15751583!   by making use of truncation due to int cast
    1576     k = get_topography_top_index_ji(j, i, 's') + bio_cell_level  !< Vertical cell center closest to 1.1m
     1584    k = INT( get_topography_top_index_ji(j, i, 's') + bio_cell_level )  !< Vertical cell center closest to 1.1m
    15771585
    15781586!
     
    15891597       ta = bio_fill_value
    15901598       IF ( ALLOCATED( pt_av ) )  THEN
    1591           ta = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
     1599          ta = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * ( k + .5_wp ) ) - degc_to_k
    15921600       ENDIF
    15931601
     
    22092217    REAL(wp) ::  clon           !< clo for neutral conditions   (clo)
    22102218    REAL(wp) ::  ireq_minimal   !< Minimal required clothing insulation (clo)
    2211 !     REAL(wp) ::  clo_fanger     !< clo for fanger subroutine, unused
    22122219    REAL(wp) ::  pmv_w          !< Fangers predicted mean vote for winter clothing
    22132220    REAL(wp) ::  pmv_s          !< Fangers predicted mean vote for summer clothing
     
    22162223    REAL(wp) ::  d_std          !< factor to threshold for sultriness
    22172224    REAL(wp) ::  pmvs           !< pred. mean vote considering sultrieness
    2218     REAL(wp) ::  top            !< Gagge's operative temperatures (degree_C)
    22192225
    22202226    INTEGER(iwp) :: ncount      !< running index
     
    22442250!--    First guess: winter clothing insulation: cold stress
    22452251       clo = wclo
    2246        CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, top )
     2252       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
    22472253       pmv_w = pmva
    22482254
     
    22512257!--       Case summer clothing insulation: heat load ?
    22522258          clo = sclo
    2253           CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
    2254              top )
     2259          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
    22552260          pmv_s = pmva
    22562261          IF ( pmva <= 0._wp )  THEN
     
    22592264!--          Between winter and summer set values
    22602265             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
    2261                 pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
     2266                pmv_s, wclo, pmv_w, eps, pmva, ncount, clo )
    22622267             IF ( ncount < 0_iwp )  THEN
    22632268                nerr = -1_iwp
     
    22672272             clo = 0.5_wp
    22682273             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
    2269                            pmva, top )
     2274                           pmva )
    22702275          ENDIF
    22712276       ELSE IF ( pmva < -0.11_wp )  THEN
    22722277          clo = 1.75_wp
    2273           CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
    2274              top )
     2278          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
    22752279       ENDIF
    22762280    ELSE
     
    22782282!--    First guess: summer clothing insulation: heat load
    22792283       clo = sclo
    2280        CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, top )
     2284       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
    22812285       pmv_s = pmva
    22822286
     
    22852289!--       Case winter clothing insulation: cold stress ?
    22862290          clo = wclo
    2287           CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
    2288              top )
     2291          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
    22892292          pmv_w = pmva
    22902293
     
    22942297!--          between winter and summer set values
    22952298             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
    2296                                pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
     2299                               pmv_s, wclo, pmv_w, eps, pmva, ncount, clo )
    22972300             IF ( ncount < 0_iwp )  THEN
    22982301                nerr = -1_iwp
     
    23022305             clo = 1.75_wp
    23032306             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
    2304                            pmva, top )
     2307                           pmva )
    23052308          ENDIF
    23062309       ELSE IF ( pmva > 0.06_wp )  THEN
    23072310          clo = 0.5_wp
    2308           CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
    2309              top )
     2311          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
    23102312       ENDIF
    23112313
     
    23332335!
    23342336!--    Required clothing insulation (ireq) is exclusively defined for
    2335 !--    operative temperatures (top) less 10 (C) for a
     2337!--    perceived temperatures (perct) less 10 (C) for a
    23362338!--    reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
    23372339       clon = ireq_neutral ( perct_ij, ireq_minimal, nerr )
     
    23932395!
    23942396!-- Saturation water vapour pressure
    2395     svp_ta = 6.1078_wp * EXP ( b * ta / ( c + ta ) )
     2397    svp_ta = 6.1078_wp * EXP( b * ta / ( c + ta ) )
    23962398
    23972399 END SUBROUTINE saturation_vapor_pressure
     
    24052407!------------------------------------------------------------------------------!
    24062408 SUBROUTINE iso_ridder( ta, tmrt, vp, ws, pair, actlev, eta, sclo,             &
    2407                        pmv_s, wclo, pmv_w, eps, pmva, top, nerr,               &
     2409                       pmv_s, wclo, pmv_w, eps, pmva, nerr,               &
    24082410                       clo_res )
    24092411
     
    24272429!-- Output variables of argument list:
    24282430    REAL(wp), INTENT ( OUT ) :: pmva     !< 0 (set to zero, because clo is evaluated for comfort)
    2429     REAL(wp), INTENT ( OUT ) :: top      !< Operative temperature (degC) at found root of Fanger's PMV
    24302431    REAL(wp), INTENT ( OUT ) :: clo_res  !< Resulting clothing insulation value (clo)
    24312432    INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error status / quality flag
     
    24712472          x_average = 0.5_wp * ( x_lower + x_upper )
    24722473          CALL fanger ( ta, tmrt, vp, ws, pair, x_average, actlev, eta,        &
    2473                         y_average, top )
    2474           sroot = SQRT ( y_average**2 - y_lower * y_upper )
     2474                        y_average )
     2475          sroot = SQRT( y_average**2 - y_lower * y_upper )
    24752476          IF ( ABS( sroot ) < 0.00001_wp )  THEN
    24762477             clo_res = x_average
     
    24802481          x_new = x_average + ( x_average - x_lower ) *                        &
    24812482                      ( SIGN ( 1._wp, y_lower - y_upper ) * y_average / sroot )
    2482           IF ( ABS ( x_new - x_ridder ) <= eps )  THEN
     2483          IF ( ABS( x_new - x_ridder ) <= eps )  THEN
    24832484             clo_res = x_ridder
    24842485             nerr       = j
     
    24872488          x_ridder = x_new
    24882489          CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta,         &
    2489                         y_new, top )
     2490                        y_new )
    24902491          IF ( ABS( y_new ) < 0.00001_wp )  THEN
    24912492             clo_res = x_ridder
     
    25112512             RETURN
    25122513          ENDIF
    2513           IF ( ABS ( x_upper - x_lower ) <= eps )  THEN
     2514          IF ( ABS( x_upper - x_lower ) <= eps )  THEN
    25142515             clo_res = x_ridder
    25152516             nerr    = j
     
    25742575!> The case of free convection (ws < 0.1 m/s) is dealt with ws = 0.1 m/s
    25752576!------------------------------------------------------------------------------!
    2576  SUBROUTINE fanger( ta, tmrt, pa, in_ws, pair, in_clo, actlev, eta, pmva, top )
     2577 SUBROUTINE fanger( ta, tmrt, pa, in_ws, pair, in_clo, actlev, eta, pmva )
    25772578
    25782579    IMPLICIT NONE
     
    25922593                                         !< dimensionless) according to Fanger corresponding to meteorological
    25932594                                         !< (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta)
    2594     REAL(wp), INTENT ( OUT ) ::  top     !< operative temperature (degC)
    25952595!
    25962596!-- Internal variables
     
    26312631!
    26322632!-- Heat_convection = forced convection
    2633     heat_convection = 12.1_wp * SQRT ( ws * pair / 1013.25_wp )
     2633    heat_convection = 12.1_wp * SQRT( ws * pair / 1013.25_wp )
    26342634!
    26352635!-- Activity = inner heat produktion per standardized surface
     
    26562656!-- Empiric factor for the adaption of the heat ballance equation
    26572657!-- to the psycho-physical scale (Equ. 40 in FANGER)
    2658     z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp )
     2658    z1 = ( .303_wp * EXP( -.036_wp * actlev ) + .0275_wp )
    26592659!
    26602660!-- Water vapour diffution through the skin
     
    26712671    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + degc_to_k )**4 - ( tmrt +        &
    26722672       degc_to_k )**4 )
    2673     IF ( ABS ( t_clothing - tmrt ) > 0._wp )  THEN
     2673    IF ( ABS( t_clothing - tmrt ) > 0._wp )  THEN
    26742674       hr = z5 / f_cl / ( t_clothing - tmrt )
    26752675    ELSE
     
    26822682!-- Predicted Mean Vote
    26832683    pmva = z1 * ( activity - z2 - z3 - z4 - z5 - z6 )
    2684 !
    2685 !-- Operative temperatur
    2686     top = ( hr * tmrt + heat_convection * ta ) / ( hr + heat_convection )
    26872684
    26882685 END SUBROUTINE fanger
     
    27312728    REAL(wp) ::  dpmv_2       !<
    27322729    REAL(wp) ::  pmvs         !<
    2733     REAL(wp) ::  bpmv(0:7)    !<
    2734     REAL(wp) ::  bpa_p50(0:7) !<
    2735     REAL(wp) ::  bpa(0:7)     !<
    2736     REAL(wp) ::  bapa(0:7)    !<
    2737     REAL(wp) ::  bdapa(0:7)   !<
    2738     REAL(wp) ::  bsqvel(0:7)  !<
    2739     REAL(wp) ::  bta(0:7)     !<
    2740     REAL(wp) ::  bdtmrt(0:7)  !<
    2741     REAL(wp) ::  aconst(0:7)  !<
    27422730    INTEGER(iwp) :: nreg      !<
    27432731
    2744     DATA bpmv     /                                                            &
    2745      -0.0556602_wp, -0.1528680_wp, -0.2336104_wp, -0.2789387_wp, -0.3551048_wp,&
    2746      -0.4304076_wp, -0.4884961_wp, -0.4897495_wp /
    2747     DATA bpa_p50 /                                                             &
    2748      -0.1607154_wp, -0.4177296_wp, -0.4120541_wp, -0.0886564_wp, +0.4285938_wp,&
    2749      +0.6281256_wp, +0.5067361_wp, +0.3965169_wp /
    2750     DATA bpa     /                                                             &
    2751      +0.0580284_wp, +0.0836264_wp, +0.1009919_wp, +0.1020777_wp, +0.0898681_wp,&
    2752      +0.0839116_wp, +0.0853258_wp, +0.0866589_wp /
    2753     DATA bapa    /                                                             &
    2754      -1.7838788_wp, -2.9306231_wp, -1.6350334_wp, +0.6211547_wp, +3.3918083_wp,&
    2755      +5.5521025_wp, +8.4897418_wp, +16.6265851_wp /
    2756     DATA bdapa   /                                                             &
    2757      +1.6752720_wp, +2.7379504_wp, +1.2940526_wp, -1.0985759_wp, -3.9054732_wp,&
    2758      -6.0403012_wp, -8.9437119_wp, -17.0671201_wp /
    2759     DATA bsqvel  /                                                             &
    2760      -0.0315598_wp, -0.0286272_wp, -0.0009228_wp, +0.0483344_wp, +0.0992366_wp,&
    2761      +0.1491379_wp, +0.1951452_wp, +0.2133949_wp /
    2762     DATA bta     /                                                             &
    2763      +0.0953986_wp, +0.1524760_wp, +0.0564241_wp, -0.0893253_wp, -0.2398868_wp,&
    2764      -0.3515237_wp, -0.5095144_wp, -0.9469258_wp /
    2765     DATA bdtmrt  /                                                             &
    2766      -0.0004672_wp, -0.0000514_wp, -0.0018037_wp, -0.0049440_wp, -0.0069036_wp,&
    2767      -0.0075844_wp, -0.0079602_wp, -0.0089439_wp /
    2768     DATA aconst  /                                                             &
    2769      +1.8686215_wp, +3.4260713_wp, +2.0116185_wp, -0.7777552_wp, -4.6715853_wp,&
    2770      -7.7314281_wp, -11.7602578_wp, -23.5934198_wp /
     2732!
     2733!-- Regression coefficients:
     2734    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpmv = (/                          &
     2735     -0.0556602_wp, -0.1528680_wp, -0.2336104_wp, -0.2789387_wp,               &
     2736     -0.3551048_wp, -0.4304076_wp, -0.4884961_wp, -0.4897495_wp /)
     2737
     2738    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpa_p50 = (/                       &
     2739     -0.1607154_wp, -0.4177296_wp, -0.4120541_wp, -0.0886564_wp,               &
     2740      0.4285938_wp,  0.6281256_wp,  0.5067361_wp,  0.3965169_wp /)
     2741
     2742    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpa = (/                           &
     2743      0.0580284_wp,  0.0836264_wp,  0.1009919_wp,  0.1020777_wp,               &
     2744      0.0898681_wp,  0.0839116_wp,  0.0853258_wp,  0.0866589_wp /)
     2745
     2746    REAL(wp), DIMENSION(0:7), PARAMETER ::  bapa = (/                          &
     2747     -1.7838788_wp, -2.9306231_wp, -1.6350334_wp,   0.6211547_wp,              &
     2748      3.3918083_wp,  5.5521025_wp,  8.4897418_wp,  16.6265851_wp /)
     2749
     2750    REAL(wp), DIMENSION(0:7), PARAMETER ::  bdapa = (/                         &
     2751      1.6752720_wp,  2.7379504_wp,  1.2940526_wp,  -1.0985759_wp,              &
     2752     -3.9054732_wp, -6.0403012_wp, -8.9437119_wp, -17.0671201_wp /)
     2753
     2754    REAL(wp), DIMENSION(0:7), PARAMETER ::  bsqvel = (/                        &
     2755     -0.0315598_wp, -0.0286272_wp, -0.0009228_wp,  0.0483344_wp,               &
     2756      0.0992366_wp,  0.1491379_wp,  0.1951452_wp,  0.2133949_wp /)
     2757
     2758    REAL(wp), DIMENSION(0:7), PARAMETER ::  bta = (/                           &
     2759      0.0953986_wp,  0.1524760_wp,  0.0564241_wp, -0.0893253_wp,               &
     2760     -0.2398868_wp, -0.3515237_wp, -0.5095144_wp, -0.9469258_wp /)
     2761
     2762    REAL(wp), DIMENSION(0:7), PARAMETER ::  bdtmrt = (/                        &
     2763     -0.0004672_wp, -0.0000514_wp, -0.0018037_wp, -0.0049440_wp,               &
     2764     -0.0069036_wp, -0.0075844_wp, -0.0079602_wp, -0.0089439_wp /)
     2765
     2766    REAL(wp), DIMENSION(0:7), PARAMETER ::  aconst = (/                        &
     2767      1.8686215_wp,  3.4260713_wp,   2.0116185_wp,  -0.7777552_wp,             &
     2768     -4.6715853_wp, -7.7314281_wp, -11.7602578_wp, -23.5934198_wp /)
     2769
     2770
    27712771!
    27722772!-- Test for compliance with regression range
     
    28002800!
    28012801!--    Natural logarithm of pa
    2802        apa = LOG ( pa )
     2802       apa = LOG( pa )
    28032803    ELSE
    28042804       apa = -5._wp
     
    28082808    pa_p50   = 0.5_wp * svp_ta
    28092809    IF ( pa_p50 > 0._wp  .AND.  pa > 0._wp )  THEN
    2810        dapa   = apa - LOG ( pa_p50 )
     2810       dapa   = apa - LOG( pa_p50 )
    28112811       pa_p50 = pa / pa_p50
    28122812    ELSE
     
    28172817!-- Square root of wind velocity
    28182818    IF ( ws >= 0._wp )  THEN
    2819        sqvel = SQRT ( ws )
     2819       sqvel = SQRT( ws )
    28202820    ELSE
    28212821       sqvel = 0._wp
     
    28262826!
    28272827!-- Select the valid regression coefficients
    2828     nreg = INT ( pmv )
     2828    nreg = INT( pmv )
    28292829    IF ( nreg < 0_iwp )  THEN
    28302830!
     
    28362836    IF ( weight < 0._wp )  weight = 0._wp
    28372837    IF ( nreg > 5_iwp )  THEN
    2838        ! nreg=6
    28392838       nreg  = 5_iwp
    28402839       weight   = pmv - 5._wp
     
    28452844    ENDIF
    28462845!
    2847 !-- Regression valid for 0. <= pmv <= 6.
     2846!-- Regression valid for 0. <= pmv <= 6., bounds are checked above
    28482847    dpmv_1 =                                                                   &
    28492848       + bpa(nreg) * pa                                                        &
     
    28572856       + aconst(nreg)
    28582857
    2859     dpmv_2 = 0._wp
    2860     IF ( nreg < 6_iwp )  THEN
    2861        dpmv_2 =                                                                &
    2862           + bpa(nreg+1_iwp)     * pa                                           &
    2863           + bpmv(nreg+1_iwp)    * pmv                                          &
    2864           + bapa(nreg+1_iwp)    * apa                                          &
    2865           + bta(nreg+1_iwp)     * ta                                           &
    2866           + bdtmrt(nreg+1_iwp)  * dtmrt                                        &
    2867           + bdapa(nreg+1_iwp)   * dapa                                         &
    2868           + bsqvel(nreg+1_iwp)  * sqvel                                        &
    2869           + bpa_p50(nreg+1_iwp) * pa_p50                                       &
    2870           + aconst(nreg+1_iwp)
    2871     ENDIF
     2858!    dpmv_2 = 0._wp
     2859!    IF ( nreg < 6_iwp )  THEN  !< nreg is always <= 5, see above
     2860    dpmv_2 =                                                                   &
     2861       + bpa(nreg+1_iwp)     * pa                                              &
     2862       + bpmv(nreg+1_iwp)    * pmv                                             &
     2863       + bapa(nreg+1_iwp)    * apa                                             &
     2864       + bta(nreg+1_iwp)     * ta                                              &
     2865       + bdtmrt(nreg+1_iwp)  * dtmrt                                           &
     2866       + bdapa(nreg+1_iwp)   * dapa                                            &
     2867       + bsqvel(nreg+1_iwp)  * sqvel                                           &
     2868       + bpa_p50(nreg+1_iwp) * pa_p50                                          &
     2869       + aconst(nreg+1_iwp)
     2870!    ENDIF
    28722871!
    28732872!-- Calculate pmv modification
     
    28802879       IF ( pmvs > -0.11_wp )  THEN
    28812880!
    2882 !--       Threshold from FUNCTION perct_regression for winter clothing insulation
     2881!--       Threshold from perct_regression for winter clothing insulation
    28832882          deltapmv = deltapmv + 0.11_wp
    28842883       ELSE
     
    29152914!
    29162915!-- Types of coefficients mean deviation: third order polynomial
    2917     REAL(wp), PARAMETER ::  dperctka = +7.5776086_wp
     2916    REAL(wp), PARAMETER ::  dperctka =  7.5776086_wp
    29182917    REAL(wp), PARAMETER ::  dperctkb = -0.740603_wp
    2919     REAL(wp), PARAMETER ::  dperctkc = +0.0213324_wp
     2918    REAL(wp), PARAMETER ::  dperctkc =  0.0213324_wp
    29202919    REAL(wp), PARAMETER ::  dperctkd = -0.00027797237_wp
    29212920!
    29222921!-- Types of coefficients mean deviation plus standard deviation
    29232922!-- regression coefficients: third order polynomial
    2924     REAL(wp), PARAMETER ::  dperctsa = +0.0268918_wp
    2925     REAL(wp), PARAMETER ::  dperctsb = +0.0465957_wp
     2923    REAL(wp), PARAMETER ::  dperctsa =  0.0268918_wp
     2924    REAL(wp), PARAMETER ::  dperctsb =  0.0465957_wp
    29262925    REAL(wp), PARAMETER ::  dperctsc = -0.00054709752_wp
    2927     REAL(wp), PARAMETER ::  dperctsd = +0.0000063714823_wp
     2926    REAL(wp), PARAMETER ::  dperctsd =  0.0000063714823_wp
    29282927!
    29292928!-- Factor to mean standard deviation defining SIGNificance for
     
    29322931!
    29332932!-- Initialise
    2934     sultr_res = +99._wp
     2933    sultr_res = 99._wp
    29352934    dperctm   = 0._wp
    29362935    dperctstd = 999999._wp
     
    29522951!-- Value of the FUNCTION
    29532952    sultr_res = dperctm + faktor * dperctstd
    2954     IF ( ABS ( sultr_res ) > 99._wp )  sultr_res = +99._wp
     2953    IF ( ABS( sultr_res ) > 99._wp )  sultr_res = +99._wp
    29552954
    29562955 END SUBROUTINE calc_sultr
     
    29832982    REAL(wp) ::  pmv_cross(2)
    29842983    REAL(wp) ::  reg_a(3)
    2985     REAL(wp) ::  coeff(3,5)
    2986     REAL(wp) ::  r_nenner
    2987     REAL(wp) ::  pmvc
    2988     REAL(wp) ::  dtmrt
    2989     REAL(wp) ::  sqrt_ws
    2990     INTEGER(iwp) ::  i
    2991     INTEGER(iwp) ::  i_bin
    2992 !
    2993 !-- Coefficient of the 3 regression lines
    2994 !     1:const   2:*pmvc  3:*ta      4:*sqrt_ws  5:*dtmrt
    2995     coeff(1,1:5) =                                                             &
    2996        (/ +0.161_wp, +0.130_wp, -1.125E-03_wp, +1.106E-03_wp, -4.570E-04_wp /)
    2997     coeff(2,1:5) =                                                             &
    2998        (/ 0.795_wp, 0.713_wp, -8.880E-03_wp, -1.803E-03_wp, -2.816E-03_wp/)
    2999     coeff(3,1:5) =                                                             &
    3000        (/ +0.05761_wp, +0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp /)
     2984    REAL(wp) ::  r_denominator  !< the regression equations denominator
     2985    REAL(wp) ::  dtmrt          !< delta mean radiant temperature
     2986    REAL(wp) ::  sqrt_ws        !< sqare root of wind speed
     2987    INTEGER(iwp) ::  i          !< running index
     2988    INTEGER(iwp) ::  i_bin      !< result row number
     2989
     2990!    REAL(wp) ::  coeff(3,5)  !< unsafe! array is (re-)writable!
     2991!    coeff(1,1:5) =                                                             &
     2992!       (/ +0.161_wp,   +0.130_wp, -1.125E-03_wp, +1.106E-03_wp, -4.570E-04_wp /)
     2993!    coeff(2,1:5) =                                                             &
     2994!       (/  0.795_wp,    0.713_wp, -8.880E-03_wp, -1.803E-03_wp, -2.816E-03_wp /)
     2995!    coeff(3,1:5) =                                                             &
     2996!       (/ +0.05761_wp, +0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp /)
     2997
     2998!
     2999!-- Coefficient of the 3 regression lines:
     3000!      1:const      2:*pmva    3:*ta          4:*sqrt_ws     5:*dtmrt
     3001    REAL(wp), DIMENSION(1:3,1:5), PARAMETER ::  coeff = RESHAPE( (/            &
     3002        0.161_wp,   0.130_wp, -1.125E-03_wp,  1.106E-03_wp, -4.570E-04_wp,     &
     3003        0.795_wp,   0.713_wp, -8.880E-03_wp, -1.803E-03_wp, -2.816E-03_wp,     &
     3004        0.05761_wp, 0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp      &
     3005       /), SHAPE(coeff), order=(/ 2, 1 /) )
    30013006!
    30023007!-- Initialise
    30033008    nerr           = 0_iwp
    30043009    dpmv_cold_res  = 0._wp
    3005     pmvc           = pmva
    30063010    dtmrt          = tmrt - ta
    30073011    sqrt_ws        = ws
    3008     IF ( sqrt_ws < 0.10_wp )  THEN
    3009        sqrt_ws = 0.10_wp
     3012    IF ( sqrt_ws < 0.1_wp )  THEN
     3013       sqrt_ws = 0.1_wp
    30103014    ELSE
    3011        sqrt_ws = SQRT ( sqrt_ws )
     3015       sqrt_ws = SQRT( sqrt_ws )
    30123016    ENDIF
    30133017
    30143018    delta_cold = 0._wp
    3015     pmv_cross = pmvc
    3016 
     3019    pmv_cross = pmva
     3020
     3021!
     3022!-- Determine regression constant for given meteorological conditions
    30173023    DO  i = 1, 3
    3018 !
    3019 !--    Regression constant for given meteorological variables
    30203024       reg_a(i) = coeff(i,1) + coeff(i,3) * ta + coeff(i,4) *                  &
    30213025                  sqrt_ws + coeff(i,5)*dtmrt
    3022        delta_cold(i) = reg_a(i) + coeff(i,2) * pmvc
     3026       delta_cold(i) = reg_a(i) + coeff(i,2) * pmva
    30233027    ENDDO
    30243028!
    30253029!-- Intersection points of regression lines in terms of Fanger's PMV
    30263030    DO  i = 1, 2
    3027        r_nenner = coeff(i,2) - coeff(i+1,2)
    3028        IF ( ABS ( r_nenner ) > 0.00001_wp )  THEN
    3029           pmv_cross(i) = ( reg_a(i+1) - reg_a(i) ) / r_nenner
     3031       r_denominator = coeff(i,2) - coeff(i+1,2)
     3032       IF ( ABS( r_denominator ) > 0.00001_wp )  THEN
     3033          pmv_cross(i) = ( reg_a(i+1) - reg_a(i) ) / r_denominator
    30303034       ELSE
    30313035          nerr = 1_iwp
     
    30333037       ENDIF
    30343038    ENDDO
    3035 
     3039!
     3040!-- Select result row number
    30363041    i_bin = 3
    30373042    DO  i = 1, 2
     
    30443049!-- Adjust to operative temperature scaled according
    30453050!-- to classical PMV (Fanger)
    3046     dpmv_cold_res = delta_cold(i_bin) - dpmv_adj(pmva)
     3051    dpmv_cold_res = delta_cold(i_bin) - dpmv_cold_adj(pmva)
    30473052
    30483053 END SUBROUTINE dpmv_cold
     
    30513056! Description:
    30523057! ------------
    3053 !> Calculates the summand dpmv_adj adjusting to the operative temperature
    3054 !> scaled according to classical PMV (Fanger)
    3055 !> Reference environment: v_1m = 0.10 m/s, dTMRT = 0 K, r.h. = 50 %
    3056 !------------------------------------------------------------------------------!
    3057  REAL(wp) FUNCTION dpmv_adj( pmva )
     3058!> Calculates the summand dpmv_cold_adj adjusting to the operative temperature
     3059!> scaled according to classical PMV (Fanger) for cold conditions.
     3060!> Valid for reference environment: v (1m) = 0.10 m/s, dTMRT = 0 K, r.h. = 50 %
     3061!------------------------------------------------------------------------------!
     3062 REAL(wp) FUNCTION dpmv_cold_adj( pmva )
    30583063
    30593064    IMPLICIT NONE
    30603065
    3061     REAL(wp), INTENT ( IN ) ::  pmva
    3062     INTEGER(iwp), PARAMETER ::  n_bin = 3
    3063     INTEGER(iwp), PARAMETER ::  n_ca = 0
    3064     INTEGER(iwp), PARAMETER ::  n_ce = 3
    3065     REAL(wp), dimension (n_bin, n_ca:n_ce) ::  coef
    3066 
    3067     REAL(wp)      ::  pmv
    3068     INTEGER(iwp)  ::  i, i_bin
    3069 !
    3070 !--                                 range_1        range_2        range_3
    3071     DATA (coef(i,0), i = 1, n_bin) /0.0941540_wp, -0.1506620_wp, -0.0871439_wp/
    3072     DATA (coef(i,1), i = 1, n_bin) /0.0783162_wp, -1.0612651_wp,  0.1695040_wp/
    3073     DATA (coef(i,2), i = 1, n_bin) /0.1350144_wp, -1.0049144_wp, -0.0167627_wp/
    3074     DATA (coef(i,3), i = 1, n_bin) /0.1104037_wp, -0.2005277_wp, -0.0003230_wp/
    3075 
    3076     IF ( pmva <= -2.1226_wp )  THEN
    3077        i_bin = 3_iwp
    3078     ELSE IF ( pmva <= -1.28_wp )  THEN
    3079        i_bin = 2_iwp
    3080     ELSE
    3081        i_bin = 1_iwp
    3082     ENDIF
    3083 
    3084     dpmv_adj   = coef(i_bin,n_ca)
    3085     pmv        = 1._wp
    3086 
    3087     DO  i = n_ca + 1, n_ce
    3088        pmv      = pmv * pmva
    3089        dpmv_adj = dpmv_adj + coef(i_bin,i) * pmv
     3066    REAL(wp), INTENT ( IN ) ::  pmva        !< (adjusted) Predicted Mean Vote
     3067
     3068    REAL(wp)      ::  pmv     !< pmv-part of the regression
     3069    INTEGER(iwp)  ::  i       !< running index
     3070    INTEGER(iwp)  ::  thr     !< thermal range
     3071!
     3072!-- Provide regression coefficients for three thermal ranges:
     3073!--    slightly cold  cold           very cold
     3074    REAL(wp), DIMENSION(1:3,0:3), PARAMETER ::  coef = RESHAPE( (/             &
     3075       0.0941540_wp, -0.1506620_wp, -0.0871439_wp,                             &
     3076       0.0783162_wp, -1.0612651_wp,  0.1695040_wp,                             &
     3077       0.1350144_wp, -1.0049144_wp, -0.0167627_wp,                             &
     3078       0.1104037_wp, -0.2005277_wp, -0.0003230_wp                              &
     3079       /), SHAPE(coef), order=(/ 2, 1 /) )
     3080!
     3081!-- Select thermal range
     3082    IF ( pmva <= -2.1226_wp )  THEN     !< very cold
     3083       thr = 3_iwp
     3084    ELSE IF ( pmva <= -1.28_wp )  THEN  !< cold
     3085       thr = 2_iwp
     3086    ELSE                                !< slightly cold
     3087       thr = 1_iwp
     3088    ENDIF
     3089!
     3090!-- Initialize
     3091    dpmv_cold_adj = coef(thr,0)
     3092    pmv           = 1._wp
     3093!
     3094!-- Calculate pmv adjustment (dpmv_cold_adj)
     3095    DO  i = 1, 3
     3096       pmv           = pmv * pmva
     3097       dpmv_cold_adj = dpmv_cold_adj + coef(thr,i) * pmv
    30903098    ENDDO
     3099
    30913100    RETURN
    3092  END FUNCTION dpmv_adj
     3101 END FUNCTION dpmv_cold_adj
    30933102
    30943103!------------------------------------------------------------------------------!
     
    31113120!
    31123121!-- Type declaration for internal varables
    3113     REAL(wp)                     ::  top02
     3122    REAL(wp)                     ::  perct02
    31143123!
    31153124!-- Initialise
     
    31173126!
    31183127!-- Convert perceived temperature from basis 0.1 m/s to basis 0.2 m/s
    3119     top02 = 1.8788_wp + 0.9296_wp * perct_ij
     3128    perct02 = 1.8788_wp + 0.9296_wp * perct_ij
    31203129!
    31213130!-- IREQ neutral conditions (thermal comfort)
    3122     ireq_neutral = 1.62_wp - 0.0564_wp * top02
     3131    ireq_neutral = 1.62_wp - 0.0564_wp * perct02
    31233132!
    31243133!-- Regression only defined for perct <= 10 (degC)
     
    31313140!
    31323141!-- Minimal required clothing insulation: maximal acceptable body cooling
    3133     ireq_minimal = 1.26_wp - 0.0588_wp * top02
     3142    ireq_minimal = 1.26_wp - 0.0588_wp * perct02
    31343143    IF ( nerr > 0_iwp )  THEN
    31353144       ireq_minimal = ireq_neutral
     
    31693178!
    31703179!-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if
    3171 !-- height and weight be known. In: Arch. Int. Med.. 17, 1916, S. 863?871.
     3180!-- height and weight be known. In: Arch. Int. Med.. 17, 1916, pp. 863:871.
    31723181    surf = 0.007184_wp * height**0.725_wp * weight**0.425_wp
    31733182    RETURN
     
    32763285    REAL(wp) ::  d_std
    32773286    REAL(wp) ::  pmvs
    3278     REAL(wp) ::  top
    32793287    REAL(wp) ::  a_surf
    32803288!     REAL(wp) ::  acti
     
    33243332!--          between winter and summer set values
    33253333             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta , sclo,     &
    3326                             pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
     3334                            pmv_s, wclo, pmv_w, eps, pmva, ncount, clo )
    33273335             IF ( ncount < 0_iwp )  THEN
    33283336                nerr = -1_iwp
     
    33653373!--          between winter and summer set values
    33663374             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
    3367                                pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
     3375                               pmv_s, wclo, pmv_w, eps, pmva, ncount, clo )
    33683376             IF ( ncount < 0_wp )  THEN
    33693377                nerr = -1_iwp
     
    34063414!
    34073415!--    Required clothing insulation (ireq) is exclusively defined for
    3408 !--    operative temperatures (top) less 10 (C) for a
     3416!--    perceived temperatures (ipt) less 10 (C) for a
    34093417!--    reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
    34103418       clon = ireq_neutral ( ipt, ireq_minimal, nerr )
     
    34933501!-- Determine perceived temperature by regression equation + adjustments
    34943502    CALL perct_regression( pmva, clo, ipt )
    3495 
     3503!
     3504!-- Consider cold conditions
    34963505    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
    34973506!
     
    35103519    ptc = ipt
    35113520    CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
    3512     sultrieness    = .FALSE.
    3513     d_std      = -99._wp
     3521    sultrieness = .FALSE.
     3522    d_std       = -99._wp
    35143523    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
    35153524!
    35163525!--    Adjust for warm/humid conditions according to Gagge 1986
    35173526       CALL saturation_vapor_pressure ( ta, svp_ta )
    3518        d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
    3519        pmvs   = pmva + d_pmv
     3527       d_pmv = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
     3528       pmvs  = pmva + d_pmv
    35203529       CALL perct_regression( pmvs, clo, ipt )
    35213530       IF ( sult_lim < 99._wp )  THEN
     
    36023611!
    36033612!-- Heat_convection = forced convection
    3604     heat_convection = 12.1_wp * SQRT ( ws * pair / 1013.25_wp )
     3613    heat_convection = 12.1_wp * SQRT( ws * pair / 1013.25_wp )
    36053614!
    36063615!-- Average skin temperature
     
    36183627    niter = INT( dt * 10._wp, KIND=iwp )
    36193628    IF ( niter < 1 )  niter = 1_iwp
    3620     adjustrate = 1._wp - EXP ( -1._wp * ( 10._wp / time_equil ) * dt )
     3629    adjustrate = 1._wp - EXP( -1._wp * ( 10._wp / time_equil ) * dt )
    36213630    IF ( adjustrate >= 1._wp )  adjustrate = 1._wp
    36223631    adjustrate_cloth = adjustrate * 30._wp
    36233632    t_clothing = t_cloth
    3624 
    3625     IF ( t_cloth <= -998.0_wp )  THEN  ! If initial run
     3633!
     3634!-- Set initial values for niter, adjustrates and t_clothing if this is the
     3635!-- first call
     3636    IF ( t_cloth <= -998._wp )  THEN  ! If initial run
    36263637       niter = 3_iwp
    36273638       adjustrate = 1._wp
     
    36293640       t_clothing = ta
    36303641    ENDIF
    3631 
     3642!
     3643!-- Update clothing temperature
    36323644    DO  i = 1, niter
    36333645       t_clothing = t_clothing - adjustrate_cloth * ( ( t_clothing +           &
     
    36383650!-- Empiric factor for the adaption of the heat ballance equation
    36393651!-- to the psycho-physical scale (Equ. 40 in FANGER)
    3640     z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp )
     3652    z1 = ( .303_wp * EXP( -.036_wp * actlev ) + .0275_wp )
    36413653!
    36423654!-- Water vapour diffution through the skin
     
    37423754
    37433755    CALL heat_exch( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, ht,  &
    3744             int_heat,mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa,      &
     3756            int_heat, mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa,     &
    37453757            vpts, wetsk )
    37463758
     
    37683780    REAL(wp), INTENT( IN )  ::  ht        !< Persons height           (m)
    37693781    REAL(wp), INTENT( IN )  ::  work      !< Work load                (W)
    3770     REAL(wp), INTENT( IN )  ::  eta       !< Work efficiency      (dimensionless)
     3782    REAL(wp), INTENT( IN )  ::  eta       !< Work efficiency     (dimensionless)
    37713783!
    37723784!-- Output arguments:
Note: See TracChangeset for help on using the changeset viewer.