Ignore:
Timestamp:
Nov 27, 2018 5:03:40 PM (3 years ago)
Author:
kanani
Message:

Fix for biomet output (ticket:757), merge of uv_exposure into biometeorology_mod

File:
1 edited

Legend:

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

    r3525 r3569  
    2121!
    2222! Current revisions:
    23 ! -----------------
     23! ------------------
    2424!
    2525!
     
    2727! -----------------
    2828! $Id$
     29! Consistently use bio_fill_value everywhere,
     30! move allocation and initialization of output variables to bio_check_data_output
     31! and bio_3d_data_averaging,
     32! dom_dwd_user, Schrempf:
     33! - integration of UVEM module part from r3474 (edited)
     34! - split UTCI regression into 6 parts
     35! - all data_output_3d is now explicity casted to sp
     36!
     37! 3525 2018-11-14 16:06:14Z kanani
    2938! Clean up, renaming from "biom" to "bio", summary of thermal index calculation
    3039! into one module (dom_dwd_user)
     
    4958! Authors:
    5059! --------
    51 ! @author Dominik Froehlich <dominik.froehlich@dwd.de>
    52 ! @author Jaroslav Resler <resler@cs.cas.cz>
     60! @author Dominik Froehlich <dominik.froehlich@dwd.de>, thermal indices
     61! @author Jaroslav Resler <resler@cs.cas.cz>, mean radiant temperature
     62! @author Michael Schrempf <schrempf@muk.uni-hannover.de>, uv exposure
    5363!
    5464!
    5565! Description:
    5666! ------------
    57 !> Human thermal comfort module calculating thermal perception of a sample
     67!> Biometeorology module consisting of two parts:
     68!> 1.: Human thermal comfort module calculating thermal perception of a sample
    5869!> human being under the current meteorological conditions.
     70!> 2.: Calculation of vitamin-D weighted UV exposure
    5971!>
    6072!> @todo Alphabetical sorting of "USE ..." lists, "ONLY" list, variable declarations
    6173!>       (per subroutine: first all CHARACTERs, then INTEGERs, LOGICALs, REALs, )
    6274!> @todo Comments start with capital letter --> "!-- Include..."
    63 !> @todo Variable and routine names strictly in lowercase letters and in English
     75!> @todo uv_vitd3dose-->new output type necessary (cumulative)
     76!> @todo consider upwelling radiation in UV
    6477!>
    6578!> @note nothing now
     
    7689
    7790    USE basic_constants_and_equations_mod,                                     &
    78        ONLY:  c_p, degc_to_k, l_v, magnus, sigma_sb
     91       ONLY:  c_p, degc_to_k, l_v, magnus, sigma_sb, pi
    7992
    8093    USE control_parameters,                                                    &
     
    8396              simulated_time, surface_pressure
    8497
     98    USE date_and_time_mod,                                                                                                        &
     99        ONLY:  calc_date_and_time, day_of_year, time_utc
     100
    85101    USE grid_variables,                                                        &
    86102       ONLY:  ddx, dx, ddy, dy
     
    92108    USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM
    93109
     110    USE netcdf_data_input_mod,                                                                                                    &
     111        ONLY:  netcdf_data_input_uvem, uvem_projarea_f, uvem_radiance_f,                                                          &
     112               uvem_irradiance_f, uvem_integration_f, building_obstruction_f
     113!
    94114!-- Import radiation model to obtain input for mean radiant temperature
    95115    USE radiation_model_mod,                                                   &
     
    106126    PRIVATE
    107127
     128!
    108129!-- Declare all global variables within the module (alphabetical order)
    109130    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_grid  !< tmrt results (°C)
    110131    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct      !< PT results   (°C)
    111     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_grid  !< UTCI results (°C)
    112     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_grid   !< PET results  (°C)
     132    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci       !< UTCI results (°C)
     133    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet        !< PET results  (°C)
     134!
    113135!-- Grids for averaged thermal indices
    114136    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  mrt_av_grid   !< time average mean
    115137    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct_av      !< PT results (aver. input)   (°C)
    116     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_av_grid  !< UTCI results (aver. input) (°C)
    117     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_av_grid   !< PET results (aver. input)  (°C)
     138    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_av       !< UTCI results (aver. input) (°C)
     139    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_av        !< PET results (aver. input)  (°C)
    118140
    119141
     
    123145    REAL ( wp ), PARAMETER ::  human_absorb = 0.7_wp  !< SW absorbtivity of a human body (Fanger 1972)
    124146    REAL ( wp ), PARAMETER ::  human_emiss = 0.97_wp  !< LW emissivity of a human body after (Fanger 1972)
     147    REAL ( wp ), PARAMETER ::  bio_fill_value = -9999._wp  !< set module fill value, replace by global fill value as soon as available
     148!
    125149!--
    126 
    127150    LOGICAL ::  aver_perct = .FALSE.  !< switch: do perct averaging in this module? (if .FALSE. this is done globally)
    128151    LOGICAL ::  aver_q     = .FALSE.  !< switch: do e  averaging in this module?
     
    135158    LOGICAL ::  average_trigger_pet   = .FALSE.  !< update averaged input on call to bio_pet?
    136159
     160    LOGICAL ::  thermal_comfort = .TRUE.  !< Turn all thermal indices on or off
    137161    LOGICAL ::  bio_perct     = .TRUE.   !< Turn index PT (instant. input) on or off
    138162    LOGICAL ::  bio_perct_av  = .TRUE.   !< Turn index PT (averaged input) on or off
     
    142166    LOGICAL ::  bio_utci_av   = .TRUE.   !< Turn index UTCI (averaged input) on or off
    143167
     168!
     169!-- UVEM parameters from here
     170!
     171!-- Declare all global variables within the module (alphabetical order)
     172    INTEGER(iwp) ::  ai                      = 0  !< loop index in azimuth direction
     173    INTEGER(iwp) ::  bi                      = 0  !< loop index of bit location within an 8bit-integer (one Byte)
     174    INTEGER(iwp) ::  clothing                = 1  !< clothing (0=unclothed, 1=Arms,Hands,Face free, 3=Hand,Face free)
     175    INTEGER(iwp) ::  iq                      = 0  !< loop index of irradiance quantity
     176    INTEGER(iwp) ::  pobi                    = 0  !< loop index of the position of corresponding byte within ibset byte vektor
     177    INTEGER(iwp) ::  obstruction_direct_beam = 0  !< Obstruction information for direct beam   
     178    INTEGER(iwp) ::  zi                      = 0  !< loop index in zenith direction
     179
     180    INTEGER(KIND=1), DIMENSION(0:44)  ::  obstruction_temp1 = 0  !< temporary obstruction information stored with ibset
     181    INTEGER(iwp),    DIMENSION(0:359) ::  obstruction_temp2 = 0  !< restored temporary obstruction information from ibset file
     182
     183    INTEGER(iwp), DIMENSION(0:35,0:9) ::  obstruction       = 1  !< final 2D obstruction information array
     184
     185    LOGICAL ::  consider_obstructions = .TRUE.   !< namelist parameter (see documentation)
     186    LOGICAL ::  sun_in_south          = .FALSE.  !< namelist parameter (see documentation)
     187    LOGICAL ::  turn_to_sun           = .TRUE.   !< namelist parameter (see documentation)
     188    LOGICAL ::  uv_exposure           = .FALSE.  !< namelist parameter (see documentation)
     189
     190    REAL(wp) ::  diffuse_exposure            =   0.0_wp  !< calculated exposure by diffuse radiation
     191    REAL(wp) ::  direct_exposure             =   0.0_wp  !< calculated exposure by direct solar beam   
     192    REAL(wp) ::  orientation_angle           =   0.0_wp  !< orientation of front/face of the human model   
     193    REAL(wp) ::  projection_area_direct_beam =   0.0_wp  !< projection area for direct solar beam
     194    REAL(wp) ::  saa                         = 180.0_wp  !< solar azimuth angle
     195    REAL(wp) ::  startpos_human              =   0.0_wp  !< start value for azimuth interpolation of human geometry array
     196    REAL(wp) ::  startpos_saa_float          =   0.0_wp  !< start value for azimuth interpolation of radiance array
     197    REAL(wp) ::  sza                         =  20.0_wp  !< solar zenith angle
     198    REAL(wp) ::  xfactor                     =   0.0_wp  !< relative x-position used for interpolation
     199    REAL(wp) ::  yfactor                     =   0.0_wp  !< relative y-position used for interpolation
     200
     201    REAL(wp), DIMENSION(0:2)  ::  irradiance =   0.0_wp  !< iradiance values extracted from irradiance lookup table 
     202
     203    REAL(wp), DIMENSION(0:2,0:90) ::  irradiance_lookup_table      = 0.0_wp  !< irradiance lookup table
     204    REAL(wp), DIMENSION(0:35,0:9) ::  integration_array            = 0.0_wp  !< solid angle factors for hemispherical integration
     205    REAL(wp), DIMENSION(0:35,0:9) ::  projection_area              = 0.0_wp  !< projection areas of a human (all directions)
     206    REAL(wp), DIMENSION(0:35,0:9) ::  projection_area_lookup_table = 0.0_wp  !< human geometry lookup table (projection areas)
     207    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_direct_temp  = 0.0_wp  !< temporary projection area for direct solar beam
     208    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_temp         = 0.0_wp  !< temporary projection area for all directions
     209    REAL(wp), DIMENSION(0:35,0:9) ::  radiance_array               = 0.0_wp  !< radiance extracted from radiance_lookup_table 
     210    REAL(wp), DIMENSION(0:71,0:9) ::  radiance_array_temp          = 0.0_wp  !< temporary radiance data
     211
     212    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vitd3_exposure     !< result variable for instantaneous vitamin-D weighted exposures
     213    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vitd3_exposure_av  !< result variable for summation of vitamin-D weighted exposures
     214
     215    REAL(wp), DIMENSION(0:35,0:9,0:90) ::  radiance_lookup_table   = 0.0_wp  !< radiance lookup table
    144216
    145217!
     
    150222    bio_check_parameters, bio_data_output_3d, bio_data_output_2d,              &
    151223    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
     224    bio_init, bio_parin, bio_perct, bio_perct_av, bio_pet,                     &
     225    bio_pet_av, bio_utci, bio_utci_av, thermal_comfort, time_bio_results,      &
     226!
     227!-- UVEM PUBLIC variables and methods
     228    uvem_calc_exposure, uv_exposure
    154229
    155230!
     
    160235       MODULE PROCEDURE bio_3d_data_averaging
    161236    END INTERFACE bio_3d_data_averaging
    162 
     237!
    163238!-- Calculate mtr from rtm fluxes and assign into 2D grid
    164239    INTERFACE bio_calculate_mrt_grid
    165240       MODULE PROCEDURE bio_calculate_mrt_grid
    166241    END INTERFACE bio_calculate_mrt_grid
    167 
     242!
    168243!-- Calculate static thermal indices PT, UTCI and/or PET
    169244    INTERFACE bio_calculate_thermal_index_maps
    170245       MODULE PROCEDURE bio_calculate_thermal_index_maps
    171246    END INTERFACE bio_calculate_thermal_index_maps
    172 
     247!
    173248!-- Calculate the dynamic index iPT (to be caled by the agent model)
    174249    INTERFACE bio_calc_ipt
    175250       MODULE PROCEDURE bio_calc_ipt
    176251    END INTERFACE bio_calc_ipt
    177 
     252!
    178253!-- Data output checks for 2D/3D data to be done in check_parameters
    179      INTERFACE bio_check_data_output
    180         MODULE PROCEDURE bio_check_data_output
    181      END INTERFACE bio_check_data_output
    182 
     254    INTERFACE bio_check_data_output
     255       MODULE PROCEDURE bio_check_data_output
     256    END INTERFACE bio_check_data_output
     257!
    183258!-- Input parameter checks to be done in check_parameters
    184259    INTERFACE bio_check_parameters
    185260       MODULE PROCEDURE bio_check_parameters
    186261    END INTERFACE bio_check_parameters
    187 
     262!
    188263!-- Data output of 2D quantities
    189264    INTERFACE bio_data_output_2d
    190265       MODULE PROCEDURE bio_data_output_2d
    191266    END INTERFACE bio_data_output_2d
    192 
     267!
    193268!-- no 3D data, thus, no averaging of 3D data, removed
    194269    INTERFACE bio_data_output_3d
    195270       MODULE PROCEDURE bio_data_output_3d
    196271    END INTERFACE bio_data_output_3d
    197 
     272!
    198273!-- Definition of data output quantities
    199274    INTERFACE bio_define_netcdf_grid
    200275       MODULE PROCEDURE bio_define_netcdf_grid
    201276    END INTERFACE bio_define_netcdf_grid
    202 
     277!
    203278!-- Obtains all relevant input values to estimate local thermal comfort/stress
    204279    INTERFACE bio_get_thermal_index_input_ij
    205280       MODULE PROCEDURE bio_get_thermal_index_input_ij
    206281    END INTERFACE bio_get_thermal_index_input_ij
    207 
     282!
    208283!-- Output of information to the header file
    209284    INTERFACE bio_header
    210285       MODULE PROCEDURE bio_header
    211286    END INTERFACE bio_header
    212 
     287!
    213288!-- Initialization actions
    214289    INTERFACE bio_init
    215290       MODULE PROCEDURE bio_init
    216291    END INTERFACE bio_init
    217 
    218 !-- Initialization of arrays
    219     INTERFACE bio_init_arrays
    220        MODULE PROCEDURE bio_init_arrays
    221     END INTERFACE bio_init_arrays
    222 
     292!
    223293!-- Reading of NAMELIST parameters
    224294    INTERFACE bio_parin
     
    226296    END INTERFACE bio_parin
    227297
     298!
     299!-- Calculate UV exposure grid
     300    INTERFACE uvem_calc_exposure
     301       MODULE PROCEDURE uvem_calc_exposure
     302    END INTERFACE uvem_calc_exposure
    228303
    229304 CONTAINS
     
    235310!> Sum up and time-average biom input quantities as well as allocate
    236311!> the array necessary for storing the average.
     312!> There is a considerable difference to the 3d_data_averaging subroutines
     313!> used by other modules:
     314!> For the thermal indices, the module needs to average the input conditions
     315!> not the result!
    237316!------------------------------------------------------------------------------!
    238317 SUBROUTINE bio_3d_data_averaging( mode, variable )
     
    260339          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
    261340
    262 !--          Indices in unknown order as depending on input file, determine
    263 !            first index to average und update only once
    264              IF ( .NOT. average_trigger_perct .AND. .NOT. average_trigger_utci &
    265                 .AND. .NOT. average_trigger_pet ) THEN
     341!
     342!--          Averaging, as well as the allocation of the required grids must be
     343!            done only once, independent from for how many thermal indices
     344!            averaged output is desired.
     345!            Therefore wee need to memorize which index is the one that controls
     346!            the averaging (what must be the first thermal index called).
     347!            Indices are in unknown order as depending on the input file,
     348!            determine first index to average und update only once
     349
     350!--          Only proceed here if this was not done for any index before. This
     351!            is done only once during the whole model run.
     352             IF ( .NOT. average_trigger_perct .AND.                            &
     353                  .NOT. average_trigger_utci  .AND.                            &
     354                  .NOT. average_trigger_pet ) THEN
     355!
     356!--             Allocate the required grids
     357                IF ( .NOT. ALLOCATED( perct_av ) ) THEN
     358                   ALLOCATE( perct_av (nys:nyn,nxl:nxr) )
     359                ENDIF
     360                perct_av = REAL( bio_fill_value, KIND = wp )
     361
     362                IF ( .NOT. ALLOCATED( utci_av ) ) THEN
     363                   ALLOCATE( utci_av (nys:nyn,nxl:nxr) )
     364                ENDIF
     365                utci_av = REAL( bio_fill_value, KIND = wp )
     366
     367                IF ( .NOT. ALLOCATED( pet_av ) ) THEN
     368                   ALLOCATE( pet_av (nys:nyn,nxl:nxr) )
     369                ENDIF
     370                pet_av = REAL( bio_fill_value, KIND = wp )
     371!
     372!--             Memorize the first index called to control averaging
    266373                IF ( TRIM( variable ) == 'bio_perct*' ) THEN
    267374                    average_trigger_perct = .TRUE.
     
    274381                ENDIF
    275382             ENDIF
    276 
    277 !--          Only continue if updateindex
     383!
     384!--          Only continue if var is the index, that controls averaging.
     385!            Break immediatelly (doing nothing) for the other indices.
    278386             IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') RETURN
    279387             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*')   RETURN
    280388             IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'bio_pet*')    RETURN
    281 
    282 !--          Set averaging switch to .TRUE. if not done by other module before!
     389!
     390!--          Now memorize which of the input grids are not averaged by other
     391!            modules. Set averaging switch to .TRUE. in that case.
    283392             IF ( .NOT. ALLOCATED( pt_av ) )  THEN
    284393                ALLOCATE( pt_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
    285394                aver_perct = .TRUE.
     395                pt_av = 0.0_wp
    286396             ENDIF
    287              pt_av = 0.0_wp
    288397
    289398             IF ( .NOT. ALLOCATED( q_av ) )  THEN
    290399                ALLOCATE( q_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
    291400                aver_q = .TRUE.
     401                q_av = 0.0_wp
    292402             ENDIF
    293              q_av = 0.0_wp
    294403
    295404             IF ( .NOT. ALLOCATED( u_av ) )  THEN
    296405                ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    297406                aver_u = .TRUE.
     407                u_av = 0.0_wp
    298408             ENDIF
    299              u_av = 0.0_wp
    300409
    301410             IF ( .NOT. ALLOCATED( v_av ) )  THEN
    302411                ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    303412                aver_v = .TRUE.
     413                v_av = 0.0_wp
    304414             ENDIF
    305              v_av = 0.0_wp
    306415
    307416             IF ( .NOT. ALLOCATED( w_av ) )  THEN
    308417                ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    309418                aver_w = .TRUE.
     419                w_av = 0.0_wp
    310420             ENDIF
    311              w_av = 0.0_wp
     421
     422          CASE ( 'uvem_vitd3dose*' )
     423             IF ( .NOT. ALLOCATED( vitd3_exposure_av ) )  THEN
     424                ALLOCATE( vitd3_exposure_av(nysg:nyng,nxlg:nxrg) )
     425             ENDIF
     426             vitd3_exposure_av = 0.0_wp
    312427
    313428          CASE DEFAULT
     
    335450
    336451          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
    337 
    338 !--          Only continue if updateindex
     452!
     453!--          Only continue if updateindex, see above
    339454             IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') &
    340455                RETURN
     
    393508                ENDDO
    394509             ENDIF
    395 
    396            CASE DEFAULT
     510!
     511!--       This is a cumulated dose. No mode == 'average' for this quantity.
     512          CASE ( 'uvem_vitd3dose*' )
     513             IF ( ALLOCATED( vitd3_exposure_av ) ) THEN
     514                DO  i = nxlg, nxrg
     515                   DO  j = nysg, nyng
     516                      vitd3_exposure_av(j,i) = vitd3_exposure_av(j,i) + vitd3_exposure(j,i)
     517                   ENDDO
     518                ENDDO
     519             ENDIF
     520
     521          CASE DEFAULT
    397522             CONTINUE
    398523
     
    409534
    410535          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
    411 
    412 !--          Only continue if update index
    413              IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') &
    414                 RETURN
    415              IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*')  &
    416                 RETURN
    417              IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'bio_pet*')   &
    418                 RETURN
     536!
     537!--          Only continue if update index, see above
     538             IF ( average_trigger_perct .AND.                                 &
     539                TRIM( variable ) /= 'bio_perct*' ) RETURN
     540             IF ( average_trigger_utci .AND.                                   &
     541                TRIM( variable ) /= 'bio_utci*' ) RETURN
     542             IF ( average_trigger_pet  .AND.                                   &
     543                TRIM( variable ) /= 'bio_pet*' ) RETURN
    419544
    420545             IF ( ALLOCATED( pt_av ) .AND. aver_perct ) THEN
     
    422547                   DO  j = nys, nyn
    423548                      DO  k = nzb, nzt+1
    424                          pt_av(k,j,i) = pt_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     549                         pt_av(k,j,i) = pt_av(k,j,i) /                         &
     550                            REAL( average_count_3d, KIND=wp )
    425551                      ENDDO
    426552                   ENDDO
     
    432558                   DO  j = nys, nyn
    433559                      DO  k = nzb, nzt+1
    434                          q_av(k,j,i) = q_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     560                         q_av(k,j,i) = q_av(k,j,i) /                           &
     561                            REAL( average_count_3d, KIND=wp )
    435562                      ENDDO
    436563                   ENDDO
     
    442569                   DO  j = nysg, nyng
    443570                      DO  k = nzb, nzt+1
    444                          u_av(k,j,i) = u_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     571                         u_av(k,j,i) = u_av(k,j,i) /                           &
     572                            REAL( average_count_3d, KIND=wp )
    445573                      ENDDO
    446574                   ENDDO
     
    452580                   DO  j = nysg, nyng
    453581                      DO  k = nzb, nzt+1
    454                          v_av(k,j,i) = v_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     582                         v_av(k,j,i) = v_av(k,j,i) /                           &
     583                            REAL( average_count_3d, KIND=wp )
    455584                      ENDDO
    456585                   ENDDO
     
    462591                   DO  j = nysg, nyng
    463592                      DO  k = nzb, nzt+1
    464                          w_av(k,j,i) = w_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     593                         w_av(k,j,i) = w_av(k,j,i) /                           &
     594                            REAL( average_count_3d, KIND=wp )
    465595                      ENDDO
    466596                   ENDDO
    467597                ENDDO
    468598             ENDIF
    469 
    470 !--          Udate thermal indices with derived averages
     599!
     600!--          Udate all thermal index grids with updated averaged input
    471601             CALL bio_calculate_thermal_index_maps ( .TRUE. )
     602
     603!
     604!--     No averaging for UVEM since we are calculating a dose (only sum is
     605!--     calculated and saved to av.nc file)
    472606
    473607        END SELECT
     
    485619!> Check data output for biometeorology model
    486620!------------------------------------------------------------------------------!
    487  SUBROUTINE bio_check_data_output( var, unit )
     621 SUBROUTINE bio_check_data_output( var, unit, i, ilen, k )
    488622
    489623    USE control_parameters,                                                    &
     
    495629    CHARACTER (LEN=*) ::  var      !< The variable in question
    496630
     631    INTEGER(iwp) ::  i      !<
     632    INTEGER(iwp) ::  ilen   !<   
     633    INTEGER(iwp) ::  k      !<
    497634
    498635    SELECT CASE ( TRIM( var ) )
    499      
    500        CASE( 'bio_mrt', 'bio_pet*', 'bio_perct*', 'bio_utci*' )
     636!
     637!-- Allocate a temporary array with the desired output dimensions.
     638       CASE ( 'bio_mrt')
    501639          unit = 'degree_C'
     640          IF ( .NOT. ALLOCATED( tmrt_grid ) )  THEN
     641             ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
     642          ENDIF
     643          tmrt_grid = REAL( bio_fill_value, KIND = wp )
     644
     645       CASE ( 'bio_perct*' )
     646          unit = 'degree_C'
     647          IF ( .NOT. ALLOCATED( perct ) )  THEN
     648             ALLOCATE( perct (nys:nyn,nxl:nxr) )
     649          ENDIF
     650          perct = REAL( bio_fill_value, KIND = wp )
     651
     652       CASE ( 'bio_utci*' )
     653          unit = 'degree_C'
     654          IF ( .NOT. ALLOCATED( utci ) )  THEN
     655             ALLOCATE( utci (nys:nyn,nxl:nxr) )
     656          ENDIF
     657          utci = REAL( bio_fill_value, KIND = wp )
     658
     659       CASE ( 'bio_pet*' )
     660          unit = 'degree_C'
     661          IF ( .NOT. ALLOCATED( pet ) )  THEN
     662             ALLOCATE( pet (nys:nyn,nxl:nxr) )
     663          ENDIF
     664          pet = REAL( bio_fill_value, KIND = wp )
     665
     666       CASE ( 'uvem_vitd3*' )
     667          IF (  .NOT.  uv_exposure )  THEN
     668             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
     669                      'res a namelist &uvexposure_par'
     670             CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
     671          ENDIF
     672          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
     673             message_string = 'illegal value for data_output: "' //            &
     674                              TRIM( var ) // '" & only 2d-horizontal ' //      &
     675                              'cross sections are allowed for this value'
     676             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
     677          ENDIF
     678          unit = 'IU/s'
     679          IF ( .NOT. ALLOCATED( vitd3_exposure ) )  THEN
     680             ALLOCATE( vitd3_exposure(nysg:nyng,nxlg:nxrg) )
     681          ENDIF
     682          vitd3_exposure = 0.0_wp
     683
     684       CASE ( 'uvem_vitd3dose*' )
     685          IF (  .NOT.  uv_exposure )  THEN
     686             message_string = 'output of "' // TRIM( var ) // '" requi' //     &
     687                      'res  a namelist &uvexposure_par'
     688             CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
     689          ENDIF
     690          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
     691             message_string = 'illegal value for data_output: "' //            &
     692                              TRIM( var ) // '" & only 2d-horizontal ' //      &
     693                              'cross sections are allowed for this value'
     694             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
     695          ENDIF
     696          unit = 'IU/av-h'
     697          IF ( .NOT. ALLOCATED( vitd3_exposure_av ) )  THEN
     698             ALLOCATE( vitd3_exposure_av(nysg:nyng,nxlg:nxrg) )
     699          ENDIF
     700          vitd3_exposure_av = 0.0_wp
    502701
    503702       CASE DEFAULT
     
    506705    END SELECT
    507706
    508     IF ( unit /= 'illegal' )  THEN
     707    IF ( thermal_comfort .AND. unit == 'degree_C' )  THEN
    509708!
    510709!--    Futher checks if output belongs to biometeorology
     710!      Break if required modules "radiation" and "humidity" are not running.
    511711       IF ( .NOT.  radiation ) THEN
    512712          message_string = 'output of "' // TRIM( var ) // '" require'         &
     
    541741    IMPLICIT NONE
    542742
    543 
     743!
    544744!--    Disabled as long as radiation model not available
    545745
    546        IF ( .NOT. humidity )  THEN
     746       IF ( thermal_comfort  .AND.  .NOT. humidity )  THEN
    547747          message_string = 'The estimation of thermal comfort requires '    // &
    548748                           'air humidity information, but humidity module ' // &
     
    562762!------------------------------------------------------------------------------!
    563763 SUBROUTINE bio_data_output_2d( av, variable, found, grid, local_pf,           &
    564                                       two_d, nzb_do, nzt_do, fill_value )
    565 
    566     USE indices,                                                               &
    567        ONLY: nxl, nxl, nxr, nxr, nyn, nyn, nys, nys, nzb, nzt
     764                                two_d, nzb_do, nzt_do )
     765
    568766
    569767    USE kinds
     
    571769
    572770    IMPLICIT NONE
    573 
     771!
    574772!-- Input variables
    575773    CHARACTER (LEN=*), INTENT(IN) ::  variable    !< Char identifier to select var for output
     
    577775    INTEGER(iwp), INTENT(IN)      ::  nzb_do      !< Unused. 2D. nz bottom to nz top
    578776    INTEGER(iwp), INTENT(IN)      ::  nzt_do      !< Unused.
    579     REAL(wp), INTENT(in)          ::  fill_value  !< Fill value for unassigned locations
    580 
     777!
    581778!-- Output variables
    582779    CHARACTER (LEN=*), INTENT(OUT) ::  grid   !< Grid type (always "zu1" for biom)
     
    584781    LOGICAL, INTENT(OUT)           ::  two_d  !< Flag parameter that indicates 2D variables, horizontal cross sections, must be .TRUE.
    585782    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< Temp. result grid to return
    586 
     783!
    587784!-- Internal variables
    588     CHARACTER (LEN=:), allocatable ::  variable_short  !< Trimmed version of char variable
    589785    INTEGER(iwp) ::  i        !< Running index, x-dir
    590786    INTEGER(iwp) ::  j        !< Running index, y-dir
     
    593789
    594790
    595     variable_short = TRIM( variable )
    596     IF ( variable_short(1:4) /= 'bio_' ) THEN
    597        found = .FALSE.
    598        grid  = 'none'
    599     ENDIF
    600 
    601791    found = .TRUE.
    602     local_pf = fill_value
    603 
    604     SELECT CASE ( variable_short )
     792    local_pf = bio_fill_value
     793
     794    SELECT CASE ( TRIM( variable ) )
    605795
    606796
    607797        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 )
     798           grid = 'zu1'
     799           two_d = .FALSE.  !< can be calculated for several levels
     800           local_pf = REAL( bio_fill_value, KIND = wp )
     801           DO  l = 1, nmrtbl
     802              i = mrtbl(ix,l)
     803              j = mrtbl(iy,l)
     804              k = mrtbl(iz,l)
     805              IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. j > nyn .OR.   &
     806                 i < nxl .OR. i > nxr ) CYCLE
     807              IF ( av == 0 )  THEN
     808                 IF ( mrt_include_sw )  THEN
     809                    local_pf(i,j,k) = ((human_absorb * mrtinsw(l) +            &
     810                    human_emiss * mrtinlw(l)) /                                &
     811                    (human_emiss * sigma_sb)) ** .25_wp - degc_to_k
     812                 ELSE
     813                    local_pf(i,j,k) = (human_emiss * mrtinlw(l) /              &
     814                                       sigma_sb) ** .25_wp - degc_to_k
     815                 ENDIF
     816              ELSE
     817                 local_pf(i,j,k) = mrt_av_grid(l)
     818              ENDIF
     819           ENDDO
     820
     821
     822        CASE ( 'bio_perct*_xy' )        ! 2d-array
     823           grid = 'zu1'
     824           two_d = .TRUE.
     825           IF ( av == 0 )  THEN
     826              DO  i = nxl, nxr
     827                 DO  j = nys, nyn
     828                    local_pf(i,j,nzb+1) = perct(j,i)
     829                 ENDDO
     830              ENDDO
     831           ELSE
     832              DO  i = nxl, nxr
     833                 DO  j = nys, nyn
     834                    local_pf(i,j,nzb+1) = perct_av(j,i)
     835                 ENDDO
     836              ENDDO
     837           END IF
     838
     839
     840        CASE ( 'bio_utci*_xy' )        ! 2d-array
     841           grid = 'zu1'
     842           two_d = .TRUE.
     843           IF ( av == 0 )  THEN
     844              DO  i = nxl, nxr
     845                 DO  j = nys, nyn
     846                    local_pf(i,j,nzb+1) = utci(j,i)
     847                 ENDDO
     848              ENDDO
     849           ELSE
     850              DO  i = nxl, nxr
     851                 DO  j = nys, nyn
     852                    local_pf(i,j,nzb+1) = utci_av(j,i)
     853                 ENDDO
     854              ENDDO
     855           END IF
     856
     857
     858        CASE ( 'bio_pet*_xy' )        ! 2d-array
     859           grid = 'zu1'
     860           two_d = .TRUE.
     861           IF ( av == 0 )  THEN
     862              DO  i = nxl, nxr
     863                 DO  j = nys, nyn
     864                    local_pf(i,j,nzb+1) = pet(j,i)
     865                 ENDDO
     866              ENDDO
     867           ELSE
     868              DO  i = nxl, nxr
     869                 DO  j = nys, nyn
     870                    local_pf(i,j,nzb+1) = pet_av(j,i)
     871                 ENDDO
     872              ENDDO
     873           END IF
     874
     875           !
     876!--    Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein.
     877!--    However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged.
     878       CASE ( 'uvem_vitd3*_xy' )        ! 2d-array
     879          IF ( av == 0 )  THEN
     880             DO  i = nxl, nxr
     881                DO  j = nys, nyn
     882                   local_pf(i,j,nzb+1) = vitd3_exposure(j,i)
     883                ENDDO
     884             ENDDO
     885          ENDIF
     886
     887          two_d = .TRUE.
     888          grid = 'zu1'
     889
     890       CASE ( 'uvem_vitd3dose*_xy' )        ! 2d-array
     891          IF ( av == 1 )  THEN
     892             DO  i = nxl, nxr
     893                DO  j = nys, nyn
     894                   local_pf(i,j,nzb+1) = vitd3_exposure_av(j,i)
     895                ENDDO
     896             ENDDO
     897          ENDIF
     898
     899          two_d = .TRUE.
     900          grid = 'zu1'
     901
     902
     903       CASE DEFAULT
     904          found = .FALSE.
     905          grid  = 'none'
     906
     907    END SELECT
     908
     909
     910 END SUBROUTINE bio_data_output_2d
     911
     912
     913!------------------------------------------------------------------------------!
     914!
     915! Description:
     916! ------------
     917!> Subroutine defining 3D output variables (dummy, always 2d!)
     918!> data_output_3d 709ff
     919!------------------------------------------------------------------------------!
     920 SUBROUTINE bio_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
     921
     922    USE indices
     923
     924    USE kinds
     925
     926
     927    IMPLICIT NONE
     928!
     929!-- Input variables
     930    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< Char identifier to select var for output
     931    INTEGER(iwp), INTENT(IN) ::  av       !< Use averaged data? 0 = no, 1 = yes?
     932    INTEGER(iwp), INTENT(IN) ::  nzb_do   !< Unused. 2D. nz bottom to nz top
     933    INTEGER(iwp), INTENT(IN) ::  nzt_do   !< Unused.
     934!
     935!-- Output variables
     936    LOGICAL, INTENT(OUT) ::  found   !< Output found?
     937    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< Temp. result grid to return
     938!
     939!-- Internal variables
     940    INTEGER(iwp) ::  l    !< Running index, radiation grid
     941    INTEGER(iwp) ::  i    !< Running index, x-dir
     942    INTEGER(iwp) ::  j    !< Running index, y-dir
     943    INTEGER(iwp) ::  k    !< Running index, z-dir
     944
     945!     REAL(wp) ::  mrt  !< Buffer for mean radiant temperature
     946
     947    found = .TRUE.
     948
     949    SELECT CASE ( TRIM( variable ) )
     950
     951        CASE ( 'bio_mrt' )
     952            local_pf = REAL( bio_fill_value, KIND = sp )
    611953            DO  l = 1, nmrtbl
    612954               i = mrtbl(ix,l)
     
    617959               IF ( av == 0 )  THEN
    618960                  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
     961                     local_pf(i,j,k) = REAL(((human_absorb * mrtinsw(l) +      &
     962                                    human_emiss * mrtinlw(l)) /                &
     963                                    (human_emiss * sigma_sb)) ** .25_wp -      &
     964                                    degc_to_k, kind=sp )
    622965                  ELSE
    623                      local_pf(i,j,k) = (human_emiss * mrtinlw(l) /             &
    624                                         sigma_sb) ** .25_wp - degc_to_k
     966                     local_pf(i,j,k) = REAL((human_emiss * mrtinlw(l) /        &
     967                                    sigma_sb) ** .25_wp - degc_to_k, kind=sp )  !< why not (human_emiss * sigma_sb) as above?
    625968                  ENDIF
    626969               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.
    635             IF ( av == 0 )  THEN
    636               DO  i = nxl, nxr
    637                  DO  j = nys, nyn
    638                     local_pf(i,j,nzb+1) = perct(j,i)
    639                  ENDDO
    640               ENDDO
    641             ELSE
    642               DO  i = nxl, nxr
    643                  DO  j = nys, nyn
    644                     local_pf(i,j,nzb+1) = perct_av(j,i)
    645                  ENDDO
    646               ENDDO
    647             END IF
    648 
    649 
    650         CASE ( 'bio_utci*_xy' )        ! 2d-array
    651             grid = 'zu1'
    652             two_d = .TRUE.
    653             IF ( av == 0 )  THEN
    654               DO  i = nxl, nxr
    655                  DO  j = nys, nyn
    656                     local_pf(i,j,nzb+1) = utci_grid(j,i)
    657                  ENDDO
    658               ENDDO
    659             ELSE
    660               DO  i = nxl, nxr
    661                  DO  j = nys, nyn
    662                     local_pf(i,j,nzb+1) = utci_av_grid(j,i)
    663                  ENDDO
    664               ENDDO
    665             END IF
    666 
    667 
    668         CASE ( 'bio_pet*_xy' )        ! 2d-array
    669             grid = 'zu1'
    670             two_d = .TRUE.
    671             IF ( av == 0 )  THEN
    672               DO  i = nxl, nxr
    673                  DO  j = nys, nyn
    674                     local_pf(i,j,nzb+1) = pet_grid(j,i)
    675                  ENDDO
    676               ENDDO
    677             ELSE
    678               DO  i = nxl, nxr
    679                  DO  j = nys, nyn
    680                     local_pf(i,j,nzb+1) = pet_av_grid(j,i)
    681                  ENDDO
    682               ENDDO
    683             END IF
    684 
    685 
    686        CASE DEFAULT
    687           found = .FALSE.
    688           grid  = 'none'
    689 
    690     END SELECT
    691 
    692 
    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)
     970                  local_pf(i,j,k) = REAL(mrt_av_grid(l), kind=sp)
    762971               ENDIF
    763972            ENDDO
     
    780989
    781990    IMPLICIT NONE
    782 
     991!
    783992!-- Input variables
    784993    CHARACTER (LEN=*), INTENT(IN)  ::  var      !< Name of output variable
    785 
     994!
    786995!-- Output variables
    787996    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
     
    790999
    7911000    LOGICAL, INTENT(OUT)           ::  found    !< Flag if output var is found
    792 
     1001!
    7931002!-- Local variables
    7941003    LOGICAL      :: is2d  !< Var is 2d?
     
    8101019       grid_y = 'y'
    8111020       grid_z = 'zu'
    812        IF ( is2d ) grid_z = 'zu1'
     1021       IF ( is2d  .AND.  var(1:7) /= 'bio_mrt' ) grid_z = 'zu1'
     1022    ENDIF
     1023
     1024    IF ( is2d  .AND.  var(1:4) == 'uvem' ) THEN
     1025       grid_x = 'x'
     1026       grid_y = 'y'
     1027       grid_z = 'zu1'
    8131028    ENDIF
    8141029
     
    8241039
    8251040    IMPLICIT NONE
    826 
     1041!
    8271042!-- Input variables
    8281043    INTEGER(iwp), INTENT(IN) ::  io           !< Unit of the output file
    829 
     1044!
    8301045!-- Internal variables
    8311046    CHARACTER (LEN=86) ::  output_height_chr  !< String for output height
     
    8571072        ONLY: message_string
    8581073
     1074    USE netcdf_data_input_mod,                                                &
     1075        ONLY:  netcdf_data_input_uvem, uvem_projarea_f, uvem_radiance_f,      &
     1076               uvem_irradiance_f, uvem_integration_f, building_obstruction_f
     1077
    8591078    IMPLICIT NONE
    860 
     1079!
    8611080!-- Internal vriables
    8621081    REAL ( wp )  :: height  !< current height in meters
    863 
    864     INTEGER ( iwp )  :: i  !< iteration index
    865 
     1082!
    8661083!-- Determine cell level corresponding to 1.1 m above ground level
    8671084!   (gravimetric center of sample human)
     
    8761093
    8771094    IF ( .NOT. radiation_interactions ) THEN
    878        message_string = 'The mrt calculation requires ' // &
    879                         'enabled radiation_interactions but it ' // &
     1095       message_string = 'The mrt calculation requires ' //                     &
     1096                        'enabled radiation_interactions but it ' //            &
    8801097                        'is disabled!'
    8811098       CALL message( 'check_parameters', 'PAHU03', 0, 0, -1, 6, 0 )
    8821099    ENDIF
    8831100
     1101!
     1102!-- Init UVEM and load lookup tables
     1103    CALL netcdf_data_input_uvem
     1104
    8841105 END SUBROUTINE bio_init
    885 
    886 
    887 !------------------------------------------------------------------------------!
    888 ! Description:
    889 ! ------------
    890 !> Allocate biom arrays and define pointers if required
    891 !> init_3d_model 1047ff
    892 !------------------------------------------------------------------------------!
    893  SUBROUTINE bio_init_arrays
    894 
    895     IMPLICIT NONE
    896 
    897 !-- Allocate a temporary array with the desired output dimensions.
    898 !   Initialization omitted for performance, will be overwritten anyway
    899     IF ( .NOT. ALLOCATED( tmrt_grid ) ) THEN
    900       ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
    901     ENDIF
    902 
    903     IF ( bio_perct ) THEN
    904       IF ( .NOT. ALLOCATED( perct ) ) THEN
    905         ALLOCATE( perct (nys:nyn,nxl:nxr) )
    906       ENDIF
    907     ENDIF
    908 
    909     IF ( bio_utci ) THEN
    910       IF ( .NOT. ALLOCATED( utci_grid ) ) THEN
    911         ALLOCATE( utci_grid (nys:nyn,nxl:nxr) )
    912       ENDIF
    913     ENDIF
    914 
    915     IF ( bio_pet ) THEN
    916       IF ( .NOT. ALLOCATED( pet_grid ) ) THEN
    917         ALLOCATE( pet_grid (nys:nyn,nxl:nxr) )
    918       ENDIF
    919     END IF
    920 
    921     IF ( bio_perct_av ) THEN
    922       IF ( .NOT. ALLOCATED( perct_av ) ) THEN
    923         ALLOCATE( perct_av (nys:nyn,nxl:nxr) )
    924       ENDIF
    925     ENDIF
    926 
    927     IF ( bio_utci_av ) THEN
    928       IF ( .NOT. ALLOCATED( utci_av_grid ) ) THEN
    929         ALLOCATE( utci_av_grid (nys:nyn,nxl:nxr) )
    930       ENDIF
    931     ENDIF
    932 
    933     IF ( bio_pet_av ) THEN
    934       IF ( .NOT. ALLOCATED( pet_av_grid ) ) THEN
    935         ALLOCATE( pet_av_grid (nys:nyn,nxl:nxr) )
    936       ENDIF
    937 
    938     END IF
    939 
    940  END SUBROUTINE bio_init_arrays
    9411106
    9421107
     
    9591124                                          bio_perct_av,                        &
    9601125                                          bio_utci,                            &
    961                                           bio_utci_av
     1126                                          bio_utci_av,                         &
     1127                                          thermal_comfort,                     &
     1128!
     1129!-- UVEM namelist parameters
     1130                                          clothing,                            &
     1131                                          consider_obstructions,               &
     1132                                          orientation_angle,                   &
     1133                                          sun_in_south,                        &
     1134                                          turn_to_sun,                         &
     1135                                          uv_exposure
    9621136
    9631137
     
    10031177
    10041178    LOGICAL, INTENT(IN)         ::  av    !< use averaged input?
     1179!
    10051180!-- Internal variables
    10061181    INTEGER(iwp)                ::  i     !< Running index, x-dir, radiation coordinates
     
    10091184    INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
    10101185
    1011 
     1186!
    10121187!-- Calculate biometeorology MRT from local radiation
    1013 !-- fluxes calculated by RTM and assign into 2D grid
    1014     tmrt_grid = -999.
     1188!  fluxes calculated by RTM and assign into 2D grid
     1189    tmrt_grid = bio_fill_value
    10151190    DO  l = 1, nmrtbl
    10161191       i = mrtbl(ix,l)
     
    10421217
    10431218    IMPLICIT NONE
    1044 
     1219!
    10451220!-- Input variables
    10461221    LOGICAL,      INTENT ( IN ) ::  average_input  !< Determine averaged input conditions?
    10471222    INTEGER(iwp), INTENT ( IN ) ::  i     !< Running index, x-dir
    10481223    INTEGER(iwp), INTENT ( IN ) ::  j     !< Running index, y-dir
    1049 
     1224!
    10501225!-- Output parameters
    10511226    REAL(wp), INTENT ( OUT )    ::  tmrt  !< Mean radiant temperature        (°C)
     
    10541229    REAL(wp), INTENT ( OUT )    ::  ws    !< Wind speed    (local level)     (m/s)
    10551230    REAL(wp), INTENT ( OUT )    ::  pair  !< Air pressure                    (hPa)
    1056 
     1231!
    10571232!-- Internal variables
    10581233    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
     1234!     INTEGER(iwp)                ::  ir    !< Running index, x-dir, radiation coordinates
     1235!     INTEGER(iwp)                ::  jr    !< Running index, y-dir, radiation coordinates
     1236!     INTEGER(iwp)                ::  kr    !< Running index, y-dir, radiation coordinates
    10621237    INTEGER(iwp)                ::  k_wind  !< Running index, z-dir, wind speed only
    1063     INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
     1238!     INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
    10641239
    10651240    REAL(wp)                    ::  vp_sat  !< Saturation vapor pressure     (hPa)
    10661241
    1067 
     1242!
    10681243!-- Determine cell level closest to 1.1m above ground
    10691244!   by making use of truncation due to int cast
     
    10761251       k_wind = k + 1_iwp
    10771252    ENDIF
    1078 
     1253!
    10791254!-- Determine local values:
    10801255    IF ( average_input ) THEN
     1256!
    10811257!--    Calculate ta from Tp assuming dry adiabatic laps rate
    10821258       ta = pt_av(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
    10831259
    1084        vp = -999._wp
     1260       vp = bio_fill_value
    10851261       IF ( humidity .AND. ALLOCATED( q_av ) ) THEN
    10861262          vp = q_av(k, j, i)
     
    10911267          0.5_wp * ABS( w_av(k_wind, j, i) + w_av(k_wind+1, j, i) ) )
    10921268    ELSE
     1269!
    10931270!-- Calculate ta from Tp assuming dry adiabatic laps rate
    10941271       ta = pt(k, j, i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
    10951272
    1096        vp = -999._wp
     1273       vp = bio_fill_value
    10971274       IF ( humidity ) THEN
    10981275          vp = q(k, j, i)
     
    11041281
    11051282    ENDIF
    1106 
     1283!
    11071284!-- Local air pressure
    11081285    pair = surface_pressure
     
    11161293       IF ( vp > vp_sat ) vp = vp_sat
    11171294    ENDIF
    1118 
     1295!
    11191296!-- local mtr value at [i,j]
    1120     tmrt = -999.  !< this can be a valid result (e.g. for inside some ostacle)
     1297    tmrt = bio_fill_value  !< this can be a valid result (e.g. for inside some ostacle)
    11211298    IF ( radiation ) THEN
     1299!
    11221300!--    Use MRT from RTM precalculated in tmrt_grid
    11231301       tmrt = tmrt_grid(j,i)
     
    11361314
    11371315    IMPLICIT NONE
    1138 
     1316!
    11391317!-- Input attributes
    11401318    LOGICAL, INTENT ( IN ) ::  av  !< Calculate based on averaged input conditions?
    1141 
     1319!
    11421320!-- Internal variables
    1143     INTEGER(iwp) ::  i, j      !< Running index
    1144 
    1145     REAL(wp) ::  clo           !< Clothing index                (no dimension)
    1146     REAL(wp) ::  ta            !< Air temperature                  (°C)
    1147     REAL(wp) ::  vp            !< Vapour pressure                  (hPa)
    1148     REAL(wp) ::  ws            !< Wind speed    (local level)      (m/s)
    1149     REAL(wp) ::  pair          !< Air pressure                     (hPa)
    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 
     1321    INTEGER(iwp) ::  i, j     !< Running index
     1322
     1323    REAL(wp) ::  clo          !< Clothing index                (no dimension)
     1324    REAL(wp) ::  ta           !< Air temperature                  (°C)
     1325    REAL(wp) ::  vp           !< Vapour pressure                  (hPa)
     1326    REAL(wp) ::  ws           !< Wind speed    (local level)      (m/s)
     1327    REAL(wp) ::  pair         !< Air pressure                     (hPa)
     1328    REAL(wp) ::  perct_ij     !< Perceived temperature            (°C)
     1329    REAL(wp) ::  utci_ij      !< Universal thermal climate index  (°C)
     1330    REAL(wp) ::  pet_ij       !< Physiologically equivalent temperature  (°C)
     1331    REAL(wp) ::  tmrt_ij      !< Mean radiant temperature         (°C)
     1332
     1333!
    11571334!-- fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...)
    1158     CALL bio_calculate_mrt_grid ( av )
    1159 
     1335    IF ( simulated_time > 0.0_wp ) THEN
     1336       CALL bio_calculate_mrt_grid ( av )
     1337    ENDIF
    11601338
    11611339    DO i = nxl, nxr
    11621340       DO j = nys, nyn
     1341!
    11631342!--       Determine local input conditions
    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 
     1343          tmrt_ij = bio_fill_value
     1344          vp      = bio_fill_value
     1345!
     1346!--       Determine input only if
     1347          IF ( simulated_time > 0.0_wp ) THEN
     1348             CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp,           &
     1349                ws, pair, tmrt_ij )
     1350          END IF
     1351!
    11681352!--       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
    1172              utci_tmp = -999._wp
    1173           ELSE
     1353          pet_ij   = bio_fill_value   !< set fail value, e.g. valid for within
     1354          perct_ij = bio_fill_value   !< some obstacle
     1355          utci_ij  = bio_fill_value
     1356          IF ( .NOT. ( tmrt_ij <= -998._wp .OR. vp <= -998._wp ) ) THEN
     1357!
    11741358!--          Calculate static thermal indices based on local tmrt
    1175              clo = -999._wp
    1176              
     1359             clo = bio_fill_value
     1360
    11771361             IF ( bio_perct ) THEN
     1362!
    11781363!--          Estimate local perceived temperature
    1179                 CALL calculate_perct_static( ta, vp, ws, tmrt_tmp, pair, clo,  &
    1180                    perct_tmp )
     1364                CALL calculate_perct_static( ta, vp, ws, tmrt_ij, pair, clo,   &
     1365                   perct_ij )
    11811366             ENDIF
    1182              
     1367
    11831368             IF ( bio_utci ) THEN
     1369!
    11841370!--          Estimate local universal thermal climate index
    1185                 CALL calculate_utci_static( ta, vp, ws, tmrt_tmp,              &
    1186                    bio_output_height, utci_tmp )
     1371                CALL calculate_utci_static( ta, vp, ws, tmrt_ij,               &
     1372                   bio_output_height, utci_ij )
    11871373             ENDIF
    1188              
     1374
    11891375             IF ( bio_pet ) THEN
     1376!
    11901377!--          Estimate local physiologically equivalent temperature
    1191                 CALL calculate_pet_static( ta, vp, ws, tmrt_tmp, pair, pet_tmp )
     1378                CALL calculate_pet_static( ta, vp, ws, tmrt_ij, pair, pet_ij )
    11921379             ENDIF
    11931380          END IF
     
    11951382
    11961383          IF ( av ) THEN
     1384!
    11971385!--          Write results for selected averaged indices
    11981386             IF ( bio_perct_av )  THEN
    1199                 perct_av(j, i) = perct_tmp
     1387                perct_av(j, i) = perct_ij
    12001388             END IF
    12011389             IF ( bio_utci_av ) THEN
    1202                 utci_av_grid(j, i) = utci_tmp
     1390                utci_av(j, i) = utci_ij
    12031391             END IF
    12041392             IF ( bio_pet_av ) THEN
    1205                 pet_av_grid(j, i)  = pet_tmp
     1393                pet_av(j, i)  = pet_ij
    12061394             END IF
    12071395          ELSE
     1396!
    12081397!--          Write result for selected indices
    12091398             IF ( bio_perct )  THEN
    1210                 perct(j, i) = perct_tmp
     1399                perct(j, i) = perct_ij
    12111400             END IF
    12121401             IF ( bio_utci ) THEN
    1213                 utci_grid(j, i) = utci_tmp
     1402                utci(j, i) = utci_ij
    12141403             END IF
    12151404             IF ( bio_pet ) THEN
    1216                 pet_grid(j, i)  = pet_tmp
     1405                pet(j, i)  = pet_ij
    12171406             END IF
    12181407          END IF
     
    12321421
    12331422    IMPLICIT NONE
    1234 
     1423!
    12351424!-- Input parameters
    12361425    REAL(wp), INTENT ( IN )  ::  ta   !< Air temperature                  (°C)
     
    12461435                                         !  (without metabolism!)         (W)
    12471436    INTEGER(iwp), INTENT ( IN ) ::  sex  !< Sex of agent (1 = male, 2 = female)
    1248 
     1437!
    12491438!-- Both, input and output parameters
    12501439    Real(wp), INTENT ( INOUT )  ::  energy_storage    !< Energy storage   (W/m²)
     
    12531442    Real(wp), INTENT ( INOUT )  ::  actlev  !< Individuals activity level
    12541443                                            !  per unit surface area      (W/m²)
     1444!
    12551445!-- Output parameters
    12561446    REAL(wp), INTENT ( OUT ) ::  ipt    !< Instationary perceived temp.   (°C)
    1257 
     1447!
    12581448!-- If clo equals the initial value, this is the initial call
    12591449    IF ( clo <= -998._wp ) THEN
     1450!
    12601451!--    Initialize instationary perceived temperature with personalized
    12611452!      PT as an initial guess, set actlev and clo
     
    12641455          ipt )
    12651456    ELSE
     1457!
    12661458!--    Estimate local instatinoary perceived temperature
    12671459       CALL ipt_cycle ( ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo,     &
     
    12871479!> www.utci.org
    12881480!------------------------------------------------------------------------------!
    1289  SUBROUTINE calculate_utci_static( ta_in, vp, ws_hag, tmrt, hag, utci )
     1481 SUBROUTINE calculate_utci_static( ta_in, vp, ws_hag, tmrt, hag, utci_ij )
    12901482
    12911483    IMPLICIT NONE
    1292 
     1484!
    12931485!-- 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)
     1486    REAL(WP), INTENT ( IN )  ::  ta_in    !< Local air temperature (°C)
     1487    REAL(WP), INTENT ( IN )  ::  vp       !< Loacl vapour pressure (hPa)
     1488    REAL(WP), INTENT ( IN )  ::  ws_hag   !< Incident wind speed (m/s)
     1489    REAL(WP), INTENT ( IN )  ::  tmrt     !< Local mean radiant temperature (°C)
     1490    REAL(WP), INTENT ( IN )  ::  hag      !< Height of wind speed input (m)
     1491!
    12991492!-- 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 
     1493    REAL(wp), INTENT ( OUT ) ::  utci_ij  !< Universal Thermal Climate Index (°C)
     1494
     1495    REAL(WP) ::  ta           !< air temperature modified by offset (°C)
     1496    REAL(WP) ::  pa           !< air pressure in kPa      (kPa)
     1497    REAL(WP) ::  d_tmrt       !< delta-tmrt               (°C)
     1498    REAL(WP) ::  va           !< wind speed at 10 m above ground level    (m/s)
     1499    REAL(WP) ::  offset       !< utci deviation by ta cond. exceeded      (°C)
     1500    REAL(WP) ::  part_ta      !< Air temperature related part of the regression
     1501    REAL(WP) ::  ta2          !< 2 times ta
     1502    REAL(WP) ::  ta3          !< 3 times ta
     1503    REAL(WP) ::  ta4          !< 4 times ta
     1504    REAL(WP) ::  ta5          !< 5 times ta
     1505    REAL(WP) ::  ta6          !< 6 times ta
     1506    REAL(WP) ::  part_va      !< Vapour pressure related part of the regression
     1507    REAL(WP) ::  va2          !< 2 times va
     1508    REAL(WP) ::  va3          !< 3 times va
     1509    REAL(WP) ::  va4          !< 4 times va
     1510    REAL(WP) ::  va5          !< 5 times va
     1511    REAL(WP) ::  va6          !< 6 times va
     1512    REAL(WP) ::  part_d_tmrt  !< Mean radiant temp. related part of the reg.
     1513    REAL(WP) ::  d_tmrt2      !< 2 times d_tmrt
     1514    REAL(WP) ::  d_tmrt3      !< 3 times d_tmrt
     1515    REAL(WP) ::  d_tmrt4      !< 4 times d_tmrt
     1516    REAL(WP) ::  d_tmrt5      !< 5 times d_tmrt
     1517    REAL(WP) ::  d_tmrt6      !< 6 times d_tmrt
     1518    REAL(WP) ::  part_pa      !< Air pressure related part of the regression
     1519    REAL(WP) ::  pa2          !< 2 times pa
     1520    REAL(WP) ::  pa3          !< 3 times pa
     1521    REAL(WP) ::  pa4          !< 4 times pa
     1522    REAL(WP) ::  pa5          !< 5 times pa
     1523    REAL(WP) ::  pa6          !< 6 times pa
     1524    REAL(WP) ::  part_pa2     !< Air pressure^2 related part of the regression
     1525    REAL(WP) ::  part_pa3     !< Air pressure^3 related part of the regression
     1526    REAL(WP) ::  part_pa46    !< Air pressure^4-6 related part of the regression
     1527
     1528!
    13091529!-- Initialize
    13101530    offset = 0._wp
    13111531    ta = ta_in
    13121532    d_tmrt = tmrt - ta_in
    1313 
     1533!
    13141534!-- Use vapour pressure in kpa
    13151535    pa = vp / 10.0_wp
    1316 
     1536!
    13171537!-- Wind altitude correction from hag to 10m after Broede et al. (2012), eq.3
    13181538!   z(0) is set to 0.01 according to UTCI profile definition
    13191539    va = ws_hag *  log ( 10.0_wp / 0.01_wp ) / log ( hag / 0.01_wp )
    1320 
     1540!
    13211541!-- Check if input values in range after Broede et al. (2012)
    13221542    IF ( ( d_tmrt > 70._wp ) .OR. ( d_tmrt < -30._wp ) .OR.                    &
    13231543       ( vp >= 50._wp ) ) THEN
    1324        utci = -999._wp
     1544       utci_ij = bio_fill_value
    13251545       RETURN
    13261546    ENDIF
    1327 
     1547!
    13281548!-- Apply eq. 2 in Broede et al. (2012) for ta out of bounds
    13291549    IF ( ta > 50._wp ) THEN
     
    13351555       ta = -50._wp
    13361556    ENDIF
    1337 
     1557!
    13381558!-- For routine application. For wind speeds and relative
    13391559!   humidity values below 0.5 m/s or 5%, respectively, the
     
    13421562    IF ( va > 17._wp ) va = 17._wp
    13431563
     1564!
     1565!-- Pre-calculate multiples of input parameters to save time later
     1566
     1567    ta2 = ta  * ta
     1568    ta3 = ta2 * ta
     1569    ta4 = ta3 * ta
     1570    ta5 = ta4 * ta
     1571    ta6 = ta5 * ta
     1572
     1573    va2 = va  * va
     1574    va3 = va2 * va
     1575    va4 = va3 * va
     1576    va5 = va4 * va
     1577    va6 = va5 * va
     1578
     1579    d_tmrt2 = d_tmrt  * d_tmrt
     1580    d_tmrt3 = d_tmrt2 * d_tmrt
     1581    d_tmrt4 = d_tmrt3 * d_tmrt
     1582    d_tmrt5 = d_tmrt4 * d_tmrt
     1583    d_tmrt6 = d_tmrt5 * d_tmrt
     1584
     1585    pa2 = pa  * pa
     1586    pa3 = pa2 * pa
     1587    pa4 = pa3 * pa
     1588    pa5 = pa4 * pa
     1589    pa6 = pa5 * pa
     1590
     1591!
     1592!-- Pre-calculate parts of the regression equation
     1593    part_ta = (  6.07562052e-01_wp )       +                                   &
     1594              ( -2.27712343e-02_wp ) * ta  +                                   &
     1595              (  8.06470249e-04_wp ) * ta2 +                                   &
     1596              ( -1.54271372e-04_wp ) * ta3 +                                   &
     1597              ( -3.24651735e-06_wp ) * ta4 +                                   &
     1598              (  7.32602852e-08_wp ) * ta5 +                                   &
     1599              (  1.35959073e-09_wp ) * ta6
     1600
     1601    part_va = ( -2.25836520e+00_wp ) * va +                                    &
     1602        (  8.80326035e-02_wp ) * ta  * va +                                    &
     1603        (  2.16844454e-03_wp ) * ta2 * va +                                    &
     1604        ( -1.53347087e-05_wp ) * ta3 * va +                                    &
     1605        ( -5.72983704e-07_wp ) * ta4 * va +                                    &
     1606        ( -2.55090145e-09_wp ) * ta5 * va +                                    &
     1607        ( -7.51269505e-01_wp ) *       va2 +                                   &
     1608        ( -4.08350271e-03_wp ) * ta  * va2 +                                   &
     1609        ( -5.21670675e-05_wp ) * ta2 * va2 +                                   &
     1610        (  1.94544667e-06_wp ) * ta3 * va2 +                                   &
     1611        (  1.14099531e-08_wp ) * ta4 * va2 +                                   &
     1612        (  1.58137256e-01_wp ) *       va3 +                                   &
     1613        ( -6.57263143e-05_wp ) * ta  * va3 +                                   &
     1614        (  2.22697524e-07_wp ) * ta2 * va3 +                                   &
     1615        ( -4.16117031e-08_wp ) * ta3 * va3 +                                   &
     1616        ( -1.27762753e-02_wp ) *       va4 +                                   &
     1617        (  9.66891875e-06_wp ) * ta  * va4 +                                   &
     1618        (  2.52785852e-09_wp ) * ta2 * va4 +                                   &
     1619        (  4.56306672e-04_wp ) *       va5 +                                   &
     1620        ( -1.74202546e-07_wp ) * ta  * va5 +                                   &
     1621        ( -5.91491269e-06_wp ) * va6
     1622
     1623    part_d_tmrt = (  3.98374029e-01_wp ) *       d_tmrt +                      &
     1624            (  1.83945314e-04_wp ) * ta  *       d_tmrt +                      &
     1625            ( -1.73754510e-04_wp ) * ta2 *       d_tmrt +                      &
     1626            ( -7.60781159e-07_wp ) * ta3 *       d_tmrt +                      &
     1627            (  3.77830287e-08_wp ) * ta4 *       d_tmrt +                      &
     1628            (  5.43079673e-10_wp ) * ta5 *       d_tmrt +                      &
     1629            ( -2.00518269e-02_wp ) *       va  * d_tmrt +                      &
     1630            (  8.92859837e-04_wp ) * ta  * va  * d_tmrt +                      &
     1631            (  3.45433048e-06_wp ) * ta2 * va  * d_tmrt +                      &
     1632            ( -3.77925774e-07_wp ) * ta3 * va  * d_tmrt +                      &
     1633            ( -1.69699377e-09_wp ) * ta4 * va  * d_tmrt +                      &
     1634            (  1.69992415e-04_wp ) *       va2 * d_tmrt +                      &
     1635            ( -4.99204314e-05_wp ) * ta  * va2 * d_tmrt +                      &
     1636            (  2.47417178e-07_wp ) * ta2 * va2 * d_tmrt +                      &
     1637            (  1.07596466e-08_wp ) * ta3 * va2 * d_tmrt +                      &
     1638            (  8.49242932e-05_wp ) *       va3 * d_tmrt +                      &
     1639            (  1.35191328e-06_wp ) * ta  * va3 * d_tmrt +                      &
     1640            ( -6.21531254e-09_wp ) * ta2 * va3 * d_tmrt +                      &
     1641            ( -4.99410301e-06_wp ) * va4 *       d_tmrt +                      &
     1642            ( -1.89489258e-08_wp ) * ta  * va4 * d_tmrt +                      &
     1643            (  8.15300114e-08_wp ) *       va5 * d_tmrt +                      &
     1644            (  7.55043090e-04_wp ) *             d_tmrt2 +                     &
     1645            ( -5.65095215e-05_wp ) * ta  *       d_tmrt2 +                     &
     1646            ( -4.52166564e-07_wp ) * ta2 *       d_tmrt2 +                     &
     1647            (  2.46688878e-08_wp ) * ta3 *       d_tmrt2 +                     &
     1648            (  2.42674348e-10_wp ) * ta4 *       d_tmrt2 +                     &
     1649            (  1.54547250e-04_wp ) *       va  * d_tmrt2 +                     &
     1650            (  5.24110970e-06_wp ) * ta  * va  * d_tmrt2 +                     &
     1651            ( -8.75874982e-08_wp ) * ta2 * va  * d_tmrt2 +                     &
     1652            ( -1.50743064e-09_wp ) * ta3 * va  * d_tmrt2 +                     &
     1653            ( -1.56236307e-05_wp ) *       va2 * d_tmrt2 +                     &
     1654            ( -1.33895614e-07_wp ) * ta  * va2 * d_tmrt2 +                     &
     1655            (  2.49709824e-09_wp ) * ta2 * va2 * d_tmrt2 +                     &
     1656            (  6.51711721e-07_wp ) *       va3 * d_tmrt2 +                     &
     1657            (  1.94960053e-09_wp ) * ta  * va3 * d_tmrt2 +                     &
     1658            ( -1.00361113e-08_wp ) *       va4 * d_tmrt2 +                     &
     1659            ( -1.21206673e-05_wp ) *             d_tmrt3 +                     &
     1660            ( -2.18203660e-07_wp ) * ta  *       d_tmrt3 +                     &
     1661            (  7.51269482e-09_wp ) * ta2 *       d_tmrt3 +                     &
     1662            (  9.79063848e-11_wp ) * ta3 *       d_tmrt3 +                     &
     1663            (  1.25006734e-06_wp ) *       va  * d_tmrt3 +                     &
     1664            ( -1.81584736e-09_wp ) * ta  * va  * d_tmrt3 +                     &
     1665            ( -3.52197671e-10_wp ) * ta2 * va  * d_tmrt3 +                     &
     1666            ( -3.36514630e-08_wp ) *       va2 * d_tmrt3 +                     &
     1667            (  1.35908359e-10_wp ) * ta  * va2 * d_tmrt3 +                     &
     1668            (  4.17032620e-10_wp ) *       va3 * d_tmrt3 +                     &
     1669            ( -1.30369025e-09_wp ) *             d_tmrt4 +                     &
     1670            (  4.13908461e-10_wp ) * ta  *       d_tmrt4 +                     &
     1671            (  9.22652254e-12_wp ) * ta2 *       d_tmrt4 +                     &
     1672            ( -5.08220384e-09_wp ) *       va  * d_tmrt4 +                     &
     1673            ( -2.24730961e-11_wp ) * ta  * va  * d_tmrt4 +                     &
     1674            (  1.17139133e-10_wp ) *       va2 * d_tmrt4 +                     &
     1675            (  6.62154879e-10_wp ) *             d_tmrt5 +                     &
     1676            (  4.03863260e-13_wp ) * ta  *       d_tmrt5 +                     &
     1677            (  1.95087203e-12_wp ) *       va  * d_tmrt5 +                     &
     1678            ( -4.73602469e-12_wp ) *             d_tmrt6
     1679
     1680    part_pa = (  5.12733497e+00_wp ) *                pa +                     &
     1681       ( -3.12788561e-01_wp ) * ta  *                 pa +                     &
     1682       ( -1.96701861e-02_wp ) * ta2 *                 pa +                     &
     1683       (  9.99690870e-04_wp ) * ta3 *                 pa +                     &
     1684       (  9.51738512e-06_wp ) * ta4 *                 pa +                     &
     1685       ( -4.66426341e-07_wp ) * ta5 *                 pa +                     &
     1686       (  5.48050612e-01_wp ) *       va  *           pa +                     &
     1687       ( -3.30552823e-03_wp ) * ta  * va  *           pa +                     &
     1688       ( -1.64119440e-03_wp ) * ta2 * va  *           pa +                     &
     1689       ( -5.16670694e-06_wp ) * ta3 * va  *           pa +                     &
     1690       (  9.52692432e-07_wp ) * ta4 * va  *           pa +                     &
     1691       ( -4.29223622e-02_wp ) *       va2 *           pa +                     &
     1692       (  5.00845667e-03_wp ) * ta  * va2 *           pa +                     &
     1693       (  1.00601257e-06_wp ) * ta2 * va2 *           pa +                     &
     1694       ( -1.81748644e-06_wp ) * ta3 * va2 *           pa +                     &
     1695       ( -1.25813502e-03_wp ) *       va3 *           pa +                     &
     1696       ( -1.79330391e-04_wp ) * ta  * va3 *           pa +                     &
     1697       (  2.34994441e-06_wp ) * ta2 * va3 *           pa +                     &
     1698       (  1.29735808e-04_wp ) *       va4 *           pa +                     &
     1699       (  1.29064870e-06_wp ) * ta  * va4 *           pa +                     &
     1700       ( -2.28558686e-06_wp ) *       va5 *           pa +                     &
     1701       ( -3.69476348e-02_wp ) *             d_tmrt  * pa +                     &
     1702       (  1.62325322e-03_wp ) * ta  *       d_tmrt  * pa +                     &
     1703       ( -3.14279680e-05_wp ) * ta2 *       d_tmrt  * pa +                     &
     1704       (  2.59835559e-06_wp ) * ta3 *       d_tmrt  * pa +                     &
     1705       ( -4.77136523e-08_wp ) * ta4 *       d_tmrt  * pa +                     &
     1706       (  8.64203390e-03_wp ) *       va  * d_tmrt  * pa +                     &
     1707       ( -6.87405181e-04_wp ) * ta  * va  * d_tmrt  * pa +                     &
     1708       ( -9.13863872e-06_wp ) * ta2 * va  * d_tmrt  * pa +                     &
     1709       (  5.15916806e-07_wp ) * ta3 * va  * d_tmrt  * pa +                     &
     1710       ( -3.59217476e-05_wp ) *       va2 * d_tmrt  * pa +                     &
     1711       (  3.28696511e-05_wp ) * ta  * va2 * d_tmrt  * pa +                     &
     1712       ( -7.10542454e-07_wp ) * ta2 * va2 * d_tmrt  * pa +                     &
     1713       ( -1.24382300e-05_wp ) *       va3 * d_tmrt  * pa +                     &
     1714       ( -7.38584400e-09_wp ) * ta  * va3 * d_tmrt  * pa +                     &
     1715       (  2.20609296e-07_wp ) *       va4 * d_tmrt  * pa +                     &
     1716       ( -7.32469180e-04_wp ) *             d_tmrt2 * pa +                     &
     1717       ( -1.87381964e-05_wp ) * ta  *       d_tmrt2 * pa +                     &
     1718       (  4.80925239e-06_wp ) * ta2 *       d_tmrt2 * pa +                     &
     1719       ( -8.75492040e-08_wp ) * ta3 *       d_tmrt2 * pa +                     &
     1720       (  2.77862930e-05_wp ) *       va  * d_tmrt2 * pa +                     &
     1721       ( -5.06004592e-06_wp ) * ta  * va  * d_tmrt2 * pa +                     &
     1722       (  1.14325367e-07_wp ) * ta2 * va  * d_tmrt2 * pa +                     &
     1723       (  2.53016723e-06_wp ) *       va2 * d_tmrt2 * pa +                     &
     1724       ( -1.72857035e-08_wp ) * ta  * va2 * d_tmrt2 * pa +                     &
     1725       ( -3.95079398e-08_wp ) *       va3 * d_tmrt2 * pa +                     &
     1726       ( -3.59413173e-07_wp ) *             d_tmrt3 * pa +                     &
     1727       (  7.04388046e-07_wp ) * ta  *       d_tmrt3 * pa +                     &
     1728       ( -1.89309167e-08_wp ) * ta2 *       d_tmrt3 * pa +                     &
     1729       ( -4.79768731e-07_wp ) *       va  * d_tmrt3 * pa +                     &
     1730       (  7.96079978e-09_wp ) * ta  * va  * d_tmrt3 * pa +                     &
     1731       (  1.62897058e-09_wp ) *       va2 * d_tmrt3 * pa +                     &
     1732       (  3.94367674e-08_wp ) *             d_tmrt4 * pa +                     &
     1733       ( -1.18566247e-09_wp ) * ta *        d_tmrt4 * pa +                     &
     1734       (  3.34678041e-10_wp ) *       va  * d_tmrt4 * pa +                     &
     1735       ( -1.15606447e-10_wp ) *             d_tmrt5 * pa
     1736
     1737    part_pa2 = ( -2.80626406e+00_wp ) *               pa2 +                    &
     1738       (  5.48712484e-01_wp ) * ta  *                 pa2 +                    &
     1739       ( -3.99428410e-03_wp ) * ta2 *                 pa2 +                    &
     1740       ( -9.54009191e-04_wp ) * ta3 *                 pa2 +                    &
     1741       (  1.93090978e-05_wp ) * ta4 *                 pa2 +                    &
     1742       ( -3.08806365e-01_wp ) *       va *            pa2 +                    &
     1743       (  1.16952364e-02_wp ) * ta  * va *            pa2 +                    &
     1744       (  4.95271903e-04_wp ) * ta2 * va *            pa2 +                    &
     1745       ( -1.90710882e-05_wp ) * ta3 * va *            pa2 +                    &
     1746       (  2.10787756e-03_wp ) *       va2 *           pa2 +                    &
     1747       ( -6.98445738e-04_wp ) * ta  * va2 *           pa2 +                    &
     1748       (  2.30109073e-05_wp ) * ta2 * va2 *           pa2 +                    &
     1749       (  4.17856590e-04_wp ) *       va3 *           pa2 +                    &
     1750       ( -1.27043871e-05_wp ) * ta  * va3 *           pa2 +                    &
     1751       ( -3.04620472e-06_wp ) *       va4 *           pa2 +                    &
     1752       (  5.14507424e-02_wp ) *             d_tmrt  * pa2 +                    &
     1753       ( -4.32510997e-03_wp ) * ta  *       d_tmrt  * pa2 +                    &
     1754       (  8.99281156e-05_wp ) * ta2 *       d_tmrt  * pa2 +                    &
     1755       ( -7.14663943e-07_wp ) * ta3 *       d_tmrt  * pa2 +                    &
     1756       ( -2.66016305e-04_wp ) *       va  * d_tmrt  * pa2 +                    &
     1757       (  2.63789586e-04_wp ) * ta  * va  * d_tmrt  * pa2 +                    &
     1758       ( -7.01199003e-06_wp ) * ta2 * va  * d_tmrt  * pa2 +                    &
     1759       ( -1.06823306e-04_wp ) *       va2 * d_tmrt  * pa2 +                    &
     1760       (  3.61341136e-06_wp ) * ta  * va2 * d_tmrt  * pa2 +                    &
     1761       (  2.29748967e-07_wp ) *       va3 * d_tmrt  * pa2 +                    &
     1762       (  3.04788893e-04_wp ) *             d_tmrt2 * pa2 +                    &
     1763       ( -6.42070836e-05_wp ) * ta  *       d_tmrt2 * pa2 +                    &
     1764       (  1.16257971e-06_wp ) * ta2 *       d_tmrt2 * pa2 +                    &
     1765       (  7.68023384e-06_wp ) *       va  * d_tmrt2 * pa2 +                    &
     1766       ( -5.47446896e-07_wp ) * ta  * va  * d_tmrt2 * pa2 +                    &
     1767       ( -3.59937910e-08_wp ) *       va2 * d_tmrt2 * pa2 +                    &
     1768       ( -4.36497725e-06_wp ) *             d_tmrt3 * pa2 +                    &
     1769       (  1.68737969e-07_wp ) * ta  *       d_tmrt3 * pa2 +                    &
     1770       (  2.67489271e-08_wp ) *       va  * d_tmrt3 * pa2 +                    &
     1771       (  3.23926897e-09_wp ) *             d_tmrt4 * pa2
     1772
     1773    part_pa3 = ( -3.53874123e-02_wp ) *               pa3 +                    &
     1774       ( -2.21201190e-01_wp ) * ta  *                 pa3 +                    &
     1775       (  1.55126038e-02_wp ) * ta2 *                 pa3 +                    &
     1776       ( -2.63917279e-04_wp ) * ta3 *                 pa3 +                    &
     1777       (  4.53433455e-02_wp ) *       va  *           pa3 +                    &
     1778       ( -4.32943862e-03_wp ) * ta  * va  *           pa3 +                    &
     1779       (  1.45389826e-04_wp ) * ta2 * va  *           pa3 +                    &
     1780       (  2.17508610e-04_wp ) *       va2 *           pa3 +                    &
     1781       ( -6.66724702e-05_wp ) * ta  * va2 *           pa3 +                    &
     1782       (  3.33217140e-05_wp ) *       va3 *           pa3 +                    &
     1783       ( -2.26921615e-03_wp ) *             d_tmrt  * pa3 +                    &
     1784       (  3.80261982e-04_wp ) * ta  *       d_tmrt  * pa3 +                    &
     1785       ( -5.45314314e-09_wp ) * ta2 *       d_tmrt  * pa3 +                    &
     1786       ( -7.96355448e-04_wp ) *       va  * d_tmrt  * pa3 +                    &
     1787       (  2.53458034e-05_wp ) * ta  * va  * d_tmrt  * pa3 +                    &
     1788       ( -6.31223658e-06_wp ) *       va2 * d_tmrt  * pa3 +                    &
     1789       (  3.02122035e-04_wp ) *             d_tmrt2 * pa3 +                    &
     1790       ( -4.77403547e-06_wp ) * ta  *       d_tmrt2 * pa3 +                    &
     1791       (  1.73825715e-06_wp ) *       va  * d_tmrt2 * pa3 +                    &
     1792       ( -4.09087898e-07_wp ) *             d_tmrt3 * pa3
     1793
     1794    part_pa46 = (  6.14155345e-01_wp ) *              pa4 +                    &
     1795       ( -6.16755931e-02_wp ) * ta  *                 pa4 +                    &
     1796       (  1.33374846e-03_wp ) * ta2 *                 pa4 +                    &
     1797       (  3.55375387e-03_wp ) *       va  *           pa4 +                    &
     1798       ( -5.13027851e-04_wp ) * ta  * va *            pa4 +                    &
     1799       (  1.02449757e-04_wp ) *       va2 *           pa4 +                    &
     1800       ( -1.48526421e-03_wp ) *             d_tmrt  * pa4 +                    &
     1801       ( -4.11469183e-05_wp ) * ta  *       d_tmrt  * pa4 +                    &
     1802       ( -6.80434415e-06_wp ) *       va  * d_tmrt  * pa4 +                    &
     1803       ( -9.77675906e-06_wp ) *             d_tmrt2 * pa4 +                    &
     1804       (  8.82773108e-02_wp ) *                       pa5 +                    &
     1805       ( -3.01859306e-03_wp ) * ta  *                 pa5 +                    &
     1806       (  1.04452989e-03_wp ) *       va  *           pa5 +                    &
     1807       (  2.47090539e-04_wp ) *             d_tmrt  * pa5 +                    &
     1808       (  1.48348065e-03_wp ) *                       pa6
     1809!
    13441810!-- 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 
     1811    utci_ij = ta + part_ta + part_va + part_d_tmrt + part_pa + part_pa2 +      &
     1812        part_pa3 + part_pa46
     1813!
    15581814!-- Consider offset in result
    1559     utci = utci + offset
     1815    utci_ij = utci_ij + offset
    15601816
    15611817 END SUBROUTINE calculate_utci_static
     
    15701826!> Value of perct is the Perceived Temperature, degree centigrade
    15711827!------------------------------------------------------------------------------!
    1572  SUBROUTINE calculate_perct_static( ta, vp, ws, tmrt, pair, clo, perct )
     1828 SUBROUTINE calculate_perct_static( ta, vp, ws, tmrt, pair, clo, perct_ij )
    15731829
    15741830    IMPLICIT NONE
    1575 
     1831!
    15761832!-- Type of input of the argument list
    15771833    REAL(wp), INTENT ( IN )  :: ta   !< Local air temperature (degC)
     
    15801836    REAL(wp), INTENT ( IN )  :: ws   !< Local wind velocitry (m/s)
    15811837    REAL(wp), INTENT ( IN )  :: pair !< Local barometric air pressure (hPa)
    1582 
     1838!
    15831839!-- 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 
     1840    REAL(wp), INTENT ( OUT ) :: perct_ij  !< Perceived temperature (degC)
     1841    REAL(wp), INTENT ( OUT ) :: clo       !< Clothing index (dimensionless)
     1842!
    15871843!-- Parameters for standard "Klima-Michel"
    15881844    REAL(wp), PARAMETER :: eta = 0._wp  !< Mechanical work efficiency for walking on flat ground (compare to Fanger (1972) pp 24f)
    15891845    REAL(wp), PARAMETER :: actlev = 134.6862_wp !< Workload by activity per standardized surface (A_Du)
    1590 
     1846!
    15911847!-- Type of program variables
    15921848    REAL(wp), PARAMETER :: eps = 0.0005  !< Accuracy in clothing insulation (clo) for evaluation the root of Fanger's PMV (pmva=0)
     
    16141870
    16151871    LOGICAL :: sultrieness
    1616 
     1872!
    16171873!-- Initialise
    1618     perct = 9999.0_wp
     1874    perct_ij = bio_fill_value
    16191875
    16201876    nerr     = 0_iwp
    16211877    ncount   = 0_iwp
    16221878    sultrieness  = .FALSE.
     1879!
    16231880!-- Tresholds: clothing insulation (account for model inaccuracies)
     1881!
    16241882!   summer clothing
    16251883    sclo     = 0.44453_wp
     1884!
    16261885!   winter clothing
    16271886    wclo     = 1.76267_wp
    1628 
     1887!
    16291888!-- decision: firstly calculate for winter or summer clothing
    16301889    IF ( ta <= 10._wp ) THEN
    1631 
     1890!
    16321891!--    First guess: winter clothing insulation: cold stress
    16331892       clo = wclo
     
    16361895
    16371896       IF ( pmva > 0._wp ) THEN
    1638 
     1897!
    16391898!--       Case summer clothing insulation: heat load ?
    16401899          clo = sclo
     
    16431902          pmv_s = pmva
    16441903          IF ( pmva <= 0._wp ) THEN
     1904!
    16451905!--          Case: comfort achievable by varying clothing insulation
    1646 !--          Between winter and summer set values
     1906!            Between winter and summer set values
    16471907             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
    16481908                pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
     
    16621922       ENDIF
    16631923    ELSE
    1664 
     1924!
    16651925!--    First guess: summer clothing insulation: heat load
    16661926       clo = sclo
     
    16691929
    16701930       IF ( pmva < 0._wp ) THEN
    1671 
     1931!
    16721932!--       Case winter clothing insulation: cold stress ?
    16731933          clo = wclo
     
    16771937
    16781938          IF ( pmva >= 0._wp ) THEN
    1679 
     1939!
    16801940!--          Case: comfort achievable by varying clothing insulation
    16811941!            between winter and summer set values
     
    16981958
    16991959    ENDIF
    1700 
     1960!
    17011961!-- Determine perceived temperature by regression equation + adjustments
    17021962    pmvs = pmva
    1703     CALL perct_regression ( pmva, clo, perct )
    1704     ptc = perct
     1963    CALL perct_regression ( pmva, clo, perct_ij )
     1964    ptc = perct_ij
    17051965    IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
     1966!
    17061967!--    Adjust for cold conditions according to Gagge 1986
    17071968       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
     
    17121973          pmvs   = -0.11_wp
    17131974       ENDIF
    1714        CALL perct_regression ( pmvs, clo, perct )
     1975       CALL perct_regression ( pmvs, clo, perct_ij )
    17151976    ENDIF
    17161977!     clo_fanger = clo
    17171978    clon = clo
    1718     IF ( clo > 0.5_wp .AND. perct <= 8.73_wp ) THEN
     1979    IF ( clo > 0.5_wp .AND. perct_ij <= 8.73_wp ) THEN
     1980!
    17191981!--    Required clothing insulation (ireq) is exclusively defined for
    17201982!      operative temperatures (top) less 10 (C) for a
    17211983!      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 )
     1984       clon = ireq_neutral ( perct_ij, ireq_minimal, nerr )
    17231985       clo = clon
    17241986    ENDIF
     
    17271989    d_std = -99._wp
    17281990    IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
     1991!
    17291992!--    Adjust for warm/humid conditions according to Gagge 1986
    17301993       CALL saturation_vapor_pressure ( ta, svp_ta )
    17311994       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
    17321995       pmvs   = pmva + d_pmv
    1733        CALL perct_regression ( pmvs, clo, perct )
     1996       CALL perct_regression ( pmvs, clo, perct_ij )
    17341997       IF ( sult_lim < 99._wp ) THEN
    1735           IF ( (perct - ptc) > sult_lim ) sultrieness = .TRUE.
     1998          IF ( (perct_ij - ptc) > sult_lim ) sultrieness = .TRUE.
     1999!
    17362000!--       Set factor to threshold for sultriness
    17372001          IF ( dgtcstd /= 0_iwp ) THEN
    1738              d_std = ( ( perct - ptc ) - dgtcm ) / dgtcstd
     2002             d_std = ( ( perct_ij - ptc ) - dgtcm ) / dgtcstd
    17392003          ENDIF
    17402004       ENDIF
     
    17622026
    17632027    IF ( ta < 0._wp ) THEN
     2028!
    17642029!--    ta  < 0 (degC): saturation water vapour pressure over ice
    17652030       b = 17.84362_wp
    17662031       c = 245.425_wp
    17672032    ELSE
     2033!
    17682034!--    ta >= 0 (degC): saturation water vapour pressure over water
    17692035       b = 17.08085_wp
    17702036       c = 234.175_wp
    17712037    ENDIF
    1772 
     2038!
    17732039!-- Saturation water vapour pressure
    17742040    svp_ta = 6.1078_wp * EXP ( b * ta / ( c + ta ) )
     
    17882054
    17892055    IMPLICIT NONE
    1790 
     2056!
    17912057!-- Input variables of argument list:
    17922058    REAL(wp), INTENT ( IN )  :: ta     !< Ambient temperature (degC)
     
    18032069    REAL(wp), INTENT ( IN )  :: pmv_w    !< Fanger's PMV corresponding to wclo
    18042070    REAL(wp), INTENT ( IN )  :: pmv_s    !< Fanger's PMV corresponding to sclo
    1805 
     2071!
    18062072! Output variables of argument list:
    18072073    REAL(wp), INTENT ( OUT ) :: pmva     !< 0 (set to zero, because clo is evaluated for comfort)
     
    18142080!           nerr = -2: error = maximum iterations (max_iteration) exceeded
    18152081!           nerr = -3: error = root not bracketed between sclo and wclo
    1816 
     2082!
    18172083!-- Type of program variables
    18182084    INTEGER(iwp), PARAMETER  ::  max_iteration = 15_iwp       !< max number of iterations
     
    18312097    REAL(wp) ::  sroot       !< sqrt of PMV-guess
    18322098    INTEGER(iwp) ::  j       !< running index
    1833 
     2099!
    18342100!-- Initialise
    18352101    nerr    = 0_iwp
    1836 
     2102!
    18372103!-- Set pmva = 0 (comfort): Root of PMV depending on clothing insulation
     2104    x_ridder    = bio_fill_value
    18382105    pmva        = 0._wp
    18392106    clo_lower   = sclo
     
    18842151             y_lower = y_new
    18852152          ELSE
     2153!
    18862154!--          Never get here in x_ridder: singularity in y
    18872155             nerr       = -1_iwp
     
    18952163          ENDIF
    18962164       ENDDO
     2165!
    18972166!--    x_ridder exceed maximum iterations
    18982167       nerr       = -2_iwp
     
    19042173       x_ridder = clo_upper
    19052174    ELSE
     2175!
    19062176!--    x_ridder not bracketed by u_clo and o_clo
    19072177       nerr = -3_iwp
     
    19192189!> individual and reference environment.
    19202190!------------------------------------------------------------------------------!
    1921  SUBROUTINE perct_regression( pmv, clo, perct )
     2191 SUBROUTINE perct_regression( pmv, clo, perct_ij )
    19222192
    19232193    IMPLICIT NONE
     
    19262196    REAL(wp), INTENT ( IN ) ::  clo   !< clothing insulation index (clo)
    19272197
    1928     REAL(wp), INTENT ( OUT ) ::  perct   !< perct (degC) corresponding to given PMV / clo
     2198    REAL(wp), INTENT ( OUT ) ::  perct_ij   !< perct (degC) corresponding to given PMV / clo
    19292199
    19302200    IF ( pmv <= -0.11_wp ) THEN
    1931        perct = 5.805_wp + 12.6784_wp * pmv
     2201       perct_ij = 5.805_wp + 12.6784_wp * pmv
    19322202    ELSE
    19332203       IF ( pmv >= + 0.01_wp ) THEN
    1934           perct = 16.826_wp + 6.163_wp * pmv
     2204          perct_ij = 16.826_wp + 6.163_wp * pmv
    19352205       ELSE
    1936           perct = 21.258_wp - 9.558_wp * clo
     2206          perct_ij = 21.258_wp - 9.558_wp * clo
    19372207       ENDIF
    19382208    ENDIF
     
    19532223
    19542224    IMPLICIT NONE
    1955 
     2225!
    19562226!-- Input variables of argument list:
    19572227    REAL(wp), INTENT ( IN ) ::  ta       !< Ambient air temperature (degC)
     
    19632233    REAL(wp), INTENT ( IN ) ::  actlev   !< Individuals activity level per unit surface area (W/m2)
    19642234    REAL(wp), INTENT ( IN ) ::  eta      !< Individuals mechanical work efficiency (dimensionless)
    1965 
     2235!
    19662236!-- Output variables of argument list:
    19672237    REAL(wp), INTENT ( OUT ) ::  pmva    !< Actual Predicted Mean Vote (PMV,
     
    19692239!            (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta)
    19702240    REAL(wp), INTENT ( OUT ) ::  top     !< operative temperature (degC)
    1971 
     2241!
    19722242!-- Internal variables
    19732243    REAL(wp) ::  f_cl         !< Increase in surface due to clothing    (factor)
     
    19922262    REAL(wp) ::  z6           !< Heat loss through forced convection
    19932263    INTEGER(iwp) :: i         !< running index
    1994 
     2264!
    19952265!-- Clo must be > 0. to avoid div. by 0!
    19962266    clo = in_clo
    19972267    IF ( clo <= 0._wp ) clo = .001_wp
    1998 
     2268!
    19992269!-- f_cl = Increase in surface due to clothing
    20002270    f_cl = 1._wp + .15_wp * clo
    2001 
     2271!
    20022272!-- Case of free convection (ws < 0.1 m/s ) not considered
    20032273    ws = in_ws
     
    20052275       ws = .1_wp
    20062276    ENDIF
    2007 
     2277!
    20082278!-- Heat_convection = forced convection
    20092279    heat_convection = 12.1_wp * SQRT ( ws * pair / 1013.25_wp )
    2010 
     2280!
    20112281!-- Activity = inner heat produktion per standardized surface
    20122282    activity = actlev * ( 1._wp - eta )
    2013 
     2283!
    20142284!-- T_skin_aver = average skin temperature
    20152285    t_skin_aver = 35.7_wp - .0275_wp * activity
    2016 
     2286!
    20172287!-- Calculation of constants for evaluation below
    20182288    bc = .155_wp * clo * 3.96_wp * 10._wp**( -8 ) * f_cl
     
    20212291    dc = ( 1._wp + ec * cc ) / bc
    20222292    gc = ( t_skin_aver + bc * ( tmrt + degc_to_k )**4 + ec * cc * ta ) / bc
    2023 
     2293!
    20242294!-- Calculation of clothing surface temperature (t_clothing) based on
    20252295!   Newton-approximation with air temperature as initial guess
     
    20292299          * dc - gc ) / ( 4._wp * ( t_clothing + degc_to_k )**3 + dc )
    20302300    ENDDO
    2031 
     2301!
    20322302!-- Empiric factor for the adaption of the heat ballance equation
    20332303!   to the psycho-physical scale (Equ. 40 in FANGER)
    20342304    z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp )
    2035 
     2305!
    20362306!-- Water vapour diffution through the skin
    20372307    z2 = .31_wp * ( 57.3_wp - .07_wp * activity-pa )
    2038 
     2308!
    20392309!-- Sweat evaporation from the skin surface
    20402310    z3 = .42_wp * ( activity - 58._wp )
    2041 
     2311!
    20422312!-- Loss of latent heat through respiration
    20432313    z4 = .0017_wp * actlev * ( 58.7_wp - pa ) + .0014_wp * actlev *            &
    20442314      ( 34._wp - ta )
    2045 
     2315!
    20462316!-- Loss of radiational heat
    20472317    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + degc_to_k )**4 - ( tmrt +        &
     
    20522322       hr = 0._wp
    20532323    ENDIF
    2054 
     2324!
    20552325!-- Heat loss through forced convection cc*(t_clothing-TT)
    20562326    z6 = cc * ( t_clothing - ta )
    2057 
     2327!
    20582328!-- Predicted Mean Vote
    20592329    pmva = z1 * ( activity - z2 - z3 - z4 - z5 - z6 )
    2060 
     2330!
    20612331!-- Operative temperatur
    20622332    top = ( hr * tmrt + heat_convection * ta ) / ( hr + heat_convection )
     
    20732343
    20742344    IMPLICIT NONE
    2075 
     2345!
    20762346!-- Input variables of argument list:
    20772347    REAL(wp),     INTENT ( IN )  :: pmva     !< Actual Predicted Mean Vote (PMV) according to Fanger
     
    20812351    REAL(wp),     INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC) at screen level
    20822352    REAL(wp),     INTENT ( IN )  :: ws       !< Wind speed (m/s) 1 m above ground
    2083 
     2353!
    20842354!-- Output variables of argument list:
    20852355    INTEGER(iwp), INTENT ( OUT ) :: nerr     !< Error status / quality flag
     
    20882358!            -3 = rel. humidity set to 5 % or 95 %, respectively
    20892359!            -4 = deltapmv set to avoid pmvs < 0
    2090 
     2360!
    20912361!-- Internal variable types:
    20922362    REAL(wp) ::  pmv          !<
     
    21422412     +1.8686215_wp, +3.4260713_wp, +2.0116185_wp, -0.7777552_wp, -4.6715853_wp,&
    21432413     -7.7314281_wp, -11.7602578_wp, -23.5934198_wp /
    2144 
     2414!
    21452415!-- Test for compliance with regression range
    21462416    IF ( pmva < -1.0_wp .OR. pmva > 7.0_wp ) THEN
     
    21492419       nerr = 0_iwp
    21502420    ENDIF
    2151 
     2421!
    21522422!-- Initialise classic PMV
    21532423    pmv  = pmva
    2154 
     2424!
    21552425!-- Water vapour pressure of air
    21562426    p10  = 0.05_wp * svp_ta
     
    21612431       nerr = -3_iwp
    21622432       IF ( vp < p10 ) THEN
     2433!
    21632434!--       Due to conditions of regression: r.H. >= 5 %
    21642435          pa = p10
    21652436       ELSE
     2437!
    21662438!--       Due to conditions of regression: r.H. <= 95 %
    21672439          pa = p95
     
    21692441    ENDIF
    21702442    IF ( pa > 0._wp ) THEN
     2443!
    21712444!--    Natural logarithm of pa
    21722445       apa = LOG ( pa )
     
    21742447       apa = -5._wp
    21752448    ENDIF
    2176 
     2449!
    21772450!-- Ratio actual water vapour pressure to that of a r.H. of 50 %
    21782451    pa_p50   = 0.5_wp * svp_ta
     
    21842457       pa_p50 = 0._wp
    21852458    ENDIF
     2459!
    21862460!-- Square root of wind velocity
    21872461    IF ( ws >= 0._wp ) THEN
     
    21902464       sqvel = 0._wp
    21912465    ENDIF
     2466!
    21922467!-- Air temperature
    21932468!    ta = ta
    21942469!-- Difference mean radiation to air temperature
    21952470    dtmrt = tmrt - ta
    2196 
     2471!
    21972472!-- Select the valid regression coefficients
    21982473    nreg = INT ( pmv )
    21992474    IF ( nreg < 0_iwp ) THEN
     2475!
    22002476!--    value of the FUNCTION in the case pmv <= -1
    22012477       deltapmv = 0._wp
     
    22132489       ENDIF
    22142490    ENDIF
    2215 
     2491!
    22162492!-- Regression valid for 0. <= pmv <= 6.
    22172493    dpmv_1 =                                                                   &
     
    22392515          + aconst ( nreg + 1_iwp )
    22402516    ENDIF
    2241 
     2517!
    22422518!-- Calculate pmv modification
    22432519    deltapmv = ( 1._wp - gew ) * dpmv_1 + gew * dpmv_2
    22442520    pmvs = pmva + deltapmv
    22452521    IF ( ( pmvs ) < 0._wp ) THEN
     2522!
    22462523!--    Prevent negative pmv* due to problems with clothing insulation
    22472524       nerr = -4_iwp
    22482525       IF ( pmvs > -0.11_wp ) THEN
     2526!
    22492527!--       Threshold from FUNCTION perct_regression for winter clothing insulation
    22502528          deltapmv = deltapmv + 0.11_wp
    22512529       ELSE
     2530!
    22522531!--       Set pmvs to "0" for compliance with summer clothing insulation
    22532532          deltapmv = -1._wp * pmva
     
    22662545!> perct.
    22672546!------------------------------------------------------------------------------!
    2268  SUBROUTINE calc_sultr( perct, dperctm, dperctstd, sultr_res )
     2547 SUBROUTINE calc_sultr( perct_ij, dperctm, dperctstd, sultr_res )
    22692548
    22702549    IMPLICIT NONE
    2271 
     2550!
    22722551!-- Input of the argument list:
    2273     REAL(wp), INTENT ( IN )  ::  perct      !< Classical perceived temperature: Base is Fanger's PMV
    2274 
     2552    REAL(wp), INTENT ( IN )  ::  perct_ij      !< Classical perceived temperature: Base is Fanger's PMV
     2553!
    22752554!-- Additional output variables of argument list:
    22762555    REAL(wp), INTENT ( OUT ) ::  dperctm    !< Mean deviation perct (classical gt) to gt* (rational gt
     
    22792558!             determining the significance to perceive sultriness
    22802559    REAL(wp), INTENT ( OUT ) ::  sultr_res
    2281 
     2560!
    22822561!-- Types of coefficients mean deviation: third order polynomial
    22832562    REAL(wp), PARAMETER ::  dperctka = +7.5776086_wp
     
    22852564    REAL(wp), PARAMETER ::  dperctkc = +0.0213324_wp
    22862565    REAL(wp), PARAMETER ::  dperctkd = -0.00027797237_wp
    2287 
     2566!
    22882567!-- Types of coefficients mean deviation plus standard deviation
    22892568!   regression coefficients: third order polynomial
     
    22922571    REAL(wp), PARAMETER ::  dperctsc = -0.00054709752_wp
    22932572    REAL(wp), PARAMETER ::  dperctsd = +0.0000063714823_wp
    2294 
     2573!
    22952574!-- Factor to mean standard deviation defining SIGNificance for
    22962575!   sultriness
    22972576    REAL(wp), PARAMETER :: faktor = 1._wp
    2298 
     2577!
    22992578!-- Initialise
    23002579    sultr_res = +99._wp
     
    23022581    dperctstd = 999999._wp
    23032582
    2304     IF ( perct < 16.826_wp .OR. perct > 56._wp ) THEN
    2305 !--    Unallowed classical PMV/perct
     2583    IF ( perct_ij < 16.826_wp .OR. perct_ij > 56._wp ) THEN
     2584!
     2585!--    Unallowed value of classical perct!
    23062586       RETURN
    23072587    ENDIF
    2308 
     2588!
    23092589!-- Mean deviation dependent on perct
    2310     dperctm = dperctka + dperctkb * perct + dperctkc * perct**2._wp + dperctkd * perct**3._wp
    2311 
     2590    dperctm = dperctka + dperctkb * perct_ij + dperctkc * perct_ij**2._wp +    &
     2591       dperctkd * perct_ij**3._wp
     2592!
    23122593!-- Mean deviation plus its standard deviation
    2313     dperctstd = dperctsa + dperctsb * perct + dperctsc * perct**2._wp + dperctsd * perct**3._wp
    2314 
     2594    dperctstd = dperctsa + dperctsb * perct_ij + dperctsc * perct_ij**2._wp +  &
     2595       dperctsd * perct_ij**3._wp
     2596!
    23152597!-- Value of the FUNCTION
    23162598    sultr_res = dperctm + faktor * dperctstd
     
    23302612
    23312613    IMPLICIT NONE
    2332 
     2614!
    23332615!-- Type of input arguments
    23342616    REAL(wp), INTENT ( IN ) ::  pmva   !< Fanger's classical predicted mean vote
     
    23362618    REAL(wp), INTENT ( IN ) ::  ws     !< Relative wind velocity 1 m above ground (m/s)
    23372619    REAL(wp), INTENT ( IN ) ::  tmrt   !< Mean radiant temperature (degC)
    2338 
     2620!
    23392621!-- Type of output argument
    23402622    INTEGER(iwp), INTENT ( OUT ) ::  nerr !< Error indicator: 0 = o.k., +1 = denominator for intersection = 0
    23412623    REAL(wp),     INTENT ( OUT ) ::  dpmv_cold_res    !< Increment to adjust pmva according to the results of Gagge's
    23422624!               2 node model depending on the input
    2343 
     2625!
    23442626!-- Type of program variables
    23452627    REAL(wp) ::  delta_cold(3)
     
    23522634    REAL(wp) ::  sqrt_ws
    23532635    INTEGER(iwp) ::  i
    2354     INTEGER(iwp) ::  j
     2636!     INTEGER(iwp) ::  j
    23552637    INTEGER(iwp) ::  i_bin
    2356 
     2638!
    23572639!-- Coefficient of the 3 regression lines
    2358     !     1:const   2:*pmvc  3:*ta      4:*sqrt_ws  5:*dtmrt
     2640!     1:const   2:*pmvc  3:*ta      4:*sqrt_ws  5:*dtmrt
    23592641    coeff(1,1:5) =                                                             &
    23602642       (/ +0.161_wp, +0.130_wp, -1.125E-03_wp, +1.106E-03_wp, -4.570E-04_wp /)
     
    23632645    coeff(3,1:5) =                                                             &
    23642646       (/ +0.05761_wp, +0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp /)
    2365 
     2647!
    23662648!-- Initialise
    23672649    nerr       = 0_iwp
     
    23842666
    23852667    DO i = 1, 3
     2668!
    23862669!--    Regression constant for given meteorological variables
    2387        reg_a(i) = coeff(i, 1) + coeff(i, 3) * ta + coeff(i, 4) *             &
     2670       reg_a(i) = coeff(i, 1) + coeff(i, 3) * ta + coeff(i, 4) *               &
    23882671                  sqrt_ws + coeff(i,5)*dtmrt
    23892672       delta_cold(i) = reg_a(i) + coeff(i, 2) * pmvc
    23902673    ENDDO
    2391 
     2674!
    23922675!-- Intersection points of regression lines in terms of Fanger's PMV
    23932676    DO i = 1, 2
     
    24082691       ENDIF
    24092692    ENDDO
     2693!
    24102694!-- Adjust to operative temperature scaled according
    24112695!   to classical PMV (Fanger)
     
    24342718    INTEGER(iwp)  ::  i, i_bin
    24352719
    2436 !                             range_1     range_2     range_3
     2720!                                   range_1        range_2        range_3
    24372721    DATA (coef(i, 0), i = 1, n_bin) /0.0941540_wp, -0.1506620_wp, -0.0871439_wp/
    24382722    DATA (coef(i, 1), i = 1, n_bin) /0.0783162_wp, -1.0612651_wp,  0.1695040_wp/
     
    24672751!> for perct < 10 (degC)
    24682752!------------------------------------------------------------------------------!
    2469  REAL(wp) FUNCTION ireq_neutral( perct, ireq_minimal, nerr )
     2753 REAL(wp) FUNCTION ireq_neutral( perct_ij, ireq_minimal, nerr )
    24702754
    24712755    IMPLICIT NONE
    2472 
     2756!
    24732757!-- Type declaration of arguments
    2474     REAL(wp),     INTENT ( IN )  ::  perct
     2758    REAL(wp),     INTENT ( IN )  ::  perct_ij
    24752759    REAL(wp),     INTENT ( OUT ) ::  ireq_minimal
    24762760    INTEGER(iwp), INTENT ( OUT ) ::  nerr
    2477 
     2761!
    24782762!-- Type declaration for internal varables
    24792763    REAL(wp)                     ::  top02
    2480 
     2764!
    24812765!-- Initialise
    24822766    nerr = 0_iwp
    2483 
     2767!
    24842768!-- 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 
     2769    top02 = 1.8788_wp + 0.9296_wp * perct_ij
     2770!
    24872771!-- IREQ neutral conditions (thermal comfort)
    24882772    ireq_neutral = 1.62_wp - 0.0564_wp * top02
    2489 
     2773!
    24902774!-- Regression only defined for perct <= 10 (degC)
    24912775    IF ( ireq_neutral < 0.5_wp ) THEN
     
    24952779       ireq_neutral = 0.5_wp
    24962780    ENDIF
    2497 
     2781!
    24982782!-- Minimal required clothing insulation: maximal acceptable body cooling
    24992783    ireq_minimal = 1.26_wp - 0.0588_wp * top02
     
    25232807
    25242808    height = height_cm * 100._wp
    2525 
     2809!
    25262810!-- According to Gehan-George, for children
    25272811    IF ( age < 19_iwp ) THEN
     
    25332817       RETURN
    25342818    END IF
    2535 
     2819!
    25362820!-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if
    25372821!   height and weight be known. In: Arch. Int. Med.. 17, 1916, S. 863?871.
     
    25982882
    25992883 SUBROUTINE ipt_init (age, weight, height, sex, work, actlev, clo,             &
    2600      ta, vp, ws, tmrt, pair, dt, storage, t_clothing,                         &
     2884     ta, vp, ws, tmrt, pair, dt, storage, t_clothing,                          &
    26012885     ipt )
    26022886
    26032887    IMPLICIT NONE
    2604 
     2888!
    26052889!-- Input parameters
    26062890    REAL(wp), INTENT(in) ::  age        !< Persons age          (years)
     
    26152899    REAL(wp), INTENT(in) ::  dt         !< Timestep             (s)
    26162900    INTEGER(iwp), INTENT(in)  :: sex    !< Persons sex (1 = male, 2 = female)
    2617 
     2901!
    26182902!-- Output parameters
    26192903    REAL(wp), INTENT(out) ::  actlev
     
    26222906    REAL(wp), INTENT(out) ::  t_clothing
    26232907    REAL(wp), INTENT(out) ::  ipt
    2624 
     2908!
    26252909!-- Internal variables
    26262910    REAL(wp), PARAMETER :: eps = 0.0005_wp
     
    26442928    REAL(wp) ::  top
    26452929    REAL(wp) ::  a_surf
    2646     REAL(wp) ::  acti
     2930!     REAL(wp) ::  acti
    26472931    INTEGER(iwp) ::  ncount
    26482932    INTEGER(iwp) ::  nerr_cold
     
    26532937    storage = 0._wp
    26542938    CALL persdat ( age, weight, height, sex, work, a_surf, actlev )
    2655 
     2939!
    26562940!-- Initialise
    2657     t_clothing = -999.0_wp
    2658     ipt        = -999.0_wp
     2941    t_clothing = bio_fill_value
     2942    ipt        = bio_fill_value
    26592943    nerr       = 0_wp
    26602944    ncount     = 0_wp
    26612945    sultrieness    = .FALSE.
     2946!
    26622947!-- Tresholds: clothing insulation (account for model inaccuracies)
    26632948!-- Summer clothing
     
    26652950!-- Winter clothing
    26662951    wclo     = 1.76267_wp
    2667 
     2952!
    26682953!-- Decision: firstly calculate for winter or summer clothing
    26692954    IF ( ta <= 10._wp ) THEN
    2670 
     2955!
    26712956!--    First guess: winter clothing insulation: cold stress
    26722957       clo = wclo
    2673        t_clothing = -999.0_wp  ! force initial run
     2958       t_clothing = bio_fill_value  ! force initial run
    26742959       CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,         &
    26752960          t_clothing, storage, dt, pmva )
     
    26772962
    26782963       IF ( pmva > 0._wp ) THEN
    2679 
     2964!
    26802965!--       Case summer clothing insulation: heat load ?           
    26812966          clo = sclo
    2682           t_clothing = -999.0_wp  ! force initial run
     2967          t_clothing = bio_fill_value  ! force initial run
    26832968          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
    26842969             t_clothing, storage, dt, pmva )
    26852970          pmv_s = pmva
    26862971          IF ( pmva <= 0._wp ) THEN
     2972!
    26872973!--          Case: comfort achievable by varying clothing insulation
    26882974!            between winter and summer set values
     
    26952981          ELSE IF ( pmva > 0.06_wp ) THEN
    26962982             clo = 0.5_wp
    2697              t_clothing = -999.0_wp
     2983             t_clothing = bio_fill_value
    26982984             CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,   &
    26992985                t_clothing, storage, dt, pmva )
     
    27012987       ELSE IF ( pmva < -0.11_wp ) THEN
    27022988          clo = 1.75_wp
    2703           t_clothing = -999.0_wp
     2989          t_clothing = bio_fill_value
    27042990          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
    27052991             t_clothing, storage, dt, pmva )
     
    27072993
    27082994    ELSE
    2709 
     2995!
    27102996!--    First guess: summer clothing insulation: heat load
    27112997       clo = sclo
    2712        t_clothing = -999.0_wp
     2998       t_clothing = bio_fill_value
    27132999       CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,         &
    27143000          t_clothing, storage, dt, pmva )
     
    27163002
    27173003       IF ( pmva < 0._wp ) THEN
    2718 
     3004!
    27193005!--       Case winter clothing insulation: cold stress ?
    27203006          clo = wclo
    2721           t_clothing = -999.0_wp
     3007          t_clothing = bio_fill_value
    27223008          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
    27233009             t_clothing, storage, dt, pmva )
     
    27253011
    27263012          IF ( pmva >= 0._wp ) THEN
    2727 
     3013!
    27283014!--          Case: comfort achievable by varying clothing insulation
    27293015!            between winter and summer set values
     
    27363022          ELSE IF ( pmva < -0.11_wp ) THEN
    27373023             clo = 1.75_wp
    2738              t_clothing = -999.0_wp
     3024             t_clothing = bio_fill_value
    27393025             CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,   &
    27403026                t_clothing, storage, dt, pmva )
     
    27423028       ELSE IF ( pmva > 0.06_wp ) THEN
    27433029          clo = 0.5_wp
    2744           t_clothing = -999.0_wp
     3030          t_clothing = bio_fill_value
    27453031          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,      &
    27463032             t_clothing, storage, dt, pmva )
     
    27483034
    27493035    ENDIF
    2750 
     3036!
    27513037!-- Determine perceived temperature by regression equation + adjustments
    27523038    pmvs = pmva
     
    27543040    ptc = ipt
    27553041    IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
     3042!
    27563043!--    Adjust for cold conditions according to Gagge 1986
    27573044       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
     
    27673054    clon = clo
    27683055    IF ( clo > 0.5_wp .AND. ipt <= 8.73_wp ) THEN
     3056!
    27693057!--    Required clothing insulation (ireq) is exclusively defined for
    27703058!      operative temperatures (top) less 10 (C) for a
     
    27773065    d_std      = -99._wp
    27783066    IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
     3067!
    27793068!--    Adjust for warm/humid conditions according to Gagge 1986
    27803069       CALL saturation_vapor_pressure ( ta, svp_ta )
     
    28083097
    28093098    IMPLICIT NONE
    2810 
     3099!
    28113100!-- Type of input of the argument list
    28123101    REAL(wp), INTENT ( IN )  ::  ta      !< Air temperature             (°C)
     
    28193108    REAL(wp), INTENT ( IN )  ::  actlev  !< Internal heat production    (W)
    28203109    REAL(wp), INTENT ( IN )  ::  work    !< Mechanical work load        (W)
    2821 
     3110!
    28223111!-- In and output parameters
    28233112    REAL(wp), INTENT (INOUT) ::  storage     !< Heat storage            (W)
    28243113    REAL(wp), INTENT (INOUT) ::  t_clothing  !< Clothig temperature     (°C)
    2825 
     3114!
    28263115!-- Type of output of the argument list
    28273116    REAL(wp), INTENT ( OUT ) ::  ipt  !< Instationary perceived temperature (°C)
    2828 
     3117!
    28293118!-- Type of internal variables
    28303119    REAL(wp) ::  d_pmv
     
    28413130
    28423131    LOGICAL ::  sultrieness
    2843 
     3132!
    28443133!-- Initialise
    2845     ipt = -999.0_wp
     3134    ipt = bio_fill_value
    28463135
    28473136    nerr     = 0_iwp
    28483137    sultrieness  = .FALSE.
    2849 
     3138!
    28503139!-- Determine pmv_adjusted for current conditions
    28513140    CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work,            &
    28523141       t_clothing, storage, dt, pmva )
    2853 
     3142!
    28543143!-- Determine perceived temperature by regression equation + adjustments
    28553144    CALL perct_regression ( pmva, clo, ipt )
    28563145
    28573146    IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
     3147!
    28583148!--    Adjust for cold conditions according to Gagge 1986
    28593149       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
     
    28663156       CALL perct_regression ( pmvs, clo, ipt )
    28673157    ENDIF
    2868 
     3158!
    28693159!-- Consider sultriness if appropriate
    28703160    ptc = ipt
     
    28733163    d_std      = -99._wp
    28743164    IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
     3165!
    28753166!--    Adjust for warm/humid conditions according to Gagge 1986
    28763167       CALL saturation_vapor_pressure ( ta, svp_ta )
     
    28983189
    28993190    IMPLICIT NONE
    2900 
     3191!
    29013192!--  Input argument types
    29023193    REAL(wp), INTENT ( IN )  ::  ta       !< Air temperature          (°C)
     
    29093200    REAL(wp), INTENT ( IN )  ::  activity !< Work load                (W/m²)
    29103201    REAL(wp), INTENT ( IN )  ::  in_clo   !< Clothing index (clo)     (no dim)
    2911 
     3202!
    29123203!-- Output argument types
    29133204    REAL(wp), INTENT ( OUT ) ::  pmva  !< actual Predicted Mean Vote (no dim)
     
    29153206    REAL(wp), INTENT (INOUT) ::  s  !< storage var. of energy balance (W/m2)
    29163207    REAL(wp), INTENT (INOUT) ::  t_cloth  !< clothing temperature (°C)
    2917 
     3208!
    29183209!-- Internal variables
    29193210    REAL(wp), PARAMETER  ::  time_equil = 7200._wp
     
    29283219    REAL(wp) ::  gc           !< preliminary result storage
    29293220    REAL(wp) ::  t_clothing   !< clothing temperature                   (°C)
    2930     REAL(wp) ::  hr           !< radiational heat resistence
     3221!     REAL(wp) ::  hr           !< radiational heat resistence
    29313222    REAL(wp) ::  clo          !< clothing insulation index              (clo)
    29323223    REAL(wp) ::  ws           !< wind speed                             (m/s)
     
    29463237    INTEGER(iwp) ::  niter  !< Running index
    29473238
    2948 
     3239!
    29493240!-- Clo must be > 0. to avoid div. by 0!
    29503241    clo = in_clo
    29513242    IF ( clo < 001._wp ) clo = .001_wp
    2952 
     3243!
    29533244!-- Increase in surface due to clothing
    29543245    f_cl = 1._wp + .15_wp * clo
    2955 
     3246!
    29563247!-- Case of free convection (ws < 0.1 m/s ) not considered
    29573248    ws = in_ws
     
    29593250       ws = .1_wp
    29603251    ENDIF
    2961 
     3252!
    29623253!-- Heat_convection = forced convection
    29633254    heat_convection = 12.1_wp * SQRT ( ws * pair / 1013.25_wp )
    2964 
     3255!
    29653256!-- Average skin temperature
    29663257    t_skin_aver = 35.7_wp - .0275_wp * activity
    2967 
     3258!
    29683259!-- Calculation of constants for evaluation below
    29693260    bc = .155_wp * clo * 3.96_wp * 10._wp**( -8._wp ) * f_cl
     
    29723263    dc = ( 1._wp + ec * cc ) / bc
    29733264    gc = ( t_skin_aver + bc * ( tmrt + 273.2_wp )**4._wp + ec * cc * ta ) / bc
    2974 
     3265!
    29753266!-- Calculation of clothing surface temperature (t_clothing) based on
    29763267!   newton-approximation with air temperature as initial guess
    2977     niter = dt
     3268    niter = INT( dt * 10._wp, KIND=iwp )
     3269    IF ( niter < 1 ) niter = 1_iwp
    29783270    adjustrate = 1._wp - EXP ( -1._wp * ( 10._wp / time_equil ) * dt )
    29793271    IF ( adjustrate >= 1._wp ) adjustrate = 1._wp
     
    29823274
    29833275    IF ( t_cloth <= -998.0_wp ) THEN  ! If initial run
    2984        niter = 3_wp
     3276       niter = 3_iwp
    29853277       adjustrate = 1._wp
    29863278       adjustrate_cloth = 1._wp
     
    29933285          dc - gc ) / ( 4._wp * ( t_clothing + 273.2_wp )**3._wp + dc )
    29943286    ENDDO
    2995 
     3287!
    29963288!-- Empiric factor for the adaption of the heat ballance equation
    29973289!   to the psycho-physical scale (Equ. 40 in FANGER)
    29983290    z1 = ( .303_wp * EXP ( -.036_wp * actlev ) + .0275_wp )
    2999 
     3291!
    30003292!-- Water vapour diffution through the skin
    30013293    z2 = .31_wp * ( 57.3_wp - .07_wp * activity-pa )
    3002 
     3294!
    30033295!-- Sweat evaporation from the skin surface
    30043296    z3 = .42_wp * ( activity - 58._wp )
    3005 
     3297!
    30063298!-- Loss of latent heat through respiration
    30073299    z4 = .0017_wp * actlev * ( 58.7_wp - pa ) + .0014_wp * actlev *            &
    30083300      ( 34._wp - ta )
    3009 
     3301!
    30103302!-- Loss of radiational heat
    30113303    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + 273.2_wp )**4 - ( tmrt +         &
    30123304       273.2_wp )**4 )
    3013 
     3305!
    30143306!-- Heat loss through forced convection
    30153307    z6 = cc * ( t_clothing - ta )
    3016 
     3308!
    30173309!-- Write together as energy ballance
    30183310    en = activity - z2 - z3 - z4 - z5 - z6
    3019 
     3311!
    30203312!-- Manage storage
    30213313    d_s = adjustrate * en + ( 1._wp - adjustrate ) * s
    3022 
     3314!
    30233315!-- Predicted Mean Vote
    30243316    pmva = z1 * d_s
    3025 
     3317!
    30263318!-- Update storage
    30273319    s = d_s
     
    30413333!------------------------------------------------------------------------------!
    30423334
    3043  SUBROUTINE calculate_pet_static( ta, vpa, v, tmrt, pair, pet )
     3335 SUBROUTINE calculate_pet_static( ta, vpa, v, tmrt, pair, pet_ij )
    30443336
    30453337    IMPLICIT NONE
    3046 
     3338!
    30473339!-- Input arguments:
    30483340    REAL(wp), INTENT( IN ) ::  ta    !< Air temperature             (°C)
     
    30513343    REAL(wp), INTENT( IN ) ::  vpa   !< Vapor pressure              (hPa)
    30523344    REAL(wp), INTENT( IN ) ::  pair  !< Air pressure                (hPa)
    3053 
     3345!
    30543346!-- Output arguments:
    3055     REAL(wp), INTENT ( OUT ) ::  pet  !< PET                         (°C)
    3056 
     3347    REAL(wp), INTENT ( OUT ) ::  pet_ij  !< PET                     (°C)
     3348!
    30573349!-- Internal variables:
    30583350    REAL(wp) ::  acl        !< clothing area                        (m²)
     
    30713363    REAL(wp) ::  tcl        !< Clothing temperature                 (°C)
    30723364    REAL(wp) ::  wetsk      !< Fraction of wet skin                 (factor)
    3073 
     3365!
    30743366!-- Variables:
    30753367    REAL(wp) :: int_heat    !< Internal heat        (W)
    3076 
     3368!
    30773369!-- MEMI configuration
    30783370    REAL(wp) :: age         !< Persons age          (a)
     
    30853377!     INTEGER(iwp) :: pos     !< Posture: 1 = standing, 2 = sitting
    30863378!     INTEGER(iwp) :: sex     !< Sex: 1 = male, 2 = female
    3087 
     3379!
    30883380!-- Configuration, keep standard parameters!
    30893381    age   = 35._wp
     
    30943386    clo   =  0.9_wp
    30953387    fcl   =  1.15_wp
    3096 
     3388!
    30973389!-- Call subfunctions
    30983390    CALL in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta,    &
     
    31033395            vpts, wetsk )
    31043396
    3105     CALL pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair, rdcl,&
    3106             rdsk, rtv, ta, tcl, tsk, pet, vpts, wetsk )
     3397    CALL pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair,      &
     3398            rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk )
    31073399
    31083400
     
    31173409 SUBROUTINE in_body ( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, &
    31183410    vpa, work )
    3119 
     3411!
    31203412!-- Input arguments:
    31213413    REAL(wp), INTENT( IN )  ::  pair      !< air pressure             (hPa)
     
    31273419    REAL(wp), INTENT( IN )  ::  work      !< Work load                (W)
    31283420    REAL(wp), INTENT( IN )  ::  eta       !< Work efficiency      (dimensionless)
    3129 
     3421!
    31303422!-- Output arguments:
    31313423    REAL(wp), INTENT( OUT ) ::  ere       !< energy ballance          (W)
     
    31333425    REAL(wp), INTENT( OUT ) ::  int_heat  !< internal heat production (W)
    31343426    REAL(wp), INTENT( OUT ) ::  rtv       !< respiratory volume
    3135 
     3427!
    31363428!-- Constants:
    31373429!     REAL(wp), PARAMETER :: cair = 1010._wp        !< replaced by c_p
    31383430!     REAL(wp), PARAMETER :: evap = 2.42_wp * 10._wp **6._wp  !< replaced by l_v
    3139 
     3431!
    31403432!-- Internal variables:
    31413433    REAL(wp) ::  eres                     !< Sensible respiratory heat flux (W)
     
    31513443
    31523444    int_heat = met * (1._wp - eta)
    3153 
     3445!
    31543446!-- Sensible respiration energy
    31553447    tex  = 0.47_wp * ta + 21.0_wp
    31563448    rtv  = 1.44_wp * 10._wp ** (-6._wp) * met
    31573449    eres = c_p * (ta - tex) * rtv
    3158 
     3450!
    31593451!-- Latent respiration energy
    31603452    vpex = 6.11_wp * 10._wp ** ( 7.45_wp * tex / ( 235._wp + tex ) )
    31613453    erel = 0.623_wp * l_v / pair * ( vpa - vpex ) * rtv
    3162 
     3454!
    31633455!-- Sum of the results
    31643456    ere = eres + erel
     
    31763468        vpts, wetsk )
    31773469
    3178 
     3470!
    31793471!-- Input arguments:
    31803472    REAL(wp), INTENT( IN )  ::  ere    !< Energy ballance          (W)
     
    31903482    REAL(wp), INTENT( IN )  ::  clo    !< clothing insulation      (clo)
    31913483    REAL(wp), INTENT( IN )  ::  fcl    !< factor for surface area increase by clothing
    3192 
     3484!
    31933485!-- Output arguments:
    31943486    REAL(wp), INTENT( OUT ) ::  acl    !< Clothing surface area        (m²)
     
    32043496    REAL(wp), INTENT( OUT ) ::  vpts   !< Sat. vapor pressure over skin (hPa)
    32053497    REAL(wp), INTENT( OUT ) ::  wetsk  !< Fraction of wet skin (dimensionless)
    3206 
     3498!
    32073499!-- Cconstants:
    32083500!     REAL(wp), PARAMETER :: cair = 1010._wp      !< replaced by c_p
     
    32143506    REAL(wp), PARAMETER :: po   = 1013.25_wp      !< Air pressure at sea level (hPa)
    32153507    REAL(wp), PARAMETER :: rob  =    1.06_wp      !<
    3216 
     3508!
    32173509!-- Internal variables
    32183510    REAL(wp) ::  c(0:10)        !< Core temperature array           (°C)
     
    32393531    REAL(wp) ::  rsum           !< Radiational loss or gain         (W/m²)
    32403532    REAL(wp) ::  sw             !<
    3241     REAL(wp) ::  swf            !<
     3533!     REAL(wp) ::  swf            !< female factor, currently unused
    32423534    REAL(wp) ::  swm            !<
    32433535    REAL(wp) ::  tbody          !<
     
    32793571
    32803572    di   = r2 - r1
    3281 
     3573!
    32823574!-- Skin temperatur
    32833575    DO j = 1, 7
     
    32953587          htcl  = 6.28_wp * ht * y * di / ( rcl * LOG( r2 / r1 ) * acl )
    32963588          tsk   = 1._wp / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl
    3297 
     3589!
    32983590!--       Radiation saldo
    32993591          aeff  = adu * feff
     
    33033595            ( ( tmrt + degc_to_k )** 4._wp - ( tcl + degc_to_k )** 4._wp )
    33043596          rsum  = rbare + rclo
    3305 
     3597!
    33063598!--       Convection
    33073599          cbare = hc * ( ta - tsk ) * adu * ( 1._wp - facl )
    33083600          cclo  = hc * ( ta - tcl ) * acl
    33093601          csum  = cbare + cclo
    3310 
     3602!
    33113603!--       Core temperature
    33123604          c(0)  = int_heat + ere
     
    33373629             tcore(4) = c(0) / ( 5.28_wp * adu + c(1) * 1._wp / 40._wp ) + tsk
    33383630          END IF
    3339 
     3631!
    33403632!--       Transpiration
    33413633          tbody = 0.1_wp * tsk + 0.9_wp * tcore(j)
     
    33603652          IF ( eswdif > 0._wp ) esw = eswphy   !< Limit is sweat production
    33613653          IF ( esw  > 0._wp )   esw = 0._wp    !< Sweat can't be evaporated, no more cooling effect
    3362 
     3654!
    33633655!--       Diffusion
    33643656          rdsk = 0.79_wp * 10._wp ** 7._wp
     
    33663658          ed   = l_v / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( vpa -     &
    33673659             vpts )
    3368 
     3660!
    33693661!--       Max vb
    33703662          vb1 = 34._wp - tsk
     
    33743666          IF ( vb1 < 0._wp ) vb1 = 0._wp
    33753667          vb = ( 6.3_wp + 75._wp * vb2 ) / ( 1._wp + 0.5_wp * vb1 )
    3376 
     3668!
    33773669!--       Energy ballence
    33783670          enbal = int_heat + ed + ere + esw + csum + rsum + food
    3379 
     3671!
    33803672!--       Clothing temperature
    33813673          xx = 0.001_wp
     
    34293721       IF ( ( j == 4_iwp ) .AND. ( vb < 89._wp ) ) CYCLE
    34303722       IF ( vb > 90._wp ) vb = 90._wp
    3431 
     3723!
    34323724!--    Loses by water
    34333725       ws = sw * 3600._wp * 1000._wp
     
    34483740!------------------------------------------------------------------------------!
    34493741 SUBROUTINE pet_iteration ( acl, adu, aeff, esw, facl, feff, int_heat, pair,   &
    3450         rdcl, rdsk, rtv, ta, tcl, tsk, pet, vpts, wetsk )
    3451 
     3742        rdcl, rdsk, rtv, ta, tcl, tsk, pet_ij, vpts, wetsk )
     3743!
    34523744!-- Input arguments:
    34533745    REAL(wp), INTENT( IN ) ::  acl   !< clothing surface area        (m²)
     
    34663758    REAL(wp), INTENT( IN ) ::  vpts  !< sat. vapor pressure over skin (hPa)
    34673759    REAL(wp), INTENT( IN ) ::  wetsk !< fraction of wet skin (dimensionless)
    3468 
     3760!
    34693761!-- Output arguments:
    3470     REAL(wp), INTENT( OUT ) ::  aeff  !< effective surface area       (m²)
    3471     REAL(wp), INTENT( OUT ) ::  pet   !< PET                          (°C)
    3472 
     3762    REAL(wp), INTENT( OUT ) ::  aeff     !< effective surface area       (m²)
     3763    REAL(wp), INTENT( OUT ) ::  pet_ij   !< PET                          (°C)
     3764!
    34733765!-- Cconstants:
    3474 !     REAL(wp), PARAMETER :: cair = 1010._wp        !< replaced by c_p
    34753766    REAL(wp), PARAMETER :: emcl =    0.95_wp      !< Longwave emission coef. of cloth
    34763767    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
    34783768    REAL(wp), PARAMETER :: po   = 1013.25_wp      !< Air pressure at sea level (hPa)
    3479 
     3769!
    34803770!-- Internal variables
    34813771    REAL ( wp ) ::  cbare             !< Convection through bare skin
     
    34993789    INTEGER ( iwp ) ::  i             !< running index
    35003790
    3501     pet = ta
     3791    pet_ij = ta
    35023792    enbal2 = 0._wp
    35033793
     
    35063796          hc = 2.67_wp + 6.5_wp * 0.1_wp ** 0.67_wp
    35073797          hc = hc * ( pair / po ) ** 0.55_wp
    3508 
     3798!
    35093799!--       Radiation
    35103800          aeff  = adu * feff
    35113801          rbare = aeff * ( 1._wp - facl ) * emsk * sigma_sb *                  &
    3512               ( ( pet + degc_to_k ) ** 4._wp - ( tsk + degc_to_k ) ** 4._wp )
     3802              ( ( pet_ij + degc_to_k ) ** 4._wp - ( tsk + degc_to_k ) ** 4._wp )
    35133803          rclo  = feff * acl * emcl * sigma_sb *                               &
    3514               ( ( pet + degc_to_k ) ** 4._wp - ( tcl + degc_to_k ) ** 4._wp )
     3804              ( ( pet_ij + degc_to_k ) ** 4._wp - ( tcl + degc_to_k ) ** 4._wp )
    35153805          rsum  = rbare + rclo
    3516 
     3806!
    35173807!--       Covection
    3518           cbare = hc * ( pet - tsk ) * adu * ( 1._wp - facl )
    3519           cclo  = hc * ( pet - tcl ) * acl
     3808          cbare = hc * ( pet_ij - tsk ) * adu * ( 1._wp - facl )
     3809          cclo  = hc * ( pet_ij - tcl ) * acl
    35203810          csum  = cbare + cclo
    3521 
     3811!
    35223812!--       Diffusion
    35233813          ed = l_v / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( 12._wp -    &
    35243814             vpts )
    3525 
     3815!
    35263816!--       Respiration
    3527           tex  = 0.47_wp * pet + 21._wp
    3528           eres = c_p * ( pet - tex ) * rtv
     3817          tex  = 0.47_wp * pet_ij + 21._wp
     3818          eres = c_p * ( pet_ij - tex ) * rtv
    35293819          vpex = 6.11_wp * 10._wp ** ( 7.45_wp * tex / ( 235._wp + tex ) )
    35303820          erel = 0.623_wp * l_v / pair * ( 12._wp - vpex ) * rtv
    35313821          ere  = eres + erel
    3532 
     3822!
    35333823!--       Energy ballance
    35343824          enbal = int_heat + ed + ere + esw + csum + rsum
    3535 
     3825!
    35363826!--       Iteration concerning ta
     3827          xx = 0.001_wp
    35373828          IF ( count1 == 0_iwp )  xx = 1._wp
    35383829          IF ( count1 == 1_iwp )  xx = 0.1_wp
    35393830          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
     3831!           IF ( count1 == 3_iwp )  xx = 0.001_wp
     3832          IF ( enbal > 0._wp )  pet_ij = pet_ij - xx
     3833          IF ( enbal < 0._wp )  pet_ij = pet_ij + xx
    35433834          IF ( ( enbal <= 0._wp ) .AND. ( enbal2 > 0._wp ) ) EXIT
    35443835          IF ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) EXIT
     
    35493840 END SUBROUTINE pet_iteration
    35503841
     3842!
     3843!-- UVEM specific subroutines
     3844
     3845!---------------------------------------------------------------------------------------------------------------------!
     3846! Description:
     3847! ------------
     3848!> Module-specific routine for new module
     3849!---------------------------------------------------------------------------------------------------------------------!
     3850 SUBROUTINE uvem_solar_position
     3851   
     3852    USE date_and_time_mod,                                                                                            &
     3853       ONLY:  calc_date_and_time, day_of_year, time_utc
     3854   
     3855    USE control_parameters,                                                                                           &
     3856       ONLY:  latitude, longitude   
     3857
     3858    IMPLICIT NONE
     3859   
     3860   
     3861    REAL(wp) ::  alpha       = 0.0_wp   !< solar azimuth angle in radiant   
     3862    REAL(wp) ::  doy_r       = 0.0_wp   !< real format of day_of_year           
     3863    REAL(wp) ::  declination = 0.0_wp   !< declination
     3864    REAL(wp) ::  dtor        = 0.0_wp   !< factor to convert degree to radiant
     3865    REAL(wp) ::  js          = 0.0_wp   !< parameter for solar position calculation
     3866    REAL(wp) ::  lat         = 52.39_wp !< latitude
     3867    REAL(wp) ::  lon         = 9.7_wp   !< longitude       
     3868    REAL(wp) ::  thetar      = 0.0_wp   !< angle for solar zenith angle calculation
     3869    REAL(wp) ::  thetasr     = 0.0_wp   !< angle for solar azimuth angle calculation   
     3870    REAL(wp) ::  zgl         = 0.0_wp   !< calculated exposure by direct beam   
     3871    REAL(wp) ::  woz         = 0.0_wp   !< calculated exposure by diffuse radiation
     3872    REAL(wp) ::  wsp         = 0.0_wp   !< calculated exposure by direct beam   
     3873   
     3874
     3875    CALL calc_date_and_time
     3876    doy_r = real(day_of_year)   
     3877    dtor = pi / 180.0_wp
     3878    lat = latitude
     3879    lon = longitude
     3880!
     3881!-- calculation of js, necessary for calculation of equation of time (zgl) :
     3882    js=  72.0_wp * ( doy_r + ( time_utc / 86400.0_wp ) ) / 73.0_wp
     3883!
     3884!-- calculation of equation of time (zgl):
     3885    zgl = 0.0066_wp + 7.3525_wp * cos( ( js + 85.9_wp ) * dtor ) + 9.9359_wp *                                        &
     3886    cos( ( 2.0_wp * js + 108.9_wp ) * dtor ) + 0.3387_wp * cos( ( 3 * js + 105.2_wp ) * dtor )
     3887!
     3888!-- calculation of apparent solar time woz:
     3889    woz = ( ( time_utc / 3600.0_wp ) - ( 4.0_wp * ( 15.0_wp - lon ) ) / 60.0_wp ) + ( zgl / 60.0_wp )
     3890!
     3891!-- calculation of hour angle (wsp):
     3892    wsp = ( woz - 12.0_wp ) * 15.0_wp
     3893!
     3894!-- calculation of declination:
     3895    declination = 0.3948_wp - 23.2559_wp * cos( ( js + 9.1_wp ) * dtor ) -                                            &
     3896    0.3915_wp * cos( ( 2.0_wp * js + 5.4_wp ) * dtor ) - 0.1764_wp * cos( ( 3.0_wp * js + 26.0_wp ) * dtor )
     3897!
     3898!-- calculation of solar zenith angle
     3899    thetar  = acos( sin( lat * dtor) * sin( declination * dtor ) + cos( wsp * dtor ) *                                &
     3900    cos( lat * dtor ) * cos( declination * dtor ) )
     3901    thetasr = asin( sin( lat * dtor) * sin( declination * dtor ) + cos( wsp * dtor ) *                                &
     3902    cos( lat * dtor ) * cos( declination * dtor ) )
     3903    sza = thetar / dtor
     3904!
     3905!-- calculation of solar azimuth angle
     3906    IF (woz .LE. 12.0_wp) alpha = pi - acos( ( sin(thetasr) * sin( lat * dtor ) -                                     &
     3907    sin( declination * dtor ) ) / ( cos(thetasr) * cos( lat * dtor ) ) )   
     3908    IF (woz .GT. 12.0_wp) alpha = pi + acos( ( sin(thetasr) * sin( lat * dtor ) -                                     &
     3909    sin( declination * dtor ) ) / ( cos(thetasr) * cos( lat * dtor ) ) )   
     3910    saa = alpha / dtor
     3911
     3912 END SUBROUTINE uvem_solar_position
     3913
     3914
     3915!------------------------------------------------------------------------------!
     3916! Description:
     3917! ------------
     3918!> Module-specific routine for new module
     3919!---------------------------------------------------------------------------------------------------------------------!
     3920 SUBROUTINE uvem_calc_exposure
     3921
     3922    USE indices,                                                                                                      &
     3923        ONLY:  nxlg, nxrg, nyng, nysg, nys, nyn, nxl, nxr
     3924   
     3925   
     3926    IMPLICIT NONE   
     3927   
     3928    INTEGER(iwp) ::  i     !< loop index in x direction
     3929    INTEGER(iwp) ::  j     !< loop index in y direction
     3930    INTEGER(iwp) ::  szai  !< loop index for different sza values
     3931
     3932    CALL uvem_solar_position
     3933     
     3934    IF (sza  .GE.  90) THEN
     3935       vitd3_exposure(:,:) = 0.0_wp
     3936    ELSE
     3937       
     3938       DO  ai = 0, 35
     3939          DO  zi = 0, 9
     3940                projection_area_lookup_table(ai,zi) = uvem_projarea_f%var(clothing,zi,ai)
     3941          ENDDO
     3942       ENDDO
     3943       DO  ai = 0, 35
     3944          DO  zi = 0, 9
     3945                integration_array(ai,zi) = uvem_integration_f%var(zi,ai)
     3946          ENDDO
     3947       ENDDO
     3948       DO  ai = 0, 2
     3949          DO  zi = 0, 90
     3950                irradiance_lookup_table(ai,zi) = uvem_irradiance_f%var(zi,ai)
     3951          ENDDO
     3952       ENDDO
     3953       DO  ai = 0, 35
     3954          DO  zi = 0, 9
     3955             DO  szai = 0, 90
     3956                radiance_lookup_table(ai,zi,szai) = uvem_radiance_f%var(szai,zi,ai)
     3957             ENDDO
     3958          ENDDO
     3959       ENDDO
     3960       
     3961       
     3962       
     3963!--    rotate 3D-Model human to desired direction  -----------------------------
     3964       projection_area_temp( 0:35,:) = projection_area_lookup_table
     3965       projection_area_temp(36:71,:) = projection_area_lookup_table               
     3966       IF (  .NOT.  turn_to_sun ) startpos_human = orientation_angle / 10.0_wp
     3967       IF (       turn_to_sun ) startpos_human = saa / 10.0_wp       
     3968       DO  ai = 0, 35
     3969          xfactor = ( startpos_human ) - INT( startpos_human )
     3970          DO  zi = 0, 9
     3971             projection_area(ai,zi) = ( projection_area_temp( 36 - INT( startpos_human ) - 1 + ai , zi) *             &
     3972                                      ( xfactor ) )                                                                   &
     3973                                      +( projection_area_temp( 36 - INT( startpos_human ) + ai , zi) *                &
     3974                                      ( 1.0_wp - xfactor ) )
     3975          ENDDO
     3976       ENDDO           
     3977!             
     3978!           
     3979!--    interpolate to accurate Solar Zenith Angle  ------------------         
     3980       DO  ai = 0, 35
     3981          xfactor = (sza)-INT(sza)
     3982          DO  zi = 0, 9
     3983             radiance_array(ai,zi) = ( radiance_lookup_table(ai, zi, INT(sza) ) * ( 1.0_wp - xfactor) ) +             &
     3984             ( radiance_lookup_table(ai,zi,INT(sza) + 1) * xfactor )
     3985          ENDDO
     3986       ENDDO
     3987       DO  iq = 0, 2
     3988          irradiance(iq) = ( irradiance_lookup_table(iq, INT(sza) ) * ( 1.0_wp - xfactor)) +                          &
     3989          (irradiance_lookup_table(iq, INT(sza) + 1) * xfactor )
     3990       ENDDO   
     3991!         
     3992!--    interpolate to accurate Solar Azimuth Angle ------------------
     3993       IF ( sun_in_south )  THEN
     3994          startpos_saa_float = 180.0_wp / 10.0_wp
     3995       ELSE
     3996          startpos_saa_float = saa / 10.0_wp
     3997       ENDIF
     3998       radiance_array_temp( 0:35,:) = radiance_array
     3999       radiance_array_temp(36:71,:) = radiance_array
     4000       xfactor = (startpos_saa_float) - INT(startpos_saa_float)
     4001       DO  ai = 0, 35
     4002          DO  zi = 0, 9
     4003             radiance_array(ai,zi) = ( radiance_array_temp( 36 - INT( startpos_saa_float ) - 1 + ai , zi ) *          &
     4004                                     ( xfactor ) )                                                                    &
     4005                                     + ( radiance_array_temp( 36 - INT( startpos_saa_float ) + ai , zi )              &
     4006                                     * ( 1.0_wp - xfactor ) )
     4007          ENDDO
     4008       ENDDO
     4009!       
     4010!     
     4011!--    calculate Projectionarea for direct beam -----------------------------'
     4012       projection_area_direct_temp( 0:35,:) = projection_area
     4013       projection_area_direct_temp(36:71,:) = projection_area
     4014       yfactor = ( sza / 10.0_wp ) - INT( sza / 10.0_wp )
     4015       xfactor = ( startpos_saa_float ) - INT( startpos_saa_float )
     4016       projection_area_direct_beam = ( projection_area_direct_temp( INT(startpos_saa_float)    ,INT(sza/10.0_wp)  ) * &
     4017                                     ( 1.0_wp - xfactor ) * ( 1.0_wp - yfactor ) ) +                                  &
     4018                                     ( projection_area_direct_temp( INT(startpos_saa_float) + 1,INT(sza/10.0_wp)  ) * &
     4019                                     (          xfactor ) * ( 1.0_wp - yfactor ) ) +                                  &
     4020                                     ( projection_area_direct_temp( INT(startpos_saa_float)    ,INT(sza/10.0_wp)+1) * &
     4021                                     ( 1.0_wp - xfactor ) * (          yfactor ) ) +                                  &
     4022                                     ( projection_area_direct_temp( INT(startpos_saa_float) + 1,INT(sza/10.0_wp)+1) * &
     4023                                     (          xfactor ) * (          yfactor ) )
     4024!                                               
     4025!                                               
     4026!                                               
     4027       DO  i = nxl, nxr    !nxlg, nxrg
     4028          DO  j = nys, nyn    !nysg, nyng
     4029!                   
     4030! !--        extract obstruction from IBSET-Integer_Array ------------------'
     4031             IF (consider_obstructions ) THEN
     4032                obstruction_temp1 = building_obstruction_f%var_3d(:,j,i)
     4033                IF (obstruction_temp1(0)  .NE.  9) THEN
     4034                   DO  pobi = 0, 44
     4035                      DO  bi = 0, 7
     4036                         IF ( btest( obstruction_temp1(pobi), bi )  .EQV.  .TRUE.) THEN
     4037                            obstruction_temp2( ( pobi * 8 ) + bi ) = 1
     4038                         ELSE
     4039                            obstruction_temp2( ( pobi * 8 ) + bi ) = 0
     4040                         ENDIF
     4041                      ENDDO
     4042                   ENDDO       
     4043                   DO  zi = 0, 9                                         
     4044                      obstruction(:,zi) = obstruction_temp2( zi * 36 :( zi * 36) + 35 )
     4045                   ENDDO
     4046                ELSE
     4047                   obstruction(:,:) = 0
     4048                ENDIF
     4049             ENDIF
     4050!             
     4051! !--        calculated human exposure ------------------' 
     4052             diffuse_exposure = SUM( radiance_array * projection_area * integration_array * obstruction )     
     4053         
     4054             obstruction_direct_beam = obstruction( nint(startpos_saa_float), nint( sza / 10.0_wp ) )
     4055             IF (sza  .GE.  89.99_wp) THEN
     4056                sza = 89.99999_wp
     4057             ENDIF
     4058!             
     4059!--          calculate direct normal irradiance (direct beam) ------------------'
     4060             direct_exposure = ( irradiance(1) / cos( pi * sza / 180.0_wp ) ) * &
     4061             projection_area_direct_beam * obstruction_direct_beam 
     4062               
     4063             vitd3_exposure(j,i) = ( diffuse_exposure + direct_exposure ) / 1000.0_wp * 70.97_wp
     4064             ! unit = international units vitamin D per second             
     4065          ENDDO
     4066       ENDDO
     4067    ENDIF
     4068
     4069 END SUBROUTINE uvem_calc_exposure
    35514070
    35524071 END MODULE biometeorology_mod
Note: See TracChangeset for help on using the changeset viewer.