Ignore:
Timestamp:
Nov 14, 2018 4:06:14 PM (3 years ago)
Author:
kanani
Message:

Changes related to clean-up of biometeorology (by dom_dwd_user)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/biometeorology_mod.f90

    r3479 r3525  
    2727! -----------------
    2828! $Id$
     29! Clean up, renaming from "biom" to "bio", summary of thermal index calculation
     30! into one module (dom_dwd_user)
     31!
     32! 3479 2018-11-01 16:00:30Z gronemeier
    2933! - reworked check for output quantities
    3034! - assign new palm-error numbers
     
    7276
    7377    USE basic_constants_and_equations_mod,                                     &
    74        ONLY:  degc_to_k, magnus, sigma_sb
    75 
    76     USE biometeorology_ipt_mod
    77 
    78     USE biometeorology_pet_mod
    79 
    80     USE biometeorology_pt_mod,                                                 &
    81        ONLY: calculate_pt_static
    82 
    83     USE biometeorology_utci_mod
     78       ONLY:  c_p, degc_to_k, l_v, magnus, sigma_sb
    8479
    8580    USE control_parameters,                                                    &
     
    9287
    9388    USE indices,                                                               &
    94        ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr
     89       ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg,    &
     90              nysg, nyng
    9591
    9692    USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM
     
    112108!-- Declare all global variables within the module (alphabetical order)
    113109    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_grid  !< tmrt results (°C)
    114     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_grid    !< PT results   (°C)
     110    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct      !< PT results   (°C)
    115111    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_grid  !< UTCI results (°C)
    116112    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_grid   !< PET results  (°C)
    117 !-- Grids for averaged thermal indices       
    118     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pt_av_grid    !< PT results (aver. input)   (°C)
     113!-- Grids for averaged thermal indices
     114    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  mrt_av_grid   !< time average mean
     115    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct_av      !< PT results (aver. input)   (°C)
    119116    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_av_grid  !< UTCI results (aver. input) (°C)
    120117    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_av_grid   !< PET results (aver. input)  (°C)
    121 !-- Grids for radiation_model
    122     REAL(wp), DIMENSION(:), ALLOCATABLE ::  biom_mrt     !< biometeorology mean radiant temperature for each MRT box
    123     REAL(wp), DIMENSION(:), ALLOCATABLE ::  biom_mrt_av  !< time average mean
    124 
    125     INTEGER( iwp ) ::  biom_cell_level     !< cell level biom calculates for
    126     REAL ( wp )    ::  biom_output_height  !< height output is calculated in m
    127     REAL ( wp )    ::  time_biom_results   !< the time, the last set of biom results have been calculated for
     118
     119
     120    INTEGER( iwp ) ::  bio_cell_level     !< cell level biom calculates for
     121    REAL ( wp )    ::  bio_output_height  !< height output is calculated in m
     122    REAL ( wp )    ::  time_bio_results   !< the time, the last set of biom results have been calculated for
    128123    REAL ( wp ), PARAMETER ::  human_absorb = 0.7_wp  !< SW absorbtivity of a human body (Fanger 1972)
    129124    REAL ( wp ), PARAMETER ::  human_emiss = 0.97_wp  !< LW emissivity of a human body after (Fanger 1972)
    130125!--
    131126
    132     LOGICAL ::  aver_pt  = .FALSE.  !< switch: do pt averaging in this module? (if .FALSE. this is done globally)
    133     LOGICAL ::  aver_q   = .FALSE.  !< switch: do e  averaging in this module?
    134     LOGICAL ::  aver_u   = .FALSE.  !< switch: do u  averaging in this module?
    135     LOGICAL ::  aver_v   = .FALSE.  !< switch: do v  averaging in this module?
    136     LOGICAL ::  aver_w   = .FALSE.  !< switch: do w  averaging in this module?
    137     LOGICAL ::  average_trigger_pt   = .FALSE.  !< update averaged input on call to biom_pt?
    138     LOGICAL ::  average_trigger_utci = .FALSE.  !< update averaged input on call to biom_utci?
    139     LOGICAL ::  average_trigger_pet  = .FALSE.  !< update averaged input on call to biom_pet?
    140 
    141     LOGICAL ::  mrt_from_rtm   = .TRUE.   !< switch: use mrt calculated by RTM for calculation of thermal indices
    142     LOGICAL ::  biom_pt        = .TRUE.   !< Turn index PT (instant. input) on or off
    143     LOGICAL ::  biom_pt_av     = .TRUE.   !< Turn index PT (averaged input) on or off
    144     LOGICAL ::  biom_pet       = .TRUE.   !< Turn index PET (instant. input) on or off
    145     LOGICAL ::  biom_pet_av    = .TRUE.   !< Turn index PET (averaged input) on or off
    146     LOGICAL ::  biom_utci      = .TRUE.   !< Turn index UTCI (instant. input) on or off
    147     LOGICAL ::  biom_utci_av   = .TRUE.   !< Turn index UTCI (averaged input) on or off
    148 
    149 !-- Add INTERFACES that must be available to other modules (alphabetical order)
    150 
    151     PUBLIC biom_3d_data_averaging, biom_check_data_output,                     &
    152     biom_calculate_static_grid, biom_calc_ipt,                                 &
    153     biom_check_parameters, biom_data_output_3d, biom_data_output_2d,           &
    154     biom_define_netcdf_grid, biom_determine_input_at, biom_header, biom_init,  &
    155     biom_init_arrays, biom_parin, biom_pt, biom_pt_av, biom_pet, biom_pet_av,  &
    156     biom_utci, biom_utci_av, time_biom_results
     127    LOGICAL ::  aver_perct = .FALSE.  !< switch: do perct averaging in this module? (if .FALSE. this is done globally)
     128    LOGICAL ::  aver_q     = .FALSE.  !< switch: do e  averaging in this module?
     129    LOGICAL ::  aver_u     = .FALSE.  !< switch: do u  averaging in this module?
     130    LOGICAL ::  aver_v     = .FALSE.  !< switch: do v  averaging in this module?
     131    LOGICAL ::  aver_w     = .FALSE.  !< switch: do w  averaging in this module?
     132    LOGICAL ::  aver_mrt   = .FALSE.  !< switch: do mrt averaging in this module?
     133    LOGICAL ::  average_trigger_perct = .FALSE.  !< update averaged input on call to bio_perct?
     134    LOGICAL ::  average_trigger_utci  = .FALSE.  !< update averaged input on call to bio_utci?
     135    LOGICAL ::  average_trigger_pet   = .FALSE.  !< update averaged input on call to bio_pet?
     136
     137    LOGICAL ::  bio_perct     = .TRUE.   !< Turn index PT (instant. input) on or off
     138    LOGICAL ::  bio_perct_av  = .TRUE.   !< Turn index PT (averaged input) on or off
     139    LOGICAL ::  bio_pet       = .TRUE.   !< Turn index PET (instant. input) on or off
     140    LOGICAL ::  bio_pet_av    = .TRUE.   !< Turn index PET (averaged input) on or off
     141    LOGICAL ::  bio_utci      = .TRUE.   !< Turn index UTCI (instant. input) on or off
     142    LOGICAL ::  bio_utci_av   = .TRUE.   !< Turn index UTCI (averaged input) on or off
     143
     144
     145!
     146!-- INTERFACES that must be available to other modules (alphabetical order)
     147
     148    PUBLIC bio_3d_data_averaging, bio_check_data_output,                       &
     149    bio_calculate_mrt_grid, bio_calculate_thermal_index_maps, bio_calc_ipt,    &
     150    bio_check_parameters, bio_data_output_3d, bio_data_output_2d,              &
     151    bio_define_netcdf_grid, bio_get_thermal_index_input_ij, bio_header,        &
     152    bio_init, bio_init_arrays, bio_parin, bio_perct, bio_perct_av, bio_pet,    &
     153    bio_pet_av, bio_utci, bio_utci_av, time_bio_results
    157154
    158155!
     
    160157!
    161158!-- 3D averaging for HTCM _INPUT_ variables
    162     INTERFACE biom_3d_data_averaging
    163        MODULE PROCEDURE biom_3d_data_averaging
    164     END INTERFACE biom_3d_data_averaging
     159    INTERFACE bio_3d_data_averaging
     160       MODULE PROCEDURE bio_3d_data_averaging
     161    END INTERFACE bio_3d_data_averaging
     162
     163!-- Calculate mtr from rtm fluxes and assign into 2D grid
     164    INTERFACE bio_calculate_mrt_grid
     165       MODULE PROCEDURE bio_calculate_mrt_grid
     166    END INTERFACE bio_calculate_mrt_grid
    165167
    166168!-- Calculate static thermal indices PT, UTCI and/or PET
    167     INTERFACE biom_calculate_static_grid
    168        MODULE PROCEDURE biom_calculate_static_grid
    169     END INTERFACE biom_calculate_static_grid
     169    INTERFACE bio_calculate_thermal_index_maps
     170       MODULE PROCEDURE bio_calculate_thermal_index_maps
     171    END INTERFACE bio_calculate_thermal_index_maps
    170172
    171173!-- Calculate the dynamic index iPT (to be caled by the agent model)
    172     INTERFACE biom_calc_ipt
    173        MODULE PROCEDURE biom_calc_ipt
    174     END INTERFACE biom_calc_ipt
     174    INTERFACE bio_calc_ipt
     175       MODULE PROCEDURE bio_calc_ipt
     176    END INTERFACE bio_calc_ipt
    175177
    176178!-- Data output checks for 2D/3D data to be done in check_parameters
    177      INTERFACE biom_check_data_output
    178         MODULE PROCEDURE biom_check_data_output
    179      END INTERFACE biom_check_data_output
     179     INTERFACE bio_check_data_output
     180        MODULE PROCEDURE bio_check_data_output
     181     END INTERFACE bio_check_data_output
    180182
    181183!-- Input parameter checks to be done in check_parameters
    182     INTERFACE biom_check_parameters
    183        MODULE PROCEDURE biom_check_parameters
    184     END INTERFACE biom_check_parameters
     184    INTERFACE bio_check_parameters
     185       MODULE PROCEDURE bio_check_parameters
     186    END INTERFACE bio_check_parameters
    185187
    186188!-- Data output of 2D quantities
    187     INTERFACE biom_data_output_2d
    188        MODULE PROCEDURE biom_data_output_2d
    189     END INTERFACE biom_data_output_2d
     189    INTERFACE bio_data_output_2d
     190       MODULE PROCEDURE bio_data_output_2d
     191    END INTERFACE bio_data_output_2d
    190192
    191193!-- no 3D data, thus, no averaging of 3D data, removed
    192     INTERFACE biom_data_output_3d
    193        MODULE PROCEDURE biom_data_output_3d
    194     END INTERFACE biom_data_output_3d
     194    INTERFACE bio_data_output_3d
     195       MODULE PROCEDURE bio_data_output_3d
     196    END INTERFACE bio_data_output_3d
    195197
    196198!-- Definition of data output quantities
    197     INTERFACE biom_define_netcdf_grid
    198        MODULE PROCEDURE biom_define_netcdf_grid
    199     END INTERFACE biom_define_netcdf_grid
     199    INTERFACE bio_define_netcdf_grid
     200       MODULE PROCEDURE bio_define_netcdf_grid
     201    END INTERFACE bio_define_netcdf_grid
    200202
    201203!-- Obtains all relevant input values to estimate local thermal comfort/stress
    202     INTERFACE biom_determine_input_at
    203        MODULE PROCEDURE biom_determine_input_at
    204     END INTERFACE biom_determine_input_at
     204    INTERFACE bio_get_thermal_index_input_ij
     205       MODULE PROCEDURE bio_get_thermal_index_input_ij
     206    END INTERFACE bio_get_thermal_index_input_ij
    205207
    206208!-- Output of information to the header file
    207     INTERFACE biom_header
    208        MODULE PROCEDURE biom_header
    209     END INTERFACE biom_header
     209    INTERFACE bio_header
     210       MODULE PROCEDURE bio_header
     211    END INTERFACE bio_header
    210212
    211213!-- Initialization actions
    212     INTERFACE biom_init
    213        MODULE PROCEDURE biom_init
    214     END INTERFACE biom_init
     214    INTERFACE bio_init
     215       MODULE PROCEDURE bio_init
     216    END INTERFACE bio_init
    215217
    216218!-- Initialization of arrays
    217     INTERFACE biom_init_arrays
    218        MODULE PROCEDURE biom_init_arrays
    219     END INTERFACE biom_init_arrays
     219    INTERFACE bio_init_arrays
     220       MODULE PROCEDURE bio_init_arrays
     221    END INTERFACE bio_init_arrays
    220222
    221223!-- Reading of NAMELIST parameters
    222     INTERFACE biom_parin
    223        MODULE PROCEDURE biom_parin
    224     END INTERFACE biom_parin
    225 
    226  
     224    INTERFACE bio_parin
     225       MODULE PROCEDURE bio_parin
     226    END INTERFACE bio_parin
     227
     228
    227229 CONTAINS
    228  
    229  
     230
     231
    230232!------------------------------------------------------------------------------!
    231233! Description:
     
    234236!> the array necessary for storing the average.
    235237!------------------------------------------------------------------------------!
    236  SUBROUTINE biom_3d_data_averaging( mode, variable )
     238 SUBROUTINE bio_3d_data_averaging( mode, variable )
    237239
    238240    IMPLICIT NONE
    239241
    240     CHARACTER (LEN=*) ::  mode     !<
    241     CHARACTER (LEN=*) ::  variable !<
    242 
    243     INTEGER(iwp) ::  i  !<
    244     INTEGER(iwp) ::  j  !<
    245     INTEGER(iwp) ::  k  !<
     242    CHARACTER (LEN=*) ::  mode     !< averaging mode: allocate, sum, or average
     243    CHARACTER (LEN=*) ::  variable !< The variable in question
     244
     245    INTEGER(iwp) ::  i        !< Running index, x-dir
     246    INTEGER(iwp) ::  j        !< Running index, y-dir
     247    INTEGER(iwp) ::  k        !< Running index, z-dir
    246248
    247249
     
    250252       SELECT CASE ( TRIM( variable ) )
    251253
    252           CASE ( 'biom_mrt' )
    253                 IF ( .NOT. ALLOCATED( biom_mrt_av ) )  THEN
    254                    ALLOCATE( biom_mrt_av(nmrtbl) )
     254          CASE ( 'bio_mrt' )
     255                IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
     256                   ALLOCATE( mrt_av_grid(nmrtbl) )
    255257                ENDIF
    256                 biom_mrt_av = 0.0_wp
    257 
    258           CASE ( 'biom_pt', 'biom_utci', 'biom_pet' )
     258                mrt_av_grid = 0.0_wp
     259
     260          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
    259261
    260262!--          Indices in unknown order as depending on input file, determine
    261263!            first index to average und update only once
    262              IF ( .NOT. average_trigger_pt .AND. .NOT. average_trigger_utci    &
     264             IF ( .NOT. average_trigger_perct .AND. .NOT. average_trigger_utci &
    263265                .AND. .NOT. average_trigger_pet ) THEN
    264                 IF ( TRIM( variable ) == 'biom_pt' ) THEN
    265                     average_trigger_pt = .TRUE.
     266                IF ( TRIM( variable ) == 'bio_perct*' ) THEN
     267                    average_trigger_perct = .TRUE.
    266268                ENDIF
    267                 IF ( TRIM( variable ) == 'biom_utci' ) THEN
     269                IF ( TRIM( variable ) == 'bio_utci*' ) THEN
    268270                    average_trigger_utci = .TRUE.
    269271                ENDIF
    270                 IF ( TRIM( variable ) == 'biom_pet' ) THEN
     272                IF ( TRIM( variable ) == 'bio_pet*' ) THEN
    271273                    average_trigger_pet = .TRUE.
    272274                ENDIF
     
    274276
    275277!--          Only continue if updateindex
    276              IF ( average_trigger_pt   .AND. TRIM( variable ) /= 'biom_pt')  RETURN
    277              IF ( average_trigger_utci .AND. TRIM( variable ) /= 'biom_utci') RETURN
    278              IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'biom_pet')  RETURN
     278             IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') RETURN
     279             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*')  RETURN
     280             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'bio_pet*')    RETURN
    279281
    280282!--          Set averaging switch to .TRUE. if not done by other module before!
    281283             IF ( .NOT. ALLOCATED( pt_av ) )  THEN
    282284                ALLOCATE( pt_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
    283                 aver_pt = .TRUE.
     285                aver_perct = .TRUE.
    284286             ENDIF
    285287             pt_av = 0.0_wp
     
    292294
    293295             IF ( .NOT. ALLOCATED( u_av ) )  THEN
    294                 ALLOCATE( u_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     296                ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    295297                aver_u = .TRUE.
    296298             ENDIF
     
    298300
    299301             IF ( .NOT. ALLOCATED( v_av ) )  THEN
    300                 ALLOCATE( v_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     302                ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    301303                aver_v = .TRUE.
    302304             ENDIF
     
    304306
    305307             IF ( .NOT. ALLOCATED( w_av ) )  THEN
    306                 ALLOCATE( w_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
     308                ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    307309                aver_w = .TRUE.
    308310             ENDIF
     
    318320       SELECT CASE ( TRIM( variable ) )
    319321
    320           CASE ( 'biom_mrt' )
    321              IF ( ALLOCATED( biom_mrt_av ) )  THEN
    322 
    323                  IF ( nmrtbl > 0 )  THEN
    324                     IF ( mrt_include_sw )  THEN
    325                        biom_mrt_av(:) = biom_mrt_av(:) + &
    326                           ((human_absorb*mrtinsw(:) + human_emiss*mrtinlw(:))  &
    327                           / (human_emiss*sigma_sb)) ** .25_wp - degc_to_k
    328                     ELSE
    329                        biom_mrt_av(:) = biom_mrt_av(:) + &
    330                           (human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp      &
    331                           - degc_to_k
    332                     ENDIF
    333                  ENDIF
     322          CASE ( 'bio_mrt' )
     323             IF ( ALLOCATED( mrt_av_grid ) )  THEN
     324
     325                IF ( mrt_include_sw )  THEN
     326                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
     327                      (( human_absorb * mrtinsw(:) + human_emiss * mrtinlw(:)) &
     328                      / (human_emiss * sigma_sb)) ** .25_wp - degc_to_k
     329                ELSE
     330                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
     331                      (human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp          &
     332                      - degc_to_k
     333                ENDIF
    334334             ENDIF
    335335
    336           CASE ( 'biom_pt', 'biom_utci', 'biom_pet' )
     336          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
    337337
    338338!--          Only continue if updateindex
    339              IF ( average_trigger_pt   .AND. TRIM( variable ) /= 'biom_pt')    &
     339             IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') &
    340340                RETURN
    341              IF ( average_trigger_utci .AND. TRIM( variable ) /= 'biom_utci')  &
     341             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*')  &
    342342                RETURN
    343              IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'biom_pet')   &
     343             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'bio_pet*')   &
    344344                RETURN
    345345
    346              IF ( ALLOCATED( pt_av ) .AND. aver_pt ) THEN
     346             IF ( ALLOCATED( pt_av ) .AND. aver_perct ) THEN
    347347                DO  i = nxl, nxr
    348348                   DO  j = nys, nyn
     
    365365
    366366             IF ( ALLOCATED( u_av )  .AND. aver_u ) THEN
    367                 DO  i = nxl, nxr
    368                    DO  j = nys, nyn
     367                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
     368                   DO  j = nysg, nyng
    369369                      DO  k = nzb, nzt+1
    370370                         u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
     
    375375
    376376             IF ( ALLOCATED( v_av )  .AND. aver_v ) THEN
    377                 DO  i = nxl, nxr
    378                    DO  j = nys, nyn
     377                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
     378                   DO  j = nysg, nyng
    379379                      DO  k = nzb, nzt+1
    380380                         v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
     
    385385
    386386             IF ( ALLOCATED( w_av )  .AND. aver_w ) THEN
    387                 DO  i = nxl, nxr
    388                    DO  j = nys, nyn
     387                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
     388                   DO  j = nysg, nyng
    389389                      DO  k = nzb, nzt+1
    390390                         w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
     
    403403       SELECT CASE ( TRIM( variable ) )
    404404
    405           CASE ( 'biom_mrt' )
    406              IF ( ALLOCATED( biom_mrt_av ) )  THEN
    407                 biom_mrt_av(:) = biom_mrt_av(:) / REAL( average_count_3d, KIND=wp )
     405          CASE ( 'bio_mrt' )
     406             IF ( ALLOCATED( mrt_av_grid ) )  THEN
     407                mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp )
    408408             ENDIF
    409409
    410           CASE ( 'biom_pt', 'biom_utci', 'biom_pet' )
     410          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
    411411
    412412!--          Only continue if update index
    413              IF ( average_trigger_pt   .AND. TRIM( variable ) /= 'biom_pt')    &
     413             IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') &
    414414                RETURN
    415              IF ( average_trigger_utci .AND. TRIM( variable ) /= 'biom_utci')  &
     415             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*')  &
    416416                RETURN
    417              IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'biom_pet')   &
     417             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'bio_pet*')   &
    418418                RETURN
    419419
    420              IF ( ALLOCATED( pt_av ) .AND. aver_pt ) THEN
     420             IF ( ALLOCATED( pt_av ) .AND. aver_perct ) THEN
    421421                DO  i = nxl, nxr
    422422                   DO  j = nys, nyn
     
    439439
    440440             IF ( ALLOCATED( u_av ) .AND. aver_u ) THEN
    441                 DO  i = nxl, nxr
    442                    DO  j = nys, nyn
     441                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
     442                   DO  j = nysg, nyng
    443443                      DO  k = nzb, nzt+1
    444444                         u_av(k,j,i) = u_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     
    449449
    450450             IF ( ALLOCATED( v_av ) .AND. aver_v ) THEN
    451                 DO  i = nxl, nxr
    452                    DO  j = nys, nyn
     451                DO  i = nxlg, nxrg
     452                   DO  j = nysg, nyng
    453453                      DO  k = nzb, nzt+1
    454454                         v_av(k,j,i) = v_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     
    459459
    460460             IF ( ALLOCATED( w_av ) .AND. aver_w ) THEN
    461                 DO  i = nxl, nxr
    462                    DO  j = nys, nyn
     461                DO  i = nxlg, nxrg
     462                   DO  j = nysg, nyng
    463463                      DO  k = nzb, nzt+1
    464464                         w_av(k,j,i) = w_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     
    469469
    470470!--          Udate thermal indices with derived averages
    471              CALL biom_calculate_static_grid ( .TRUE. )
     471             CALL bio_calculate_thermal_index_maps ( .TRUE. )
    472472
    473473        END SELECT
     
    476476
    477477
    478  END SUBROUTINE biom_3d_data_averaging
     478 END SUBROUTINE bio_3d_data_averaging
    479479
    480480
     
    485485!> Check data output for biometeorology model
    486486!------------------------------------------------------------------------------!
    487  SUBROUTINE biom_check_data_output( var, unit )
     487 SUBROUTINE bio_check_data_output( var, unit )
    488488
    489489    USE control_parameters,                                                    &
     
    498498    SELECT CASE ( TRIM( var ) )
    499499     
    500        CASE( 'biom_tmrt', 'biom_mrt', 'biom_pet', 'biom_pt', 'biom_utci' )
     500       CASE( 'bio_mrt', 'bio_pet*', 'bio_perct*', 'bio_utci*' )
    501501          unit = 'degree_C'
    502502
     
    509509!
    510510!--    Futher checks if output belongs to biometeorology
    511        IF ( .NOT. biometeorology .OR. .NOT. radiation ) THEN
    512          message_string = 'output of "' // TRIM( var ) // '" req' //  &
    513                'uires biometeorology = .TRUE. and radiation = .TRUE. ' &
    514                // 'in initialisation_parameters'
    515          CALL message( 'biom_check_data_output', 'PA0561', 1, 2, 0, 6, 0 )
     511       IF ( .NOT. radiation ) THEN
     512          message_string = 'output of "' // TRIM( var ) // '" require'         &
     513                           // 's radiation = .TRUE.'
     514          CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
     515          unit = 'illegal'
    516516       ENDIF
    517517       IF ( mrt_nlevels == 0 ) THEN
    518           message_string = 'output of "' // TRIM( var ) // '" require'&
     518          message_string = 'output of "' // TRIM( var ) // '" require'         &
    519519                           // 's mrt_nlevels > 0'
    520           CALL message( 'biom_check_data_output', 'PA0562', 1, 2, 0, 6, 0 )
     520          CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
     521          unit = 'illegal'
    521522       ENDIF
    522523
    523     ENDIF
    524 
    525  END SUBROUTINE biom_check_data_output
     524
     525    ENDIF
     526
     527 END SUBROUTINE bio_check_data_output
    526528
    527529!------------------------------------------------------------------------------!
     
    531533!> check_parameters 1302
    532534!------------------------------------------------------------------------------!
    533  SUBROUTINE biom_check_parameters
     535 SUBROUTINE bio_check_parameters
    534536
    535537    USE control_parameters,                                                    &
     
    541543
    542544!--    Disabled as long as radiation model not available
    543        IF ( .NOT. radiation )  THEN
    544           message_string = 'The human thermal comfort module does require ' // &
    545                            'radiation information in terms of the mean ' //    &
    546                            'radiant temperature, but radiation is not enabled!'
    547           CALL message( 'check_parameters', 'PAHU01', 0, 0, 0, 6, 0 )
    548        ENDIF
    549545
    550546       IF ( .NOT. humidity )  THEN
    551           message_string = 'The human thermal comfort module does require ' // &
     547          message_string = 'The estimation of thermal comfort requires '    // &
    552548                           'air humidity information, but humidity module ' // &
    553549                           'is disabled!'
    554           CALL message( 'check_parameters', 'PAHU02', 0, 0, 0, 6, 0 )
     550          CALL message( 'check_parameters', 'PA0561', 0, 0, 0, 6, 0 )
    555551       ENDIF
    556552
    557  END SUBROUTINE biom_check_parameters
    558 
    559 
    560 !------------------------------------------------------------------------------!
    561 !
    562 ! Description:
    563 ! ------------
    564 !> Subroutine defining 3D output variables (dummy, always 2d!)
    565 !> data_output_3d 709ff
    566 !------------------------------------------------------------------------------!
    567  SUBROUTINE biom_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
    568 
    569     USE indices
    570 
    571     USE kinds
    572 
    573 
    574     IMPLICIT NONE
    575 
    576 !-- Input variables
    577     CHARACTER (LEN=*), INTENT(IN) ::  variable   !< Char identifier to select var for output
    578     INTEGER(iwp), INTENT(IN) ::  av       !< Use averaged data? 0 = no, 1 = yes?
    579     INTEGER(iwp), INTENT(IN) ::  nzb_do   !< Unused. 2D. nz bottom to nz top
    580     INTEGER(iwp), INTENT(IN) ::  nzt_do   !< Unused.
    581 
    582 !-- Output variables
    583     LOGICAL, INTENT(OUT) ::  found   !< Output found?
    584     REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< Temp. result grid to return
    585 
    586 !-- Internal variables
    587     INTEGER(iwp) ::  l    !< Running index, radiation grid
    588     INTEGER(iwp) ::  i    !< Running index, x-dir
    589     INTEGER(iwp) ::  j    !< Running index, y-dir
    590     INTEGER(iwp) ::  k    !< Running index, z-dir
    591 
    592     CHARACTER (LEN=:), allocatable ::  variable_short  !< Trimmed version of char variable
    593 
    594     REAL(wp), PARAMETER ::  fill_value = -999._wp
    595     REAL(wp) ::  mrt  !< Buffer for mean radiant temperature
    596 
    597     found = .TRUE.
    598     variable_short = TRIM( variable )
    599 
    600     IF ( variable_short(1:4) /= 'biom' ) THEN
    601 !--   Nothing to do, set found to FALSE and return immediatelly
    602       found = .FALSE.
    603       RETURN
    604     ENDIF
    605 
    606     SELECT CASE ( variable_short )
    607 
    608         CASE ( 'biom_mrt' )
    609 
    610             local_pf = REAL( fill_value, KIND = wp )
    611             DO  l = 1, nmrtbl
    612                 i = mrtbl(ix,l)
    613                 j = mrtbl(iy,l)
    614                 k = mrtbl(iz,l)
    615                 IF ( mrt_include_sw )  THEN
    616                     mrt = ((human_absorb*mrtinsw(l) + human_emiss*mrtinlw(l)) &
    617                            / (human_emiss*sigma_sb)) ** .25_wp
    618                 ELSE
    619                     mrt = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp
    620                 ENDIF
    621                 local_pf(i,j,k) = mrt
    622             ENDDO
    623 
    624         CASE ( 'biom_tmrt' )        ! 2d-array
    625               DO  i = nxl, nxr
    626                  DO  j = nys, nyn
    627                     local_pf(i,j,nzb_do) = REAL( tmrt_grid(j,i), sp )
    628                  ENDDO
    629               ENDDO
    630 
    631         CASE ( 'biom_pt' )        ! 2d-array
    632             IF ( av == 0 )  THEN
    633               DO  i = nxl, nxr
    634                  DO  j = nys, nyn
    635                     local_pf(i,j,nzb_do) = REAL( pt_grid(j,i), sp )
    636                  ENDDO
    637               ENDDO
    638             ELSE
    639               DO  i = nxl, nxr
    640                  DO  j = nys, nyn
    641                     local_pf(i,j,nzb_do) = REAL( pt_av_grid(j,i), sp )
    642                  ENDDO
    643               ENDDO
    644             END IF
    645 
    646         CASE ( 'biom_utci' )        ! 2d-array
    647             IF ( av == 0 )  THEN
    648               DO  i = nxl, nxr
    649                  DO  j = nys, nyn
    650                     local_pf(i,j,nzb_do) = REAL( utci_grid(j,i), sp )
    651                  ENDDO
    652               ENDDO
    653             ELSE
    654               DO  i = nxl, nxr
    655                  DO  j = nys, nyn
    656                     local_pf(i,j,nzb_do) = REAL( utci_av_grid(j,i), sp )
    657                  ENDDO
    658               ENDDO
    659             END IF
    660 
    661         CASE ( 'biom_pet' )        ! 2d-array
    662             IF ( av == 0 )  THEN
    663               DO  i = nxl, nxr
    664                  DO  j = nys, nyn
    665                     local_pf(i,j,nzb_do) = REAL( pet_grid(j,i), sp )
    666                  ENDDO
    667               ENDDO
    668             ELSE
    669               DO  i = nxl, nxr
    670                  DO  j = nys, nyn
    671                     local_pf(i,j,nzb_do) = REAL( pet_av_grid(j,i), sp )
    672                  ENDDO
    673               ENDDO
    674             END IF
    675 
    676        CASE DEFAULT
    677           found = .FALSE.
    678 
    679     END SELECT
    680 
    681  END SUBROUTINE biom_data_output_3d
     553 END SUBROUTINE bio_check_parameters
     554
    682555
    683556!------------------------------------------------------------------------------!
     
    688561!> data_output_2d 1188ff
    689562!------------------------------------------------------------------------------!
    690  SUBROUTINE biom_data_output_2d( av, variable, found, grid, local_pf,          &
     563 SUBROUTINE bio_data_output_2d( av, variable, found, grid, local_pf,           &
    691564                                      two_d, nzb_do, nzt_do, fill_value )
    692565
     
    717590    INTEGER(iwp) ::  j        !< Running index, y-dir
    718591    INTEGER(iwp) ::  k        !< Running index, z-dir
    719 
    720 
    721     found = .TRUE.
     592    INTEGER(iwp) ::  l        !< Running index, radiation grid
     593
     594
    722595    variable_short = TRIM( variable )
    723     IF ( variable_short(1:4) == 'biom' ) THEN
    724        two_d = .TRUE.
    725        grid = 'zu1'
    726     ELSE
     596    IF ( variable_short(1:4) /= 'bio_' ) THEN
    727597       found = .FALSE.
    728598       grid  = 'none'
    729        RETURN
    730     ENDIF
    731 
     599    ENDIF
     600
     601    found = .TRUE.
    732602    local_pf = fill_value
    733603
     
    735605
    736606
    737         CASE ( 'biom_tmrt_xy' )        ! 2d-array
    738               DO  i = nxl, nxr
    739                  DO  j = nys, nyn
    740                     local_pf(i,j,1) = tmrt_grid(j,i)
    741                  ENDDO
    742               ENDDO
    743 
    744 
    745         CASE ( 'biom_pt_xy' )        ! 2d-array
     607        CASE ( 'bio_mrt_xy' )
     608            grid = 'zu1'
     609            two_d = .FALSE.  !< can be calculated for several levels
     610            local_pf = REAL( fill_value, KIND = wp )
     611            DO  l = 1, nmrtbl
     612               i = mrtbl(ix,l)
     613               j = mrtbl(iy,l)
     614               k = mrtbl(iz,l)
     615               IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. j > nyn .OR.  &
     616                  i < nxl .OR. i > nxr ) CYCLE
     617               IF ( av == 0 )  THEN
     618                  IF ( mrt_include_sw )  THEN
     619                     local_pf(i,j,k) = ((human_absorb * mrtinsw(l) +           &
     620                     human_emiss * mrtinlw(l)) /                               &
     621                     (human_emiss * sigma_sb)) ** .25_wp - degc_to_k
     622                  ELSE
     623                     local_pf(i,j,k) = (human_emiss * mrtinlw(l) /             &
     624                                        sigma_sb) ** .25_wp - degc_to_k
     625                  ENDIF
     626               ELSE
     627                  local_pf(i,j,k) = mrt_av_grid(l)
     628               ENDIF
     629            ENDDO
     630
     631
     632        CASE ( 'bio_perct*_xy' )        ! 2d-array
     633            grid = 'zu1'
     634            two_d = .TRUE.
    746635            IF ( av == 0 )  THEN
    747636              DO  i = nxl, nxr
    748637                 DO  j = nys, nyn
    749                     local_pf(i,j,nzb+1) = pt_grid(j,i)
     638                    local_pf(i,j,nzb+1) = perct(j,i)
    750639                 ENDDO
    751640              ENDDO
     
    753642              DO  i = nxl, nxr
    754643                 DO  j = nys, nyn
    755                     local_pf(i,j,nzb+1) = pt_av_grid(j,i)
     644                    local_pf(i,j,nzb+1) = perct_av(j,i)
    756645                 ENDDO
    757646              ENDDO
     
    759648
    760649
    761         CASE ( 'biom_utci_xy' )        ! 2d-array
     650        CASE ( 'bio_utci*_xy' )        ! 2d-array
     651            grid = 'zu1'
     652            two_d = .TRUE.
    762653            IF ( av == 0 )  THEN
    763654              DO  i = nxl, nxr
     
    775666
    776667
    777         CASE ( 'biom_pet_xy' )        ! 2d-array
     668        CASE ( 'bio_pet*_xy' )        ! 2d-array
     669            grid = 'zu1'
     670            two_d = .TRUE.
    778671            IF ( av == 0 )  THEN
    779672              DO  i = nxl, nxr
     
    798691
    799692
    800  END SUBROUTINE biom_data_output_2d
    801 
     693 END SUBROUTINE bio_data_output_2d
     694
     695
     696!------------------------------------------------------------------------------!
     697!
     698! Description:
     699! ------------
     700!> Subroutine defining 3D output variables (dummy, always 2d!)
     701!> data_output_3d 709ff
     702!------------------------------------------------------------------------------!
     703 SUBROUTINE bio_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
     704
     705    USE indices
     706
     707    USE kinds
     708
     709
     710    IMPLICIT NONE
     711
     712!-- Input variables
     713    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< Char identifier to select var for output
     714    INTEGER(iwp), INTENT(IN) ::  av       !< Use averaged data? 0 = no, 1 = yes?
     715    INTEGER(iwp), INTENT(IN) ::  nzb_do   !< Unused. 2D. nz bottom to nz top
     716    INTEGER(iwp), INTENT(IN) ::  nzt_do   !< Unused.
     717
     718!-- Output variables
     719    LOGICAL, INTENT(OUT) ::  found   !< Output found?
     720    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< Temp. result grid to return
     721
     722!-- Internal variables
     723    INTEGER(iwp) ::  l    !< Running index, radiation grid
     724    INTEGER(iwp) ::  i    !< Running index, x-dir
     725    INTEGER(iwp) ::  j    !< Running index, y-dir
     726    INTEGER(iwp) ::  k    !< Running index, z-dir
     727
     728    CHARACTER (LEN=:), allocatable ::  variable_short  !< Trimmed version of char variable
     729
     730    REAL(wp), PARAMETER ::  fill_value = -999._wp
     731    REAL(wp) ::  mrt  !< Buffer for mean radiant temperature
     732
     733    found = .TRUE.
     734    variable_short = TRIM( variable )
     735
     736    IF ( variable_short(1:4) /= 'bio_' ) THEN
     737!--   Nothing to do, set found to FALSE and return immediatelly
     738      found = .FALSE.
     739      RETURN
     740    ENDIF
     741
     742    SELECT CASE ( variable_short )
     743
     744        CASE ( 'bio_mrt' )
     745            local_pf = REAL( fill_value, KIND = wp )
     746            DO  l = 1, nmrtbl
     747               i = mrtbl(ix,l)
     748               j = mrtbl(iy,l)
     749               k = mrtbl(iz,l)
     750               IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. j > nyn .OR.  &
     751                  i < nxl .OR. i > nxr ) CYCLE
     752               IF ( av == 0 )  THEN
     753                  IF ( mrt_include_sw )  THEN
     754                     local_pf(i,j,k) = ((human_absorb * mrtinsw(l) + human_emiss * mrtinlw(l)) /   &
     755                                        (human_emiss * sigma_sb)) ** .25_wp - degc_to_k
     756                  ELSE
     757                     local_pf(i,j,k) = (human_emiss * mrtinlw(l) /             &
     758                                        sigma_sb) ** .25_wp - degc_to_k
     759                  ENDIF
     760               ELSE
     761                  local_pf(i,j,k) = mrt_av_grid(l)
     762               ENDIF
     763            ENDDO
     764
     765       CASE DEFAULT
     766          found = .FALSE.
     767
     768    END SELECT
     769
     770 END SUBROUTINE bio_data_output_3d
    802771
    803772!------------------------------------------------------------------------------!
     
    808777!> netcdf_interface_mod 918ff
    809778!------------------------------------------------------------------------------!
    810  SUBROUTINE biom_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
     779 SUBROUTINE bio_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
    811780
    812781    IMPLICIT NONE
     
    836805    is2d = ( var(l-1:l) == 'xy' )
    837806
    838 
    839     IF ( var(1:4) == 'biom' ) THEN
     807    IF ( var(1:4) == 'bio_' ) THEN
    840808       found  = .TRUE.
    841809       grid_x = 'x'
     
    845813    ENDIF
    846814
    847  END SUBROUTINE biom_define_netcdf_grid
     815 END SUBROUTINE bio_define_netcdf_grid
    848816
    849817!------------------------------------------------------------------------------!
     
    853821!> header 982
    854822!------------------------------------------------------------------------------!
    855  SUBROUTINE biom_header( io )
     823 SUBROUTINE bio_header( io )
    856824
    857825    IMPLICIT NONE
     
    863831    CHARACTER (LEN=86) ::  output_height_chr  !< String for output height
    864832
    865     WRITE( output_height_chr, '(F8.1,7X)' )  biom_output_height
     833    WRITE( output_height_chr, '(F8.1,7X)' )  bio_output_height
    866834!
    867835!-- Write biom header
    868836    WRITE( io, 1 )
    869837    WRITE( io, 2 )  TRIM( output_height_chr )
    870     WRITE( io, 3 )  TRIM( ACHAR( biom_cell_level ) )
     838    WRITE( io, 3 )  TRIM( ACHAR( bio_cell_level ) )
    871839
    8728401   FORMAT (//' Human thermal comfort module information:'/                    &
     
    8758433   FORMAT ('    --> This corresponds to cell level : ', A )
    876844
    877  END SUBROUTINE biom_header
     845 END SUBROUTINE bio_header
    878846
    879847
     
    884852!> init_3d_model 1987ff
    885853!------------------------------------------------------------------------------!
    886  SUBROUTINE biom_init
    887 
    888     USE control_parameters,                                                 &
     854 SUBROUTINE bio_init
     855
     856    USE control_parameters,                                                    &
    889857        ONLY: message_string
    890858
     
    899867!   (gravimetric center of sample human)
    900868
    901     time_biom_results = 0.0_wp
    902     biom_cell_level = 0_iwp
    903     biom_output_height = 0.5_wp * dz(1)
     869    time_bio_results = 0.0_wp
     870    bio_cell_level = 0_iwp
     871    bio_output_height = 0.5_wp * dz(1)
    904872    height = 0.0_wp
    905873
    906     biom_cell_level = INT ( 1.099_wp / dz(1) )
    907     biom_output_height = biom_output_height + biom_cell_level * dz(1)
    908 
    909     IF ( .NOT. radiation_interactions .AND. mrt_from_rtm ) THEN
    910        message_string = 'The mrt_from_rtm switch require ' // &
     874    bio_cell_level = INT ( 1.099_wp / dz(1) )
     875    bio_output_height = bio_output_height + bio_cell_level * dz(1)
     876
     877    IF ( .NOT. radiation_interactions ) THEN
     878       message_string = 'The mrt calculation requires ' // &
    911879                        'enabled radiation_interactions but it ' // &
    912                         'is disabled! Set mrt_from_rtm to .F.'
     880                        'is disabled!'
    913881       CALL message( 'check_parameters', 'PAHU03', 0, 0, -1, 6, 0 )
    914        mrt_from_rtm = .FALSE.
    915     ENDIF
    916 
    917  END SUBROUTINE biom_init
     882    ENDIF
     883
     884 END SUBROUTINE bio_init
    918885
    919886
     
    924891!> init_3d_model 1047ff
    925892!------------------------------------------------------------------------------!
    926  SUBROUTINE biom_init_arrays
     893 SUBROUTINE bio_init_arrays
    927894
    928895    IMPLICIT NONE
     
    934901    ENDIF
    935902
    936     IF ( biom_pt ) THEN
    937       IF ( .NOT. ALLOCATED( pt_grid ) ) THEN
    938         ALLOCATE( pt_grid (nys:nyn,nxl:nxr) )
     903    IF ( bio_perct ) THEN
     904      IF ( .NOT. ALLOCATED( perct ) ) THEN
     905        ALLOCATE( perct (nys:nyn,nxl:nxr) )
    939906      ENDIF
    940907    ENDIF
    941908
    942     IF ( biom_utci ) THEN
     909    IF ( bio_utci ) THEN
    943910      IF ( .NOT. ALLOCATED( utci_grid ) ) THEN
    944911        ALLOCATE( utci_grid (nys:nyn,nxl:nxr) )
     
    946913    ENDIF
    947914
    948     IF ( biom_pet ) THEN
     915    IF ( bio_pet ) THEN
    949916      IF ( .NOT. ALLOCATED( pet_grid ) ) THEN
    950917        ALLOCATE( pet_grid (nys:nyn,nxl:nxr) )
     
    952919    END IF
    953920
    954     IF ( biom_pt_av ) THEN
    955       IF ( .NOT. ALLOCATED( pt_av_grid ) ) THEN
    956         ALLOCATE( pt_av_grid (nys:nyn,nxl:nxr) )
     921    IF ( bio_perct_av ) THEN
     922      IF ( .NOT. ALLOCATED( perct_av ) ) THEN
     923        ALLOCATE( perct_av (nys:nyn,nxl:nxr) )
    957924      ENDIF
    958925    ENDIF
    959926
    960     IF ( biom_utci_av ) THEN
     927    IF ( bio_utci_av ) THEN
    961928      IF ( .NOT. ALLOCATED( utci_av_grid ) ) THEN
    962929        ALLOCATE( utci_av_grid (nys:nyn,nxl:nxr) )
     
    964931    ENDIF
    965932
    966     IF ( biom_pet_av ) THEN
     933    IF ( bio_pet_av ) THEN
    967934      IF ( .NOT. ALLOCATED( pet_av_grid ) ) THEN
    968935        ALLOCATE( pet_av_grid (nys:nyn,nxl:nxr) )
     
    971938    END IF
    972939
    973  END SUBROUTINE biom_init_arrays
     940 END SUBROUTINE bio_init_arrays
    974941
    975942
     
    979946!> Parin for &biometeorology_parameters for reading biomet parameters
    980947!------------------------------------------------------------------------------!
    981  SUBROUTINE biom_parin
     948 SUBROUTINE bio_parin
    982949
    983950    IMPLICIT NONE
     
    987954    CHARACTER (LEN=80) ::  line  !< Dummy string for current line in parameter file
    988955
    989     NAMELIST /biometeorology_parameters/  mrt_from_rtm,                        &
    990                                           biom_pet,                            &
    991                                           biom_pet_av,                         &
    992                                           biom_pt,                             &
    993                                           biom_pt_av,                          &
    994                                           biom_utci,                           &
    995                                           biom_utci_av
     956    NAMELIST /biometeorology_parameters/  bio_pet,                             &
     957                                          bio_pet_av,                          &
     958                                          bio_perct,                           &
     959                                          bio_perct_av,                        &
     960                                          bio_utci,                            &
     961                                          bio_utci_av
    996962
    997963
     
    1025991
    1026992
    1027  END SUBROUTINE biom_parin
     993 END SUBROUTINE bio_parin
    1028994
    1029995!------------------------------------------------------------------------------!
    1030996! Description:
    1031997! ------------
    1032 !> Calculates the mean radiant temperature (tmrt) based on the Six-directions
    1033 !> method according to VDI 3787 2.
    1034 !------------------------------------------------------------------------------!
    1035  SUBROUTINE calculate_tmrt_6_directions( SW_N, SW_E, SW_S, SW_W,               &
    1036     SW_U, SW_D, LW_N, LW_E, LW_S, LW_W, LW_U, LW_D, tmrt )
     998!> Calculate biometeorology MRT for all 2D grid
     999!------------------------------------------------------------------------------!
     1000 SUBROUTINE bio_calculate_mrt_grid ( av )
    10371001
    10381002    IMPLICIT NONE
    10391003
    1040 !-- Type of input of the argument list
    1041 !   Short- (SW_) and longwave (LW_) radiation fluxes from the six directions
    1042 !   North (N), East (E), South (S), West (W), up (U) and down (D)
    1043     REAL(wp), INTENT ( IN )  ::  SW_N  !< Sw radflux density from N    (W/m²)
    1044     REAL(wp), INTENT ( IN )  ::  SW_E  !< Sw radflux density from E    (W/m²)
    1045     REAL(wp), INTENT ( IN )  ::  SW_S  !< Sw radflux density from S    (W/m²)
    1046     REAL(wp), INTENT ( IN )  ::  SW_W  !< Sw radflux density from W    (W/m²)
    1047     REAL(wp), INTENT ( IN )  ::  SW_U  !< Sw radflux density from U    (W/m²)
    1048     REAL(wp), INTENT ( IN )  ::  SW_D  !< Sw radflux density from D    (W/m²)
    1049     REAL(wp), INTENT ( IN )  ::  LW_N  !< Lw radflux density from N    (W/m²)
    1050     REAL(wp), INTENT ( IN )  ::  LW_E  !< Lw radflux density from E    (W/m²)
    1051     REAL(wp), INTENT ( IN )  ::  LW_S  !< Lw radflux density from S    (W/m²)
    1052     REAL(wp), INTENT ( IN )  ::  LW_W  !< Lw radflux density from W    (W/m²)
    1053     REAL(wp), INTENT ( IN )  ::  LW_U  !< Lw radflux density from U    (W/m²)
    1054     REAL(wp), INTENT ( IN )  ::  LW_D  !< Lw radflux density from D    (W/m²)
    1055 
    1056 !-- Type of output of the argument list
    1057     REAL(wp), INTENT ( OUT ) ::  tmrt  !< Mean radiant temperature (°C)
    1058 
    1059 !-- Directional weighting factors
    1060     REAL(wp), PARAMETER      ::  weight_h  = 0.22_wp
    1061     REAL(wp), PARAMETER      ::  weight_v  = 0.06_wp
    1062 
    1063     REAL(wp) ::  nrfd           !<  Net radiation flux density (W/m²)
    1064 
    1065 !-- Initialise
    1066     nrfd = 0._wp
    1067 
    1068 !-- Compute mean radiation flux density absorbed by rotational symetric human
    1069     nrfd = ( weight_h * ( human_absorb * SW_N + human_emiss * LW_N ) ) +       &
    1070        ( weight_h * ( human_absorb * SW_E + human_emiss * LW_E ) ) +           &
    1071        ( weight_h * ( human_absorb * SW_S + human_emiss * LW_S ) ) +           &
    1072        ( weight_h * ( human_absorb * SW_W + human_emiss * LW_W ) ) +           &
    1073        ( weight_v * ( human_absorb * SW_U + human_emiss * LW_U ) ) +           &
    1074        ( weight_v * ( human_absorb * SW_D + human_emiss * LW_D ) )
    1075 
    1076 !-- Compute mean radiant temperature
    1077     tmrt = ( nrfd / (human_emiss * sigma_sb) )**0.25_wp - degc_to_k
    1078 
    1079  END SUBROUTINE calculate_tmrt_6_directions
    1080 
    1081 !------------------------------------------------------------------------------!
    1082 ! Description:
    1083 ! ------------
    1084 !> Very crude approximation of mean radiant temperature based on upwards and
    1085 !> downwards radiation fluxes only (other directions curently not available,
    1086 !> replace as soon as possible!)
    1087 !------------------------------------------------------------------------------!
    1088  SUBROUTINE calculate_tmrt_2_directions( sw_u, sw_d, lw_u, lw_d, ta, tmrt )
    1089 
    1090     IMPLICIT NONE
    1091 
    1092 !-- Type of input of the argument list
    1093     REAL(wp), INTENT ( IN )  ::  sw_u  !< Shortwave radiation flux density from upper direction (W/m²)
    1094     REAL(wp), INTENT ( IN )  ::  sw_d  !< Shortwave radiation flux density from lower direction (W/m²)
    1095     REAL(wp), INTENT ( IN )  ::  lw_u  !< Longwave radiation flux density from upper direction (W/m²)
    1096     REAL(wp), INTENT ( IN )  ::  lw_d  !< Longwave radiation flux density from lower direction (W/m²)
    1097     REAL(wp), INTENT ( IN )  ::  ta    !< Air temperature (°C)
    1098 
    1099 !-- Type of output of the argument list
    1100     REAL(wp), INTENT ( OUT ) ::  tmrt  !< mean radiant temperature, (°C)
    1101 
    1102 !-- Directional weighting factors and parameters
    1103     REAL(wp), PARAMETER ::  weight_h  = 0.22_wp     !< Weight for horizontal radiational gain after Fanger (1972)
    1104     REAL(wp), PARAMETER ::  weight_v  = 0.06_wp     !< Weight for vertical radiational gain after Fanger (1972)
    1105 
    1106 !-- Other internal variables
    1107     REAL(wp) ::  sw_in
    1108     REAL(wp) ::  sw_out
    1109     REAL(wp) ::  lw_in
    1110     REAL(wp) ::  lw_out
    1111     REAL(wp) ::  nrfd     !< Net radiation flux density (W/m²)
    1112     REAL(wp) ::  lw_air   !< Longwave emission by surrounding air volume (W/m²)
    1113     REAL(wp) ::  sw_side  !< Shortwave immission from the sides (W/m²)
    1114 
    1115     INTEGER(iwp) ::  no_input  !< Count missing input radiation fluxes
    1116 
    1117 !-- initialise
    1118     sw_in    = sw_u
    1119     sw_out   = sw_d
    1120     lw_in    = lw_u
    1121     lw_out   = lw_d
    1122     nrfd     = 0._wp
    1123     no_input = 0_iwp
    1124 
    1125 !-- test for missing input data
    1126     IF ( sw_in <= -998._wp .OR. sw_out <= -998._wp .OR. lw_in <= -998._wp .OR. &
    1127        lw_out <= -998._wp .OR. ta <= -998._wp ) THEN
    1128        IF ( sw_in <= -998._wp ) THEN
    1129           sw_in = 0.
    1130           no_input = no_input + 1
     1004    LOGICAL, INTENT(IN)         ::  av    !< use averaged input?
     1005!-- Internal variables
     1006    INTEGER(iwp)                ::  i     !< Running index, x-dir, radiation coordinates
     1007    INTEGER(iwp)                ::  j     !< Running index, y-dir, radiation coordinates
     1008    INTEGER(iwp)                ::  k     !< Running index, y-dir, radiation coordinates
     1009    INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
     1010
     1011
     1012!-- Calculate biometeorology MRT from local radiation
     1013!-- fluxes calculated by RTM and assign into 2D grid
     1014    tmrt_grid = -999.
     1015    DO  l = 1, nmrtbl
     1016       i = mrtbl(ix,l)
     1017       j = mrtbl(iy,l)
     1018       k = mrtbl(iz,l)
     1019       IF ( k - get_topography_top_index_ji( j, i, 's' ) == bio_cell_level +   &
     1020             1_iwp) THEN
     1021          IF ( mrt_include_sw )  THEN
     1022              tmrt_grid(j,i) = ((human_absorb*mrtinsw(l) +                     &
     1023                                human_emiss*mrtinlw(l))  /                     &
     1024                               (human_emiss*sigma_sb)) ** .25_wp - degc_to_k
     1025          ELSE
     1026              tmrt_grid(j,i) = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp   &
     1027                                 - degc_to_k
     1028          ENDIF
    11311029       ENDIF
    1132        IF ( sw_out <= -998._wp ) THEN
    1133           sw_out = 0.
    1134           no_input = no_input + 1
    1135        ENDIF
    1136        IF ( lw_in <= -998._wp ) THEN
    1137           lw_in = 0.
    1138           no_input = no_input + 1
    1139        ENDIF
    1140        IF ( lw_out <= -998._wp ) THEN
    1141           lw_out = 0.
    1142           no_input = no_input + 1
    1143        ENDIF
    1144 
    1145 !-- Accept two missing radiation flux directions, fail otherwise as error might become too large
    1146        IF ( ta <= -998._wp .OR. no_input >= 3 ) THEN
    1147           tmrt = -999._wp
    1148           RETURN
    1149        ENDIF
    1150     ENDIF
    1151 
    1152     sw_side = sw_in * 0.125_wp  ! distribute half of upper sw_in to the 4 sides
    1153     lw_air = ( sigma_sb * 0.95_wp * ( ta + degc_to_k )**4 )
    1154 
    1155 !-- Compute mean radiation flux density absorbed by rotational symetric human
    1156     nrfd = ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) +         &
    1157        ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) +             &
    1158        ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) +             &
    1159        ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) +             &
    1160        ( weight_v * ( human_absorb * (sw_in * 0.5_wp) + human_emiss * lw_in ) ) +     &
    1161        ( weight_v * ( human_absorb * sw_out + human_emiss * lw_out ) )
    1162 
    1163 !-- Compute mean radiant temperature
    1164     tmrt = ( nrfd / (human_emiss * sigma_sb) )**0.25_wp - degc_to_k
    1165 
    1166  END SUBROUTINE calculate_tmrt_2_directions
    1167 
    1168 !------------------------------------------------------------------------------!
    1169 ! Description:
    1170 ! ------------
    1171 !> Calculate static thermal indices for given meteorological conditions
    1172 !------------------------------------------------------------------------------!
    1173  SUBROUTINE calculate_static_thermal_indices ( ta, vp, ws, pair, tmrt,         &
    1174     pt_static, utci_static, pet_static )
    1175 
    1176     IMPLICIT NONE
    1177 
    1178 !-- Input parameters
    1179     REAL(wp), INTENT ( IN )  ::  ta           !< Air temperature                  (°C)
    1180     REAL(wp), INTENT ( IN )  ::  vp           !< Vapour pressure                  (hPa)
    1181     REAL(wp), INTENT ( IN )  ::  ws           !< Wind speed    (local level)      (m/s)
    1182     REAL(wp), INTENT ( IN )  ::  pair         !< Air pressure                     (hPa)
    1183     REAL(wp), INTENT ( IN )  ::  tmrt         !< Mean radiant temperature         (°C)
    1184 !-- Output parameters
    1185     REAL(wp), INTENT ( OUT ) ::  pt_static    !< Perceived temperature            (°C)
    1186     REAL(wp), INTENT ( OUT ) ::  utci_static  !< Universal thermal climate index  (°C)
    1187     REAL(wp), INTENT ( OUT ) ::  pet_static   !< Physiologically equivalent temp. (°C)
    1188 !-- Temporary field, not used here
    1189     REAL(wp)                 ::  clo          !< Clothing index (no dim.)
    1190 
    1191     clo = -999._wp
    1192 
    1193     IF ( biom_pt ) THEN
    1194 !-- Estimate local perceived temperature
    1195        CALL calculate_pt_static( ta, vp, ws, tmrt, pair, clo, pt_static )
    1196     ENDIF
    1197 
    1198     IF ( biom_utci ) THEN
    1199 !-- Estimate local universal thermal climate index
    1200        CALL calculate_utci_static( ta, vp, ws, tmrt, biom_output_height,       &
    1201           utci_static )
    1202     ENDIF
    1203 
    1204     IF ( biom_pet ) THEN
    1205 !-- Estimate local physiologically equivalent temperature
    1206        CALL calculate_pet_static( ta, vp, ws, tmrt, pair, pet_static )
    1207     ENDIF
    1208 
    1209  END SUBROUTINE calculate_static_thermal_indices
     1030    ENDDO
     1031
     1032END SUBROUTINE bio_calculate_mrt_grid
    12101033
    12111034
     
    12151038!> Calculate static thermal indices for 2D grid point i, j
    12161039!------------------------------------------------------------------------------!
    1217  SUBROUTINE biom_determine_input_at( average_input, i, j, ta, vp, ws, pair,    &
    1218     tmrt )
     1040 SUBROUTINE bio_get_thermal_index_input_ij( average_input, i, j, ta, vp, ws,   &
     1041    pair, tmrt )
    12191042
    12201043    IMPLICIT NONE
     
    12341057!-- Internal variables
    12351058    INTEGER(iwp)                ::  k     !< Running index, z-dir
     1059    INTEGER(iwp)                ::  ir    !< Running index, x-dir, radiation coordinates
     1060    INTEGER(iwp)                ::  jr    !< Running index, y-dir, radiation coordinates
     1061    INTEGER(iwp)                ::  kr    !< Running index, y-dir, radiation coordinates
    12361062    INTEGER(iwp)                ::  k_wind  !< Running index, z-dir, wind speed only
     1063    INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
    12371064
    12381065    REAL(wp)                    ::  vp_sat  !< Saturation vapor pressure     (hPa)
     
    12411068!-- Determine cell level closest to 1.1m above ground
    12421069!   by making use of truncation due to int cast
    1243     k = get_topography_top_index_ji(j, i, 's') + biom_cell_level  !< Vertical cell center closest to 1.1m
     1070    k = get_topography_top_index_ji(j, i, 's') + bio_cell_level  !< Vertical cell center closest to 1.1m
     1071
     1072!
     1073!-- Avoid non-representative horizontal u and v of 0.0 m/s too close to ground
    12441074    k_wind = k
    1245     IF( k_wind < 1_iwp ) THEN  ! Avoid horizontal u and v of 0.0 m/s close to ground
    1246        k_wind = 1_iwp
     1075    IF ( bio_cell_level < 1_iwp ) THEN
     1076       k_wind = k + 1_iwp
    12471077    ENDIF
    12481078
     
    12521082       ta = pt_av(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
    12531083
    1254        vp = 0.034_wp
     1084       vp = -999._wp
    12551085       IF ( humidity .AND. ALLOCATED( q_av ) ) THEN
    12561086          vp = q_av(k, j, i)
     
    12641094       ta = pt(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
    12651095
    1266        vp = q(k, j, i)
     1096       vp = -999._wp
     1097       IF ( humidity ) THEN
     1098          vp = q(k, j, i)
     1099       ENDIF
    12671100
    12681101       ws = ( 0.5_wp * ABS( u(k_wind, j, i) + u(k_wind, j, i+1) )  +           &
     
    12781111!-- The magnus formula is limited to temperatures up to 333.15 K to
    12791112!   avoid negative values of vp_sat
    1280     vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) )
    1281     vp  = vp * pair / ( vp + 0.622_wp )
    1282     IF ( vp > vp_sat ) vp = vp_sat
    1283 
    1284     tmrt = ta
     1113    IF ( vp > -998._wp ) THEN
     1114       vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) )
     1115       vp  = vp * pair / ( vp + 0.622_wp )
     1116       IF ( vp > vp_sat ) vp = vp_sat
     1117    ENDIF
     1118
     1119!-- local mtr value at [i,j]
     1120    tmrt = -999.  !< this can be a valid result (e.g. for inside some ostacle)
    12851121    IF ( radiation ) THEN
    1286        IF ( mrt_from_rtm ) THEN
    1287           tmrt = tmrt_grid(j, i)
    1288        ELSE
    1289           CALL calculate_tmrt_2_directions (rad_sw_in(0, j, i),                &
    1290              rad_sw_out(0, j, i), rad_lw_in(0, j, i), rad_lw_out(0, j, i), ta, &
    1291              tmrt )
    1292        ENDIF
    1293     ENDIF
    1294 
    1295  END SUBROUTINE biom_determine_input_at
     1122!--    Use MRT from RTM precalculated in tmrt_grid
     1123       tmrt = tmrt_grid(j,i)
     1124    ENDIF
     1125
     1126 END SUBROUTINE bio_get_thermal_index_input_ij
    12961127
    12971128
     
    13021133!> time_integration.f90: 1065ff
    13031134!------------------------------------------------------------------------------!
    1304  SUBROUTINE biom_calculate_static_grid ( average_input )
     1135 SUBROUTINE bio_calculate_thermal_index_maps ( av )
    13051136
    13061137    IMPLICIT NONE
    13071138
    13081139!-- Input attributes
    1309     LOGICAL, INTENT ( IN ) ::  average_input  !< Calculate based on averaged input? conditions?
     1140    LOGICAL, INTENT ( IN ) ::  av  !< Calculate based on averaged input conditions?
    13101141
    13111142!-- Internal variables
    1312     INTEGER(iwp) ::  i, j, k, l   !< Running index
    1313 
     1143    INTEGER(iwp) ::  i, j      !< Running index
     1144
     1145    REAL(wp) ::  clo           !< Clothing index                (no dimension)
    13141146    REAL(wp) ::  ta            !< Air temperature                  (°C)
    13151147    REAL(wp) ::  vp            !< Vapour pressure                  (hPa)
    13161148    REAL(wp) ::  ws            !< Wind speed    (local level)      (m/s)
    13171149    REAL(wp) ::  pair          !< Air pressure                     (hPa)
    1318     REAL(wp) ::  tmrt_tmp      !< Mean radiant temperature
    1319     REAL(wp) ::  pt_tmp        !< Perceived temperature
    1320     REAL(wp) ::  utci_tmp      !< Universal thermal climate index
    1321     REAL(wp) ::  pet_tmp       !< Physiologically equivalent temperature
    1322 
    1323 
    1324     CALL biom_init_arrays
    1325 
    1326     IF ( mrt_from_rtm ) THEN
    1327        tmrt_grid = -999._wp
    1328        DO  l = 1, nmrtbl
    1329            i = mrtbl(ix,l)
    1330            j = mrtbl(iy,l)
    1331            k = mrtbl(iz,l)
    1332            IF ( k - get_topography_top_index_ji( j, i, 's' ) == 1  ) THEN
    1333               IF ( mrt_include_sw )  THEN
    1334                   tmrt_tmp = ((human_absorb*mrtinsw(l) + human_emiss*mrtinlw(l)) &
    1335                              / (human_emiss*sigma_sb)) ** .25_wp
    1336               ELSE
    1337                   tmrt_tmp = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp
    1338               ENDIF
    1339               tmrt_grid(j, i) = tmrt_tmp - degc_to_k
    1340            ENDIF
    1341        ENDDO
    1342     ENDIF
     1150    REAL(wp) ::  perct_tmp     !< Perceived temperature            (°C)
     1151    REAL(wp) ::  utci_tmp      !< Universal thermal climate index  (°C)
     1152    REAL(wp) ::  pet_tmp       !< Physiologically equivalent temperature  (°C)
     1153    REAL(wp) ::  tmrt_tmp      !< Mean radiant temperature         (°C)
     1154
     1155    CALL bio_init_arrays
     1156
     1157!-- fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...)
     1158    CALL bio_calculate_mrt_grid ( av )
     1159
    13431160
    13441161    DO i = nxl, nxr
    13451162       DO j = nys, nyn
    13461163!--       Determine local input conditions
    1347           CALL biom_determine_input_at ( average_input, i, j, ta, vp, ws,      &
    1348                pair, tmrt_tmp )
    1349           tmrt_grid(j, i) = tmrt_tmp
    1350 
    1351 !--       Only proceed if tmrt is available
    1352           IF ( tmrt_tmp <= -998._wp ) THEN
    1353              pt_tmp   = -999._wp
     1164          CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp,              &
     1165               ws, pair, tmrt_tmp )
     1166!           tmrt_grid(j, i) = tmrt_tmp  !< already set by bio_calculate_mrt_grid!
     1167
     1168!--       Only proceed if input is available
     1169          IF ( tmrt_tmp <= -998._wp .OR. vp <= -998._wp ) THEN
     1170             pet_tmp = -999._wp         !< set fail value, e.g. valid for within
     1171             perct_tmp  = -999._wp      !< some obstacle
    13541172             utci_tmp = -999._wp
    1355              pet_tmp  = -999._wp
    1356              CYCLE
     1173          ELSE
     1174!--          Calculate static thermal indices based on local tmrt
     1175             clo = -999._wp
     1176             
     1177             IF ( bio_perct ) THEN
     1178!--          Estimate local perceived temperature
     1179                CALL calculate_perct_static( ta, vp, ws, tmrt_tmp, pair, clo,  &
     1180                   perct_tmp )
     1181             ENDIF
     1182             
     1183             IF ( bio_utci ) THEN
     1184!--          Estimate local universal thermal climate index
     1185                CALL calculate_utci_static( ta, vp, ws, tmrt_tmp,              &
     1186                   bio_output_height, utci_tmp )
     1187             ENDIF
     1188             
     1189             IF ( bio_pet ) THEN
     1190!--          Estimate local physiologically equivalent temperature
     1191                CALL calculate_pet_static( ta, vp, ws, tmrt_tmp, pair, pet_tmp )
     1192             ENDIF
    13571193          END IF
    13581194
    1359 !--       Calculate static thermal indices based on local tmrt
    1360           CALL calculate_static_thermal_indices ( ta, vp, ws,                  &
    1361              pair, tmrt_tmp, pt_tmp, utci_tmp, pet_tmp )
    1362 
    1363           IF ( average_input ) THEN
     1195
     1196          IF ( av ) THEN
    13641197!--          Write results for selected averaged indices
    1365              IF ( biom_pt_av )  THEN
    1366                 pt_av_grid(j, i)   = pt_tmp
     1198             IF ( bio_perct_av )  THEN
     1199                perct_av(j, i) = perct_tmp
    13671200             END IF
    1368              IF ( biom_utci_av ) THEN
     1201             IF ( bio_utci_av ) THEN
    13691202                utci_av_grid(j, i) = utci_tmp
    13701203             END IF
    1371              IF ( biom_pet_av ) THEN
     1204             IF ( bio_pet_av ) THEN
    13721205                pet_av_grid(j, i)  = pet_tmp
    13731206             END IF
    13741207          ELSE
    13751208!--          Write result for selected indices
    1376              IF ( biom_pt )  THEN
    1377                 pt_grid(j, i)   = pt_tmp
     1209             IF ( bio_perct )  THEN
     1210                perct(j, i) = perct_tmp
    13781211             END IF
    1379              IF ( biom_utci ) THEN
     1212             IF ( bio_utci ) THEN
    13801213                utci_grid(j, i) = utci_tmp
    13811214             END IF
    1382              IF ( biom_pet ) THEN
     1215             IF ( bio_pet ) THEN
    13831216                pet_grid(j, i)  = pet_tmp
    13841217             END IF
     
    13881221    END DO
    13891222
    1390  END SUBROUTINE biom_calculate_static_grid
     1223 END SUBROUTINE bio_calculate_thermal_index_maps
    13911224
    13921225!------------------------------------------------------------------------------!
     
    13951228!> Calculate dynamic thermal indices (currently only iPT, but expandable)
    13961229!------------------------------------------------------------------------------!
    1397  SUBROUTINE biom_calc_ipt( ta, vp, ws, pair, tmrt, dt, energy_storage,         &
     1230 SUBROUTINE bio_calc_ipt( ta, vp, ws, pair, tmrt, dt, energy_storage,          &
    13981231    t_clo, clo, actlev, age, weight, height, work, sex, ipt )
    13991232
     
    14361269    ENDIF
    14371270
    1438  END SUBROUTINE biom_calc_ipt
     1271 END SUBROUTINE bio_calc_ipt
     1272
     1273
     1274!------------------------------------------------------------------------------!
     1275! Description:
     1276! ------------
     1277!> SUBROUTINE for calculating UTCI Temperature (UTCI)
     1278!> computed by a 6th order approximation
     1279!
     1280!> UTCI regression equation after
     1281!> Bröde P, Fiala D, Blazejczyk K, Holmér I, Jendritzky G, Kampmann B, Tinz B,
     1282!> Havenith G (2012) Deriving the operational procedure for the Universal Thermal
     1283!> Climate Index (UTCI). International Journal of Biometeorology 56 (3):481-494.
     1284!> doi:10.1007/s00484-011-0454-1
     1285!
     1286!> original source available at:
     1287!> www.utci.org
     1288!------------------------------------------------------------------------------!
     1289 SUBROUTINE calculate_utci_static( ta_in, vp, ws_hag, tmrt, hag, utci )
     1290
     1291    IMPLICIT NONE
     1292
     1293!-- Type of input of the argument list
     1294    REAL(WP), INTENT ( IN )  ::  ta_in  !< Local air temperature (°C)
     1295    REAL(WP), INTENT ( IN )  ::  vp     !< Loacl vapour pressure (hPa)
     1296    REAL(WP), INTENT ( IN )  ::  ws_hag !< Incident wind speed (m/s)
     1297    REAL(WP), INTENT ( IN )  ::  tmrt   !< Local mean radiant temperature (°C)
     1298    REAL(WP), INTENT ( IN )  ::  hag    !< Height of wind speed input (m)
     1299!-- Type of output of the argument list
     1300    REAL(wp), INTENT ( OUT ) ::  utci   !< Universal Thermal Climate Index (°C)
     1301
     1302!-- Make sure precission is sufficient for regression equation
     1303    REAL(WP) ::  ta         !< internal air temperature (°C)
     1304    REAL(WP) ::  pa         !< air pressure in kPa      (kPa)
     1305    REAL(WP) ::  d_tmrt     !< delta-tmrt               (°C)
     1306    REAL(WP) ::  va         !< wind speed at 10 m above ground level    (m/s)
     1307    REAL(WP) ::  offset     !< utci deviation by ta cond. exceeded      (°C)
     1308
     1309!-- Initialize
     1310    offset = 0._wp
     1311    ta = ta_in
     1312    d_tmrt = tmrt - ta_in
     1313
     1314!-- Use vapour pressure in kpa
     1315    pa = vp / 10.0_wp
     1316
     1317!-- Wind altitude correction from hag to 10m after Broede et al. (2012), eq.3
     1318!   z(0) is set to 0.01 according to UTCI profile definition
     1319    va = ws_hag *  log ( 10.0_wp / 0.01_wp ) / log ( hag / 0.01_wp )
     1320
     1321!-- Check if input values in range after Broede et al. (2012)
     1322    IF ( ( d_tmrt > 70._wp ) .OR. ( d_tmrt < -30._wp ) .OR.                    &
     1323       ( vp >= 50._wp ) ) THEN
     1324       utci = -999._wp
     1325       RETURN
     1326    ENDIF
     1327
     1328!-- Apply eq. 2 in Broede et al. (2012) for ta out of bounds
     1329    IF ( ta > 50._wp ) THEN
     1330       offset = ta - 50._wp
     1331       ta = 50._wp
     1332    ENDIF
     1333    IF ( ta < -50._wp ) THEN
     1334       offset = ta + 50._wp
     1335       ta = -50._wp
     1336    ENDIF
     1337
     1338!-- For routine application. For wind speeds and relative
     1339!   humidity values below 0.5 m/s or 5%, respectively, the
     1340!   user is advised to use the lower bounds for the calculations.
     1341    IF ( va < 0.5_wp ) va = 0.5_wp
     1342    IF ( va > 17._wp ) va = 17._wp
     1343
     1344!-- Calculate 6th order polynomial as approximation
     1345    utci = ta +                                                                &
     1346       ( 6.07562052e-01 )   +                                                  &
     1347       ( -2.27712343e-02 ) * ta +                                              &
     1348       ( 8.06470249e-04 )  * ta * ta +                                         &
     1349       ( -1.54271372e-04 ) * ta * ta * ta +                                    &
     1350       ( -3.24651735e-06 ) * ta * ta * ta * ta +                               &
     1351       ( 7.32602852e-08 )  * ta * ta * ta * ta * ta +                          &
     1352       ( 1.35959073e-09 )  * ta * ta * ta * ta * ta * ta +                     &
     1353       ( -2.25836520e+00 ) * va +                                              &
     1354       ( 8.80326035e-02 )  * ta * va +                                         &
     1355       ( 2.16844454e-03 )  * ta * ta * va +                                    &
     1356       ( -1.53347087e-05 ) * ta * ta * ta * va +                               &
     1357       ( -5.72983704e-07 ) * ta * ta * ta * ta * va +                          &
     1358       ( -2.55090145e-09 ) * ta * ta * ta * ta * ta * va +                     &
     1359       ( -7.51269505e-01 ) * va * va +                                         &
     1360       ( -4.08350271e-03 ) * ta * va * va +                                    &
     1361       ( -5.21670675e-05 ) * ta * ta * va * va +                               &
     1362       ( 1.94544667e-06 )  * ta * ta * ta * va * va +                          &
     1363       ( 1.14099531e-08 )  * ta * ta * ta * ta * va * va +                     &
     1364       ( 1.58137256e-01 )  * va * va * va +                                    &
     1365       ( -6.57263143e-05 ) * ta * va * va * va +                               &
     1366       ( 2.22697524e-07 )  * ta * ta * va * va * va +                          &
     1367       ( -4.16117031e-08 ) * ta * ta * ta * va * va * va +                     &
     1368       ( -1.27762753e-02 ) * va * va * va * va +                               &
     1369       ( 9.66891875e-06 )  * ta * va * va * va * va +                          &
     1370       ( 2.52785852e-09 )  * ta * ta * va * va * va * va +                     &
     1371       ( 4.56306672e-04 )  * va * va * va * va * va +                          &
     1372       ( -1.74202546e-07 ) * ta * va * va * va * va * va +                     &
     1373       ( -5.91491269e-06 ) * va * va * va * va * va * va +                     &
     1374       ( 3.98374029e-01 )  * d_tmrt +                                          &
     1375       ( 1.83945314e-04 )  * ta * d_tmrt +                                     &
     1376       ( -1.73754510e-04 ) * ta * ta * d_tmrt +                                &
     1377       ( -7.60781159e-07 ) * ta * ta * ta * d_tmrt +                           &
     1378       ( 3.77830287e-08 )  * ta * ta * ta * ta * d_tmrt +                      &
     1379       ( 5.43079673e-10 )  * ta * ta * ta * ta * ta * d_tmrt +                 &
     1380       ( -2.00518269e-02 ) * va * d_tmrt +                                     &
     1381       ( 8.92859837e-04 )  * ta * va * d_tmrt +                                &
     1382       ( 3.45433048e-06 )  * ta * ta * va * d_tmrt +                           &
     1383       ( -3.77925774e-07 ) * ta * ta * ta * va * d_tmrt +                      &
     1384       ( -1.69699377e-09 ) * ta * ta * ta * ta * va * d_tmrt +                 &
     1385       ( 1.69992415e-04 )  * va * va * d_tmrt +                                &
     1386       ( -4.99204314e-05 ) * ta * va * va * d_tmrt +                           &
     1387       ( 2.47417178e-07 )  * ta * ta * va * va * d_tmrt +                      &
     1388       ( 1.07596466e-08 )  * ta * ta * ta * va * va * d_tmrt +                 &
     1389       ( 8.49242932e-05 )  * va * va * va * d_tmrt +                           &
     1390       ( 1.35191328e-06 )  * ta * va * va * va * d_tmrt +                      &
     1391       ( -6.21531254e-09 ) * ta * ta * va * va * va * d_tmrt +                 &
     1392       ( -4.99410301e-06 ) * va * va * va * va * d_tmrt +                      &
     1393       ( -1.89489258e-08 ) * ta * va * va * va * va * d_tmrt +                 &
     1394       ( 8.15300114e-08 )  * va * va * va * va * va * d_tmrt +                 &
     1395       ( 7.55043090e-04 )  * d_tmrt * d_tmrt +                                 &
     1396       ( -5.65095215e-05 ) * ta * d_tmrt * d_tmrt +                            &
     1397       ( -4.52166564e-07 ) * ta * ta * d_tmrt * d_tmrt +                       &
     1398       ( 2.46688878e-08 )  * ta * ta * ta * d_tmrt * d_tmrt +                  &
     1399       ( 2.42674348e-10 )  * ta * ta * ta * ta * d_tmrt * d_tmrt +             &
     1400       ( 1.54547250e-04 )  * va * d_tmrt * d_tmrt +                            &
     1401       ( 5.24110970e-06 )  * ta * va * d_tmrt * d_tmrt +                       &
     1402       ( -8.75874982e-08 ) * ta * ta * va * d_tmrt * d_tmrt +                  &
     1403       ( -1.50743064e-09 ) * ta * ta * ta * va * d_tmrt * d_tmrt +             &
     1404       ( -1.56236307e-05 ) * va * va * d_tmrt * d_tmrt +                       &
     1405       ( -1.33895614e-07 ) * ta * va * va * d_tmrt * d_tmrt +                  &
     1406       ( 2.49709824e-09 )  * ta * ta * va * va * d_tmrt * d_tmrt +             &
     1407       ( 6.51711721e-07 )  * va * va * va * d_tmrt * d_tmrt +                  &
     1408       ( 1.94960053e-09 )  * ta * va * va * va * d_tmrt * d_tmrt +             &
     1409       ( -1.00361113e-08 ) * va * va * va * va * d_tmrt * d_tmrt +             &
     1410       ( -1.21206673e-05 ) * d_tmrt * d_tmrt * d_tmrt +                        &
     1411       ( -2.18203660e-07 ) * ta * d_tmrt * d_tmrt * d_tmrt +                   &
     1412       ( 7.51269482e-09 )  * ta * ta * d_tmrt * d_tmrt * d_tmrt +              &
     1413       ( 9.79063848e-11 )  * ta * ta * ta * d_tmrt * d_tmrt * d_tmrt +         &
     1414       ( 1.25006734e-06 )  * va * d_tmrt * d_tmrt * d_tmrt +                   &
     1415       ( -1.81584736e-09 ) * ta * va * d_tmrt * d_tmrt * d_tmrt +              &
     1416       ( -3.52197671e-10 ) * ta * ta * va * d_tmrt * d_tmrt * d_tmrt +         &
     1417       ( -3.36514630e-08 ) * va * va * d_tmrt * d_tmrt * d_tmrt +              &
     1418       ( 1.35908359e-10 )  * ta * va * va * d_tmrt * d_tmrt * d_tmrt +         &
     1419       ( 4.17032620e-10 )  * va * va * va * d_tmrt * d_tmrt * d_tmrt +         &
     1420       ( -1.30369025e-09 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt +               &
     1421       ( 4.13908461e-10 )  * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt +          &
     1422       ( 9.22652254e-12 )  * ta * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt +     &
     1423       ( -5.08220384e-09 ) * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt +          &
     1424       ( -2.24730961e-11 ) * ta * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt +     &
     1425       ( 1.17139133e-10 )  * va * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt +     &
     1426       ( 6.62154879e-10 )  * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt +      &
     1427       ( 4.03863260e-13 )  * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt + &
     1428       ( 1.95087203e-12 )  * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt + &
     1429       ( -4.73602469e-12 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt *      &
     1430       d_tmrt +                                                                &
     1431       ( 5.12733497e+00 )  * pa +                                              &
     1432       ( -3.12788561e-01 ) * ta * pa +                                         &
     1433       ( -1.96701861e-02 ) * ta * ta * pa +                                    &
     1434       ( 9.99690870e-04 )  * ta * ta * ta * pa +                               &
     1435       ( 9.51738512e-06 )  * ta * ta * ta * ta * pa +                          &
     1436       ( -4.66426341e-07 ) * ta * ta * ta * ta * ta * pa +                     &
     1437       ( 5.48050612e-01 )  * va * pa +                                         &
     1438       ( -3.30552823e-03 ) * ta * va * pa +                                    &
     1439       ( -1.64119440e-03 ) * ta * ta * va * pa +                               &
     1440       ( -5.16670694e-06 ) * ta * ta * ta * va * pa +                          &
     1441       ( 9.52692432e-07 )  * ta * ta * ta * ta * va * pa +                     &
     1442       ( -4.29223622e-02 ) * va * va * pa +                                    &
     1443       ( 5.00845667e-03 )  * ta * va * va * pa +                               &
     1444       ( 1.00601257e-06 )  * ta * ta * va * va * pa +                          &
     1445       ( -1.81748644e-06 ) * ta * ta * ta * va * va * pa +                     &
     1446       ( -1.25813502e-03 ) * va * va * va * pa +                               &
     1447       ( -1.79330391e-04 ) * ta * va * va * va * pa +                          &
     1448       ( 2.34994441e-06 )  * ta * ta * va * va * va * pa +                     &
     1449       ( 1.29735808e-04 )  * va * va * va * va * pa +                          &
     1450       ( 1.29064870e-06 )  * ta * va * va * va * va * pa +                     &
     1451       ( -2.28558686e-06 ) * va * va * va * va * va * pa +                     &
     1452       ( -3.69476348e-02 ) * d_tmrt * pa +                                     &
     1453       ( 1.62325322e-03 )  * ta * d_tmrt * pa +                                &
     1454       ( -3.14279680e-05 ) * ta * ta * d_tmrt * pa +                           &
     1455       ( 2.59835559e-06 )  * ta * ta * ta * d_tmrt * pa +                      &
     1456       ( -4.77136523e-08 ) * ta * ta * ta * ta * d_tmrt * pa +                 &
     1457       ( 8.64203390e-03 )  * va * d_tmrt * pa +                                &
     1458       ( -6.87405181e-04 ) * ta * va * d_tmrt * pa +                           &
     1459       ( -9.13863872e-06 ) * ta * ta * va * d_tmrt * pa +                      &
     1460       ( 5.15916806e-07 )  * ta * ta * ta * va * d_tmrt * pa +                 &
     1461       ( -3.59217476e-05 ) * va * va * d_tmrt * pa +                           &
     1462       ( 3.28696511e-05 )  * ta * va * va * d_tmrt * pa +                      &
     1463       ( -7.10542454e-07 ) * ta * ta * va * va * d_tmrt * pa +                 &
     1464       ( -1.24382300e-05 ) * va * va * va * d_tmrt * pa +                      &
     1465       ( -7.38584400e-09 ) * ta * va * va * va * d_tmrt * pa +                 &
     1466       ( 2.20609296e-07 )  * va * va * va * va * d_tmrt * pa +                 &
     1467       ( -7.32469180e-04 ) * d_tmrt * d_tmrt * pa +                            &
     1468       ( -1.87381964e-05 ) * ta * d_tmrt * d_tmrt * pa +                       &
     1469       ( 4.80925239e-06 )  * ta * ta * d_tmrt * d_tmrt * pa +                  &
     1470       ( -8.75492040e-08 ) * ta * ta * ta * d_tmrt * d_tmrt * pa +             &
     1471       ( 2.77862930e-05 )  * va * d_tmrt * d_tmrt * pa +                       &
     1472       ( -5.06004592e-06 ) * ta * va * d_tmrt * d_tmrt * pa +                  &
     1473       ( 1.14325367e-07 )  * ta * ta * va * d_tmrt * d_tmrt * pa +             &
     1474       ( 2.53016723e-06 )  * va * va * d_tmrt * d_tmrt * pa +                  &
     1475       ( -1.72857035e-08 ) * ta * va * va * d_tmrt * d_tmrt * pa +             &
     1476       ( -3.95079398e-08 ) * va * va * va * d_tmrt * d_tmrt * pa +             &
     1477       ( -3.59413173e-07 ) * d_tmrt * d_tmrt * d_tmrt * pa +                   &
     1478       ( 7.04388046e-07 )  * ta * d_tmrt * d_tmrt * d_tmrt * pa +              &
     1479       ( -1.89309167e-08 ) * ta * ta * d_tmrt * d_tmrt * d_tmrt * pa +         &
     1480       ( -4.79768731e-07 ) * va * d_tmrt * d_tmrt * d_tmrt * pa +              &
     1481       ( 7.96079978e-09 )  * ta * va * d_tmrt * d_tmrt * d_tmrt * pa +         &
     1482       ( 1.62897058e-09 )  * va * va * d_tmrt * d_tmrt * d_tmrt * pa +         &
     1483       ( 3.94367674e-08 )  * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa +          &
     1484       ( -1.18566247e-09 ) * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa +     &
     1485       ( 3.34678041e-10 )  * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa +     &
     1486       ( -1.15606447e-10 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa + &
     1487       ( -2.80626406e+00 ) * pa * pa +                                         &
     1488       ( 5.48712484e-01 )  * ta * pa * pa +                                    &
     1489       ( -3.99428410e-03 ) * ta * ta * pa * pa +                               &
     1490       ( -9.54009191e-04 ) * ta * ta * ta * pa * pa +                          &
     1491       ( 1.93090978e-05 )  * ta * ta * ta * ta * pa * pa +                     &
     1492       ( -3.08806365e-01 ) * va * pa * pa +                                    &
     1493       ( 1.16952364e-02 )  * ta * va * pa * pa +                               &
     1494       ( 4.95271903e-04 )  * ta * ta * va * pa * pa +                          &
     1495       ( -1.90710882e-05 ) * ta * ta * ta * va * pa * pa +                     &
     1496       ( 2.10787756e-03 )  * va * va * pa * pa +                               &
     1497       ( -6.98445738e-04 ) * ta * va * va * pa * pa +                          &
     1498       ( 2.30109073e-05 )  * ta * ta * va * va * pa * pa +                     &
     1499       ( 4.17856590e-04 )  * va * va * va * pa * pa +                          &
     1500       ( -1.27043871e-05 ) * ta * va * va * va * pa * pa +                     &
     1501       ( -3.04620472e-06 ) * va * va * va * va * pa * pa +                     &
     1502       ( 5.14507424e-02 )  * d_tmrt * pa * pa +                                &
     1503       ( -4.32510997e-03 ) * ta * d_tmrt * pa * pa +                           &
     1504       ( 8.99281156e-05 )  * ta * ta * d_tmrt * pa * pa +                      &
     1505       ( -7.14663943e-07 ) * ta * ta * ta * d_tmrt * pa * pa +                 &
     1506       ( -2.66016305e-04 ) * va * d_tmrt * pa * pa +                           &
     1507       ( 2.63789586e-04 )  * ta * va * d_tmrt * pa * pa +                      &
     1508       ( -7.01199003e-06 ) * ta * ta * va * d_tmrt * pa * pa +                 &
     1509       ( -1.06823306e-04 ) * va * va * d_tmrt * pa * pa +                      &
     1510       ( 3.61341136e-06 )  * ta * va * va * d_tmrt * pa * pa +                 &
     1511       ( 2.29748967e-07 )  * va * va * va * d_tmrt * pa * pa +                 &
     1512       ( 3.04788893e-04 )  * d_tmrt * d_tmrt * pa * pa +                       &
     1513       ( -6.42070836e-05 ) * ta * d_tmrt * d_tmrt * pa * pa +                  &
     1514       ( 1.16257971e-06 )  * ta * ta * d_tmrt * d_tmrt * pa * pa +             &
     1515       ( 7.68023384e-06 )  * va * d_tmrt * d_tmrt * pa * pa +                  &
     1516       ( -5.47446896e-07 ) * ta * va * d_tmrt * d_tmrt * pa * pa +             &
     1517       ( -3.59937910e-08 ) * va * va * d_tmrt * d_tmrt * pa * pa +             &
     1518       ( -4.36497725e-06 ) * d_tmrt * d_tmrt * d_tmrt * pa * pa +              &
     1519       ( 1.68737969e-07 )  * ta * d_tmrt * d_tmrt * d_tmrt * pa * pa +         &
     1520       ( 2.67489271e-08 )  * va * d_tmrt * d_tmrt * d_tmrt * pa * pa +         &
     1521       ( 3.23926897e-09 )  * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa * pa +     &
     1522       ( -3.53874123e-02 ) * pa * pa * pa +                                    &
     1523       ( -2.21201190e-01 ) * ta * pa * pa * pa +                               &
     1524       ( 1.55126038e-02 )  * ta * ta * pa * pa * pa +                          &
     1525       ( -2.63917279e-04 ) * ta * ta * ta * pa * pa * pa +                     &
     1526       ( 4.53433455e-02 )  * va * pa * pa * pa +                               &
     1527       ( -4.32943862e-03 ) * ta * va * pa * pa * pa +                          &
     1528       ( 1.45389826e-04 )  * ta * ta * va * pa * pa * pa +                     &
     1529       ( 2.17508610e-04 )  * va * va * pa * pa * pa +                          &
     1530       ( -6.66724702e-05 ) * ta * va * va * pa * pa * pa +                     &
     1531       ( 3.33217140e-05 )  * va * va * va * pa * pa * pa +                     &
     1532       ( -2.26921615e-03 ) * d_tmrt * pa * pa * pa +                           &
     1533       ( 3.80261982e-04 )  * ta * d_tmrt * pa * pa * pa +                      &
     1534       ( -5.45314314e-09 ) * ta * ta * d_tmrt * pa * pa * pa +                 &
     1535       ( -7.96355448e-04 ) * va * d_tmrt * pa * pa * pa +                      &
     1536       ( 2.53458034e-05 )  * ta * va * d_tmrt * pa * pa * pa +                 &
     1537       ( -6.31223658e-06 ) * va * va * d_tmrt * pa * pa * pa +                 &
     1538       ( 3.02122035e-04 )  * d_tmrt * d_tmrt * pa * pa * pa +                  &
     1539       ( -4.77403547e-06 ) * ta * d_tmrt * d_tmrt * pa * pa * pa +             &
     1540       ( 1.73825715e-06 )  * va * d_tmrt * d_tmrt * pa * pa * pa +             &
     1541       ( -4.09087898e-07 ) * d_tmrt * d_tmrt * d_tmrt * pa * pa * pa +         &
     1542       ( 6.14155345e-01 )  * pa * pa * pa * pa +                               &
     1543       ( -6.16755931e-02 ) * ta * pa * pa * pa * pa +                          &
     1544       ( 1.33374846e-03 )  * ta * ta * pa * pa * pa * pa +                     &
     1545       ( 3.55375387e-03 )  * va * pa * pa * pa * pa +                          &
     1546       ( -5.13027851e-04 ) * ta * va * pa * pa * pa * pa +                     &
     1547       ( 1.02449757e-04 )  * va * va * pa * pa * pa * pa +                     &
     1548       ( -1.48526421e-03 ) * d_tmrt * pa * pa * pa * pa +                      &
     1549       ( -4.11469183e-05 ) * ta * d_tmrt * pa * pa * pa * pa +                 &
     1550       ( -6.80434415e-06 ) * va * d_tmrt * pa * pa * pa * pa +                 &
     1551       ( -9.77675906e-06 ) * d_tmrt * d_tmrt * pa * pa * pa * pa +             &
     1552       ( 8.82773108e-02 )  * pa * pa * pa * pa * pa +                          &
     1553       ( -3.01859306e-03 ) * ta * pa * pa * pa * pa * pa +                     &
     1554       ( 1.04452989e-03 )  * va * pa * pa * pa * pa * pa +                     &
     1555       ( 2.47090539e-04 )  * d_tmrt * pa * pa * pa * pa * pa +                 &
     1556       ( 1.48348065e-03 )  * pa * pa * pa * pa * pa * pa
     1557
     1558!-- Consider offset in result
     1559    utci = utci + offset
     1560
     1561 END SUBROUTINE calculate_utci_static
     1562
     1563
     1564
     1565
     1566!------------------------------------------------------------------------------!
     1567! Description:
     1568! ------------
     1569!> calculate_perct_static: Estimation of perceived temperature (PT, degC)
     1570!> Value of perct is the Perceived Temperature, degree centigrade
     1571!------------------------------------------------------------------------------!
     1572 SUBROUTINE calculate_perct_static( ta, vp, ws, tmrt, pair, clo, perct )
     1573
     1574    IMPLICIT NONE
     1575
     1576!-- Type of input of the argument list
     1577    REAL(wp), INTENT ( IN )  :: ta   !< Local air temperature (degC)
     1578    REAL(wp), INTENT ( IN )  :: vp   !< Local vapour pressure (hPa)
     1579    REAL(wp), INTENT ( IN )  :: tmrt !< Local mean radiant temperature (degC)
     1580    REAL(wp), INTENT ( IN )  :: ws   !< Local wind velocitry (m/s)
     1581    REAL(wp), INTENT ( IN )  :: pair !< Local barometric air pressure (hPa)
     1582
     1583!-- Type of output of the argument list
     1584    REAL(wp), INTENT ( OUT ) :: perct   !< Perceived temperature (degC)
     1585    REAL(wp), INTENT ( OUT ) :: clo     !< Clothing index (dimensionless)
     1586
     1587!-- Parameters for standard "Klima-Michel"
     1588    REAL(wp), PARAMETER :: eta = 0._wp  !< Mechanical work efficiency for walking on flat ground (compare to Fanger (1972) pp 24f)
     1589    REAL(wp), PARAMETER :: actlev = 134.6862_wp !< Workload by activity per standardized surface (A_Du)
     1590
     1591!-- Type of program variables
     1592    REAL(wp), PARAMETER :: eps = 0.0005  !< Accuracy in clothing insulation (clo) for evaluation the root of Fanger's PMV (pmva=0)
     1593    REAL(wp) ::  sclo           !< summer clothing insulation
     1594    REAL(wp) ::  wclo           !< winter clothing insulation
     1595    REAL(wp) ::  d_pmv          !< PMV deviation (dimensionless --> PMV)
     1596    REAL(wp) ::  svp_ta         !< saturation vapor pressure    (hPa)
     1597    REAL(wp) ::  sult_lim       !< threshold for sultrieness    (hPa)
     1598    REAL(wp) ::  dgtcm          !< Mean deviation dependent on perct
     1599    REAL(wp) ::  dgtcstd        !< Mean deviation plus its standard deviation
     1600    REAL(wp) ::  clon           !< clo for neutral conditions   (clo)
     1601    REAL(wp) ::  ireq_minimal   !< Minimal required clothing insulation (clo)
     1602!     REAL(wp) ::  clo_fanger     !< clo for fanger subroutine, unused
     1603    REAL(wp) ::  pmv_w          !< Fangers predicted mean vote for winter clothing
     1604    REAL(wp) ::  pmv_s          !< Fangers predicted mean vote for summer clothing
     1605    REAL(wp) ::  pmva           !< adjusted predicted mean vote
     1606    REAL(wp) ::  ptc            !< perceived temp. for cold conditions (°C)
     1607    REAL(wp) ::  d_std          !< factor to threshold for sultriness
     1608    REAL(wp) ::  pmvs           !< pred. mean vote considering sultrieness
     1609    REAL(wp) ::  top            !< Gagge's operative temperatures (°C)
     1610
     1611    INTEGER(iwp) :: ncount      !< running index
     1612    INTEGER(iwp) :: nerr_cold   !< error number (cold conditions)
     1613    INTEGER(iwp) :: nerr        !< error number
     1614
     1615    LOGICAL :: sultrieness
     1616
     1617!-- Initialise
     1618    perct = 9999.0_wp
     1619
     1620    nerr     = 0_iwp
     1621    ncount   = 0_iwp
     1622    sultrieness  = .FALSE.
     1623!-- Tresholds: clothing insulation (account for model inaccuracies)
     1624!   summer clothing
     1625    sclo     = 0.44453_wp
     1626!   winter clothing
     1627    wclo     = 1.76267_wp
     1628
     1629!-- decision: firstly calculate for winter or summer clothing
     1630    IF ( ta <= 10._wp ) THEN
     1631
     1632!--    First guess: winter clothing insulation: cold stress
     1633       clo = wclo
     1634       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, top )
     1635       pmv_w = pmva
     1636
     1637       IF ( pmva > 0._wp ) THEN
     1638
     1639!--       Case summer clothing insulation: heat load ?
     1640          clo = sclo
     1641          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
     1642             top )
     1643          pmv_s = pmva
     1644          IF ( pmva <= 0._wp ) THEN
     1645!--          Case: comfort achievable by varying clothing insulation
     1646!--          Between winter and summer set values
     1647             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
     1648                pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
     1649             IF ( ncount < 0_iwp ) THEN
     1650                nerr = -1_iwp
     1651                RETURN
     1652             ENDIF
     1653          ELSE IF ( pmva > 0.06_wp ) THEN
     1654             clo = 0.5_wp
     1655             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
     1656                           pmva, top )
     1657          ENDIF
     1658       ELSE IF ( pmva < -0.11_wp ) THEN
     1659          clo = 1.75_wp
     1660          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
     1661             top )
     1662       ENDIF
     1663    ELSE
     1664
     1665!--    First guess: summer clothing insulation: heat load
     1666       clo = sclo
     1667       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, top )
     1668       pmv_s = pmva
     1669
     1670       IF ( pmva < 0._wp ) THEN
     1671
     1672!--       Case winter clothing insulation: cold stress ?
     1673          clo = wclo
     1674          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
     1675             top )
     1676          pmv_w = pmva
     1677
     1678          IF ( pmva >= 0._wp ) THEN
     1679
     1680!--          Case: comfort achievable by varying clothing insulation
     1681!            between winter and summer set values
     1682             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
     1683                               pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
     1684             IF ( ncount < 0_iwp ) THEN
     1685                nerr = -1_iwp
     1686                RETURN
     1687             ENDIF
     1688          ELSE IF ( pmva < -0.11_wp ) THEN
     1689             clo = 1.75_wp
     1690             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
     1691                           pmva, top )
     1692          ENDIF
     1693       ELSE IF ( pmva > 0.06_wp ) THEN
     1694          clo = 0.5_wp
     1695          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
     1696             top )
     1697       ENDIF
     1698
     1699    ENDIF
     1700
     1701!-- Determine perceived temperature by regression equation + adjustments
     1702    pmvs = pmva
     1703    CALL perct_regression ( pmva, clo, perct )
     1704    ptc = perct
     1705    IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
     1706!--    Adjust for cold conditions according to Gagge 1986
     1707       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
     1708       IF ( nerr_cold > 0_iwp ) nerr = -5_iwp
     1709       pmvs = pmva - d_pmv
     1710       IF ( pmvs > -0.11_wp ) THEN
     1711          d_pmv  = 0._wp
     1712          pmvs   = -0.11_wp
     1713       ENDIF
     1714       CALL perct_regression ( pmvs, clo, perct )
     1715    ENDIF
     1716!     clo_fanger = clo
     1717    clon = clo
     1718    IF ( clo > 0.5_wp .AND. perct <= 8.73_wp ) THEN
     1719!--    Required clothing insulation (ireq) is exclusively defined for
     1720!      operative temperatures (top) less 10 (C) for a
     1721!      reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
     1722       clon = ireq_neutral ( perct, ireq_minimal, nerr )
     1723       clo = clon
     1724    ENDIF
     1725    CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
     1726    sultrieness    = .FALSE.
     1727    d_std = -99._wp
     1728    IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
     1729!--    Adjust for warm/humid conditions according to Gagge 1986
     1730       CALL saturation_vapor_pressure ( ta, svp_ta )
     1731       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
     1732       pmvs   = pmva + d_pmv
     1733       CALL perct_regression ( pmvs, clo, perct )
     1734       IF ( sult_lim < 99._wp ) THEN
     1735          IF ( (perct - ptc) > sult_lim ) sultrieness = .TRUE.
     1736!--       Set factor to threshold for sultriness
     1737          IF ( dgtcstd /= 0_iwp ) THEN
     1738             d_std = ( ( perct - ptc ) - dgtcm ) / dgtcstd
     1739          ENDIF
     1740       ENDIF
     1741    ENDIF
     1742
     1743 END SUBROUTINE calculate_perct_static
     1744
     1745!------------------------------------------------------------------------------!
     1746! Description:
     1747! ------------
     1748!> The SUBROUTINE calculates the saturation water vapour pressure
     1749!> (hPa = hecto Pascal) for a given temperature ta (degC).
     1750!> For example, ta can be the air temperature or the dew point temperature.
     1751!------------------------------------------------------------------------------!
     1752 SUBROUTINE saturation_vapor_pressure( ta, svp_ta )
     1753
     1754    IMPLICIT NONE
     1755
     1756    REAL(wp), INTENT ( IN )  ::  ta     !< ambient air temperature (degC)
     1757    REAL(wp), INTENT ( OUT ) ::  svp_ta !< saturation water vapour pressure (hPa)
     1758
     1759    REAL(wp)      ::  b
     1760    REAL(wp)      ::  c
     1761
     1762
     1763    IF ( ta < 0._wp ) THEN
     1764!--    ta  < 0 (degC): saturation water vapour pressure over ice
     1765       b = 17.84362_wp
     1766       c = 245.425_wp
     1767    ELSE
     1768!--    ta >= 0 (degC): saturation water vapour pressure over water
     1769       b = 17.08085_wp
     1770       c = 234.175_wp
     1771    ENDIF
     1772
     1773!-- Saturation water vapour pressure
     1774    svp_ta = 6.1078_wp * EXP ( b * ta / ( c + ta ) )
     1775
     1776 END SUBROUTINE saturation_vapor_pressure
     1777
     1778!------------------------------------------------------------------------------!
     1779! Description:
     1780! ------------
     1781!> Find the clothing insulation value clo_res (clo) to make Fanger's Predicted
     1782!> Mean Vote (PMV) equal comfort (pmva=0) for actual meteorological conditions
     1783!> (ta,tmrt, vp, ws, pair) and values of individual's activity level
     1784!------------------------------------------------------------------------------!
     1785 SUBROUTINE iso_ridder( ta, tmrt, vp, ws, pair, actlev, eta, sclo,             &
     1786                       pmv_s, wclo, pmv_w, eps, pmva, top, nerr,               &
     1787                       clo_res )
     1788
     1789    IMPLICIT NONE
     1790
     1791!-- Input variables of argument list:
     1792    REAL(wp), INTENT ( IN )  :: ta     !< Ambient temperature (degC)
     1793    REAL(wp), INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC)
     1794    REAL(wp), INTENT ( IN )  :: vp     !< Water vapour pressure (hPa)
     1795    REAL(wp), INTENT ( IN )  :: ws    !< Wind speed (m/s) 1 m above ground
     1796    REAL(wp), INTENT ( IN )  :: pair       !< Barometric pressure (hPa)
     1797    REAL(wp), INTENT ( IN )  :: actlev   !< Individuals activity level per unit surface area (W/m2)
     1798    REAL(wp), INTENT ( IN )  :: eta      !< Individuals work efficiency (dimensionless)
     1799    REAL(wp), INTENT ( IN )  :: sclo     !< Lower threshold of bracketing clothing insulation (clo)
     1800    REAL(wp), INTENT ( IN )  :: wclo     !< Upper threshold of bracketing clothing insulation (clo)
     1801    REAL(wp), INTENT ( IN )  :: eps      !< (0.05) accuracy in clothing insulation (clo) for
     1802!                                          evaluation the root of Fanger's PMV (pmva=0)
     1803    REAL(wp), INTENT ( IN )  :: pmv_w    !< Fanger's PMV corresponding to wclo
     1804    REAL(wp), INTENT ( IN )  :: pmv_s    !< Fanger's PMV corresponding to sclo
     1805
     1806! Output variables of argument list:
     1807    REAL(wp), INTENT ( OUT ) :: pmva     !< 0 (set to zero, because clo is evaluated for comfort)
     1808    REAL(wp), INTENT ( OUT ) :: top      !< Operative temperature (degC) at found root of Fanger's PMV
     1809    REAL(wp), INTENT ( OUT ) :: clo_res  !< Resulting clothing insulation value (clo)
     1810    INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error status / quality flag
     1811!           nerr >= 0, o.k., and nerr is the number of iterations for
     1812!                              convergence
     1813!           nerr = -1: error = malfunction of Ridder's convergence method
     1814!           nerr = -2: error = maximum iterations (max_iteration) exceeded
     1815!           nerr = -3: error = root not bracketed between sclo and wclo
     1816
     1817!-- Type of program variables
     1818    INTEGER(iwp), PARAMETER  ::  max_iteration = 15_iwp       !< max number of iterations
     1819    REAL(wp),     PARAMETER  ::  guess_0       = -1.11e30_wp  !< initial guess
     1820    REAL(wp) ::  x_ridder    !< current guess for clothing insulation   (clo)
     1821    REAL(wp) ::  clo_lower   !< lower limit of clothing insulation      (clo)
     1822    REAL(wp) ::  clo_upper   !< upper limit of clothing insulation      (clo)
     1823    REAL(wp) ::  x_lower     !< lower guess for clothing insulation     (clo)
     1824    REAL(wp) ::  x_upper     !< upper guess for clothing insulation     (clo)
     1825    REAL(wp) ::  x_average   !< average of x_lower and x_upper          (clo)
     1826    REAL(wp) ::  x_new       !< preliminary result for clothing insulation (clo)
     1827    REAL(wp) ::  y_lower     !< predicted mean vote for summer clothing
     1828    REAL(wp) ::  y_upper     !< predicted mean vote for winter clothing
     1829    REAL(wp) ::  y_average   !< average of y_lower and y_upper
     1830    REAL(wp) ::  y_new       !< preliminary result for pred. mean vote
     1831    REAL(wp) ::  sroot       !< sqrt of PMV-guess
     1832    INTEGER(iwp) ::  j       !< running index
     1833
     1834!-- Initialise
     1835    nerr    = 0_iwp
     1836
     1837!-- Set pmva = 0 (comfort): Root of PMV depending on clothing insulation
     1838    pmva        = 0._wp
     1839    clo_lower   = sclo
     1840    y_lower     = pmv_s
     1841    clo_upper   = wclo
     1842    y_upper     = pmv_w
     1843    IF ( ( y_lower > 0._wp .AND. y_upper < 0._wp ) .OR.                        &
     1844         ( y_lower < 0._wp .AND. y_upper > 0._wp ) ) THEN
     1845       x_lower  = clo_lower
     1846       x_upper  = clo_upper
     1847       x_ridder = guess_0
     1848
     1849       DO j = 1_iwp, max_iteration
     1850          x_average = 0.5_wp * ( x_lower + x_upper )
     1851          CALL fanger ( ta, tmrt, vp, ws, pair, x_average, actlev, eta,        &
     1852                        y_average, top )
     1853          sroot = SQRT ( y_average**2 - y_lower * y_upper )
     1854          IF ( sroot == 0._wp ) THEN
     1855             clo_res = x_average
     1856             nerr = j
     1857             RETURN
     1858          ENDIF
     1859          x_new = x_average + ( x_average - x_lower ) *                        &
     1860                      ( SIGN ( 1._wp, y_lower - y_upper ) * y_average / sroot )
     1861          IF ( ABS ( x_new - x_ridder ) <= eps ) THEN
     1862             clo_res = x_ridder
     1863             nerr       = j
     1864             RETURN
     1865          ENDIF
     1866          x_ridder = x_new
     1867          CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta,         &
     1868                        y_new, top )
     1869          IF ( y_new == 0._wp ) THEN
     1870             clo_res = x_ridder
     1871             nerr       = j
     1872             RETURN
     1873          ENDIF
     1874          IF ( SIGN ( y_average, y_new ) /= y_average ) THEN
     1875             x_lower = x_average
     1876             y_lower = y_average
     1877             x_upper  = x_ridder
     1878             y_upper  = y_new
     1879          ELSE IF ( SIGN ( y_lower, y_new ) /= y_lower ) THEN
     1880             x_upper  = x_ridder
     1881             y_upper  = y_new
     1882          ELSE IF ( SIGN ( y_upper, y_new ) /= y_upper ) THEN
     1883             x_lower = x_ridder
     1884             y_lower = y_new
     1885          ELSE
     1886!--          Never get here in x_ridder: singularity in y
     1887             nerr       = -1_iwp
     1888             clo_res = x_ridder
     1889             RETURN
     1890          ENDIF
     1891          IF ( ABS ( x_upper - x_lower ) <= eps ) THEN
     1892             clo_res = x_ridder
     1893             nerr       = j
     1894             RETURN
     1895          ENDIF
     1896       ENDDO
     1897!--    x_ridder exceed maximum iterations
     1898       nerr       = -2_iwp
     1899       clo_res = y_new
     1900       RETURN
     1901    ELSE IF ( y_lower == 0. ) THEN
     1902       x_ridder = clo_lower
     1903    ELSE IF ( y_upper == 0. ) THEN
     1904       x_ridder = clo_upper
     1905    ELSE
     1906!--    x_ridder not bracketed by u_clo and o_clo
     1907       nerr = -3_iwp
     1908       clo_res = x_ridder
     1909       RETURN
     1910    ENDIF
     1911
     1912 END SUBROUTINE iso_ridder
     1913
     1914!------------------------------------------------------------------------------!
     1915! Description:
     1916! ------------
     1917!> Regression relations between perceived temperature (perct) and (adjusted)
     1918!> PMV. The regression presumes the Klima-Michel settings for reference
     1919!> individual and reference environment.
     1920!------------------------------------------------------------------------------!
     1921 SUBROUTINE perct_regression( pmv, clo, perct )
     1922
     1923    IMPLICIT NONE
     1924
     1925    REAL(wp), INTENT ( IN ) ::  pmv   !< Fangers predicted mean vote (dimensionless)
     1926    REAL(wp), INTENT ( IN ) ::  clo   !< clothing insulation index (clo)
     1927
     1928    REAL(wp), INTENT ( OUT ) ::  perct   !< perct (degC) corresponding to given PMV / clo
     1929
     1930    IF ( pmv <= -0.11_wp ) THEN
     1931       perct = 5.805_wp + 12.6784_wp * pmv
     1932    ELSE
     1933       IF ( pmv >= + 0.01_wp ) THEN
     1934          perct = 16.826_wp + 6.163_wp * pmv
     1935       ELSE
     1936          perct = 21.258_wp - 9.558_wp * clo
     1937       ENDIF
     1938    ENDIF
     1939
     1940 END SUBROUTINE perct_regression
     1941
     1942!------------------------------------------------------------------------------!
     1943! Description:
     1944! ------------
     1945!> FANGER.F90
     1946!>
     1947!> SI-VERSION: ACTLEV W m-2, DAMPFDRUCK hPa
     1948!> Berechnet das aktuelle Predicted Mean Vote nach Fanger
     1949!>
     1950!> The case of free convection (ws < 0.1 m/s) is dealt with ws = 0.1 m/s
     1951!------------------------------------------------------------------------------!
     1952 SUBROUTINE fanger( ta, tmrt, pa, in_ws, pair, in_clo, actlev, eta, pmva, top )
     1953
     1954    IMPLICIT NONE
     1955
     1956!-- Input variables of argument list:
     1957    REAL(wp), INTENT ( IN ) ::  ta       !< Ambient air temperature (degC)
     1958    REAL(wp), INTENT ( IN ) ::  tmrt     !< Mean radiant temperature (degC)
     1959    REAL(wp), INTENT ( IN ) ::  pa       !< Water vapour pressure (hPa)
     1960    REAL(wp), INTENT ( IN ) ::  pair     !< Barometric pressure (hPa) at site
     1961    REAL(wp), INTENT ( IN ) ::  in_ws    !< Wind speed (m/s) 1 m above ground
     1962    REAL(wp), INTENT ( IN ) ::  in_clo   !< Clothing insulation (clo)
     1963    REAL(wp), INTENT ( IN ) ::  actlev   !< Individuals activity level per unit surface area (W/m2)
     1964    REAL(wp), INTENT ( IN ) ::  eta      !< Individuals mechanical work efficiency (dimensionless)
     1965
     1966!-- Output variables of argument list:
     1967    REAL(wp), INTENT ( OUT ) ::  pmva    !< Actual Predicted Mean Vote (PMV,
     1968!            dimensionless) according to Fanger corresponding to meteorological
     1969!            (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta)
     1970    REAL(wp), INTENT ( OUT ) ::  top     !< operative temperature (degC)
     1971
     1972!-- Internal variables
     1973    REAL(wp) ::  f_cl         !< Increase in surface due to clothing    (factor)
     1974    REAL(wp) ::  heat_convection  !< energy loss by autocnvection       (W)
     1975    REAL(wp) ::  activity     !< persons activity  (must stay == actlev, W)
     1976    REAL(wp) ::  t_skin_aver  !< average skin temperature               (°C)
     1977    REAL(wp) ::  bc           !< preliminary result storage
     1978    REAL(wp) ::  cc           !< preliminary result storage
     1979    REAL(wp) ::  dc           !< preliminary result storage
     1980    REAL(wp) ::  ec           !< preliminary result storage
     1981    REAL(wp) ::  gc           !< preliminary result storage
     1982    REAL(wp) ::  t_clothing   !< clothing temperature                   (°C)
     1983    REAL(wp) ::  hr           !< radiational heat resistence
     1984    REAL(wp) ::  clo          !< clothing insulation index              (clo)
     1985    REAL(wp) ::  ws           !< wind speed                             (m/s)
     1986    REAL(wp) ::  z1           !< Empiric factor for the adaption of the heat
     1987!             ballance equation to the psycho-physical scale (Equ. 40 in FANGER)
     1988    REAL(wp) ::  z2           !< Water vapour diffution through the skin
     1989    REAL(wp) ::  z3           !< Sweat evaporation from the skin surface
     1990    REAL(wp) ::  z4           !< Loss of latent heat through respiration
     1991    REAL(wp) ::  z5           !< Loss of radiational heat
     1992    REAL(wp) ::  z6           !< Heat loss through forced convection
     1993    INTEGER(iwp) :: i         !< running index
     1994
     1995!-- Clo must be > 0. to avoid div. by 0!
     1996    clo = in_clo
     1997    IF ( clo <= 0._wp ) clo = .001_wp
     1998
     1999!-- f_cl = Increase in surface due to clothing
     2000    f_cl = 1._wp + .15_wp * clo
     2001
     2002!-- Case of free convection (ws < 0.1 m/s ) not considered
     2003    ws = in_ws
     2004    IF ( ws < .1_wp ) THEN
     2005       ws = .1_wp
     2006    ENDIF
     2007
     2008!-- Heat_convection = forced convection
     2009    heat_convection = 12.1_wp * SQRT ( ws * pair / 1013.25_wp )
     2010
     2011!-- Activity = inner heat produktion per standardized surface
     2012    activity = actlev * ( 1._wp - eta )
     2013
     2014!-- T_skin_aver = average skin temperature
     2015    t_skin_aver = 35.7_wp - .0275_wp * activity
     2016
     2017!-- Calculation of constants for evaluation below
     2018    bc = .155_wp * clo * 3.96_wp * 10._wp**( -8 ) * f_cl
     2019    cc = f_cl * heat_convection
     2020    ec = .155_wp * clo
     2021    dc = ( 1._wp + ec * cc ) / bc
     2022    gc = ( t_skin_aver + bc * ( tmrt + degc_to_k )**4 + ec * cc * ta ) / bc
     2023
     2024!-- Calculation of clothing surface temperature (t_clothing) based on
     2025!   Newton-approximation with air temperature as initial guess
     2026    t_clothing = ta
     2027    DO i = 1, 3
     2028       t_clothing = t_clothing - ( ( t_clothing + degc_to_k )**4 + t_clothing  &
     2029          * dc - gc ) / ( 4._wp * ( t_clothing + degc_to_k )**3 + dc )
     2030    ENDDO
     2031
     2032!-- Empiric factor for the adaption of the heat ballance equation
     2033!   to the psycho-physical scale (Equ. 40 in FANGER)
     2034    z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp )
     2035
     2036!-- Water vapour diffution through the skin
     2037    z2 = .31_wp * ( 57.3_wp - .07_wp * activity-pa )
     2038
     2039!-- Sweat evaporation from the skin surface
     2040    z3 = .42_wp * ( activity - 58._wp )
     2041
     2042!-- Loss of latent heat through respiration
     2043    z4 = .0017_wp * actlev * ( 58.7_wp - pa ) + .0014_wp * actlev *            &
     2044      ( 34._wp - ta )
     2045
     2046!-- Loss of radiational heat
     2047    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + degc_to_k )**4 - ( tmrt +        &
     2048       degc_to_k )**4 )
     2049    IF ( ABS ( t_clothing - tmrt ) > 0._wp ) THEN
     2050       hr = z5 / f_cl / ( t_clothing - tmrt )
     2051    ELSE
     2052       hr = 0._wp
     2053    ENDIF
     2054
     2055!-- Heat loss through forced convection cc*(t_clothing-TT)
     2056    z6 = cc * ( t_clothing - ta )
     2057
     2058!-- Predicted Mean Vote
     2059    pmva = z1 * ( activity - z2 - z3 - z4 - z5 - z6 )
     2060
     2061!-- Operative temperatur
     2062    top = ( hr * tmrt + heat_convection * ta ) / ( hr + heat_convection )
     2063
     2064 END SUBROUTINE fanger
     2065
     2066!------------------------------------------------------------------------------!
     2067! Description:
     2068! ------------
     2069!> For pmva > 0 and clo =0.5 the increment (deltapmv) is calculated
     2070!> that converts pmva into Gagge's et al. (1986) PMV*.
     2071!------------------------------------------------------------------------------!
     2072 REAL(wp) FUNCTION deltapmv( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
     2073
     2074    IMPLICIT NONE
     2075
     2076!-- Input variables of argument list:
     2077    REAL(wp),     INTENT ( IN )  :: pmva     !< Actual Predicted Mean Vote (PMV) according to Fanger
     2078    REAL(wp),     INTENT ( IN )  :: ta       !< Ambient temperature (degC) at screen level
     2079    REAL(wp),     INTENT ( IN )  :: vp       !< Water vapour pressure (hPa) at screen level
     2080    REAL(wp),     INTENT ( IN )  :: svp_ta   !< Saturation water vapour pressure (hPa) at ta
     2081    REAL(wp),     INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC) at screen level
     2082    REAL(wp),     INTENT ( IN )  :: ws       !< Wind speed (m/s) 1 m above ground
     2083
     2084!-- Output variables of argument list:
     2085    INTEGER(iwp), INTENT ( OUT ) :: nerr     !< Error status / quality flag
     2086!             0 = o.k.
     2087!            -2 = pmva outside valid regression range
     2088!            -3 = rel. humidity set to 5 % or 95 %, respectively
     2089!            -4 = deltapmv set to avoid pmvs < 0
     2090
     2091!-- Internal variable types:
     2092    REAL(wp) ::  pmv          !<
     2093    REAL(wp) ::  pa_p50       !<
     2094    REAL(wp) ::  pa           !<
     2095    REAL(wp) ::  apa          !<
     2096    REAL(wp) ::  dapa         !<
     2097    REAL(wp) ::  sqvel        !<
     2098    REAL(wp) ::  dtmrt        !<
     2099    REAL(wp) ::  p10          !<
     2100    REAL(wp) ::  p95          !<
     2101    REAL(wp) ::  gew          !<
     2102    REAL(wp) ::  gew2         !<
     2103    REAL(wp) ::  dpmv_1       !<
     2104    REAL(wp) ::  dpmv_2       !<
     2105    REAL(wp) ::  pmvs         !<
     2106    REAL(wp) ::  bpmv(0:7)    !<
     2107    REAL(wp) ::  bpa_p50(0:7) !<
     2108    REAL(wp) ::  bpa(0:7)     !<
     2109    REAL(wp) ::  bapa(0:7)    !<
     2110    REAL(wp) ::  bdapa(0:7)   !<
     2111    REAL(wp) ::  bsqvel(0:7)  !<
     2112    REAL(wp) ::  bta(0:7)     !<
     2113    REAL(wp) ::  bdtmrt(0:7)  !<
     2114    REAL(wp) ::  aconst(0:7)  !<
     2115    INTEGER(iwp) :: nreg      !<
     2116
     2117    DATA bpmv     /                                                            &
     2118     -0.0556602_wp, -0.1528680_wp, -0.2336104_wp, -0.2789387_wp, -0.3551048_wp,&
     2119     -0.4304076_wp, -0.4884961_wp, -0.4897495_wp /
     2120    DATA bpa_p50 /                                                             &
     2121     -0.1607154_wp, -0.4177296_wp, -0.4120541_wp, -0.0886564_wp, +0.4285938_wp,&
     2122     +0.6281256_wp, +0.5067361_wp, +0.3965169_wp /
     2123    DATA bpa     /                                                             &
     2124     +0.0580284_wp, +0.0836264_wp, +0.1009919_wp, +0.1020777_wp, +0.0898681_wp,&
     2125     +0.0839116_wp, +0.0853258_wp, +0.0866589_wp /
     2126    DATA bapa    /                                                             &
     2127     -1.7838788_wp, -2.9306231_wp, -1.6350334_wp, +0.6211547_wp, +3.3918083_wp,&
     2128     +5.5521025_wp, +8.4897418_wp, +16.6265851_wp /
     2129    DATA bdapa   /                                                             &
     2130     +1.6752720_wp, +2.7379504_wp, +1.2940526_wp, -1.0985759_wp, -3.9054732_wp,&
     2131     -6.0403012_wp, -8.9437119_wp, -17.0671201_wp /
     2132    DATA bsqvel  /                                                             &
     2133     -0.0315598_wp, -0.0286272_wp, -0.0009228_wp, +0.0483344_wp, +0.0992366_wp,&
     2134     +0.1491379_wp, +0.1951452_wp, +0.2133949_wp /
     2135    DATA bta     /                                                             &
     2136     +0.0953986_wp, +0.1524760_wp, +0.0564241_wp, -0.0893253_wp, -0.2398868_wp,&
     2137     -0.3515237_wp, -0.5095144_wp, -0.9469258_wp /
     2138    DATA bdtmrt  /                                                             &
     2139     -0.0004672_wp, -0.0000514_wp, -0.0018037_wp, -0.0049440_wp, -0.0069036_wp,&
     2140     -0.0075844_wp, -0.0079602_wp, -0.0089439_wp /
     2141    DATA aconst  /                                                             &
     2142     +1.8686215_wp, +3.4260713_wp, +2.0116185_wp, -0.7777552_wp, -4.6715853_wp,&
     2143     -7.7314281_wp, -11.7602578_wp, -23.5934198_wp /
     2144
     2145!-- Test for compliance with regression range
     2146    IF ( pmva < -1.0_wp .OR. pmva > 7.0_wp ) THEN
     2147       nerr = -2_iwp
     2148    ELSE
     2149       nerr = 0_iwp
     2150    ENDIF
     2151
     2152!-- Initialise classic PMV
     2153    pmv  = pmva
     2154
     2155!-- Water vapour pressure of air
     2156    p10  = 0.05_wp * svp_ta
     2157    p95  = 1.00_wp * svp_ta
     2158    IF ( vp >= p10 .AND. vp <= p95 ) THEN
     2159       pa = vp
     2160    ELSE
     2161       nerr = -3_iwp
     2162       IF ( vp < p10 ) THEN
     2163!--       Due to conditions of regression: r.H. >= 5 %
     2164          pa = p10
     2165       ELSE
     2166!--       Due to conditions of regression: r.H. <= 95 %
     2167          pa = p95
     2168       ENDIF
     2169    ENDIF
     2170    IF ( pa > 0._wp ) THEN
     2171!--    Natural logarithm of pa
     2172       apa = LOG ( pa )
     2173    ELSE
     2174       apa = -5._wp
     2175    ENDIF
     2176
     2177!-- Ratio actual water vapour pressure to that of a r.H. of 50 %
     2178    pa_p50   = 0.5_wp * svp_ta
     2179    IF ( pa_p50 > 0._wp .AND. pa > 0._wp ) THEN
     2180       dapa   = apa - LOG ( pa_p50 )
     2181       pa_p50 = pa / pa_p50
     2182    ELSE
     2183       dapa   = -5._wp
     2184       pa_p50 = 0._wp
     2185    ENDIF
     2186!-- Square root of wind velocity
     2187    IF ( ws >= 0._wp ) THEN
     2188       sqvel = SQRT ( ws )
     2189    ELSE
     2190       sqvel = 0._wp
     2191    ENDIF
     2192!-- Air temperature
     2193!    ta = ta
     2194!-- Difference mean radiation to air temperature
     2195    dtmrt = tmrt - ta
     2196
     2197!-- Select the valid regression coefficients
     2198    nreg = INT ( pmv )
     2199    IF ( nreg < 0_iwp ) THEN
     2200!--    value of the FUNCTION in the case pmv <= -1
     2201       deltapmv = 0._wp
     2202       RETURN
     2203    ENDIF
     2204    gew = MOD ( pmv, 1._wp )
     2205    IF ( gew < 0._wp ) gew = 0._wp
     2206    IF ( nreg > 5_iwp ) THEN
     2207       ! nreg=6
     2208       nreg  = 5_iwp
     2209       gew   = pmv - 5._wp
     2210       gew2  = pmv - 6._wp
     2211       IF ( gew2 > 0_iwp ) THEN
     2212          gew = ( gew - gew2 ) / gew
     2213       ENDIF
     2214    ENDIF
     2215
     2216!-- Regression valid for 0. <= pmv <= 6.
     2217    dpmv_1 =                                                                   &
     2218       + bpa ( nreg ) * pa                                                     &
     2219       + bpmv ( nreg ) * pmv                                                   &
     2220       + bapa ( nreg ) * apa                                                   &
     2221       + bta ( nreg ) * ta                                                     &
     2222       + bdtmrt ( nreg ) * dtmrt                                               &
     2223       + bdapa ( nreg ) * dapa                                                 &
     2224       + bsqvel ( nreg ) * sqvel                                               &
     2225       + bpa_p50 ( nreg ) * pa_p50                                             &
     2226       + aconst ( nreg )
     2227
     2228    dpmv_2 = 0._wp
     2229    IF ( nreg < 6_iwp ) THEN
     2230       dpmv_2 =                                                                &
     2231          + bpa ( nreg + 1_iwp )     * pa                                      &
     2232          + bpmv ( nreg + 1_iwp )    * pmv                                     &
     2233          + bapa ( nreg + 1_iwp )    * apa                                     &
     2234          + bta ( nreg + 1_iwp )     * ta                                      &
     2235          + bdtmrt ( nreg + 1_iwp )  * dtmrt                                   &
     2236          + bdapa ( nreg + 1_iwp )   * dapa                                    &
     2237          + bsqvel ( nreg + 1_iwp )  * sqvel                                   &
     2238          + bpa_p50 ( nreg + 1_iwp ) * pa_p50                                  &
     2239          + aconst ( nreg + 1_iwp )
     2240    ENDIF
     2241
     2242!-- Calculate pmv modification
     2243    deltapmv = ( 1._wp - gew ) * dpmv_1 + gew * dpmv_2
     2244    pmvs = pmva + deltapmv
     2245    IF ( ( pmvs ) < 0._wp ) THEN
     2246!--    Prevent negative pmv* due to problems with clothing insulation
     2247       nerr = -4_iwp
     2248       IF ( pmvs > -0.11_wp ) THEN
     2249!--       Threshold from FUNCTION perct_regression for winter clothing insulation
     2250          deltapmv = deltapmv + 0.11_wp
     2251       ELSE
     2252!--       Set pmvs to "0" for compliance with summer clothing insulation
     2253          deltapmv = -1._wp * pmva
     2254       ENDIF
     2255    ENDIF
     2256
     2257 END FUNCTION deltapmv
     2258
     2259!------------------------------------------------------------------------------!
     2260! Description:
     2261! ------------
     2262!> The subroutine "calc_sultr" returns a threshold value to perceived
     2263!> temperature allowing to decide whether the actual perceived temperature
     2264!> is linked to perecption of sultriness. The threshold values depends
     2265!> on the Fanger's classical PMV, expressed here as perceived temperature
     2266!> perct.
     2267!------------------------------------------------------------------------------!
     2268 SUBROUTINE calc_sultr( perct, dperctm, dperctstd, sultr_res )
     2269
     2270    IMPLICIT NONE
     2271
     2272!-- Input of the argument list:
     2273    REAL(wp), INTENT ( IN )  ::  perct      !< Classical perceived temperature: Base is Fanger's PMV
     2274
     2275!-- Additional output variables of argument list:
     2276    REAL(wp), INTENT ( OUT ) ::  dperctm    !< Mean deviation perct (classical gt) to gt* (rational gt
     2277!             calculated based on Gagge's rational PMV*)
     2278    REAL(wp), INTENT ( OUT ) ::  dperctstd  !< dperctm plus its standard deviation times a factor
     2279!             determining the significance to perceive sultriness
     2280    REAL(wp), INTENT ( OUT ) ::  sultr_res
     2281
     2282!-- Types of coefficients mean deviation: third order polynomial
     2283    REAL(wp), PARAMETER ::  dperctka = +7.5776086_wp
     2284    REAL(wp), PARAMETER ::  dperctkb = -0.740603_wp
     2285    REAL(wp), PARAMETER ::  dperctkc = +0.0213324_wp
     2286    REAL(wp), PARAMETER ::  dperctkd = -0.00027797237_wp
     2287
     2288!-- Types of coefficients mean deviation plus standard deviation
     2289!   regression coefficients: third order polynomial
     2290    REAL(wp), PARAMETER ::  dperctsa = +0.0268918_wp
     2291    REAL(wp), PARAMETER ::  dperctsb = +0.0465957_wp
     2292    REAL(wp), PARAMETER ::  dperctsc = -0.00054709752_wp
     2293    REAL(wp), PARAMETER ::  dperctsd = +0.0000063714823_wp
     2294
     2295!-- Factor to mean standard deviation defining SIGNificance for
     2296!   sultriness
     2297    REAL(wp), PARAMETER :: faktor = 1._wp
     2298
     2299!-- Initialise
     2300    sultr_res = +99._wp
     2301    dperctm   = 0._wp
     2302    dperctstd = 999999._wp
     2303
     2304    IF ( perct < 16.826_wp .OR. perct > 56._wp ) THEN
     2305!--    Unallowed classical PMV/perct
     2306       RETURN
     2307    ENDIF
     2308
     2309!-- Mean deviation dependent on perct
     2310    dperctm = dperctka + dperctkb * perct + dperctkc * perct**2._wp + dperctkd * perct**3._wp
     2311
     2312!-- Mean deviation plus its standard deviation
     2313    dperctstd = dperctsa + dperctsb * perct + dperctsc * perct**2._wp + dperctsd * perct**3._wp
     2314
     2315!-- Value of the FUNCTION
     2316    sultr_res = dperctm + faktor * dperctstd
     2317    IF ( ABS ( sultr_res ) > 99._wp ) sultr_res = +99._wp
     2318
     2319 END SUBROUTINE calc_sultr
     2320
     2321!------------------------------------------------------------------------------!
     2322! Description:
     2323! ------------
     2324!> Multiple linear regression to calculate an increment delta_cold,
     2325!> to adjust Fanger's classical PMV (pmva) by Gagge's 2 node model,
     2326!> applying Fanger's convective heat transfer coefficient, hcf.
     2327!> Wind velocitiy of the reference environment is 0.10 m/s
     2328!------------------------------------------------------------------------------!
     2329 SUBROUTINE dpmv_cold( pmva, ta, ws, tmrt, nerr, dpmv_cold_res )
     2330
     2331    IMPLICIT NONE
     2332
     2333!-- Type of input arguments
     2334    REAL(wp), INTENT ( IN ) ::  pmva   !< Fanger's classical predicted mean vote
     2335    REAL(wp), INTENT ( IN ) ::  ta     !< Air temperature 2 m above ground (degC)
     2336    REAL(wp), INTENT ( IN ) ::  ws     !< Relative wind velocity 1 m above ground (m/s)
     2337    REAL(wp), INTENT ( IN ) ::  tmrt   !< Mean radiant temperature (degC)
     2338
     2339!-- Type of output argument
     2340    INTEGER(iwp), INTENT ( OUT ) ::  nerr !< Error indicator: 0 = o.k., +1 = denominator for intersection = 0
     2341    REAL(wp),     INTENT ( OUT ) ::  dpmv_cold_res    !< Increment to adjust pmva according to the results of Gagge's
     2342!               2 node model depending on the input
     2343
     2344!-- Type of program variables
     2345    REAL(wp) ::  delta_cold(3)
     2346    REAL(wp) ::  pmv_cross(2)
     2347    REAL(wp) ::  reg_a(3)
     2348    REAL(wp) ::  coeff(3,5)
     2349    REAL(wp) ::  r_nenner
     2350    REAL(wp) ::  pmvc
     2351    REAL(wp) ::  dtmrt
     2352    REAL(wp) ::  sqrt_ws
     2353    INTEGER(iwp) ::  i
     2354    INTEGER(iwp) ::  j
     2355    INTEGER(iwp) ::  i_bin
     2356
     2357!-- Coefficient of the 3 regression lines
     2358    !     1:const   2:*pmvc  3:*ta      4:*sqrt_ws  5:*dtmrt
     2359    coeff(1,1:5) =                                                             &
     2360       (/ +0.161_wp, +0.130_wp, -1.125E-03_wp, +1.106E-03_wp, -4.570E-04_wp /)
     2361    coeff(2,1:5) =                                                             &
     2362       (/ 0.795_wp, 0.713_wp, -8.880E-03_wp, -1.803E-03_wp, -2.816E-03_wp/)
     2363    coeff(3,1:5) =                                                             &
     2364       (/ +0.05761_wp, +0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp /)
     2365
     2366!-- Initialise
     2367    nerr       = 0_iwp
     2368    dpmv_cold_res  = 0._wp
     2369    pmvc       = pmva
     2370    dtmrt      = tmrt - ta
     2371    sqrt_ws   = ws
     2372    IF ( sqrt_ws < 0.10_wp ) THEN
     2373       sqrt_ws = 0.10_wp
     2374    ELSE
     2375       sqrt_ws = SQRT ( sqrt_ws )
     2376    ENDIF
     2377
     2378    DO i = 1, 3
     2379       delta_cold (i) = 0._wp
     2380       IF ( i < 3 ) THEN
     2381          pmv_cross (i) = pmvc
     2382       ENDIF
     2383    ENDDO
     2384
     2385    DO i = 1, 3
     2386!--    Regression constant for given meteorological variables
     2387       reg_a(i) = coeff(i, 1) + coeff(i, 3) * ta + coeff(i, 4) *             &
     2388                  sqrt_ws + coeff(i,5)*dtmrt
     2389       delta_cold(i) = reg_a(i) + coeff(i, 2) * pmvc
     2390    ENDDO
     2391
     2392!-- Intersection points of regression lines in terms of Fanger's PMV
     2393    DO i = 1, 2
     2394       r_nenner = coeff (i, 2) - coeff (i + 1, 2)
     2395       IF ( ABS ( r_nenner ) > 0.00001_wp ) THEN
     2396          pmv_cross(i) = ( reg_a (i + 1) - reg_a (i) ) / r_nenner
     2397       ELSE
     2398          nerr = 1_iwp
     2399          RETURN
     2400       ENDIF
     2401    ENDDO
     2402
     2403    i_bin = 3
     2404    DO i = 1, 2
     2405       IF ( pmva > pmv_cross (i) ) THEN
     2406          i_bin = i
     2407          EXIT
     2408       ENDIF
     2409    ENDDO
     2410!-- Adjust to operative temperature scaled according
     2411!   to classical PMV (Fanger)
     2412    dpmv_cold_res = delta_cold(i_bin) - dpmv_adj(pmva)
     2413
     2414 END SUBROUTINE dpmv_cold
     2415
     2416!------------------------------------------------------------------------------!
     2417! Description:
     2418! ------------
     2419!> Calculates the summand dpmv_adj adjusting to the operative temperature
     2420!> scaled according to classical PMV (Fanger)
     2421!> Reference environment: v_1m = 0.10 m/s, dTMRT = 0 K, r.h. = 50 %
     2422!------------------------------------------------------------------------------!
     2423 REAL(wp) FUNCTION dpmv_adj( pmva )
     2424
     2425    IMPLICIT NONE
     2426
     2427    REAL(wp), INTENT ( IN ) ::  pmva
     2428    INTEGER(iwp), PARAMETER ::  n_bin = 3
     2429    INTEGER(iwp), PARAMETER ::  n_ca = 0
     2430    INTEGER(iwp), PARAMETER ::  n_ce = 3
     2431    REAL(wp), dimension (n_bin, n_ca:n_ce) ::  coef
     2432
     2433    REAL(wp)      ::  pmv
     2434    INTEGER(iwp)  ::  i, i_bin
     2435
     2436!                             range_1     range_2     range_3
     2437    DATA (coef(i, 0), i = 1, n_bin) /0.0941540_wp, -0.1506620_wp, -0.0871439_wp/
     2438    DATA (coef(i, 1), i = 1, n_bin) /0.0783162_wp, -1.0612651_wp,  0.1695040_wp/
     2439    DATA (coef(i, 2), i = 1, n_bin) /0.1350144_wp, -1.0049144_wp, -0.0167627_wp/
     2440    DATA (coef(i, 3), i = 1, n_bin) /0.1104037_wp, -0.2005277_wp, -0.0003230_wp/
     2441
     2442    IF ( pmva <= -2.1226_wp ) THEN
     2443       i_bin = 3_iwp
     2444    ELSE IF ( pmva <= -1.28_wp ) THEN
     2445       i_bin = 2_iwp
     2446    ELSE
     2447       i_bin = 1_iwp
     2448    ENDIF
     2449
     2450    dpmv_adj   = coef( i_bin, n_ca )
     2451    pmv        = 1._wp
     2452
     2453    DO i = n_ca + 1, n_ce
     2454       pmv      = pmv * pmva
     2455       dpmv_adj = dpmv_adj + coef(i_bin, i) * pmv
     2456    ENDDO
     2457    RETURN
     2458 END FUNCTION dpmv_adj
     2459
     2460!------------------------------------------------------------------------------!
     2461! Description:
     2462! ------------
     2463!> Based on perceived temperature (perct) as input, ireq_neutral determines
     2464!> the required clothing insulation (clo) for thermally neutral conditions
     2465!> (neither body cooling nor body heating). It is related to the Klima-
     2466!> Michel activity level (134.682 W/m2). IREQ_neutral is only defined
     2467!> for perct < 10 (degC)
     2468!------------------------------------------------------------------------------!
     2469 REAL(wp) FUNCTION ireq_neutral( perct, ireq_minimal, nerr )
     2470
     2471    IMPLICIT NONE
     2472
     2473!-- Type declaration of arguments
     2474    REAL(wp),     INTENT ( IN )  ::  perct
     2475    REAL(wp),     INTENT ( OUT ) ::  ireq_minimal
     2476    INTEGER(iwp), INTENT ( OUT ) ::  nerr
     2477
     2478!-- Type declaration for internal varables
     2479    REAL(wp)                     ::  top02
     2480
     2481!-- Initialise
     2482    nerr = 0_iwp
     2483
     2484!-- Convert perceived temperature from basis 0.1 m/s to basis 0.2 m/s
     2485    top02 = 1.8788_wp + 0.9296_wp * perct
     2486
     2487!-- IREQ neutral conditions (thermal comfort)
     2488    ireq_neutral = 1.62_wp - 0.0564_wp * top02
     2489
     2490!-- Regression only defined for perct <= 10 (degC)
     2491    IF ( ireq_neutral < 0.5_wp ) THEN
     2492       IF ( ireq_neutral < 0.48_wp ) THEN
     2493          nerr = 1_iwp
     2494       ENDIF
     2495       ireq_neutral = 0.5_wp
     2496    ENDIF
     2497
     2498!-- Minimal required clothing insulation: maximal acceptable body cooling
     2499    ireq_minimal = 1.26_wp - 0.0588_wp * top02
     2500    IF ( nerr > 0_iwp ) THEN
     2501       ireq_minimal = ireq_neutral
     2502    ENDIF
     2503
     2504    RETURN
     2505 END FUNCTION ireq_neutral
     2506
     2507
     2508!------------------------------------------------------------------------------!
     2509! Description:
     2510! ------------
     2511!> The SUBROUTINE surface area calculates the surface area of the individual
     2512!> according to its height (m), weight (kg), and age (y)
     2513!------------------------------------------------------------------------------!
     2514 SUBROUTINE surface_area ( height_cm, weight, age, surf )
     2515
     2516    IMPLICIT NONE
     2517
     2518    REAL(wp)    , INTENT(in)  ::  weight
     2519    REAL(wp)    , INTENT(in)  ::  height_cm
     2520    INTEGER(iwp), INTENT(in)  ::  age
     2521    REAL(wp)    , INTENT(out) ::  surf
     2522    REAL(wp)                  ::  height
     2523
     2524    height = height_cm * 100._wp
     2525
     2526!-- According to Gehan-George, for children
     2527    IF ( age < 19_iwp ) THEN
     2528       IF ( age < 5_iwp ) THEN
     2529          surf = 0.02667_wp * height**0.42246_wp * weight**0.51456_wp
     2530          RETURN
     2531       END IF
     2532       surf = 0.03050_wp * height**0.35129_wp * weight**0.54375_wp
     2533       RETURN
     2534    END IF
     2535
     2536!-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if
     2537!   height and weight be known. In: Arch. Int. Med.. 17, 1916, S. 863?871.
     2538    surf = 0.007184_wp * height**0.725_wp * weight**0.425_wp
     2539    RETURN
     2540
     2541 END SUBROUTINE surface_area
     2542
     2543!------------------------------------------------------------------------------!
     2544! Description:
     2545! ------------
     2546!> The SUBROUTINE persdat calculates
     2547!>  - the total internal energy production = metabolic + workload,
     2548!>  - the total internal energy production for a standardized surface (actlev)
     2549!>  - the DuBois - area (a_surf [m2])
     2550!> from
     2551!>  - the persons age (years),
     2552!>  - weight (kg),
     2553!>  - height (m),
     2554!>  - sex (1 = male, 2 = female),
     2555!>  - work load (W)
     2556!> for a sample human.
     2557!------------------------------------------------------------------------------!
     2558 SUBROUTINE persdat ( age, weight, height, sex, work, a_surf, actlev )
     2559!
     2560    IMPLICIT NONE
     2561
     2562    REAL(wp), INTENT(in) ::  age
     2563    REAL(wp), INTENT(in) ::  weight
     2564    REAL(wp), INTENT(in) ::  height
     2565    REAL(wp), INTENT(in) ::  work
     2566    INTEGER(iwp), INTENT(in) ::  sex
     2567    REAL(wp), INTENT(out) ::  actlev
     2568    REAL(wp) ::  a_surf
     2569    REAL(wp) ::  energy_prod
     2570    REAL(wp) ::  s
     2571    REAL(wp) ::  factor
     2572    REAL(wp) ::  basic_heat_prod
     2573
     2574    CALL surface_area ( height, weight, INT( age ), a_surf )
     2575    s = height * 100._wp / ( weight ** ( 1._wp / 3._wp ) )
     2576    factor = 1._wp + .004_wp  * ( 30._wp - age )
     2577    basic_heat_prod = 0.
     2578    IF ( sex == 1_iwp ) THEN
     2579       basic_heat_prod = 3.45_wp * weight ** ( 3._wp / 4._wp ) * ( factor +    &
     2580                     .01_wp  * ( s - 43.4_wp ) )
     2581    ELSE IF ( sex == 2_iwp ) THEN
     2582       basic_heat_prod = 3.19_wp * weight ** ( 3._wp / 4._wp ) * ( factor +    &
     2583                    .018_wp * ( s - 42.1_wp ) )
     2584    END IF
     2585
     2586    energy_prod = work + basic_heat_prod
     2587    actlev = energy_prod / a_surf
     2588
     2589 END SUBROUTINE persdat
     2590
     2591
     2592!------------------------------------------------------------------------------!
     2593! Description:
     2594! ------------
     2595!> SUBROUTINE ipt_init
     2596!> initializes the instationary perceived temperature
     2597!------------------------------------------------------------------------------!
     2598
     2599 SUBROUTINE ipt_init (age, weight, height, sex, work, actlev, clo,             &
     2600     ta, vp, ws, tmrt, pair, dt, storage, t_clothing,                         &
     2601     ipt )
     2602
     2603    IMPLICIT NONE
     2604
     2605!-- Input parameters
     2606    REAL(wp), INTENT(in) ::  age        !< Persons age          (years)
     2607    REAL(wp), INTENT(in) ::  weight     !< Persons weight       (kg)
     2608    REAL(wp), INTENT(in) ::  height     !< Persons height       (m)
     2609    REAL(wp), INTENT(in) ::  work       !< Current workload     (W)
     2610    REAL(wp), INTENT(in) ::  ta         !< Air Temperature      (°C)
     2611    REAL(wp), INTENT(in) ::  vp         !< Vapor pressure       (hPa)
     2612    REAL(wp), INTENT(in) ::  ws         !< Wind speed in approx. 1.1m (m/s)
     2613    REAL(wp), INTENT(in) ::  tmrt       !< Mean radiant temperature   (°C)
     2614    REAL(wp), INTENT(in) ::  pair       !< Air pressure         (hPa)
     2615    REAL(wp), INTENT(in) ::  dt         !< Timestep             (s)
     2616    INTEGER(iwp), INTENT(in)  :: sex    !< Persons sex (1 = male, 2 = female)
     2617
     2618!-- Output parameters
     2619    REAL(wp), INTENT(out) ::  actlev
     2620    REAL(wp), INTENT(out) ::  clo
     2621    REAL(wp), INTENT(out) ::  storage
     2622    REAL(wp), INTENT(out) ::  t_clothing
     2623    REAL(wp), INTENT(out) ::  ipt
     2624
     2625!-- Internal variables
     2626    REAL(wp), PARAMETER :: eps = 0.0005_wp
     2627    REAL(wp), PARAMETER :: eta = 0._wp
     2628    REAL(wp) ::  sclo
     2629    REAL(wp) ::  wclo
     2630    REAL(wp) ::  d_pmv
     2631    REAL(wp) ::  svp_ta
     2632    REAL(wp) ::  sult_lim
     2633    REAL(wp) ::  dgtcm
     2634    REAL(wp) ::  dgtcstd
     2635    REAL(wp) ::  clon
     2636    REAL(wp) ::  ireq_minimal
     2637!     REAL(wp) ::  clo_fanger
     2638    REAL(wp) ::  pmv_w
     2639    REAL(wp) ::  pmv_s
     2640    REAL(wp) ::  pmva
     2641    REAL(wp) ::  ptc
     2642    REAL(wp) ::  d_std
     2643    REAL(wp) ::  pmvs
     2644    REAL(wp) ::  top
     2645    REAL(wp) ::  a_surf
     2646    REAL(wp) ::  acti
     2647    INTEGER(iwp) ::  ncount
     2648    INTEGER(iwp) ::  nerr_cold
     2649    INTEGER(iwp) ::  nerr
     2650
     2651    LOGICAL ::  sultrieness
     2652
     2653    storage = 0._wp
     2654    CALL persdat ( age, weight, height, sex, work, a_surf, actlev )
     2655
     2656!-- Initialise
     2657    t_clothing = -999.0_wp
     2658    ipt        = -999.0_wp
     2659    nerr       = 0_wp
     2660    ncount     = 0_wp
     2661    sultrieness    = .FALSE.
     2662!-- Tresholds: clothing insulation (account for model inaccuracies)
     2663!-- Summer clothing
     2664    sclo     = 0.44453_wp
     2665!-- Winter clothing
     2666    wclo     = 1.76267_wp
     2667
     2668!-- Decision: firstly calculate for winter or summer clothing
     2669    IF ( ta <= 10._wp ) THEN
     2670
     2671!--    First guess: winter clothing insulation: cold stress
     2672       clo = wclo
     2673       t_clothing = -999.0_wp  ! force initial run
     2674       CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,         &
     2675          t_clothing, storage, dt, pmva )
     2676       pmv_w = pmva
     2677
     2678       IF ( pmva > 0._wp ) THEN
     2679
     2680!--       Case summer clothing insulation: heat load ?           
     2681          clo = sclo
     2682          t_clothing = -999.0_wp  ! force initial run
     2683          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
     2684             t_clothing, storage, dt, pmva )
     2685          pmv_s = pmva
     2686          IF ( pmva <= 0._wp ) THEN
     2687!--          Case: comfort achievable by varying clothing insulation
     2688!            between winter and summer set values
     2689             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta , sclo,     &
     2690                            pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
     2691             IF ( ncount < 0_iwp ) THEN
     2692                nerr = -1_iwp
     2693                RETURN
     2694             ENDIF
     2695          ELSE IF ( pmva > 0.06_wp ) THEN
     2696             clo = 0.5_wp
     2697             t_clothing = -999.0_wp
     2698             CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,   &
     2699                t_clothing, storage, dt, pmva )
     2700          ENDIF
     2701       ELSE IF ( pmva < -0.11_wp ) THEN
     2702          clo = 1.75_wp
     2703          t_clothing = -999.0_wp
     2704          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
     2705             t_clothing, storage, dt, pmva )
     2706       ENDIF
     2707
     2708    ELSE
     2709
     2710!--    First guess: summer clothing insulation: heat load
     2711       clo = sclo
     2712       t_clothing = -999.0_wp
     2713       CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,         &
     2714          t_clothing, storage, dt, pmva )
     2715       pmv_s = pmva
     2716
     2717       IF ( pmva < 0._wp ) THEN
     2718
     2719!--       Case winter clothing insulation: cold stress ?
     2720          clo = wclo
     2721          t_clothing = -999.0_wp
     2722          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
     2723             t_clothing, storage, dt, pmva )
     2724          pmv_w = pmva
     2725
     2726          IF ( pmva >= 0._wp ) THEN
     2727
     2728!--          Case: comfort achievable by varying clothing insulation
     2729!            between winter and summer set values
     2730             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
     2731                               pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
     2732             IF ( ncount < 0_wp ) THEN
     2733                nerr = -1_iwp
     2734                RETURN
     2735             ENDIF
     2736          ELSE IF ( pmva < -0.11_wp ) THEN
     2737             clo = 1.75_wp
     2738             t_clothing = -999.0_wp
     2739             CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,   &
     2740                t_clothing, storage, dt, pmva )
     2741          ENDIF
     2742       ELSE IF ( pmva > 0.06_wp ) THEN
     2743          clo = 0.5_wp
     2744          t_clothing = -999.0_wp
     2745          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
     2746             t_clothing, storage, dt, pmva )
     2747       ENDIF
     2748
     2749    ENDIF
     2750
     2751!-- Determine perceived temperature by regression equation + adjustments
     2752    pmvs = pmva
     2753    CALL perct_regression ( pmva, clo, ipt )
     2754    ptc = ipt
     2755    IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
     2756!--    Adjust for cold conditions according to Gagge 1986
     2757       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
     2758       IF ( nerr_cold > 0_iwp ) nerr = -5_iwp
     2759       pmvs = pmva - d_pmv
     2760       IF ( pmvs > -0.11_wp ) THEN
     2761          d_pmv  = 0._wp
     2762          pmvs   = -0.11_wp
     2763       ENDIF
     2764       CALL perct_regression ( pmvs, clo, ipt )
     2765    ENDIF
     2766!     clo_fanger = clo
     2767    clon = clo
     2768    IF ( clo > 0.5_wp .AND. ipt <= 8.73_wp ) THEN
     2769!--    Required clothing insulation (ireq) is exclusively defined for
     2770!      operative temperatures (top) less 10 (C) for a
     2771!      reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
     2772       clon = ireq_neutral ( ipt, ireq_minimal, nerr )
     2773       clo = clon
     2774    ENDIF
     2775    CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
     2776    sultrieness    = .FALSE.
     2777    d_std      = -99._wp
     2778    IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
     2779!--    Adjust for warm/humid conditions according to Gagge 1986
     2780       CALL saturation_vapor_pressure ( ta, svp_ta )
     2781       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
     2782       pmvs   = pmva + d_pmv
     2783       CALL perct_regression ( pmvs, clo, ipt )
     2784       IF ( sult_lim < 99._wp ) THEN
     2785          IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE.
     2786       ENDIF
     2787    ENDIF
     2788
     2789 
     2790 END SUBROUTINE ipt_init
     2791 
     2792!------------------------------------------------------------------------------!
     2793! Description:
     2794! ------------
     2795!> SUBROUTINE ipt_cycle
     2796!> Calculates one timestep for the instationary version of perceived
     2797!> temperature (iPT, °C) for
     2798!>  - standard measured/predicted meteorological values and TMRT
     2799!>    as input;
     2800!>  - regressions for determination of PT;
     2801!>  - adjustment to Gagge's PMV* (2-node-model, 1986) as base of PT
     2802!>    under warm/humid conditions (Icl= 0.50 clo) and under cold
     2803!>    conditions (Icl= 1.75 clo)
     2804!>
     2805!------------------------------------------------------------------------------!
     2806 SUBROUTINE ipt_cycle( ta, vp, ws, tmrt, pair, dt, storage, t_clothing, clo,   &
     2807     actlev, work, ipt )
     2808
     2809    IMPLICIT NONE
     2810
     2811!-- Type of input of the argument list
     2812    REAL(wp), INTENT ( IN )  ::  ta      !< Air temperature             (°C)
     2813    REAL(wp), INTENT ( IN )  ::  vp      !< Vapor pressure              (hPa)
     2814    REAL(wp), INTENT ( IN )  ::  tmrt    !< Mean radiant temperature    (°C)
     2815    REAL(wp), INTENT ( IN )  ::  ws      !< Wind speed                  (m/s)
     2816    REAL(wp), INTENT ( IN )  ::  pair    !< Air pressure                (hPa)
     2817    REAL(wp), INTENT ( IN )  ::  dt      !< Timestep                    (s)
     2818    REAL(wp), INTENT ( IN )  ::  clo     !< Clothing index              (no dim)
     2819    REAL(wp), INTENT ( IN )  ::  actlev  !< Internal heat production    (W)
     2820    REAL(wp), INTENT ( IN )  ::  work    !< Mechanical work load        (W)
     2821
     2822!-- In and output parameters
     2823    REAL(wp), INTENT (INOUT) ::  storage     !< Heat storage            (W)
     2824    REAL(wp), INTENT (INOUT) ::  t_clothing  !< Clothig temperature     (°C)
     2825
     2826!-- Type of output of the argument list
     2827    REAL(wp), INTENT ( OUT ) ::  ipt  !< Instationary perceived temperature (°C)
     2828
     2829!-- Type of internal variables
     2830    REAL(wp) ::  d_pmv
     2831    REAL(wp) ::  svp_ta
     2832    REAL(wp) ::  sult_lim
     2833    REAL(wp) ::  dgtcm
     2834    REAL(wp) ::  dgtcstd
     2835    REAL(wp) ::  pmva
     2836    REAL(wp) ::  ptc
     2837    REAL(wp) ::  d_std
     2838    REAL(wp) ::  pmvs
     2839    INTEGER(iwp) ::  nerr_cold
     2840    INTEGER(iwp) ::  nerr
     2841
     2842    LOGICAL ::  sultrieness
     2843
     2844!-- Initialise
     2845    ipt = -999.0_wp
     2846
     2847    nerr     = 0_iwp
     2848    sultrieness  = .FALSE.
     2849
     2850!-- Determine pmv_adjusted for current conditions
     2851    CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,            &
     2852       t_clothing, storage, dt, pmva )
     2853
     2854!-- Determine perceived temperature by regression equation + adjustments
     2855    CALL perct_regression ( pmva, clo, ipt )
     2856
     2857    IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
     2858!--    Adjust for cold conditions according to Gagge 1986
     2859       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
     2860       IF ( nerr_cold > 0_iwp ) nerr = -5_iwp
     2861       pmvs = pmva - d_pmv
     2862       IF ( pmvs > -0.11_wp ) THEN
     2863          d_pmv  = 0._wp
     2864          pmvs   = -0.11_wp
     2865       ENDIF
     2866       CALL perct_regression ( pmvs, clo, ipt )
     2867    ENDIF
     2868
     2869!-- Consider sultriness if appropriate
     2870    ptc = ipt
     2871    CALL calc_sultr ( ptc, dgtcm, dgtcstd, sult_lim )
     2872    sultrieness    = .FALSE.
     2873    d_std      = -99._wp
     2874    IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
     2875!--    Adjust for warm/humid conditions according to Gagge 1986
     2876       CALL saturation_vapor_pressure ( ta, svp_ta )
     2877       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
     2878       pmvs   = pmva + d_pmv
     2879       CALL perct_regression ( pmvs, clo, ipt )
     2880       IF ( sult_lim < 99._wp ) THEN
     2881          IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE.
     2882       ENDIF
     2883    ENDIF
     2884
     2885 END SUBROUTINE ipt_cycle
     2886
     2887!------------------------------------------------------------------------------!
     2888! Description:
     2889! ------------
     2890!> SUBROUTINE fanger_s calculates the
     2891!> actual Predicted Mean Vote (dimensionless) according
     2892!> to Fanger corresponding to meteorological (ta,tmrt,pa,ws,pair)
     2893!> and individual variables (clo, actlev, eta) considering a storage
     2894!> and clothing temperature for a given timestep.
     2895!------------------------------------------------------------------------------!
     2896 SUBROUTINE fanger_s_acti( ta, tmrt, pa, in_ws, pair, in_clo, actlev,          &
     2897    activity, t_cloth, s, dt, pmva )
     2898
     2899    IMPLICIT NONE
     2900
     2901!--  Input argument types
     2902    REAL(wp), INTENT ( IN )  ::  ta       !< Air temperature          (°C)
     2903    REAL(wp), INTENT ( IN )  ::  tmrt     !< Mean radiant temperature (°C)
     2904    REAL(wp), INTENT ( IN )  ::  pa       !< Vapour pressure          (hPa)
     2905    REAL(wp), INTENT ( IN )  ::  pair     !< Air pressure             (hPa)
     2906    REAL(wp), INTENT ( IN )  ::  in_ws   !< Wind speed               (m/s)
     2907    REAL(wp), INTENT ( IN )  ::  actlev   !< Metabolic + work energy  (W/m²)
     2908    REAL(wp), INTENT ( IN )  ::  dt       !< Timestep                 (s)
     2909    REAL(wp), INTENT ( IN )  ::  activity !< Work load                (W/m²)
     2910    REAL(wp), INTENT ( IN )  ::  in_clo   !< Clothing index (clo)     (no dim)
     2911
     2912!-- Output argument types
     2913    REAL(wp), INTENT ( OUT ) ::  pmva  !< actual Predicted Mean Vote (no dim)
     2914
     2915    REAL(wp), INTENT (INOUT) ::  s  !< storage var. of energy balance (W/m2)
     2916    REAL(wp), INTENT (INOUT) ::  t_cloth  !< clothing temperature (°C)
     2917
     2918!-- Internal variables
     2919    REAL(wp), PARAMETER  ::  time_equil = 7200._wp
     2920
     2921    REAL(wp) ::  f_cl         !< Increase in surface due to clothing    (factor)
     2922    REAL(wp) ::  heat_convection  !< energy loss by autocnvection       (W)
     2923    REAL(wp) ::  t_skin_aver  !< average skin temperature               (°C)
     2924    REAL(wp) ::  bc           !< preliminary result storage
     2925    REAL(wp) ::  cc           !< preliminary result storage
     2926    REAL(wp) ::  dc           !< preliminary result storage
     2927    REAL(wp) ::  ec           !< preliminary result storage
     2928    REAL(wp) ::  gc           !< preliminary result storage
     2929    REAL(wp) ::  t_clothing   !< clothing temperature                   (°C)
     2930    REAL(wp) ::  hr           !< radiational heat resistence
     2931    REAL(wp) ::  clo          !< clothing insulation index              (clo)
     2932    REAL(wp) ::  ws           !< wind speed                             (m/s)
     2933    REAL(wp) ::  z1           !< Empiric factor for the adaption of the heat
     2934!             ballance equation to the psycho-physical scale (Equ. 40 in FANGER)
     2935    REAL(wp) ::  z2           !< Water vapour diffution through the skin
     2936    REAL(wp) ::  z3           !< Sweat evaporation from the skin surface
     2937    REAL(wp) ::  z4           !< Loss of latent heat through respiration
     2938    REAL(wp) ::  z5           !< Loss of radiational heat
     2939    REAL(wp) ::  z6           !< Heat loss through forced convection
     2940    REAL(wp) ::  en           !< Energy ballance                        (W)
     2941    REAL(wp) ::  d_s          !< Storage delta                          (W)
     2942    REAL(wp) ::  adjustrate   !< Max storage adjustment rate
     2943    REAL(wp) ::  adjustrate_cloth  !< max clothing temp. adjustment rate
     2944
     2945    INTEGER(iwp) :: i         !< running index
     2946    INTEGER(iwp) ::  niter  !< Running index
     2947
     2948
     2949!-- Clo must be > 0. to avoid div. by 0!
     2950    clo = in_clo
     2951    IF ( clo < 001._wp ) clo = .001_wp
     2952
     2953!-- Increase in surface due to clothing
     2954    f_cl = 1._wp + .15_wp * clo
     2955
     2956!-- Case of free convection (ws < 0.1 m/s ) not considered
     2957    ws = in_ws
     2958    IF ( ws < .1_wp ) THEN
     2959       ws = .1_wp
     2960    ENDIF
     2961
     2962!-- Heat_convection = forced convection
     2963    heat_convection = 12.1_wp * SQRT ( ws * pair / 1013.25_wp )
     2964
     2965!-- Average skin temperature
     2966    t_skin_aver = 35.7_wp - .0275_wp * activity
     2967
     2968!-- Calculation of constants for evaluation below
     2969    bc = .155_wp * clo * 3.96_wp * 10._wp**( -8._wp ) * f_cl
     2970    cc = f_cl * heat_convection
     2971    ec = .155_wp * clo
     2972    dc = ( 1._wp + ec * cc ) / bc
     2973    gc = ( t_skin_aver + bc * ( tmrt + 273.2_wp )**4._wp + ec * cc * ta ) / bc
     2974
     2975!-- Calculation of clothing surface temperature (t_clothing) based on
     2976!   newton-approximation with air temperature as initial guess
     2977    niter = dt
     2978    adjustrate = 1._wp - EXP ( -1._wp * ( 10._wp / time_equil ) * dt )
     2979    IF ( adjustrate >= 1._wp ) adjustrate = 1._wp
     2980    adjustrate_cloth = adjustrate * 30._wp
     2981    t_clothing = t_cloth
     2982
     2983    IF ( t_cloth <= -998.0_wp ) THEN  ! If initial run
     2984       niter = 3_wp
     2985       adjustrate = 1._wp
     2986       adjustrate_cloth = 1._wp
     2987       t_clothing = ta
     2988    ENDIF
     2989
     2990    DO i = 1, niter
     2991       t_clothing = t_clothing - adjustrate_cloth * ( ( t_clothing +           &
     2992          273.2_wp )**4._wp + t_clothing *                                     &
     2993          dc - gc ) / ( 4._wp * ( t_clothing + 273.2_wp )**3._wp + dc )
     2994    ENDDO
     2995
     2996!-- Empiric factor for the adaption of the heat ballance equation
     2997!   to the psycho-physical scale (Equ. 40 in FANGER)
     2998    z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp )
     2999
     3000!-- Water vapour diffution through the skin
     3001    z2 = .31_wp * ( 57.3_wp - .07_wp * activity-pa )
     3002
     3003!-- Sweat evaporation from the skin surface
     3004    z3 = .42_wp * ( activity - 58._wp )
     3005
     3006!-- Loss of latent heat through respiration
     3007    z4 = .0017_wp * actlev * ( 58.7_wp - pa ) + .0014_wp * actlev *            &
     3008      ( 34._wp - ta )
     3009
     3010!-- Loss of radiational heat
     3011    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + 273.2_wp )**4 - ( tmrt +         &
     3012       273.2_wp )**4 )
     3013
     3014!-- Heat loss through forced convection
     3015    z6 = cc * ( t_clothing - ta )
     3016
     3017!-- Write together as energy ballance
     3018    en = activity - z2 - z3 - z4 - z5 - z6
     3019
     3020!-- Manage storage
     3021    d_s = adjustrate * en + ( 1._wp - adjustrate ) * s
     3022
     3023!-- Predicted Mean Vote
     3024    pmva = z1 * d_s
     3025
     3026!-- Update storage
     3027    s = d_s
     3028    t_cloth = t_clothing
     3029
     3030 END SUBROUTINE fanger_s_acti
     3031
     3032
     3033
     3034!------------------------------------------------------------------------------!
     3035!
     3036! Description:
     3037! ------------
     3038!> Physiologically Equivalent Temperature (PET),
     3039!> stationary (calculated based on MEMI),
     3040!> Subroutine based on PETBER vers. 1.5.1996 by P. Hoeppe
     3041!------------------------------------------------------------------------------!
     3042
     3043 SUBROUTINE calculate_pet_static( ta, vpa, v, tmrt, pair, pet )
     3044
     3045    IMPLICIT NONE
     3046
     3047!-- Input arguments:
     3048    REAL(wp), INTENT( IN ) ::  ta    !< Air temperature             (°C)
     3049    REAL(wp), INTENT( IN ) ::  tmrt  !< Mean radiant temperature    (°C)
     3050    REAL(wp), INTENT( IN ) ::  v     !< Wind speed                  (m/s)
     3051    REAL(wp), INTENT( IN ) ::  vpa   !< Vapor pressure              (hPa)
     3052    REAL(wp), INTENT( IN ) ::  pair  !< Air pressure                (hPa)
     3053
     3054!-- Output arguments:
     3055    REAL(wp), INTENT ( OUT ) ::  pet  !< PET                         (°C)
     3056
     3057!-- Internal variables:
     3058    REAL(wp) ::  acl        !< clothing area                        (m²)
     3059    REAL(wp) ::  adu        !< Du Bois area                         (m²)
     3060    REAL(wp) ::  aeff       !< effective area                       (m²)
     3061    REAL(wp) ::  ere        !< energy ballance                      (W)
     3062    REAL(wp) ::  erel       !< latent energy ballance               (W)
     3063    REAL(wp) ::  esw        !< Energy-loss through sweat evap.      (W)
     3064    REAL(wp) ::  facl       !< Surface area extension through clothing (factor)
     3065    REAL(wp) ::  feff       !< Surface modification by posture      (factor)
     3066    REAL(wp) ::  rdcl       !< Diffusion resistence of clothing     (factor)
     3067    REAL(wp) ::  rdsk       !< Diffusion resistence of skin         (factor)
     3068    REAL(wp) ::  rtv
     3069    REAL(wp) ::  vpts       !< Sat. vapor pressure over skin        (hPa)
     3070    REAL(wp) ::  tsk        !< Skin temperature                     (°C)
     3071    REAL(wp) ::  tcl        !< Clothing temperature                 (°C)
     3072    REAL(wp) ::  wetsk      !< Fraction of wet skin                 (factor)
     3073
     3074!-- Variables:
     3075    REAL(wp) :: int_heat    !< Internal heat        (W)
     3076
     3077!-- MEMI configuration
     3078    REAL(wp) :: age         !< Persons age          (a)
     3079    REAL(wp) :: mbody       !< Persons body mass    (kg)
     3080    REAL(wp) :: ht          !< Persons height       (m)
     3081    REAL(wp) :: work        !< Work load            (W)
     3082    REAL(wp) :: eta         !< Work efficiency      (dimensionless)
     3083    REAL(wp) :: clo         !< Clothing insulation index (clo)
     3084    REAL(wp) :: fcl         !< Surface area modification by clothing (factor)
     3085!     INTEGER(iwp) :: pos     !< Posture: 1 = standing, 2 = sitting
     3086!     INTEGER(iwp) :: sex     !< Sex: 1 = male, 2 = female
     3087
     3088!-- Configuration, keep standard parameters!
     3089    age   = 35._wp
     3090    mbody = 75._wp
     3091    ht    =  1.75_wp
     3092    work  = 80._wp
     3093    eta   =  0._wp
     3094    clo   =  0.9_wp
     3095    fcl   =  1.15_wp
     3096
     3097!-- Call subfunctions
     3098    CALL in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta,    &
     3099            vpa, work )
     3100
     3101    CALL heat_exch ( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, ht, &
     3102            int_heat,mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa,      &
     3103            vpts, wetsk )
     3104
     3105    CALL pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair, rdcl,&
     3106            rdsk, rtv, ta, tcl, tsk, pet, vpts, wetsk )
     3107
     3108
     3109 END SUBROUTINE calculate_pet_static
     3110
     3111
     3112!------------------------------------------------------------------------------!
     3113! Description:
     3114! ------------
     3115!> Calculate internal energy ballance
     3116!------------------------------------------------------------------------------!
     3117 SUBROUTINE in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, &
     3118    vpa, work )
     3119
     3120!-- Input arguments:
     3121    REAL(wp), INTENT( IN )  ::  pair      !< air pressure             (hPa)
     3122    REAL(wp), INTENT( IN )  ::  ta        !< air temperature          (°C)
     3123    REAL(wp), INTENT( IN )  ::  vpa       !< vapor pressure           (hPa)
     3124    REAL(wp), INTENT( IN )  ::  age       !< Persons age              (a)
     3125    REAL(wp), INTENT( IN )  ::  mbody     !< Persons body mass        (kg)
     3126    REAL(wp), INTENT( IN )  ::  ht        !< Persons height           (m)
     3127    REAL(wp), INTENT( IN )  ::  work      !< Work load                (W)
     3128    REAL(wp), INTENT( IN )  ::  eta       !< Work efficiency      (dimensionless)
     3129
     3130!-- Output arguments:
     3131    REAL(wp), INTENT( OUT ) ::  ere       !< energy ballance          (W)
     3132    REAL(wp), INTENT( OUT ) ::  erel      !< latent energy ballance   (W)
     3133    REAL(wp), INTENT( OUT ) ::  int_heat  !< internal heat production (W)
     3134    REAL(wp), INTENT( OUT ) ::  rtv       !< respiratory volume
     3135
     3136!-- Constants:
     3137!     REAL(wp), PARAMETER :: cair = 1010._wp        !< replaced by c_p
     3138!     REAL(wp), PARAMETER :: evap = 2.42_wp * 10._wp **6._wp  !< replaced by l_v
     3139
     3140!-- Internal variables:
     3141    REAL(wp) ::  eres                     !< Sensible respiratory heat flux (W)
     3142    REAL(wp) ::  met
     3143    REAL(wp) ::  tex
     3144    REAL(wp) ::  vpex
     3145
     3146    met = 3.45_wp * mbody ** ( 3._wp / 4._wp ) * (1._wp + 0.004_wp *           &
     3147          ( 30._wp - age) + 0.010_wp * ( ( ht * 100._wp /                      &
     3148          ( mbody ** ( 1._wp / 3._wp ) ) ) - 43.4_wp ) )
     3149
     3150    met = work + met
     3151
     3152    int_heat = met * (1._wp - eta)
     3153
     3154!-- Sensible respiration energy
     3155    tex  = 0.47_wp * ta + 21.0_wp
     3156    rtv  = 1.44_wp * 10._wp ** (-6._wp) * met
     3157    eres = c_p * (ta - tex) * rtv
     3158
     3159!-- Latent respiration energy
     3160    vpex = 6.11_wp * 10._wp ** ( 7.45_wp * tex / ( 235._wp + tex ) )
     3161    erel = 0.623_wp * l_v / pair * ( vpa - vpex ) * rtv
     3162
     3163!-- Sum of the results
     3164    ere = eres + erel
     3165
     3166 END SUBROUTINE in_body
     3167
     3168
     3169!------------------------------------------------------------------------------!
     3170! Description:
     3171! ------------
     3172!> Calculate heat gain or loss
     3173!------------------------------------------------------------------------------!
     3174 SUBROUTINE heat_exch ( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff,  &
     3175        ht, int_heat, mbody, pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa,     &
     3176        vpts, wetsk )
     3177
     3178
     3179!-- Input arguments:
     3180    REAL(wp), INTENT( IN )  ::  ere    !< Energy ballance          (W)
     3181    REAL(wp), INTENT( IN )  ::  erel   !< Latent energy ballance   (W)
     3182    REAL(wp), INTENT( IN )  ::  int_heat  !< internal heat production (W)
     3183    REAL(wp), INTENT( IN )  ::  pair   !< Air pressure             (hPa)
     3184    REAL(wp), INTENT( IN )  ::  ta     !< Air temperature          (°C)
     3185    REAL(wp), INTENT( IN )  ::  tmrt   !< Mean radiant temperature (°C)
     3186    REAL(wp), INTENT( IN )  ::  v      !< Wind speed               (m/s)
     3187    REAL(wp), INTENT( IN )  ::  vpa    !< Vapor pressure           (hPa)
     3188    REAL(wp), INTENT( IN )  ::  mbody  !< body mass                (kg)
     3189    REAL(wp), INTENT( IN )  ::  ht     !< height                   (m)
     3190    REAL(wp), INTENT( IN )  ::  clo    !< clothing insulation      (clo)
     3191    REAL(wp), INTENT( IN )  ::  fcl    !< factor for surface area increase by clothing
     3192
     3193!-- Output arguments:
     3194    REAL(wp), INTENT( OUT ) ::  acl    !< Clothing surface area        (m²)
     3195    REAL(wp), INTENT( OUT ) ::  adu    !< Du-Bois area                 (m²)
     3196    REAL(wp), INTENT( OUT ) ::  aeff   !< Effective surface area       (m²)
     3197    REAL(wp), INTENT( OUT ) ::  esw    !< Energy-loss through sweat evap. (W)
     3198    REAL(wp), INTENT( OUT ) ::  facl   !< Surface area extension through clothing (factor)
     3199    REAL(wp), INTENT( OUT ) ::  feff   !< Surface modification by posture (factor)
     3200    REAL(wp), INTENT( OUT ) ::  rdcl   !< Diffusion resistence of clothing (factor)
     3201    REAL(wp), INTENT( OUT ) ::  rdsk   !< Diffusion resistence of skin (factor)
     3202    REAL(wp), INTENT( OUT ) ::  tcl    !< Clothing temperature         (°C)
     3203    REAL(wp), INTENT( OUT ) ::  tsk    !< Skin temperature             (°C)
     3204    REAL(wp), INTENT( OUT ) ::  vpts   !< Sat. vapor pressure over skin (hPa)
     3205    REAL(wp), INTENT( OUT ) ::  wetsk  !< Fraction of wet skin (dimensionless)
     3206
     3207!-- Cconstants:
     3208!     REAL(wp), PARAMETER :: cair = 1010._wp      !< replaced by c_p
     3209    REAL(wp), PARAMETER :: cb   = 3640._wp        !<
     3210    REAL(wp), PARAMETER :: emcl =    0.95_wp      !< Longwave emission coef. of cloth
     3211    REAL(wp), PARAMETER :: emsk =    0.99_wp      !< Longwave emission coef. of skin
     3212!     REAL(wp), PARAMETER :: evap = 2.42_wp * 10._wp **6._wp  !< replaced by l_v
     3213    REAL(wp), PARAMETER :: food =    0._wp        !< Heat gain by food        (W)
     3214    REAL(wp), PARAMETER :: po   = 1013.25_wp      !< Air pressure at sea level (hPa)
     3215    REAL(wp), PARAMETER :: rob  =    1.06_wp      !<
     3216
     3217!-- Internal variables
     3218    REAL(wp) ::  c(0:10)        !< Core temperature array           (°C)
     3219    REAL(wp) ::  cbare          !< Convection through bare skin
     3220    REAL(wp) ::  cclo           !< Convection through clothing
     3221    REAL(wp) ::  csum           !< Convection in total
     3222    REAL(wp) ::  di             !< difference between r1 and r2
     3223    REAL(wp) ::  ed             !< energy transfer by diffusion     (W)
     3224    REAL(wp) ::  enbal          !< energy ballance                  (W)
     3225    REAL(wp) ::  enbal2         !< energy ballance (storage, last cycle)
     3226    REAL(wp) ::  eswdif         !< difference between sweat production and evaporation potential
     3227    REAL(wp) ::  eswphy         !< sweat created by physiology
     3228    REAL(wp) ::  eswpot         !< potential sweat evaporation
     3229    REAL(wp) ::  fec            !<
     3230    REAL(wp) ::  hc             !<
     3231    REAL(wp) ::  he             !<
     3232    REAL(wp) ::  htcl           !<
     3233    REAL(wp) ::  r1             !<
     3234    REAL(wp) ::  r2             !<
     3235    REAL(wp) ::  rbare          !< Radiational loss of bare skin    (W/m²)
     3236    REAL(wp) ::  rcl            !<
     3237    REAL(wp) ::  rclo           !< Radiational loss of clothing     (W/m²)
     3238    REAL(wp) ::  rclo2          !< Longwave radiation gain or loss  (W/m²)
     3239    REAL(wp) ::  rsum           !< Radiational loss or gain         (W/m²)
     3240    REAL(wp) ::  sw             !<
     3241    REAL(wp) ::  swf            !<
     3242    REAL(wp) ::  swm            !<
     3243    REAL(wp) ::  tbody          !<
     3244    REAL(wp) ::  tcore(1:7)     !<
     3245    REAL(wp) ::  vb             !<
     3246    REAL(wp) ::  vb1            !<
     3247    REAL(wp) ::  vb2            !<
     3248    REAL(wp) ::  wd             !<
     3249    REAL(wp) ::  wr             !<
     3250    REAL(wp) ::  ws             !<
     3251    REAL(wp) ::  wsum           !<
     3252    REAL(wp) ::  xx             !< modification step                (K)
     3253    REAL(wp) ::  y              !< fraction of bare skin
     3254    INTEGER(iwp) ::  count1     !< running index
     3255    INTEGER(iwp) ::  count3     !< running index
     3256    INTEGER(iwp) ::  j          !< running index
     3257    INTEGER(iwp) ::  i          !< running index
     3258    LOGICAL ::  skipIncreaseCount   !< iteration control flag
     3259
     3260    wetsk = 0._wp
     3261    adu = 0.203_wp * mbody ** 0.425_wp * ht ** 0.725_wp
     3262
     3263    hc = 2.67_wp + ( 6.5_wp * v ** 0.67_wp)
     3264    hc   = hc * (pair /po) ** 0.55_wp
     3265    feff = 0.725_wp                     !< Posture: 0.725 for stading
     3266    facl = (- 2.36_wp + 173.51_wp * clo - 100.76_wp * clo * clo + 19.28_wp     &
     3267          * (clo ** 3._wp)) / 100._wp
     3268
     3269    IF ( facl > 1._wp )   facl = 1._wp
     3270    rcl = ( clo / 6.45_wp ) / facl
     3271    IF ( clo >= 2._wp )  y  = 1._wp
     3272
     3273    IF ( ( clo > 0.6_wp ) .AND. ( clo < 2._wp ) )  y = ( ht - 0.2_wp ) / ht
     3274    IF ( ( clo <= 0.6_wp ) .AND. ( clo > 0.3_wp ) ) y = 0.5_wp
     3275    IF ( ( clo <= 0.3_wp ) .AND. ( clo > 0._wp ) )  y = 0.1_wp
     3276
     3277    r2   = adu * (fcl - 1._wp + facl) / (2._wp * 3.14_wp * ht * y)
     3278    r1   = facl * adu / (2._wp * 3.14_wp * ht * y)
     3279
     3280    di   = r2 - r1
     3281
     3282!-- Skin temperatur
     3283    DO j = 1, 7
     3284
     3285       tsk    = 34._wp
     3286       count1 = 0_iwp
     3287       tcl    = ( ta + tmrt + tsk ) / 3._wp
     3288       count3 = 1_iwp
     3289       enbal2 = 0._wp
     3290
     3291       DO i = 1, 100  ! allow for 100 iterations max
     3292          acl   = adu * facl + adu * ( fcl - 1._wp )
     3293          rclo2 = emcl * sigma_sb * ( ( tcl + degc_to_k )**4._wp -             &
     3294            ( tmrt + degc_to_k )** 4._wp ) * feff
     3295          htcl  = 6.28_wp * ht * y * di / ( rcl * LOG( r2 / r1 ) * acl )
     3296          tsk   = 1._wp / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl
     3297
     3298!--       Radiation saldo
     3299          aeff  = adu * feff
     3300          rbare = aeff * ( 1._wp - facl ) * emsk * sigma_sb *                  &
     3301            ( ( tmrt + degc_to_k )** 4._wp - ( tsk + degc_to_k )** 4._wp )
     3302          rclo  = feff * acl * emcl * sigma_sb *                               &
     3303            ( ( tmrt + degc_to_k )** 4._wp - ( tcl + degc_to_k )** 4._wp )
     3304          rsum  = rbare + rclo
     3305
     3306!--       Convection
     3307          cbare = hc * ( ta - tsk ) * adu * ( 1._wp - facl )
     3308          cclo  = hc * ( ta - tcl ) * acl
     3309          csum  = cbare + cclo
     3310
     3311!--       Core temperature
     3312          c(0)  = int_heat + ere
     3313          c(1)  = adu * rob * cb
     3314          c(2)  = 18._wp - 0.5_wp * tsk
     3315          c(3)  = 5.28_wp * adu * c(2)
     3316          c(4)  = 0.0208_wp * c(1)
     3317          c(5)  = 0.76075_wp * c(1)
     3318          c(6)  = c(3) - c(5) - tsk * c(4)
     3319          c(7)  = - c(0) * c(2) - tsk * c(3) + tsk * c(5)
     3320          c(8)  = c(6) * c(6) - 4._wp * c(4) * c(7)
     3321          c(9)  = 5.28_wp * adu - c(5) - c(4) * tsk
     3322          c(10) = c(9) * c(9) - 4._wp * c(4) *                                 &
     3323             ( c(5) * tsk - c(0) - 5.28_wp * adu * tsk )
     3324
     3325          IF ( tsk == 36._wp ) tsk = 36.01_wp
     3326          tcore(7) = c(0) / ( 5.28_wp * adu + c(1) * 6.3_wp / 3600._wp ) + tsk
     3327          tcore(3) = c(0) / ( 5.28_wp * adu + ( c(1) * 6.3_wp / 3600._wp ) /   &
     3328            ( 1._wp + 0.5_wp * ( 34._wp - tsk ) ) ) + tsk
     3329          IF ( c(10) >= 0._wp ) THEN
     3330             tcore(6) = ( - c(9) - c(10)**0.5_wp ) / ( 2._wp * c(4) )
     3331             tcore(1) = ( - c(9) + c(10)**0.5_wp ) / ( 2._wp * c(4) )
     3332          END IF
     3333
     3334          IF ( c(8) >= 0._wp ) THEN
     3335             tcore(2) = ( - c(6) + ABS( c(8) ) ** 0.5_wp ) / ( 2._wp * c(4) )
     3336             tcore(5) = ( - c(6) - ABS( c(8) ) ** 0.5_wp ) / ( 2._wp * c(4) )
     3337             tcore(4) = c(0) / ( 5.28_wp * adu + c(1) * 1._wp / 40._wp ) + tsk
     3338          END IF
     3339
     3340!--       Transpiration
     3341          tbody = 0.1_wp * tsk + 0.9_wp * tcore(j)
     3342          swm   = 304.94_wp * ( tbody - 36.6_wp ) * adu / 3600000._wp
     3343          vpts  = 6.11_wp * 10._wp**( 7.45_wp * tsk / ( 235._wp + tsk ) )
     3344
     3345          IF ( tbody <= 36.6_wp ) swm = 0._wp  !< no need for sweating
     3346
     3347          sw = swm
     3348          eswphy = - sw * l_v
     3349          he     = 0.633_wp * hc / ( pair * c_p )
     3350          fec    = 1._wp / ( 1._wp + 0.92_wp * hc * rcl )
     3351          eswpot = he * ( vpa - vpts ) * adu * l_v * fec
     3352          wetsk  = eswphy / eswpot
     3353
     3354          IF ( wetsk > 1._wp ) wetsk = 1._wp
     3355!
     3356!--       Sweat production > evaporation?
     3357          eswdif = eswphy - eswpot
     3358
     3359          IF ( eswdif <= 0._wp ) esw = eswpot  !< Limit is evaporation
     3360          IF ( eswdif > 0._wp ) esw = eswphy   !< Limit is sweat production
     3361          IF ( esw  > 0._wp )   esw = 0._wp    !< Sweat can't be evaporated, no more cooling effect
     3362
     3363!--       Diffusion
     3364          rdsk = 0.79_wp * 10._wp ** 7._wp
     3365          rdcl = 0._wp
     3366          ed   = l_v / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( vpa -     &
     3367             vpts )
     3368
     3369!--       Max vb
     3370          vb1 = 34._wp - tsk
     3371          vb2 = tcore(j) - 36.6_wp
     3372
     3373          IF ( vb2 < 0._wp ) vb2 = 0._wp
     3374          IF ( vb1 < 0._wp ) vb1 = 0._wp
     3375          vb = ( 6.3_wp + 75._wp * vb2 ) / ( 1._wp + 0.5_wp * vb1 )
     3376
     3377!--       Energy ballence
     3378          enbal = int_heat + ed + ere + esw + csum + rsum + food
     3379
     3380!--       Clothing temperature
     3381          xx = 0.001_wp
     3382          IF ( count1 == 0_iwp ) xx = 1._wp
     3383          IF ( count1 == 1_iwp ) xx = 0.1_wp
     3384          IF ( count1 == 2_iwp ) xx = 0.01_wp
     3385          IF ( count1 == 3_iwp ) xx = 0.001_wp
     3386
     3387          IF ( enbal > 0._wp ) tcl = tcl + xx
     3388          IF ( enbal < 0._wp ) tcl = tcl - xx
     3389
     3390          skipIncreaseCount = .FALSE.
     3391          IF ( ( (enbal <= 0._wp ) .AND. (enbal2 > 0._wp ) ) .OR.              &
     3392             ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) ) THEN
     3393             skipIncreaseCount = .TRUE.
     3394          ELSE
     3395             enbal2 = enbal
     3396             count3 = count3 + 1_iwp
     3397          END IF
     3398
     3399          IF ( ( count3 > 200_iwp ) .OR. skipIncreaseCount ) THEN
     3400             IF ( count1 < 3_iwp ) THEN
     3401                count1 = count1 + 1_iwp
     3402                enbal2 = 0._wp
     3403             ELSE
     3404                EXIT
     3405             END IF
     3406          END IF
     3407       END DO
     3408
     3409       IF ( count1 == 3_iwp ) THEN
     3410          SELECT CASE ( j )
     3411             CASE ( 2, 5)
     3412                IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND.                     &
     3413                   ( tsk <= 34.050_wp ) ) ) CYCLE
     3414             CASE ( 6, 1 )
     3415                IF ( c(10) < 0._wp ) CYCLE
     3416                IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND.                     &
     3417                   ( tsk > 33.850_wp ) ) ) CYCLE
     3418             CASE ( 3 )
     3419                IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND.                      &
     3420                   ( tsk <= 34.000_wp ) ) ) CYCLE
     3421             CASE ( 7 )
     3422                IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND.                      &
     3423                   ( tsk > 34.000_wp ) ) ) CYCLE
     3424             CASE default
     3425          END SELECT
     3426       END IF
     3427
     3428       IF ( ( j /= 4_iwp ) .AND. ( vb >= 91._wp ) ) CYCLE
     3429       IF ( ( j == 4_iwp ) .AND. ( vb < 89._wp ) ) CYCLE
     3430       IF ( vb > 90._wp ) vb = 90._wp
     3431
     3432!--    Loses by water
     3433       ws = sw * 3600._wp * 1000._wp
     3434       IF ( ws > 2000._wp ) ws = 2000._wp
     3435       wd = ed / l_v * 3600._wp * ( -1000._wp )
     3436       wr = erel / l_v * 3600._wp * ( -1000._wp )
     3437
     3438       wsum = ws + wr + wd
     3439
     3440       RETURN
     3441    END DO
     3442 END SUBROUTINE heat_exch
     3443
     3444!------------------------------------------------------------------------------!
     3445! Description:
     3446! ------------
     3447!> Calculate PET
     3448!------------------------------------------------------------------------------!
     3449 SUBROUTINE pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair,   &
     3450        rdcl, rdsk, rtv, ta, tcl, tsk, pet, vpts, wetsk )
     3451
     3452!-- Input arguments:
     3453    REAL(wp), INTENT( IN ) ::  acl   !< clothing surface area        (m²)
     3454    REAL(wp), INTENT( IN ) ::  adu   !< Du-Bois area                 (m²)
     3455    REAL(wp), INTENT( IN ) ::  esw   !< energy-loss through sweat evap. (W)
     3456    REAL(wp), INTENT( IN ) ::  facl  !< surface area extension through clothing (factor)
     3457    REAL(wp), INTENT( IN ) ::  feff  !< surface modification by posture (factor)
     3458    REAL(wp), INTENT( IN ) ::  int_heat  !< internal heat production (W)
     3459    REAL(wp), INTENT( IN ) ::  pair  !< air pressure                 (hPa)
     3460    REAL(wp), INTENT( IN ) ::  rdcl  !< diffusion resistence of clothing (factor)
     3461    REAL(wp), INTENT( IN ) ::  rdsk  !< diffusion resistence of skin (factor)
     3462    REAL(wp), INTENT( IN ) ::  rtv   !< respiratory volume
     3463    REAL(wp), INTENT( IN ) ::  ta    !< air temperature              (°C)
     3464    REAL(wp), INTENT( IN ) ::  tcl   !< clothing temperature         (°C)
     3465    REAL(wp), INTENT( IN ) ::  tsk   !< skin temperature             (°C)
     3466    REAL(wp), INTENT( IN ) ::  vpts  !< sat. vapor pressure over skin (hPa)
     3467    REAL(wp), INTENT( IN ) ::  wetsk !< fraction of wet skin (dimensionless)
     3468
     3469!-- Output arguments:
     3470    REAL(wp), INTENT( OUT ) ::  aeff  !< effective surface area       (m²)
     3471    REAL(wp), INTENT( OUT ) ::  pet   !< PET                          (°C)
     3472
     3473!-- Cconstants:
     3474!     REAL(wp), PARAMETER :: cair = 1010._wp        !< replaced by c_p
     3475    REAL(wp), PARAMETER :: emcl =    0.95_wp      !< Longwave emission coef. of cloth
     3476    REAL(wp), PARAMETER :: emsk =    0.99_wp      !< Longwave emission coef. of skin
     3477!     REAL(wp), PARAMETER :: evap = 2.42_wp * 10._wp **6._wp  !< replaced by l_v
     3478    REAL(wp), PARAMETER :: po   = 1013.25_wp      !< Air pressure at sea level (hPa)
     3479
     3480!-- Internal variables
     3481    REAL ( wp ) ::  cbare             !< Convection through bare skin
     3482    REAL ( wp ) ::  cclo              !< Convection through clothing
     3483    REAL ( wp ) ::  csum              !< Convection in total
     3484    REAL ( wp ) ::  ed                !< Diffusion                      (W)
     3485    REAL ( wp ) ::  enbal             !< Energy ballance                (W)
     3486    REAL ( wp ) ::  enbal2            !< Energy ballance (last iteration cycle)
     3487    REAL ( wp ) ::  ere               !< Energy ballance result         (W)
     3488    REAL ( wp ) ::  erel              !< Latent energy ballance         (W)
     3489    REAL ( wp ) ::  eres              !< Sensible respiratory heat flux (W)
     3490    REAL ( wp ) ::  hc                !<
     3491    REAL ( wp ) ::  rbare             !< Radiational loss of bare skin  (W/m²)
     3492    REAL ( wp ) ::  rclo              !< Radiational loss of clothing   (W/m²)
     3493    REAL ( wp ) ::  rsum              !< Radiational loss or gain       (W/m²)
     3494    REAL ( wp ) ::  tex               !< Temperat. of exhaled air       (°C)
     3495    REAL ( wp ) ::  vpex              !< Vapor pressure of exhaled air  (hPa)
     3496    REAL ( wp ) ::  xx                !< Delta PET per iteration        (K)
     3497
     3498    INTEGER ( iwp ) ::  count1        !< running index
     3499    INTEGER ( iwp ) ::  i             !< running index
     3500
     3501    pet = ta
     3502    enbal2 = 0._wp
     3503
     3504    DO count1 = 0, 3
     3505       DO i = 1, 125  ! 500 / 4
     3506          hc = 2.67_wp + 6.5_wp * 0.1_wp ** 0.67_wp
     3507          hc = hc * ( pair / po ) ** 0.55_wp
     3508
     3509!--       Radiation
     3510          aeff  = adu * feff
     3511          rbare = aeff * ( 1._wp - facl ) * emsk * sigma_sb *                  &
     3512              ( ( pet + degc_to_k ) ** 4._wp - ( tsk + degc_to_k ) ** 4._wp )
     3513          rclo  = feff * acl * emcl * sigma_sb *                               &
     3514              ( ( pet + degc_to_k ) ** 4._wp - ( tcl + degc_to_k ) ** 4._wp )
     3515          rsum  = rbare + rclo
     3516
     3517!--       Covection
     3518          cbare = hc * ( pet - tsk ) * adu * ( 1._wp - facl )
     3519          cclo  = hc * ( pet - tcl ) * acl
     3520          csum  = cbare + cclo
     3521
     3522!--       Diffusion
     3523          ed = l_v / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( 12._wp -    &
     3524             vpts )
     3525
     3526!--       Respiration
     3527          tex  = 0.47_wp * pet + 21._wp
     3528          eres = c_p * ( pet - tex ) * rtv
     3529          vpex = 6.11_wp * 10._wp ** ( 7.45_wp * tex / ( 235._wp + tex ) )
     3530          erel = 0.623_wp * l_v / pair * ( 12._wp - vpex ) * rtv
     3531          ere  = eres + erel
     3532
     3533!--       Energy ballance
     3534          enbal = int_heat + ed + ere + esw + csum + rsum
     3535
     3536!--       Iteration concerning ta
     3537          IF ( count1 == 0_iwp )  xx = 1._wp
     3538          IF ( count1 == 1_iwp )  xx = 0.1_wp
     3539          IF ( count1 == 2_iwp )  xx = 0.01_wp
     3540          IF ( count1 == 3_iwp )  xx = 0.001_wp
     3541          IF ( enbal > 0._wp )  pet = pet - xx
     3542          IF ( enbal < 0._wp )  pet = pet + xx
     3543          IF ( ( enbal <= 0._wp ) .AND. ( enbal2 > 0._wp ) ) EXIT
     3544          IF ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) EXIT
     3545
     3546          enbal2 = enbal
     3547       END DO
     3548    END DO
     3549 END SUBROUTINE pet_iteration
     3550
    14393551
    14403552 END MODULE biometeorology_mod
Note: See TracChangeset for help on using the changeset viewer.