Index: palm/trunk/SOURCE/netcdf_data_input_mod.f90
===================================================================
--- palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 4177)
+++ palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 4178)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Implement input of external radiation forcing. Therefore, provide public
+! subroutines and variables. 
+! 
+! 4150 2019-08-08 20:00:47Z suehring
 ! Some variables are given the public attribute, in order to call netcdf input
 ! from single routines
@@ -608,5 +612,13 @@
        LOGICAL ::  from_file = .FALSE. !< flag indicating whether an input variable is available and read from file or default values are used
     END TYPE int_2d_32bit
-
+!
+!-- Define data type to read 1D real variables
+    TYPE real_1d
+       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
+
+       REAL(wp) ::  fill = -9999.9_wp                  !< fill value
+       
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  var     !< respective variable
+    END TYPE real_1d   
 !
 !-- Define data type to read 2D real variables
@@ -834,9 +846,12 @@
     CHARACTER(LEN=100) ::  input_file_vm      = 'PIDS_VM'      !< Name of file which comprises virtual measurement data
     
-    CHARACTER(LEN=25), ALLOCATABLE, DIMENSION(:)    ::  string_values  !< output of string variables read from netcdf input files
+    CHARACTER(LEN=25), ALLOCATABLE, DIMENSION(:) ::  string_values  !< output of string variables read from netcdf input files
+    CHARACTER(LEN=50), DIMENSION(:), ALLOCATABLE ::  vars_pids      !< variable in input file
 
     INTEGER(iwp)                                     ::  id_emis        !< NetCDF id of input file for chemistry emissions: TBD: It has to be removed
 
     INTEGER(iwp) ::  nc_stat         !< return value of nf90 function call
+    INTEGER(iwp) ::  num_var_pids    !< number of variables in file
+    INTEGER(iwp) ::  pids_id         !< file id
 
     LOGICAL ::  input_pids_static  = .FALSE.   !< Flag indicating whether Palm-input-data-standard file containing static information exists
@@ -949,5 +964,6 @@
 !
 !-- Public data structures
-    PUBLIC real_2d
+    PUBLIC real_1d,                                                            &
+           real_2d
 !
 !-- Public variables
@@ -956,12 +972,17 @@
            chem_emis, chem_emis_att, chem_emis_att_type, chem_emis_val_type,   &
            coord_ref_sys,                                                      &
-           init_3d, init_model, input_file_atts, input_file_static,            &
+           init_3d, init_model, input_file_atts,                               &
+           input_file_dynamic,                                                 &
+           input_file_static,                                                  &
            input_pids_static,                                                  &
            input_pids_dynamic, input_pids_vm, input_file_vm,                   &
            leaf_area_density_f, nest_offl,                                     &
+           num_var_pids,                                                       &
            pavement_pars_f, pavement_subsurface_pars_f, pavement_type_f,       &
+           pids_id,                                                            &
            root_area_density_lad_f, root_area_density_lsm_f, soil_pars_f,      &
            soil_type_f, street_crossing_f, street_type_f, surface_fraction_f,  &
            terrain_height_f, vegetation_pars_f, vegetation_type_f,             &
+           vars_pids,                                                          &
            water_pars_f, water_type_f
 !
Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4177)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4178)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! External radiation forcing implemented.
+! 
+! 4168 2019-08-16 13:50:17Z suehring
 ! Replace function get_topography_top_index by topo_top_ind
 ! 
@@ -667,4 +670,5 @@
         ONLY:  cloud_droplets, coupling_char,                                  &
                debug_output, debug_output_timestep, debug_string,              &
+               dt_3d,                                                          &
                dz, dt_spinup, end_time,                                        &
                humidity,                                                       &
@@ -704,6 +708,24 @@
 
     USE netcdf_data_input_mod,                                                 &
-        ONLY:  albedo_type_f, albedo_pars_f, building_type_f, pavement_type_f, &
-               vegetation_type_f, water_type_f
+        ONLY:  albedo_type_f,                                                  &
+               albedo_pars_f,                                                  &
+               building_type_f,                                                &
+               pavement_type_f,                                                &
+               vegetation_type_f,                                              &
+               water_type_f,                                                   &
+               char_fill,                                                      &
+               check_existence,                                                &
+               close_input_file,                                               &
+               get_attribute,                                                  &
+               get_variable,                                                   &
+               inquire_num_variables,                                          &
+               inquire_variable_names,                                         &
+               input_file_dynamic,                                             &
+               netcdf_data_input_get_dimension_length,                         &
+               num_var_pids,                                                   &
+               pids_id,                                                        &
+               open_read_file,                                                 &
+               real_1d,                                                        &
+               vars_pids
 
     USE plant_canopy_model_mod,                                                &
@@ -1236,4 +1258,10 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE            :: albedo_surf        !< albedo of the surface
     REAL(wp), DIMENSION(:), ALLOCATABLE            :: emiss_surf         !< emissivity of the wall surface
+!
+!-- External radiation
+    TYPE( real_1d ) ::  rad_lw_in_f     !< external incoming longwave radiation, from observation or model 
+    TYPE( real_1d ) ::  rad_sw_in_f     !< external incoming shortwave radiation, from observation or model 
+    TYPE( real_1d ) ::  rad_sw_in_dif_f !< external incoming shortwave radiation, diffuse part, from observation or model
+    TYPE( real_1d ) ::  time_rad_f      !< time dimension for external radiation, from observation or model 
 
 
@@ -1414,4 +1442,13 @@
           CASE ( 'rrtmg' )
              CALL radiation_rrtmg
+             
+          CASE ( 'external' )
+!
+!--          During spinup apply clear-sky model
+             IF ( time_since_reference_point < 0.0_wp )  THEN
+                CALL radiation_clearsky
+             ELSE
+                CALL radiation_external
+             ENDIF
 
           CASE DEFAULT
@@ -1804,5 +1841,6 @@
        IF ( radiation_scheme /= 'constant'   .AND.                             &
             radiation_scheme /= 'clear-sky'  .AND.                             &
-            radiation_scheme /= 'rrtmg' )  THEN
+            radiation_scheme /= 'rrtmg'      .AND.                             &
+            radiation_scheme /= 'external' )  THEN
           message_string = 'unknown radiation_scheme = '//                     &
                            TRIM( radiation_scheme )
@@ -1899,4 +1937,5 @@
        INTEGER(iwp) ::  l         !< running index for orientation of vertical surfaces 
        INTEGER(iwp) ::  m         !< running index for surface elements
+       INTEGER(iwp) ::  ntime = 0 !< number of available external radiation timesteps 
 #if defined( __rrtmg )
        INTEGER(iwp) ::  ind_type  !< running index for subgrid-surface tiles
@@ -2113,7 +2152,6 @@
 
        IF ( radiation_scheme == 'clear-sky'  .OR.                              &
-            radiation_scheme == 'constant')  THEN
-
-
+            radiation_scheme == 'constant'   .OR.                              &
+            radiation_scheme == 'external' )  THEN
 !
 !--       Allocate arrays for incoming/outgoing short/longwave radiation
@@ -2896,5 +2934,129 @@
 #endif
        ENDIF
-
+!
+!--    Initializaion actions exclusively required for external 
+!--    radiation forcing
+       IF ( radiation_scheme == 'external' )  THEN
+!
+!--       Open the radiation input file
+#if defined( __netcdf )
+          CALL open_read_file( TRIM( input_file_dynamic ) //                   &
+                               TRIM( coupling_char ),                          &
+                               pids_id )
+                               
+          CALL inquire_num_variables( pids_id, num_var_pids )
+!         
+!--       Allocate memory to store variable names and read them
+          ALLOCATE( vars_pids(1:num_var_pids) )
+          CALL inquire_variable_names( pids_id, vars_pids )
+          
+!         
+!--       Input time dimension.
+          IF ( check_existence( vars_pids, 'time_rad' ) )  THEN
+             CALL netcdf_data_input_get_dimension_length( pids_id,             &
+                                                          ntime, 'time_rad' )
+          
+             ALLOCATE( time_rad_f%var(0:ntime-1) )
+!                                                                                 
+!--          Read variable                    
+             CALL get_variable( pids_id, 'time_rad', time_rad_f%var )
+                                
+             time_rad_f%from_file = .TRUE.
+          ENDIF           
+!         
+!--       Input shortwave downwelling.
+          IF ( check_existence( vars_pids, 'rad_sw_in' ) )  THEN
+             ALLOCATE( rad_sw_in_f%var(0:ntime-1) )
+!         
+!--          Read _FillValue attribute
+             CALL get_attribute( pids_id, char_fill, rad_sw_in_f%fill,         &
+                                 .FALSE., 'rad_sw_in' )
+!                                                                                 
+!--          Read variable                    
+             CALL get_variable( pids_id, 'rad_sw_in', rad_sw_in_f%var )
+                                
+             rad_sw_in_f%from_file = .TRUE.
+          ENDIF
+!         
+!--       Input longwave downwelling.
+          IF ( check_existence( vars_pids, 'rad_lw_in' ) )  THEN
+             ALLOCATE( rad_lw_in_f%var(0:ntime-1) )
+!         
+!--          Read _FillValue attribute
+             CALL get_attribute( pids_id, char_fill, rad_lw_in_f%fill,         &
+                                 .FALSE., 'rad_lw_in' )
+!                                                                                 
+!--          Read variable                    
+             CALL get_variable( pids_id, 'rad_lw_in', rad_lw_in_f%var )
+                                
+             rad_lw_in_f%from_file = .TRUE.
+          ENDIF
+!         
+!--       Input shortwave downwelling, diffuse part.
+          IF ( check_existence( vars_pids, 'rad_sw_in_dif' ) )  THEN
+             ALLOCATE( rad_sw_in_dif_f%var(0:ntime-1) )
+!         
+!--          Read _FillValue attribute
+             CALL get_attribute( pids_id, char_fill, rad_sw_in_dif_f%fill,     &
+                                 .FALSE., 'rad_sw_in_dif' )
+!                                                                                 
+!--          Read variable                    
+             CALL get_variable( pids_id, 'rad_sw_in_dif', rad_sw_in_dif_f%var )
+                                
+             rad_sw_in_dif_f%from_file = .TRUE.
+          ENDIF
+!         
+!--       Finally, close the input file.
+          CALL close_input_file( pids_id )
+#endif
+!
+!--       Make some consistency checks.
+          IF ( .NOT. rad_sw_in_f%from_file  .OR.                               &
+               .NOT. rad_lw_in_f%from_file )  THEN
+             message_string = 'In case of external radiation forcing ' //      &
+                              'both, rad_sw_in and rad_lw_in are required.'
+             CALL message( 'radiation_init', 'PA0195', 1, 2, 0, 6, 0 )
+          ENDIF
+          
+          IF ( .NOT. time_rad_f%from_file )  THEN
+             message_string = 'In case of external radiation forcing ' //      &
+                              'dimension time_rad is required.'
+             CALL message( 'radiation_init', 'PA0196', 1, 2, 0, 6, 0 )
+          ENDIF
+          
+          IF ( ANY( rad_sw_in_f%var ==  rad_sw_in_f%fill ) )  THEN
+             message_string = 'External radiation forcing: rad_sw_in must ' // &
+                              'not contain any fill values'
+             CALL message( 'radiation_init', 'PA0197', 1, 2, 0, 6, 0 )
+          ENDIF
+          
+          IF ( ANY( rad_lw_in_f%var ==  rad_lw_in_f%fill ) )  THEN
+             message_string = 'External radiation forcing: rad_lw_in must ' // &
+                              'not contain any fill values'
+             CALL message( 'radiation_init', 'PA0198', 1, 2, 0, 6, 0 )
+          ENDIF
+          
+          IF ( rad_sw_in_dif_f%from_file )  THEN
+             IF ( ANY( rad_sw_in_dif_f%var ==  rad_sw_in_dif_f%fill ) )  THEN
+                message_string = 'External radiation forcing: ' //             &
+                                 'rad_lw_in_dif must not contain any fill ' // &
+                                 'values'
+                CALL message( 'radiation_init', 'PA0199', 1, 2, 0, 6, 0 )
+             ENDIF
+          ENDIF
+          
+          IF ( time_rad_f%var(0) /= 0.0_wp )  THEN
+             message_string = 'External radiation forcing: first point in ' // &
+                              'time is /= 0.0.'
+             CALL message( 'radiation_init', 'PA0313', 1, 2, 0, 6, 0 )
+          ENDIF
+          
+          IF ( end_time - spinup_time > time_rad_f%var(ntime-1) )  THEN
+             message_string = 'External radiation forcing does not cover ' //  &
+                              'the entire simulation time.'
+             CALL message( 'radiation_init', 'PA0314', 1, 2, 0, 6, 0 )
+          ENDIF
+          
+       ENDIF
 !
 !--    Perform user actions if required
@@ -2913,10 +3075,19 @@
           CASE ( 'constant' )
              CALL radiation_constant
+             
+          CASE ( 'external' )
+!
+!--          During spinup apply clear-sky model
+             IF ( time_since_reference_point < 0.0_wp )  THEN
+                CALL radiation_clearsky
+             ELSE
+                CALL radiation_external
+             ENDIF
 
           CASE DEFAULT
 
        END SELECT
-
-! readjust date and time to its initial value
+!
+!--    Readjust date and time to its initial value
        CALL init_date_and_time
 
@@ -2958,4 +3129,181 @@
 
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> A simple clear sky radiation model
+!------------------------------------------------------------------------------!
+    SUBROUTINE radiation_external
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  l   !< running index for surface orientation
+       INTEGER(iwp) ::  t   !< index of current timestep
+       INTEGER(iwp) ::  tm  !< index of previous timestep 
+       
+       REAL(wp) ::  fac_dt     !< interpolation factor
+       REAL(wp) ::  lw_in      !< downwelling longwave radiation, interpolated value     
+       REAL(wp) ::  sw_in      !< downwelling shortwave radiation, interpolated value       
+       REAL(wp) ::  sw_in_dif  !< downwelling diffuse shortwave radiation, interpolated value   
+
+       TYPE(surf_type), POINTER ::  surf !< pointer on respective surface type, used to generalize routine   
+
+!
+!--    Calculate current zenith angle
+       CALL calc_zenith
+!
+!--    Interpolate external radiation on current timestep 
+       IF ( time_since_reference_point  <= 0.0_wp )  THEN
+          t      = 0
+          tm     = 0
+          fac_dt = 0
+       ELSE
+          t = 0
+          DO WHILE ( time_rad_f%var(t) <= time_since_reference_point )
+             t = t + 1
+          ENDDO
+          
+          tm = MAX( t-1, 0 )
+          
+          fac_dt = ( time_since_reference_point - time_rad_f%var(tm) + dt_3d ) &
+                 / ( time_rad_f%var(t)  - time_rad_f%var(tm) )
+          fac_dt = MIN( 1.0_wp, fac_dt )
+       ENDIF
+             
+       sw_in = ( 1.0_wp - fac_dt ) * rad_sw_in_f%var(tm)                       &
+                          + fac_dt * rad_sw_in_f%var(t)
+                                         
+       lw_in = ( 1.0_wp - fac_dt ) * rad_lw_in_f%var(tm)                       &
+                          + fac_dt * rad_lw_in_f%var(t)
+!
+!--    Limit shortwave incoming radiation to positive values, in order 
+!--    to overcome possible observation errors.
+       sw_in = MAX( 0.0_wp, sw_in )
+       sw_in = MERGE( sw_in, 0.0_wp, sun_up )
+!
+!--    Call clear-sky calculation for each surface orientation. 
+!--    First, horizontal surfaces
+       surf => surf_lsm_h
+       CALL radiation_external_surf
+       surf => surf_usm_h
+       CALL radiation_external_surf
+!
+!--    Vertical surfaces
+       DO  l = 0, 3
+          surf => surf_lsm_v(l)
+          CALL radiation_external_surf
+          surf => surf_usm_v(l)
+          CALL radiation_external_surf
+       ENDDO
+!
+!--    In case of average radiation check if external diffuse shortwave 
+!--    radiation is available and, if required, interpolate it onto the
+!--    respective arrays.
+       IF ( average_radiation )  THEN
+          IF ( rad_sw_in_dif_f%from_file )  THEN 
+             sw_in_dif= ( 1.0_wp - fac_dt ) * rad_sw_in_dif_f%var(tm)          &
+                                   + fac_dt * rad_sw_in_dif_f%var(t)
+                                         
+             IF ( ALLOCATED( rad_sw_in_diff ) )  rad_sw_in_diff = sw_in_dif
+             IF ( ALLOCATED( rad_sw_in_dir  ) )  rad_sw_in_dir  = sw_in - sw_in_dif
+!
+!--          Diffuse longwave radiation equals the total downwelling longwaver
+!--          radiation
+             IF ( ALLOCATED( rad_lw_in_diff ) )  rad_lw_in_diff = lw_in
+          ENDIF
+       ENDIF
+       
+       CONTAINS
+
+          SUBROUTINE radiation_external_surf
+
+             USE control_parameters
+          
+             IMPLICIT NONE
+
+             INTEGER(iwp) ::  i         !< grid index along x-dimension   
+             INTEGER(iwp) ::  j         !< grid index along y-dimension  
+             INTEGER(iwp) ::  k         !< grid index along z-dimension    
+             INTEGER(iwp) ::  m         !< running index for surface elements       
+
+             IF ( surf%ns < 1 )  RETURN
+                          
+             surf%rad_sw_in = sw_in                                         
+             surf%rad_lw_in = lw_in
+
+             IF ( average_radiation )  THEN
+             
+                surf%rad_sw_out = albedo_urb * surf%rad_sw_in
+                
+                surf%rad_lw_out = emissivity_urb * sigma_sb * t_rad_urb**4     &
+                                 + ( 1.0_wp - emissivity_urb ) * surf%rad_lw_in
+                
+                surf%rad_net = surf%rad_sw_in - surf%rad_sw_out                &
+                             + surf%rad_lw_in - surf%rad_lw_out
+                
+                surf%rad_lw_out_change_0 = 4.0_wp * emissivity_urb * sigma_sb  &
+                                                  * t_rad_urb**3
+!
+!--          Calculate upwelling radiation fluxes and net radiation for each 
+!--          surface element depending on its properties.
+             ELSE
+
+                DO  m = 1, surf%ns
+                   i = surf%i(m)
+                   j = surf%j(m)
+                   k = surf%k(m)
+!
+!--                Weighted average according to surface fraction. 
+                   surf%rad_sw_out(m) = ( surf%frac(ind_veg_wall,m)  *         &
+                                          surf%albedo(ind_veg_wall,m)          &
+                                        + surf%frac(ind_pav_green,m) *         &
+                                          surf%albedo(ind_pav_green,m)         &
+                                        + surf%frac(ind_wat_win,m)   *         &
+                                          surf%albedo(ind_wat_win,m) )         &
+                                        * surf%rad_sw_in(m)
+
+                   surf%rad_lw_out(m) = ( surf%frac(ind_veg_wall,m)  *         &
+                                          surf%emissivity(ind_veg_wall,m)      &
+                                        + surf%frac(ind_pav_green,m) *         &
+                                          surf%emissivity(ind_pav_green,m)     &
+                                        + surf%frac(ind_wat_win,m)   *         &
+                                          surf%emissivity(ind_wat_win,m)       &
+                                        )                                      &
+                                        * sigma_sb                             &
+                                        * ( surf%pt_surface(m) * exner(k) )**4
+
+                   surf%rad_lw_out_change_0(m) =                               &
+                                      ( surf%frac(ind_veg_wall,m)  *           &
+                                        surf%emissivity(ind_veg_wall,m)        &
+                                      + surf%frac(ind_pav_green,m) *           &
+                                        surf%emissivity(ind_pav_green,m)       &
+                                      + surf%frac(ind_wat_win,m)   *           &
+                                        surf%emissivity(ind_wat_win,m)         &
+                                      ) * 4.0_wp * sigma_sb                    &
+                                      * ( surf%pt_surface(m) * exner(k) )**3
+
+                   surf%rad_net(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m)    &
+                                   + surf%rad_lw_in(m) - surf%rad_lw_out(m)
+
+                ENDDO
+
+             ENDIF
+                
+!
+!--          Fill out values in radiation arrays
+             DO  m = 1, surf%ns
+                i = surf%i(m)
+                j = surf%j(m)
+                rad_sw_out(0,j,i) = surf%rad_sw_out(m)
+                rad_lw_out(0,j,i) = surf%rad_lw_out(m)
+             ENDDO
+             
+             rad_sw_in(0,:,:) = sw_in
+             rad_lw_in(0,:,:) = lw_in
+ 
+          END SUBROUTINE radiation_external_surf
+
+    END SUBROUTINE radiation_external    
+    
 !------------------------------------------------------------------------------!
 ! Description:
@@ -3061,4 +3409,6 @@
                 
                 surf%rad_lw_in  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k+1))**4
+                
+                print*, "clear sky", emissivity_atm_clsky * sigma_sb * (pt1 * exner(k+1))**4
 
                 surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4   &
@@ -3136,7 +3486,7 @@
                 i = surf%i(m)
                 j = surf%j(m)
-                rad_sw_in(0,j,i) = surf%rad_sw_in(m)
+                rad_sw_in(0,j,i)  = surf%rad_sw_in(m)
                 rad_sw_out(0,j,i) = surf%rad_sw_out(m)
-                rad_lw_in(0,j,i) = surf%rad_lw_in(m)
+                rad_lw_in(0,j,i)  = surf%rad_lw_in(m)
                 rad_lw_out(0,j,i) = surf%rad_lw_out(m)
              ENDDO
@@ -3355,4 +3705,6 @@
           IF ( .NOT. lw_radiation )  WRITE( io, 10 )
           IF ( .NOT. sw_radiation )  WRITE( io, 11 )
+       ELSEIF ( radiation_scheme == "external" )  THEN
+          WRITE( io, 14 )
        ENDIF 
 
@@ -3390,4 +3742,5 @@
 13 FORMAT (/'    Albedo is set individually for each xy-location, according ', &
                  'to given surface type.')
+14 FORMAT ('    --> External radiation forcing is used')
 
 
@@ -5328,8 +5681,10 @@
          ENDIF
      ENDIF
-
-!--  if radiation scheme is not RRTMG, split diffusion and direct part of the solar downward radiation
-!--  comming from radiation model and store it in 2D arrays
-     IF (  radiation_scheme /= 'rrtmg'  )  CALL calc_diffusion_radiation
+!
+!--  Split downwelling shortwave radiation into a diffuse and a direct part. 
+!--  Note, if radiation scheme is RRTMG or diffuse radiation is externally
+!--  prescribed, this is not required. 
+     IF (  radiation_scheme /= 'rrtmg'  .AND.                                  &
+           .NOT. rad_sw_in_dif_f%from_file )  CALL calc_diffusion_radiation
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -6118,14 +6473,17 @@
     SUBROUTINE calc_diffusion_radiation 
     
-        REAL(wp), PARAMETER                          :: lowest_solarUp = 0.1_wp  !< limit the sun elevation to protect stability of the calculation
-        INTEGER(iwp)                                 :: i, j
-        REAL(wp)                                     ::  year_angle              !< angle
-        REAL(wp)                                     ::  etr                     !< extraterestrial radiation
-        REAL(wp)                                     ::  corrected_solarUp       !< corrected solar up radiation
-        REAL(wp)                                     ::  horizontalETR           !< horizontal extraterestrial radiation
-        REAL(wp)                                     ::  clearnessIndex          !< clearness index
-        REAL(wp)                                     ::  diff_frac               !< diffusion fraction of the radiation
-
+        INTEGER(iwp)         ::  i                       !< grid index x-direction
+        INTEGER(iwp)         ::  j                       !< grid index y-direction
         
+        REAL(wp)             ::  year_angle              !< angle
+        REAL(wp)             ::  etr                     !< extraterestrial radiation
+        REAL(wp)             ::  corrected_solarUp       !< corrected solar up radiation
+        REAL(wp)             ::  horizontalETR           !< horizontal extraterestrial radiation
+        REAL(wp)             ::  clearnessIndex          !< clearness index
+        REAL(wp)             ::  diff_frac               !< diffusion fraction of the radiation
+        
+        REAL(wp), PARAMETER  ::  lowest_solarUp = 0.1_wp  !< limit the sun elevation to protect stability of the calculation
+print*, "in diffusion_radiation"
+!  
 !--     Calculate current day and time based on the initial values and simulation time
         year_angle = ( (day_of_year_init * 86400) + time_utc_init              &
@@ -6138,10 +6496,9 @@
                           0.000719_wp * cos(2.0_wp * year_angle) +             &
                           0.000077_wp * sin(2.0_wp * year_angle))
-        
+!       
 !--    
 !--     Under a very low angle, we keep extraterestrial radiation at
 !--     the last small value, therefore the clearness index will be pushed
-!--     towards 0 while keeping full continuity.
-!--    
+!--     towards 0 while keeping full continuity.    
         IF ( cos_zenith <= lowest_solarUp )  THEN
             corrected_solarUp = lowest_solarUp
@@ -11451,5 +11808,6 @@
           IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 ALLOCATE( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
              ELSE
@@ -11459,5 +11817,6 @@
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 READ ( 13 )  tmp_3d2
                 rad_lw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
@@ -11473,5 +11832,6 @@
           IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 ALLOCATE( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
              ELSE
@@ -11481,5 +11841,6 @@
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 READ ( 13 )  tmp_3d2
                 rad_lw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
@@ -11495,5 +11856,6 @@
           IF ( .NOT. ALLOCATED( rad_lw_out ) )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 ALLOCATE( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
              ELSE
@@ -11503,5 +11865,6 @@
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 READ ( 13 )  tmp_3d2
                 rad_lw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =  &
@@ -11517,5 +11880,6 @@
           IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 ALLOCATE( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
              ELSE
@@ -11525,5 +11889,6 @@
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 READ ( 13 )  tmp_3d2
                 rad_lw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) &
@@ -11571,5 +11936,6 @@
           IF ( .NOT. ALLOCATED( rad_sw_in ) )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 ALLOCATE( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
              ELSE
@@ -11579,5 +11945,6 @@
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 READ ( 13 )  tmp_3d2
                 rad_sw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
@@ -11593,5 +11960,6 @@
           IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 ALLOCATE( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
              ELSE
@@ -11601,5 +11969,6 @@
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 READ ( 13 )  tmp_3d2
                 rad_sw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
@@ -11615,5 +11984,6 @@
           IF ( .NOT. ALLOCATED( rad_sw_out ) )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 ALLOCATE( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
              ELSE
@@ -11623,5 +11993,6 @@
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 READ ( 13 )  tmp_3d2
                 rad_sw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =  &
@@ -11637,5 +12008,6 @@
           IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 ALLOCATE( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
              ELSE
@@ -11645,5 +12017,6 @@
           IF ( k == 1 )  THEN
              IF ( radiation_scheme == 'clear-sky'  .OR.                    &
-                  radiation_scheme == 'constant')  THEN
+                  radiation_scheme == 'constant'   .OR.                    &
+                  radiation_scheme == 'external' )  THEN
                 READ ( 13 )  tmp_3d2
                 rad_sw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) &
