From 87d6403177d4c92da7f7aa58e717104d3db9819b Mon Sep 17 00:00:00 2001
From: Nick Wogan <nicholaswogan@gmail.com>
Date: Thu, 26 Oct 2023 20:07:23 -0700
Subject: [PATCH] added function to make string of opacities

---
 src/radtran/clima_radtran.f90              |  59 +++++++++
 src/radtran/clima_radtran_types.f90        | 135 ++++++++++++++++++++-
 src/radtran/clima_radtran_types_create.f90 |  22 +++-
 tests/test_radtran.f90                     |   2 +
 4 files changed, 212 insertions(+), 6 deletions(-)

diff --git a/src/radtran/clima_radtran.f90 b/src/radtran/clima_radtran.f90
index b2a4479..b082ca4 100644
--- a/src/radtran/clima_radtran.f90
+++ b/src/radtran/clima_radtran.f90
@@ -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
@@ -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 !!!
   !!!!!!!!!!!!!!!!!
diff --git a/src/radtran/clima_radtran_types.f90 b/src/radtran/clima_radtran_types.f90
index 279cd93..75e8c7a 100644
--- a/src/radtran/clima_radtran_types.f90
+++ b/src/radtran/clima_radtran_types.f90
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
\ No newline at end of file
diff --git a/src/radtran/clima_radtran_types_create.f90 b/src/radtran/clima_radtran_types_create.f90
index 25df632..b48b0af 100644
--- a/src/radtran/clima_radtran_types_create.f90
+++ b/src/radtran/clima_radtran_types_create.f90
@@ -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
@@ -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 !!!
     !!!!!!!!!!!!!!!!!!!!!!!
@@ -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
@@ -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
 
@@ -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
@@ -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)
diff --git a/tests/test_radtran.f90 b/tests/test_radtran.f90
index 14d0efe..4bb6bf1 100644
--- a/tests/test_radtran.f90
+++ b/tests/test_radtran.f90
@@ -144,6 +144,8 @@ subroutine test_Radtran()
     write(1) rad%photons_sol
     close(1)
 
+    call rad%print_opacities()
+
   end subroutine
   
 end program
\ No newline at end of file