source: palm/trunk/SOURCE/biometeorology_mod.f90

Last change on this file was 4843, checked in by raasch, 5 months ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

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