Index: palm/trunk/SOURCE/biometeorology_mod.f90
===================================================================
--- palm/trunk/SOURCE/biometeorology_mod.f90	(revision 3739)
+++ palm/trunk/SOURCE/biometeorology_mod.f90	(revision 3740)
@@ -27,4 +27,8 @@
 ! -----------------
 ! $Id$
+! - Added safety-meassure to catch the case that 'bio_mrt_av' is stated after
+! 'bio_<index>' in the output section of the p3d file.
+! 
+! 3739 2019-02-13 08:05:17Z dom_dwd_user
 ! - Auto-adjusting thermal_comfort flag if not set by user, but thermal_indices
 ! set as output quantities.
@@ -395,8 +399,16 @@
 
           CASE ( 'bio_mrt' )
-                IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
+!
+!--             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
                    ALLOCATE( mrt_av_grid(nmrtbl) )
                 ENDIF
                 mrt_av_grid = 0.0_wp
+
 
           CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
@@ -410,5 +422,5 @@
 !--          Indices are in unknown order as depending on the input file,
 !--          determine first index to average und update only once
-
+!
 !--          Only proceed here if this was not done for any index before. This
 !--          is done only once during the whole model run.
@@ -430,5 +442,5 @@
 !
 !--          Only continue if var is the index, that controls averaging.
-!            Break immediatelly (doing nothing) for the other indices.
+!--          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
@@ -436,5 +448,5 @@
 !
 !--          Now memorize which of the input grids are not averaged by other 
-!            modules. Set averaging switch to .TRUE. in that case.
+!--          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) )
@@ -489,5 +501,84 @@
 
           CASE ( 'bio_mrt' )
-             IF ( ALLOCATED( mrt_av_grid ) )  THEN
+!
+!--          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. 
+             IF ( .NOT. do_average_mrt .AND. ALLOCATED( mrt_av_grid ) )  THEN
+
+                IF ( mrt_include_sw )  THEN
+                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
+                      ( ( human_absorb * mrtinsw(:) +                          &
+                      mrtinlw(:) ) /                                           &
+                      ( 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
+                ENDIF
+             ENDIF
+
+          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
+!
+!--          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
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+             IF ( ALLOCATED( q_av )  .AND. do_average_q ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb, nzt+1
+                         q_av(k,j,i) = q_av(k,j,i) + q(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+             IF ( ALLOCATED( u_av )  .AND. do_average_u ) THEN
+                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
+                   DO  j = nysg, nyng
+                      DO  k = nzb, nzt+1
+                         u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+             IF ( ALLOCATED( v_av )  .AND. do_average_v ) THEN
+                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
+                   DO  j = nysg, nyng
+                      DO  k = nzb, nzt+1
+                         v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+             IF ( ALLOCATED( w_av )  .AND. do_average_w ) THEN
+                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
+                   DO  j = nysg, nyng
+                      DO  k = nzb, nzt+1
+                         w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+
+             IF ( ALLOCATED( mrt_av_grid ) .AND. do_average_mrt )  THEN
 
                 IF ( mrt_include_sw )  THEN
@@ -502,78 +593,4 @@
                 ENDIF
              ENDIF
-
-          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
-!
-!--          Only continue if updateindex, 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
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      DO  k = nzb, nzt+1
-                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-             IF ( ALLOCATED( q_av )  .AND. do_average_q ) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      DO  k = nzb, nzt+1
-                         q_av(k,j,i) = q_av(k,j,i) + q(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-             IF ( ALLOCATED( u_av )  .AND. do_average_u ) THEN
-                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-             IF ( ALLOCATED( v_av )  .AND. do_average_v ) THEN
-                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-             IF ( ALLOCATED( w_av )  .AND. do_average_w ) THEN
-                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
-                   DO  j = nysg, nyng
-                      DO  k = nzb, nzt+1
-                         w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
-                      ENDDO
-                   ENDDO
-                ENDDO
-             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
-                ELSE
-                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
-                      ( mrtinlw(:) /                             &  ! ( human_emiss * mrtinlw(:) /
-                      ( human_emiss * sigma_sb ) ) ** .25_wp - degc_to_k
-                ENDIF
-             ENDIF
 !
 !--       This is a cumulated dose. No mode == 'average' for this quantity.
@@ -597,5 +614,9 @@
 
           CASE ( 'bio_mrt' )
-             IF ( ALLOCATED( mrt_av_grid ) )  THEN
+!
+!--          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. 
+             IF ( .NOT. do_average_mrt .AND. ALLOCATED( mrt_av_grid ) )  THEN
                 mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp )
              ENDIF
@@ -815,8 +836,10 @@
     END SELECT
 
+!
+!--    Further checks if thermal comfort output is desired.
     IF ( thermal_comfort .AND. unit == 'degree_C' )  THEN
-!
-!--    Futher checks if output belongs to biometeorology
-!      Break if required modules "radiation" and "humidity" are not running.
+
+!
+!--    Break if required modules "radiation" and "humidity" are not avalable.
        IF ( .NOT.  radiation ) THEN
           message_string = 'output of "' // TRIM( var ) // '" require'         &
@@ -848,5 +871,5 @@
 ! ------------
 !> Check parameters routine for biom module
-!> check_parameters 1302
+!> Currently unused but might come in handy for future checks?
 !------------------------------------------------------------------------------!
  SUBROUTINE bio_check_parameters
@@ -854,5 +877,4 @@
     USE control_parameters,                                                    &
         ONLY:  message_string
-
 
     IMPLICIT NONE
@@ -917,9 +939,9 @@
                  IF ( mrt_include_sw )  THEN
                     local_pf(i,j,k) = ( ( human_absorb * mrtinsw(l) +          &
-                                    mrtinlw(l) ) /               &  ! human_emiss * mrtinlw(l) ) /
+                                    mrtinlw(l) ) /                             &
                                     ( human_emiss * sigma_sb ) ) ** .25_wp -   &
                                     degc_to_k
                  ELSE
-                    local_pf(i,j,k) = ( mrtinlw(l)  /          &  ! ( (human_emiss * mrtinlw(l) ) /
+                    local_pf(i,j,k) = ( mrtinlw(l)  /                          &
                                     ( human_emiss * sigma_sb ) ) ** .25_wp -   &
                                     degc_to_k
@@ -1070,9 +1092,9 @@
                   IF ( mrt_include_sw )  THEN
                      local_pf(i,j,k) = REAL( ( ( human_absorb * mrtinsw(l) +   &
-                                    mrtinlw(l) ) /               &  ! human_emiss * mrtinlw(l) ) /
+                                    mrtinlw(l) ) /                             &
                                     ( human_emiss * sigma_sb ) ) ** .25_wp -   &
                                     degc_to_k, KIND = sp )
                   ELSE
-                     local_pf(i,j,k) = REAL( ( mrtinlw(l) /  &  ! REAL( ( ( human_emiss * mrtinlw(l) ) /
+                     local_pf(i,j,k) = REAL( ( mrtinlw(l) /                    &
                                     ( human_emiss * sigma_sb ) ) ** .25_wp -   &
                                     degc_to_k, KIND = sp )
@@ -1461,5 +1483,9 @@
 
 
+!
+!-- We need to differentiate if averaged input is desired (av == .TRUE.) or not.
     IF ( av ) THEN
+!
+!-- Make sure tmrt_av_grid is present and initialize with the fill value
        IF ( .NOT. ALLOCATED( tmrt_av_grid ) )  THEN
           ALLOCATE( tmrt_av_grid (nys:nyn,nxl:nxr) )
@@ -1467,16 +1493,27 @@
        tmrt_av_grid = REAL( bio_fill_value, KIND = wp )
 
-       DO  l = 1, nmrtbl
-          i = mrtbl(ix,l)
-          j = mrtbl(iy,l)
-          k = mrtbl(iz,l)
-          IF ( k - get_topography_top_index_ji( j, i, 's' ) ==                 &
-                bio_cell_level + 1_iwp) THEN
-
-             tmrt_av_grid(j,i) = mrt_av_grid(l)
-
-          ENDIF
-       ENDDO
-
+!
+!-- mrt_av_grid should always be allcoated here, but better make sure ist actually is.
+       IF ( ALLOCATED( mrt_av_grid ) )  THEN
+!
+!-- Iterate over the radiation grid (radiation coordinates) and fill the 
+!-- tmrt_av_grid (x, y coordinates) where appropriate:
+!-- tmrt_av_grid is written for all i / j if level (k) matches output height.
+          DO  l = 1, nmrtbl
+             i = mrtbl(ix,l)
+             j = mrtbl(iy,l)
+             k = mrtbl(iz,l)
+             IF ( k - get_topography_top_index_ji( j, i, 's' ) ==              &
+                   bio_cell_level + 1_iwp) THEN
+!
+!-- Averaging was done before, so we can just copy the result here
+                tmrt_av_grid(j,i) = mrt_av_grid(l)
+
+             ENDIF
+          ENDDO
+       ENDIF
+
+!
+!-- In case instantaneous input is desired, mrt values will be re-calculated.
     ELSE
 !
@@ -1497,9 +1534,9 @@
              IF ( mrt_include_sw )  THEN
                  tmrt_grid(j,i) = ( ( human_absorb * mrtinsw(l) +              &
-                                  mrtinlw(l) )  /                &  ! human_emiss * mrtinlw(l) )  /
+                                  mrtinlw(l) )  /                              &
                                   ( human_emiss * sigma_sb ) ) ** .25_wp -     &
                                   degc_to_k
              ELSE
-                 tmrt_grid(j,i) = ( mrtinlw(l)  /            &  ! ( ( human_emiss * mrtinlw(l) )  /
+                 tmrt_grid(j,i) = ( mrtinlw(l)  /                              &
                                   ( human_emiss * sigma_sb ) ) ** .25_wp -     &
                                   degc_to_k
@@ -1536,9 +1573,5 @@
 !-- Internal variables
     INTEGER(iwp)                ::  k     !< Running index, z-dir
-!     INTEGER(iwp)                ::  ir    !< Running index, x-dir, radiation coordinates
-!     INTEGER(iwp)                ::  jr    !< Running index, y-dir, radiation coordinates
-!     INTEGER(iwp)                ::  kr    !< Running index, y-dir, radiation coordinates
     INTEGER(iwp)                ::  k_wind  !< Running index, z-dir, wind speed only
-!     INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
 
     REAL(wp)                    ::  vp_sat  !< Saturation vapor pressure     (hPa)
