Changeset 3749
- Timestamp:
- Feb 19, 2019 7:18:38 AM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/biometeorology_mod.f90 (modified) (41 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/biometeorology_mod.f90
r3742 r3749 27 27 ! ----------------- 28 28 ! $Id$ 29 ! - Added addittional safety meassures to bio_calculate_thermal_index_maps. 30 ! - Replaced several REAL (un-)equality comparisons. 31 ! 32 ! 3742 2019-02-14 11:25:22Z dom_dwd_user 29 33 ! - Allocation of the input _av grids was moved to the "sum" section of 30 34 ! bio_3d_data_averaging to make sure averaging is only done once! … … 523 527 ENDIF 524 528 529 ! 530 !-- u_av, v_av and w_av are always allocated 525 531 IF ( .NOT. ALLOCATED( u_av ) ) THEN 526 532 ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) … … 866 872 SUBROUTINE bio_check_parameters 867 873 868 USE control_parameters, &869 ONLY: message_string870 874 871 875 IMPLICIT NONE … … 1455 1459 END SUBROUTINE bio_wrd_local 1456 1460 1457 1458 1461 !------------------------------------------------------------------------------! 1459 1462 ! Description: … … 1584 1587 ! 1585 1588 !-- Calculate ta from Tp assuming dry adiabatic laps rate 1586 ta = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * ( k + .5_wp ) ) - degc_to_k 1589 ta = bio_fill_value 1590 IF ( ALLOCATED( pt_av ) ) THEN 1591 ta = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * ( k + .5_wp ) ) - degc_to_k 1592 ENDIF 1587 1593 1588 1594 vp = bio_fill_value … … 1591 1597 ENDIF 1592 1598 1593 ws = ( 0.5_wp * ABS( u_av(k_wind, j, i) + u_av(k_wind, j, i+1) ) + & 1594 0.5_wp * ABS( v_av(k_wind, j, i) + v_av(k_wind, j+1, i) ) + & 1595 0.5_wp * ABS( w_av(k_wind, j, i) + w_av(k_wind+1, j, i) ) ) 1599 ws = bio_fill_value 1600 IF ( ALLOCATED( u_av ) .AND. ALLOCATED( v_av ) .AND. & 1601 ALLOCATED( w_av ) ) THEN 1602 ws = ( 0.5_wp * ABS( u_av(k_wind,j,i) + u_av(k_wind,j,i+1) ) + & 1603 0.5_wp * ABS( v_av(k_wind,j,i) + v_av(k_wind,j+1,i) ) + & 1604 0.5_wp * ABS( w_av(k_wind,j,i) + w_av(k_wind+1,j,i) ) ) 1605 ENDIF 1596 1606 ELSE 1597 1607 ! … … 1604 1614 ENDIF 1605 1615 1606 ws = ( 0.5_wp * ABS( u(k_wind, j, i) + u(k_wind, j, i+1) ) +&1607 0.5_wp * ABS( v(k_wind, j, i) + v(k_wind, j+1, i) ) +&1608 0.5_wp * ABS( w(k_wind, j, i) + w(k_wind+1, j,i) ) )1616 ws = ( 0.5_wp * ABS( u(k_wind,j,i) + u(k_wind,j,i+1) ) + & 1617 0.5_wp * ABS( v(k_wind,j,i) + v(k_wind,j+1,i) ) + & 1618 0.5_wp * ABS( w(k_wind,j,i) + w(k_wind+1,j,i) ) ) 1609 1619 1610 1620 ENDIF … … 1641 1651 !> time_integration.f90: 1065ff 1642 1652 !------------------------------------------------------------------------------! 1643 SUBROUTINE bio_calculate_thermal_index_maps 1653 SUBROUTINE bio_calculate_thermal_index_maps( av ) 1644 1654 1645 1655 IMPLICIT NONE … … 1662 1672 1663 1673 ! 1664 !-- fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...) 1665 CALL bio_calculate_mrt_grid ( av ) 1666 1667 DO i = nxl, nxr 1668 DO j = nys, nyn 1669 ! 1670 !-- Determine local input conditions 1671 tmrt_ij = bio_fill_value 1672 vp = bio_fill_value 1673 ! 1674 !-- Determine input only if 1675 CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp, & 1676 ws, pair, tmrt_ij ) 1677 ! 1678 !-- Only proceed if input is available 1679 pet_ij = bio_fill_value !< set fail value, e.g. valid for within 1680 perct_ij = bio_fill_value !< some obstacle 1681 utci_ij = bio_fill_value 1682 IF ( .NOT. ( tmrt_ij <= -998._wp .OR. vp <= -998._wp ) ) THEN 1683 ! 1684 !-- Calculate static thermal indices based on local tmrt 1685 clo = bio_fill_value 1686 1687 IF ( do_calculate_perct .OR. do_calculate_perct_av ) THEN 1688 ! 1689 !-- Estimate local perceived temperature 1690 CALL calculate_perct_static( ta, vp, ws, tmrt_ij, pair, clo, & 1691 perct_ij ) 1674 !-- Check if some thermal index is desired. Don't do anything if, e.g. only 1675 !-- bio_mrt is desired. 1676 IF ( do_calculate_perct .OR. do_calculate_perct_av .OR. & 1677 do_calculate_utci .OR. do_calculate_utci_av .OR. & 1678 do_calculate_pet .OR. do_calculate_pet_av ) THEN 1679 1680 ! 1681 !-- fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...) 1682 CALL bio_calculate_mrt_grid ( av ) 1683 1684 DO i = nxl, nxr 1685 DO j = nys, nyn 1686 ! 1687 !-- Determine local input conditions 1688 tmrt_ij = bio_fill_value 1689 vp = bio_fill_value 1690 ! 1691 !-- Determine local meteorological conditions 1692 CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp, & 1693 ws, pair, tmrt_ij ) 1694 ! 1695 !-- Only proceed if input is available 1696 pet_ij = bio_fill_value !< set fail value, e.g. valid for 1697 perct_ij = bio_fill_value !< within some obstacle 1698 utci_ij = bio_fill_value 1699 IF ( .NOT. ( tmrt_ij <= -998._wp .OR. vp <= -998._wp .OR. & 1700 ws <= -998._wp .OR. ta <= -998._wp ) ) THEN 1701 ! 1702 !-- Calculate static thermal indices based on local tmrt 1703 clo = bio_fill_value 1704 1705 IF ( do_calculate_perct .OR. do_calculate_perct_av ) THEN 1706 ! 1707 !-- Estimate local perceived temperature 1708 CALL calculate_perct_static( ta, vp, ws, tmrt_ij, pair, & 1709 clo, perct_ij ) 1710 ENDIF 1711 1712 IF ( do_calculate_utci .OR. do_calculate_utci_av ) THEN 1713 ! 1714 !-- Estimate local universal thermal climate index 1715 CALL calculate_utci_static( ta, vp, ws, tmrt_ij, & 1716 bio_output_height, utci_ij ) 1717 ENDIF 1718 1719 IF ( do_calculate_pet .OR. do_calculate_pet_av ) THEN 1720 ! 1721 !-- Estimate local physiologically equivalent temperature 1722 CALL calculate_pet_static( ta, vp, ws, tmrt_ij, pair, & 1723 pet_ij ) 1724 ENDIF 1692 1725 ENDIF 1693 1726 1694 IF ( do_calculate_utci .OR. do_calculate_utci_av ) THEN 1695 ! 1696 !-- Estimate local universal thermal climate index 1697 CALL calculate_utci_static( ta, vp, ws, tmrt_ij, & 1698 bio_output_height, utci_ij ) 1727 1728 IF ( av ) THEN 1729 ! 1730 !-- Write results for selected averaged indices 1731 IF ( do_calculate_perct_av ) THEN 1732 perct_av(j, i) = perct_ij 1733 ENDIF 1734 IF ( do_calculate_utci_av ) THEN 1735 utci_av(j, i) = utci_ij 1736 ENDIF 1737 IF ( do_calculate_pet_av ) THEN 1738 pet_av(j, i) = pet_ij 1739 ENDIF 1740 ELSE 1741 ! 1742 !-- Write result for selected indices 1743 IF ( do_calculate_perct ) THEN 1744 perct(j, i) = perct_ij 1745 ENDIF 1746 IF ( do_calculate_utci ) THEN 1747 utci(j, i) = utci_ij 1748 ENDIF 1749 IF ( do_calculate_pet ) THEN 1750 pet(j, i) = pet_ij 1751 ENDIF 1699 1752 ENDIF 1700 1753 1701 IF ( do_calculate_pet .OR. do_calculate_pet_av ) THEN 1702 ! 1703 !-- Estimate local physiologically equivalent temperature 1704 CALL calculate_pet_static( ta, vp, ws, tmrt_ij, pair, pet_ij ) 1705 ENDIF 1706 ENDIF 1707 1708 1709 IF ( av ) THEN 1710 ! 1711 !-- Write results for selected averaged indices 1712 IF ( do_calculate_perct_av ) THEN 1713 perct_av(j, i) = perct_ij 1714 ENDIF 1715 IF ( do_calculate_utci_av ) THEN 1716 utci_av(j, i) = utci_ij 1717 ENDIF 1718 IF ( do_calculate_pet_av ) THEN 1719 pet_av(j, i) = pet_ij 1720 ENDIF 1721 ELSE 1722 ! 1723 !-- Write result for selected indices 1724 IF ( do_calculate_perct ) THEN 1725 perct(j, i) = perct_ij 1726 ENDIF 1727 IF ( do_calculate_utci ) THEN 1728 utci(j, i) = utci_ij 1729 ENDIF 1730 IF ( do_calculate_pet ) THEN 1731 pet(j, i) = pet_ij 1732 ENDIF 1733 ENDIF 1734 1754 ENDDO 1735 1755 ENDDO 1736 END DO1756 ENDIF 1737 1757 1738 1758 END SUBROUTINE bio_calculate_thermal_index_maps … … 1782 1802 !-- Initialize instationary perceived temperature with personalized 1783 1803 ! PT as an initial guess, set actlev and clo 1784 CALL ipt_init 1804 CALL ipt_init( age, weight, height, sex, work, actlev, clo, & 1785 1805 ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo, & 1786 1806 ipt ) … … 2294 2314 !-- Determine perceived temperature by regression equation + adjustments 2295 2315 pmvs = pmva 2296 CALL perct_regression 2316 CALL perct_regression( pmva, clo, perct_ij ) 2297 2317 ptc = perct_ij 2298 2318 IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN … … 2306 2326 pmvs = -0.11_wp 2307 2327 ENDIF 2308 CALL perct_regression 2328 CALL perct_regression( pmvs, clo, perct_ij ) 2309 2329 ENDIF 2310 2330 ! clo_fanger = clo … … 2318 2338 clo = clon 2319 2339 ENDIF 2320 CALL calc_sultr 2340 CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim ) 2321 2341 sultrieness = .FALSE. 2322 2342 d_std = -99._wp … … 2327 2347 d_pmv = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr ) 2328 2348 pmvs = pmva + d_pmv 2329 CALL perct_regression 2349 CALL perct_regression( pmvs, clo, perct_ij ) 2330 2350 IF ( sult_lim < 99._wp ) THEN 2331 2351 IF ( (perct_ij - ptc) > sult_lim ) sultrieness = .TRUE. 2332 2352 ! 2333 2353 !-- Set factor to threshold for sultriness 2334 IF ( dgtcstd /= 0_iwp ) THEN2354 IF ( ABS( dgtcstd ) > 0.00001_wp ) THEN 2335 2355 d_std = ( ( perct_ij - ptc ) - dgtcm ) / dgtcstd 2336 2356 ENDIF … … 2453 2473 y_average, top ) 2454 2474 sroot = SQRT ( y_average**2 - y_lower * y_upper ) 2455 IF ( sroot == 0._wp ) THEN2475 IF ( ABS( sroot ) < 0.00001_wp ) THEN 2456 2476 clo_res = x_average 2457 2477 nerr = j … … 2468 2488 CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta, & 2469 2489 y_new, top ) 2470 IF ( y_new == 0._wp ) THEN2490 IF ( ABS( y_new ) < 0.00001_wp ) THEN 2471 2491 clo_res = x_ridder 2472 2492 nerr = j 2473 2493 RETURN 2474 2494 ENDIF 2475 IF ( SIGN ( y_average, y_new ) /= y_average) THEN2495 IF ( ABS( SIGN( y_average, y_new ) - y_average ) > 0.00001_wp ) THEN 2476 2496 x_lower = x_average 2477 2497 y_lower = y_average 2478 2498 x_upper = x_ridder 2479 2499 y_upper = y_new 2480 ELSE IF ( SIGN ( y_lower, y_new ) /= y_lower) THEN2500 ELSE IF ( ABS( SIGN( y_lower, y_new ) - y_lower ) > 0.00001_wp ) THEN 2481 2501 x_upper = x_ridder 2482 2502 y_upper = y_new 2483 ELSE IF ( SIGN ( y_upper, y_new ) /= y_upper) THEN2503 ELSE IF ( ABS( SIGN( y_upper, y_new ) - y_upper ) > 0.00001_wp ) THEN 2484 2504 x_lower = x_ridder 2485 2505 y_lower = y_new … … 2487 2507 ! 2488 2508 !-- Never get here in x_ridder: singularity in y 2489 nerr 2509 nerr = -1_iwp 2490 2510 clo_res = x_ridder 2491 2511 RETURN … … 2493 2513 IF ( ABS ( x_upper - x_lower ) <= eps ) THEN 2494 2514 clo_res = x_ridder 2495 nerr 2515 nerr = j 2496 2516 RETURN 2497 2517 ENDIF … … 2502 2522 clo_res = y_new 2503 2523 RETURN 2504 ELSE IF ( y_lower == 0.) THEN2524 ELSE IF ( ABS( y_lower ) < 0.00001_wp ) THEN 2505 2525 x_ridder = clo_lower 2506 ELSE IF ( y_upper == 0.) THEN2526 ELSE IF ( ABS( y_upper ) < 0.00001_wp ) THEN 2507 2527 x_ridder = clo_upper 2508 2528 ELSE … … 2992 3012 ENDIF 2993 3013 2994 DO i = 1, 3 2995 delta_cold(i) = 0._wp 2996 IF ( i < 3 ) THEN 2997 pmv_cross(i) = pmvc 2998 ENDIF 2999 ENDDO 3014 delta_cold = 0._wp 3015 pmv_cross = pmvc 3000 3016 3001 3017 DO i = 1, 3 … … 3004 3020 reg_a(i) = coeff(i,1) + coeff(i,3) * ta + coeff(i,4) * & 3005 3021 sqrt_ws + coeff(i,5)*dtmrt 3006 delta_cold(i ) = reg_a(i) + coeff(i,2) * pmvc3022 delta_cold(i-1) = reg_a(i) + coeff(i,2) * pmvc 3007 3023 ENDDO 3008 3024 ! … … 3011 3027 r_nenner = coeff(i,2) - coeff(i+1,2) 3012 3028 IF ( ABS ( r_nenner ) > 0.00001_wp ) THEN 3013 pmv_cross(i ) = ( reg_a(i+1) - reg_a(i) ) / r_nenner3029 pmv_cross(i-1) = ( reg_a(i+1) - reg_a(i) ) / r_nenner 3014 3030 ELSE 3015 3031 nerr = 1_iwp … … 3020 3036 i_bin = 3 3021 3037 DO i = 1, 2 3022 IF ( pmva > pmv_cross(i ) ) THEN3038 IF ( pmva > pmv_cross(i-1) ) THEN 3023 3039 i_bin = i 3024 3040 EXIT … … 3028 3044 !-- Adjust to operative temperature scaled according 3029 3045 !-- to classical PMV (Fanger) 3030 dpmv_cold_res = delta_cold(i_bin ) - dpmv_adj(pmva)3046 dpmv_cold_res = delta_cold(i_bin-1) - dpmv_adj(pmva) 3031 3047 3032 3048 END SUBROUTINE dpmv_cold … … 3130 3146 !> according to its height (m), weight (kg), and age (y) 3131 3147 !------------------------------------------------------------------------------! 3132 SUBROUTINE surface_area 3148 SUBROUTINE surface_area( height_cm, weight, age, surf ) 3133 3149 3134 3150 IMPLICIT NONE … … 3174 3190 !> for a sample human. 3175 3191 !------------------------------------------------------------------------------! 3176 SUBROUTINE persdat 3192 SUBROUTINE persdat( age, weight, height, sex, work, a_surf, actlev ) 3177 3193 3178 3194 IMPLICIT NONE … … 3190 3206 REAL(wp) :: basic_heat_prod 3191 3207 3192 CALL surface_area 3208 CALL surface_area( height, weight, INT( age ), a_surf ) 3193 3209 s = height * 100._wp / ( weight**( 1._wp / 3._wp ) ) 3194 3210 factor = 1._wp + .004_wp * ( 30._wp - age ) … … 3215 3231 !------------------------------------------------------------------------------! 3216 3232 3217 SUBROUTINE ipt_init (age, weight, height, sex, work, actlev, clo, &3233 SUBROUTINE ipt_init( age, weight, height, sex, work, actlev, clo, & 3218 3234 ta, vp, ws, tmrt, pair, dt, storage, t_clothing, & 3219 3235 ipt ) … … 3270 3286 3271 3287 storage = 0._wp 3272 CALL persdat 3288 CALL persdat( age, weight, height, sex, work, a_surf, actlev ) 3273 3289 ! 3274 3290 !-- Initialise … … 3371 3387 !-- Determine perceived temperature by regression equation + adjustments 3372 3388 pmvs = pmva 3373 CALL perct_regression 3389 CALL perct_regression( pmva, clo, ipt ) 3374 3390 ptc = ipt 3375 3391 IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN … … 3383 3399 pmvs = -0.11_wp 3384 3400 ENDIF 3385 CALL perct_regression 3401 CALL perct_regression( pmvs, clo, ipt ) 3386 3402 ENDIF 3387 3403 ! clo_fanger = clo … … 3395 3411 clo = clon 3396 3412 ENDIF 3397 CALL calc_sultr 3413 CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim ) 3398 3414 sultrieness = .FALSE. 3399 3415 d_std = -99._wp … … 3404 3420 d_pmv = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr ) 3405 3421 pmvs = pmva + d_pmv 3406 CALL perct_regression 3422 CALL perct_regression( pmvs, clo, ipt ) 3407 3423 IF ( sult_lim < 99._wp ) THEN 3408 3424 IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE. … … 3476 3492 ! 3477 3493 !-- Determine perceived temperature by regression equation + adjustments 3478 CALL perct_regression 3494 CALL perct_regression( pmva, clo, ipt ) 3479 3495 3480 3496 IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN … … 3488 3504 pmvs = -0.11_wp 3489 3505 ENDIF 3490 CALL perct_regression 3506 CALL perct_regression( pmvs, clo, ipt ) 3491 3507 ENDIF 3492 3508 ! 3493 3509 !-- Consider sultriness if appropriate 3494 3510 ptc = ipt 3495 CALL calc_sultr 3511 CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim ) 3496 3512 sultrieness = .FALSE. 3497 3513 d_std = -99._wp … … 3502 3518 d_pmv = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr ) 3503 3519 pmvs = pmva + d_pmv 3504 CALL perct_regression 3520 CALL perct_regression( pmvs, clo, ipt ) 3505 3521 IF ( sult_lim < 99._wp ) THEN 3506 3522 IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE. … … 3722 3738 ! 3723 3739 !-- Call subfunctions 3724 CALL in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta,&3740 CALL in_body( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, & 3725 3741 vpa, work ) 3726 3742 3727 CALL heat_exch ( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, ht,&3743 CALL heat_exch( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, ht, & 3728 3744 int_heat,mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa, & 3729 3745 vpts, wetsk ) 3730 3746 3731 CALL pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair,&3747 CALL pet_iteration( acl, adu, aeff, esw, facl, feff, int_heat, pair, & 3732 3748 rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk ) 3733 3749 … … 3741 3757 !> Calculate internal energy ballance 3742 3758 !------------------------------------------------------------------------------! 3743 SUBROUTINE in_body 3759 SUBROUTINE in_body( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, & 3744 3760 vpa, work ) 3745 3761 ! … … 3794 3810 !> Calculate heat gain or loss 3795 3811 !------------------------------------------------------------------------------! 3796 SUBROUTINE heat_exch ( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff,&3812 SUBROUTINE heat_exch( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, & 3797 3813 ht, int_heat, mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa, & 3798 3814 vpts, wetsk ) … … 3953 3969 ( c(5) * tsk - c(0) - 5.28_wp * adu * tsk ) 3954 3970 3955 IF ( tsk == 36._wp ) tsk = 36.01_wp3971 IF ( ABS( tsk - 36._wp ) < 0.00001_wp ) tsk = 36.01_wp 3956 3972 tcore(7) = c(0) / ( 5.28_wp * adu + c(1) * 6.3_wp / 3600._wp ) + tsk 3957 3973 tcore(3) = c(0) / ( 5.28_wp * adu + ( c(1) * 6.3_wp / 3600._wp ) / & … … 4077 4093 !> Calculate PET 4078 4094 !------------------------------------------------------------------------------! 4079 SUBROUTINE pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair,&4095 SUBROUTINE pet_iteration( acl, adu, aeff, esw, facl, feff, int_heat, pair, & 4080 4096 rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk ) 4081 4097 !
Note: See TracChangeset
for help on using the changeset viewer.
