From d241df6258cc16bdb9338cddb3d24101e23b4815 Mon Sep 17 00:00:00 2001 From: Nick Wogan Date: Thu, 26 Sep 2024 14:13:10 -0700 Subject: [PATCH] moved to new xsections --- .gitignore | 22 +- CMakeLists.txt | 1 + clima/cython/AdiabatClimate.pyx | 2 +- clima/cython/_clima.pyx | 1 + src/clima_useful.f90 | 16 + src/dependencies/CMakeLists.txt | 4 +- src/radtran/clima_radtran_types_create.f90 | 328 ++++++++------------- tests/CMakeLists.txt | 15 + tests/test_adiabat.f90 | 4 +- tests/test_climate.f90 | 2 +- tests/test_radtran.f90 | 2 +- 11 files changed, 161 insertions(+), 236 deletions(-) diff --git a/.gitignore b/.gitignore index 1ad9409..dd65195 100644 --- a/.gitignore +++ b/.gitignore @@ -1,26 +1,8 @@ .DS_Store .ipynb_checkpoints __pycache__ -*.so -*.o -Photo.run -*.mod build -*.a -a.out -test.run -*.dSYM -.vscode - - -src/Stringifor -build_dir -src/dependencies/lib -src/dependencies/include -src/dependencies/examples -src/dependencies/cmake -src/dependencies/modules - -random.ipynb +data +*.so diff --git a/CMakeLists.txt b/CMakeLists.txt index 99d8451..1584d50 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,6 +2,7 @@ cmake_minimum_required(VERSION "3.14") cmake_policy(SET CMP0148 OLD) # To suppress a warning emerging from scikit-build project(Clima LANGUAGES Fortran C VERSION "0.4.7") +set(PHOTOCHEM_CLIMA_DATA_VERSION "0.2.0") set(CMAKE_Fortran_MODULE_DIRECTORY "${CMAKE_BINARY_DIR}/modules") diff --git a/clima/cython/AdiabatClimate.pyx b/clima/cython/AdiabatClimate.pyx index 2fc3899..72112c8 100644 --- a/clima/cython/AdiabatClimate.pyx +++ b/clima/cython/AdiabatClimate.pyx @@ -48,7 +48,7 @@ cdef class AdiabatClimate: self._init_called = True if data_dir == None: - data_dir_ = os.path.dirname(os.path.realpath(__file__))+'/data' + data_dir_ = photochem_clima_data.DATA_DIR else: data_dir_ = data_dir diff --git a/clima/cython/_clima.pyx b/clima/cython/_clima.pyx index db0f3df..fda8541 100644 --- a/clima/cython/_clima.pyx +++ b/clima/cython/_clima.pyx @@ -5,6 +5,7 @@ from cpython.object cimport PyObject_GenericSetAttr import numpy as np import ctypes as ct import os +import photochem_clima_data DEF S_STR_LEN = 20; DEF ERR_LEN = 1024; diff --git a/src/clima_useful.f90 b/src/clima_useful.f90 index 95ab4dd..c7eea4f 100644 --- a/src/clima_useful.f90 +++ b/src/clima_useful.f90 @@ -1,12 +1,21 @@ module clima_useful use clima_const, only: dp use minpack_module, only: fcn_hybrj + use h5fortran, only: hdf5_file implicit none private + public :: hdf5_file_closer public :: MinpackHybrd1Vars, MinpackHybrj public :: linear_solve + ! Helps close hdf5 files + type :: hdf5_file_closer + type(hdf5_file), pointer :: h => null() + contains + final :: hdf5_file_closer_final + end type + type :: MinpackHybrd1Vars integer :: n real(dp), allocatable :: x(:) @@ -64,6 +73,13 @@ subroutine dgesv(n, nrhs, a, lda, ipiv, b, ldb, info) contains + subroutine hdf5_file_closer_final(self) + type(hdf5_file_closer), intent(inout) :: self + if (associated(self%h)) then + if (self%h%is_open) call self%h%close() + endif + end subroutine + function create_MinpackHybrd1Vars(n, tol) result(v) integer, intent(in) :: n real(dp), optional, intent(in) :: tol diff --git a/src/dependencies/CMakeLists.txt b/src/dependencies/CMakeLists.txt index 65d39b8..03d92c6 100644 --- a/src/dependencies/CMakeLists.txt +++ b/src/dependencies/CMakeLists.txt @@ -32,9 +32,9 @@ CPMAddPackage( # futils CPMAddPackage( NAME futils - VERSION 0.1.13 + VERSION 0.1.14 GITHUB_REPOSITORY "Nicholaswogan/futils" - GIT_TAG "v0.1.13" + GIT_TAG "v0.1.14" EXCLUDE_FROM_ALL ON ) diff --git a/src/radtran/clima_radtran_types_create.f90 b/src/radtran/clima_radtran_types_create.f90 index d2c8ed0..1338df7 100644 --- a/src/radtran/clima_radtran_types_create.f90 +++ b/src/radtran/clima_radtran_types_create.f90 @@ -488,7 +488,7 @@ module function create_OpticalProperties(datadir, optype, species_names, & ! need to go see what photolysis data is avaliable j = 0 do i = 1,size(species_names) - filename = datadir//"/xsections/"//trim(species_names(i))//"/"//trim(species_names(i))//"_xs.txt" + filename = datadir//"/xsections/"//trim(species_names(i))//".h5" inquire(file=filename, exist=file_exists) if (file_exists) j = j + 1 enddo @@ -498,7 +498,7 @@ module function create_OpticalProperties(datadir, optype, species_names, & allocate(tmp_str_list(op%npxs)) j = 1 do i = 1,size(species_names) - filename = datadir//"/xsections/"//trim(species_names(i))//"/"//trim(species_names(i))//"_xs.txt" + filename = datadir//"/xsections/"//trim(species_names(i))//".h5" inquire(file=filename, exist=file_exists) if (file_exists) then tmp_str_list(j) = species_names(i) @@ -522,7 +522,7 @@ module function create_OpticalProperties(datadir, optype, species_names, & '"photolysis-xs" is not in the list of species.' return endif - filename = datadir//"/xsections/"//trim(tmp_str_list(i))//"/"//trim(tmp_str_list(i))//"_xs.txt" + filename = datadir//"/xsections/"//trim(tmp_str_list(i))//".h5" op%pxs(i) = create_PhotolysisXsection(filename, trim(tmp_str_list(i)), ind1, op%wavl, err) if (allocated(err)) return @@ -547,7 +547,7 @@ module function create_OpticalProperties(datadir, optype, species_names, & return endif - filename = datadir//"/aerosol_xsections/"//sop%particle_xs(i)%dat//"/mie_"//sop%particle_xs(i)%dat//'.dat' + filename = datadir//"/aerosol_xsections/"//sop%particle_xs(i)%dat//"/mie_"//sop%particle_xs(i)%dat//'.h5' op%part(i) = create_ParticleXsection(filename, ind1, sop%particle_xs(i)%dat, op%wavl, err) if (allocated(err)) return @@ -627,7 +627,9 @@ subroutine read_wavl(filename, optype, wavl, err) end subroutine function create_ParticleXsection(filename, p_ind, dat_name, wavl, err) result(part) - use futils, only: addpnt, inter2 + use h5fortran + use clima_useful, only: hdf5_file_closer + use futils, only: interp_discrete_to_bins character(*), intent(in) :: filename integer, intent(in) :: p_ind character(*), intent(in) :: dat_name @@ -636,157 +638,103 @@ function create_ParticleXsection(filename, p_ind, dat_name, wavl, err) result(pa type(ParticleXsection) :: part - real(dp), allocatable :: wavl_tmp(:) - real(dp), allocatable :: w0_tmp(:,:), qext_tmp(:,:), g_tmp(:,:) - real(dp), allocatable :: w0_file(:,:), qext_file(:,:), g_file(:,:) - real(dp), allocatable :: temp_data(:), temp_wavelength(:) - - integer :: nw_tmp - real(dp) :: dum - integer :: i, j, io, ierr + type(hdf5_file), target :: h + integer(HSIZE_T), allocatable :: dims(:) + real(dp), allocatable :: wv_tmp(:) + real(dp), allocatable :: w0_tmp(:,:), qext_tmp(:,:), g0_tmp(:,:) + real(dp), allocatable :: w0_file(:,:), qext_file(:,:), g0_file(:,:) + integer :: i, 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 + if (.not.is_hdf5(filename)) then err = "Was unable to open mie data file "//trim(filename) - close(2) return endif - read(2, iostat=io) nw_tmp - if (io /= 0) then - err = "Problem reading mie data file "//trim(filename) - close(2) - return - endif - allocate(wavl_tmp(nw_tmp)) - read(2, iostat=io) wavl_tmp - if (io /= 0) then - err = "Problem reading mie data file "//trim(filename) - close(2) - return - endif - - read(2, iostat=io) part%nrad - if (io /= 0) then - err = "Problem reading mie data file "//trim(filename) - close(2) - return - endif - allocate(part%radii(part%nrad)) - read(2, iostat=io) part%radii - if (io /= 0) then - err = "Problem reading mie data file "//trim(filename) - close(2) - return - endif + block + type(hdf5_file_closer) :: h_closer + + ! Open file + call h%open(filename,'r') + h_closer%h => h + + ! Wavelengths + call check_h5_dataset(h, 'wavelengths', 1, H5T_FLOAT_F, filename, err) + if (allocated(err)) return + call h%shape("wavelengths", dims) + allocate(wv_tmp(dims(1))) + call h%read("wavelengths", wv_tmp) + + ! Radii + call check_h5_dataset(h, 'radii', 1, H5T_FLOAT_F, filename, err) + if (allocated(err)) return + call h%shape("radii", dims) + part%nrad = dims(1) + allocate(part%radii(dims(1))) + call h%read("radii", part%radii) + ! convert from micron to cm part%radii = part%radii/1.0e4_dp part%r_min = minval(part%radii) part%r_max = maxval(part%radii) - - allocate(w0_tmp(nw_tmp,part%nrad)) - allocate(qext_tmp(nw_tmp,part%nrad)) - allocate(g_tmp(nw_tmp,part%nrad)) - - read(2, iostat=io) w0_tmp - if (io /= 0) then - err = "Problem reading mie data file "//trim(filename) - close(2) - return - endif - read(2, iostat=io) qext_tmp - if (io /= 0) then - err = "Problem reading mie data file "//trim(filename) - close(2) + + ! w0 + call check_h5_dataset(h, 'w0', 2, H5T_FLOAT_F, filename, err) + if (allocated(err)) return + call h%shape("w0", dims) + allocate(w0_tmp(dims(1),dims(2))) + call h%read("w0", w0_tmp) + + ! qext + call check_h5_dataset(h, 'qext', 2, H5T_FLOAT_F, filename, err) + if (allocated(err)) return + call h%shape("qext", dims) + allocate(qext_tmp(dims(1),dims(2))) + call h%read("qext", qext_tmp) + + ! g0 + call check_h5_dataset(h, 'g0', 2, H5T_FLOAT_F, filename, err) + if (allocated(err)) return + call h%shape("g0", dims) + allocate(g0_tmp(dims(1),dims(2))) + call h%read("g0", g0_tmp) + + if (size(w0_tmp,1) /= part%nrad .or. size(w0_tmp,2) /= size(wv_tmp)) then + err = '"w0" has the wrong shape in "'//filename//'"' return endif - read(2, iostat=io) g_tmp - if (io /= 0) then - err = "Problem reading mie data file "//trim(filename) - close(2) + if (size(qext_tmp,1) /= part%nrad .or. size(qext_tmp,2) /= size(wv_tmp)) then + err = '"qext" has the wrong shape in "'//filename//'"' return endif - - ! this next read should be usuccessful - read(2, iostat=io) dum - if (io == 0) then - err = "Problem reading mie data file "//trim(filename)// & - ". Should have reached the end of the file, but did not." - close(2) + if (size(g0_tmp,1) /= part%nrad .or. size(g0_tmp,2) /= size(wv_tmp)) then + err = '"g0" has the wrong shape in "'//filename//'"' return endif - close(2) - - ! now lets interpolate a bunch + + end block + + ! Now lets interpolate a bunch allocate(w0_file(part%nrad,size(wavl)-1)) allocate(qext_file(part%nrad,size(wavl)-1)) - allocate(g_file(part%nrad,size(wavl)-1)) - - allocate(temp_data(nw_tmp+2)) ! for interpolation - allocate(temp_wavelength(nw_tmp+2)) - + allocate(g0_file(part%nrad,size(wavl)-1)) + do i = 1, part%nrad - - ! w0 (single scattering albedo) - temp_data(1:nw_tmp) = w0_tmp(1:nw_tmp,i) - temp_wavelength(1:nw_tmp) = wavl_tmp(1:nw_tmp) - j = nw_tmp - call addpnt(temp_wavelength, temp_data, nw_tmp+2, j, 0.0_dp, w0_tmp(1,i), ierr) - call addpnt(temp_wavelength, temp_data, nw_tmp+2, j, huge(1.0_dp), w0_tmp(nw_tmp,i), ierr) - if (ierr /= 0) then - err = "Problems interpolating mie data from file "//trim(filename)// & - " to the wavelength grid" - return - endif - call inter2(size(wavl), wavl, w0_file(i,:), & - nw_tmp+2, temp_wavelength, temp_data, ierr) - if (ierr /= 0) then - err = "Problems interpolating mie data from file "//trim(filename)// & - " to the wavelength grid" - return - endif - ! qext (The extinction efficiency) - temp_data(1:nw_tmp) = qext_tmp(1:nw_tmp,i) - temp_wavelength(1:nw_tmp) = wavl_tmp(1:nw_tmp) - j = nw_tmp - call addpnt(temp_wavelength, temp_data, nw_tmp+2, j, 0.0_dp, qext_tmp(1,i), ierr) - call addpnt(temp_wavelength, temp_data, nw_tmp+2, j, huge(1.0_dp), qext_tmp(nw_tmp,i), ierr) - if (ierr /= 0) then - err = "Problems interpolating mie data from file "//trim(filename)// & - " to the wavelength grid" - return - endif - call inter2(size(wavl), wavl, qext_file(i,:), & - nw_tmp+2, temp_wavelength, temp_data, ierr) - if (ierr /= 0) then - err = "Problems interpolating mie data from file "//trim(filename)// & - " to the wavelength grid" - return - endif - - ! g (The scattering anisotropy or asymmetry factor) - temp_data(1:nw_tmp) = g_tmp(1:nw_tmp,i) - temp_wavelength(1:nw_tmp) = wavl_tmp(1:nw_tmp) - j = nw_tmp - call addpnt(temp_wavelength, temp_data, nw_tmp+2, j, 0.0_dp, g_tmp(1,i), ierr) - call addpnt(temp_wavelength, temp_data, nw_tmp+2, j, huge(1.0_dp), g_tmp(nw_tmp,i), ierr) - if (ierr /= 0) then - err = "Problems interpolating mie data from file "//trim(filename)// & - " to the wavelength grid" - return - endif - call inter2(size(wavl), wavl, g_file(i,:), & - nw_tmp+2, temp_wavelength, temp_data, ierr) - if (ierr /= 0) then - err = "Problems interpolating mie data from file "//trim(filename)// & - " to the wavelength grid" - return - endif - + ! w0 + call interp_discrete_to_bins(wavl, wv_tmp, w0_tmp(i,:), w0_file(i,:), 'Constant', err=err) + if (allocated(err)) return + + ! ext + call interp_discrete_to_bins(wavl, wv_tmp, qext_tmp(i,:), qext_file(i,:), 'Constant', err=err) + if (allocated(err)) return + + ! g0 + call interp_discrete_to_bins(wavl, wv_tmp, g0_tmp(i,:), g0_file(i,:), 'Constant', err=err) + if (allocated(err)) return + enddo allocate(part%w0(size(wavl) - 1)) @@ -803,7 +751,7 @@ function create_ParticleXsection(filename, p_ind, dat_name, wavl, err) result(pa err = 'Failed to initialize interpolator for "'//filename//'"' return endif - call part%gt(i)%initialize(part%radii, g_file(:,i), ierr) + call part%gt(i)%initialize(part%radii, g0_file(:,i), ierr) if (ierr /= 0) then err = 'Failed to initialize interpolator for "'//filename//'"' return @@ -1366,8 +1314,10 @@ subroutine check_h5_dataset(h, dataset, ndims, dtype, prefix, err) end subroutine function create_PhotolysisXsection(xsfilename, sp, sp_ind, wavl, err) result(xs) - use clima_const, only: s_str_len, log10tiny - use futils, only: inter2, addpnt + use h5fortran + use clima_useful, only: hdf5_file_closer + use clima_const, only: log10tiny + use futils, only: interp_discrete_to_bins character(*), intent(in) :: xsfilename character(*), intent(in) :: sp integer, intent(in) :: sp_ind @@ -1376,93 +1326,53 @@ function create_PhotolysisXsection(xsfilename, sp, sp_ind, wavl, err) result(xs) type(Xsection) :: xs - integer, parameter :: maxcols = 200 - character(len=10000) :: line - character(len=100) :: tmp(maxcols) - real(dp), parameter :: rdelta = 1.0e-4_dp - real(dp), allocatable :: wav_f(:), log10_xs_0d(:) - - integer :: k, kk, io, ierr + type(hdf5_file), target :: h + integer(HSIZE_T), allocatable :: dims(:) + real(dp), allocatable :: wv_tmp(:), xs_tmp(:) xs%xs_type = PhotolysisXsection allocate(xs%sp_ind(1)) xs%sp_ind(1) = sp_ind - - open(101, file=xsfilename,status='old',iostat=io) - if (io /= 0) then + xs%dim = 0 + + if (.not.is_hdf5(xsfilename)) then err = 'Species "'//sp//'" does not have photolysis xsection data' return endif - read(101,*) - read(101,'(A)') line - do k=1,maxcols - read(line,*,iostat=io) tmp(1:k) - if (io /= 0) exit - enddo - if (k == maxcols+1) then - err = 'More cross section temperature data than allowed for '//sp - return - endif + ! Block to close hdf5 file when we leave this function + block + type(hdf5_file_closer) :: h_closer - if (k - 2 == 1) then - xs%dim = 0 - else - xs%dim = 1 - endif + ! Open file + call h%open(xsfilename,'r') + h_closer%h => h - if (xs%dim == 0) then + ! Wavelengths + call check_h5_dataset(h, 'wavelengths', 1, H5T_FLOAT_F, xsfilename, err) + if (allocated(err)) return + call h%shape("wavelengths", dims) + allocate(wv_tmp(dims(1))) + call h%read("wavelengths", wv_tmp) - ! count lines - k = 0 - io = 0 - do while(io == 0) - read(101,*,iostat=io) - if (io == 0) k = k + 1 - enddo + ! Photoabsorption + call check_h5_dataset(h, 'photoabsorption', 1, H5T_FLOAT_F, xsfilename, err) + if (allocated(err)) return + call h%shape("photoabsorption", dims) + allocate(xs_tmp(dims(1))) + call h%read("photoabsorption", xs_tmp) - allocate(wav_f(k+4)) - allocate(log10_xs_0d(k+4)) - log10_xs_0d = 0.0_dp - - rewind(101) - read(101,*) - read(101,*) - do k = 1,size(wav_f)-4 - read(101,*,iostat=io) wav_f(k), log10_xs_0d(k) - if (io/=0) then - err = 'Problem readling "'//xsfilename//'"' - return - endif - enddo - close(101) + xs_tmp = max(xs_tmp,tiny(1.0_dp)) + xs_tmp = log10(xs_tmp) - log10_xs_0d = log10_xs_0d + tiny(1.0_dp) - log10_xs_0d = log10(log10_xs_0d) + ! Interpolate + allocate(xs%xs_0d(size(wavl)-1)) + call interp_discrete_to_bins(wavl, wv_tmp, xs_tmp, xs%xs_0d, 'FillValue', log10tiny, err) + if (allocated(err)) return - kk = size(wav_f) - k = size(wav_f)-4 - call addpnt(wav_f, log10_xs_0d, kk, k, wav_f(1)*(1.0_dp-rdelta), log10tiny,ierr) - call addpnt(wav_f, log10_xs_0d, kk, k, 0.0_dp, log10tiny,ierr) - call addpnt(wav_f, log10_xs_0d, kk, k, wav_f(k)*(1.0_dp+rdelta), log10tiny,ierr) - call addpnt(wav_f, log10_xs_0d, kk, k, huge(rdelta),log10tiny,ierr) - if (ierr /= 0) then - err = 'Problem interpolating data in "'//trim(xsfilename)//'"' - return - endif - allocate(xs%xs_0d(size(wavl)-1)) - call inter2(size(wavl), wavl, xs%xs_0d, & - kk, wav_f, log10_xs_0d, ierr) - if (ierr /= 0) then - err = 'Problem interpolating data in "'//trim(xsfilename)//'"' - return - endif - xs%xs_0d = 10.0_dp**xs%xs_0d + xs%xs_0d = 10.0_dp**xs%xs_0d - elseif (xs%dim == 1) then - err = 'Photolysis cross sections for "'//sp//'" can not depend on T' - return - endif + end block end function diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index e1cd52f..375f2a2 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,3 +1,18 @@ +# First download data, if needed +if (NOT EXISTS "${CMAKE_SOURCE_DIR}/data") + file(DOWNLOAD + "https://github.com/Nicholaswogan/photochem_clima_data/archive/v${PHOTOCHEM_CLIMA_DATA_VERSION}.tar.gz" + "photochem_clima_data_v${PHOTOCHEM_CLIMA_DATA_VERSION}.tar.gz" + ) + file(ARCHIVE_EXTRACT + INPUT "photochem_clima_data_v${PHOTOCHEM_CLIMA_DATA_VERSION}.tar.gz" + DESTINATION ${CMAKE_BINARY_DIR} + ) + file(RENAME + "${CMAKE_BINARY_DIR}/photochem_clima_data-${PHOTOCHEM_CLIMA_DATA_VERSION}/photochem_clima_data/data" + "${CMAKE_SOURCE_DIR}/data" + ) +endif() set(CLIMA_TESTS test_radtran test_climate test_adiabat) diff --git a/tests/test_adiabat.f90 b/tests/test_adiabat.f90 index cebe871..9ce4583 100644 --- a/tests/test_adiabat.f90 +++ b/tests/test_adiabat.f90 @@ -15,7 +15,7 @@ program test c = AdiabatClimate('../templates/AdiabatClimate/species.yaml', & '../templates/AdiabatClimate/Earth/settings.yaml', & '../templates/ModernEarth/Sun_now.txt', & - '../clima/data', & + '../data', & err=err) if (allocated(err)) then print*,err @@ -143,7 +143,7 @@ subroutine test_RCE() c = AdiabatClimate('../templates/AdiabatClimate/species.yaml', & '../tests/settings_RCE_test.yaml', & '../templates/ModernEarth/Sun_now.txt', & - '../clima/data', & + '../data', & err=err) if (allocated(err)) then print*,err diff --git a/tests/test_climate.f90 b/tests/test_climate.f90 index 4efa393..e479a83 100644 --- a/tests/test_climate.f90 +++ b/tests/test_climate.f90 @@ -20,7 +20,7 @@ subroutine test_Climate_evolve() "../templates/ModernEarth/settings.yaml", & "../templates/ModernEarth/Sun_now.txt", & "../templates/ModernEarth/atmosphere.txt", & - "../clima/data", & + "../data", & err) if (allocated(err)) then print*,err diff --git a/tests/test_radtran.f90 b/tests/test_radtran.f90 index 3127245..c360cd9 100644 --- a/tests/test_radtran.f90 +++ b/tests/test_radtran.f90 @@ -35,7 +35,7 @@ subroutine test_Radtran() num_zeniths = 8 surface_albedo = 0.15_dp rad = Radtran("../templates/ModernEarth/settings.yaml",& - "../templates/ModernEarth/Sun_now.txt", num_zeniths, surface_albedo, nz, "../clima/data", err) + "../templates/ModernEarth/Sun_now.txt", num_zeniths, surface_albedo, nz, "../data", err) if (allocated(err)) then print*,err stop 1