Home

Context Navigation

biometeorology_mod.f90 @ 4676

Last change on this file since 4676 was 4633, checked in by suehring, 5 years ago

Biomet: bugfix in check for humidity

Property svn:keywords set to Id

Property svn:mergeinfo set to (toggle deleted branches)

/palm/branches/chemistry/SOURCE/biometeorology_mod.f90	2047-3190,3218-3297
/palm/branches/resler/SOURCE/biometeorology_mod.f90	2023-4492
/palm/branches/salsa/SOURCE/biometeorology_mod.f90	2503-3581
/palm/trunk/SOURCE/biometeorology_mod.f90	merged	eligible
/palm/branches/forwind/SOURCE/biometeorology_mod.f90	1564-1913
/palm/branches/fricke/SOURCE/biometeorology_mod.f90	942-977
/palm/branches/hoffmann/SOURCE/biometeorology_mod.f90	989-1052
/palm/branches/letzel/masked_output/SOURCE/biometeorology_mod.f90	296-409
/palm/branches/palm4u/SOURCE/biometeorology_mod.f90	2540-2692
/palm/branches/rans/SOURCE/biometeorology_mod.f90	2078-3128
/palm/branches/suehring/biometeorology_mod.f90	423-666

File size: 193.6 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: