Changeset 3742


Ignore:
Timestamp:
Feb 14, 2019 11:25:22 AM (5 years ago)
Author:
dom_dwd_user
Message:

biometeorology_mod.f90:
(C) Allocation of the input _av grids was moved to the "sum" section of bio_3d_data_averaging to make sure averaging is only done once!
(B) Moved call of bio_calculate_thermal_index_maps from biometeorology module to time_integration (ll 1712) to make sure averaged input is updated before calculating.

time_integration.f90:
(B) Moved call of bio_calculate_thermal_index_maps from biometeorology module to time_integration (ll 1712) to make sure averaged input is updated before calculating.

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/biometeorology_mod.f90

    r3740 r3742  
    2727! -----------------
    2828! $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
    2935! - Added safety-meassure to catch the case that 'bio_mrt_av' is stated after
    3036! 'bio_<index>' in the output section of the p3d file.
     
    125131
    126132    USE arrays_3d,                                                             &
    127        ONLY:  pt, p, u, v, w, q
     133        ONLY:  pt, p, u, v, w, q
    128134
    129135    USE averaging,                                                             &
    130        ONLY:  pt_av, q_av, u_av, v_av, w_av
     136        ONLY:  pt_av, q_av, u_av, v_av, w_av
    131137
    132138    USE basic_constants_and_equations_mod,                                     &
    133        ONLY:  c_p, degc_to_k, l_v, magnus, sigma_sb, pi
     139        ONLY:  c_p, degc_to_k, l_v, magnus, sigma_sb, pi
    134140
    135141    USE control_parameters,                                                    &
    136        ONLY:  average_count_3d, biometeorology, dz, dz_stretch_factor,         &
    137               dz_stretch_level, humidity, initializing_actions, nz_do3d,       &
    138               surface_pressure
     142        ONLY:  average_count_3d, biometeorology, dz, dz_stretch_factor,        &
     143               dz_stretch_level, humidity, initializing_actions, nz_do3d,      &
     144               surface_pressure
    139145
    140146    USE date_and_time_mod,                                                     &
     
    142148
    143149    USE grid_variables,                                                        &
    144        ONLY:  ddx, dx, ddy, dy
     150        ONLY:  ddx, dx, ddy, dy
    145151
    146152    USE indices,                                                               &
    147        ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg,    &
    148               nysg, nyng
     153        ONLY:  nxl, nxr, nys, nyn, nzb, nzt, nys, nyn, nxl, nxr, nxlg, nxrg,   &
     154               nysg, nyng
    149155
    150156    USE kinds  !< Set precision of INTEGER and REAL arrays according to PALM
     
    156162!-- Import radiation model to obtain input for mean radiant temperature
    157163    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_out
     164        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
    162168
    163169    USE surface_mod,                                                           &
    164        ONLY: get_topography_top_index_ji
     170        ONLY: get_topography_top_index_ji
    165171
    166172    IMPLICIT NONE
     
    190196!
    191197!--
     198    LOGICAL ::  thermal_comfort  = .FALSE.  !< Enables or disables the entire thermal comfort part
    192199    LOGICAL ::  do_average_theta = .FALSE.  !< switch: do theta averaging in this module? (if .FALSE. this is done globally)
    193200    LOGICAL ::  do_average_q     = .FALSE.  !< switch: do e averaging in this module?
     
    196203    LOGICAL ::  do_average_w     = .FALSE.  !< switch: do w averaging in this module?
    197204    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
    209214
    210215!
     
    399404
    400405          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
    409408                   ALLOCATE( mrt_av_grid(nmrtbl) )
    410409                ENDIF
    411410                mrt_av_grid = 0.0_wp
     411                do_average_mrt = .FALSE.  !< overwrite if that was enabled somehow
    412412
    413413
     
    425425!--          Only proceed here if this was not done for any index before. This
    426426!--          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 ) THEN
     427             IF ( .NOT. average_trigger_perct  .AND.                           &
     428                  .NOT. average_trigger_utci   .AND.                           &
     429                  .NOT. average_trigger_pet )  THEN
    430430!
    431431!--             Memorize the first index called to control averaging
    432                 IF ( TRIM( variable ) == 'bio_perct*' ) THEN
     432                IF ( TRIM( variable ) == 'bio_perct*' )  THEN
    433433                    average_trigger_perct = .TRUE.
    434434                ENDIF
    435                 IF ( TRIM( variable ) == 'bio_utci*' ) THEN
     435                IF ( TRIM( variable ) == 'bio_utci*' )  THEN
    436436                    average_trigger_utci = .TRUE.
    437437                ENDIF
    438                 IF ( TRIM( variable ) == 'bio_pet*' ) THEN
     438                IF ( TRIM( variable ) == 'bio_pet*' )  THEN
    439439                    average_trigger_pet = .TRUE.
    440440                ENDIF
    441441             ENDIF
    442442!
    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
    486446
    487447          CASE ( 'uvem_vitd3dose*' )
     
    505465!--          that case do_average_mrt will be .TRUE. leading to a double-
    506466!--          averaging.
    507              IF ( .NOT. do_average_mrt .AND. ALLOCATED( mrt_av_grid ) )  THEN
     467             IF ( .NOT. do_average_mrt  .AND. ALLOCATED( mrt_av_grid ) )  THEN
    508468
    509469                IF ( mrt_include_sw )  THEN
     
    511471                      ( ( human_absorb * mrtinsw(:) +                          &
    512472                      mrtinlw(:) ) /                                           &
    513                       ( human_emiss * sigma_sb ) ) ** .25_wp - degc_to_k
     473                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
    514474                ELSE
    515475                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
    516476                      ( mrtinlw(:) /                                           &
    517                       ( human_emiss * sigma_sb ) ) ** .25_wp - degc_to_k
     477                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
    518478                ENDIF
    519479             ENDIF
     
    523483!--          Only continue if the current index is the one to trigger the input
    524484!--          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
    533501                DO  i = nxl, nxr
    534502                   DO  j = nys, nyn
     
    540508             ENDIF
    541509
    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
    543516                DO  i = nxl, nxr
    544517                   DO  j = nys, nyn
     
    550523             ENDIF
    551524
    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
    553531                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
    554532                   DO  j = nysg, nyng
     
    560538             ENDIF
    561539
    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
    563546                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
    564547                   DO  j = nysg, nyng
     
    570553             ENDIF
    571554
    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
    573561                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
    574562                   DO  j = nysg, nyng
     
    580568             ENDIF
    581569
    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
    583576
    584577                IF ( mrt_include_sw )  THEN
    585578                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
    586579                      ( ( human_absorb * mrtinsw(:) +                          &
    587                       mrtinlw(:) ) /                             &  ! human_emiss * mrtinlw(:) /
    588                       ( human_emiss * sigma_sb ) ) ** .25_wp - degc_to_k
     580                      mrtinlw(:) ) /                                           &
     581                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
    589582                ELSE
    590583                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
    591                       ( mrtinlw(:) /                             &  ! ( human_emiss * mrtinlw(:) /
    592                       ( human_emiss * sigma_sb ) ) ** .25_wp - degc_to_k
     584                      ( mrtinlw(:) /                                           &
     585                      ( human_emiss * sigma_sb ) )**.25_wp - degc_to_k
    593586                ENDIF
    594587             ENDIF
     
    596589!--       This is a cumulated dose. No mode == 'average' for this quantity.
    597590          CASE ( 'uvem_vitd3dose*' )
    598              IF ( ALLOCATED( vitd3_exposure_av ) ) THEN
     591             IF ( ALLOCATED( vitd3_exposure_av ) )  THEN
    599592                DO  i = nxlg, nxrg
    600593                   DO  j = nysg, nyng
     
    618611!--          that case do_average_mrt will be .TRUE. leading to a double-
    619612!--          averaging.
    620              IF ( .NOT. do_average_mrt .AND. ALLOCATED( mrt_av_grid ) )  THEN
     613             IF ( .NOT. do_average_mrt  .AND. ALLOCATED( mrt_av_grid ) )  THEN
    621614                mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp )
    622615             ENDIF
     
    625618!
    626619!--          Only continue if update index, see above
    627              IF ( average_trigger_perct .AND.                                  &
    628                 TRIM( variable ) /= 'bio_perct*' ) RETURN
    629              IF ( average_trigger_utci .AND.                                   &
    630                 TRIM( variable ) /= 'bio_utci*' ) RETURN
    631              IF ( average_trigger_pet  .AND.                                   &
    632                 TRIM( variable ) /= 'bio_pet*' ) RETURN
    633 
    634              IF ( ALLOCATED( pt_av ) .AND. do_average_theta ) THEN
     620             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
    635628                DO  i = nxl, nxr
    636629                   DO  j = nys, nyn
     
    643636             ENDIF
    644637
    645              IF ( ALLOCATED( q_av ) .AND. do_average_q ) THEN
     638             IF ( ALLOCATED( q_av )  .AND.  do_average_q ) THEN
    646639                DO  i = nxl, nxr
    647640                   DO  j = nys, nyn
     
    654647             ENDIF
    655648
    656              IF ( ALLOCATED( u_av ) .AND. do_average_u ) THEN
     649             IF ( ALLOCATED( u_av )  .AND.  do_average_u ) THEN
    657650                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
    658651                   DO  j = nysg, nyng
     
    665658             ENDIF
    666659
    667              IF ( ALLOCATED( v_av ) .AND. do_average_v ) THEN
     660             IF ( ALLOCATED( v_av )  .AND.  do_average_v ) THEN
    668661                DO  i = nxlg, nxrg
    669662                   DO  j = nysg, nyng
     
    676669             ENDIF
    677670
    678              IF ( ALLOCATED( w_av ) .AND. do_average_w ) THEN
     671             IF ( ALLOCATED( w_av )  .AND.  do_average_w ) THEN
    679672                DO  i = nxlg, nxrg
    680673                   DO  j = nysg, nyng
     
    687680             ENDIF
    688681
    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
    696686
    697687!
     
    730720    SELECT CASE ( TRIM( var ) )
    731721!
    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).
    737728       CASE ( 'bio_mrt' )
    738729          unit = 'degree_C'
     
    746737          unit = 'degree_C'
    747738          thermal_comfort = .TRUE.
    748           IF ( j == 0 ) THEN                !< if instantaneous input
     739          IF ( j == 0 )  THEN                !< if instantaneous input
    749740             do_calculate_perct = .TRUE.
    750741             IF ( .NOT. ALLOCATED( perct ) )  THEN
     
    763754          unit = 'degree_C'
    764755          thermal_comfort = .TRUE.
    765           IF ( j == 0 ) THEN
     756          IF ( j == 0 )  THEN
    766757             do_calculate_utci = .TRUE.
    767758             IF ( .NOT. ALLOCATED( utci ) )  THEN
     
    780771          unit = 'degree_C'
    781772          thermal_comfort = .TRUE.
    782           IF ( j == 0 ) THEN
     773          IF ( j == 0 )  THEN
    783774             do_calculate_pet = .TRUE.
    784775             IF ( .NOT. ALLOCATED( pet ) )  THEN
     
    838829!
    839830!--    Further checks if thermal comfort output is desired.
    840     IF ( thermal_comfort .AND. unit == 'degree_C' )  THEN
     831    IF ( thermal_comfort  .AND. unit == 'degree_C' )  THEN
    841832
    842833!
    843834!--    Break if required modules "radiation" and "humidity" are not avalable.
    844        IF ( .NOT.  radiation ) THEN
     835       IF ( .NOT.  radiation )  THEN
    845836          message_string = 'output of "' // TRIM( var ) // '" require'         &
    846837                           // 's radiation = .TRUE.'
     
    855846          unit = 'illegal'
    856847       ENDIF
    857        IF ( mrt_nlevels == 0 ) THEN
     848       IF ( mrt_nlevels == 0 )  THEN
    858849          message_string = 'output of "' // TRIM( var ) // '" require'         &
    859850                           // 's mrt_nlevels > 0'
     
    934925              j = mrtbl(iy,l)
    935926              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 ) CYCLE
     927              IF ( k < nzb_do  .OR.  k > nzt_do  .OR.  j < nys  .OR.           &
     928                 j > nyn  .OR.  i < nxl  .OR.  i > nxr ) CYCLE
    938929              IF ( av == 0 )  THEN
    939930                 IF ( mrt_include_sw )  THEN
    940931                    local_pf(i,j,k) = ( ( human_absorb * mrtinsw(l) +          &
    941932                                    mrtinlw(l) ) /                             &
    942                                     ( human_emiss * sigma_sb ) ) ** .25_wp -   &
     933                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
    943934                                    degc_to_k
    944935                 ELSE
    945936                    local_pf(i,j,k) = ( mrtinlw(l)  /                          &
    946                                     ( human_emiss * sigma_sb ) ) ** .25_wp -   &
     937                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
    947938                                    degc_to_k
    948939                 ENDIF
     
    968959                 ENDDO
    969960              ENDDO
    970            END IF
     961           ENDIF
    971962
    972963
     
    986977                 ENDDO
    987978              ENDDO
    988            END IF
     979           ENDIF
    989980
    990981
     
    1004995                 ENDDO
    1005996              ENDDO
    1006            END IF
     997           ENDIF
    1007998
    1008999!
     
    10871078               j = mrtbl(iy,l)
    10881079               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 ) CYCLE
     1080               IF ( k < nzb_do  .OR.  k > nzt_do  .OR.  j < nys  .OR.          &
     1081                  j > nyn  .OR.  i < nxl  .OR.  i > nxr ) CYCLE
    10911082               IF ( av == 0 )  THEN
    10921083                  IF ( mrt_include_sw )  THEN
    10931084                     local_pf(i,j,k) = REAL( ( ( human_absorb * mrtinsw(l) +   &
    10941085                                    mrtinlw(l) ) /                             &
    1095                                     ( human_emiss * sigma_sb ) ) ** .25_wp -   &
     1086                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
    10961087                                    degc_to_k, KIND = sp )
    10971088                  ELSE
    10981089                     local_pf(i,j,k) = REAL( ( mrtinlw(l) /                    &
    1099                                     ( human_emiss * sigma_sb ) ) ** .25_wp -   &
     1090                                    ( human_emiss * sigma_sb ) )**.25_wp -     &
    11001091                                    degc_to_k, KIND = sp )
    11011092                  ENDIF
     
    11471138    is2d = ( var(l-1:l) == 'xy' )
    11481139
    1149     IF ( var(1:4) == 'bio_' ) THEN
     1140    IF ( var(1:4) == 'bio_' )  THEN
    11501141       found  = .TRUE.
    11511142       grid_x = 'x'
    11521143       grid_y = 'y'
    11531144       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' ) THEN
     1145       IF ( is2d  .AND.  var(1:7) /= 'bio_mrt' )  grid_z = 'zu1'
     1146    ENDIF
     1147
     1148    IF ( is2d  .AND.  var(1:4) == 'uvem' )  THEN
    11581149       grid_x = 'x'
    11591150       grid_y = 'y'
     
    12411232        ONLY: message_string
    12421233
    1243     IF ( .NOT. radiation_interactions ) THEN
     1234    IF ( .NOT. radiation_interactions )  THEN
    12441235       message_string = 'The mrt calculation requires ' //                     &
    12451236                        'enabled radiation_interactions but it ' //            &
     
    14851476!
    14861477!-- We need to differentiate if averaged input is desired (av == .TRUE.) or not.
    1487     IF ( av ) THEN
     1478    IF ( av )  THEN
    14881479!
    14891480!-- Make sure tmrt_av_grid is present and initialize with the fill value
     
    15051496             k = mrtbl(iz,l)
    15061497             IF ( k - get_topography_top_index_ji( j, i, 's' ) ==              &
    1507                    bio_cell_level + 1_iwp) THEN
     1498                   bio_cell_level + 1_iwp)  THEN
    15081499!
    15091500!-- Averaging was done before, so we can just copy the result here
     
    15311522          k = mrtbl(iz,l)
    15321523          IF ( k - get_topography_top_index_ji( j, i, 's' ) ==                 &
    1533                 bio_cell_level + 1_iwp) THEN
     1524                bio_cell_level + 1_iwp)  THEN
    15341525             IF ( mrt_include_sw )  THEN
    15351526                 tmrt_grid(j,i) = ( ( human_absorb * mrtinsw(l) +              &
    15361527                                  mrtinlw(l) )  /                              &
    1537                                   ( human_emiss * sigma_sb ) ) ** .25_wp -     &
     1528                                  ( human_emiss * sigma_sb ) )**.25_wp -       &
    15381529                                  degc_to_k
    15391530             ELSE
    15401531                 tmrt_grid(j,i) = ( mrtinlw(l)  /                              &
    1541                                   ( human_emiss * sigma_sb ) ) ** .25_wp -     &
     1532                                  ( human_emiss * sigma_sb ) )**.25_wp -       &
    15421533                                  degc_to_k
    15431534             ENDIF
     
    15851576!-- Avoid non-representative horizontal u and v of 0.0 m/s too close to ground
    15861577    k_wind = k
    1587     IF ( bio_cell_level < 1_iwp ) THEN
     1578    IF ( bio_cell_level < 1_iwp )  THEN
    15881579       k_wind = k + 1_iwp
    15891580    ENDIF
    15901581!
    15911582!-- Determine local values:
    1592     IF ( average_input ) THEN
     1583    IF ( average_input )  THEN
    15931584!
    15941585!--    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_k
     1586       ta = pt_av(k,j,i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
    15961587
    15971588       vp = bio_fill_value
    1598        IF ( humidity .AND. ALLOCATED( q_av ) ) THEN
    1599           vp = q_av(k, j, i)
     1589       IF ( humidity  .AND.  ALLOCATED( q_av ) ) THEN
     1590          vp = q_av(k,j,i)
    16001591       ENDIF
    16011592
     
    16061597!
    16071598!-- 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_k
     1599       ta = pt(k,j,i) - ( 0.0098_wp * dz(1) * (  k + .5_wp ) ) - degc_to_k
    16091600
    16101601       vp = bio_fill_value
    1611        IF ( humidity ) THEN
    1612           vp = q(k, j, i)
     1602       IF ( humidity )  THEN
     1603          vp = q(k,j,i)
    16131604       ENDIF
    16141605
     
    16251616!-- The magnus formula is limited to temperatures up to 333.15 K to
    16261617!   avoid negative values of vp_sat
    1627     IF ( vp > -998._wp ) THEN
     1618    IF ( vp > -998._wp )  THEN
    16281619       vp_sat = 0.01_wp * magnus( MIN( ta + degc_to_k, 333.15_wp ) )
    16291620       vp  = vp * pair / ( vp + 0.622_wp )
    1630        IF ( vp > vp_sat ) vp = vp_sat
     1621       IF ( vp > vp_sat )  vp = vp_sat
    16311622    ENDIF
    16321623!
    16331624!-- local mtr value at [i,j]
    16341625    tmrt = bio_fill_value  !< this can be a valid result (e.g. for inside some ostacle)
    1635     IF ( .NOT. average_input ) THEN
     1626    IF ( .NOT. average_input )  THEN
    16361627!
    16371628!--    Use MRT from RTM precalculated in tmrt_grid
     
    16741665    CALL bio_calculate_mrt_grid ( av )
    16751666
    1676     DO i = nxl, nxr
    1677        DO j = nys, nyn
     1667    DO  i = nxl, nxr
     1668       DO  j = nys, nyn
    16781669!
    16791670!--       Determine local input conditions
     
    16891680          perct_ij = bio_fill_value   !< some obstacle
    16901681          utci_ij  = bio_fill_value
    1691           IF ( .NOT. ( tmrt_ij <= -998._wp .OR. vp <= -998._wp ) ) THEN
     1682          IF ( .NOT. ( tmrt_ij <= -998._wp  .OR.  vp <= -998._wp ) ) THEN
    16921683!
    16931684!--          Calculate static thermal indices based on local tmrt
    16941685             clo = bio_fill_value
    16951686
    1696              IF ( do_calculate_perct .OR. do_calculate_perct_av ) THEN
     1687             IF ( do_calculate_perct  .OR.  do_calculate_perct_av ) THEN
    16971688!
    16981689!--          Estimate local perceived temperature
     
    17011692             ENDIF
    17021693
    1703              IF ( do_calculate_utci  .OR. do_calculate_utci_av ) THEN
     1694             IF ( do_calculate_utci  .OR.  do_calculate_utci_av ) THEN
    17041695!
    17051696!--          Estimate local universal thermal climate index
     
    17081699             ENDIF
    17091700
    1710              IF ( do_calculate_pet .OR. do_calculate_pet_av ) THEN
     1701             IF ( do_calculate_pet  .OR.  do_calculate_pet_av ) THEN
    17111702!
    17121703!--          Estimate local physiologically equivalent temperature
    17131704                CALL calculate_pet_static( ta, vp, ws, tmrt_ij, pair, pet_ij )
    17141705             ENDIF
    1715           END IF
    1716 
    1717 
    1718           IF ( av ) THEN
     1706          ENDIF
     1707
     1708
     1709          IF ( av )  THEN
    17191710!
    17201711!--          Write results for selected averaged indices
    17211712             IF ( do_calculate_perct_av )  THEN
    17221713                perct_av(j, i) = perct_ij
    1723              END IF
    1724              IF ( do_calculate_utci_av ) THEN
     1714             ENDIF
     1715             IF ( do_calculate_utci_av )  THEN
    17251716                utci_av(j, i) = utci_ij
    1726              END IF
    1727              IF ( do_calculate_pet_av ) THEN
     1717             ENDIF
     1718             IF ( do_calculate_pet_av )  THEN
    17281719                pet_av(j, i)  = pet_ij
    1729              END IF
     1720             ENDIF
    17301721          ELSE
    17311722!
     
    17331724             IF ( do_calculate_perct )  THEN
    17341725                perct(j, i) = perct_ij
    1735              END IF
    1736              IF ( do_calculate_utci ) THEN
     1726             ENDIF
     1727             IF ( do_calculate_utci )  THEN
    17371728                utci(j, i) = utci_ij
    1738              END IF
    1739              IF ( do_calculate_pet ) THEN
     1729             ENDIF
     1730             IF ( do_calculate_pet )  THEN
    17401731                pet(j, i)  = pet_ij
    1741              END IF
    1742           END IF
    1743 
    1744        END DO
    1745     END DO
     1732             ENDIF
     1733          ENDIF
     1734
     1735       ENDDO
     1736    ENDDO
    17461737
    17471738 END SUBROUTINE bio_calculate_thermal_index_maps
     
    17821773!
    17831774!-- return immediatelly if nothing to do!
    1784     IF ( .NOT. thermal_comfort ) THEN
     1775    IF ( .NOT. thermal_comfort )  THEN
    17851776        RETURN
    17861777    ENDIF
    17871778!
    17881779!-- If clo equals the initial value, this is the initial call
    1789     IF ( clo <= -998._wp ) THEN
     1780    IF ( clo <= -998._wp )  THEN
    17901781!
    17911782!--    Initialize instationary perceived temperature with personalized
     
    18811872!
    18821873!-- 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 ) ) THEN
     1874    IF ( ( d_tmrt > 70._wp )  .OR.  ( d_tmrt < -30._wp )  .OR.                 &
     1875       ( vp >= 50._wp ) )  THEN
    18851876       utci_ij = bio_fill_value
    18861877       RETURN
     
    18881879!
    18891880!-- Apply eq. 2 in Broede et al. (2012) for ta out of bounds
    1890     IF ( ta > 50._wp ) THEN
     1881    IF ( ta > 50._wp )  THEN
    18911882       offset = ta - 50._wp
    18921883       ta = 50._wp
    18931884    ENDIF
    1894     IF ( ta < -50._wp ) THEN
     1885    IF ( ta < -50._wp )  THEN
    18951886       offset = ta + 50._wp
    18961887       ta = -50._wp
     
    19001891!-- humidity values below 0.5 m/s or 5%, respectively, the
    19011892!-- user is advised to use the lower bounds for the calculations.
    1902     IF ( va < 0.5_wp ) va = 0.5_wp
    1903     IF ( va > 17._wp ) va = 17._wp
     1893    IF ( va < 0.5_wp )  va = 0.5_wp
     1894    IF ( va > 17._wp )  va = 17._wp
    19041895
    19051896!
     
    22292220!
    22302221!-- decision: firstly calculate for winter or summer clothing
    2231     IF ( ta <= 10._wp ) THEN
     2222    IF ( ta <= 10._wp )  THEN
    22322223!
    22332224!--    First guess: winter clothing insulation: cold stress
     
    22362227       pmv_w = pmva
    22372228
    2238        IF ( pmva > 0._wp ) THEN
     2229       IF ( pmva > 0._wp )  THEN
    22392230!
    22402231!--       Case summer clothing insulation: heat load ?
     
    22432234             top )
    22442235          pmv_s = pmva
    2245           IF ( pmva <= 0._wp ) THEN
     2236          IF ( pmva <= 0._wp )  THEN
    22462237!
    22472238!--          Case: comfort achievable by varying clothing insulation
     
    22492240             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
    22502241                pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
    2251              IF ( ncount < 0_iwp ) THEN
     2242             IF ( ncount < 0_iwp )  THEN
    22522243                nerr = -1_iwp
    22532244                RETURN
    22542245             ENDIF
    2255           ELSE IF ( pmva > 0.06_wp ) THEN
     2246          ELSE IF ( pmva > 0.06_wp )  THEN
    22562247             clo = 0.5_wp
    22572248             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
    22582249                           pmva, top )
    22592250          ENDIF
    2260        ELSE IF ( pmva < -0.11_wp ) THEN
     2251       ELSE IF ( pmva < -0.11_wp )  THEN
    22612252          clo = 1.75_wp
    22622253          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
     
    22702261       pmv_s = pmva
    22712262
    2272        IF ( pmva < 0._wp ) THEN
     2263       IF ( pmva < 0._wp )  THEN
    22732264!
    22742265!--       Case winter clothing insulation: cold stress ?
     
    22782269          pmv_w = pmva
    22792270
    2280           IF ( pmva >= 0._wp ) THEN
     2271          IF ( pmva >= 0._wp )  THEN
    22812272!
    22822273!--          Case: comfort achievable by varying clothing insulation
     
    22842275             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
    22852276                               pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
    2286              IF ( ncount < 0_iwp ) THEN
     2277             IF ( ncount < 0_iwp )  THEN
    22872278                nerr = -1_iwp
    22882279                RETURN
    22892280             ENDIF
    2290           ELSE IF ( pmva < -0.11_wp ) THEN
     2281          ELSE IF ( pmva < -0.11_wp )  THEN
    22912282             clo = 1.75_wp
    22922283             CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta,           &
    22932284                           pmva, top )
    22942285          ENDIF
    2295        ELSE IF ( pmva > 0.06_wp ) THEN
     2286       ELSE IF ( pmva > 0.06_wp )  THEN
    22962287          clo = 0.5_wp
    22972288          CALL fanger ( ta, tmrt, vp, ws, pair, clo, actlev, eta, pmva,        &
     
    23052296    CALL perct_regression ( pmva, clo, perct_ij )
    23062297    ptc = perct_ij
    2307     IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
     2298    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp ) THEN
    23082299!
    23092300!--    Adjust for cold conditions according to Gagge 1986
    23102301       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
    2311        IF ( nerr_cold > 0_iwp ) nerr = -5_iwp
     2302       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
    23122303       pmvs = pmva - d_pmv
    2313        IF ( pmvs > -0.11_wp ) THEN
     2304       IF ( pmvs > -0.11_wp )  THEN
    23142305          d_pmv  = 0._wp
    23152306          pmvs   = -0.11_wp
     
    23192310!     clo_fanger = clo
    23202311    clon = clo
    2321     IF ( clo > 0.5_wp .AND. perct_ij <= 8.73_wp ) THEN
     2312    IF ( clo > 0.5_wp  .AND.  perct_ij <= 8.73_wp ) THEN
    23222313!
    23232314!--    Required clothing insulation (ireq) is exclusively defined for
     
    23302321    sultrieness    = .FALSE.
    23312322    d_std = -99._wp
    2332     IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
     2323    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp ) THEN
    23332324!
    23342325!--    Adjust for warm/humid conditions according to Gagge 1986
     
    23372328       pmvs   = pmva + d_pmv
    23382329       CALL perct_regression ( pmvs, clo, perct_ij )
    2339        IF ( sult_lim < 99._wp ) THEN
    2340           IF ( (perct_ij - ptc) > sult_lim ) sultrieness = .TRUE.
     2330       IF ( sult_lim < 99._wp )  THEN
     2331          IF ( (perct_ij - ptc) > sult_lim )  sultrieness = .TRUE.
    23412332!
    23422333!--       Set factor to threshold for sultriness
    2343           IF ( dgtcstd /= 0_iwp ) THEN
     2334          IF ( dgtcstd /= 0_iwp )  THEN
    23442335             d_std = ( ( perct_ij - ptc ) - dgtcm ) / dgtcstd
    23452336          ENDIF
     
    23692360
    23702361
    2371     IF ( ta < 0._wp ) THEN
     2362    IF ( ta < 0._wp )  THEN
    23722363!
    23732364!--    ta  < 0 (degC): water vapour pressure over ice
     
    24512442    clo_upper   = wclo
    24522443    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 ) ) THEN
     2444    IF ( ( y_lower > 0._wp  .AND.  y_upper < 0._wp )  .OR.                     &
     2445         ( y_lower < 0._wp  .AND.  y_upper > 0._wp ) ) THEN
    24552446       x_lower  = clo_lower
    24562447       x_upper  = clo_upper
    24572448       x_ridder = guess_0
    24582449
    2459        DO j = 1_iwp, max_iteration
     2450       DO  j = 1_iwp, max_iteration
    24602451          x_average = 0.5_wp * ( x_lower + x_upper )
    24612452          CALL fanger ( ta, tmrt, vp, ws, pair, x_average, actlev, eta,        &
    24622453                        y_average, top )
    24632454          sroot = SQRT ( y_average**2 - y_lower * y_upper )
    2464           IF ( sroot == 0._wp ) THEN
     2455          IF ( sroot == 0._wp )  THEN
    24652456             clo_res = x_average
    24662457             nerr = j
     
    24692460          x_new = x_average + ( x_average - x_lower ) *                        &
    24702461                      ( SIGN ( 1._wp, y_lower - y_upper ) * y_average / sroot )
    2471           IF ( ABS ( x_new - x_ridder ) <= eps ) THEN
     2462          IF ( ABS ( x_new - x_ridder ) <= eps )  THEN
    24722463             clo_res = x_ridder
    24732464             nerr       = j
     
    24772468          CALL fanger ( ta, tmrt, vp, ws, pair, x_ridder, actlev, eta,         &
    24782469                        y_new, top )
    2479           IF ( y_new == 0._wp ) THEN
     2470          IF ( y_new == 0._wp )  THEN
    24802471             clo_res = x_ridder
    24812472             nerr       = j
    24822473             RETURN
    24832474          ENDIF
    2484           IF ( SIGN ( y_average, y_new ) /= y_average ) THEN
     2475          IF ( SIGN ( y_average, y_new ) /= y_average )  THEN
    24852476             x_lower = x_average
    24862477             y_lower = y_average
    24872478             x_upper  = x_ridder
    24882479             y_upper  = y_new
    2489           ELSE IF ( SIGN ( y_lower, y_new ) /= y_lower ) THEN
     2480          ELSE IF ( SIGN ( y_lower, y_new ) /= y_lower )  THEN
    24902481             x_upper  = x_ridder
    24912482             y_upper  = y_new
    2492           ELSE IF ( SIGN ( y_upper, y_new ) /= y_upper ) THEN
     2483          ELSE IF ( SIGN ( y_upper, y_new ) /= y_upper )  THEN
    24932484             x_lower = x_ridder
    24942485             y_lower = y_new
     
    25002491             RETURN
    25012492          ENDIF
    2502           IF ( ABS ( x_upper - x_lower ) <= eps ) THEN
     2493          IF ( ABS ( x_upper - x_lower ) <= eps )  THEN
    25032494             clo_res = x_ridder
    25042495             nerr       = j
     
    25112502       clo_res = y_new
    25122503       RETURN
    2513     ELSE IF ( y_lower == 0. ) THEN
     2504    ELSE IF ( y_lower == 0. )  THEN
    25142505       x_ridder = clo_lower
    2515     ELSE IF ( y_upper == 0. ) THEN
     2506    ELSE IF ( y_upper == 0. )  THEN
    25162507       x_ridder = clo_upper
    25172508    ELSE
     
    25412532    REAL(wp), INTENT ( OUT ) ::  perct_ij   !< perct (degC) corresponding to given PMV / clo
    25422533
    2543     IF ( pmv <= -0.11_wp ) THEN
     2534    IF ( pmv <= -0.11_wp )  THEN
    25442535       perct_ij = 5.805_wp + 12.6784_wp * pmv
    25452536    ELSE
    2546        IF ( pmv >= + 0.01_wp ) THEN
     2537       IF ( pmv >= + 0.01_wp )  THEN
    25472538          perct_ij = 16.826_wp + 6.163_wp * pmv
    25482539       ELSE
     
    26082599!-- Clo must be > 0. to avoid div. by 0!
    26092600    clo = in_clo
    2610     IF ( clo <= 0._wp ) clo = .001_wp
     2601    IF ( clo <= 0._wp )  clo = .001_wp
    26112602!
    26122603!-- f_cl = Increase in surface due to clothing
     
    26152606!-- Case of free convection (ws < 0.1 m/s ) not considered
    26162607    ws = in_ws
    2617     IF ( ws < .1_wp ) THEN
     2608    IF ( ws < .1_wp )  THEN
    26182609       ws = .1_wp
    26192610    ENDIF
     
    26382629!-- Newton-approximation with air temperature as initial guess
    26392630    t_clothing = ta
    2640     DO i = 1, 3
     2631    DO  i = 1, 3
    26412632       t_clothing = t_clothing - ( ( t_clothing + degc_to_k )**4 + t_clothing  &
    26422633          * dc - gc ) / ( 4._wp * ( t_clothing + degc_to_k )**3 + dc )
     
    26602651    z5 = 3.96e-8_wp * f_cl * ( ( t_clothing + degc_to_k )**4 - ( tmrt +        &
    26612652       degc_to_k )**4 )
    2662     IF ( ABS ( t_clothing - tmrt ) > 0._wp ) THEN
     2653    IF ( ABS ( t_clothing - tmrt ) > 0._wp )  THEN
    26632654       hr = z5 / f_cl / ( t_clothing - tmrt )
    26642655    ELSE
     
    27602751!
    27612752!-- Test for compliance with regression range
    2762     IF ( pmva < -1.0_wp .OR. pmva > 7.0_wp ) THEN
     2753    IF ( pmva < -1.0_wp  .OR.  pmva > 7.0_wp ) THEN
    27632754       nerr = -2_iwp
    27642755    ELSE
     
    27722763    p10  = 0.05_wp * svp_ta
    27732764    p95  = 1.00_wp * svp_ta
    2774     IF ( vp >= p10 .AND. vp <= p95 ) THEN
     2765    IF ( vp >= p10  .AND.  vp <= p95 ) THEN
    27752766       pa = vp
    27762767    ELSE
    27772768       nerr = -3_iwp
    2778        IF ( vp < p10 ) THEN
     2769       IF ( vp < p10 )  THEN
    27792770!
    27802771!--       Due to conditions of regression: r.H. >= 5 %
     
    27862777       ENDIF
    27872778    ENDIF
    2788     IF ( pa > 0._wp ) THEN
     2779    IF ( pa > 0._wp )  THEN
    27892780!
    27902781!--    Natural logarithm of pa
     
    27962787!-- Ratio actual water vapour pressure to that of a r.H. of 50 %
    27972788    pa_p50   = 0.5_wp * svp_ta
    2798     IF ( pa_p50 > 0._wp .AND. pa > 0._wp ) THEN
     2789    IF ( pa_p50 > 0._wp  .AND.  pa > 0._wp ) THEN
    27992790       dapa   = apa - LOG ( pa_p50 )
    28002791       pa_p50 = pa / pa_p50
     
    28052796!
    28062797!-- Square root of wind velocity
    2807     IF ( ws >= 0._wp ) THEN
     2798    IF ( ws >= 0._wp )  THEN
    28082799       sqvel = SQRT ( ws )
    28092800    ELSE
     
    28162807!-- Select the valid regression coefficients
    28172808    nreg = INT ( pmv )
    2818     IF ( nreg < 0_iwp ) THEN
     2809    IF ( nreg < 0_iwp )  THEN
    28192810!
    28202811!--    value of the FUNCTION in the case pmv <= -1
     
    28232814    ENDIF
    28242815    weight = MOD ( pmv, 1._wp )
    2825     IF ( weight < 0._wp ) weight = 0._wp
    2826     IF ( nreg > 5_iwp ) THEN
     2816    IF ( weight < 0._wp )  weight = 0._wp
     2817    IF ( nreg > 5_iwp )  THEN
    28272818       ! nreg=6
    28282819       nreg  = 5_iwp
    28292820       weight   = pmv - 5._wp
    28302821       weight2  = pmv - 6._wp
    2831        IF ( weight2 > 0_iwp ) THEN
     2822       IF ( weight2 > 0_iwp )  THEN
    28322823          weight = ( weight - weight2 ) / weight
    28332824       ENDIF
     
    28362827!-- Regression valid for 0. <= pmv <= 6.
    28372828    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)
    28472838
    28482839    dpmv_2 = 0._wp
    2849     IF ( nreg < 6_iwp ) THEN
     2840    IF ( nreg < 6_iwp )  THEN
    28502841       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)
    28602851    ENDIF
    28612852!
     
    28632854    deltapmv = ( 1._wp - weight ) * dpmv_1 + weight * dpmv_2
    28642855    pmvs = pmva + deltapmv
    2865     IF ( ( pmvs ) < 0._wp ) THEN
     2856    IF ( ( pmvs ) < 0._wp )  THEN
    28662857!
    28672858!--    Prevent negative pmv* due to problems with clothing insulation
    28682859       nerr = -4_iwp
    2869        IF ( pmvs > -0.11_wp ) THEN
     2860       IF ( pmvs > -0.11_wp )  THEN
    28702861!
    28712862!--       Threshold from FUNCTION perct_regression for winter clothing insulation
     
    29252916    dperctstd = 999999._wp
    29262917
    2927     IF ( perct_ij < 16.826_wp .OR. perct_ij > 56._wp ) THEN
     2918    IF ( perct_ij < 16.826_wp  .OR.  perct_ij > 56._wp ) THEN
    29282919!
    29292920!--    Unallowed value of classical perct!
     
    29412932!-- Value of the FUNCTION
    29422933    sultr_res = dperctm + faktor * dperctstd
    2943     IF ( ABS ( sultr_res ) > 99._wp ) sultr_res = +99._wp
     2934    IF ( ABS ( sultr_res ) > 99._wp )  sultr_res = +99._wp
    29442935
    29452936 END SUBROUTINE calc_sultr
     
    29952986    dtmrt          = tmrt - ta
    29962987    sqrt_ws        = ws
    2997     IF ( sqrt_ws < 0.10_wp ) THEN
     2988    IF ( sqrt_ws < 0.10_wp )  THEN
    29982989       sqrt_ws = 0.10_wp
    29992990    ELSE
     
    30012992    ENDIF
    30022993
    3003     DO i = 1, 3
    3004        delta_cold (i) = 0._wp
    3005        IF ( i < 3 ) THEN
    3006           pmv_cross (i) = pmvc
     2994    DO  i = 1, 3
     2995       delta_cold(i) = 0._wp
     2996       IF ( i < 3 )  THEN
     2997          pmv_cross(i) = pmvc
    30072998       ENDIF
    30082999    ENDDO
    30093000
    3010     DO i = 1, 3
     3001    DO  i = 1, 3
    30113002!
    30123003!--    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) *                  &
    30143005                  sqrt_ws + coeff(i,5)*dtmrt
    3015        delta_cold(i) = reg_a(i) + coeff(i, 2) * pmvc
     3006       delta_cold(i) = reg_a(i) + coeff(i,2) * pmvc
    30163007    ENDDO
    30173008!
    30183009!-- Intersection points of regression lines in terms of Fanger's PMV
    3019     DO i = 1, 2
    3020        r_nenner = coeff (i, 2) - coeff (i + 1, 2)
    3021        IF ( ABS ( r_nenner ) > 0.00001_wp ) THEN
    3022           pmv_cross(i) = ( reg_a (i + 1) - reg_a (i) ) / r_nenner
     3010    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
    30233014       ELSE
    30243015          nerr = 1_iwp
     
    30283019
    30293020    i_bin = 3
    3030     DO i = 1, 2
    3031        IF ( pmva > pmv_cross (i) ) THEN
     3021    DO  i = 1, 2
     3022       IF ( pmva > pmv_cross(i) ) THEN
    30323023          i_bin = i
    30333024          EXIT
     
    30623053!
    30633054!--                                 range_1        range_2        range_3
    3064     DATA (coef(i, 0), i = 1, n_bin) /0.0941540_wp, -0.1506620_wp, -0.0871439_wp/
    3065     DATA (coef(i, 1), i = 1, n_bin) /0.0783162_wp, -1.0612651_wp,  0.1695040_wp/
    3066     DATA (coef(i, 2), i = 1, n_bin) /0.1350144_wp, -1.0049144_wp, -0.0167627_wp/
    3067     DATA (coef(i, 3), i = 1, n_bin) /0.1104037_wp, -0.2005277_wp, -0.0003230_wp/
    3068 
    3069     IF ( pmva <= -2.1226_wp ) THEN
     3055    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
    30703061       i_bin = 3_iwp
    3071     ELSE IF ( pmva <= -1.28_wp ) THEN
     3062    ELSE IF ( pmva <= -1.28_wp )  THEN
    30723063       i_bin = 2_iwp
    30733064    ELSE
     
    30753066    ENDIF
    30763067
    3077     dpmv_adj   = coef( i_bin, n_ca )
     3068    dpmv_adj   = coef(i_bin,n_ca)
    30783069    pmv        = 1._wp
    30793070
    3080     DO i = n_ca + 1, n_ce
     3071    DO  i = n_ca + 1, n_ce
    30813072       pmv      = pmv * pmva
    3082        dpmv_adj = dpmv_adj + coef(i_bin, i) * pmv
     3073       dpmv_adj = dpmv_adj + coef(i_bin,i) * pmv
    30833074    ENDDO
    30843075    RETURN
     
    31163107!
    31173108!-- Regression only defined for perct <= 10 (degC)
    3118     IF ( ireq_neutral < 0.5_wp ) THEN
    3119        IF ( ireq_neutral < 0.48_wp ) THEN
     3109    IF ( ireq_neutral < 0.5_wp )  THEN
     3110       IF ( ireq_neutral < 0.48_wp )  THEN
    31203111          nerr = 1_iwp
    31213112       ENDIF
     
    31253116!-- Minimal required clothing insulation: maximal acceptable body cooling
    31263117    ireq_minimal = 1.26_wp - 0.0588_wp * top02
    3127     IF ( nerr > 0_iwp ) THEN
     3118    IF ( nerr > 0_iwp )  THEN
    31283119       ireq_minimal = ireq_neutral
    31293120    ENDIF
     
    31523143!
    31533144!-- According to Gehan-George, for children
    3154     IF ( age < 19_iwp ) THEN
    3155        IF ( age < 5_iwp ) THEN
     3145    IF ( age < 19_iwp )  THEN
     3146       IF ( age < 5_iwp )  THEN
    31563147          surf = 0.02667_wp * height**0.42246_wp * weight**0.51456_wp
    31573148          RETURN
    3158        END IF
     3149       ENDIF
    31593150       surf = 0.03050_wp * height**0.35129_wp * weight**0.54375_wp
    31603151       RETURN
    3161     END IF
     3152    ENDIF
    31623153!
    31633154!-- DuBois D, DuBois EF: A formula to estimate the approximate surface area if
     
    32003191
    32013192    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 ) )
    32033194    factor = 1._wp + .004_wp  * ( 30._wp - age )
    32043195    basic_heat_prod = 0.
    3205     IF ( sex == 1_iwp ) THEN
    3206        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 +      &
    32073198                     .01_wp  * ( s - 43.4_wp ) )
    3208     ELSE IF ( sex == 2_iwp ) THEN
    3209        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 +      &
    32103201                    .018_wp * ( s - 42.1_wp ) )
    3211     END IF
     3202    ENDIF
    32123203
    32133204    energy_prod = work + basic_heat_prod
     
    32953286!
    32963287!-- Decision: firstly calculate for winter or summer clothing
    3297     IF ( ta <= 10._wp ) THEN
     3288    IF ( ta <= 10._wp )  THEN
    32983289!
    32993290!--    First guess: winter clothing insulation: cold stress
     
    33043295       pmv_w = pmva
    33053296
    3306        IF ( pmva > 0._wp ) THEN
     3297       IF ( pmva > 0._wp )  THEN
    33073298!
    33083299!--       Case summer clothing insulation: heat load ?           
     
    33123303             t_clothing, storage, dt, pmva )
    33133304          pmv_s = pmva
    3314           IF ( pmva <= 0._wp ) THEN
     3305          IF ( pmva <= 0._wp )  THEN
    33153306!
    33163307!--          Case: comfort achievable by varying clothing insulation
     
    33183309             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta , sclo,     &
    33193310                            pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
    3320              IF ( ncount < 0_iwp ) THEN
     3311             IF ( ncount < 0_iwp )  THEN
    33213312                nerr = -1_iwp
    33223313                RETURN
    33233314             ENDIF
    3324           ELSE IF ( pmva > 0.06_wp ) THEN
     3315          ELSE IF ( pmva > 0.06_wp )  THEN
    33253316             clo = 0.5_wp
    33263317             t_clothing = bio_fill_value
     
    33283319                t_clothing, storage, dt, pmva )
    33293320          ENDIF
    3330        ELSE IF ( pmva < -0.11_wp ) THEN
     3321       ELSE IF ( pmva < -0.11_wp )  THEN
    33313322          clo = 1.75_wp
    33323323          t_clothing = bio_fill_value
     
    33443335       pmv_s = pmva
    33453336
    3346        IF ( pmva < 0._wp ) THEN
     3337       IF ( pmva < 0._wp )  THEN
    33473338!
    33483339!--       Case winter clothing insulation: cold stress ?
     
    33533344          pmv_w = pmva
    33543345
    3355           IF ( pmva >= 0._wp ) THEN
     3346          IF ( pmva >= 0._wp )  THEN
    33563347!
    33573348!--          Case: comfort achievable by varying clothing insulation
     
    33593350             CALL iso_ridder ( ta, tmrt, vp, ws, pair, actlev, eta, sclo,      &
    33603351                               pmv_s, wclo, pmv_w, eps, pmva, top, ncount, clo )
    3361              IF ( ncount < 0_wp ) THEN
     3352             IF ( ncount < 0_wp )  THEN
    33623353                nerr = -1_iwp
    33633354                RETURN
    33643355             ENDIF
    3365           ELSE IF ( pmva < -0.11_wp ) THEN
     3356          ELSE IF ( pmva < -0.11_wp )  THEN
    33663357             clo = 1.75_wp
    33673358             t_clothing = bio_fill_value
     
    33693360                t_clothing, storage, dt, pmva )
    33703361          ENDIF
    3371        ELSE IF ( pmva > 0.06_wp ) THEN
     3362       ELSE IF ( pmva > 0.06_wp )  THEN
    33723363          clo = 0.5_wp
    33733364          t_clothing = bio_fill_value
     
    33823373    CALL perct_regression ( pmva, clo, ipt )
    33833374    ptc = ipt
    3384     IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
     3375    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp ) THEN
    33853376!
    33863377!--    Adjust for cold conditions according to Gagge 1986
    33873378       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
    3388        IF ( nerr_cold > 0_iwp ) nerr = -5_iwp
     3379       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
    33893380       pmvs = pmva - d_pmv
    3390        IF ( pmvs > -0.11_wp ) THEN
     3381       IF ( pmvs > -0.11_wp )  THEN
    33913382          d_pmv  = 0._wp
    33923383          pmvs   = -0.11_wp
     
    33963387!     clo_fanger = clo
    33973388    clon = clo
    3398     IF ( clo > 0.5_wp .AND. ipt <= 8.73_wp ) THEN
     3389    IF ( clo > 0.5_wp  .AND.  ipt <= 8.73_wp ) THEN
    33993390!
    34003391!--    Required clothing insulation (ireq) is exclusively defined for
     
    34073398    sultrieness    = .FALSE.
    34083399    d_std      = -99._wp
    3409     IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
     3400    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp ) THEN
    34103401!
    34113402!--    Adjust for warm/humid conditions according to Gagge 1986
     
    34143405       pmvs   = pmva + d_pmv
    34153406       CALL perct_regression ( pmvs, clo, ipt )
    3416        IF ( sult_lim < 99._wp ) THEN
    3417           IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE.
     3407       IF ( sult_lim < 99._wp )  THEN
     3408          IF ( (ipt - ptc) > sult_lim )  sultrieness = .TRUE.
    34183409       ENDIF
    34193410    ENDIF
     
    34873478    CALL perct_regression ( pmva, clo, ipt )
    34883479
    3489     IF ( clo >= 1.75_wp .AND. pmva <= -0.11_wp ) THEN
     3480    IF ( clo >= 1.75_wp  .AND.  pmva <= -0.11_wp ) THEN
    34903481!
    34913482!--    Adjust for cold conditions according to Gagge 1986
    34923483       CALL dpmv_cold ( pmva, ta, ws, tmrt, nerr_cold, d_pmv )
    3493        IF ( nerr_cold > 0_iwp ) nerr = -5_iwp
     3484       IF ( nerr_cold > 0_iwp )  nerr = -5_iwp
    34943485       pmvs = pmva - d_pmv
    3495        IF ( pmvs > -0.11_wp ) THEN
     3486       IF ( pmvs > -0.11_wp )  THEN
    34963487          d_pmv  = 0._wp
    34973488          pmvs   = -0.11_wp
     
    35053496    sultrieness    = .FALSE.
    35063497    d_std      = -99._wp
    3507     IF ( pmva > 0.06_wp .AND. clo <= 0.5_wp ) THEN
     3498    IF ( pmva > 0.06_wp  .AND.  clo <= 0.5_wp ) THEN
    35083499!
    35093500!--    Adjust for warm/humid conditions according to Gagge 1986
     
    35123503       pmvs   = pmva + d_pmv
    35133504       CALL perct_regression ( pmvs, clo, ipt )
    3514        IF ( sult_lim < 99._wp ) THEN
    3515           IF ( (ipt - ptc) > sult_lim ) sultrieness = .TRUE.
     3505       IF ( sult_lim < 99._wp )  THEN
     3506          IF ( (ipt - ptc) > sult_lim )  sultrieness = .TRUE.
    35163507       ENDIF
    35173508    ENDIF
     
    35833574!-- Clo must be > 0. to avoid div. by 0!
    35843575    clo = in_clo
    3585     IF ( clo < 001._wp ) clo = .001_wp
     3576    IF ( clo < 001._wp )  clo = .001_wp
    35863577!
    35873578!-- Increase in surface due to clothing
     
    35903581!-- Case of free convection (ws < 0.1 m/s ) not considered
    35913582    ws = in_ws
    3592     IF ( ws < .1_wp ) THEN
     3583    IF ( ws < .1_wp )  THEN
    35933584       ws = .1_wp
    35943585    ENDIF
     
    36103601!-- newton-approximation with air temperature as initial guess
    36113602    niter = INT( dt * 10._wp, KIND=iwp )
    3612     IF ( niter < 1 ) niter = 1_iwp
     3603    IF ( niter < 1 )  niter = 1_iwp
    36133604    adjustrate = 1._wp - EXP ( -1._wp * ( 10._wp / time_equil ) * dt )
    3614     IF ( adjustrate >= 1._wp ) adjustrate = 1._wp
     3605    IF ( adjustrate >= 1._wp )  adjustrate = 1._wp
    36153606    adjustrate_cloth = adjustrate * 30._wp
    36163607    t_clothing = t_cloth
    36173608
    3618     IF ( t_cloth <= -998.0_wp ) THEN  ! If initial run
     3609    IF ( t_cloth <= -998.0_wp )  THEN  ! If initial run
    36193610       niter = 3_iwp
    36203611       adjustrate = 1._wp
     
    36233614    ENDIF
    36243615
    3625     DO i = 1, niter
     3616    DO  i = 1, niter
    36263617       t_clothing = t_clothing - adjustrate_cloth * ( ( t_clothing +           &
    36273618          273.2_wp )**4._wp + t_clothing *                                     &
     
    37753766    REAL(wp) ::  vpex
    37763767
    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 *             &
    37783771          ( 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 ) )
    37813773    met = work + met
    3782 
    37833774    int_heat = met * (1._wp - eta)
    37843775!
    37853776!-- Sensible respiration energy
    37863777    tex  = 0.47_wp * ta + 21.0_wp
    3787     rtv  = 1.44_wp * 10._wp ** (-6._wp) * met
     3778    rtv  = 1.44_wp * 10._wp**(-6._wp) * met
    37883779    eres = c_p * (ta - tex) * rtv
    37893780!
    37903781!-- 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 ) )
    37923783    erel = 0.623_wp * l_v / pair * ( vpa - vpex ) * rtv
    37933784!
     
    38893880    LOGICAL ::  skipIncreaseCount   !< iteration control flag
    38903881
    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)
    38963894    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
    39013902    rcl = ( clo / 6.45_wp ) / facl
    39023903    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 )
    39113909    di   = r2 - r1
    3912 !
    3913 !-- Skin temperatur
    3914     DO j = 1, 7
     3910
     3911!
     3912!-- Estimate skin temperatur
     3913    DO  j = 1, 7
    39153914
    39163915       tsk    = 34._wp
     
    39203919       enbal2 = 0._wp
    39213920
    3922        DO i = 1, 100  ! allow for 100 iterations max
     3921       DO  i = 1, 100  ! allow for 100 iterations max
    39233922          acl   = adu * facl + adu * ( fcl - 1._wp )
    39243923          rclo2 = emcl * sigma_sb * ( ( tcl + degc_to_k )**4._wp -             &
    3925             ( tmrt + degc_to_k )** 4._wp ) * feff
     3924            ( tmrt + degc_to_k )**4._wp ) * feff
    39263925          htcl  = 6.28_wp * ht * y * di / ( rcl * LOG( r2 / r1 ) * acl )
    39273926          tsk   = 1._wp / htcl * ( hc * ( tcl - ta ) + rclo2 ) + tcl
     
    39303929          aeff  = adu * feff
    39313930          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 )
    39333932          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 )
    39353934          rsum  = rbare + rclo
    39363935!
     
    39523951          c(9)  = 5.28_wp * adu - c(5) - c(4) * tsk
    39533952          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_wp
     3953                  ( c(5) * tsk - c(0) - 5.28_wp * adu * tsk )
     3954
     3955          IF ( tsk == 36._wp )  tsk = 36.01_wp
    39573956          tcore(7) = c(0) / ( 5.28_wp * adu + c(1) * 6.3_wp / 3600._wp ) + tsk
    39583957          tcore(3) = c(0) / ( 5.28_wp * adu + ( c(1) * 6.3_wp / 3600._wp ) /   &
    39593958            ( 1._wp + 0.5_wp * ( 34._wp - tsk ) ) ) + tsk
    3960           IF ( c(10) >= 0._wp ) THEN
     3959          IF ( c(10) >= 0._wp )  THEN
    39613960             tcore(6) = ( - c(9) - c(10)**0.5_wp ) / ( 2._wp * c(4) )
    39623961             tcore(1) = ( - c(9) + c(10)**0.5_wp ) / ( 2._wp * c(4) )
    3963           END IF
    3964 
    3965           IF ( c(8) >= 0._wp ) THEN
    3966              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) )
    39683967             tcore(4) = c(0) / ( 5.28_wp * adu + c(1) * 1._wp / 40._wp ) + tsk
    3969           END IF
     3968          ENDIF
    39703969!
    39713970!--       Transpiration
     
    39743973          vpts  = 6.11_wp * 10._wp**( 7.45_wp * tsk / ( 235._wp + tsk ) )
    39753974
    3976           IF ( tbody <= 36.6_wp ) swm = 0._wp  !< no need for sweating
     3975          IF ( tbody <= 36.6_wp )  swm = 0._wp  !< no need for sweating
    39773976
    39783977          sw = swm
     
    39833982          wetsk  = eswphy / eswpot
    39843983
    3985           IF ( wetsk > 1._wp ) wetsk = 1._wp
     3984          IF ( wetsk > 1._wp )  wetsk = 1._wp
    39863985!
    39873986!--       Sweat production > evaporation?
    39883987          eswdif = eswphy - eswpot
    39893988
    3990           IF ( eswdif <= 0._wp ) esw = eswpot  !< Limit is evaporation
    3991           IF ( eswdif > 0._wp ) esw = eswphy   !< Limit is sweat production
    3992           IF ( esw  > 0._wp )   esw = 0._wp    !< Sweat can't be evaporated, no more cooling effect
     3989          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
    39933992!
    39943993!--       Diffusion
    3995           rdsk = 0.79_wp * 10._wp ** 7._wp
     3994          rdsk = 0.79_wp * 10._wp**7._wp
    39963995          rdcl = 0._wp
    39973996          ed   = l_v / ( rdsk + rdcl ) * adu * ( 1._wp - wetsk ) * ( vpa -     &
     
    40024001          vb2 = tcore(j) - 36.6_wp
    40034002
    4004           IF ( vb2 < 0._wp ) vb2 = 0._wp
    4005           IF ( vb1 < 0._wp ) vb1 = 0._wp
     4003          IF ( vb2 < 0._wp )  vb2 = 0._wp
     4004          IF ( vb1 < 0._wp )  vb1 = 0._wp
    40064005          vb = ( 6.3_wp + 75._wp * vb2 ) / ( 1._wp + 0.5_wp * vb1 )
    40074006!
     
    40114010!--       Clothing temperature
    40124011          xx = 0.001_wp
    4013           IF ( count1 == 0_iwp ) xx = 1._wp
    4014           IF ( count1 == 1_iwp ) xx = 0.1_wp
    4015           IF ( count1 == 2_iwp ) xx = 0.01_wp
    4016           IF ( count1 == 3_iwp ) xx = 0.001_wp
    4017 
    4018           IF ( enbal > 0._wp ) tcl = tcl + xx
    4019           IF ( enbal < 0._wp ) tcl = tcl - xx
     4012          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
    40204019
    40214020          skipIncreaseCount = .FALSE.
    4022           IF ( ( (enbal <= 0._wp ) .AND. (enbal2 > 0._wp ) ) .OR.              &
    4023              ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) ) THEN
     4021          IF ( ( (enbal <= 0._wp )  .AND.  (enbal2 > 0._wp ) )  .OR.           &
     4022             ( ( enbal >= 0._wp )  .AND.  ( enbal2 < 0._wp ) ) ) THEN
    40244023             skipIncreaseCount = .TRUE.
    40254024          ELSE
    40264025             enbal2 = enbal
    40274026             count3 = count3 + 1_iwp
    4028           END IF
    4029 
    4030           IF ( ( count3 > 200_iwp ) .OR. skipIncreaseCount ) THEN
    4031              IF ( count1 < 3_iwp ) THEN
     4027          ENDIF
     4028
     4029          IF ( ( count3 > 200_iwp )  .OR.  skipIncreaseCount ) THEN
     4030             IF ( count1 < 3_iwp )  THEN
    40324031                count1 = count1 + 1_iwp
    40334032                enbal2 = 0._wp
    40344033             ELSE
    40354034                EXIT
    4036              END IF
    4037           END IF
    4038        END DO
    4039 
    4040        IF ( count1 == 3_iwp ) THEN
     4035             ENDIF
     4036          ENDIF
     4037       ENDDO
     4038
     4039       IF ( count1 == 3_iwp )  THEN
    40414040          SELECT CASE ( j )
    40424041             CASE ( 2, 5)
    4043                 IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND.                     &
    4044                    ( tsk <= 34.050_wp ) ) ) CYCLE
     4042                IF ( .NOT. ( ( tcore(j) >= 36.6_wp )  .AND.                    &
     4043                   ( tsk <= 34.050_wp ) ) )  CYCLE
    40454044             CASE ( 6, 1 )
    40464045                IF ( c(10) < 0._wp ) CYCLE
    4047                 IF ( .NOT. ( ( tcore(j) >= 36.6_wp ) .AND.                     &
    4048                    ( tsk > 33.850_wp ) ) ) CYCLE
     4046                IF ( .NOT. ( ( tcore(j) >= 36.6_wp )  .AND.                    &
     4047                   ( tsk > 33.850_wp ) ) )  CYCLE
    40494048             CASE ( 3 )
    4050                 IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND.                      &
    4051                    ( tsk <= 34.000_wp ) ) ) CYCLE
     4049                IF ( .NOT. ( ( tcore(j) < 36.6_wp )  .AND.                     &
     4050                   ( tsk <= 34.000_wp ) ) )  CYCLE
    40524051             CASE ( 7 )
    4053                 IF ( .NOT. ( ( tcore(j) < 36.6_wp ) .AND.                      &
    4054                    ( tsk > 34.000_wp ) ) ) CYCLE
     4052                IF ( .NOT. ( ( tcore(j) < 36.6_wp )  .AND.                     &
     4053                   ( tsk > 34.000_wp ) ) )  CYCLE
    40554054             CASE default
    40564055          END SELECT
    4057        END IF
    4058 
    4059        IF ( ( j /= 4_iwp ) .AND. ( vb >= 91._wp ) ) CYCLE
    4060        IF ( ( j == 4_iwp ) .AND. ( vb < 89._wp ) ) CYCLE
     4056       ENDIF
     4057
     4058       IF ( ( j /= 4_iwp )  .AND.  ( vb >= 91._wp ) ) CYCLE
     4059       IF ( ( j == 4_iwp )  .AND.  ( vb < 89._wp ) ) CYCLE
    40614060       IF ( vb > 90._wp ) vb = 90._wp
    40624061!
    40634062!--    Loses by water
    40644063       ws = sw * 3600._wp * 1000._wp
    4065        IF ( ws > 2000._wp ) ws = 2000._wp
     4064       IF ( ws > 2000._wp )  ws = 2000._wp
    40664065       wd = ed / l_v * 3600._wp * ( -1000._wp )
    40674066       wr = erel / l_v * 3600._wp * ( -1000._wp )
     
    40704069
    40714070       RETURN
    4072     END DO
     4071    ENDDO
    40734072 END SUBROUTINE heat_exch
    40744073
     
    41314130    enbal2 = 0._wp
    41324131
    4133     DO count1 = 0, 3
    4134        DO i = 1, 125  ! 500 / 4
    4135           hc = 2.67_wp + 6.5_wp * 0.1_wp ** 0.67_wp
    4136           hc = hc * ( pair / po ) ** 0.55_wp
     4132    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
    41374136!
    41384137!--       Radiation
    41394138          aeff  = adu * feff
    41404139          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 )
    41424141          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 )
    41444143          rsum  = rbare + rclo
    41454144!
     
    41564155          tex  = 0.47_wp * pet_ij + 21._wp
    41574156          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 ) )
    41594158          erel = 0.623_wp * l_v / pair * ( 12._wp - vpex ) * rtv
    41604159          ere  = eres + erel
     
    41714170          IF ( enbal > 0._wp )  pet_ij = pet_ij - xx
    41724171          IF ( enbal < 0._wp )  pet_ij = pet_ij + xx
    4173           IF ( ( enbal <= 0._wp ) .AND. ( enbal2 > 0._wp ) ) EXIT
    4174           IF ( ( enbal >= 0._wp ) .AND. ( enbal2 < 0._wp ) ) EXIT
     4172          IF ( ( enbal <= 0._wp )  .AND.  ( enbal2 > 0._wp ) ) EXIT
     4173          IF ( ( enbal >= 0._wp )  .AND.  ( enbal2 < 0._wp ) ) EXIT
    41754174
    41764175          enbal2 = enbal
    4177        END DO
    4178     END DO
     4176       ENDDO
     4177    ENDDO
    41794178 END SUBROUTINE pet_iteration
    41804179
     
    42714270    CALL uvem_solar_position
    42724271     
    4273     IF (sza  .GE.  90) THEN
     4272    IF (sza  .GE.  90)  THEN
    42744273       vitd3_exposure(:,:) = 0.0_wp
    42754274    ELSE
     
    43684367!                   
    43694368! !--        extract obstruction from IBSET-Integer_Array ------------------'
    4370              IF (consider_obstructions ) THEN
     4369             IF (consider_obstructions )  THEN
    43714370                obstruction_temp1 = building_obstruction_f%var_3d(:,j,i)
    4372                 IF (obstruction_temp1(0)  .NE.  9) THEN
     4371                IF (obstruction_temp1(0)  .NE.  9)  THEN
    43734372                   DO  pobi = 0, 44
    43744373                      DO  bi = 0, 7
    4375                          IF ( btest( obstruction_temp1(pobi), bi )  .EQV.  .TRUE.) THEN
     4374                         IF ( btest( obstruction_temp1(pobi), bi )  .EQV.  .TRUE.)  THEN
    43764375                            obstruction_temp2( ( pobi * 8 ) + bi ) = 1
    43774376                         ELSE
     
    43924391         
    43934392             obstruction_direct_beam = obstruction( nint(startpos_saa_float), nint( sza / 10.0_wp ) )
    4394              IF (sza  .GE.  89.99_wp) THEN
     4393             IF (sza  .GE.  89.99_wp)  THEN
    43954394                sza = 89.99999_wp
    43964395             ENDIF
  • palm/trunk/SOURCE/time_integration.f90

    r3739 r3742  
    2525! -----------------
    2626! $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
    2731! Removed everything related to "time_bio_results" as this is never used.
    2832!
     
    16351639       IF ( time_dopr_listing >= dt_dopr_listing )  THEN
    16361640          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,    &
    16381642                                                           dt_3d ) )
    16391643       ENDIF
     
    16411645!
    16421646!--    Graphic output for PROFIL
    1643        IF (        time_dopr >= dt_dopr                                     &
     1647       IF (        time_dopr >= dt_dopr                                        &
    16441648            .AND.  time_since_reference_point >= skip_time_dopr )  THEN
    16451649          IF ( dopr_n /= 0 )  CALL data_output_profiles
     
    16581662!--    Output of spectra (formatted for use with PROFIL), in case of no
    16591663!--    time averaging, spectra has to be calculated before
    1660        IF (         time_dosp >= dt_dosp                                     &
     1664       IF (         time_dosp >= dt_dosp                                       &
    16611665             .AND.  time_since_reference_point >= skip_time_dosp )  THEN
    16621666          IF ( average_count_sp == 0 )  CALL calc_spectra
     
    16671671!
    16681672!--    2d-data output (cross-sections)
    1669        IF (         time_do2d_xy >= dt_do2d_xy                               &
     1673       IF (         time_do2d_xy >= dt_do2d_xy                                 &
    16701674             .AND.  time_since_reference_point >= skip_time_do2d_xy )  THEN
    16711675          CALL data_output_2d( 'xy', 0 )
    16721676          time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
    16731677       ENDIF
    1674        IF (         time_do2d_xz >= dt_do2d_xz                               &
     1678       IF (         time_do2d_xz >= dt_do2d_xz                                 &
    16751679             .AND.  time_since_reference_point >= skip_time_do2d_xz )  THEN
    16761680          CALL data_output_2d( 'xz', 0 )
    16771681          time_do2d_xz = MOD( time_do2d_xz, MAX( dt_do2d_xz, dt_3d ) )
    16781682       ENDIF
    1679        IF (         time_do2d_yz >= dt_do2d_yz                               &
     1683       IF (         time_do2d_yz >= dt_do2d_yz                                 &
    16801684             .AND.  time_since_reference_point >= skip_time_do2d_yz )  THEN
    16811685          CALL data_output_2d( 'yz', 0 )
     
    16851689!
    16861690!--    3d-data output (volume data)
    1687        IF (         time_do3d >= dt_do3d                                     &
     1691       IF (         time_do3d >= dt_do3d                                       &
    16881692             .AND.  time_since_reference_point >= skip_time_do3d )  THEN
    16891693          CALL data_output_3d( 0 )
     
    16941698!--    Masked data output
    16951699       DO  mid = 1, masks
    1696           IF (         time_domask(mid) >= dt_domask(mid)                          &
     1700          IF (         time_domask(mid) >= dt_domask(mid)                      &
    16971701                .AND.  time_since_reference_point >= skip_time_domask(mid) )  THEN
    16981702             CALL data_output_mask( 0 )
    1699              time_domask(mid) = MOD( time_domask(mid),                             &
     1703             time_domask(mid) = MOD( time_domask(mid),                         &
    17001704                                     MAX( dt_domask(mid), dt_3d ) )
    17011705          ENDIF
     
    17041708!
    17051709!--    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                            &
    17071711             .AND.  time_since_reference_point >= skip_time_data_output_av )  THEN
    17081712          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
    17091718          CALL data_output_2d( 'xy', 1 )
    17101719          CALL data_output_2d( 'xz', 1 )
     
    17191728!--    Output of surface data, instantaneous and averaged data
    17201729       IF ( surface_output )  THEN
    1721           IF (         time_dosurf >= dt_dosurf                                    &
     1730          IF (         time_dosurf >= dt_dosurf                                &
    17221731                .AND.  time_since_reference_point >= skip_time_dosurf )  THEN
    17231732             CALL surface_data_output( 0 )
    17241733             time_dosurf = MOD( time_dosurf, MAX( dt_dosurf, dt_3d ) )
    17251734          ENDIF
    1726           IF (         time_dosurf_av >= dt_dosurf_av                              &
     1735          IF (         time_dosurf_av >= dt_dosurf_av                          &
    17271736                .AND.  time_since_reference_point >= skip_time_dosurf_av )  THEN
    17281737             CALL surface_data_output( 1 )
Note: See TracChangeset for help on using the changeset viewer.