Ignore:
Timestamp:
Feb 19, 2019 7:18:38 AM (2 years ago)
Author:
dom_dwd_user
Message:

biometeorology_mod.f90:
(B) Added addittional safety meassures to bio_calculate_thermal_index_maps.
(B) Replaced several REAL (un-)equality comparisons.

File:
1 edited

Legend:

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

    r3742 r3749  
    2727! -----------------
    2828! $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
    2933! - Allocation of the input _av grids was moved to the "sum" section of
    3034! bio_3d_data_averaging to make sure averaging is only done once!
     
    523527             ENDIF
    524528
     529!
     530!-- u_av, v_av and w_av are always allocated
    525531             IF ( .NOT. ALLOCATED( u_av ) )  THEN
    526532                ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     
    866872 SUBROUTINE bio_check_parameters
    867873
    868     USE control_parameters,                                                    &
    869         ONLY:  message_string
    870874
    871875    IMPLICIT NONE
     
    14551459 END SUBROUTINE bio_wrd_local
    14561460
    1457 
    14581461!------------------------------------------------------------------------------!
    14591462! Description:
     
    15841587!
    15851588!--    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
    15871593
    15881594       vp = bio_fill_value
     
    15911597       ENDIF
    15921598
    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
    15961606    ELSE
    15971607!
     
    16041614       ENDIF
    16051615
    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) ) )
    16091619
    16101620    ENDIF
     
    16411651!> time_integration.f90: 1065ff
    16421652!------------------------------------------------------------------------------!
    1643  SUBROUTINE bio_calculate_thermal_index_maps ( av )
     1653 SUBROUTINE bio_calculate_thermal_index_maps( av )
    16441654
    16451655    IMPLICIT NONE
     
    16621672
    16631673!
    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
    16921725             ENDIF
    16931726
    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
    16991752             ENDIF
    17001753
    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
    17351755       ENDDO
    1736     ENDDO
     1756    ENDIF
    17371757
    17381758 END SUBROUTINE bio_calculate_thermal_index_maps
     
    17821802!--    Initialize instationary perceived temperature with personalized
    17831803!      PT as an initial guess, set actlev and clo
    1784        CALL ipt_init ( age, weight, height, sex, work, actlev, clo,            &
     1804       CALL ipt_init( age, weight, height, sex, work, actlev, clo,            &
    17851805          ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo,                   &
    17861806          ipt )
     
    22942314!-- Determine perceived temperature by regression equation + adjustments
    22952315    pmvs = pmva
    2296     CALL perct_regression ( pmva, clo, perct_ij )
     2316    CALL perct_regression( pmva, clo, perct_ij )
    22972317    ptc = perct_ij
    22982318    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
     
    23062326          pmvs   = -0.11_wp
    23072327       ENDIF
    2308        CALL perct_regression ( pmvs, clo, perct_ij )
     2328       CALL perct_regression( pmvs, clo, perct_ij )
    23092329    ENDIF
    23102330!     clo_fanger = clo
     
    23182338       clo = clon
    23192339    ENDIF
    2320     CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
     2340    CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
    23212341    sultrieness    = .FALSE.
    23222342    d_std = -99._wp
     
    23272347       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
    23282348       pmvs   = pmva + d_pmv
    2329        CALL perct_regression ( pmvs, clo, perct_ij )
     2349       CALL perct_regression( pmvs, clo, perct_ij )
    23302350       IF ( sult_lim < 99._wp )  THEN
    23312351          IF ( (perct_ij - ptc) > sult_lim )  sultrieness = .TRUE.
    23322352!
    23332353!--       Set factor to threshold for sultriness
    2334           IF ( dgtcstd /= 0_iwp )  THEN
     2354          IF ( ABS( dgtcstd ) > 0.00001_wp )  THEN
    23352355             d_std = ( ( perct_ij - ptc ) - dgtcm ) / dgtcstd
    23362356          ENDIF
     
    24532473                        y_average, top )
    24542474          sroot = SQRT ( y_average**2 - y_lower * y_upper )
    2455           IF ( sroot == 0._wp )  THEN
     2475          IF ( ABS( sroot ) < 0.00001_wp )  THEN
    24562476             clo_res = x_average
    24572477             nerr = j
     
    24682488          CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta,         &
    24692489                        y_new, top )
    2470           IF ( y_new == 0._wp )  THEN
     2490          IF ( ABS( y_new ) < 0.00001_wp )  THEN
    24712491             clo_res = x_ridder
    24722492             nerr       = j
    24732493             RETURN
    24742494          ENDIF
    2475           IF ( SIGN ( y_average, y_new ) /= y_average )  THEN
     2495          IF ( ABS( SIGN( y_average, y_new ) - y_average ) > 0.00001_wp )  THEN
    24762496             x_lower = x_average
    24772497             y_lower = y_average
    24782498             x_upper  = x_ridder
    24792499             y_upper  = y_new
    2480           ELSE IF ( SIGN ( y_lower, y_new ) /= y_lower )  THEN
     2500          ELSE IF ( ABS( SIGN( y_lower, y_new ) - y_lower ) > 0.00001_wp )  THEN
    24812501             x_upper  = x_ridder
    24822502             y_upper  = y_new
    2483           ELSE IF ( SIGN ( y_upper, y_new ) /= y_upper )  THEN
     2503          ELSE IF ( ABS( SIGN( y_upper, y_new ) - y_upper ) > 0.00001_wp )  THEN
    24842504             x_lower = x_ridder
    24852505             y_lower = y_new
     
    24872507!
    24882508!--          Never get here in x_ridder: singularity in y
    2489              nerr       = -1_iwp
     2509             nerr    = -1_iwp
    24902510             clo_res = x_ridder
    24912511             RETURN
     
    24932513          IF ( ABS ( x_upper - x_lower ) <= eps )  THEN
    24942514             clo_res = x_ridder
    2495              nerr       = j
     2515             nerr    = j
    24962516             RETURN
    24972517          ENDIF
     
    25022522       clo_res = y_new
    25032523       RETURN
    2504     ELSE IF ( y_lower == 0. )  THEN
     2524    ELSE IF ( ABS( y_lower ) < 0.00001_wp )  THEN
    25052525       x_ridder = clo_lower
    2506     ELSE IF ( y_upper == 0. )  THEN
     2526    ELSE IF ( ABS( y_upper ) < 0.00001_wp )  THEN
    25072527       x_ridder = clo_upper
    25082528    ELSE
     
    29923012    ENDIF
    29933013
    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
    30003016
    30013017    DO  i = 1, 3
     
    30043020       reg_a(i) = coeff(i,1) + coeff(i,3) * ta + coeff(i,4) *                  &
    30053021                  sqrt_ws + coeff(i,5)*dtmrt
    3006        delta_cold(i) = reg_a(i) + coeff(i,2) * pmvc
     3022       delta_cold(i-1) = reg_a(i) + coeff(i,2) * pmvc
    30073023    ENDDO
    30083024!
     
    30113027       r_nenner = coeff(i,2) - coeff(i+1,2)
    30123028       IF ( ABS ( r_nenner ) > 0.00001_wp )  THEN
    3013           pmv_cross(i) = ( reg_a(i+1) - reg_a(i) ) / r_nenner
     3029          pmv_cross(i-1) = ( reg_a(i+1) - reg_a(i) ) / r_nenner
    30143030       ELSE
    30153031          nerr = 1_iwp
     
    30203036    i_bin = 3
    30213037    DO  i = 1, 2
    3022        IF ( pmva > pmv_cross(i) )  THEN
     3038       IF ( pmva > pmv_cross(i-1) )  THEN
    30233039          i_bin = i
    30243040          EXIT
     
    30283044!-- Adjust to operative temperature scaled according
    30293045!-- 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)
    30313047
    30323048 END SUBROUTINE dpmv_cold
     
    31303146!> according to its height (m), weight (kg), and age (y)
    31313147!------------------------------------------------------------------------------!
    3132  SUBROUTINE surface_area ( height_cm, weight, age, surf )
     3148 SUBROUTINE surface_area( height_cm, weight, age, surf )
    31333149
    31343150    IMPLICIT NONE
     
    31743190!> for a sample human.
    31753191!------------------------------------------------------------------------------!
    3176  SUBROUTINE persdat ( age, weight, height, sex, work, a_surf, actlev )
     3192 SUBROUTINE persdat( age, weight, height, sex, work, a_surf, actlev )
    31773193
    31783194    IMPLICIT NONE
     
    31903206    REAL(wp) ::  basic_heat_prod
    31913207
    3192     CALL surface_area ( height, weight, INT( age ), a_surf )
     3208    CALL surface_area( height, weight, INT( age ), a_surf )
    31933209    s = height * 100._wp / ( weight**( 1._wp / 3._wp ) )
    31943210    factor = 1._wp + .004_wp  * ( 30._wp - age )
     
    32153231!------------------------------------------------------------------------------!
    32163232
    3217  SUBROUTINE ipt_init (age, weight, height, sex, work, actlev, clo,             &
     3233 SUBROUTINE ipt_init( age, weight, height, sex, work, actlev, clo,             &
    32183234     ta, vp, ws, tmrt, pair, dt, storage, t_clothing,                          &
    32193235     ipt )
     
    32703286
    32713287    storage = 0._wp
    3272     CALL persdat ( age, weight, height, sex, work, a_surf, actlev )
     3288    CALL persdat( age, weight, height, sex, work, a_surf, actlev )
    32733289!
    32743290!-- Initialise
     
    33713387!-- Determine perceived temperature by regression equation + adjustments
    33723388    pmvs = pmva
    3373     CALL perct_regression ( pmva, clo, ipt )
     3389    CALL perct_regression( pmva, clo, ipt )
    33743390    ptc = ipt
    33753391    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
     
    33833399          pmvs   = -0.11_wp
    33843400       ENDIF
    3385        CALL perct_regression ( pmvs, clo, ipt )
     3401       CALL perct_regression( pmvs, clo, ipt )
    33863402    ENDIF
    33873403!     clo_fanger = clo
     
    33953411       clo = clon
    33963412    ENDIF
    3397     CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
     3413    CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
    33983414    sultrieness    = .FALSE.
    33993415    d_std      = -99._wp
     
    34043420       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
    34053421       pmvs   = pmva + d_pmv
    3406        CALL perct_regression ( pmvs, clo, ipt )
     3422       CALL perct_regression( pmvs, clo, ipt )
    34073423       IF ( sult_lim < 99._wp )  THEN
    34083424          IF ( (ipt - ptc) > sult_lim )  sultrieness = .TRUE.
     
    34763492!
    34773493!-- Determine perceived temperature by regression equation + adjustments
    3478     CALL perct_regression ( pmva, clo, ipt )
     3494    CALL perct_regression( pmva, clo, ipt )
    34793495
    34803496    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
     
    34883504          pmvs   = -0.11_wp
    34893505       ENDIF
    3490        CALL perct_regression ( pmvs, clo, ipt )
     3506       CALL perct_regression( pmvs, clo, ipt )
    34913507    ENDIF
    34923508!
    34933509!-- Consider sultriness if appropriate
    34943510    ptc = ipt
    3495     CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
     3511    CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
    34963512    sultrieness    = .FALSE.
    34973513    d_std      = -99._wp
     
    35023518       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
    35033519       pmvs   = pmva + d_pmv
    3504        CALL perct_regression ( pmvs, clo, ipt )
     3520       CALL perct_regression( pmvs, clo, ipt )
    35053521       IF ( sult_lim < 99._wp )  THEN
    35063522          IF ( (ipt - ptc) > sult_lim )  sultrieness = .TRUE.
     
    37223738!
    37233739!-- 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,     &
    37253741            vpa, work )
    37263742
    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, &
    37283744            int_heat,mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa,      &
    37293745            vpts, wetsk )
    37303746
    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,       &
    37323748            rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk )
    37333749
     
    37413757!> Calculate internal energy ballance
    37423758!------------------------------------------------------------------------------!
    3743  SUBROUTINE in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, &
     3759 SUBROUTINE in_body( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, &
    37443760    vpa, work )
    37453761!
     
    37943810!> Calculate heat gain or loss
    37953811!------------------------------------------------------------------------------!
    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,   &
    37973813        ht, int_heat, mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa,     &
    37983814        vpts, wetsk )
     
    39533969                  ( c(5) * tsk - c(0) - 5.28_wp * adu * tsk )
    39543970
    3955           IF ( tsk == 36._wp )  tsk = 36.01_wp
     3971          IF ( ABS( tsk - 36._wp ) < 0.00001_wp )  tsk = 36.01_wp
    39563972          tcore(7) = c(0) / ( 5.28_wp * adu + c(1) * 6.3_wp / 3600._wp ) + tsk
    39573973          tcore(3) = c(0) / ( 5.28_wp * adu + ( c(1) * 6.3_wp / 3600._wp ) /   &
     
    40774093!> Calculate PET
    40784094!------------------------------------------------------------------------------!
    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,    &
    40804096        rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk )
    40814097!
Note: See TracChangeset for help on using the changeset viewer.