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

Last change on this file since 4733 was 4633, checked in by suehring, 4 years ago

Biomet: bugfix in check for humidity

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