Saved mean intensity during radiative transfer

clima/cython/ClimaRadtranWrk.pyx

@@ -53,6 +53,17 @@ cdef class ClimaRadtranWrk:
       rwrk_pxd.climaradtranwrk_fdn_n_get(self._ptr, &dim1, <double *>arr.data)
       return arr
 
+  property amean:
+    """ndarray[double,ndim=2], size (nz+1,nw). Mean intensity in photons/cm^2/s.
+    Only used in solar radiative transfer.
+    """
+    def __get__(self):
+      cdef int dim1, dim2
+      rwrk_pxd.climaradtranwrk_amean_get_size(self._ptr, &dim1, &dim2)
+      cdef ndarray arr = np.empty((dim1, dim2), np.double, order="F")
+      rwrk_pxd.climaradtranwrk_amean_get(self._ptr, &dim1, &dim2, <double *>arr.data)
+      return arr
+
   property tau_band:
     "ndarray[double,ndim=2], size (nz,nw). Band optical thickness."
     def __get__(self):

clima/cython/ClimaRadtranWrk_pxd.pxd

@@ -17,5 +17,8 @@ cdef extern void climaradtranwrk_fup_n_get(ClimaRadtranWrk *ptr, int *dim1, doub
 cdef extern void climaradtranwrk_fdn_n_get_size(ClimaRadtranWrk *ptr, int *dim1)
 cdef extern void climaradtranwrk_fdn_n_get(ClimaRadtranWrk *ptr, int *dim1, double *arr)
 
+cdef extern void climaradtranwrk_amean_get_size(ClimaRadtranWrk *ptr, int *dim1, int *dim2)
+cdef extern void climaradtranwrk_amean_get(ClimaRadtranWrk *ptr, int *dim1, int *dim2, double *arr)
+
 cdef extern void climaradtranwrk_tau_band_get_size(ClimaRadtranWrk *ptr, int *dim1, int *dim2)
 cdef extern void climaradtranwrk_tau_band_get(ClimaRadtranWrk *ptr, int *dim1, int *dim2, double *arr)

clima/fortran/ClimaRadtranWrk.f90

@@ -82,6 +82,26 @@ subroutine climaradtranwrk_fdn_n_get(ptr, dim1, arr) bind(c)
   arr = w%fdn_n
 end subroutine
 
+subroutine climaradtranwrk_amean_get_size(ptr, dim1, dim2) bind(c)
+  use clima_radtran, only: ClimaRadtranWrk
+  type(c_ptr), value, intent(in) :: ptr
+  integer(c_int), intent(out) :: dim1, dim2
+  type(ClimaRadtranWrk), pointer :: w
+  call c_f_pointer(ptr, w)
+  dim1 = size(w%amean,1)
+  dim2 = size(w%amean,2)
+end subroutine
+
+subroutine climaradtranwrk_amean_get(ptr, dim1, dim2, arr) bind(c)
+  use clima_radtran, only: ClimaRadtranWrk
+  type(c_ptr), value, intent(in) :: ptr
+  integer(c_int), intent(in) :: dim1, dim2
+  real(c_double), intent(out) :: arr(dim1,dim2)
+  type(ClimaRadtranWrk), pointer :: w
+  call c_f_pointer(ptr, w)
+  arr = w%amean
+end subroutine
+
 subroutine climaradtranwrk_tau_band_get_size(ptr, dim1, dim2) bind(c)
   use clima_radtran, only: ClimaRadtranWrk
   type(c_ptr), value, intent(in) :: ptr

src/radtran/clima_radtran.f90

@@ -19,6 +19,9 @@ module clima_radtran
     !> (nz+1) mW/m2 at the edges of the vertical grid 
     !> (integral of fup_a and fdn_a over wavelength grid)
     real(dp), allocatable :: fup_n(:), fdn_n(:)
+    !> (nz+1,nw) Mean intensity in photons/cm^2/s. Only used
+    !> in solar radiative transfer.
+    real(dp), allocatable :: amean(:,:)
     !> Band optical thickness (nz,nw)
     real(dp), allocatable :: tau_band(:,:)
 
@@ -183,6 +186,8 @@ function create_Radtran_2(species_names, particle_names, s, star_f, &
     allocate(rad%wrk_ir%fdn_a(nz+1, rad%ir%nw))
     allocate(rad%wrk_ir%fup_n(nz+1))
     allocate(rad%wrk_ir%fdn_n(nz+1))
+    allocate(rad%wrk_ir%amean(nz+1, rad%ir%nw))
+    rad%wrk_ir%amean = 0.0_dp ! not used
     allocate(rad%wrk_ir%tau_band(nz,rad%ir%nw))
 
     ! Solar work arrays
@@ -192,6 +197,7 @@ function create_Radtran_2(species_names, particle_names, s, star_f, &
     allocate(rad%wrk_sol%fdn_a(nz+1, rad%sol%nw))
     allocate(rad%wrk_sol%fup_n(nz+1))
     allocate(rad%wrk_sol%fdn_n(nz+1))
+    allocate(rad%wrk_sol%amean(nz+1, rad%sol%nw))
     allocate(rad%wrk_sol%tau_band(nz,rad%sol%nw))
 
     ! total flux
@@ -238,7 +244,7 @@ subroutine Radtran_radiate(self, T_surface, T, P, densities, dz, pdensities, rad
                    P, T_surface, T, densities, dz, &
                    pdensities, radii, &
                    wrk_ir%rx, wrk_ir%rz, &
-                   wrk_ir%fup_a, wrk_ir%fdn_a, wrk_ir%fup_n, wrk_ir%fdn_n, wrk_ir%tau_band)
+                   wrk_ir%fup_a, wrk_ir%fdn_a, wrk_ir%fup_n, wrk_ir%fdn_n, wrk_ir%amean, wrk_ir%tau_band)
     if (ierr /= 0) then
       err = 'Input particle radii are outside the data range.'
       return
@@ -252,7 +258,7 @@ subroutine Radtran_radiate(self, T_surface, T, P, densities, dz, pdensities, rad
                    P, T_surface, T, densities, dz, &
                    pdensities, radii, &
                    wrk_sol%rx, wrk_sol%rz, &
-                   wrk_sol%fup_a, wrk_sol%fdn_a, wrk_sol%fup_n, wrk_sol%fdn_n, wrk_sol%tau_band)
+                   wrk_sol%fup_a, wrk_sol%fdn_a, wrk_sol%fup_n, wrk_sol%fdn_n, wrk_sol%amean, wrk_sol%tau_band)
     if (ierr /= 0) then
       err = 'Input particle radii are outside the data range.'
       return

src/radtran/clima_radtran_radiate.f90

@@ -18,13 +18,13 @@ function radiate(op, &
                    P, T_surface, T, densities, dz, &
                    pdensities, radii, &
                    rw, rz, &
-                   fup_a, fdn_a, fup_n, fdn_n, tau_band) result(ierr)
+                   fup_a, fdn_a, fup_n, fdn_n, amean, tau_band) result(ierr)
     use clima_radtran_types, only: RadiateXSWrk, RadiateZWrk, Ksettings
     use clima_radtran_types, only: OpticalProperties, Kcoefficients
     use clima_radtran_types, only: FarUVOpticalProperties, SolarOpticalProperties, IROpticalProperties
     use clima_radtran_types, only: k_RandomOverlap, k_RandomOverlapResortRebin, k_AdaptiveEquivalentExtinction
     use clima_eqns, only: planck_fcn, ten2power
-    use clima_const, only: k_boltz, pi
+    use clima_const, only: k_boltz, pi, plank, c_light
 
     type(OpticalProperties), target, intent(inout) :: op !! Optical properties
 
@@ -55,14 +55,15 @@ function radiate(op, &
                                                     !! at the edges of the vertical grid
     real(dp), intent(out) :: fup_n(:), fdn_n(:) !! (nz+1) mW/m2 at the edges of the vertical grid 
                                                 !! (integral of fup_a and fdn_a over wavelength grid)
+    real(dp), intent(out) :: amean(:,:) !! (nz+1,nw) Mean intensity in photons/cm^2/s (Only for Solar)
     real(dp), intent(out) :: tau_band(:,:) !! (nz,nw) The optical depth of each layer
     integer :: ierr !! if ierr /= 0 on return, then there was an error
 
     type(Ksettings), pointer :: kset
     integer :: nz, ng
     integer :: i, j, k, l, n, jj
     integer :: ierr_0
-    real(dp) :: avg_freq
+    real(dp) :: avg_freq, avg_wavl
 
     ! array of indexes for recursive
     ! correlated-k
@@ -273,6 +274,7 @@ function radiate(op, &
 
       rz%fup2 = 0.0_dp
       rz%fdn2 = 0.0_dp
+      rz%amean2 = 0.0_dp
 
       if (op%nk /= 0) then
 
@@ -284,10 +286,15 @@ function radiate(op, &
           do j = 1,size(zenith_u)
             rz%fup1 = 0.0_dp
             rz%fdn1 = 0.0_dp
+            rz%amean1 = 0.0_dp
             rz%tau_band = 0.0_dp
             call k_loops(op, zenith_u(j), albedo, emissivity, cols, rw%ks, rz, iks, 1)
             rz%fup2(:) = rz%fup2(:) + rz%fup1(:)*zenith_weights(j)
             rz%fdn2(:) = rz%fdn2(:) + rz%fdn1(:)*zenith_weights(j)
+            if (op%op_type == FarUVOpticalProperties .or. &
+                op%op_type == SolarOpticalProperties) then
+              rz%amean2 = rz%amean2 + rz%amean1*zenith_weights(j)
+            endif
           enddo
         elseif (kset%k_method == k_RandomOverlapResortRebin) then
           ! Random Overlap with Resorting and Rebinning.
@@ -310,6 +317,13 @@ function radiate(op, &
         fup_a(i,l) = rz%fup2(n)
         fdn_a(i,l) = rz%fdn2(n)
       enddo
+      if (op%op_type == FarUVOpticalProperties .or. &
+          op%op_type == SolarOpticalProperties) then
+        do i = 1,nz+1
+          n = nz+2-i
+          amean(i,l) = rz%amean2(n)
+        enddo
+      endif
       do i = 1,nz
         n = nz+1-i
         tau_band(i,l) = rz%tau_band(n) ! band optical thickness
@@ -330,6 +344,14 @@ function radiate(op, &
       do l = 1,op%nw
         fup_a(:,l) = fup_a(:,l)*photons_sol(l)*diurnal_fac
         fdn_a(:,l) = fdn_a(:,l)*photons_sol(l)*diurnal_fac
+        amean(:,l) = amean(:,l)*photons_sol(l)*diurnal_fac
+
+        ! Convert from from mW/m^2/Hz to mW/m^2/nm
+        avg_freq = 0.5_dp*(op%freq(l) + op%freq(l+1))
+        avg_wavl = 1.0e9_dp*c_light/avg_freq ! nm
+        amean(:,l) = amean(:,l)*(avg_freq/avg_wavl)
+        ! Convert from mW/m^2/nm to photons/cm^2/s
+        amean(:,l) = amean(:,l)*(avg_wavl/(plank*c_light*1.0e16_dp))*(op%wavl(l+1) - op%wavl(l))
       enddo
     endif
 
@@ -390,6 +412,7 @@ subroutine k_aee(op, kset, zenith_u, zenith_weights, surface_albedo, surface_emi
     do ii = 1,size(zenith_u)
     rz%fup1 = 0.0_dp
     rz%fdn1 = 0.0_dp
+    rz%amean1 = 0.0_dp
     rz%tau_band = 0.0_dp
 
     do i = 1,op%k(1)%ngauss
@@ -426,6 +449,10 @@ subroutine k_aee(op, kset, zenith_u, zenith_weights, surface_albedo, surface_emi
       ! k-coeff weights (kset%wbin(i))
       rz%fup1(:) = rz%fup1(:) + rz%fup(:)*op%k(1)%weights(i)
       rz%fdn1(:) = rz%fdn1(:) + rz%fdn(:)*op%k(1)%weights(i)
+      if (op%op_type == FarUVOpticalProperties .or. &
+          op%op_type == SolarOpticalProperties) then
+        rz%amean1 = rz%amean1 + rz%amean(:)*op%k(1)%weights(i)
+      endif
 
       ! Save the optical thickness in the band.
       rz%tau_band(:) = rz%tau_band(:) + rz%tau(:)*op%k(1)%weights(i)
@@ -434,6 +461,10 @@ subroutine k_aee(op, kset, zenith_u, zenith_weights, surface_albedo, surface_emi
 
     rz%fup2 = rz%fup2 + rz%fup1*zenith_weights(ii)
     rz%fdn2 = rz%fdn2 + rz%fdn1*zenith_weights(ii)
+    if (op%op_type == FarUVOpticalProperties .or. &
+          op%op_type == SolarOpticalProperties) then
+      rz%amean2 = rz%amean2 + rz%amean1*zenith_weights(ii)
+    endif
     enddo
 
   end subroutine
@@ -521,6 +552,7 @@ subroutine k_rorr(op, kset, zenith_u, zenith_weights, surface_albedo, surface_em
     do ii = 1,size(zenith_u)
     rz%fup1 = 0.0_dp
     rz%fdn1 = 0.0_dp
+    rz%amean1 = 0.0_dp
     rz%tau_band = 0.0_dp
 
     ! loop over mixed k-coeffs
@@ -554,6 +586,10 @@ subroutine k_rorr(op, kset, zenith_u, zenith_weights, surface_albedo, surface_em
       ! k-coeff weights (kset%wbin(i))
       rz%fup1 = rz%fup1 + rz%fup*kset%wbin(i)
       rz%fdn1 = rz%fdn1 + rz%fdn*kset%wbin(i)
+      if (op%op_type == FarUVOpticalProperties .or. &
+          op%op_type == SolarOpticalProperties) then
+        rz%amean1 = rz%amean1 + rz%amean*kset%wbin(i)
+      endif
 
       ! Save the optical thickness in the band.
       rz%tau_band(:) = rz%tau_band(:) + rz%tau(:)*kset%wbin(i)
@@ -562,6 +598,10 @@ subroutine k_rorr(op, kset, zenith_u, zenith_weights, surface_albedo, surface_em
 
     rz%fup2 = rz%fup2 + rz%fup1*zenith_weights(ii)
     rz%fdn2 = rz%fdn2 + rz%fdn1*zenith_weights(ii)
+    if (op%op_type == FarUVOpticalProperties .or. &
+        op%op_type == SolarOpticalProperties) then
+      rz%amean2 = rz%amean2 + rz%amean1*zenith_weights(ii)
+    endif
     enddo
 
   end subroutine
@@ -629,6 +669,10 @@ recursive subroutine k_loops(op, u0, surface_albedo, surface_emissivity, cols, k
 
         rz%fup1 = rz%fup1 + rz%fup*gauss_weight
         rz%fdn1 = rz%fdn1 + rz%fdn*gauss_weight
+        if (op%op_type == FarUVOpticalProperties .or. &
+          op%op_type == SolarOpticalProperties) then
+          rz%amean1 = rz%amean1 + rz%amean*gauss_weight
+        endif
 
         ! Save the optical thickness in the band.
         rz%tau_band(:) = rz%tau_band(:) + rz%tau(:)*gauss_weight

src/radtran/clima_radtran_types.f90

@@ -213,6 +213,8 @@ module function create_RadiateXSWrk(op, nz) result(rw)
     real(dp), allocatable :: tau_band(:) !! band optical thickness.
 
     real(dp), allocatable :: amean(:)
+    real(dp), allocatable :: amean1(:)
+    real(dp), allocatable :: amean2(:)
     real(dp), allocatable :: fup1(:)
     real(dp), allocatable :: fdn1(:)
     real(dp), allocatable :: fup2(:)

src/radtran/clima_radtran_types_create.f90

@@ -146,6 +146,8 @@ module function create_RadiateZWrk(nz, npart) result(rz)
     allocate(rz%tau_band(nz))
 
     allocate(rz%amean(nz+1))
+    allocate(rz%amean1(nz+1))
+    allocate(rz%amean2(nz+1))
     allocate(rz%fup1(nz+1))
     allocate(rz%fdn1(nz+1))
     allocate(rz%fup2(nz+1))