Index: /palm/trunk/SOURCE/biometeorology_mod.f90
===================================================================
--- /palm/trunk/SOURCE/biometeorology_mod.f90	(revision 3748)
+++ /palm/trunk/SOURCE/biometeorology_mod.f90	(revision 3749)
@@ -27,4 +27,8 @@
 ! -----------------
 ! $Id$
+! - Added addittional safety meassures to bio_calculate_thermal_index_maps.
+! - Replaced several REAL (un-)equality comparisons.
+! 
+! 3742 2019-02-14 11:25:22Z dom_dwd_user
 ! - Allocation of the input _av grids was moved to the "sum" section of
 ! bio_3d_data_averaging to make sure averaging is only done once!
@@ -523,4 +527,6 @@
              ENDIF
 
+!
+!-- u_av, v_av and w_av are always allocated
              IF ( .NOT. ALLOCATED( u_av ) )  THEN
                 ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
@@ -866,6 +872,4 @@
  SUBROUTINE bio_check_parameters
 
-    USE control_parameters,                                                    &
-        ONLY:  message_string
 
     IMPLICIT NONE
@@ -1455,5 +1459,4 @@
  END SUBROUTINE bio_wrd_local
 
-
 !------------------------------------------------------------------------------!
 ! Description:
@@ -1584,5 +1587,8 @@
 !
 !--    Calculate ta from Tp assuming dry adiabatic laps rate
-       ta = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
+       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
+       ENDIF
 
        vp = bio_fill_value
@@ -1591,7 +1597,11 @@
        ENDIF
 
-       ws = ( 0.5_wp * ABS( u_av(k_wind, j, i) + u_av(k_wind, j, i+1) )  +     &
-          0.5_wp * ABS( v_av(k_wind, j, i) + v_av(k_wind, j+1, i) )  +         &
-          0.5_wp * ABS( w_av(k_wind, j, i) + w_av(k_wind+1, j, i) ) )
+       ws = bio_fill_value
+       IF ( ALLOCATED( u_av )  .AND.  ALLOCATED( v_av )  .AND.                 &
+          ALLOCATED( w_av ) )  THEN
+             ws = ( 0.5_wp * ABS( u_av(k_wind,j,i) + u_av(k_wind,j,i+1) )  +   &
+             0.5_wp * ABS( v_av(k_wind,j,i) + v_av(k_wind,j+1,i) )  +          &
+             0.5_wp * ABS( w_av(k_wind,j,i) + w_av(k_wind+1,j,i) ) )
+       ENDIF
     ELSE
 !
@@ -1604,7 +1614,7 @@
        ENDIF
 
-       ws = ( 0.5_wp * ABS( u(k_wind, j, i) + u(k_wind, j, i+1) )  +           &
-          0.5_wp * ABS( v(k_wind, j, i) + v(k_wind, j+1, i) )  +               &
-          0.5_wp * ABS( w(k_wind, j, i) + w(k_wind+1, j, i) ) )
+       ws = ( 0.5_wp * ABS( u(k_wind,j,i) + u(k_wind,j,i+1) )  +               &
+          0.5_wp * ABS( v(k_wind,j,i) + v(k_wind,j+1,i) )  +                   &
+          0.5_wp * ABS( w(k_wind,j,i) + w(k_wind+1,j,i) ) )
 
     ENDIF
@@ -1641,5 +1651,5 @@
 !> time_integration.f90: 1065ff
 !------------------------------------------------------------------------------!
- SUBROUTINE bio_calculate_thermal_index_maps ( av )
+ SUBROUTINE bio_calculate_thermal_index_maps( av )
 
     IMPLICIT NONE
@@ -1662,77 +1672,87 @@
 
 !
-!-- fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...)
-    CALL bio_calculate_mrt_grid ( av )
-
-    DO  i = nxl, nxr
-       DO  j = nys, nyn
-!
-!--       Determine local input conditions
-          tmrt_ij = bio_fill_value
-          vp      = bio_fill_value
-!
-!--       Determine input only if 
-          CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp,              &
-                                                ws, pair, tmrt_ij )
-!
-!--       Only proceed if input is available
-          pet_ij   = bio_fill_value   !< set fail value, e.g. valid for within
-          perct_ij = bio_fill_value   !< some obstacle
-          utci_ij  = bio_fill_value
-          IF ( .NOT. ( tmrt_ij <= -998._wp  .OR.  vp <= -998._wp ) )  THEN
-!
-!--          Calculate static thermal indices based on local tmrt
-             clo = bio_fill_value
-
-             IF ( do_calculate_perct  .OR.  do_calculate_perct_av )  THEN
-!
-!--          Estimate local perceived temperature
-                CALL calculate_perct_static( ta, vp, ws, tmrt_ij, pair, clo,   &
-                   perct_ij )
+!-- Check if some thermal index is desired. Don't do anything if, e.g. only
+!-- bio_mrt is desired.
+    IF ( do_calculate_perct  .OR.  do_calculate_perct_av  .OR.                 &
+       do_calculate_utci  .OR.  do_calculate_utci_av  .OR.                     &
+       do_calculate_pet  .OR.  do_calculate_pet_av )  THEN
+
+!
+!--    fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...)
+       CALL bio_calculate_mrt_grid ( av )
+
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+!
+!--          Determine local input conditions
+             tmrt_ij = bio_fill_value
+             vp      = bio_fill_value
+!
+!--          Determine local meteorological conditions
+             CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp,           &
+                                                   ws, pair, tmrt_ij )
+!
+!--          Only proceed if input is available
+             pet_ij   = bio_fill_value   !< set fail value, e.g. valid for
+             perct_ij = bio_fill_value   !< within some obstacle
+             utci_ij  = bio_fill_value
+             IF ( .NOT. ( tmrt_ij <= -998._wp  .OR.  vp <= -998._wp  .OR.      &
+                ws <= -998._wp  .OR.  ta <= -998._wp ) )  THEN
+!
+!--             Calculate static thermal indices based on local tmrt
+                clo = bio_fill_value
+
+                IF ( do_calculate_perct  .OR.  do_calculate_perct_av )  THEN
+!
+!--                Estimate local perceived temperature
+                   CALL calculate_perct_static( ta, vp, ws, tmrt_ij, pair,     &
+                      clo, perct_ij )
+                ENDIF
+
+                IF ( do_calculate_utci  .OR.  do_calculate_utci_av )  THEN
+!
+!--                Estimate local universal thermal climate index
+                   CALL calculate_utci_static( ta, vp, ws, tmrt_ij,            &
+                      bio_output_height, utci_ij )
+                ENDIF
+
+                IF ( do_calculate_pet  .OR.  do_calculate_pet_av )  THEN
+!
+!--                Estimate local physiologically equivalent temperature
+                   CALL calculate_pet_static( ta, vp, ws, tmrt_ij, pair,       &
+                      pet_ij )
+                ENDIF
              ENDIF
 
-             IF ( do_calculate_utci  .OR.  do_calculate_utci_av )  THEN
-!
-!--          Estimate local universal thermal climate index
-                CALL calculate_utci_static( ta, vp, ws, tmrt_ij,               &
-                   bio_output_height, utci_ij )
+
+             IF ( av )  THEN
+!
+!--             Write results for selected averaged indices
+                IF ( do_calculate_perct_av )  THEN
+                   perct_av(j, i) = perct_ij
+                ENDIF
+                IF ( do_calculate_utci_av )  THEN
+                   utci_av(j, i) = utci_ij
+                ENDIF
+                IF ( do_calculate_pet_av )  THEN
+                   pet_av(j, i)  = pet_ij
+                ENDIF
+             ELSE
+!
+!--             Write result for selected indices
+                IF ( do_calculate_perct )  THEN
+                   perct(j, i) = perct_ij
+                ENDIF
+                IF ( do_calculate_utci )  THEN
+                   utci(j, i) = utci_ij
+                ENDIF
+                IF ( do_calculate_pet )  THEN
+                   pet(j, i)  = pet_ij
+                ENDIF
              ENDIF
 
-             IF ( do_calculate_pet  .OR.  do_calculate_pet_av )  THEN
-!
-!--          Estimate local physiologically equivalent temperature
-                CALL calculate_pet_static( ta, vp, ws, tmrt_ij, pair, pet_ij )
-             ENDIF
-          ENDIF
-
-
-          IF ( av )  THEN
-!
-!--          Write results for selected averaged indices
-             IF ( do_calculate_perct_av )  THEN
-                perct_av(j, i) = perct_ij
-             ENDIF
-             IF ( do_calculate_utci_av )  THEN
-                utci_av(j, i) = utci_ij
-             ENDIF
-             IF ( do_calculate_pet_av )  THEN
-                pet_av(j, i)  = pet_ij
-             ENDIF
-          ELSE
-!
-!--          Write result for selected indices
-             IF ( do_calculate_perct )  THEN
-                perct(j, i) = perct_ij
-             ENDIF
-             IF ( do_calculate_utci )  THEN
-                utci(j, i) = utci_ij
-             ENDIF
-             IF ( do_calculate_pet )  THEN
-                pet(j, i)  = pet_ij
-             ENDIF
-          ENDIF
-
+          ENDDO
        ENDDO
-    ENDDO
+    ENDIF
 
  END SUBROUTINE bio_calculate_thermal_index_maps
@@ -1782,5 +1802,5 @@
 !--    Initialize instationary perceived temperature with personalized 
 !      PT as an initial guess, set actlev and clo
-       CALL ipt_init ( age, weight, height, sex, work, actlev, clo,            &
+       CALL ipt_init( age, weight, height, sex, work, actlev, clo,            &
           ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo,                   &
           ipt )
@@ -2294,5 +2314,5 @@
 !-- Determine perceived temperature by regression equation + adjustments
     pmvs = pmva
-    CALL perct_regression ( pmva, clo, perct_ij )
+    CALL perct_regression( pmva, clo, perct_ij )
     ptc = perct_ij
     IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
@@ -2306,5 +2326,5 @@
           pmvs   = -0.11_wp
        ENDIF
-       CALL perct_regression ( pmvs, clo, perct_ij )
+       CALL perct_regression( pmvs, clo, perct_ij )
     ENDIF
 !     clo_fanger = clo
@@ -2318,5 +2338,5 @@
        clo = clon
     ENDIF
-    CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
+    CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
     sultrieness    = .FALSE.
     d_std = -99._wp
@@ -2327,10 +2347,10 @@
        d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
        pmvs   = pmva + d_pmv
-       CALL perct_regression ( pmvs, clo, perct_ij )
+       CALL perct_regression( pmvs, clo, perct_ij )
        IF ( sult_lim < 99._wp )  THEN
           IF ( (perct_ij - ptc) > sult_lim )  sultrieness = .TRUE.
 !
 !--       Set factor to threshold for sultriness
-          IF ( dgtcstd /= 0_iwp )  THEN
+          IF ( ABS( dgtcstd ) > 0.00001_wp )  THEN
              d_std = ( ( perct_ij - ptc ) - dgtcm ) / dgtcstd
           ENDIF
@@ -2453,5 +2473,5 @@
                         y_average, top )
           sroot = SQRT ( y_average**2 - y_lower * y_upper )
-          IF ( sroot == 0._wp )  THEN
+          IF ( ABS( sroot ) < 0.00001_wp )  THEN
              clo_res = x_average
              nerr = j
@@ -2468,18 +2488,18 @@
           CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta,         &
                         y_new, top )
-          IF ( y_new == 0._wp )  THEN
+          IF ( ABS( y_new ) < 0.00001_wp )  THEN
              clo_res = x_ridder
              nerr       = j
              RETURN
           ENDIF
-          IF ( SIGN ( y_average, y_new ) /= y_average )  THEN
+          IF ( ABS( SIGN( y_average, y_new ) - y_average ) > 0.00001_wp )  THEN
              x_lower = x_average
              y_lower = y_average
              x_upper  = x_ridder
              y_upper  = y_new
-          ELSE IF ( SIGN ( y_lower, y_new ) /= y_lower )  THEN
+          ELSE IF ( ABS( SIGN( y_lower, y_new ) - y_lower ) > 0.00001_wp )  THEN
              x_upper  = x_ridder
              y_upper  = y_new
-          ELSE IF ( SIGN ( y_upper, y_new ) /= y_upper )  THEN
+          ELSE IF ( ABS( SIGN( y_upper, y_new ) - y_upper ) > 0.00001_wp )  THEN
              x_lower = x_ridder
              y_lower = y_new
@@ -2487,5 +2507,5 @@
 !
 !--          Never get here in x_ridder: singularity in y
-             nerr       = -1_iwp
+             nerr    = -1_iwp
              clo_res = x_ridder
              RETURN
@@ -2493,5 +2513,5 @@
           IF ( ABS ( x_upper - x_lower ) <= eps )  THEN
              clo_res = x_ridder
-             nerr       = j
+             nerr    = j
              RETURN
           ENDIF
@@ -2502,7 +2522,7 @@
        clo_res = y_new
        RETURN
-    ELSE IF ( y_lower == 0. )  THEN
+    ELSE IF ( ABS( y_lower ) < 0.00001_wp )  THEN
        x_ridder = clo_lower
-    ELSE IF ( y_upper == 0. )  THEN
+    ELSE IF ( ABS( y_upper ) < 0.00001_wp )  THEN
        x_ridder = clo_upper
     ELSE
@@ -2992,10 +3012,6 @@
     ENDIF
 
-    DO  i = 1, 3
-       delta_cold(i) = 0._wp
-       IF ( i < 3 )  THEN
-          pmv_cross(i) = pmvc
-       ENDIF
-    ENDDO
+    delta_cold = 0._wp
+    pmv_cross = pmvc
 
     DO  i = 1, 3
@@ -3004,5 +3020,5 @@
        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-1) = reg_a(i) + coeff(i,2) * pmvc
     ENDDO
 !
@@ -3011,5 +3027,5 @@
        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
+          pmv_cross(i-1) = ( reg_a(i+1) - reg_a(i) ) / r_nenner
        ELSE
           nerr = 1_iwp
@@ -3020,5 +3036,5 @@
     i_bin = 3
     DO  i = 1, 2
-       IF ( pmva > pmv_cross(i) )  THEN
+       IF ( pmva > pmv_cross(i-1) )  THEN
           i_bin = i
           EXIT
@@ -3028,5 +3044,5 @@
 !-- 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-1) - dpmv_adj(pmva)
 
  END SUBROUTINE dpmv_cold
@@ -3130,5 +3146,5 @@
 !> according to its height (m), weight (kg), and age (y)
 !------------------------------------------------------------------------------!
- SUBROUTINE surface_area ( height_cm, weight, age, surf )
+ SUBROUTINE surface_area( height_cm, weight, age, surf )
 
     IMPLICIT NONE
@@ -3174,5 +3190,5 @@
 !> for a sample human.
 !------------------------------------------------------------------------------!
- SUBROUTINE persdat ( age, weight, height, sex, work, a_surf, actlev )
+ SUBROUTINE persdat( age, weight, height, sex, work, a_surf, actlev )
 
     IMPLICIT NONE
@@ -3190,5 +3206,5 @@
     REAL(wp) ::  basic_heat_prod
 
-    CALL surface_area ( height, weight, INT( age ), a_surf )
+    CALL surface_area( height, weight, INT( age ), a_surf )
     s = height * 100._wp / ( weight**( 1._wp / 3._wp ) )
     factor = 1._wp + .004_wp  * ( 30._wp - age )
@@ -3215,5 +3231,5 @@
 !------------------------------------------------------------------------------!
 
- SUBROUTINE ipt_init (age, weight, height, sex, work, actlev, clo,             &
+ SUBROUTINE ipt_init( age, weight, height, sex, work, actlev, clo,             &
      ta, vp, ws, tmrt, pair, dt, storage, t_clothing,                          &
      ipt )
@@ -3270,5 +3286,5 @@
 
     storage = 0._wp
-    CALL persdat ( age, weight, height, sex, work, a_surf, actlev )
+    CALL persdat( age, weight, height, sex, work, a_surf, actlev )
 !
 !-- Initialise
@@ -3371,5 +3387,5 @@
 !-- Determine perceived temperature by regression equation + adjustments
     pmvs = pmva
-    CALL perct_regression ( pmva, clo, ipt )
+    CALL perct_regression( pmva, clo, ipt )
     ptc = ipt
     IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
@@ -3383,5 +3399,5 @@
           pmvs   = -0.11_wp
        ENDIF
-       CALL perct_regression ( pmvs, clo, ipt )
+       CALL perct_regression( pmvs, clo, ipt )
     ENDIF
 !     clo_fanger = clo
@@ -3395,5 +3411,5 @@
        clo = clon
     ENDIF
-    CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
+    CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
     sultrieness    = .FALSE.
     d_std      = -99._wp
@@ -3404,5 +3420,5 @@
        d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
        pmvs   = pmva + d_pmv
-       CALL perct_regression ( pmvs, clo, ipt )
+       CALL perct_regression( pmvs, clo, ipt )
        IF ( sult_lim < 99._wp )  THEN
           IF ( (ipt - ptc) > sult_lim )  sultrieness = .TRUE.
@@ -3476,5 +3492,5 @@
 !
 !-- Determine perceived temperature by regression equation + adjustments
-    CALL perct_regression ( pmva, clo, ipt )
+    CALL perct_regression( pmva, clo, ipt )
 
     IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
@@ -3488,10 +3504,10 @@
           pmvs   = -0.11_wp
        ENDIF
-       CALL perct_regression ( pmvs, clo, ipt )
+       CALL perct_regression( pmvs, clo, ipt )
     ENDIF
 !
 !-- Consider sultriness if appropriate
     ptc = ipt
-    CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
+    CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
     sultrieness    = .FALSE.
     d_std      = -99._wp
@@ -3502,5 +3518,5 @@
        d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
        pmvs   = pmva + d_pmv
-       CALL perct_regression ( pmvs, clo, ipt )
+       CALL perct_regression( pmvs, clo, ipt )
        IF ( sult_lim < 99._wp )  THEN
           IF ( (ipt - ptc) > sult_lim )  sultrieness = .TRUE.
@@ -3722,12 +3738,12 @@
 !
 !-- Call subfunctions
-    CALL in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta,    &
+    CALL in_body( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta,     &
             vpa, work )
 
-    CALL heat_exch ( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, ht, &
+    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,      &
             vpts, wetsk )
 
-    CALL pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair,      &
+    CALL pet_iteration( acl, adu, aeff, esw, facl, feff, int_heat, pair,       &
             rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk )
 
@@ -3741,5 +3757,5 @@
 !> Calculate internal energy ballance
 !------------------------------------------------------------------------------!
- SUBROUTINE in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, &
+ SUBROUTINE in_body( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, &
     vpa, work )
 !
@@ -3794,5 +3810,5 @@
 !> Calculate heat gain or loss
 !------------------------------------------------------------------------------!
- SUBROUTINE heat_exch ( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff,  &
+ SUBROUTINE 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,     &
         vpts, wetsk )
@@ -3953,5 +3969,5 @@
                   ( c(5) * tsk - c(0) - 5.28_wp * adu * tsk )
 
-          IF ( tsk == 36._wp )  tsk = 36.01_wp
+          IF ( ABS( tsk - 36._wp ) < 0.00001_wp )  tsk = 36.01_wp
           tcore(7) = c(0) / ( 5.28_wp * adu + c(1) * 6.3_wp / 3600._wp ) + tsk
           tcore(3) = c(0) / ( 5.28_wp * adu + ( c(1) * 6.3_wp / 3600._wp ) /   &
@@ -4077,5 +4093,5 @@
 !> Calculate PET
 !------------------------------------------------------------------------------!
- SUBROUTINE pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair,   &
+ SUBROUTINE pet_iteration( acl, adu, aeff, esw, facl, feff, int_heat, pair,    &
         rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk )
 !
