Changeset 3753

Timestamp:

Feb 19, 2019 2:48:54 PM (6 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:

• palm/trunk/SOURCE/biometeorology_mod.f90 (modified) (62 diffs)

palm/trunk/SOURCE/biometeorology_mod.f90

@@ -27 +27 @@
 ! -----------------
 ! $Id$
+! - Added automatic setting of mrt_nlevels in case it was not part of 
+! radiation_parameters namelist (or set to 0 accidentially).
+! - Minor speed improvoemnts in perceived temperature calculations.
+! - Perceived temperature regression arrays now declared as PARAMETERs.
+! 
+! 3750 2019-02-19 07:29:39Z dom_dwd_user
 ! - Added addittional safety meassures to bio_calculate_thermal_index_maps.
 ! - Replaced several REAL (un-)equality comparisons.
@@ -836 +842 @@
 !--    Further checks if thermal comfort output is desired.
     IF ( thermal_comfort  .AND.  unit == 'degree_C' )  THEN
-
-!
-!--    Break if required modules "radiation" and "humidity" are not avalable.
+!
+!--    Break if required modules "radiation" is not avalable.
        IF ( .NOT.  radiation )  THEN
           message_string = 'output of "' // TRIM( var ) // '" require'         &
@@ -845 +850 @@
           unit = 'illegal'
        ENDIF
-       IF ( .NOT.  humidity )  THEN
-          message_string = 'The estimation of thermal comfort requires '    // &
-                           'air humidity information, but humidity module ' // &
-                           'is disabled!'
-          CALL message( 'check_parameters', 'PA0561', 1, 2, 0, 6, 0 )
-          unit = 'illegal'
-       ENDIF
-       IF ( mrt_nlevels == 0 )  THEN
-          message_string = 'output of "' // TRIM( var ) // '" require'         &
-                           // 's mrt_nlevels > 0'
-          CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
-          unit = 'illegal'
+!
+!--    All "thermal_comfort" outputs except from 'bio_mrt' will also need 
+!--    humidity input. Check also for that.
+       IF ( TRIM( var ) /= 'bio_mrt' )  THEN
+          IF ( .NOT.  humidity )  THEN
+             message_string = 'The estimation of thermal comfort '    //       &
+                              'requires air humidity information, but ' //     &
+                              'humidity module is disabled!'
+             CALL message( 'check_parameters', 'PA0561', 1, 2, 0, 6, 0 )
+             unit = 'illegal'
+          ENDIF
        ENDIF
 
@@ -1214 +1218 @@
     height = 0.0_wp
 
-    bio_cell_level = INT ( 1.099_wp / dz(1) )
+    bio_cell_level = INT( 1.099_wp / dz(1) )
     bio_output_height = bio_output_height + bio_cell_level * dz(1)
-
+!
+!-- Set radiation level if not done by user
+    IF ( mrt_nlevels == 0 )  THEN
+       mrt_nlevels = bio_cell_level + 1_iwp
+    ENDIF
 !
 !-- Init UVEM and load lookup tables
@@ -1574 +1582 @@
 !-- Determine cell level closest to 1.1m above ground
 !   by making use of truncation due to int cast
-    k = get_topography_top_index_ji(j, i, 's') + bio_cell_level  !< Vertical cell center closest to 1.1m
+    k = INT( get_topography_top_index_ji(j, i, 's') + bio_cell_level )  !< Vertical cell center closest to 1.1m
 
 !
@@ -1589 +1597 @@
        ta = bio_fill_value
        IF ( ALLOCATED( pt_av ) )  THEN
-          ta = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
+          ta = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * ( k + .5_wp ) ) - degc_to_k
        ENDIF
 
@@ -2209 +2217 @@
     REAL(wp) ::  clon           !< clo for neutral conditions   (clo)
     REAL(wp) ::  ireq_minimal   !< Minimal required clothing insulation (clo)
-!     REAL(wp) ::  clo_fanger     !< clo for fanger subroutine, unused
     REAL(wp) ::  pmv_w          !< Fangers predicted mean vote for winter clothing
     REAL(wp) ::  pmv_s          !< Fangers predicted mean vote for summer clothing
@@ -2216 +2223 @@
     REAL(wp) ::  d_std          !< factor to threshold for sultriness
     REAL(wp) ::  pmvs           !< pred. mean vote considering sultrieness
-    REAL(wp) ::  top            !< Gagge's operative temperatures (degree_C)
 
     INTEGER(iwp) :: ncount      !< running index
@@ -2244 +2250 @@
 !--    First guess: winter clothing insulation: cold stress
        clo = wclo
-       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, top )
+       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
        pmv_w = pmva
 
@@ -2251 +2257 @@
 !--       Case summer clothing insulation: heat load ?
           clo = sclo
-          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
-             top )
+          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
           pmv_s = pmva
           IF ( pmva <= 0._wp )  THEN
@@ -2259 +2264 @@
 !--          Between winter and summer set values
              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
-                pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
+                pmv_s, wclo, pmv_w, eps, pmva, ncount, clo )
              IF ( ncount < 0_iwp )  THEN
                 nerr = -1_iwp
@@ -2267 +2272 @@
              clo = 0.5_wp
              CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
-                           pmva, top )
+                           pmva )
           ENDIF
        ELSE IF ( pmva < -0.11_wp )  THEN
           clo = 1.75_wp
-          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
-             top )
+          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
        ENDIF
     ELSE
@@ -2278 +2282 @@
 !--    First guess: summer clothing insulation: heat load
        clo = sclo
-       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, top )
+       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
        pmv_s = pmva
 
@@ -2285 +2289 @@
 !--       Case winter clothing insulation: cold stress ?
           clo = wclo
-          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
-             top )
+          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
           pmv_w = pmva
 
@@ -2294 +2297 @@
 !--          between winter and summer set values
              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
-                               pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
+                               pmv_s, wclo, pmv_w, eps, pmva, ncount, clo )
              IF ( ncount < 0_iwp )  THEN
                 nerr = -1_iwp
@@ -2302 +2305 @@
              clo = 1.75_wp
              CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
-                           pmva, top )
+                           pmva )
           ENDIF
        ELSE IF ( pmva > 0.06_wp )  THEN
           clo = 0.5_wp
-          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
-             top )
+          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
        ENDIF
 
@@ -2333 +2335 @@
 !
 !--    Required clothing insulation (ireq) is exclusively defined for 
-!--    operative temperatures (top) less 10 (C) for a 
+!--    perceived temperatures (perct) less 10 (C) for a 
 !--    reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
        clon = ireq_neutral ( perct_ij, ireq_minimal, nerr )
@@ -2393 +2395 @@
 !
 !-- Saturation water vapour pressure
-    svp_ta = 6.1078_wp * EXP ( b * ta / ( c + ta ) )
+    svp_ta = 6.1078_wp * EXP( b * ta / ( c + ta ) )
 
  END SUBROUTINE saturation_vapor_pressure
@@ -2405 +2407 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE iso_ridder( ta, tmrt, vp, ws, pair, actlev, eta, sclo,             &
-                       pmv_s, wclo, pmv_w, eps, pmva, top, nerr,               &
+                       pmv_s, wclo, pmv_w, eps, pmva, nerr,               &
                        clo_res )
 
@@ -2427 +2429 @@
 !-- Output variables of argument list:
     REAL(wp), INTENT ( OUT ) :: pmva     !< 0 (set to zero, because clo is evaluated for comfort)
-    REAL(wp), INTENT ( OUT ) :: top      !< Operative temperature (degC) at found root of Fanger's PMV
     REAL(wp), INTENT ( OUT ) :: clo_res  !< Resulting clothing insulation value (clo)
     INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error status / quality flag
@@ -2471 +2472 @@
           x_average = 0.5_wp * ( x_lower + x_upper )
           CALL fanger ( ta, tmrt, vp, ws, pair, x_average, actlev, eta,        &
-                        y_average, top )
-          sroot = SQRT ( y_average**2 - y_lower * y_upper )
+                        y_average )
+          sroot = SQRT( y_average**2 - y_lower * y_upper )
           IF ( ABS( sroot ) < 0.00001_wp )  THEN
              clo_res = x_average
@@ -2480 +2481 @@
           x_new = x_average + ( x_average - x_lower ) *                        &
                       ( SIGN ( 1._wp, y_lower - y_upper ) * y_average / sroot )
-          IF ( ABS ( x_new - x_ridder ) <= eps )  THEN
+          IF ( ABS( x_new - x_ridder ) <= eps )  THEN
              clo_res = x_ridder
              nerr       = j
@@ -2487 +2488 @@
           x_ridder = x_new
           CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta,         &
-                        y_new, top )
+                        y_new )
           IF ( ABS( y_new ) < 0.00001_wp )  THEN
              clo_res = x_ridder
@@ -2511 +2512 @@
              RETURN
           ENDIF
-          IF ( ABS ( x_upper - x_lower ) <= eps )  THEN
+          IF ( ABS( x_upper - x_lower ) <= eps )  THEN
              clo_res = x_ridder
              nerr    = j
@@ -2574 +2575 @@
 !> The case of free convection (ws < 0.1 m/s) is dealt with ws = 0.1 m/s
 !------------------------------------------------------------------------------!
- SUBROUTINE fanger( ta, tmrt, pa, in_ws, pair, in_clo, actlev, eta, pmva, top )
+ SUBROUTINE fanger( ta, tmrt, pa, in_ws, pair, in_clo, actlev, eta, pmva )
 
     IMPLICIT NONE
@@ -2592 +2593 @@
                                          !< dimensionless) according to Fanger corresponding to meteorological 
                                          !< (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta)
-    REAL(wp), INTENT ( OUT ) ::  top     !< operative temperature (degC)
 !
 !-- Internal variables
@@ -2631 +2631 @@
 !
 !-- Heat_convection = forced convection
-    heat_convection = 12.1_wp * SQRT ( ws * pair / 1013.25_wp )
+    heat_convection = 12.1_wp * SQRT( ws * pair / 1013.25_wp )
 !
 !-- Activity = inner heat produktion per standardized surface
@@ -2656 +2656 @@
 !-- Empiric factor for the adaption of the heat ballance equation
 !-- to the psycho-physical scale (Equ. 40 in FANGER)
-    z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp )
+    z1 = ( .303_wp * EXP( -.036_wp * actlev ) + .0275_wp )
 !
 !-- Water vapour diffution through the skin
@@ -2671 +2671 @@
     z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + degc_to_k )**4 - ( tmrt +        &
        degc_to_k )**4 )
-    IF ( ABS ( t_clothing - tmrt ) > 0._wp )  THEN
+    IF ( ABS( t_clothing - tmrt ) > 0._wp )  THEN
        hr = z5 / f_cl / ( t_clothing - tmrt )
     ELSE
@@ -2682 +2682 @@
 !-- Predicted Mean Vote
     pmva = z1 * ( activity - z2 - z3 - z4 - z5 - z6 )
-!
-!-- Operative temperatur
-    top = ( hr * tmrt + heat_convection * ta ) / ( hr + heat_convection )
 
  END SUBROUTINE fanger
@@ -2731 +2728 @@
     REAL(wp) ::  dpmv_2       !< 
     REAL(wp) ::  pmvs         !< 
-    REAL(wp) ::  bpmv(0:7)    !< 
-    REAL(wp) ::  bpa_p50(0:7) !< 
-    REAL(wp) ::  bpa(0:7)     !< 
-    REAL(wp) ::  bapa(0:7)    !< 
-    REAL(wp) ::  bdapa(0:7)   !< 
-    REAL(wp) ::  bsqvel(0:7)  !< 
-    REAL(wp) ::  bta(0:7)     !< 
-    REAL(wp) ::  bdtmrt(0:7)  !< 
-    REAL(wp) ::  aconst(0:7)  !< 
     INTEGER(iwp) :: nreg      !< 
 
-    DATA bpmv     /                                                            &
-     -0.0556602_wp, -0.1528680_wp, -0.2336104_wp, -0.2789387_wp, -0.3551048_wp,&
-     -0.4304076_wp, -0.4884961_wp, -0.4897495_wp /
-    DATA bpa_p50 /                                                             &
-     -0.1607154_wp, -0.4177296_wp, -0.4120541_wp, -0.0886564_wp, +0.4285938_wp,&
-     +0.6281256_wp, +0.5067361_wp, +0.3965169_wp /
-    DATA bpa     /                                                             &
-     +0.0580284_wp, +0.0836264_wp, +0.1009919_wp, +0.1020777_wp, +0.0898681_wp,&
-     +0.0839116_wp, +0.0853258_wp, +0.0866589_wp /
-    DATA bapa    /                                                             &
-     -1.7838788_wp, -2.9306231_wp, -1.6350334_wp, +0.6211547_wp, +3.3918083_wp,&
-     +5.5521025_wp, +8.4897418_wp, +16.6265851_wp /
-    DATA bdapa   /                                                             &
-     +1.6752720_wp, +2.7379504_wp, +1.2940526_wp, -1.0985759_wp, -3.9054732_wp,&
-     -6.0403012_wp, -8.9437119_wp, -17.0671201_wp /
-    DATA bsqvel  /                                                             &
-     -0.0315598_wp, -0.0286272_wp, -0.0009228_wp, +0.0483344_wp, +0.0992366_wp,&
-     +0.1491379_wp, +0.1951452_wp, +0.2133949_wp /
-    DATA bta     /                                                             &
-     +0.0953986_wp, +0.1524760_wp, +0.0564241_wp, -0.0893253_wp, -0.2398868_wp,&
-     -0.3515237_wp, -0.5095144_wp, -0.9469258_wp /
-    DATA bdtmrt  /                                                             &
-     -0.0004672_wp, -0.0000514_wp, -0.0018037_wp, -0.0049440_wp, -0.0069036_wp,&
-     -0.0075844_wp, -0.0079602_wp, -0.0089439_wp /
-    DATA aconst  /                                                             &
-     +1.8686215_wp, +3.4260713_wp, +2.0116185_wp, -0.7777552_wp, -4.6715853_wp,&
-     -7.7314281_wp, -11.7602578_wp, -23.5934198_wp /
+!
+!-- Regression coefficients:
+    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpmv = (/                          &
+     -0.0556602_wp, -0.1528680_wp, -0.2336104_wp, -0.2789387_wp,               &
+     -0.3551048_wp, -0.4304076_wp, -0.4884961_wp, -0.4897495_wp /)
+
+    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpa_p50 = (/                       &
+     -0.1607154_wp, -0.4177296_wp, -0.4120541_wp, -0.0886564_wp,               &
+      0.4285938_wp,  0.6281256_wp,  0.5067361_wp,  0.3965169_wp /)
+
+    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpa = (/                           &
+      0.0580284_wp,  0.0836264_wp,  0.1009919_wp,  0.1020777_wp,               &
+      0.0898681_wp,  0.0839116_wp,  0.0853258_wp,  0.0866589_wp /)
+
+    REAL(wp), DIMENSION(0:7), PARAMETER ::  bapa = (/                          &
+     -1.7838788_wp, -2.9306231_wp, -1.6350334_wp,   0.6211547_wp,              &
+      3.3918083_wp,  5.5521025_wp,  8.4897418_wp,  16.6265851_wp /)
+
+    REAL(wp), DIMENSION(0:7), PARAMETER ::  bdapa = (/                         &
+      1.6752720_wp,  2.7379504_wp,  1.2940526_wp,  -1.0985759_wp,              &
+     -3.9054732_wp, -6.0403012_wp, -8.9437119_wp, -17.0671201_wp /)
+
+    REAL(wp), DIMENSION(0:7), PARAMETER ::  bsqvel = (/                        &
+     -0.0315598_wp, -0.0286272_wp, -0.0009228_wp,  0.0483344_wp,               &
+      0.0992366_wp,  0.1491379_wp,  0.1951452_wp,  0.2133949_wp /)
+
+    REAL(wp), DIMENSION(0:7), PARAMETER ::  bta = (/                           &
+      0.0953986_wp,  0.1524760_wp,  0.0564241_wp, -0.0893253_wp,               &
+     -0.2398868_wp, -0.3515237_wp, -0.5095144_wp, -0.9469258_wp /)
+
+    REAL(wp), DIMENSION(0:7), PARAMETER ::  bdtmrt = (/                        &
+     -0.0004672_wp, -0.0000514_wp, -0.0018037_wp, -0.0049440_wp,               &
+     -0.0069036_wp, -0.0075844_wp, -0.0079602_wp, -0.0089439_wp /)
+
+    REAL(wp), DIMENSION(0:7), PARAMETER ::  aconst = (/                        &
+      1.8686215_wp,  3.4260713_wp,   2.0116185_wp,  -0.7777552_wp,             &
+     -4.6715853_wp, -7.7314281_wp, -11.7602578_wp, -23.5934198_wp /)
+
+
 !
 !-- Test for compliance with regression range
@@ -2800 +2800 @@
 !
 !--    Natural logarithm of pa
-       apa = LOG ( pa )
+       apa = LOG( pa )
     ELSE
        apa = -5._wp
@@ -2808 +2808 @@
     pa_p50   = 0.5_wp * svp_ta
     IF ( pa_p50 > 0._wp  .AND.  pa > 0._wp )  THEN
-       dapa   = apa - LOG ( pa_p50 )
+       dapa   = apa - LOG( pa_p50 )
        pa_p50 = pa / pa_p50
     ELSE
@@ -2817 +2817 @@
 !-- Square root of wind velocity
     IF ( ws >= 0._wp )  THEN
-       sqvel = SQRT ( ws )
+       sqvel = SQRT( ws )
     ELSE
        sqvel = 0._wp
@@ -2826 +2826 @@
 !
 !-- Select the valid regression coefficients
-    nreg = INT ( pmv )
+    nreg = INT( pmv )
     IF ( nreg < 0_iwp )  THEN
 !
@@ -2836 +2836 @@
     IF ( weight < 0._wp )  weight = 0._wp
     IF ( nreg > 5_iwp )  THEN
-       ! nreg=6
        nreg  = 5_iwp
        weight   = pmv - 5._wp
@@ -2845 +2844 @@
     ENDIF
 !
-!-- Regression valid for 0. <= pmv <= 6.
+!-- Regression valid for 0. <= pmv <= 6., bounds are checked above
     dpmv_1 =                                                                   &
        + bpa(nreg) * pa                                                        &
@@ -2857 +2856 @@
        + aconst(nreg)
 
-    dpmv_2 = 0._wp
-    IF ( nreg < 6_iwp )  THEN
-       dpmv_2 =                                                                &
-          + bpa(nreg+1_iwp)     * pa                                           &
-          + bpmv(nreg+1_iwp)    * pmv                                          &
-          + bapa(nreg+1_iwp)    * apa                                          &
-          + bta(nreg+1_iwp)     * ta                                           &
-          + bdtmrt(nreg+1_iwp)  * dtmrt                                        &
-          + bdapa(nreg+1_iwp)   * dapa                                         &
-          + bsqvel(nreg+1_iwp)  * sqvel                                        &
-          + bpa_p50(nreg+1_iwp) * pa_p50                                       &
-          + aconst(nreg+1_iwp)
-    ENDIF
+!    dpmv_2 = 0._wp
+!    IF ( nreg < 6_iwp )  THEN  !< nreg is always <= 5, see above
+    dpmv_2 =                                                                   &
+       + bpa(nreg+1_iwp)     * pa                                              &
+       + bpmv(nreg+1_iwp)    * pmv                                             &
+       + bapa(nreg+1_iwp)    * apa                                             &
+       + bta(nreg+1_iwp)     * ta                                              &
+       + bdtmrt(nreg+1_iwp)  * dtmrt                                           &
+       + bdapa(nreg+1_iwp)   * dapa                                            &
+       + bsqvel(nreg+1_iwp)  * sqvel                                           &
+       + bpa_p50(nreg+1_iwp) * pa_p50                                          &
+       + aconst(nreg+1_iwp)
+!    ENDIF
 !
 !-- Calculate pmv modification
@@ -2880 +2879 @@
        IF ( pmvs > -0.11_wp )  THEN
 !
-!--       Threshold from FUNCTION perct_regression for winter clothing insulation
+!--       Threshold from perct_regression for winter clothing insulation
           deltapmv = deltapmv + 0.11_wp
        ELSE
@@ -2915 +2914 @@
 !
 !-- Types of coefficients mean deviation: third order polynomial
-    REAL(wp), PARAMETER ::  dperctka = +7.5776086_wp
+    REAL(wp), PARAMETER ::  dperctka =  7.5776086_wp
     REAL(wp), PARAMETER ::  dperctkb = -0.740603_wp
-    REAL(wp), PARAMETER ::  dperctkc = +0.0213324_wp
+    REAL(wp), PARAMETER ::  dperctkc =  0.0213324_wp
     REAL(wp), PARAMETER ::  dperctkd = -0.00027797237_wp
 !
 !-- Types of coefficients mean deviation plus standard deviation
 !-- regression coefficients: third order polynomial
-    REAL(wp), PARAMETER ::  dperctsa = +0.0268918_wp
-    REAL(wp), PARAMETER ::  dperctsb = +0.0465957_wp
+    REAL(wp), PARAMETER ::  dperctsa =  0.0268918_wp
+    REAL(wp), PARAMETER ::  dperctsb =  0.0465957_wp
     REAL(wp), PARAMETER ::  dperctsc = -0.00054709752_wp
-    REAL(wp), PARAMETER ::  dperctsd = +0.0000063714823_wp
+    REAL(wp), PARAMETER ::  dperctsd =  0.0000063714823_wp
 !
 !-- Factor to mean standard deviation defining SIGNificance for 
@@ -2932 +2931 @@
 !
 !-- Initialise
-    sultr_res = +99._wp
+    sultr_res = 99._wp
     dperctm   = 0._wp
     dperctstd = 999999._wp
@@ -2952 +2951 @@
 !-- Value of the FUNCTION
     sultr_res = dperctm + faktor * dperctstd
-    IF ( ABS ( sultr_res ) > 99._wp )  sultr_res = +99._wp
+    IF ( ABS( sultr_res ) > 99._wp )  sultr_res = +99._wp
 
  END SUBROUTINE calc_sultr
@@ -2983 +2982 @@
     REAL(wp) ::  pmv_cross(2)
     REAL(wp) ::  reg_a(3)
-    REAL(wp) ::  coeff(3,5)
-    REAL(wp) ::  r_nenner
-    REAL(wp) ::  pmvc
-    REAL(wp) ::  dtmrt
-    REAL(wp) ::  sqrt_ws
-    INTEGER(iwp) ::  i
-    INTEGER(iwp) ::  i_bin
-!
-!-- Coefficient of the 3 regression lines
-!     1:const   2:*pmvc  3:*ta      4:*sqrt_ws  5:*dtmrt
-    coeff(1,1:5) =                                                             &
-       (/ +0.161_wp, +0.130_wp, -1.125E-03_wp, +1.106E-03_wp, -4.570E-04_wp /)
-    coeff(2,1:5) =                                                             &
-       (/ 0.795_wp, 0.713_wp, -8.880E-03_wp, -1.803E-03_wp, -2.816E-03_wp/)
-    coeff(3,1:5) =                                                             &
-       (/ +0.05761_wp, +0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp /)
+    REAL(wp) ::  r_denominator  !< the regression equations denominator
+    REAL(wp) ::  dtmrt          !< delta mean radiant temperature
+    REAL(wp) ::  sqrt_ws        !< sqare root of wind speed
+    INTEGER(iwp) ::  i          !< running index
+    INTEGER(iwp) ::  i_bin      !< result row number
+
+!    REAL(wp) ::  coeff(3,5)  !< unsafe! array is (re-)writable!
+!    coeff(1,1:5) =                                                             &
+!       (/ +0.161_wp,   +0.130_wp, -1.125E-03_wp, +1.106E-03_wp, -4.570E-04_wp /)
+!    coeff(2,1:5) =                                                             &
+!       (/  0.795_wp,    0.713_wp, -8.880E-03_wp, -1.803E-03_wp, -2.816E-03_wp /)
+!    coeff(3,1:5) =                                                             &
+!       (/ +0.05761_wp, +0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp /)
+
+!
+!-- Coefficient of the 3 regression lines:
+!      1:const      2:*pmva    3:*ta          4:*sqrt_ws     5:*dtmrt
+    REAL(wp), DIMENSION(1:3,1:5), PARAMETER ::  coeff = RESHAPE( (/            &
+        0.161_wp,   0.130_wp, -1.125E-03_wp,  1.106E-03_wp, -4.570E-04_wp,     &
+        0.795_wp,   0.713_wp, -8.880E-03_wp, -1.803E-03_wp, -2.816E-03_wp,     &
+        0.05761_wp, 0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp      &
+       /), SHAPE(coeff), order=(/ 2, 1 /) )
 !
 !-- Initialise
     nerr           = 0_iwp
     dpmv_cold_res  = 0._wp
-    pmvc           = pmva
     dtmrt          = tmrt - ta
     sqrt_ws        = ws
-    IF ( sqrt_ws < 0.10_wp )  THEN
-       sqrt_ws = 0.10_wp
+    IF ( sqrt_ws < 0.1_wp )  THEN
+       sqrt_ws = 0.1_wp
     ELSE
-       sqrt_ws = SQRT ( sqrt_ws )
+       sqrt_ws = SQRT( sqrt_ws )
     ENDIF
 
     delta_cold = 0._wp
-    pmv_cross = pmvc
-
+    pmv_cross = pmva
+
+!
+!-- Determine regression constant for given meteorological conditions
     DO  i = 1, 3
-!
-!--    Regression constant for given meteorological variables
        reg_a(i) = coeff(i,1) + coeff(i,3) * ta + coeff(i,4) *                  &
                   sqrt_ws + coeff(i,5)*dtmrt 
-       delta_cold(i) = reg_a(i) + coeff(i,2) * pmvc
+       delta_cold(i) = reg_a(i) + coeff(i,2) * pmva
     ENDDO
 !
 !-- Intersection points of regression lines in terms of Fanger's PMV
     DO  i = 1, 2
-       r_nenner = coeff(i,2) - coeff(i+1,2)
-       IF ( ABS ( r_nenner ) > 0.00001_wp )  THEN
-          pmv_cross(i) = ( reg_a(i+1) - reg_a(i) ) / r_nenner
+       r_denominator = coeff(i,2) - coeff(i+1,2)
+       IF ( ABS( r_denominator ) > 0.00001_wp )  THEN
+          pmv_cross(i) = ( reg_a(i+1) - reg_a(i) ) / r_denominator
        ELSE
           nerr = 1_iwp
@@ -3033 +3037 @@
        ENDIF
     ENDDO
-
+!
+!-- Select result row number
     i_bin = 3
     DO  i = 1, 2
@@ -3044 +3049 @@
 !-- Adjust to operative temperature scaled according 
 !-- to classical PMV (Fanger)
-    dpmv_cold_res = delta_cold(i_bin) - dpmv_adj(pmva)
+    dpmv_cold_res = delta_cold(i_bin) - dpmv_cold_adj(pmva)
 
  END SUBROUTINE dpmv_cold
@@ -3051 +3056 @@
 ! Description:
 ! ------------
-!> Calculates the summand dpmv_adj adjusting to the operative temperature 
-!> scaled according to classical PMV (Fanger)
-!> Reference environment: v_1m = 0.10 m/s, dTMRT = 0 K, r.h. = 50 %
-!------------------------------------------------------------------------------!
- REAL(wp) FUNCTION dpmv_adj( pmva )
+!> Calculates the summand dpmv_cold_adj adjusting to the operative temperature 
+!> scaled according to classical PMV (Fanger) for cold conditions.
+!> Valid for reference environment: v (1m) = 0.10 m/s, dTMRT = 0 K, r.h. = 50 %
+!------------------------------------------------------------------------------!
+ REAL(wp) FUNCTION dpmv_cold_adj( pmva )
 
     IMPLICIT NONE
 
-    REAL(wp), INTENT ( IN ) ::  pmva
-    INTEGER(iwp), PARAMETER ::  n_bin = 3
-    INTEGER(iwp), PARAMETER ::  n_ca = 0
-    INTEGER(iwp), PARAMETER ::  n_ce = 3
-    REAL(wp), dimension (n_bin, n_ca:n_ce) ::  coef
-
-    REAL(wp)      ::  pmv
-    INTEGER(iwp)  ::  i, i_bin
-!
-!--                                 range_1        range_2        range_3
-    DATA (coef(i,0), i = 1, n_bin) /0.0941540_wp, -0.1506620_wp, -0.0871439_wp/
-    DATA (coef(i,1), i = 1, n_bin) /0.0783162_wp, -1.0612651_wp,  0.1695040_wp/
-    DATA (coef(i,2), i = 1, n_bin) /0.1350144_wp, -1.0049144_wp, -0.0167627_wp/
-    DATA (coef(i,3), i = 1, n_bin) /0.1104037_wp, -0.2005277_wp, -0.0003230_wp/
-
-    IF ( pmva <= -2.1226_wp )  THEN
-       i_bin = 3_iwp
-    ELSE IF ( pmva <= -1.28_wp )  THEN
-       i_bin = 2_iwp
-    ELSE
-       i_bin = 1_iwp
-    ENDIF
-
-    dpmv_adj   = coef(i_bin,n_ca)
-    pmv        = 1._wp
-
-    DO  i = n_ca + 1, n_ce
-       pmv      = pmv * pmva
-       dpmv_adj = dpmv_adj + coef(i_bin,i) * pmv
+    REAL(wp), INTENT ( IN ) ::  pmva        !< (adjusted) Predicted Mean Vote
+
+    REAL(wp)      ::  pmv     !< pmv-part of the regression
+    INTEGER(iwp)  ::  i       !< running index
+    INTEGER(iwp)  ::  thr     !< thermal range
+!
+!-- Provide regression coefficients for three thermal ranges:
+!--    slightly cold  cold           very cold
+    REAL(wp), DIMENSION(1:3,0:3), PARAMETER ::  coef = RESHAPE( (/             &
+       0.0941540_wp, -0.1506620_wp, -0.0871439_wp,                             &
+       0.0783162_wp, -1.0612651_wp,  0.1695040_wp,                             &
+       0.1350144_wp, -1.0049144_wp, -0.0167627_wp,                             &
+       0.1104037_wp, -0.2005277_wp, -0.0003230_wp                              &
+       /), SHAPE(coef), order=(/ 2, 1 /) )
+!
+!-- Select thermal range
+    IF ( pmva <= -2.1226_wp )  THEN     !< very cold
+       thr = 3_iwp
+    ELSE IF ( pmva <= -1.28_wp )  THEN  !< cold
+       thr = 2_iwp
+    ELSE                                !< slightly cold
+       thr = 1_iwp
+    ENDIF
+!
+!-- Initialize
+    dpmv_cold_adj = coef(thr,0)
+    pmv           = 1._wp
+!
+!-- Calculate pmv adjustment (dpmv_cold_adj)
+    DO  i = 1, 3
+       pmv           = pmv * pmva
+       dpmv_cold_adj = dpmv_cold_adj + coef(thr,i) * pmv
     ENDDO
+
     RETURN
- END FUNCTION dpmv_adj
+ END FUNCTION dpmv_cold_adj
 
 !------------------------------------------------------------------------------!
@@ -3111 +3120 @@
 !
 !-- Type declaration for internal varables
-    REAL(wp)                     ::  top02
+    REAL(wp)                     ::  perct02
 !
 !-- Initialise
@@ -3117 +3126 @@
 !
 !-- Convert perceived temperature from basis 0.1 m/s to basis 0.2 m/s
-    top02 = 1.8788_wp + 0.9296_wp * perct_ij
+    perct02 = 1.8788_wp + 0.9296_wp * perct_ij
 !
 !-- IREQ neutral conditions (thermal comfort)
-    ireq_neutral = 1.62_wp - 0.0564_wp * top02
+    ireq_neutral = 1.62_wp - 0.0564_wp * perct02
 !
 !-- Regression only defined for perct <= 10 (degC)
@@ -3131 +3140 @@
 !
 !-- Minimal required clothing insulation: maximal acceptable body cooling
-    ireq_minimal = 1.26_wp - 0.0588_wp * top02
+    ireq_minimal = 1.26_wp - 0.0588_wp * perct02
     IF ( nerr > 0_iwp )  THEN
        ireq_minimal = ireq_neutral
@@ -3169 +3178 @@
 !
 !-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if 
-!-- height and weight be known. In: Arch. Int. Med.. 17, 1916, S. 863?871.
+!-- height and weight be known. In: Arch. Int. Med.. 17, 1916, pp. 863:871.
     surf = 0.007184_wp * height**0.725_wp * weight**0.425_wp
     RETURN
@@ -3276 +3285 @@
     REAL(wp) ::  d_std
     REAL(wp) ::  pmvs
-    REAL(wp) ::  top
     REAL(wp) ::  a_surf
 !     REAL(wp) ::  acti
@@ -3324 +3332 @@
 !--          between winter and summer set values
              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta , sclo,     &
-                            pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
+                            pmv_s, wclo, pmv_w, eps, pmva, ncount, clo )
              IF ( ncount < 0_iwp )  THEN
                 nerr = -1_iwp
@@ -3365 +3373 @@
 !--          between winter and summer set values
              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
-                               pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
+                               pmv_s, wclo, pmv_w, eps, pmva, ncount, clo )
              IF ( ncount < 0_wp )  THEN
                 nerr = -1_iwp
@@ -3406 +3414 @@
 !
 !--    Required clothing insulation (ireq) is exclusively defined for 
-!--    operative temperatures (top) less 10 (C) for a 
+!--    perceived temperatures (ipt) less 10 (C) for a 
 !--    reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
        clon = ireq_neutral ( ipt, ireq_minimal, nerr )
@@ -3493 +3501 @@
 !-- Determine perceived temperature by regression equation + adjustments
     CALL perct_regression( pmva, clo, ipt )
-
+!
+!-- Consider cold conditions
     IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
 !
@@ -3510 +3519 @@
     ptc = ipt
     CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
-    sultrieness    = .FALSE.
-    d_std      = -99._wp
+    sultrieness = .FALSE.
+    d_std       = -99._wp
     IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
 !
 !--    Adjust for warm/humid conditions according to Gagge 1986
        CALL saturation_vapor_pressure ( ta, svp_ta )
-       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
-       pmvs   = pmva + d_pmv
+       d_pmv = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
+       pmvs  = pmva + d_pmv
        CALL perct_regression( pmvs, clo, ipt )
        IF ( sult_lim < 99._wp )  THEN
@@ -3602 +3611 @@
 !
 !-- Heat_convection = forced convection
-    heat_convection = 12.1_wp * SQRT ( ws * pair / 1013.25_wp )
+    heat_convection = 12.1_wp * SQRT( ws * pair / 1013.25_wp )
 !
 !-- Average skin temperature
@@ -3618 +3627 @@
     niter = INT( dt * 10._wp, KIND=iwp )
     IF ( niter < 1 )  niter = 1_iwp
-    adjustrate = 1._wp - EXP ( -1._wp * ( 10._wp / time_equil ) * dt )
+    adjustrate = 1._wp - EXP( -1._wp * ( 10._wp / time_equil ) * dt )
     IF ( adjustrate >= 1._wp )  adjustrate = 1._wp
     adjustrate_cloth = adjustrate * 30._wp
     t_clothing = t_cloth
-
-    IF ( t_cloth <= -998.0_wp )  THEN  ! If initial run
+!
+!-- Set initial values for niter, adjustrates and t_clothing if this is the
+!-- first call
+    IF ( t_cloth <= -998._wp )  THEN  ! If initial run
        niter = 3_iwp
        adjustrate = 1._wp
@@ -3629 +3640 @@
        t_clothing = ta
     ENDIF
-
+!
+!-- Update clothing temperature
     DO  i = 1, niter
        t_clothing = t_clothing - adjustrate_cloth * ( ( t_clothing +           &
@@ -3638 +3650 @@
 !-- Empiric factor for the adaption of the heat ballance equation
 !-- to the psycho-physical scale (Equ. 40 in FANGER)
-    z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp )
+    z1 = ( .303_wp * EXP( -.036_wp * actlev ) + .0275_wp )
 !
 !-- Water vapour diffution through the skin
@@ -3742 +3754 @@
 
     CALL heat_exch( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, ht,  &
-            int_heat,mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa,      &
+            int_heat, mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa,     &
             vpts, wetsk )
 
@@ -3768 +3780 @@
     REAL(wp), INTENT( IN )  ::  ht        !< Persons height           (m)
     REAL(wp), INTENT( IN )  ::  work      !< Work load                (W)
-    REAL(wp), INTENT( IN )  ::  eta       !< Work efficiency      (dimensionless)
+    REAL(wp), INTENT( IN )  ::  eta       !< Work efficiency     (dimensionless)
 !
 !-- Output arguments: