Index: /palm/trunk/SOURCE/biometeorology_mod.f90
===================================================================
--- /palm/trunk/SOURCE/biometeorology_mod.f90	(revision 3692)
+++ /palm/trunk/SOURCE/biometeorology_mod.f90	(revision 3693)
@@ -27,4 +27,8 @@
 ! -----------------
 ! $Id$
+! Added usage of time_averaged mean radiant temperature, together with calculation, 
+! grid and restart routines. General cleanup and commenting.
+! 
+! 3685 2019-01-21 01:02:11Z knoop
 ! Some interface calls moved to module_interface + cleanup
 ! 
@@ -113,5 +117,5 @@
               surface_pressure
 
-    USE date_and_time_mod,                                                                                                        &
+    USE date_and_time_mod,                                                     &
         ONLY:  calc_date_and_time, day_of_year, time_utc
 
@@ -125,6 +129,6 @@
     USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM
 
-    USE netcdf_data_input_mod,                                                                                                    &
-        ONLY:  netcdf_data_input_uvem, uvem_projarea_f, uvem_radiance_f,                                                          &
+    USE netcdf_data_input_mod,                                                 &
+        ONLY:  netcdf_data_input_uvem, uvem_projarea_f, uvem_radiance_f,       &
                uvem_irradiance_f, uvem_integration_f, building_obstruction_f
 !
@@ -152,4 +156,5 @@
 !-- Grids for averaged thermal indices
     REAL(wp), DIMENSION(:), ALLOCATABLE   ::  mrt_av_grid   !< time average mean
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_av_grid  !< tmrt results (degree_C)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct_av      !< PT results (aver. input)   (degree_C)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_av       !< UTCI results (aver. input) (degree_C)
@@ -165,14 +170,15 @@
 !
 !--
-    LOGICAL ::  aver_perct = .FALSE.  !< switch: do perct averaging in this module? (if .FALSE. this is done globally)
-    LOGICAL ::  aver_q     = .FALSE.  !< switch: do e  averaging in this module?
-    LOGICAL ::  aver_u     = .FALSE.  !< switch: do u  averaging in this module?
-    LOGICAL ::  aver_v     = .FALSE.  !< switch: do v  averaging in this module?
-    LOGICAL ::  aver_w     = .FALSE.  !< switch: do w  averaging in this module?
+    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?
+    LOGICAL ::  do_average_u     = .FALSE.  !< switch: do u averaging in this module?
+    LOGICAL ::  do_average_v     = .FALSE.  !< switch: do v averaging in this module?
+    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 = .TRUE.  !< Turn all thermal indices on or off
+    LOGICAL ::  thermal_comfort = .FALSE.  !< Turn all thermal indices on or off
     LOGICAL ::  bio_perct     = .TRUE.   !< Turn index PT (instant. input) on or off
     LOGICAL ::  bio_perct_av  = .TRUE.   !< Turn index PT (averaged input) on or off
@@ -378,13 +384,13 @@
 !
 !--          Averaging, as well as the allocation of the required grids must be
-!            done only once, independent from for how many thermal indices 
-!            averaged output is desired.
-!            Therefore wee need to memorize which index is the one that controls
-!            the averaging (what must be the first thermal index called).
-!            Indices are in unknown order as depending on the input file, 
-!            determine first index to average und update only once
+!--          done only once, independent from for how many thermal indices 
+!--          averaged output is desired.
+!--          Therefore wee need to memorize which index is the one that controls
+!--          the averaging (what must be the first thermal index called).
+!--          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.
+!--          is done only once during the whole model run.
              IF ( .NOT. average_trigger_perct .AND.                            &
                   .NOT. average_trigger_utci  .AND.                            &
@@ -413,5 +419,5 @@
              IF ( .NOT. ALLOCATED( pt_av ) )  THEN
                 ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                aver_perct = .TRUE.
+                do_average_theta = .TRUE.
                 pt_av = 0.0_wp
              ENDIF
@@ -419,5 +425,5 @@
              IF ( .NOT. ALLOCATED( q_av ) )  THEN
                 ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                aver_q = .TRUE.
+                do_average_q = .TRUE.
                 q_av = 0.0_wp
              ENDIF
@@ -425,5 +431,5 @@
              IF ( .NOT. ALLOCATED( u_av ) )  THEN
                 ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                aver_u = .TRUE.
+                do_average_u = .TRUE.
                 u_av = 0.0_wp
              ENDIF
@@ -431,5 +437,5 @@
              IF ( .NOT. ALLOCATED( v_av ) )  THEN
                 ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                aver_v = .TRUE.
+                do_average_v = .TRUE.
                 v_av = 0.0_wp
              ENDIF
@@ -437,6 +443,12 @@
              IF ( .NOT. ALLOCATED( w_av ) )  THEN
                 ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-                aver_w = .TRUE.
+                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
 
@@ -480,5 +492,5 @@
                 RETURN
 
-             IF ( ALLOCATED( pt_av ) .AND. aver_perct ) THEN
+             IF ( ALLOCATED( pt_av ) .AND. do_average_theta ) THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -490,5 +502,5 @@
              ENDIF
 
-             IF ( ALLOCATED( q_av )  .AND. aver_q ) THEN
+             IF ( ALLOCATED( q_av )  .AND. do_average_q ) THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -500,5 +512,5 @@
              ENDIF
 
-             IF ( ALLOCATED( u_av )  .AND. aver_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
@@ -510,5 +522,5 @@
              ENDIF
 
-             IF ( ALLOCATED( v_av )  .AND. aver_v ) THEN
+             IF ( ALLOCATED( v_av )  .AND. do_average_v ) THEN
                 DO  i = nxlg, nxrg       !< yes, ghost points are required here!
                    DO  j = nysg, nyng
@@ -520,5 +532,5 @@
              ENDIF
 
-             IF ( ALLOCATED( w_av )  .AND. aver_w ) THEN
+             IF ( ALLOCATED( w_av )  .AND. do_average_w ) THEN
                 DO  i = nxlg, nxrg       !< yes, ghost points are required here!
                    DO  j = nysg, nyng
@@ -528,4 +540,17 @@
                    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(:) + human_emiss * mrtinlw(:)) &
+                      / (human_emiss * sigma_sb)) ** .25_wp - degc_to_k
+                ELSE
+                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
+                      (human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp          &
+                      - degc_to_k
+                ENDIF
              ENDIF
 !
@@ -564,5 +589,5 @@
                 TRIM( variable ) /= 'bio_pet*' ) RETURN
 
-             IF ( ALLOCATED( pt_av ) .AND. aver_perct ) THEN
+             IF ( ALLOCATED( pt_av ) .AND. do_average_theta ) THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -575,5 +600,5 @@
              ENDIF
 
-             IF ( ALLOCATED( q_av ) .AND. aver_q ) THEN
+             IF ( ALLOCATED( q_av ) .AND. do_average_q ) THEN
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
@@ -586,5 +611,5 @@
              ENDIF
 
-             IF ( ALLOCATED( u_av ) .AND. aver_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
@@ -597,5 +622,5 @@
              ENDIF
 
-             IF ( ALLOCATED( v_av ) .AND. aver_v ) THEN
+             IF ( ALLOCATED( v_av ) .AND. do_average_v ) THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -608,5 +633,5 @@
              ENDIF
 
-             IF ( ALLOCATED( w_av ) .AND. aver_w ) THEN
+             IF ( ALLOCATED( w_av ) .AND. do_average_w ) THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -618,4 +643,9 @@
                 ENDDO
              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
@@ -793,5 +823,4 @@
 
 !------------------------------------------------------------------------------!
-!
 ! Description:
 ! ------------
@@ -817,5 +846,6 @@
     CHARACTER (LEN=*), INTENT(OUT) ::  grid   !< Grid type (always "zu1" for biom)
     LOGICAL, INTENT(OUT)           ::  found  !< Output found?
-    LOGICAL, INTENT(OUT)           ::  two_d  !< Flag parameter that indicates 2D variables, horizontal cross sections, must be .TRUE.
+    LOGICAL, INTENT(OUT)           ::  two_d  !< Flag parameter that indicates 2D variables,
+                                              !< horizontal cross sections, must be .TRUE. for thermal indices and uv
     REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< Temp. result grid to return
 !
@@ -911,5 +941,5 @@
            END IF
 
-           !
+!
 !--    Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein. 
 !--    However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged. 
@@ -950,5 +980,4 @@
 
 !------------------------------------------------------------------------------!
-!
 ! Description:
 ! ------------
@@ -1211,4 +1240,161 @@
 ! Description:
 ! ------------
+!> Soubroutine reads global biometeorology configuration from restart file(s)
+!------------------------------------------------------------------------------!
+ SUBROUTINE bio_rrd_global( found )
+
+    USE control_parameters,                                                    &
+        ONLY:  length, restart_string
+
+
+    IMPLICIT NONE
+
+    LOGICAL, INTENT(OUT) ::  found      !< variable found? yes = .T., no = .F.
+
+    found = .TRUE.
+
+
+    SELECT CASE ( restart_string(1:length) )
+
+!
+!--    read control flags to determine if input grids need to be averaged
+       CASE ( 'do_average_theta' )
+          READ ( 13 )  do_average_theta
+
+       CASE ( 'do_average_q' )
+          READ ( 13 )  do_average_q
+
+       CASE ( 'do_average_u' )
+          READ ( 13 )  do_average_u
+
+       CASE ( 'do_average_v' )
+          READ ( 13 )  do_average_v
+
+       CASE ( 'do_average_w' )
+          READ ( 13 )  do_average_w
+
+       CASE ( 'do_average_mrt' )
+          READ ( 13 )  do_average_mrt
+
+!
+!--    read control flags to determine which thermal index needs to trigger averaging
+       CASE ( 'average_trigger_perct' )
+          READ ( 13 )  average_trigger_perct
+
+       CASE ( 'average_trigger_utci' )
+          READ ( 13 )  average_trigger_utci
+
+       CASE ( 'average_trigger_pet' )
+          READ ( 13 )  average_trigger_pet
+
+
+       CASE DEFAULT
+
+          found = .FALSE.
+
+    END SELECT
+
+
+ END SUBROUTINE bio_rrd_global
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Soubroutine reads local biometeorology configuration from restart file(s)
+!------------------------------------------------------------------------------!
+ SUBROUTINE bio_rrd_local( found )
+
+
+    USE control_parameters,                                                    &
+        ONLY:  length, restart_string
+
+
+    IMPLICIT NONE
+
+
+    LOGICAL, INTENT(OUT) ::  found      !< variable found? yes = .T., no = .F.
+
+    found = .TRUE.
+
+
+    SELECT CASE ( restart_string(1:length) )
+
+       CASE ( 'nmrtbl' )
+          READ ( 13 )  bio_nmrtbl
+
+       CASE ( 'mrt_av_grid' )
+          IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
+             ALLOCATE( mrt_av_grid(bio_nmrtbl) )
+          ENDIF
+          READ ( 13 )  mrt_av_grid
+
+
+       CASE DEFAULT
+
+          found = .FALSE.
+
+    END SELECT
+
+
+ END SUBROUTINE bio_rrd_local
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Write global restart data for the biometeorology module.
+!------------------------------------------------------------------------------!
+ SUBROUTINE bio_wrd_global
+
+    IMPLICIT NONE
+
+    CALL wrd_write_string( 'do_average_theta' )
+    WRITE ( 14 )  do_average_theta
+    CALL wrd_write_string( 'do_average_q' )
+    WRITE ( 14 )  do_average_q
+    CALL wrd_write_string( 'do_average_u' )
+    WRITE ( 14 )  do_average_u
+    CALL wrd_write_string( 'do_average_v' )
+    WRITE ( 14 )  do_average_v
+    CALL wrd_write_string( 'do_average_w' )
+    WRITE ( 14 )  do_average_w
+    CALL wrd_write_string( 'do_average_mrt' )
+    WRITE ( 14 )  do_average_mrt
+    CALL wrd_write_string( 'average_trigger_perct' )
+    WRITE ( 14 )  average_trigger_perct
+    CALL wrd_write_string( 'average_trigger_utci' )
+    WRITE ( 14 )  average_trigger_utci
+    CALL wrd_write_string( 'average_trigger_pet' )
+    WRITE ( 14 )  average_trigger_pet
+
+ END SUBROUTINE bio_wrd_global
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Write local restart data for the biometeorology module.
+!------------------------------------------------------------------------------!
+ SUBROUTINE bio_wrd_local
+
+    IMPLICIT NONE
+
+!
+!-- First nmrtbl has to be written/read, because it is the dimension of mrt_av_grid
+    CALL wrd_write_string( 'nmrtbl' )
+    WRITE ( 14 )  nmrtbl
+
+    IF ( ALLOCATED( mrt_av_grid ) )  THEN
+       CALL wrd_write_string( 'mrt_av_grid' )
+       WRITE ( 14 )  mrt_av_grid
+    ENDIF
+
+
+ END SUBROUTINE bio_wrd_local
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
 !> Calculate biometeorology MRT for all 2D grid
 !------------------------------------------------------------------------------!
@@ -1225,29 +1411,50 @@
     INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
 
+
+    IF ( av ) THEN
+       IF ( .NOT. ALLOCATED( tmrt_av_grid ) )  THEN
+          ALLOCATE( tmrt_av_grid (nys:nyn,nxl:nxr) )
+       ENDIF
+
+       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
+
+    ELSE
+
 !
 !-- Calculate biometeorology MRT from local radiation fluxes calculated by RTM and assign
 !-- into 2D grid. Depending on selected output quantities, tmrt_grid might not have been
 !-- allocated in bio_check_data_output yet.
-    IF ( .NOT. ALLOCATED( tmrt_grid ) )  THEN
-       ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
-    ENDIF
-    tmrt_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
-          IF ( mrt_include_sw )  THEN
-              tmrt_grid(j,i) = ((human_absorb*mrtinsw(l) +                     &
-                                human_emiss*mrtinlw(l))  /                     &
-                               (human_emiss*sigma_sb)) ** .25_wp - degc_to_k
-          ELSE
-              tmrt_grid(j,i) = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp   &
-                                 - degc_to_k
+       IF ( .NOT. ALLOCATED( tmrt_grid ) )  THEN
+          ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
+       ENDIF
+       tmrt_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
+             IF ( mrt_include_sw )  THEN
+                 tmrt_grid(j,i) = ((human_absorb*mrtinsw(l) +                     &
+                                   human_emiss*mrtinlw(l))  /                     &
+                                  (human_emiss*sigma_sb)) ** .25_wp - degc_to_k
+             ELSE
+                 tmrt_grid(j,i) = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp   &
+                                    - degc_to_k
+             ENDIF
           ENDIF
-       ENDIF
-    ENDDO
+       ENDDO
+    ENDIF
 
 END SUBROUTINE bio_calculate_mrt_grid
@@ -1342,8 +1549,10 @@
 !-- local mtr value at [i,j]
     tmrt = bio_fill_value  !< this can be a valid result (e.g. for inside some ostacle)
-    IF ( radiation ) THEN
+    IF ( .NOT. average_input ) THEN
 !
 !--    Use MRT from RTM precalculated in tmrt_grid
        tmrt = tmrt_grid(j,i)
+    ELSE
+       tmrt = tmrt_av_grid(j,i)
     ENDIF
 
@@ -1488,4 +1697,9 @@
     REAL(wp), INTENT ( OUT ) ::  ipt    !< Instationary perceived temp.   (degree_C)
 !
+!-- return immediatelly if nothing to do!
+    IF ( .NOT. thermal_comfort ) THEN
+        RETURN
+    ENDIF
+!
 !-- If clo equals the initial value, this is the initial call
     IF ( clo <= -998._wp ) THEN
@@ -1505,153 +1719,4 @@
  END SUBROUTINE bio_calc_ipt
 
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Soubroutine reads global biometeorology configuration from restart file(s)
-!------------------------------------------------------------------------------!
- SUBROUTINE bio_rrd_global( found )
-
-    USE control_parameters,                                                    &
-        ONLY:  length, restart_string
-
-
-    IMPLICIT NONE
-
-    LOGICAL, INTENT(OUT) ::  found      !< variable found? yes = .T., no = .F.
-
-    found = .TRUE.
-
-
-    SELECT CASE ( restart_string(1:length) )
-
-!
-!--    read control flags to determine if input grids need to be averaged
-       CASE ( 'aver_perct' )
-          READ ( 13 )  aver_perct
-
-       CASE ( 'aver_q' )
-          READ ( 13 )  aver_q
-
-       CASE ( 'aver_u' )
-          READ ( 13 )  aver_u
-
-       CASE ( 'aver_v' )
-          READ ( 13 )  aver_v
-
-       CASE ( 'aver_w' )
-          READ ( 13 )  aver_w
-!
-!--    read control flags to determine which thermal index needs to trigger averaging
-       CASE ( 'average_trigger_perct' )
-          READ ( 13 )  average_trigger_perct
-
-       CASE ( 'average_trigger_utci' )
-          READ ( 13 )  average_trigger_utci
-
-       CASE ( 'average_trigger_pet' )
-          READ ( 13 )  average_trigger_pet
-
-
-       CASE DEFAULT
-
-          found = .FALSE.
-
-    END SELECT
-
-
- END SUBROUTINE bio_rrd_global
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Soubroutine reads local biometeorology configuration from restart file(s)
-!------------------------------------------------------------------------------!
- SUBROUTINE bio_rrd_local( found )
-
-
-    USE control_parameters,                                                    &
-        ONLY:  length, restart_string
-
-
-    IMPLICIT NONE
-
-
-    LOGICAL, INTENT(OUT) ::  found      !< variable found? yes = .T., no = .F.
-
-    found = .TRUE.
-
-
-    SELECT CASE ( restart_string(1:length) )
-
-       CASE ( 'nmrtbl' )
-          READ ( 13 )  bio_nmrtbl
-
-       CASE ( 'mrt_av_grid' )
-          IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
-             ALLOCATE( mrt_av_grid(bio_nmrtbl) )
-          ENDIF
-          READ ( 13 )  mrt_av_grid
-
-
-       CASE DEFAULT
-
-          found = .FALSE.
-
-    END SELECT
-
-
- END SUBROUTINE bio_rrd_local
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Write global restart data for the biometeorology module.
-!------------------------------------------------------------------------------!
- SUBROUTINE bio_wrd_global
-
-    IMPLICIT NONE
-
-    CALL wrd_write_string( 'aver_perct' )
-    WRITE ( 14 )  aver_perct
-    CALL wrd_write_string( 'aver_q' )
-    WRITE ( 14 )  aver_q
-    CALL wrd_write_string( 'aver_u' )
-    WRITE ( 14 )  aver_u
-    CALL wrd_write_string( 'aver_v' )
-    WRITE ( 14 )  aver_v
-    CALL wrd_write_string( 'aver_w' )
-    WRITE ( 14 )  aver_w
-    CALL wrd_write_string( 'average_trigger_perct' )
-    WRITE ( 14 )  average_trigger_perct
-    CALL wrd_write_string( 'average_trigger_utci' )
-    WRITE ( 14 )  average_trigger_utci
-    CALL wrd_write_string( 'average_trigger_pet' )
-    WRITE ( 14 )  average_trigger_pet
-
- END SUBROUTINE bio_wrd_global
-
-
-!------------------------------------------------------------------------------!
-! Description:
-! ------------
-!> Write local restart data for the biometeorology module.
-!------------------------------------------------------------------------------!
- SUBROUTINE bio_wrd_local
-
-    IMPLICIT NONE
-
-!
-!-- First nmrtbl has to be written/read, because it is the dimension of mrt_av_grid
-    CALL wrd_write_string( 'nmrtbl' )
-    WRITE ( 14 )  nmrtbl
-
-    IF ( ALLOCATED( mrt_av_grid ) )  THEN
-       CALL wrd_write_string( 'mrt_av_grid' )
-       WRITE ( 14 )  mrt_av_grid
-    ENDIF
-
-
- END SUBROUTINE bio_wrd_local
 
 
@@ -1661,5 +1726,5 @@
 !> SUBROUTINE for calculating UTCI Temperature (UTCI)
 !> computed by a 6th order approximation
-!
+!>
 !> UTCI regression equation after
 !> BrÃ¶de P, Fiala D, Blazejczyk K, HolmÃ©r I, Jendritzky G, Kampmann B, Tinz B, 
@@ -1667,5 +1732,5 @@
 !> Climate Index (UTCI). International Journal of Biometeorology 56 (3):481-494. 
 !> doi:10.1007/s00484-011-0454-1
-!
+!>
 !> original source available at:
 !> www.utci.org
@@ -1728,5 +1793,5 @@
 !
 !-- Wind altitude correction from hag to 10m after Broede et al. (2012), eq.3
-!   z(0) is set to 0.01 according to UTCI profile definition
+!-- z(0) is set to 0.01 according to UTCI profile definition
     va = ws_hag *  log ( 10.0_wp / 0.01_wp ) / log ( hag / 0.01_wp )
 !
@@ -1749,6 +1814,6 @@
 !
 !-- For routine application. For wind speeds and relative
-!   humidity values below 0.5 m/s or 5%, respectively, the
-!   user is advised to use the lower bounds for the calculations.
+!-- 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
@@ -1988,5 +2053,5 @@
        (  1.33374846e-03_wp ) * ta2 *                 pa4 +                    &
        (  3.55375387e-03_wp ) *       va  *           pa4 +                    &
-       ( -5.13027851e-04_wp ) * ta  * va *            pa4 +                    &
+       ( -5.13027851e-04_wp ) * ta  * va  *           pa4 +                    &
        (  1.02449757e-04_wp ) *       va2 *           pa4 +                    &
        ( -1.48526421e-03_wp ) *             d_tmrt  * pa4 +                    &
@@ -2034,6 +2099,7 @@
 !
 !-- Parameters for standard "Klima-Michel"
-    REAL(wp), PARAMETER :: eta = 0._wp  !< Mechanical work efficiency for walking on flat ground (compare to Fanger (1972) pp 24f)
-    REAL(wp), PARAMETER :: actlev = 134.6862_wp !< Workload by activity per standardized surface (A_Du)
+    REAL(wp), PARAMETER :: eta = 0._wp  !< Mechanical work efficiency for walking on flat ground
+                                        !< (compare to Fanger (1972) pp 24f)
+    REAL(wp), PARAMETER :: actlev = 134.6862_wp  !< Workload by activity per standardized surface (A_Du)
 !
 !-- Type of program variables
@@ -2072,8 +2138,8 @@
 !-- Tresholds: clothing insulation (account for model inaccuracies)
 !
-!   summer clothing
+!-- summer clothing
     sclo     = 0.44453_wp
 !
-!   winter clothing
+!-- winter clothing
     wclo     = 1.76267_wp
 !
@@ -2096,5 +2162,5 @@
 !
 !--          Case: comfort achievable by varying clothing insulation 
-!            Between winter and summer set values
+!--          Between winter and summer set values
              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
                 pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
@@ -2131,5 +2197,5 @@
 !
 !--          Case: comfort achievable by varying clothing insulation 
-!            between winter and summer set values
+!--          between winter and summer set values
              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
                                pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
@@ -2172,6 +2238,6 @@
 !
 !--    Required clothing insulation (ireq) is exclusively defined for 
-!      operative temperatures (top) less 10 (C) for a 
-!      reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
+!--    operative temperatures (top) less 10 (C) for a 
+!--    reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
        clon = ireq_neutral ( perct_ij, ireq_minimal, nerr )
        clo = clon
@@ -2202,7 +2268,9 @@
 ! Description:
 ! ------------
-!> The SUBROUTINE calculates the saturation water vapour pressure 
+!> The SUBROUTINE calculates the (saturation) water vapour pressure 
 !> (hPa = hecto Pascal) for a given temperature ta (degC).
-!> For example, ta can be the air temperature or the dew point temperature.
+!> 'ta' can be the air temperature or the dew point temperature. The first will
+!> result in the current vapor pressure (hPa), the latter will calulate the
+!> saturation vapor pressure (hPa).
 !------------------------------------------------------------------------------!
  SUBROUTINE saturation_vapor_pressure( ta, svp_ta )
@@ -2211,5 +2279,5 @@
 
     REAL(wp), INTENT ( IN )  ::  ta     !< ambient air temperature (degC)
-    REAL(wp), INTENT ( OUT ) ::  svp_ta !< saturation water vapour pressure (hPa)
+    REAL(wp), INTENT ( OUT ) ::  svp_ta !< water vapour pressure (hPa)
 
     REAL(wp)      ::  b
@@ -2219,10 +2287,10 @@
     IF ( ta < 0._wp ) THEN
 !
-!--    ta  < 0 (degC): saturation water vapour pressure over ice
+!--    ta  < 0 (degC): water vapour pressure over ice
        b = 17.84362_wp
        c = 245.425_wp
     ELSE
 !
-!--    ta >= 0 (degC): saturation water vapour pressure over water
+!--    ta >= 0 (degC): water vapour pressure over water
        b = 17.08085_wp
        c = 234.175_wp
@@ -2248,9 +2316,9 @@
 !
 !-- Input variables of argument list:
-    REAL(wp), INTENT ( IN )  :: ta     !< Ambient temperature (degC)
+    REAL(wp), INTENT ( IN )  :: ta       !< Ambient temperature (degC)
     REAL(wp), INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC)
-    REAL(wp), INTENT ( IN )  :: vp     !< Water vapour pressure (hPa)
-    REAL(wp), INTENT ( IN )  :: ws    !< Wind speed (m/s) 1 m above ground
-    REAL(wp), INTENT ( IN )  :: pair       !< Barometric pressure (hPa)
+    REAL(wp), INTENT ( IN )  :: vp       !< Water vapour pressure (hPa)
+    REAL(wp), INTENT ( IN )  :: ws       !< Wind speed (m/s) 1 m above ground
+    REAL(wp), INTENT ( IN )  :: pair     !< Barometric air pressure (hPa)
     REAL(wp), INTENT ( IN )  :: actlev   !< Individuals activity level per unit surface area (W/m2)
     REAL(wp), INTENT ( IN )  :: eta      !< Individuals work efficiency (dimensionless)
@@ -2262,14 +2330,13 @@
     REAL(wp), INTENT ( IN )  :: pmv_s    !< Fanger's PMV corresponding to sclo
 !
-! Output variables of argument list:
+!-- Output variables of argument list:
     REAL(wp), INTENT ( OUT ) :: pmva     !< 0 (set to zero, because clo is evaluated for comfort)
     REAL(wp), INTENT ( OUT ) :: top      !< Operative temperature (degC) at found root of Fanger's PMV
     REAL(wp), INTENT ( OUT ) :: clo_res  !< Resulting clothing insulation value (clo)
     INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error status / quality flag
-!           nerr >= 0, o.k., and nerr is the number of iterations for
-!                              convergence
-!           nerr = -1: error = malfunction of Ridder's convergence method
-!           nerr = -2: error = maximum iterations (max_iteration) exceeded 
-!           nerr = -3: error = root not bracketed between sclo and wclo
+                                         !< nerr >= 0, o.k., and nerr is the number of iterations for convergence
+                                         !< nerr = -1: error = malfunction of Ridder's convergence method
+                                         !< nerr = -2: error = maximum iterations (max_iteration) exceeded 
+                                         !< nerr = -3: error = root not bracketed between sclo and wclo
 !
 !-- Type of program variables
@@ -2428,6 +2495,6 @@
 !-- Output variables of argument list:
     REAL(wp), INTENT ( OUT ) ::  pmva    !< Actual Predicted Mean Vote (PMV, 
-!            dimensionless) according to Fanger corresponding to meteorological 
-!            (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta)
+                                         !< dimensionless) according to Fanger corresponding to meteorological 
+                                         !< (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta)
     REAL(wp), INTENT ( OUT ) ::  top     !< operative temperature (degC)
 !
@@ -2447,5 +2514,5 @@
     REAL(wp) ::  ws           !< wind speed                             (m/s)
     REAL(wp) ::  z1           !< Empiric factor for the adaption of the heat 
-!             ballance equation to the psycho-physical scale (Equ. 40 in FANGER)
+                              !< ballance equation to the psycho-physical scale (Equ. 40 in FANGER)
     REAL(wp) ::  z2           !< Water vapour diffution through the skin
     REAL(wp) ::  z3           !< Sweat evaporation from the skin surface
@@ -2485,5 +2552,5 @@
 !
 !-- Calculation of clothing surface temperature (t_clothing) based on
-!   Newton-approximation with air temperature as initial guess
+!-- Newton-approximation with air temperature as initial guess
     t_clothing = ta
     DO i = 1, 3
@@ -2493,5 +2560,5 @@
 !
 !-- Empiric factor for the adaption of the heat ballance equation
-!   to the psycho-physical scale (Equ. 40 in FANGER)
+!-- to the psycho-physical scale (Equ. 40 in FANGER)
     z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp )
 !
@@ -2535,4 +2602,5 @@
 
     IMPLICIT NONE
+
 !
 !-- Input variables of argument list:
@@ -2543,24 +2611,26 @@
     REAL(wp),     INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC) at screen level
     REAL(wp),     INTENT ( IN )  :: ws       !< Wind speed (m/s) 1 m above ground
+
 !
 !-- Output variables of argument list:
     INTEGER(iwp), INTENT ( OUT ) :: nerr     !< Error status / quality flag
-!             0 = o.k.
-!            -2 = pmva outside valid regression range
-!            -3 = rel. humidity set to 5 % or 95 %, respectively
-!            -4 = deltapmv set to avoid pmvs < 0
-!
-!-- Internal variable types:
-    REAL(wp) ::  pmv          !< 
-    REAL(wp) ::  pa_p50       !< 
-    REAL(wp) ::  pa           !< 
-    REAL(wp) ::  apa          !< 
-    REAL(wp) ::  dapa         !< 
-    REAL(wp) ::  sqvel        !< 
-    REAL(wp) ::  dtmrt        !< 
-    REAL(wp) ::  p10          !< 
-    REAL(wp) ::  p95          !< 
-    REAL(wp) ::  gew          !< 
-    REAL(wp) ::  gew2         !< 
+                                             !<  0 = o.k.
+                                             !< -2 = pmva outside valid regression range
+                                             !< -3 = rel. humidity set to 5 % or 95 %, respectively
+                                             !< -4 = deltapmv set to avoid pmvs < 0
+
+!
+!-- Internal variables:
+    REAL(wp) ::  pmv          !< temp storage og predicted mean vote
+    REAL(wp) ::  pa_p50       !< ratio actual water vapour pressure to that of relative humidity of 50 %
+    REAL(wp) ::  pa           !< vapor pressure (hPa) with hard bounds
+    REAL(wp) ::  apa          !< natural logarithm of pa (with hard lower border)
+    REAL(wp) ::  dapa         !< difference of apa and pa_p50
+    REAL(wp) ::  sqvel        !< square root of local wind velocity
+    REAL(wp) ::  dtmrt        !< difference mean radiation to air temperature
+    REAL(wp) ::  p10          !< lower bound for pa
+    REAL(wp) ::  p95          !< upper bound for pa
+    REAL(wp) ::  weight       !< 
+    REAL(wp) ::  weight2      !< 
     REAL(wp) ::  dpmv_1       !< 
     REAL(wp) ::  dpmv_2       !< 
@@ -2657,6 +2727,4 @@
     ENDIF
 !
-!-- Air temperature
-!    ta = ta
 !-- Difference mean radiation to air temperature
     dtmrt = tmrt - ta
@@ -2670,13 +2738,13 @@
        RETURN
     ENDIF
-    gew = MOD ( pmv, 1._wp )
-    IF ( gew < 0._wp ) gew = 0._wp
+    weight = MOD ( pmv, 1._wp )
+    IF ( weight < 0._wp ) weight = 0._wp
     IF ( nreg > 5_iwp ) THEN
        ! nreg=6
        nreg  = 5_iwp
-       gew   = pmv - 5._wp
-       gew2  = pmv - 6._wp
-       IF ( gew2 > 0_iwp ) THEN
-          gew = ( gew - gew2 ) / gew
+       weight   = pmv - 5._wp
+       weight2  = pmv - 6._wp
+       IF ( weight2 > 0_iwp ) THEN
+          weight = ( weight - weight2 ) / weight
        ENDIF
     ENDIF
@@ -2709,5 +2777,5 @@
 !
 !-- Calculate pmv modification
-    deltapmv = ( 1._wp - gew ) * dpmv_1 + gew * dpmv_2
+    deltapmv = ( 1._wp - weight ) * dpmv_1 + weight * dpmv_2
     pmvs = pmva + deltapmv
     IF ( ( pmvs ) < 0._wp ) THEN
@@ -2746,7 +2814,7 @@
 !-- Additional output variables of argument list:
     REAL(wp), INTENT ( OUT ) ::  dperctm    !< Mean deviation perct (classical gt) to gt* (rational gt 
-!             calculated based on Gagge's rational PMV*)
+                                            !< calculated based on Gagge's rational PMV*)
     REAL(wp), INTENT ( OUT ) ::  dperctstd  !< dperctm plus its standard deviation times a factor 
-!             determining the significance to perceive sultriness
+                                            !< determining the significance to perceive sultriness
     REAL(wp), INTENT ( OUT ) ::  sultr_res
 !
@@ -2758,5 +2826,5 @@
 !
 !-- Types of coefficients mean deviation plus standard deviation
-!   regression coefficients: third order polynomial
+!-- regression coefficients: third order polynomial
     REAL(wp), PARAMETER ::  dperctsa = +0.0268918_wp
     REAL(wp), PARAMETER ::  dperctsb = +0.0465957_wp
@@ -2765,5 +2833,5 @@
 !
 !-- Factor to mean standard deviation defining SIGNificance for 
-!   sultriness
+!-- sultriness
     REAL(wp), PARAMETER :: faktor = 1._wp
 !
@@ -2814,5 +2882,5 @@
     INTEGER(iwp), INTENT ( OUT ) ::  nerr !< Error indicator: 0 = o.k., +1 = denominator for intersection = 0
     REAL(wp),     INTENT ( OUT ) ::  dpmv_cold_res    !< Increment to adjust pmva according to the results of Gagge's 
-!               2 node model depending on the input
+                                                      !< 2 node model depending on the input
 !
 !-- Type of program variables
@@ -2826,5 +2894,4 @@
     REAL(wp) ::  sqrt_ws
     INTEGER(iwp) ::  i
-!     INTEGER(iwp) ::  j
     INTEGER(iwp) ::  i_bin
 !
@@ -2839,9 +2906,9 @@
 !
 !-- Initialise
-    nerr       = 0_iwp
+    nerr           = 0_iwp
     dpmv_cold_res  = 0._wp
-    pmvc       = pmva
-    dtmrt      = tmrt - ta
-    sqrt_ws   = ws
+    pmvc           = pmva
+    dtmrt          = tmrt - ta
+    sqrt_ws        = ws
     IF ( sqrt_ws < 0.10_wp ) THEN
        sqrt_ws = 0.10_wp
@@ -2885,5 +2952,5 @@
 !
 !-- Adjust to operative temperature scaled according 
-!   to classical PMV (Fanger)
+!-- to classical PMV (Fanger)
     dpmv_cold_res = delta_cold(i_bin) - dpmv_adj(pmva)
 
@@ -2909,6 +2976,6 @@
     REAL(wp)      ::  pmv
     INTEGER(iwp)  ::  i, i_bin
-
-!                                   range_1        range_2        range_3
+!
+!--                                 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/
@@ -3011,5 +3078,5 @@
 !
 !-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if 
-!   height and weight be known. In: Arch. Int. Med.. 17, 1916, S. 863?871.
+!-- height and weight be known. In: Arch. Int. Med.. 17, 1916, S. 863?871.
     surf = 0.007184_wp * height**0.725_wp * weight**0.425_wp
     RETURN
@@ -3033,5 +3100,5 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE persdat ( age, weight, height, sex, work, a_surf, actlev )
-!
+
     IMPLICIT NONE
 
@@ -3164,5 +3231,5 @@
 !
 !--          Case: comfort achievable by varying clothing insulation 
-!            between winter and summer set values
+!--          between winter and summer set values
              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta , sclo,     &
                             pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
@@ -3205,5 +3272,5 @@
 !
 !--          Case: comfort achievable by varying clothing insulation 
-!            between winter and summer set values
+!--          between winter and summer set values
              CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
                                pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
@@ -3248,6 +3315,6 @@
 !
 !--    Required clothing insulation (ireq) is exclusively defined for 
-!      operative temperatures (top) less 10 (C) for a 
-!      reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
+!--    operative temperatures (top) less 10 (C) for a 
+!--    reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
        clon = ireq_neutral ( ipt, ireq_minimal, nerr )
        clo = clon
@@ -3387,5 +3454,5 @@
     REAL(wp), INTENT ( IN )  ::  pa       !< Vapour pressure          (hPa)
     REAL(wp), INTENT ( IN )  ::  pair     !< Air pressure             (hPa)
-    REAL(wp), INTENT ( IN )  ::  in_ws   !< Wind speed               (m/s)
+    REAL(wp), INTENT ( IN )  ::  in_ws    !< Wind speed               (m/s)
     REAL(wp), INTENT ( IN )  ::  actlev   !< Metabolic + work energy  (W/mÂ²)
     REAL(wp), INTENT ( IN )  ::  dt       !< Timestep                 (s)
@@ -3415,5 +3482,5 @@
     REAL(wp) ::  ws           !< wind speed                             (m/s)
     REAL(wp) ::  z1           !< Empiric factor for the adaption of the heat 
-!             ballance equation to the psycho-physical scale (Equ. 40 in FANGER)
+                              !< ballance equation to the psycho-physical scale (Equ. 40 in FANGER)
     REAL(wp) ::  z2           !< Water vapour diffution through the skin
     REAL(wp) ::  z3           !< Sweat evaporation from the skin surface
@@ -3427,5 +3494,5 @@
 
     INTEGER(iwp) :: i         !< running index
-    INTEGER(iwp) ::  niter  !< Running index
+    INTEGER(iwp) ::  niter    !< Running index
 
 !
@@ -3457,5 +3524,5 @@
 !
 !-- Calculation of clothing surface temperature (t_clothing) based on
-!   newton-approximation with air temperature as initial guess
+!-- newton-approximation with air temperature as initial guess
     niter = INT( dt * 10._wp, KIND=iwp )
     IF ( niter < 1 ) niter = 1_iwp
@@ -3479,5 +3546,5 @@
 !
 !-- Empiric factor for the adaption of the heat ballance equation
-!   to the psycho-physical scale (Equ. 40 in FANGER)
+!-- to the psycho-physical scale (Equ. 40 in FANGER)
     z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp )
 !
@@ -3617,8 +3684,4 @@
     REAL(wp), INTENT( OUT ) ::  int_heat  !< internal heat production (W)
     REAL(wp), INTENT( OUT ) ::  rtv       !< respiratory volume
-!
-!-- Constants:
-!     REAL(wp), PARAMETER :: cair = 1010._wp        !< replaced by c_p
-!     REAL(wp), PARAMETER :: evap = 2.42_wp * 10._wp **6._wp  !< replaced by l_v
 !
 !-- Internal variables:
