Ignore:
Timestamp:
Jan 31, 2020 4:14:57 PM (20 months ago)
Author:
pavelkrc
Message:

Merge branch resler (updated documentation, optional flux tracing, bugfix)

Location:
palm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/radiation_model_mod.f90

    r4360 r4392  
    1515! PALM. If not, see <http://www.gnu.org/licenses/>.
    1616!
    17 ! Copyright 2015-2019 Institute of Computer Science of the
     17! Copyright 2015-2020 Institute of Computer Science of the
    1818!                     Czech Academy of Sciences, Prague
    1919! Copyright 2015-2019 Czech Technical University in Prague
     
    2828! -----------------
    2929! $Id$
     30! - Add debug tracing of large radiative fluxes (option trace_fluxes_above)
     31! - Print exact counts of SVF and CSF if debut_output is enabled
     32! - Update descriptions of RTM 3.0 and related comments
     33!
     34! 4360 2020-01-07 11:25:50Z suehring
    3035! Renamed pc_heating_rate, pc_transpiration_rate, pc_transpiration_rate to
    3136! pcm_heating_rate, pcm_latent_rate, pcm_transpiration_rate
     
    214219! Description:
    215220! ------------
    216 !> Radiation models and interfaces
    217 !> @todo Replace dz(1) appropriatly to account for grid stretching
     221!> Radiation models and interfaces:
     222!> constant, simple and RRTMG models, interface to external radiation model
     223!> Radiative Transfer Model (RTM) version 3.0 for modelling of radiation
     224!> interactions within urban canopy or other surface layer in complex terrain
     225!> Integrations of RTM with other PALM-4U modules:
     226!> integration with RRTMG, USM, LSM, PCM, BIO modules
     227!>
    218228!> @todo move variable definitions used in radiation_init only to the subroutine
    219229!>       as they are no longer required after initialization.
    220230!> @todo Output of full column vertical profiles used in RRTMG
    221231!> @todo Output of other rrtm arrays (such as volume mixing ratios)
    222 !> @todo Check for mis-used NINT() calls in raytrace_2d
    223 !>       RESULT: Original was correct (carefully verified formula), the change
    224 !>               to INT broke raytracing      -- P. Krc
    225232!> @todo Optimize radiation_tendency routines
    226 !> @todo Consider rotated model domains (rotation_angle/=0.0)
    227233!>
    228234!> @note Many variables have a leading dummy dimension (0:0) in order to
     
    430436                skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
    431437                sky_trans,                       & !< sky transmissivity
    432                 time_radiation = 0.0_wp            !< time since last call of radiation code
     438                time_radiation = 0.0_wp,         & !< time since last call of radiation code
     439                trace_fluxes_above = -1.0_wp       !< NAMELIST option for debug tracing of large radiative fluxes (W/m2;W/m3)
    433440
    434441    INTEGER(iwp) ::  day_of_year   !< day of the current year
     
    15321539! Description:
    15331540! ------------
    1534 !> Initialization of the radiation model
     1541!> Initialization of the radiation model and Radiative Transfer Model
    15351542!------------------------------------------------------------------------------!
    15361543    SUBROUTINE radiation_init
     
    38483855! Description:
    38493856! ------------
    3850 !> Parin for &radiation_parameters for radiation model
     3857!> Parin for &radiation_parameters for radiation model and RTM
    38513858!------------------------------------------------------------------------------!
    38523859    SUBROUTINE radiation_parin
     
    38683875                                  raytrace_discrete_azims,                      &
    38693876                                  raytrace_discrete_elevs, raytrace_mpi_rma,    &
     3877                                  trace_fluxes_above,                           &
    38703878                                  skip_time_do_radiation, surface_reflections,  &
    38713879                                  svfnorm_report_thresh, sw_radiation,          &
     
    38843892                                  raytrace_discrete_azims,                      &
    38853893                                  raytrace_discrete_elevs, raytrace_mpi_rma,    &
     3894                                  trace_fluxes_above,                           &
    38863895                                  skip_time_do_radiation, surface_reflections,  &
    38873896                                  svfnorm_report_thresh, sw_radiation,          &
     
    56985707! Description:
    56995708! ------------
    5700 !> This subroutine calculates interaction of the solar radiation
     5709!> Radiative Transfer Model (RTM) version 3.0 for modelling of radiation
     5710!> interactions within urban canopy or inside of surface layer in complex terrain.
     5711!> This subroutine calculates interaction of the solar SW and LW radiation
    57015712!> with urban and land surfaces and updates all surface heatfluxes.
    5702 !> It calculates also the required parameters for RRTMG lower BC.
     5713!> It also calculates interactions of SW and LW radiation with resolved
     5714!> plant canopy and calculates the corresponding plant canopy heat fluxes.
     5715!> The subroutine also models spatial and temporal distribution of Mean
     5716!> Radiant Temperature (MRT). The resulting values are provided to other
     5717!> PALM-4U modules (RRTMG, USM, LSM, PCM and BIO).
    57035718!>
    5704 !> For more info. see Resler et al. 2017
     5719!> The new version 3.0 was radically rewriten from version 1.0.
     5720!> The most significant changes include new angular discretization scheme,
     5721!> redesigned and significantly optimized raytracing scheme, new processes
     5722!> included in modelling (e.g. intetrations of LW radiation with PC),
     5723!> integrated calculation of Mean Radiant Temperature (MRT), and improved
     5724!> and enhanced output and debug capabilities. This new version significantly
     5725!> improves effectivity of the paralelization and the scalability of the model
     5726!> and allows simulation of extensive domain with appropriate HPC resources.
    57055727!>
    5706 !> The new version 2.0 was radically rewriten, the discretization scheme
    5707 !> has been changed. This new version significantly improves effectivity
    5708 !> of the paralelization and the scalability of the model.
     5728!> More info about RTM v.1.0. see:
     5729!> Resler et al., GMD. 2017, https://doi.org/10.5194/gmd-10-3635-2017
     5730!> Info about RTM v. 3.0 see:
     5731!> Krc et al. 2020 (to appear in GMD),
     5732!> Maronga et al., GMDD 2019,  https://doi.org/10.5194/gmd-2019-103
     5733!>
     5734
     5735
    57095736!------------------------------------------------------------------------------!
    57105737
     
    59005927     ENDDO
    59015928
     5929     IF ( trace_fluxes_above >= 0._wp )  THEN
     5930        CALL radiation_print_debug_surf( 'surfoutll before initial pass', surfoutll )
     5931        CALL radiation_print_debug_horz( 'rad_lw_in_diff before initial pass', rad_lw_in_diff )
     5932        CALL radiation_print_debug_horz( 'rad_sw_in_diff before initial pass', rad_sw_in_diff )
     5933        CALL radiation_print_debug_horz( 'rad_sw_in_dir before initial pass', rad_sw_in_dir )
     5934     ENDIF
     5935
    59025936#if defined( __parallel )
    59035937!--     might be optimized and gather only values relevant for current processor
     
    60396073     ENDIF
    60406074
     6075     IF ( trace_fluxes_above >= 0._wp )  THEN
     6076        CALL radiation_print_debug_surf( 'surfinl after initial pass', surfinl )
     6077        CALL radiation_print_debug_surf( 'surfinlwdif after initial pass', surfinlwdif )
     6078        CALL radiation_print_debug_surf( 'surfinswdif after initial pass', surfinswdif )
     6079        CALL radiation_print_debug_surf( 'surfinswdir after initial pass', surfinswdir )
     6080        IF ( npcbl > 0 )  THEN
     6081           CALL radiation_print_debug_pcb( 'pcbinlw after initial pass', pcbinlw )
     6082           CALL radiation_print_debug_pcb( 'pcbinswdif after initial pass', pcbinswdif )
     6083           CALL radiation_print_debug_pcb( 'pcbinswdir after initial pass', pcbinswdir )
     6084        ENDIF
     6085     ENDIF
     6086
    60416087     IF ( plant_lw_interact )  THEN
    60426088!
     
    60546100     ENDIF
    60556101
     6102     IF ( trace_fluxes_above >= 0._wp )  THEN
     6103        CALL radiation_print_debug_surf( 'surfinl after PC emiss', surfinl )
     6104     ENDIF
     6105
    60566106     surfins = surfinswdir + surfinswdif
    60576107     surfinl = surfinl + surfinlwdif
     
    60816131!--      for non-transparent surfaces, longwave albedo is 1 - emissivity
    60826132         surfoutll = (1._wp - emiss_surf) * surfinl
     6133
     6134         IF ( trace_fluxes_above >= 0._wp )  THEN
     6135            CALL radiation_print_debug_surf( 'surfoutll before reflective pass', surfoutll, refstep )
     6136            CALL radiation_print_debug_surf( 'surfoutsl before reflective pass', surfoutsl, refstep )
     6137         ENDIF
    60836138
    60846139#if defined( __parallel )
     
    61486203            mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc)
    61496204         ENDDO
     6205
     6206         IF ( trace_fluxes_above >= 0._wp )  THEN
     6207            CALL radiation_print_debug_surf( 'surfinl after reflected pass', surfinl, refstep )
     6208            CALL radiation_print_debug_surf( 'surfins after reflected pass', surfins, refstep )
     6209            CALL radiation_print_debug_pcb( 'pcbinlw after reflected pass', pcbinlw, refstep )
     6210            CALL radiation_print_debug_pcb( 'pcbinsw after reflected pass', pcbinsw, refstep )
     6211         ENDIF
    61506212
    61516213         surfinsw = surfinsw  + surfins
     
    65826644! Description:
    65836645! ------------
    6584 !> This subroutine splits direct and diffusion dw radiation
     6646!> This subroutine splits direct and diffusion dw radiation for RTM processing.
    65856647!> It sould not be called in case the radiation model already does it
    65866648!> It follows Boland, Ridley & Brown (2008)
     
    66426704    END SUBROUTINE calc_diffusion_radiation
    66436705
     6706!------------------------------------------------------------------------------!
     6707! Description:
     6708! ------------
     6709!> Print consecutive radiative extremes if requested to trace early radiation
     6710!> interaction instabilities.
     6711!------------------------------------------------------------------------------!
     6712    SUBROUTINE radiation_print_debug_surf( description, values, step )
     6713
     6714       CHARACTER (LEN=*), INTENT(in)      ::  description
     6715       REAL(wp), DIMENSION(:), INTENT(in) ::  values
     6716       INTEGER(iwp), INTENT(in), OPTIONAL ::  step
     6717
     6718       CHARACTER (LEN=50)                 ::  location
     6719       CHARACTER (LEN=1024)               ::  debug_string
     6720       INTEGER                            ::  isurf
     6721       REAL(wp)                           ::  x
     6722
     6723       isurf = MAXLOC( values, DIM=1 )
     6724       x = values(isurf)
     6725       IF ( x < trace_fluxes_above )  RETURN
     6726
     6727       IF ( PRESENT( step ) )  THEN
     6728          WRITE( location, '(A," #",I0)' ) description, step
     6729       ELSE
     6730          location = description
     6731       ENDIF
     6732
     6733       WRITE( debug_string, '("Maximum of ",A50," = ",F12.1," at coords ' //  &
     6734                            'i=",I4,", j=",I4,", k=",I4,", d=",I1,". '    //  &
     6735                            'Alb=",F7.3,", emis=",F7.3)'                    ) &
     6736              location, x, surfl(ix,isurf), surfl(iy,isurf),                  &
     6737              surfl(iz,isurf), surfl(id,isurf), albedo_surf(isurf),           &
     6738              emiss_surf(isurf)
     6739       CALL debug_message( debug_string, 'info' )
     6740
     6741    END SUBROUTINE
     6742
     6743    SUBROUTINE radiation_print_debug_pcb( description, values, step )
     6744
     6745       CHARACTER (LEN=*), INTENT(in)      ::  description
     6746       REAL(wp), DIMENSION(:), INTENT(in) ::  values
     6747       INTEGER(iwp), INTENT(in), OPTIONAL ::  step
     6748
     6749       CHARACTER (LEN=50)                 ::  location
     6750       CHARACTER (LEN=1024)               ::  debug_string
     6751       INTEGER                            ::  ipcb
     6752       REAL(wp)                           ::  x
     6753
     6754       IF ( npcbl <= 0 )  RETURN
     6755       ipcb = MAXLOC( values, DIM=1 )
     6756       x = values(ipcb) / (dx*dy*dz(1))
     6757       IF ( x < trace_fluxes_above )  RETURN
     6758
     6759       IF ( PRESENT( step ) )  THEN
     6760          WRITE( location, '(A," #",I0)' ) description, step
     6761       ELSE
     6762          location = description
     6763       ENDIF
     6764
     6765       WRITE( debug_string, '("Maximum of ",A50," = ",F12.1," at coords ' //  &
     6766                            'i=",I4,", j=",I4,", k=",I4)'                   ) &
     6767              location, x, pcbl(ix,ipcb), pcbl(iy,ipcb), pcbl(iz,ipcb)
     6768       CALL debug_message( debug_string, 'info' )
     6769
     6770    END SUBROUTINE
     6771
     6772    SUBROUTINE radiation_print_debug_horz( description, values, step )
     6773
     6774       CHARACTER (LEN=*), INTENT(in)        ::  description
     6775       REAL(wp), DIMENSION(:,:), INTENT(in) ::  values
     6776       INTEGER(iwp), INTENT(in), OPTIONAL   ::  step
     6777
     6778       CHARACTER (LEN=50)                   ::  location
     6779       CHARACTER (LEN=1024)                 ::  debug_string
     6780       INTEGER, DIMENSION(2)                ::  ji
     6781       REAL(wp)                             ::  x
     6782
     6783       ji = MAXLOC( values )
     6784       x = values(ji(1),ji(2))
     6785       IF ( x < trace_fluxes_above )  RETURN
     6786
     6787       IF ( PRESENT( step ) )  THEN
     6788          WRITE( location, '(A," #",I0)' ) description, step
     6789       ELSE
     6790          location = description
     6791       ENDIF
     6792
     6793       WRITE( debug_string, '("Maximum of ",A50," = ",F12.1," at coords ' //  &
     6794                            'i=",I4,", j=",I4)'                             ) &
     6795              location, x, ji(2), ji(1)
     6796       CALL debug_message( debug_string, 'info' )
     6797
     6798    END SUBROUTINE
     6799
    66446800 END SUBROUTINE radiation_interaction
    66456801   
     
    66476803! Description:
    66486804! ------------
    6649 !> This subroutine initializes structures needed for radiative transfer
    6650 !> model. This model calculates transformation processes of the
     6805!> This subroutine initializes structures needed for Radiative Transfer
     6806!> Model (RTM). This model calculates transformation processes of the
    66516807!> radiation inside urban and land canopy layer. The module includes also
    66526808!> the interaction of the radiation with the resolved plant canopy.
    6653 !>
    6654 !> For more info. see Resler et al. 2017
    6655 !>
    6656 !> The new version 2.0 was radically rewriten, the discretization scheme
    6657 !> has been changed. This new version significantly improves effectivity
    6658 !> of the paralelization and the scalability of the model.
    66596809!>
    66606810!------------------------------------------------------------------------------!
     
    71827332!> Calculates shape view factors (SVF), plant sink canopy factors (PCSF),
    71837333!> sky-view factors, discretized path for direct solar radiation, MRT factors
    7184 !> and other preprocessed data needed for radiation_interaction.
     7334!> and other preprocessed data needed for radiation_interaction inside RTM.
     7335!> This subroutine is called only one at the beginning of the simulation.
     7336!> The resulting factors can be stored to files and reused with other
     7337!> simulations utilizing the same surface and plant canopy structure.
    71857338!------------------------------------------------------------------------------!
    71867339    SUBROUTINE radiation_calc_svf
     
    79448097
    79458098        IF ( rad_angular_discretization )  THEN
    7946            IF ( debug_output )  CALL debug_message( 'Load svf from the structure array to plain arrays', 'info' )
     8099           IF ( debug_output )  THEN
     8100              WRITE( debug_string, '("Load ",I0," SVFs from the structure array to plain arrays")' ) nsvfl
     8101              CALL debug_message( debug_string, 'info' )
     8102           ENDIF
    79478103           ALLOCATE( svf(ndsvf,nsvfl) )
    79488104           ALLOCATE( svfsurf(idsvf,nsvfl) )
     
    79588114
    79598115           !< load svf from the structure array to plain arrays
    7960            IF ( debug_output )  CALL debug_message( 'Load svf from the structure array to plain arrays', 'info' )
     8116           IF ( debug_output )  THEN
     8117              WRITE( debug_string, '("Load ",I0," SVFs from the structure array to plain arrays")' ) nsvfl
     8118              CALL debug_message( debug_string, 'info' )
     8119           ENDIF
    79618120           ALLOCATE( svf(ndsvf,nsvfl) )
    79628121           ALLOCATE( svfsurf(idsvf,nsvfl) )
     
    81898348            DEALLOCATE( pcsflt_l )
    81908349            DEALLOCATE( kpcsflt_l )
    8191             IF ( debug_output )  CALL debug_message( 'End of aggregate csf', 'info' )
     8350            IF ( debug_output )  THEN
     8351               WRITE( debug_string, '("Finished aggregating ",I0," CSFs.")') ncsfl
     8352               CALL debug_message( debug_string, 'info' )
     8353            ENDIF
    81928354           
    81938355        ENDIF
     
    82128374! ------------
    82138375!> Raytracing for detecting obstacles and calculating compound canopy sink
    8214 !> factors. (A simple obstacle detection would only need to process faces in
    8215 !> 3 dimensions without any ordering.)
     8376!> factors for RTM. (A simple obstacle detection would only need to process
     8377!> faces in 3 dimensions without any ordering.)
    82168378!> Assumtions:
    82178379!> -----------
     
    84108572! ------------
    84118573!> A new, more efficient version of ray tracing algorithm that processes a whole
    8412 !> arc instead of a single ray.
     8574!> arc instead of a single ray (new in RTM version 2.5).
    84138575!>
    84148576!> In all comments, horizon means tangent of horizon angle, i.e.
     
    89259087! ------------
    89269088!> Calculates apparent solar positions for all timesteps and stores discretized
    8927 !> positions.
     9089!> positions for RTM.
    89289090!------------------------------------------------------------------------------!
    89299091   SUBROUTINE radiation_presimulate_solar_pos
     
    90139175! Description:
    90149176! ------------
    9015 !> Determines whether two faces are oriented towards each other. Since the
     9177!> Determines whether two faces are oriented towards each other in RTM. Since the
    90169178!> surfaces follow the gird box surfaces, it checks first whether the two surfaces
    90179179!> are directed in the same direction, then it checks if the two surfaces are
     
    90749236! Description:
    90759237! ------------
    9076 !> Reads svf, svfsurf, csf, csfsurf and mrt factors data from saved file
    9077 !> SVF means sky view factors and CSF means canopy sink factors
     9238!> Reads svf, svfsurf, csf, csfsurf and mrt factors data from saved file.
     9239!> This allows to skip their calculation during of RTM init phase.
     9240!> SVF means sky view factors and CSF means canopy sink factors.
    90789241!------------------------------------------------------------------------------!
    90799242    SUBROUTINE radiation_read_svf
     
    92409403! ------------
    92419404!> Subroutine stores svf, svfsurf, csf, csfsurf and mrt data to a file.
     9405!> The stored factors can be reused in future simulation with the same
     9406!> geometry structure of the surfaces and resolved plant canopy.
    92429407!------------------------------------------------------------------------------!
    92439408    SUBROUTINE radiation_write_svf
     
    93049469! Description:
    93059470! ------------
    9306 !> Block of auxiliary subroutines:
     9471!> Block of auxiliary subroutines for RTM:
    93079472!> 1. quicksort and corresponding comparison
    93089473!> 2. merge_and_grow_csf for implementation of "dynamical growing"
     
    94409605! Description:
    94419606! ------------
    9442 !> Grows the CSF array exponentially after it is full. During that, the ray
    9443 !> canopy sink factors with common source face and target plant canopy grid
    9444 !> cell are merged together so that the size doesn't grow out of control.
     9607!> Grows the CSF array in RTM exponentially when it is full. During that,
     9608!> the ray canopy sink factors with common source face and target plant canopy
     9609!> grid cell are merged together so that the size doesn't grow out of control.
    94459610!------------------------------------------------------------------------------!
    94469611    SUBROUTINE merge_and_grow_csf(newsize)
Note: See TracChangeset for help on using the changeset viewer.