diff --git a/clima/cython/ClimaRadtranWrk.pyx b/clima/cython/ClimaRadtranWrk.pyx index 950e672..f843b88 100644 --- a/clima/cython/ClimaRadtranWrk.pyx +++ b/clima/cython/ClimaRadtranWrk.pyx @@ -53,6 +53,17 @@ cdef class ClimaRadtranWrk: rwrk_pxd.climaradtranwrk_fdn_n_get(self._ptr, &dim1, 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, arr.data) + return arr + property tau_band: "ndarray[double,ndim=2], size (nz,nw). Band optical thickness." def __get__(self): diff --git a/clima/cython/ClimaRadtranWrk_pxd.pxd b/clima/cython/ClimaRadtranWrk_pxd.pxd index bca72ab..0ccfd3c 100644 --- a/clima/cython/ClimaRadtranWrk_pxd.pxd +++ b/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) \ No newline at end of file diff --git a/clima/fortran/ClimaRadtranWrk.f90 b/clima/fortran/ClimaRadtranWrk.f90 index 286aa90..086deff 100644 --- a/clima/fortran/ClimaRadtranWrk.f90 +++ b/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 diff --git a/src/radtran/clima_radtran.f90 b/src/radtran/clima_radtran.f90 index d409dea..fe59cff 100644 --- a/src/radtran/clima_radtran.f90 +++ b/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 diff --git a/src/radtran/clima_radtran_radiate.f90 b/src/radtran/clima_radtran_radiate.f90 index 172f54d..5c4d9c9 100644 --- a/src/radtran/clima_radtran_radiate.f90 +++ b/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,6 +55,7 @@ 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 @@ -62,7 +63,7 @@ function radiate(op, & 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 diff --git a/src/radtran/clima_radtran_types.f90 b/src/radtran/clima_radtran_types.f90 index 03028fc..32c250e 100644 --- a/src/radtran/clima_radtran_types.f90 +++ b/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(:) diff --git a/src/radtran/clima_radtran_types_create.f90 b/src/radtran/clima_radtran_types_create.f90 index 747b260..d2c8ed0 100644 --- a/src/radtran/clima_radtran_types_create.f90 +++ b/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))