Home

Context Navigation

biometeorology_mod.f90 @ 3785

Last change on this file since 3785 was 3753, checked in by dom_dwd_user, 6 years ago

biometeorology_mod.f90:
(C) Added automatic setting of mrt_nlevels in case it was not part of radiation_parameters namelist (or set to 0 accidentially).
(C) Minor speed improvoemnts in perceived temperature calculations.
(C) Perceived temperature regression arrays now declared as PARAMETERs.

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-3320,3337-3474
/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: 180.4 KB

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

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

Context Navigation

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

Download in other formats: