Ignore:
Timestamp:
Feb 13, 2019 12:35:12 PM (2 years ago)
Author:
dom_dwd_user
Message:

biometeorology_mod.f90:
(B) Added safety-meassure to catch the case that 'bio_mrt_av' is stated after 'bio_<index>' in the output section of the p3d file. This could lead to double-averaging and thus to silly results otherwise.

File:
1 edited

Legend:

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

    r3739 r3740  
    2727! -----------------
    2828! $Id$
     29! - Added safety-meassure to catch the case that 'bio_mrt_av' is stated after
     30! 'bio_<index>' in the output section of the p3d file.
     31!
     32! 3739 2019-02-13 08:05:17Z dom_dwd_user
    2933! - Auto-adjusting thermal_comfort flag if not set by user, but thermal_indices
    3034! set as output quantities.
     
    395399
    396400          CASE ( 'bio_mrt' )
    397                 IF ( .NOT. ALLOCATED( mrt_av_grid ) )  THEN
     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
    398409                   ALLOCATE( mrt_av_grid(nmrtbl) )
    399410                ENDIF
    400411                mrt_av_grid = 0.0_wp
     412
    401413
    402414          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
     
    410422!--          Indices are in unknown order as depending on the input file,
    411423!--          determine first index to average und update only once
    412 
     424!
    413425!--          Only proceed here if this was not done for any index before. This
    414426!--          is done only once during the whole model run.
     
    430442!
    431443!--          Only continue if var is the index, that controls averaging.
    432 !            Break immediatelly (doing nothing) for the other indices.
     444!--          Break immediatelly (doing nothing) for the other indices.
    433445             IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') RETURN
    434446             IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*')   RETURN
     
    436448!
    437449!--          Now memorize which of the input grids are not averaged by other
    438 !            modules. Set averaging switch to .TRUE. in that case.
     450!--          modules. Set averaging switch to .TRUE. in that case.
    439451             IF ( .NOT. ALLOCATED( pt_av ) )  THEN
    440452                ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     
    489501
    490502          CASE ( 'bio_mrt' )
    491              IF ( ALLOCATED( mrt_av_grid ) )  THEN
     503!
     504!--          Consider the case 'bio_mrt' is called after some thermal index. In
     505!--          that case do_average_mrt will be .TRUE. leading to a double-
     506!--          averaging.
     507             IF ( .NOT. do_average_mrt .AND. ALLOCATED( mrt_av_grid ) )  THEN
     508
     509                IF ( mrt_include_sw )  THEN
     510                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
     511                      ( ( human_absorb * mrtinsw(:) +                          &
     512                      mrtinlw(:) ) /                                           &
     513                      ( human_emiss * sigma_sb ) ) ** .25_wp - degc_to_k
     514                ELSE
     515                   mrt_av_grid(:) = mrt_av_grid(:) +                           &
     516                      ( mrtinlw(:) /                                           &
     517                      ( human_emiss * sigma_sb ) ) ** .25_wp - degc_to_k
     518                ENDIF
     519             ENDIF
     520
     521          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
     522!
     523!--          Only continue if the current index is the one to trigger the input
     524!--          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
     533                DO  i = nxl, nxr
     534                   DO  j = nys, nyn
     535                      DO  k = nzb, nzt+1
     536                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
     537                      ENDDO
     538                   ENDDO
     539                ENDDO
     540             ENDIF
     541
     542             IF ( ALLOCATED( q_av )  .AND. do_average_q ) THEN
     543                DO  i = nxl, nxr
     544                   DO  j = nys, nyn
     545                      DO  k = nzb, nzt+1
     546                         q_av(k,j,i) = q_av(k,j,i) + q(k,j,i)
     547                      ENDDO
     548                   ENDDO
     549                ENDDO
     550             ENDIF
     551
     552             IF ( ALLOCATED( u_av )  .AND. do_average_u ) THEN
     553                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
     554                   DO  j = nysg, nyng
     555                      DO  k = nzb, nzt+1
     556                         u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
     557                      ENDDO
     558                   ENDDO
     559                ENDDO
     560             ENDIF
     561
     562             IF ( ALLOCATED( v_av )  .AND. do_average_v ) THEN
     563                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
     564                   DO  j = nysg, nyng
     565                      DO  k = nzb, nzt+1
     566                         v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
     567                      ENDDO
     568                   ENDDO
     569                ENDDO
     570             ENDIF
     571
     572             IF ( ALLOCATED( w_av )  .AND. do_average_w ) THEN
     573                DO  i = nxlg, nxrg       !< yes, ghost points are required here!
     574                   DO  j = nysg, nyng
     575                      DO  k = nzb, nzt+1
     576                         w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
     577                      ENDDO
     578                   ENDDO
     579                ENDDO
     580             ENDIF
     581
     582             IF ( ALLOCATED( mrt_av_grid ) .AND. do_average_mrt )  THEN
    492583
    493584                IF ( mrt_include_sw )  THEN
     
    502593                ENDIF
    503594             ENDIF
    504 
    505           CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
    506 !
    507 !--          Only continue if updateindex, see above
    508              IF ( average_trigger_perct .AND. TRIM( variable ) /= 'bio_perct*') &
    509                 RETURN
    510              IF ( average_trigger_utci .AND. TRIM( variable ) /= 'bio_utci*')  &
    511                 RETURN
    512              IF ( average_trigger_pet  .AND. TRIM( variable ) /= 'bio_pet*')   &
    513                 RETURN
    514 
    515              IF ( ALLOCATED( pt_av ) .AND. do_average_theta ) THEN
    516                 DO  i = nxl, nxr
    517                    DO  j = nys, nyn
    518                       DO  k = nzb, nzt+1
    519                          pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
    520                       ENDDO
    521                    ENDDO
    522                 ENDDO
    523              ENDIF
    524 
    525              IF ( ALLOCATED( q_av )  .AND. do_average_q ) THEN
    526                 DO  i = nxl, nxr
    527                    DO  j = nys, nyn
    528                       DO  k = nzb, nzt+1
    529                          q_av(k,j,i) = q_av(k,j,i) + q(k,j,i)
    530                       ENDDO
    531                    ENDDO
    532                 ENDDO
    533              ENDIF
    534 
    535              IF ( ALLOCATED( u_av )  .AND. do_average_u ) THEN
    536                 DO  i = nxlg, nxrg       !< yes, ghost points are required here!
    537                    DO  j = nysg, nyng
    538                       DO  k = nzb, nzt+1
    539                          u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
    540                       ENDDO
    541                    ENDDO
    542                 ENDDO
    543              ENDIF
    544 
    545              IF ( ALLOCATED( v_av )  .AND. do_average_v ) THEN
    546                 DO  i = nxlg, nxrg       !< yes, ghost points are required here!
    547                    DO  j = nysg, nyng
    548                       DO  k = nzb, nzt+1
    549                          v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
    550                       ENDDO
    551                    ENDDO
    552                 ENDDO
    553              ENDIF
    554 
    555              IF ( ALLOCATED( w_av )  .AND. do_average_w ) THEN
    556                 DO  i = nxlg, nxrg       !< yes, ghost points are required here!
    557                    DO  j = nysg, nyng
    558                       DO  k = nzb, nzt+1
    559                          w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
    560                       ENDDO
    561                    ENDDO
    562                 ENDDO
    563              ENDIF
    564 
    565              IF ( ALLOCATED( mrt_av_grid ) .AND. do_average_mrt )  THEN
    566 
    567                 IF ( mrt_include_sw )  THEN
    568                    mrt_av_grid(:) = mrt_av_grid(:) +                           &
    569                       ( ( human_absorb * mrtinsw(:) +                          &
    570                       mrtinlw(:) ) /                             &  ! human_emiss * mrtinlw(:) /
    571                       ( human_emiss * sigma_sb ) ) ** .25_wp - degc_to_k
    572                 ELSE
    573                    mrt_av_grid(:) = mrt_av_grid(:) +                           &
    574                       ( mrtinlw(:) /                             &  ! ( human_emiss * mrtinlw(:) /
    575                       ( human_emiss * sigma_sb ) ) ** .25_wp - degc_to_k
    576                 ENDIF
    577              ENDIF
    578595!
    579596!--       This is a cumulated dose. No mode == 'average' for this quantity.
     
    597614
    598615          CASE ( 'bio_mrt' )
    599              IF ( ALLOCATED( mrt_av_grid ) )  THEN
     616!
     617!--          Consider the case 'bio_mrt' is called after some thermal index. In
     618!--          that case do_average_mrt will be .TRUE. leading to a double-
     619!--          averaging.
     620             IF ( .NOT. do_average_mrt .AND. ALLOCATED( mrt_av_grid ) )  THEN
    600621                mrt_av_grid(:) = mrt_av_grid(:) / REAL( average_count_3d, KIND=wp )
    601622             ENDIF
     
    815836    END SELECT
    816837
     838!
     839!--    Further checks if thermal comfort output is desired.
    817840    IF ( thermal_comfort .AND. unit == 'degree_C' )  THEN
    818 !
    819 !--    Futher checks if output belongs to biometeorology
    820 !      Break if required modules "radiation" and "humidity" are not running.
     841
     842!
     843!--    Break if required modules "radiation" and "humidity" are not avalable.
    821844       IF ( .NOT.  radiation ) THEN
    822845          message_string = 'output of "' // TRIM( var ) // '" require'         &
     
    848871! ------------
    849872!> Check parameters routine for biom module
    850 !> check_parameters 1302
     873!> Currently unused but might come in handy for future checks?
    851874!------------------------------------------------------------------------------!
    852875 SUBROUTINE bio_check_parameters
     
    854877    USE control_parameters,                                                    &
    855878        ONLY:  message_string
    856 
    857879
    858880    IMPLICIT NONE
     
    917939                 IF ( mrt_include_sw )  THEN
    918940                    local_pf(i,j,k) = ( ( human_absorb * mrtinsw(l) +          &
    919                                     mrtinlw(l) ) /               &  ! human_emiss * mrtinlw(l) ) /
     941                                    mrtinlw(l) ) /                             &
    920942                                    ( human_emiss * sigma_sb ) ) ** .25_wp -   &
    921943                                    degc_to_k
    922944                 ELSE
    923                     local_pf(i,j,k) = ( mrtinlw(l)  /          &  ! ( (human_emiss * mrtinlw(l) ) /
     945                    local_pf(i,j,k) = ( mrtinlw(l)  /                          &
    924946                                    ( human_emiss * sigma_sb ) ) ** .25_wp -   &
    925947                                    degc_to_k
     
    10701092                  IF ( mrt_include_sw )  THEN
    10711093                     local_pf(i,j,k) = REAL( ( ( human_absorb * mrtinsw(l) +   &
    1072                                     mrtinlw(l) ) /               &  ! human_emiss * mrtinlw(l) ) /
     1094                                    mrtinlw(l) ) /                             &
    10731095                                    ( human_emiss * sigma_sb ) ) ** .25_wp -   &
    10741096                                    degc_to_k, KIND = sp )
    10751097                  ELSE
    1076                      local_pf(i,j,k) = REAL( ( mrtinlw(l) /  &  ! REAL( ( ( human_emiss * mrtinlw(l) ) /
     1098                     local_pf(i,j,k) = REAL( ( mrtinlw(l) /                    &
    10771099                                    ( human_emiss * sigma_sb ) ) ** .25_wp -   &
    10781100                                    degc_to_k, KIND = sp )
     
    14611483
    14621484
     1485!
     1486!-- We need to differentiate if averaged input is desired (av == .TRUE.) or not.
    14631487    IF ( av ) THEN
     1488!
     1489!-- Make sure tmrt_av_grid is present and initialize with the fill value
    14641490       IF ( .NOT. ALLOCATED( tmrt_av_grid ) )  THEN
    14651491          ALLOCATE( tmrt_av_grid (nys:nyn,nxl:nxr) )
     
    14671493       tmrt_av_grid = REAL( bio_fill_value, KIND = wp )
    14681494
    1469        DO  l = 1, nmrtbl
    1470           i = mrtbl(ix,l)
    1471           j = mrtbl(iy,l)
    1472           k = mrtbl(iz,l)
    1473           IF ( k - get_topography_top_index_ji( j, i, 's' ) ==                 &
    1474                 bio_cell_level + 1_iwp) THEN
    1475 
    1476              tmrt_av_grid(j,i) = mrt_av_grid(l)
    1477 
    1478           ENDIF
    1479        ENDDO
    1480 
     1495!
     1496!-- mrt_av_grid should always be allcoated here, but better make sure ist actually is.
     1497       IF ( ALLOCATED( mrt_av_grid ) )  THEN
     1498!
     1499!-- Iterate over the radiation grid (radiation coordinates) and fill the
     1500!-- tmrt_av_grid (x, y coordinates) where appropriate:
     1501!-- tmrt_av_grid is written for all i / j if level (k) matches output height.
     1502          DO  l = 1, nmrtbl
     1503             i = mrtbl(ix,l)
     1504             j = mrtbl(iy,l)
     1505             k = mrtbl(iz,l)
     1506             IF ( k - get_topography_top_index_ji( j, i, 's' ) ==              &
     1507                   bio_cell_level + 1_iwp) THEN
     1508!
     1509!-- Averaging was done before, so we can just copy the result here
     1510                tmrt_av_grid(j,i) = mrt_av_grid(l)
     1511
     1512             ENDIF
     1513          ENDDO
     1514       ENDIF
     1515
     1516!
     1517!-- In case instantaneous input is desired, mrt values will be re-calculated.
    14811518    ELSE
    14821519!
     
    14971534             IF ( mrt_include_sw )  THEN
    14981535                 tmrt_grid(j,i) = ( ( human_absorb * mrtinsw(l) +              &
    1499                                   mrtinlw(l) )  /                &  ! human_emiss * mrtinlw(l) )  /
     1536                                  mrtinlw(l) )  /                              &
    15001537                                  ( human_emiss * sigma_sb ) ) ** .25_wp -     &
    15011538                                  degc_to_k
    15021539             ELSE
    1503                  tmrt_grid(j,i) = ( mrtinlw(l)  /            &  ! ( ( human_emiss * mrtinlw(l) )  /
     1540                 tmrt_grid(j,i) = ( mrtinlw(l)  /                              &
    15041541                                  ( human_emiss * sigma_sb ) ) ** .25_wp -     &
    15051542                                  degc_to_k
     
    15361573!-- Internal variables
    15371574    INTEGER(iwp)                ::  k     !< Running index, z-dir
    1538 !     INTEGER(iwp)                ::  ir    !< Running index, x-dir, radiation coordinates
    1539 !     INTEGER(iwp)                ::  jr    !< Running index, y-dir, radiation coordinates
    1540 !     INTEGER(iwp)                ::  kr    !< Running index, y-dir, radiation coordinates
    15411575    INTEGER(iwp)                ::  k_wind  !< Running index, z-dir, wind speed only
    1542 !     INTEGER(iwp)                ::  l     !< Running index, radiation coordinates
    15431576
    15441577    REAL(wp)                    ::  vp_sat  !< Saturation vapor pressure     (hPa)
Note: See TracChangeset for help on using the changeset viewer.