updated out2atmosphere_txt with optional decimals

CMakeLists.txt

@@ -1,7 +1,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.5.3")
+project(Clima LANGUAGES Fortran C VERSION "0.5.4")
 set(PHOTOCHEM_CLIMA_DATA_VERSION "0.2.0")
 
 set(CMAKE_Fortran_MODULE_DIRECTORY "${CMAKE_BINARY_DIR}/modules")

clima/cython/AdiabatClimate.pyx

@@ -441,16 +441,18 @@ cdef class AdiabatClimate:
     if len(err.strip()) > 0:
       raise ClimaException(err.decode("utf-8").strip())
 
-  def out2atmosphere_txt(self, str filename, ndarray[double, ndim=1] eddy, cbool overwrite = False, cbool clip = True):
+  def out2atmosphere_txt(self, str filename, ndarray[double, ndim=1] eddy, int number_of_decimals=5, cbool overwrite = False, cbool clip = True):
     """Saves state of the atmosphere to a file.
 
     Parameters
     ----------
     filename : str
         Output filename
-    eddy: ndarray[double,ndim=1]
+    eddy : ndarray[double,ndim=1]
         Array of eddy diffusions (cm^2/s) to write to the output file. This is
         useful for coupling to the photochemical model.
+    number_of_decimals : int
+        Number of decimals
     overwrite : bool, optional
         If true, then output file can be overwritten, by default False
     clip : bool, optional
@@ -461,7 +463,7 @@ cdef class AdiabatClimate:
     cdef bytes filename_b = pystring2cstring(filename)
     cdef char *filename_c = filename_b
     cdef char err[ERR_LEN+1]
-    wa_pxd.adiabatclimate_out2atmosphere_txt_wrapper(self._ptr, filename_c, &nz, <double *>eddy.data, &overwrite, &clip, err)  
+    wa_pxd.adiabatclimate_out2atmosphere_txt_wrapper(self._ptr, filename_c, &nz, <double *>eddy.data, &number_of_decimals, &overwrite, &clip, err)  
     if len(err.strip()) > 0:
       raise ClimaException(err.decode("utf-8").strip())
 

clima/cython/AdiabatClimate_pxd.pxd

@@ -73,7 +73,7 @@ cdef extern void adiabatclimate_set_ocean_solubility_fcn_wrapper(AdiabatClimate
 
 cdef extern void adiabatclimate_to_regular_grid_wrapper(AdiabatClimate *ptr, char *err)
 
-cdef extern void adiabatclimate_out2atmosphere_txt_wrapper(AdiabatClimate *ptr, char *filename, int *nz, double *eddy, bool *overwrite, bool *clip, char *err)
+cdef extern void adiabatclimate_out2atmosphere_txt_wrapper(AdiabatClimate *ptr, char *filename, int *nz, double *eddy, int *number_of_decimals, bool *overwrite, bool *clip, char *err)
 
 cdef extern void adiabatclimate_heat_redistribution_parameters_wrapper(AdiabatClimate *ptr, double *tau_LW, double *k_term, double *f_term, char *err)
 

clima/fortran/AdiabatClimate.f90

@@ -433,12 +433,13 @@ subroutine adiabatclimate_to_regular_grid_wrapper(ptr, err) bind(c)
 
 end subroutine
 
-subroutine adiabatclimate_out2atmosphere_txt_wrapper(ptr, filename, nz, eddy, overwrite, clip, err) bind(c)
+subroutine adiabatclimate_out2atmosphere_txt_wrapper(ptr, filename, nz, eddy, number_of_decimals, overwrite, clip, err) bind(c)
   use clima, only: AdiabatClimate
   type(c_ptr), value, intent(in) :: ptr
   character(kind=c_char), intent(in) :: filename(*)
   integer(c_int), intent(in) :: nz
   real(c_double), intent(in) :: eddy(nz)
+  integer(c_int), intent(in) :: number_of_decimals
   logical(c_bool), intent(in) :: overwrite, clip
   character(kind=c_char), intent(out) :: err(err_len+1)
 
@@ -454,7 +455,7 @@ subroutine adiabatclimate_out2atmosphere_txt_wrapper(ptr, filename, nz, eddy, ov
   overwrite_f = overwrite
   clip_f = clip
 
-  call c%out2atmosphere_txt(filename_f, eddy, overwrite_f, clip_f, err_f)
+  call c%out2atmosphere_txt(filename_f, eddy, number_of_decimals, overwrite_f, clip_f, err_f)
   err(1) = c_null_char
   if (allocated(err_f)) then
     call copy_string_ftoc(err_f, err)

src/adiabat/clima_adiabat.f90

@@ -1049,17 +1049,22 @@ subroutine AdiabatClimate_to_regular_grid(self, err)
 
   end subroutine
 
-  subroutine AdiabatClimate_out2atmosphere_txt(self, filename, eddy, overwrite, clip, err)
+  subroutine AdiabatClimate_out2atmosphere_txt(self, filename, eddy, number_of_decimals, overwrite, clip, err)
     use futils, only: FileCloser
     class(AdiabatClimate), target, intent(inout) :: self
     character(len=*), intent(in) :: filename
     real(dp), intent(in) :: eddy(:)
+    integer, intent(in) :: number_of_decimals
     logical, intent(in) :: overwrite, clip
     character(:), allocatable, intent(out) :: err
-    type(FileCloser) :: file
 
     character(len=100) :: tmp
     integer :: io, i, j
+    integer :: number_of_spaces
+    character(len=10) :: number_of_decimals_str, number_of_spaces_str
+    character(:), allocatable :: fmt_label, fmt_number
+    real(dp) :: clip_value
+    type(FileCloser) :: file
 
     call self%to_regular_grid(err)
     if (allocated(err)) return
@@ -1068,6 +1073,12 @@ subroutine AdiabatClimate_out2atmosphere_txt(self, filename, eddy, overwrite, cl
       err = '"eddy" has the wrong size'
       return
     endif
+
+    if (clip) then
+      clip_value = 1.0e-40_dp
+    else
+      clip_value = - huge(1.0_dp)
+    endif
 
     if (overwrite) then
       open(2, file=filename, form='formatted', status='replace', iostat=io)
@@ -1084,35 +1095,58 @@ subroutine AdiabatClimate_out2atmosphere_txt(self, filename, eddy, overwrite, cl
         return
       endif
     endif
+
+    ! number of decimals must be reasonable
+    if (number_of_decimals < 2 .or. number_of_decimals > 17) then
+      err = '"number_of_decimals" should be between 1 and 17.'
+      return
+    endif
+    number_of_spaces = number_of_decimals + 9
+    ! make sure number of spaces works with the length of species names
+    do i = 1,self%sp%ng
+      number_of_spaces = max(number_of_spaces,len_trim(self%species_names(i)) + 3)
+    enddo
+    write(number_of_decimals_str,'(i10)') number_of_decimals
+    write(number_of_spaces_str,'(i10)') number_of_spaces
+
+    fmt_label = "(a"//trim(adjustl(number_of_spaces_str))//")"
+    fmt_number = "(es"//trim(adjustl(number_of_spaces_str))//"."//trim(adjustl(number_of_decimals_str))//"e3)"
 
     tmp = 'alt'
-    write(unit=2,fmt="(3x,a27)",advance='no') tmp
+    write(unit=2,fmt=fmt_label,advance='no') tmp
     tmp = 'press'
-    write(unit=2,fmt="(a27)",advance='no') tmp
+    write(unit=2,fmt=fmt_label,advance='no') tmp
     tmp = 'den'
-    write(unit=2,fmt="(a27)",advance='no') tmp
+    write(unit=2,fmt=fmt_label,advance='no') tmp
     tmp = 'temp'
-    write(unit=2,fmt="(a27)",advance='no') tmp
+    write(unit=2,fmt=fmt_label,advance='no') tmp
     tmp = 'eddy'
-    write(unit=2,fmt="(a27)",advance='no') tmp
+    write(unit=2,fmt=fmt_label,advance='no') tmp
     do j = 1,self%sp%ng
       tmp = self%species_names(j)
-      write(unit=2,fmt="(a27)",advance='no') tmp
+      write(unit=2,fmt=fmt_label,advance='no') tmp
     enddo
 
     do i = 1,self%nz
       write(2,*)
-      write(unit=2,fmt="(es27.17e3)",advance='no') self%z(i)/1.e5_dp
-      write(unit=2,fmt="(es27.17e3)",advance='no') self%P(i)/1.e6_dp
-      write(unit=2,fmt="(es27.17e3)",advance='no') sum(self%densities(i,:))
-      write(unit=2,fmt="(es27.17e3)",advance='no') self%T(i)
-      write(unit=2,fmt="(es27.17e3)",advance='no') eddy(i)
+      write(tmp,fmt=fmt_number) self%z(i)/1.e5_dp
+      write(unit=2,fmt=fmt_label,advance='no') adjustl(tmp)
+
+      write(tmp,fmt=fmt_number) self%P(i)/1.e6_dp
+      write(unit=2,fmt=fmt_label,advance='no') adjustl(tmp)
+
+      write(tmp,fmt=fmt_number) sum(self%densities(i,:))
+      write(unit=2,fmt=fmt_label,advance='no') adjustl(tmp)
+
+      write(tmp,fmt=fmt_number) self%T(i)
+      write(unit=2,fmt=fmt_label,advance='no') adjustl(tmp)
+
+      write(tmp,fmt=fmt_number) eddy(i)
+      write(unit=2,fmt=fmt_label,advance='no') adjustl(tmp)
+
       do j = 1,self%sp%ng
-        if (clip) then
-          write(unit=2,fmt="(es27.17e3)",advance='no') max(self%f_i(i,j),1.0e-40_dp)
-        else
-          write(unit=2,fmt="(es27.17e3)",advance='no') self%f_i(i,j)
-        endif
+        write(tmp,fmt=fmt_number) max(self%f_i(i,j), clip_value)
+        write(unit=2,fmt=fmt_label,advance='no') adjustl(tmp)
       enddo
     enddo
 

tests/test_adiabat.f90

@@ -7,7 +7,7 @@ program test
   type(AdiabatClimate) :: c
   character(:), allocatable :: err
   real(dp) :: T, OLR, ISR, OLR1, ISR1
-  real(dp), allocatable :: P_i_surf(:), N_i_surf(:), f_i(:,:)
+  real(dp), allocatable :: P_i_surf(:), N_i_surf(:), f_i(:,:), eddy(:)
   integer :: i
   logical :: converged
   procedure(ocean_solubility_fcn), pointer :: ocean_fcn_ptr
@@ -36,6 +36,15 @@ program test
   endif
   print*,T, c%T_trop
 
+  allocate(eddy(size(c%P)))
+  eddy = 1.0e6_dp
+  call c%out2atmosphere_txt('test.txt', eddy, 5, .true., .true., err)
+  if (allocated(err)) then
+    print*,err
+    stop 1
+  endif
+  deallocate(eddy)
+
   ! Test surface_temperature_column
   c%tidally_locked_dayside = .false.
   N_i_surf = [15.0e3_dp, 400e-6_dp*23.0_dp, 1.0*36.0_dp, 1.0e-10_dp, 1.0e-10_dp, 1.0e-10_dp, 1.0e-10_dp]