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

Last change on this file since 3750 was 3750, checked in by dom_dwd_user, 2 years ago

biometeorology_mod.f90:
(C) Fixed several ugly/bad style issues. Now compiles with -Wall and -pedantic without errors or notices.

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