Changeset 4392 for palm/trunk/SOURCE/radiation_model_mod.f90

Timestamp:

Jan 31, 2020 4:14:57 PM (5 years ago)

Author:

pavelkrc

Message:

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

Location:

palm/trunk

Files:

• . (modified) (1 prop)
• SOURCE (modified) (1 prop)
• SOURCE/radiation_model_mod.f90 (modified) (29 diffs, 1 prop)

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -15 +15 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 2015-2019 Institute of Computer Science of the 
+! Copyright 2015-2020 Institute of Computer Science of the
 !                     Czech Academy of Sciences, Prague
 ! Copyright 2015-2019 Czech Technical University in Prague
@@ -28 +28 @@
 ! -----------------
 ! $Id$
+! - Add debug tracing of large radiative fluxes (option trace_fluxes_above)
+! - Print exact counts of SVF and CSF if debut_output is enabled
+! - Update descriptions of RTM 3.0 and related comments
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Renamed pc_heating_rate, pc_transpiration_rate, pc_transpiration_rate to
 ! pcm_heating_rate, pcm_latent_rate, pcm_transpiration_rate
@@ -214 +219 @@
 ! Description:
 ! ------------
-!> Radiation models and interfaces
-!> @todo Replace dz(1) appropriatly to account for grid stretching
+!> Radiation models and interfaces:
+!> constant, simple and RRTMG models, interface to external radiation model
+!> Radiative Transfer Model (RTM) version 3.0 for modelling of radiation
+!> interactions within urban canopy or other surface layer in complex terrain
+!> Integrations of RTM with other PALM-4U modules:
+!> integration with RRTMG, USM, LSM, PCM, BIO modules
+!>
 !> @todo move variable definitions used in radiation_init only to the subroutine
 !>       as they are no longer required after initialization.
 !> @todo Output of full column vertical profiles used in RRTMG
 !> @todo Output of other rrtm arrays (such as volume mixing ratios)
-!> @todo Check for mis-used NINT() calls in raytrace_2d
-!>       RESULT: Original was correct (carefully verified formula), the change
-!>               to INT broke raytracing      -- P. Krc
 !> @todo Optimize radiation_tendency routines
-!> @todo Consider rotated model domains (rotation_angle/=0.0)
 !>
 !> @note Many variables have a leading dummy dimension (0:0) in order to
@@ -430 +436 @@
                 skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
                 sky_trans,                       & !< sky transmissivity
-                time_radiation = 0.0_wp            !< time since last call of radiation code
+                time_radiation = 0.0_wp,         & !< time since last call of radiation code
+                trace_fluxes_above = -1.0_wp       !< NAMELIST option for debug tracing of large radiative fluxes (W/m2;W/m3)
 
     INTEGER(iwp) ::  day_of_year   !< day of the current year
@@ -1532 +1539 @@
 ! Description:
 ! ------------
-!> Initialization of the radiation model
+!> Initialization of the radiation model and Radiative Transfer Model
 !------------------------------------------------------------------------------!
     SUBROUTINE radiation_init
@@ -3848 +3855 @@
 ! Description:
 ! ------------
-!> Parin for &radiation_parameters for radiation model
+!> Parin for &radiation_parameters for radiation model and RTM
 !------------------------------------------------------------------------------!
     SUBROUTINE radiation_parin
@@ -3868 +3875 @@
                                   raytrace_discrete_azims,                      &
                                   raytrace_discrete_elevs, raytrace_mpi_rma,    &
+                                  trace_fluxes_above,                           &
                                   skip_time_do_radiation, surface_reflections,  &
                                   svfnorm_report_thresh, sw_radiation,          &
@@ -3884 +3892 @@
                                   raytrace_discrete_azims,                      &
                                   raytrace_discrete_elevs, raytrace_mpi_rma,    &
+                                  trace_fluxes_above,                           &
                                   skip_time_do_radiation, surface_reflections,  &
                                   svfnorm_report_thresh, sw_radiation,          &
@@ -5698 +5707 @@
 ! Description:
 ! ------------
-!> This subroutine calculates interaction of the solar radiation
+!> Radiative Transfer Model (RTM) version 3.0 for modelling of radiation
+!> interactions within urban canopy or inside of surface layer in complex terrain.
+!> This subroutine calculates interaction of the solar SW and LW radiation
 !> with urban and land surfaces and updates all surface heatfluxes.
-!> It calculates also the required parameters for RRTMG lower BC.
+!> It also calculates interactions of SW and LW radiation with resolved
+!> plant canopy and calculates the corresponding plant canopy heat fluxes.
+!> The subroutine also models spatial and temporal distribution of Mean
+!> Radiant Temperature (MRT). The resulting values are provided to other
+!> PALM-4U modules (RRTMG, USM, LSM, PCM and BIO).
 !>
-!> For more info. see Resler et al. 2017
+!> The new version 3.0 was radically rewriten from version 1.0.
+!> The most significant changes include new angular discretization scheme,
+!> redesigned and significantly optimized raytracing scheme, new processes
+!> included in modelling (e.g. intetrations of LW radiation with PC),
+!> integrated calculation of Mean Radiant Temperature (MRT), and improved
+!> and enhanced output and debug capabilities. This new version significantly
+!> improves effectivity of the paralelization and the scalability of the model
+!> and allows simulation of extensive domain with appropriate HPC resources.
 !>
-!> The new version 2.0 was radically rewriten, the discretization scheme
-!> has been changed. This new version significantly improves effectivity
-!> of the paralelization and the scalability of the model.
+!> More info about RTM v.1.0. see:
+!> Resler et al., GMD. 2017, https://doi.org/10.5194/gmd-10-3635-2017
+!> Info about RTM v. 3.0 see:
+!> Krc et al. 2020 (to appear in GMD),
+!> Maronga et al., GMDD 2019,  https://doi.org/10.5194/gmd-2019-103
+!>
+
+
 !------------------------------------------------------------------------------!
 
@@ -5900 +5927 @@
      ENDDO
 
+     IF ( trace_fluxes_above >= 0._wp )  THEN
+        CALL radiation_print_debug_surf( 'surfoutll before initial pass', surfoutll )
+        CALL radiation_print_debug_horz( 'rad_lw_in_diff before initial pass', rad_lw_in_diff )
+        CALL radiation_print_debug_horz( 'rad_sw_in_diff before initial pass', rad_sw_in_diff )
+        CALL radiation_print_debug_horz( 'rad_sw_in_dir before initial pass', rad_sw_in_dir )
+     ENDIF
+
 #if defined( __parallel )
 !--     might be optimized and gather only values relevant for current processor
@@ -6039 +6073 @@
      ENDIF
 
+     IF ( trace_fluxes_above >= 0._wp )  THEN
+        CALL radiation_print_debug_surf( 'surfinl after initial pass', surfinl )
+        CALL radiation_print_debug_surf( 'surfinlwdif after initial pass', surfinlwdif )
+        CALL radiation_print_debug_surf( 'surfinswdif after initial pass', surfinswdif )
+        CALL radiation_print_debug_surf( 'surfinswdir after initial pass', surfinswdir )
+        IF ( npcbl > 0 )  THEN
+           CALL radiation_print_debug_pcb( 'pcbinlw after initial pass', pcbinlw )
+           CALL radiation_print_debug_pcb( 'pcbinswdif after initial pass', pcbinswdif )
+           CALL radiation_print_debug_pcb( 'pcbinswdir after initial pass', pcbinswdir )
+        ENDIF
+     ENDIF
+
      IF ( plant_lw_interact )  THEN
 !
@@ -6054 +6100 @@
      ENDIF
 
+     IF ( trace_fluxes_above >= 0._wp )  THEN
+        CALL radiation_print_debug_surf( 'surfinl after PC emiss', surfinl )
+     ENDIF
+
      surfins = surfinswdir + surfinswdif
      surfinl = surfinl + surfinlwdif
@@ -6081 +6131 @@
 !--      for non-transparent surfaces, longwave albedo is 1 - emissivity
          surfoutll = (1._wp - emiss_surf) * surfinl
+
+         IF ( trace_fluxes_above >= 0._wp )  THEN
+            CALL radiation_print_debug_surf( 'surfoutll before reflective pass', surfoutll, refstep )
+            CALL radiation_print_debug_surf( 'surfoutsl before reflective pass', surfoutsl, refstep )
+         ENDIF
 
 #if defined( __parallel )
@@ -6148 +6203 @@
             mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc)
          ENDDO
+
+         IF ( trace_fluxes_above >= 0._wp )  THEN
+            CALL radiation_print_debug_surf( 'surfinl after reflected pass', surfinl, refstep )
+            CALL radiation_print_debug_surf( 'surfins after reflected pass', surfins, refstep )
+            CALL radiation_print_debug_pcb( 'pcbinlw after reflected pass', pcbinlw, refstep )
+            CALL radiation_print_debug_pcb( 'pcbinsw after reflected pass', pcbinsw, refstep )
+         ENDIF
 
          surfinsw = surfinsw  + surfins
@@ -6582 +6644 @@
 ! Description:
 ! ------------
-!> This subroutine splits direct and diffusion dw radiation
+!> This subroutine splits direct and diffusion dw radiation for RTM processing.
 !> It sould not be called in case the radiation model already does it
 !> It follows Boland, Ridley & Brown (2008)
@@ -6642 +6704 @@
     END SUBROUTINE calc_diffusion_radiation
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Print consecutive radiative extremes if requested to trace early radiation
+!> interaction instabilities.
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_print_debug_surf( description, values, step )
+
+       CHARACTER (LEN=*), INTENT(in)      ::  description
+       REAL(wp), DIMENSION(:), INTENT(in) ::  values
+       INTEGER(iwp), INTENT(in), OPTIONAL ::  step
+
+       CHARACTER (LEN=50)                 ::  location
+       CHARACTER (LEN=1024)               ::  debug_string
+       INTEGER                            ::  isurf
+       REAL(wp)                           ::  x
+
+       isurf = MAXLOC( values, DIM=1 )
+       x = values(isurf)
+       IF ( x < trace_fluxes_above )  RETURN
+
+       IF ( PRESENT( step ) )  THEN
+          WRITE( location, '(A," #",I0)' ) description, step
+       ELSE
+          location = description
+       ENDIF
+
+       WRITE( debug_string, '("Maximum of ",A50," = ",F12.1," at coords ' //  &
+                            'i=",I4,", j=",I4,", k=",I4,", d=",I1,". '    //  &
+                            'Alb=",F7.3,", emis=",F7.3)'                    ) &
+              location, x, surfl(ix,isurf), surfl(iy,isurf),                  &
+              surfl(iz,isurf), surfl(id,isurf), albedo_surf(isurf),           &
+              emiss_surf(isurf)
+       CALL debug_message( debug_string, 'info' )
+
+    END SUBROUTINE
+
+    SUBROUTINE radiation_print_debug_pcb( description, values, step )
+
+       CHARACTER (LEN=*), INTENT(in)      ::  description
+       REAL(wp), DIMENSION(:), INTENT(in) ::  values
+       INTEGER(iwp), INTENT(in), OPTIONAL ::  step
+
+       CHARACTER (LEN=50)                 ::  location
+       CHARACTER (LEN=1024)               ::  debug_string
+       INTEGER                            ::  ipcb
+       REAL(wp)                           ::  x
+
+       IF ( npcbl <= 0 )  RETURN
+       ipcb = MAXLOC( values, DIM=1 )
+       x = values(ipcb) / (dx*dy*dz(1))
+       IF ( x < trace_fluxes_above )  RETURN
+
+       IF ( PRESENT( step ) )  THEN
+          WRITE( location, '(A," #",I0)' ) description, step
+       ELSE
+          location = description
+       ENDIF
+
+       WRITE( debug_string, '("Maximum of ",A50," = ",F12.1," at coords ' //  &
+                            'i=",I4,", j=",I4,", k=",I4)'                   ) &
+              location, x, pcbl(ix,ipcb), pcbl(iy,ipcb), pcbl(iz,ipcb)
+       CALL debug_message( debug_string, 'info' )
+
+    END SUBROUTINE
+
+    SUBROUTINE radiation_print_debug_horz( description, values, step )
+
+       CHARACTER (LEN=*), INTENT(in)        ::  description
+       REAL(wp), DIMENSION(:,:), INTENT(in) ::  values
+       INTEGER(iwp), INTENT(in), OPTIONAL   ::  step
+
+       CHARACTER (LEN=50)                   ::  location
+       CHARACTER (LEN=1024)                 ::  debug_string
+       INTEGER, DIMENSION(2)                ::  ji
+       REAL(wp)                             ::  x
+
+       ji = MAXLOC( values )
+       x = values(ji(1),ji(2))
+       IF ( x < trace_fluxes_above )  RETURN
+
+       IF ( PRESENT( step ) )  THEN
+          WRITE( location, '(A," #",I0)' ) description, step
+       ELSE
+          location = description
+       ENDIF
+
+       WRITE( debug_string, '("Maximum of ",A50," = ",F12.1," at coords ' //  &
+                            'i=",I4,", j=",I4)'                             ) &
+              location, x, ji(2), ji(1)
+       CALL debug_message( debug_string, 'info' )
+
+    END SUBROUTINE
+
  END SUBROUTINE radiation_interaction
     
@@ -6647 +6803 @@
 ! Description:
 ! ------------
-!> This subroutine initializes structures needed for radiative transfer
-!> model. This model calculates transformation processes of the
+!> This subroutine initializes structures needed for Radiative Transfer
+!> Model (RTM). This model calculates transformation processes of the
 !> radiation inside urban and land canopy layer. The module includes also
 !> the interaction of the radiation with the resolved plant canopy.
-!>
-!> For more info. see Resler et al. 2017
-!>
-!> The new version 2.0 was radically rewriten, the discretization scheme
-!> has been changed. This new version significantly improves effectivity
-!> of the paralelization and the scalability of the model.
 !>
 !------------------------------------------------------------------------------!
@@ -7182 +7332 @@
 !> Calculates shape view factors (SVF), plant sink canopy factors (PCSF),
 !> sky-view factors, discretized path for direct solar radiation, MRT factors
-!> and other preprocessed data needed for radiation_interaction.
+!> and other preprocessed data needed for radiation_interaction inside RTM.
+!> This subroutine is called only one at the beginning of the simulation.
+!> The resulting factors can be stored to files and reused with other
+!> simulations utilizing the same surface and plant canopy structure.
 !------------------------------------------------------------------------------!
     SUBROUTINE radiation_calc_svf
@@ -7944 +8097 @@
 
         IF ( rad_angular_discretization )  THEN
-           IF ( debug_output )  CALL debug_message( 'Load svf from the structure array to plain arrays', 'info' )
+           IF ( debug_output )  THEN
+              WRITE( debug_string, '("Load ",I0," SVFs from the structure array to plain arrays")' ) nsvfl
+              CALL debug_message( debug_string, 'info' )
+           ENDIF
            ALLOCATE( svf(ndsvf,nsvfl) )
            ALLOCATE( svfsurf(idsvf,nsvfl) )
@@ -7958 +8114 @@
 
            !< load svf from the structure array to plain arrays
-           IF ( debug_output )  CALL debug_message( 'Load svf from the structure array to plain arrays', 'info' )
+           IF ( debug_output )  THEN
+              WRITE( debug_string, '("Load ",I0," SVFs from the structure array to plain arrays")' ) nsvfl
+              CALL debug_message( debug_string, 'info' )
+           ENDIF
            ALLOCATE( svf(ndsvf,nsvfl) )
            ALLOCATE( svfsurf(idsvf,nsvfl) )
@@ -8189 +8348 @@
             DEALLOCATE( pcsflt_l )
             DEALLOCATE( kpcsflt_l )
-            IF ( debug_output )  CALL debug_message( 'End of aggregate csf', 'info' )
+            IF ( debug_output )  THEN
+               WRITE( debug_string, '("Finished aggregating ",I0," CSFs.")') ncsfl
+               CALL debug_message( debug_string, 'info' )
+            ENDIF
             
         ENDIF
@@ -8212 +8374 @@
 ! ------------
 !> Raytracing for detecting obstacles and calculating compound canopy sink
-!> factors. (A simple obstacle detection would only need to process faces in
-!> 3 dimensions without any ordering.)
+!> factors for RTM. (A simple obstacle detection would only need to process
+!> faces in 3 dimensions without any ordering.)
 !> Assumtions:
 !> -----------
@@ -8410 +8572 @@
 ! ------------
 !> A new, more efficient version of ray tracing algorithm that processes a whole
-!> arc instead of a single ray.
+!> arc instead of a single ray (new in RTM version 2.5).
 !>
 !> In all comments, horizon means tangent of horizon angle, i.e.
@@ -8925 +9087 @@
 ! ------------
 !> Calculates apparent solar positions for all timesteps and stores discretized
-!> positions.
+!> positions for RTM.
 !------------------------------------------------------------------------------!
    SUBROUTINE radiation_presimulate_solar_pos
@@ -9013 +9175 @@
 ! Description:
 ! ------------
-!> Determines whether two faces are oriented towards each other. Since the
+!> Determines whether two faces are oriented towards each other in RTM. Since the
 !> surfaces follow the gird box surfaces, it checks first whether the two surfaces
 !> are directed in the same direction, then it checks if the two surfaces are
@@ -9074 +9236 @@
 ! Description:
 ! ------------
-!> Reads svf, svfsurf, csf, csfsurf and mrt factors data from saved file
-!> SVF means sky view factors and CSF means canopy sink factors
+!> Reads svf, svfsurf, csf, csfsurf and mrt factors data from saved file.
+!> This allows to skip their calculation during of RTM init phase.
+!> SVF means sky view factors and CSF means canopy sink factors.
 !------------------------------------------------------------------------------!
     SUBROUTINE radiation_read_svf
@@ -9240 +9403 @@
 ! ------------
 !> Subroutine stores svf, svfsurf, csf, csfsurf and mrt data to a file.
+!> The stored factors can be reused in future simulation with the same
+!> geometry structure of the surfaces and resolved plant canopy.
 !------------------------------------------------------------------------------!
     SUBROUTINE radiation_write_svf
@@ -9304 +9469 @@
 ! Description:
 ! ------------
-!> Block of auxiliary subroutines:
+!> Block of auxiliary subroutines for RTM:
 !> 1. quicksort and corresponding comparison
 !> 2. merge_and_grow_csf for implementation of "dynamical growing"
@@ -9440 +9605 @@
 ! Description:
 ! ------------
-!> Grows the CSF array exponentially after it is full. During that, the ray
-!> canopy sink factors with common source face and target plant canopy grid
-!> cell are merged together so that the size doesn't grow out of control.
+!> Grows the CSF array in RTM exponentially when it is full. During that,
+!> the ray canopy sink factors with common source face and target plant canopy
+!> grid cell are merged together so that the size doesn't grow out of control.
 !------------------------------------------------------------------------------!
     SUBROUTINE merge_and_grow_csf(newsize)