Changeset 3742
- Timestamp:
- Feb 14, 2019 11:25:22 AM (6 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 2 edited
-
biometeorology_mod.f90 (modified) (130 diffs)
-
time_integration.f90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/biometeorology_mod.f90
r3740 r3742 27 27 ! ----------------- 28 28 ! $Id$ 29 ! - Allocation of the input _av grids was moved to the "sum" section of 30 ! bio_3d_data_averaging to make sure averaging is only done once! 31 ! - Moved call of bio_calculate_thermal_index_maps from biometeorology module to 32 ! time_integration to make sure averaged input is updated before calculating. 33 ! 34 ! 3740 2019-02-13 12:35:12Z dom_dwd_user 29 35 ! - Added safety-meassure to catch the case that 'bio_mrt_av' is stated after 30 36 ! 'bio_<index>' in the output section of the p3d file. … … 125 131 126 132 USE arrays_3d, & 127 ONLY: pt, p, u, v, w, q133 ONLY: pt, p, u, v, w, q 128 134 129 135 USE averaging, & 130 ONLY: pt_av, q_av, u_av, v_av, w_av136 ONLY: pt_av, q_av, u_av, v_av, w_av 131 137 132 138 USE basic_constants_and_equations_mod, & 133 ONLY: c_p, degc_to_k, l_v, magnus, sigma_sb, pi139 ONLY: c_p, degc_to_k, l_v, magnus, sigma_sb, pi 134 140 135 141 USE control_parameters, & 136 ONLY: average_count_3d, biometeorology, dz, dz_stretch_factor,&137 dz_stretch_level, humidity, initializing_actions, nz_do3d,&138 surface_pressure142 ONLY: average_count_3d, biometeorology, dz, dz_stretch_factor, & 143 dz_stretch_level, humidity, initializing_actions, nz_do3d, & 144 surface_pressure 139 145 140 146 USE date_and_time_mod, & … … 142 148 143 149 USE grid_variables, & 144 ONLY: ddx, dx, ddy, dy150 ONLY: ddx, dx, ddy, dy 145 151 146 152 USE indices, & 147 ONLY: nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg,&148 nysg, nyng153 ONLY: nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg, & 154 nysg, nyng 149 155 150 156 USE kinds !< Set precision of INTEGER and REAL arrays according to PALM … … 156 162 !-- Import radiation model to obtain input for mean radiant temperature 157 163 USE radiation_model_mod, & 158 ONLY: ix, iy, iz, id, mrt_nlevels, mrt_include_sw,&159 mrtinsw, mrtinlw, mrtbl, nmrtbl, radiation,&160 radiation_interactions, rad_sw_in,&161 rad_sw_out, rad_lw_in, rad_lw_out164 ONLY: ix, iy, iz, id, mrt_nlevels, mrt_include_sw, & 165 mrtinsw, mrtinlw, mrtbl, nmrtbl, radiation, & 166 radiation_interactions, rad_sw_in, & 167 rad_sw_out, rad_lw_in, rad_lw_out 162 168 163 169 USE surface_mod, & 164 ONLY:get_topography_top_index_ji170 ONLY: get_topography_top_index_ji 165 171 166 172 IMPLICIT NONE … … 190 196 ! 191 197 !-- 198 LOGICAL :: thermal_comfort = .FALSE. !< Enables or disables the entire thermal comfort part 192 199 LOGICAL :: do_average_theta = .FALSE. !< switch: do theta averaging in this module? (if .FALSE. this is done globally) 193 200 LOGICAL :: do_average_q = .FALSE. !< switch: do e averaging in this module? … … 196 203 LOGICAL :: do_average_w = .FALSE. !< switch: do w averaging in this module? 197 204 LOGICAL :: do_average_mrt = .FALSE. !< switch: do mrt averaging in this module? 198 LOGICAL :: average_trigger_perct = .FALSE. !< update averaged input on call to bio_perct? 199 LOGICAL :: average_trigger_utci = .FALSE. !< update averaged input on call to bio_utci? 200 LOGICAL :: average_trigger_pet = .FALSE. !< update averaged input on call to bio_pet? 201 202 LOGICAL :: thermal_comfort = .FALSE. !< Enables or disables thermal comfort part 203 LOGICAL :: do_calculate_perct = .FALSE. !< Turn index PT (instant. input) on or off 204 LOGICAL :: do_calculate_perct_av = .FALSE. !< Turn index PT (averaged input) on or off 205 LOGICAL :: do_calculate_pet = .FALSE. !< Turn index PET (instant. input) on or off 206 LOGICAL :: do_calculate_pet_av = .FALSE. !< Turn index PET (averaged input) on or off 207 LOGICAL :: do_calculate_utci = .FALSE. !< Turn index UTCI (instant. input) on or off 208 LOGICAL :: do_calculate_utci_av = .FALSE. !< Turn index UTCI (averaged input) on or off 205 LOGICAL :: average_trigger_perct = .FALSE. !< update averaged input on call to bio_perct? 206 LOGICAL :: average_trigger_utci = .FALSE. !< update averaged input on call to bio_utci? 207 LOGICAL :: average_trigger_pet = .FALSE. !< update averaged input on call to bio_pet? 208 LOGICAL :: do_calculate_perct = .FALSE. !< Turn index PT (instant. input) on or off 209 LOGICAL :: do_calculate_perct_av = .FALSE. !< Turn index PT (averaged input) on or off 210 LOGICAL :: do_calculate_pet = .FALSE. !< Turn index PET (instant. input) on or off 211 LOGICAL :: do_calculate_pet_av = .FALSE. !< Turn index PET (averaged input) on or off 212 LOGICAL :: do_calculate_utci = .FALSE. !< Turn index UTCI (instant. input) on or off 213 LOGICAL :: do_calculate_utci_av = .FALSE. !< Turn index UTCI (averaged input) on or off 209 214 210 215 ! … … 399 404 400 405 CASE ( 'bio_mrt' ) 401 ! 402 !-- Consider the case 'bio_mrt' is called after some thermal index. 403 !-- In that case do_average_mrt will be .TRUE. leading to a double 404 !-- averaging. The other averaged input quantities are core 405 !-- components, that will hopefully always be treated before the 406 !-- bio methods are called. 407 IF ( .NOT. do_average_mrt .AND. & 408 .NOT. ALLOCATED( mrt_av_grid ) ) THEN 406 407 IF ( .NOT. ALLOCATED( mrt_av_grid ) ) THEN 409 408 ALLOCATE( mrt_av_grid(nmrtbl) ) 410 409 ENDIF 411 410 mrt_av_grid = 0.0_wp 411 do_average_mrt = .FALSE. !< overwrite if that was enabled somehow 412 412 413 413 … … 425 425 !-- Only proceed here if this was not done for any index before. This 426 426 !-- is done only once during the whole model run. 427 IF ( .NOT. average_trigger_perct .AND.&428 .NOT. average_trigger_utci .AND.&429 .NOT. average_trigger_pet ) THEN427 IF ( .NOT. average_trigger_perct .AND. & 428 .NOT. average_trigger_utci .AND. & 429 .NOT. average_trigger_pet ) THEN 430 430 ! 431 431 !-- Memorize the first index called to control averaging 432 IF ( TRIM( variable ) == 'bio_perct*' ) THEN432 IF ( TRIM( variable ) == 'bio_perct*' ) THEN 433 433 average_trigger_perct = .TRUE. 434 434 ENDIF 435 IF ( TRIM( variable ) == 'bio_utci*' ) THEN435 IF ( TRIM( variable ) == 'bio_utci*' ) THEN 436 436 average_trigger_utci = .TRUE. 437 437 ENDIF 438 IF ( TRIM( variable ) == 'bio_pet*' ) THEN438 IF ( TRIM( variable ) == 'bio_pet*' ) THEN 439 439 average_trigger_pet = .TRUE. 440 440 ENDIF 441 441 ENDIF 442 442 ! 443 !-- Only continue if var is the index, that controls averaging. 444 !-- Break immediatelly (doing nothing) for the other indices. 445 IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') RETURN 446 IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*') RETURN 447 IF ( average_trigger_pet .AND. TRIM( variable ) /= 'bio_pet*') RETURN 448 ! 449 !-- Now memorize which of the input grids are not averaged by other 450 !-- modules. Set averaging switch to .TRUE. in that case. 451 IF ( .NOT. ALLOCATED( pt_av ) ) THEN 452 ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 453 do_average_theta = .TRUE. 454 pt_av = 0.0_wp 455 ENDIF 456 457 IF ( .NOT. ALLOCATED( q_av ) ) THEN 458 ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 459 do_average_q = .TRUE. 460 q_av = 0.0_wp 461 ENDIF 462 463 IF ( .NOT. ALLOCATED( u_av ) ) THEN 464 ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 465 do_average_u = .TRUE. 466 u_av = 0.0_wp 467 ENDIF 468 469 IF ( .NOT. ALLOCATED( v_av ) ) THEN 470 ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 471 do_average_v = .TRUE. 472 v_av = 0.0_wp 473 ENDIF 474 475 IF ( .NOT. ALLOCATED( w_av ) ) THEN 476 ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 477 do_average_w = .TRUE. 478 w_av = 0.0_wp 479 ENDIF 480 481 IF ( .NOT. ALLOCATED( mrt_av_grid ) ) THEN 482 ALLOCATE( mrt_av_grid(nmrtbl) ) 483 do_average_mrt = .TRUE. 484 mrt_av_grid = 0.0_wp 485 ENDIF 443 !-- Allocation of the input _av grids was moved to the "sum" section to 444 !-- make sure averaging is only done once! 445 486 446 487 447 CASE ( 'uvem_vitd3dose*' ) … … 505 465 !-- that case do_average_mrt will be .TRUE. leading to a double- 506 466 !-- averaging. 507 IF ( .NOT. do_average_mrt .AND.ALLOCATED( mrt_av_grid ) ) THEN467 IF ( .NOT. do_average_mrt .AND. ALLOCATED( mrt_av_grid ) ) THEN 508 468 509 469 IF ( mrt_include_sw ) THEN … … 511 471 ( ( human_absorb * mrtinsw(:) + & 512 472 mrtinlw(:) ) / & 513 ( human_emiss * sigma_sb ) ) **.25_wp - degc_to_k473 ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k 514 474 ELSE 515 475 mrt_av_grid(:) = mrt_av_grid(:) + & 516 476 ( mrtinlw(:) / & 517 ( human_emiss * sigma_sb ) ) **.25_wp - degc_to_k477 ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k 518 478 ENDIF 519 479 ENDIF … … 523 483 !-- Only continue if the current index is the one to trigger the input 524 484 !-- averaging, see above 525 IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') & 526 RETURN 527 IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*') & 528 RETURN 529 IF ( average_trigger_pet .AND. TRIM( variable ) /= 'bio_pet*') & 530 RETURN 531 532 IF ( ALLOCATED( pt_av ) .AND. do_average_theta ) THEN 485 IF ( average_trigger_perct .AND. TRIM( variable ) /= & 486 'bio_perct*') RETURN 487 IF ( average_trigger_utci .AND. TRIM( variable ) /= & 488 'bio_utci*') RETURN 489 IF ( average_trigger_pet .AND. TRIM( variable ) /= & 490 'bio_pet*') RETURN 491 ! 492 !-- Now memorize which of the input grids are not averaged by other 493 !-- modules. Set averaging switch to .TRUE. and allocate the respective 494 !-- grid in that case. 495 IF ( .NOT. ALLOCATED( pt_av ) ) THEN !< if not averaged by other module 496 ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 497 do_average_theta = .TRUE. !< memorize, that bio is responsible 498 pt_av = 0.0_wp 499 ENDIF 500 IF ( ALLOCATED( pt_av ) .AND. do_average_theta ) THEN 533 501 DO i = nxl, nxr 534 502 DO j = nys, nyn … … 540 508 ENDIF 541 509 542 IF ( ALLOCATED( q_av ) .AND. do_average_q ) THEN 510 IF ( .NOT. ALLOCATED( q_av ) ) THEN 511 ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 512 do_average_q = .TRUE. 513 q_av = 0.0_wp 514 ENDIF 515 IF ( ALLOCATED( q_av ) .AND. do_average_q ) THEN 543 516 DO i = nxl, nxr 544 517 DO j = nys, nyn … … 550 523 ENDIF 551 524 552 IF ( ALLOCATED( u_av ) .AND. do_average_u ) THEN 525 IF ( .NOT. ALLOCATED( u_av ) ) THEN 526 ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 527 do_average_u = .TRUE. 528 u_av = 0.0_wp 529 ENDIF 530 IF ( ALLOCATED( u_av ) .AND. do_average_u ) THEN 553 531 DO i = nxlg, nxrg !< yes, ghost points are required here! 554 532 DO j = nysg, nyng … … 560 538 ENDIF 561 539 562 IF ( ALLOCATED( v_av ) .AND. do_average_v ) THEN 540 IF ( .NOT. ALLOCATED( v_av ) ) THEN 541 ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 542 do_average_v = .TRUE. 543 v_av = 0.0_wp 544 ENDIF 545 IF ( ALLOCATED( v_av ) .AND. do_average_v ) THEN 563 546 DO i = nxlg, nxrg !< yes, ghost points are required here! 564 547 DO j = nysg, nyng … … 570 553 ENDIF 571 554 572 IF ( ALLOCATED( w_av ) .AND. do_average_w ) THEN 555 IF ( .NOT. ALLOCATED( w_av ) ) THEN 556 ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 557 do_average_w = .TRUE. 558 w_av = 0.0_wp 559 ENDIF 560 IF ( ALLOCATED( w_av ) .AND. do_average_w ) THEN 573 561 DO i = nxlg, nxrg !< yes, ghost points are required here! 574 562 DO j = nysg, nyng … … 580 568 ENDIF 581 569 582 IF ( ALLOCATED( mrt_av_grid ) .AND. do_average_mrt ) THEN 570 IF ( .NOT. ALLOCATED( mrt_av_grid ) ) THEN 571 ALLOCATE( mrt_av_grid(nmrtbl) ) 572 do_average_mrt = .TRUE. 573 mrt_av_grid = 0.0_wp 574 ENDIF 575 IF ( ALLOCATED( mrt_av_grid ) .AND. do_average_mrt ) THEN 583 576 584 577 IF ( mrt_include_sw ) THEN 585 578 mrt_av_grid(:) = mrt_av_grid(:) + & 586 579 ( ( human_absorb * mrtinsw(:) + & 587 mrtinlw(:) ) / & ! human_emiss * mrtinlw(:) /588 ( human_emiss * sigma_sb ) ) **.25_wp - degc_to_k580 mrtinlw(:) ) / & 581 ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k 589 582 ELSE 590 583 mrt_av_grid(:) = mrt_av_grid(:) + & 591 ( mrtinlw(:) / & ! ( human_emiss * mrtinlw(:) /592 ( human_emiss * sigma_sb ) ) **.25_wp - degc_to_k584 ( mrtinlw(:) / & 585 ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k 593 586 ENDIF 594 587 ENDIF … … 596 589 !-- This is a cumulated dose. No mode == 'average' for this quantity. 597 590 CASE ( 'uvem_vitd3dose*' ) 598 IF ( ALLOCATED( vitd3_exposure_av ) ) THEN591 IF ( ALLOCATED( vitd3_exposure_av ) ) THEN 599 592 DO i = nxlg, nxrg 600 593 DO j = nysg, nyng … … 618 611 !-- that case do_average_mrt will be .TRUE. leading to a double- 619 612 !-- averaging. 620 IF ( .NOT. do_average_mrt .AND.ALLOCATED( mrt_av_grid ) ) THEN613 IF ( .NOT. do_average_mrt .AND. ALLOCATED( mrt_av_grid ) ) THEN 621 614 mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp ) 622 615 ENDIF … … 625 618 ! 626 619 !-- Only continue if update index, see above 627 IF ( average_trigger_perct .AND.&628 TRIM( variable ) /= 'bio_perct*' ) RETURN629 IF ( average_trigger_utci .AND.&630 TRIM( variable ) /= 'bio_utci*' ) RETURN631 IF ( average_trigger_pet .AND.&632 TRIM( variable ) /= 'bio_pet*' ) RETURN633 634 IF ( ALLOCATED( pt_av ) .AND. do_average_theta )THEN620 IF ( average_trigger_perct .AND. & 621 TRIM( variable ) /= 'bio_perct*' ) RETURN 622 IF ( average_trigger_utci .AND. & 623 TRIM( variable ) /= 'bio_utci*' ) RETURN 624 IF ( average_trigger_pet .AND. & 625 TRIM( variable ) /= 'bio_pet*' ) RETURN 626 627 IF ( ALLOCATED( pt_av ) .AND. do_average_theta ) THEN 635 628 DO i = nxl, nxr 636 629 DO j = nys, nyn … … 643 636 ENDIF 644 637 645 IF ( ALLOCATED( q_av ) .AND. do_average_q )THEN638 IF ( ALLOCATED( q_av ) .AND. do_average_q ) THEN 646 639 DO i = nxl, nxr 647 640 DO j = nys, nyn … … 654 647 ENDIF 655 648 656 IF ( ALLOCATED( u_av ) .AND. do_average_u )THEN649 IF ( ALLOCATED( u_av ) .AND. do_average_u ) THEN 657 650 DO i = nxlg, nxrg !< yes, ghost points are required here! 658 651 DO j = nysg, nyng … … 665 658 ENDIF 666 659 667 IF ( ALLOCATED( v_av ) .AND. do_average_v )THEN660 IF ( ALLOCATED( v_av ) .AND. do_average_v ) THEN 668 661 DO i = nxlg, nxrg 669 662 DO j = nysg, nyng … … 676 669 ENDIF 677 670 678 IF ( ALLOCATED( w_av ) .AND. do_average_w )THEN671 IF ( ALLOCATED( w_av ) .AND. do_average_w ) THEN 679 672 DO i = nxlg, nxrg 680 673 DO j = nysg, nyng … … 687 680 ENDIF 688 681 689 IF ( ALLOCATED( mrt_av_grid ) .AND. do_average_mrt ) THEN 690 mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp ) 691 ENDIF 692 693 ! 694 !-- Udate all thermal index grids with updated averaged input 695 CALL bio_calculate_thermal_index_maps ( .TRUE. ) 682 IF ( ALLOCATED( mrt_av_grid ) .AND. do_average_mrt ) THEN 683 mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, & 684 KIND=wp ) 685 ENDIF 696 686 697 687 ! … … 730 720 SELECT CASE ( TRIM( var ) ) 731 721 ! 732 !-- Allocate a temporary array with the desired output dimensions. 733 !-- Arrays for time-averaged thermal indices are also allocated here because they are not running 734 !-- through the standard averaging procedure in bio_3d_data_averaging as the values of the 735 !-- averaged thermal indices are derived in a single step based on priorly averaged arrays 736 !-- (see bio_calculate_thermal_index_maps). 722 !-- Allocate a temporary array with the desired output dimensions. 723 !-- Arrays for time-averaged thermal indices are also allocated here because 724 !-- they are not running through the standard averaging procedure in 725 !-- bio_3d_data_averaging as the values of the averaged thermal indices are 726 !-- derived in a single step based on priorly averaged arrays (see 727 !-- bio_calculate_thermal_index_maps). 737 728 CASE ( 'bio_mrt' ) 738 729 unit = 'degree_C' … … 746 737 unit = 'degree_C' 747 738 thermal_comfort = .TRUE. 748 IF ( j == 0 ) THEN !< if instantaneous input739 IF ( j == 0 ) THEN !< if instantaneous input 749 740 do_calculate_perct = .TRUE. 750 741 IF ( .NOT. ALLOCATED( perct ) ) THEN … … 763 754 unit = 'degree_C' 764 755 thermal_comfort = .TRUE. 765 IF ( j == 0 ) THEN756 IF ( j == 0 ) THEN 766 757 do_calculate_utci = .TRUE. 767 758 IF ( .NOT. ALLOCATED( utci ) ) THEN … … 780 771 unit = 'degree_C' 781 772 thermal_comfort = .TRUE. 782 IF ( j == 0 ) THEN773 IF ( j == 0 ) THEN 783 774 do_calculate_pet = .TRUE. 784 775 IF ( .NOT. ALLOCATED( pet ) ) THEN … … 838 829 ! 839 830 !-- Further checks if thermal comfort output is desired. 840 IF ( thermal_comfort .AND.unit == 'degree_C' ) THEN831 IF ( thermal_comfort .AND. unit == 'degree_C' ) THEN 841 832 842 833 ! 843 834 !-- Break if required modules "radiation" and "humidity" are not avalable. 844 IF ( .NOT. radiation ) THEN835 IF ( .NOT. radiation ) THEN 845 836 message_string = 'output of "' // TRIM( var ) // '" require' & 846 837 // 's radiation = .TRUE.' … … 855 846 unit = 'illegal' 856 847 ENDIF 857 IF ( mrt_nlevels == 0 ) THEN848 IF ( mrt_nlevels == 0 ) THEN 858 849 message_string = 'output of "' // TRIM( var ) // '" require' & 859 850 // 's mrt_nlevels > 0' … … 934 925 j = mrtbl(iy,l) 935 926 k = mrtbl(iz,l) 936 IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. j > nyn .OR.&937 i < nxl .OR. i > nxr )CYCLE927 IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. & 928 j > nyn .OR. i < nxl .OR. i > nxr ) CYCLE 938 929 IF ( av == 0 ) THEN 939 930 IF ( mrt_include_sw ) THEN 940 931 local_pf(i,j,k) = ( ( human_absorb * mrtinsw(l) + & 941 932 mrtinlw(l) ) / & 942 ( human_emiss * sigma_sb ) ) ** .25_wp -&933 ( human_emiss * sigma_sb ) )**.25_wp - & 943 934 degc_to_k 944 935 ELSE 945 936 local_pf(i,j,k) = ( mrtinlw(l) / & 946 ( human_emiss * sigma_sb ) ) ** .25_wp -&937 ( human_emiss * sigma_sb ) )**.25_wp - & 947 938 degc_to_k 948 939 ENDIF … … 968 959 ENDDO 969 960 ENDDO 970 END 961 ENDIF 971 962 972 963 … … 986 977 ENDDO 987 978 ENDDO 988 END 979 ENDIF 989 980 990 981 … … 1004 995 ENDDO 1005 996 ENDDO 1006 END 997 ENDIF 1007 998 1008 999 ! … … 1087 1078 j = mrtbl(iy,l) 1088 1079 k = mrtbl(iz,l) 1089 IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. j > nyn .OR.&1090 i < nxl .OR. i > nxr )CYCLE1080 IF ( k < nzb_do .OR. k > nzt_do .OR. j < nys .OR. & 1081 j > nyn .OR. i < nxl .OR. i > nxr ) CYCLE 1091 1082 IF ( av == 0 ) THEN 1092 1083 IF ( mrt_include_sw ) THEN 1093 1084 local_pf(i,j,k) = REAL( ( ( human_absorb * mrtinsw(l) + & 1094 1085 mrtinlw(l) ) / & 1095 ( human_emiss * sigma_sb ) ) ** .25_wp -&1086 ( human_emiss * sigma_sb ) )**.25_wp - & 1096 1087 degc_to_k, KIND = sp ) 1097 1088 ELSE 1098 1089 local_pf(i,j,k) = REAL( ( mrtinlw(l) / & 1099 ( human_emiss * sigma_sb ) ) ** .25_wp -&1090 ( human_emiss * sigma_sb ) )**.25_wp - & 1100 1091 degc_to_k, KIND = sp ) 1101 1092 ENDIF … … 1147 1138 is2d = ( var(l-1:l) == 'xy' ) 1148 1139 1149 IF ( var(1:4) == 'bio_' ) THEN1140 IF ( var(1:4) == 'bio_' ) THEN 1150 1141 found = .TRUE. 1151 1142 grid_x = 'x' 1152 1143 grid_y = 'y' 1153 1144 grid_z = 'zu' 1154 IF ( is2d .AND. var(1:7) /= 'bio_mrt' ) grid_z = 'zu1'1155 ENDIF 1156 1157 IF ( is2d .AND. var(1:4) == 'uvem' ) THEN1145 IF ( is2d .AND. var(1:7) /= 'bio_mrt' ) grid_z = 'zu1' 1146 ENDIF 1147 1148 IF ( is2d .AND. var(1:4) == 'uvem' ) THEN 1158 1149 grid_x = 'x' 1159 1150 grid_y = 'y' … … 1241 1232 ONLY: message_string 1242 1233 1243 IF ( .NOT. radiation_interactions ) THEN1234 IF ( .NOT. radiation_interactions ) THEN 1244 1235 message_string = 'The mrt calculation requires ' // & 1245 1236 'enabled radiation_interactions but it ' // & … … 1485 1476 ! 1486 1477 !-- We need to differentiate if averaged input is desired (av == .TRUE.) or not. 1487 IF ( av ) THEN1478 IF ( av ) THEN 1488 1479 ! 1489 1480 !-- Make sure tmrt_av_grid is present and initialize with the fill value … … 1505 1496 k = mrtbl(iz,l) 1506 1497 IF ( k - get_topography_top_index_ji( j, i, 's' ) == & 1507 bio_cell_level + 1_iwp) THEN1498 bio_cell_level + 1_iwp) THEN 1508 1499 ! 1509 1500 !-- Averaging was done before, so we can just copy the result here … … 1531 1522 k = mrtbl(iz,l) 1532 1523 IF ( k - get_topography_top_index_ji( j, i, 's' ) == & 1533 bio_cell_level + 1_iwp) THEN1524 bio_cell_level + 1_iwp) THEN 1534 1525 IF ( mrt_include_sw ) THEN 1535 1526 tmrt_grid(j,i) = ( ( human_absorb * mrtinsw(l) + & 1536 1527 mrtinlw(l) ) / & 1537 ( human_emiss * sigma_sb ) ) ** .25_wp -&1528 ( human_emiss * sigma_sb ) )**.25_wp - & 1538 1529 degc_to_k 1539 1530 ELSE 1540 1531 tmrt_grid(j,i) = ( mrtinlw(l) / & 1541 ( human_emiss * sigma_sb ) ) ** .25_wp -&1532 ( human_emiss * sigma_sb ) )**.25_wp - & 1542 1533 degc_to_k 1543 1534 ENDIF … … 1585 1576 !-- Avoid non-representative horizontal u and v of 0.0 m/s too close to ground 1586 1577 k_wind = k 1587 IF ( bio_cell_level < 1_iwp ) THEN1578 IF ( bio_cell_level < 1_iwp ) THEN 1588 1579 k_wind = k + 1_iwp 1589 1580 ENDIF 1590 1581 ! 1591 1582 !-- Determine local values: 1592 IF ( average_input ) THEN1583 IF ( average_input ) THEN 1593 1584 ! 1594 1585 !-- Calculate ta from Tp assuming dry adiabatic laps rate 1595 ta = pt_av(k, j,i) - ( 0.0098_wp * dz(1) * ( k + .5_wp ) ) - degc_to_k1586 ta = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * ( k + .5_wp ) ) - degc_to_k 1596 1587 1597 1588 vp = bio_fill_value 1598 IF ( humidity .AND. ALLOCATED( q_av ) )THEN1599 vp = q_av(k, j,i)1589 IF ( humidity .AND. ALLOCATED( q_av ) ) THEN 1590 vp = q_av(k,j,i) 1600 1591 ENDIF 1601 1592 … … 1606 1597 ! 1607 1598 !-- Calculate ta from Tp assuming dry adiabatic laps rate 1608 ta = pt(k, j,i) - ( 0.0098_wp * dz(1) * ( k + .5_wp ) ) - degc_to_k1599 ta = pt(k,j,i) - ( 0.0098_wp * dz(1) * ( k + .5_wp ) ) - degc_to_k 1609 1600 1610 1601 vp = bio_fill_value 1611 IF ( humidity ) THEN1612 vp = q(k, j,i)1602 IF ( humidity ) THEN 1603 vp = q(k,j,i) 1613 1604 ENDIF 1614 1605 … … 1625 1616 !-- The magnus formula is limited to temperatures up to 333.15 K to 1626 1617 ! avoid negative values of vp_sat 1627 IF ( vp > -998._wp ) THEN1618 IF ( vp > -998._wp ) THEN 1628 1619 vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) ) 1629 1620 vp = vp * pair / ( vp + 0.622_wp ) 1630 IF ( vp > vp_sat ) vp = vp_sat1621 IF ( vp > vp_sat ) vp = vp_sat 1631 1622 ENDIF 1632 1623 ! 1633 1624 !-- local mtr value at [i,j] 1634 1625 tmrt = bio_fill_value !< this can be a valid result (e.g. for inside some ostacle) 1635 IF ( .NOT. average_input ) THEN1626 IF ( .NOT. average_input ) THEN 1636 1627 ! 1637 1628 !-- Use MRT from RTM precalculated in tmrt_grid … … 1674 1665 CALL bio_calculate_mrt_grid ( av ) 1675 1666 1676 DO i = nxl, nxr1677 DO j = nys, nyn1667 DO i = nxl, nxr 1668 DO j = nys, nyn 1678 1669 ! 1679 1670 !-- Determine local input conditions … … 1689 1680 perct_ij = bio_fill_value !< some obstacle 1690 1681 utci_ij = bio_fill_value 1691 IF ( .NOT. ( tmrt_ij <= -998._wp .OR. vp <= -998._wp ) )THEN1682 IF ( .NOT. ( tmrt_ij <= -998._wp .OR. vp <= -998._wp ) ) THEN 1692 1683 ! 1693 1684 !-- Calculate static thermal indices based on local tmrt 1694 1685 clo = bio_fill_value 1695 1686 1696 IF ( do_calculate_perct .OR. do_calculate_perct_av )THEN1687 IF ( do_calculate_perct .OR. do_calculate_perct_av ) THEN 1697 1688 ! 1698 1689 !-- Estimate local perceived temperature … … 1701 1692 ENDIF 1702 1693 1703 IF ( do_calculate_utci .OR. do_calculate_utci_av )THEN1694 IF ( do_calculate_utci .OR. do_calculate_utci_av ) THEN 1704 1695 ! 1705 1696 !-- Estimate local universal thermal climate index … … 1708 1699 ENDIF 1709 1700 1710 IF ( do_calculate_pet .OR. do_calculate_pet_av )THEN1701 IF ( do_calculate_pet .OR. do_calculate_pet_av ) THEN 1711 1702 ! 1712 1703 !-- Estimate local physiologically equivalent temperature 1713 1704 CALL calculate_pet_static( ta, vp, ws, tmrt_ij, pair, pet_ij ) 1714 1705 ENDIF 1715 END 1716 1717 1718 IF ( av ) THEN1706 ENDIF 1707 1708 1709 IF ( av ) THEN 1719 1710 ! 1720 1711 !-- Write results for selected averaged indices 1721 1712 IF ( do_calculate_perct_av ) THEN 1722 1713 perct_av(j, i) = perct_ij 1723 END 1724 IF ( do_calculate_utci_av ) THEN1714 ENDIF 1715 IF ( do_calculate_utci_av ) THEN 1725 1716 utci_av(j, i) = utci_ij 1726 END 1727 IF ( do_calculate_pet_av ) THEN1717 ENDIF 1718 IF ( do_calculate_pet_av ) THEN 1728 1719 pet_av(j, i) = pet_ij 1729 END 1720 ENDIF 1730 1721 ELSE 1731 1722 ! … … 1733 1724 IF ( do_calculate_perct ) THEN 1734 1725 perct(j, i) = perct_ij 1735 END 1736 IF ( do_calculate_utci ) THEN1726 ENDIF 1727 IF ( do_calculate_utci ) THEN 1737 1728 utci(j, i) = utci_ij 1738 END 1739 IF ( do_calculate_pet ) THEN1729 ENDIF 1730 IF ( do_calculate_pet ) THEN 1740 1731 pet(j, i) = pet_ij 1741 END 1742 END 1743 1744 END 1745 END 1732 ENDIF 1733 ENDIF 1734 1735 ENDDO 1736 ENDDO 1746 1737 1747 1738 END SUBROUTINE bio_calculate_thermal_index_maps … … 1782 1773 ! 1783 1774 !-- return immediatelly if nothing to do! 1784 IF ( .NOT. thermal_comfort ) THEN1775 IF ( .NOT. thermal_comfort ) THEN 1785 1776 RETURN 1786 1777 ENDIF 1787 1778 ! 1788 1779 !-- If clo equals the initial value, this is the initial call 1789 IF ( clo <= -998._wp ) THEN1780 IF ( clo <= -998._wp ) THEN 1790 1781 ! 1791 1782 !-- Initialize instationary perceived temperature with personalized … … 1881 1872 ! 1882 1873 !-- Check if input values in range after Broede et al. (2012) 1883 IF ( ( d_tmrt > 70._wp ) .OR. ( d_tmrt < -30._wp ) .OR.&1884 ( vp >= 50._wp ) ) THEN1874 IF ( ( d_tmrt > 70._wp ) .OR. ( d_tmrt < -30._wp ) .OR. & 1875 ( vp >= 50._wp ) ) THEN 1885 1876 utci_ij = bio_fill_value 1886 1877 RETURN … … 1888 1879 ! 1889 1880 !-- Apply eq. 2 in Broede et al. (2012) for ta out of bounds 1890 IF ( ta > 50._wp ) THEN1881 IF ( ta > 50._wp ) THEN 1891 1882 offset = ta - 50._wp 1892 1883 ta = 50._wp 1893 1884 ENDIF 1894 IF ( ta < -50._wp ) THEN1885 IF ( ta < -50._wp ) THEN 1895 1886 offset = ta + 50._wp 1896 1887 ta = -50._wp … … 1900 1891 !-- humidity values below 0.5 m/s or 5%, respectively, the 1901 1892 !-- user is advised to use the lower bounds for the calculations. 1902 IF ( va < 0.5_wp ) va = 0.5_wp1903 IF ( va > 17._wp ) va = 17._wp1893 IF ( va < 0.5_wp ) va = 0.5_wp 1894 IF ( va > 17._wp ) va = 17._wp 1904 1895 1905 1896 ! … … 2229 2220 ! 2230 2221 !-- decision: firstly calculate for winter or summer clothing 2231 IF ( ta <= 10._wp ) THEN2222 IF ( ta <= 10._wp ) THEN 2232 2223 ! 2233 2224 !-- First guess: winter clothing insulation: cold stress … … 2236 2227 pmv_w = pmva 2237 2228 2238 IF ( pmva > 0._wp ) THEN2229 IF ( pmva > 0._wp ) THEN 2239 2230 ! 2240 2231 !-- Case summer clothing insulation: heat load ? … … 2243 2234 top ) 2244 2235 pmv_s = pmva 2245 IF ( pmva <= 0._wp ) THEN2236 IF ( pmva <= 0._wp ) THEN 2246 2237 ! 2247 2238 !-- Case: comfort achievable by varying clothing insulation … … 2249 2240 CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo, & 2250 2241 pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo ) 2251 IF ( ncount < 0_iwp ) THEN2242 IF ( ncount < 0_iwp ) THEN 2252 2243 nerr = -1_iwp 2253 2244 RETURN 2254 2245 ENDIF 2255 ELSE IF ( pmva > 0.06_wp ) THEN2246 ELSE IF ( pmva > 0.06_wp ) THEN 2256 2247 clo = 0.5_wp 2257 2248 CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, & 2258 2249 pmva, top ) 2259 2250 ENDIF 2260 ELSE IF ( pmva < -0.11_wp ) THEN2251 ELSE IF ( pmva < -0.11_wp ) THEN 2261 2252 clo = 1.75_wp 2262 2253 CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, & … … 2270 2261 pmv_s = pmva 2271 2262 2272 IF ( pmva < 0._wp ) THEN2263 IF ( pmva < 0._wp ) THEN 2273 2264 ! 2274 2265 !-- Case winter clothing insulation: cold stress ? … … 2278 2269 pmv_w = pmva 2279 2270 2280 IF ( pmva >= 0._wp ) THEN2271 IF ( pmva >= 0._wp ) THEN 2281 2272 ! 2282 2273 !-- Case: comfort achievable by varying clothing insulation … … 2284 2275 CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo, & 2285 2276 pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo ) 2286 IF ( ncount < 0_iwp ) THEN2277 IF ( ncount < 0_iwp ) THEN 2287 2278 nerr = -1_iwp 2288 2279 RETURN 2289 2280 ENDIF 2290 ELSE IF ( pmva < -0.11_wp ) THEN2281 ELSE IF ( pmva < -0.11_wp ) THEN 2291 2282 clo = 1.75_wp 2292 2283 CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, & 2293 2284 pmva, top ) 2294 2285 ENDIF 2295 ELSE IF ( pmva > 0.06_wp ) THEN2286 ELSE IF ( pmva > 0.06_wp ) THEN 2296 2287 clo = 0.5_wp 2297 2288 CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva, & … … 2305 2296 CALL perct_regression ( pmva, clo, perct_ij ) 2306 2297 ptc = perct_ij 2307 IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp )THEN2298 IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN 2308 2299 ! 2309 2300 !-- Adjust for cold conditions according to Gagge 1986 2310 2301 CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv ) 2311 IF ( nerr_cold > 0_iwp ) nerr = -5_iwp2302 IF ( nerr_cold > 0_iwp ) nerr = -5_iwp 2312 2303 pmvs = pmva - d_pmv 2313 IF ( pmvs > -0.11_wp ) THEN2304 IF ( pmvs > -0.11_wp ) THEN 2314 2305 d_pmv = 0._wp 2315 2306 pmvs = -0.11_wp … … 2319 2310 ! clo_fanger = clo 2320 2311 clon = clo 2321 IF ( clo > 0.5_wp .AND. perct_ij <= 8.73_wp )THEN2312 IF ( clo > 0.5_wp .AND. perct_ij <= 8.73_wp ) THEN 2322 2313 ! 2323 2314 !-- Required clothing insulation (ireq) is exclusively defined for … … 2330 2321 sultrieness = .FALSE. 2331 2322 d_std = -99._wp 2332 IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp )THEN2323 IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN 2333 2324 ! 2334 2325 !-- Adjust for warm/humid conditions according to Gagge 1986 … … 2337 2328 pmvs = pmva + d_pmv 2338 2329 CALL perct_regression ( pmvs, clo, perct_ij ) 2339 IF ( sult_lim < 99._wp ) THEN2340 IF ( (perct_ij - ptc) > sult_lim ) sultrieness = .TRUE.2330 IF ( sult_lim < 99._wp ) THEN 2331 IF ( (perct_ij - ptc) > sult_lim ) sultrieness = .TRUE. 2341 2332 ! 2342 2333 !-- Set factor to threshold for sultriness 2343 IF ( dgtcstd /= 0_iwp ) THEN2334 IF ( dgtcstd /= 0_iwp ) THEN 2344 2335 d_std = ( ( perct_ij - ptc ) - dgtcm ) / dgtcstd 2345 2336 ENDIF … … 2369 2360 2370 2361 2371 IF ( ta < 0._wp ) THEN2362 IF ( ta < 0._wp ) THEN 2372 2363 ! 2373 2364 !-- ta < 0 (degC): water vapour pressure over ice … … 2451 2442 clo_upper = wclo 2452 2443 y_upper = pmv_w 2453 IF ( ( y_lower > 0._wp .AND. y_upper < 0._wp ) .OR.&2454 ( y_lower < 0._wp .AND. y_upper > 0._wp ) )THEN2444 IF ( ( y_lower > 0._wp .AND. y_upper < 0._wp ) .OR. & 2445 ( y_lower < 0._wp .AND. y_upper > 0._wp ) ) THEN 2455 2446 x_lower = clo_lower 2456 2447 x_upper = clo_upper 2457 2448 x_ridder = guess_0 2458 2449 2459 DO j = 1_iwp, max_iteration2450 DO j = 1_iwp, max_iteration 2460 2451 x_average = 0.5_wp * ( x_lower + x_upper ) 2461 2452 CALL fanger ( ta, tmrt, vp, ws, pair, x_average, actlev, eta, & 2462 2453 y_average, top ) 2463 2454 sroot = SQRT ( y_average**2 - y_lower * y_upper ) 2464 IF ( sroot == 0._wp ) THEN2455 IF ( sroot == 0._wp ) THEN 2465 2456 clo_res = x_average 2466 2457 nerr = j … … 2469 2460 x_new = x_average + ( x_average - x_lower ) * & 2470 2461 ( SIGN ( 1._wp, y_lower - y_upper ) * y_average / sroot ) 2471 IF ( ABS ( x_new - x_ridder ) <= eps ) THEN2462 IF ( ABS ( x_new - x_ridder ) <= eps ) THEN 2472 2463 clo_res = x_ridder 2473 2464 nerr = j … … 2477 2468 CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta, & 2478 2469 y_new, top ) 2479 IF ( y_new == 0._wp ) THEN2470 IF ( y_new == 0._wp ) THEN 2480 2471 clo_res = x_ridder 2481 2472 nerr = j 2482 2473 RETURN 2483 2474 ENDIF 2484 IF ( SIGN ( y_average, y_new ) /= y_average ) THEN2475 IF ( SIGN ( y_average, y_new ) /= y_average ) THEN 2485 2476 x_lower = x_average 2486 2477 y_lower = y_average 2487 2478 x_upper = x_ridder 2488 2479 y_upper = y_new 2489 ELSE IF ( SIGN ( y_lower, y_new ) /= y_lower ) THEN2480 ELSE IF ( SIGN ( y_lower, y_new ) /= y_lower ) THEN 2490 2481 x_upper = x_ridder 2491 2482 y_upper = y_new 2492 ELSE IF ( SIGN ( y_upper, y_new ) /= y_upper ) THEN2483 ELSE IF ( SIGN ( y_upper, y_new ) /= y_upper ) THEN 2493 2484 x_lower = x_ridder 2494 2485 y_lower = y_new … … 2500 2491 RETURN 2501 2492 ENDIF 2502 IF ( ABS ( x_upper - x_lower ) <= eps ) THEN2493 IF ( ABS ( x_upper - x_lower ) <= eps ) THEN 2503 2494 clo_res = x_ridder 2504 2495 nerr = j … … 2511 2502 clo_res = y_new 2512 2503 RETURN 2513 ELSE IF ( y_lower == 0. ) THEN2504 ELSE IF ( y_lower == 0. ) THEN 2514 2505 x_ridder = clo_lower 2515 ELSE IF ( y_upper == 0. ) THEN2506 ELSE IF ( y_upper == 0. ) THEN 2516 2507 x_ridder = clo_upper 2517 2508 ELSE … … 2541 2532 REAL(wp), INTENT ( OUT ) :: perct_ij !< perct (degC) corresponding to given PMV / clo 2542 2533 2543 IF ( pmv <= -0.11_wp ) THEN2534 IF ( pmv <= -0.11_wp ) THEN 2544 2535 perct_ij = 5.805_wp + 12.6784_wp * pmv 2545 2536 ELSE 2546 IF ( pmv >= + 0.01_wp ) THEN2537 IF ( pmv >= + 0.01_wp ) THEN 2547 2538 perct_ij = 16.826_wp + 6.163_wp * pmv 2548 2539 ELSE … … 2608 2599 !-- Clo must be > 0. to avoid div. by 0! 2609 2600 clo = in_clo 2610 IF ( clo <= 0._wp ) clo = .001_wp2601 IF ( clo <= 0._wp ) clo = .001_wp 2611 2602 ! 2612 2603 !-- f_cl = Increase in surface due to clothing … … 2615 2606 !-- Case of free convection (ws < 0.1 m/s ) not considered 2616 2607 ws = in_ws 2617 IF ( ws < .1_wp ) THEN2608 IF ( ws < .1_wp ) THEN 2618 2609 ws = .1_wp 2619 2610 ENDIF … … 2638 2629 !-- Newton-approximation with air temperature as initial guess 2639 2630 t_clothing = ta 2640 DO i = 1, 32631 DO i = 1, 3 2641 2632 t_clothing = t_clothing - ( ( t_clothing + degc_to_k )**4 + t_clothing & 2642 2633 * dc - gc ) / ( 4._wp * ( t_clothing + degc_to_k )**3 + dc ) … … 2660 2651 z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + degc_to_k )**4 - ( tmrt + & 2661 2652 degc_to_k )**4 ) 2662 IF ( ABS ( t_clothing - tmrt ) > 0._wp ) THEN2653 IF ( ABS ( t_clothing - tmrt ) > 0._wp ) THEN 2663 2654 hr = z5 / f_cl / ( t_clothing - tmrt ) 2664 2655 ELSE … … 2760 2751 ! 2761 2752 !-- Test for compliance with regression range 2762 IF ( pmva < -1.0_wp .OR. pmva > 7.0_wp )THEN2753 IF ( pmva < -1.0_wp .OR. pmva > 7.0_wp ) THEN 2763 2754 nerr = -2_iwp 2764 2755 ELSE … … 2772 2763 p10 = 0.05_wp * svp_ta 2773 2764 p95 = 1.00_wp * svp_ta 2774 IF ( vp >= p10 .AND. vp <= p95 )THEN2765 IF ( vp >= p10 .AND. vp <= p95 ) THEN 2775 2766 pa = vp 2776 2767 ELSE 2777 2768 nerr = -3_iwp 2778 IF ( vp < p10 ) THEN2769 IF ( vp < p10 ) THEN 2779 2770 ! 2780 2771 !-- Due to conditions of regression: r.H. >= 5 % … … 2786 2777 ENDIF 2787 2778 ENDIF 2788 IF ( pa > 0._wp ) THEN2779 IF ( pa > 0._wp ) THEN 2789 2780 ! 2790 2781 !-- Natural logarithm of pa … … 2796 2787 !-- Ratio actual water vapour pressure to that of a r.H. of 50 % 2797 2788 pa_p50 = 0.5_wp * svp_ta 2798 IF ( pa_p50 > 0._wp .AND. pa > 0._wp )THEN2789 IF ( pa_p50 > 0._wp .AND. pa > 0._wp ) THEN 2799 2790 dapa = apa - LOG ( pa_p50 ) 2800 2791 pa_p50 = pa / pa_p50 … … 2805 2796 ! 2806 2797 !-- Square root of wind velocity 2807 IF ( ws >= 0._wp ) THEN2798 IF ( ws >= 0._wp ) THEN 2808 2799 sqvel = SQRT ( ws ) 2809 2800 ELSE … … 2816 2807 !-- Select the valid regression coefficients 2817 2808 nreg = INT ( pmv ) 2818 IF ( nreg < 0_iwp ) THEN2809 IF ( nreg < 0_iwp ) THEN 2819 2810 ! 2820 2811 !-- value of the FUNCTION in the case pmv <= -1 … … 2823 2814 ENDIF 2824 2815 weight = MOD ( pmv, 1._wp ) 2825 IF ( weight < 0._wp ) weight = 0._wp2826 IF ( nreg > 5_iwp ) THEN2816 IF ( weight < 0._wp ) weight = 0._wp 2817 IF ( nreg > 5_iwp ) THEN 2827 2818 ! nreg=6 2828 2819 nreg = 5_iwp 2829 2820 weight = pmv - 5._wp 2830 2821 weight2 = pmv - 6._wp 2831 IF ( weight2 > 0_iwp ) THEN2822 IF ( weight2 > 0_iwp ) THEN 2832 2823 weight = ( weight - weight2 ) / weight 2833 2824 ENDIF … … 2836 2827 !-- Regression valid for 0. <= pmv <= 6. 2837 2828 dpmv_1 = & 2838 + bpa ( nreg ) * pa&2839 + bpmv ( nreg ) * pmv&2840 + bapa ( nreg ) * apa&2841 + bta ( nreg ) * ta&2842 + bdtmrt ( nreg ) * dtmrt&2843 + bdapa ( nreg ) * dapa&2844 + bsqvel ( nreg ) * sqvel&2845 + bpa_p50 ( nreg ) * pa_p50&2846 + aconst ( nreg)2829 + bpa(nreg) * pa & 2830 + bpmv(nreg) * pmv & 2831 + bapa(nreg) * apa & 2832 + bta(nreg) * ta & 2833 + bdtmrt(nreg) * dtmrt & 2834 + bdapa(nreg) * dapa & 2835 + bsqvel(nreg) * sqvel & 2836 + bpa_p50(nreg) * pa_p50 & 2837 + aconst(nreg) 2847 2838 2848 2839 dpmv_2 = 0._wp 2849 IF ( nreg < 6_iwp ) THEN2840 IF ( nreg < 6_iwp ) THEN 2850 2841 dpmv_2 = & 2851 + bpa ( nreg + 1_iwp ) * pa&2852 + bpmv ( nreg + 1_iwp ) * pmv&2853 + bapa ( nreg + 1_iwp ) * apa&2854 + bta ( nreg + 1_iwp ) * ta&2855 + bdtmrt ( nreg + 1_iwp ) * dtmrt&2856 + bdapa ( nreg + 1_iwp ) * dapa&2857 + bsqvel ( nreg + 1_iwp ) * sqvel&2858 + bpa_p50 ( nreg + 1_iwp ) * pa_p50&2859 + aconst ( nreg + 1_iwp)2842 + bpa(nreg+1_iwp) * pa & 2843 + bpmv(nreg+1_iwp) * pmv & 2844 + bapa(nreg+1_iwp) * apa & 2845 + bta(nreg+1_iwp) * ta & 2846 + bdtmrt(nreg+1_iwp) * dtmrt & 2847 + bdapa(nreg+1_iwp) * dapa & 2848 + bsqvel(nreg+1_iwp) * sqvel & 2849 + bpa_p50(nreg+1_iwp) * pa_p50 & 2850 + aconst(nreg+1_iwp) 2860 2851 ENDIF 2861 2852 ! … … 2863 2854 deltapmv = ( 1._wp - weight ) * dpmv_1 + weight * dpmv_2 2864 2855 pmvs = pmva + deltapmv 2865 IF ( ( pmvs ) < 0._wp ) THEN2856 IF ( ( pmvs ) < 0._wp ) THEN 2866 2857 ! 2867 2858 !-- Prevent negative pmv* due to problems with clothing insulation 2868 2859 nerr = -4_iwp 2869 IF ( pmvs > -0.11_wp ) THEN2860 IF ( pmvs > -0.11_wp ) THEN 2870 2861 ! 2871 2862 !-- Threshold from FUNCTION perct_regression for winter clothing insulation … … 2925 2916 dperctstd = 999999._wp 2926 2917 2927 IF ( perct_ij < 16.826_wp .OR. perct_ij > 56._wp )THEN2918 IF ( perct_ij < 16.826_wp .OR. perct_ij > 56._wp ) THEN 2928 2919 ! 2929 2920 !-- Unallowed value of classical perct! … … 2941 2932 !-- Value of the FUNCTION 2942 2933 sultr_res = dperctm + faktor * dperctstd 2943 IF ( ABS ( sultr_res ) > 99._wp ) sultr_res = +99._wp2934 IF ( ABS ( sultr_res ) > 99._wp ) sultr_res = +99._wp 2944 2935 2945 2936 END SUBROUTINE calc_sultr … … 2995 2986 dtmrt = tmrt - ta 2996 2987 sqrt_ws = ws 2997 IF ( sqrt_ws < 0.10_wp ) THEN2988 IF ( sqrt_ws < 0.10_wp ) THEN 2998 2989 sqrt_ws = 0.10_wp 2999 2990 ELSE … … 3001 2992 ENDIF 3002 2993 3003 DO i = 1, 33004 delta_cold 3005 IF ( i < 3 ) THEN3006 pmv_cross 2994 DO i = 1, 3 2995 delta_cold(i) = 0._wp 2996 IF ( i < 3 ) THEN 2997 pmv_cross(i) = pmvc 3007 2998 ENDIF 3008 2999 ENDDO 3009 3000 3010 DO i = 1, 33001 DO i = 1, 3 3011 3002 ! 3012 3003 !-- Regression constant for given meteorological variables 3013 reg_a(i) = coeff(i, 1) + coeff(i, 3) * ta + coeff(i, 4) *&3004 reg_a(i) = coeff(i,1) + coeff(i,3) * ta + coeff(i,4) * & 3014 3005 sqrt_ws + coeff(i,5)*dtmrt 3015 delta_cold(i) = reg_a(i) + coeff(i, 3006 delta_cold(i) = reg_a(i) + coeff(i,2) * pmvc 3016 3007 ENDDO 3017 3008 ! 3018 3009 !-- Intersection points of regression lines in terms of Fanger's PMV 3019 DO i = 1, 23020 r_nenner = coeff (i, 2) - coeff (i + 1,2)3021 IF ( ABS ( r_nenner ) > 0.00001_wp ) THEN3022 pmv_cross(i) = ( reg_a (i + 1) - reg_a(i) ) / r_nenner3010 DO i = 1, 2 3011 r_nenner = coeff(i,2) - coeff(i+1,2) 3012 IF ( ABS ( r_nenner ) > 0.00001_wp ) THEN 3013 pmv_cross(i) = ( reg_a(i+1) - reg_a(i) ) / r_nenner 3023 3014 ELSE 3024 3015 nerr = 1_iwp … … 3028 3019 3029 3020 i_bin = 3 3030 DO i = 1, 23031 IF ( pmva > pmv_cross (i) )THEN3021 DO i = 1, 2 3022 IF ( pmva > pmv_cross(i) ) THEN 3032 3023 i_bin = i 3033 3024 EXIT … … 3062 3053 ! 3063 3054 !-- range_1 range_2 range_3 3064 DATA (coef(i, 3065 DATA (coef(i, 3066 DATA (coef(i, 3067 DATA (coef(i, 3068 3069 IF ( pmva <= -2.1226_wp ) THEN3055 DATA (coef(i,0), i = 1, n_bin) /0.0941540_wp, -0.1506620_wp, -0.0871439_wp/ 3056 DATA (coef(i,1), i = 1, n_bin) /0.0783162_wp, -1.0612651_wp, 0.1695040_wp/ 3057 DATA (coef(i,2), i = 1, n_bin) /0.1350144_wp, -1.0049144_wp, -0.0167627_wp/ 3058 DATA (coef(i,3), i = 1, n_bin) /0.1104037_wp, -0.2005277_wp, -0.0003230_wp/ 3059 3060 IF ( pmva <= -2.1226_wp ) THEN 3070 3061 i_bin = 3_iwp 3071 ELSE IF ( pmva <= -1.28_wp ) THEN3062 ELSE IF ( pmva <= -1.28_wp ) THEN 3072 3063 i_bin = 2_iwp 3073 3064 ELSE … … 3075 3066 ENDIF 3076 3067 3077 dpmv_adj = coef( i_bin, n_ca)3068 dpmv_adj = coef(i_bin,n_ca) 3078 3069 pmv = 1._wp 3079 3070 3080 DO i = n_ca + 1, n_ce3071 DO i = n_ca + 1, n_ce 3081 3072 pmv = pmv * pmva 3082 dpmv_adj = dpmv_adj + coef(i_bin, 3073 dpmv_adj = dpmv_adj + coef(i_bin,i) * pmv 3083 3074 ENDDO 3084 3075 RETURN … … 3116 3107 ! 3117 3108 !-- Regression only defined for perct <= 10 (degC) 3118 IF ( ireq_neutral < 0.5_wp ) THEN3119 IF ( ireq_neutral < 0.48_wp ) THEN3109 IF ( ireq_neutral < 0.5_wp ) THEN 3110 IF ( ireq_neutral < 0.48_wp ) THEN 3120 3111 nerr = 1_iwp 3121 3112 ENDIF … … 3125 3116 !-- Minimal required clothing insulation: maximal acceptable body cooling 3126 3117 ireq_minimal = 1.26_wp - 0.0588_wp * top02 3127 IF ( nerr > 0_iwp ) THEN3118 IF ( nerr > 0_iwp ) THEN 3128 3119 ireq_minimal = ireq_neutral 3129 3120 ENDIF … … 3152 3143 ! 3153 3144 !-- According to Gehan-George, for children 3154 IF ( age < 19_iwp ) THEN3155 IF ( age < 5_iwp ) THEN3145 IF ( age < 19_iwp ) THEN 3146 IF ( age < 5_iwp ) THEN 3156 3147 surf = 0.02667_wp * height**0.42246_wp * weight**0.51456_wp 3157 3148 RETURN 3158 END 3149 ENDIF 3159 3150 surf = 0.03050_wp * height**0.35129_wp * weight**0.54375_wp 3160 3151 RETURN 3161 END 3152 ENDIF 3162 3153 ! 3163 3154 !-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if … … 3200 3191 3201 3192 CALL surface_area ( height, weight, INT( age ), a_surf ) 3202 s = height * 100._wp / ( weight **( 1._wp / 3._wp ) )3193 s = height * 100._wp / ( weight**( 1._wp / 3._wp ) ) 3203 3194 factor = 1._wp + .004_wp * ( 30._wp - age ) 3204 3195 basic_heat_prod = 0. 3205 IF ( sex == 1_iwp ) THEN3206 basic_heat_prod = 3.45_wp * weight ** ( 3._wp / 4._wp ) * ( factor +&3196 IF ( sex == 1_iwp ) THEN 3197 basic_heat_prod = 3.45_wp * weight**( 3._wp / 4._wp ) * ( factor + & 3207 3198 .01_wp * ( s - 43.4_wp ) ) 3208 ELSE IF ( sex == 2_iwp ) THEN3209 basic_heat_prod = 3.19_wp * weight ** ( 3._wp / 4._wp ) * ( factor +&3199 ELSE IF ( sex == 2_iwp ) THEN 3200 basic_heat_prod = 3.19_wp * weight**( 3._wp / 4._wp ) * ( factor + & 3210 3201 .018_wp * ( s - 42.1_wp ) ) 3211 END 3202 ENDIF 3212 3203 3213 3204 energy_prod = work + basic_heat_prod … … 3295 3286 ! 3296 3287 !-- Decision: firstly calculate for winter or summer clothing 3297 IF ( ta <= 10._wp ) THEN3288 IF ( ta <= 10._wp ) THEN 3298 3289 ! 3299 3290 !-- First guess: winter clothing insulation: cold stress … … 3304 3295 pmv_w = pmva 3305 3296 3306 IF ( pmva > 0._wp ) THEN3297 IF ( pmva > 0._wp ) THEN 3307 3298 ! 3308 3299 !-- Case summer clothing insulation: heat load ? … … 3312 3303 t_clothing, storage, dt, pmva ) 3313 3304 pmv_s = pmva 3314 IF ( pmva <= 0._wp ) THEN3305 IF ( pmva <= 0._wp ) THEN 3315 3306 ! 3316 3307 !-- Case: comfort achievable by varying clothing insulation … … 3318 3309 CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta , sclo, & 3319 3310 pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo ) 3320 IF ( ncount < 0_iwp ) THEN3311 IF ( ncount < 0_iwp ) THEN 3321 3312 nerr = -1_iwp 3322 3313 RETURN 3323 3314 ENDIF 3324 ELSE IF ( pmva > 0.06_wp ) THEN3315 ELSE IF ( pmva > 0.06_wp ) THEN 3325 3316 clo = 0.5_wp 3326 3317 t_clothing = bio_fill_value … … 3328 3319 t_clothing, storage, dt, pmva ) 3329 3320 ENDIF 3330 ELSE IF ( pmva < -0.11_wp ) THEN3321 ELSE IF ( pmva < -0.11_wp ) THEN 3331 3322 clo = 1.75_wp 3332 3323 t_clothing = bio_fill_value … … 3344 3335 pmv_s = pmva 3345 3336 3346 IF ( pmva < 0._wp ) THEN3337 IF ( pmva < 0._wp ) THEN 3347 3338 ! 3348 3339 !-- Case winter clothing insulation: cold stress ? … … 3353 3344 pmv_w = pmva 3354 3345 3355 IF ( pmva >= 0._wp ) THEN3346 IF ( pmva >= 0._wp ) THEN 3356 3347 ! 3357 3348 !-- Case: comfort achievable by varying clothing insulation … … 3359 3350 CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo, & 3360 3351 pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo ) 3361 IF ( ncount < 0_wp ) THEN3352 IF ( ncount < 0_wp ) THEN 3362 3353 nerr = -1_iwp 3363 3354 RETURN 3364 3355 ENDIF 3365 ELSE IF ( pmva < -0.11_wp ) THEN3356 ELSE IF ( pmva < -0.11_wp ) THEN 3366 3357 clo = 1.75_wp 3367 3358 t_clothing = bio_fill_value … … 3369 3360 t_clothing, storage, dt, pmva ) 3370 3361 ENDIF 3371 ELSE IF ( pmva > 0.06_wp ) THEN3362 ELSE IF ( pmva > 0.06_wp ) THEN 3372 3363 clo = 0.5_wp 3373 3364 t_clothing = bio_fill_value … … 3382 3373 CALL perct_regression ( pmva, clo, ipt ) 3383 3374 ptc = ipt 3384 IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp )THEN3375 IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN 3385 3376 ! 3386 3377 !-- Adjust for cold conditions according to Gagge 1986 3387 3378 CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv ) 3388 IF ( nerr_cold > 0_iwp ) nerr = -5_iwp3379 IF ( nerr_cold > 0_iwp ) nerr = -5_iwp 3389 3380 pmvs = pmva - d_pmv 3390 IF ( pmvs > -0.11_wp ) THEN3381 IF ( pmvs > -0.11_wp ) THEN 3391 3382 d_pmv = 0._wp 3392 3383 pmvs = -0.11_wp … … 3396 3387 ! clo_fanger = clo 3397 3388 clon = clo 3398 IF ( clo > 0.5_wp .AND. ipt <= 8.73_wp )THEN3389 IF ( clo > 0.5_wp .AND. ipt <= 8.73_wp ) THEN 3399 3390 ! 3400 3391 !-- Required clothing insulation (ireq) is exclusively defined for … … 3407 3398 sultrieness = .FALSE. 3408 3399 d_std = -99._wp 3409 IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp )THEN3400 IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN 3410 3401 ! 3411 3402 !-- Adjust for warm/humid conditions according to Gagge 1986 … … 3414 3405 pmvs = pmva + d_pmv 3415 3406 CALL perct_regression ( pmvs, clo, ipt ) 3416 IF ( sult_lim < 99._wp ) THEN3417 IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE.3407 IF ( sult_lim < 99._wp ) THEN 3408 IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE. 3418 3409 ENDIF 3419 3410 ENDIF … … 3487 3478 CALL perct_regression ( pmva, clo, ipt ) 3488 3479 3489 IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp )THEN3480 IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN 3490 3481 ! 3491 3482 !-- Adjust for cold conditions according to Gagge 1986 3492 3483 CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv ) 3493 IF ( nerr_cold > 0_iwp ) nerr = -5_iwp3484 IF ( nerr_cold > 0_iwp ) nerr = -5_iwp 3494 3485 pmvs = pmva - d_pmv 3495 IF ( pmvs > -0.11_wp ) THEN3486 IF ( pmvs > -0.11_wp ) THEN 3496 3487 d_pmv = 0._wp 3497 3488 pmvs = -0.11_wp … … 3505 3496 sultrieness = .FALSE. 3506 3497 d_std = -99._wp 3507 IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp )THEN3498 IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN 3508 3499 ! 3509 3500 !-- Adjust for warm/humid conditions according to Gagge 1986 … … 3512 3503 pmvs = pmva + d_pmv 3513 3504 CALL perct_regression ( pmvs, clo, ipt ) 3514 IF ( sult_lim < 99._wp ) THEN3515 IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE.3505 IF ( sult_lim < 99._wp ) THEN 3506 IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE. 3516 3507 ENDIF 3517 3508 ENDIF … … 3583 3574 !-- Clo must be > 0. to avoid div. by 0! 3584 3575 clo = in_clo 3585 IF ( clo < 001._wp ) clo = .001_wp3576 IF ( clo < 001._wp ) clo = .001_wp 3586 3577 ! 3587 3578 !-- Increase in surface due to clothing … … 3590 3581 !-- Case of free convection (ws < 0.1 m/s ) not considered 3591 3582 ws = in_ws 3592 IF ( ws < .1_wp ) THEN3583 IF ( ws < .1_wp ) THEN 3593 3584 ws = .1_wp 3594 3585 ENDIF … … 3610 3601 !-- newton-approximation with air temperature as initial guess 3611 3602 niter = INT( dt * 10._wp, KIND=iwp ) 3612 IF ( niter < 1 ) niter = 1_iwp3603 IF ( niter < 1 ) niter = 1_iwp 3613 3604 adjustrate = 1._wp - EXP ( -1._wp * ( 10._wp / time_equil ) * dt ) 3614 IF ( adjustrate >= 1._wp ) adjustrate = 1._wp3605 IF ( adjustrate >= 1._wp ) adjustrate = 1._wp 3615 3606 adjustrate_cloth = adjustrate * 30._wp 3616 3607 t_clothing = t_cloth 3617 3608 3618 IF ( t_cloth <= -998.0_wp ) THEN ! If initial run3609 IF ( t_cloth <= -998.0_wp ) THEN ! If initial run 3619 3610 niter = 3_iwp 3620 3611 adjustrate = 1._wp … … 3623 3614 ENDIF 3624 3615 3625 DO i = 1, niter3616 DO i = 1, niter 3626 3617 t_clothing = t_clothing - adjustrate_cloth * ( ( t_clothing + & 3627 3618 273.2_wp )**4._wp + t_clothing * & … … 3775 3766 REAL(wp) :: vpex 3776 3767 3777 met = 3.45_wp * mbody ** ( 3._wp / 4._wp ) * (1._wp + 0.004_wp * & 3768 ! 3769 !-- Metabolic heat production 3770 met = 3.45_wp * mbody**( 3._wp / 4._wp ) * (1._wp + 0.004_wp * & 3778 3771 ( 30._wp - age) + 0.010_wp * ( ( ht * 100._wp / & 3779 ( mbody ** ( 1._wp / 3._wp ) ) ) - 43.4_wp ) ) 3780 3772 ( mbody**( 1._wp / 3._wp ) ) ) - 43.4_wp ) ) 3781 3773 met = work + met 3782 3783 3774 int_heat = met * (1._wp - eta) 3784 3775 ! 3785 3776 !-- Sensible respiration energy 3786 3777 tex = 0.47_wp * ta + 21.0_wp 3787 rtv = 1.44_wp * 10._wp **(-6._wp) * met3778 rtv = 1.44_wp * 10._wp**(-6._wp) * met 3788 3779 eres = c_p * (ta - tex) * rtv 3789 3780 ! 3790 3781 !-- Latent respiration energy 3791 vpex = 6.11_wp * 10._wp **( 7.45_wp * tex / ( 235._wp + tex ) )3782 vpex = 6.11_wp * 10._wp**( 7.45_wp * tex / ( 235._wp + tex ) ) 3792 3783 erel = 0.623_wp * l_v / pair * ( vpa - vpex ) * rtv 3793 3784 ! … … 3889 3880 LOGICAL :: skipIncreaseCount !< iteration control flag 3890 3881 3891 wetsk = 0._wp 3892 adu = 0.203_wp * mbody ** 0.425_wp * ht ** 0.725_wp 3893 3894 hc = 2.67_wp + ( 6.5_wp * v ** 0.67_wp) 3895 hc = hc * (pair /po) ** 0.55_wp 3882 ! 3883 !-- Initialize 3884 wetsk = 0._wp !< skin is dry everywhere on init (no non-evaporated sweat) 3885 ! 3886 !-- Set Du Bois Area for the sample person 3887 adu = 0.203_wp * mbody**0.425_wp * ht**0.725_wp 3888 ! 3889 !-- Initialize convective heat considering local air preassure 3890 hc = 2.67_wp + ( 6.5_wp * v**0.67_wp ) 3891 hc = hc * ( pair / po )**0.55_wp 3892 ! 3893 !-- Set surface modification by posture (the person will always stand) 3896 3894 feff = 0.725_wp !< Posture: 0.725 for stading 3897 facl = (- 2.36_wp + 173.51_wp * clo - 100.76_wp * clo * clo + 19.28_wp & 3898 * (clo ** 3._wp)) / 100._wp 3899 3900 IF ( facl > 1._wp ) facl = 1._wp 3895 ! 3896 !-- Set surface modification by clothing 3897 facl = ( - 2.36_wp + 173.51_wp * clo - 100.76_wp * clo * clo + 19.28_wp & 3898 * ( clo**3._wp ) ) / 100._wp 3899 IF ( facl > 1._wp ) facl = 1._wp 3900 ! 3901 !-- Initialize heat resistences 3901 3902 rcl = ( clo / 6.45_wp ) / facl 3902 3903 IF ( clo >= 2._wp ) y = 1._wp 3903 3904 IF ( ( clo > 0.6_wp ) .AND. ( clo < 2._wp ) ) y = ( ht - 0.2_wp ) / ht 3905 IF ( ( clo <= 0.6_wp ) .AND. ( clo > 0.3_wp ) ) y = 0.5_wp 3906 IF ( ( clo <= 0.3_wp ) .AND. ( clo > 0._wp ) ) y = 0.1_wp 3907 3908 r2 = adu * (fcl - 1._wp + facl) / (2._wp * 3.14_wp * ht * y) 3909 r1 = facl * adu / (2._wp * 3.14_wp * ht * y) 3910 3904 IF ( ( clo > 0.6_wp ) .AND. ( clo < 2._wp ) ) y = ( ht - 0.2_wp ) / ht 3905 IF ( ( clo <= 0.6_wp ) .AND. ( clo > 0.3_wp ) ) y = 0.5_wp 3906 IF ( ( clo <= 0.3_wp ) .AND. ( clo > 0._wp ) ) y = 0.1_wp 3907 r2 = adu * ( fcl - 1._wp + facl ) / ( 2._wp * 3.14_wp * ht * y ) 3908 r1 = facl * adu / ( 2._wp * 3.14_wp * ht * y ) 3911 3909 di = r2 - r1 3912 ! 3913 !-- Skin temperatur 3914 DO j = 1, 7 3910 3911 ! 3912 !-- Estimate skin temperatur 3913 DO j = 1, 7 3915 3914 3916 3915 tsk = 34._wp … … 3920 3919 enbal2 = 0._wp 3921 3920 3922 DO i = 1, 100 ! allow for 100 iterations max3921 DO i = 1, 100 ! allow for 100 iterations max 3923 3922 acl = adu * facl + adu * ( fcl - 1._wp ) 3924 3923 rclo2 = emcl * sigma_sb * ( ( tcl + degc_to_k )**4._wp - & 3925 ( tmrt + degc_to_k )** 3924 ( tmrt + degc_to_k )**4._wp ) * feff 3926 3925 htcl = 6.28_wp * ht * y * di / ( rcl * LOG( r2 / r1 ) * acl ) 3927 3926 tsk = 1._wp / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl … … 3930 3929 aeff = adu * feff 3931 3930 rbare = aeff * ( 1._wp - facl ) * emsk * sigma_sb * & 3932 ( ( tmrt + degc_to_k )** 4._wp - ( tsk + degc_to_k )**4._wp )3931 ( ( tmrt + degc_to_k )**4._wp - ( tsk + degc_to_k )**4._wp ) 3933 3932 rclo = feff * acl * emcl * sigma_sb * & 3934 ( ( tmrt + degc_to_k )** 4._wp - ( tcl + degc_to_k )**4._wp )3933 ( ( tmrt + degc_to_k )**4._wp - ( tcl + degc_to_k )**4._wp ) 3935 3934 rsum = rbare + rclo 3936 3935 ! … … 3952 3951 c(9) = 5.28_wp * adu - c(5) - c(4) * tsk 3953 3952 c(10) = c(9) * c(9) - 4._wp * c(4) * & 3954 ( c(5) * tsk - c(0) - 5.28_wp * adu * tsk )3955 3956 IF ( tsk == 36._wp ) tsk = 36.01_wp3953 ( c(5) * tsk - c(0) - 5.28_wp * adu * tsk ) 3954 3955 IF ( tsk == 36._wp ) tsk = 36.01_wp 3957 3956 tcore(7) = c(0) / ( 5.28_wp * adu + c(1) * 6.3_wp / 3600._wp ) + tsk 3958 3957 tcore(3) = c(0) / ( 5.28_wp * adu + ( c(1) * 6.3_wp / 3600._wp ) / & 3959 3958 ( 1._wp + 0.5_wp * ( 34._wp - tsk ) ) ) + tsk 3960 IF ( c(10) >= 0._wp ) THEN3959 IF ( c(10) >= 0._wp ) THEN 3961 3960 tcore(6) = ( - c(9) - c(10)**0.5_wp ) / ( 2._wp * c(4) ) 3962 3961 tcore(1) = ( - c(9) + c(10)**0.5_wp ) / ( 2._wp * c(4) ) 3963 END 3964 3965 IF ( c(8) >= 0._wp ) THEN3966 tcore(2) = ( - c(6) + ABS( c(8) ) **0.5_wp ) / ( 2._wp * c(4) )3967 tcore(5) = ( - c(6) - ABS( c(8) ) **0.5_wp ) / ( 2._wp * c(4) )3962 ENDIF 3963 3964 IF ( c(8) >= 0._wp ) THEN 3965 tcore(2) = ( - c(6) + ABS( c(8) )**0.5_wp ) / ( 2._wp * c(4) ) 3966 tcore(5) = ( - c(6) - ABS( c(8) )**0.5_wp ) / ( 2._wp * c(4) ) 3968 3967 tcore(4) = c(0) / ( 5.28_wp * adu + c(1) * 1._wp / 40._wp ) + tsk 3969 END 3968 ENDIF 3970 3969 ! 3971 3970 !-- Transpiration … … 3974 3973 vpts = 6.11_wp * 10._wp**( 7.45_wp * tsk / ( 235._wp + tsk ) ) 3975 3974 3976 IF ( tbody <= 36.6_wp ) swm = 0._wp !< no need for sweating3975 IF ( tbody <= 36.6_wp ) swm = 0._wp !< no need for sweating 3977 3976 3978 3977 sw = swm … … 3983 3982 wetsk = eswphy / eswpot 3984 3983 3985 IF ( wetsk > 1._wp ) wetsk = 1._wp3984 IF ( wetsk > 1._wp ) wetsk = 1._wp 3986 3985 ! 3987 3986 !-- Sweat production > evaporation? 3988 3987 eswdif = eswphy - eswpot 3989 3988 3990 IF ( eswdif <= 0._wp ) esw = eswpot !< Limit is evaporation3991 IF ( eswdif > 0._wp ) esw = eswphy!< Limit is sweat production3992 IF ( esw > 0._wp ) esw = 0._wp!< Sweat can't be evaporated, no more cooling effect3989 IF ( eswdif <= 0._wp ) esw = eswpot !< Limit is evaporation 3990 IF ( eswdif > 0._wp ) esw = eswphy !< Limit is sweat production 3991 IF ( esw > 0._wp ) esw = 0._wp !< Sweat can't be evaporated, no more cooling effect 3993 3992 ! 3994 3993 !-- Diffusion 3995 rdsk = 0.79_wp * 10._wp **7._wp3994 rdsk = 0.79_wp * 10._wp**7._wp 3996 3995 rdcl = 0._wp 3997 3996 ed = l_v / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( vpa - & … … 4002 4001 vb2 = tcore(j) - 36.6_wp 4003 4002 4004 IF ( vb2 < 0._wp ) vb2 = 0._wp4005 IF ( vb1 < 0._wp ) vb1 = 0._wp4003 IF ( vb2 < 0._wp ) vb2 = 0._wp 4004 IF ( vb1 < 0._wp ) vb1 = 0._wp 4006 4005 vb = ( 6.3_wp + 75._wp * vb2 ) / ( 1._wp + 0.5_wp * vb1 ) 4007 4006 ! … … 4011 4010 !-- Clothing temperature 4012 4011 xx = 0.001_wp 4013 IF ( count1 == 0_iwp ) xx = 1._wp4014 IF ( count1 == 1_iwp ) xx = 0.1_wp4015 IF ( count1 == 2_iwp ) xx = 0.01_wp4016 IF ( count1 == 3_iwp ) xx = 0.001_wp4017 4018 IF ( enbal > 0._wp ) tcl = tcl + xx4019 IF ( enbal < 0._wp ) tcl = tcl - xx4012 IF ( count1 == 0_iwp ) xx = 1._wp 4013 IF ( count1 == 1_iwp ) xx = 0.1_wp 4014 IF ( count1 == 2_iwp ) xx = 0.01_wp 4015 IF ( count1 == 3_iwp ) xx = 0.001_wp 4016 4017 IF ( enbal > 0._wp ) tcl = tcl + xx 4018 IF ( enbal < 0._wp ) tcl = tcl - xx 4020 4019 4021 4020 skipIncreaseCount = .FALSE. 4022 IF ( ( (enbal <= 0._wp ) .AND. (enbal2 > 0._wp ) ) .OR.&4023 ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) )THEN4021 IF ( ( (enbal <= 0._wp ) .AND. (enbal2 > 0._wp ) ) .OR. & 4022 ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) ) THEN 4024 4023 skipIncreaseCount = .TRUE. 4025 4024 ELSE 4026 4025 enbal2 = enbal 4027 4026 count3 = count3 + 1_iwp 4028 END 4029 4030 IF ( ( count3 > 200_iwp ) .OR. skipIncreaseCount )THEN4031 IF ( count1 < 3_iwp ) THEN4027 ENDIF 4028 4029 IF ( ( count3 > 200_iwp ) .OR. skipIncreaseCount ) THEN 4030 IF ( count1 < 3_iwp ) THEN 4032 4031 count1 = count1 + 1_iwp 4033 4032 enbal2 = 0._wp 4034 4033 ELSE 4035 4034 EXIT 4036 END 4037 END 4038 END 4039 4040 IF ( count1 == 3_iwp ) THEN4035 ENDIF 4036 ENDIF 4037 ENDDO 4038 4039 IF ( count1 == 3_iwp ) THEN 4041 4040 SELECT CASE ( j ) 4042 4041 CASE ( 2, 5) 4043 IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND.&4044 ( tsk <= 34.050_wp ) ) ) CYCLE4042 IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND. & 4043 ( tsk <= 34.050_wp ) ) ) CYCLE 4045 4044 CASE ( 6, 1 ) 4046 4045 IF ( c(10) < 0._wp ) CYCLE 4047 IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND.&4048 ( tsk > 33.850_wp ) ) ) CYCLE4046 IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND. & 4047 ( tsk > 33.850_wp ) ) ) CYCLE 4049 4048 CASE ( 3 ) 4050 IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND.&4051 ( tsk <= 34.000_wp ) ) ) CYCLE4049 IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND. & 4050 ( tsk <= 34.000_wp ) ) ) CYCLE 4052 4051 CASE ( 7 ) 4053 IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND.&4054 ( tsk > 34.000_wp ) ) ) CYCLE4052 IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND. & 4053 ( tsk > 34.000_wp ) ) ) CYCLE 4055 4054 CASE default 4056 4055 END SELECT 4057 END 4058 4059 IF ( ( j /= 4_iwp ) .AND. ( vb >= 91._wp ) )CYCLE4060 IF ( ( j == 4_iwp ) .AND. ( vb < 89._wp ) )CYCLE4056 ENDIF 4057 4058 IF ( ( j /= 4_iwp ) .AND. ( vb >= 91._wp ) ) CYCLE 4059 IF ( ( j == 4_iwp ) .AND. ( vb < 89._wp ) ) CYCLE 4061 4060 IF ( vb > 90._wp ) vb = 90._wp 4062 4061 ! 4063 4062 !-- Loses by water 4064 4063 ws = sw * 3600._wp * 1000._wp 4065 IF ( ws > 2000._wp ) ws = 2000._wp4064 IF ( ws > 2000._wp ) ws = 2000._wp 4066 4065 wd = ed / l_v * 3600._wp * ( -1000._wp ) 4067 4066 wr = erel / l_v * 3600._wp * ( -1000._wp ) … … 4070 4069 4071 4070 RETURN 4072 END 4071 ENDDO 4073 4072 END SUBROUTINE heat_exch 4074 4073 … … 4131 4130 enbal2 = 0._wp 4132 4131 4133 DO count1 = 0, 34134 DO i = 1, 125 ! 500 / 44135 hc = 2.67_wp + 6.5_wp * 0.1_wp **0.67_wp4136 hc = hc * ( pair / po ) **0.55_wp4132 DO count1 = 0, 3 4133 DO i = 1, 125 ! 500 / 4 4134 hc = 2.67_wp + 6.5_wp * 0.1_wp**0.67_wp 4135 hc = hc * ( pair / po )**0.55_wp 4137 4136 ! 4138 4137 !-- Radiation 4139 4138 aeff = adu * feff 4140 4139 rbare = aeff * ( 1._wp - facl ) * emsk * sigma_sb * & 4141 ( ( pet_ij + degc_to_k ) ** 4._wp - ( tsk + degc_to_k ) **4._wp )4140 ( ( pet_ij + degc_to_k )**4._wp - ( tsk + degc_to_k )**4._wp ) 4142 4141 rclo = feff * acl * emcl * sigma_sb * & 4143 ( ( pet_ij + degc_to_k ) ** 4._wp - ( tcl + degc_to_k ) **4._wp )4142 ( ( pet_ij + degc_to_k )**4._wp - ( tcl + degc_to_k )**4._wp ) 4144 4143 rsum = rbare + rclo 4145 4144 ! … … 4156 4155 tex = 0.47_wp * pet_ij + 21._wp 4157 4156 eres = c_p * ( pet_ij - tex ) * rtv 4158 vpex = 6.11_wp * 10._wp **( 7.45_wp * tex / ( 235._wp + tex ) )4157 vpex = 6.11_wp * 10._wp**( 7.45_wp * tex / ( 235._wp + tex ) ) 4159 4158 erel = 0.623_wp * l_v / pair * ( 12._wp - vpex ) * rtv 4160 4159 ere = eres + erel … … 4171 4170 IF ( enbal > 0._wp ) pet_ij = pet_ij - xx 4172 4171 IF ( enbal < 0._wp ) pet_ij = pet_ij + xx 4173 IF ( ( enbal <= 0._wp ) .AND. ( enbal2 > 0._wp ) )EXIT4174 IF ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) )EXIT4172 IF ( ( enbal <= 0._wp ) .AND. ( enbal2 > 0._wp ) ) EXIT 4173 IF ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) EXIT 4175 4174 4176 4175 enbal2 = enbal 4177 END 4178 END 4176 ENDDO 4177 ENDDO 4179 4178 END SUBROUTINE pet_iteration 4180 4179 … … 4271 4270 CALL uvem_solar_position 4272 4271 4273 IF (sza .GE. 90) THEN4272 IF (sza .GE. 90) THEN 4274 4273 vitd3_exposure(:,:) = 0.0_wp 4275 4274 ELSE … … 4368 4367 ! 4369 4368 ! !-- extract obstruction from IBSET-Integer_Array ------------------' 4370 IF (consider_obstructions ) THEN4369 IF (consider_obstructions ) THEN 4371 4370 obstruction_temp1 = building_obstruction_f%var_3d(:,j,i) 4372 IF (obstruction_temp1(0) .NE. 9) THEN4371 IF (obstruction_temp1(0) .NE. 9) THEN 4373 4372 DO pobi = 0, 44 4374 4373 DO bi = 0, 7 4375 IF ( btest( obstruction_temp1(pobi), bi ) .EQV. .TRUE.) THEN4374 IF ( btest( obstruction_temp1(pobi), bi ) .EQV. .TRUE.) THEN 4376 4375 obstruction_temp2( ( pobi * 8 ) + bi ) = 1 4377 4376 ELSE … … 4392 4391 4393 4392 obstruction_direct_beam = obstruction( nint(startpos_saa_float), nint( sza / 10.0_wp ) ) 4394 IF (sza .GE. 89.99_wp) THEN4393 IF (sza .GE. 89.99_wp) THEN 4395 4394 sza = 89.99999_wp 4396 4395 ENDIF -
palm/trunk/SOURCE/time_integration.f90
r3739 r3742 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Moved call of bio_calculate_thermal_index_maps from biometeorology module to 28 ! time_integration to make sure averaged input is updated before calculating. 29 ! 30 ! 3739 2019-02-13 08:05:17Z dom_dwd_user 27 31 ! Removed everything related to "time_bio_results" as this is never used. 28 32 ! … … 1635 1639 IF ( time_dopr_listing >= dt_dopr_listing ) THEN 1636 1640 CALL print_1d 1637 time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing, &1641 time_dopr_listing = MOD( time_dopr_listing, MAX( dt_dopr_listing, & 1638 1642 dt_3d ) ) 1639 1643 ENDIF … … 1641 1645 ! 1642 1646 !-- Graphic output for PROFIL 1643 IF ( time_dopr >= dt_dopr &1647 IF ( time_dopr >= dt_dopr & 1644 1648 .AND. time_since_reference_point >= skip_time_dopr ) THEN 1645 1649 IF ( dopr_n /= 0 ) CALL data_output_profiles … … 1658 1662 !-- Output of spectra (formatted for use with PROFIL), in case of no 1659 1663 !-- time averaging, spectra has to be calculated before 1660 IF ( time_dosp >= dt_dosp &1664 IF ( time_dosp >= dt_dosp & 1661 1665 .AND. time_since_reference_point >= skip_time_dosp ) THEN 1662 1666 IF ( average_count_sp == 0 ) CALL calc_spectra … … 1667 1671 ! 1668 1672 !-- 2d-data output (cross-sections) 1669 IF ( time_do2d_xy >= dt_do2d_xy &1673 IF ( time_do2d_xy >= dt_do2d_xy & 1670 1674 .AND. time_since_reference_point >= skip_time_do2d_xy ) THEN 1671 1675 CALL data_output_2d( 'xy', 0 ) 1672 1676 time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) ) 1673 1677 ENDIF 1674 IF ( time_do2d_xz >= dt_do2d_xz &1678 IF ( time_do2d_xz >= dt_do2d_xz & 1675 1679 .AND. time_since_reference_point >= skip_time_do2d_xz ) THEN 1676 1680 CALL data_output_2d( 'xz', 0 ) 1677 1681 time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) ) 1678 1682 ENDIF 1679 IF ( time_do2d_yz >= dt_do2d_yz &1683 IF ( time_do2d_yz >= dt_do2d_yz & 1680 1684 .AND. time_since_reference_point >= skip_time_do2d_yz ) THEN 1681 1685 CALL data_output_2d( 'yz', 0 ) … … 1685 1689 ! 1686 1690 !-- 3d-data output (volume data) 1687 IF ( time_do3d >= dt_do3d &1691 IF ( time_do3d >= dt_do3d & 1688 1692 .AND. time_since_reference_point >= skip_time_do3d ) THEN 1689 1693 CALL data_output_3d( 0 ) … … 1694 1698 !-- Masked data output 1695 1699 DO mid = 1, masks 1696 IF ( time_domask(mid) >= dt_domask(mid) 1700 IF ( time_domask(mid) >= dt_domask(mid) & 1697 1701 .AND. time_since_reference_point >= skip_time_domask(mid) ) THEN 1698 1702 CALL data_output_mask( 0 ) 1699 time_domask(mid) = MOD( time_domask(mid), 1703 time_domask(mid) = MOD( time_domask(mid), & 1700 1704 MAX( dt_domask(mid), dt_3d ) ) 1701 1705 ENDIF … … 1704 1708 ! 1705 1709 !-- Output of time-averaged 2d/3d/masked data 1706 IF ( time_do_av >= dt_data_output_av 1710 IF ( time_do_av >= dt_data_output_av & 1707 1711 .AND. time_since_reference_point >= skip_time_data_output_av ) THEN 1708 1712 CALL average_3d_data 1713 ! 1714 !-- Udate thermal comfort indices based on updated averaged input 1715 IF ( biometeorology .AND. thermal_comfort ) THEN 1716 CALL bio_calculate_thermal_index_maps ( .TRUE. ) 1717 ENDIF 1709 1718 CALL data_output_2d( 'xy', 1 ) 1710 1719 CALL data_output_2d( 'xz', 1 ) … … 1719 1728 !-- Output of surface data, instantaneous and averaged data 1720 1729 IF ( surface_output ) THEN 1721 IF ( time_dosurf >= dt_dosurf 1730 IF ( time_dosurf >= dt_dosurf & 1722 1731 .AND. time_since_reference_point >= skip_time_dosurf ) THEN 1723 1732 CALL surface_data_output( 0 ) 1724 1733 time_dosurf = MOD( time_dosurf, MAX( dt_dosurf, dt_3d ) ) 1725 1734 ENDIF 1726 IF ( time_dosurf_av >= dt_dosurf_av 1735 IF ( time_dosurf_av >= dt_dosurf_av & 1727 1736 .AND. time_since_reference_point >= skip_time_dosurf_av ) THEN 1728 1737 CALL surface_data_output( 1 )
Note: See TracChangeset
for help on using the changeset viewer.
