Index: /palm/trunk/SOURCE/chem_photolysis_mod.f90
===================================================================
--- /palm/trunk/SOURCE/chem_photolysis_mod.f90	(revision 3823)
+++ /palm/trunk/SOURCE/chem_photolysis_mod.f90	(revision 3824)
@@ -234,5 +234,5 @@
 
        USE radiation_model_mod,                                                   &
-           ONLY:  calc_zenith, zenith
+           ONLY:  calc_zenith, cos_zenith
 
        IMPLICIT NONE
@@ -247,6 +247,6 @@
        CALL calc_zenith
 
-       IF ( zenith(0) > 0.0_wp ) THEN
-          coszi = 1. / zenith(0)
+       IF ( cos_zenith > 0.0_wp ) THEN
+          coszi = 1. / cos_zenith
 
           DO iphot = 1,nphot
@@ -254,5 +254,5 @@
                 IF ( TRIM( names_s(iav) ) == TRIM( phot_names(iphot) ) ) then
                    phot_frequen(iphot)%freq(nzb+1:nzt,:,:) =    &
-                                par_l(iav) * zenith(0)**par_m(iav) *  EXP( -par_n(iav) * coszi ) 
+                                par_l(iav) * cos_zenith**par_m(iav) *  EXP( -par_n(iav) * coszi ) 
                 ENDIF
              ENDDO
Index: /palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3823)
+++ /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3824)
@@ -2799,5 +2799,5 @@
 
     USE radiation_model_mod,                                                &
-         ONLY:  zenith
+         ONLY:  cos_zenith
 
 
@@ -3419,5 +3419,5 @@
                 !
                 !-- Get Rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
                      relh, lai, sai, nwet, lu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
                      r_aero_lu , Rb )
@@ -3566,5 +3566,5 @@
                 !
                 !-- Get Rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
                      relh, lai, sai, nwet, lup_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
                      r_aero_lu , Rb )
@@ -3713,5 +3713,5 @@
                 !
                 !-- Get Rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
                      relh, lai, sai, nwet, luw_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
                      r_aero_lu , Rb )
@@ -4109,5 +4109,5 @@
                 !
                 !-- Get Rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
                      relh, lai, sai, nwet, luu_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
                      r_aero_lu , Rb )
@@ -4249,5 +4249,5 @@
                 !
                 !-- Get Rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
                      relh, lai, sai, nwet, lug_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
                      r_aero_lu , Rb )
@@ -4395,5 +4395,5 @@
                 !
                 !-- Get Rc_tot
-                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, zenith(0), &
+                CALL drydepos_gas_depac( spc_names(lsp), day_of_year, latitude, ts, ustar_lu, glrad, cos_zenith, &
                      relh, lai, sai, nwet, lud_dep, 2, Rc_tot, ccomp_tot(lsp), hyp(nzb), catm, diffc, &
                      r_aero_lu , Rb )
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3823)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3824)
@@ -28,4 +28,8 @@
 ! -----------------
 ! $Id$
+! Change zenith(0:0) and others to scalar.
+! Code review.
+! 
+! 3814 2019-03-26 08:40:31Z pavelkrc
 ! Rename exported nzu, nzp and related variables due to name conflict
 ! 
@@ -747,7 +751,7 @@
 
 
-    REAL(wp), DIMENSION(0:0) ::  zenith,         & !< cosine of solar zenith angle
-                                 sun_dir_lat,    & !< solar directional vector in latitudes
-                                 sun_dir_lon       !< solar directional vector in longitudes
+    REAL(wp) ::  cos_zenith        !< cosine of solar zenith angle, also z-coordinate of solar unit vector
+    REAL(wp) ::  sun_dir_lat       !< y-coordinate of solar unit vector
+    REAL(wp) ::  sun_dir_lon       !< x-coordinate of solar unit vector
 
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_net_av       !< average of net radiation (rad_net) at surface
@@ -1294,7 +1298,7 @@
            rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
            rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,              &
-           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb, solar_constant, &
+           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, solar_constant,           &
            skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
-           zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,       &
+           cos_zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,   &
            idir, jdir, kdir, id, iz, iy, ix,                                   &
            iup_u, inorth_u, isouth_u, ieast_u, iwest_u,                        &
@@ -2834,5 +2838,5 @@
 !
 !--    Calculate sky transmissivity
-       sky_trans = 0.6_wp + 0.2_wp * zenith(0)
+       sky_trans = 0.6_wp + 0.2_wp * cos_zenith
 
 !
@@ -2909,5 +2913,5 @@
                 k = nz_urban_t
 
-                surf%rad_sw_in  = solar_constant * sky_trans * zenith(0)
+                surf%rad_sw_in  = solar_constant * sky_trans * cos_zenith
                 surf%rad_sw_out = albedo_urb * surf%rad_sw_in
                 
@@ -2933,5 +2937,5 @@
                    k = surf%k(m)
 
-                   surf%rad_sw_in(m) = solar_constant * sky_trans * zenith(0)
+                   surf%rad_sw_in(m) = solar_constant * sky_trans * cos_zenith
 
 !
@@ -3379,4 +3383,5 @@
 
        REAL(wp), DIMENSION(0:nzt+1) :: pt_av, q_av, ql_av
+       REAL(wp), DIMENSION(0:0)     :: zenith   !< to provide indexed array
 !
 !--    Just dummy arguments
@@ -3395,4 +3400,5 @@
 !--    Calculate current (cosine of) zenith angle and whether the sun is up
        CALL calc_zenith     
+       zenith(0) = cos_zenith
 !
 !--    Calculate surface albedo. In case average radiation is applied, 
@@ -4172,7 +4178,7 @@
 !
 !--    Calculate cosine of solar zenith angle
-       zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)   &
+       cos_zenith = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)   &
                                             * COS(hour_angle)
-       zenith(0) = MAX(0.0_wp,zenith(0))
+       cos_zenith = MAX(0.0_wp,cos_zenith)
 
 !
@@ -4182,9 +4188,9 @@
 !
 !--       Direction in longitudes equals to sin(solar_azimuth) * sin(zenith)
-          sun_dir_lon(0) = -SIN(hour_angle) * COS(declination)
+          sun_dir_lon = -SIN(hour_angle) * COS(declination)
 
 !
 !--       Direction in latitues equals to cos(solar_azimuth) * sin(zenith)
-          sun_dir_lat(0) = SIN(declination) * COS(lat) - COS(hour_angle) &
+          sun_dir_lat = SIN(declination) * COS(lat) - COS(hour_angle) &
                               * COS(declination) * SIN(lat)
        ENDIF
@@ -4192,5 +4198,5 @@
 !
 !--    Check if the sun is up (otheriwse shortwave calculations can be skipped)
-       IF ( zenith(0) > 0.0_wp )  THEN
+       IF ( cos_zenith > 0.0_wp )  THEN
           sun_up = .TRUE.
        ELSE
@@ -4227,21 +4233,21 @@
                  IF ( surf%albedo_type(ind_type,m) == 1 )  THEN
                     surf%rrtm_aldir(ind_type,m) = 0.026_wp /                    &
-                                                ( zenith(0)**1.7_wp + 0.065_wp )&
-                                     + 0.15_wp * ( zenith(0) - 0.1_wp )         &
-                                               * ( zenith(0) - 0.5_wp )         &
-                                               * ( zenith(0) - 1.0_wp )
+                                                ( cos_zenith**1.7_wp + 0.065_wp )&
+                                     + 0.15_wp * ( cos_zenith - 0.1_wp )         &
+                                               * ( cos_zenith - 0.5_wp )         &
+                                               * ( cos_zenith - 1.0_wp )
                     surf%rrtm_asdir(ind_type,m) = surf%rrtm_aldir(ind_type,m)
 !
 !--              Snow
                  ELSEIF ( surf%albedo_type(ind_type,m) == 16 )  THEN
-                    IF ( zenith(0) < 0.5_wp )  THEN
+                    IF ( cos_zenith < 0.5_wp )  THEN
                        surf%rrtm_aldir(ind_type,m) =                           &
                                  0.5_wp * ( 1.0_wp - surf%aldif(ind_type,m) )  &
                                         * ( 3.0_wp / ( 1.0_wp + 4.0_wp         &
-                                        * zenith(0) ) ) - 1.0_wp
+                                        * cos_zenith ) ) - 1.0_wp
                        surf%rrtm_asdir(ind_type,m) =                           &
                                  0.5_wp * ( 1.0_wp - surf%asdif(ind_type,m) )  &
                                         * ( 3.0_wp / ( 1.0_wp + 4.0_wp         &
-                                        * zenith(0) ) ) - 1.0_wp
+                                        * cos_zenith ) ) - 1.0_wp
 
                        surf%rrtm_aldir(ind_type,m) =                           &
@@ -4282,8 +4288,8 @@
                           surf%rrtm_aldir(ind_type,m) =                        &
                                 surf%aldif(ind_type,m) * 1.4_wp /              &
-                                           ( 1.0_wp + 0.8_wp * zenith(0) )
+                                           ( 1.0_wp + 0.8_wp * cos_zenith )
                           surf%rrtm_asdir(ind_type,m) =                        &
                                 surf%asdif(ind_type,m) * 1.4_wp /              &
-                                           ( 1.0_wp + 0.8_wp * zenith(0) )
+                                           ( 1.0_wp + 0.8_wp * cos_zenith )
 !
 !--                    Surface types with weak zenith dependence
@@ -4291,8 +4297,8 @@
                           surf%rrtm_aldir(ind_type,m) =                        &
                                 surf%aldif(ind_type,m) * 1.1_wp /              &
-                                           ( 1.0_wp + 0.2_wp * zenith(0) )
+                                           ( 1.0_wp + 0.2_wp * cos_zenith )
                           surf%rrtm_asdir(ind_type,m) =                        &
                                 surf%asdif(ind_type,m) * 1.1_wp /              &
-                                           ( 1.0_wp + 0.2_wp * zenith(0) )
+                                           ( 1.0_wp + 0.2_wp * cos_zenith )
 
                        CASE DEFAULT
@@ -5018,5 +5024,5 @@
      mrot(2, :) = (/ 0._wp,  COS(alpha), SIN(alpha) /)
      mrot(3, :) = (/ 0._wp, -SIN(alpha), COS(alpha) /)
-     sunorig = (/ zenith(0), sun_dir_lat, sun_dir_lon /)
+     sunorig = (/ cos_zenith, sun_dir_lat, sun_dir_lon /)
      sunorig = MATMUL(mrot, sunorig)
      DO d = 0, nsurf_type
@@ -5024,5 +5030,5 @@
      ENDDO
 
-     IF ( zenith(0) > 0 )  THEN
+     IF ( cos_zenith > 0 )  THEN
 !--      now we will "squash" the sunorig vector by grid box size in
 !--      each dimension, so that this new direction vector will allow us
@@ -5186,9 +5192,9 @@
 
      !-- direct radiation
-     IF ( zenith(0) > 0 )  THEN
+     IF ( cos_zenith > 0 )  THEN
         !--Identify solar direction vector (discretized number) 1)
         !--
-        j = FLOOR(ACOS(zenith(0)) / pi * raytrace_discrete_elevs)
-        i = MODULO(NINT(ATAN2(sun_dir_lon(0), sun_dir_lat(0))               &
+        j = FLOOR(ACOS(cos_zenith) / pi * raytrace_discrete_elevs)
+        i = MODULO(NINT(ATAN2(sun_dir_lon, sun_dir_lat)               &
                         / (2._wp*pi) * raytrace_discrete_azims-.5_wp, iwp), &
                    raytrace_discrete_azims)
@@ -5199,5 +5205,5 @@
            i = surfl(ix, isurf)
            surfinswdir(isurf) = rad_sw_in_dir(j,i) *                        &
-                costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / zenith(0)
+                costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / cos_zenith
         ENDDO
 !
@@ -5207,5 +5213,5 @@
            i = mrtbl(ix, imrt)
            mrtinsw(imrt) = mrtinsw(imrt) + mrtdsit(imrt, isd) * rad_sw_in_dir(j,i) &
-                                     / zenith(0) / 4._wp ! normal to sphere
+                                     / cos_zenith / 4._wp ! normal to sphere
         ENDDO
      ENDIF
@@ -5217,5 +5223,7 @@
         mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc)
      ENDDO
-
+!
+!--  Absorption in each local plant canopy grid box from the first atmospheric
+!--  pass of radiation
      IF ( npcbl > 0 )  THEN
 
@@ -5223,6 +5231,5 @@
          pcbinswdif(:) = 0._wp
          pcbinlw(:) = 0._wp
-!
-!--      pcsf first pass
+
          DO icsf = 1, ncsfl
              ipcgb = csfsurf(1, icsf)
@@ -5234,8 +5241,8 @@
              IF ( isurfsrc == -1 )  THEN
 !
-!--             Diffuse rad from sky.
+!--             Diffuse radiation from sky
                 pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i)
 !
-!--             Absorbed diffuse LW from sky minus emitted to sky
+!--             Absorbed diffuse LW radiation from sky minus emitted to sky
                 IF ( plant_lw_interact )  THEN
                    pcbinlw(ipcgb) = csf(1,icsf)                                  &
@@ -5244,6 +5251,6 @@
                 ENDIF
 !
-!--             Direct rad
-                IF ( zenith(0) > 0 )  THEN
+!--             Direct solar radiation
+                IF ( cos_zenith > 0 )  THEN
 !--                Estimate directed box absorption
                    pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i))
@@ -5312,5 +5319,6 @@
 
          surfoutsl = albedo_surf * surfins
-!--         for non-transparent surfaces, longwave albedo is 1 - emissivity
+!
+!--      for non-transparent surfaces, longwave albedo is 1 - emissivity
          surfoutll = (1._wp - emiss_surf) * surfinl
 
@@ -5318,5 +5326,5 @@
          CALL MPI_AllGatherv(surfoutsl, nsurfl, MPI_REAL, &
              surfouts, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
-         IF ( ierr /= 0 ) THEN
+         IF ( ierr /= 0 )  THEN
              WRITE(9,*) 'Error MPI_AllGatherv2:', ierr, SIZE(surfoutsl), nsurfl, &
                         SIZE(surfouts), nsurfs, surfstart
@@ -5326,5 +5334,5 @@
          CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, &
              surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
-         IF ( ierr /= 0 ) THEN
+         IF ( ierr /= 0 )  THEN
              WRITE(9,*) 'Error MPI_AllGatherv3:', ierr, SIZE(surfoutll), nsurfl, &
                         SIZE(surfoutl), nsurfs, surfstart
@@ -5336,10 +5344,10 @@
          surfoutl = surfoutll
 #endif
-
-!--         reset for next pass input
+!
+!--      Reset for the input from next reflective pass
          surfins = 0._wp
          surfinl = 0._wp
-
-!--         reflected radiation
+!
+!--      Reflected radiation
          DO isvf = 1, nsvfl
              isurf = svfsurf(1, isvf)
@@ -5825,5 +5833,5 @@
 !> This subroutine splits direct and diffusion dw radiation
 !> It sould not be called in case the radiation model already does it
-!> It follows <CITATION>
+!> It follows Boland, Ridley & Brown (2008)
 !------------------------------------------------------------------------------!
     SUBROUTINE calc_diffusion_radiation 
@@ -5855,8 +5863,8 @@
 !--     towards 0 while keeping full continuity.
 !--    
-        IF ( zenith(0) <= lowest_solarUp )  THEN
+        IF ( cos_zenith <= lowest_solarUp )  THEN
             corrected_solarUp = lowest_solarUp
         ELSE
-            corrected_solarUp = zenith(0)
+            corrected_solarUp = cos_zenith
         ENDIF
         
@@ -8225,11 +8233,11 @@
 !--      Update apparent solar position based on modified t_s_r_p
          CALL calc_zenith
-         IF ( zenith(0) > 0 )  THEN
+         IF ( cos_zenith > 0 )  THEN
 !--         
 !--         Identify solar direction vector (discretized number) 1)
-            i = MODULO(NINT(ATAN2(sun_dir_lon(0), sun_dir_lat(0))               &
+            i = MODULO(NINT(ATAN2(sun_dir_lon, sun_dir_lat)               &
                             / (2._wp*pi) * raytrace_discrete_azims-.5_wp, iwp), &
                        raytrace_discrete_azims)
-            j = FLOOR(ACOS(zenith(0)) / pi * raytrace_discrete_elevs)
+            j = FLOOR(ACOS(cos_zenith) / pi * raytrace_discrete_elevs)
             IF ( dsidir_rev(j, i) == -1 )  THEN
                ndsidir = ndsidir + 1
Index: /palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3823)
+++ /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3824)
@@ -28,4 +28,8 @@
 ! -----------------
 ! $Id$
+! Remove unused imports
+!
+! 
+! 3814 2019-03-26 08:40:31Z pavelkrc
 ! unused subroutine commented out
 ! 
@@ -438,5 +442,5 @@
 
     USE basic_constants_and_equations_mod,                                     &
-        ONLY:  c_p, g, kappa, pi, r_d, rho_l, l_v
+        ONLY:  c_p, g, kappa, pi, r_d, rho_l, l_v, sigma_sb
 
     USE control_parameters,                                                    &
@@ -473,7 +477,6 @@
         
     USE radiation_model_mod,                                                   &
-        ONLY:  albedo_type, radiation_interaction, calc_zenith, zenith,        &
+        ONLY:  albedo_type, radiation_interaction,                             &
                radiation, rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out,        &
-               sigma_sb, sun_direction, sun_dir_lat, sun_dir_lon,              &
                force_radiation_call, iup_u, inorth_u, isouth_u, ieast_u,       &
                iwest_u, iup_l, inorth_l, isouth_l, ieast_l, iwest_l, id,       &
