Home

Context Navigation

biometeorology_mod.f90 @ 3615

Last change on this file since 3615 was 3614, checked in by raasch, 6 years ago

unused variables removed, abort renamed inifor_abort to avoid intrinsic problem in Fortran

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

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

Context Navigation

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

Download in other formats: