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

Timestamp:

Apr 10, 2020 9:49:43 AM (5 years ago)

Author:

pavelkrc

Message:

Merge brach resler

Location:

palm/trunk

Files:

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

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -28 +28 @@
 ! -----------------
 ! $Id$
+! Avoid unstable direct normal radiation near horizon
+! 
+! 4481 2020-03-31 18:55:54Z maronga
 ! use statement for exchange horiz added
 ! 
@@ -454 +457 @@
                 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 +5810 @@
      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 +6024 @@
         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 +6044 @@
            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 +6052 @@
            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 +6095 @@
                    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 +7877 @@
 
 !--
-!--     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