Skip to content

Commit

Permalink
added function to make string of opacities
Browse files Browse the repository at this point in the history
  • Loading branch information
Nicholaswogan committed Oct 27, 2023
1 parent 3751150 commit 87d6403
Show file tree
Hide file tree
Showing 4 changed files with 212 additions and 6 deletions.
59 changes: 59 additions & 0 deletions src/radtran/clima_radtran.f90
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,8 @@ module clima_radtran
procedure :: TOA_fluxes => Radtran_TOA_fluxes
procedure :: skin_temperature => Radtran_skin_temperature
procedure :: equilibrium_temperature => Radtran_equilibrium_temperature
procedure :: opacities2yaml => Radtran_opacities2yaml
procedure :: print_opacities => Radtran_print_opacities
end type

interface Radtran
Expand Down Expand Up @@ -452,6 +454,63 @@ function Radtran_equilibrium_temperature(self, bond_albedo) result(T_eq)
T_eq = equilibrium_temperature(stellar_radiation, bond_albedo)
end function

function Radtran_opacities2yaml(self) result(out)
class(Radtran), target, intent(inout) :: self
character(:), allocatable :: out

character(:), allocatable :: line

out = ''

line = ''
line = line//'optical-properties:'
out = out//line

out = out//'\n'
line = ' '
line = line//'ir:'
out = out//line

out = out//'\n'
out = out//self%ir%opacities2yaml()

out = out//'\n'
line = ' '
line = line//'solar:'
out = out//line

out = out//'\n'
out = out//self%sol%opacities2yaml()

end function

subroutine Radtran_print_opacities(self)
use iso_fortran_env, only: output_unit
class(Radtran), target, intent(inout) :: self
character(:), allocatable :: out

integer :: i
logical :: skip_next

out = self%opacities2yaml()

skip_next = .false.
do i = 1,len(out)-1
if (out(i:i+1) == '\n') then
write(output_unit,'(a)')
skip_next = .true.
else
if (.not. skip_next) then
write(output_unit,'(a)',advance='no') out(i:i)
else
skip_next = .false.
endif
endif
enddo
write(output_unit,'(a)') out(len(out):len(out))

end subroutine

!!!!!!!!!!!!!!!!!
!!! Utilities !!!
!!!!!!!!!!!!!!!!!
Expand Down
135 changes: 132 additions & 3 deletions src/radtran/clima_radtran_types.f90
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@

module clima_radtran_types
use iso_c_binding
use clima_const, only: dp
use clima_const, only: dp, s_str_len
use linear_interpolation_module, only: linear_interp_1d, linear_interp_2d
implicit none
public
Expand Down Expand Up @@ -51,6 +51,7 @@ module clima_radtran_types
end type

type :: ParticleXsection
character(:), allocatable :: dat_name
integer :: p_ind
integer :: nrad
real(dp), allocatable :: radii(:) ! cm
Expand All @@ -61,6 +62,7 @@ module clima_radtran_types
end type

type :: WaterContinuum
character(:), allocatable :: model
integer :: LH2O ! index of H2O
integer :: ntemp ! number of temperatures
real(dp), allocatable :: temp(:) ! (ntemp) Kelvin
Expand All @@ -76,8 +78,10 @@ module clima_radtran_types
end enum

type :: Ksettings
integer, allocatable :: new_num_k_bins
! approach to combining k-distributions
integer :: k_method
character(:), allocatable :: k_method_name ! name
integer :: k_method ! enum (see above)
! if k_method == k_RandomOverlapResortRebin
! then these are the weights we are re-binning to.
integer :: nbin = -1
Expand All @@ -94,6 +98,10 @@ module clima_radtran_types
integer :: nw
real(dp), allocatable :: wavl(:)
real(dp), allocatable :: freq(:)

! Copy of species and particle names
character(s_str_len), allocatable :: species_names(:)
character(s_str_len), allocatable :: particle_names(:)

! K-distributions (e.g. H2O)
integer :: ngauss_max = -1
Expand All @@ -120,7 +128,9 @@ module clima_radtran_types
! Water Continuum absorption. Can only be a thing if H2O is a species
! and if there are other gases in the atmosphere.
type(WaterContinuum), allocatable :: cont


contains
procedure :: opacities2yaml => OpticalProperties_opacities2yaml
end type

interface
Expand Down Expand Up @@ -234,5 +244,124 @@ module subroutine read_stellar_flux(star_file, nw, wavl, photon_scale_factor, ph
character(:), allocatable, intent(out) :: err
end subroutine
end interface

contains

function OpticalProperties_opacities2yaml(self) result(out)
use clima_const, only:
class(OpticalProperties), intent(inout) :: self
character(:), allocatable :: out

character(:), allocatable :: line
character(s_str_len) :: tmp_str
integer :: i

out = ''

line = ' '
line = line//'k-method: '//self%kset%k_method_name
out = out//line

if (self%kset%k_method == k_RandomOverlapResortRebin) then
out = out//'\n'
write(tmp_str,'(i0)') self%kset%nbin
line = ' '
line = line//'number-of-bins: '//trim(tmp_str)
out = out//line
endif

if (allocated(self%kset%new_num_k_bins)) then
out = out//'\n'
write(tmp_str,'(i0)') self%kset%new_num_k_bins
line = ' '
line = line//'new-num-k-bins: '//trim(tmp_str)
out = out//line
endif

out = out//'\n'
line = ' '
line = line//'opacities:'
out = out//line

if (allocated(self%k)) then
out = out//'\n'
line = ' '
line = line//'k-distributions: ['
do i = 1,self%nk
line = line//trim(self%species_names(self%k(i)%sp_ind))
if (i /= self%nk) then
line = line//', '
endif
enddo
line = line//']'
out = out//line
endif

if (allocated(self%cia)) then
out = out//'\n'
line = ' '
line = line//'CIA: ['
do i = 1,self%ncia
line = line//trim(self%species_names(self%cia(i)%sp_ind(1)))// &
'-'//trim(self%species_names(self%cia(i)%sp_ind(2)))
if (i /= self%ncia) then
line = line//', '
endif
enddo
line = line//']'
out = out//line
endif

if (allocated(self%ray)) then
out = out//'\n'
line = ' '
line = line//'rayleigh: ['
do i = 1,self%nray
line = line//trim(self%species_names(self%ray(i)%sp_ind(1)))
if (i /= self%nray) then
line = line//', '
endif
enddo
line = line//']'
out = out//line
endif

if (allocated(self%pxs)) then
out = out//'\n'
line = ' '
line = line//'photolysis-xs: ['
do i = 1,self%npxs
line = line//trim(self%species_names(self%pxs(i)%sp_ind(1)))
if (i /= self%npxs) then
line = line//', '
endif
enddo
line = line//']'
out = out//line
endif

if (allocated(self%cont)) then
out = out//'\n'
line = ' '
line = line//'water-continuum: '//self%cont%model
out = out//line
endif

if (allocated(self%part)) then
out = out//'\n'
line = ' '
line = line//'particle-xs: ['
do i = 1,self%npart
line = line//'{name: '//trim(self%particle_names(self%part(i)%p_ind))
line = line//', data: '//trim(self%part(i)%dat_name)//'}'
if (i /= self%npart) then
line = line//', '
endif
enddo
line = line//']'
out = out//line
endif

end function

end module
22 changes: 19 additions & 3 deletions src/radtran/clima_radtran_types_create.f90
Original file line number Diff line number Diff line change
Expand Up @@ -164,6 +164,13 @@ function create_Ksettings(sop) result(kset)
type(Ksettings) :: kset

real(dp), allocatable :: tmp(:) ! dummy variable.

if (allocated(sop%new_num_k_bins)) then
allocate(kset%new_num_k_bins)
kset%new_num_k_bins = sop%new_num_k_bins
endif

kset%k_method_name = sop%k_method

! method for mixing k-distributions.
if (sop%k_method == "RandomOverlapResortRebin") then
Expand Down Expand Up @@ -218,7 +225,12 @@ module function create_OpticalProperties(datadir, optype, species_names, &
op%nw = size(op%wavl) - 1
allocate(op%freq(size(op%wavl)))
op%freq = c_light/(op%wavl*1.0e-9_dp)


allocate(op%species_names(size(species_names)))
op%species_names = species_names
allocate(op%particle_names(size(particle_names)))
op%particle_names = particle_names

!!!!!!!!!!!!!!!!!!!!!!!
!!! k-distributions !!!
!!!!!!!!!!!!!!!!!!!!!!!
Expand Down Expand Up @@ -540,7 +552,7 @@ module function create_OpticalProperties(datadir, optype, species_names, &
endif

filename = datadir//"/aerosol_xsections/"//sop%particle_xs(i)%dat//"/mie_"//sop%particle_xs(i)%dat//'.dat'
op%part(i) = create_ParticleXsection(filename, ind1, op%wavl, err)
op%part(i) = create_ParticleXsection(filename, ind1, sop%particle_xs(i)%dat, op%wavl, err)
if (allocated(err)) return

enddo
Expand Down Expand Up @@ -618,10 +630,11 @@ subroutine read_wavl(filename, optype, wavl, err)

end subroutine

function create_ParticleXsection(filename, p_ind, wavl, err) result(part)
function create_ParticleXsection(filename, p_ind, dat_name, wavl, err) result(part)
use futils, only: addpnt, inter2
character(*), intent(in) :: filename
integer, intent(in) :: p_ind
character(*), intent(in) :: dat_name
real(dp), intent(in) :: wavl(:)
character(:), allocatable, intent(out) :: err

Expand All @@ -637,6 +650,7 @@ function create_ParticleXsection(filename, p_ind, wavl, err) result(part)
integer :: i, j, io, ierr

part%p_ind = p_ind
part%dat_name = trim(dat_name)

open(2,file=filename,form="unformatted",status='old',iostat=io)
if (io /= 0) then
Expand Down Expand Up @@ -822,6 +836,8 @@ function create_WaterContinuum(model, filename, species_names, wavl, err) result
real(dp), parameter :: rdelta = 1.0e-4_dp

type(hdf5_file) :: h

cont%model = trim(model)

! Look for H2O
ind = findloc(species_names,"H2O",1)
Expand Down
2 changes: 2 additions & 0 deletions tests/test_radtran.f90
Original file line number Diff line number Diff line change
Expand Up @@ -144,6 +144,8 @@ subroutine test_Radtran()
write(1) rad%photons_sol
close(1)

call rad%print_opacities()

end subroutine

end program

0 comments on commit 87d6403

Please sign in to comment.