Changeset 3753
 Timestamp:
 Feb 19, 2019 2:48:54 PM (2 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/biometeorology_mod.f90
r3750 r3753 27 27 !  28 28 ! $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 20190219 07:29:39Z dom_dwd_user 29 35 !  Added addittional safety meassures to bio_calculate_thermal_index_maps. 30 36 !  Replaced several REAL (un)equality comparisons. … … 836 842 ! Further checks if thermal comfort output is desired. 837 843 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. 841 846 IF ( .NOT. radiation ) THEN 842 847 message_string = 'output of "' // TRIM( var ) // '" require' & … … 845 850 unit = 'illegal' 846 851 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 859 863 ENDIF 860 864 … … 1214 1218 height = 0.0_wp 1215 1219 1216 bio_cell_level = INT 1220 bio_cell_level = INT( 1.099_wp / dz(1) ) 1217 1221 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 1219 1227 ! 1220 1228 ! Init UVEM and load lookup tables … … 1574 1582 ! Determine cell level closest to 1.1m above ground 1575 1583 ! 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.1m1584 k = INT( get_topography_top_index_ji(j, i, 's') + bio_cell_level ) !< Vertical cell center closest to 1.1m 1577 1585 1578 1586 ! … … 1589 1597 ta = bio_fill_value 1590 1598 IF ( ALLOCATED( pt_av ) ) THEN 1591 ta = pt_av(k,j,i)  ( 0.0098_wp * dz(1) * ( 1599 ta = pt_av(k,j,i)  ( 0.0098_wp * dz(1) * ( k + .5_wp ) )  degc_to_k 1592 1600 ENDIF 1593 1601 … … 2209 2217 REAL(wp) :: clon !< clo for neutral conditions (clo) 2210 2218 REAL(wp) :: ireq_minimal !< Minimal required clothing insulation (clo) 2211 ! REAL(wp) :: clo_fanger !< clo for fanger subroutine, unused2212 2219 REAL(wp) :: pmv_w !< Fangers predicted mean vote for winter clothing 2213 2220 REAL(wp) :: pmv_s !< Fangers predicted mean vote for summer clothing … … 2216 2223 REAL(wp) :: d_std !< factor to threshold for sultriness 2217 2224 REAL(wp) :: pmvs !< pred. mean vote considering sultrieness 2218 REAL(wp) :: top !< Gagge's operative temperatures (degree_C)2219 2225 2220 2226 INTEGER(iwp) :: ncount !< running index … … 2244 2250 ! First guess: winter clothing insulation: cold stress 2245 2251 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 ) 2247 2253 pmv_w = pmva 2248 2254 … … 2251 2257 ! Case summer clothing insulation: heat load ? 2252 2258 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 ) 2255 2260 pmv_s = pmva 2256 2261 IF ( pmva <= 0._wp ) THEN … … 2259 2264 ! Between winter and summer set values 2260 2265 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 ) 2262 2267 IF ( ncount < 0_iwp ) THEN 2263 2268 nerr = 1_iwp … … 2267 2272 clo = 0.5_wp 2268 2273 CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, & 2269 pmva , top)2274 pmva ) 2270 2275 ENDIF 2271 2276 ELSE IF ( pmva < 0.11_wp ) THEN 2272 2277 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 ) 2275 2279 ENDIF 2276 2280 ELSE … … 2278 2282 ! First guess: summer clothing insulation: heat load 2279 2283 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 ) 2281 2285 pmv_s = pmva 2282 2286 … … 2285 2289 ! Case winter clothing insulation: cold stress ? 2286 2290 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 ) 2289 2292 pmv_w = pmva 2290 2293 … … 2294 2297 ! between winter and summer set values 2295 2298 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 ) 2297 2300 IF ( ncount < 0_iwp ) THEN 2298 2301 nerr = 1_iwp … … 2302 2305 clo = 1.75_wp 2303 2306 CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, & 2304 pmva , top)2307 pmva ) 2305 2308 ENDIF 2306 2309 ELSE IF ( pmva > 0.06_wp ) THEN 2307 2310 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 ) 2310 2312 ENDIF 2311 2313 … … 2333 2335 ! 2334 2336 ! Required clothing insulation (ireq) is exclusively defined for 2335 ! operative temperatures (top) less 10 (C) for a2337 ! perceived temperatures (perct) less 10 (C) for a 2336 2338 ! reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s 2337 2339 clon = ireq_neutral ( perct_ij, ireq_minimal, nerr ) … … 2393 2395 ! 2394 2396 ! Saturation water vapour pressure 2395 svp_ta = 6.1078_wp * EXP 2397 svp_ta = 6.1078_wp * EXP( b * ta / ( c + ta ) ) 2396 2398 2397 2399 END SUBROUTINE saturation_vapor_pressure … … 2405 2407 !! 2406 2408 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, & 2408 2410 clo_res ) 2409 2411 … … 2427 2429 ! Output variables of argument list: 2428 2430 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 PMV2430 2431 REAL(wp), INTENT ( OUT ) :: clo_res !< Resulting clothing insulation value (clo) 2431 2432 INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error status / quality flag … … 2471 2472 x_average = 0.5_wp * ( x_lower + x_upper ) 2472 2473 CALL fanger ( ta, tmrt, vp, ws, pair, x_average, actlev, eta, & 2473 y_average , top)2474 sroot = SQRT 2474 y_average ) 2475 sroot = SQRT( y_average**2  y_lower * y_upper ) 2475 2476 IF ( ABS( sroot ) < 0.00001_wp ) THEN 2476 2477 clo_res = x_average … … 2480 2481 x_new = x_average + ( x_average  x_lower ) * & 2481 2482 ( SIGN ( 1._wp, y_lower  y_upper ) * y_average / sroot ) 2482 IF ( ABS 2483 IF ( ABS( x_new  x_ridder ) <= eps ) THEN 2483 2484 clo_res = x_ridder 2484 2485 nerr = j … … 2487 2488 x_ridder = x_new 2488 2489 CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta, & 2489 y_new , top)2490 y_new ) 2490 2491 IF ( ABS( y_new ) < 0.00001_wp ) THEN 2491 2492 clo_res = x_ridder … … 2511 2512 RETURN 2512 2513 ENDIF 2513 IF ( ABS 2514 IF ( ABS( x_upper  x_lower ) <= eps ) THEN 2514 2515 clo_res = x_ridder 2515 2516 nerr = j … … 2574 2575 !> The case of free convection (ws < 0.1 m/s) is dealt with ws = 0.1 m/s 2575 2576 !! 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 ) 2577 2578 2578 2579 IMPLICIT NONE … … 2592 2593 !< dimensionless) according to Fanger corresponding to meteorological 2593 2594 !< (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta) 2594 REAL(wp), INTENT ( OUT ) :: top !< operative temperature (degC)2595 2595 ! 2596 2596 ! Internal variables … … 2631 2631 ! 2632 2632 ! Heat_convection = forced convection 2633 heat_convection = 12.1_wp * SQRT 2633 heat_convection = 12.1_wp * SQRT( ws * pair / 1013.25_wp ) 2634 2634 ! 2635 2635 ! Activity = inner heat produktion per standardized surface … … 2656 2656 ! Empiric factor for the adaption of the heat ballance equation 2657 2657 ! to the psychophysical scale (Equ. 40 in FANGER) 2658 z1 = ( .303_wp * EXP 2658 z1 = ( .303_wp * EXP( .036_wp * actlev ) + .0275_wp ) 2659 2659 ! 2660 2660 ! Water vapour diffution through the skin … … 2671 2671 z5 = 3.96e8_wp * f_cl * ( ( t_clothing + degc_to_k )**4  ( tmrt + & 2672 2672 degc_to_k )**4 ) 2673 IF ( ABS 2673 IF ( ABS( t_clothing  tmrt ) > 0._wp ) THEN 2674 2674 hr = z5 / f_cl / ( t_clothing  tmrt ) 2675 2675 ELSE … … 2682 2682 ! Predicted Mean Vote 2683 2683 pmva = z1 * ( activity  z2  z3  z4  z5  z6 ) 2684 !2685 ! Operative temperatur2686 top = ( hr * tmrt + heat_convection * ta ) / ( hr + heat_convection )2687 2684 2688 2685 END SUBROUTINE fanger … … 2731 2728 REAL(wp) :: dpmv_2 !< 2732 2729 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) !<2742 2730 INTEGER(iwp) :: nreg !< 2743 2731 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 2771 2771 ! 2772 2772 ! Test for compliance with regression range … … 2800 2800 ! 2801 2801 ! Natural logarithm of pa 2802 apa = LOG 2802 apa = LOG( pa ) 2803 2803 ELSE 2804 2804 apa = 5._wp … … 2808 2808 pa_p50 = 0.5_wp * svp_ta 2809 2809 IF ( pa_p50 > 0._wp .AND. pa > 0._wp ) THEN 2810 dapa = apa  LOG 2810 dapa = apa  LOG( pa_p50 ) 2811 2811 pa_p50 = pa / pa_p50 2812 2812 ELSE … … 2817 2817 ! Square root of wind velocity 2818 2818 IF ( ws >= 0._wp ) THEN 2819 sqvel = SQRT 2819 sqvel = SQRT( ws ) 2820 2820 ELSE 2821 2821 sqvel = 0._wp … … 2826 2826 ! 2827 2827 ! Select the valid regression coefficients 2828 nreg = INT 2828 nreg = INT( pmv ) 2829 2829 IF ( nreg < 0_iwp ) THEN 2830 2830 ! … … 2836 2836 IF ( weight < 0._wp ) weight = 0._wp 2837 2837 IF ( nreg > 5_iwp ) THEN 2838 ! nreg=62839 2838 nreg = 5_iwp 2840 2839 weight = pmv  5._wp … … 2845 2844 ENDIF 2846 2845 ! 2847 ! Regression valid for 0. <= pmv <= 6. 2846 ! Regression valid for 0. <= pmv <= 6., bounds are checked above 2848 2847 dpmv_1 = & 2849 2848 + bpa(nreg) * pa & … … 2857 2856 + aconst(nreg) 2858 2857 2859 dpmv_2 = 0._wp2860 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 2871 ENDIF2858 ! 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 2872 2871 ! 2873 2872 ! Calculate pmv modification … … 2880 2879 IF ( pmvs > 0.11_wp ) THEN 2881 2880 ! 2882 ! Threshold from FUNCTIONperct_regression for winter clothing insulation2881 ! Threshold from perct_regression for winter clothing insulation 2883 2882 deltapmv = deltapmv + 0.11_wp 2884 2883 ELSE … … 2915 2914 ! 2916 2915 ! Types of coefficients mean deviation: third order polynomial 2917 REAL(wp), PARAMETER :: dperctka = +7.5776086_wp2916 REAL(wp), PARAMETER :: dperctka = 7.5776086_wp 2918 2917 REAL(wp), PARAMETER :: dperctkb = 0.740603_wp 2919 REAL(wp), PARAMETER :: dperctkc = +0.0213324_wp2918 REAL(wp), PARAMETER :: dperctkc = 0.0213324_wp 2920 2919 REAL(wp), PARAMETER :: dperctkd = 0.00027797237_wp 2921 2920 ! 2922 2921 ! Types of coefficients mean deviation plus standard deviation 2923 2922 ! regression coefficients: third order polynomial 2924 REAL(wp), PARAMETER :: dperctsa = +0.0268918_wp2925 REAL(wp), PARAMETER :: dperctsb = +0.0465957_wp2923 REAL(wp), PARAMETER :: dperctsa = 0.0268918_wp 2924 REAL(wp), PARAMETER :: dperctsb = 0.0465957_wp 2926 2925 REAL(wp), PARAMETER :: dperctsc = 0.00054709752_wp 2927 REAL(wp), PARAMETER :: dperctsd = +0.0000063714823_wp2926 REAL(wp), PARAMETER :: dperctsd = 0.0000063714823_wp 2928 2927 ! 2929 2928 ! Factor to mean standard deviation defining SIGNificance for … … 2932 2931 ! 2933 2932 ! Initialise 2934 sultr_res = +99._wp2933 sultr_res = 99._wp 2935 2934 dperctm = 0._wp 2936 2935 dperctstd = 999999._wp … … 2952 2951 ! Value of the FUNCTION 2953 2952 sultr_res = dperctm + faktor * dperctstd 2954 IF ( ABS 2953 IF ( ABS( sultr_res ) > 99._wp ) sultr_res = +99._wp 2955 2954 2956 2955 END SUBROUTINE calc_sultr … … 2983 2982 REAL(wp) :: pmv_cross(2) 2984 2983 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.125E03_wp, +1.106E03_wp, 4.570E04_wp /) 2997 coeff(2,1:5) = & 2998 (/ 0.795_wp, 0.713_wp, 8.880E03_wp, 1.803E03_wp, 2.816E03_wp/) 2999 coeff(3,1:5) = & 3000 (/ +0.05761_wp, +0.458_wp, 1.829E02_wp, 5.577E03_wp, 1.970E03_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.125E03_wp, +1.106E03_wp, 4.570E04_wp /) 2993 ! coeff(2,1:5) = & 2994 ! (/ 0.795_wp, 0.713_wp, 8.880E03_wp, 1.803E03_wp, 2.816E03_wp /) 2995 ! coeff(3,1:5) = & 2996 ! (/ +0.05761_wp, +0.458_wp, 1.829E02_wp, 5.577E03_wp, 1.970E03_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.125E03_wp, 1.106E03_wp, 4.570E04_wp, & 3003 0.795_wp, 0.713_wp, 8.880E03_wp, 1.803E03_wp, 2.816E03_wp, & 3004 0.05761_wp, 0.458_wp, 1.829E02_wp, 5.577E03_wp, 1.970E03_wp & 3005 /), SHAPE(coeff), order=(/ 2, 1 /) ) 3001 3006 ! 3002 3007 ! Initialise 3003 3008 nerr = 0_iwp 3004 3009 dpmv_cold_res = 0._wp 3005 pmvc = pmva3006 3010 dtmrt = tmrt  ta 3007 3011 sqrt_ws = ws 3008 IF ( sqrt_ws < 0.1 0_wp ) THEN3009 sqrt_ws = 0.1 0_wp3012 IF ( sqrt_ws < 0.1_wp ) THEN 3013 sqrt_ws = 0.1_wp 3010 3014 ELSE 3011 sqrt_ws = SQRT 3015 sqrt_ws = SQRT( sqrt_ws ) 3012 3016 ENDIF 3013 3017 3014 3018 delta_cold = 0._wp 3015 pmv_cross = pmvc 3016 3019 pmv_cross = pmva 3020 3021 ! 3022 ! Determine regression constant for given meteorological conditions 3017 3023 DO i = 1, 3 3018 !3019 ! Regression constant for given meteorological variables3020 3024 reg_a(i) = coeff(i,1) + coeff(i,3) * ta + coeff(i,4) * & 3021 3025 sqrt_ws + coeff(i,5)*dtmrt 3022 delta_cold(i) = reg_a(i) + coeff(i,2) * pmv c3026 delta_cold(i) = reg_a(i) + coeff(i,2) * pmva 3023 3027 ENDDO 3024 3028 ! 3025 3029 ! Intersection points of regression lines in terms of Fanger's PMV 3026 3030 DO i = 1, 2 3027 r_ nenner = coeff(i,2)  coeff(i+1,2)3028 IF ( ABS ( r_nenner ) > 0.00001_wp ) THEN3029 pmv_cross(i) = ( reg_a(i+1)  reg_a(i) ) / r_ nenner3031 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 3030 3034 ELSE 3031 3035 nerr = 1_iwp … … 3033 3037 ENDIF 3034 3038 ENDDO 3035 3039 ! 3040 ! Select result row number 3036 3041 i_bin = 3 3037 3042 DO i = 1, 2 … … 3044 3049 ! Adjust to operative temperature scaled according 3045 3050 ! 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) 3047 3052 3048 3053 END SUBROUTINE dpmv_cold … … 3051 3056 ! Description: 3052 3057 !  3053 !> Calculates the summand dpmv_ adj adjusting to the operative temperature3054 !> 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 ) 3058 3063 3059 3064 IMPLICIT NONE 3060 3065 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 !< pmvpart 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 3090 3098 ENDDO 3099 3091 3100 RETURN 3092 END FUNCTION dpmv_ adj3101 END FUNCTION dpmv_cold_adj 3093 3102 3094 3103 !! … … 3111 3120 ! 3112 3121 ! Type declaration for internal varables 3113 REAL(wp) :: top023122 REAL(wp) :: perct02 3114 3123 ! 3115 3124 ! Initialise … … 3117 3126 ! 3118 3127 ! Convert perceived temperature from basis 0.1 m/s to basis 0.2 m/s 3119 top02 = 1.8788_wp + 0.9296_wp * perct_ij3128 perct02 = 1.8788_wp + 0.9296_wp * perct_ij 3120 3129 ! 3121 3130 ! IREQ neutral conditions (thermal comfort) 3122 ireq_neutral = 1.62_wp  0.0564_wp * top023131 ireq_neutral = 1.62_wp  0.0564_wp * perct02 3123 3132 ! 3124 3133 ! Regression only defined for perct <= 10 (degC) … … 3131 3140 ! 3132 3141 ! Minimal required clothing insulation: maximal acceptable body cooling 3133 ireq_minimal = 1.26_wp  0.0588_wp * top023142 ireq_minimal = 1.26_wp  0.0588_wp * perct02 3134 3143 IF ( nerr > 0_iwp ) THEN 3135 3144 ireq_minimal = ireq_neutral … … 3169 3178 ! 3170 3179 ! 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. 3172 3181 surf = 0.007184_wp * height**0.725_wp * weight**0.425_wp 3173 3182 RETURN … … 3276 3285 REAL(wp) :: d_std 3277 3286 REAL(wp) :: pmvs 3278 REAL(wp) :: top3279 3287 REAL(wp) :: a_surf 3280 3288 ! REAL(wp) :: acti … … 3324 3332 ! between winter and summer set values 3325 3333 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 ) 3327 3335 IF ( ncount < 0_iwp ) THEN 3328 3336 nerr = 1_iwp … … 3365 3373 ! between winter and summer set values 3366 3374 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 ) 3368 3376 IF ( ncount < 0_wp ) THEN 3369 3377 nerr = 1_iwp … … 3406 3414 ! 3407 3415 ! Required clothing insulation (ireq) is exclusively defined for 3408 ! operative temperatures (top) less 10 (C) for a3416 ! perceived temperatures (ipt) less 10 (C) for a 3409 3417 ! reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s 3410 3418 clon = ireq_neutral ( ipt, ireq_minimal, nerr ) … … 3493 3501 ! Determine perceived temperature by regression equation + adjustments 3494 3502 CALL perct_regression( pmva, clo, ipt ) 3495 3503 ! 3504 ! Consider cold conditions 3496 3505 IF ( clo >= 1.75_wp .AND. pmva <= 0.11_wp ) THEN 3497 3506 ! … … 3510 3519 ptc = ipt 3511 3520 CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim ) 3512 sultrieness 3513 d_std = 99._wp3521 sultrieness = .FALSE. 3522 d_std = 99._wp 3514 3523 IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN 3515 3524 ! 3516 3525 ! Adjust for warm/humid conditions according to Gagge 1986 3517 3526 CALL saturation_vapor_pressure ( ta, svp_ta ) 3518 d_pmv 3519 pmvs 3527 d_pmv = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr ) 3528 pmvs = pmva + d_pmv 3520 3529 CALL perct_regression( pmvs, clo, ipt ) 3521 3530 IF ( sult_lim < 99._wp ) THEN … … 3602 3611 ! 3603 3612 ! Heat_convection = forced convection 3604 heat_convection = 12.1_wp * SQRT 3613 heat_convection = 12.1_wp * SQRT( ws * pair / 1013.25_wp ) 3605 3614 ! 3606 3615 ! Average skin temperature … … 3618 3627 niter = INT( dt * 10._wp, KIND=iwp ) 3619 3628 IF ( niter < 1 ) niter = 1_iwp 3620 adjustrate = 1._wp  EXP 3629 adjustrate = 1._wp  EXP( 1._wp * ( 10._wp / time_equil ) * dt ) 3621 3630 IF ( adjustrate >= 1._wp ) adjustrate = 1._wp 3622 3631 adjustrate_cloth = adjustrate * 30._wp 3623 3632 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 3626 3637 niter = 3_iwp 3627 3638 adjustrate = 1._wp … … 3629 3640 t_clothing = ta 3630 3641 ENDIF 3631 3642 ! 3643 ! Update clothing temperature 3632 3644 DO i = 1, niter 3633 3645 t_clothing = t_clothing  adjustrate_cloth * ( ( t_clothing + & … … 3638 3650 ! Empiric factor for the adaption of the heat ballance equation 3639 3651 ! to the psychophysical scale (Equ. 40 in FANGER) 3640 z1 = ( .303_wp * EXP 3652 z1 = ( .303_wp * EXP( .036_wp * actlev ) + .0275_wp ) 3641 3653 ! 3642 3654 ! Water vapour diffution through the skin … … 3742 3754 3743 3755 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, & 3745 3757 vpts, wetsk ) 3746 3758 … … 3768 3780 REAL(wp), INTENT( IN ) :: ht !< Persons height (m) 3769 3781 REAL(wp), INTENT( IN ) :: work !< Work load (W) 3770 REAL(wp), INTENT( IN ) :: eta !< Work efficiency 3782 REAL(wp), INTENT( IN ) :: eta !< Work efficiency (dimensionless) 3771 3783 ! 3772 3784 ! Output arguments:
Note: See TracChangeset
for help on using the changeset viewer.