Changeset 3525 for palm/trunk/SOURCE/biometeorology_mod.f90
 Timestamp:
 Nov 14, 2018 4:06:14 PM (4 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/biometeorology_mod.f90
r3479 r3525 27 27 !  28 28 ! $Id$ 29 ! Clean up, renaming from "biom" to "bio", summary of thermal index calculation 30 ! into one module (dom_dwd_user) 31 ! 32 ! 3479 20181101 16:00:30Z gronemeier 29 33 !  reworked check for output quantities 30 34 !  assign new palmerror numbers … … 72 76 73 77 USE basic_constants_and_equations_mod, & 74 ONLY: degc_to_k, magnus, sigma_sb 75 76 USE biometeorology_ipt_mod 77 78 USE biometeorology_pet_mod 79 80 USE biometeorology_pt_mod, & 81 ONLY: calculate_pt_static 82 83 USE biometeorology_utci_mod 78 ONLY: c_p, degc_to_k, l_v, magnus, sigma_sb 84 79 85 80 USE control_parameters, & … … 92 87 93 88 USE indices, & 94 ONLY: nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr 89 ONLY: nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg, & 90 nysg, nyng 95 91 96 92 USE kinds !< Set precision of INTEGER and REAL arrays according to PALM … … 112 108 ! Declare all global variables within the module (alphabetical order) 113 109 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tmrt_grid !< tmrt results (Â°C) 114 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: p t_grid!< PT results (Â°C)110 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: perct !< PT results (Â°C) 115 111 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: utci_grid !< UTCI results (Â°C) 116 112 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pet_grid !< PET results (Â°C) 117 ! Grids for averaged thermal indices 118 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt_av_grid !< PT results (aver. input) (Â°C) 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) 119 116 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: utci_av_grid !< UTCI results (aver. input) (Â°C) 120 117 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pet_av_grid !< PET results (aver. input) (Â°C) 121 ! Grids for radiation_model 122 REAL(wp), DIMENSION(:), ALLOCATABLE :: biom_mrt !< biometeorology mean radiant temperature for each MRT box 123 REAL(wp), DIMENSION(:), ALLOCATABLE :: biom_mrt_av !< time average mean 124 125 INTEGER( iwp ) :: biom_cell_level !< cell level biom calculates for 126 REAL ( wp ) :: biom_output_height !< height output is calculated in m 127 REAL ( wp ) :: time_biom_results !< the time, the last set of biom results have been calculated for 118 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 128 123 REAL ( wp ), PARAMETER :: human_absorb = 0.7_wp !< SW absorbtivity of a human body (Fanger 1972) 129 124 REAL ( wp ), PARAMETER :: human_emiss = 0.97_wp !< LW emissivity of a human body after (Fanger 1972) 130 125 ! 131 126 132 LOGICAL :: aver_pt = .FALSE. !< switch: do pt averaging in this module? (if .FALSE. this is done globally) 133 LOGICAL :: aver_q = .FALSE. !< switch: do e averaging in this module? 134 LOGICAL :: aver_u = .FALSE. !< switch: do u averaging in this module? 135 LOGICAL :: aver_v = .FALSE. !< switch: do v averaging in this module? 136 LOGICAL :: aver_w = .FALSE. !< switch: do w averaging in this module? 137 LOGICAL :: average_trigger_pt = .FALSE. !< update averaged input on call to biom_pt? 138 LOGICAL :: average_trigger_utci = .FALSE. !< update averaged input on call to biom_utci? 139 LOGICAL :: average_trigger_pet = .FALSE. !< update averaged input on call to biom_pet? 140 141 LOGICAL :: mrt_from_rtm = .TRUE. !< switch: use mrt calculated by RTM for calculation of thermal indices 142 LOGICAL :: biom_pt = .TRUE. !< Turn index PT (instant. input) on or off 143 LOGICAL :: biom_pt_av = .TRUE. !< Turn index PT (averaged input) on or off 144 LOGICAL :: biom_pet = .TRUE. !< Turn index PET (instant. input) on or off 145 LOGICAL :: biom_pet_av = .TRUE. !< Turn index PET (averaged input) on or off 146 LOGICAL :: biom_utci = .TRUE. !< Turn index UTCI (instant. input) on or off 147 LOGICAL :: biom_utci_av = .TRUE. !< Turn index UTCI (averaged input) on or off 148 149 ! Add INTERFACES that must be available to other modules (alphabetical order) 150 151 PUBLIC biom_3d_data_averaging, biom_check_data_output, & 152 biom_calculate_static_grid, biom_calc_ipt, & 153 biom_check_parameters, biom_data_output_3d, biom_data_output_2d, & 154 biom_define_netcdf_grid, biom_determine_input_at, biom_header, biom_init, & 155 biom_init_arrays, biom_parin, biom_pt, biom_pt_av, biom_pet, biom_pet_av, & 156 biom_utci, biom_utci_av, time_biom_results 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? 136 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 143 144 145 ! 146 ! INTERFACES that must be available to other modules (alphabetical order) 147 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 157 154 158 155 ! … … 160 157 ! 161 158 ! 3D averaging for HTCM _INPUT_ variables 162 INTERFACE biom_3d_data_averaging 163 MODULE PROCEDURE biom_3d_data_averaging 164 END INTERFACE biom_3d_data_averaging 159 INTERFACE bio_3d_data_averaging 160 MODULE PROCEDURE bio_3d_data_averaging 161 END INTERFACE bio_3d_data_averaging 162 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 165 167 166 168 ! Calculate static thermal indices PT, UTCI and/or PET 167 INTERFACE bio m_calculate_static_grid168 MODULE PROCEDURE bio m_calculate_static_grid169 END INTERFACE bio m_calculate_static_grid169 INTERFACE bio_calculate_thermal_index_maps 170 MODULE PROCEDURE bio_calculate_thermal_index_maps 171 END INTERFACE bio_calculate_thermal_index_maps 170 172 171 173 ! Calculate the dynamic index iPT (to be caled by the agent model) 172 INTERFACE bio m_calc_ipt173 MODULE PROCEDURE bio m_calc_ipt174 END INTERFACE bio m_calc_ipt174 INTERFACE bio_calc_ipt 175 MODULE PROCEDURE bio_calc_ipt 176 END INTERFACE bio_calc_ipt 175 177 176 178 ! Data output checks for 2D/3D data to be done in check_parameters 177 INTERFACE bio m_check_data_output178 MODULE PROCEDURE bio m_check_data_output179 END INTERFACE bio m_check_data_output179 INTERFACE bio_check_data_output 180 MODULE PROCEDURE bio_check_data_output 181 END INTERFACE bio_check_data_output 180 182 181 183 ! Input parameter checks to be done in check_parameters 182 INTERFACE bio m_check_parameters183 MODULE PROCEDURE bio m_check_parameters184 END INTERFACE bio m_check_parameters184 INTERFACE bio_check_parameters 185 MODULE PROCEDURE bio_check_parameters 186 END INTERFACE bio_check_parameters 185 187 186 188 ! Data output of 2D quantities 187 INTERFACE bio m_data_output_2d188 MODULE PROCEDURE bio m_data_output_2d189 END INTERFACE bio m_data_output_2d189 INTERFACE bio_data_output_2d 190 MODULE PROCEDURE bio_data_output_2d 191 END INTERFACE bio_data_output_2d 190 192 191 193 ! no 3D data, thus, no averaging of 3D data, removed 192 INTERFACE bio m_data_output_3d193 MODULE PROCEDURE bio m_data_output_3d194 END INTERFACE bio m_data_output_3d194 INTERFACE bio_data_output_3d 195 MODULE PROCEDURE bio_data_output_3d 196 END INTERFACE bio_data_output_3d 195 197 196 198 ! Definition of data output quantities 197 INTERFACE bio m_define_netcdf_grid198 MODULE PROCEDURE bio m_define_netcdf_grid199 END INTERFACE bio m_define_netcdf_grid199 INTERFACE bio_define_netcdf_grid 200 MODULE PROCEDURE bio_define_netcdf_grid 201 END INTERFACE bio_define_netcdf_grid 200 202 201 203 ! Obtains all relevant input values to estimate local thermal comfort/stress 202 INTERFACE bio m_determine_input_at203 MODULE PROCEDURE bio m_determine_input_at204 END INTERFACE bio m_determine_input_at204 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 205 207 206 208 ! Output of information to the header file 207 INTERFACE bio m_header208 MODULE PROCEDURE bio m_header209 END INTERFACE bio m_header209 INTERFACE bio_header 210 MODULE PROCEDURE bio_header 211 END INTERFACE bio_header 210 212 211 213 ! Initialization actions 212 INTERFACE bio m_init213 MODULE PROCEDURE bio m_init214 END INTERFACE bio m_init214 INTERFACE bio_init 215 MODULE PROCEDURE bio_init 216 END INTERFACE bio_init 215 217 216 218 ! Initialization of arrays 217 INTERFACE bio m_init_arrays218 MODULE PROCEDURE bio m_init_arrays219 END INTERFACE bio m_init_arrays219 INTERFACE bio_init_arrays 220 MODULE PROCEDURE bio_init_arrays 221 END INTERFACE bio_init_arrays 220 222 221 223 ! Reading of NAMELIST parameters 222 INTERFACE bio m_parin223 MODULE PROCEDURE bio m_parin224 END INTERFACE bio m_parin225 226 224 INTERFACE bio_parin 225 MODULE PROCEDURE bio_parin 226 END INTERFACE bio_parin 227 228 227 229 CONTAINS 228 229 230 231 230 232 !! 231 233 ! Description: … … 234 236 !> the array necessary for storing the average. 235 237 !! 236 SUBROUTINE bio m_3d_data_averaging( mode, variable )238 SUBROUTINE bio_3d_data_averaging( mode, variable ) 237 239 238 240 IMPLICIT NONE 239 241 240 CHARACTER (LEN=*) :: mode !< 241 CHARACTER (LEN=*) :: variable !< 242 243 INTEGER(iwp) :: i !<244 INTEGER(iwp) :: j !<245 INTEGER(iwp) :: k !<242 CHARACTER (LEN=*) :: mode !< averaging mode: allocate, sum, or average 243 CHARACTER (LEN=*) :: variable !< The variable in question 244 245 INTEGER(iwp) :: i !< Running index, xdir 246 INTEGER(iwp) :: j !< Running index, ydir 247 INTEGER(iwp) :: k !< Running index, zdir 246 248 247 249 … … 250 252 SELECT CASE ( TRIM( variable ) ) 251 253 252 CASE ( 'bio m_mrt' )253 IF ( .NOT. ALLOCATED( biom_mrt_av) ) THEN254 ALLOCATE( biom_mrt_av(nmrtbl) )254 CASE ( 'bio_mrt' ) 255 IF ( .NOT. ALLOCATED( mrt_av_grid ) ) THEN 256 ALLOCATE( mrt_av_grid(nmrtbl) ) 255 257 ENDIF 256 biom_mrt_av= 0.0_wp257 258 CASE ( 'bio m_pt', 'biom_utci', 'biom_pet' )258 mrt_av_grid = 0.0_wp 259 260 CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' ) 259 261 260 262 ! Indices in unknown order as depending on input file, determine 261 263 ! first index to average und update only once 262 IF ( .NOT. average_trigger_p t .AND. .NOT. average_trigger_utci&264 IF ( .NOT. average_trigger_perct .AND. .NOT. average_trigger_utci & 263 265 .AND. .NOT. average_trigger_pet ) THEN 264 IF ( TRIM( variable ) == 'bio m_pt' ) THEN265 average_trigger_p t = .TRUE.266 IF ( TRIM( variable ) == 'bio_perct*' ) THEN 267 average_trigger_perct = .TRUE. 266 268 ENDIF 267 IF ( TRIM( variable ) == 'bio m_utci' ) THEN269 IF ( TRIM( variable ) == 'bio_utci*' ) THEN 268 270 average_trigger_utci = .TRUE. 269 271 ENDIF 270 IF ( TRIM( variable ) == 'bio m_pet' ) THEN272 IF ( TRIM( variable ) == 'bio_pet*' ) THEN 271 273 average_trigger_pet = .TRUE. 272 274 ENDIF … … 274 276 275 277 ! Only continue if updateindex 276 IF ( average_trigger_p t .AND. TRIM( variable ) /= 'biom_pt')RETURN277 IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio m_utci')RETURN278 IF ( average_trigger_pet .AND. TRIM( variable ) /= 'bio m_pet')RETURN278 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 279 281 280 282 ! Set averaging switch to .TRUE. if not done by other module before! 281 283 IF ( .NOT. ALLOCATED( pt_av ) ) THEN 282 284 ALLOCATE( pt_av(nzb:nzt+1,nys:nyn,nxl:nxr) ) 283 aver_p t = .TRUE.285 aver_perct = .TRUE. 284 286 ENDIF 285 287 pt_av = 0.0_wp … … 292 294 293 295 IF ( .NOT. ALLOCATED( u_av ) ) THEN 294 ALLOCATE( u_av(nzb:nzt+1,nys :nyn,nxl:nxr) )296 ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 295 297 aver_u = .TRUE. 296 298 ENDIF … … 298 300 299 301 IF ( .NOT. ALLOCATED( v_av ) ) THEN 300 ALLOCATE( v_av(nzb:nzt+1,nys :nyn,nxl:nxr) )302 ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 301 303 aver_v = .TRUE. 302 304 ENDIF … … 304 306 305 307 IF ( .NOT. ALLOCATED( w_av ) ) THEN 306 ALLOCATE( w_av(nzb:nzt+1,nys :nyn,nxl:nxr) )308 ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 307 309 aver_w = .TRUE. 308 310 ENDIF … … 318 320 SELECT CASE ( TRIM( variable ) ) 319 321 320 CASE ( 'biom_mrt' ) 321 IF ( ALLOCATED( biom_mrt_av ) ) THEN 322 323 IF ( nmrtbl > 0 ) THEN 324 IF ( mrt_include_sw ) THEN 325 biom_mrt_av(:) = biom_mrt_av(:) + & 326 ((human_absorb*mrtinsw(:) + human_emiss*mrtinlw(:)) & 327 / (human_emiss*sigma_sb)) ** .25_wp  degc_to_k 328 ELSE 329 biom_mrt_av(:) = biom_mrt_av(:) + & 330 (human_emiss * mrtinlw(:) / sigma_sb) ** .25_wp & 331  degc_to_k 332 ENDIF 333 ENDIF 322 CASE ( 'bio_mrt' ) 323 IF ( ALLOCATED( mrt_av_grid ) ) THEN 324 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 334 334 ENDIF 335 335 336 CASE ( 'bio m_pt', 'biom_utci', 'biom_pet' )336 CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' ) 337 337 338 338 ! Only continue if updateindex 339 IF ( average_trigger_p t .AND. TRIM( variable ) /= 'biom_pt')&339 IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') & 340 340 RETURN 341 IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio m_utci') &341 IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*') & 342 342 RETURN 343 IF ( average_trigger_pet .AND. TRIM( variable ) /= 'bio m_pet') &343 IF ( average_trigger_pet .AND. TRIM( variable ) /= 'bio_pet*') & 344 344 RETURN 345 345 346 IF ( ALLOCATED( pt_av ) .AND. aver_p t ) THEN346 IF ( ALLOCATED( pt_av ) .AND. aver_perct ) THEN 347 347 DO i = nxl, nxr 348 348 DO j = nys, nyn … … 365 365 366 366 IF ( ALLOCATED( u_av ) .AND. aver_u ) THEN 367 DO i = nxl , nxr368 DO j = nys , nyn367 DO i = nxlg, nxrg !< yes, ghost points are required here! 368 DO j = nysg, nyng 369 369 DO k = nzb, nzt+1 370 370 u_av(k,j,i) = u_av(k,j,i) + u(k,j,i) … … 375 375 376 376 IF ( ALLOCATED( v_av ) .AND. aver_v ) THEN 377 DO i = nxl , nxr378 DO j = nys , nyn377 DO i = nxlg, nxrg !< yes, ghost points are required here! 378 DO j = nysg, nyng 379 379 DO k = nzb, nzt+1 380 380 v_av(k,j,i) = v_av(k,j,i) + v(k,j,i) … … 385 385 386 386 IF ( ALLOCATED( w_av ) .AND. aver_w ) THEN 387 DO i = nxl , nxr388 DO j = nys , nyn387 DO i = nxlg, nxrg !< yes, ghost points are required here! 388 DO j = nysg, nyng 389 389 DO k = nzb, nzt+1 390 390 w_av(k,j,i) = w_av(k,j,i) + w(k,j,i) … … 403 403 SELECT CASE ( TRIM( variable ) ) 404 404 405 CASE ( 'bio m_mrt' )406 IF ( ALLOCATED( biom_mrt_av) ) THEN407 biom_mrt_av(:) = biom_mrt_av(:) / REAL( average_count_3d, KIND=wp )405 CASE ( 'bio_mrt' ) 406 IF ( ALLOCATED( mrt_av_grid ) ) THEN 407 mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp ) 408 408 ENDIF 409 409 410 CASE ( 'bio m_pt', 'biom_utci', 'biom_pet' )410 CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' ) 411 411 412 412 ! Only continue if update index 413 IF ( average_trigger_p t .AND. TRIM( variable ) /= 'biom_pt')&413 IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') & 414 414 RETURN 415 IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio m_utci') &415 IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*') & 416 416 RETURN 417 IF ( average_trigger_pet .AND. TRIM( variable ) /= 'bio m_pet') &417 IF ( average_trigger_pet .AND. TRIM( variable ) /= 'bio_pet*') & 418 418 RETURN 419 419 420 IF ( ALLOCATED( pt_av ) .AND. aver_p t ) THEN420 IF ( ALLOCATED( pt_av ) .AND. aver_perct ) THEN 421 421 DO i = nxl, nxr 422 422 DO j = nys, nyn … … 439 439 440 440 IF ( ALLOCATED( u_av ) .AND. aver_u ) THEN 441 DO i = nxl , nxr442 DO j = nys , nyn441 DO i = nxlg, nxrg !< yes, ghost points are required here! 442 DO j = nysg, nyng 443 443 DO k = nzb, nzt+1 444 444 u_av(k,j,i) = u_av(k,j,i) / REAL( average_count_3d, KIND=wp ) … … 449 449 450 450 IF ( ALLOCATED( v_av ) .AND. aver_v ) THEN 451 DO i = nxl , nxr452 DO j = nys , nyn451 DO i = nxlg, nxrg 452 DO j = nysg, nyng 453 453 DO k = nzb, nzt+1 454 454 v_av(k,j,i) = v_av(k,j,i) / REAL( average_count_3d, KIND=wp ) … … 459 459 460 460 IF ( ALLOCATED( w_av ) .AND. aver_w ) THEN 461 DO i = nxl , nxr462 DO j = nys , nyn461 DO i = nxlg, nxrg 462 DO j = nysg, nyng 463 463 DO k = nzb, nzt+1 464 464 w_av(k,j,i) = w_av(k,j,i) / REAL( average_count_3d, KIND=wp ) … … 469 469 470 470 ! Udate thermal indices with derived averages 471 CALL bio m_calculate_static_grid( .TRUE. )471 CALL bio_calculate_thermal_index_maps ( .TRUE. ) 472 472 473 473 END SELECT … … 476 476 477 477 478 END SUBROUTINE bio m_3d_data_averaging478 END SUBROUTINE bio_3d_data_averaging 479 479 480 480 … … 485 485 !> Check data output for biometeorology model 486 486 !! 487 SUBROUTINE bio m_check_data_output( var, unit )487 SUBROUTINE bio_check_data_output( var, unit ) 488 488 489 489 USE control_parameters, & … … 498 498 SELECT CASE ( TRIM( var ) ) 499 499 500 CASE( 'bio m_tmrt', 'biom_mrt', 'biom_pet', 'biom_pt', 'biom_utci' )500 CASE( 'bio_mrt', 'bio_pet*', 'bio_perct*', 'bio_utci*' ) 501 501 unit = 'degree_C' 502 502 … … 509 509 ! 510 510 ! Futher checks if output belongs to biometeorology 511 IF ( .NOT. biometeorology .OR. .NOT.radiation ) THEN512 message_string = 'output of "' // TRIM( var ) // '" req' //&513 'uires biometeorology = .TRUE. and radiation = .TRUE. ' &514 // 'in initialisation_parameters'515 CALL message( 'biom_check_data_output', 'PA0561', 1, 2, 0, 6, 0 )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' 516 516 ENDIF 517 517 IF ( mrt_nlevels == 0 ) THEN 518 message_string = 'output of "' // TRIM( var ) // '" require' &518 message_string = 'output of "' // TRIM( var ) // '" require' & 519 519 // 's mrt_nlevels > 0' 520 CALL message( 'biom_check_data_output', 'PA0562', 1, 2, 0, 6, 0 ) 520 CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 ) 521 unit = 'illegal' 521 522 ENDIF 522 523 523 ENDIF 524 525 END SUBROUTINE biom_check_data_output 524 525 ENDIF 526 527 END SUBROUTINE bio_check_data_output 526 528 527 529 !! … … 531 533 !> check_parameters 1302 532 534 !! 533 SUBROUTINE bio m_check_parameters535 SUBROUTINE bio_check_parameters 534 536 535 537 USE control_parameters, & … … 541 543 542 544 ! Disabled as long as radiation model not available 543 IF ( .NOT. radiation ) THEN544 message_string = 'The human thermal comfort module does require ' // &545 'radiation information in terms of the mean ' // &546 'radiant temperature, but radiation is not enabled!'547 CALL message( 'check_parameters', 'PAHU01', 0, 0, 0, 6, 0 )548 ENDIF549 545 550 546 IF ( .NOT. humidity ) THEN 551 message_string = 'The human thermal comfort module does require '// &547 message_string = 'The estimation of thermal comfort requires ' // & 552 548 'air humidity information, but humidity module ' // & 553 549 'is disabled!' 554 CALL message( 'check_parameters', 'PA HU02', 0, 0, 0, 6, 0 )550 CALL message( 'check_parameters', 'PA0561', 0, 0, 0, 6, 0 ) 555 551 ENDIF 556 552 557 END SUBROUTINE biom_check_parameters 558 559 560 !! 561 ! 562 ! Description: 563 !  564 !> Subroutine defining 3D output variables (dummy, always 2d!) 565 !> data_output_3d 709ff 566 !! 567 SUBROUTINE biom_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do ) 568 569 USE indices 570 571 USE kinds 572 573 574 IMPLICIT NONE 575 576 ! Input variables 577 CHARACTER (LEN=*), INTENT(IN) :: variable !< Char identifier to select var for output 578 INTEGER(iwp), INTENT(IN) :: av !< Use averaged data? 0 = no, 1 = yes? 579 INTEGER(iwp), INTENT(IN) :: nzb_do !< Unused. 2D. nz bottom to nz top 580 INTEGER(iwp), INTENT(IN) :: nzt_do !< Unused. 581 582 ! Output variables 583 LOGICAL, INTENT(OUT) :: found !< Output found? 584 REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< Temp. result grid to return 585 586 ! Internal variables 587 INTEGER(iwp) :: l !< Running index, radiation grid 588 INTEGER(iwp) :: i !< Running index, xdir 589 INTEGER(iwp) :: j !< Running index, ydir 590 INTEGER(iwp) :: k !< Running index, zdir 591 592 CHARACTER (LEN=:), allocatable :: variable_short !< Trimmed version of char variable 593 594 REAL(wp), PARAMETER :: fill_value = 999._wp 595 REAL(wp) :: mrt !< Buffer for mean radiant temperature 596 597 found = .TRUE. 598 variable_short = TRIM( variable ) 599 600 IF ( variable_short(1:4) /= 'biom' ) THEN 601 ! Nothing to do, set found to FALSE and return immediatelly 602 found = .FALSE. 603 RETURN 604 ENDIF 605 606 SELECT CASE ( variable_short ) 607 608 CASE ( 'biom_mrt' ) 609 610 local_pf = REAL( fill_value, KIND = wp ) 611 DO l = 1, nmrtbl 612 i = mrtbl(ix,l) 613 j = mrtbl(iy,l) 614 k = mrtbl(iz,l) 615 IF ( mrt_include_sw ) THEN 616 mrt = ((human_absorb*mrtinsw(l) + human_emiss*mrtinlw(l)) & 617 / (human_emiss*sigma_sb)) ** .25_wp 618 ELSE 619 mrt = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp 620 ENDIF 621 local_pf(i,j,k) = mrt 622 ENDDO 623 624 CASE ( 'biom_tmrt' ) ! 2darray 625 DO i = nxl, nxr 626 DO j = nys, nyn 627 local_pf(i,j,nzb_do) = REAL( tmrt_grid(j,i), sp ) 628 ENDDO 629 ENDDO 630 631 CASE ( 'biom_pt' ) ! 2darray 632 IF ( av == 0 ) THEN 633 DO i = nxl, nxr 634 DO j = nys, nyn 635 local_pf(i,j,nzb_do) = REAL( pt_grid(j,i), sp ) 636 ENDDO 637 ENDDO 638 ELSE 639 DO i = nxl, nxr 640 DO j = nys, nyn 641 local_pf(i,j,nzb_do) = REAL( pt_av_grid(j,i), sp ) 642 ENDDO 643 ENDDO 644 END IF 645 646 CASE ( 'biom_utci' ) ! 2darray 647 IF ( av == 0 ) THEN 648 DO i = nxl, nxr 649 DO j = nys, nyn 650 local_pf(i,j,nzb_do) = REAL( utci_grid(j,i), sp ) 651 ENDDO 652 ENDDO 653 ELSE 654 DO i = nxl, nxr 655 DO j = nys, nyn 656 local_pf(i,j,nzb_do) = REAL( utci_av_grid(j,i), sp ) 657 ENDDO 658 ENDDO 659 END IF 660 661 CASE ( 'biom_pet' ) ! 2darray 662 IF ( av == 0 ) THEN 663 DO i = nxl, nxr 664 DO j = nys, nyn 665 local_pf(i,j,nzb_do) = REAL( pet_grid(j,i), sp ) 666 ENDDO 667 ENDDO 668 ELSE 669 DO i = nxl, nxr 670 DO j = nys, nyn 671 local_pf(i,j,nzb_do) = REAL( pet_av_grid(j,i), sp ) 672 ENDDO 673 ENDDO 674 END IF 675 676 CASE DEFAULT 677 found = .FALSE. 678 679 END SELECT 680 681 END SUBROUTINE biom_data_output_3d 553 END SUBROUTINE bio_check_parameters 554 682 555 683 556 !! … … 688 561 !> data_output_2d 1188ff 689 562 !! 690 SUBROUTINE bio m_data_output_2d( av, variable, found, grid, local_pf,&563 SUBROUTINE bio_data_output_2d( av, variable, found, grid, local_pf, & 691 564 two_d, nzb_do, nzt_do, fill_value ) 692 565 … … 717 590 INTEGER(iwp) :: j !< Running index, ydir 718 591 INTEGER(iwp) :: k !< Running index, zdir 719 720 721 found = .TRUE. 592 INTEGER(iwp) :: l !< Running index, radiation grid 593 594 722 595 variable_short = TRIM( variable ) 723 IF ( variable_short(1:4) == 'biom' ) THEN 724 two_d = .TRUE. 725 grid = 'zu1' 726 ELSE 596 IF ( variable_short(1:4) /= 'bio_' ) THEN 727 597 found = .FALSE. 728 598 grid = 'none' 729 RETURN730 ENDIF 731 599 ENDIF 600 601 found = .TRUE. 732 602 local_pf = fill_value 733 603 … … 735 605 736 606 737 CASE ( 'biom_tmrt_xy' ) ! 2darray 738 DO i = nxl, nxr 739 DO j = nys, nyn 740 local_pf(i,j,1) = tmrt_grid(j,i) 741 ENDDO 742 ENDDO 743 744 745 CASE ( 'biom_pt_xy' ) ! 2darray 607 CASE ( 'bio_mrt_xy' ) 608 grid = 'zu1' 609 two_d = .FALSE. !< can be calculated for several levels 610 local_pf = REAL( fill_value, KIND = wp ) 611 DO l = 1, nmrtbl 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 629 ENDDO 630 631 632 CASE ( 'bio_perct*_xy' ) ! 2darray 633 grid = 'zu1' 634 two_d = .TRUE. 746 635 IF ( av == 0 ) THEN 747 636 DO i = nxl, nxr 748 637 DO j = nys, nyn 749 local_pf(i,j,nzb+1) = p t_grid(j,i)638 local_pf(i,j,nzb+1) = perct(j,i) 750 639 ENDDO 751 640 ENDDO … … 753 642 DO i = nxl, nxr 754 643 DO j = nys, nyn 755 local_pf(i,j,nzb+1) = p t_av_grid(j,i)644 local_pf(i,j,nzb+1) = perct_av(j,i) 756 645 ENDDO 757 646 ENDDO … … 759 648 760 649 761 CASE ( 'biom_utci_xy' ) ! 2darray 650 CASE ( 'bio_utci*_xy' ) ! 2darray 651 grid = 'zu1' 652 two_d = .TRUE. 762 653 IF ( av == 0 ) THEN 763 654 DO i = nxl, nxr … … 775 666 776 667 777 CASE ( 'biom_pet_xy' ) ! 2darray 668 CASE ( 'bio_pet*_xy' ) ! 2darray 669 grid = 'zu1' 670 two_d = .TRUE. 778 671 IF ( av == 0 ) THEN 779 672 DO i = nxl, nxr … … 798 691 799 692 800 END SUBROUTINE biom_data_output_2d 801 693 END SUBROUTINE bio_data_output_2d 694 695 696 !! 697 ! 698 ! Description: 699 !  700 !> Subroutine defining 3D output variables (dummy, always 2d!) 701 !> data_output_3d 709ff 702 !! 703 SUBROUTINE bio_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do ) 704 705 USE indices 706 707 USE kinds 708 709 710 IMPLICIT NONE 711 712 ! Input variables 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. 717 718 ! Output variables 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 721 722 ! Internal variables 723 INTEGER(iwp) :: l !< Running index, radiation grid 724 INTEGER(iwp) :: i !< Running index, xdir 725 INTEGER(iwp) :: j !< Running index, ydir 726 INTEGER(iwp) :: k !< Running index, zdir 727 728 CHARACTER (LEN=:), allocatable :: variable_short !< Trimmed version of char variable 729 730 REAL(wp), PARAMETER :: fill_value = 999._wp 731 REAL(wp) :: mrt !< Buffer for mean radiant temperature 732 733 found = .TRUE. 734 variable_short = TRIM( variable ) 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 740 ENDIF 741 742 SELECT CASE ( variable_short ) 743 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 764 765 CASE DEFAULT 766 found = .FALSE. 767 768 END SELECT 769 770 END SUBROUTINE bio_data_output_3d 802 771 803 772 !! … … 808 777 !> netcdf_interface_mod 918ff 809 778 !! 810 SUBROUTINE bio m_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )779 SUBROUTINE bio_define_netcdf_grid( var, found, grid_x, grid_y, grid_z ) 811 780 812 781 IMPLICIT NONE … … 836 805 is2d = ( var(l1:l) == 'xy' ) 837 806 838 839 IF ( var(1:4) == 'biom' ) THEN 807 IF ( var(1:4) == 'bio_' ) THEN 840 808 found = .TRUE. 841 809 grid_x = 'x' … … 845 813 ENDIF 846 814 847 END SUBROUTINE bio m_define_netcdf_grid815 END SUBROUTINE bio_define_netcdf_grid 848 816 849 817 !! … … 853 821 !> header 982 854 822 !! 855 SUBROUTINE bio m_header( io )823 SUBROUTINE bio_header( io ) 856 824 857 825 IMPLICIT NONE … … 863 831 CHARACTER (LEN=86) :: output_height_chr !< String for output height 864 832 865 WRITE( output_height_chr, '(F8.1,7X)' ) bio m_output_height833 WRITE( output_height_chr, '(F8.1,7X)' ) bio_output_height 866 834 ! 867 835 ! Write biom header 868 836 WRITE( io, 1 ) 869 837 WRITE( io, 2 ) TRIM( output_height_chr ) 870 WRITE( io, 3 ) TRIM( ACHAR( bio m_cell_level ) )838 WRITE( io, 3 ) TRIM( ACHAR( bio_cell_level ) ) 871 839 872 840 1 FORMAT (//' Human thermal comfort module information:'/ & … … 875 843 3 FORMAT (' > This corresponds to cell level : ', A ) 876 844 877 END SUBROUTINE bio m_header845 END SUBROUTINE bio_header 878 846 879 847 … … 884 852 !> init_3d_model 1987ff 885 853 !! 886 SUBROUTINE bio m_init887 888 USE control_parameters, &854 SUBROUTINE bio_init 855 856 USE control_parameters, & 889 857 ONLY: message_string 890 858 … … 899 867 ! (gravimetric center of sample human) 900 868 901 time_bio m_results = 0.0_wp902 bio m_cell_level = 0_iwp903 bio m_output_height = 0.5_wp * dz(1)869 time_bio_results = 0.0_wp 870 bio_cell_level = 0_iwp 871 bio_output_height = 0.5_wp * dz(1) 904 872 height = 0.0_wp 905 873 906 bio m_cell_level = INT ( 1.099_wp / dz(1) )907 bio m_output_height = biom_output_height + biom_cell_level * dz(1)908 909 IF ( .NOT. radiation_interactions .AND. mrt_from_rtm) THEN910 message_string = 'The mrt _from_rtm switch require' // &874 bio_cell_level = INT ( 1.099_wp / dz(1) ) 875 bio_output_height = bio_output_height + bio_cell_level * dz(1) 876 877 IF ( .NOT. radiation_interactions ) THEN 878 message_string = 'The mrt calculation requires ' // & 911 879 'enabled radiation_interactions but it ' // & 912 'is disabled! Set mrt_from_rtm to .F.'880 'is disabled!' 913 881 CALL message( 'check_parameters', 'PAHU03', 0, 0, 1, 6, 0 ) 914 mrt_from_rtm = .FALSE. 915 ENDIF 916 917 END SUBROUTINE biom_init 882 ENDIF 883 884 END SUBROUTINE bio_init 918 885 919 886 … … 924 891 !> init_3d_model 1047ff 925 892 !! 926 SUBROUTINE bio m_init_arrays893 SUBROUTINE bio_init_arrays 927 894 928 895 IMPLICIT NONE … … 934 901 ENDIF 935 902 936 IF ( bio m_pt ) THEN937 IF ( .NOT. ALLOCATED( p t_grid) ) THEN938 ALLOCATE( p t_grid(nys:nyn,nxl:nxr) )903 IF ( bio_perct ) THEN 904 IF ( .NOT. ALLOCATED( perct ) ) THEN 905 ALLOCATE( perct (nys:nyn,nxl:nxr) ) 939 906 ENDIF 940 907 ENDIF 941 908 942 IF ( bio m_utci ) THEN909 IF ( bio_utci ) THEN 943 910 IF ( .NOT. ALLOCATED( utci_grid ) ) THEN 944 911 ALLOCATE( utci_grid (nys:nyn,nxl:nxr) ) … … 946 913 ENDIF 947 914 948 IF ( bio m_pet ) THEN915 IF ( bio_pet ) THEN 949 916 IF ( .NOT. ALLOCATED( pet_grid ) ) THEN 950 917 ALLOCATE( pet_grid (nys:nyn,nxl:nxr) ) … … 952 919 END IF 953 920 954 IF ( bio m_pt_av ) THEN955 IF ( .NOT. ALLOCATED( p t_av_grid) ) THEN956 ALLOCATE( p t_av_grid(nys:nyn,nxl:nxr) )921 IF ( bio_perct_av ) THEN 922 IF ( .NOT. ALLOCATED( perct_av ) ) THEN 923 ALLOCATE( perct_av (nys:nyn,nxl:nxr) ) 957 924 ENDIF 958 925 ENDIF 959 926 960 IF ( bio m_utci_av ) THEN927 IF ( bio_utci_av ) THEN 961 928 IF ( .NOT. ALLOCATED( utci_av_grid ) ) THEN 962 929 ALLOCATE( utci_av_grid (nys:nyn,nxl:nxr) ) … … 964 931 ENDIF 965 932 966 IF ( bio m_pet_av ) THEN933 IF ( bio_pet_av ) THEN 967 934 IF ( .NOT. ALLOCATED( pet_av_grid ) ) THEN 968 935 ALLOCATE( pet_av_grid (nys:nyn,nxl:nxr) ) … … 971 938 END IF 972 939 973 END SUBROUTINE bio m_init_arrays940 END SUBROUTINE bio_init_arrays 974 941 975 942 … … 979 946 !> Parin for &biometeorology_parameters for reading biomet parameters 980 947 !! 981 SUBROUTINE bio m_parin948 SUBROUTINE bio_parin 982 949 983 950 IMPLICIT NONE … … 987 954 CHARACTER (LEN=80) :: line !< Dummy string for current line in parameter file 988 955 989 NAMELIST /biometeorology_parameters/ mrt_from_rtm, & 990 biom_pet, & 991 biom_pet_av, & 992 biom_pt, & 993 biom_pt_av, & 994 biom_utci, & 995 biom_utci_av 956 NAMELIST /biometeorology_parameters/ bio_pet, & 957 bio_pet_av, & 958 bio_perct, & 959 bio_perct_av, & 960 bio_utci, & 961 bio_utci_av 996 962 997 963 … … 1025 991 1026 992 1027 END SUBROUTINE bio m_parin993 END SUBROUTINE bio_parin 1028 994 1029 995 !! 1030 996 ! Description: 1031 997 !  1032 !> Calculates the mean radiant temperature (tmrt) based on the Sixdirections 1033 !> method according to VDI 3787 2. 1034 !! 1035 SUBROUTINE calculate_tmrt_6_directions( SW_N, SW_E, SW_S, SW_W, & 1036 SW_U, SW_D, LW_N, LW_E, LW_S, LW_W, LW_U, LW_D, tmrt ) 998 !> Calculate biometeorology MRT for all 2D grid 999 !! 1000 SUBROUTINE bio_calculate_mrt_grid ( av ) 1037 1001 1038 1002 IMPLICIT NONE 1039 1003 1040 ! Type of input of the argument list 1041 ! Short (SW_) and longwave (LW_) radiation fluxes from the six directions 1042 ! North (N), East (E), South (S), West (W), up (U) and down (D) 1043 REAL(wp), INTENT ( IN ) :: SW_N !< Sw radflux density from N (W/mÂ²) 1044 REAL(wp), INTENT ( IN ) :: SW_E !< Sw radflux density from E (W/mÂ²) 1045 REAL(wp), INTENT ( IN ) :: SW_S !< Sw radflux density from S (W/mÂ²) 1046 REAL(wp), INTENT ( IN ) :: SW_W !< Sw radflux density from W (W/mÂ²) 1047 REAL(wp), INTENT ( IN ) :: SW_U !< Sw radflux density from U (W/mÂ²) 1048 REAL(wp), INTENT ( IN ) :: SW_D !< Sw radflux density from D (W/mÂ²) 1049 REAL(wp), INTENT ( IN ) :: LW_N !< Lw radflux density from N (W/mÂ²) 1050 REAL(wp), INTENT ( IN ) :: LW_E !< Lw radflux density from E (W/mÂ²) 1051 REAL(wp), INTENT ( IN ) :: LW_S !< Lw radflux density from S (W/mÂ²) 1052 REAL(wp), INTENT ( IN ) :: LW_W !< Lw radflux density from W (W/mÂ²) 1053 REAL(wp), INTENT ( IN ) :: LW_U !< Lw radflux density from U (W/mÂ²) 1054 REAL(wp), INTENT ( IN ) :: LW_D !< Lw radflux density from D (W/mÂ²) 1055 1056 ! Type of output of the argument list 1057 REAL(wp), INTENT ( OUT ) :: tmrt !< Mean radiant temperature (Â°C) 1058 1059 ! Directional weighting factors 1060 REAL(wp), PARAMETER :: weight_h = 0.22_wp 1061 REAL(wp), PARAMETER :: weight_v = 0.06_wp 1062 1063 REAL(wp) :: nrfd !< Net radiation flux density (W/mÂ²) 1064 1065 ! Initialise 1066 nrfd = 0._wp 1067 1068 ! Compute mean radiation flux density absorbed by rotational symetric human 1069 nrfd = ( weight_h * ( human_absorb * SW_N + human_emiss * LW_N ) ) + & 1070 ( weight_h * ( human_absorb * SW_E + human_emiss * LW_E ) ) + & 1071 ( weight_h * ( human_absorb * SW_S + human_emiss * LW_S ) ) + & 1072 ( weight_h * ( human_absorb * SW_W + human_emiss * LW_W ) ) + & 1073 ( weight_v * ( human_absorb * SW_U + human_emiss * LW_U ) ) + & 1074 ( weight_v * ( human_absorb * SW_D + human_emiss * LW_D ) ) 1075 1076 ! Compute mean radiant temperature 1077 tmrt = ( nrfd / (human_emiss * sigma_sb) )**0.25_wp  degc_to_k 1078 1079 END SUBROUTINE calculate_tmrt_6_directions 1080 1081 !! 1082 ! Description: 1083 !  1084 !> Very crude approximation of mean radiant temperature based on upwards and 1085 !> downwards radiation fluxes only (other directions curently not available, 1086 !> replace as soon as possible!) 1087 !! 1088 SUBROUTINE calculate_tmrt_2_directions( sw_u, sw_d, lw_u, lw_d, ta, tmrt ) 1089 1090 IMPLICIT NONE 1091 1092 ! Type of input of the argument list 1093 REAL(wp), INTENT ( IN ) :: sw_u !< Shortwave radiation flux density from upper direction (W/mÂ²) 1094 REAL(wp), INTENT ( IN ) :: sw_d !< Shortwave radiation flux density from lower direction (W/mÂ²) 1095 REAL(wp), INTENT ( IN ) :: lw_u !< Longwave radiation flux density from upper direction (W/mÂ²) 1096 REAL(wp), INTENT ( IN ) :: lw_d !< Longwave radiation flux density from lower direction (W/mÂ²) 1097 REAL(wp), INTENT ( IN ) :: ta !< Air temperature (Â°C) 1098 1099 ! Type of output of the argument list 1100 REAL(wp), INTENT ( OUT ) :: tmrt !< mean radiant temperature, (Â°C) 1101 1102 ! Directional weighting factors and parameters 1103 REAL(wp), PARAMETER :: weight_h = 0.22_wp !< Weight for horizontal radiational gain after Fanger (1972) 1104 REAL(wp), PARAMETER :: weight_v = 0.06_wp !< Weight for vertical radiational gain after Fanger (1972) 1105 1106 ! Other internal variables 1107 REAL(wp) :: sw_in 1108 REAL(wp) :: sw_out 1109 REAL(wp) :: lw_in 1110 REAL(wp) :: lw_out 1111 REAL(wp) :: nrfd !< Net radiation flux density (W/mÂ²) 1112 REAL(wp) :: lw_air !< Longwave emission by surrounding air volume (W/mÂ²) 1113 REAL(wp) :: sw_side !< Shortwave immission from the sides (W/mÂ²) 1114 1115 INTEGER(iwp) :: no_input !< Count missing input radiation fluxes 1116 1117 ! initialise 1118 sw_in = sw_u 1119 sw_out = sw_d 1120 lw_in = lw_u 1121 lw_out = lw_d 1122 nrfd = 0._wp 1123 no_input = 0_iwp 1124 1125 ! test for missing input data 1126 IF ( sw_in <= 998._wp .OR. sw_out <= 998._wp .OR. lw_in <= 998._wp .OR. & 1127 lw_out <= 998._wp .OR. ta <= 998._wp ) THEN 1128 IF ( sw_in <= 998._wp ) THEN 1129 sw_in = 0. 1130 no_input = no_input + 1 1004 LOGICAL, INTENT(IN) :: av !< use averaged input? 1005 ! Internal variables 1006 INTEGER(iwp) :: i !< Running index, xdir, radiation coordinates 1007 INTEGER(iwp) :: j !< Running index, ydir, radiation coordinates 1008 INTEGER(iwp) :: k !< Running index, ydir, radiation coordinates 1009 INTEGER(iwp) :: l !< Running index, radiation coordinates 1010 1011 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_emiss*sigma_sb)) ** .25_wp  degc_to_k 1025 ELSE 1026 tmrt_grid(j,i) = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp & 1027  degc_to_k 1028 ENDIF 1131 1029 ENDIF 1132 IF ( sw_out <= 998._wp ) THEN 1133 sw_out = 0. 1134 no_input = no_input + 1 1135 ENDIF 1136 IF ( lw_in <= 998._wp ) THEN 1137 lw_in = 0. 1138 no_input = no_input + 1 1139 ENDIF 1140 IF ( lw_out <= 998._wp ) THEN 1141 lw_out = 0. 1142 no_input = no_input + 1 1143 ENDIF 1144 1145 ! Accept two missing radiation flux directions, fail otherwise as error might become too large 1146 IF ( ta <= 998._wp .OR. no_input >= 3 ) THEN 1147 tmrt = 999._wp 1148 RETURN 1149 ENDIF 1150 ENDIF 1151 1152 sw_side = sw_in * 0.125_wp ! distribute half of upper sw_in to the 4 sides 1153 lw_air = ( sigma_sb * 0.95_wp * ( ta + degc_to_k )**4 ) 1154 1155 ! Compute mean radiation flux density absorbed by rotational symetric human 1156 nrfd = ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) + & 1157 ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) + & 1158 ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) + & 1159 ( weight_h * ( human_absorb * sw_side + human_emiss * lw_air ) ) + & 1160 ( weight_v * ( human_absorb * (sw_in * 0.5_wp) + human_emiss * lw_in ) ) + & 1161 ( weight_v * ( human_absorb * sw_out + human_emiss * lw_out ) ) 1162 1163 ! Compute mean radiant temperature 1164 tmrt = ( nrfd / (human_emiss * sigma_sb) )**0.25_wp  degc_to_k 1165 1166 END SUBROUTINE calculate_tmrt_2_directions 1167 1168 !! 1169 ! Description: 1170 !  1171 !> Calculate static thermal indices for given meteorological conditions 1172 !! 1173 SUBROUTINE calculate_static_thermal_indices ( ta, vp, ws, pair, tmrt, & 1174 pt_static, utci_static, pet_static ) 1175 1176 IMPLICIT NONE 1177 1178 ! Input parameters 1179 REAL(wp), INTENT ( IN ) :: ta !< Air temperature (Â°C) 1180 REAL(wp), INTENT ( IN ) :: vp !< Vapour pressure (hPa) 1181 REAL(wp), INTENT ( IN ) :: ws !< Wind speed (local level) (m/s) 1182 REAL(wp), INTENT ( IN ) :: pair !< Air pressure (hPa) 1183 REAL(wp), INTENT ( IN ) :: tmrt !< Mean radiant temperature (Â°C) 1184 ! Output parameters 1185 REAL(wp), INTENT ( OUT ) :: pt_static !< Perceived temperature (Â°C) 1186 REAL(wp), INTENT ( OUT ) :: utci_static !< Universal thermal climate index (Â°C) 1187 REAL(wp), INTENT ( OUT ) :: pet_static !< Physiologically equivalent temp. (Â°C) 1188 ! Temporary field, not used here 1189 REAL(wp) :: clo !< Clothing index (no dim.) 1190 1191 clo = 999._wp 1192 1193 IF ( biom_pt ) THEN 1194 ! Estimate local perceived temperature 1195 CALL calculate_pt_static( ta, vp, ws, tmrt, pair, clo, pt_static ) 1196 ENDIF 1197 1198 IF ( biom_utci ) THEN 1199 ! Estimate local universal thermal climate index 1200 CALL calculate_utci_static( ta, vp, ws, tmrt, biom_output_height, & 1201 utci_static ) 1202 ENDIF 1203 1204 IF ( biom_pet ) THEN 1205 ! Estimate local physiologically equivalent temperature 1206 CALL calculate_pet_static( ta, vp, ws, tmrt, pair, pet_static ) 1207 ENDIF 1208 1209 END SUBROUTINE calculate_static_thermal_indices 1030 ENDDO 1031 1032 END SUBROUTINE bio_calculate_mrt_grid 1210 1033 1211 1034 … … 1215 1038 !> Calculate static thermal indices for 2D grid point i, j 1216 1039 !! 1217 SUBROUTINE bio m_determine_input_at( average_input, i, j, ta, vp, ws, pair,&1218 tmrt )1040 SUBROUTINE bio_get_thermal_index_input_ij( average_input, i, j, ta, vp, ws, & 1041 pair, tmrt ) 1219 1042 1220 1043 IMPLICIT NONE … … 1234 1057 ! Internal variables 1235 1058 INTEGER(iwp) :: k !< Running index, zdir 1059 INTEGER(iwp) :: ir !< Running index, xdir, radiation coordinates 1060 INTEGER(iwp) :: jr !< Running index, ydir, radiation coordinates 1061 INTEGER(iwp) :: kr !< Running index, ydir, radiation coordinates 1236 1062 INTEGER(iwp) :: k_wind !< Running index, zdir, wind speed only 1063 INTEGER(iwp) :: l !< Running index, radiation coordinates 1237 1064 1238 1065 REAL(wp) :: vp_sat !< Saturation vapor pressure (hPa) … … 1241 1068 ! Determine cell level closest to 1.1m above ground 1242 1069 ! by making use of truncation due to int cast 1243 k = get_topography_top_index_ji(j, i, 's') + biom_cell_level !< Vertical cell center closest to 1.1m 1070 k = get_topography_top_index_ji(j, i, 's') + bio_cell_level !< Vertical cell center closest to 1.1m 1071 1072 ! 1073 ! Avoid nonrepresentative horizontal u and v of 0.0 m/s too close to ground 1244 1074 k_wind = k 1245 IF ( k_wind < 1_iwp ) THEN ! Avoid horizontal u and v of 0.0 m/s close to ground1246 k_wind = 1_iwp1075 IF ( bio_cell_level < 1_iwp ) THEN 1076 k_wind = k + 1_iwp 1247 1077 ENDIF 1248 1078 … … 1252 1082 ta = pt_av(k, j, i)  ( 0.0098_wp * dz(1) * ( k + .5_wp ) )  degc_to_k 1253 1083 1254 vp = 0.034_wp1084 vp = 999._wp 1255 1085 IF ( humidity .AND. ALLOCATED( q_av ) ) THEN 1256 1086 vp = q_av(k, j, i) … … 1264 1094 ta = pt(k, j, i)  ( 0.0098_wp * dz(1) * ( k + .5_wp ) )  degc_to_k 1265 1095 1266 vp = q(k, j, i) 1096 vp = 999._wp 1097 IF ( humidity ) THEN 1098 vp = q(k, j, i) 1099 ENDIF 1267 1100 1268 1101 ws = ( 0.5_wp * ABS( u(k_wind, j, i) + u(k_wind, j, i+1) ) + & … … 1278 1111 ! The magnus formula is limited to temperatures up to 333.15 K to 1279 1112 ! avoid negative values of vp_sat 1280 vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) ) 1281 vp = vp * pair / ( vp + 0.622_wp ) 1282 IF ( vp > vp_sat ) vp = vp_sat 1283 1284 tmrt = ta 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 1118 1119 ! local mtr value at [i,j] 1120 tmrt = 999. !< this can be a valid result (e.g. for inside some ostacle) 1285 1121 IF ( radiation ) THEN 1286 IF ( mrt_from_rtm ) THEN 1287 tmrt = tmrt_grid(j, i) 1288 ELSE 1289 CALL calculate_tmrt_2_directions (rad_sw_in(0, j, i), & 1290 rad_sw_out(0, j, i), rad_lw_in(0, j, i), rad_lw_out(0, j, i), ta, & 1291 tmrt ) 1292 ENDIF 1293 ENDIF 1294 1295 END SUBROUTINE biom_determine_input_at 1122 ! Use MRT from RTM precalculated in tmrt_grid 1123 tmrt = tmrt_grid(j,i) 1124 ENDIF 1125 1126 END SUBROUTINE bio_get_thermal_index_input_ij 1296 1127 1297 1128 … … 1302 1133 !> time_integration.f90: 1065ff 1303 1134 !! 1304 SUBROUTINE bio m_calculate_static_grid ( average_input)1135 SUBROUTINE bio_calculate_thermal_index_maps ( av ) 1305 1136 1306 1137 IMPLICIT NONE 1307 1138 1308 1139 ! Input attributes 1309 LOGICAL, INTENT ( IN ) :: av erage_input !< Calculate based on averaged input?conditions?1140 LOGICAL, INTENT ( IN ) :: av !< Calculate based on averaged input conditions? 1310 1141 1311 1142 ! Internal variables 1312 INTEGER(iwp) :: i, j, k, l !< Running index 1313 1143 INTEGER(iwp) :: i, j !< Running index 1144 1145 REAL(wp) :: clo !< Clothing index (no dimension) 1314 1146 REAL(wp) :: ta !< Air temperature (Â°C) 1315 1147 REAL(wp) :: vp !< Vapour pressure (hPa) 1316 1148 REAL(wp) :: ws !< Wind speed (local level) (m/s) 1317 1149 REAL(wp) :: pair !< Air pressure (hPa) 1318 REAL(wp) :: tmrt_tmp !< Mean radiant temperature 1319 REAL(wp) :: pt_tmp !< Perceived temperature 1320 REAL(wp) :: utci_tmp !< Universal thermal climate index 1321 REAL(wp) :: pet_tmp !< Physiologically equivalent temperature 1322 1323 1324 CALL biom_init_arrays 1325 1326 IF ( mrt_from_rtm ) THEN 1327 tmrt_grid = 999._wp 1328 DO l = 1, nmrtbl 1329 i = mrtbl(ix,l) 1330 j = mrtbl(iy,l) 1331 k = mrtbl(iz,l) 1332 IF ( k  get_topography_top_index_ji( j, i, 's' ) == 1 ) THEN 1333 IF ( mrt_include_sw ) THEN 1334 tmrt_tmp = ((human_absorb*mrtinsw(l) + human_emiss*mrtinlw(l)) & 1335 / (human_emiss*sigma_sb)) ** .25_wp 1336 ELSE 1337 tmrt_tmp = (human_emiss*mrtinlw(l) / sigma_sb) ** .25_wp 1338 ENDIF 1339 tmrt_grid(j, i) = tmrt_tmp  degc_to_k 1340 ENDIF 1341 ENDDO 1342 ENDIF 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) 1154 1155 CALL bio_init_arrays 1156 1157 ! fill out the MRT 2D grid from appropriate source (RTM, RRTMG,...) 1158 CALL bio_calculate_mrt_grid ( av ) 1159 1343 1160 1344 1161 DO i = nxl, nxr 1345 1162 DO j = nys, nyn 1346 1163 ! Determine local input conditions 1347 CALL biom_determine_input_at ( average_input, i, j, ta, vp, ws, & 1348 pair, tmrt_tmp ) 1349 tmrt_grid(j, i) = tmrt_tmp 1350 1351 ! Only proceed if tmrt is available 1352 IF ( tmrt_tmp <= 998._wp ) THEN 1353 pt_tmp = 999._wp 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! 1167 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 1354 1172 utci_tmp = 999._wp 1355 pet_tmp = 999._wp 1356 CYCLE 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 1357 1193 END IF 1358 1194 1359 ! Calculate static thermal indices based on local tmrt 1360 CALL calculate_static_thermal_indices ( ta, vp, ws, & 1361 pair, tmrt_tmp, pt_tmp, utci_tmp, pet_tmp ) 1362 1363 IF ( average_input ) THEN 1195 1196 IF ( av ) THEN 1364 1197 ! Write results for selected averaged indices 1365 IF ( bio m_pt_av ) THEN1366 p t_av_grid(j, i) = pt_tmp1198 IF ( bio_perct_av ) THEN 1199 perct_av(j, i) = perct_tmp 1367 1200 END IF 1368 IF ( bio m_utci_av ) THEN1201 IF ( bio_utci_av ) THEN 1369 1202 utci_av_grid(j, i) = utci_tmp 1370 1203 END IF 1371 IF ( bio m_pet_av ) THEN1204 IF ( bio_pet_av ) THEN 1372 1205 pet_av_grid(j, i) = pet_tmp 1373 1206 END IF 1374 1207 ELSE 1375 1208 ! Write result for selected indices 1376 IF ( bio m_pt ) THEN1377 p t_grid(j, i) = pt_tmp1209 IF ( bio_perct ) THEN 1210 perct(j, i) = perct_tmp 1378 1211 END IF 1379 IF ( bio m_utci ) THEN1212 IF ( bio_utci ) THEN 1380 1213 utci_grid(j, i) = utci_tmp 1381 1214 END IF 1382 IF ( bio m_pet ) THEN1215 IF ( bio_pet ) THEN 1383 1216 pet_grid(j, i) = pet_tmp 1384 1217 END IF … … 1388 1221 END DO 1389 1222 1390 END SUBROUTINE bio m_calculate_static_grid1223 END SUBROUTINE bio_calculate_thermal_index_maps 1391 1224 1392 1225 !! … … 1395 1228 !> Calculate dynamic thermal indices (currently only iPT, but expandable) 1396 1229 !! 1397 SUBROUTINE bio m_calc_ipt( ta, vp, ws, pair, tmrt, dt, energy_storage,&1230 SUBROUTINE bio_calc_ipt( ta, vp, ws, pair, tmrt, dt, energy_storage, & 1398 1231 t_clo, clo, actlev, age, weight, height, work, sex, ipt ) 1399 1232 … … 1436 1269 ENDIF 1437 1270 1438 END SUBROUTINE biom_calc_ipt 1271 END SUBROUTINE bio_calc_ipt 1272 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):481494. 1284 !> doi:10.1007/s0048401104541 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 !< deltatmrt (Â°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.07562052e01 ) + & 1347 ( 2.27712343e02 ) * ta + & 1348 ( 8.06470249e04 ) * ta * ta + & 1349 ( 1.54271372e04 ) * ta * ta * ta + & 1350 ( 3.24651735e06 ) * ta * ta * ta * ta + & 1351 ( 7.32602852e08 ) * ta * ta * ta * ta * ta + & 1352 ( 1.35959073e09 ) * ta * ta * ta * ta * ta * ta + & 1353 ( 2.25836520e+00 ) * va + & 1354 ( 8.80326035e02 ) * ta * va + & 1355 ( 2.16844454e03 ) * ta * ta * va + & 1356 ( 1.53347087e05 ) * ta * ta * ta * va + & 1357 ( 5.72983704e07 ) * ta * ta * ta * ta * va + & 1358 ( 2.55090145e09 ) * ta * ta * ta * ta * ta * va + & 1359 ( 7.51269505e01 ) * va * va + & 1360 ( 4.08350271e03 ) * ta * va * va + & 1361 ( 5.21670675e05 ) * ta * ta * va * va + & 1362 ( 1.94544667e06 ) * ta * ta * ta * va * va + & 1363 ( 1.14099531e08 ) * ta * ta * ta * ta * va * va + & 1364 ( 1.58137256e01 ) * va * va * va + & 1365 ( 6.57263143e05 ) * ta * va * va * va + & 1366 ( 2.22697524e07 ) * ta * ta * va * va * va + & 1367 ( 4.16117031e08 ) * ta * ta * ta * va * va * va + & 1368 ( 1.27762753e02 ) * va * va * va * va + & 1369 ( 9.66891875e06 ) * ta * va * va * va * va + & 1370 ( 2.52785852e09 ) * ta * ta * va * va * va * va + & 1371 ( 4.56306672e04 ) * va * va * va * va * va + & 1372 ( 1.74202546e07 ) * ta * va * va * va * va * va + & 1373 ( 5.91491269e06 ) * va * va * va * va * va * va + & 1374 ( 3.98374029e01 ) * d_tmrt + & 1375 ( 1.83945314e04 ) * ta * d_tmrt + & 1376 ( 1.73754510e04 ) * ta * ta * d_tmrt + & 1377 ( 7.60781159e07 ) * ta * ta * ta * d_tmrt + & 1378 ( 3.77830287e08 ) * ta * ta * ta * ta * d_tmrt + & 1379 ( 5.43079673e10 ) * ta * ta * ta * ta * ta * d_tmrt + & 1380 ( 2.00518269e02 ) * va * d_tmrt + & 1381 ( 8.92859837e04 ) * ta * va * d_tmrt + & 1382 ( 3.45433048e06 ) * ta * ta * va * d_tmrt + & 1383 ( 3.77925774e07 ) * ta * ta * ta * va * d_tmrt + & 1384 ( 1.69699377e09 ) * ta * ta * ta * ta * va * d_tmrt + & 1385 ( 1.69992415e04 ) * va * va * d_tmrt + & 1386 ( 4.99204314e05 ) * ta * va * va * d_tmrt + & 1387 ( 2.47417178e07 ) * ta * ta * va * va * d_tmrt + & 1388 ( 1.07596466e08 ) * ta * ta * ta * va * va * d_tmrt + & 1389 ( 8.49242932e05 ) * va * va * va * d_tmrt + & 1390 ( 1.35191328e06 ) * ta * va * va * va * d_tmrt + & 1391 ( 6.21531254e09 ) * ta * ta * va * va * va * d_tmrt + & 1392 ( 4.99410301e06 ) * va * va * va * va * d_tmrt + & 1393 ( 1.89489258e08 ) * ta * va * va * va * va * d_tmrt + & 1394 ( 8.15300114e08 ) * va * va * va * va * va * d_tmrt + & 1395 ( 7.55043090e04 ) * d_tmrt * d_tmrt + & 1396 ( 5.65095215e05 ) * ta * d_tmrt * d_tmrt + & 1397 ( 4.52166564e07 ) * ta * ta * d_tmrt * d_tmrt + & 1398 ( 2.46688878e08 ) * ta * ta * ta * d_tmrt * d_tmrt + & 1399 ( 2.42674348e10 ) * ta * ta * ta * ta * d_tmrt * d_tmrt + & 1400 ( 1.54547250e04 ) * va * d_tmrt * d_tmrt + & 1401 ( 5.24110970e06 ) * ta * va * d_tmrt * d_tmrt + & 1402 ( 8.75874982e08 ) * ta * ta * va * d_tmrt * d_tmrt + & 1403 ( 1.50743064e09 ) * ta * ta * ta * va * d_tmrt * d_tmrt + & 1404 ( 1.56236307e05 ) * va * va * d_tmrt * d_tmrt + & 1405 ( 1.33895614e07 ) * ta * va * va * d_tmrt * d_tmrt + & 1406 ( 2.49709824e09 ) * ta * ta * va * va * d_tmrt * d_tmrt + & 1407 ( 6.51711721e07 ) * va * va * va * d_tmrt * d_tmrt + & 1408 ( 1.94960053e09 ) * ta * va * va * va * d_tmrt * d_tmrt + & 1409 ( 1.00361113e08 ) * va * va * va * va * d_tmrt * d_tmrt + & 1410 ( 1.21206673e05 ) * d_tmrt * d_tmrt * d_tmrt + & 1411 ( 2.18203660e07 ) * ta * d_tmrt * d_tmrt * d_tmrt + & 1412 ( 7.51269482e09 ) * ta * ta * d_tmrt * d_tmrt * d_tmrt + & 1413 ( 9.79063848e11 ) * ta * ta * ta * d_tmrt * d_tmrt * d_tmrt + & 1414 ( 1.25006734e06 ) * va * d_tmrt * d_tmrt * d_tmrt + & 1415 ( 1.81584736e09 ) * ta * va * d_tmrt * d_tmrt * d_tmrt + & 1416 ( 3.52197671e10 ) * ta * ta * va * d_tmrt * d_tmrt * d_tmrt + & 1417 ( 3.36514630e08 ) * va * va * d_tmrt * d_tmrt * d_tmrt + & 1418 ( 1.35908359e10 ) * ta * va * va * d_tmrt * d_tmrt * d_tmrt + & 1419 ( 4.17032620e10 ) * va * va * va * d_tmrt * d_tmrt * d_tmrt + & 1420 ( 1.30369025e09 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt + & 1421 ( 4.13908461e10 ) * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt + & 1422 ( 9.22652254e12 ) * ta * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt + & 1423 ( 5.08220384e09 ) * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt + & 1424 ( 2.24730961e11 ) * ta * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt + & 1425 ( 1.17139133e10 ) * va * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt + & 1426 ( 6.62154879e10 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt + & 1427 ( 4.03863260e13 ) * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt + & 1428 ( 1.95087203e12 ) * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt + & 1429 ( 4.73602469e12 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt * & 1430 d_tmrt + & 1431 ( 5.12733497e+00 ) * pa + & 1432 ( 3.12788561e01 ) * ta * pa + & 1433 ( 1.96701861e02 ) * ta * ta * pa + & 1434 ( 9.99690870e04 ) * ta * ta * ta * pa + & 1435 ( 9.51738512e06 ) * ta * ta * ta * ta * pa + & 1436 ( 4.66426341e07 ) * ta * ta * ta * ta * ta * pa + & 1437 ( 5.48050612e01 ) * va * pa + & 1438 ( 3.30552823e03 ) * ta * va * pa + & 1439 ( 1.64119440e03 ) * ta * ta * va * pa + & 1440 ( 5.16670694e06 ) * ta * ta * ta * va * pa + & 1441 ( 9.52692432e07 ) * ta * ta * ta * ta * va * pa + & 1442 ( 4.29223622e02 ) * va * va * pa + & 1443 ( 5.00845667e03 ) * ta * va * va * pa + & 1444 ( 1.00601257e06 ) * ta * ta * va * va * pa + & 1445 ( 1.81748644e06 ) * ta * ta * ta * va * va * pa + & 1446 ( 1.25813502e03 ) * va * va * va * pa + & 1447 ( 1.79330391e04 ) * ta * va * va * va * pa + & 1448 ( 2.34994441e06 ) * ta * ta * va * va * va * pa + & 1449 ( 1.29735808e04 ) * va * va * va * va * pa + & 1450 ( 1.29064870e06 ) * ta * va * va * va * va * pa + & 1451 ( 2.28558686e06 ) * va * va * va * va * va * pa + & 1452 ( 3.69476348e02 ) * d_tmrt * pa + & 1453 ( 1.62325322e03 ) * ta * d_tmrt * pa + & 1454 ( 3.14279680e05 ) * ta * ta * d_tmrt * pa + & 1455 ( 2.59835559e06 ) * ta * ta * ta * d_tmrt * pa + & 1456 ( 4.77136523e08 ) * ta * ta * ta * ta * d_tmrt * pa + & 1457 ( 8.64203390e03 ) * va * d_tmrt * pa + & 1458 ( 6.87405181e04 ) * ta * va * d_tmrt * pa + & 1459 ( 9.13863872e06 ) * ta * ta * va * d_tmrt * pa + & 1460 ( 5.15916806e07 ) * ta * ta * ta * va * d_tmrt * pa + & 1461 ( 3.59217476e05 ) * va * va * d_tmrt * pa + & 1462 ( 3.28696511e05 ) * ta * va * va * d_tmrt * pa + & 1463 ( 7.10542454e07 ) * ta * ta * va * va * d_tmrt * pa + & 1464 ( 1.24382300e05 ) * va * va * va * d_tmrt * pa + & 1465 ( 7.38584400e09 ) * ta * va * va * va * d_tmrt * pa + & 1466 ( 2.20609296e07 ) * va * va * va * va * d_tmrt * pa + & 1467 ( 7.32469180e04 ) * d_tmrt * d_tmrt * pa + & 1468 ( 1.87381964e05 ) * ta * d_tmrt * d_tmrt * pa + & 1469 ( 4.80925239e06 ) * ta * ta * d_tmrt * d_tmrt * pa + & 1470 ( 8.75492040e08 ) * ta * ta * ta * d_tmrt * d_tmrt * pa + & 1471 ( 2.77862930e05 ) * va * d_tmrt * d_tmrt * pa + & 1472 ( 5.06004592e06 ) * ta * va * d_tmrt * d_tmrt * pa + & 1473 ( 1.14325367e07 ) * ta * ta * va * d_tmrt * d_tmrt * pa + & 1474 ( 2.53016723e06 ) * va * va * d_tmrt * d_tmrt * pa + & 1475 ( 1.72857035e08 ) * ta * va * va * d_tmrt * d_tmrt * pa + & 1476 ( 3.95079398e08 ) * va * va * va * d_tmrt * d_tmrt * pa + & 1477 ( 3.59413173e07 ) * d_tmrt * d_tmrt * d_tmrt * pa + & 1478 ( 7.04388046e07 ) * ta * d_tmrt * d_tmrt * d_tmrt * pa + & 1479 ( 1.89309167e08 ) * ta * ta * d_tmrt * d_tmrt * d_tmrt * pa + & 1480 ( 4.79768731e07 ) * va * d_tmrt * d_tmrt * d_tmrt * pa + & 1481 ( 7.96079978e09 ) * ta * va * d_tmrt * d_tmrt * d_tmrt * pa + & 1482 ( 1.62897058e09 ) * va * va * d_tmrt * d_tmrt * d_tmrt * pa + & 1483 ( 3.94367674e08 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa + & 1484 ( 1.18566247e09 ) * ta * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa + & 1485 ( 3.34678041e10 ) * va * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa + & 1486 ( 1.15606447e10 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa + & 1487 ( 2.80626406e+00 ) * pa * pa + & 1488 ( 5.48712484e01 ) * ta * pa * pa + & 1489 ( 3.99428410e03 ) * ta * ta * pa * pa + & 1490 ( 9.54009191e04 ) * ta * ta * ta * pa * pa + & 1491 ( 1.93090978e05 ) * ta * ta * ta * ta * pa * pa + & 1492 ( 3.08806365e01 ) * va * pa * pa + & 1493 ( 1.16952364e02 ) * ta * va * pa * pa + & 1494 ( 4.95271903e04 ) * ta * ta * va * pa * pa + & 1495 ( 1.90710882e05 ) * ta * ta * ta * va * pa * pa + & 1496 ( 2.10787756e03 ) * va * va * pa * pa + & 1497 ( 6.98445738e04 ) * ta * va * va * pa * pa + & 1498 ( 2.30109073e05 ) * ta * ta * va * va * pa * pa + & 1499 ( 4.17856590e04 ) * va * va * va * pa * pa + & 1500 ( 1.27043871e05 ) * ta * va * va * va * pa * pa + & 1501 ( 3.04620472e06 ) * va * va * va * va * pa * pa + & 1502 ( 5.14507424e02 ) * d_tmrt * pa * pa + & 1503 ( 4.32510997e03 ) * ta * d_tmrt * pa * pa + & 1504 ( 8.99281156e05 ) * ta * ta * d_tmrt * pa * pa + & 1505 ( 7.14663943e07 ) * ta * ta * ta * d_tmrt * pa * pa + & 1506 ( 2.66016305e04 ) * va * d_tmrt * pa * pa + & 1507 ( 2.63789586e04 ) * ta * va * d_tmrt * pa * pa + & 1508 ( 7.01199003e06 ) * ta * ta * va * d_tmrt * pa * pa + & 1509 ( 1.06823306e04 ) * va * va * d_tmrt * pa * pa + & 1510 ( 3.61341136e06 ) * ta * va * va * d_tmrt * pa * pa + & 1511 ( 2.29748967e07 ) * va * va * va * d_tmrt * pa * pa + & 1512 ( 3.04788893e04 ) * d_tmrt * d_tmrt * pa * pa + & 1513 ( 6.42070836e05 ) * ta * d_tmrt * d_tmrt * pa * pa + & 1514 ( 1.16257971e06 ) * ta * ta * d_tmrt * d_tmrt * pa * pa + & 1515 ( 7.68023384e06 ) * va * d_tmrt * d_tmrt * pa * pa + & 1516 ( 5.47446896e07 ) * ta * va * d_tmrt * d_tmrt * pa * pa + & 1517 ( 3.59937910e08 ) * va * va * d_tmrt * d_tmrt * pa * pa + & 1518 ( 4.36497725e06 ) * d_tmrt * d_tmrt * d_tmrt * pa * pa + & 1519 ( 1.68737969e07 ) * ta * d_tmrt * d_tmrt * d_tmrt * pa * pa + & 1520 ( 2.67489271e08 ) * va * d_tmrt * d_tmrt * d_tmrt * pa * pa + & 1521 ( 3.23926897e09 ) * d_tmrt * d_tmrt * d_tmrt * d_tmrt * pa * pa + & 1522 ( 3.53874123e02 ) * pa * pa * pa + & 1523 ( 2.21201190e01 ) * ta * pa * pa * pa + & 1524 ( 1.55126038e02 ) * ta * ta * pa * pa * pa + & 1525 ( 2.63917279e04 ) * ta * ta * ta * pa * pa * pa + & 1526 ( 4.53433455e02 ) * va * pa * pa * pa + & 1527 ( 4.32943862e03 ) * ta * va * pa * pa * pa + & 1528 ( 1.45389826e04 ) * ta * ta * va * pa * pa * pa + & 1529 ( 2.17508610e04 ) * va * va * pa * pa * pa + & 1530 ( 6.66724702e05 ) * ta * va * va * pa * pa * pa + & 1531 ( 3.33217140e05 ) * va * va * va * pa * pa * pa + & 1532 ( 2.26921615e03 ) * d_tmrt * pa * pa * pa + & 1533 ( 3.80261982e04 ) * ta * d_tmrt * pa * pa * pa + & 1534 ( 5.45314314e09 ) * ta * ta * d_tmrt * pa * pa * pa + & 1535 ( 7.96355448e04 ) * va * d_tmrt * pa * pa * pa + & 1536 ( 2.53458034e05 ) * ta * va * d_tmrt * pa * pa * pa + & 1537 ( 6.31223658e06 ) * va * va * d_tmrt * pa * pa * pa + & 1538 ( 3.02122035e04 ) * d_tmrt * d_tmrt * pa * pa * pa + & 1539 ( 4.77403547e06 ) * ta * d_tmrt * d_tmrt * pa * pa * pa + & 1540 ( 1.73825715e06 ) * va * d_tmrt * d_tmrt * pa * pa * pa + & 1541 ( 4.09087898e07 ) * d_tmrt * d_tmrt * d_tmrt * pa * pa * pa + & 1542 ( 6.14155345e01 ) * pa * pa * pa * pa + & 1543 ( 6.16755931e02 ) * ta * pa * pa * pa * pa + & 1544 ( 1.33374846e03 ) * ta * ta * pa * pa * pa * pa + & 1545 ( 3.55375387e03 ) * va * pa * pa * pa * pa + & 1546 ( 5.13027851e04 ) * ta * va * pa * pa * pa * pa + & 1547 ( 1.02449757e04 ) * va * va * pa * pa * pa * pa + & 1548 ( 1.48526421e03 ) * d_tmrt * pa * pa * pa * pa + & 1549 ( 4.11469183e05 ) * ta * d_tmrt * pa * pa * pa * pa + & 1550 ( 6.80434415e06 ) * va * d_tmrt * pa * pa * pa * pa + & 1551 ( 9.77675906e06 ) * d_tmrt * d_tmrt * pa * pa * pa * pa + & 1552 ( 8.82773108e02 ) * pa * pa * pa * pa * pa + & 1553 ( 3.01859306e03 ) * ta * pa * pa * pa * pa * pa + & 1554 ( 1.04452989e03 ) * va * pa * pa * pa * pa * pa + & 1555 ( 2.47090539e04 ) * d_tmrt * pa * pa * pa * pa * pa + & 1556 ( 1.48348065e03 ) * 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 "KlimaMichel" 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 PMVguess 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 KlimaMichel 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 !> SIVERSION: ACTLEV W m2, 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 psychophysical 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 ! Newtonapproximation 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 psychophysical 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 * activitypa ) 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.96e8_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_clothingTT) 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.125E03_wp, +1.106E03_wp, 4.570E04_wp /) 2361 coeff(2,1:5) = & 2362 (/ 0.795_wp, 0.713_wp, 8.880E03_wp, 1.803E03_wp, 2.816E03_wp/) 2363 coeff(3,1:5) = & 2364 (/ +0.05761_wp, +0.458_wp, 1.829E02_wp, 5.577E03_wp, 1.970E03_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 GehanGeorge, 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* (2nodemodel, 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 psychophysical 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 ! newtonapproximation 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 psychophysical 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 * activitypa ) 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.96e8_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 !< Energyloss 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 !< DuBois area (mÂ²) 3196 REAL(wp), INTENT( OUT ) :: aeff !< Effective surface area (mÂ²) 3197 REAL(wp), INTENT( OUT ) :: esw !< Energyloss 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 !< DuBois area (mÂ²) 3455 REAL(wp), INTENT( IN ) :: esw !< energyloss 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 1439 3551 1440 3552 END MODULE biometeorology_mod
Note: See TracChangeset
for help on using the changeset viewer.