source: palm/trunk/SOURCE/biometeorology_mod.f90 @ 4182

Last change on this file since 4182 was 4182, checked in by scharf, 2 years ago
  • corrected "Former revisions" section
  • minor formatting in "Former revisions" section
  • added "Author" section
  • Property svn:keywords set to Id
  • Property svn:mergeinfo set to (toggle deleted branches)
    /palm/branches/chemistry/SOURCE/biometeorology_mod.f902047-3190,​3218-3297
    /palm/branches/resler/SOURCE/biometeorology_mod.f902023-3320,​3337-3474
    /palm/branches/salsa/SOURCE/biometeorology_mod.f902503-3581
    /palm/trunk/SOURCE/biometeorology_mod.f90mergedeligible
    /palm/branches/forwind/SOURCE/biometeorology_mod.f901564-1913
    /palm/branches/fricke/SOURCE/biometeorology_mod.f90942-977
    /palm/branches/hoffmann/SOURCE/biometeorology_mod.f90989-1052
    /palm/branches/letzel/masked_output/SOURCE/biometeorology_mod.f90296-409
    /palm/branches/palm4u/SOURCE/biometeorology_mod.f902540-2692
    /palm/branches/rans/SOURCE/biometeorology_mod.f902078-3128
    /palm/branches/suehring/biometeorology_mod.f90423-666
File size: 181.3 KB
RevLine 
[3448]1!> @file biometeorology_mod.f90
2!--------------------------------------------------------------------------------!
[3321]3! This file is part of PALM-4U.
4!
5! PALM-4U is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM-4U is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
[3650]17! Copyright 2018-2019 Deutscher Wetterdienst (DWD)
18! Copyright 2018-2019 Institute of Computer Science, Academy of Sciences, Prague
19! Copyright 2018-2019 Leibniz Universitaet Hannover
[3448]20!--------------------------------------------------------------------------------!
[3321]21!
[3337]22! Current revisions:
[3569]23! ------------------
[3337]24!
25!
26! Former revisions:
[3321]27! -----------------
28! $Id: biometeorology_mod.f90 4182 2019-08-22 15:20:23Z scharf $
[4182]29! Corrected "Former revisions" section
30!
31! 4168 2019-08-16 13:50:17Z suehring
[4168]32! Replace function get_topography_top_index by topo_top_ind
33!
34! 4144 2019-08-06 09:11:47Z raasch
[4144]35! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
36!
37! 4127 2019-07-30 14:47:10Z suehring
[4127]38! Output for bio_mrt added (merge from branch resler)
39!
40! 4126 2019-07-30 11:09:11Z gronemeier
[4126]41! renamed vitd3_exposure_av into vitd3_dose,
42! renamed uvem_calc_exposure into bio_calculate_uv_exposure
43!
44! 3885 2019-04-11 11:29:34Z kanani
[3885]45! Changes related to global restructuring of location messages and introduction
46! of additional debug messages
47!
48! 3753 2019-02-19 14:48:54Z dom_dwd_user
[3753]49! - Added automatic setting of mrt_nlevels in case it was not part of
50! radiation_parameters namelist (or set to 0 accidentially).
51! - Minor speed improvoemnts in perceived temperature calculations.
52! - Perceived temperature regression arrays now declared as PARAMETERs.
53!
54! 3750 2019-02-19 07:29:39Z dom_dwd_user
[3749]55! - Added addittional safety meassures to bio_calculate_thermal_index_maps.
56! - Replaced several REAL (un-)equality comparisons.
57!
58! 3742 2019-02-14 11:25:22Z dom_dwd_user
[3742]59! - Allocation of the input _av grids was moved to the "sum" section of
60! bio_3d_data_averaging to make sure averaging is only done once!
61! - Moved call of bio_calculate_thermal_index_maps from biometeorology module to
62! time_integration to make sure averaged input is updated before calculating.
63!
64! 3740 2019-02-13 12:35:12Z dom_dwd_user
[3740]65! - Added safety-meassure to catch the case that 'bio_mrt_av' is stated after
66! 'bio_<index>' in the output section of the p3d file.
67!
68! 3739 2019-02-13 08:05:17Z dom_dwd_user
[3739]69! - Auto-adjusting thermal_comfort flag if not set by user, but thermal_indices
70! set as output quantities.
71! - Renamed flags "bio_<index>" to "do_calculate_<index>" for better readability
72! - Removed everything related to "time_bio_results" as this is never used.
73! - Moved humidity warning to check_data_output
74! - Fixed bug in mrt calculation introduced with my commit yesterday.
75!
76! 3735 2019-02-12 09:52:40Z dom_dwd_user
[3735]77! - Fixed auto-setting of thermal index calculation flags by output
78!  as originally proposed by resler.
79! - removed bio_pet and outher configuration variables.
80! - Updated namelist.
81!
82! 3711 2019-01-31 13:44:26Z knoop
[3711]83! Introduced interface routine bio_init_checks + small error message changes
84!
85! 3693 2019-01-23 15:20:53Z dom_dwd_user
[3693]86! Added usage of time_averaged mean radiant temperature, together with calculation,
87! grid and restart routines. General cleanup and commenting.
88!
89! 3685 2019-01-21 01:02:11Z knoop
[3685]90! Some interface calls moved to module_interface + cleanup
91!
92! 3650 2019-01-04 13:01:33Z kanani
[4182]93! Bugfixes and additions for enabling restarts with biometeorology
[3650]94!
[4182]95! 3448 2018-10-29 18:14:31Z kanani
96! Initial revision
[3569]97!
[4182]98!
99!
100! Authors:
101! --------
102! @author Dominik Froehlich <dominik.froehlich@dwd.de>, thermal indices
103! @author Jaroslav Resler <resler@cs.cas.cz>, mean radiant temperature
104! @author Michael Schrempf <schrempf@muk.uni-hannover.de>, uv exposure
105!
106!
[3448]107! Description:
108! ------------
[3569]109!> Biometeorology module consisting of two parts:
110!> 1.: Human thermal comfort module calculating thermal perception of a sample
[3448]111!> human being under the current meteorological conditions.
[3569]112!> 2.: Calculation of vitamin-D weighted UV exposure
[3448]113!>
114!> @todo Alphabetical sorting of "USE ..." lists, "ONLY" list, variable declarations
115!>       (per subroutine: first all CHARACTERs, then INTEGERs, LOGICALs, REALs, )
116!> @todo Comments start with capital letter --> "!-- Include..."
[3569]117!> @todo uv_vitd3dose-->new output type necessary (cumulative)
118!> @todo consider upwelling radiation in UV
[3448]119!>
120!> @note nothing now
121!>
122!> @bug  no known bugs by now
[3321]123!------------------------------------------------------------------------------!
[3448]124 MODULE biometeorology_mod
[3321]125
126    USE arrays_3d,                                                             &
[3742]127        ONLY:  pt, p, u, v, w, q
[3321]128
[3448]129    USE averaging,                                                             &
[3742]130        ONLY:  pt_av, q_av, u_av, v_av, w_av
[3448]131
[3361]132    USE basic_constants_and_equations_mod,                                     &
[3742]133        ONLY:  c_p, degc_to_k, l_v, magnus, sigma_sb, pi
[3361]134
[3448]135    USE control_parameters,                                                    &
[3885]136        ONLY:  average_count_3d, biometeorology,                               &
137               debug_output,                                                   &
138               dz, dz_stretch_factor,                                          &
[3742]139               dz_stretch_level, humidity, initializing_actions, nz_do3d,      &
140               surface_pressure
[3321]141
[3693]142    USE date_and_time_mod,                                                     &
[3569]143        ONLY:  calc_date_and_time, day_of_year, time_utc
144
[3448]145    USE grid_variables,                                                        &
[3742]146        ONLY:  ddx, dx, ddy, dy
[3321]147
[3448]148    USE indices,                                                               &
[3742]149        ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg,   &
[4168]150               nysg, nyng, topo_top_ind
[3321]151
[3448]152    USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM
[3321]153
[3693]154    USE netcdf_data_input_mod,                                                 &
155        ONLY:  netcdf_data_input_uvem, uvem_projarea_f, uvem_radiance_f,       &
[3569]156               uvem_irradiance_f, uvem_integration_f, building_obstruction_f
157!
[3448]158!-- Import radiation model to obtain input for mean radiant temperature
159    USE radiation_model_mod,                                                   &
[3742]160        ONLY:  ix, iy, iz, id, mrt_nlevels, mrt_include_sw,                    &
161               mrtinsw, mrtinlw, mrtbl, nmrtbl, radiation,                     &
162               radiation_interactions, rad_sw_in,                              &
163               rad_sw_out, rad_lw_in, rad_lw_out
[3321]164
[3448]165    IMPLICIT NONE
[3321]166
[3448]167    PRIVATE
[3321]168
[3569]169!
[3448]170!-- Declare all global variables within the module (alphabetical order)
[3593]171    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_grid  !< tmrt results (degree_C)
172    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct      !< PT results   (degree_C)
173    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci       !< UTCI results (degree_C)
174    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet        !< PET results  (degree_C)
[3569]175!
[3525]176!-- Grids for averaged thermal indices
177    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  mrt_av_grid   !< time average mean
[3693]178    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_av_grid  !< tmrt results (degree_C)
[3593]179    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct_av      !< PT results (aver. input)   (degree_C)
180    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_av       !< UTCI results (aver. input) (degree_C)
181    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_av        !< PET results (aver. input)  (degree_C)
[3321]182
[3525]183
184    INTEGER( iwp ) ::  bio_cell_level     !< cell level biom calculates for
185    REAL ( wp )    ::  bio_output_height  !< height output is calculated in m
[3448]186    REAL ( wp ), PARAMETER ::  human_absorb = 0.7_wp  !< SW absorbtivity of a human body (Fanger 1972)
187    REAL ( wp ), PARAMETER ::  human_emiss = 0.97_wp  !< LW emissivity of a human body after (Fanger 1972)
[3569]188    REAL ( wp ), PARAMETER ::  bio_fill_value = -9999._wp  !< set module fill value, replace by global fill value as soon as available
189!
[3448]190!--
[3742]191    LOGICAL ::  thermal_comfort  = .FALSE.  !< Enables or disables the entire thermal comfort part
[3693]192    LOGICAL ::  do_average_theta = .FALSE.  !< switch: do theta averaging in this module? (if .FALSE. this is done globally)
193    LOGICAL ::  do_average_q     = .FALSE.  !< switch: do e averaging in this module?
194    LOGICAL ::  do_average_u     = .FALSE.  !< switch: do u averaging in this module?
195    LOGICAL ::  do_average_v     = .FALSE.  !< switch: do v averaging in this module?
196    LOGICAL ::  do_average_w     = .FALSE.  !< switch: do w averaging in this module?
197    LOGICAL ::  do_average_mrt   = .FALSE.  !< switch: do mrt averaging in this module?
[3742]198    LOGICAL ::  average_trigger_perct  = .FALSE.  !< update averaged input on call to bio_perct?
199    LOGICAL ::  average_trigger_utci   = .FALSE.  !< update averaged input on call to bio_utci?
200    LOGICAL ::  average_trigger_pet    = .FALSE.  !< update averaged input on call to bio_pet?
[4127]201    LOGICAL ::  average_trigger_mrt    = .FALSE.  !< update averaged input on call to bio_pet?
[3742]202    LOGICAL ::  do_calculate_perct     = .FALSE.  !< Turn index PT (instant. input) on or off
203    LOGICAL ::  do_calculate_perct_av  = .FALSE.  !< Turn index PT (averaged input) on or off
204    LOGICAL ::  do_calculate_pet       = .FALSE.  !< Turn index PET (instant. input) on or off
205    LOGICAL ::  do_calculate_pet_av    = .FALSE.  !< Turn index PET (averaged input) on or off
206    LOGICAL ::  do_calculate_utci      = .FALSE.  !< Turn index UTCI (instant. input) on or off
207    LOGICAL ::  do_calculate_utci_av   = .FALSE.  !< Turn index UTCI (averaged input) on or off
[4127]208    LOGICAL ::  do_calculate_mrt2d     = .FALSE.  !< Turn index MRT 2D (averaged or inst) on or off
[3321]209
[3569]210!
211!-- UVEM parameters from here
212!
213!-- Declare all global variables within the module (alphabetical order)
[3650]214    INTEGER(iwp) ::  bio_nmrtbl
[3569]215    INTEGER(iwp) ::  ai                      = 0  !< loop index in azimuth direction
216    INTEGER(iwp) ::  bi                      = 0  !< loop index of bit location within an 8bit-integer (one Byte)
217    INTEGER(iwp) ::  clothing                = 1  !< clothing (0=unclothed, 1=Arms,Hands,Face free, 3=Hand,Face free)
218    INTEGER(iwp) ::  iq                      = 0  !< loop index of irradiance quantity
219    INTEGER(iwp) ::  pobi                    = 0  !< loop index of the position of corresponding byte within ibset byte vektor
220    INTEGER(iwp) ::  obstruction_direct_beam = 0  !< Obstruction information for direct beam   
221    INTEGER(iwp) ::  zi                      = 0  !< loop index in zenith direction
[3321]222
[3569]223    INTEGER(KIND=1), DIMENSION(0:44)  ::  obstruction_temp1 = 0  !< temporary obstruction information stored with ibset
224    INTEGER(iwp),    DIMENSION(0:359) ::  obstruction_temp2 = 0  !< restored temporary obstruction information from ibset file
225
226    INTEGER(iwp), DIMENSION(0:35,0:9) ::  obstruction       = 1  !< final 2D obstruction information array
227
228    LOGICAL ::  consider_obstructions = .TRUE.   !< namelist parameter (see documentation)
229    LOGICAL ::  sun_in_south          = .FALSE.  !< namelist parameter (see documentation)
230    LOGICAL ::  turn_to_sun           = .TRUE.   !< namelist parameter (see documentation)
231    LOGICAL ::  uv_exposure           = .FALSE.  !< namelist parameter (see documentation)
232
233    REAL(wp) ::  diffuse_exposure            =   0.0_wp  !< calculated exposure by diffuse radiation
234    REAL(wp) ::  direct_exposure             =   0.0_wp  !< calculated exposure by direct solar beam   
235    REAL(wp) ::  orientation_angle           =   0.0_wp  !< orientation of front/face of the human model   
236    REAL(wp) ::  projection_area_direct_beam =   0.0_wp  !< projection area for direct solar beam
237    REAL(wp) ::  saa                         = 180.0_wp  !< solar azimuth angle
238    REAL(wp) ::  startpos_human              =   0.0_wp  !< start value for azimuth interpolation of human geometry array
239    REAL(wp) ::  startpos_saa_float          =   0.0_wp  !< start value for azimuth interpolation of radiance array
240    REAL(wp) ::  sza                         =  20.0_wp  !< solar zenith angle
241    REAL(wp) ::  xfactor                     =   0.0_wp  !< relative x-position used for interpolation
242    REAL(wp) ::  yfactor                     =   0.0_wp  !< relative y-position used for interpolation
243
244    REAL(wp), DIMENSION(0:2)  ::  irradiance =   0.0_wp  !< iradiance values extracted from irradiance lookup table 
245
246    REAL(wp), DIMENSION(0:2,0:90) ::  irradiance_lookup_table      = 0.0_wp  !< irradiance lookup table
247    REAL(wp), DIMENSION(0:35,0:9) ::  integration_array            = 0.0_wp  !< solid angle factors for hemispherical integration
248    REAL(wp), DIMENSION(0:35,0:9) ::  projection_area              = 0.0_wp  !< projection areas of a human (all directions)
249    REAL(wp), DIMENSION(0:35,0:9) ::  projection_area_lookup_table = 0.0_wp  !< human geometry lookup table (projection areas)
250    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_direct_temp  = 0.0_wp  !< temporary projection area for direct solar beam
251    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_temp         = 0.0_wp  !< temporary projection area for all directions
252    REAL(wp), DIMENSION(0:35,0:9) ::  radiance_array               = 0.0_wp  !< radiance extracted from radiance_lookup_table 
253    REAL(wp), DIMENSION(0:71,0:9) ::  radiance_array_temp          = 0.0_wp  !< temporary radiance data
254
[4126]255    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vitd3_exposure  !< result variable for instantaneous vitamin-D weighted exposures
256    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vitd3_dose      !< result variable for summation of vitamin-D weighted exposures
[3569]257
258    REAL(wp), DIMENSION(0:35,0:9,0:90) ::  radiance_lookup_table   = 0.0_wp  !< radiance lookup table
259
[3525]260!
261!-- INTERFACES that must be available to other modules (alphabetical order)
[3321]262
[3525]263    PUBLIC bio_3d_data_averaging, bio_check_data_output,                       &
264    bio_calculate_mrt_grid, bio_calculate_thermal_index_maps, bio_calc_ipt,    &
265    bio_check_parameters, bio_data_output_3d, bio_data_output_2d,              &
266    bio_define_netcdf_grid, bio_get_thermal_index_input_ij, bio_header,        &
[3739]267    bio_init, bio_init_checks, bio_parin, thermal_comfort,                     &
[3650]268    bio_nmrtbl, bio_wrd_local, bio_rrd_local, bio_wrd_global, bio_rrd_global
[3569]269!
270!-- UVEM PUBLIC variables and methods
[4126]271    PUBLIC bio_calculate_uv_exposure, uv_exposure
[3525]272
[3448]273!
274!-- PALM interfaces:
275!
276!-- 3D averaging for HTCM _INPUT_ variables
[3525]277    INTERFACE bio_3d_data_averaging
278       MODULE PROCEDURE bio_3d_data_averaging
279    END INTERFACE bio_3d_data_averaging
[3569]280!
[3525]281!-- Calculate mtr from rtm fluxes and assign into 2D grid
282    INTERFACE bio_calculate_mrt_grid
283       MODULE PROCEDURE bio_calculate_mrt_grid
284    END INTERFACE bio_calculate_mrt_grid
[3569]285!
[3448]286!-- Calculate static thermal indices PT, UTCI and/or PET
[3525]287    INTERFACE bio_calculate_thermal_index_maps
288       MODULE PROCEDURE bio_calculate_thermal_index_maps
289    END INTERFACE bio_calculate_thermal_index_maps
[3569]290!
[3448]291!-- Calculate the dynamic index iPT (to be caled by the agent model)
[3525]292    INTERFACE bio_calc_ipt
293       MODULE PROCEDURE bio_calc_ipt
294    END INTERFACE bio_calc_ipt
[3569]295!
[3448]296!-- Data output checks for 2D/3D data to be done in check_parameters
[3569]297    INTERFACE bio_check_data_output
298       MODULE PROCEDURE bio_check_data_output
299    END INTERFACE bio_check_data_output
300!
[3448]301!-- Input parameter checks to be done in check_parameters
[3525]302    INTERFACE bio_check_parameters
303       MODULE PROCEDURE bio_check_parameters
304    END INTERFACE bio_check_parameters
[3569]305!
[3448]306!-- Data output of 2D quantities
[3525]307    INTERFACE bio_data_output_2d
308       MODULE PROCEDURE bio_data_output_2d
309    END INTERFACE bio_data_output_2d
[3569]310!
[3448]311!-- no 3D data, thus, no averaging of 3D data, removed
[3525]312    INTERFACE bio_data_output_3d
313       MODULE PROCEDURE bio_data_output_3d
314    END INTERFACE bio_data_output_3d
[3569]315!
[3448]316!-- Definition of data output quantities
[3525]317    INTERFACE bio_define_netcdf_grid
318       MODULE PROCEDURE bio_define_netcdf_grid
319    END INTERFACE bio_define_netcdf_grid
[3569]320!
[3448]321!-- Obtains all relevant input values to estimate local thermal comfort/stress
[3525]322    INTERFACE bio_get_thermal_index_input_ij
323       MODULE PROCEDURE bio_get_thermal_index_input_ij
324    END INTERFACE bio_get_thermal_index_input_ij
[3569]325!
[3448]326!-- Output of information to the header file
[3525]327    INTERFACE bio_header
328       MODULE PROCEDURE bio_header
329    END INTERFACE bio_header
[3569]330!
[3448]331!-- Initialization actions
[3525]332    INTERFACE bio_init
333       MODULE PROCEDURE bio_init
334    END INTERFACE bio_init
[3569]335!
[3711]336!-- Initialization checks
337    INTERFACE bio_init_checks
338       MODULE PROCEDURE bio_init_checks
339    END INTERFACE bio_init_checks
340!
[3448]341!-- Reading of NAMELIST parameters
[3525]342    INTERFACE bio_parin
343       MODULE PROCEDURE bio_parin
344    END INTERFACE bio_parin
[3569]345!
[3650]346!-- Read global restart parameters
347    INTERFACE bio_rrd_global
348       MODULE PROCEDURE bio_rrd_global
349    END INTERFACE bio_rrd_global
350!
351!-- Read local restart parameters
352    INTERFACE bio_rrd_local
353       MODULE PROCEDURE bio_rrd_local
354    END INTERFACE bio_rrd_local
355!
356!-- Write global restart parameters
357    INTERFACE bio_wrd_global
358       MODULE PROCEDURE bio_wrd_global
359    END INTERFACE bio_wrd_global
360!
361!-- Write local restart parameters
362    INTERFACE bio_wrd_local
363       MODULE PROCEDURE bio_wrd_local
364    END INTERFACE bio_wrd_local
365!
[3569]366!-- Calculate UV exposure grid
[4126]367    INTERFACE bio_calculate_uv_exposure
368       MODULE PROCEDURE bio_calculate_uv_exposure
369    END INTERFACE bio_calculate_uv_exposure
[3525]370
[3448]371 CONTAINS
[3525]372
373
[3321]374!------------------------------------------------------------------------------!
375! Description:
376! ------------
[3448]377!> Sum up and time-average biom input quantities as well as allocate
378!> the array necessary for storing the average.
[3569]379!> There is a considerable difference to the 3d_data_averaging subroutines
380!> used by other modules:
381!> For the thermal indices, the module needs to average the input conditions
382!> not the result!
[3321]383!------------------------------------------------------------------------------!
[3525]384 SUBROUTINE bio_3d_data_averaging( mode, variable )
[3321]385
386    IMPLICIT NONE
387
[3525]388    CHARACTER (LEN=*) ::  mode     !< averaging mode: allocate, sum, or average
389    CHARACTER (LEN=*) ::  variable !< The variable in question
[3321]390
[3525]391    INTEGER(iwp) ::  i        !< Running index, x-dir
392    INTEGER(iwp) ::  j        !< Running index, y-dir
393    INTEGER(iwp) ::  k        !< Running index, z-dir
[3321]394
395
396    IF ( mode == 'allocate' )  THEN
397
398       SELECT CASE ( TRIM( variable ) )
[3448]399
[3525]400          CASE ( 'bio_mrt' )
[3742]401
402                IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
[3525]403                   ALLOCATE( mrt_av_grid(nmrtbl) )
[3321]404                ENDIF
[3525]405                mrt_av_grid = 0.0_wp
[3742]406                do_average_mrt = .FALSE.  !< overwrite if that was enabled somehow
[3321]407
[3740]408
[4127]409          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*', 'bio_mrt*' )
[3448]410
[3569]411!
412!--          Averaging, as well as the allocation of the required grids must be
[3693]413!--          done only once, independent from for how many thermal indices
414!--          averaged output is desired.
415!--          Therefore wee need to memorize which index is the one that controls
416!--          the averaging (what must be the first thermal index called).
417!--          Indices are in unknown order as depending on the input file,
418!--          determine first index to average und update only once
[3740]419!
[3569]420!--          Only proceed here if this was not done for any index before. This
[3693]421!--          is done only once during the whole model run.
[3742]422             IF ( .NOT. average_trigger_perct  .AND.                           &
423                  .NOT. average_trigger_utci   .AND.                           &
[4127]424                  .NOT. average_trigger_pet    .AND.                           &
425                  .NOT. average_trigger_mrt )  THEN
[3569]426!
427!--             Memorize the first index called to control averaging
[3742]428                IF ( TRIM( variable ) == 'bio_perct*' )  THEN
[3525]429                    average_trigger_perct = .TRUE.
[3321]430                ENDIF
[3742]431                IF ( TRIM( variable ) == 'bio_utci*' )  THEN
[3448]432                    average_trigger_utci = .TRUE.
433                ENDIF
[3742]434                IF ( TRIM( variable ) == 'bio_pet*' )  THEN
[3448]435                    average_trigger_pet = .TRUE.
436                ENDIF
[4127]437                IF ( TRIM( variable ) == 'bio_mrt*' )  THEN
438                    average_trigger_mrt = .TRUE.
439                ENDIF
[3448]440             ENDIF
[3569]441!
[3742]442!--          Allocation of the input _av grids was moved to the "sum" section to
443!--          make sure averaging is only done once!
[3448]444
445
[3569]446          CASE ( 'uvem_vitd3dose*' )
[4126]447             IF ( .NOT. ALLOCATED( vitd3_dose ) )  THEN
448                ALLOCATE( vitd3_dose(nysg:nyng,nxlg:nxrg) )
[3569]449             ENDIF
[4126]450             vitd3_dose = 0.0_wp
[3569]451
[3321]452          CASE DEFAULT
453             CONTINUE
454
455       END SELECT
456
457    ELSEIF ( mode == 'sum' )  THEN
458
459       SELECT CASE ( TRIM( variable ) )
460
[3525]461          CASE ( 'bio_mrt' )
[3740]462!
463!--          Consider the case 'bio_mrt' is called after some thermal index. In
464!--          that case do_average_mrt will be .TRUE. leading to a double-
465!--          averaging.
[3742]466             IF ( .NOT. do_average_mrt  .AND.  ALLOCATED( mrt_av_grid ) )  THEN
[3321]467
[3525]468                IF ( mrt_include_sw )  THEN
469                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
[3739]470                      ( ( human_absorb * mrtinsw(:) +                          &
[3740]471                      mrtinlw(:) ) /                                           &
[3742]472                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
[3525]473                ELSE
474                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
[3740]475                      ( mrtinlw(:) /                                           &
[3742]476                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
[3525]477                ENDIF
[3321]478             ENDIF
479
[4127]480          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*', 'bio_mrt*' )
[3569]481!
[3740]482!--          Only continue if the current index is the one to trigger the input
483!--          averaging, see above
[3742]484             IF ( average_trigger_perct  .AND.  TRIM( variable ) /=            &
485                'bio_perct*')  RETURN
486             IF ( average_trigger_utci   .AND.  TRIM( variable ) /=            &
487                'bio_utci*')   RETURN
488             IF ( average_trigger_pet    .AND.  TRIM( variable ) /=            &
489                'bio_pet*')    RETURN
[4127]490             IF ( average_trigger_mrt    .AND.  TRIM( variable ) /=            &
491                'bio_mrt*')    RETURN
[3742]492!
493!--          Now memorize which of the input grids are not averaged by other
494!--          modules. Set averaging switch to .TRUE. and allocate the respective
495!--          grid in that case.
496             IF ( .NOT. ALLOCATED( pt_av ) )  THEN  !< if not averaged by other module
497                ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
498                do_average_theta = .TRUE.  !< memorize, that bio is responsible
499                pt_av = 0.0_wp
500             ENDIF
501             IF ( ALLOCATED( pt_av )  .AND.  do_average_theta )  THEN
[3448]502                DO  i = nxl, nxr
503                   DO  j = nys, nyn
504                      DO  k = nzb, nzt+1
505                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
506                      ENDDO
507                   ENDDO
508                ENDDO
509             ENDIF
[3321]510
[3742]511             IF ( .NOT. ALLOCATED( q_av ) )  THEN
512                ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
513                do_average_q = .TRUE.
514                q_av = 0.0_wp
515             ENDIF
516             IF ( ALLOCATED( q_av )  .AND.  do_average_q )  THEN
[3448]517                DO  i = nxl, nxr
518                   DO  j = nys, nyn
519                      DO  k = nzb, nzt+1
520                         q_av(k,j,i) = q_av(k,j,i) + q(k,j,i)
521                      ENDDO
522                   ENDDO
[3321]523                ENDDO
524             ENDIF
525
[3749]526!
527!-- u_av, v_av and w_av are always allocated
[3742]528             IF ( .NOT. ALLOCATED( u_av ) )  THEN
529                ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
530                do_average_u = .TRUE.
531                u_av = 0.0_wp
532             ENDIF
533             IF ( ALLOCATED( u_av )  .AND.  do_average_u )  THEN
[3525]534                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
535                   DO  j = nysg, nyng
[3448]536                      DO  k = nzb, nzt+1
537                         u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
538                      ENDDO
539                   ENDDO
540                ENDDO
541             ENDIF
542
[3742]543             IF ( .NOT. ALLOCATED( v_av ) )  THEN
544                ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
545                do_average_v = .TRUE.
546                v_av = 0.0_wp
547             ENDIF
548             IF ( ALLOCATED( v_av )  .AND.  do_average_v )  THEN
[3525]549                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
550                   DO  j = nysg, nyng
[3448]551                      DO  k = nzb, nzt+1
552                         v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
553                      ENDDO
554                   ENDDO
555                ENDDO
556             ENDIF
557
[3742]558             IF ( .NOT. ALLOCATED( w_av ) )  THEN
559                ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
560                do_average_w = .TRUE.
561                w_av = 0.0_wp
562             ENDIF
563             IF ( ALLOCATED( w_av )  .AND.  do_average_w )  THEN
[3525]564                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
565                   DO  j = nysg, nyng
[3448]566                      DO  k = nzb, nzt+1
567                         w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
568                      ENDDO
569                   ENDDO
570                ENDDO
571             ENDIF
[3693]572
[3742]573             IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
574                ALLOCATE( mrt_av_grid(nmrtbl) )
575                do_average_mrt = .TRUE.
576                mrt_av_grid = 0.0_wp
577             ENDIF
578             IF ( ALLOCATED( mrt_av_grid )  .AND.  do_average_mrt )  THEN
[3693]579
580                IF ( mrt_include_sw )  THEN
581                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
[3739]582                      ( ( human_absorb * mrtinsw(:) +                          &
[3742]583                      mrtinlw(:) ) /                                           &
584                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
[3693]585                ELSE
586                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
[3742]587                      ( mrtinlw(:) /                                           &
588                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
[3693]589                ENDIF
590             ENDIF
[3569]591!
592!--       This is a cumulated dose. No mode == 'average' for this quantity.
593          CASE ( 'uvem_vitd3dose*' )
[4126]594             IF ( ALLOCATED( vitd3_dose ) )  THEN
[3569]595                DO  i = nxlg, nxrg
596                   DO  j = nysg, nyng
[4126]597                      vitd3_dose(j,i) = vitd3_dose(j,i) + vitd3_exposure(j,i)
[3569]598                   ENDDO
599                ENDDO
600             ENDIF
[3448]601
[3569]602          CASE DEFAULT
[3321]603             CONTINUE
604
605       END SELECT
606
607    ELSEIF ( mode == 'average' )  THEN
608
609       SELECT CASE ( TRIM( variable ) )
610
[3525]611          CASE ( 'bio_mrt' )
[3740]612!
613!--          Consider the case 'bio_mrt' is called after some thermal index. In
614!--          that case do_average_mrt will be .TRUE. leading to a double-
615!--          averaging.
[3742]616             IF ( .NOT. do_average_mrt  .AND.  ALLOCATED( mrt_av_grid ) )  THEN
[3525]617                mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp )
[3321]618             ENDIF
619
[4127]620          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*', 'bio_mrt*' )
[3569]621!
622!--          Only continue if update index, see above
[3742]623             IF ( average_trigger_perct  .AND.                                 &
624                TRIM( variable ) /= 'bio_perct*' )  RETURN
625             IF ( average_trigger_utci  .AND.                                  &
626                TRIM( variable ) /= 'bio_utci*' )  RETURN
627             IF ( average_trigger_pet   .AND.                                  &
628                TRIM( variable ) /= 'bio_pet*' )  RETURN
[4127]629             IF ( average_trigger_mrt   .AND.                                  &
630                TRIM( variable ) /= 'bio_mrt*' )  RETURN
[3448]631
[3742]632             IF ( ALLOCATED( pt_av )  .AND.  do_average_theta )  THEN
[3448]633                DO  i = nxl, nxr
634                   DO  j = nys, nyn
635                      DO  k = nzb, nzt+1
[3569]636                         pt_av(k,j,i) = pt_av(k,j,i) /                         &
637                            REAL( average_count_3d, KIND=wp )
[3448]638                      ENDDO
639                   ENDDO
640                ENDDO
[3321]641             ENDIF
642
[3742]643             IF ( ALLOCATED( q_av )  .AND.  do_average_q )  THEN
[3448]644                DO  i = nxl, nxr
645                   DO  j = nys, nyn
646                      DO  k = nzb, nzt+1
[3569]647                         q_av(k,j,i) = q_av(k,j,i) /                           &
648                            REAL( average_count_3d, KIND=wp )
[3448]649                      ENDDO
650                   ENDDO
651                ENDDO
652             ENDIF
[3321]653
[3742]654             IF ( ALLOCATED( u_av )  .AND.  do_average_u )  THEN
[3525]655                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
656                   DO  j = nysg, nyng
[3448]657                      DO  k = nzb, nzt+1
[3569]658                         u_av(k,j,i) = u_av(k,j,i) /                           &
659                            REAL( average_count_3d, KIND=wp )
[3448]660                      ENDDO
661                   ENDDO
662                ENDDO
663             ENDIF
664
[3742]665             IF ( ALLOCATED( v_av )  .AND.  do_average_v )  THEN
[3525]666                DO  i = nxlg, nxrg
667                   DO  j = nysg, nyng
[3448]668                      DO  k = nzb, nzt+1
[3569]669                         v_av(k,j,i) = v_av(k,j,i) /                           &
670                            REAL( average_count_3d, KIND=wp )
[3448]671                      ENDDO
672                   ENDDO
673                ENDDO
674             ENDIF
675
[3742]676             IF ( ALLOCATED( w_av )  .AND.  do_average_w )  THEN
[3525]677                DO  i = nxlg, nxrg
678                   DO  j = nysg, nyng
[3448]679                      DO  k = nzb, nzt+1
[3569]680                         w_av(k,j,i) = w_av(k,j,i) /                           &
681                            REAL( average_count_3d, KIND=wp )
[3448]682                      ENDDO
683                   ENDDO
684                ENDDO
685             ENDIF
[3693]686
[3742]687             IF ( ALLOCATED( mrt_av_grid )  .AND.  do_average_mrt )  THEN
688                mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d,      &
689                KIND=wp )
[3693]690             ENDIF
691
[3569]692!
693!--     No averaging for UVEM since we are calculating a dose (only sum is
694!--     calculated and saved to av.nc file)
695
[3448]696        END SELECT
697
[3321]698    ENDIF
699
700
[3525]701 END SUBROUTINE bio_3d_data_averaging
[3321]702
[3448]703
704
[3321]705!------------------------------------------------------------------------------!
706! Description:
707! ------------
[3448]708!> Check data output for biometeorology model
[3321]709!------------------------------------------------------------------------------!
[3735]710 SUBROUTINE bio_check_data_output( var, unit, i, j, ilen, k )
[3321]711
[3479]712    USE control_parameters,                                                    &
713        ONLY: data_output, message_string
[3321]714
[3479]715    IMPLICIT NONE
[3321]716
[3479]717    CHARACTER (LEN=*) ::  unit     !< The unit for the variable var
718    CHARACTER (LEN=*) ::  var      !< The variable in question
[3321]719
[3739]720    INTEGER(iwp), INTENT(IN) ::  i     !< Current element of data_output
721    INTEGER(iwp), INTENT(IN) ::  j     !< Average quantity? 0 = no, 1 = yes
722    INTEGER(iwp), INTENT(IN) ::  ilen  !< Length of current entry in data_output
723    INTEGER(iwp), INTENT(IN) ::  k     !< Output is xy mode? 0 = no, 1 = yes
[3321]724
[3479]725    SELECT CASE ( TRIM( var ) )
[3569]726!
[3742]727!--    Allocate a temporary array with the desired output dimensions.
728!--    Arrays for time-averaged thermal indices are also allocated here because
729!--    they are not running through the standard averaging procedure in
730!--    bio_3d_data_averaging as the values of the averaged thermal indices are
731!--    derived in a single step based on priorly averaged arrays (see
732!--    bio_calculate_thermal_index_maps).
[4127]733       CASE ( 'bio_mrt', 'bio_mrt*' )
[3479]734          unit = 'degree_C'
[3739]735          thermal_comfort = .TRUE.  !< enable thermal_comfort if user forgot to do so
[3569]736          IF ( .NOT. ALLOCATED( tmrt_grid ) )  THEN
737             ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
[3650]738             tmrt_grid = REAL( bio_fill_value, KIND = wp )
[3569]739          ENDIF
[4127]740          IF ( TRIM( var ) == 'bio_mrt*' )  THEN
741             do_calculate_mrt2d = .TRUE.
742          END IF
[3321]743
[3569]744       CASE ( 'bio_perct*' )
745          unit = 'degree_C'
[3739]746          thermal_comfort = .TRUE.
[3742]747          IF ( j == 0 )  THEN                !< if instantaneous input
[3739]748             do_calculate_perct = .TRUE.
[3735]749             IF ( .NOT. ALLOCATED( perct ) )  THEN
750                ALLOCATE( perct (nys:nyn,nxl:nxr) )
751                perct = REAL( bio_fill_value, KIND = wp )
752             ENDIF
753          ELSE                              !< if averaged input
[3739]754             do_calculate_perct_av = .TRUE.
[3735]755             IF ( .NOT. ALLOCATED( perct_av ) )  THEN
756                ALLOCATE( perct_av (nys:nyn,nxl:nxr) )
757                perct_av = REAL( bio_fill_value, KIND = wp )
758             ENDIF
[3569]759          ENDIF
760
761       CASE ( 'bio_utci*' )
762          unit = 'degree_C'
[3739]763          thermal_comfort = .TRUE.
[3742]764          IF ( j == 0 )  THEN
[3739]765             do_calculate_utci = .TRUE.
[3735]766             IF ( .NOT. ALLOCATED( utci ) )  THEN
767                ALLOCATE( utci (nys:nyn,nxl:nxr) )
768                utci = REAL( bio_fill_value, KIND = wp )
769             ENDIF
770          ELSE
[3739]771             do_calculate_utci_av = .TRUE.
[3735]772             IF ( .NOT. ALLOCATED( utci_av ) )  THEN
773                ALLOCATE( utci_av (nys:nyn,nxl:nxr) )
774                utci_av = REAL( bio_fill_value, KIND = wp )
775             ENDIF
[3569]776          ENDIF
777
778       CASE ( 'bio_pet*' )
779          unit = 'degree_C'
[3739]780          thermal_comfort = .TRUE.
[3742]781          IF ( j == 0 )  THEN
[3739]782             do_calculate_pet = .TRUE.
[3735]783             IF ( .NOT. ALLOCATED( pet ) )  THEN
784                ALLOCATE( pet (nys:nyn,nxl:nxr) )
785                pet = REAL( bio_fill_value, KIND = wp )
786             ENDIF
787          ELSE
[3739]788             do_calculate_pet_av = .TRUE.
[3735]789             IF ( .NOT. ALLOCATED( pet_av ) )  THEN
790                ALLOCATE( pet_av (nys:nyn,nxl:nxr) )
791                pet_av = REAL( bio_fill_value, KIND = wp )
792             ENDIF
[3569]793          ENDIF
794
[3650]795
[3569]796       CASE ( 'uvem_vitd3*' )
[4126]797!           IF (  .NOT.  uv_exposure )  THEN
798!              message_string = 'output of "' // TRIM( var ) // '" requi' //     &
799!                       'res a namelist &uvexposure_par'
800!              CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
801!           ENDIF
[3569]802          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
803             message_string = 'illegal value for data_output: "' //            &
804                              TRIM( var ) // '" & only 2d-horizontal ' //      &
805                              'cross sections are allowed for this value'
806             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
807          ENDIF
808          unit = 'IU/s'
809          IF ( .NOT. ALLOCATED( vitd3_exposure ) )  THEN
810             ALLOCATE( vitd3_exposure(nysg:nyng,nxlg:nxrg) )
811          ENDIF
812          vitd3_exposure = 0.0_wp
813
814       CASE ( 'uvem_vitd3dose*' )
[4126]815!           IF (  .NOT.  uv_exposure )  THEN
816!              message_string = 'output of "' // TRIM( var ) // '" requi' //     &
817!                       'res  a namelist &uvexposure_par'
818!              CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
819!           ENDIF
[3569]820          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
821             message_string = 'illegal value for data_output: "' //            &
822                              TRIM( var ) // '" & only 2d-horizontal ' //      &
823                              'cross sections are allowed for this value'
824             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
825          ENDIF
826          unit = 'IU/av-h'
[4126]827          IF ( .NOT. ALLOCATED( vitd3_dose ) )  THEN
828             ALLOCATE( vitd3_dose(nysg:nyng,nxlg:nxrg) )
[3569]829          ENDIF
[4126]830          vitd3_dose = 0.0_wp
[3569]831
[3479]832       CASE DEFAULT
833          unit = 'illegal'
834
835    END SELECT
836
[3740]837!
838!--    Further checks if thermal comfort output is desired.
[3742]839    IF ( thermal_comfort  .AND.  unit == 'degree_C' )  THEN
[3479]840!
[3753]841!--    Break if required modules "radiation" is not avalable.
[3742]842       IF ( .NOT.  radiation )  THEN
[3525]843          message_string = 'output of "' // TRIM( var ) // '" require'         &
844                           // 's radiation = .TRUE.'
845          CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
846          unit = 'illegal'
[3479]847       ENDIF
[3753]848!
849!--    All "thermal_comfort" outputs except from 'bio_mrt' will also need
850!--    humidity input. Check also for that.
851       IF ( TRIM( var ) /= 'bio_mrt' )  THEN
852          IF ( .NOT.  humidity )  THEN
853             message_string = 'The estimation of thermal comfort '    //       &
854                              'requires air humidity information, but ' //     &
855                              'humidity module is disabled!'
856             CALL message( 'check_parameters', 'PA0561', 1, 2, 0, 6, 0 )
857             unit = 'illegal'
858          ENDIF
[3739]859       ENDIF
[3479]860
[3525]861
[3479]862    ENDIF
863
[3525]864 END SUBROUTINE bio_check_data_output
[3321]865
[3448]866!------------------------------------------------------------------------------!
867! Description:
868! ------------
869!> Check parameters routine for biom module
[3740]870!> Currently unused but might come in handy for future checks?
[3448]871!------------------------------------------------------------------------------!
[3525]872 SUBROUTINE bio_check_parameters
[3321]873
874
[3448]875    IMPLICIT NONE
[3321]876
[3448]877
878
[3525]879 END SUBROUTINE bio_check_parameters
[3448]880
881
[3321]882!------------------------------------------------------------------------------!
883! Description:
884! ------------
[3525]885!> Subroutine defining 2D output variables
886!> data_output_2d 1188ff
[3321]887!------------------------------------------------------------------------------!
[3525]888 SUBROUTINE bio_data_output_2d( av, variable, found, grid, local_pf,           &
[3569]889                                two_d, nzb_do, nzt_do )
[3321]890
891
892    USE kinds
893
894
895    IMPLICIT NONE
[3569]896!
[3448]897!-- Input variables
[3525]898    CHARACTER (LEN=*), INTENT(IN) ::  variable    !< Char identifier to select var for output
899    INTEGER(iwp), INTENT(IN)      ::  av          !< Use averaged data? 0 = no, 1 = yes?
900    INTEGER(iwp), INTENT(IN)      ::  nzb_do      !< Unused. 2D. nz bottom to nz top
901    INTEGER(iwp), INTENT(IN)      ::  nzt_do      !< Unused.
[3569]902!
[3448]903!-- Output variables
[3525]904    CHARACTER (LEN=*), INTENT(OUT) ::  grid   !< Grid type (always "zu1" for biom)
905    LOGICAL, INTENT(OUT)           ::  found  !< Output found?
[3693]906    LOGICAL, INTENT(OUT)           ::  two_d  !< Flag parameter that indicates 2D variables,
907                                              !< horizontal cross sections, must be .TRUE. for thermal indices and uv
[3525]908    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< Temp. result grid to return
[3569]909!
[3448]910!-- Internal variables
[3525]911    INTEGER(iwp) ::  i        !< Running index, x-dir
912    INTEGER(iwp) ::  j        !< Running index, y-dir
913    INTEGER(iwp) ::  k        !< Running index, z-dir
914    INTEGER(iwp) ::  l        !< Running index, radiation grid
[3321]915
916
[3525]917    found = .TRUE.
[3569]918    local_pf = bio_fill_value
[3525]919
[3569]920    SELECT CASE ( TRIM( variable ) )
[3321]921
922
[3525]923        CASE ( 'bio_mrt_xy' )
[3569]924           grid = 'zu1'
925           two_d = .FALSE.  !< can be calculated for several levels
926           local_pf = REAL( bio_fill_value, KIND = wp )
927           DO  l = 1, nmrtbl
928              i = mrtbl(ix,l)
929              j = mrtbl(iy,l)
930              k = mrtbl(iz,l)
[3742]931              IF ( k < nzb_do  .OR.  k > nzt_do  .OR.  j < nys  .OR.           &
932                 j > nyn  .OR.  i < nxl  .OR.  i > nxr )  CYCLE
[3569]933              IF ( av == 0 )  THEN
934                 IF ( mrt_include_sw )  THEN
[3739]935                    local_pf(i,j,k) = ( ( human_absorb * mrtinsw(l) +          &
[3740]936                                    mrtinlw(l) ) /                             &
[3742]937                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
[3739]938                                    degc_to_k
[3569]939                 ELSE
[3740]940                    local_pf(i,j,k) = ( mrtinlw(l)  /                          &
[3742]941                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
[3739]942                                    degc_to_k
[3569]943                 ENDIF
944              ELSE
945                 local_pf(i,j,k) = mrt_av_grid(l)
946              ENDIF
947           ENDDO
[3321]948
[4127]949        CASE ( 'bio_mrt*_xy' )        ! 2d-array
950           grid = 'zu1'
951           two_d = .TRUE.
952           IF ( av == 0 )  THEN
953              DO  i = nxl, nxr
954                 DO  j = nys, nyn
955                    local_pf(i,j,nzb+1) = tmrt_grid(j,i)
956                 ENDDO
957              ENDDO
958           ELSE
959              DO  i = nxl, nxr
960                 DO  j = nys, nyn
961                    local_pf(i,j,nzb+1) = tmrt_av_grid(j,i)
962                 ENDDO
963              ENDDO
964           ENDIF
[3448]965
[4127]966
[3525]967        CASE ( 'bio_perct*_xy' )        ! 2d-array
[3569]968           grid = 'zu1'
969           two_d = .TRUE.
970           IF ( av == 0 )  THEN
[3448]971              DO  i = nxl, nxr
972                 DO  j = nys, nyn
[3525]973                    local_pf(i,j,nzb+1) = perct(j,i)
[3448]974                 ENDDO
975              ENDDO
[3569]976           ELSE
[3448]977              DO  i = nxl, nxr
978                 DO  j = nys, nyn
[3525]979                    local_pf(i,j,nzb+1) = perct_av(j,i)
[3448]980                 ENDDO
981              ENDDO
[3742]982           ENDIF
[3321]983
[3525]984
985        CASE ( 'bio_utci*_xy' )        ! 2d-array
[3569]986           grid = 'zu1'
987           two_d = .TRUE.
988           IF ( av == 0 )  THEN
[3448]989              DO  i = nxl, nxr
990                 DO  j = nys, nyn
[3569]991                    local_pf(i,j,nzb+1) = utci(j,i)
[3448]992                 ENDDO
993              ENDDO
[3569]994           ELSE
[3448]995              DO  i = nxl, nxr
996                 DO  j = nys, nyn
[3569]997                    local_pf(i,j,nzb+1) = utci_av(j,i)
[3448]998                 ENDDO
999              ENDDO
[3742]1000           ENDIF
[3321]1001
[3525]1002
1003        CASE ( 'bio_pet*_xy' )        ! 2d-array
[3569]1004           grid = 'zu1'
1005           two_d = .TRUE.
1006           IF ( av == 0 )  THEN
[3448]1007              DO  i = nxl, nxr
1008                 DO  j = nys, nyn
[3569]1009                    local_pf(i,j,nzb+1) = pet(j,i)
[3448]1010                 ENDDO
1011              ENDDO
[3569]1012           ELSE
[3448]1013              DO  i = nxl, nxr
1014                 DO  j = nys, nyn
[3569]1015                    local_pf(i,j,nzb+1) = pet_av(j,i)
[3448]1016                 ENDDO
1017              ENDDO
[3742]1018           ENDIF
[3321]1019
[3693]1020!
[3569]1021!--    Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein.
1022!--    However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged.
1023       CASE ( 'uvem_vitd3*_xy' )        ! 2d-array
1024          IF ( av == 0 )  THEN
1025             DO  i = nxl, nxr
1026                DO  j = nys, nyn
1027                   local_pf(i,j,nzb+1) = vitd3_exposure(j,i)
1028                ENDDO
1029             ENDDO
1030          ENDIF
[3525]1031
[3569]1032          two_d = .TRUE.
1033          grid = 'zu1'
1034
1035       CASE ( 'uvem_vitd3dose*_xy' )        ! 2d-array
1036          IF ( av == 1 )  THEN
1037             DO  i = nxl, nxr
1038                DO  j = nys, nyn
[4126]1039                   local_pf(i,j,nzb+1) = vitd3_dose(j,i)
[3569]1040                ENDDO
1041             ENDDO
1042          ENDIF
1043
1044          two_d = .TRUE.
1045          grid = 'zu1'
1046
1047
[3448]1048       CASE DEFAULT
1049          found = .FALSE.
[3525]1050          grid  = 'none'
[3321]1051
[3448]1052    END SELECT
[3321]1053
[3448]1054
[3525]1055 END SUBROUTINE bio_data_output_2d
1056
1057
[3448]1058!------------------------------------------------------------------------------!
1059! Description:
1060! ------------
[3525]1061!> Subroutine defining 3D output variables (dummy, always 2d!)
1062!> data_output_3d 709ff
[3448]1063!------------------------------------------------------------------------------!
[3525]1064 SUBROUTINE bio_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
[3448]1065
[3525]1066    USE indices
[3448]1067
1068    USE kinds
1069
1070
1071    IMPLICIT NONE
[3569]1072!
[3448]1073!-- Input variables
[3525]1074    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< Char identifier to select var for output
1075    INTEGER(iwp), INTENT(IN) ::  av       !< Use averaged data? 0 = no, 1 = yes?
1076    INTEGER(iwp), INTENT(IN) ::  nzb_do   !< Unused. 2D. nz bottom to nz top
1077    INTEGER(iwp), INTENT(IN) ::  nzt_do   !< Unused.
[3569]1078!
[3448]1079!-- Output variables
[3525]1080    LOGICAL, INTENT(OUT) ::  found   !< Output found?
1081    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< Temp. result grid to return
[3569]1082!
[3448]1083!-- Internal variables
[3525]1084    INTEGER(iwp) ::  l    !< Running index, radiation grid
1085    INTEGER(iwp) ::  i    !< Running index, x-dir
1086    INTEGER(iwp) ::  j    !< Running index, y-dir
1087    INTEGER(iwp) ::  k    !< Running index, z-dir
1088
[3569]1089!     REAL(wp) ::  mrt  !< Buffer for mean radiant temperature
[3448]1090
1091    found = .TRUE.
[3525]1092
[3569]1093    SELECT CASE ( TRIM( variable ) )
[3448]1094
[3525]1095        CASE ( 'bio_mrt' )
[3569]1096            local_pf = REAL( bio_fill_value, KIND = sp )
[3525]1097            DO  l = 1, nmrtbl
1098               i = mrtbl(ix,l)
1099               j = mrtbl(iy,l)
1100               k = mrtbl(iz,l)
[3742]1101               IF ( k < nzb_do  .OR.  k > nzt_do  .OR.  j < nys  .OR.          &
1102                  j > nyn  .OR.  i < nxl  .OR.  i > nxr )  CYCLE
[3525]1103               IF ( av == 0 )  THEN
1104                  IF ( mrt_include_sw )  THEN
[3735]1105                     local_pf(i,j,k) = REAL( ( ( human_absorb * mrtinsw(l) +   &
[3740]1106                                    mrtinlw(l) ) /                             &
[3742]1107                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
[3735]1108                                    degc_to_k, KIND = sp )
[3525]1109                  ELSE
[3740]1110                     local_pf(i,j,k) = REAL( ( mrtinlw(l) /                    &
[3742]1111                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
[3735]1112                                    degc_to_k, KIND = sp )
[3525]1113                  ENDIF
1114               ELSE
[3735]1115                  local_pf(i,j,k) = REAL( mrt_av_grid(l), KIND = sp )
[3525]1116               ENDIF
1117            ENDDO
[3448]1118
[3321]1119       CASE DEFAULT
1120          found = .FALSE.
1121
1122    END SELECT
1123
[3525]1124 END SUBROUTINE bio_data_output_3d
[3321]1125
[3448]1126!------------------------------------------------------------------------------!
1127! Description:
1128! ------------
1129!> Subroutine defining appropriate grid for netcdf variables.
1130!> It is called out from subroutine netcdf_interface_mod.
1131!> netcdf_interface_mod 918ff
1132!------------------------------------------------------------------------------!
[3525]1133 SUBROUTINE bio_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
[3321]1134
[3448]1135    IMPLICIT NONE
[3569]1136!
[3448]1137!-- Input variables
1138    CHARACTER (LEN=*), INTENT(IN)  ::  var      !< Name of output variable
[3569]1139!
[3448]1140!-- Output variables
1141    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
1142    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !< y grid of output variable
1143    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !< z grid of output variable
1144
1145    LOGICAL, INTENT(OUT)           ::  found    !< Flag if output var is found
[3569]1146!
[3448]1147!-- Local variables
1148    LOGICAL      :: is2d  !< Var is 2d?
1149
1150    INTEGER(iwp) :: l     !< Length of the var array
1151
1152
1153    found  = .FALSE.
1154    grid_x = 'none'
1155    grid_y = 'none'
1156    grid_z = 'none'
1157
1158    l = MAX( 2, LEN_TRIM( var ) )
1159    is2d = ( var(l-1:l) == 'xy' )
1160
[3742]1161    IF ( var(1:4) == 'bio_' )  THEN
[3448]1162       found  = .TRUE.
1163       grid_x = 'x'
1164       grid_y = 'y'
1165       grid_z = 'zu'
[3742]1166       IF ( is2d  .AND.  var(1:7) /= 'bio_mrt' )  grid_z = 'zu1'
[3448]1167    ENDIF
1168
[3742]1169    IF ( is2d  .AND.  var(1:4) == 'uvem' )  THEN
[3569]1170       grid_x = 'x'
1171       grid_y = 'y'
1172       grid_z = 'zu1'
1173    ENDIF
1174
[3525]1175 END SUBROUTINE bio_define_netcdf_grid
[3448]1176
[3321]1177!------------------------------------------------------------------------------!
1178! Description:
1179! ------------
[3448]1180!> Header output for biom module
1181!> header 982
1182!------------------------------------------------------------------------------!
[3525]1183 SUBROUTINE bio_header( io )
[3448]1184
1185    IMPLICIT NONE
[3569]1186!
[3448]1187!-- Input variables
1188    INTEGER(iwp), INTENT(IN) ::  io           !< Unit of the output file
[3569]1189!
[3448]1190!-- Internal variables
1191    CHARACTER (LEN=86) ::  output_height_chr  !< String for output height
1192
[3525]1193    WRITE( output_height_chr, '(F8.1,7X)' )  bio_output_height
[3321]1194!
[3448]1195!-- Write biom header
1196    WRITE( io, 1 )
1197    WRITE( io, 2 )  TRIM( output_height_chr )
[3525]1198    WRITE( io, 3 )  TRIM( ACHAR( bio_cell_level ) )
[3448]1199
12001   FORMAT (//' Human thermal comfort module information:'/                    &
1201              ' ------------------------------'/)
12022   FORMAT ('    --> All indices calculated for a height of (m): ', A )
12033   FORMAT ('    --> This corresponds to cell level : ', A )
1204
[3525]1205 END SUBROUTINE bio_header
[3448]1206
1207
[3321]1208!------------------------------------------------------------------------------!
[3448]1209! Description:
1210! ------------
1211!> Initialization of the HTCM
1212!> init_3d_model 1987ff
1213!------------------------------------------------------------------------------!
[3525]1214 SUBROUTINE bio_init
[3321]1215
[3614]1216    USE netcdf_data_input_mod,                                                 &
1217        ONLY:  netcdf_data_input_uvem
[3569]1218
[3448]1219    IMPLICIT NONE
[3569]1220!
[3448]1221!-- Internal vriables
1222    REAL ( wp )  :: height  !< current height in meters
[3685]1223
[3885]1224    IF ( debug_output )  CALL debug_message( 'bio_init', 'start' )
[3569]1225!
[3448]1226!-- Determine cell level corresponding to 1.1 m above ground level
1227!   (gravimetric center of sample human)
[3321]1228
[3525]1229    bio_cell_level = 0_iwp
1230    bio_output_height = 0.5_wp * dz(1)
[3448]1231    height = 0.0_wp
[3321]1232
[3753]1233    bio_cell_level = INT( 1.099_wp / dz(1) )
[3525]1234    bio_output_height = bio_output_height + bio_cell_level * dz(1)
[3569]1235!
[3753]1236!-- Set radiation level if not done by user
1237    IF ( mrt_nlevels == 0 )  THEN
1238       mrt_nlevels = bio_cell_level + 1_iwp
1239    ENDIF
1240!
[3569]1241!-- Init UVEM and load lookup tables
[3650]1242    IF ( uv_exposure )  CALL netcdf_data_input_uvem
[3735]1243
[3885]1244    IF ( debug_output )  CALL debug_message( 'bio_init', 'end' )
[3569]1245
[3525]1246 END SUBROUTINE bio_init
[3321]1247
1248
[3448]1249!------------------------------------------------------------------------------!
1250! Description:
1251! ------------
[3711]1252!> Checks done after the Initialization
1253!------------------------------------------------------------------------------!
1254 SUBROUTINE bio_init_checks
1255
1256    USE control_parameters,                                                    &
1257        ONLY: message_string
1258
[4126]1259    IF ( (.NOT. radiation_interactions) .AND. ( thermal_comfort ) )  THEN
[3711]1260       message_string = 'The mrt calculation requires ' //                     &
1261                        'enabled radiation_interactions but it ' //            &
1262                        'is disabled!'
1263       CALL message( 'bio_init_checks', 'PAHU03', 1, 2, 0, 6, 0 )
1264    ENDIF
1265
1266
1267 END SUBROUTINE bio_init_checks
1268
1269
1270!------------------------------------------------------------------------------!
1271! Description:
1272! ------------
[3448]1273!> Parin for &biometeorology_parameters for reading biomet parameters
[3321]1274!------------------------------------------------------------------------------!
[3525]1275 SUBROUTINE bio_parin
[3321]1276
[3448]1277    IMPLICIT NONE
[3321]1278
[3448]1279!
1280!-- Internal variables
1281    CHARACTER (LEN=80) ::  line  !< Dummy string for current line in parameter file
[3321]1282
[3735]1283    NAMELIST /biometeorology_parameters/  thermal_comfort,                     &
1284
[3569]1285!
1286!-- UVEM namelist parameters
1287                                          clothing,                            &
1288                                          consider_obstructions,               &
1289                                          orientation_angle,                   &
1290                                          sun_in_south,                        &
1291                                          turn_to_sun,                         &
1292                                          uv_exposure
[3321]1293
1294
[3448]1295!-- Try to find biometeorology_parameters namelist
1296    REWIND ( 11 )
1297    line = ' '
1298    DO WHILE ( INDEX( line, '&biometeorology_parameters' ) == 0 )
1299       READ ( 11, '(A)', END = 20 )  line
1300    ENDDO
1301    BACKSPACE ( 11 )
[3321]1302
[3448]1303!
1304!-- Read biometeorology_parameters namelist
1305    READ ( 11, biometeorology_parameters, ERR = 10, END = 20 )
[3321]1306
[3448]1307!
1308!-- Set flag that indicates that the biomet_module is switched on
1309    biometeorology = .TRUE.
[3321]1310
[3448]1311    GOTO 20
[3321]1312
[3448]1313!
1314!-- In case of error
1315 10 BACKSPACE( 11 )
1316    READ( 11 , '(A)') line
1317    CALL parin_fail_message( 'biometeorology_parameters', line )
[3321]1318
[3448]1319!
1320!-- Complete
1321 20 CONTINUE
[3321]1322
1323
[3525]1324 END SUBROUTINE bio_parin
[3321]1325
1326!------------------------------------------------------------------------------!
1327! Description:
1328! ------------
[3693]1329!> Soubroutine reads global biometeorology configuration from restart file(s)
1330!------------------------------------------------------------------------------!
1331 SUBROUTINE bio_rrd_global( found )
1332
1333    USE control_parameters,                                                    &
1334        ONLY:  length, restart_string
1335
1336
1337    IMPLICIT NONE
1338
1339    LOGICAL, INTENT(OUT) ::  found      !< variable found? yes = .T., no = .F.
1340
1341    found = .TRUE.
1342
1343
1344    SELECT CASE ( restart_string(1:length) )
1345
1346!
1347!--    read control flags to determine if input grids need to be averaged
1348       CASE ( 'do_average_theta' )
1349          READ ( 13 )  do_average_theta
1350
1351       CASE ( 'do_average_q' )
1352          READ ( 13 )  do_average_q
1353
1354       CASE ( 'do_average_u' )
1355          READ ( 13 )  do_average_u
1356
1357       CASE ( 'do_average_v' )
1358          READ ( 13 )  do_average_v
1359
1360       CASE ( 'do_average_w' )
1361          READ ( 13 )  do_average_w
1362
1363       CASE ( 'do_average_mrt' )
1364          READ ( 13 )  do_average_mrt
1365
1366!
1367!--    read control flags to determine which thermal index needs to trigger averaging
1368       CASE ( 'average_trigger_perct' )
1369          READ ( 13 )  average_trigger_perct
1370
1371       CASE ( 'average_trigger_utci' )
1372          READ ( 13 )  average_trigger_utci
1373
1374       CASE ( 'average_trigger_pet' )
1375          READ ( 13 )  average_trigger_pet
1376
[4127]1377       CASE ( 'average_trigger_mrt' )
1378          READ ( 13 )  average_trigger_mrt
[3693]1379
[4127]1380
[3693]1381       CASE DEFAULT
1382
1383          found = .FALSE.
1384
1385    END SELECT
1386
1387
1388 END SUBROUTINE bio_rrd_global
1389
1390
1391!------------------------------------------------------------------------------!
1392! Description:
1393! ------------
1394!> Soubroutine reads local biometeorology configuration from restart file(s)
1395!------------------------------------------------------------------------------!
1396 SUBROUTINE bio_rrd_local( found )
1397
1398
1399    USE control_parameters,                                                    &
1400        ONLY:  length, restart_string
1401
1402
1403    IMPLICIT NONE
1404
1405
1406    LOGICAL, INTENT(OUT) ::  found      !< variable found? yes = .T., no = .F.
1407
1408    found = .TRUE.
1409
1410
1411    SELECT CASE ( restart_string(1:length) )
1412
1413       CASE ( 'nmrtbl' )
1414          READ ( 13 )  bio_nmrtbl
1415
1416       CASE ( 'mrt_av_grid' )
1417          IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
1418             ALLOCATE( mrt_av_grid(bio_nmrtbl) )
1419          ENDIF
1420          READ ( 13 )  mrt_av_grid
1421
1422
1423       CASE DEFAULT
1424
1425          found = .FALSE.
1426
1427    END SELECT
1428
1429
1430 END SUBROUTINE bio_rrd_local
1431
1432!------------------------------------------------------------------------------!
1433! Description:
1434! ------------
1435!> Write global restart data for the biometeorology module.
1436!------------------------------------------------------------------------------!
1437 SUBROUTINE bio_wrd_global
1438
1439    IMPLICIT NONE
1440
1441    CALL wrd_write_string( 'do_average_theta' )
1442    WRITE ( 14 )  do_average_theta
1443    CALL wrd_write_string( 'do_average_q' )
1444    WRITE ( 14 )  do_average_q
1445    CALL wrd_write_string( 'do_average_u' )
1446    WRITE ( 14 )  do_average_u
1447    CALL wrd_write_string( 'do_average_v' )
1448    WRITE ( 14 )  do_average_v
1449    CALL wrd_write_string( 'do_average_w' )
1450    WRITE ( 14 )  do_average_w
1451    CALL wrd_write_string( 'do_average_mrt' )
1452    WRITE ( 14 )  do_average_mrt
1453    CALL wrd_write_string( 'average_trigger_perct' )
1454    WRITE ( 14 )  average_trigger_perct
1455    CALL wrd_write_string( 'average_trigger_utci' )
1456    WRITE ( 14 )  average_trigger_utci
1457    CALL wrd_write_string( 'average_trigger_pet' )
1458    WRITE ( 14 )  average_trigger_pet
[4127]1459    CALL wrd_write_string( 'average_trigger_mrt' )
1460    WRITE ( 14 )  average_trigger_mrt
[3693]1461
1462 END SUBROUTINE bio_wrd_global
1463
1464
1465!------------------------------------------------------------------------------!
1466! Description:
1467! ------------
1468!> Write local restart data for the biometeorology module.
1469!------------------------------------------------------------------------------!
1470 SUBROUTINE bio_wrd_local
1471
1472    IMPLICIT NONE
1473
1474!
1475!-- First nmrtbl has to be written/read, because it is the dimension of mrt_av_grid
1476    CALL wrd_write_string( 'nmrtbl' )
1477    WRITE ( 14 )  nmrtbl
1478
1479    IF ( ALLOCATED( mrt_av_grid ) )  THEN
1480       CALL wrd_write_string( 'mrt_av_grid' )
1481       WRITE ( 14 )  mrt_av_grid
1482    ENDIF
1483
1484
1485 END SUBROUTINE bio_wrd_local
1486
1487!------------------------------------------------------------------------------!
1488! Description:
1489! ------------
[3525]1490!> Calculate biometeorology MRT for all 2D grid
[3321]1491!------------------------------------------------------------------------------!
[3525]1492 SUBROUTINE bio_calculate_mrt_grid ( av )
[3321]1493
[3448]1494    IMPLICIT NONE
[3321]1495
[3525]1496    LOGICAL, INTENT(IN)         ::  av    !< use averaged input?
[3569]1497!
[3525]1498!-- Internal variables
1499    INTEGER(iwp)                ::  i     !< Running index, x-dir, radiation coordinates
1500    INTEGER(iwp)                ::  j     !< Running index, y-dir, radiation coordinates
1501    INTEGER(iwp)                ::  k     !< Running index, y-dir, radiation coordinates
1502    INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
[3321]1503
[3693]1504
[3740]1505!
1506!-- We need to differentiate if averaged input is desired (av == .TRUE.) or not.
[3742]1507    IF ( av )  THEN
[3740]1508!
1509!-- Make sure tmrt_av_grid is present and initialize with the fill value
[3693]1510       IF ( .NOT. ALLOCATED( tmrt_av_grid ) )  THEN
1511          ALLOCATE( tmrt_av_grid (nys:nyn,nxl:nxr) )
1512       ENDIF
[3735]1513       tmrt_av_grid = REAL( bio_fill_value, KIND = wp )
[3693]1514
[3740]1515!
1516!-- mrt_av_grid should always be allcoated here, but better make sure ist actually is.
1517       IF ( ALLOCATED( mrt_av_grid ) )  THEN
1518!
1519!-- Iterate over the radiation grid (radiation coordinates) and fill the
1520!-- tmrt_av_grid (x, y coordinates) where appropriate:
1521!-- tmrt_av_grid is written for all i / j if level (k) matches output height.
1522          DO  l = 1, nmrtbl
1523             i = mrtbl(ix,l)
1524             j = mrtbl(iy,l)
1525             k = mrtbl(iz,l)
[4168]1526             IF ( k - topo_top_ind(j,i,0) == bio_cell_level + 1_iwp)  THEN
[3740]1527!
1528!-- Averaging was done before, so we can just copy the result here
1529                tmrt_av_grid(j,i) = mrt_av_grid(l)
[3693]1530
[3740]1531             ENDIF
1532          ENDDO
1533       ENDIF
[3693]1534
[3740]1535!
1536!-- In case instantaneous input is desired, mrt values will be re-calculated.
[3693]1537    ELSE
[3569]1538!
[3593]1539!-- Calculate biometeorology MRT from local radiation fluxes calculated by RTM and assign
1540!-- into 2D grid. Depending on selected output quantities, tmrt_grid might not have been
1541!-- allocated in bio_check_data_output yet.
[3693]1542       IF ( .NOT. ALLOCATED( tmrt_grid ) )  THEN
1543          ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
1544       ENDIF
1545       tmrt_grid = REAL( bio_fill_value, KIND = wp )
[3593]1546
[3693]1547       DO  l = 1, nmrtbl
1548          i = mrtbl(ix,l)
1549          j = mrtbl(iy,l)
1550          k = mrtbl(iz,l)
[4168]1551          IF ( k - topo_top_ind(j,i,0) == bio_cell_level + 1_iwp)  THEN
[3693]1552             IF ( mrt_include_sw )  THEN
[3735]1553                 tmrt_grid(j,i) = ( ( human_absorb * mrtinsw(l) +              &
[3740]1554                                  mrtinlw(l) )  /                              &
[3742]1555                                  ( human_emiss * sigma_sb ) )**.25_wp -       &
[3735]1556                                  degc_to_k
[3693]1557             ELSE
[3740]1558                 tmrt_grid(j,i) = ( mrtinlw(l)  /                              &
[3742]1559                                  ( human_emiss * sigma_sb ) )**.25_wp -       &
[3735]1560                                  degc_to_k
[3693]1561             ENDIF
[3525]1562          ENDIF
[3693]1563       ENDDO
1564    ENDIF
[3321]1565
[3525]1566END SUBROUTINE bio_calculate_mrt_grid
[3321]1567
1568
[3448]1569!------------------------------------------------------------------------------!
1570! Description:
1571! ------------
1572!> Calculate static thermal indices for 2D grid point i, j
1573!------------------------------------------------------------------------------!
[3525]1574 SUBROUTINE bio_get_thermal_index_input_ij( average_input, i, j, ta, vp, ws,   &
[4168]1575                                            pair, tmrt )
[3321]1576
[3448]1577    IMPLICIT NONE
[3569]1578!
[3448]1579!-- Input variables
1580    LOGICAL,      INTENT ( IN ) ::  average_input  !< Determine averaged input conditions?
1581    INTEGER(iwp), INTENT ( IN ) ::  i     !< Running index, x-dir
1582    INTEGER(iwp), INTENT ( IN ) ::  j     !< Running index, y-dir
[3569]1583!
[3448]1584!-- Output parameters
[3593]1585    REAL(wp), INTENT ( OUT )    ::  tmrt  !< Mean radiant temperature        (degree_C)
1586    REAL(wp), INTENT ( OUT )    ::  ta    !< Air temperature                 (degree_C)
[3448]1587    REAL(wp), INTENT ( OUT )    ::  vp    !< Vapour pressure                 (hPa)
1588    REAL(wp), INTENT ( OUT )    ::  ws    !< Wind speed    (local level)     (m/s)
1589    REAL(wp), INTENT ( OUT )    ::  pair  !< Air pressure                    (hPa)
[3569]1590!
[3448]1591!-- Internal variables
1592    INTEGER(iwp)                ::  k     !< Running index, z-dir
1593    INTEGER(iwp)                ::  k_wind  !< Running index, z-dir, wind speed only
[3321]1594
[3448]1595    REAL(wp)                    ::  vp_sat  !< Saturation vapor pressure     (hPa)
[3321]1596
[3569]1597!
[3448]1598!-- Determine cell level closest to 1.1m above ground
1599!   by making use of truncation due to int cast
[4168]1600    k = INT( topo_top_ind(j,i,0) + bio_cell_level )  !< Vertical cell center closest to 1.1m
[3525]1601
1602!
1603!-- Avoid non-representative horizontal u and v of 0.0 m/s too close to ground
[3448]1604    k_wind = k
[3742]1605    IF ( bio_cell_level < 1_iwp )  THEN
[3525]1606       k_wind = k + 1_iwp
[3448]1607    ENDIF
[3569]1608!
[3448]1609!-- Determine local values:
[3742]1610    IF ( average_input )  THEN
[3569]1611!
[3448]1612!--    Calculate ta from Tp assuming dry adiabatic laps rate
[3749]1613       ta = bio_fill_value
1614       IF ( ALLOCATED( pt_av ) )  THEN
[3753]1615          ta = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * ( k + .5_wp ) ) - degc_to_k
[3749]1616       ENDIF
[3321]1617
[3569]1618       vp = bio_fill_value
[3742]1619       IF ( humidity  .AND.  ALLOCATED( q_av ) )  THEN
1620          vp = q_av(k,j,i)
[3448]1621       ENDIF
[3321]1622
[3749]1623       ws = bio_fill_value
1624       IF ( ALLOCATED( u_av )  .AND.  ALLOCATED( v_av )  .AND.                 &
1625          ALLOCATED( w_av ) )  THEN
1626             ws = ( 0.5_wp * ABS( u_av(k_wind,j,i) + u_av(k_wind,j,i+1) )  +   &
1627             0.5_wp * ABS( v_av(k_wind,j,i) + v_av(k_wind,j+1,i) )  +          &
1628             0.5_wp * ABS( w_av(k_wind,j,i) + w_av(k_wind+1,j,i) ) )
1629       ENDIF
[3448]1630    ELSE
[3569]1631!
[3448]1632!-- Calculate ta from Tp assuming dry adiabatic laps rate
[3742]1633       ta = pt(k,j,i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
[3321]1634
[3569]1635       vp = bio_fill_value
[3742]1636       IF ( humidity )  THEN
1637          vp = q(k,j,i)
[3525]1638       ENDIF
[3321]1639
[3749]1640       ws = ( 0.5_wp * ABS( u(k_wind,j,i) + u(k_wind,j,i+1) )  +               &
1641          0.5_wp * ABS( v(k_wind,j,i) + v(k_wind,j+1,i) )  +                   &
1642          0.5_wp * ABS( w(k_wind,j,i) + w(k_wind+1,j,i) ) )
[3321]1643
[3448]1644    ENDIF
[3569]1645!
[3448]1646!-- Local air pressure
1647    pair = surface_pressure
1648!
1649!-- Calculate water vapour pressure at saturation and convert to hPa
1650!-- The magnus formula is limited to temperatures up to 333.15 K to
1651!   avoid negative values of vp_sat
[3742]1652    IF ( vp > -998._wp )  THEN
[3525]1653       vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) )
1654       vp  = vp * pair / ( vp + 0.622_wp )
[3742]1655       IF ( vp > vp_sat )  vp = vp_sat
[3525]1656    ENDIF
[3569]1657!
[3525]1658!-- local mtr value at [i,j]
[3569]1659    tmrt = bio_fill_value  !< this can be a valid result (e.g. for inside some ostacle)
[3742]1660    IF ( .NOT. average_input )  THEN
[3569]1661!
[3525]1662!--    Use MRT from RTM precalculated in tmrt_grid
1663       tmrt = tmrt_grid(j,i)
[3693]1664    ELSE
1665       tmrt = tmrt_av_grid(j,i)
[3448]1666    ENDIF
1667
[3525]1668 END SUBROUTINE bio_get_thermal_index_input_ij
[3448]1669
1670
[3321]1671!------------------------------------------------------------------------------!
1672! Description:
1673! ------------
[3448]1674!> Calculate static thermal indices for any point within a 2D grid
1675!> time_integration.f90: 1065ff
[3321]1676!------------------------------------------------------------------------------!
[3749]1677 SUBROUTINE bio_calculate_thermal_index_maps( av )
[3321]1678
[3448]1679    IMPLICIT NONE
[3569]1680!
[3448]1681!-- Input attributes
[3525]1682    LOGICAL, INTENT ( IN ) ::  av  !< Calculate based on averaged input conditions?
[3569]1683!
[3448]1684!-- Internal variables
[3569]1685    INTEGER(iwp) ::  i, j     !< Running index
[3321]1686
[3569]1687    REAL(wp) ::  clo          !< Clothing index                (no dimension)
[3593]1688    REAL(wp) ::  ta           !< Air temperature                  (degree_C)
[3569]1689    REAL(wp) ::  vp           !< Vapour pressure                  (hPa)
1690    REAL(wp) ::  ws           !< Wind speed    (local level)      (m/s)
1691    REAL(wp) ::  pair         !< Air pressure                     (hPa)
[3593]1692    REAL(wp) ::  perct_ij     !< Perceived temperature            (degree_C)
1693    REAL(wp) ::  utci_ij      !< Universal thermal climate index  (degree_C)
1694    REAL(wp) ::  pet_ij       !< Physiologically equivalent temperature  (degree_C)
1695    REAL(wp) ::  tmrt_ij      !< Mean radiant temperature         (degree_C)
[3321]1696
[3569]1697!
[3749]1698!-- Check if some thermal index is desired. Don't do anything if, e.g. only
1699!-- bio_mrt is desired.
1700    IF ( do_calculate_perct  .OR.  do_calculate_perct_av  .OR.                 &
1701       do_calculate_utci  .OR.  do_calculate_utci_av  .OR.                     &
[4127]1702       do_calculate_pet  .OR.  do_calculate_pet_av  .OR.                       &
1703       do_calculate_mrt2d )  THEN
[3321]1704
[3569]1705!
[3749]1706!--    fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...)
1707       CALL bio_calculate_mrt_grid ( av )
1708
1709       DO  i = nxl, nxr
1710          DO  j = nys, nyn
[3569]1711!
[3749]1712!--          Determine local input conditions
1713             tmrt_ij = bio_fill_value
1714             vp      = bio_fill_value
[3569]1715!
[3749]1716!--          Determine local meteorological conditions
1717             CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp,           &
1718                                                   ws, pair, tmrt_ij )
[3569]1719!
[3749]1720!--          Only proceed if input is available
1721             pet_ij   = bio_fill_value   !< set fail value, e.g. valid for
1722             perct_ij = bio_fill_value   !< within some obstacle
1723             utci_ij  = bio_fill_value
1724             IF ( .NOT. ( tmrt_ij <= -998._wp  .OR.  vp <= -998._wp  .OR.      &
1725                ws <= -998._wp  .OR.  ta <= -998._wp ) )  THEN
1726!
1727!--             Calculate static thermal indices based on local tmrt
1728                clo = bio_fill_value
[3569]1729
[3749]1730                IF ( do_calculate_perct  .OR.  do_calculate_perct_av )  THEN
[3569]1731!
[3749]1732!--                Estimate local perceived temperature
1733                   CALL calculate_perct_static( ta, vp, ws, tmrt_ij, pair,     &
1734                      clo, perct_ij )
1735                ENDIF
[3569]1736
[3749]1737                IF ( do_calculate_utci  .OR.  do_calculate_utci_av )  THEN
[3569]1738!
[3749]1739!--                Estimate local universal thermal climate index
1740                   CALL calculate_utci_static( ta, vp, ws, tmrt_ij,            &
1741                      bio_output_height, utci_ij )
1742                ENDIF
[3569]1743
[3749]1744                IF ( do_calculate_pet  .OR.  do_calculate_pet_av )  THEN
[3569]1745!
[3749]1746!--                Estimate local physiologically equivalent temperature
1747                   CALL calculate_pet_static( ta, vp, ws, tmrt_ij, pair,       &
1748                      pet_ij )
1749                ENDIF
[3525]1750             ENDIF
[3321]1751
1752
[3749]1753             IF ( av )  THEN
[3569]1754!
[3749]1755!--             Write results for selected averaged indices
1756                IF ( do_calculate_perct_av )  THEN
1757                   perct_av(j, i) = perct_ij
1758                ENDIF
1759                IF ( do_calculate_utci_av )  THEN
1760                   utci_av(j, i) = utci_ij
1761                ENDIF
1762                IF ( do_calculate_pet_av )  THEN
1763                   pet_av(j, i)  = pet_ij
1764                ENDIF
1765             ELSE
[3569]1766!
[3749]1767!--             Write result for selected indices
1768                IF ( do_calculate_perct )  THEN
1769                   perct(j, i) = perct_ij
1770                ENDIF
1771                IF ( do_calculate_utci )  THEN
1772                   utci(j, i) = utci_ij
1773                ENDIF
1774                IF ( do_calculate_pet )  THEN
1775                   pet(j, i)  = pet_ij
1776                ENDIF
[3742]1777             ENDIF
[3321]1778
[3749]1779          ENDDO
[3742]1780       ENDDO
[3749]1781    ENDIF
[3321]1782
[3525]1783 END SUBROUTINE bio_calculate_thermal_index_maps
[3321]1784
[3448]1785!------------------------------------------------------------------------------!
1786! Description:
1787! ------------
1788!> Calculate dynamic thermal indices (currently only iPT, but expandable)
1789!------------------------------------------------------------------------------!
[3525]1790 SUBROUTINE bio_calc_ipt( ta, vp, ws, pair, tmrt, dt, energy_storage,          &
[3448]1791    t_clo, clo, actlev, age, weight, height, work, sex, ipt )
[3321]1792
[3448]1793    IMPLICIT NONE
[3569]1794!
[3448]1795!-- Input parameters
[3593]1796    REAL(wp), INTENT ( IN )  ::  ta   !< Air temperature                  (degree_C)
[3448]1797    REAL(wp), INTENT ( IN )  ::  vp   !< Vapour pressure                  (hPa)
1798    REAL(wp), INTENT ( IN )  ::  ws   !< Wind speed    (local level)      (m/s)
1799    REAL(wp), INTENT ( IN )  ::  pair !< Air pressure                     (hPa)
[3593]1800    REAL(wp), INTENT ( IN )  ::  tmrt !< Mean radiant temperature         (degree_C)
[3448]1801    REAL(wp), INTENT ( IN )  ::  dt   !< Time past since last calculation (s)
1802    REAL(wp), INTENT ( IN )  ::  age  !< Age of agent                     (y)
1803    REAL(wp), INTENT ( IN )  ::  weight  !< Weight of agent               (Kg)
1804    REAL(wp), INTENT ( IN )  ::  height  !< Height of agent               (m)
1805    REAL(wp), INTENT ( IN )  ::  work    !< Mechanical workload of agent
1806                                         !  (without metabolism!)         (W)
1807    INTEGER(iwp), INTENT ( IN ) ::  sex  !< Sex of agent (1 = male, 2 = female)
[3569]1808!
[3448]1809!-- Both, input and output parameters
1810    Real(wp), INTENT ( INOUT )  ::  energy_storage    !< Energy storage   (W/m²)
[3593]1811    Real(wp), INTENT ( INOUT )  ::  t_clo   !< Clothing temperature       (degree_C)
[3448]1812    Real(wp), INTENT ( INOUT )  ::  clo     !< Current clothing in sulation
1813    Real(wp), INTENT ( INOUT )  ::  actlev  !< Individuals activity level
1814                                            !  per unit surface area      (W/m²)
[3569]1815!
[3448]1816!-- Output parameters
[3593]1817    REAL(wp), INTENT ( OUT ) ::  ipt    !< Instationary perceived temp.   (degree_C)
[3569]1818!
[3693]1819!-- return immediatelly if nothing to do!
[3742]1820    IF ( .NOT. thermal_comfort )  THEN
[3693]1821        RETURN
1822    ENDIF
1823!
[3448]1824!-- If clo equals the initial value, this is the initial call
[3742]1825    IF ( clo <= -998._wp )  THEN
[3569]1826!
[3448]1827!--    Initialize instationary perceived temperature with personalized
1828!      PT as an initial guess, set actlev and clo
[3749]1829       CALL ipt_init( age, weight, height, sex, work, actlev, clo,            &
[3448]1830          ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo,                   &
1831          ipt )
1832    ELSE
[3569]1833!
[3448]1834!--    Estimate local instatinoary perceived temperature
1835       CALL ipt_cycle ( ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo,     &
1836          clo, actlev, work, ipt )
1837    ENDIF
[3321]1838
[3525]1839 END SUBROUTINE bio_calc_ipt
[3448]1840
[3525]1841
[3650]1842
[3525]1843!------------------------------------------------------------------------------!
1844! Description:
1845! ------------
1846!> SUBROUTINE for calculating UTCI Temperature (UTCI)
1847!> computed by a 6th order approximation
[3693]1848!>
[3525]1849!> UTCI regression equation after
1850!> Bröde P, Fiala D, Blazejczyk K, Holmér I, Jendritzky G, Kampmann B, Tinz B,
1851!> Havenith G (2012) Deriving the operational procedure for the Universal Thermal
1852!> Climate Index (UTCI). International Journal of Biometeorology 56 (3):481-494.
1853!> doi:10.1007/s00484-011-0454-1
[3693]1854!>
[3525]1855!> original source available at:
1856!> www.utci.org
1857!------------------------------------------------------------------------------!
[3569]1858 SUBROUTINE calculate_utci_static( ta_in, vp, ws_hag, tmrt, hag, utci_ij )
[3525]1859
1860    IMPLICIT NONE
[3569]1861!
[3525]1862!-- Type of input of the argument list
[3593]1863    REAL(WP), INTENT ( IN )  ::  ta_in    !< Local air temperature (degree_C)
[3569]1864    REAL(WP), INTENT ( IN )  ::  vp       !< Loacl vapour pressure (hPa)
1865    REAL(WP), INTENT ( IN )  ::  ws_hag   !< Incident wind speed (m/s)
[3593]1866    REAL(WP), INTENT ( IN )  ::  tmrt     !< Local mean radiant temperature (degree_C)
[3569]1867    REAL(WP), INTENT ( IN )  ::  hag      !< Height of wind speed input (m)
1868!
[3525]1869!-- Type of output of the argument list
[3593]1870    REAL(wp), INTENT ( OUT ) ::  utci_ij  !< Universal Thermal Climate Index (degree_C)
[3525]1871
[3593]1872    REAL(WP) ::  ta           !< air temperature modified by offset (degree_C)
[3569]1873    REAL(WP) ::  pa           !< air pressure in kPa      (kPa)
[3593]1874    REAL(WP) ::  d_tmrt       !< delta-tmrt               (degree_C)
[3569]1875    REAL(WP) ::  va           !< wind speed at 10 m above ground level    (m/s)
[3593]1876    REAL(WP) ::  offset       !< utci deviation by ta cond. exceeded      (degree_C)
[3569]1877    REAL(WP) ::  part_ta      !< Air temperature related part of the regression
1878    REAL(WP) ::  ta2          !< 2 times ta
1879    REAL(WP) ::  ta3          !< 3 times ta
1880    REAL(WP) ::  ta4          !< 4 times ta
1881    REAL(WP) ::  ta5          !< 5 times ta
1882    REAL(WP) ::  ta6          !< 6 times ta
1883    REAL(WP) ::  part_va      !< Vapour pressure related part of the regression
1884    REAL(WP) ::  va2          !< 2 times va
1885    REAL(WP) ::  va3          !< 3 times va
1886    REAL(WP) ::  va4          !< 4 times va
1887    REAL(WP) ::  va5          !< 5 times va
1888    REAL(WP) ::  va6          !< 6 times va
1889    REAL(WP) ::  part_d_tmrt  !< Mean radiant temp. related part of the reg.
1890    REAL(WP) ::  d_tmrt2      !< 2 times d_tmrt
1891    REAL(WP) ::  d_tmrt3      !< 3 times d_tmrt
1892    REAL(WP) ::  d_tmrt4      !< 4 times d_tmrt
1893    REAL(WP) ::  d_tmrt5      !< 5 times d_tmrt
1894    REAL(WP) ::  d_tmrt6      !< 6 times d_tmrt
1895    REAL(WP) ::  part_pa      !< Air pressure related part of the regression
1896    REAL(WP) ::  pa2          !< 2 times pa
1897    REAL(WP) ::  pa3          !< 3 times pa
1898    REAL(WP) ::  pa4          !< 4 times pa
1899    REAL(WP) ::  pa5          !< 5 times pa
1900    REAL(WP) ::  pa6          !< 6 times pa
1901    REAL(WP) ::  part_pa2     !< Air pressure^2 related part of the regression
1902    REAL(WP) ::  part_pa3     !< Air pressure^3 related part of the regression
1903    REAL(WP) ::  part_pa46    !< Air pressure^4-6 related part of the regression
[3525]1904
[3569]1905!
[3525]1906!-- Initialize
1907    offset = 0._wp
1908    ta = ta_in
1909    d_tmrt = tmrt - ta_in
[3569]1910!
[3525]1911!-- Use vapour pressure in kpa
1912    pa = vp / 10.0_wp
[3569]1913!
[3525]1914!-- Wind altitude correction from hag to 10m after Broede et al. (2012), eq.3
[3693]1915!-- z(0) is set to 0.01 according to UTCI profile definition
[3525]1916    va = ws_hag *  log ( 10.0_wp / 0.01_wp ) / log ( hag / 0.01_wp )
[3569]1917!
[3525]1918!-- Check if input values in range after Broede et al. (2012)
[3742]1919    IF ( ( d_tmrt > 70._wp )  .OR.  ( d_tmrt < -30._wp )  .OR.                 &
1920       ( vp >= 50._wp ) )  THEN
[3569]1921       utci_ij = bio_fill_value
[3525]1922       RETURN
1923    ENDIF
[3569]1924!
[3525]1925!-- Apply eq. 2 in Broede et al. (2012) for ta out of bounds
[3742]1926    IF ( ta > 50._wp )  THEN
[3525]1927       offset = ta - 50._wp
1928       ta = 50._wp
1929    ENDIF
[3742]1930    IF ( ta < -50._wp )  THEN
[3525]1931       offset = ta + 50._wp
1932       ta = -50._wp
1933    ENDIF
[3569]1934!
[3525]1935!-- For routine application. For wind speeds and relative
[3693]1936!-- humidity values below 0.5 m/s or 5%, respectively, the
1937!-- user is advised to use the lower bounds for the calculations.
[3742]1938    IF ( va < 0.5_wp )  va = 0.5_wp
1939    IF ( va > 17._wp )  va = 17._wp
[3525]1940
[3569]1941!
1942!-- Pre-calculate multiples of input parameters to save time later
1943
1944    ta2 = ta  * ta
1945    ta3 = ta2 * ta
1946    ta4 = ta3 * ta
1947    ta5 = ta4 * ta
1948    ta6 = ta5 * ta
1949
1950    va2 = va  * va
1951    va3 = va2 * va
1952    va4 = va3 * va
1953    va5 = va4 * va
1954    va6 = va5 * va
1955
1956    d_tmrt2 = d_tmrt  * d_tmrt
1957    d_tmrt3 = d_tmrt2 * d_tmrt
1958    d_tmrt4 = d_tmrt3 * d_tmrt
1959    d_tmrt5 = d_tmrt4 * d_tmrt
1960    d_tmrt6 = d_tmrt5 * d_tmrt
1961
1962    pa2 = pa  * pa
1963    pa3 = pa2 * pa
1964    pa4 = pa3 * pa
1965    pa5 = pa4 * pa
1966    pa6 = pa5 * pa
1967
1968!
1969!-- Pre-calculate parts of the regression equation
1970    part_ta = (  6.07562052e-01_wp )       +                                   &
1971              ( -2.27712343e-02_wp ) * ta  +                                   &
1972              (  8.06470249e-04_wp ) * ta2 +                                   &
1973              ( -1.54271372e-04_wp ) * ta3 +                                   &
1974              ( -3.24651735e-06_wp ) * ta4 +                                   &
1975              (  7.32602852e-08_wp ) * ta5 +                                   &
1976              (  1.35959073e-09_wp ) * ta6
1977
1978    part_va = ( -2.25836520e+00_wp ) * va +                                    &
1979        (  8.80326035e-02_wp ) * ta  * va +                                    &
1980        (  2.16844454e-03_wp ) * ta2 * va +                                    &
1981        ( -1.53347087e-05_wp ) * ta3 * va +                                    &
1982        ( -5.72983704e-07_wp ) * ta4 * va +                                    &
1983        ( -2.55090145e-09_wp ) * ta5 * va +                                    &
1984        ( -7.51269505e-01_wp ) *       va2 +                                   &
1985        ( -4.08350271e-03_wp ) * ta  * va2 +                                   &
1986        ( -5.21670675e-05_wp ) * ta2 * va2 +                                   &
1987        (  1.94544667e-06_wp ) * ta3 * va2 +                                   &
1988        (  1.14099531e-08_wp ) * ta4 * va2 +                                   &
1989        (  1.58137256e-01_wp ) *       va3 +                                   &
1990        ( -6.57263143e-05_wp ) * ta  * va3 +                                   &
1991        (  2.22697524e-07_wp ) * ta2 * va3 +                                   &
1992        ( -4.16117031e-08_wp ) * ta3 * va3 +                                   &
1993        ( -1.27762753e-02_wp ) *       va4 +                                   &
1994        (  9.66891875e-06_wp ) * ta  * va4 +                                   &
1995        (  2.52785852e-09_wp ) * ta2 * va4 +                                   &
1996        (  4.56306672e-04_wp ) *       va5 +                                   &
1997        ( -1.74202546e-07_wp ) * ta  * va5 +                                   &
1998        ( -5.91491269e-06_wp ) * va6
1999
2000    part_d_tmrt = (  3.98374029e-01_wp ) *       d_tmrt +                      &
2001            (  1.83945314e-04_wp ) * ta  *       d_tmrt +                      &
2002            ( -1.73754510e-04_wp ) * ta2 *       d_tmrt +                      &
2003            ( -7.60781159e-07_wp ) * ta3 *       d_tmrt +                      &
2004            (  3.77830287e-08_wp ) * ta4 *       d_tmrt +                      &
2005            (  5.43079673e-10_wp ) * ta5 *       d_tmrt +                      &
2006            ( -2.00518269e-02_wp ) *       va  * d_tmrt +                      &
2007            (  8.92859837e-04_wp ) * ta  * va  * d_tmrt +                      &
2008            (  3.45433048e-06_wp ) * ta2 * va  * d_tmrt +                      &
2009            ( -3.77925774e-07_wp ) * ta3 * va  * d_tmrt +                      &
2010            ( -1.69699377e-09_wp ) * ta4 * va  * d_tmrt +                      &
2011            (  1.69992415e-04_wp ) *       va2 * d_tmrt +                      &
2012            ( -4.99204314e-05_wp ) * ta  * va2 * d_tmrt +                      &
2013            (  2.47417178e-07_wp ) * ta2 * va2 * d_tmrt +                      &
2014            (  1.07596466e-08_wp ) * ta3 * va2 * d_tmrt +                      &
2015            (  8.49242932e-05_wp ) *       va3 * d_tmrt +                      &
2016            (  1.35191328e-06_wp ) * ta  * va3 * d_tmrt +                      &
2017            ( -6.21531254e-09_wp ) * ta2 * va3 * d_tmrt +                      &
2018            ( -4.99410301e-06_wp ) * va4 *       d_tmrt +                      &
2019            ( -1.89489258e-08_wp ) * ta  * va4 * d_tmrt +                      &
2020            (  8.15300114e-08_wp ) *       va5 * d_tmrt +                      &
2021            (  7.55043090e-04_wp ) *             d_tmrt2 +                     &
2022            ( -5.65095215e-05_wp ) * ta  *       d_tmrt2 +                     &
2023            ( -4.52166564e-07_wp ) * ta2 *       d_tmrt2 +                     &
2024            (  2.46688878e-08_wp ) * ta3 *       d_tmrt2 +                     &
2025            (  2.42674348e-10_wp ) * ta4 *       d_tmrt2 +                     &
2026            (  1.54547250e-04_wp ) *       va  * d_tmrt2 +                     &
2027            (  5.24110970e-06_wp ) * ta  * va  * d_tmrt2 +                     &
2028            ( -8.75874982e-08_wp ) * ta2 * va  * d_tmrt2 +                     &
2029            ( -1.50743064e-09_wp ) * ta3 * va  * d_tmrt2 +                     &
2030            ( -1.56236307e-05_wp ) *       va2 * d_tmrt2 +                     &
2031            ( -1.33895614e-07_wp ) * ta  * va2 * d_tmrt2 +                     &
2032            (  2.49709824e-09_wp ) * ta2 * va2 * d_tmrt2 +                     &
2033            (  6.51711721e-07_wp ) *       va3 * d_tmrt2 +                     &
2034            (  1.94960053e-09_wp ) * ta  * va3 * d_tmrt2 +                     &
2035            ( -1.00361113e-08_wp ) *       va4 * d_tmrt2 +                     &
2036            ( -1.21206673e-05_wp ) *             d_tmrt3 +                     &
2037            ( -2.18203660e-07_wp ) * ta  *       d_tmrt3 +                     &
2038            (  7.51269482e-09_wp ) * ta2 *       d_tmrt3 +                     &
2039            (  9.79063848e-11_wp ) * ta3 *       d_tmrt3 +                     &
2040            (  1.25006734e-06_wp ) *       va  * d_tmrt3 +                     &
2041            ( -1.81584736e-09_wp ) * ta  * va  * d_tmrt3 +                     &
2042            ( -3.52197671e-10_wp ) * ta2 * va  * d_tmrt3 +                     &
2043            ( -3.36514630e-08_wp ) *       va2 * d_tmrt3 +                     &
2044            (  1.35908359e-10_wp ) * ta  * va2 * d_tmrt3 +                     &
2045            (  4.17032620e-10_wp ) *       va3 * d_tmrt3 +                     &
2046            ( -1.30369025e-09_wp ) *             d_tmrt4 +                     &
2047            (  4.13908461e-10_wp ) * ta  *       d_tmrt4 +                     &
2048            (  9.22652254e-12_wp ) * ta2 *       d_tmrt4 +                     &
2049            ( -5.08220384e-09_wp ) *       va  * d_tmrt4 +                     &
2050            ( -2.24730961e-11_wp ) * ta  * va  * d_tmrt4 +                     &
2051            (  1.17139133e-10_wp ) *       va2 * d_tmrt4 +                     &
2052            (  6.62154879e-10_wp ) *             d_tmrt5 +                     &
2053            (  4.03863260e-13_wp ) * ta  *       d_tmrt5 +                     &
2054            (  1.95087203e-12_wp ) *       va  * d_tmrt5 +                     &
2055            ( -4.73602469e-12_wp ) *             d_tmrt6
2056
2057    part_pa = (  5.12733497e+00_wp ) *                pa +                     &
2058       ( -3.12788561e-01_wp ) * ta  *                 pa +                     &
2059       ( -1.96701861e-02_wp ) * ta2 *                 pa +                     &
2060       (  9.99690870e-04_wp ) * ta3 *                 pa +                     &
2061       (  9.51738512e-06_wp ) * ta4 *                 pa +                     &
2062       ( -4.66426341e-07_wp ) * ta5 *                 pa +                     &
2063       (  5.48050612e-01_wp ) *       va  *           pa +                     &
2064       ( -3.30552823e-03_wp ) * ta  * va  *           pa +                     &
2065       ( -1.64119440e-03_wp ) * ta2 * va  *           pa +                     &
2066       ( -5.16670694e-06_wp ) * ta3 * va  *           pa +                     &
2067       (  9.52692432e-07_wp ) * ta4 * va  *           pa +                     &
2068       ( -4.29223622e-02_wp ) *       va2 *           pa +                     &
2069       (  5.00845667e-03_wp ) * ta  * va2 *           pa +                     &
2070       (  1.00601257e-06_wp ) * ta2 * va2 *           pa +                     &
2071       ( -1.81748644e-06_wp ) * ta3 * va2 *           pa +                     &
2072       ( -1.25813502e-03_wp ) *       va3 *           pa +                     &
2073       ( -1.79330391e-04_wp ) * ta  * va3 *           pa +                     &
2074       (  2.34994441e-06_wp ) * ta2 * va3 *           pa +                     &
2075       (  1.29735808e-04_wp ) *       va4 *           pa +                     &
2076       (  1.29064870e-06_wp ) * ta  * va4 *           pa +                     &
2077       ( -2.28558686e-06_wp ) *       va5 *           pa +                     &
2078       ( -3.69476348e-02_wp ) *             d_tmrt  * pa +                     &
2079       (  1.62325322e-03_wp ) * ta  *       d_tmrt  * pa +                     &
2080       ( -3.14279680e-05_wp ) * ta2 *       d_tmrt  * pa +                     &
2081       (  2.59835559e-06_wp ) * ta3 *       d_tmrt  * pa +                     &
2082       ( -4.77136523e-08_wp ) * ta4 *       d_tmrt  * pa +                     &
2083       (  8.64203390e-03_wp ) *       va  * d_tmrt  * pa +                     &
2084       ( -6.87405181e-04_wp ) * ta  * va  * d_tmrt  * pa +                     &
2085       ( -9.13863872e-06_wp ) * ta2 * va  * d_tmrt  * pa +                     &
2086       (  5.15916806e-07_wp ) * ta3 * va  * d_tmrt  * pa +                     &
2087       ( -3.59217476e-05_wp ) *       va2 * d_tmrt  * pa +                     &
2088       (  3.28696511e-05_wp ) * ta  * va2 * d_tmrt  * pa +                     &
2089       ( -7.10542454e-07_wp ) * ta2 * va2 * d_tmrt  * pa +                     &
2090       ( -1.24382300e-05_wp ) *       va3 * d_tmrt  * pa +                     &
2091       ( -7.38584400e-09_wp ) * ta  * va3 * d_tmrt  * pa +                     &
2092       (  2.20609296e-07_wp ) *       va4 * d_tmrt  * pa +                     &
2093       ( -7.32469180e-04_wp ) *             d_tmrt2 * pa +                     &
2094       ( -1.87381964e-05_wp ) * ta  *       d_tmrt2 * pa +                     &
2095       (  4.80925239e-06_wp ) * ta2 *       d_tmrt2 * pa +                     &
2096       ( -8.75492040e-08_wp ) * ta3 *       d_tmrt2 * pa +                     &
2097       (  2.77862930e-05_wp ) *       va  * d_tmrt2 * pa +                     &
2098       ( -5.06004592e-06_wp ) * ta  * va  * d_tmrt2 * pa +                     &
2099       (  1.14325367e-07_wp ) * ta2 * va  * d_tmrt2 * pa +                     &
2100       (  2.53016723e-06_wp ) *       va2 * d_tmrt2 * pa +                     &
2101       ( -1.72857035e-08_wp ) * ta  * va2 * d_tmrt2 * pa +                     &
2102       ( -3.95079398e-08_wp ) *       va3 * d_tmrt2 * pa +                     &
2103       ( -3.59413173e-07_wp ) *             d_tmrt3 * pa +                     &
2104       (  7.04388046e-07_wp ) * ta  *       d_tmrt3 * pa +                     &
2105       ( -1.89309167e-08_wp ) * ta2 *       d_tmrt3 * pa +                     &
2106       ( -4.79768731e-07_wp ) *       va  * d_tmrt3 * pa +                     &
2107       (  7.96079978e-09_wp ) * ta  * va  * d_tmrt3 * pa +                     &
2108       (  1.62897058e-09_wp ) *       va2 * d_tmrt3 * pa +                     &
2109       (  3.94367674e-08_wp ) *             d_tmrt4 * pa +                     &
2110       ( -1.18566247e-09_wp ) * ta *        d_tmrt4 * pa +                     &
2111       (  3.34678041e-10_wp ) *       va  * d_tmrt4 * pa +                     &
2112       ( -1.15606447e-10_wp ) *             d_tmrt5 * pa
2113
2114    part_pa2 = ( -2.80626406e+00_wp ) *               pa2 +                    &
2115       (  5.48712484e-01_wp ) * ta  *                 pa2 +                    &
2116       ( -3.99428410e-03_wp ) * ta2 *                 pa2 +                    &
2117       ( -9.54009191e-04_wp ) * ta3 *                 pa2 +                    &
2118       (  1.93090978e-05_wp ) * ta4 *                 pa2 +                    &
2119       ( -3.08806365e-01_wp ) *       va *            pa2 +                    &
2120       (  1.16952364e-02_wp ) * ta  * va *            pa2 +                    &
2121       (  4.95271903e-04_wp ) * ta2 * va *            pa2 +                    &
2122       ( -1.90710882e-05_wp ) * ta3 * va *            pa2 +                    &
2123       (  2.10787756e-03_wp ) *       va2 *           pa2 +                    &
2124       ( -6.98445738e-04_wp ) * ta  * va2 *           pa2 +                    &
2125       (  2.30109073e-05_wp ) * ta2 * va2 *           pa2 +                    &
2126       (  4.17856590e-04_wp ) *       va3 *           pa2 +                    &
2127       ( -1.27043871e-05_wp ) * ta  * va3 *           pa2 +                    &
2128       ( -3.04620472e-06_wp ) *       va4 *           pa2 +                    &
2129       (  5.14507424e-02_wp ) *             d_tmrt  * pa2 +                    &
2130       ( -4.32510997e-03_wp ) * ta  *       d_tmrt  * pa2 +                    &
2131       (  8.99281156e-05_wp ) * ta2 *       d_tmrt  * pa2 +                    &
2132       ( -7.14663943e-07_wp ) * ta3 *       d_tmrt  * pa2 +                    &
2133       ( -2.66016305e-04_wp ) *       va  * d_tmrt  * pa2 +                    &
2134       (  2.63789586e-04_wp ) * ta  * va  * d_tmrt  * pa2 +                    &
2135       ( -7.01199003e-06_wp ) * ta2 * va  * d_tmrt  * pa2 +                    &
2136       ( -1.06823306e-04_wp ) *       va2 * d_tmrt  * pa2 +                    &
2137       (  3.61341136e-06_wp ) * ta  * va2 * d_tmrt  * pa2 +                    &
2138       (  2.29748967e-07_wp ) *       va3 * d_tmrt  * pa2 +                    &
2139       (  3.04788893e-04_wp ) *             d_tmrt2 * pa2 +                    &
2140       ( -6.42070836e-05_wp ) * ta  *       d_tmrt2 * pa2 +                    &
2141       (  1.16257971e-06_wp ) * ta2 *       d_tmrt2 * pa2 +                    &
2142       (  7.68023384e-06_wp ) *       va  * d_tmrt2 * pa2 +                    &
2143       ( -5.47446896e-07_wp ) * ta  * va  * d_tmrt2 * pa2 +                    &
2144       ( -3.59937910e-08_wp ) *       va2 * d_tmrt2 * pa2 +                    &
2145       ( -4.36497725e-06_wp ) *             d_tmrt3 * pa2 +                    &
2146       (  1.68737969e-07_wp ) * ta  *       d_tmrt3 * pa2 +                    &
2147       (  2.67489271e-08_wp ) *       va  * d_tmrt3 * pa2 +                    &
2148       (  3.23926897e-09_wp ) *             d_tmrt4 * pa2
2149
2150    part_pa3 = ( -3.53874123e-02_wp ) *               pa3 +                    &
2151       ( -2.21201190e-01_wp ) * ta  *                 pa3 +                    &
2152       (  1.55126038e-02_wp ) * ta2 *                 pa3 +                    &
2153       ( -2.63917279e-04_wp ) * ta3 *                 pa3 +                    &
2154       (  4.53433455e-02_wp ) *       va  *           pa3 +                    &
2155       ( -4.32943862e-03_wp ) * ta  * va  *           pa3 +                    &
2156       (  1.45389826e-04_wp ) * ta2 * va  *           pa3 +                    &
2157       (  2.17508610e-04_wp ) *       va2 *           pa3 +                    &
2158       ( -6.66724702e-05_wp ) * ta  * va2 *           pa3 +                    &
2159       (  3.33217140e-05_wp ) *       va3 *           pa3 +                    &
2160       ( -2.26921615e-03_wp ) *             d_tmrt  * pa3 +                    &
2161       (  3.80261982e-04_wp ) * ta  *       d_tmrt  * pa3 +                    &
2162       ( -5.45314314e-09_wp ) * ta2 *       d_tmrt  * pa3 +                    &
2163       ( -7.96355448e-04_wp ) *       va  * d_tmrt  * pa3 +                    &
2164       (  2.53458034e-05_wp ) * ta  * va  * d_tmrt  * pa3 +                    &
2165       ( -6.31223658e-06_wp ) *       va2 * d_tmrt  * pa3 +                    &
2166       (  3.02122035e-04_wp ) *             d_tmrt2 * pa3 +                    &
2167       ( -4.77403547e-06_wp ) * ta  *       d_tmrt2 * pa3 +                    &
2168       (  1.73825715e-06_wp ) *       va  * d_tmrt2 * pa3 +                    &
2169       ( -4.09087898e-07_wp ) *             d_tmrt3 * pa3
2170
2171    part_pa46 = (  6.14155345e-01_wp ) *              pa4 +                    &
2172       ( -6.16755931e-02_wp ) * ta  *                 pa4 +                    &
2173       (  1.33374846e-03_wp ) * ta2 *                 pa4 +                    &
2174       (  3.55375387e-03_wp ) *       va  *           pa4 +                    &
[3693]2175       ( -5.13027851e-04_wp ) * ta  * va  *           pa4 +                    &
[3569]2176       (  1.02449757e-04_wp ) *       va2 *           pa4 +                    &
2177       ( -1.48526421e-03_wp ) *             d_tmrt  * pa4 +                    &
2178       ( -4.11469183e-05_wp ) * ta  *       d_tmrt  * pa4 +                    &
2179       ( -6.80434415e-06_wp ) *       va  * d_tmrt  * pa4 +                    &
2180       ( -9.77675906e-06_wp ) *             d_tmrt2 * pa4 +                    &
2181       (  8.82773108e-02_wp ) *                       pa5 +                    &
2182       ( -3.01859306e-03_wp ) * ta  *                 pa5 +                    &
2183       (  1.04452989e-03_wp ) *       va  *           pa5 +                    &
2184       (  2.47090539e-04_wp ) *             d_tmrt  * pa5 +                    &
2185       (  1.48348065e-03_wp ) *                       pa6
2186!
[3525]2187!-- Calculate 6th order polynomial as approximation
[3569]2188    utci_ij = ta + part_ta + part_va + part_d_tmrt + part_pa + part_pa2 +      &
2189        part_pa3 + part_pa46
2190!
[3525]2191!-- Consider offset in result
[3569]2192    utci_ij = utci_ij + offset
[3525]2193
2194 END SUBROUTINE calculate_utci_static
2195
2196
2197
2198
2199!------------------------------------------------------------------------------!
2200! Description:
2201! ------------
[3593]2202!> calculate_perct_static: Estimation of perceived temperature (PT, degree_C)
[3525]2203!> Value of perct is the Perceived Temperature, degree centigrade
2204!------------------------------------------------------------------------------!
[3569]2205 SUBROUTINE calculate_perct_static( ta, vp, ws, tmrt, pair, clo, perct_ij )
[3525]2206
2207    IMPLICIT NONE
[3569]2208!
[3525]2209!-- Type of input of the argument list
2210    REAL(wp), INTENT ( IN )  :: ta   !< Local air temperature (degC)
2211    REAL(wp), INTENT ( IN )  :: vp   !< Local vapour pressure (hPa)
2212    REAL(wp), INTENT ( IN )  :: tmrt !< Local mean radiant temperature (degC)
2213    REAL(wp), INTENT ( IN )  :: ws   !< Local wind velocitry (m/s)
2214    REAL(wp), INTENT ( IN )  :: pair !< Local barometric air pressure (hPa)
[3569]2215!
[3525]2216!-- Type of output of the argument list
[3569]2217    REAL(wp), INTENT ( OUT ) :: perct_ij  !< Perceived temperature (degC)
2218    REAL(wp), INTENT ( OUT ) :: clo       !< Clothing index (dimensionless)
2219!
[3525]2220!-- Parameters for standard "Klima-Michel"
[3693]2221    REAL(wp), PARAMETER :: eta = 0._wp  !< Mechanical work efficiency for walking on flat ground
2222                                        !< (compare to Fanger (1972) pp 24f)
2223    REAL(wp), PARAMETER :: actlev = 134.6862_wp  !< Workload by activity per standardized surface (A_Du)
[3569]2224!
[3525]2225!-- Type of program variables
2226    REAL(wp), PARAMETER :: eps = 0.0005  !< Accuracy in clothing insulation (clo) for evaluation the root of Fanger's PMV (pmva=0)
2227    REAL(wp) ::  sclo           !< summer clothing insulation
2228    REAL(wp) ::  wclo           !< winter clothing insulation
2229    REAL(wp) ::  d_pmv          !< PMV deviation (dimensionless --> PMV)
2230    REAL(wp) ::  svp_ta         !< saturation vapor pressure    (hPa)
2231    REAL(wp) ::  sult_lim       !< threshold for sultrieness    (hPa)
2232    REAL(wp) ::  dgtcm          !< Mean deviation dependent on perct
2233    REAL(wp) ::  dgtcstd        !< Mean deviation plus its standard deviation
2234    REAL(wp) ::  clon           !< clo for neutral conditions   (clo)
2235    REAL(wp) ::  ireq_minimal   !< Minimal required clothing insulation (clo)
2236    REAL(wp) ::  pmv_w          !< Fangers predicted mean vote for winter clothing
2237    REAL(wp) ::  pmv_s          !< Fangers predicted mean vote for summer clothing
2238    REAL(wp) ::  pmva           !< adjusted predicted mean vote
[3593]2239    REAL(wp) ::  ptc            !< perceived temp. for cold conditions (degree_C)
[3525]2240    REAL(wp) ::  d_std          !< factor to threshold for sultriness
2241    REAL(wp) ::  pmvs           !< pred. mean vote considering sultrieness
2242
2243    INTEGER(iwp) :: ncount      !< running index
2244    INTEGER(iwp) :: nerr_cold   !< error number (cold conditions)
2245    INTEGER(iwp) :: nerr        !< error number
2246
2247    LOGICAL :: sultrieness
[3569]2248!
[3525]2249!-- Initialise
[3569]2250    perct_ij = bio_fill_value
[3525]2251
2252    nerr     = 0_iwp
2253    ncount   = 0_iwp
2254    sultrieness  = .FALSE.
[3569]2255!
[3525]2256!-- Tresholds: clothing insulation (account for model inaccuracies)
[3569]2257!
[3693]2258!-- summer clothing
[3525]2259    sclo     = 0.44453_wp
[3569]2260!
[3693]2261!-- winter clothing
[3525]2262    wclo     = 1.76267_wp
[3569]2263!
[3525]2264!-- decision: firstly calculate for winter or summer clothing
[3742]2265    IF ( ta <= 10._wp )  THEN
[3569]2266!
[3525]2267!--    First guess: winter clothing insulation: cold stress
2268       clo = wclo
[3753]2269       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
[3525]2270       pmv_w = pmva
2271
[3742]2272       IF ( pmva > 0._wp )  THEN
[3569]2273!
[3525]2274!--       Case summer clothing insulation: heat load ?
2275          clo = sclo
[3753]2276          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
[3525]2277          pmv_s = pmva
[3742]2278          IF ( pmva <= 0._wp )  THEN
[3569]2279!
[3525]2280!--          Case: comfort achievable by varying clothing insulation
[3693]2281!--          Between winter and summer set values
[3525]2282             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
[3753]2283                pmv_s, wclo, pmv_w, eps, pmva, ncount, clo )
[3742]2284             IF ( ncount < 0_iwp )  THEN
[3525]2285                nerr = -1_iwp
2286                RETURN
2287             ENDIF
[3742]2288          ELSE IF ( pmva > 0.06_wp )  THEN
[3525]2289             clo = 0.5_wp
2290             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
[3753]2291                           pmva )
[3525]2292          ENDIF
[3742]2293       ELSE IF ( pmva < -0.11_wp )  THEN
[3525]2294          clo = 1.75_wp
[3753]2295          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
[3525]2296       ENDIF
2297    ELSE
[3569]2298!
[3525]2299!--    First guess: summer clothing insulation: heat load
2300       clo = sclo
[3753]2301       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
[3525]2302       pmv_s = pmva
2303
[3742]2304       IF ( pmva < 0._wp )  THEN
[3569]2305!
[3525]2306!--       Case winter clothing insulation: cold stress ?
2307          clo = wclo
[3753]2308          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
[3525]2309          pmv_w = pmva
2310
[3742]2311          IF ( pmva >= 0._wp )  THEN
[3569]2312!
[3525]2313!--          Case: comfort achievable by varying clothing insulation
[3693]2314!--          between winter and summer set values
[3525]2315             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
[3753]2316                               pmv_s, wclo, pmv_w, eps, pmva, ncount, clo )
[3742]2317             IF ( ncount < 0_iwp )  THEN
[3525]2318                nerr = -1_iwp
2319                RETURN
2320             ENDIF
[3742]2321          ELSE IF ( pmva < -0.11_wp )  THEN
[3525]2322             clo = 1.75_wp
2323             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
[3753]2324                           pmva )
[3525]2325          ENDIF
[3742]2326       ELSE IF ( pmva > 0.06_wp )  THEN
[3525]2327          clo = 0.5_wp
[3753]2328          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
[3525]2329       ENDIF
2330
2331    ENDIF
[3569]2332!
[3525]2333!-- Determine perceived temperature by regression equation + adjustments
2334    pmvs = pmva
[3749]2335    CALL perct_regression( pmva, clo, perct_ij )
[3569]2336    ptc = perct_ij
[3742]2337    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
[3569]2338!
[3525]2339!--    Adjust for cold conditions according to Gagge 1986
2340       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
[3742]2341       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
[3525]2342       pmvs = pmva - d_pmv
[3742]2343       IF ( pmvs > -0.11_wp )  THEN
[3525]2344          d_pmv  = 0._wp
2345          pmvs   = -0.11_wp
2346       ENDIF
[3749]2347       CALL perct_regression( pmvs, clo, perct_ij )
[3525]2348    ENDIF
2349!     clo_fanger = clo
2350    clon = clo
[3742]2351    IF ( clo > 0.5_wp  .AND.  perct_ij <= 8.73_wp )  THEN
[3569]2352!
[3525]2353!--    Required clothing insulation (ireq) is exclusively defined for
[3753]2354!--    perceived temperatures (perct) less 10 (C) for a
[3693]2355!--    reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
[3569]2356       clon = ireq_neutral ( perct_ij, ireq_minimal, nerr )
[3525]2357       clo = clon
2358    ENDIF
[3749]2359    CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
[3525]2360    sultrieness    = .FALSE.
2361    d_std = -99._wp
[3742]2362    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
[3569]2363!
[3525]2364!--    Adjust for warm/humid conditions according to Gagge 1986
2365       CALL saturation_vapor_pressure ( ta, svp_ta )
2366       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
2367       pmvs   = pmva + d_pmv
[3749]2368       CALL perct_regression( pmvs, clo, perct_ij )
[3742]2369       IF ( sult_lim < 99._wp )  THEN
2370          IF ( (perct_ij - ptc) > sult_lim )  sultrieness = .TRUE.
[3569]2371!
[3525]2372!--       Set factor to threshold for sultriness
[3749]2373          IF ( ABS( dgtcstd ) > 0.00001_wp )  THEN
[3569]2374             d_std = ( ( perct_ij - ptc ) - dgtcm ) / dgtcstd
[3525]2375          ENDIF
2376       ENDIF
2377    ENDIF
2378
2379 END SUBROUTINE calculate_perct_static
2380
2381!------------------------------------------------------------------------------!
2382! Description:
2383! ------------
[3693]2384!> The SUBROUTINE calculates the (saturation) water vapour pressure
[3525]2385!> (hPa = hecto Pascal) for a given temperature ta (degC).
[3693]2386!> 'ta' can be the air temperature or the dew point temperature. The first will
2387!> result in the current vapor pressure (hPa), the latter will calulate the
2388!> saturation vapor pressure (hPa).
[3525]2389!------------------------------------------------------------------------------!
2390 SUBROUTINE saturation_vapor_pressure( ta, svp_ta )
2391
2392    IMPLICIT NONE
2393
2394    REAL(wp), INTENT ( IN )  ::  ta     !< ambient air temperature (degC)
[3693]2395    REAL(wp), INTENT ( OUT ) ::  svp_ta !< water vapour pressure (hPa)
[3525]2396
2397    REAL(wp)      ::  b
2398    REAL(wp)      ::  c
2399
2400
[3742]2401    IF ( ta < 0._wp )  THEN
[3569]2402!
[3693]2403!--    ta  < 0 (degC): water vapour pressure over ice
[3525]2404       b = 17.84362_wp
2405       c = 245.425_wp
2406    ELSE
[3569]2407!
[3693]2408!--    ta >= 0 (degC): water vapour pressure over water
[3525]2409       b = 17.08085_wp
2410       c = 234.175_wp
2411    ENDIF
[3569]2412!
[3525]2413!-- Saturation water vapour pressure
[3753]2414    svp_ta = 6.1078_wp * EXP( b * ta / ( c + ta ) )
[3525]2415
2416 END SUBROUTINE saturation_vapor_pressure
2417
2418!------------------------------------------------------------------------------!
2419! Description:
2420! ------------
2421!> Find the clothing insulation value clo_res (clo) to make Fanger's Predicted
2422!> Mean Vote (PMV) equal comfort (pmva=0) for actual meteorological conditions
2423!> (ta,tmrt, vp, ws, pair) and values of individual's activity level
2424!------------------------------------------------------------------------------!
2425 SUBROUTINE iso_ridder( ta, tmrt, vp, ws, pair, actlev, eta, sclo,             &
[3753]2426                       pmv_s, wclo, pmv_w, eps, pmva, nerr,               &
[3525]2427                       clo_res )
2428
2429    IMPLICIT NONE
[3569]2430!
[3525]2431!-- Input variables of argument list:
[3693]2432    REAL(wp), INTENT ( IN )  :: ta       !< Ambient temperature (degC)
[3525]2433    REAL(wp), INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC)
[3693]2434    REAL(wp), INTENT ( IN )  :: vp       !< Water vapour pressure (hPa)
2435    REAL(wp), INTENT ( IN )  :: ws       !< Wind speed (m/s) 1 m above ground
2436    REAL(wp), INTENT ( IN )  :: pair     !< Barometric air pressure (hPa)
[3525]2437    REAL(wp), INTENT ( IN )  :: actlev   !< Individuals activity level per unit surface area (W/m2)
2438    REAL(wp), INTENT ( IN )  :: eta      !< Individuals work efficiency (dimensionless)
2439    REAL(wp), INTENT ( IN )  :: sclo     !< Lower threshold of bracketing clothing insulation (clo)
2440    REAL(wp), INTENT ( IN )  :: wclo     !< Upper threshold of bracketing clothing insulation (clo)
2441    REAL(wp), INTENT ( IN )  :: eps      !< (0.05) accuracy in clothing insulation (clo) for
2442!                                          evaluation the root of Fanger's PMV (pmva=0)
2443    REAL(wp), INTENT ( IN )  :: pmv_w    !< Fanger's PMV corresponding to wclo
2444    REAL(wp), INTENT ( IN )  :: pmv_s    !< Fanger's PMV corresponding to sclo
[3569]2445!
[3693]2446!-- Output variables of argument list:
[3525]2447    REAL(wp), INTENT ( OUT ) :: pmva     !< 0 (set to zero, because clo is evaluated for comfort)
2448    REAL(wp), INTENT ( OUT ) :: clo_res  !< Resulting clothing insulation value (clo)
2449    INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error status / quality flag
[3693]2450                                         !< nerr >= 0, o.k., and nerr is the number of iterations for convergence
2451                                         !< nerr = -1: error = malfunction of Ridder's convergence method
2452                                         !< nerr = -2: error = maximum iterations (max_iteration) exceeded
2453                                         !< nerr = -3: error = root not bracketed between sclo and wclo
[3569]2454!
[3525]2455!-- Type of program variables
2456    INTEGER(iwp), PARAMETER  ::  max_iteration = 15_iwp       !< max number of iterations
2457    REAL(wp),     PARAMETER  ::  guess_0       = -1.11e30_wp  !< initial guess
2458    REAL(wp) ::  x_ridder    !< current guess for clothing insulation   (clo)
2459    REAL(wp) ::  clo_lower   !< lower limit of clothing insulation      (clo)
2460    REAL(wp) ::  clo_upper   !< upper limit of clothing insulation      (clo)
2461    REAL(wp) ::  x_lower     !< lower guess for clothing insulation     (clo)
2462    REAL(wp) ::  x_upper     !< upper guess for clothing insulation     (clo)
2463    REAL(wp) ::  x_average   !< average of x_lower and x_upper          (clo)
2464    REAL(wp) ::  x_new       !< preliminary result for clothing insulation (clo)
2465    REAL(wp) ::  y_lower     !< predicted mean vote for summer clothing
2466    REAL(wp) ::  y_upper     !< predicted mean vote for winter clothing
2467    REAL(wp) ::  y_average   !< average of y_lower and y_upper
2468    REAL(wp) ::  y_new       !< preliminary result for pred. mean vote
2469    REAL(wp) ::  sroot       !< sqrt of PMV-guess
2470    INTEGER(iwp) ::  j       !< running index
[3569]2471!
[3525]2472!-- Initialise
2473    nerr    = 0_iwp
[3569]2474!
[3525]2475!-- Set pmva = 0 (comfort): Root of PMV depending on clothing insulation
[3569]2476    x_ridder    = bio_fill_value
[3525]2477    pmva        = 0._wp
2478    clo_lower   = sclo
2479    y_lower     = pmv_s
2480    clo_upper   = wclo
2481    y_upper     = pmv_w
[3742]2482    IF ( ( y_lower > 0._wp  .AND.  y_upper < 0._wp )  .OR.                     &
2483         ( y_lower < 0._wp  .AND.  y_upper > 0._wp ) )  THEN
[3525]2484       x_lower  = clo_lower
2485       x_upper  = clo_upper
2486       x_ridder = guess_0
2487
[3742]2488       DO  j = 1_iwp, max_iteration
[3525]2489          x_average = 0.5_wp * ( x_lower + x_upper )
2490          CALL fanger ( ta, tmrt, vp, ws, pair, x_average, actlev, eta,        &
[3753]2491                        y_average )
2492          sroot = SQRT( y_average**2 - y_lower * y_upper )
[3749]2493          IF ( ABS( sroot ) < 0.00001_wp )  THEN
[3525]2494             clo_res = x_average
2495             nerr = j
2496             RETURN
2497          ENDIF
2498          x_new = x_average + ( x_average - x_lower ) *                        &
2499                      ( SIGN ( 1._wp, y_lower - y_upper ) * y_average / sroot )
[3753]2500          IF ( ABS( x_new - x_ridder ) <= eps )  THEN
[3525]2501             clo_res = x_ridder
2502             nerr       = j
2503             RETURN
2504          ENDIF
2505          x_ridder = x_new
2506          CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta,         &
[3753]2507                        y_new )
[3749]2508          IF ( ABS( y_new ) < 0.00001_wp )  THEN
[3525]2509             clo_res = x_ridder
2510             nerr       = j
2511             RETURN
2512          ENDIF
[3749]2513          IF ( ABS( SIGN( y_average, y_new ) - y_average ) > 0.00001_wp )  THEN
[3525]2514             x_lower = x_average
2515             y_lower = y_average
2516             x_upper  = x_ridder
2517             y_upper  = y_new
[3749]2518          ELSE IF ( ABS( SIGN( y_lower, y_new ) - y_lower ) > 0.00001_wp )  THEN
[3525]2519             x_upper  = x_ridder
2520             y_upper  = y_new
[3749]2521          ELSE IF ( ABS( SIGN( y_upper, y_new ) - y_upper ) > 0.00001_wp )  THEN
[3525]2522             x_lower = x_ridder
2523             y_lower = y_new
2524          ELSE
[3569]2525!
[3525]2526!--          Never get here in x_ridder: singularity in y
[3749]2527             nerr    = -1_iwp
[3525]2528             clo_res = x_ridder
2529             RETURN
2530          ENDIF
[3753]2531          IF ( ABS( x_upper - x_lower ) <= eps )  THEN
[3525]2532             clo_res = x_ridder
[3749]2533             nerr    = j
[3525]2534             RETURN
2535          ENDIF
2536       ENDDO
[3569]2537!
[3525]2538!--    x_ridder exceed maximum iterations
2539       nerr       = -2_iwp
2540       clo_res = y_new
2541       RETURN
[3749]2542    ELSE IF ( ABS( y_lower ) < 0.00001_wp )  THEN
[3525]2543       x_ridder = clo_lower
[3749]2544    ELSE IF ( ABS( y_upper ) < 0.00001_wp )  THEN
[3525]2545       x_ridder = clo_upper
2546    ELSE
[3569]2547!
[3525]2548!--    x_ridder not bracketed by u_clo and o_clo
2549       nerr = -3_iwp
2550       clo_res = x_ridder
2551       RETURN
2552    ENDIF
2553
2554 END SUBROUTINE iso_ridder
2555
2556!------------------------------------------------------------------------------!
2557! Description:
2558! ------------
2559!> Regression relations between perceived temperature (perct) and (adjusted)
2560!> PMV. The regression presumes the Klima-Michel settings for reference
2561!> individual and reference environment.
2562!------------------------------------------------------------------------------!
[3569]2563 SUBROUTINE perct_regression( pmv, clo, perct_ij )
[3525]2564
2565    IMPLICIT NONE
2566
2567    REAL(wp), INTENT ( IN ) ::  pmv   !< Fangers predicted mean vote (dimensionless)
2568    REAL(wp), INTENT ( IN ) ::  clo   !< clothing insulation index (clo)
2569
[3569]2570    REAL(wp), INTENT ( OUT ) ::  perct_ij   !< perct (degC) corresponding to given PMV / clo
[3525]2571
[3742]2572    IF ( pmv <= -0.11_wp )  THEN
[3569]2573       perct_ij = 5.805_wp + 12.6784_wp * pmv
[3525]2574    ELSE
[3742]2575       IF ( pmv >= + 0.01_wp )  THEN
[3569]2576          perct_ij = 16.826_wp + 6.163_wp * pmv
[3525]2577       ELSE
[3569]2578          perct_ij = 21.258_wp - 9.558_wp * clo
[3525]2579       ENDIF
2580    ENDIF
2581
2582 END SUBROUTINE perct_regression
2583
2584!------------------------------------------------------------------------------!
2585! Description:
2586! ------------
2587!> FANGER.F90
2588!>
2589!> SI-VERSION: ACTLEV W m-2, DAMPFDRUCK hPa
2590!> Berechnet das aktuelle Predicted Mean Vote nach Fanger
2591!>
2592!> The case of free convection (ws < 0.1 m/s) is dealt with ws = 0.1 m/s
2593!------------------------------------------------------------------------------!
[3753]2594 SUBROUTINE fanger( ta, tmrt, pa, in_ws, pair, in_clo, actlev, eta, pmva )
[3525]2595
2596    IMPLICIT NONE
[3569]2597!
[3525]2598!-- Input variables of argument list:
2599    REAL(wp), INTENT ( IN ) ::  ta       !< Ambient air temperature (degC)
2600    REAL(wp), INTENT ( IN ) ::  tmrt     !< Mean radiant temperature (degC)
2601    REAL(wp), INTENT ( IN ) ::  pa       !< Water vapour pressure (hPa)
2602    REAL(wp), INTENT ( IN ) ::  pair     !< Barometric pressure (hPa) at site
2603    REAL(wp), INTENT ( IN ) ::  in_ws    !< Wind speed (m/s) 1 m above ground
2604    REAL(wp), INTENT ( IN ) ::  in_clo   !< Clothing insulation (clo)
2605    REAL(wp), INTENT ( IN ) ::  actlev   !< Individuals activity level per unit surface area (W/m2)
2606    REAL(wp), INTENT ( IN ) ::  eta      !< Individuals mechanical work efficiency (dimensionless)
[3569]2607!
[3525]2608!-- Output variables of argument list:
2609    REAL(wp), INTENT ( OUT ) ::  pmva    !< Actual Predicted Mean Vote (PMV,
[3693]2610                                         !< dimensionless) according to Fanger corresponding to meteorological
2611                                         !< (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta)
[3569]2612!
[3525]2613!-- Internal variables
2614    REAL(wp) ::  f_cl         !< Increase in surface due to clothing    (factor)
2615    REAL(wp) ::  heat_convection  !< energy loss by autocnvection       (W)
2616    REAL(wp) ::  activity     !< persons activity  (must stay == actlev, W)
[3593]2617    REAL(wp) ::  t_skin_aver  !< average skin temperature               (degree_C)
[3525]2618    REAL(wp) ::  bc           !< preliminary result storage
2619    REAL(wp) ::  cc           !< preliminary result storage
2620    REAL(wp) ::  dc           !< preliminary result storage
2621    REAL(wp) ::  ec           !< preliminary result storage
2622    REAL(wp) ::  gc           !< preliminary result storage
[3593]2623    REAL(wp) ::  t_clothing   !< clothing temperature                   (degree_C)
[3525]2624    REAL(wp) ::  hr           !< radiational heat resistence
2625    REAL(wp) ::  clo          !< clothing insulation index              (clo)
2626    REAL(wp) ::  ws           !< wind speed                             (m/s)
2627    REAL(wp) ::  z1           !< Empiric factor for the adaption of the heat
[3693]2628                              !< ballance equation to the psycho-physical scale (Equ. 40 in FANGER)
[3525]2629    REAL(wp) ::  z2           !< Water vapour diffution through the skin
2630    REAL(wp) ::  z3           !< Sweat evaporation from the skin surface
2631    REAL(wp) ::  z4           !< Loss of latent heat through respiration
2632    REAL(wp) ::  z5           !< Loss of radiational heat
2633    REAL(wp) ::  z6           !< Heat loss through forced convection
2634    INTEGER(iwp) :: i         !< running index
[3569]2635!
[3525]2636!-- Clo must be > 0. to avoid div. by 0!
2637    clo = in_clo
[3742]2638    IF ( clo <= 0._wp )  clo = .001_wp
[3569]2639!
[3525]2640!-- f_cl = Increase in surface due to clothing
2641    f_cl = 1._wp + .15_wp * clo
[3569]2642!
[3525]2643!-- Case of free convection (ws < 0.1 m/s ) not considered
2644    ws = in_ws
[3742]2645    IF ( ws < .1_wp )  THEN
[3525]2646       ws = .1_wp
2647    ENDIF
[3569]2648!
[3525]2649!-- Heat_convection = forced convection
[3753]2650    heat_convection = 12.1_wp * SQRT( ws * pair / 1013.25_wp )
[3569]2651!
[3525]2652!-- Activity = inner heat produktion per standardized surface
2653    activity = actlev * ( 1._wp - eta )
[3569]2654!
[3525]2655!-- T_skin_aver = average skin temperature
2656    t_skin_aver = 35.7_wp - .0275_wp * activity
[3569]2657!
[3525]2658!-- Calculation of constants for evaluation below
2659    bc = .155_wp * clo * 3.96_wp * 10._wp**( -8 ) * f_cl
2660    cc = f_cl * heat_convection
2661    ec = .155_wp * clo
2662    dc = ( 1._wp + ec * cc ) / bc
2663    gc = ( t_skin_aver + bc * ( tmrt + degc_to_k )**4 + ec * cc * ta ) / bc
[3569]2664!
[3525]2665!-- Calculation of clothing surface temperature (t_clothing) based on
[3693]2666!-- Newton-approximation with air temperature as initial guess
[3525]2667    t_clothing = ta
[3742]2668    DO  i = 1, 3
[3525]2669       t_clothing = t_clothing - ( ( t_clothing + degc_to_k )**4 + t_clothing  &
2670          * dc - gc ) / ( 4._wp * ( t_clothing + degc_to_k )**3 + dc )
2671    ENDDO
[3569]2672!
[3525]2673!-- Empiric factor for the adaption of the heat ballance equation
[3693]2674!-- to the psycho-physical scale (Equ. 40 in FANGER)
[3753]2675    z1 = ( .303_wp * EXP( -.036_wp * actlev ) + .0275_wp )
[3569]2676!
[3525]2677!-- Water vapour diffution through the skin
2678    z2 = .31_wp * ( 57.3_wp - .07_wp * activity-pa )
[3569]2679!
[3525]2680!-- Sweat evaporation from the skin surface
2681    z3 = .42_wp * ( activity - 58._wp )
[3569]2682!
[3525]2683!-- Loss of latent heat through respiration
2684    z4 = .0017_wp * actlev * ( 58.7_wp - pa ) + .0014_wp * actlev *            &
2685      ( 34._wp - ta )
[3569]2686!
[3525]2687!-- Loss of radiational heat
2688    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + degc_to_k )**4 - ( tmrt +        &
2689       degc_to_k )**4 )
[3753]2690    IF ( ABS( t_clothing - tmrt ) > 0._wp )  THEN
[3525]2691       hr = z5 / f_cl / ( t_clothing - tmrt )
2692    ELSE
2693       hr = 0._wp
2694    ENDIF
[3569]2695!
[3525]2696!-- Heat loss through forced convection cc*(t_clothing-TT)
2697    z6 = cc * ( t_clothing - ta )
[3569]2698!
[3525]2699!-- Predicted Mean Vote
2700    pmva = z1 * ( activity - z2 - z3 - z4 - z5 - z6 )
2701
2702 END SUBROUTINE fanger
2703
2704!------------------------------------------------------------------------------!
2705! Description:
2706! ------------
2707!> For pmva > 0 and clo =0.5 the increment (deltapmv) is calculated
2708!> that converts pmva into Gagge's et al. (1986) PMV*.
2709!------------------------------------------------------------------------------!
2710 REAL(wp) FUNCTION deltapmv( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
2711
2712    IMPLICIT NONE
[3693]2713
[3569]2714!
[3525]2715!-- Input variables of argument list:
2716    REAL(wp),     INTENT ( IN )  :: pmva     !< Actual Predicted Mean Vote (PMV) according to Fanger
2717    REAL(wp),     INTENT ( IN )  :: ta       !< Ambient temperature (degC) at screen level
2718    REAL(wp),     INTENT ( IN )  :: vp       !< Water vapour pressure (hPa) at screen level
2719    REAL(wp),     INTENT ( IN )  :: svp_ta   !< Saturation water vapour pressure (hPa) at ta
2720    REAL(wp),     INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC) at screen level
2721    REAL(wp),     INTENT ( IN )  :: ws       !< Wind speed (m/s) 1 m above ground
[3693]2722
[3569]2723!
[3525]2724!-- Output variables of argument list:
2725    INTEGER(iwp), INTENT ( OUT ) :: nerr     !< Error status / quality flag
[3693]2726                                             !<  0 = o.k.
2727                                             !< -2 = pmva outside valid regression range
2728                                             !< -3 = rel. humidity set to 5 % or 95 %, respectively
2729                                             !< -4 = deltapmv set to avoid pmvs < 0
2730
[3569]2731!
[3693]2732!-- Internal variables:
2733    REAL(wp) ::  pmv          !< temp storage og predicted mean vote
2734    REAL(wp) ::  pa_p50       !< ratio actual water vapour pressure to that of relative humidity of 50 %
2735    REAL(wp) ::  pa           !< vapor pressure (hPa) with hard bounds
2736    REAL(wp) ::  apa          !< natural logarithm of pa (with hard lower border)
2737    REAL(wp) ::  dapa         !< difference of apa and pa_p50
2738    REAL(wp) ::  sqvel        !< square root of local wind velocity
2739    REAL(wp) ::  dtmrt        !< difference mean radiation to air temperature
2740    REAL(wp) ::  p10          !< lower bound for pa
2741    REAL(wp) ::  p95          !< upper bound for pa
2742    REAL(wp) ::  weight       !<
2743    REAL(wp) ::  weight2      !<
[3525]2744    REAL(wp) ::  dpmv_1       !<
2745    REAL(wp) ::  dpmv_2       !<
2746    REAL(wp) ::  pmvs         !<
2747    INTEGER(iwp) :: nreg      !<
2748
[3569]2749!
[3753]2750!-- Regression coefficients:
2751    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpmv = (/                          &
2752     -0.0556602_wp, -0.1528680_wp, -0.2336104_wp, -0.2789387_wp,               &
2753     -0.3551048_wp, -0.4304076_wp, -0.4884961_wp, -0.4897495_wp /)
2754
2755    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpa_p50 = (/                       &
2756     -0.1607154_wp, -0.4177296_wp, -0.4120541_wp, -0.0886564_wp,               &
2757      0.4285938_wp,  0.6281256_wp,  0.5067361_wp,  0.3965169_wp /)
2758
2759    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpa = (/                           &
2760      0.0580284_wp,  0.0836264_wp,  0.1009919_wp,  0.1020777_wp,               &
2761      0.0898681_wp,  0.0839116_wp,  0.0853258_wp,  0.0866589_wp /)
2762
2763    REAL(wp), DIMENSION(0:7), PARAMETER ::  bapa = (/                          &
2764     -1.7838788_wp, -2.9306231_wp, -1.6350334_wp,   0.6211547_wp,              &
2765      3.3918083_wp,  5.5521025_wp,  8.4897418_wp,  16.6265851_wp /)
2766
2767    REAL(wp), DIMENSION(0:7), PARAMETER ::  bdapa = (/                         &
2768      1.6752720_wp,  2.7379504_wp,  1.2940526_wp,  -1.0985759_wp,              &
2769     -3.9054732_wp, -6.0403012_wp, -8.9437119_wp, -17.0671201_wp /)
2770
2771    REAL(wp), DIMENSION(0:7), PARAMETER ::  bsqvel = (/                        &
2772     -0.0315598_wp, -0.0286272_wp, -0.0009228_wp,  0.0483344_wp,               &
2773      0.0992366_wp,  0.1491379_wp,  0.1951452_wp,  0.2133949_wp /)
2774
2775    REAL(wp), DIMENSION(0:7), PARAMETER ::  bta = (/                           &
2776      0.0953986_wp,  0.1524760_wp,  0.0564241_wp, -0.0893253_wp,               &
2777     -0.2398868_wp, -0.3515237_wp, -0.5095144_wp, -0.9469258_wp /)
2778
2779    REAL(wp), DIMENSION(0:7), PARAMETER ::  bdtmrt = (/                        &
2780     -0.0004672_wp, -0.0000514_wp, -0.0018037_wp, -0.0049440_wp,               &
2781     -0.0069036_wp, -0.0075844_wp, -0.0079602_wp, -0.0089439_wp /)
2782
2783    REAL(wp), DIMENSION(0:7), PARAMETER ::  aconst = (/                        &
2784      1.8686215_wp,  3.4260713_wp,   2.0116185_wp,  -0.7777552_wp,             &
2785     -4.6715853_wp, -7.7314281_wp, -11.7602578_wp, -23.5934198_wp /)
2786
2787
2788!
[3525]2789!-- Test for compliance with regression range
[3742]2790    IF ( pmva < -1.0_wp  .OR.  pmva > 7.0_wp )  THEN
[3525]2791       nerr = -2_iwp
2792    ELSE
2793       nerr = 0_iwp
2794    ENDIF
[3569]2795!
[3525]2796!-- Initialise classic PMV
2797    pmv  = pmva
[3569]2798!
[3525]2799!-- Water vapour pressure of air
2800    p10  = 0.05_wp * svp_ta
2801    p95  = 1.00_wp * svp_ta
[3742]2802    IF ( vp >= p10  .AND.  vp <= p95 )  THEN
[3525]2803       pa = vp
2804    ELSE
2805       nerr = -3_iwp
[3742]2806       IF ( vp < p10 )  THEN
[3569]2807!
[3525]2808!--       Due to conditions of regression: r.H. >= 5 %
2809          pa = p10
2810       ELSE
[3569]2811!
[3525]2812!--       Due to conditions of regression: r.H. <= 95 %
2813          pa = p95
2814       ENDIF
2815    ENDIF
[3742]2816    IF ( pa > 0._wp )  THEN
[3569]2817!
[3525]2818!--    Natural logarithm of pa
[3753]2819       apa = LOG( pa )
[3525]2820    ELSE
2821       apa = -5._wp
2822    ENDIF
[3569]2823!
[3525]2824!-- Ratio actual water vapour pressure to that of a r.H. of 50 %
2825    pa_p50   = 0.5_wp * svp_ta
[3742]2826    IF ( pa_p50 > 0._wp  .AND.  pa > 0._wp )  THEN
[3753]2827       dapa   = apa - LOG( pa_p50 )
[3525]2828       pa_p50 = pa / pa_p50
2829    ELSE
2830       dapa   = -5._wp
2831       pa_p50 = 0._wp
2832    ENDIF
[3569]2833!
[3525]2834!-- Square root of wind velocity
[3742]2835    IF ( ws >= 0._wp )  THEN
[3753]2836       sqvel = SQRT( ws )
[3525]2837    ELSE
2838       sqvel = 0._wp
2839    ENDIF
[3569]2840!
[3525]2841!-- Difference mean radiation to air temperature
2842    dtmrt = tmrt - ta
[3569]2843!
[3525]2844!-- Select the valid regression coefficients
[3753]2845    nreg = INT( pmv )
[3742]2846    IF ( nreg < 0_iwp )  THEN
[3569]2847!
[3525]2848!--    value of the FUNCTION in the case pmv <= -1
2849       deltapmv = 0._wp
2850       RETURN
2851    ENDIF
[3693]2852    weight = MOD ( pmv, 1._wp )
[3742]2853    IF ( weight < 0._wp )  weight = 0._wp
2854    IF ( nreg > 5_iwp )  THEN
[3525]2855       nreg  = 5_iwp
[3693]2856       weight   = pmv - 5._wp
2857       weight2  = pmv - 6._wp
[3742]2858       IF ( weight2 > 0_iwp )  THEN
[3693]2859          weight = ( weight - weight2 ) / weight
[3525]2860       ENDIF
2861    ENDIF
[3569]2862!
[3753]2863!-- Regression valid for 0. <= pmv <= 6., bounds are checked above
[3525]2864    dpmv_1 =                                                                   &
[3742]2865       + bpa(nreg) * pa                                                        &
2866       + bpmv(nreg) * pmv                                                      &
2867       + bapa(nreg) * apa                                                      &
2868       + bta(nreg) * ta                                                        &
2869       + bdtmrt(nreg) * dtmrt                                                  &
2870       + bdapa(nreg) * dapa                                                    &
2871       + bsqvel(nreg) * sqvel                                                  &
2872       + bpa_p50(nreg) * pa_p50                                                &
2873       + aconst(nreg)
[3525]2874
[3753]2875!    dpmv_2 = 0._wp
2876!    IF ( nreg < 6_iwp )  THEN  !< nreg is always <= 5, see above
2877    dpmv_2 =                                                                   &
2878       + bpa(nreg+1_iwp)     * pa                                              &
2879       + bpmv(nreg+1_iwp)    * pmv                                             &
2880       + bapa(nreg+1_iwp)    * apa                                             &
2881       + bta(nreg+1_iwp)     * ta                                              &
2882       + bdtmrt(nreg+1_iwp)  * dtmrt                                           &
2883       + bdapa(nreg+1_iwp)   * dapa                                            &
2884       + bsqvel(nreg+1_iwp)  * sqvel                                           &
2885       + bpa_p50(nreg+1_iwp) * pa_p50                                          &
2886       + aconst(nreg+1_iwp)
2887!    ENDIF
[3569]2888!
[3525]2889!-- Calculate pmv modification
[3693]2890    deltapmv = ( 1._wp - weight ) * dpmv_1 + weight * dpmv_2
[3525]2891    pmvs = pmva + deltapmv
[3742]2892    IF ( ( pmvs ) < 0._wp )  THEN
[3569]2893!
[3525]2894!--    Prevent negative pmv* due to problems with clothing insulation
2895       nerr = -4_iwp
[3742]2896       IF ( pmvs > -0.11_wp )  THEN
[3569]2897!
[3753]2898!--       Threshold from perct_regression for winter clothing insulation
[3525]2899          deltapmv = deltapmv + 0.11_wp
2900       ELSE
[3569]2901!
[3525]2902!--       Set pmvs to "0" for compliance with summer clothing insulation
2903          deltapmv = -1._wp * pmva
2904       ENDIF
2905    ENDIF
2906
2907 END FUNCTION deltapmv
2908
2909!------------------------------------------------------------------------------!
2910! Description:
2911! ------------
2912!> The subroutine "calc_sultr" returns a threshold value to perceived
2913!> temperature allowing to decide whether the actual perceived temperature
2914!> is linked to perecption of sultriness. The threshold values depends
2915!> on the Fanger's classical PMV, expressed here as perceived temperature
2916!> perct.
2917!------------------------------------------------------------------------------!
[3569]2918 SUBROUTINE calc_sultr( perct_ij, dperctm, dperctstd, sultr_res )
[3525]2919
2920    IMPLICIT NONE
[3569]2921!
[3525]2922!-- Input of the argument list:
[3569]2923    REAL(wp), INTENT ( IN )  ::  perct_ij      !< Classical perceived temperature: Base is Fanger's PMV
2924!
[3525]2925!-- Additional output variables of argument list:
2926    REAL(wp), INTENT ( OUT ) ::  dperctm    !< Mean deviation perct (classical gt) to gt* (rational gt
[3693]2927                                            !< calculated based on Gagge's rational PMV*)
[3525]2928    REAL(wp), INTENT ( OUT ) ::  dperctstd  !< dperctm plus its standard deviation times a factor
[3693]2929                                            !< determining the significance to perceive sultriness
[3525]2930    REAL(wp), INTENT ( OUT ) ::  sultr_res
[3569]2931!
[3525]2932!-- Types of coefficients mean deviation: third order polynomial
[3753]2933    REAL(wp), PARAMETER ::  dperctka =  7.5776086_wp
[3525]2934    REAL(wp), PARAMETER ::  dperctkb = -0.740603_wp
[3753]2935    REAL(wp), PARAMETER ::  dperctkc =  0.0213324_wp
[3525]2936    REAL(wp), PARAMETER ::  dperctkd = -0.00027797237_wp
[3569]2937!
[3525]2938!-- Types of coefficients mean deviation plus standard deviation
[3693]2939!-- regression coefficients: third order polynomial
[3753]2940    REAL(wp), PARAMETER ::  dperctsa =  0.0268918_wp
2941    REAL(wp), PARAMETER ::  dperctsb =  0.0465957_wp
[3525]2942    REAL(wp), PARAMETER ::  dperctsc = -0.00054709752_wp
[3753]2943    REAL(wp), PARAMETER ::  dperctsd =  0.0000063714823_wp
[3569]2944!
[3525]2945!-- Factor to mean standard deviation defining SIGNificance for
[3693]2946!-- sultriness
[3525]2947    REAL(wp), PARAMETER :: faktor = 1._wp
[3569]2948!
[3525]2949!-- Initialise
[3753]2950    sultr_res = 99._wp