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

Last change on this file since 4797 was 4768, checked in by suehring, 4 years ago

Enable 3D data output also with 64-bit precision

  • 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-4492
    /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: 193.6 KB
Line 
1!> @file biometeorology_mod.f90
2!--------------------------------------------------------------------------------------------------!
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 terms of the GNU
6! General Public License as published by the Free Software Foundation, either version 3 of the
7! License, or (at your option) any later version.
8!
9! PALM-4U is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even
10! the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 2018-2020 Deutscher Wetterdienst (DWD)
17! Copyright 2018-2020 Institute of Computer Science, Academy of Sciences, Prague
18! Copyright 2018-2020 Leibniz Universitaet Hannover
19!--------------------------------------------------------------------------------------------------!
20!
21! Current revisions:
22! ------------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: biometeorology_mod.f90 4768 2020-11-02 19:11:23Z raasch $
28! Enable 3D data output also with 64-bit precision
29!
30! 4633 2020-08-05 14:21:14Z suehring
31! Bugfix in check for humidity
32!
33! 4590 2020-07-06 14:34:59Z suehring
34! Enable restart via mpi-IO. Therefore, allocated array mrt_av_grid as 3D array instead of an
35! 1D array.
36!
37! 4577 2020-06-25 09:53:58Z raasch
38! further re-formatting concerning Fortran parameter variables
39!
40! 4540 2020-05-18 15:23:29Z raasch
41! file re-formatted to follow the PALM coding standard
42!
43! 4535 2020-05-15 12:07:23Z raasch
44! bugfix for restart data format query
45!
46! 4517 2020-05-03 14:29:30Z raasch
47! added restart with MPI-IO for reading local arrays
48!
49! 4495 2020-04-13 20:11:20Z raasch
50! restart data handling with MPI-IO added
51!
52! 4493 2020-04-10 09:49:43Z pavelkrc
53! Revise bad formatting
54!
55! 4286 2019-10-30 16:01:14Z resler
56! implement new palm_date_time_mod
57!
58! 4223 2019-09-10 09:20:47Z gronemeier
59! Corrected "Former revisions" section
60!
61! 4168 2019-08-16 13:50:17Z suehring
62! Replace function get_topography_top_index by topo_top_ind
63!
64! 4144 2019-08-06 09:11:47Z raasch
65! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
66!
67! 4127 2019-07-30 14:47:10Z suehring
68! Output for bio_mrt added (merge from branch resler)
69!
70! 4126 2019-07-30 11:09:11Z gronemeier
71! renamed vitd3_exposure_av into vitd3_dose,
72! renamed uvem_calc_exposure into bio_calculate_uv_exposure
73!
74! 3885 2019-04-11 11:29:34Z kanani
75! Changes related to global restructuring of location messages and introduction of additional debug
76! messages
77!
78! 3753 2019-02-19 14:48:54Z dom_dwd_user
79! - Added automatic setting of mrt_nlevels in case it was not part of radiation_parameters namelist
80!   (or set to 0 accidentially).
81! - Minor speed improvoemnts in perceived temperature calculations.
82! - Perceived temperature regression arrays now declared as PARAMETERs.
83!
84! 3750 2019-02-19 07:29:39Z dom_dwd_user
85! - Added addittional safety meassures to bio_calculate_thermal_index_maps.
86! - Replaced several REAL (un-)equality comparisons.
87!
88! 3742 2019-02-14 11:25:22Z dom_dwd_user
89! - Allocation of the input _av grids was moved to the "sum" section of bio_3d_data_averaging to
90!   make sure averaging is only done once!
91! - Moved call of bio_calculate_thermal_index_maps from biometeorology module to time_integration to
92!   make sure averaged input is updated before calculating.
93!
94! 3740 2019-02-13 12:35:12Z dom_dwd_user
95! - Added safety-meassure to catch the case that 'bio_mrt_av' is stated after 'bio_<index>' in the
96!   output section of the p3d file.
97!
98! 3739 2019-02-13 08:05:17Z dom_dwd_user
99! - Auto-adjusting thermal_comfort flag if not set by user, but thermal_indices set as output
100!   quantities.
101! - Renamed flags "bio_<index>" to "do_calculate_<index>" for better readability
102! - Removed everything related to "time_bio_results" as this is never used.
103! - Moved humidity warning to check_data_output.
104! - Fixed bug in mrt calculation introduced with my commit yesterday.
105!
106! 3735 2019-02-12 09:52:40Z dom_dwd_user
107! - Fixed auto-setting of thermal index calculation flags by output as originally proposed by
108!   resler.
109! - removed bio_pet and outher configuration variables.
110! - Updated namelist.
111!
112! 3711 2019-01-31 13:44:26Z knoop
113! Introduced interface routine bio_init_checks + small error message changes
114!
115! 3693 2019-01-23 15:20:53Z dom_dwd_user
116! Added usage of time_averaged mean radiant temperature, together with calculation, grid and restart
117! routines. General cleanup and commenting.
118!
119! 3685 2019-01-21 01:02:11Z knoop
120! Some interface calls moved to module_interface + cleanup
121!
122! 3650 2019-01-04 13:01:33Z kanani
123! Bugfixes and additions for enabling restarts with biometeorology
124!
125! 3448 2018-10-29 18:14:31Z kanani
126! Initial revision
127!
128!
129!
130! Authors:
131! --------
132! @author Dominik Froehlich <dominik.froehlich@dwd.de>, thermal indices
133! @author Jaroslav Resler <resler@cs.cas.cz>, mean radiant temperature
134! @author Michael Schrempf <schrempf@muk.uni-hannover.de>, uv exposure
135!
136!
137! Description:
138! ------------
139!> Biometeorology module consisting of two parts:
140!> 1.: Human thermal comfort module calculating thermal perception of a sample human being under the
141!> current meteorological conditions.
142!> 2.: Calculation of vitamin-D weighted UV exposure
143!>
144!> @todo Alphabetical sorting of "USE ..." lists, "ONLY" list, variable declarations
145!>       (per subroutine: first all CHARACTERs, then INTEGERs, LOGICALs, REALs, )
146!> @todo Comments start with capital letter --> "!-- Include..."
147!> @todo uv_vitd3dose-->new output type necessary (cumulative)
148!> @todo consider upwelling radiation in UV
149!>
150!> @note nothing now
151!>
152!> @bug  no known bugs by now
153!--------------------------------------------------------------------------------------------------!
154 MODULE biometeorology_mod
155
156    USE arrays_3d,                                                                                 &
157        ONLY:  pt, p, u, v, w, q
158
159    USE averaging,                                                                                 &
160        ONLY:  pt_av, q_av, u_av, v_av, w_av
161
162    USE basic_constants_and_equations_mod,                                                         &
163        ONLY: degc_to_k, c_p, l_v, magnus, pi, sigma_sb
164
165    USE control_parameters,                                                                        &
166        ONLY:  average_count_3d, biometeorology,                                                   &
167               debug_output,                                                                       &
168               dz, dz_stretch_factor,                                                              &
169               dz_stretch_level, humidity, initializing_actions, nz_do3d,                          &
170               restart_data_format_output, surface_pressure
171
172    USE grid_variables,                                                                            &
173        ONLY:  ddx, dx, ddy, dy
174
175    USE indices,                                                                                   &
176        ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg, nysg, nyng,           &
177               topo_top_ind
178
179    USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM
180
181    USE netcdf_data_input_mod,                                                                     &
182        ONLY: building_obstruction_f, netcdf_data_input_uvem, uvem_integration_f,                  &
183              uvem_irradiance_f, uvem_projarea_f, uvem_radiance_f
184
185    USE palm_date_time_mod,                                                                        &
186        ONLY:  get_date_time
187!
188!-- Import radiation model to obtain input for mean radiant temperature
189    USE radiation_model_mod,                                                                       &
190        ONLY:  id, ix, iy, iz, mrt_include_sw, mrt_nlevels,                                        &
191               mrtbl, mrtinlw, mrtinsw, nmrtbl, radiation,                                         &
192               rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out, radiation_interactions
193
194    USE restart_data_mpi_io_mod,                                                                   &
195        ONLY:  rrd_mpi_io,                                                                         &
196               rd_mpi_io_check_array,                                                              &
197               wrd_mpi_io
198
199
200    IMPLICIT NONE
201
202!
203!-- Declare all global variables within the module (alphabetical order)
204    REAL(wp), PARAMETER ::  bio_fill_value = -9999.0_wp  !< set module fill value, replace by global fill value as soon as available
205    REAL(wp), PARAMETER ::  human_absorb = 0.7_wp  !< SW absorbtivity of a human body (Fanger 1972)
206    REAL(wp), PARAMETER ::  human_emiss = 0.97_wp  !< LW emissivity of a human body after (Fanger 1972)
207
208    INTEGER(iwp) ::  bio_cell_level     !< cell level biom calculates for
209
210    LOGICAL ::  thermal_comfort        = .FALSE.  !< Enables or disables the entire thermal comfort part
211    LOGICAL ::  do_average_theta       = .FALSE.  !< switch: do theta averaging in this module? (if .FALSE. this is done globally)
212    LOGICAL ::  do_average_q           = .FALSE.  !< switch: do e averaging in this module?
213    LOGICAL ::  do_average_u           = .FALSE.  !< switch: do u averaging in this module?
214    LOGICAL ::  do_average_v           = .FALSE.  !< switch: do v averaging in this module?
215    LOGICAL ::  do_average_w           = .FALSE.  !< switch: do w averaging in this module?
216    LOGICAL ::  do_average_mrt         = .FALSE.  !< switch: do mrt averaging in this module?
217    LOGICAL ::  average_trigger_perct  = .FALSE.  !< update averaged input on call to bio_perct?
218    LOGICAL ::  average_trigger_utci   = .FALSE.  !< update averaged input on call to bio_utci?
219    LOGICAL ::  average_trigger_pet    = .FALSE.  !< update averaged input on call to bio_pet?
220    LOGICAL ::  average_trigger_mrt    = .FALSE.  !< update averaged input on call to bio_pet?
221    LOGICAL ::  do_calculate_perct     = .FALSE.  !< Turn index PT (instant. input) on or off
222    LOGICAL ::  do_calculate_perct_av  = .FALSE.  !< Turn index PT (averaged input) on or off
223    LOGICAL ::  do_calculate_pet       = .FALSE.  !< Turn index PET (instant. input) on or off
224    LOGICAL ::  do_calculate_pet_av    = .FALSE.  !< Turn index PET (averaged input) on or off
225    LOGICAL ::  do_calculate_utci      = .FALSE.  !< Turn index UTCI (instant. input) on or off
226    LOGICAL ::  do_calculate_utci_av   = .FALSE.  !< Turn index UTCI (averaged input) on or off
227    LOGICAL ::  do_calculate_mrt2d     = .FALSE.  !< Turn index MRT 2D (averaged or inst) on or off
228
229    REAL(wp)    ::  bio_output_height  !< height output is calculated in m
230
231    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct      !< PT results   (degree_C)
232    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet        !< PET results  (degree_C)
233    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_grid  !< tmrt results (degree_C)
234    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci       !< UTCI results (degree_C)
235!
236!-- Grids for averaged thermal indices
237    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  perct_av      !< PT results (aver. input)   (degree_C)
238    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pet_av        !< PET results (aver. input)  (degree_C)
239    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  tmrt_av_grid  !< tmrt results (degree_C)
240    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  utci_av       !< UTCI results (aver. input) (degree_C)
241
242    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  mrt_av_grid   !< time average mean radiant temperature
243
244!
245!-- UVEM parameters from here
246    INTEGER(iwp) ::  ai                      = 0  !< loop index in azimuth direction
247    INTEGER(iwp) ::  bi                      = 0  !< loop index of bit location within an 8bit-integer (one Byte)
248    INTEGER(iwp) ::  bio_nmrtbl
249    INTEGER(iwp) ::  clothing                = 1  !< clothing (0=unclothed, 1=Arms,Hands,Face free, 3=Hand,Face free)
250    INTEGER(iwp) ::  iq                      = 0  !< loop index of irradiance quantity
251    INTEGER(iwp) ::  pobi                    = 0  !< loop index of the position of corresponding byte within ibset byte vektor
252    INTEGER(iwp) ::  obstruction_direct_beam = 0  !< Obstruction information for direct beam
253    INTEGER(iwp) ::  zi                      = 0  !< loop index in zenith direction
254
255    INTEGER(KIND=1), DIMENSION(0:44)  ::  obstruction_temp1 = 0  !< temporary obstruction information stored with ibset
256    INTEGER(iwp),    DIMENSION(0:359) ::  obstruction_temp2 = 0  !< restored temporary obstruction information from ibset file
257
258    INTEGER(iwp), DIMENSION(0:35,0:9) ::  obstruction       = 1  !< final 2D obstruction information array
259
260    LOGICAL ::  consider_obstructions = .TRUE.   !< namelist parameter (see documentation)
261    LOGICAL ::  sun_in_south          = .FALSE.  !< namelist parameter (see documentation)
262    LOGICAL ::  turn_to_sun           = .TRUE.   !< namelist parameter (see documentation)
263    LOGICAL ::  uv_exposure           = .FALSE.  !< namelist parameter (see documentation)
264
265    REAL(wp) ::  diffuse_exposure            =   0.0_wp  !< calculated exposure by diffuse radiation
266    REAL(wp) ::  direct_exposure             =   0.0_wp  !< calculated exposure by direct solar beam
267    REAL(wp) ::  orientation_angle           =   0.0_wp  !< orientation of front/face of the human model
268    REAL(wp) ::  projection_area_direct_beam =   0.0_wp  !< projection area for direct solar beam
269    REAL(wp) ::  saa                         = 180.0_wp  !< solar azimuth angle
270    REAL(wp) ::  startpos_human              =   0.0_wp  !< start value for azimuth interpolation of human geometry array
271    REAL(wp) ::  startpos_saa_float          =   0.0_wp  !< start value for azimuth interpolation of radiance array
272    REAL(wp) ::  sza                         =  20.0_wp  !< solar zenith angle
273    REAL(wp) ::  xfactor                     =   0.0_wp  !< relative x-position used for interpolation
274    REAL(wp) ::  yfactor                     =   0.0_wp  !< relative y-position used for interpolation
275
276    REAL(wp), DIMENSION(0:2)  ::  irradiance =   0.0_wp  !< iradiance values extracted from irradiance lookup table
277
278    REAL(wp), DIMENSION(0:2,0:90) ::  irradiance_lookup_table      = 0.0_wp  !< irradiance lookup table
279    REAL(wp), DIMENSION(0:35,0:9) ::  integration_array            = 0.0_wp  !< solid angle factors for hemispherical integration
280    REAL(wp), DIMENSION(0:35,0:9) ::  projection_area              = 0.0_wp  !< projection areas of a human (all directions)
281    REAL(wp), DIMENSION(0:35,0:9) ::  projection_area_lookup_table = 0.0_wp  !< human geometry lookup table (projection areas)
282    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_direct_temp  = 0.0_wp  !< temporary projection area for direct solar beam
283    REAL(wp), DIMENSION(0:71,0:9) ::  projection_area_temp         = 0.0_wp  !< temporary projection area for all directions
284    REAL(wp), DIMENSION(0:35,0:9) ::  radiance_array               = 0.0_wp  !< radiance extracted from radiance_lookup_table
285    REAL(wp), DIMENSION(0:71,0:9) ::  radiance_array_temp          = 0.0_wp  !< temporary radiance data
286
287    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vitd3_exposure  !< result variable for instantaneous vitamin-D weighted exposures
288    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  vitd3_dose      !< result variable for summation of vitamin-D weighted exposures
289
290    REAL(wp), DIMENSION(0:35,0:9,0:90) ::  radiance_lookup_table   = 0.0_wp  !< radiance lookup table
291
292
293    PRIVATE
294
295!
296!-- INTERFACES that must be available to other modules (alphabetical order)
297    PUBLIC bio_3d_data_averaging,  bio_calculate_mrt_grid, bio_calculate_thermal_index_maps,       &
298           bio_calc_ipt, bio_check_data_output, bio_check_parameters,                              &
299           bio_data_output_2d, bio_data_output_3d,  bio_define_netcdf_grid,                        &
300           bio_get_thermal_index_input_ij, bio_header, bio_init, bio_init_checks, bio_nmrtbl,      &
301           bio_parin, bio_rrd_global, bio_rrd_local, bio_wrd_global, bio_wrd_local, thermal_comfort
302!
303!-- UVEM PUBLIC variables and methods
304    PUBLIC bio_calculate_uv_exposure, uv_exposure
305
306!
307!-- PALM interfaces:
308!
309!-- 3D averaging for HTCM _INPUT_ variables
310    INTERFACE bio_3d_data_averaging
311       MODULE PROCEDURE bio_3d_data_averaging
312    END INTERFACE bio_3d_data_averaging
313!
314!-- Calculate mtr from rtm fluxes and assign into 2D grid
315    INTERFACE bio_calculate_mrt_grid
316       MODULE PROCEDURE bio_calculate_mrt_grid
317    END INTERFACE bio_calculate_mrt_grid
318!
319!-- Calculate static thermal indices PT, UTCI and/or PET
320    INTERFACE bio_calculate_thermal_index_maps
321       MODULE PROCEDURE bio_calculate_thermal_index_maps
322    END INTERFACE bio_calculate_thermal_index_maps
323!
324!-- Calculate the dynamic index iPT (to be caled by the agent model)
325    INTERFACE bio_calc_ipt
326       MODULE PROCEDURE bio_calc_ipt
327    END INTERFACE bio_calc_ipt
328!
329!-- Data output checks for 2D/3D data to be done in check_parameters
330    INTERFACE bio_check_data_output
331       MODULE PROCEDURE bio_check_data_output
332    END INTERFACE bio_check_data_output
333!
334!-- Input parameter checks to be done in check_parameters
335    INTERFACE bio_check_parameters
336       MODULE PROCEDURE bio_check_parameters
337    END INTERFACE bio_check_parameters
338!
339!-- Data output of 2D quantities
340    INTERFACE bio_data_output_2d
341       MODULE PROCEDURE bio_data_output_2d
342    END INTERFACE bio_data_output_2d
343!
344!-- no 3D data, thus, no averaging of 3D data, removed
345    INTERFACE bio_data_output_3d
346       MODULE PROCEDURE bio_data_output_3d
347    END INTERFACE bio_data_output_3d
348!
349!-- Definition of data output quantities
350    INTERFACE bio_define_netcdf_grid
351       MODULE PROCEDURE bio_define_netcdf_grid
352    END INTERFACE bio_define_netcdf_grid
353!
354!-- Obtains all relevant input values to estimate local thermal comfort/stress
355    INTERFACE bio_get_thermal_index_input_ij
356       MODULE PROCEDURE bio_get_thermal_index_input_ij
357    END INTERFACE bio_get_thermal_index_input_ij
358!
359!-- Output of information to the header file
360    INTERFACE bio_header
361       MODULE PROCEDURE bio_header
362    END INTERFACE bio_header
363!
364!-- Initialization actions
365    INTERFACE bio_init
366       MODULE PROCEDURE bio_init
367    END INTERFACE bio_init
368!
369!-- Initialization checks
370    INTERFACE bio_init_checks
371       MODULE PROCEDURE bio_init_checks
372    END INTERFACE bio_init_checks
373!
374!-- Reading of NAMELIST parameters
375    INTERFACE bio_parin
376       MODULE PROCEDURE bio_parin
377    END INTERFACE bio_parin
378!
379!-- Read global restart parameters
380    INTERFACE bio_rrd_global
381       MODULE PROCEDURE bio_rrd_global_ftn
382       MODULE PROCEDURE bio_rrd_global_mpi
383    END INTERFACE bio_rrd_global
384!
385!-- Read local restart parameters
386    INTERFACE bio_rrd_local
387       MODULE PROCEDURE bio_rrd_local_ftn
388       MODULE PROCEDURE bio_rrd_local_mpi
389    END INTERFACE bio_rrd_local
390!
391!-- Write global restart parameters
392    INTERFACE bio_wrd_global
393       MODULE PROCEDURE bio_wrd_global
394    END INTERFACE bio_wrd_global
395!
396!-- Write local restart parameters
397    INTERFACE bio_wrd_local
398       MODULE PROCEDURE bio_wrd_local
399    END INTERFACE bio_wrd_local
400!
401!-- Calculate UV exposure grid
402    INTERFACE bio_calculate_uv_exposure
403       MODULE PROCEDURE bio_calculate_uv_exposure
404    END INTERFACE bio_calculate_uv_exposure
405
406 CONTAINS
407
408
409!--------------------------------------------------------------------------------------------------!
410! Description:
411! ------------
412!> Sum up and time-average biom input quantities as well as allocate the array necessary for storing
413!> the average.
414!> There is a considerable difference to the 3d_data_averaging subroutines used by other modules:
415!> For the thermal indices, the module needs to average the input conditions, not the result!
416!--------------------------------------------------------------------------------------------------!
417 SUBROUTINE bio_3d_data_averaging( mode, variable )
418
419    IMPLICIT NONE
420
421    CHARACTER (LEN=*) ::  mode     !< Averaging mode: allocate, sum, or average
422    CHARACTER (LEN=*) ::  variable !< The variable in question
423
424    INTEGER(iwp) ::  i        !< Running index, x-direction
425    INTEGER(iwp) ::  j        !< Running index, y-direction
426    INTEGER(iwp) ::  k        !< Running index, z-direction
427    INTEGER(iwp) ::  l        !< index used to link radiation arrays to 3d grid arrays
428
429
430    IF ( mode == 'allocate' )  THEN
431
432       SELECT CASE ( TRIM( variable ) )
433
434          CASE ( 'bio_mrt' )
435
436                IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
437                   ALLOCATE( mrt_av_grid(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
438                ENDIF
439                mrt_av_grid = 0.0_wp
440                do_average_mrt = .FALSE.  !< Overwrite if that was enabled somehow
441
442
443          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*', 'bio_mrt*' )
444
445!
446!--          Averaging, as well as the allocation of the required grids must be done only once,
447!--          independent of for how many thermal indices averaged output is desired.
448!--          Therefore we need to memorize which index is the one that controls the averaging
449!--          (what must be the first thermal index called).
450!--          Indices are in unknown order as depending on the input file, determine first index to
451!--          average und update only once.
452!
453!--          Only proceed here if this was not done for any index before. This is done only once
454!--          during the whole model run.
455             IF ( .NOT. average_trigger_perct  .AND.                                               &
456                  .NOT. average_trigger_utci   .AND.                                               &
457                  .NOT. average_trigger_pet    .AND.                                               &
458                  .NOT. average_trigger_mrt )  THEN
459!
460!--             Memorize the first index called to control averaging
461                IF ( TRIM( variable ) == 'bio_perct*' )  THEN
462                    average_trigger_perct = .TRUE.
463                ENDIF
464                IF ( TRIM( variable ) == 'bio_utci*' )  THEN
465                    average_trigger_utci = .TRUE.
466                ENDIF
467                IF ( TRIM( variable ) == 'bio_pet*' )  THEN
468                    average_trigger_pet = .TRUE.
469                ENDIF
470                IF ( TRIM( variable ) == 'bio_mrt*' )  THEN
471                    average_trigger_mrt = .TRUE.
472                ENDIF
473             ENDIF
474!
475!--          Allocation of the input _av grids was moved to the "sum" section to make sure averaging
476!--          is only done once!
477
478          CASE ( 'uvem_vitd3dose*' )
479
480             IF ( .NOT. ALLOCATED( vitd3_dose ) )  THEN
481                ALLOCATE( vitd3_dose(nysg:nyng,nxlg:nxrg) )
482             ENDIF
483             vitd3_dose = 0.0_wp
484
485          CASE DEFAULT
486             CONTINUE
487
488       END SELECT
489
490    ELSEIF ( mode == 'sum' )  THEN
491
492       SELECT CASE ( TRIM( variable ) )
493
494          CASE ( 'bio_mrt' )
495!
496!--          Consider the case 'bio_mrt' is called after some thermal index. In that case
497!            do_average_mrt will be .TRUE. leading to a double-averaging.
498             IF ( .NOT. do_average_mrt  .AND.  ALLOCATED( mrt_av_grid ) )  THEN
499
500                IF ( mrt_include_sw )  THEN
501                   DO  l = 1, nmrtbl
502                      i = mrtbl(ix,l)
503                      j = mrtbl(iy,l)
504                      k = mrtbl(iz,l)
505                      mrt_av_grid(k,j,i) = mrt_av_grid(k,j,i)  +                                   &
506                                    ( ( human_absorb * mrtinsw(l) +                                &
507                                    mrtinlw(l) ) / ( human_emiss * sigma_sb ) )**0.25_wp - degc_to_k
508                   ENDDO
509                ELSE
510                   DO  l = 1, nmrtbl
511                      i = mrtbl(ix,l)
512                      j = mrtbl(iy,l)
513                      k = mrtbl(iz,l)
514                      mrt_av_grid(k,j,i) = mrt_av_grid(k,j,i)  +                                   &
515                                    ( mrtinlw(l) / ( human_emiss * sigma_sb ) )**0.25_wp - degc_to_k
516                   ENDDO
517                ENDIF
518             ENDIF
519
520          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*', 'bio_mrt*' )
521!
522!--          Only continue if the current index is the one to trigger the input averaging, see
523!--          above.
524             IF ( average_trigger_perct  .AND.  TRIM( variable ) /= 'bio_perct*')    RETURN
525             IF ( average_trigger_utci   .AND.  TRIM( variable ) /= 'bio_utci*' )    RETURN
526             IF ( average_trigger_pet    .AND.  TRIM( variable ) /= 'bio_pet*'  )    RETURN
527             IF ( average_trigger_mrt    .AND.  TRIM( variable ) /= 'bio_mrt*'  )    RETURN
528!
529!--          Now memorize which of the input grids are not averaged by other modules. Set averaging
530!--          switch to .TRUE. and allocate the respective grid in that case.
531             IF ( .NOT. ALLOCATED( pt_av ) )  THEN  !< if not averaged by other module
532                ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
533                do_average_theta = .TRUE.  !< memorize, that bio is responsible
534                pt_av = 0.0_wp
535             ENDIF
536             IF ( ALLOCATED( pt_av )  .AND.  do_average_theta )  THEN
537                DO  i = nxl, nxr
538                   DO  j = nys, nyn
539                      DO  k = nzb, nzt+1
540                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
541                      ENDDO
542                   ENDDO
543                ENDDO
544             ENDIF
545
546             IF ( .NOT. ALLOCATED( q_av ) )  THEN
547                ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
548                do_average_q = .TRUE.
549                q_av = 0.0_wp
550             ENDIF
551             IF ( ALLOCATED( q_av )  .AND.  do_average_q )  THEN
552                DO  i = nxl, nxr
553                   DO  j = nys, nyn
554                      DO  k = nzb, nzt+1
555                         q_av(k,j,i) = q_av(k,j,i) + q(k,j,i)
556                      ENDDO
557                   ENDDO
558                ENDDO
559             ENDIF
560
561!
562!--          u_av, v_av and w_av are always allocated
563             IF ( .NOT. ALLOCATED( u_av ) )  THEN
564                ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
565                do_average_u = .TRUE.
566                u_av = 0.0_wp
567             ENDIF
568             IF ( ALLOCATED( u_av )  .AND.  do_average_u )  THEN
569                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
570                   DO  j = nysg, nyng
571                      DO  k = nzb, nzt+1
572                         u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
573                      ENDDO
574                   ENDDO
575                ENDDO
576             ENDIF
577
578             IF ( .NOT. ALLOCATED( v_av ) )  THEN
579                ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
580                do_average_v = .TRUE.
581                v_av = 0.0_wp
582             ENDIF
583             IF ( ALLOCATED( v_av )  .AND.  do_average_v )  THEN
584                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
585                   DO  j = nysg, nyng
586                      DO  k = nzb, nzt+1
587                         v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
588                      ENDDO
589                   ENDDO
590                ENDDO
591             ENDIF
592
593             IF ( .NOT. ALLOCATED( w_av ) )  THEN
594                ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
595                do_average_w = .TRUE.
596                w_av = 0.0_wp
597             ENDIF
598             IF ( ALLOCATED( w_av )  .AND.  do_average_w )  THEN
599                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
600                   DO  j = nysg, nyng
601                      DO  k = nzb, nzt+1
602                         w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
603                      ENDDO
604                   ENDDO
605                ENDDO
606             ENDIF
607
608             IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
609                ALLOCATE( mrt_av_grid(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
610                do_average_mrt = .TRUE.
611                mrt_av_grid = 0.0_wp
612             ENDIF
613             IF ( ALLOCATED( mrt_av_grid )  .AND.  do_average_mrt )  THEN
614
615                IF ( mrt_include_sw )  THEN
616                   DO  l = 1, nmrtbl
617                      i = mrtbl(ix,l)
618                      j = mrtbl(iy,l)
619                      k = mrtbl(iz,l)
620                      mrt_av_grid(k,j,i) = mrt_av_grid(k,j,i)  +                                   &
621                                    ( ( human_absorb * mrtinsw(l) +                                &
622                                    mrtinlw(l) ) /                                                 &
623                                      ( human_emiss  * sigma_sb ) )**0.25_wp - degc_to_k
624                   ENDDO
625                ELSE
626                   DO  l = 1, nmrtbl
627                      i = mrtbl(ix,l)
628                      j = mrtbl(iy,l)
629                      k = mrtbl(iz,l)
630                      mrt_av_grid(k,j,i) = mrt_av_grid(k,j,i) +                                    &
631                                    ( mrtinlw(l) /                                                 &
632                                    ( human_emiss * sigma_sb ) )**0.25_wp - degc_to_k
633                   ENDDO
634                ENDIF
635             ENDIF
636!
637!--       This is a cumulated dose. No mode == 'average' for this quantity.
638          CASE ( 'uvem_vitd3dose*' )
639             IF ( ALLOCATED( vitd3_dose ) )  THEN
640                DO  i = nxlg, nxrg
641                   DO  j = nysg, nyng
642                      vitd3_dose(j,i) = vitd3_dose(j,i) + vitd3_exposure(j,i)
643                   ENDDO
644                ENDDO
645             ENDIF
646
647          CASE DEFAULT
648             CONTINUE
649
650       END SELECT
651
652    ELSEIF ( mode == 'average' )  THEN
653
654       SELECT CASE ( TRIM( variable ) )
655
656          CASE ( 'bio_mrt' )
657!
658!--          Consider the case 'bio_mrt' is called after some thermal index. In that case
659!--          do_average_mrt will be .TRUE. leading to a double-averaging.
660             IF ( .NOT. do_average_mrt  .AND.  ALLOCATED( mrt_av_grid ) )  THEN
661                DO  i = nxl, nxr
662                   DO  j = nys, nyn
663                      DO  k = nzb, nzt+1
664                         mrt_av_grid(k,j,i) = mrt_av_grid(k,j,i) / REAL( average_count_3d, KIND=wp )
665                      ENDDO
666                   ENDDO
667                ENDDO
668             ENDIF
669
670          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*', 'bio_mrt*' )
671!
672!--          Only continue if update index, see above
673             IF ( average_trigger_perct  .AND.                                                     &
674                TRIM( variable ) /= 'bio_perct*' )  RETURN
675             IF ( average_trigger_utci  .AND.                                                      &
676                TRIM( variable ) /= 'bio_utci*' )  RETURN
677             IF ( average_trigger_pet   .AND.                                                      &
678                TRIM( variable ) /= 'bio_pet*' )  RETURN
679             IF ( average_trigger_mrt   .AND.                                                      &
680                TRIM( variable ) /= 'bio_mrt*' )  RETURN
681
682             IF ( ALLOCATED( pt_av )  .AND.  do_average_theta )  THEN
683                DO  i = nxl, nxr
684                   DO  j = nys, nyn
685                      DO  k = nzb, nzt+1
686                         pt_av(k,j,i) = pt_av(k,j,i) /  REAL( average_count_3d, KIND = wp )
687                      ENDDO
688                   ENDDO
689                ENDDO
690             ENDIF
691
692             IF ( ALLOCATED( q_av )  .AND.  do_average_q )  THEN
693                DO  i = nxl, nxr
694                   DO  j = nys, nyn
695                      DO  k = nzb, nzt+1
696                         q_av(k,j,i) = q_av(k,j,i) / REAL( average_count_3d, KIND = wp )
697                      ENDDO
698                   ENDDO
699                ENDDO
700             ENDIF
701
702             IF ( ALLOCATED( u_av )  .AND.  do_average_u )  THEN
703                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
704                   DO  j = nysg, nyng
705                      DO  k = nzb, nzt+1
706                         u_av(k,j,i) = u_av(k,j,i) / REAL( average_count_3d, KIND = wp )
707                      ENDDO
708                   ENDDO
709                ENDDO
710             ENDIF
711
712             IF ( ALLOCATED( v_av )  .AND.  do_average_v )  THEN
713                DO  i = nxlg, nxrg
714                   DO  j = nysg, nyng
715                      DO  k = nzb, nzt+1
716                         v_av(k,j,i) = v_av(k,j,i) / REAL( average_count_3d, KIND = wp )
717                      ENDDO
718                   ENDDO
719                ENDDO
720             ENDIF
721
722             IF ( ALLOCATED( w_av )  .AND.  do_average_w )  THEN
723                DO  i = nxlg, nxrg
724                   DO  j = nysg, nyng
725                      DO  k = nzb, nzt+1
726                         w_av(k,j,i) = w_av(k,j,i) / REAL( average_count_3d, KIND = wp )
727                      ENDDO
728                   ENDDO
729                ENDDO
730             ENDIF
731
732             IF ( ALLOCATED( mrt_av_grid )  .AND.  do_average_mrt )  THEN
733                DO  i = nxl, nxr
734                   DO  j = nys, nyn
735                      DO  k = nzb, nzt+1
736                         mrt_av_grid(k,j,i) = mrt_av_grid(k,j,i) / REAL( average_count_3d, KIND=wp )
737                      ENDDO
738                   ENDDO
739                ENDDO
740             ENDIF
741
742!
743!--     No averaging for UVEM SINce we are calculating a dose (only sum is calculated and saved to
744!--     av.nc file)
745        END SELECT
746
747    ENDIF
748
749
750 END SUBROUTINE bio_3d_data_averaging
751
752
753
754!--------------------------------------------------------------------------------------------------!
755! Description:
756! ------------
757!> Check data output for biometeorology model
758!--------------------------------------------------------------------------------------------------!
759 SUBROUTINE bio_check_data_output( var, unit, i, j, ilen, k )
760
761    USE control_parameters,                                                                        &
762        ONLY: data_output, message_string
763
764    IMPLICIT NONE
765
766    CHARACTER (LEN=*) ::  unit     !< The unit for the variable var
767    CHARACTER (LEN=*) ::  var      !< The variable in question
768
769    INTEGER(iwp), INTENT(IN) ::  i     !< Current element of data_output
770    INTEGER(iwp), INTENT(IN) ::  j     !< Average quantity? 0 = no, 1 = yes
771    INTEGER(iwp), INTENT(IN) ::  ilen  !< Length of current entry in data_output
772    INTEGER(iwp), INTENT(IN) ::  k     !< Output is xy mode? 0 = no, 1 = yes
773
774    SELECT CASE ( TRIM( var ) )
775!
776!--    Allocate a temporary array with the desired output dimensions.
777!--    Arrays for time-averaged thermal indices are also allocated here because they are not running
778!--    through the standard averaging procedure in bio_3d_data_averaging as the values of the
779!--    averaged thermal indices are derived in a SINgle step based on priorly averaged arrays (see
780!--    bio_calculate_thermal_index_maps).
781       CASE ( 'bio_mrt', 'bio_mrt*' )
782          unit = 'degree_C'
783          thermal_comfort = .TRUE.  !< enable thermal_comfort if user forgot to do so
784          IF ( .NOT. ALLOCATED( tmrt_grid ) )  THEN
785             ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
786             tmrt_grid = REAL( bio_fill_value, KIND = wp )
787          ENDIF
788          IF ( TRIM( var ) == 'bio_mrt*' )  THEN
789             do_calculate_mrt2d = .TRUE.
790          END IF
791
792       CASE ( 'bio_perct*' )
793          unit = 'degree_C'
794          thermal_comfort = .TRUE.
795          IF ( j == 0 )  THEN                !< if instantaneous input
796             do_calculate_perct = .TRUE.
797             IF ( .NOT. ALLOCATED( perct ) )  THEN
798                ALLOCATE( perct (nys:nyn,nxl:nxr) )
799                perct = REAL( bio_fill_value, KIND = wp )
800             ENDIF
801          ELSE                               !< if averaged input
802             do_calculate_perct_av = .TRUE.
803             IF ( .NOT. ALLOCATED( perct_av ) )  THEN
804                ALLOCATE( perct_av (nys:nyn,nxl:nxr) )
805                perct_av = REAL( bio_fill_value, KIND = wp )
806             ENDIF
807          ENDIF
808
809       CASE ( 'bio_utci*' )
810          unit = 'degree_C'
811          thermal_comfort = .TRUE.
812          IF ( j == 0 )  THEN
813             do_calculate_utci = .TRUE.
814             IF ( .NOT. ALLOCATED( utci ) )  THEN
815                ALLOCATE( utci (nys:nyn,nxl:nxr) )
816                utci = REAL( bio_fill_value, KIND = wp )
817             ENDIF
818          ELSE
819             do_calculate_utci_av = .TRUE.
820             IF ( .NOT. ALLOCATED( utci_av ) )  THEN
821                ALLOCATE( utci_av (nys:nyn,nxl:nxr) )
822                utci_av = REAL( bio_fill_value, KIND = wp )
823             ENDIF
824          ENDIF
825
826       CASE ( 'bio_pet*' )
827          unit = 'degree_C'
828          thermal_comfort = .TRUE.
829          IF ( j == 0 )  THEN
830             do_calculate_pet = .TRUE.
831             IF ( .NOT. ALLOCATED( pet ) )  THEN
832                ALLOCATE( pet (nys:nyn,nxl:nxr) )
833                pet = REAL( bio_fill_value, KIND = wp )
834             ENDIF
835          ELSE
836             do_calculate_pet_av = .TRUE.
837             IF ( .NOT. ALLOCATED( pet_av ) )  THEN
838                ALLOCATE( pet_av (nys:nyn,nxl:nxr) )
839                pet_av = REAL( bio_fill_value, KIND = wp )
840             ENDIF
841          ENDIF
842
843
844       CASE ( 'uvem_vitd3*' )
845!           IF ( .NOT. uv_exposure )  THEN
846!              message_string = 'output of "' // TRIM( var ) // '" requi' //                       &
847!                       'res a namelist &uvexposure_par'
848!              CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
849!           ENDIF
850          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
851             message_string = 'illegal value for data_output: "' //                                &
852                              TRIM( var ) // '" & only 2d-horizontal ' //                          &
853                              'cross sections are allowed for this value'
854             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
855          ENDIF
856          unit = 'IU/s'
857          IF ( .NOT. ALLOCATED( vitd3_exposure ) )  THEN
858             ALLOCATE( vitd3_exposure(nysg:nyng,nxlg:nxrg) )
859          ENDIF
860          vitd3_exposure = 0.0_wp
861
862       CASE ( 'uvem_vitd3dose*' )
863!           IF (  .NOT.  uv_exposure )  THEN
864!              message_string = 'output of "' // TRIM( var ) // '" requi' //     &
865!                       'res  a namelist &uvexposure_par'
866!              CALL message( 'uvem_check_data_output', 'UV0001', 1, 2, 0, 6, 0 )
867!           ENDIF
868          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
869             message_string = 'illegal value for data_output: "' //                                &
870                              TRIM( var ) // '" & only 2d-horizontal ' //                          &
871                              'cross sections are allowed for this value'
872             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
873          ENDIF
874          unit = 'IU/av-h'
875          IF ( .NOT. ALLOCATED( vitd3_dose ) )  THEN
876             ALLOCATE( vitd3_dose(nysg:nyng,nxlg:nxrg) )
877          ENDIF
878          vitd3_dose = 0.0_wp
879
880       CASE DEFAULT
881          unit = 'illegal'
882
883    END SELECT
884
885!
886!-- Further checks if thermal comfort output is desired.
887    IF ( thermal_comfort  .AND.  unit(1:6) == 'degree' )  THEN
888!
889!--    Break if required modules "radiation" is not available.
890       IF ( .NOT.  radiation )  THEN
891          message_string = 'output of "' // TRIM( var ) // '" require' // 's radiation = .TRUE.'
892          CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
893          unit = 'illegal'
894       ENDIF
895!
896!--    All "thermal_comfort" outputs except from 'bio_mrt' will also need  humidity input. Check
897!--    also for that.
898       IF ( TRIM( var ) /= 'bio_mrt' )  THEN
899          IF ( .NOT.  humidity )  THEN
900             message_string = 'The estimation of thermal comfort '    //                           &
901                              'requires air humidity information, but ' //                         &
902                              'humidity module is disabled!'
903             CALL message( 'check_parameters', 'PA0561', 1, 2, 0, 6, 0 )
904             unit = 'illegal'
905          ENDIF
906       ENDIF
907
908
909    ENDIF
910
911 END SUBROUTINE bio_check_data_output
912
913!--------------------------------------------------------------------------------------------------!
914! Description:
915! ------------
916!> Check parameters routine for biom module
917!> Currently unused but might come in handy for future checks?
918!--------------------------------------------------------------------------------------------------!
919 SUBROUTINE bio_check_parameters
920
921
922    IMPLICIT NONE
923
924
925 END SUBROUTINE bio_check_parameters
926
927
928!--------------------------------------------------------------------------------------------------!
929! Description:
930! ------------
931!> Subroutine defining 2D output variables
932!> data_output_2d 1188ff
933!--------------------------------------------------------------------------------------------------!
934 SUBROUTINE bio_data_output_2d( av, variable, found, grid, local_pf, two_d, nzb_do, nzt_do)
935
936
937    USE kinds
938
939
940    IMPLICIT NONE
941!
942!-- Input variables
943    CHARACTER (LEN=*), INTENT(IN) ::  variable    !< Char identifier to select var for output
944    INTEGER(iwp), INTENT(IN)      ::  av          !< Use averaged data? 0 = no, 1 = yes?
945    INTEGER(iwp), INTENT(IN)      ::  nzb_do      !< Unused. 2D. nz bottom to nz top
946    INTEGER(iwp), INTENT(IN)      ::  nzt_do      !< Unused.
947!
948!-- Output variables
949    CHARACTER (LEN=*), INTENT(OUT) ::  grid   !< Grid type (always "zu1" for biom)
950
951    LOGICAL, INTENT(OUT)           ::  found  !< Output found?
952    LOGICAL, INTENT(OUT)           ::  two_d  !< Flag parameter that indicates 2D variables,
953                                              !< horizontal cross sections, must be .TRUE. for thermal indices and uv
954    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< Temp. result grid to return
955!
956!-- Internal variables
957    INTEGER(iwp) ::  i        !< Running index, x-dir
958    INTEGER(iwp) ::  j        !< Running index, y-dir
959    INTEGER(iwp) ::  k        !< Running index, z-dir
960    INTEGER(iwp) ::  l        !< Running index, radiation grid
961
962
963    found = .TRUE.
964    local_pf = bio_fill_value
965
966    SELECT CASE ( TRIM( variable ) )
967
968
969        CASE ( 'bio_mrt_xy' )
970           grid = 'zu1'
971           two_d = .FALSE.  !< can be calculated for several levels
972           local_pf = REAL( bio_fill_value, KIND = wp )
973           DO  l = 1, nmrtbl
974              i = mrtbl(ix,l)
975              j = mrtbl(iy,l)
976              k = mrtbl(iz,l)
977              IF ( k < nzb_do  .OR.  k > nzt_do  .OR.  j < nys  .OR.                               &
978                 j > nyn  .OR.  i < nxl  .OR.  i > nxr )  CYCLE
979              IF ( av == 0 )  THEN
980                 IF ( mrt_include_sw )  THEN
981                    local_pf(i,j,k) = ( ( human_absorb * mrtinsw(l) +                              &
982                                      mrtinlw(l) ) /                                               &
983                                      ( human_emiss * sigma_sb ) )**0.25_wp - degc_to_k
984                 ELSE
985                    local_pf(i,j,k) = ( mrtinlw(l) /                                               &
986                                      ( human_emiss * sigma_sb ) )**0.25_wp - degc_to_k
987                 ENDIF
988              ELSE
989                 local_pf(i,j,k) = mrt_av_grid(k,j,i)
990              ENDIF
991           ENDDO
992
993        CASE ( 'bio_mrt*_xy' )        ! 2d-array
994           grid = 'zu1'
995           two_d = .TRUE.
996           IF ( av == 0 )  THEN
997              DO  i = nxl, nxr
998                 DO  j = nys, nyn
999                    local_pf(i,j,nzb+1) = tmrt_grid(j,i)
1000                 ENDDO
1001              ENDDO
1002           ELSE
1003              DO  i = nxl, nxr
1004                 DO  j = nys, nyn
1005                    local_pf(i,j,nzb+1) = tmrt_av_grid(j,i)
1006                 ENDDO
1007              ENDDO
1008           ENDIF
1009
1010
1011        CASE ( 'bio_perct*_xy' )        ! 2d-array
1012           grid = 'zu1'
1013           two_d = .TRUE.
1014           IF ( av == 0 )  THEN
1015              DO  i = nxl, nxr
1016                 DO  j = nys, nyn
1017                    local_pf(i,j,nzb+1) = perct(j,i)
1018                 ENDDO
1019              ENDDO
1020           ELSE
1021              DO  i = nxl, nxr
1022                 DO  j = nys, nyn
1023                    local_pf(i,j,nzb+1) = perct_av(j,i)
1024                 ENDDO
1025              ENDDO
1026           ENDIF
1027
1028
1029        CASE ( 'bio_utci*_xy' )        ! 2d-array
1030           grid = 'zu1'
1031           two_d = .TRUE.
1032           IF ( av == 0 )  THEN
1033              DO  i = nxl, nxr
1034                 DO  j = nys, nyn
1035                    local_pf(i,j,nzb+1) = utci(j,i)
1036                 ENDDO
1037              ENDDO
1038           ELSE
1039              DO  i = nxl, nxr
1040                 DO  j = nys, nyn
1041                    local_pf(i,j,nzb+1) = utci_av(j,i)
1042                 ENDDO
1043              ENDDO
1044           ENDIF
1045
1046
1047        CASE ( 'bio_pet*_xy' )        ! 2d-array
1048           grid = 'zu1'
1049           two_d = .TRUE.
1050           IF ( av == 0 )  THEN
1051              DO  i = nxl, nxr
1052                 DO  j = nys, nyn
1053                    local_pf(i,j,nzb+1) = pet(j,i)
1054                 ENDDO
1055              ENDDO
1056           ELSE
1057              DO  i = nxl, nxr
1058                 DO  j = nys, nyn
1059                    local_pf(i,j,nzb+1) = pet_av(j,i)
1060                 ENDDO
1061              ENDDO
1062           ENDIF
1063
1064!
1065!--    Before data is transfered to local_pf, transfer is in 2D dummy variable and exchange ghost
1066!--    points therein. However, at this point this is only required for instantaneous arrays,
1067!--    time-averaged quantities are already exchanged.
1068       CASE ( 'uvem_vitd3*_xy' )        ! 2d-array
1069          IF ( av == 0 )  THEN
1070             DO  i = nxl, nxr
1071                DO  j = nys, nyn
1072                   local_pf(i,j,nzb+1) = vitd3_exposure(j,i)
1073                ENDDO
1074             ENDDO
1075          ENDIF
1076
1077          two_d = .TRUE.
1078          grid = 'zu1'
1079
1080       CASE ( 'uvem_vitd3dose*_xy' )        ! 2d-array
1081          IF ( av == 1 )  THEN
1082             DO  i = nxl, nxr
1083                DO  j = nys, nyn
1084                   local_pf(i,j,nzb+1) = vitd3_dose(j,i)
1085                ENDDO
1086             ENDDO
1087          ENDIF
1088
1089          two_d = .TRUE.
1090          grid = 'zu1'
1091
1092
1093       CASE DEFAULT
1094          found = .FALSE.
1095          grid  = 'none'
1096
1097    END SELECT
1098
1099
1100 END SUBROUTINE bio_data_output_2d
1101
1102
1103!--------------------------------------------------------------------------------------------------!
1104! Description:
1105! ------------
1106!> Subroutine defining 3D output variables (dummy, always 2d!)
1107!> data_output_3d 709ff
1108!--------------------------------------------------------------------------------------------------!
1109 SUBROUTINE bio_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
1110
1111    USE indices
1112
1113    USE kinds
1114
1115
1116    IMPLICIT NONE
1117!
1118!-- Input variables
1119    CHARACTER (LEN=*), INTENT(IN) ::  variable   !< Char identifier to select var for output
1120
1121    INTEGER(iwp), INTENT(IN) ::  av       !< Use averaged data? 0 = no, 1 = yes?
1122    INTEGER(iwp), INTENT(IN) ::  nzb_do   !< Unused. 2D. nz bottom to nz top
1123    INTEGER(iwp), INTENT(IN) ::  nzt_do   !< Unused.
1124!
1125!-- Output variables
1126    LOGICAL, INTENT(OUT) ::  found   !< Output found?
1127
1128    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< Temp. result grid to return
1129!
1130!-- Internal variables
1131    INTEGER(iwp) ::  l    !< Running index, radiation grid
1132    INTEGER(iwp) ::  i    !< Running index, x-dir
1133    INTEGER(iwp) ::  j    !< Running index, y-dir
1134    INTEGER(iwp) ::  k    !< Running index, z-dir
1135
1136!    REAL(wp) ::  mrt  !< Buffer for mean radiant temperature
1137
1138    found = .TRUE.
1139
1140    SELECT CASE ( TRIM( variable ) )
1141
1142        CASE ( 'bio_mrt' )
1143            local_pf = REAL( bio_fill_value, KIND = sp )
1144            DO  l = 1, nmrtbl
1145               i = mrtbl(ix,l)
1146               j = mrtbl(iy,l)
1147               k = mrtbl(iz,l)
1148               IF ( k < nzb_do  .OR.  k > nzt_do  .OR.  j < nys  .OR.                              &
1149                  j > nyn  .OR.  i < nxl  .OR.  i > nxr )  CYCLE
1150               IF ( av == 0 )  THEN
1151                  IF ( mrt_include_sw )  THEN
1152                     local_pf(i,j,k) = REAL( ( ( human_absorb * mrtinsw(l) +                       &
1153                                       mrtinlw(l) ) /                                              &
1154                                       ( human_emiss * sigma_sb ) )**0.25_wp - degc_to_k,          &
1155                                         KIND = sp )
1156                  ELSE
1157                     local_pf(i,j,k) = REAL( ( mrtinlw(l) /                                        &
1158                                       ( human_emiss * sigma_sb ) )**0.25_wp - degc_to_k,          &
1159                                         KIND = sp )
1160                  ENDIF
1161               ELSE
1162                  local_pf(i,j,k) = REAL( mrt_av_grid(k,j,i), KIND = sp )
1163               ENDIF
1164            ENDDO
1165
1166       CASE DEFAULT
1167          found = .FALSE.
1168
1169    END SELECT
1170
1171 END SUBROUTINE bio_data_output_3d
1172
1173!--------------------------------------------------------------------------------------------------!
1174! Description:
1175! ------------
1176!> Subroutine defining appropriate grid for netcdf variables.
1177!> It is called out from subroutine netcdf_interface_mod.
1178!> netcdf_interface_mod 918ff
1179!--------------------------------------------------------------------------------------------------!
1180 SUBROUTINE bio_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
1181
1182    IMPLICIT NONE
1183!
1184!-- Input variables
1185    CHARACTER (LEN=*), INTENT(IN)  ::  var      !< Name of output variable
1186!
1187!-- Output variables
1188    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
1189    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !< y grid of output variable
1190    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !< z grid of output variable
1191
1192    LOGICAL, INTENT(OUT)           ::  found    !< Flag if output var is found
1193!
1194!-- Local variables
1195    INTEGER(iwp) :: l     !< Length of the var array
1196
1197    LOGICAL      :: is2d  !< Var is 2d?
1198
1199    found  = .FALSE.
1200    grid_x = 'none'
1201    grid_y = 'none'
1202    grid_z = 'none'
1203
1204    l = MAX( 2, LEN_TRIM( var ) )
1205    is2d = ( var(l-1:l) == 'xy' )
1206
1207    IF ( var(1:4) == 'bio_' )  THEN
1208       found  = .TRUE.
1209       grid_x = 'x'
1210       grid_y = 'y'
1211       grid_z = 'zu'
1212       IF ( is2d  .AND.  var(1:7) /= 'bio_mrt' )  grid_z = 'zu1'
1213    ENDIF
1214
1215    IF ( is2d  .AND.  var(1:4) == 'uvem' )  THEN
1216       grid_x = 'x'
1217       grid_y = 'y'
1218       grid_z = 'zu1'
1219    ENDIF
1220
1221 END SUBROUTINE bio_define_netcdf_grid
1222
1223!--------------------------------------------------------------------------------------------------!
1224! Description:
1225! ------------
1226!> Header output for biom module
1227!> header 982
1228!--------------------------------------------------------------------------------------------------!
1229 SUBROUTINE bio_header( io )
1230
1231    IMPLICIT NONE
1232!
1233!-- Input variables
1234    INTEGER(iwp), INTENT(IN) ::  io           !< Unit of the output file
1235!
1236!-- Internal variables
1237    CHARACTER (LEN=86) ::  output_height_chr  !< String for output height
1238
1239    WRITE( output_height_chr, '(F8.1,7X)' )  bio_output_height
1240!
1241!-- Write biom header
1242    WRITE( io, 1 )
1243    WRITE( io, 2 )  TRIM( output_height_chr )
1244    WRITE( io, 3 )  TRIM( ACHAR( bio_cell_level ) )
1245
12461   FORMAT (//' Human thermal comfort module information:'/                                        &
1247              ' ------------------------------'/)
12482   FORMAT ('    --> All indices calculated for a height of (m): ', A )
12493   FORMAT ('    --> This corresponds to cell level : ', A )
1250
1251 END SUBROUTINE bio_header
1252
1253
1254!--------------------------------------------------------------------------------------------------!
1255! Description:
1256! ------------
1257!> Initialization of the HTCM
1258!> init_3d_model 1987ff
1259!--------------------------------------------------------------------------------------------------!
1260 SUBROUTINE bio_init
1261
1262    USE netcdf_data_input_mod,                                                                     &
1263        ONLY:  netcdf_data_input_uvem
1264
1265    IMPLICIT NONE
1266!
1267!-- Internal vriables
1268    REAL ( wp )  :: height  !< current height in meters
1269
1270    IF ( debug_output )  CALL debug_message( 'bio_init', 'start' )
1271!
1272!-- Determine cell level corresponding to 1.1 m above ground level (gravimetric center of sample
1273!-- human)
1274
1275    bio_cell_level = 0_iwp
1276    bio_output_height = 0.5_wp * dz(1)
1277    height = 0.0_wp
1278
1279    bio_cell_level = INT( 1.099_wp / dz(1) )
1280    bio_output_height = bio_output_height + bio_cell_level * dz(1)
1281!
1282!-- Set radiation level if not done by user
1283    IF ( mrt_nlevels == 0 )  THEN
1284       mrt_nlevels = bio_cell_level + 1_iwp
1285    ENDIF
1286!
1287!-- Init UVEM and load lookup tables
1288    IF ( uv_exposure )  CALL netcdf_data_input_uvem
1289
1290    IF ( debug_output )  CALL debug_message( 'bio_init', 'end' )
1291
1292 END SUBROUTINE bio_init
1293
1294
1295!--------------------------------------------------------------------------------------------------!
1296! Description:
1297! ------------
1298!> Checks done after the Initialization
1299!--------------------------------------------------------------------------------------------------!
1300 SUBROUTINE bio_init_checks
1301
1302    USE control_parameters,                                                                        &
1303        ONLY: message_string
1304
1305    IF ( (.NOT. radiation_interactions) .AND. ( thermal_comfort ) )  THEN
1306       message_string = 'The mrt calculation requires ' //                                         &
1307                        'enabled radiation_interactions but it ' //                                &
1308                        'is disabled!'
1309       CALL message( 'bio_init_checks', 'PAHU03', 1, 2, 0, 6, 0 )
1310    ENDIF
1311
1312
1313 END SUBROUTINE bio_init_checks
1314
1315
1316!--------------------------------------------------------------------------------------------------!
1317! Description:
1318! ------------
1319!> Parin for &biometeorology_parameters for reading biomet parameters
1320!--------------------------------------------------------------------------------------------------!
1321 SUBROUTINE bio_parin
1322
1323    IMPLICIT NONE
1324
1325!
1326!-- Internal variables
1327    CHARACTER (LEN=80) ::  line  !< Dummy string for current line in parameter file
1328
1329    NAMELIST /biometeorology_parameters/  clothing,                                                &
1330                                          consider_obstructions,                                   &
1331                                          orientation_angle,                                       &
1332                                          sun_in_south,                                            &
1333                                          thermal_comfort,                                         &
1334                                          turn_to_sun,                                             &
1335                                          uv_exposure
1336
1337
1338!-- Try to find biometeorology_parameters namelist
1339    REWIND ( 11 )
1340    line = ' '
1341    DO WHILE ( INDEX( line, '&biometeorology_parameters' ) == 0 )
1342       READ ( 11, '(A)', END = 20 )  line
1343    ENDDO
1344    BACKSPACE ( 11 )
1345
1346!
1347!-- Read biometeorology_parameters namelist
1348    READ ( 11, biometeorology_parameters, ERR = 10, END = 20 )
1349
1350!
1351!-- Set flag that indicates that the biomet_module is switched on
1352    biometeorology = .TRUE.
1353
1354    GOTO 20
1355
1356!
1357!-- In case of error
1358 10 BACKSPACE( 11 )
1359    READ( 11 , '(A)') line
1360    CALL parin_fail_message( 'biometeorology_parameters', line )
1361
1362!
1363!-- Complete
1364 20 CONTINUE
1365
1366
1367 END SUBROUTINE bio_parin
1368
1369!--------------------------------------------------------------------------------------------------!
1370! Description:
1371! ------------
1372!> Read module-specific global restart data (Fortran binary format).
1373!--------------------------------------------------------------------------------------------------!
1374 SUBROUTINE bio_rrd_global_ftn( found )
1375
1376    USE control_parameters,                                                                        &
1377        ONLY:  length, restart_string
1378
1379
1380    IMPLICIT NONE
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!
1390!--    Read control flags to determine if input grids need to be averaged.
1391       CASE ( 'do_average_theta' )
1392          READ ( 13 )  do_average_theta
1393
1394       CASE ( 'do_average_q' )
1395          READ ( 13 )  do_average_q
1396
1397       CASE ( 'do_average_u' )
1398          READ ( 13 )  do_average_u
1399
1400       CASE ( 'do_average_v' )
1401          READ ( 13 )  do_average_v
1402
1403       CASE ( 'do_average_w' )
1404          READ ( 13 )  do_average_w
1405
1406       CASE ( 'do_average_mrt' )
1407          READ ( 13 )  do_average_mrt
1408
1409!
1410!--    Read control flags to determine which thermal index needs to trigger averaging.
1411       CASE ( 'average_trigger_perct' )
1412          READ ( 13 )  average_trigger_perct
1413
1414       CASE ( 'average_trigger_utci' )
1415          READ ( 13 )  average_trigger_utci
1416
1417       CASE ( 'average_trigger_pet' )
1418          READ ( 13 )  average_trigger_pet
1419
1420       CASE ( 'average_trigger_mrt' )
1421          READ ( 13 )  average_trigger_mrt
1422
1423
1424       CASE DEFAULT
1425
1426          found = .FALSE.
1427
1428    END SELECT
1429
1430
1431 END SUBROUTINE bio_rrd_global_ftn
1432
1433
1434!--------------------------------------------------------------------------------------------------!
1435! Description:
1436! ------------
1437!> Read module-specific global restart data (MPI-IO).
1438!--------------------------------------------------------------------------------------------------!
1439 SUBROUTINE bio_rrd_global_mpi
1440
1441
1442!
1443!--  Read control flags to determine if input grids need to be averaged
1444     CALL rrd_mpi_io( 'do_average_theta', do_average_theta )
1445     CALL rrd_mpi_io( 'do_average_q', do_average_q )
1446     CALL rrd_mpi_io( 'do_average_u', do_average_u )
1447     CALL rrd_mpi_io( 'do_average_v', do_average_v )
1448     CALL rrd_mpi_io( 'do_average_w', do_average_w )
1449     CALL rrd_mpi_io( 'do_average_mrt', do_average_mrt )
1450!
1451!--  Rad control flags to determine which thermal index needs to trigger averaging
1452     CALL rrd_mpi_io( 'average_trigger_perct', average_trigger_perct )
1453     CALL rrd_mpi_io( 'average_trigger_utci', average_trigger_utci )
1454     CALL rrd_mpi_io( 'average_trigger_pet', average_trigger_pet )
1455     CALL rrd_mpi_io( 'average_trigger_mrt', average_trigger_mrt )
1456
1457 END SUBROUTINE bio_rrd_global_mpi
1458
1459
1460!--------------------------------------------------------------------------------------------------!
1461! Description:
1462! ------------
1463!> Read module-specific local restart data arrays (Fortran binary format).
1464!--------------------------------------------------------------------------------------------------!
1465 SUBROUTINE bio_rrd_local_ftn( found )
1466
1467
1468    USE control_parameters,                                                                        &
1469        ONLY:  length, restart_string
1470
1471
1472    IMPLICIT NONE
1473
1474
1475    LOGICAL, INTENT(OUT) ::  found      !< variable found? yes = .TRUE., no = .FALSE.
1476
1477    found = .TRUE.
1478
1479
1480    SELECT CASE ( restart_string(1:length) )
1481
1482       CASE ( 'mrt_av_grid' )
1483          IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
1484             ALLOCATE( mrt_av_grid(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1485             mrt_av_grid = 0.0_wp
1486          ENDIF
1487          READ ( 13 )  mrt_av_grid
1488
1489       CASE DEFAULT
1490
1491          found = .FALSE.
1492
1493    END SELECT
1494
1495
1496 END SUBROUTINE bio_rrd_local_ftn
1497
1498
1499!------------------------------------------------------------------------------!
1500! Description:
1501! ------------
1502!> Read module-specific local restart data arrays (Fortran binary format).
1503!------------------------------------------------------------------------------!
1504 SUBROUTINE bio_rrd_local_mpi
1505
1506    USE control_parameters
1507
1508    USE indices
1509
1510    USE kinds
1511
1512
1513    IMPLICIT NONE
1514
1515    LOGICAL      ::  array_found  !<
1516
1517    CALL rd_mpi_io_check_array( 'mrt_av_grid' , found = array_found )
1518    IF ( array_found )  THEN
1519       IF ( .NOT. ALLOCATED( mrt_av_grid ) )  ALLOCATE( mrt_av_grid(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1520       CALL rrd_mpi_io( 'mrt_av_grid', mrt_av_grid )
1521    ENDIF
1522
1523 END SUBROUTINE bio_rrd_local_mpi
1524
1525
1526!--------------------------------------------------------------------------------------------------!
1527! Description:
1528! ------------
1529!> Write global restart data for the biometeorology module.
1530!--------------------------------------------------------------------------------------------------!
1531 SUBROUTINE bio_wrd_global
1532
1533    IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
1534
1535       CALL wrd_write_string( 'do_average_theta' )
1536       WRITE ( 14 )  do_average_theta
1537       CALL wrd_write_string( 'do_average_q' )
1538       WRITE ( 14 )  do_average_q
1539       CALL wrd_write_string( 'do_average_u' )
1540       WRITE ( 14 )  do_average_u
1541       CALL wrd_write_string( 'do_average_v' )
1542       WRITE ( 14 )  do_average_v
1543       CALL wrd_write_string( 'do_average_w' )
1544       WRITE ( 14 )  do_average_w
1545       CALL wrd_write_string( 'do_average_mrt' )
1546       WRITE ( 14 )  do_average_mrt
1547       CALL wrd_write_string( 'average_trigger_perct' )
1548       WRITE ( 14 )  average_trigger_perct
1549       CALL wrd_write_string( 'average_trigger_utci' )
1550       WRITE ( 14 )  average_trigger_utci
1551       CALL wrd_write_string( 'average_trigger_pet' )
1552       WRITE ( 14 )  average_trigger_pet
1553       CALL wrd_write_string( 'average_trigger_mrt' )
1554       WRITE ( 14 )  average_trigger_mrt
1555
1556    ELSEIF ( TRIM( restart_data_format_output(1:3) ) == 'mpi' )  THEN
1557
1558       CALL wrd_mpi_io( 'do_average_theta', do_average_theta )
1559       CALL wrd_mpi_io( 'do_average_q', do_average_q )
1560       CALL wrd_mpi_io( 'do_average_u', do_average_u )
1561       CALL wrd_mpi_io( 'do_average_v', do_average_v )
1562       CALL wrd_mpi_io( 'do_average_w', do_average_w )
1563       CALL wrd_mpi_io( 'do_average_mrt', do_average_mrt )
1564       CALL wrd_mpi_io( 'average_trigger_perct', average_trigger_perct )
1565       CALL wrd_mpi_io( 'average_trigger_utci', average_trigger_utci )
1566       CALL wrd_mpi_io( 'average_trigger_pet', average_trigger_pet )
1567       CALL wrd_mpi_io( 'average_trigger_mrt', average_trigger_mrt )
1568
1569    ENDIF
1570
1571 END SUBROUTINE bio_wrd_global
1572
1573
1574!--------------------------------------------------------------------------------------------------!
1575! Description:
1576! ------------
1577!> Write local restart data for the biometeorology module.
1578!--------------------------------------------------------------------------------------------------!
1579 SUBROUTINE bio_wrd_local
1580
1581    IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
1582
1583       IF ( ALLOCATED( mrt_av_grid ) )  THEN
1584          CALL wrd_write_string( 'mrt_av_grid' )
1585          WRITE ( 14 )  mrt_av_grid
1586       ENDIF
1587
1588    ELSEIF ( TRIM( restart_data_format_output(1:3) ) == 'mpi' )  THEN
1589
1590       IF ( ALLOCATED( mrt_av_grid ) )  CALL wrd_mpi_io( 'mrt_av_grid', mrt_av_grid )
1591
1592    ENDIF
1593
1594 END SUBROUTINE bio_wrd_local
1595
1596!--------------------------------------------------------------------------------------------------!
1597! Description:
1598! ------------
1599!> Calculate biometeorology MRT for all 2D grid
1600!--------------------------------------------------------------------------------------------------!
1601 SUBROUTINE bio_calculate_mrt_grid ( av )
1602
1603    IMPLICIT NONE
1604
1605    LOGICAL, INTENT(IN)         ::  av    !< use averaged input?
1606!
1607!-- Internal variables
1608    INTEGER(iwp)                ::  i     !< Running index, x-dir, radiation coordinates
1609    INTEGER(iwp)                ::  j     !< Running index, y-dir, radiation coordinates
1610    INTEGER(iwp)                ::  k     !< Running index, y-dir, radiation coordinates
1611    INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
1612
1613
1614!
1615!-- We need to differentiate if averaged input is desired (av == .TRUE.) or not.
1616    IF ( av )  THEN
1617!
1618!--    Make sure tmrt_av_grid is present and initialize with the fill value
1619       IF ( .NOT. ALLOCATED( tmrt_av_grid ) )  THEN
1620          ALLOCATE( tmrt_av_grid(nys:nyn,nxl:nxr) )
1621       ENDIF
1622       tmrt_av_grid = REAL( bio_fill_value, KIND = wp )
1623
1624!
1625!--    mrt_av_grid should always be allcoated here, but better make sure ist actually is.
1626       IF ( ALLOCATED( mrt_av_grid ) )  THEN
1627!
1628!--       Iterate over the radiation grid (radiation coordinates) and fill the tmrt_av_grid
1629!--       (x, y coordinates) where appropriate: tmrt_av_grid is written for all i / j if level (k)
1630!--       matches output height.
1631          DO  l = 1, nmrtbl
1632             i = mrtbl(ix,l)
1633             j = mrtbl(iy,l)
1634             k = mrtbl(iz,l)
1635             IF ( k - topo_top_ind(j,i,0) == bio_cell_level + 1_iwp)  THEN
1636!
1637!--             Averaging was done before, so we can just copy the result here.
1638                tmrt_av_grid(j,i) = mrt_av_grid(k,j,i)
1639
1640             ENDIF
1641          ENDDO
1642       ENDIF
1643
1644!
1645!-- In case instantaneous input is desired, mrt values will be re-calculated.
1646    ELSE
1647!
1648!--    Calculate biometeorology MRT from local radiation fluxes calculated by RTM and assign into 2D
1649!--    grid. Depending on selected output quantities, tmrt_grid might not have been allocated in
1650!--    bio_check_data_output yet.
1651       IF ( .NOT. ALLOCATED( tmrt_grid ) )  THEN
1652          ALLOCATE( tmrt_grid (nys:nyn,nxl:nxr) )
1653       ENDIF
1654       tmrt_grid = REAL( bio_fill_value, KIND = wp )
1655
1656       DO  l = 1, nmrtbl
1657          i = mrtbl(ix,l)
1658          j = mrtbl(iy,l)
1659          k = mrtbl(iz,l)
1660          IF ( k - topo_top_ind(j,i,0) == bio_cell_level + 1_iwp)  THEN
1661             IF ( mrt_include_sw )  THEN
1662                tmrt_grid(j,i) = ( ( human_absorb * mrtinsw(l) +                                   &
1663                                 mrtinlw(l) )  /                                                   &
1664                                 ( human_emiss * sigma_sb ) )**0.25_wp -                           &
1665                                 degc_to_k
1666             ELSE
1667                tmrt_grid(j,i) = ( mrtinlw(l)  /                                                   &
1668                                 ( human_emiss * sigma_sb ) )**0.25_wp -                           &
1669                                 degc_to_k
1670             ENDIF
1671          ENDIF
1672       ENDDO
1673    ENDIF
1674
1675 END SUBROUTINE bio_calculate_mrt_grid
1676
1677
1678!--------------------------------------------------------------------------------------------------!
1679! Description:
1680! ------------
1681!> Calculate static thermal indices for 2D grid point i, j
1682!--------------------------------------------------------------------------------------------------!
1683 SUBROUTINE bio_get_thermal_index_input_ij( average_input, i, j, ta, vp, ws, pair, tmrt )
1684
1685    IMPLICIT NONE
1686!
1687!-- Input variables
1688    LOGICAL,      INTENT ( IN ) ::  average_input  !< Determine averaged input conditions?
1689    INTEGER(iwp), INTENT ( IN ) ::  i     !< Running index, x-dir
1690    INTEGER(iwp), INTENT ( IN ) ::  j     !< Running index, y-dir
1691!
1692!-- Output parameters
1693    REAL(wp), INTENT ( OUT )    ::  pair  !< Air pressure                    (hPa)
1694    REAL(wp), INTENT ( OUT )    ::  ta    !< Air temperature                 (degree_C)
1695    REAL(wp), INTENT ( OUT )    ::  tmrt  !< Mean radiant temperature        (degree_C)
1696    REAL(wp), INTENT ( OUT )    ::  vp    !< Vapour pressure                 (hPa)
1697    REAL(wp), INTENT ( OUT )    ::  ws    !< Wind speed    (local level)     (m/s)
1698!
1699!-- Internal variables
1700    INTEGER(iwp)                ::  k     !< Running index, z-dir
1701    INTEGER(iwp)                ::  k_wind  !< Running index, z-dir, wind speed only
1702
1703    REAL(wp)                    ::  vp_sat  !< Saturation vapor pressure     (hPa)
1704
1705!
1706!-- Determine cell level closest to 1.1m above ground by making use of truncation due to int cast.
1707    k = INT( topo_top_ind(j,i,0) + bio_cell_level )  !< Vertical cell center closest to 1.1m
1708
1709!
1710!-- Avoid non-representative horizontal u and v of 0.0 m/s too close to ground
1711    k_wind = k
1712    IF ( bio_cell_level < 1_iwp )  THEN
1713       k_wind = k + 1_iwp
1714    ENDIF
1715!
1716!-- Determine local values:
1717    IF ( average_input )  THEN
1718!
1719!--    Calculate ta from Tp assuming dry adiabatic laps rate
1720       ta = bio_fill_value
1721       IF ( ALLOCATED( pt_av ) )  THEN
1722          ta = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * ( k + 0.5_wp ) ) - degc_to_k
1723       ENDIF
1724
1725       vp = bio_fill_value
1726       IF ( humidity  .AND.  ALLOCATED( q_av ) )  THEN
1727          vp = q_av(k,j,i)
1728       ENDIF
1729
1730       ws = bio_fill_value
1731       IF ( ALLOCATED( u_av )  .AND.  ALLOCATED( v_av )  .AND.                                     &
1732          ALLOCATED( w_av ) )  THEN
1733             ws = ( 0.5_wp * ABS( u_av(k_wind,j,i) + u_av(k_wind,j,i+1) ) +                        &
1734                    0.5_wp * ABS( v_av(k_wind,j,i) + v_av(k_wind,j+1,i) ) +                        &
1735                    0.5_wp * ABS( w_av(k_wind,j,i) + w_av(k_wind+1,j,i) ) )
1736       ENDIF
1737    ELSE
1738!
1739!--    Calculate ta from Tp assuming dry adiabatic laps rate
1740       ta = pt(k,j,i) - ( 0.0098_wp * dz(1) * (  k + 0.5_wp ) ) - degc_to_k
1741
1742       vp = bio_fill_value
1743       IF ( humidity )  THEN
1744          vp = q(k,j,i)
1745       ENDIF
1746
1747       ws = ( 0.5_wp * ABS( u(k_wind,j,i) + u(k_wind,j,i+1) )  +                                   &
1748              0.5_wp * ABS( v(k_wind,j,i) + v(k_wind,j+1,i) )  +                                   &
1749              0.5_wp * ABS( w(k_wind,j,i) + w(k_wind+1,j,i) ) )
1750
1751    ENDIF
1752!
1753!-- Local air pressure
1754    pair = surface_pressure
1755!
1756!-- Calculate water vapour pressure at saturation and convert to hPa.
1757!-- The magnus formula is limited to temperatures up to 333.15 K to avoid negative values of vp_sat.
1758    IF ( vp > -998.0_wp )  THEN
1759       vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) )
1760       vp  = vp * pair / ( vp + 0.622_wp )
1761       IF ( vp > vp_sat )  vp = vp_sat
1762    ENDIF
1763!
1764!-- Local mtr value at [i,j]
1765    tmrt = bio_fill_value  !< this can be a valid result (e.g. for inside some ostacle)
1766    IF ( .NOT. average_input )  THEN
1767!
1768!--    Use MRT from RTM precalculated in tmrt_grid
1769       tmrt = tmrt_grid(j,i)
1770    ELSE
1771       tmrt = tmrt_av_grid(j,i)
1772    ENDIF
1773
1774 END SUBROUTINE bio_get_thermal_index_input_ij
1775
1776
1777!--------------------------------------------------------------------------------------------------!
1778! Description:
1779! ------------
1780!> Calculate static thermal indices for any point within a 2D grid time_integration.f90: 1065ff
1781!--------------------------------------------------------------------------------------------------!
1782 SUBROUTINE bio_calculate_thermal_index_maps( av )
1783
1784    IMPLICIT NONE
1785!
1786!-- Input attributes
1787    LOGICAL, INTENT ( IN ) ::  av  !< Calculate based on averaged input conditions?
1788!
1789!-- Internal variables
1790    INTEGER(iwp) ::  i, j     !< Running index
1791
1792    REAL(wp) ::  clo          !< Clothing index                (no dimension)
1793    REAL(wp) ::  pair         !< Air pressure                     (hPa)
1794    REAL(wp) ::  perct_ij     !< Perceived temperature            (degree_C)
1795    REAL(wp) ::  pet_ij       !< Physiologically equivalent temperature  (degree_C)
1796    REAL(wp) ::  ta           !< Air temperature                  (degree_C)
1797    REAL(wp) ::  tmrt_ij      !< Mean radiant temperature         (degree_C)
1798    REAL(wp) ::  utci_ij      !< Universal thermal climate index  (degree_C)
1799    REAL(wp) ::  vp           !< Vapour pressure                  (hPa)
1800    REAL(wp) ::  ws           !< Wind speed    (local level)      (m/s)
1801
1802!
1803!-- Check if some thermal index is desired. Don't do anything if, e.g. only bio_mrt is desired.
1804    IF ( do_calculate_perct    .OR.  do_calculate_perct_av  .OR.  do_calculate_utci    .OR.        &
1805         do_calculate_utci_av  .OR.  do_calculate_pet       .OR.  do_calculate_pet_av  .OR.        &
1806         do_calculate_mrt2d )  THEN
1807
1808!
1809!--    fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...)
1810       CALL bio_calculate_mrt_grid ( av )
1811
1812       DO  i = nxl, nxr
1813          DO  j = nys, nyn
1814!
1815!--          Determine local input conditions
1816             tmrt_ij = bio_fill_value
1817             vp      = bio_fill_value
1818!
1819!--          Determine local meteorological conditions
1820             CALL bio_get_thermal_index_input_ij ( av, i, j, ta, vp, ws, pair, tmrt_ij )
1821!
1822!--          Only proceed if input is available
1823             pet_ij   = bio_fill_value   !< set fail value, e.g. valid for
1824             perct_ij = bio_fill_value   !< within some obstacle
1825             utci_ij  = bio_fill_value
1826             IF ( .NOT. ( tmrt_ij <= -998.0_wp  .OR.  vp <= -998.0_wp  .OR.   ws <= -998.0_wp  .OR.&
1827                  ta <= -998.0_wp ) )  THEN
1828!
1829!--             Calculate static thermal indices based on local tmrt
1830                clo = bio_fill_value
1831
1832                IF ( do_calculate_perct  .OR.  do_calculate_perct_av )  THEN
1833!
1834!--                Estimate local perceived temperature
1835                   CALL calculate_perct_static( ta, vp, ws, tmrt_ij, pair, clo, perct_ij )
1836                ENDIF
1837
1838                IF ( do_calculate_utci  .OR.  do_calculate_utci_av )  THEN
1839!
1840!--                Estimate local universal thermal climate index
1841                   CALL calculate_utci_static( ta, vp, ws, tmrt_ij, bio_output_height, utci_ij )
1842                ENDIF
1843
1844                IF ( do_calculate_pet  .OR.  do_calculate_pet_av )  THEN
1845!
1846!--                Estimate local physiologically equivalent temperature
1847                   CALL calculate_pet_static( ta, vp, ws, tmrt_ij, pair, pet_ij )
1848                ENDIF
1849             ENDIF
1850
1851
1852             IF ( av )  THEN
1853!
1854!--             Write results for selected averaged indices
1855                IF ( do_calculate_perct_av )  THEN
1856                   perct_av(j, i) = perct_ij
1857                ENDIF
1858                IF ( do_calculate_utci_av )  THEN
1859                   utci_av(j, i) = utci_ij
1860                ENDIF
1861                IF ( do_calculate_pet_av )  THEN
1862                   pet_av(j, i)  = pet_ij
1863                ENDIF
1864             ELSE
1865!
1866!--             Write result for selected indices
1867                IF ( do_calculate_perct )  THEN
1868                   perct(j, i) = perct_ij
1869                ENDIF
1870                IF ( do_calculate_utci )  THEN
1871                   utci(j, i) = utci_ij
1872                ENDIF
1873                IF ( do_calculate_pet )  THEN
1874                   pet(j, i)  = pet_ij
1875                ENDIF
1876             ENDIF
1877
1878          ENDDO
1879       ENDDO
1880    ENDIF
1881
1882 END SUBROUTINE bio_calculate_thermal_index_maps
1883
1884!--------------------------------------------------------------------------------------------------!
1885! Description:
1886! ------------
1887!> Calculate dynamic thermal indices (currently only iPT, but expandable)
1888!--------------------------------------------------------------------------------------------------!
1889 SUBROUTINE bio_calc_ipt( ta, vp, ws, pair, tmrt, dt, energy_storage, t_clo, clo, actlev, age,     &
1890                          weight, height, work, sex, ipt )
1891
1892    IMPLICIT NONE
1893!
1894!-- Input parameters
1895    INTEGER(iwp), INTENT ( IN ) ::  sex  !< Sex of agent (1 = male, 2 = female)
1896
1897    REAL(wp), INTENT ( IN )  ::  age     !< Age of agent                     (y)
1898    REAL(wp), INTENT ( IN )  ::  dt      !< Time past SINce last calculation (s)
1899    REAL(wp), INTENT ( IN )  ::  height  !< Height of agent                  (m)
1900    REAL(wp), INTENT ( IN )  ::  pair    !< Air pressure                     (hPa)
1901    REAL(wp), INTENT ( IN )  ::  ta      !< Air temperature                  (degree_C)
1902    REAL(wp), INTENT ( IN )  ::  tmrt    !< Mean radiant temperature         (degree_C)
1903    REAL(wp), INTENT ( IN )  ::  vp      !< Vapour pressure                  (hPa)
1904    REAL(wp), INTENT ( IN )  ::  weight  !< Weight of agent                  (Kg)
1905    REAL(wp), INTENT ( IN )  ::  work    !< Mechanical workload of agent  (without metabolism!) (W)
1906    REAL(wp), INTENT ( IN )  ::  ws      !< Wind speed  (local level)        (m/s)
1907
1908!
1909!-- Both, input and output parameters
1910    Real(wp), INTENT ( INOUT )  ::  actlev            !< Individuals activity level
1911                                                      !< per unit surface area      (W/m²)
1912    Real(wp), INTENT ( INOUT )  ::  clo               !< Current clothing in sulation
1913    Real(wp), INTENT ( INOUT )  ::  energy_storage    !< Energy storage   (W/m²)
1914    Real(wp), INTENT ( INOUT )  ::  t_clo             !< Clothing temperature       (degree_C)
1915!
1916!-- Output parameters
1917    REAL(wp), INTENT ( OUT ) ::  ipt    !< Instationary perceived temp.   (degree_C)
1918!
1919!-- Return immediatelly if nothing to do!
1920    IF ( .NOT. thermal_comfort )  THEN
1921        RETURN
1922    ENDIF
1923!
1924!-- If clo equals the initial value, this is the initial call
1925    IF ( clo <= -998.0_wp )  THEN
1926!
1927!--    Initialize instationary perceived temperature with personalized PT as an initial guess, set
1928!--    actlev and clo
1929       CALL ipt_init( age, weight, height, sex, work, actlev, clo, ta, vp, ws, tmrt, pair, dt,     &
1930                      energy_storage, t_clo, ipt )
1931    ELSE
1932!
1933!--    Estimate local instatinoary perceived temperature
1934       CALL ipt_cycle ( ta, vp, ws, tmrt, pair, dt, energy_storage, t_clo, clo, actlev, work, ipt )
1935    ENDIF
1936
1937 END SUBROUTINE bio_calc_ipt
1938
1939
1940
1941!--------------------------------------------------------------------------------------------------!
1942! Description:
1943! ------------
1944!> SUBROUTINE for calculating UTCI Temperature (UTCI)
1945!> computed by a 6th order approximation
1946!>
1947!> UTCI regression equation according to
1948!> Bröde P, Fiala D, Blazejczyk K, Holmér I, Jendritzky G, Kampmann B, Tinz B, Havenith G (2012)
1949!> Deriving the operational procedure for the Universal Thermal Climate Index (UTCI). International
1950!> Journal of Biometeorology 56 (3):481-494. doi:10.1007/s00484-011-0454-1
1951!>
1952!> original source available at:
1953!> www.utci.org
1954!--------------------------------------------------------------------------------------------------!
1955 SUBROUTINE calculate_utci_static( ta_in, vp, ws_hag, tmrt, hag, utci_ij )
1956
1957    IMPLICIT NONE
1958!
1959!-- Type of input of the argument list
1960    REAL(WP), INTENT ( IN )  ::  hag      !< Height of wind speed input (m)
1961    REAL(WP), INTENT ( IN )  ::  ta_in    !< Local air temperature (degree_C)
1962    REAL(WP), INTENT ( IN )  ::  tmrt     !< Local mean radiant temperature (degree_C)
1963    REAL(WP), INTENT ( IN )  ::  vp       !< Loacl vapour pressure (hPa)
1964    REAL(WP), INTENT ( IN )  ::  ws_hag   !< Incident wind speed (m/s)
1965!
1966!-- Type of output of the argument list
1967    REAL(WP) ::  d_tmrt       !< delta-tmrt               (degree_C)
1968    REAL(WP) ::  d_tmrt2      !< 2 times d_tmrt
1969    REAL(WP) ::  d_tmrt3      !< 3 times d_tmrt
1970    REAL(WP) ::  d_tmrt4      !< 4 times d_tmrt
1971    REAL(WP) ::  d_tmrt5      !< 5 times d_tmrt
1972    REAL(WP) ::  d_tmrt6      !< 6 times d_tmrt
1973    REAL(WP) ::  offset       !< utci deviation by ta cond. exceeded      (degree_C)
1974    REAL(WP) ::  pa           !< air pressure in kPa      (kPa)
1975    REAL(WP) ::  pa2          !< 2 times pa
1976    REAL(WP) ::  pa3          !< 3 times pa
1977    REAL(WP) ::  pa4          !< 4 times pa
1978    REAL(WP) ::  pa5          !< 5 times pa
1979    REAL(WP) ::  pa6          !< 6 times pa
1980    REAL(WP) ::  part_d_tmrt  !< Mean radiant temp. related part of the reg.
1981    REAL(WP) ::  part_pa      !< Air pressure related part of the regression
1982    REAL(WP) ::  part_pa2     !< Air pressure^2 related part of the regression
1983    REAL(WP) ::  part_pa3     !< Air pressure^3 related part of the regression
1984    REAL(WP) ::  part_pa46    !< Air pressure^4-6 related part of the regression
1985    REAL(WP) ::  part_ta      !< Air temperature related part of the regression
1986    REAL(WP) ::  part_va      !< Vapour pressure related part of the regression
1987    REAL(WP) ::  ta           !< air temperature modified by offset (degree_C)
1988    REAL(WP) ::  ta2          !< 2 times ta
1989    REAL(WP) ::  ta3          !< 3 times ta
1990    REAL(WP) ::  ta4          !< 4 times ta
1991    REAL(WP) ::  ta5          !< 5 times ta
1992    REAL(WP) ::  ta6          !< 6 times ta
1993    REAL(WP) ::  va           !< wind speed at 10 m above ground level    (m/s)
1994    REAL(WP) ::  va2          !< 2 times va
1995    REAL(WP) ::  va3          !< 3 times va
1996    REAL(WP) ::  va4          !< 4 times va
1997    REAL(WP) ::  va5          !< 5 times va
1998    REAL(WP) ::  va6          !< 6 times va
1999
2000
2001    REAL(wp), INTENT ( OUT ) ::  utci_ij  !< Universal Thermal Climate Index (degree_C)
2002
2003!
2004!-- Initialize
2005    offset = 0.0_wp
2006    ta = ta_in
2007    d_tmrt = tmrt - ta_in
2008!
2009!-- Use vapour pressure in kpa
2010    pa = vp / 10.0_wp
2011!
2012!-- Wind altitude correction from hag to 10m after Broede et al. (2012), eq.3
2013!-- z(0) is set to 0.01 according to UTCI profile definition
2014    va = ws_hag *  log ( 10.0_wp / 0.01_wp ) / log ( hag / 0.01_wp )
2015!
2016!-- Check if input values in range after Broede et al. (2012)
2017    IF ( ( d_tmrt > 70.0_wp )  .OR.  ( d_tmrt < -30.0_wp )  .OR.  ( vp >= 50.0_wp ) )  THEN
2018       utci_ij = bio_fill_value
2019       RETURN
2020    ENDIF
2021!
2022!-- Apply eq. 2 in Broede et al. (2012) for ta out of bounds
2023    IF ( ta > 50.0_wp )  THEN
2024       offset = ta - 50.0_wp
2025       ta = 50.0_wp
2026    ENDIF
2027    IF ( ta < -50.0_wp )  THEN
2028       offset = ta + 50.0_wp
2029       ta = -50.0_wp
2030    ENDIF
2031!
2032!-- For routine application. For wind speeds and relative humidity values below 0.5 m/s or 5%,
2033!-- respectively, the user is advised to use the lower bounds for the calculations.
2034    IF ( va < 0.5_wp )  va = 0.5_wp
2035    IF ( va > 17.0_wp )  va = 17.0_wp
2036
2037!
2038!-- Pre-calculate multiples of input parameters to save time later
2039    ta2 = ta  * ta
2040    ta3 = ta2 * ta
2041    ta4 = ta3 * ta
2042    ta5 = ta4 * ta
2043    ta6 = ta5 * ta
2044
2045    va2 = va  * va
2046    va3 = va2 * va
2047    va4 = va3 * va
2048    va5 = va4 * va
2049    va6 = va5 * va
2050
2051    d_tmrt2 = d_tmrt  * d_tmrt
2052    d_tmrt3 = d_tmrt2 * d_tmrt
2053    d_tmrt4 = d_tmrt3 * d_tmrt
2054    d_tmrt5 = d_tmrt4 * d_tmrt
2055    d_tmrt6 = d_tmrt5 * d_tmrt
2056
2057    pa2 = pa  * pa
2058    pa3 = pa2 * pa
2059    pa4 = pa3 * pa
2060    pa5 = pa4 * pa
2061    pa6 = pa5 * pa
2062
2063!
2064!-- Pre-calculate parts of the regression equation
2065    part_ta = (  6.07562052e-01_wp )       +                                                       &
2066              ( -2.27712343e-02_wp ) * ta  +                                                       &
2067              (  8.06470249e-04_wp ) * ta2 +                                                       &
2068              ( -1.54271372e-04_wp ) * ta3 +                                                       &
2069              ( -3.24651735e-06_wp ) * ta4 +                                                       &
2070              (  7.32602852e-08_wp ) * ta5 +                                                       &
2071              (  1.35959073e-09_wp ) * ta6
2072
2073    part_va = ( -2.25836520e+00_wp ) *       va  +                                                 &
2074              (  8.80326035e-02_wp ) * ta  * va  +                                                 &
2075              (  2.16844454e-03_wp ) * ta2 * va  +                                                 &
2076              ( -1.53347087e-05_wp ) * ta3 * va  +                                                 &
2077              ( -5.72983704e-07_wp ) * ta4 * va  +                                                 &
2078              ( -2.55090145e-09_wp ) * ta5 * va  +                                                 &
2079              ( -7.51269505e-01_wp ) *       va2 +                                                 &
2080              ( -4.08350271e-03_wp ) * ta  * va2 +                                                 &
2081              ( -5.21670675e-05_wp ) * ta2 * va2 +                                                 &
2082              (  1.94544667e-06_wp ) * ta3 * va2 +                                                 &
2083              (  1.14099531e-08_wp ) * ta4 * va2 +                                                 &
2084              (  1.58137256e-01_wp ) *       va3 +                                                 &
2085              ( -6.57263143e-05_wp ) * ta  * va3 +                                                 &
2086              (  2.22697524e-07_wp ) * ta2 * va3 +                                                 &
2087              ( -4.16117031e-08_wp ) * ta3 * va3 +                                                 &
2088              ( -1.27762753e-02_wp ) *       va4 +                                                 &
2089              (  9.66891875e-06_wp ) * ta  * va4 +                                                 &
2090              (  2.52785852e-09_wp ) * ta2 * va4 +                                                 &
2091              (  4.56306672e-04_wp ) *       va5 +                                                 &
2092              ( -1.74202546e-07_wp ) * ta  * va5 +                                                 &
2093              ( -5.91491269e-06_wp ) * va6
2094
2095    part_d_tmrt = (  3.98374029e-01_wp ) *             d_tmrt  +                                   &
2096                  (  1.83945314e-04_wp ) * ta  *       d_tmrt  +                                   &
2097                  ( -1.73754510e-04_wp ) * ta2 *       d_tmrt  +                                   &
2098                  ( -7.60781159e-07_wp ) * ta3 *       d_tmrt  +                                   &
2099                  (  3.77830287e-08_wp ) * ta4 *       d_tmrt  +                                   &
2100                  (  5.43079673e-10_wp ) * ta5 *       d_tmrt  +                                   &
2101                  ( -2.00518269e-02_wp ) *       va  * d_tmrt  +                                   &
2102                  (  8.92859837e-04_wp ) * ta  * va  * d_tmrt  +                                   &
2103                  (  3.45433048e-06_wp ) * ta2 * va  * d_tmrt  +                                   &
2104                  ( -3.77925774e-07_wp ) * ta3 * va  * d_tmrt  +                                   &
2105                  ( -1.69699377e-09_wp ) * ta4 * va  * d_tmrt  +                                   &
2106                  (  1.69992415e-04_wp ) *       va2 * d_tmrt  +                                   &
2107                  ( -4.99204314e-05_wp ) * ta  * va2 * d_tmrt  +                                   &
2108                  (  2.47417178e-07_wp ) * ta2 * va2 * d_tmrt  +                                   &
2109                  (  1.07596466e-08_wp ) * ta3 * va2 * d_tmrt  +                                   &
2110                  (  8.49242932e-05_wp ) *       va3 * d_tmrt  +                                   &
2111                  (  1.35191328e-06_wp ) * ta  * va3 * d_tmrt  +                                   &
2112                  ( -6.21531254e-09_wp ) * ta2 * va3 * d_tmrt  +                                   &
2113                  ( -4.99410301e-06_wp ) * va4 *       d_tmrt  +                                   &
2114                  ( -1.89489258e-08_wp ) * ta  * va4 * d_tmrt  +                                   &
2115                  (  8.15300114e-08_wp ) *       va5 * d_tmrt  +                                   &
2116                  (  7.55043090e-04_wp ) *             d_tmrt2 +                                   &
2117                  ( -5.65095215e-05_wp ) * ta  *       d_tmrt2 +                                   &
2118                  ( -4.52166564e-07_wp ) * ta2 *       d_tmrt2 +                                   &
2119                  (  2.46688878e-08_wp ) * ta3 *       d_tmrt2 +                                   &
2120                  (  2.42674348e-10_wp ) * ta4 *       d_tmrt2 +                                   &
2121                  (  1.54547250e-04_wp ) *       va  * d_tmrt2 +                                   &
2122                  (  5.24110970e-06_wp ) * ta  * va  * d_tmrt2 +                                   &
2123                  ( -8.75874982e-08_wp ) * ta2 * va  * d_tmrt2 +                                   &
2124                  ( -1.50743064e-09_wp ) * ta3 * va  * d_tmrt2 +                                   &
2125                  ( -1.56236307e-05_wp ) *       va2 * d_tmrt2 +                                   &
2126                  ( -1.33895614e-07_wp ) * ta  * va2 * d_tmrt2 +                                   &
2127                  (  2.49709824e-09_wp ) * ta2 * va2 * d_tmrt2 +                                   &
2128                  (  6.51711721e-07_wp ) *       va3 * d_tmrt2 +                                   &
2129                  (  1.94960053e-09_wp ) * ta  * va3 * d_tmrt2 +                                   &
2130                  ( -1.00361113e-08_wp ) *       va4 * d_tmrt2 +                                   &
2131                  ( -1.21206673e-05_wp ) *             d_tmrt3 +                                   &
2132                  ( -2.18203660e-07_wp ) * ta  *       d_tmrt3 +                                   &
2133                  (  7.51269482e-09_wp ) * ta2 *       d_tmrt3 +                                   &
2134                  (  9.79063848e-11_wp ) * ta3 *       d_tmrt3 +                                   &
2135                  (  1.25006734e-06_wp ) *       va  * d_tmrt3 +                                   &
2136                  ( -1.81584736e-09_wp ) * ta  * va  * d_tmrt3 +                                   &
2137                  ( -3.52197671e-10_wp ) * ta2 * va  * d_tmrt3 +                                   &
2138                  ( -3.36514630e-08_wp ) *       va2 * d_tmrt3 +                                   &
2139                  (  1.35908359e-10_wp ) * ta  * va2 * d_tmrt3 +                                   &
2140                  (  4.17032620e-10_wp ) *       va3 * d_tmrt3 +                                   &
2141                  ( -1.30369025e-09_wp ) *             d_tmrt4 +                                   &
2142                  (  4.13908461e-10_wp ) * ta  *       d_tmrt4 +                                   &
2143                  (  9.22652254e-12_wp ) * ta2 *       d_tmrt4 +                                   &
2144                  ( -5.08220384e-09_wp ) *       va  * d_tmrt4 +                                   &
2145                  ( -2.24730961e-11_wp ) * ta  * va  * d_tmrt4 +                                   &
2146                  (  1.17139133e-10_wp ) *       va2 * d_tmrt4 +                                   &
2147                  (  6.62154879e-10_wp ) *             d_tmrt5 +                                   &
2148                  (  4.03863260e-13_wp ) * ta  *       d_tmrt5 +                                   &
2149                  (  1.95087203e-12_wp ) *       va  * d_tmrt5 +                                   &
2150                  ( -4.73602469e-12_wp ) *             d_tmrt6
2151
2152    part_pa = (  5.12733497e+00_wp ) *                       pa +                                  &
2153              ( -3.12788561e-01_wp ) * ta  *                 pa +                                  &
2154              ( -1.96701861e-02_wp ) * ta2 *                 pa +                                  &
2155              (  9.99690870e-04_wp ) * ta3 *                 pa +                                  &
2156              (  9.51738512e-06_wp ) * ta4 *                 pa +                                  &
2157              ( -4.66426341e-07_wp ) * ta5 *                 pa +                                  &
2158              (  5.48050612e-01_wp ) *       va  *           pa +                                  &
2159              ( -3.30552823e-03_wp ) * ta  * va  *           pa +                                  &
2160              ( -1.64119440e-03_wp ) * ta2 * va  *           pa +                                  &
2161              ( -5.16670694e-06_wp ) * ta3 * va  *           pa +                                  &
2162              (  9.52692432e-07_wp ) * ta4 * va  *           pa +                                  &
2163              ( -4.29223622e-02_wp ) *       va2 *           pa +                                  &
2164              (  5.00845667e-03_wp ) * ta  * va2 *           pa +                                  &
2165              (  1.00601257e-06_wp ) * ta2 * va2 *           pa +                                  &
2166              ( -1.81748644e-06_wp ) * ta3 * va2 *           pa +                                  &
2167              ( -1.25813502e-03_wp ) *       va3 *           pa +                                  &
2168              ( -1.79330391e-04_wp ) * ta  * va3 *           pa +                                  &
2169              (  2.34994441e-06_wp ) * ta2 * va3 *           pa +                                  &
2170              (  1.29735808e-04_wp ) *       va4 *           pa +                                  &
2171              (  1.29064870e-06_wp ) * ta  * va4 *           pa +                                  &
2172              ( -2.28558686e-06_wp ) *       va5 *           pa +                                  &
2173              ( -3.69476348e-02_wp ) *             d_tmrt  * pa +                                  &
2174              (  1.62325322e-03_wp ) * ta  *       d_tmrt  * pa +                                  &
2175              ( -3.14279680e-05_wp ) * ta2 *       d_tmrt  * pa +                                  &
2176              (  2.59835559e-06_wp ) * ta3 *       d_tmrt  * pa +                                  &
2177              ( -4.77136523e-08_wp ) * ta4 *       d_tmrt  * pa +                                  &
2178              (  8.64203390e-03_wp ) *       va  * d_tmrt  * pa +                                  &
2179              ( -6.87405181e-04_wp ) * ta  * va  * d_tmrt  * pa +                                  &
2180              ( -9.13863872e-06_wp ) * ta2 * va  * d_tmrt  * pa +                                  &
2181              (  5.15916806e-07_wp ) * ta3 * va  * d_tmrt  * pa +                                  &
2182              ( -3.59217476e-05_wp ) *       va2 * d_tmrt  * pa +                                  &
2183              (  3.28696511e-05_wp ) * ta  * va2 * d_tmrt  * pa +                                  &
2184              ( -7.10542454e-07_wp ) * ta2 * va2 * d_tmrt  * pa +                                  &
2185              ( -1.24382300e-05_wp ) *       va3 * d_tmrt  * pa +                                  &
2186              ( -7.38584400e-09_wp ) * ta  * va3 * d_tmrt  * pa +                                  &
2187              (  2.20609296e-07_wp ) *       va4 * d_tmrt  * pa +                                  &
2188              ( -7.32469180e-04_wp ) *             d_tmrt2 * pa +                                  &
2189              ( -1.87381964e-05_wp ) * ta  *       d_tmrt2 * pa +                                  &
2190              (  4.80925239e-06_wp ) * ta2 *       d_tmrt2 * pa +                                  &
2191              ( -8.75492040e-08_wp ) * ta3 *       d_tmrt2 * pa +                                  &
2192              (  2.77862930e-05_wp ) *       va  * d_tmrt2 * pa +                                  &
2193              ( -5.06004592e-06_wp ) * ta  * va  * d_tmrt2 * pa +                                  &
2194              (  1.14325367e-07_wp ) * ta2 * va  * d_tmrt2 * pa +                                  &
2195              (  2.53016723e-06_wp ) *       va2 * d_tmrt2 * pa +                                  &
2196              ( -1.72857035e-08_wp ) * ta  * va2 * d_tmrt2 * pa +                                  &
2197              ( -3.95079398e-08_wp ) *       va3 * d_tmrt2 * pa +                                  &
2198              ( -3.59413173e-07_wp ) *             d_tmrt3 * pa +                                  &
2199              (  7.04388046e-07_wp ) * ta  *       d_tmrt3 * pa +                                  &
2200              ( -1.89309167e-08_wp ) * ta2 *       d_tmrt3 * pa +                                  &
2201              ( -4.79768731e-07_wp ) *       va  * d_tmrt3 * pa +                                  &
2202              (  7.96079978e-09_wp ) * ta  * va  * d_tmrt3 * pa +                                  &
2203              (  1.62897058e-09_wp ) *       va2 * d_tmrt3 * pa +                                  &
2204              (  3.94367674e-08_wp ) *             d_tmrt4 * pa +                                  &
2205              ( -1.18566247e-09_wp ) * ta *        d_tmrt4 * pa +                                  &
2206              (  3.34678041e-10_wp ) *       va  * d_tmrt4 * pa +                                  &
2207              ( -1.15606447e-10_wp ) *             d_tmrt5 * pa
2208
2209    part_pa2 = ( -2.80626406e+00_wp ) *                       pa2 +                                &
2210               (  5.48712484e-01_wp ) * ta  *                 pa2 +                                &
2211               ( -3.99428410e-03_wp ) * ta2 *                 pa2 +                                &
2212               ( -9.54009191e-04_wp ) * ta3 *                 pa2 +                                &
2213               (  1.93090978e-05_wp ) * ta4 *                 pa2 +                                &
2214               ( -3.08806365e-01_wp ) *       va *            pa2 +                                &
2215               (  1.16952364e-02_wp ) * ta  * va *            pa2 +                                &
2216               (  4.95271903e-04_wp ) * ta2 * va *            pa2 +                                &
2217               ( -1.90710882e-05_wp ) * ta3 * va *            pa2 +                                &
2218               (  2.10787756e-03_wp ) *       va2 *           pa2 +                                &
2219               ( -6.98445738e-04_wp ) * ta  * va2 *           pa2 +                                &
2220               (  2.30109073e-05_wp ) * ta2 * va2 *           pa2 +                                &
2221               (  4.17856590e-04_wp ) *       va3 *           pa2 +                                &
2222               ( -1.27043871e-05_wp ) * ta  * va3 *           pa2 +                                &
2223               ( -3.04620472e-06_wp ) *       va4 *           pa2 +                                &
2224               (  5.14507424e-02_wp ) *             d_tmrt  * pa2 +                                &
2225               ( -4.32510997e-03_wp ) * ta  *       d_tmrt  * pa2 +                                &
2226               (  8.99281156e-05_wp ) * ta2 *       d_tmrt  * pa2 +                                &
2227               ( -7.14663943e-07_wp ) * ta3 *       d_tmrt  * pa2 +                                &
2228               ( -2.66016305e-04_wp ) *       va  * d_tmrt  * pa2 +                                &
2229               (  2.63789586e-04_wp ) * ta  * va  * d_tmrt  * pa2 +                                &
2230               ( -7.01199003e-06_wp ) * ta2 * va  * d_tmrt  * pa2 +                                &
2231               ( -1.06823306e-04_wp ) *       va2 * d_tmrt  * pa2 +                                &
2232               (  3.61341136e-06_wp ) * ta  * va2 * d_tmrt  * pa2 +                                &
2233               (  2.29748967e-07_wp ) *       va3 * d_tmrt  * pa2 +                                &
2234               (  3.04788893e-04_wp ) *             d_tmrt2 * pa2 +                                &
2235               ( -6.42070836e-05_wp ) * ta  *       d_tmrt2 * pa2 +                                &
2236               (  1.16257971e-06_wp ) * ta2 *       d_tmrt2 * pa2 +                                &
2237               (  7.68023384e-06_wp ) *       va  * d_tmrt2 * pa2 +                                &
2238               ( -5.47446896e-07_wp ) * ta  * va  * d_tmrt2 * pa2 +                                &
2239               ( -3.59937910e-08_wp ) *       va2 * d_tmrt2 * pa2 +                                &
2240               ( -4.36497725e-06_wp ) *             d_tmrt3 * pa2 +                                &
2241               (  1.68737969e-07_wp ) * ta  *       d_tmrt3 * pa2 +                                &
2242               (  2.67489271e-08_wp ) *       va  * d_tmrt3 * pa2 +                                &
2243               (  3.23926897e-09_wp ) *             d_tmrt4 * pa2
2244
2245    part_pa3 = ( -3.53874123e-02_wp ) *                       pa3 +                                &
2246               ( -2.21201190e-01_wp ) * ta  *                 pa3 +                                &
2247               (  1.55126038e-02_wp ) * ta2 *                 pa3 +                                &
2248               ( -2.63917279e-04_wp ) * ta3 *                 pa3 +                                &
2249               (  4.53433455e-02_wp ) *       va  *           pa3 +                                &
2250               ( -4.32943862e-03_wp ) * ta  * va  *           pa3 +                                &
2251               (  1.45389826e-04_wp ) * ta2 * va  *           pa3 +                                &
2252               (  2.17508610e-04_wp ) *       va2 *           pa3 +                                &
2253               ( -6.66724702e-05_wp ) * ta  * va2 *           pa3 +                                &
2254               (  3.33217140e-05_wp ) *       va3 *           pa3 +                                &
2255               ( -2.26921615e-03_wp ) *             d_tmrt  * pa3 +                                &
2256               (  3.80261982e-04_wp ) * ta  *       d_tmrt  * pa3 +                                &
2257               ( -5.45314314e-09_wp ) * ta2 *       d_tmrt  * pa3 +                                &
2258               ( -7.96355448e-04_wp ) *       va  * d_tmrt  * pa3 +                                &
2259               (  2.53458034e-05_wp ) * ta  * va  * d_tmrt  * pa3 +                                &
2260               ( -6.31223658e-06_wp ) *       va2 * d_tmrt  * pa3 +                                &
2261               (  3.02122035e-04_wp ) *             d_tmrt2 * pa3 +                                &
2262               ( -4.77403547e-06_wp ) * ta  *       d_tmrt2 * pa3 +                                &
2263               (  1.73825715e-06_wp ) *       va  * d_tmrt2 * pa3 +                                &
2264               ( -4.09087898e-07_wp ) *             d_tmrt3 * pa3
2265
2266    part_pa46 = (  6.14155345e-01_wp ) *                       pa4 +                               &
2267                ( -6.16755931e-02_wp ) * ta  *                 pa4 +                               &
2268                (  1.33374846e-03_wp ) * ta2 *                 pa4 +                               &
2269                (  3.55375387e-03_wp ) *       va  *           pa4 +                               &
2270                ( -5.13027851e-04_wp ) * ta  * va  *           pa4 +                               &
2271                (  1.02449757e-04_wp ) *       va2 *           pa4 +                               &
2272                ( -1.48526421e-03_wp ) *             d_tmrt  * pa4 +                               &
2273                ( -4.11469183e-05_wp ) * ta  *       d_tmrt  * pa4 +                               &
2274                ( -6.80434415e-06_wp ) *       va  * d_tmrt  * pa4 +                               &
2275                ( -9.77675906e-06_wp ) *             d_tmrt2 * pa4 +                               &
2276                (  8.82773108e-02_wp ) *                       pa5 +                               &
2277                ( -3.01859306e-03_wp ) * ta  *                 pa5 +                               &
2278                (  1.04452989e-03_wp ) *       va  *           pa5 +                               &
2279                (  2.47090539e-04_wp ) *             d_tmrt  * pa5 +                               &
2280                (  1.48348065e-03_wp ) *                       pa6
2281!
2282!-- Calculate 6th order polynomial as approximation
2283    utci_ij = ta + part_ta + part_va + part_d_tmrt + part_pa + part_pa2 + part_pa3 + part_pa46
2284!
2285!-- Consider offset in result
2286    utci_ij = utci_ij + offset
2287
2288 END SUBROUTINE calculate_utci_static
2289
2290
2291
2292
2293!--------------------------------------------------------------------------------------------------!
2294! Description:
2295! ------------
2296!> Calculate_perct_static: Estimation of perceived temperature (PT, degree_C)
2297!> Value of perct is the Perceived Temperature, degree centigrade
2298!--------------------------------------------------------------------------------------------------!
2299 SUBROUTINE calculate_perct_static( ta, vp, ws, tmrt, pair, clo, perct_ij )
2300
2301    IMPLICIT NONE
2302!
2303!-- Type of input of the argument list
2304    REAL(wp), INTENT ( IN )  :: pair !< Local barometric air pressure (hPa)
2305    REAL(wp), INTENT ( IN )  :: ta   !< Local air temperature (degC)
2306    REAL(wp), INTENT ( IN )  :: tmrt !< Local mean radiant temperature (degC)
2307    REAL(wp), INTENT ( IN )  :: vp   !< Local vapour pressure (hPa)
2308    REAL(wp), INTENT ( IN )  :: ws   !< Local wind velocitry (m/s)
2309!
2310!-- Type of output of the argument list
2311    REAL(wp), INTENT ( OUT ) :: clo       !< Clothing index (dimensionless)
2312    REAL(wp), INTENT ( OUT ) :: perct_ij  !< Perceived temperature (degC)
2313!
2314!-- Parameters for standard "Klima-Michel"
2315    REAL(wp), PARAMETER :: actlev = 134.6862_wp  !< Workload by activity per standardized surface (A_Du)
2316    REAL(wp), PARAMETER :: eta = 0.0_wp          !< Mechanical work efficiency for walking on flat
2317                                                 !< ground (compare to Fanger (1972) pp 24f)
2318!
2319!-- Type of program variables
2320    REAL(wp), PARAMETER :: eps = 0.0005  !< Accuracy in clothing insulation (clo) for evaluation the root of Fanger's PMV (pmva=0)
2321
2322    INTEGER(iwp) :: ncount      !< running index
2323    INTEGER(iwp) :: nerr_cold   !< error number (cold conditions)
2324    INTEGER(iwp) :: nerr        !< error number
2325
2326    LOGICAL :: sultrieness
2327
2328    REAL(wp) ::  clon           !< clo for neutral conditions   (clo)
2329    REAL(wp) ::  d_pmv          !< PMV deviation (dimensionless --> PMV)
2330    REAL(wp) ::  dgtcm          !< Mean deviation dependent on perct
2331    REAL(wp) ::  dgtcstd        !< Mean deviation plus its standard deviation
2332    REAL(wp) ::  d_std          !< factor to threshold for sultriness
2333    REAL(wp) ::  ireq_minimal   !< Minimal required clothing insulation (clo)
2334    REAL(wp) ::  pmv_s          !< Fangers predicted mean vote for summer clothing
2335    REAL(wp) ::  pmv_w          !< Fangers predicted mean vote for winter clothing
2336    REAL(wp) ::  pmva           !< adjusted predicted mean vote
2337    REAL(wp) ::  pmvs           !< pred. mean vote considering sultrieness
2338    REAL(wp) ::  ptc            !< perceived temp. for cold conditions (degree_C)
2339    REAL(wp) ::  sclo           !< summer clothing insulation
2340    REAL(wp) ::  svp_ta         !< saturation vapor pressure    (hPa)
2341    REAL(wp) ::  sult_lim       !< threshold for sultrieness    (hPa)
2342    REAL(wp) ::  wclo           !< winter clothing insulation
2343
2344!
2345!-- Initialise
2346    perct_ij = bio_fill_value
2347
2348    nerr     = 0_iwp
2349    ncount   = 0_iwp
2350    sultrieness  = .FALSE.
2351!
2352!-- Tresholds: clothing insulation (account for model inaccuracies)
2353!-- Summer clothing
2354    sclo     = 0.44453_wp
2355!
2356!-- Winter clothing
2357    wclo     = 1.76267_wp
2358!
2359!-- Eecision: first calculate for winter or summer clothing
2360    IF ( ta <= 10.0_wp )  THEN
2361!
2362!--    First guess: winter clothing insulation: cold stress
2363       clo = wclo
2364       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
2365       pmv_w = pmva
2366
2367       IF ( pmva > 0.0_wp )  THEN
2368!
2369!--       Case summer clothing insulation: heat load ?
2370          clo = sclo
2371          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
2372          pmv_s = pmva
2373          IF ( pmva <= 0.0_wp )  THEN
2374!
2375!--          Case: comfort achievable by varying clothing insulation between winter and summer set
2376!--                values
2377             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo, pmv_s, wclo, pmv_w, eps, &
2378                               pmva, ncount, clo )
2379             IF ( ncount < 0_iwp )  THEN
2380                nerr = -1_iwp
2381                RETURN
2382             ENDIF
2383          ELSE IF ( pmva > 0.06_wp )  THEN
2384             clo = 0.5_wp
2385             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,  pmva )
2386          ENDIF
2387       ELSE IF ( pmva < - 0.11_wp )  THEN
2388          clo = 1.75_wp
2389          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
2390       ENDIF
2391    ELSE
2392!
2393!--    First guess: summer clothing insulation: heat load
2394       clo = sclo
2395       CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
2396       pmv_s = pmva
2397
2398       IF ( pmva < 0.0_wp )  THEN
2399!
2400!--       Case winter clothing insulation: cold stress ?
2401          clo = wclo
2402          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
2403          pmv_w = pmva
2404
2405          IF ( pmva >= 0.0_wp )  THEN
2406!
2407!--          Case: comfort achievable by varying clothing insulation between winter and summer set
2408!--                values
2409             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo, pmv_s, wclo, pmv_w, eps, &
2410                               pmva, ncount, clo )
2411             IF ( ncount < 0_iwp )  THEN
2412                nerr = -1_iwp
2413                RETURN
2414             ENDIF
2415          ELSE IF ( pmva < - 0.11_wp )  THEN
2416             clo = 1.75_wp
2417             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
2418          ENDIF
2419       ELSE IF ( pmva > 0.06_wp )  THEN
2420          clo = 0.5_wp
2421          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva )
2422       ENDIF
2423
2424    ENDIF
2425!
2426!-- Determine perceived temperature by regression equation + adjustments
2427    pmvs = pmva
2428    CALL perct_regression( pmva, clo, perct_ij )
2429    ptc = perct_ij
2430    IF ( clo >= 1.75_wp  .AND.  pmva <= - 0.11_wp )  THEN
2431!
2432!--    Adjust for cold conditions according to Gagge 1986
2433       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
2434       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
2435       pmvs = pmva - d_pmv
2436       IF ( pmvs > - 0.11_wp )  THEN
2437          d_pmv  = 0.0_wp
2438          pmvs   = - 0.11_wp
2439       ENDIF
2440       CALL perct_regression( pmvs, clo, perct_ij )
2441    ENDIF
2442!     clo_fanger = clo
2443    clon = clo
2444    IF ( clo > 0.5_wp  .AND.  perct_ij <= 8.73_wp )  THEN
2445!
2446!--    Required clothing insulation (ireq) is exclusively defined for perceived temperatures (perct)
2447!--    less 10 (C) for a reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s.
2448       clon = ireq_neutral ( perct_ij, ireq_minimal, nerr )
2449       clo = clon
2450    ENDIF
2451    CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
2452    sultrieness    = .FALSE.
2453    d_std = -99.0_wp
2454    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
2455!
2456!--    Adjust for warm/humid conditions according to Gagge 1986
2457       CALL saturation_vapor_pressure ( ta, svp_ta )
2458       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
2459       pmvs   = pmva + d_pmv
2460       CALL perct_regression( pmvs, clo, perct_ij )
2461       IF ( sult_lim < 99.0_wp )  THEN
2462          IF ( (perct_ij - ptc) > sult_lim )  sultrieness = .TRUE.
2463!
2464!--       Set factor to threshold for sultriness
2465          IF ( ABS( dgtcstd ) > 0.00001_wp )  THEN
2466             d_std = ( ( perct_ij - ptc ) - dgtcm ) / dgtcstd
2467          ENDIF
2468       ENDIF
2469    ENDIF
2470
2471 END SUBROUTINE calculate_perct_static
2472
2473!--------------------------------------------------------------------------------------------------!
2474! Description:
2475! ------------
2476!> The SUBROUTINE calculates the (saturation) water vapour pressure (hPa = hecto Pascal) for a given
2477!> temperature ta (degC).
2478!>'ta' can be the air temperature or the dew point temperature. The first will result in the current
2479!> vapor pressure (hPa), the latter will calulate the saturation vapor pressure (hPa).
2480!--------------------------------------------------------------------------------------------------!
2481 SUBROUTINE saturation_vapor_pressure( ta, svp_ta )
2482
2483    IMPLICIT NONE
2484
2485    REAL(wp), INTENT ( IN )  ::  ta     !< ambient air temperature (degC)
2486    REAL(wp), INTENT ( OUT ) ::  svp_ta !< water vapour pressure (hPa)
2487
2488    REAL(wp)      ::  b
2489    REAL(wp)      ::  c
2490
2491
2492    IF ( ta < 0.0_wp )  THEN
2493!
2494!--    ta  < 0 (degC): water vapour pressure over ice
2495       b = 17.84362_wp
2496       c = 245.425_wp
2497    ELSE
2498!
2499!--    ta >= 0 (degC): water vapour pressure over water
2500       b = 17.08085_wp
2501       c = 234.175_wp
2502    ENDIF
2503!
2504!-- Saturation water vapour pressure
2505    svp_ta = 6.1078_wp * EXP( b * ta / ( c + ta ) )
2506
2507 END SUBROUTINE saturation_vapor_pressure
2508
2509!--------------------------------------------------------------------------------------------------!
2510! Description:
2511! ------------
2512!> Find the clothing insulation value clo_res (clo) to make Fanger's Predicted Mean Vote (PMV) equal
2513!> comfort (pmva=0) for actual meteorological conditions (ta,tmrt, vp, ws, pair) and values of
2514!> individual's activity level.
2515!--------------------------------------------------------------------------------------------------!
2516 SUBROUTINE iso_ridder( ta, tmrt, vp, ws, pair, actlev, eta, sclo, pmv_s, wclo, pmv_w, eps, pmva,  &
2517                        nerr, clo_res )
2518
2519    IMPLICIT NONE
2520!
2521!-- Input variables of argument list:
2522    REAL(wp), INTENT ( IN )  :: actlev   !< Individuals activity level per unit surface area (W/m2)
2523    REAL(wp), INTENT ( IN )  :: eps      !< (0.05) accuracy in clothing insulation (clo) for evaluation the root of Fanger's PMV (pmva=0)
2524    REAL(wp), INTENT ( IN )  :: eta      !< Individuals work efficiency (dimensionless)
2525    REAL(wp), INTENT ( IN )  :: pair     !< Barometric air pressure (hPa)
2526    REAL(wp), INTENT ( IN )  :: pmv_s    !< Fanger's PMV corresponding to sclo
2527    REAL(wp), INTENT ( IN )  :: pmv_w    !< Fanger's PMV corresponding to wclo
2528    REAL(wp), INTENT ( IN )  :: sclo     !< Lower threshold of bracketing clothing insulation (clo)
2529    REAL(wp), INTENT ( IN )  :: ta       !< Ambient temperature (degC)
2530    REAL(wp), INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC)
2531    REAL(wp), INTENT ( IN )  :: vp       !< Water vapour pressure (hPa)
2532    REAL(wp), INTENT ( IN )  :: wclo     !< Upper threshold of bracketing clothing insulation (clo)
2533    REAL(wp), INTENT ( IN )  :: ws       !< Wind speed (m/s) 1 m above ground
2534!
2535!-- Output variables of argument list:
2536    INTEGER(iwp), INTENT ( OUT ) :: nerr !< Error status / quality flag
2537                                         !< nerr >= 0, o.k., and nerr is the number of iterations for convergence
2538                                         !< nerr = -1: error = malfunction of Ridder's convergence method
2539                                         !< nerr = -2: error = maximum iterations (max_iteration) exceeded
2540                                         !< nerr = -3: error = root not bracketed between sclo and wclo
2541
2542    REAL(wp), INTENT ( OUT ) :: clo_res  !< Resulting clothing insulation value (clo)
2543    REAL(wp), INTENT ( OUT ) :: pmva     !< 0 (set to zero, because clo is evaluated for comfort)
2544!
2545!-- Type of program variables
2546    INTEGER(iwp), PARAMETER  ::  max_iteration = 15_iwp       !< max number of iterations
2547
2548    REAL(wp),     PARAMETER  ::  guess_0       = -1.11e30_wp  !< initial guess
2549
2550    INTEGER(iwp) ::  j       !< running index
2551
2552    REAL(wp) ::  clo_lower   !< lower limit of clothing insulation      (clo)
2553    REAL(wp) ::  clo_upper   !< upper limit of clothing insulation      (clo)
2554    REAL(wp) ::  sroot       !< sqrt of PMV-guess
2555    REAL(wp) ::  x_average   !< average of x_lower and x_upper          (clo)
2556    REAL(wp) ::  x_lower     !< lower guess for clothing insulation     (clo)
2557    REAL(wp) ::  x_new       !< preliminary result for clothing insulation (clo)
2558    REAL(wp) ::  x_ridder    !< current guess for clothing insulation   (clo)
2559    REAL(wp) ::  x_upper     !< upper guess for clothing insulation     (clo)
2560    REAL(wp) ::  y_average   !< average of y_lower and y_upper
2561    REAL(wp) ::  y_new       !< preliminary result for pred. mean vote
2562    REAL(wp) ::  y_lower     !< predicted mean vote for summer clothing
2563    REAL(wp) ::  y_upper     !< predicted mean vote for winter clothing
2564!
2565!-- Initialise
2566    nerr    = 0_iwp
2567!
2568!-- Set pmva = 0 (comfort): Root of PMV depending on clothing insulation
2569    x_ridder    = bio_fill_value
2570    pmva        = 0.0_wp
2571    clo_lower   = sclo
2572    y_lower     = pmv_s
2573    clo_upper   = wclo
2574    y_upper     = pmv_w
2575    IF ( ( y_lower > 0.0_wp .AND. y_upper < 0.0_wp )  .OR.                                         &
2576         ( y_lower < 0.0_wp .AND. y_upper > 0.0_wp ) )  THEN
2577       x_lower  = clo_lower
2578       x_upper  = clo_upper
2579       x_ridder = guess_0
2580
2581       DO  j = 1_iwp, max_iteration
2582          x_average = 0.5_wp * ( x_lower + x_upper )
2583          CALL fanger ( ta, tmrt, vp, ws, pair, x_average, actlev, eta, y_average )
2584          sroot = SQRT( y_average**2 - y_lower * y_upper )
2585          IF ( ABS( sroot ) < 0.00001_wp )  THEN
2586             clo_res = x_average
2587             nerr = j
2588             RETURN
2589          ENDIF
2590          x_new = x_average + ( x_average - x_lower ) *                                            &
2591                  ( SIGN ( 1.0_wp, y_lower - y_upper ) * y_average / sroot )
2592          IF ( ABS( x_new - x_ridder ) <= eps )  THEN
2593             clo_res = x_ridder
2594             nerr       = j
2595             RETURN
2596          ENDIF
2597          x_ridder = x_new
2598          CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta, y_new )
2599          IF ( ABS( y_new ) < 0.00001_wp )  THEN
2600             clo_res = x_ridder
2601             nerr       = j
2602             RETURN
2603          ENDIF
2604          IF ( ABS( SIGN( y_average, y_new ) - y_average ) > 0.00001_wp )  THEN
2605             x_lower = x_average
2606             y_lower = y_average
2607             x_upper  = x_ridder
2608             y_upper  = y_new
2609          ELSE IF ( ABS( SIGN( y_lower, y_new ) - y_lower ) > 0.00001_wp )  THEN
2610             x_upper  = x_ridder
2611             y_upper  = y_new
2612          ELSE IF ( ABS( SIGN( y_upper, y_new ) - y_upper ) > 0.00001_wp )  THEN
2613             x_lower = x_ridder
2614             y_lower = y_new
2615          ELSE
2616!
2617!--          Never get here in x_ridder: SINgularity in y
2618             nerr    = -1_iwp
2619             clo_res = x_ridder
2620             RETURN
2621          ENDIF
2622          IF ( ABS( x_upper - x_lower ) <= eps )  THEN
2623             clo_res = x_ridder
2624             nerr    = j
2625             RETURN
2626          ENDIF
2627       ENDDO
2628!
2629!--    x_ridder exceed maximum iterations
2630       nerr       = -2_iwp
2631       clo_res = y_new
2632       RETURN
2633    ELSE IF ( ABS( y_lower ) < 0.00001_wp )  THEN
2634       x_ridder = clo_lower
2635    ELSE IF ( ABS( y_upper ) < 0.00001_wp )  THEN
2636       x_ridder = clo_upper
2637    ELSE
2638!
2639!--    x_ridder not bracketed by u_clo and o_clo
2640       nerr = -3_iwp
2641       clo_res = x_ridder
2642       RETURN
2643    ENDIF
2644
2645 END SUBROUTINE iso_ridder
2646
2647!--------------------------------------------------------------------------------------------------!
2648! Description:
2649! ------------
2650!> Regression relations between perceived temperature (perct) and (adjusted) PMV. The regression
2651!> presumes the Klima-Michel settings for reference individual and reference environment.
2652!--------------------------------------------------------------------------------------------------!
2653 SUBROUTINE perct_regression( pmv, clo, perct_ij )
2654
2655    IMPLICIT NONE
2656
2657    REAL(wp), INTENT ( IN ) ::  clo   !< clothing insulation index (clo)
2658    REAL(wp), INTENT ( IN ) ::  pmv   !< Fangers predicted mean vote (dimensionless)
2659
2660    REAL(wp), INTENT ( OUT ) ::  perct_ij   !< perct (degC) corresponding to given PMV / clo
2661
2662    IF ( pmv <= - 0.11_wp )  THEN
2663       perct_ij = 5.805_wp + 12.6784_wp * pmv
2664    ELSE
2665       IF ( pmv >= + 0.01_wp )  THEN
2666          perct_ij = 16.826_wp + 6.163_wp * pmv
2667       ELSE
2668          perct_ij = 21.258_wp - 9.558_wp * clo
2669       ENDIF
2670    ENDIF
2671
2672 END SUBROUTINE perct_regression
2673
2674!--------------------------------------------------------------------------------------------------!
2675! Description:
2676! ------------
2677!> FANGER.F90
2678!>
2679!> SI-VERSION: ACTLEV W m-2, VAPOUR PRESSURE hPa
2680!> Calculates the current Predicted Mean Vote according to Fanger.
2681!> The case of free convection (ws < 0.1 m/s) is dealt with ws = 0.1 m/s
2682!--------------------------------------------------------------------------------------------------!
2683 SUBROUTINE fanger( ta, tmrt, pa, in_ws, pair, in_clo, actlev, eta, pmva )
2684
2685    IMPLICIT NONE
2686!
2687!-- Input variables of argument list:
2688    REAL(wp), INTENT ( IN ) ::  actlev   !< Individuals activity level per unit surface area (W/m2)
2689    REAL(wp), INTENT ( IN ) ::  eta      !< Individuals mechanical work efficiency (dimensionless)
2690    REAL(wp), INTENT ( IN ) ::  in_clo   !< Clothing insulation (clo)
2691    REAL(wp), INTENT ( IN ) ::  in_ws    !< Wind speed (m/s) 1 m above ground
2692    REAL(wp), INTENT ( IN ) ::  pa       !< Water vapour pressure (hPa)
2693    REAL(wp), INTENT ( IN ) ::  pair     !< Barometric pressure (hPa) at site
2694    REAL(wp), INTENT ( IN ) ::  ta       !< Ambient air temperature (degC)
2695    REAL(wp), INTENT ( IN ) ::  tmrt     !< Mean radiant temperature (degC)
2696
2697!
2698!-- Output variables of argument list:
2699    REAL(wp), INTENT ( OUT ) ::  pmva    !< Actual Predicted Mean Vote (PMV,
2700                                         !< dimensionless) according to Fanger corresponding to meteorological
2701                                         !< (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta)
2702!
2703!-- Internal variables
2704    INTEGER(iwp) :: i         !< running index
2705
2706    REAL(wp) ::  activity     !< persons activity  (must stay == actlev, W)
2707    REAL(wp) ::  bc           !< preliminary result storage
2708    REAL(wp) ::  cc           !< preliminary result storage
2709    REAL(wp) ::  clo          !< clothing insulation index              (clo)
2710    REAL(wp) ::  dc           !< preliminary result storage
2711    REAL(wp) ::  ec           !< preliminary result storage
2712    REAL(wp) ::  f_cl         !< Increase in surface due to clothing    (factor)
2713    REAL(wp) ::  gc           !< preliminary result storage
2714    REAL(wp) ::  heat_convection  !< energy loss by autocnvection       (W)
2715    REAL(wp) ::  hr           !< radiational heat resistence
2716    REAL(wp) ::  t_clothing   !< clothing temperature                   (degree_C)
2717    REAL(wp) ::  t_skin_aver  !< average skin temperature               (degree_C)
2718    REAL(wp) ::  ws           !< wind speed                             (m/s)
2719    REAL(wp) ::  z1           !< Empiric factor for the adaption of the heat
2720                              !< ballance equation to the psycho-physical scale (Equ. 40 in FANGER)
2721    REAL(wp) ::  z2           !< Water vapour diffution through the skin
2722    REAL(wp) ::  z3           !< Sweat evaporation from the skin surface
2723    REAL(wp) ::  z4           !< Loss of latent heat through respiration
2724    REAL(wp) ::  z5           !< Loss of radiational heat
2725    REAL(wp) ::  z6           !< Heat loss through forced convection
2726
2727!
2728!-- Clo must be > 0. to avoid div. by 0!
2729    clo = in_clo
2730    IF ( clo <= 0.0_wp )  clo = .001_wp
2731!
2732!-- f_cl = increase in surface due to clothing
2733    f_cl = 1.0_wp + 0.15_wp * clo
2734!
2735!-- Case of free convection (ws < 0.1 m/s ) not considered
2736    ws = in_ws
2737    IF ( ws < 0.1_wp )  THEN
2738       ws = 0.1_wp
2739    ENDIF
2740!
2741!-- Heat_convection = forced convection
2742    heat_convection = 12.1_wp * SQRT( ws * pair / 1013.25_wp )
2743!
2744!-- Activity = inner heat production per standardized surface
2745    activity = actlev * ( 1.0_wp - eta )
2746!
2747!-- t_skin_aver = average skin temperature
2748    t_skin_aver = 35.7_wp - 0.0275_wp * activity
2749!
2750!-- Calculation of constants for evaluation below
2751    bc = 0.155_wp * clo * 3.96_wp * 10.0_wp**( -8 ) * f_cl
2752    cc = f_cl * heat_convection
2753    ec = 0.155_wp * clo
2754    dc = ( 1.0_wp + ec * cc ) / bc
2755    gc = ( t_skin_aver + bc * ( tmrt + degc_to_k )**4 + ec * cc * ta ) / bc
2756!
2757!-- Calculation of clothing surface temperature (t_clothing) based on Newton-approximation with air
2758!-- temperature as initial guess.
2759    t_clothing = ta
2760    DO  i = 1, 3
2761       t_clothing = t_clothing - ( ( t_clothing + degc_to_k )**4 + t_clothing * dc - gc ) /        &
2762                    ( 4.0_wp * ( t_clothing + degc_to_k )**3 + dc )
2763    ENDDO
2764!
2765!-- Empiric factor for the adaption of the heat ballance equation to the psycho-physical scale (Equ.
2766!-- 40 in FANGER)
2767    z1 = ( 0.303_wp * EXP( - 0.036_wp * actlev ) + 0.0275_wp )
2768!
2769!-- Water vapour diffution through the skin
2770    z2 = 0.31_wp * ( 57.3_wp - 0.07_wp * activity-pa )
2771!
2772!-- Sweat evaporation from the skin surface
2773    z3 = 0.42_wp * ( activity - 58.0_wp )
2774!
2775!-- Loss of latent heat through respiration
2776    z4 = 0.0017_wp * actlev * ( 58.7_wp - pa ) + 0.0014_wp * actlev *                              &
2777         ( 34.0_wp - ta )
2778!
2779!-- Loss of radiational heat
2780    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + degc_to_k )**4 - ( tmrt + degc_to_k )**4 )
2781    IF ( ABS( t_clothing - tmrt ) > 0.0_wp )  THEN
2782       hr = z5 / f_cl / ( t_clothing - tmrt )
2783    ELSE
2784       hr = 0.0_wp
2785    ENDIF
2786!
2787!-- Heat loss through forced convection cc*(t_clothing-TT)
2788    z6 = cc * ( t_clothing - ta )
2789!
2790!-- Predicted Mean Vote
2791    pmva = z1 * ( activity - z2 - z3 - z4 - z5 - z6 )
2792
2793 END SUBROUTINE fanger
2794
2795!--------------------------------------------------------------------------------------------------!
2796! Description:
2797! ------------
2798!> For pmva > 0 and clo =0.5 the increment (deltapmv) is calculated that converts pmva into Gagge's
2799!> et al. (1986) PMV*.
2800!--------------------------------------------------------------------------------------------------!
2801 REAL(wp) FUNCTION deltapmv( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
2802
2803    IMPLICIT NONE
2804
2805!
2806!-- Input variables of argument list:
2807    REAL(wp),     INTENT ( IN )  :: pmva     !< Actual Predicted Mean Vote (PMV) according to Fanger
2808    REAL(wp),     INTENT ( IN )  :: svp_ta   !< Saturation water vapour pressure (hPa) at ta
2809    REAL(wp),     INTENT ( IN )  :: ta       !< Ambient temperature (degC) at screen level
2810    REAL(wp),     INTENT ( IN )  :: tmrt     !< Mean radiant temperature (degC) at screen level
2811    REAL(wp),     INTENT ( IN )  :: vp       !< Water vapour pressure (hPa) at screen level
2812    REAL(wp),     INTENT ( IN )  :: ws       !< Wind speed (m/s) 1 m above ground
2813
2814!
2815!-- Output variables of argument list:
2816    INTEGER(iwp), INTENT ( OUT ) :: nerr     !< Error status / quality flag
2817                                             !<  0 = o.k.
2818                                             !< -2 = pmva outside valid regression range
2819                                             !< -3 = rel. humidity set to 5 % or 95 %, respectively
2820                                             !< -4 = deltapmv set to avoid pmvs < 0
2821
2822!
2823!-- Internal variables:
2824    INTEGER(iwp) :: nreg      !<
2825
2826    REAL(wp) ::  apa          !< natural logarithm of pa (with hard lower border)
2827    REAL(wp) ::  dapa         !< difference of apa and pa_p50
2828    REAL(wp) ::  dpmv_1       !<
2829    REAL(wp) ::  dpmv_2       !<
2830    REAL(wp) ::  dtmrt        !< difference mean radiation to air temperature
2831    REAL(wp) ::  pa           !< vapor pressure (hPa) with hard bounds
2832    REAL(wp) ::  pa_p50       !< ratio actual water vapour pressure to that of relative humidity of
2833                              !< 50 %
2834    REAL(wp) ::  pmv          !< temp storage og predicted mean vote
2835    REAL(wp) ::  pmvs         !<
2836    REAL(wp) ::  p10          !< lower bound for pa
2837    REAL(wp) ::  p95          !< upper bound for pa
2838    REAL(wp) ::  sqvel        !< square root of local wind velocity
2839    REAL(wp) ::  weight       !<
2840    REAL(wp) ::  weight2      !<
2841
2842!
2843!-- Regression coefficients:
2844    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpmv = (/                                              &
2845     - 0.0556602_wp, - 0.1528680_wp, - 0.2336104_wp, - 0.2789387_wp,                               &
2846     - 0.3551048_wp, - 0.4304076_wp, - 0.4884961_wp, - 0.4897495_wp /)
2847
2848    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpa_p50 = (/                                           &
2849     - 0.1607154_wp, - 0.4177296_wp, - 0.4120541_wp, - 0.0886564_wp,                               &
2850       0.4285938_wp,   0.6281256_wp,   0.5067361_wp,   0.3965169_wp /)
2851
2852    REAL(wp), DIMENSION(0:7), PARAMETER ::  bpa = (/                                               &
2853       0.0580284_wp,   0.0836264_wp,   0.1009919_wp,   0.1020777_wp,                               &
2854       0.0898681_wp,   0.0839116_wp,   0.0853258_wp,   0.0866589_wp /)
2855
2856    REAL(wp), DIMENSION(0:7), PARAMETER ::  bapa = (/                                              &
2857     - 1.7838788_wp, - 2.9306231_wp, - 1.6350334_wp,    0.6211547_wp,                              &
2858       3.3918083_wp,   5.5521025_wp,   8.4897418_wp,   16.6265851_wp /)
2859
2860    REAL(wp), DIMENSION(0:7), PARAMETER ::  bdapa = (/                                             &
2861       1.6752720_wp,   2.7379504_wp,   1.2940526_wp, -  1.0985759_wp,                              &
2862     - 3.9054732_wp, - 6.0403012_wp, - 8.9437119_wp, - 17.0671201_wp /)
2863
2864    REAL(wp), DIMENSION(0:7), PARAMETER ::  bsqvel = (/                                            &
2865     - 0.0315598_wp, - 0.0286272_wp, - 0.0009228_wp,   0.0483344_wp,                               &
2866       0.0992366_wp,   0.1491379_wp,   0.1951452_wp,   0.2133949_wp /)
2867
2868    REAL(wp), DIMENSION(0:7), PARAMETER ::  bta = (/                                               &
2869       0.0953986_wp,   0.1524760_wp,   0.0564241_wp, - 0.0893253_wp,                               &
2870     - 0.2398868_wp, - 0.3515237_wp, - 0.5095144_wp, - 0.9469258_wp /)
2871
2872    REAL(wp), DIMENSION(0:7), PARAMETER ::  bdtmrt = (/                                            &
2873     - 0.0004672_wp, - 0.0000514_wp, - 0.0018037_wp, - 0.0049440_wp,                               &
2874     - 0.0069036_wp, - 0.0075844_wp, - 0.0079602_wp, - 0.0089439_wp /)
2875
2876    REAL(wp), DIMENSION(0:7), PARAMETER ::  aconst = (/                                            &
2877       1.8686215_wp,   3.4260713_wp,    2.0116185_wp, -  0.7777552_wp,                             &
2878     - 4.6715853_wp, - 7.7314281_wp, - 11.7602578_wp, - 23.5934198_wp /)
2879
2880
2881!
2882!-- Test for compliance with regression range
2883    IF ( pmva < -1.0_wp  .OR.  pmva > 7.0_wp )  THEN
2884       nerr = -2_iwp
2885    ELSE
2886       nerr = 0_iwp
2887    ENDIF
2888!
2889!-- Initialise classic PMV
2890    pmv  = pmva
2891!
2892!-- Water vapour pressure of air
2893    p10  = 0.05_wp * svp_ta
2894    p95  = 1.00_wp * svp_ta
2895    IF ( vp >= p10  .AND.  vp <= p95 )  THEN
2896       pa = vp
2897    ELSE
2898       nerr = -3_iwp
2899       IF ( vp < p10 )  THEN
2900!
2901!--       Due to conditions of regression: r.H. >= 5 %
2902          pa = p10
2903       ELSE
2904!
2905!--       Due to conditions of regression: r.H. <= 95 %
2906          pa = p95
2907       ENDIF
2908    ENDIF
2909    IF ( pa > 0.0_wp )  THEN
2910!
2911!--    Natural logarithm of pa
2912       apa = LOG( pa )
2913    ELSE
2914       apa = -5.0_wp
2915    ENDIF
2916!
2917!-- Ratio actual water vapour pressure to that of a r.H. of 50 %
2918    pa_p50   = 0.5_wp * svp_ta
2919    IF ( pa_p50 > 0.0_wp  .AND.  pa > 0.0_wp )  THEN
2920       dapa   = apa - LOG( pa_p50 )
2921       pa_p50 = pa / pa_p50
2922    ELSE
2923       dapa   = -5.0_wp
2924       pa_p50 = 0.0_wp
2925    ENDIF
2926!
2927!-- Square root of wind velocity
2928    IF ( ws >= 0.0_wp )  THEN
2929       sqvel = SQRT( ws )
2930    ELSE
2931       sqvel = 0.0_wp
2932    ENDIF
2933!
2934!-- Difference mean radiation to air temperature
2935    dtmrt = tmrt - ta
2936!
2937!-- Select the valid regression coefficients
2938    nreg = INT( pmv )
2939    IF ( nreg < 0_iwp )  THEN
2940!
2941!--    Value of the FUNCTION in the case pmv <= -1
2942       deltapmv = 0.0_wp
2943       RETURN
2944    ENDIF
2945    weight = MOD ( pmv, 1.0_wp )
2946    IF ( weight < 0.0_wp )  weight = 0.0_wp
2947    IF ( nreg > 5_iwp )  THEN
2948       nreg  = 5_iwp
2949       weight   = pmv - 5.0_wp
2950       weight2  = pmv - 6.0_wp
2951       IF ( weight2 > 0_iwp )  THEN
2952          weight = ( weight - weight2 ) / weight
2953       ENDIF
2954    ENDIF
2955!
2956!-- Regression valid for 0. <= pmv <= 6., bounds are checked above
2957    dpmv_1 =                                                                                       &
2958             + bpa(nreg)     * pa                                                                  &
2959             + bpmv(nreg)    * pmv                                                                 &
2960             + bapa(nreg)    * apa                                                                 &
2961             + bta(nreg)     * ta                                                                  &
2962             + bdtmrt(nreg)  * dtmrt                                                               &
2963             + bdapa(nreg)   * dapa                                                                &
2964             + bsqvel(nreg)  * sqvel                                                               &
2965             + bpa_p50(nreg) * pa_p50                                                              &
2966             + aconst(nreg)
2967
2968!    dpmv_2 = 0.0_wp
2969!    IF ( nreg < 6_iwp )  THEN  !< nreg is always <= 5, see above
2970    dpmv_2 =                                                                                       &
2971             + bpa(nreg+1_iwp)     * pa                                                            &
2972             + bpmv(nreg+1_iwp)    * pmv                                                           &
2973             + bapa(nreg+1_iwp)    * apa                                                           &
2974             + bta(nreg+1_iwp)     * ta                                                            &
2975             + bdtmrt(nreg+1_iwp)  * dtmrt                                                         &
2976             + bdapa(nreg+1_iwp)   * dapa                                                          &
2977             + bsqvel(nreg+1_iwp)  * sqvel                                                         &
2978             + bpa_p50(nreg+1_iwp) * pa_p50                                                        &
2979             + aconst(nreg+1_iwp)
2980!    ENDIF
2981!
2982!-- Calculate pmv modification
2983    deltapmv = ( 1.0_wp - weight ) * dpmv_1 + weight * dpmv_2
2984    pmvs = pmva + deltapmv
2985    IF ( ( pmvs ) < 0.0_wp )  THEN
2986!
2987!--    Prevent negative pmv* due to problems with clothing insulation
2988       nerr = -4_iwp
2989       IF ( pmvs > -0.11_wp )  THEN
2990!
2991!--       Threshold from perct_regression for winter clothing insulation
2992          deltapmv = deltapmv + 0.11_wp
2993       ELSE
2994!
2995!--       Set pmvs to "0" for compliance with summer clothing insulation
2996          deltapmv = -1.0_wp * pmva
2997       ENDIF
2998    ENDIF
2999
3000 END FUNCTION deltapmv
3001
3002!--------------------------------------------------------------------------------------------------!
3003! Description:
3004! ------------
3005!> The subroutine "calc_sultr" returns a threshold value to perceived temperature allowing to decide
3006!> whether the actual perceived temperature is linked to perecption of sultriness. The threshold
3007!> values depends on the Fanger's classical PMV, expressed here as perceived temperature perct.
3008!--------------------------------------------------------------------------------------------------!
3009 SUBROUTINE calc_sultr( perct_ij, dperctm, dperctstd, sultr_res )
3010
3011    IMPLICIT NONE
3012!
3013!-- Input of the argument list:
3014    REAL(wp), INTENT ( IN )  ::  perct_ij   !< Classical perceived temperature: Base is Fanger's PMV
3015!
3016!-- Additional output variables of argument list:
3017    REAL(wp), INTENT ( OUT ) ::  dperctm    !< Mean deviation perct (classical gt) to gt* (rational
3018                                            !< gt calculated based on Gagge's rational PMV*)
3019    REAL(wp), INTENT ( OUT ) ::  dperctstd  !< dperctm plus its standard deviation times a factor
3020                                            !< determining the significance to perceive sultriness
3021    REAL(wp), INTENT ( OUT ) ::  sultr_res
3022!
3023!-- Types of coefficients mean deviation: third order polynomial
3024    REAL(wp), PARAMETER ::  dperctka =   7.5776086_wp
3025    REAL(wp), PARAMETER ::  dperctkb = - 0.740603_wp
3026    REAL(wp), PARAMETER ::  dperctkc =   0.0213324_wp
3027    REAL(wp), PARAMETER ::  dperctkd = - 0.00027797237_wp
3028!
3029!-- Types of coefficients mean deviation plus standard deviation
3030!-- regression coefficients: third order polynomial
3031    REAL(wp), PARAMETER ::  dperctsa =   0.0268918_wp
3032    REAL(wp), PARAMETER ::  dperctsb =   0.0465957_wp
3033    REAL(wp), PARAMETER ::  dperctsc = - 0.00054709752_wp
3034    REAL(wp), PARAMETER ::  dperctsd =   0.0000063714823_wp
3035!
3036!-- Factor to mean standard deviation defining SIGNificance for
3037!-- sultriness
3038    REAL(wp), PARAMETER :: faktor = 1.0_wp
3039!
3040!-- Initialise
3041    sultr_res = 99.0_wp
3042    dperctm   = 0.0_wp
3043    dperctstd = 999999.0_wp
3044
3045    IF ( perct_ij < 16.826_wp  .OR.  perct_ij > 56.0_wp )  THEN
3046!
3047!--    Unallowed value of classical perct!
3048       RETURN
3049    ENDIF
3050!
3051!-- Mean deviation dependent on perct
3052    dperctm = dperctka + dperctkb * perct_ij + dperctkc * perct_ij**2.0_wp + dperctkd *            &
3053              perct_ij**3.0_wp
3054!
3055!-- Mean deviation plus its standard deviation
3056    dperctstd = dperctsa + dperctsb * perct_ij + dperctsc * perct_ij**2.0_wp + dperctsd *          &
3057                perct_ij**3.0_wp
3058!
3059!-- Value of the FUNCTION
3060    sultr_res = dperctm + faktor * dperctstd
3061    IF ( ABS( sultr_res ) > 99.0_wp )  sultr_res = +99.0_wp
3062
3063 END SUBROUTINE calc_sultr
3064
3065!--------------------------------------------------------------------------------------------------!
3066! Description:
3067! ------------
3068!> Multiple linear regression to calculate an increment delta_cold, to adjust Fanger's classical PMV
3069!> (pmva) by Gagge's 2 node model, applying Fanger's convective heat transfer coefficient, hcf.
3070!> Wind velocitiy of the reference environment is 0.10 m/s
3071!--------------------------------------------------------------------------------------------------!
3072 SUBROUTINE dpmv_cold( pmva, ta, ws, tmrt, nerr, dpmv_cold_res )
3073
3074    IMPLICIT NONE
3075!
3076!-- Type of input arguments
3077    REAL(wp), INTENT ( IN ) ::  pmva   !< Fanger's classical predicted mean vote
3078    REAL(wp), INTENT ( IN ) ::  ta     !< Air temperature 2 m above ground (degC)
3079    REAL(wp), INTENT ( IN ) ::  tmrt   !< Mean radiant temperature (degC)
3080    REAL(wp), INTENT ( IN ) ::  ws     !< Relative wind velocity 1 m above ground (m/s)
3081!
3082!-- Type of output argument
3083    INTEGER(iwp), INTENT ( OUT ) ::  nerr !< Error indicator: 0 = o.k., +1 = denominator for
3084                                          !< intersection = 0
3085
3086    REAL(wp),     INTENT ( OUT ) ::  dpmv_cold_res    !< Increment to adjust pmva according to the
3087                                                      !< results of Gagge's 2 node model depending on the input
3088!
3089!-- Type of program variables
3090    INTEGER(iwp) ::  i          !< running index
3091    INTEGER(iwp) ::  i_bin      !< result row number
3092
3093    REAL(wp) ::  delta_cold(3)
3094    REAL(wp) ::  dtmrt          !< delta mean radiant temperature
3095    REAL(wp) ::  pmv_cross(2)
3096    REAL(wp) ::  reg_a(3)
3097    REAL(wp) ::  r_denominator  !< the regression equations denominator
3098    REAL(wp) ::  sqrt_ws        !< sqare root of wind speed
3099
3100!    REAL(wp) ::  coeff(3,5)  !< unsafe! array is (re-)writable!
3101!    coeff(1,1:5) =                                                             &
3102!       (/ +0.161_wp,   +0.130_wp, -1.125E-03_wp, +1.106E-03_wp, -4.570E-04_wp /)
3103!    coeff(2,1:5) =                                                             &
3104!       (/  0.795_wp,    0.713_wp, -8.880E-03_wp, -1.803E-03_wp, -2.816E-03_wp /)
3105!    coeff(3,1:5) =                                                             &
3106!       (/ +0.05761_wp, +0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp /)
3107
3108!
3109!-- Coefficient of the 3 regression lines:
3110!      1:const      2:*pmva    3:*ta          4:*sqrt_ws     5:*dtmrt
3111    REAL(wp), DIMENSION(1:3,1:5), PARAMETER ::  coeff = RESHAPE( (/                                &
3112        0.161_wp,   0.130_wp, -1.125E-03_wp,  1.106E-03_wp, -4.570E-04_wp,                         &
3113        0.795_wp,   0.713_wp, -8.880E-03_wp, -1.803E-03_wp, -2.816E-03_wp,                         &
3114        0.05761_wp, 0.458_wp, -1.829E-02_wp, -5.577E-03_wp, -1.970E-03_wp                          &
3115       /), SHAPE( coeff ), order=(/ 2, 1 /)                    )
3116!
3117!-- Initialise
3118    nerr           = 0_iwp
3119    dpmv_cold_res  = 0.0_wp
3120    dtmrt          = tmrt - ta
3121    sqrt_ws        = ws
3122    IF ( sqrt_ws < 0.1_wp )  THEN
3123       sqrt_ws = 0.1_wp
3124    ELSE
3125       sqrt_ws = SQRT( sqrt_ws )
3126    ENDIF
3127
3128    delta_cold = 0.0_wp
3129    pmv_cross = pmva
3130
3131!
3132!-- Determine regression constant for given meteorological conditions
3133    DO  i = 1, 3
3134       reg_a(i) = coeff(i,1) + coeff(i,3) * ta + coeff(i,4) * sqrt_ws + coeff(i,5)*dtmrt
3135       delta_cold(i) = reg_a(i) + coeff(i,2) * pmva
3136    ENDDO
3137!
3138!-- Intersection points of regression lines in terms of Fanger's PMV
3139    DO  i = 1, 2
3140       r_denominator = coeff(i,2) - coeff(i+1,2)
3141       IF ( ABS( r_denominator ) > 0.00001_wp )  THEN
3142          pmv_cross(i) = ( reg_a(i+1) - reg_a(i) ) / r_denominator
3143       ELSE
3144          nerr = 1_iwp
3145          RETURN
3146       ENDIF
3147    ENDDO
3148!
3149!-- Select result row number
3150    i_bin = 3
3151    DO  i = 1, 2
3152       IF ( pmva > pmv_cross(i) )  THEN
3153          i_bin = i
3154          EXIT
3155       ENDIF
3156    ENDDO
3157!
3158!-- Adjust to operative temperature scaled according to classical PMV (Fanger)
3159    dpmv_cold_res = delta_cold(i_bin) - dpmv_cold_adj(pmva)
3160
3161 END SUBROUTINE dpmv_cold
3162
3163!--------------------------------------------------------------------------------------------------!
3164! Description:
3165! ------------
3166!> Calculates the summand dpmv_cold_adj adjusting to the operative temperature scaled according to
3167!> classical PMV (Fanger) for cold conditions. Valid for reference environment: v (1m) = 0.10 m/s,
3168!> dTMRT = 0 K, r.h. = 50 %
3169!--------------------------------------------------------------------------------------------------!
3170 REAL(wp) FUNCTION dpmv_cold_adj( pmva )
3171
3172    IMPLICIT NONE
3173
3174    INTEGER(iwp)  ::  i       !< running index
3175    INTEGER(iwp)  ::  thr     !< thermal range
3176
3177    REAL(wp), INTENT ( IN ) ::  pmva        !< (adjusted) Predicted Mean Vote
3178
3179    REAL(wp)  ::  pmv     !< pmv-part of the regression
3180
3181!
3182!-- Provide regression coefficients for three thermal ranges:
3183!--                                                    slightly cold  cold           very cold
3184    REAL(wp), DIMENSION(1:3,0:3), PARAMETER ::  coef = RESHAPE( (/                                 &
3185                                                       0.0941540_wp, -0.1506620_wp, -0.0871439_wp, &
3186                                                       0.0783162_wp, -1.0612651_wp,  0.1695040_wp, &
3187                                                       0.1350144_wp, -1.0049144_wp, -0.0167627_wp, &
3188                                                       0.1104037_wp, -0.2005277_wp, -0.0003230_wp  &
3189                                                                 /), SHAPE(coef), order=(/ 1, 2 /) )
3190!
3191!-- Select thermal range
3192    IF ( pmva <= -2.1226_wp )  THEN     !< very cold
3193       thr = 3_iwp
3194    ELSE IF ( pmva <= -1.28_wp )  THEN  !< cold
3195       thr = 2_iwp
3196    ELSE                                !< slightly cold
3197       thr = 1_iwp
3198    ENDIF
3199!
3200!-- Initialize
3201    dpmv_cold_adj = coef(thr,0)
3202    pmv           = 1.0_wp
3203!
3204!-- Calculate pmv adjustment (dpmv_cold_adj)
3205    DO  i = 1, 3
3206       pmv           = pmv * pmva
3207       dpmv_cold_adj = dpmv_cold_adj + coef(thr,i) * pmv
3208    ENDDO
3209
3210    RETURN
3211 END FUNCTION dpmv_cold_adj
3212
3213!--------------------------------------------------------------------------------------------------!
3214! Description:
3215! ------------
3216!> Based on perceived temperature (perct) as input, ireq_neutral determines the required clothing
3217!> insulation (clo) for thermally neutral conditions (neither body cooling nor body heating). It is
3218!> related to the Klima-Michel activity level (134.682 W/m2). IREQ_neutral is only defined for perct
3219!> < 10 (degC)
3220!--------------------------------------------------------------------------------------------------!
3221 REAL(wp) FUNCTION ireq_neutral( perct_ij, ireq_minimal, nerr )
3222
3223    IMPLICIT NONE
3224!
3225!-- Type declaration of arguments
3226    INTEGER(iwp), INTENT ( OUT ) ::  nerr
3227
3228    REAL(wp),     INTENT ( IN )  ::  perct_ij
3229    REAL(wp),     INTENT ( OUT ) ::  ireq_minimal
3230!
3231!-- Type declaration for internal varables
3232    REAL(wp)                     ::  perct02
3233!
3234!-- Initialise
3235    nerr = 0_iwp
3236!
3237!-- Convert perceived temperature from basis 0.1 m/s to basis 0.2 m/s
3238    perct02 = 1.8788_wp + 0.9296_wp * perct_ij
3239!
3240!-- IREQ neutral conditions (thermal comfort)
3241    ireq_neutral = 1.62_wp - 0.0564_wp * perct02
3242!
3243!-- Regression only defined for perct <= 10 (degC)
3244    IF ( ireq_neutral < 0.5_wp )  THEN
3245       IF ( ireq_neutral < 0.48_wp )  THEN
3246          nerr = 1_iwp
3247       ENDIF
3248       ireq_neutral = 0.5_wp
3249    ENDIF
3250!
3251!-- Minimal required clothing insulation: maximal acceptable body cooling
3252    ireq_minimal = 1.26_wp - 0.0588_wp * perct02
3253    IF ( nerr > 0_iwp )  THEN
3254       ireq_minimal = ireq_neutral
3255    ENDIF
3256
3257    RETURN
3258 END FUNCTION ireq_neutral
3259
3260
3261!--------------------------------------------------------------------------------------------------!
3262! Description:
3263! ------------
3264!> The SUBROUTINE surface area calculates the surface area of the individual according to its height
3265!> (m), weight (kg), and age (y)
3266!--------------------------------------------------------------------------------------------------!
3267 SUBROUTINE surface_area( height_cm, weight, age, surf )
3268
3269    IMPLICIT NONE
3270
3271    INTEGER(iwp), INTENT(in)  ::  age
3272
3273    REAL(wp)    , INTENT(in)  ::  height_cm
3274    REAL(wp)    , INTENT(in)  ::  weight
3275
3276    REAL(wp)    , INTENT(out) ::  surf
3277
3278    REAL(wp)                  ::  height
3279
3280    height = height_cm * 100.0_wp
3281!
3282!-- According to Gehan-George, for children
3283    IF ( age < 19_iwp )  THEN
3284       IF ( age < 5_iwp )  THEN
3285          surf = 0.02667_wp * height**0.42246_wp * weight**0.51456_wp
3286          RETURN
3287       ENDIF
3288       surf = 0.03050_wp * height**0.35129_wp * weight**0.54375_wp
3289       RETURN
3290    ENDIF
3291!
3292!-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if height and weight be
3293!>  known. In: Arch. Int. Med.. 17, 1916, pp. 863:871.
3294    surf = 0.007184_wp * height**0.725_wp * weight**0.425_wp
3295    RETURN
3296
3297 END SUBROUTINE surface_area
3298
3299!--------------------------------------------------------------------------------------------------!
3300! Description:
3301! ------------
3302!> The SUBROUTINE persdat calculates
3303!>  - the total internal energy production = metabolic + workload,
3304!>  - the total internal energy production for a standardized surface (actlev)
3305!>  - the DuBois - area (a_surf [m2])
3306!> from
3307!>  - the persons age (years),
3308!>  - weight (kg),
3309!>  - height (m),
3310!>  - sex (1 = male, 2 = female),
3311!>  - work load (W)
3312!> for a sample human.
3313!--------------------------------------------------------------------------------------------------!
3314 SUBROUTINE persdat( age, weight, height, sex, work, a_surf, actlev )
3315
3316    IMPLICIT NONE
3317
3318    INTEGER(iwp), INTENT(in) ::  sex
3319
3320
3321    REAL(wp), INTENT(in) ::  age
3322    REAL(wp), INTENT(in) ::  height
3323    REAL(wp), INTENT(in) ::  weight
3324    REAL(wp), INTENT(in) ::  work
3325
3326    REAL(wp), INTENT(out) ::  actlev
3327
3328    REAL(wp) ::  a_surf
3329    REAL(wp) ::  basic_heat_prod
3330    REAL(wp) ::  energy_prod
3331    REAL(wp) ::  factor
3332    REAL(wp) ::  s
3333
3334
3335    CALL surface_area( height, weight, INT( age ), a_surf )
3336    s = height * 100.0_wp / ( weight**( 1.0_wp / 3.0_wp ) )
3337    factor = 1.0_wp + .004_wp  * ( 30.0_wp - age )
3338    basic_heat_prod = 0.0_wp
3339    IF ( sex == 1_iwp )  THEN
3340       basic_heat_prod = 3.45_wp * weight**( 3.0_wp / 4.0_wp ) * ( factor + 0.01_wp                &
3341                         * ( s - 43.4_wp ) )
3342    ELSE IF ( sex == 2_iwp )  THEN
3343       basic_heat_prod = 3.19_wp * weight**( 3.0_wp / 4.0_wp ) * ( factor + 0.018_wp               &
3344                         * ( s - 42.1_wp ) )
3345    ENDIF
3346
3347    energy_prod = work + basic_heat_prod
3348    actlev = energy_prod / a_surf
3349
3350 END SUBROUTINE persdat
3351
3352
3353!--------------------------------------------------------------------------------------------------!
3354! Description:
3355! ------------
3356!> SUBROUTINE ipt_init
3357!> initializes the instationary perceived temperature
3358!--------------------------------------------------------------------------------------------------!
3359
3360 SUBROUTINE ipt_init( age, weight, height, sex, work, actlev, clo, ta, vp, ws, tmrt, pair, dt,     &
3361                      storage, t_clothing, ipt )
3362
3363    IMPLICIT NONE
3364!
3365!-- Input parameters
3366
3367    INTEGER(iwp), INTENT(in)  :: sex    !< Persons sex (1 = male, 2 = female)
3368
3369    REAL(wp), INTENT(in) ::  age        !< Persons age          (years)
3370    REAL(wp), INTENT(in) ::  dt         !< Timestep             (s)
3371    REAL(wp), INTENT(in) ::  height     !< Persons height       (m)7
3372    REAL(wp), INTENT(in) ::  pair       !< Air pressure         (hPa)
3373    REAL(wp), INTENT(in) ::  ta         !< Air Temperature      (degree_C)
3374    REAL(wp), INTENT(in) ::  tmrt       !< Mean radiant temperature   (degree_C)
3375    REAL(wp), INTENT(in) ::  vp         !< Vapor pressure       (hPa)
3376    REAL(wp), INTENT(in) ::  weight     !< Persons weight       (kg)
3377    REAL(wp), INTENT(in) ::  work       !< Current workload     (W)
3378    REAL(wp), INTENT(in) ::  ws         !< Wind speed in approx. 1.1m (m/s)
3379!
3380!-- Output parameters
3381    REAL(wp), INTENT(out) ::  actlev
3382    REAL(wp), INTENT(out) ::  clo
3383    REAL(wp), INTENT(out) ::  ipt
3384    REAL(wp), INTENT(out) ::  storage
3385    REAL(wp), INTENT(out) ::  t_clothing
3386!
3387!-- Internal variables
3388    REAL(wp), PARAMETER :: eps = 0.0005_wp
3389    REAL(wp), PARAMETER :: eta = 0.0_wp
3390
3391    INTEGER(iwp) ::  ncount
3392    INTEGER(iwp) ::  nerr_cold
3393    INTEGER(iwp) ::  nerr
3394
3395    LOGICAL ::  sultrieness
3396
3397!    REAL(wp) ::  acti
3398    REAL(wp) ::  a_surf
3399!    REAL(wp) ::  clo_fanger
3400    REAL(wp) ::  clon
3401    REAL(wp) ::  d_pmv
3402    REAL(wp) ::  d_std
3403    REAL(wp) ::  dgtcm
3404    REAL(wp) ::  dgtcstd
3405    REAL(wp) ::  ireq_minimal
3406    REAL(wp) ::  pmv_s
3407    REAL(wp) ::  pmv_w
3408    REAL(wp) ::  pmva
3409    REAL(wp) ::  pmvs
3410    REAL(wp) ::  ptc
3411    REAL(wp) ::  sclo
3412    REAL(wp) ::  sult_lim
3413    REAL(wp) ::  svp_ta
3414    REAL(wp) ::  wclo
3415
3416
3417    storage = 0.0_wp
3418    CALL persdat( age, weight, height, sex, work, a_surf, actlev )
3419!
3420!-- Initialise
3421    t_clothing = bio_fill_value
3422    ipt        = bio_fill_value
3423    nerr       = 0_wp
3424    ncount     = 0_wp
3425    sultrieness    = .FALSE.
3426!
3427!-- Tresholds: clothing insulation (account for model inaccuracies)
3428!-- Summer clothing
3429    sclo     = 0.44453_wp
3430!-- Winter clothing
3431    wclo     = 1.76267_wp
3432!
3433!-- Decision: firstly calculate for winter or summer clothing
3434    IF ( ta <= 10.0_wp )  THEN
3435!
3436!--    First guess: winter clothing insulation: cold stress
3437       clo = wclo
3438       t_clothing = bio_fill_value  ! force initial run
3439       CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage, dt,    &
3440                            pmva )
3441       pmv_w = pmva
3442
3443       IF ( pmva > 0.0_wp )  THEN
3444!
3445!--       Case summer clothing insulation: heat load ?
3446          clo = sclo
3447          t_clothing = bio_fill_value  ! force initial run
3448          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage, dt, &
3449                               pmva )
3450          pmv_s = pmva
3451          IF ( pmva <= 0.0_wp )  THEN
3452!
3453!--          Case: comfort achievable by varying clothing insulation between winter and summer set
3454!--                values
3455             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta , sclo, pmv_s, wclo, pmv_w, eps,&
3456                               pmva, ncount, clo )
3457             IF ( ncount < 0_iwp )  THEN
3458                nerr = -1_iwp
3459                RETURN
3460             ENDIF
3461          ELSE IF ( pmva > 0.06_wp )  THEN
3462             clo = 0.5_wp
3463             t_clothing = bio_fill_value
3464             CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage,  &
3465                                  dt, pmva )
3466          ENDIF
3467       ELSE IF ( pmva < - 0.11_wp )  THEN
3468          clo = 1.75_wp
3469          t_clothing = bio_fill_value
3470          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage, dt, &
3471                               pmva )
3472       ENDIF
3473
3474    ELSE
3475!
3476!--    First guess: summer clothing insulation: heat load
3477       clo = sclo
3478       t_clothing = bio_fill_value
3479       CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage, dt,    &
3480                            pmva )
3481       pmv_s = pmva
3482
3483       IF ( pmva < 0.0_wp )  THEN
3484!
3485!--       Case winter clothing insulation: cold stress ?
3486          clo = wclo
3487          t_clothing = bio_fill_value
3488          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage, dt, &
3489                               pmva )
3490          pmv_w = pmva
3491
3492          IF ( pmva >= 0.0_wp )  THEN
3493!
3494!--          Case: comfort achievable by varying clothing insulation between winter and summer set
3495!--                values
3496             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo, pmv_s, wclo, pmv_w, eps, &
3497                               pmva, ncount, clo )
3498             IF ( ncount < 0_wp )  THEN
3499                nerr = -1_iwp
3500                RETURN
3501             ENDIF
3502          ELSE IF ( pmva < - 0.11_wp )  THEN
3503             clo = 1.75_wp
3504             t_clothing = bio_fill_value
3505             CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage,  &
3506                                  dt, pmva )
3507          ENDIF
3508       ELSE IF ( pmva > 0.06_wp )  THEN
3509          clo = 0.5_wp
3510          t_clothing = bio_fill_value
3511          CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage, dt, &
3512                               pmva )
3513       ENDIF
3514
3515    ENDIF
3516!
3517!-- Determine perceived temperature by regression equation + adjustments
3518    pmvs = pmva
3519    CALL perct_regression( pmva, clo, ipt )
3520    ptc = ipt
3521    IF ( clo >= 1.75_wp  .AND.  pmva <= - 0.11_wp )  THEN
3522!
3523!--    Adjust for cold conditions according to Gagge 1986
3524       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
3525       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
3526       pmvs = pmva - d_pmv
3527       IF ( pmvs > - 0.11_wp )  THEN
3528          d_pmv  = 0.0_wp
3529          pmvs   = - 0.11_wp
3530       ENDIF
3531       CALL perct_regression( pmvs, clo, ipt )
3532    ENDIF
3533!     clo_fanger = clo
3534    clon = clo
3535    IF ( clo > 0.5_wp  .AND.  ipt <= 8.73_wp )  THEN
3536!
3537!--    Required clothing insulation (ireq) is exclusively defined for perceived temperatures (ipt)
3538!--    less 10 (C) for a reference wind of 0.2 m/s according to 8.73 (C) for 0.1 m/s
3539       clon = ireq_neutral ( ipt, ireq_minimal, nerr )
3540       clo = clon
3541    ENDIF
3542    CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
3543    sultrieness    = .FALSE.
3544    d_std      = - 99.0_wp
3545    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
3546!
3547!--    Adjust for warm/humid conditions according to Gagge 1986
3548       CALL saturation_vapor_pressure ( ta, svp_ta )
3549       d_pmv  = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
3550       pmvs   = pmva + d_pmv
3551       CALL perct_regression( pmvs, clo, ipt )
3552       IF ( sult_lim < 99.0_wp )  THEN
3553          IF ( (ipt - ptc) > sult_lim )  sultrieness = .TRUE.
3554       ENDIF
3555    ENDIF
3556
3557
3558 END SUBROUTINE ipt_init
3559
3560!--------------------------------------------------------------------------------------------------!
3561! Description:
3562! ------------
3563!> SUBROUTINE ipt_cycle
3564!> Calculates one timestep for the instationary version of perceived temperature (iPT, degree_C) for
3565!>  - standard measured/predicted meteorological values and TMRT as input;
3566!>  - regressions for determination of PT;
3567!>  - adjustment to Gagge's PMV* (2-node-model, 1986) as base of PT under warm/humid conditions
3568!>    (Icl= 0.50 clo) and under cold conditions (Icl= 1.75 clo)
3569!--------------------------------------------------------------------------------------------------!
3570 SUBROUTINE ipt_cycle( ta, vp, ws, tmrt, pair, dt, storage, t_clothing, clo, actlev, work, ipt )
3571
3572    IMPLICIT NONE
3573!
3574!-- Type of input of the argument list
3575    REAL(wp), INTENT ( IN )  ::  actlev  !< Internal heat production    (W)
3576    REAL(wp), INTENT ( IN )  ::  clo     !< Clothing index              (no dim)
3577    REAL(wp), INTENT ( IN )  ::  dt      !< Timestep                    (s)
3578    REAL(wp), INTENT ( IN )  ::  pair    !< Air pressure                (hPa)
3579    REAL(wp), INTENT ( IN )  ::  ta      !< Air temperature             (degree_C)
3580    REAL(wp), INTENT ( IN )  ::  tmrt    !< Mean radiant temperature    (degree_C)
3581    REAL(wp), INTENT ( IN )  ::  vp      !< Vapor pressure              (hPa)
3582    REAL(wp), INTENT ( IN )  ::  work    !< Mechanical work load        (W)
3583    REAL(wp), INTENT ( IN )  ::  ws      !< Wind speed                  (m/s)
3584!
3585!-- In and output parameters
3586    REAL(wp), INTENT (INOUT) ::  storage     !< Heat storage            (W)
3587    REAL(wp), INTENT (INOUT) ::  t_clothing  !< Clothig temperature     (degree_C)
3588!
3589!-- Type of output of the argument list
3590    REAL(wp), INTENT ( OUT ) ::  ipt  !< Instationary perceived temperature (degree_C)
3591!
3592!-- Type of internal variables
3593    INTEGER(iwp) ::  nerr
3594    INTEGER(iwp) ::  nerr_cold
3595
3596    LOGICAL ::  sultrieness
3597
3598    REAL(wp) ::  d_pmv
3599    REAL(wp) ::  d_std
3600    REAL(wp) ::  dgtcm
3601    REAL(wp) ::  dgtcstd
3602    REAL(wp) ::  pmva
3603    REAL(wp) ::  pmvs
3604    REAL(wp) ::  ptc
3605    REAL(wp) ::  sult_lim
3606    REAL(wp) ::  svp_ta
3607!
3608!-- Initialise
3609    ipt = bio_fill_value
3610
3611    nerr     = 0_iwp
3612    sultrieness  = .FALSE.
3613!
3614!-- Determine pmv_adjusted for current conditions
3615    CALL fanger_s_acti ( ta, tmrt, vp, ws, pair, clo, actlev, work, t_clothing, storage, dt, pmva )
3616!
3617!-- Determine perceived temperature by regression equation + adjustments
3618    CALL perct_regression( pmva, clo, ipt )
3619!
3620!-- Consider cold conditions
3621    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp )  THEN
3622!
3623!--    Adjust for cold conditions according to Gagge 1986
3624       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
3625       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
3626       pmvs = pmva - d_pmv
3627       IF ( pmvs > - 0.11_wp )  THEN
3628          d_pmv  = 0.0_wp
3629          pmvs   = - 0.11_wp
3630       ENDIF
3631       CALL perct_regression( pmvs, clo, ipt )
3632    ENDIF
3633!
3634!-- Consider sultriness if appropriate
3635    ptc = ipt
3636    CALL calc_sultr( ptc, dgtcm, dgtcstd, sult_lim )
3637    sultrieness = .FALSE.
3638    d_std       = - 99.0_wp
3639    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp )  THEN
3640!
3641!--    Adjust for warm/humid conditions according to Gagge 1986
3642       CALL saturation_vapor_pressure ( ta, svp_ta )
3643       d_pmv = deltapmv ( pmva, ta, vp, svp_ta, tmrt, ws, nerr )
3644       pmvs  = pmva + d_pmv
3645       CALL perct_regression( pmvs, clo, ipt )
3646       IF ( sult_lim < 99.0_wp )  THEN
3647          IF ( (ipt - ptc) > sult_lim )  sultrieness = .TRUE.
3648       ENDIF
3649    ENDIF
3650
3651 END SUBROUTINE ipt_cycle
3652
3653!--------------------------------------------------------------------------------------------------!
3654! Description:
3655! ------------
3656!> SUBROUTINE fanger_s calculates the actual Predicted Mean Vote (dimensionless) according to Fanger
3657!> corresponding to meteorological (ta,tmrt,pa,ws,pair) and individual variables (clo, actlev, eta)
3658!> considering a storage and clothing temperature for a given timestep.
3659!--------------------------------------------------------------------------------------------------!
3660 SUBROUTINE fanger_s_acti( ta, tmrt, pa, in_ws, pair, in_clo, actlev, activity, t_cloth, s, dt,    &
3661                           pmva )
3662
3663    IMPLICIT NONE
3664!
3665!--  Input argument types
3666    REAL(wp), INTENT ( IN )  ::  activity !< Work load                (W/m²)
3667    REAL(wp), INTENT ( IN )  ::  actlev   !< Metabolic + work energy  (W/m²)
3668    REAL(wp), INTENT ( IN )  ::  dt       !< Timestep                 (s)
3669    REAL(wp), INTENT ( IN )  ::  in_clo   !< Clothing index (clo)     (no dim)
3670    REAL(wp), INTENT ( IN )  ::  in_ws    !< Wind speed               (m/s)
3671    REAL(wp), INTENT ( IN )  ::  pa       !< Vapour pressure          (hPa)
3672    REAL(wp), INTENT ( IN )  ::  pair     !< Air pressure             (hPa)
3673    REAL(wp), INTENT ( IN )  ::  ta       !< Air temperature          (degree_C)
3674    REAL(wp), INTENT ( IN )  ::  tmrt     !< Mean radiant temperature (degree_C)
3675!
3676!-- Output argument types
3677    REAL(wp), INTENT ( OUT ) ::  pmva  !< actual Predicted Mean Vote (no dim)
3678
3679    REAL(wp), INTENT (INOUT) ::  s  !< storage var. of energy balance (W/m2)
3680    REAL(wp), INTENT (INOUT) ::  t_cloth  !< clothing temperature (degree_C)
3681!
3682!-- Internal variables
3683    REAL(wp), PARAMETER  ::  time_equil = 7200.0_wp
3684
3685    INTEGER(iwp) :: i         !< running index
3686    INTEGER(iwp) ::  niter    !< Running index
3687
3688    REAL(wp) ::  adjustrate        !< Max storage adjustment rate
3689    REAL(wp) ::  adjustrate_cloth  !< max clothing temp. adjustment rate
3690    REAL(wp) ::  bc                !< preliminary result storage
3691    REAL(wp) ::  cc                !< preliminary result storage
3692    REAL(wp) ::  clo               !< clothing insulation index              (clo)
3693    REAL(wp) ::  d_s               !< Storage delta                          (W)
3694    REAL(wp) ::  dc                !< preliminary result storage
3695    REAL(wp) ::  en                !< Energy ballance                        (W)
3696    REAL(wp) ::  ec                !< preliminary result storage
3697    REAL(wp) ::  f_cl              !< Increase in surface due to clothing    (factor)
3698    REAL(wp) ::  gc                !< preliminary result storage
3699    REAL(wp) ::  heat_convection   !< energy loss by autocnvection       (W)
3700!    REAL(wp) ::  hr                !< radiational heat resistence
3701    REAL(wp) ::  t_clothing        !< clothing temperature                   (degree_C)
3702    REAL(wp) ::  t_skin_aver       !< average skin temperature               (degree_C)
3703    REAL(wp) ::  ws                !< wind speed                             (m/s)
3704    REAL(wp) ::  z1                !< Empiric factor for the adaption of the heat
3705                                   !< ballance equation to the psycho-physical scale
3706                                   !< (Equ. 40 in FANGER)
3707    REAL(wp) ::  z2                !< Water vapour diffution through the skin
3708    REAL(wp) ::  z3                !< Sweat evaporation from the skin surface
3709    REAL(wp) ::  z4                !< Loss of latent heat through respiration
3710    REAL(wp) ::  z5                !< Loss of radiational heat
3711    REAL(wp) ::  z6                !< Heat loss through forced convection
3712
3713
3714
3715
3716!
3717!-- Clo must be > 0. to avoid div. by 0!
3718    clo = in_clo
3719    IF ( clo < 001.0_wp )  clo = 0.001_wp
3720!
3721!-- Increase in surface due to clothing
3722    f_cl = 1.0_wp + 0.15_wp * clo
3723!
3724!-- Case of free convection (ws < 0.1 m/s ) not considered
3725    ws = in_ws
3726    IF ( ws < 0.1_wp )  THEN
3727       ws = 0.1_wp
3728    ENDIF
3729!
3730!-- Heat_convection = forced convection
3731    heat_convection = 12.1_wp * SQRT( ws * pair / 1013.25_wp )
3732!
3733!-- Average skin temperature
3734    t_skin_aver = 35.7_wp - 0.0275_wp * activity
3735!
3736!-- Calculation of constants for evaluation below
3737    bc = 0.155_wp * clo * 3.96_wp * 10.0_wp**( -8.0_wp ) * f_cl
3738    cc = f_cl * heat_convection
3739    ec = 0.155_wp * clo
3740    dc = ( 1.0_wp + ec * cc ) / bc
3741    gc = ( t_skin_aver + bc * ( tmrt + 273.2_wp )**4.0_wp + ec * cc * ta ) / bc
3742!
3743!-- Calculation of clothing surface temperature (t_clothing) based on Newton-approximation with air
3744!-- temperature as initial guess
3745    niter = INT( dt * 10.0_wp, KIND=iwp )
3746    IF ( niter < 1 )  niter = 1_iwp
3747    adjustrate = 1.0_wp - EXP( -1.0_wp * ( 10.0_wp / time_equil ) * dt )
3748    IF ( adjustrate >= 1.0_wp )  adjustrate = 1.0_wp
3749    adjustrate_cloth = adjustrate * 30.0_wp
3750    t_clothing = t_cloth
3751!
3752!-- Set initial values for niter, adjustrates and t_clothing if this is the first call
3753    IF ( t_cloth <= -998.0_wp )  THEN  ! If initial run
3754       niter = 3_iwp
3755       adjustrate = 1.0_wp
3756       adjustrate_cloth = 1.0_wp
3757       t_clothing = ta
3758    ENDIF
3759!
3760!-- Update clothing temperature
3761    DO  i = 1, niter
3762       t_clothing = t_clothing - adjustrate_cloth * ( ( t_clothing + 273.2_wp )**4.0_wp  +         &
3763                    t_clothing * dc - gc ) / ( 4.0_wp * ( t_clothing + 273.2_wp )**3.0_wp + dc )
3764    ENDDO
3765!
3766!-- Empiric factor for the adaption of the heat ballance equation to the psycho-physical scale
3767!-- (Equ. 40 in FANGER)
3768    z1 = ( 0.303_wp * EXP( - 0.036_wp * actlev ) + 0.0275_wp )
3769!
3770!-- Water vapour diffution through the skin
3771    z2 = 0.31_wp * ( 57.3_wp - 0.07_wp * activity-pa )
3772!
3773!-- Sweat evaporation from the skin surface
3774    z3 = 0.42_wp * ( activity - 58.0_wp )
3775!
3776!-- Loss of latent heat through respiration
3777    z4 = 0.0017_wp * actlev * ( 58.7_wp - pa ) + 0.0014_wp * actlev * ( 34.0_wp - ta )
3778!
3779!-- Loss of radiational heat
3780    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + 273.2_wp )**4 - ( tmrt + 273.2_wp )**4 )
3781!
3782!-- Heat loss through forced convection
3783    z6 = cc * ( t_clothing - ta )
3784!
3785!-- Write together as energy ballance
3786    en = activity - z2 - z3 - z4 - z5 - z6
3787!
3788!-- Manage storage
3789    d_s = adjustrate * en + ( 1.0_wp - adjustrate ) * s
3790!
3791!-- Predicted Mean Vote
3792    pmva = z1 * d_s
3793!
3794!-- Update storage
3795    s = d_s
3796    t_cloth = t_clothing
3797
3798 END SUBROUTINE fanger_s_acti
3799
3800
3801
3802!--------------------------------------------------------------------------------------------------!
3803!
3804! Description:
3805! ------------
3806!> Physiologically Equivalent Temperature (PET),
3807!> stationary (calculated based on MEMI),
3808!> Subroutine based on PETBER vers. 1.5.1996 by P. Hoeppe
3809!--------------------------------------------------------------------------------------------------!
3810
3811 SUBROUTINE calculate_pet_static( ta, vpa, v, tmrt, pair, pet_ij )
3812
3813    IMPLICIT NONE
3814!
3815!-- Input arguments:
3816    REAL(wp), INTENT( IN ) ::  pair  !< Air pressure                (hPa)
3817    REAL(wp), INTENT( IN ) ::  ta    !< Air temperature             (degree_C)
3818    REAL(wp), INTENT( IN ) ::  tmrt  !< Mean radiant temperature    (degree_C)
3819    REAL(wp), INTENT( IN ) ::  v     !< Wind speed                  (m/s)
3820    REAL(wp), INTENT( IN ) ::  vpa   !< Vapor pressure              (hPa)
3821!
3822!-- Output arguments:
3823    REAL(wp), INTENT ( OUT ) ::  pet_ij  !< PET                     (degree_C)
3824!
3825!-- Internal variables:
3826    REAL(wp) ::  acl        !< clothing area                        (m²)
3827    REAL(wp) ::  adu        !< Du Bois area                         (m²)
3828    REAL(wp) ::  aeff       !< effective area                       (m²)
3829    REAL(wp) ::  ere        !< energy ballance                      (W)
3830    REAL(wp) ::  erel       !< latent energy ballance               (W)
3831    REAL(wp) ::  esw        !< Energy-loss through sweat evap.      (W)
3832    REAL(wp) ::  facl       !< Surface area extension through clothing (factor)
3833    REAL(wp) ::  feff       !< Surface modification by posture      (factor)
3834    REAL(wp) ::  rdcl       !< Diffusion resistence of clothing     (factor)
3835    REAL(wp) ::  rdsk       !< Diffusion resistence of skin         (factor)
3836    REAL(wp) ::  rtv
3837    REAL(wp) ::  vpts       !< Sat. vapor pressure over skin        (hPa)
3838    REAL(wp) ::  tcl        !< Clothing temperature                 (degree_C)
3839    REAL(wp) ::  tsk        !< Skin temperature                     (degree_C)
3840    REAL(wp) ::  wetsk      !< Fraction of wet skin                 (factor)
3841!
3842!-- Variables:
3843    REAL(wp) :: int_heat    !< Internal heat        (W)
3844!
3845!-- MEMI configuration
3846    REAL(wp) :: age         !< Persons age          (a)
3847    REAL(wp) :: clo         !< Clothing insulation index (clo)
3848    REAL(wp) :: eta         !< Work efficiency      (dimensionless)
3849    REAL(wp) :: fcl         !< Surface area modification by clothing (factor)
3850    REAL(wp) :: ht          !< Persons height       (m)
3851    REAL(wp) :: mbody       !< Persons body mass    (kg)
3852    REAL(wp) :: work        !< Work load            (W)
3853!    INTEGER(iwp) :: pos     !< Posture: 1 = standing, 2 = sitting
3854!    INTEGER(iwp) :: sex     !< Sex: 1 = male, 2 = female
3855!
3856!-- Configuration, keep standard parameters!
3857    age   = 35.0_wp
3858    mbody = 75.0_wp
3859    ht    =  1.75_wp
3860    work  = 80.0_wp
3861    eta   =  0.0_wp
3862    clo   =  0.9_wp
3863    fcl   =  1.15_wp
3864!
3865!-- Call subfunctions
3866    CALL in_body( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, vpa, work )
3867
3868    CALL heat_exch( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, ht, int_heat, mbody,     &
3869                    pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa, vpts, wetsk )
3870
3871    CALL pet_iteration( acl, adu, aeff, esw, facl, feff, int_heat, pair, rdcl, rdsk, rtv, ta, tcl, &
3872                        tsk, pet_ij, vpts, wetsk )
3873
3874
3875 END SUBROUTINE calculate_pet_static
3876
3877
3878!--------------------------------------------------------------------------------------------------!
3879! Description:
3880! ------------
3881!> Calculate internal energy ballance
3882!--------------------------------------------------------------------------------------------------!
3883 SUBROUTINE in_body( age, eta, ere, erel, ht, int_heat, mbody, pair, rtv, ta, vpa, work )
3884!
3885!-- Input arguments:
3886    REAL(wp), INTENT( IN )  ::  age       !< Persons age              (a)
3887    REAL(wp), INTENT( IN )  ::  eta       !< Work efficiency     (dimensionless)
3888    REAL(wp), INTENT( IN )  ::  ht        !< Persons height           (m)
3889    REAL(wp), INTENT( IN )  ::  mbody     !< Persons body mass        (kg)
3890    REAL(wp), INTENT( IN )  ::  pair      !< air pressure             (hPa)
3891    REAL(wp), INTENT( IN )  ::  ta        !< air temperature          (degree_C)
3892    REAL(wp), INTENT( IN )  ::  vpa       !< vapor pressure           (hPa)
3893    REAL(wp), INTENT( IN )  ::  work      !< Work load                (W)
3894!
3895!-- Output arguments:
3896    REAL(wp), INTENT( OUT ) ::  ere       !< energy ballance          (W)
3897    REAL(wp), INTENT( OUT ) ::  erel      !< latent energy ballance   (W)
3898    REAL(wp), INTENT( OUT ) ::  int_heat  !< internal heat production (W)
3899    REAL(wp), INTENT( OUT ) ::  rtv       !< respiratory volume
3900!
3901!-- Internal variables:
3902    REAL(wp) ::  eres                     !< Sensible respiratory heat flux (W)
3903    REAL(wp) ::  met
3904    REAL(wp) ::  tex
3905    REAL(wp) ::  vpex
3906
3907!
3908!-- Metabolic heat production
3909    met = 3.45_wp * mbody**( 3.0_wp / 4.0_wp ) * (1.0_wp + 0.004_wp *                              &
3910          ( 30.0_wp - age) + 0.010_wp * ( ( ht * 100.0_wp /                                        &
3911          ( mbody**( 1.0_wp / 3.0_wp ) ) ) - 43.4_wp ) )
3912    met = work + met
3913    int_heat = met * (1.0_wp - eta)
3914!
3915!-- Sensible respiration energy
3916    tex  = 0.47_wp * ta + 21.0_wp
3917    rtv  = 1.44_wp * 10.0_wp**(-6.0_wp) * met
3918    eres = c_p * (ta - tex) * rtv
3919!
3920!-- Latent respiration energy
3921    vpex = 6.11_wp * 10.0_wp**( 7.45_wp * tex / ( 235.0_wp + tex ) )
3922    erel = 0.623_wp * l_v / pair * ( vpa - vpex ) * rtv
3923!
3924!-- Sum of the results
3925    ere = eres + erel
3926
3927 END SUBROUTINE in_body
3928
3929
3930!--------------------------------------------------------------------------------------------------!
3931! Description:
3932! ------------
3933!> Calculate heat gain or loss
3934!--------------------------------------------------------------------------------------------------!
3935 SUBROUTINE heat_exch( acl, adu, aeff, clo, ere, erel, esw, facl, fcl, feff, ht, int_heat, mbody,  &
3936                       pair, rdcl, rdsk, ta, tcl, tmrt, tsk, v, vpa, vpts, wetsk )
3937
3938!
3939!-- Input arguments:
3940    REAL(wp), INTENT( IN )  ::  clo    !< clothing insulation      (clo)
3941    REAL(wp), INTENT( IN )  ::  fcl    !< factor for surface area increase by clothing
3942    REAL(wp), INTENT( IN )  ::  ere    !< Energy ballance          (W)
3943    REAL(wp), INTENT( IN )  ::  erel   !< Latent energy ballance   (W)
3944    REAL(wp), INTENT( IN )  ::  ht     !< height                   (m)
3945    REAL(wp), INTENT( IN )  ::  int_heat  !< internal heat production (W)
3946    REAL(wp), INTENT( IN )  ::  mbody  !< body mass                (kg)
3947    REAL(wp), INTENT( IN )  ::  pair   !< Air pressure             (hPa)
3948    REAL(wp), INTENT( IN )  ::  ta     !< Air temperature          (degree_C)
3949    REAL(wp), INTENT( IN )  ::  tmrt   !< Mean radiant temperature (degree_C)
3950    REAL(wp), INTENT( IN )  ::  v      !< Wind speed               (m/s)
3951    REAL(wp), INTENT( IN )  ::  vpa    !< Vapor pressure           (hPa)
3952!
3953!-- Output arguments:
3954    REAL(wp), INTENT( OUT ) ::  acl    !< Clothing surface area        (m²)
3955    REAL(wp), INTENT( OUT ) ::  adu    !< Du-Bois area                 (m²)
3956    REAL(wp), INTENT( OUT ) ::  aeff   !< Effective surface area       (m²)
3957    REAL(wp), INTENT( OUT ) ::  esw    !< Energy-loss through sweat evap. (W)
3958    REAL(wp), INTENT( OUT ) ::  facl   !< Surface area extension through clothing (factor)
3959    REAL(wp), INTENT( OUT ) ::  feff   !< Surface modification by posture (factor)
3960    REAL(wp), INTENT( OUT ) ::  rdcl   !< Diffusion resistence of clothing (factor)
3961    REAL(wp), INTENT( OUT ) ::  rdsk   !< Diffusion resistence of skin (factor)
3962    REAL(wp), INTENT( OUT ) ::  tcl    !< Clothing temperature         (degree_C)
3963    REAL(wp), INTENT( OUT ) ::  tsk    !< Skin temperature             (degree_C)
3964    REAL(wp), INTENT( OUT ) ::  vpts   !< Sat. vapor pressure over skin (hPa)
3965    REAL(wp), INTENT( OUT ) ::  wetsk  !< Fraction of wet skin (dimensionless)
3966!
3967!-- Cconstants:
3968!     REAL(wp), PARAMETER :: cair = 1010.0_wp      !< replaced by c_p
3969    REAL(wp), PARAMETER :: cb   = 3640.0_wp        !<
3970    REAL(wp), PARAMETER :: emcl =    0.95_wp      !< Longwave emission coef. of cloth
3971    REAL(wp), PARAMETER :: emsk =    0.99_wp      !< Longwave emission coef. of skin
3972!    REAL(wp), PARAMETER :: evap = 2.42_wp * 10.0_wp **6.0_wp  !< replaced by l_v
3973    REAL(wp), PARAMETER :: food =    0.0_wp        !< Heat gain by food        (W)
3974    REAL(wp), PARAMETER :: po   = 1013.25_wp      !< Air pressure at sea level (hPa)
3975    REAL(wp), PARAMETER :: rob  =    1.06_wp      !<
3976!
3977!-- Internal variables
3978    INTEGER(iwp) ::  count1     !< running index
3979    INTEGER(iwp) ::  count3     !< running index
3980    INTEGER(iwp) ::  j          !< running index
3981    INTEGER(iwp) ::  i          !< running index
3982
3983    LOGICAL ::  skipincreasecount   !< iteration control flag
3984
3985    REAL(wp) ::  cbare          !< Convection through bare skin
3986    REAL(wp) ::  cclo           !< Convection through clothing
3987    REAL(wp) ::  csum           !< Convection in total
3988    REAL(wp) ::  di             !< difference between r1 and r2
3989    REAL(wp) ::  ed             !< energy transfer by diffusion     (W)
3990    REAL(wp) ::  enbal          !< energy ballance                  (W)
3991    REAL(wp) ::  enbal2         !< energy ballance (storage, last cycle)
3992    REAL(wp) ::  eswdif         !< difference between sweat production and evaporation potential
3993    REAL(wp) ::  eswphy         !< sweat created by physiology
3994    REAL(wp) ::  eswpot         !< potential sweat evaporation
3995    REAL(wp) ::  fec            !<
3996    REAL(wp) ::  hc             !<
3997    REAL(wp) ::  he             !<
3998    REAL(wp) ::  htcl           !<
3999    REAL(wp) ::  r1             !<
4000    REAL(wp) ::  r2             !<
4001    REAL(wp) ::  rbare          !< Radiational loss of bare skin    (W/m²)
4002    REAL(wp) ::  rcl            !<
4003    REAL(wp) ::  rclo           !< Radiational loss of clothing     (W/m²)
4004    REAL(wp) ::  rclo2          !< Longwave radiation gain or loss  (W/m²)
4005    REAL(wp) ::  rsum           !< Radiational loss or gain         (W/m²)
4006    REAL(wp) ::  sw             !<
4007!    REAL(wp) ::  swf            !< female factor, currently unused
4008    REAL(wp) ::  swm            !<
4009    REAL(wp) ::  tbody          !<
4010    REAL(wp) ::  vb             !<
4011    REAL(wp) ::  vb1            !<
4012    REAL(wp) ::  vb2            !<
4013    REAL(wp) ::  wd             !<
4014    REAL(wp) ::  wr             !<
4015    REAL(wp) ::  ws             !<
4016    REAL(wp) ::  wsum           !<
4017    REAL(wp) ::  xx             !< modification step                (K)
4018    REAL(wp) ::  y              !< fraction of bare skin
4019
4020    REAL(wp) ::  c(0:10)        !< Core temperature array           (degree_C)
4021    REAL(wp) ::  tcore(1:7)     !<
4022
4023!
4024!-- Initialize
4025    wetsk = 0.0_wp  !< skin is dry everywhere on init (no non-evaporated sweat)
4026!
4027!-- Set Du Bois Area for the sample person
4028    adu = 0.203_wp * mbody**0.425_wp * ht**0.725_wp
4029!
4030!-- Initialize convective heat considering local air preassure
4031    hc = 2.67_wp + ( 6.5_wp * v**0.67_wp )
4032    hc = hc * ( pair / po )**0.55_wp
4033!
4034!-- Set surface modification by posture (the person will always stand)
4035    feff = 0.725_wp                     !< Posture: 0.725 for stading
4036!
4037!-- Set surface modification by clothing
4038    facl = ( - 2.36_wp + 173.51_wp * clo - 100.76_wp * clo * clo + 19.28_wp * ( clo**3.0_wp ) )    &
4039           / 100.0_wp
4040    IF ( facl > 1.0_wp )  facl = 1.0_wp
4041!
4042!-- Initialize heat resistences
4043    rcl = ( clo / 6.45_wp ) / facl
4044    IF ( clo >= 2.0_wp )  y  = 1.0_wp
4045    IF ( ( clo > 0.6_wp )   .AND.  ( clo < 2.0_wp ) )   y = ( ht - 0.2_wp ) / ht
4046    IF ( ( clo <= 0.6_wp )  .AND.  ( clo > 0.3_wp ) )  y = 0.5_wp
4047    IF ( ( clo <= 0.3_wp )  .AND.  ( clo > 0.0_wp ) )   y = 0.1_wp
4048    r2   = adu * ( fcl - 1.0_wp + facl ) / ( 2.0_wp * 3.14_wp * ht * y )
4049    r1   = facl * adu / ( 2.0_wp * 3.14_wp * ht * y )
4050    di   = r2 - r1
4051
4052!
4053!-- Estimate skin temperatur
4054    DO  j = 1, 7
4055
4056       tsk    = 34.0_wp
4057       count1 = 0_iwp
4058       tcl    = ( ta + tmrt + tsk ) / 3.0_wp
4059       count3 = 1_iwp
4060       enbal2 = 0.0_wp
4061
4062       DO  i = 1, 100  ! allow for 100 iterations max
4063          acl   = adu * facl + adu * ( fcl - 1.0_wp )
4064          rclo2 = emcl * sigma_sb * ( ( tcl + degc_to_k )**4.0_wp -                                &
4065                  ( tmrt + degc_to_k )**4.0_wp ) * feff
4066          htcl  = 6.28_wp * ht * y * di / ( rcl * LOG( r2 / r1 ) * acl )
4067          tsk   = 1.0_wp / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl
4068!
4069!--       Radiation saldo
4070          aeff  = adu * feff
4071          rbare = aeff * ( 1.0_wp - facl ) * emsk * sigma_sb *                                     &
4072                  ( ( tmrt + degc_to_k )**4.0_wp - ( tsk + degc_to_k )**4.0_wp )
4073          rclo  = feff * acl * emcl * sigma_sb *                                                   &
4074                  ( ( tmrt + degc_to_k )**4.0_wp - ( tcl + degc_to_k )**4.0_wp )
4075          rsum  = rbare + rclo
4076!
4077!--       Convection
4078          cbare = hc * ( ta - tsk ) * adu * ( 1.0_wp - facl )
4079          cclo  = hc * ( ta - tcl ) * acl
4080          csum  = cbare + cclo
4081!
4082!--       Core temperature
4083          c(0)  = int_heat + ere
4084          c(1)  = adu * rob * cb
4085          c(2)  = 18.0_wp - 0.5_wp * tsk
4086          c(3)  = 5.28_wp * adu * c(2)
4087          c(4)  = 0.0208_wp * c(1)
4088          c(5)  = 0.76075_wp * c(1)
4089          c(6)  = c(3) - c(5) - tsk * c(4)
4090          c(7)  = - c(0) * c(2) - tsk * c(3) + tsk * c(5)
4091          c(8)  = c(6) * c(6) - 4.0_wp * c(4) * c(7)
4092          c(9)  = 5.28_wp * adu - c(5) - c(4) * tsk
4093          c(10) = c(9) * c(9) - 4.0_wp * c(4) * ( c(5) * tsk - c(0) - 5.28_wp * adu * tsk )
4094
4095          IF ( ABS( tsk - 36.0_wp ) < 0.00001_wp )  tsk = 36.01_wp
4096          tcore(7) = c(0) / ( 5.28_wp * adu + c(1) * 6.3_wp / 3600.0_wp ) + tsk
4097          tcore(3) = c(0) / ( 5.28_wp * adu + ( c(1) * 6.3_wp / 3600.0_wp ) /   &
4098                     ( 1.0_wp + 0.5_wp * ( 34.0_wp - tsk ) ) ) + tsk
4099          IF ( c(10) >= 0.0_wp )  THEN
4100             tcore(6) = ( - c(9) - c(10)**0.5_wp ) / ( 2.0_wp * c(4) )
4101             tcore(1) = ( - c(9) + c(10)**0.5_wp ) / ( 2.0_wp * c(4) )
4102          ENDIF
4103
4104          IF ( c(8) >= 0.0_wp )  THEN
4105             tcore(2) = ( - c(6) + ABS( c(8) )**0.5_wp ) / ( 2.0_wp * c(4) )
4106             tcore(5) = ( - c(6) - ABS( c(8) )**0.5_wp ) / ( 2.0_wp * c(4) )
4107             tcore(4) = c(0) / ( 5.28_wp * adu + c(1) * 1.0_wp / 40.0_wp ) + tsk
4108          ENDIF
4109!
4110!--       Transpiration
4111          tbody = 0.1_wp * tsk + 0.9_wp * tcore(j)
4112          swm   = 304.94_wp * ( tbody - 36.6_wp ) * adu / 3600000.0_wp
4113          vpts  = 6.11_wp * 10.0_wp**( 7.45_wp * tsk / ( 235.0_wp + tsk ) )
4114
4115          IF ( tbody <= 36.6_wp )  swm = 0.0_wp  !< no need for sweating
4116
4117          sw = swm
4118          eswphy = - sw * l_v
4119          he     = 0.633_wp * hc / ( pair * c_p )
4120          fec    = 1.0_wp / ( 1.0_wp + 0.92_wp * hc * rcl )
4121          eswpot = he * ( vpa - vpts ) * adu * l_v * fec
4122          wetsk  = eswphy / eswpot
4123
4124          IF ( wetsk > 1.0_wp )  wetsk = 1.0_wp
4125!
4126!--       Sweat production > evaporation?
4127          eswdif = eswphy - eswpot
4128
4129          IF ( eswdif <= 0.0_wp )  esw = eswpot     !< Limit is evaporation
4130          IF ( eswdif > 0.0_wp )   esw = eswphy     !< Limit is sweat production
4131          IF ( esw  > 0.0_wp )     esw = 0.0_wp     !< Sweat can't be evaporated, no more cooling
4132                                                    !< effect
4133!
4134!--       Diffusion
4135          rdsk = 0.79_wp * 10.0_wp**7.0_wp
4136          rdcl = 0.0_wp
4137          ed   = l_v / ( rdsk + rdcl ) * adu * ( 1.0_wp - wetsk ) * ( vpa - vpts )
4138!
4139!--       Max vb
4140          vb1 = 34.0_wp - tsk
4141          vb2 = tcore(j) - 36.6_wp
4142
4143          IF ( vb2 < 0.0_wp )  vb2 = 0.0_wp
4144          IF ( vb1 < 0.0_wp )  vb1 = 0.0_wp
4145          vb = ( 6.3_wp + 75.0_wp * vb2 ) / ( 1.0_wp + 0.5_wp * vb1 )
4146!
4147!--       Energy ballence
4148          enbal = int_heat + ed + ere + esw + csum + rsum + food
4149!
4150!--       Clothing temperature
4151          xx = 0.001_wp
4152          IF ( count1 == 0_iwp )  xx = 1.0_wp
4153          IF ( count1 == 1_iwp )  xx = 0.1_wp
4154          IF ( count1 == 2_iwp )  xx = 0.01_wp
4155          IF ( count1 == 3_iwp )  xx = 0.001_wp
4156
4157          IF ( enbal > 0.0_wp )  tcl = tcl + xx
4158          IF ( enbal < 0.0_wp )  tcl = tcl - xx
4159
4160          skipincreasecount = .FALSE.
4161          IF ( ( (enbal <= 0.0_wp )  .AND.  (enbal2 > 0.0_wp ) )  .OR.                             &
4162             ( ( enbal >= 0.0_wp )   .AND.  ( enbal2 < 0.0_wp ) ) )  THEN
4163             skipincreasecount = .TRUE.
4164          ELSE
4165             enbal2 = enbal
4166             count3 = count3 + 1_iwp
4167          ENDIF
4168
4169          IF ( ( count3 > 200_iwp )  .OR.  skipincreasecount )  THEN
4170             IF ( count1 < 3_iwp )  THEN
4171                count1 = count1 + 1_iwp
4172                enbal2 = 0.0_wp
4173             ELSE
4174                EXIT
4175             ENDIF
4176          ENDIF
4177       ENDDO
4178
4179       IF ( count1 == 3_iwp )  THEN
4180          SELECT CASE ( j )
4181             CASE ( 2, 5)
4182                IF ( .NOT. ( ( tcore(j) >= 36.6_wp )  .AND.  ( tsk <= 34.050_wp ) ) )  CYCLE
4183             CASE ( 6, 1 )
4184                IF ( c(10) < 0.0_wp ) CYCLE
4185                IF ( .NOT. ( ( tcore(j) >= 36.6_wp )  .AND.  ( tsk > 33.850_wp ) ) )  CYCLE
4186             CASE ( 3 )
4187                IF ( .NOT. ( ( tcore(j) < 36.6_wp )  .AND.  ( tsk <= 34.000_wp ) ) )  CYCLE
4188             CASE ( 7 )
4189                IF ( .NOT. ( ( tcore(j) < 36.6_wp )  .AND.  ( tsk > 34.000_wp ) ) )  CYCLE
4190             CASE default
4191          END SELECT
4192       ENDIF
4193
4194       IF ( ( j /= 4_iwp )  .AND.  ( vb >= 91.0_wp ) )  CYCLE
4195       IF ( ( j == 4_iwp )  .AND.  ( vb < 89.0_wp ) )  CYCLE
4196       IF ( vb > 90.0_wp ) vb = 90.0_wp
4197!
4198!--    Loses by water
4199       ws = sw * 3600.0_wp * 1000.0_wp
4200       IF ( ws > 2000.0_wp )  ws = 2000.0_wp
4201       wd = ed / l_v * 3600.0_wp * ( -1000.0_wp )
4202       wr = erel / l_v * 3600.0_wp * ( -1000.0_wp )
4203
4204       wsum = ws + wr + wd
4205
4206       RETURN
4207    ENDDO
4208 END SUBROUTINE heat_exch
4209
4210!--------------------------------------------------------------------------------------------------!
4211! Description:
4212! ------------
4213!> Calculate PET
4214!--------------------------------------------------------------------------------------------------!
4215 SUBROUTINE pet_iteration( acl, adu, aeff, esw, facl, feff, int_heat, pair, rdcl, rdsk, rtv, ta,   &
4216                           tcl, tsk, pet_ij, vpts, wetsk )
4217!
4218!-- Input arguments:
4219    REAL(wp), INTENT( IN ) ::  acl       !< clothing surface area        (m²)
4220    REAL(wp), INTENT( IN ) ::  adu       !< Du-Bois area                 (m²)
4221    REAL(wp), INTENT( IN ) ::  esw       !< energy-loss through sweat evap. (W)
4222    REAL(wp), INTENT( IN ) ::  facl      !< surface area extension through clothing (factor)
4223    REAL(wp), INTENT( IN ) ::  feff      !< surface modification by posture (factor)
4224    REAL(wp), INTENT( IN ) ::  int_heat  !< internal heat production (W)
4225    REAL(wp), INTENT( IN ) ::  pair      !< air pressure                 (hPa)
4226    REAL(wp), INTENT( IN ) ::  rdcl      !< diffusion resistence of clothing (factor)
4227    REAL(wp), INTENT( IN ) ::  rdsk      !< diffusion resistence of skin (factor)
4228    REAL(wp), INTENT( IN ) ::  rtv       !< respiratory volume
4229    REAL(wp), INTENT( IN ) ::  ta        !< air temperature              (degree_C)
4230    REAL(wp), INTENT( IN ) ::  tcl       !< clothing temperature         (degree_C)
4231    REAL(wp), INTENT( IN ) ::  tsk       !< skin temperature             (degree_C)
4232    REAL(wp), INTENT( IN ) ::  vpts      !< sat. vapor pressure over skin (hPa)
4233    REAL(wp), INTENT( IN ) ::  wetsk     !< fraction of wet skin (dimensionless)
4234!
4235!-- Output arguments:
4236    REAL(wp), INTENT( OUT ) ::  aeff     !< effective surface area       (m²)
4237    REAL(wp), INTENT( OUT ) ::  pet_ij   !< PET                          (degree_C)
4238!
4239!-- Cconstants:
4240    REAL(wp), PARAMETER :: emcl =    0.95_wp      !< Longwave emission coef. of cloth
4241    REAL(wp), PARAMETER :: emsk =    0.99_wp      !< Longwave emission coef. of skin
4242    REAL(wp), PARAMETER :: po   = 1013.25_wp      !< Air pressure at sea level (hPa)
4243!
4244!-- Internal variables
4245    INTEGER ( iwp ) ::  count1        !< running index
4246    INTEGER ( iwp ) ::  i             !< running index
4247
4248    REAL ( wp ) ::  cbare             !< Convection through bare skin
4249    REAL ( wp ) ::  cclo              !< Convection through clothing
4250    REAL ( wp ) ::  csum              !< Convection in total
4251    REAL ( wp ) ::  ed                !< Diffusion                      (W)
4252    REAL ( wp ) ::  enbal             !< Energy ballance                (W)
4253    REAL ( wp ) ::  enbal2            !< Energy ballance (last iteration cycle)
4254    REAL ( wp ) ::  ere               !< Energy ballance result         (W)
4255    REAL ( wp ) ::  erel              !< Latent energy ballance         (W)
4256    REAL ( wp ) ::  eres              !< Sensible respiratory heat flux (W)
4257    REAL ( wp ) ::  hc                !<
4258    REAL ( wp ) ::  rbare             !< Radiational loss of bare skin  (W/m²)
4259    REAL ( wp ) ::  rclo              !< Radiational loss of clothing   (W/m²)
4260    REAL ( wp ) ::  rsum              !< Radiational loss or gain       (W/m²)
4261    REAL ( wp ) ::  tex               !< Temperat. of exhaled air       (degree_C)
4262    REAL ( wp ) ::  vpex              !< Vapor pressure of exhaled air  (hPa)
4263    REAL ( wp ) ::  xx                !< Delta PET per iteration        (K)
4264
4265
4266    pet_ij = ta
4267    enbal2 = 0.0_wp
4268
4269    DO  count1 = 0, 3
4270       DO  i = 1, 125  ! 500 / 4
4271          hc = 2.67_wp + 6.5_wp * 0.1_wp**0.67_wp
4272          hc = hc * ( pair / po )**0.55_wp
4273!
4274!--       Radiation
4275          aeff  = adu * feff
4276          rbare = aeff * ( 1.0_wp - facl ) * emsk * sigma_sb *                  &
4277                  ( ( pet_ij + degc_to_k )**4.0_wp - ( tsk + degc_to_k )**4.0_wp )
4278          rclo  = feff * acl * emcl * sigma_sb *                               &
4279                  ( ( pet_ij + degc_to_k )**4.0_wp - ( tcl + degc_to_k )**4.0_wp )
4280          rsum  = rbare + rclo
4281!
4282!--       Covection
4283          cbare = hc * ( pet_ij - tsk ) * adu * ( 1.0_wp - facl )
4284          cclo  = hc * ( pet_ij - tcl ) * acl
4285          csum  = cbare + cclo
4286!
4287!--       Diffusion
4288          ed = l_v / ( rdsk + rdcl ) * adu * ( 1.0_wp - wetsk ) * ( 12.0_wp - vpts )
4289!
4290!--       Respiration
4291          tex  = 0.47_wp * pet_ij + 21.0_wp
4292          eres = c_p * ( pet_ij - tex ) * rtv
4293          vpex = 6.11_wp * 10.0_wp**( 7.45_wp * tex / ( 235.0_wp + tex ) )
4294          erel = 0.623_wp * l_v / pair * ( 12.0_wp - vpex ) * rtv
4295          ere  = eres + erel
4296!
4297!--       Energy ballance
4298          enbal = int_heat + ed + ere + esw + csum + rsum
4299!
4300!--       Iteration concerning ta
4301          xx = 0.001_wp
4302          IF ( count1 == 0_iwp )  xx = 1.0_wp
4303          IF ( count1 == 1_iwp )  xx = 0.1_wp
4304          IF ( count1 == 2_iwp )  xx = 0.01_wp
4305!           IF ( count1 == 3_iwp )  xx = 0.001_wp
4306          IF ( enbal > 0.0_wp )  pet_ij = pet_ij - xx
4307          IF ( enbal < 0.0_wp )  pet_ij = pet_ij + xx
4308          IF ( ( enbal <= 0.0_wp )  .AND.  ( enbal2 > 0.0_wp ) )  EXIT
4309          IF ( ( enbal >= 0.0_wp )  .AND.  ( enbal2 < 0.0_wp ) )  EXIT
4310
4311          enbal2 = enbal
4312       ENDDO
4313    ENDDO
4314 END SUBROUTINE pet_iteration
4315
4316!
4317!-- UVEM specific subroutines
4318
4319!--------------------------------------------------------------------------------------------------!
4320! Description:
4321! ------------
4322!> Module-specific routine for new module
4323!--------------------------------------------------------------------------------------------------!
4324 SUBROUTINE uvem_solar_position
4325
4326    USE control_parameters,                                                                        &
4327       ONLY:  latitude, longitude, time_SINce_reference_point
4328
4329    IMPLICIT NONE
4330
4331    INTEGER(iwp) ::  day_of_year = 0       !< day of year
4332
4333    REAL(wp) ::  alpha         = 0.0_wp    !< solar azimuth angle in radiant
4334    REAL(wp) ::  declination   = 0.0_wp    !< declination
4335    REAL(wp) ::  dtor          = 0.0_wp    !< factor to convert degree to radiant
4336    REAL(wp) ::  js            = 0.0_wp    !< parameter for solar position calculation
4337    REAL(wp) ::  lat           = 52.39_wp  !< latitude
4338    REAL(wp) ::  lon           = 9.7_wp    !< longitude
4339    REAL(wp) ::  second_of_day = 0.0_wp    !< current second of the day
4340    REAL(wp) ::  thetar        = 0.0_wp    !< angle for solar zenith angle calculation
4341    REAL(wp) ::  thetasr       = 0.0_wp    !< angle for solar azimuth angle calculation
4342    REAL(wp) ::  zgl           = 0.0_wp    !< calculated exposure by direct beam
4343    REAL(wp) ::  woz           = 0.0_wp    !< calculated exposure by diffuse radiation
4344    REAL(wp) ::  wsp           = 0.0_wp    !< calculated exposure by direct beam
4345
4346
4347    CALL get_date_time( time_SINce_reference_point, day_of_year = day_of_year,                     &
4348                        second_of_day = second_of_day )
4349    dtor = pi / 180.0_wp
4350    lat = latitude
4351    lon = longitude
4352!
4353!-- Calculation of js, necessary for calculation of equation of time (zgl) :
4354    js=  72.0_wp * ( REAL( day_of_year, KIND = wp ) + ( second_of_day / 86400.0_wp ) ) / 73.0_wp
4355!
4356!-- Calculation of equation of time (zgl):
4357    zgl = 0.0066_wp + 7.3525_wp * COS( ( js + 85.9_wp ) * dtor ) + 9.9359_wp *                     &
4358    COS( ( 2.0_wp * js + 108.9_wp ) * dtor ) + 0.3387_wp * COS( ( 3 * js + 105.2_wp ) * dtor )
4359!
4360!-- Calculation of apparent solar time woz:
4361    woz = ( ( second_of_day / 3600.0_wp ) - ( 4.0_wp * ( 15.0_wp - lon ) ) / 60.0_wp ) +           &
4362          ( zgl / 60.0_wp )
4363!
4364!-- Calculation of hour angle (wsp):
4365    wsp = ( woz - 12.0_wp ) * 15.0_wp
4366!
4367!-- Calculation of declination:
4368    declination = 0.3948_wp - 23.2559_wp * COS( ( js + 9.1_wp ) * dtor ) -                         &
4369                  0.3915_wp * COS( ( 2.0_wp * js + 5.4_wp ) * dtor ) - 0.1764_wp *                 &
4370                  COS( ( 3.0_wp * js + 26.0_wp ) * dtor )
4371!
4372!-- Calculation of solar zenith angle
4373    thetar  = ACOS( SIN( lat * dtor) * SIN( declination * dtor ) + COS( wsp * dtor ) *             &
4374              COS( lat * dtor ) * COS( declination * dtor ) )
4375    thetasr = ASIN( SIN( lat * dtor) * SIN( declination * dtor ) + COS( wsp * dtor ) *             &
4376              COS( lat * dtor ) * COS( declination * dtor ) )
4377    sza = thetar / dtor
4378!
4379!-- calculation of solar azimuth angle
4380    IF (woz <= 12.0_wp) alpha = pi - ACOS( ( SIN(thetasr) * SIN( lat * dtor ) -                    &
4381                                SIN( declination * dtor ) ) / ( COS(thetasr) * COS( lat * dtor ) ) )
4382    IF (woz > 12.0_wp)  alpha = pi + ACOS( ( SIN(thetasr) * SIN( lat * dtor ) -                    &
4383                                SIN( declination * dtor ) ) / ( COS(thetasr) * COS( lat * dtor ) ) )
4384    saa = alpha / dtor
4385
4386 END SUBROUTINE uvem_solar_position
4387
4388
4389!--------------------------------------------------------------------------------------------------!
4390! Description:
4391! ------------
4392!> Module-specific routine for new module
4393!--------------------------------------------------------------------------------------------------!
4394 SUBROUTINE bio_calculate_uv_exposure
4395
4396    USE indices,                                                                                   &
4397        ONLY:  nxl, nxr, nyn, nys
4398
4399
4400    IMPLICIT NONE
4401
4402    INTEGER(iwp) ::  i     !< loop index in x direction
4403    INTEGER(iwp) ::  j     !< loop index in y direction
4404    INTEGER(iwp) ::  szai  !< loop index for different sza values
4405
4406    CALL uvem_solar_position
4407
4408    IF (sza  >=  90)  THEN
4409       vitd3_exposure(:,:) = 0.0_wp
4410    ELSE
4411
4412       DO  ai = 0, 35
4413          DO  zi = 0, 9
4414             projection_area_lookup_table(ai,zi) = uvem_projarea_f%var(clothing,zi,ai)
4415          ENDDO
4416       ENDDO
4417       DO  ai = 0, 35
4418          DO  zi = 0, 9
4419             integration_array(ai,zi) = uvem_integration_f%var(zi,ai)
4420          ENDDO
4421       ENDDO
4422       DO  ai = 0, 2
4423          DO  zi = 0, 90
4424             irradiance_lookup_table(ai,zi) = uvem_irradiance_f%var(zi,ai)
4425          ENDDO
4426       ENDDO
4427       DO  ai = 0, 35
4428          DO  zi = 0, 9
4429             DO  szai = 0, 90
4430                radiance_lookup_table(ai,zi,szai) = uvem_radiance_f%var(szai,zi,ai)
4431             ENDDO
4432          ENDDO
4433       ENDDO
4434
4435
4436!--    Rotate 3D-Model human to desired direction
4437       projection_area_temp( 0:35,:) = projection_area_lookup_table
4438       projection_area_temp(36:71,:) = projection_area_lookup_table
4439       IF (  .NOT.  turn_to_sun ) startpos_human = orientation_angle / 10.0_wp
4440       IF (       turn_to_sun ) startpos_human = saa / 10.0_wp
4441       DO  ai = 0, 35
4442          xfactor = ( startpos_human ) - INT( startpos_human )
4443          DO  zi = 0, 9
4444             projection_area(ai,zi) = ( projection_area_temp( 36 -                                 &
4445                                                               INT( startpos_human ) - 1 + ai , zi)&
4446                                        * ( xfactor ) )                                            &
4447                                      + ( projection_area_temp( 36 -                               &
4448                                                                INT( startpos_human ) + ai , zi)   &
4449                                        * ( 1.0_wp - xfactor ) )
4450          ENDDO
4451       ENDDO
4452!
4453!
4454!--    Interpolate to accurate Solar Zenith Angle
4455       DO  ai = 0, 35
4456          xfactor = ( sza )-INT( sza )
4457          DO  zi = 0, 9
4458             radiance_array(ai,zi) = ( radiance_lookup_table(ai, zi, INT( sza ) ) *                &
4459                                     ( 1.0_wp - xfactor) ) +                                       &
4460                                     ( radiance_lookup_table(ai,zi,INT( sza ) + 1) * xfactor )
4461          ENDDO
4462       ENDDO
4463       DO  iq = 0, 2
4464          irradiance(iq) = ( irradiance_lookup_table(iq, INT( sza ) ) * ( 1.0_wp - xfactor)) +     &
4465                           ( irradiance_lookup_table(iq, INT( sza ) + 1) * xfactor )
4466       ENDDO
4467!
4468!--    Interpolate to accurate Solar Azimuth Angle
4469       IF ( sun_in_south )  THEN
4470          startpos_saa_float = 180.0_wp / 10.0_wp
4471       ELSE
4472          startpos_saa_float = saa / 10.0_wp
4473       ENDIF
4474       radiance_array_temp( 0:35,:) = radiance_array
4475       radiance_array_temp(36:71,:) = radiance_array
4476       xfactor = (startpos_saa_float) - INT( startpos_saa_float )
4477       DO  ai = 0, 35
4478          DO  zi = 0, 9
4479             radiance_array(ai,zi) = ( radiance_array_temp(36 -                                    &
4480                                                           INT( startpos_saa_float ) - 1 + ai, zi) &
4481                                       * ( xfactor ) )                                             &
4482                                     + ( radiance_array_temp(36 -                                  &
4483                                                             INT( startpos_saa_float ) + ai, zi)   &
4484                                       * ( 1.0_wp - xfactor ) )
4485          ENDDO
4486       ENDDO
4487
4488!
4489!--    Calculate Projectionarea for direct beam
4490       projection_area_direct_temp( 0:35,:) = projection_area
4491       projection_area_direct_temp(36:71,:) = projection_area
4492       yfactor = ( sza / 10.0_wp ) - INT( sza / 10.0_wp )
4493       xfactor = ( startpos_saa_float ) - INT( startpos_saa_float )
4494       projection_area_direct_beam = ( projection_area_direct_temp( INT(startpos_saa_float)    ,INT(sza/10.0_wp)  ) * &
4495                                     ( 1.0_wp - xfactor ) * ( 1.0_wp - yfactor ) ) +                                  &
4496                                     ( projection_area_direct_temp( INT(startpos_saa_float) + 1,INT(sza/10.0_wp)  ) * &
4497                                     (          xfactor ) * ( 1.0_wp - yfactor ) ) +                                  &
4498                                     ( projection_area_direct_temp( INT(startpos_saa_float)    ,INT(sza/10.0_wp)+1) * &
4499                                     ( 1.0_wp - xfactor ) * (          yfactor ) ) +                                  &
4500                                     ( projection_area_direct_temp( INT(startpos_saa_float) + 1,INT(sza/10.0_wp)+1) * &
4501                                     (          xfactor ) * (          yfactor ) )
4502
4503
4504       DO  i = nxl, nxr
4505          DO  j = nys, nyn
4506!
4507!--          Extract obstruction from IBSET-Integer_Array
4508             IF (consider_obstructions )  THEN
4509                obstruction_temp1 = building_obstruction_f%var_3d(:,j,i)
4510                IF ( obstruction_temp1(0)  /=  9 )  THEN
4511                   DO  pobi = 0, 44
4512                      DO  bi = 0, 7
4513                         IF ( BTEST( obstruction_temp1(pobi), bi )  .EQV.  .TRUE.)  THEN
4514                            obstruction_temp2( ( pobi * 8 ) + bi ) = 1
4515                         ELSE
4516                            obstruction_temp2( ( pobi * 8 ) + bi ) = 0
4517                         ENDIF
4518                      ENDDO
4519                   ENDDO
4520                   DO  zi = 0, 9
4521                      obstruction(:,zi) = obstruction_temp2( zi * 36 :( zi * 36) + 35 )
4522                   ENDDO
4523                ELSE
4524                   obstruction(:,:) = 0
4525                ENDIF
4526             ENDIF
4527!
4528!--          Calculated human exposure
4529             diffuse_exposure = SUM( radiance_array * projection_area * integration_array *        &
4530                                     obstruction )
4531
4532             obstruction_direct_beam = obstruction( NINT( startpos_saa_float),                     &
4533                                                    NINT( sza / 10.0_wp ) )
4534             IF (sza  >=  89.99_wp)  THEN
4535                sza = 89.99999_wp
4536             ENDIF
4537!
4538!--          Calculate direct normal irradiance (direct beam)
4539             direct_exposure = ( irradiance(1) / COS( pi * sza / 180.0_wp ) ) *                    &
4540                               projection_area_direct_beam * obstruction_direct_beam
4541
4542             vitd3_exposure(j,i) = ( diffuse_exposure + direct_exposure ) / 1000.0_wp * 70.97_wp
4543!            unit = international units vitamin D per second
4544          ENDDO
4545       ENDDO
4546    ENDIF
4547
4548 END SUBROUTINE bio_calculate_uv_exposure
4549
4550 END MODULE biometeorology_mod
Note: See TracBrowser for help on using the repository browser.