Home

Context Navigation

biometeorology_mod.f90 @ 3550

Last change on this file since 3550 was 3525, checked in by kanani, 6 years ago

Changes related to clean-up of biometeorology (by dom_dwd_user)

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

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

Context Navigation

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

Download in other formats: