Index: palm/trunk/SOURCE/biometeorology_mod.f90
===================================================================
--- palm/trunk/SOURCE/biometeorology_mod.f90	(revision 3741)
+++ palm/trunk/SOURCE/biometeorology_mod.f90	(revision 3742)
@@ -27,4 +27,10 @@
 ! -----------------
 ! $Id$
+! - 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!
+! - Moved call of bio_calculate_thermal_index_maps from biometeorology module to
+! time_integration to make sure averaged input is updated before calculating.
+! 
+! 3740 2019-02-13 12:35:12Z dom_dwd_user
 ! - Added safety-meassure to catch the case that 'bio_mrt_av' is stated after
 ! 'bio_<index>' in the output section of the p3d file.
@@ -125,16 +131,16 @@
 
     USE arrays_3d,                                                             &
-       ONLY:  pt, p, u, v, w, q
+        ONLY:  pt, p, u, v, w, q
 
     USE averaging,                                                             &
-       ONLY:  pt_av, q_av, u_av, v_av, w_av
+        ONLY:  pt_av, q_av, u_av, v_av, w_av
 
     USE basic_constants_and_equations_mod,                                     &
-       ONLY:  c_p, degc_to_k, l_v, magnus, sigma_sb, pi
+        ONLY:  c_p, degc_to_k, l_v, magnus, sigma_sb, pi
 
     USE control_parameters,                                                    &
-       ONLY:  average_count_3d, biometeorology, dz, dz_stretch_factor,         &
-              dz_stretch_level, humidity, initializing_actions, nz_do3d,       &
-              surface_pressure
+        ONLY:  average_count_3d, biometeorology, dz, dz_stretch_factor,        &
+               dz_stretch_level, humidity, initializing_actions, nz_do3d,      &
+               surface_pressure
 
     USE date_and_time_mod,                                                     &
@@ -142,9 +148,9 @@
 
     USE grid_variables,                                                        &
-       ONLY:  ddx, dx, ddy, dy
+        ONLY:  ddx, dx, ddy, dy
 
     USE indices,                                                               &
-       ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg,    &
-              nysg, nyng
+        ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg,   &
+               nysg, nyng
 
     USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM
@@ -156,11 +162,11 @@
 !-- Import radiation model to obtain input for mean radiant temperature
     USE radiation_model_mod,                                                   &
-       ONLY: ix, iy, iz, id, mrt_nlevels, mrt_include_sw,                      &
-             mrtinsw, mrtinlw, mrtbl, nmrtbl, radiation,                       &
-             radiation_interactions, rad_sw_in,                                &
-             rad_sw_out, rad_lw_in, rad_lw_out
+        ONLY:  ix, iy, iz, id, mrt_nlevels, mrt_include_sw,                    &
+               mrtinsw, mrtinlw, mrtbl, nmrtbl, radiation,                     &
+               radiation_interactions, rad_sw_in,                              &
+               rad_sw_out, rad_lw_in, rad_lw_out
 
     USE surface_mod,                                                           &
-       ONLY: get_topography_top_index_ji
+        ONLY:  get_topography_top_index_ji
 
     IMPLICIT NONE
@@ -190,4 +196,5 @@
 !
 !--
+    LOGICAL ::  thermal_comfort  = .FALSE.  !< Enables or disables the entire thermal comfort part
     LOGICAL ::  do_average_theta = .FALSE.  !< switch: do theta averaging in this module? (if .FALSE. this is done globally)
     LOGICAL ::  do_average_q     = .FALSE.  !< switch: do e averaging in this module?
@@ -196,15 +203,13 @@
     LOGICAL ::  do_average_w     = .FALSE.  !< switch: do w averaging in this module?
     LOGICAL ::  do_average_mrt   = .FALSE.  !< switch: do mrt averaging in this module?
-    LOGICAL ::  average_trigger_perct = .FALSE.  !< update averaged input on call to bio_perct?
-    LOGICAL ::  average_trigger_utci  = .FALSE.  !< update averaged input on call to bio_utci?
-    LOGICAL ::  average_trigger_pet   = .FALSE.  !< update averaged input on call to bio_pet?
-
-    LOGICAL ::  thermal_comfort = .FALSE.  !< Enables or disables thermal comfort part
-    LOGICAL ::  do_calculate_perct     = .FALSE.   !< Turn index PT (instant. input) on or off
-    LOGICAL ::  do_calculate_perct_av  = .FALSE.   !< Turn index PT (averaged input) on or off
-    LOGICAL ::  do_calculate_pet       = .FALSE.   !< Turn index PET (instant. input) on or off
-    LOGICAL ::  do_calculate_pet_av    = .FALSE.   !< Turn index PET (averaged input) on or off
-    LOGICAL ::  do_calculate_utci      = .FALSE.   !< Turn index UTCI (instant. input) on or off
-    LOGICAL ::  do_calculate_utci_av   = .FALSE.   !< Turn index UTCI (averaged input) on or off
+    LOGICAL ::  average_trigger_perct  = .FALSE.  !< update averaged input on call to bio_perct?
+    LOGICAL ::  average_trigger_utci   = .FALSE.  !< update averaged input on call to bio_utci?
+    LOGICAL ::  average_trigger_pet    = .FALSE.  !< update averaged input on call to bio_pet?
+    LOGICAL ::  do_calculate_perct     = .FALSE.  !< Turn index PT (instant. input) on or off
+    LOGICAL ::  do_calculate_perct_av  = .FALSE.  !< Turn index PT (averaged input) on or off
+    LOGICAL ::  do_calculate_pet       = .FALSE.  !< Turn index PET (instant. input) on or off
+    LOGICAL ::  do_calculate_pet_av    = .FALSE.  !< Turn index PET (averaged input) on or off
+    LOGICAL ::  do_calculate_utci      = .FALSE.  !< Turn index UTCI (instant. input) on or off
+    LOGICAL ::  do_calculate_utci_av   = .FALSE.  !< Turn index UTCI (averaged input) on or off
 
 !
@@ -399,15 +404,10 @@
 
           CASE ( 'bio_mrt' )
-!
-!--             Consider the case 'bio_mrt' is called after some thermal index. 
-!--             In that case do_average_mrt will be .TRUE. leading to a double 
-!--             averaging.  The other averaged input quantities are core 
-!--             components, that will  hopefully always be treated before the 
-!--             bio methods are called.
-                IF ( .NOT. do_average_mrt .AND.                                &
-                      .NOT. ALLOCATED( mrt_av_grid ) )  THEN
+
+                IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
                    ALLOCATE( mrt_av_grid(nmrtbl) )
                 ENDIF
                 mrt_av_grid = 0.0_wp
+                do_average_mrt = .FALSE.  !< overwrite if that was enabled somehow
 
 
@@ -425,63 +425,23 @@
 !--          Only proceed here if this was not done for any index before. This
 !--          is done only once during the whole model run.
-             IF ( .NOT. average_trigger_perct .AND.                            &
-                  .NOT. average_trigger_utci  .AND.                            &
-                  .NOT. average_trigger_pet ) THEN
+             IF ( .NOT. average_trigger_perct  .AND.                           &
+                  .NOT. average_trigger_utci   .AND.                           &
+                  .NOT. average_trigger_pet )  THEN
 !
 !--             Memorize the first index called to control averaging
-                IF ( TRIM( variable ) == 'bio_perct*' ) THEN
+                IF ( TRIM( variable ) == 'bio_perct*' )  THEN
                     average_trigger_perct = .TRUE.
                 ENDIF
-                IF ( TRIM( variable ) == 'bio_utci*' ) THEN
+                IF ( TRIM( variable ) == 'bio_utci*' )  THEN
                     average_trigger_utci = .TRUE.
                 ENDIF
-                IF ( TRIM( variable ) == 'bio_pet*' ) THEN
+                IF ( TRIM( variable ) == 'bio_pet*' )  THEN
                     average_trigger_pet = .TRUE.
                 ENDIF
              ENDIF
 !
-!--          Only continue if var is the index, that controls averaging.
-!--          Break immediatelly (doing nothing) for the other indices.
-             IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') RETURN
-             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*')   RETURN
-             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'bio_pet*')    RETURN
-!
-!--          Now memorize which of the input grids are not averaged by other 
-!--          modules. Set averaging switch to .TRUE. in that case.
-             IF ( .NOT. ALLOCATED( pt_av ) )  THEN
-                ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                do_average_theta = .TRUE.
-                pt_av = 0.0_wp
-             ENDIF
-
-             IF ( .NOT. ALLOCATED( q_av ) )  THEN
-                ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                do_average_q = .TRUE.
-                q_av = 0.0_wp
-             ENDIF
-
-             IF ( .NOT. ALLOCATED( u_av ) )  THEN
-                ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                do_average_u = .TRUE.
-                u_av = 0.0_wp
-             ENDIF
-
-             IF ( .NOT. ALLOCATED( v_av ) )  THEN
-                ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                do_average_v = .TRUE.
-                v_av = 0.0_wp
-             ENDIF
-
-             IF ( .NOT. ALLOCATED( w_av ) )  THEN
-                ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                do_average_w = .TRUE.
-                w_av = 0.0_wp
-             ENDIF
-
-             IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
-                ALLOCATE( mrt_av_grid(nmrtbl) )
-                do_average_mrt = .TRUE.
-                mrt_av_grid = 0.0_wp
-             ENDIF
+!--          Allocation of the input _av grids was moved to the "sum" section to
+!--          make sure averaging is only done once!
+
 
           CASE ( 'uvem_vitd3dose*' )
@@ -505,5 +465,5 @@
 !--          that case do_average_mrt will be .TRUE. leading to a double-
 !--          averaging. 
-             IF ( .NOT. do_average_mrt .AND. ALLOCATED( mrt_av_grid ) )  THEN
+             IF ( .NOT. do_average_mrt  .AND.  ALLOCATED( mrt_av_grid ) )  THEN
 
                 IF ( mrt_include_sw )  THEN
@@ -511,9 +471,9 @@
                       ( ( human_absorb * mrtinsw(:) +                          &
                       mrtinlw(:) ) /                                           &
-                      ( human_emiss * sigma_sb ) ) ** .25_wp - degc_to_k
+                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
                 ELSE
                    mrt_av_grid(:) = mrt_av_grid(:) +                           &
                       ( mrtinlw(:) /                                           &
-                      ( human_emiss * sigma_sb ) ) ** .25_wp - degc_to_k
+                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
                 ENDIF
              ENDIF
@@ -523,12 +483,20 @@
 !--          Only continue if the current index is the one to trigger the input 
 !--          averaging, see above
-             IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') &
-                RETURN
-             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*')  &
-                RETURN
-             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'bio_pet*')   &
-                RETURN
-
-             IF ( ALLOCATED( pt_av ) .AND. do_average_theta ) THEN
+             IF ( average_trigger_perct  .AND.  TRIM( variable ) /=            &
+                'bio_perct*')  RETURN
+             IF ( average_trigger_utci   .AND.  TRIM( variable ) /=            &
+                'bio_utci*')   RETURN
+             IF ( average_trigger_pet    .AND.  TRIM( variable ) /=            &
+                'bio_pet*')    RETURN
+!
+!--          Now memorize which of the input grids are not averaged by other 
+!--          modules. Set averaging switch to .TRUE. and allocate the respective
+!--          grid in that case.
+             IF ( .NOT. ALLOCATED( pt_av ) )  THEN  !< if not averaged by other module
+                ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                do_average_theta = .TRUE.  !< memorize, that bio is responsible
+                pt_av = 0.0_wp
+             ENDIF
+             IF ( ALLOCATED( pt_av )  .AND.  do_average_theta )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -540,5 +508,10 @@
              ENDIF
 
-             IF ( ALLOCATED( q_av )  .AND. do_average_q ) THEN
+             IF ( .NOT. ALLOCATED( q_av ) )  THEN
+                ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                do_average_q = .TRUE.
+                q_av = 0.0_wp
+             ENDIF
+             IF ( ALLOCATED( q_av )  .AND.  do_average_q )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -550,5 +523,10 @@
              ENDIF
 
-             IF ( ALLOCATED( u_av )  .AND. do_average_u ) THEN
+             IF ( .NOT. ALLOCATED( u_av ) )  THEN
+                ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                do_average_u = .TRUE.
+                u_av = 0.0_wp
+             ENDIF
+             IF ( ALLOCATED( u_av )  .AND.  do_average_u )  THEN
                 DO  i = nxlg, nxrg       !< yes, ghost points are required here!
                    DO  j = nysg, nyng
@@ -560,5 +538,10 @@
              ENDIF
 
-             IF ( ALLOCATED( v_av )  .AND. do_average_v ) THEN
+             IF ( .NOT. ALLOCATED( v_av ) )  THEN
+                ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                do_average_v = .TRUE.
+                v_av = 0.0_wp
+             ENDIF
+             IF ( ALLOCATED( v_av )  .AND.  do_average_v )  THEN
                 DO  i = nxlg, nxrg       !< yes, ghost points are required here!
                    DO  j = nysg, nyng
@@ -570,5 +553,10 @@
              ENDIF
 
-             IF ( ALLOCATED( w_av )  .AND. do_average_w ) THEN
+             IF ( .NOT. ALLOCATED( w_av ) )  THEN
+                ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                do_average_w = .TRUE.
+                w_av = 0.0_wp
+             ENDIF
+             IF ( ALLOCATED( w_av )  .AND.  do_average_w )  THEN
                 DO  i = nxlg, nxrg       !< yes, ghost points are required here!
                    DO  j = nysg, nyng
@@ -580,15 +568,20 @@
              ENDIF
 
-             IF ( ALLOCATED( mrt_av_grid ) .AND. do_average_mrt )  THEN
+             IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
+                ALLOCATE( mrt_av_grid(nmrtbl) )
+                do_average_mrt = .TRUE.
+                mrt_av_grid = 0.0_wp
+             ENDIF
+             IF ( ALLOCATED( mrt_av_grid )  .AND.  do_average_mrt )  THEN
 
                 IF ( mrt_include_sw )  THEN
                    mrt_av_grid(:) = mrt_av_grid(:) +                           &
                       ( ( human_absorb * mrtinsw(:) +                          &
-                      mrtinlw(:) ) /                             &  ! human_emiss * mrtinlw(:) /
-                      ( human_emiss * sigma_sb ) ) ** .25_wp - degc_to_k
+                      mrtinlw(:) ) /                                           &
+                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
                 ELSE
                    mrt_av_grid(:) = mrt_av_grid(:) +                           &
-                      ( mrtinlw(:) /                             &  ! ( human_emiss * mrtinlw(:) /
-                      ( human_emiss * sigma_sb ) ) ** .25_wp - degc_to_k
+                      ( mrtinlw(:) /                                           &
+                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
                 ENDIF
              ENDIF
@@ -596,5 +589,5 @@
 !--       This is a cumulated dose. No mode == 'average' for this quantity.
           CASE ( 'uvem_vitd3dose*' )
-             IF ( ALLOCATED( vitd3_exposure_av ) ) THEN
+             IF ( ALLOCATED( vitd3_exposure_av ) )  THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -618,5 +611,5 @@
 !--          that case do_average_mrt will be .TRUE. leading to a double-
 !--          averaging. 
-             IF ( .NOT. do_average_mrt .AND. ALLOCATED( mrt_av_grid ) )  THEN
+             IF ( .NOT. do_average_mrt  .AND.  ALLOCATED( mrt_av_grid ) )  THEN
                 mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp )
              ENDIF
@@ -625,12 +618,12 @@
 !
 !--          Only continue if update index, see above
-             IF ( average_trigger_perct .AND.                                  &
-                TRIM( variable ) /= 'bio_perct*' ) RETURN
-             IF ( average_trigger_utci .AND.                                   &
-                TRIM( variable ) /= 'bio_utci*' ) RETURN
-             IF ( average_trigger_pet  .AND.                                   &
-                TRIM( variable ) /= 'bio_pet*' ) RETURN
-
-             IF ( ALLOCATED( pt_av ) .AND. do_average_theta ) THEN
+             IF ( average_trigger_perct  .AND.                                 &
+                TRIM( variable ) /= 'bio_perct*' )  RETURN
+             IF ( average_trigger_utci  .AND.                                  &
+                TRIM( variable ) /= 'bio_utci*' )  RETURN
+             IF ( average_trigger_pet   .AND.                                  &
+                TRIM( variable ) /= 'bio_pet*' )  RETURN
+
+             IF ( ALLOCATED( pt_av )  .AND.  do_average_theta )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -643,5 +636,5 @@
              ENDIF
 
-             IF ( ALLOCATED( q_av ) .AND. do_average_q ) THEN
+             IF ( ALLOCATED( q_av )  .AND.  do_average_q )  THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -654,5 +647,5 @@
              ENDIF
 
-             IF ( ALLOCATED( u_av ) .AND. do_average_u ) THEN
+             IF ( ALLOCATED( u_av )  .AND.  do_average_u )  THEN
                 DO  i = nxlg, nxrg       !< yes, ghost points are required here!
                    DO  j = nysg, nyng
@@ -665,5 +658,5 @@
              ENDIF
 
-             IF ( ALLOCATED( v_av ) .AND. do_average_v ) THEN
+             IF ( ALLOCATED( v_av )  .AND.  do_average_v )  THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -676,5 +669,5 @@
              ENDIF
 
-             IF ( ALLOCATED( w_av ) .AND. do_average_w ) THEN
+             IF ( ALLOCATED( w_av )  .AND.  do_average_w )  THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -687,11 +680,8 @@
              ENDIF
 
-             IF ( ALLOCATED( mrt_av_grid ) .AND. do_average_mrt )  THEN
-                mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp )
-             ENDIF
-
-!
-!--          Udate all thermal index grids with updated averaged input
-             CALL bio_calculate_thermal_index_maps ( .TRUE. )
+             IF ( ALLOCATED( mrt_av_grid )  .AND.  do_average_mrt )  THEN
+                mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d,      &
+                KIND=wp )
+             ENDIF
 
 !
@@ -730,9 +720,10 @@
     SELECT CASE ( TRIM( var ) )
 !
-!-- Allocate a temporary array with the desired output dimensions.
-!-- Arrays for time-averaged thermal indices are also allocated here because they are not running
-!-- through the standard averaging procedure in bio_3d_data_averaging as the values of the
-!-- averaged thermal indices are derived in a single step based on priorly averaged arrays
-!-- (see bio_calculate_thermal_index_maps).
+!--    Allocate a temporary array with the desired output dimensions.
+!--    Arrays for time-averaged thermal indices are also allocated here because
+!--    they are not running through the standard averaging procedure in
+!--    bio_3d_data_averaging as the values of the averaged thermal indices are
+!--    derived in a single step based on priorly averaged arrays (see
+!--    bio_calculate_thermal_index_maps).
        CASE ( 'bio_mrt' )
           unit = 'degree_C'
@@ -746,5 +737,5 @@
           unit = 'degree_C'
           thermal_comfort = .TRUE.
-          IF ( j == 0 ) THEN                !< if instantaneous input
+          IF ( j == 0 )  THEN                !< if instantaneous input
              do_calculate_perct = .TRUE.
              IF ( .NOT. ALLOCATED( perct ) )  THEN
@@ -763,5 +754,5 @@
           unit = 'degree_C'
           thermal_comfort = .TRUE.
-          IF ( j == 0 ) THEN
+          IF ( j == 0 )  THEN
              do_calculate_utci = .TRUE.
              IF ( .NOT. ALLOCATED( utci ) )  THEN
@@ -780,5 +771,5 @@
           unit = 'degree_C'
           thermal_comfort = .TRUE.
-          IF ( j == 0 ) THEN
+          IF ( j == 0 )  THEN
              do_calculate_pet = .TRUE.
              IF ( .NOT. ALLOCATED( pet ) )  THEN
@@ -838,9 +829,9 @@
 !
 !--    Further checks if thermal comfort output is desired.
-    IF ( thermal_comfort .AND. unit == 'degree_C' )  THEN
+    IF ( thermal_comfort  .AND.  unit == 'degree_C' )  THEN
 
 !
 !--    Break if required modules "radiation" and "humidity" are not avalable.
-       IF ( .NOT.  radiation ) THEN
+       IF ( .NOT.  radiation )  THEN
           message_string = 'output of "' // TRIM( var ) // '" require'         &
                            // 's radiation = .TRUE.'
@@ -855,5 +846,5 @@
           unit = 'illegal'
        ENDIF
-       IF ( mrt_nlevels == 0 ) THEN
+       IF ( mrt_nlevels == 0 )  THEN
           message_string = 'output of "' // TRIM( var ) // '" require'         &
                            // 's mrt_nlevels > 0'
@@ -934,15 +925,15 @@
               j = mrtbl(iy,l)
               k = mrtbl(iz,l)
-              IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. j > nyn .OR.   &
-                 i < nxl .OR. i > nxr ) CYCLE
+              IF ( k < nzb_do  .OR.  k > nzt_do  .OR.  j < nys  .OR.           &
+                 j > nyn  .OR.  i < nxl  .OR.  i > nxr )  CYCLE
               IF ( av == 0 )  THEN
                  IF ( mrt_include_sw )  THEN
                     local_pf(i,j,k) = ( ( human_absorb * mrtinsw(l) +          &
                                     mrtinlw(l) ) /                             &
-                                    ( human_emiss * sigma_sb ) ) ** .25_wp -   &
+                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
                                     degc_to_k
                  ELSE
                     local_pf(i,j,k) = ( mrtinlw(l)  /                          &
-                                    ( human_emiss * sigma_sb ) ) ** .25_wp -   &
+                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
                                     degc_to_k
                  ENDIF
@@ -968,5 +959,5 @@
                  ENDDO
               ENDDO
-           END IF
+           ENDIF
 
 
@@ -986,5 +977,5 @@
                  ENDDO
               ENDDO
-           END IF
+           ENDIF
 
 
@@ -1004,5 +995,5 @@
                  ENDDO
               ENDDO
-           END IF
+           ENDIF
 
 !
@@ -1087,15 +1078,15 @@
                j = mrtbl(iy,l)
                k = mrtbl(iz,l)
-               IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. j > nyn .OR.  &
-                  i < nxl .OR. i > nxr ) CYCLE
+               IF ( k < nzb_do  .OR.  k > nzt_do  .OR.  j < nys  .OR.          &
+                  j > nyn  .OR.  i < nxl  .OR.  i > nxr )  CYCLE
                IF ( av == 0 )  THEN
                   IF ( mrt_include_sw )  THEN
                      local_pf(i,j,k) = REAL( ( ( human_absorb * mrtinsw(l) +   &
                                     mrtinlw(l) ) /                             &
-                                    ( human_emiss * sigma_sb ) ) ** .25_wp -   &
+                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
                                     degc_to_k, KIND = sp )
                   ELSE
                      local_pf(i,j,k) = REAL( ( mrtinlw(l) /                    &
-                                    ( human_emiss * sigma_sb ) ) ** .25_wp -   &
+                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
                                     degc_to_k, KIND = sp )
                   ENDIF
@@ -1147,13 +1138,13 @@
     is2d = ( var(l-1:l) == 'xy' )
 
-    IF ( var(1:4) == 'bio_' ) THEN
+    IF ( var(1:4) == 'bio_' )  THEN
        found  = .TRUE.
        grid_x = 'x'
        grid_y = 'y'
        grid_z = 'zu'
-       IF ( is2d  .AND.  var(1:7) /= 'bio_mrt' ) grid_z = 'zu1'
-    ENDIF
-
-    IF ( is2d  .AND.  var(1:4) == 'uvem' ) THEN
+       IF ( is2d  .AND.  var(1:7) /= 'bio_mrt' )  grid_z = 'zu1'
+    ENDIF
+
+    IF ( is2d  .AND.  var(1:4) == 'uvem' )  THEN
        grid_x = 'x'
        grid_y = 'y'
@@ -1241,5 +1232,5 @@
         ONLY: message_string
 
-    IF ( .NOT. radiation_interactions ) THEN
+    IF ( .NOT. radiation_interactions )  THEN
        message_string = 'The mrt calculation requires ' //                     &
                         'enabled radiation_interactions but it ' //            &
@@ -1485,5 +1476,5 @@
 !
 !-- We need to differentiate if averaged input is desired (av == .TRUE.) or not.
-    IF ( av ) THEN
+    IF ( av )  THEN
 !
 !-- Make sure tmrt_av_grid is present and initialize with the fill value
@@ -1505,5 +1496,5 @@
              k = mrtbl(iz,l)
              IF ( k - get_topography_top_index_ji( j, i, 's' ) ==              &
-                   bio_cell_level + 1_iwp) THEN
+                   bio_cell_level + 1_iwp)  THEN
 !
 !-- Averaging was done before, so we can just copy the result here
@@ -1531,13 +1522,13 @@
           k = mrtbl(iz,l)
           IF ( k - get_topography_top_index_ji( j, i, 's' ) ==                 &
-                bio_cell_level + 1_iwp) THEN
+                bio_cell_level + 1_iwp)  THEN
              IF ( mrt_include_sw )  THEN
                  tmrt_grid(j,i) = ( ( human_absorb * mrtinsw(l) +              &
                                   mrtinlw(l) )  /                              &
-                                  ( human_emiss * sigma_sb ) ) ** .25_wp -     &
+                                  ( human_emiss * sigma_sb ) )**.25_wp -       &
                                   degc_to_k
              ELSE
                  tmrt_grid(j,i) = ( mrtinlw(l)  /                              &
-                                  ( human_emiss * sigma_sb ) ) ** .25_wp -     &
+                                  ( human_emiss * sigma_sb ) )**.25_wp -       &
                                   degc_to_k
              ENDIF
@@ -1585,17 +1576,17 @@
 !-- Avoid non-representative horizontal u and v of 0.0 m/s too close to ground
     k_wind = k
-    IF ( bio_cell_level < 1_iwp ) THEN
+    IF ( bio_cell_level < 1_iwp )  THEN
        k_wind = k + 1_iwp
     ENDIF
 !
 !-- Determine local values:
-    IF ( average_input ) THEN
+    IF ( average_input )  THEN
 !
 !--    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 = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
 
        vp = bio_fill_value
-       IF ( humidity .AND. ALLOCATED( q_av ) ) THEN
-          vp = q_av(k, j, i)
+       IF ( humidity  .AND.  ALLOCATED( q_av ) )  THEN
+          vp = q_av(k,j,i)
        ENDIF
 
@@ -1606,9 +1597,9 @@
 !
 !-- Calculate ta from Tp assuming dry adiabatic laps rate
-       ta = pt(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
+       ta = pt(k,j,i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
 
        vp = bio_fill_value
-       IF ( humidity ) THEN
-          vp = q(k, j, i)
+       IF ( humidity )  THEN
+          vp = q(k,j,i)
        ENDIF
 
@@ -1625,13 +1616,13 @@
 !-- The magnus formula is limited to temperatures up to 333.15 K to
 !   avoid negative values of vp_sat
-    IF ( vp > -998._wp ) THEN
+    IF ( vp > -998._wp )  THEN
        vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) )
        vp  = vp * pair / ( vp + 0.622_wp )
-       IF ( vp > vp_sat ) vp = vp_sat
+       IF ( vp > vp_sat )  vp = vp_sat
     ENDIF
 !
 !-- local mtr value at [i,j]
     tmrt = bio_fill_value  !< this can be a valid result (e.g. for inside some ostacle)
-    IF ( .NOT. average_input ) THEN
+    IF ( .NOT. average_input )  THEN
 !
 !--    Use MRT from RTM precalculated in tmrt_grid
@@ -1674,6 +1665,6 @@
     CALL bio_calculate_mrt_grid ( av )
 
-    DO i = nxl, nxr
-       DO j = nys, nyn
+    DO  i = nxl, nxr
+       DO  j = nys, nyn
 !
 !--       Determine local input conditions
@@ -1689,10 +1680,10 @@
           perct_ij = bio_fill_value   !< some obstacle
           utci_ij  = bio_fill_value
-          IF ( .NOT. ( tmrt_ij <= -998._wp .OR. vp <= -998._wp ) ) THEN
+          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
+             IF ( do_calculate_perct  .OR.  do_calculate_perct_av )  THEN
 !
 !--          Estimate local perceived temperature
@@ -1701,5 +1692,5 @@
              ENDIF
 
-             IF ( do_calculate_utci  .OR. do_calculate_utci_av ) THEN
+             IF ( do_calculate_utci  .OR.  do_calculate_utci_av )  THEN
 !
 !--          Estimate local universal thermal climate index
@@ -1708,24 +1699,24 @@
              ENDIF
 
-             IF ( do_calculate_pet .OR. do_calculate_pet_av ) THEN
+             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
-          END IF
-
-
-          IF ( av ) THEN
+          ENDIF
+
+
+          IF ( av )  THEN
 !
 !--          Write results for selected averaged indices
              IF ( do_calculate_perct_av )  THEN
                 perct_av(j, i) = perct_ij
-             END IF
-             IF ( do_calculate_utci_av ) THEN
+             ENDIF
+             IF ( do_calculate_utci_av )  THEN
                 utci_av(j, i) = utci_ij
-             END IF
-             IF ( do_calculate_pet_av ) THEN
+             ENDIF
+             IF ( do_calculate_pet_av )  THEN
                 pet_av(j, i)  = pet_ij
-             END IF
+             ENDIF
           ELSE
 !
@@ -1733,15 +1724,15 @@
              IF ( do_calculate_perct )  THEN
                 perct(j, i) = perct_ij
-             END IF
-             IF ( do_calculate_utci ) THEN
+             ENDIF
+             IF ( do_calculate_utci )  THEN
                 utci(j, i) = utci_ij
-             END IF
-             IF ( do_calculate_pet ) THEN
+             ENDIF
+             IF ( do_calculate_pet )  THEN
                 pet(j, i)  = pet_ij
-             END IF
-          END IF
-
-       END DO
-    END DO
+             ENDIF
+          ENDIF
+
+       ENDDO
+    ENDDO
 
  END SUBROUTINE bio_calculate_thermal_index_maps
@@ -1782,10 +1773,10 @@
 !
 !-- return immediatelly if nothing to do!
-    IF ( .NOT. thermal_comfort ) THEN
+    IF ( .NOT. thermal_comfort )  THEN
         RETURN
     ENDIF
 !
 !-- If clo equals the initial value, this is the initial call
-    IF ( clo <= -998._wp ) THEN
+    IF ( clo <= -998._wp )  THEN
 !
 !--    Initialize instationary perceived temperature with personalized 
@@ -1881,6 +1872,6 @@
 !
 !-- Check if input values in range after Broede et al. (2012)
-    IF ( ( d_tmrt > 70._wp ) .OR. ( d_tmrt < -30._wp ) .OR.                    &
-       ( vp >= 50._wp ) ) THEN
+    IF ( ( d_tmrt > 70._wp )  .OR.  ( d_tmrt < -30._wp )  .OR.                 &
+       ( vp >= 50._wp ) )  THEN
        utci_ij = bio_fill_value
        RETURN
@@ -1888,9 +1879,9 @@
 !
 !-- Apply eq. 2 in Broede et al. (2012) for ta out of bounds
-    IF ( ta > 50._wp ) THEN
+    IF ( ta > 50._wp )  THEN
        offset = ta - 50._wp
        ta = 50._wp
     ENDIF
-    IF ( ta < -50._wp ) THEN
+    IF ( ta < -50._wp )  THEN
        offset = ta + 50._wp
        ta = -50._wp
@@ -1900,6 +1891,6 @@
 !-- humidity values below 0.5 m/s or 5%, respectively, the
 !-- user is advised to use the lower bounds for the calculations.
-    IF ( va < 0.5_wp ) va = 0.5_wp
-    IF ( va > 17._wp ) va = 17._wp
+    IF ( va < 0.5_wp )  va = 0.5_wp
+    IF ( va > 17._wp )  va = 17._wp
 
 !
@@ -2229,5 +2220,5 @@
 !
 !-- decision: firstly calculate for winter or summer clothing
-    IF ( ta <= 10._wp ) THEN
+    IF ( ta <= 10._wp )  THEN
 !
 !--    First guess: winter clothing insulation: cold stress
@@ -2236,5 +2227,5 @@
        pmv_w = pmva
 
-       IF ( pmva > 0._wp ) THEN
+       IF ( pmva > 0._wp )  THEN
 !
 !--       Case summer clothing insulation: heat load ?
@@ -2243,5 +2234,5 @@
              top )
           pmv_s = pmva
-          IF ( pmva <= 0._wp ) THEN
+          IF ( pmva <= 0._wp )  THEN
 !
 !--          Case: comfort achievable by varying clothing insulation 
@@ -2249,14 +2240,14 @@
              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
                 pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
-             IF ( ncount < 0_iwp ) THEN
+             IF ( ncount < 0_iwp )  THEN
                 nerr = -1_iwp
                 RETURN
              ENDIF
-          ELSE IF ( pmva > 0.06_wp ) THEN
+          ELSE IF ( pmva > 0.06_wp )  THEN
              clo = 0.5_wp
              CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
                            pmva, top )
           ENDIF
-       ELSE IF ( pmva < -0.11_wp ) THEN
+       ELSE IF ( pmva < -0.11_wp )  THEN
           clo = 1.75_wp
           CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
@@ -2270,5 +2261,5 @@
        pmv_s = pmva
 
-       IF ( pmva < 0._wp ) THEN
+       IF ( pmva < 0._wp )  THEN
 !
 !--       Case winter clothing insulation: cold stress ?
@@ -2278,5 +2269,5 @@
           pmv_w = pmva
 
-          IF ( pmva >= 0._wp ) THEN
+          IF ( pmva >= 0._wp )  THEN
 !
 !--          Case: comfort achievable by varying clothing insulation 
@@ -2284,14 +2275,14 @@
              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
                                pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
-             IF ( ncount < 0_iwp ) THEN
+             IF ( ncount < 0_iwp )  THEN
                 nerr = -1_iwp
                 RETURN
              ENDIF
-          ELSE IF ( pmva < -0.11_wp ) THEN
+          ELSE IF ( pmva < -0.11_wp )  THEN
              clo = 1.75_wp
              CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
                            pmva, top )
           ENDIF
-       ELSE IF ( pmva > 0.06_wp ) THEN
+       ELSE IF ( pmva > 0.06_wp )  THEN
           clo = 0.5_wp
           CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
@@ -2305,11 +2296,11 @@
     CALL perct_regression ( pmva, clo, perct_ij )
     ptc = perct_ij
-    IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
+    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
 !
 !--    Adjust for cold conditions according to Gagge 1986
        CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
-       IF ( nerr_cold > 0_iwp ) nerr = -5_iwp
+       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
        pmvs = pmva - d_pmv
-       IF ( pmvs > -0.11_wp ) THEN
+       IF ( pmvs > -0.11_wp )  THEN
           d_pmv  = 0._wp
           pmvs   = -0.11_wp
@@ -2319,5 +2310,5 @@
 !     clo_fanger = clo
     clon = clo
-    IF ( clo > 0.5_wp .AND. perct_ij <= 8.73_wp ) THEN
+    IF ( clo > 0.5_wp  .AND.  perct_ij <= 8.73_wp )  THEN
 !
 !--    Required clothing insulation (ireq) is exclusively defined for 
@@ -2330,5 +2321,5 @@
     sultrieness    = .FALSE.
     d_std = -99._wp
-    IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
+    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
 !
 !--    Adjust for warm/humid conditions according to Gagge 1986
@@ -2337,9 +2328,9 @@
        pmvs   = pmva + d_pmv
        CALL perct_regression ( pmvs, clo, perct_ij )
-       IF ( sult_lim < 99._wp ) THEN
-          IF ( (perct_ij - ptc) > sult_lim ) sultrieness = .TRUE.
+       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 ( dgtcstd /= 0_iwp )  THEN
              d_std = ( ( perct_ij - ptc ) - dgtcm ) / dgtcstd
           ENDIF
@@ -2369,5 +2360,5 @@
 
 
-    IF ( ta < 0._wp ) THEN
+    IF ( ta < 0._wp )  THEN
 !
 !--    ta  < 0 (degC): water vapour pressure over ice
@@ -2451,16 +2442,16 @@
     clo_upper   = wclo
     y_upper     = pmv_w
-    IF ( ( y_lower > 0._wp .AND. y_upper < 0._wp ) .OR.                        &
-         ( y_lower < 0._wp .AND. y_upper > 0._wp ) ) THEN
+    IF ( ( y_lower > 0._wp  .AND.  y_upper < 0._wp )  .OR.                     &
+         ( y_lower < 0._wp  .AND.  y_upper > 0._wp ) )  THEN
        x_lower  = clo_lower
        x_upper  = clo_upper
        x_ridder = guess_0
 
-       DO j = 1_iwp, max_iteration
+       DO  j = 1_iwp, max_iteration
           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 )
-          IF ( sroot == 0._wp ) THEN
+          IF ( sroot == 0._wp )  THEN
              clo_res = x_average
              nerr = j
@@ -2469,5 +2460,5 @@
           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
@@ -2477,18 +2468,18 @@
           CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta,         &
                         y_new, top )
-          IF ( y_new == 0._wp ) THEN
+          IF ( y_new == 0._wp )  THEN
              clo_res = x_ridder
              nerr       = j
              RETURN
           ENDIF
-          IF ( SIGN ( y_average, y_new ) /= y_average ) THEN
+          IF ( SIGN ( y_average, y_new ) /= y_average )  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 ( SIGN ( y_lower, y_new ) /= y_lower )  THEN
              x_upper  = x_ridder
              y_upper  = y_new
-          ELSE IF ( SIGN ( y_upper, y_new ) /= y_upper ) THEN
+          ELSE IF ( SIGN ( y_upper, y_new ) /= y_upper )  THEN
              x_lower = x_ridder
              y_lower = y_new
@@ -2500,5 +2491,5 @@
              RETURN
           ENDIF
-          IF ( ABS ( x_upper - x_lower ) <= eps ) THEN
+          IF ( ABS ( x_upper - x_lower ) <= eps )  THEN
              clo_res = x_ridder
              nerr       = j
@@ -2511,7 +2502,7 @@
        clo_res = y_new
        RETURN
-    ELSE IF ( y_lower == 0. ) THEN
+    ELSE IF ( y_lower == 0. )  THEN
        x_ridder = clo_lower
-    ELSE IF ( y_upper == 0. ) THEN
+    ELSE IF ( y_upper == 0. )  THEN
        x_ridder = clo_upper
     ELSE
@@ -2541,8 +2532,8 @@
     REAL(wp), INTENT ( OUT ) ::  perct_ij   !< perct (degC) corresponding to given PMV / clo
 
-    IF ( pmv <= -0.11_wp ) THEN
+    IF ( pmv <= -0.11_wp )  THEN
        perct_ij = 5.805_wp + 12.6784_wp * pmv
     ELSE
-       IF ( pmv >= + 0.01_wp ) THEN
+       IF ( pmv >= + 0.01_wp )  THEN
           perct_ij = 16.826_wp + 6.163_wp * pmv
        ELSE
@@ -2608,5 +2599,5 @@
 !-- Clo must be > 0. to avoid div. by 0!
     clo = in_clo
-    IF ( clo <= 0._wp ) clo = .001_wp
+    IF ( clo <= 0._wp )  clo = .001_wp
 !
 !-- f_cl = Increase in surface due to clothing
@@ -2615,5 +2606,5 @@
 !-- Case of free convection (ws < 0.1 m/s ) not considered
     ws = in_ws
-    IF ( ws < .1_wp ) THEN
+    IF ( ws < .1_wp )  THEN
        ws = .1_wp
     ENDIF
@@ -2638,5 +2629,5 @@
 !-- Newton-approximation with air temperature as initial guess
     t_clothing = ta
-    DO i = 1, 3
+    DO  i = 1, 3
        t_clothing = t_clothing - ( ( t_clothing + degc_to_k )**4 + t_clothing  &
           * dc - gc ) / ( 4._wp * ( t_clothing + degc_to_k )**3 + dc )
@@ -2660,5 +2651,5 @@
     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
@@ -2760,5 +2751,5 @@
 !
 !-- Test for compliance with regression range
-    IF ( pmva < -1.0_wp .OR. pmva > 7.0_wp ) THEN
+    IF ( pmva < -1.0_wp  .OR.  pmva > 7.0_wp )  THEN
        nerr = -2_iwp
     ELSE
@@ -2772,9 +2763,9 @@
     p10  = 0.05_wp * svp_ta
     p95  = 1.00_wp * svp_ta
-    IF ( vp >= p10 .AND. vp <= p95 ) THEN
+    IF ( vp >= p10  .AND.  vp <= p95 )  THEN
        pa = vp
     ELSE
        nerr = -3_iwp
-       IF ( vp < p10 ) THEN
+       IF ( vp < p10 )  THEN
 !
 !--       Due to conditions of regression: r.H. >= 5 %
@@ -2786,5 +2777,5 @@
        ENDIF
     ENDIF
-    IF ( pa > 0._wp ) THEN
+    IF ( pa > 0._wp )  THEN
 !
 !--    Natural logarithm of pa
@@ -2796,5 +2787,5 @@
 !-- Ratio actual water vapour pressure to that of a r.H. of 50 %
     pa_p50   = 0.5_wp * svp_ta
-    IF ( pa_p50 > 0._wp .AND. pa > 0._wp ) THEN
+    IF ( pa_p50 > 0._wp  .AND.  pa > 0._wp )  THEN
        dapa   = apa - LOG ( pa_p50 )
        pa_p50 = pa / pa_p50
@@ -2805,5 +2796,5 @@
 !
 !-- Square root of wind velocity
-    IF ( ws >= 0._wp ) THEN
+    IF ( ws >= 0._wp )  THEN
        sqvel = SQRT ( ws )
     ELSE
@@ -2816,5 +2807,5 @@
 !-- Select the valid regression coefficients
     nreg = INT ( pmv )
-    IF ( nreg < 0_iwp ) THEN
+    IF ( nreg < 0_iwp )  THEN
 !
 !--    value of the FUNCTION in the case pmv <= -1
@@ -2823,11 +2814,11 @@
     ENDIF
     weight = MOD ( pmv, 1._wp )
-    IF ( weight < 0._wp ) weight = 0._wp
-    IF ( nreg > 5_iwp ) THEN
+    IF ( weight < 0._wp )  weight = 0._wp
+    IF ( nreg > 5_iwp )  THEN
        ! nreg=6
        nreg  = 5_iwp
        weight   = pmv - 5._wp
        weight2  = pmv - 6._wp
-       IF ( weight2 > 0_iwp ) THEN
+       IF ( weight2 > 0_iwp )  THEN
           weight = ( weight - weight2 ) / weight
        ENDIF
@@ -2836,26 +2827,26 @@
 !-- Regression valid for 0. <= pmv <= 6.
     dpmv_1 =                                                                   &
-       + bpa ( nreg ) * pa                                                     &
-       + bpmv ( nreg ) * pmv                                                   &
-       + bapa ( nreg ) * apa                                                   &
-       + bta ( nreg ) * ta                                                     &
-       + bdtmrt ( nreg ) * dtmrt                                               &
-       + bdapa ( nreg ) * dapa                                                 &
-       + bsqvel ( nreg ) * sqvel                                               &
-       + bpa_p50 ( nreg ) * pa_p50                                             &
-       + aconst ( nreg )
+       + bpa(nreg) * pa                                                        &
+       + bpmv(nreg) * pmv                                                      &
+       + bapa(nreg) * apa                                                      &
+       + bta(nreg) * ta                                                        &
+       + bdtmrt(nreg) * dtmrt                                                  &
+       + bdapa(nreg) * dapa                                                    &
+       + bsqvel(nreg) * sqvel                                                  &
+       + bpa_p50(nreg) * pa_p50                                                &
+       + aconst(nreg)
 
     dpmv_2 = 0._wp
-    IF ( nreg < 6_iwp ) THEN
+    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 )
+          + 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
 !
@@ -2863,9 +2854,9 @@
     deltapmv = ( 1._wp - weight ) * dpmv_1 + weight * dpmv_2
     pmvs = pmva + deltapmv
-    IF ( ( pmvs ) < 0._wp ) THEN
+    IF ( ( pmvs ) < 0._wp )  THEN
 !
 !--    Prevent negative pmv* due to problems with clothing insulation
        nerr = -4_iwp
-       IF ( pmvs > -0.11_wp ) THEN
+       IF ( pmvs > -0.11_wp )  THEN
 !
 !--       Threshold from FUNCTION perct_regression for winter clothing insulation
@@ -2925,5 +2916,5 @@
     dperctstd = 999999._wp
 
-    IF ( perct_ij < 16.826_wp .OR. perct_ij > 56._wp ) THEN
+    IF ( perct_ij < 16.826_wp  .OR.  perct_ij > 56._wp )  THEN
 !
 !--    Unallowed value of classical perct!
@@ -2941,5 +2932,5 @@
 !-- 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
@@ -2995,5 +2986,5 @@
     dtmrt          = tmrt - ta
     sqrt_ws        = ws
-    IF ( sqrt_ws < 0.10_wp ) THEN
+    IF ( sqrt_ws < 0.10_wp )  THEN
        sqrt_ws = 0.10_wp
     ELSE
@@ -3001,24 +2992,24 @@
     ENDIF
 
-    DO i = 1, 3
-       delta_cold (i) = 0._wp
-       IF ( i < 3 ) THEN
-          pmv_cross (i) = pmvc
+    DO  i = 1, 3
+       delta_cold(i) = 0._wp
+       IF ( i < 3 )  THEN
+          pmv_cross(i) = pmvc
        ENDIF
     ENDDO
 
-    DO i = 1, 3
+    DO  i = 1, 3
 !
 !--    Regression constant for given meteorological variables
-       reg_a(i) = coeff(i, 1) + coeff(i, 3) * ta + coeff(i, 4) *               &
+       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) * pmvc
     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
+    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
        ELSE
           nerr = 1_iwp
@@ -3028,6 +3019,6 @@
 
     i_bin = 3
-    DO i = 1, 2
-       IF ( pmva > pmv_cross (i) ) THEN
+    DO  i = 1, 2
+       IF ( pmva > pmv_cross(i) )  THEN
           i_bin = i
           EXIT
@@ -3062,12 +3053,12 @@
 !
 !--                                 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
+    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
+    ELSE IF ( pmva <= -1.28_wp )  THEN
        i_bin = 2_iwp
     ELSE
@@ -3075,10 +3066,10 @@
     ENDIF
 
-    dpmv_adj   = coef( i_bin, n_ca )
+    dpmv_adj   = coef(i_bin,n_ca)
     pmv        = 1._wp
 
-    DO i = n_ca + 1, n_ce
+    DO  i = n_ca + 1, n_ce
        pmv      = pmv * pmva
-       dpmv_adj = dpmv_adj + coef(i_bin, i) * pmv
+       dpmv_adj = dpmv_adj + coef(i_bin,i) * pmv
     ENDDO
     RETURN
@@ -3116,6 +3107,6 @@
 !
 !-- Regression only defined for perct <= 10 (degC)
-    IF ( ireq_neutral < 0.5_wp ) THEN
-       IF ( ireq_neutral < 0.48_wp ) THEN
+    IF ( ireq_neutral < 0.5_wp )  THEN
+       IF ( ireq_neutral < 0.48_wp )  THEN
           nerr = 1_iwp
        ENDIF
@@ -3125,5 +3116,5 @@
 !-- Minimal required clothing insulation: maximal acceptable body cooling
     ireq_minimal = 1.26_wp - 0.0588_wp * top02
-    IF ( nerr > 0_iwp ) THEN
+    IF ( nerr > 0_iwp )  THEN
        ireq_minimal = ireq_neutral
     ENDIF
@@ -3152,12 +3143,12 @@
 !
 !-- According to Gehan-George, for children
-    IF ( age < 19_iwp ) THEN
-       IF ( age < 5_iwp ) THEN
+    IF ( age < 19_iwp )  THEN
+       IF ( age < 5_iwp )  THEN
           surf = 0.02667_wp * height**0.42246_wp * weight**0.51456_wp
           RETURN
-       END IF
+       ENDIF
        surf = 0.03050_wp * height**0.35129_wp * weight**0.54375_wp
        RETURN
-    END IF
+    ENDIF
 !
 !-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if 
@@ -3200,14 +3191,14 @@
 
     CALL surface_area ( height, weight, INT( age ), a_surf )
-    s = height * 100._wp / ( weight ** ( 1._wp / 3._wp ) )
+    s = height * 100._wp / ( weight**( 1._wp / 3._wp ) )
     factor = 1._wp + .004_wp  * ( 30._wp - age )
     basic_heat_prod = 0.
-    IF ( sex == 1_iwp ) THEN
-       basic_heat_prod = 3.45_wp * weight ** ( 3._wp / 4._wp ) * ( factor +    &
+    IF ( sex == 1_iwp )  THEN
+       basic_heat_prod = 3.45_wp * weight**( 3._wp / 4._wp ) * ( factor +      &
                      .01_wp  * ( s - 43.4_wp ) )
-    ELSE IF ( sex == 2_iwp ) THEN
-       basic_heat_prod = 3.19_wp * weight ** ( 3._wp / 4._wp ) * ( factor +    &
+    ELSE IF ( sex == 2_iwp )  THEN
+       basic_heat_prod = 3.19_wp * weight**( 3._wp / 4._wp ) * ( factor +      &
                     .018_wp * ( s - 42.1_wp ) )
-    END IF
+    ENDIF
 
     energy_prod = work + basic_heat_prod
@@ -3295,5 +3286,5 @@
 !
 !-- Decision: firstly calculate for winter or summer clothing
-    IF ( ta <= 10._wp ) THEN
+    IF ( ta <= 10._wp )  THEN
 !
 !--    First guess: winter clothing insulation: cold stress
@@ -3304,5 +3295,5 @@
        pmv_w = pmva
 
-       IF ( pmva > 0._wp ) THEN
+       IF ( pmva > 0._wp )  THEN
 !
 !--       Case summer clothing insulation: heat load ?            
@@ -3312,5 +3303,5 @@
              t_clothing, storage, dt, pmva )
           pmv_s = pmva
-          IF ( pmva <= 0._wp ) THEN
+          IF ( pmva <= 0._wp )  THEN
 !
 !--          Case: comfort achievable by varying clothing insulation 
@@ -3318,9 +3309,9 @@
              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta , sclo,     &
                             pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
-             IF ( ncount < 0_iwp ) THEN
+             IF ( ncount < 0_iwp )  THEN
                 nerr = -1_iwp
                 RETURN
              ENDIF
-          ELSE IF ( pmva > 0.06_wp ) THEN
+          ELSE IF ( pmva > 0.06_wp )  THEN
              clo = 0.5_wp
              t_clothing = bio_fill_value
@@ -3328,5 +3319,5 @@
                 t_clothing, storage, dt, pmva )
           ENDIF
-       ELSE IF ( pmva < -0.11_wp ) THEN
+       ELSE IF ( pmva < -0.11_wp )  THEN
           clo = 1.75_wp
           t_clothing = bio_fill_value
@@ -3344,5 +3335,5 @@
        pmv_s = pmva
 
-       IF ( pmva < 0._wp ) THEN
+       IF ( pmva < 0._wp )  THEN
 !
 !--       Case winter clothing insulation: cold stress ?
@@ -3353,5 +3344,5 @@
           pmv_w = pmva
 
-          IF ( pmva >= 0._wp ) THEN
+          IF ( pmva >= 0._wp )  THEN
 !
 !--          Case: comfort achievable by varying clothing insulation 
@@ -3359,9 +3350,9 @@
              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
                                pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
-             IF ( ncount < 0_wp ) THEN
+             IF ( ncount < 0_wp )  THEN
                 nerr = -1_iwp
                 RETURN
              ENDIF
-          ELSE IF ( pmva < -0.11_wp ) THEN
+          ELSE IF ( pmva < -0.11_wp )  THEN
              clo = 1.75_wp
              t_clothing = bio_fill_value
@@ -3369,5 +3360,5 @@
                 t_clothing, storage, dt, pmva )
           ENDIF
-       ELSE IF ( pmva > 0.06_wp ) THEN
+       ELSE IF ( pmva > 0.06_wp )  THEN
           clo = 0.5_wp
           t_clothing = bio_fill_value
@@ -3382,11 +3373,11 @@
     CALL perct_regression ( pmva, clo, ipt )
     ptc = ipt
-    IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
+    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
 !
 !--    Adjust for cold conditions according to Gagge 1986
        CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
-       IF ( nerr_cold > 0_iwp ) nerr = -5_iwp
+       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
        pmvs = pmva - d_pmv
-       IF ( pmvs > -0.11_wp ) THEN
+       IF ( pmvs > -0.11_wp )  THEN
           d_pmv  = 0._wp
           pmvs   = -0.11_wp
@@ -3396,5 +3387,5 @@
 !     clo_fanger = clo
     clon = clo
-    IF ( clo > 0.5_wp .AND. ipt <= 8.73_wp ) THEN
+    IF ( clo > 0.5_wp  .AND.  ipt <= 8.73_wp )  THEN
 !
 !--    Required clothing insulation (ireq) is exclusively defined for 
@@ -3407,5 +3398,5 @@
     sultrieness    = .FALSE.
     d_std      = -99._wp
-    IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
+    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
 !
 !--    Adjust for warm/humid conditions according to Gagge 1986
@@ -3414,6 +3405,6 @@
        pmvs   = pmva + d_pmv
        CALL perct_regression ( pmvs, clo, ipt )
-       IF ( sult_lim < 99._wp ) THEN
-          IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE.
+       IF ( sult_lim < 99._wp )  THEN
+          IF ( (ipt - ptc) > sult_lim )  sultrieness = .TRUE.
        ENDIF
     ENDIF
@@ -3487,11 +3478,11 @@
     CALL perct_regression ( pmva, clo, ipt )
 
-    IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
+    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
 !
 !--    Adjust for cold conditions according to Gagge 1986
        CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
-       IF ( nerr_cold > 0_iwp ) nerr = -5_iwp
+       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
        pmvs = pmva - d_pmv
-       IF ( pmvs > -0.11_wp ) THEN
+       IF ( pmvs > -0.11_wp )  THEN
           d_pmv  = 0._wp
           pmvs   = -0.11_wp
@@ -3505,5 +3496,5 @@
     sultrieness    = .FALSE.
     d_std      = -99._wp
-    IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
+    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
 !
 !--    Adjust for warm/humid conditions according to Gagge 1986
@@ -3512,6 +3503,6 @@
        pmvs   = pmva + d_pmv
        CALL perct_regression ( pmvs, clo, ipt )
-       IF ( sult_lim < 99._wp ) THEN
-          IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE.
+       IF ( sult_lim < 99._wp )  THEN
+          IF ( (ipt - ptc) > sult_lim )  sultrieness = .TRUE.
        ENDIF
     ENDIF
@@ -3583,5 +3574,5 @@
 !-- Clo must be > 0. to avoid div. by 0!
     clo = in_clo
-    IF ( clo < 001._wp ) clo = .001_wp
+    IF ( clo < 001._wp )  clo = .001_wp
 !
 !-- Increase in surface due to clothing
@@ -3590,5 +3581,5 @@
 !-- Case of free convection (ws < 0.1 m/s ) not considered
     ws = in_ws
-    IF ( ws < .1_wp ) THEN
+    IF ( ws < .1_wp )  THEN
        ws = .1_wp
     ENDIF
@@ -3610,11 +3601,11 @@
 !-- newton-approximation with air temperature as initial guess
     niter = INT( dt * 10._wp, KIND=iwp )
-    IF ( niter < 1 ) niter = 1_iwp
+    IF ( niter < 1 )  niter = 1_iwp
     adjustrate = 1._wp - EXP ( -1._wp * ( 10._wp / time_equil ) * dt )
-    IF ( adjustrate >= 1._wp ) adjustrate = 1._wp
+    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
+    IF ( t_cloth <= -998.0_wp )  THEN  ! If initial run
        niter = 3_iwp
        adjustrate = 1._wp
@@ -3623,5 +3614,5 @@
     ENDIF
 
-    DO i = 1, niter
+    DO  i = 1, niter
        t_clothing = t_clothing - adjustrate_cloth * ( ( t_clothing +           &
           273.2_wp )**4._wp + t_clothing *                                     &
@@ -3775,19 +3766,19 @@
     REAL(wp) ::  vpex
 
-    met = 3.45_wp * mbody ** ( 3._wp / 4._wp ) * (1._wp + 0.004_wp *           &
+!
+!-- Metabolic heat production
+    met = 3.45_wp * mbody**( 3._wp / 4._wp ) * (1._wp + 0.004_wp *             &
           ( 30._wp - age) + 0.010_wp * ( ( ht * 100._wp /                      &
-          ( mbody ** ( 1._wp / 3._wp ) ) ) - 43.4_wp ) )
-
+          ( mbody**( 1._wp / 3._wp ) ) ) - 43.4_wp ) )
     met = work + met
-
     int_heat = met * (1._wp - eta)
 !
 !-- Sensible respiration energy
     tex  = 0.47_wp * ta + 21.0_wp
-    rtv  = 1.44_wp * 10._wp ** (-6._wp) * met
+    rtv  = 1.44_wp * 10._wp**(-6._wp) * met
     eres = c_p * (ta - tex) * rtv
 !
 !-- Latent respiration energy
-    vpex = 6.11_wp * 10._wp ** ( 7.45_wp * tex / ( 235._wp + tex ) )
+    vpex = 6.11_wp * 10._wp**( 7.45_wp * tex / ( 235._wp + tex ) )
     erel = 0.623_wp * l_v / pair * ( vpa - vpex ) * rtv
 !
@@ -3889,28 +3880,36 @@
     LOGICAL ::  skipIncreaseCount   !< iteration control flag
 
-    wetsk = 0._wp
-    adu = 0.203_wp * mbody ** 0.425_wp * ht ** 0.725_wp
-
-    hc = 2.67_wp + ( 6.5_wp * v ** 0.67_wp)
-    hc   = hc * (pair /po) ** 0.55_wp
+!
+!-- Initialize
+    wetsk = 0._wp  !< skin is dry everywhere on init (no non-evaporated sweat)
+!
+!-- Set Du Bois Area for the sample person
+    adu = 0.203_wp * mbody**0.425_wp * ht**0.725_wp
+!
+!-- Initialize convective heat considering local air preassure
+    hc = 2.67_wp + ( 6.5_wp * v**0.67_wp )
+    hc = hc * ( pair / po )**0.55_wp
+!
+!-- Set surface modification by posture (the person will always stand)
     feff = 0.725_wp                     !< Posture: 0.725 for stading
-    facl = (- 2.36_wp + 173.51_wp * clo - 100.76_wp * clo * clo + 19.28_wp     &
-          * (clo ** 3._wp)) / 100._wp
-
-    IF ( facl > 1._wp )   facl = 1._wp
+!
+!-- Set surface modification by clothing
+    facl = ( - 2.36_wp + 173.51_wp * clo - 100.76_wp * clo * clo + 19.28_wp    &
+          * ( clo**3._wp ) ) / 100._wp
+    IF ( facl > 1._wp )  facl = 1._wp
+!
+!-- Initialize heat resistences
     rcl = ( clo / 6.45_wp ) / facl
     IF ( clo >= 2._wp )  y  = 1._wp
-
-    IF ( ( clo > 0.6_wp ) .AND. ( clo < 2._wp ) )  y = ( ht - 0.2_wp ) / ht
-    IF ( ( clo <= 0.6_wp ) .AND. ( clo > 0.3_wp ) ) y = 0.5_wp
-    IF ( ( clo <= 0.3_wp ) .AND. ( clo > 0._wp ) )  y = 0.1_wp
-
-    r2   = adu * (fcl - 1._wp + facl) / (2._wp * 3.14_wp * ht * y)
-    r1   = facl * adu / (2._wp * 3.14_wp * ht * y)
-
+    IF ( ( clo > 0.6_wp )   .AND.  ( clo < 2._wp ) )   y = ( ht - 0.2_wp ) / ht
+    IF ( ( clo <= 0.6_wp )  .AND.  ( clo > 0.3_wp ) )  y = 0.5_wp
+    IF ( ( clo <= 0.3_wp )  .AND.  ( clo > 0._wp ) )   y = 0.1_wp
+    r2   = adu * ( fcl - 1._wp + facl ) / ( 2._wp * 3.14_wp * ht * y )
+    r1   = facl * adu / ( 2._wp * 3.14_wp * ht * y )
     di   = r2 - r1
-!
-!-- Skin temperatur
-    DO j = 1, 7
+
+!
+!-- Estimate skin temperatur
+    DO  j = 1, 7
 
        tsk    = 34._wp
@@ -3920,8 +3919,8 @@
        enbal2 = 0._wp
 
-       DO i = 1, 100  ! allow for 100 iterations max
+       DO  i = 1, 100  ! allow for 100 iterations max
           acl   = adu * facl + adu * ( fcl - 1._wp )
           rclo2 = emcl * sigma_sb * ( ( tcl + degc_to_k )**4._wp -             &
-            ( tmrt + degc_to_k )** 4._wp ) * feff
+            ( tmrt + degc_to_k )**4._wp ) * feff
           htcl  = 6.28_wp * ht * y * di / ( rcl * LOG( r2 / r1 ) * acl )
           tsk   = 1._wp / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl
@@ -3930,7 +3929,7 @@
           aeff  = adu * feff
           rbare = aeff * ( 1._wp - facl ) * emsk * sigma_sb *                  &
-            ( ( tmrt + degc_to_k )** 4._wp - ( tsk + degc_to_k )** 4._wp )
+            ( ( tmrt + degc_to_k )**4._wp - ( tsk + degc_to_k )**4._wp )
           rclo  = feff * acl * emcl * sigma_sb *                               &
-            ( ( tmrt + degc_to_k )** 4._wp - ( tcl + degc_to_k )** 4._wp )
+            ( ( tmrt + degc_to_k )**4._wp - ( tcl + degc_to_k )**4._wp )
           rsum  = rbare + rclo
 !
@@ -3952,20 +3951,20 @@
           c(9)  = 5.28_wp * adu - c(5) - c(4) * tsk
           c(10) = c(9) * c(9) - 4._wp * c(4) *                                 &
-             ( c(5) * tsk - c(0) - 5.28_wp * adu * tsk )
-
-          IF ( tsk == 36._wp ) tsk = 36.01_wp
+                  ( c(5) * tsk - c(0) - 5.28_wp * adu * tsk )
+
+          IF ( tsk == 36._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 ) /   &
             ( 1._wp + 0.5_wp * ( 34._wp - tsk ) ) ) + tsk
-          IF ( c(10) >= 0._wp ) THEN
+          IF ( c(10) >= 0._wp )  THEN
              tcore(6) = ( - c(9) - c(10)**0.5_wp ) / ( 2._wp * c(4) )
              tcore(1) = ( - c(9) + c(10)**0.5_wp ) / ( 2._wp * c(4) )
-          END IF
-
-          IF ( c(8) >= 0._wp ) THEN
-             tcore(2) = ( - c(6) + ABS( c(8) ) ** 0.5_wp ) / ( 2._wp * c(4) )
-             tcore(5) = ( - c(6) - ABS( c(8) ) ** 0.5_wp ) / ( 2._wp * c(4) )
+          ENDIF
+
+          IF ( c(8) >= 0._wp )  THEN
+             tcore(2) = ( - c(6) + ABS( c(8) )**0.5_wp ) / ( 2._wp * c(4) )
+             tcore(5) = ( - c(6) - ABS( c(8) )**0.5_wp ) / ( 2._wp * c(4) )
              tcore(4) = c(0) / ( 5.28_wp * adu + c(1) * 1._wp / 40._wp ) + tsk
-          END IF
+          ENDIF
 !
 !--       Transpiration
@@ -3974,5 +3973,5 @@
           vpts  = 6.11_wp * 10._wp**( 7.45_wp * tsk / ( 235._wp + tsk ) )
 
-          IF ( tbody <= 36.6_wp ) swm = 0._wp  !< no need for sweating
+          IF ( tbody <= 36.6_wp )  swm = 0._wp  !< no need for sweating
 
           sw = swm
@@ -3983,15 +3982,15 @@
           wetsk  = eswphy / eswpot
 
-          IF ( wetsk > 1._wp ) wetsk = 1._wp
+          IF ( wetsk > 1._wp )  wetsk = 1._wp
 !
 !--       Sweat production > evaporation?
           eswdif = eswphy - eswpot
 
-          IF ( eswdif <= 0._wp ) esw = eswpot  !< Limit is evaporation
-          IF ( eswdif > 0._wp ) esw = eswphy   !< Limit is sweat production
-          IF ( esw  > 0._wp )   esw = 0._wp    !< Sweat can't be evaporated, no more cooling effect
+          IF ( eswdif <= 0._wp )  esw = eswpot  !< Limit is evaporation
+          IF ( eswdif > 0._wp )   esw = eswphy  !< Limit is sweat production
+          IF ( esw  > 0._wp )     esw = 0._wp   !< Sweat can't be evaporated, no more cooling effect
 !
 !--       Diffusion
-          rdsk = 0.79_wp * 10._wp ** 7._wp
+          rdsk = 0.79_wp * 10._wp**7._wp
           rdcl = 0._wp
           ed   = l_v / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( vpa -     &
@@ -4002,6 +4001,6 @@
           vb2 = tcore(j) - 36.6_wp
 
-          IF ( vb2 < 0._wp ) vb2 = 0._wp
-          IF ( vb1 < 0._wp ) vb1 = 0._wp
+          IF ( vb2 < 0._wp )  vb2 = 0._wp
+          IF ( vb1 < 0._wp )  vb1 = 0._wp
           vb = ( 6.3_wp + 75._wp * vb2 ) / ( 1._wp + 0.5_wp * vb1 )
 !
@@ -4011,57 +4010,57 @@
 !--       Clothing temperature
           xx = 0.001_wp
-          IF ( count1 == 0_iwp ) xx = 1._wp
-          IF ( count1 == 1_iwp ) xx = 0.1_wp
-          IF ( count1 == 2_iwp ) xx = 0.01_wp
-          IF ( count1 == 3_iwp ) xx = 0.001_wp
-
-          IF ( enbal > 0._wp ) tcl = tcl + xx
-          IF ( enbal < 0._wp ) tcl = tcl - xx
+          IF ( count1 == 0_iwp )  xx = 1._wp
+          IF ( count1 == 1_iwp )  xx = 0.1_wp
+          IF ( count1 == 2_iwp )  xx = 0.01_wp
+          IF ( count1 == 3_iwp )  xx = 0.001_wp
+
+          IF ( enbal > 0._wp )  tcl = tcl + xx
+          IF ( enbal < 0._wp )  tcl = tcl - xx
 
           skipIncreaseCount = .FALSE.
-          IF ( ( (enbal <= 0._wp ) .AND. (enbal2 > 0._wp ) ) .OR.              &
-             ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) ) THEN
+          IF ( ( (enbal <= 0._wp )  .AND.  (enbal2 > 0._wp ) )  .OR.           &
+             ( ( enbal >= 0._wp )  .AND.  ( enbal2 < 0._wp ) ) )  THEN
              skipIncreaseCount = .TRUE.
           ELSE
              enbal2 = enbal
              count3 = count3 + 1_iwp
-          END IF
-
-          IF ( ( count3 > 200_iwp ) .OR. skipIncreaseCount ) THEN
-             IF ( count1 < 3_iwp ) THEN
+          ENDIF
+
+          IF ( ( count3 > 200_iwp )  .OR.  skipIncreaseCount )  THEN
+             IF ( count1 < 3_iwp )  THEN
                 count1 = count1 + 1_iwp
                 enbal2 = 0._wp
              ELSE
                 EXIT
-             END IF
-          END IF
-       END DO
-
-       IF ( count1 == 3_iwp ) THEN
+             ENDIF
+          ENDIF
+       ENDDO
+
+       IF ( count1 == 3_iwp )  THEN
           SELECT CASE ( j )
              CASE ( 2, 5) 
-                IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND.                     &
-                   ( tsk <= 34.050_wp ) ) ) CYCLE
+                IF ( .NOT. ( ( tcore(j) >= 36.6_wp )  .AND.                    &
+                   ( tsk <= 34.050_wp ) ) )  CYCLE
              CASE ( 6, 1 )
                 IF ( c(10) < 0._wp ) CYCLE
-                IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND.                     &
-                   ( tsk > 33.850_wp ) ) ) CYCLE
+                IF ( .NOT. ( ( tcore(j) >= 36.6_wp )  .AND.                    &
+                   ( tsk > 33.850_wp ) ) )  CYCLE
              CASE ( 3 )
-                IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND.                      &
-                   ( tsk <= 34.000_wp ) ) ) CYCLE
+                IF ( .NOT. ( ( tcore(j) < 36.6_wp )  .AND.                     &
+                   ( tsk <= 34.000_wp ) ) )  CYCLE
              CASE ( 7 )
-                IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND.                      &
-                   ( tsk > 34.000_wp ) ) ) CYCLE
+                IF ( .NOT. ( ( tcore(j) < 36.6_wp )  .AND.                     &
+                   ( tsk > 34.000_wp ) ) )  CYCLE
              CASE default
           END SELECT
-       END IF
-
-       IF ( ( j /= 4_iwp ) .AND. ( vb >= 91._wp ) ) CYCLE
-       IF ( ( j == 4_iwp ) .AND. ( vb < 89._wp ) ) CYCLE
+       ENDIF
+
+       IF ( ( j /= 4_iwp )  .AND.  ( vb >= 91._wp ) )  CYCLE
+       IF ( ( j == 4_iwp )  .AND.  ( vb < 89._wp ) )  CYCLE
        IF ( vb > 90._wp ) vb = 90._wp
 !
 !--    Loses by water
        ws = sw * 3600._wp * 1000._wp
-       IF ( ws > 2000._wp ) ws = 2000._wp
+       IF ( ws > 2000._wp )  ws = 2000._wp
        wd = ed / l_v * 3600._wp * ( -1000._wp )
        wr = erel / l_v * 3600._wp * ( -1000._wp )
@@ -4070,5 +4069,5 @@
 
        RETURN
-    END DO
+    ENDDO
  END SUBROUTINE heat_exch
 
@@ -4131,15 +4130,15 @@
     enbal2 = 0._wp
 
-    DO count1 = 0, 3
-       DO i = 1, 125  ! 500 / 4
-          hc = 2.67_wp + 6.5_wp * 0.1_wp ** 0.67_wp
-          hc = hc * ( pair / po ) ** 0.55_wp
+    DO  count1 = 0, 3
+       DO  i = 1, 125  ! 500 / 4
+          hc = 2.67_wp + 6.5_wp * 0.1_wp**0.67_wp
+          hc = hc * ( pair / po )**0.55_wp
 !
 !--       Radiation
           aeff  = adu * feff
           rbare = aeff * ( 1._wp - facl ) * emsk * sigma_sb *                  &
-              ( ( pet_ij + degc_to_k ) ** 4._wp - ( tsk + degc_to_k ) ** 4._wp )
+              ( ( pet_ij + degc_to_k )**4._wp - ( tsk + degc_to_k )**4._wp )
           rclo  = feff * acl * emcl * sigma_sb *                               &
-              ( ( pet_ij + degc_to_k ) ** 4._wp - ( tcl + degc_to_k ) ** 4._wp )
+              ( ( pet_ij + degc_to_k )**4._wp - ( tcl + degc_to_k )**4._wp )
           rsum  = rbare + rclo
 !
@@ -4156,5 +4155,5 @@
           tex  = 0.47_wp * pet_ij + 21._wp
           eres = c_p * ( pet_ij - tex ) * rtv
-          vpex = 6.11_wp * 10._wp ** ( 7.45_wp * tex / ( 235._wp + tex ) )
+          vpex = 6.11_wp * 10._wp**( 7.45_wp * tex / ( 235._wp + tex ) )
           erel = 0.623_wp * l_v / pair * ( 12._wp - vpex ) * rtv
           ere  = eres + erel
@@ -4171,10 +4170,10 @@
           IF ( enbal > 0._wp )  pet_ij = pet_ij - xx
           IF ( enbal < 0._wp )  pet_ij = pet_ij + xx
-          IF ( ( enbal <= 0._wp ) .AND. ( enbal2 > 0._wp ) ) EXIT
-          IF ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) EXIT
+          IF ( ( enbal <= 0._wp )  .AND.  ( enbal2 > 0._wp ) )  EXIT
+          IF ( ( enbal >= 0._wp )  .AND.  ( enbal2 < 0._wp ) )  EXIT
 
           enbal2 = enbal
-       END DO
-    END DO
+       ENDDO
+    ENDDO
  END SUBROUTINE pet_iteration
 
@@ -4271,5 +4270,5 @@
     CALL uvem_solar_position
      
-    IF (sza  .GE.  90) THEN
+    IF (sza  .GE.  90)  THEN
        vitd3_exposure(:,:) = 0.0_wp
     ELSE
@@ -4368,10 +4367,10 @@
 !                   
 ! !--        extract obstruction from IBSET-Integer_Array ------------------'
-             IF (consider_obstructions ) THEN
+             IF (consider_obstructions )  THEN
                 obstruction_temp1 = building_obstruction_f%var_3d(:,j,i)
-                IF (obstruction_temp1(0)  .NE.  9) THEN 
+                IF (obstruction_temp1(0)  .NE.  9)  THEN 
                    DO  pobi = 0, 44 
                       DO  bi = 0, 7 
-                         IF ( btest( obstruction_temp1(pobi), bi )  .EQV.  .TRUE.) THEN
+                         IF ( btest( obstruction_temp1(pobi), bi )  .EQV.  .TRUE.)  THEN
                             obstruction_temp2( ( pobi * 8 ) + bi ) = 1
                          ELSE
@@ -4392,5 +4391,5 @@
          
              obstruction_direct_beam = obstruction( nint(startpos_saa_float), nint( sza / 10.0_wp ) ) 
-             IF (sza  .GE.  89.99_wp) THEN
+             IF (sza  .GE.  89.99_wp)  THEN
                 sza = 89.99999_wp
              ENDIF
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 3741)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 3742)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! - Moved call of bio_calculate_thermal_index_maps from biometeorology module to
+! time_integration to make sure averaged input is updated before calculating.
+! 
+! 3739 2019-02-13 08:05:17Z dom_dwd_user
 ! Removed everything related to "time_bio_results" as this is never used.
 ! 
@@ -1635,5 +1639,5 @@
        IF ( time_dopr_listing >= dt_dopr_listing )  THEN
           CALL print_1d
-          time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing, &
+          time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing,    &
                                                            dt_3d ) )
        ENDIF
@@ -1641,5 +1645,5 @@
 !
 !--    Graphic output for PROFIL
-       IF (        time_dopr >= dt_dopr                                     &
+       IF (        time_dopr >= dt_dopr                                        &
             .AND.  time_since_reference_point >= skip_time_dopr )  THEN
           IF ( dopr_n /= 0 )  CALL data_output_profiles
@@ -1658,5 +1662,5 @@
 !--    Output of spectra (formatted for use with PROFIL), in case of no
 !--    time averaging, spectra has to be calculated before
-       IF (         time_dosp >= dt_dosp                                     &
+       IF (         time_dosp >= dt_dosp                                       &
              .AND.  time_since_reference_point >= skip_time_dosp )  THEN
           IF ( average_count_sp == 0 )  CALL calc_spectra
@@ -1667,15 +1671,15 @@
 !
 !--    2d-data output (cross-sections)
-       IF (         time_do2d_xy >= dt_do2d_xy                               &
+       IF (         time_do2d_xy >= dt_do2d_xy                                 &
              .AND.  time_since_reference_point >= skip_time_do2d_xy )  THEN
           CALL data_output_2d( 'xy', 0 )
           time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
        ENDIF
-       IF (         time_do2d_xz >= dt_do2d_xz                               &
+       IF (         time_do2d_xz >= dt_do2d_xz                                 &
              .AND.  time_since_reference_point >= skip_time_do2d_xz )  THEN
           CALL data_output_2d( 'xz', 0 )
           time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) )
        ENDIF
-       IF (         time_do2d_yz >= dt_do2d_yz                               &
+       IF (         time_do2d_yz >= dt_do2d_yz                                 &
              .AND.  time_since_reference_point >= skip_time_do2d_yz )  THEN
           CALL data_output_2d( 'yz', 0 )
@@ -1685,5 +1689,5 @@
 !
 !--    3d-data output (volume data)
-       IF (         time_do3d >= dt_do3d                                     & 
+       IF (         time_do3d >= dt_do3d                                       &
              .AND.  time_since_reference_point >= skip_time_do3d )  THEN
           CALL data_output_3d( 0 )
@@ -1694,8 +1698,8 @@
 !--    Masked data output
        DO  mid = 1, masks
-          IF (         time_domask(mid) >= dt_domask(mid)                          &
+          IF (         time_domask(mid) >= dt_domask(mid)                      &
                 .AND.  time_since_reference_point >= skip_time_domask(mid) )  THEN
              CALL data_output_mask( 0 )
-             time_domask(mid) = MOD( time_domask(mid),                             &
+             time_domask(mid) = MOD( time_domask(mid),                         &
                                      MAX( dt_domask(mid), dt_3d ) )
           ENDIF
@@ -1704,7 +1708,12 @@
 !
 !--    Output of time-averaged 2d/3d/masked data
-       IF (         time_do_av >= dt_data_output_av                                &
+       IF (         time_do_av >= dt_data_output_av                            &
              .AND.  time_since_reference_point >= skip_time_data_output_av )  THEN
           CALL average_3d_data
+!
+!--       Udate thermal comfort indices based on updated averaged input
+          IF ( biometeorology  .AND.  thermal_comfort )  THEN
+             CALL bio_calculate_thermal_index_maps ( .TRUE. )
+          ENDIF
           CALL data_output_2d( 'xy', 1 )
           CALL data_output_2d( 'xz', 1 )
@@ -1719,10 +1728,10 @@
 !--    Output of surface data, instantaneous and averaged data
        IF ( surface_output )  THEN
-          IF (         time_dosurf >= dt_dosurf                                    &
+          IF (         time_dosurf >= dt_dosurf                                &
                 .AND.  time_since_reference_point >= skip_time_dosurf )  THEN
              CALL surface_data_output( 0 )
              time_dosurf = MOD( time_dosurf, MAX( dt_dosurf, dt_3d ) )
           ENDIF
-          IF (         time_dosurf_av >= dt_dosurf_av                              &
+          IF (         time_dosurf_av >= dt_dosurf_av                          &
                 .AND.  time_since_reference_point >= skip_time_dosurf_av )  THEN
              CALL surface_data_output( 1 )
