Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4492)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4493)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! Avoid unstable direct normal radiation near horizon
+! 
+! 4481 2020-03-31 18:55:54Z maronga
 ! use statement for exchange horiz added
 ! 
@@ -454,5 +457,6 @@
                 sky_trans,                       & !< sky transmissivity
                 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)
+                trace_fluxes_above = -1.0_wp,    & !< NAMELIST option for debug tracing of large radiative fluxes (W/m2;W/m3)
+                min_stable_coszen = 0.0262_wp      !< 1.5 deg above horizon, eliminates most of circumsolar
 
     INTEGER(iwp) ::  day_of_year   !< day of the current year
@@ -5806,4 +5810,5 @@
      REAL(wp)                          ::  area_surf          !< total area of surfaces in all processor
      REAL(wp)                          ::  area_hor           !< total horizontal area of domain in all processor
+     REAL(wp)                          ::  sun_direct_factor  !< factor for direct normal radiation from direct horizontal
 #if defined( __parallel )
      REAL(wp), DIMENSION(1:7)          ::  combine_allreduce   !< dummy array used to combine several MPI_ALLREDUCE calls
@@ -6019,9 +6024,14 @@
         mrtinlw(imrt) = mrtsky(imrt) * rad_lw_in_diff(j,i)
      ENDDO
-
-     !-- direct radiation
+!
+!--  Direct radiation
      IF ( cos_zenith > 0 )  THEN
-        !--Identify solar direction vector (discretized number) 1)
-        !--
+!
+!--     To avoid numerical instability near horizon depending on what direct
+!--     radiation is used (slightly different zenith angle, considering
+!--     circumsolar etc.), we use a minimum value for cos_zenith
+        sun_direct_factor = 1._wp / MAX(min_stable_coszen, cos_zenith)
+!
+!--     Identify solar direction vector (discretized number) (1)
         j = FLOOR(ACOS(cos_zenith) / pi * raytrace_discrete_elevs)
         i = MODULO(NINT(ATAN2(sun_dir_lon, sun_dir_lat)               &
@@ -6034,5 +6044,5 @@
            i = surfl(ix, isurf)
            surfinswdir(isurf) = rad_sw_in_dir(j,i) *                        &
-                costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / cos_zenith
+                costheta(surfl(id, isurf)) * dsitrans(isurf, isd) * sun_direct_factor
         ENDDO
 !
@@ -6042,5 +6052,5 @@
            i = mrtbl(ix, imrt)
            mrtinsw(imrt) = mrtinsw(imrt) + mrtdsit(imrt, isd) * rad_sw_in_dir(j,i) &
-                                     / cos_zenith / 4.0_wp ! normal to sphere
+                                     * sun_direct_factor / 4.0_wp ! normal to sphere
         ENDDO
      ENDIF
@@ -6085,5 +6095,5 @@
                    pc_abs_frac = 1.0_wp - exp(pc_abs_eff * lad_s(k,j,i))
 !
-!--                isd has already been established, see 1)
+!--                isd has already been established, see (1)
                    pcbinswdir(ipcgb) = rad_sw_in_dir(j, i) * pc_box_area &
                                        * pc_abs_frac * dsitransc(ipcgb, isd)
@@ -7867,5 +7877,6 @@
 
 !--
-!--     Raytrace to canopy boxes to fill dsitransc TODO optimize
+!--     Raytrace to canopy boxes to fill dsitransc
+!--     TODO: consider replacing by DSI rays toward surfaces
         dsitransc(:,:) = 0._wp
         az0 = 0._wp
Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4492)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4493)
@@ -28,4 +28,9 @@
 ! -----------------
 ! $Id$
+! J.Resler, 2020/03/19
+! - remove reading of deprecated input parameters c_surface and lambda_surf
+! - and calculate them from parameters of the outer wall/roof layer
+! 
+! 4481 2020-03-31 18:55:54Z maronga
 ! use statement for exchange horiz added
 ! 
@@ -219,22 +224,10 @@
 !> Further work:
 !> -------------
-!> 1. Remove global arrays surfouts, surfoutl and only keep track of radiosity
-!>    from surfaces that are visible from local surfaces (i.e. there is a SVF
-!>    where target is local). To do that, radiosity will be exchanged after each
-!>    reflection step using MPI_Alltoall instead of current MPI_Allgather.
-!>
-!> 2. Temporarily large values of surface heat flux can be observed, up to
-!>    1.2 Km/s, which seem to be not realistic.
-!>
 !> @todo Output of _av variables in case of restarts
 !> @todo Revise flux conversion in energy-balance solver
-!> @todo Check optimizations for RMA operations
-!> @todo Alternatives for MPI_WIN_ALLOCATE? (causes problems with openmpi)
-!> @todo Check for load imbalances in CPU measures, e.g. for exchange_horiz_prog
-!>       factor 3 between min and max time
 !> @todo Check divisions in wtend (etc.) calculations for possible division
 !>       by zero, e.g. in case fraq(0,m) + fraq(1,m) = 0?!
 !> @todo Use unit 90 for OPEN/CLOSE of input files (FK)
-!> @todo Move plant canopy stuff into plant canopy code
+!> @todo remove reading of old csv inputs
 !------------------------------------------------------------------------------!
  MODULE urban_surface_mod
@@ -375,7 +368,4 @@
     INTEGER(iwp) ::  ind_alb_win_gfl       = 77   !< index in input list for albedo_type of window fraction ground floor level
     INTEGER(iwp) ::  ind_alb_win_r         = 115  !< index in input list for albedo_type of window fraction roof
-    INTEGER(iwp) ::  ind_c_surface         = 45   !< index in input list for heat capacity wall surface
-    INTEGER(iwp) ::  ind_c_surface_green   = 48   !< index in input list for heat capacity green surface
-    INTEGER(iwp) ::  ind_c_surface_win     = 47   !< index in input list for heat capacity window surface
     INTEGER(iwp) ::  ind_emis_wall_agfl    = 14   !< index in input list for wall emissivity, above ground floor level
     INTEGER(iwp) ::  ind_emis_wall_gfl     = 32   !< index in input list for wall emissivity, ground floor level
@@ -425,7 +415,4 @@
     INTEGER(iwp) ::  ind_lai_w_agfl        = 5    !< index in input list for LAI on wall, above ground floor level
     INTEGER(iwp) ::  ind_lai_w_gfl         = 25   !< index in input list for LAI on wall, ground floor level
-    INTEGER(iwp) ::  ind_lambda_surf       = 46   !< index in input list for thermal conductivity of wall surface
-    INTEGER(iwp) ::  ind_lambda_surf_green = 50   !< index in input list for thermal conductivity of green surface
-    INTEGER(iwp) ::  ind_lambda_surf_win   = 49   !< index in input list for thermal conductivity of window surface
     INTEGER(iwp) ::  ind_tc1_agfl          = 9    !< index in input list for thermal conductivity at first wall layer,
                                                   !< above ground floor level
@@ -3438,5 +3425,5 @@
     END SUBROUTINE usm_init_material_model
 
- 
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -3687,11 +3674,4 @@
            surf_usm_h%zw_window(nzb_wall+3,m)       = building_pars(ind_thick_4_win_r,building_type)
 
-           surf_usm_h%c_surface(m)           = building_pars(ind_c_surface,building_type)  
-           surf_usm_h%lambda_surf(m)         = building_pars(ind_lambda_surf,building_type)  
-           surf_usm_h%c_surface_green(m)     = building_pars(ind_c_surface_green,building_type)  
-           surf_usm_h%lambda_surf_green(m)   = building_pars(ind_lambda_surf_green,building_type)  
-           surf_usm_h%c_surface_window(m)    = building_pars(ind_c_surface_win,building_type)  
-           surf_usm_h%lambda_surf_window(m)  = building_pars(ind_lambda_surf_win,building_type)  
-           
            surf_usm_h%green_type_roof(m)     = building_pars(ind_green_type_roof,building_type)
 
@@ -3849,11 +3829,4 @@
               surf_usm_v(l)%zw_window(nzb_wall+3,m)       = building_pars(ind_thick_4_win,building_type)
 
-              surf_usm_v(l)%c_surface(m)           = building_pars(ind_c_surface,building_type)  
-              surf_usm_v(l)%lambda_surf(m)         = building_pars(ind_lambda_surf,building_type)
-              surf_usm_v(l)%c_surface_green(m)     = building_pars(ind_c_surface_green,building_type)  
-              surf_usm_v(l)%lambda_surf_green(m)   = building_pars(ind_lambda_surf_green,building_type)
-              surf_usm_v(l)%c_surface_window(m)    = building_pars(ind_c_surface_win,building_type)  
-              surf_usm_v(l)%lambda_surf_window(m)  = building_pars(ind_lambda_surf_win,building_type)
-
            ENDDO
         ENDDO
@@ -3953,11 +3926,4 @@
                  surf_usm_h%zw_window(nzb_wall+3,m)       = building_pars(ind_thick_4_win_r,st)
 
-                 surf_usm_h%c_surface(m)           = building_pars(ind_c_surface,st)  
-                 surf_usm_h%lambda_surf(m)         = building_pars(ind_lambda_surf,st)
-                 surf_usm_h%c_surface_green(m)     = building_pars(ind_c_surface_green,st)  
-                 surf_usm_h%lambda_surf_green(m)   = building_pars(ind_lambda_surf_green,st)
-                 surf_usm_h%c_surface_window(m)    = building_pars(ind_c_surface_win,st)  
-                 surf_usm_h%lambda_surf_window(m)  = building_pars(ind_lambda_surf_win,st)
-                 
                  surf_usm_h%green_type_roof(m)     = building_pars(ind_green_type_roof,st)
 
@@ -4122,12 +4088,4 @@
                     surf_usm_v(l)%zw_window(nzb_wall+3,m)       = building_pars(ind_thick_4_win,st)
 
-                    surf_usm_v(l)%c_surface(m)           = building_pars(ind_c_surface,st)  
-                    surf_usm_v(l)%lambda_surf(m)         = building_pars(ind_lambda_surf,st) 
-                    surf_usm_v(l)%c_surface_green(m)     = building_pars(ind_c_surface_green,st)  
-                    surf_usm_v(l)%lambda_surf_green(m)   = building_pars(ind_lambda_surf_green,st) 
-                    surf_usm_v(l)%c_surface_window(m)    = building_pars(ind_c_surface_win,st)  
-                    surf_usm_v(l)%lambda_surf_window(m)  = building_pars(ind_lambda_surf_win,st) 
-
-
                  ENDIF
               ENDDO
@@ -4275,5 +4233,5 @@
                                     building_pars_f%pars_xy(ind_tc1,j,i)         
                  surf_usm_h%lambda_h(nzb_wall+1,m) =                           &
-                                    building_pars_f%pars_xy(ind_tc1,j,i)        
+                                    building_pars_f%pars_xy(ind_tc1,j,i)
               ENDIF
               IF ( building_pars_f%pars_xy(ind_tc2,j,i) /=                     &  
@@ -4292,5 +4250,5 @@
                                      building_pars_f%pars_xy(ind_tc1,j,i)         
                  surf_usm_h%lambda_h_green(nzb_wall+1,m) =                     &
-                                     building_pars_f%pars_xy(ind_tc1,j,i)        
+                                     building_pars_f%pars_xy(ind_tc1,j,i)
               ENDIF
               IF ( building_pars_f%pars_xy(ind_tc2,j,i) /=                     & 
@@ -4309,5 +4267,5 @@
                                      building_pars_f%pars_xy(ind_tc1,j,i)         
                  surf_usm_h%lambda_h_window(nzb_wall+1,m) =                    &
-                                     building_pars_f%pars_xy(ind_tc1,j,i)        
+                                     building_pars_f%pars_xy(ind_tc1,j,i)
               ENDIF
               IF ( building_pars_f%pars_xy(ind_tc2,j,i) /=                     &     
@@ -4415,14 +4373,4 @@
                  surf_usm_h%zw_green(nzb_wall+3,m) =                           &
                                    building_pars_f%pars_xy(ind_thick_4_agfl,j,i)
-
-              IF ( building_pars_f%pars_xy(ind_c_surface,j,i) /=               &       
-                   building_pars_f%fill )                                      & 
-                 surf_usm_h%c_surface(m) =                                     &
-                                    building_pars_f%pars_xy(ind_c_surface,j,i)
-                 
-              IF ( building_pars_f%pars_xy(ind_lambda_surf,j,i) /=             &        
-                   building_pars_f%fill )                                      &
-                 surf_usm_h%lambda_surf(m) =                                   &
-                                    building_pars_f%pars_xy(ind_lambda_surf,j,i)
            ENDDO
 
@@ -4564,5 +4512,5 @@
                                           building_pars_f%pars_xy(ind_tc1,j,i)   
                        surf_usm_v(l)%lambda_h(nzb_wall+1,m) =                  &
-                                          building_pars_f%pars_xy(ind_tc1,j,i)  
+                                          building_pars_f%pars_xy(ind_tc1,j,i)
                     ENDIF
                     IF ( building_pars_f%pars_xy(ind_tc2,j,i) /=               &  
@@ -4581,5 +4529,5 @@
                                            building_pars_f%pars_xy(ind_tc1,j,i)   
                        surf_usm_v(l)%lambda_h_green(nzb_wall+1,m) =            &
-                                           building_pars_f%pars_xy(ind_tc1,j,i)  
+                                           building_pars_f%pars_xy(ind_tc1,j,i)
                     ENDIF
                     IF ( building_pars_f%pars_xy(ind_tc2,j,i) /=               & 
@@ -4598,5 +4546,5 @@
                                      building_pars_f%pars_xy(ind_tc1,j,i)         
                        surf_usm_v(l)%lambda_h_window(nzb_wall+1,m) =           &
-                                     building_pars_f%pars_xy(ind_tc1,j,i)        
+                                     building_pars_f%pars_xy(ind_tc1,j,i)
                     ENDIF
                     IF ( building_pars_f%pars_xy(ind_tc2,j,i) /=               &     
@@ -4707,14 +4655,4 @@
                                    building_pars_f%pars_xy(ind_thick_4_agfl,j,i)
                     
-                    IF ( building_pars_f%pars_xy(ind_c_surface,j,i) /=         &       
-                         building_pars_f%fill )                                & 
-                       surf_usm_v(l)%c_surface(m) =                            &
-                                     building_pars_f%pars_xy(ind_c_surface,j,i)
-                       
-                    IF ( building_pars_f%pars_xy(ind_lambda_surf,j,i) /=       &        
-                         building_pars_f%fill )                                &
-                       surf_usm_v(l)%lambda_surf(m) =                          &
-                                    building_pars_f%pars_xy(ind_lambda_surf,j,i)
-                    
               ENDDO
            ENDDO
@@ -5097,4 +5035,46 @@
         
         CALL usm_init_material_model()
+
+!--     init skin layer properties
+!--     (can be done after initialization of wall layers)
+
+        DO  m = 1, surf_usm_h%ns
+            i = surf_usm_h%i(m)
+            j = surf_usm_h%j(m)
+
+            surf_usm_h%c_surface(m)           = surf_usm_h%rho_c_wall(nzb_wall,m) *              &
+                                                surf_usm_h%dz_wall(nzb_wall,m) * 0.25_wp
+            surf_usm_h%lambda_surf(m)         = surf_usm_h%lambda_h(nzb_wall,m) *                &
+                                                surf_usm_h%ddz_wall(nzb_wall,m) * 2.0_wp
+            surf_usm_h%c_surface_green(m)     = surf_usm_h%rho_c_wall(nzb_wall,m) *              &
+                                                surf_usm_h%dz_wall(nzb_wall,m) * 0.25_wp
+            surf_usm_h%lambda_surf_green(m)   = surf_usm_h%lambda_h_green(nzb_wall,m) *          &
+                                                surf_usm_h%ddz_green(nzb_wall,m) * 2.0_wp
+            surf_usm_h%c_surface_window(m)    = surf_usm_h%rho_c_window(nzb_wall,m) *            &
+                                                surf_usm_h%dz_window(nzb_wall,m) * 0.25_wp
+            surf_usm_h%lambda_surf_window(m)  = surf_usm_h%lambda_h_window(nzb_wall,m) * &
+                                                surf_usm_h%ddz_window(nzb_wall,m) * 2.0_wp
+        ENDDO
+
+        DO  l = 0, 3
+            DO  m = 1, surf_usm_v(l)%ns
+                i = surf_usm_v(l)%i(m) + surf_usm_v(l)%ioff
+                j = surf_usm_v(l)%j(m) + surf_usm_v(l)%joff
+
+                surf_usm_v(l)%c_surface(m) = surf_usm_v(l)%rho_c_wall(nzb_wall,m) *              &
+                                             surf_usm_v(l)%dz_wall(nzb_wall,m) * 0.25_wp
+                surf_usm_v(l)%lambda_surf(m) = surf_usm_v(l)%lambda_h(nzb_wall,m) *                &
+                                               surf_usm_v(l)%ddz_wall(nzb_wall,m) * 2.0_wp
+                surf_usm_v(l)%c_surface_green(m) = surf_usm_v(l)%rho_c_green(nzb_wall,m) *              &
+                                                   surf_usm_v(l)%dz_green(nzb_wall,m) * 0.25_wp
+                surf_usm_v(l)%lambda_surf_green(m) = surf_usm_v(l)%lambda_h_green(nzb_wall,m) *           &
+                                                     surf_usm_v(l)%ddz_green(nzb_wall,m) * 2.0_wp
+                surf_usm_v(l)%c_surface_window(m)    = surf_usm_v(l)%rho_c_window(nzb_wall,m) *            &
+                                                       surf_usm_v(l)%dz_window(nzb_wall,m) * 0.25_wp
+                surf_usm_v(l)%lambda_surf_window(m)  = surf_usm_v(l)%lambda_h_window(nzb_wall,m) * &
+                                                       surf_usm_v(l)%ddz_window(nzb_wall,m) * 2.0_wp
+            ENDDO
+        ENDDO
+
 !        
 !--     init anthropogenic sources of heat
@@ -7839,7 +7819,7 @@
         DO  m = 1, surf_usm_h%ns
 !
-!--       During spinup set green and window fraction to zero and restore
-!--       at the end of the loop.
-!--       Note, this is a temporary fix and need to be removed later.  
+!--        During spinup set green and window fraction to zero and restore
+!--        at the end of the loop.
+!--        Note, this is a temporary fix and need to be removed later.
            IF ( during_spinup )  THEN
               frac_win   = surf_usm_h%frac(m,ind_wat_win)
@@ -7960,5 +7940,5 @@
        
 !
-!--          Calculate water vapour pressure at saturation 
+!--           Calculate water vapour pressure at saturation
               e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surf_green_h(m) &
                             - 273.16_wp ) / ( t_surf_green_h(m) - 35.86_wp ) )
@@ -8027,7 +8007,7 @@
                                      surf_usm_h%rad_lw_out(m)
 !
-!--     numerator of the prognostic equation
-!--     Todo: Adjust to tile approach. So far, emissivity for wall (element 0)
-!--     is used
+!--        numerator of the prognostic equation
+!--        Todo: Adjust to tile approach. So far, emissivity for wall (element 0)
+!--        is used
            coef_1 = surf_usm_h%rad_net_l(m) +                                  & 
                  ( 3.0_wp + 1.0_wp ) * surf_usm_h%emissivity(m,ind_veg_wall) * &
@@ -8055,5 +8035,5 @@
                                        f_shf_green * surf_usm_h%pt1(m) +       &
                                        lambda_surface_green * t_green_h(nzb_wall,m)
-          ENDIF
+           ENDIF
 !
 !--        denominator of the prognostic equation
@@ -8320,5 +8300,5 @@
               ENDIF
 !
-!--          Get indices of respective grid point
+!--           Get indices of respective grid point
               i = surf_usm_v(l)%i(m)
               j = surf_usm_v(l)%j(m)
@@ -8326,13 +8306,11 @@
  
 !
-!--          TODO - how to calculate lambda_surface for horizontal (??? do you mean verical ???) surfaces
-!--          (lambda_surface is set according to stratification in land surface model).
-!--          Please note, for vertical surfaces no ol is defined, since 
-!--          stratification is not considered in this case.
+!--           Please note, for vertical surfaces no ol is defined, since
+!--           stratification is not considered in this case.
               lambda_surface = surf_usm_v(l)%lambda_surf(m)
               lambda_surface_window = surf_usm_v(l)%lambda_surf_window(m)
               lambda_surface_green = surf_usm_v(l)%lambda_surf_green(m)
  
-!            pt1  = pt(k,j,i)
+!             pt1  = pt(k,j,i)
               IF ( humidity )  THEN
                  qv1 = q(k,j,i)
@@ -8351,11 +8329,11 @@
               ENDIF
  
-!--          Calculation of r_a for vertical surfaces
+!--           Calculation of r_a for vertical surfaces
 !--
-!--          heat transfer coefficient for forced convection along vertical walls
-!--          follows formulation in TUF3d model (Krayenhoff & Voogt, 2006)
+!--           heat transfer coefficient for forced convection along vertical walls
+!--           follows formulation in TUF3d model (Krayenhoff & Voogt, 2006)
 !--            
-!--          H = httc (Tsfc - Tair)
-!--          httc = rw * (11.8 + 4.2 * Ueff) - 4.0
+!--           H = httc (Tsfc - Tair)
+!--           httc = rw * (11.8 + 4.2 * Ueff) - 4.0
 !--            
 !--                rw: wall patch roughness relative to 1.0 for concrete
@@ -8369,10 +8347,10 @@
 !--                zH: height of the convective layer
 !--                wstar = (g/Tcan*Qs*zH)**(1./3.)
-!--          Effective velocity components must always 
-!--          be defined at scalar grid point. The wall normal component is 
-!--          obtained by simple linear interpolation. ( An alternative would
-!--          be an logarithmic interpolation. )
-!--          Parameter roughness_concrete (default value = 0.001) is used
-!--          to calculation of roughness relative to concrete
+!--           Effective velocity components must always
+!--           be defined at scalar grid point. The wall normal component is
+!--           obtained by simple linear interpolation. ( An alternative would
+!--           be an logarithmic interpolation. )
+!--           Parameter roughness_concrete (default value = 0.001) is used
+!--           to calculation of roughness relative to concrete
               surf_usm_v(l)%r_a(m) = rho_cp / ( surf_usm_v(l)%z0(m) /           &
                          roughness_concrete * ( 11.8_wp + 4.2_wp *              &
@@ -8383,5 +8361,5 @@
                             )  - 4.0_wp  ) 
 !
-!--          Limit aerodynamic resistance
+!--           Limit aerodynamic resistance
               IF ( surf_usm_v(l)%r_a(m) < 1.0_wp )  surf_usm_v(l)%r_a(m) = 1.0_wp   
               
@@ -8540,4 +8518,5 @@
               t_surf_green_v_p(l)%t(m) = t_surf_green_v_p(l)%t(m) + dt_3d * tsc(3) *       &
                                  surf_usm_v(l)%tt_surface_green_m(m)
+
 !
 !--           Store surface temperature. Further, in case humidity is used
