fixed multithread bug tau in radtran

CMakeLists.txt

@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION "3.14")
 
-project(Clima LANGUAGES Fortran C VERSION "0.3.9")
+project(Clima LANGUAGES Fortran C VERSION "0.3.10")
 
 set(CMAKE_Fortran_MODULE_DIRECTORY "${CMAKE_BINARY_DIR}/modules")
 

src/adiabat/clima_adiabat.f90

@@ -914,7 +914,7 @@ subroutine AdiabatClimate_heat_redistribution_parameters(self, tau_LW, k_term, f
       denominator = 0.0_dp
       do i = 1,self%rad%ir%nw
         dlambda = self%rad%ir%wavl(i+1) - self%rad%ir%wavl(i) ! Width of a wavelength bin (nm)
-        tau_lambda = sum(self%rad%wrk_ir%rx%tau_band(:,i)) ! IR optical depth
+        tau_lambda = sum(self%rad%wrk_ir%tau_band(:,i)) ! IR optical depth
 
         avg_freq = 0.5_dp*(self%rad%ir%freq(i) + self%rad%ir%freq(i+1))
         avg_lam = (c_light*1.0e9_dp/avg_freq) ! (m/s) * (nm/m) * (s) = nm

src/radtran/clima_radtran.f90

@@ -14,14 +14,14 @@ module clima_radtran
     type(RadiateXSWrk) :: rx
     type(RadiateZWrk) :: rz
 
-    !! (nz+1,nw) mW/m2/Hz in each wavelength bin
-    !! at the edges of the vertical grid
-    real(dp), allocatable :: fup_a(:,:)
-    real(dp), allocatable :: fdn_a(:,:)
-    !! (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(:)
-    real(dp), allocatable :: fdn_n(:) 
+    !> (nz+1,nw) mW/m2/Hz in each wavelength bin
+    !> at the edges of the vertical grid
+    real(dp), allocatable :: fup_a(:,:), fdn_a(:,:)
+    !> (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(:)
+    !> Band optical thickness (nz,nw)
+    real(dp), allocatable :: tau_band(:,:)
 
   end type
 
@@ -170,6 +170,7 @@ function create_RadtranIR_2(species_names, particle_names, s, nz, datadir, err)
     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%tau_band(nz,rad%ir%nw))
 
   end function
 
@@ -199,7 +200,7 @@ subroutine RadtranIR_radiate(self, T_surface, T, P, densities, dz, pdensities, r
                    P, T_surface, T, densities, dz, &
                    pdensities, radii, &
                    wrk%rx, wrk%rz, &
-                   wrk%fup_a, wrk%fdn_a, wrk%fup_n, wrk%fdn_n)
+                   wrk%fup_a, wrk%fdn_a, wrk%fup_n, wrk%fdn_n, wrk%tau_band)
     if (ierr /= 0) then
       err = 'Input particle radii are outside the data range.'
       return
@@ -326,6 +327,7 @@ 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%tau_band(nz,rad%ir%nw))
 
     ! Solar work arrays
     rad%wrk_sol%rx = RadiateXSWrk(rad%sol, nz)
@@ -334,6 +336,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%tau_band(nz,rad%sol%nw))
 
     ! total flux
     allocate(rad%f_total(nz+1))
@@ -371,7 +374,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%fup_a, wrk_ir%fdn_a, wrk_ir%fup_n, wrk_ir%fdn_n, wrk_ir%tau_band)
     if (ierr /= 0) then
       err = 'Input particle radii are outside the data range.'
       return
@@ -383,7 +386,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%fup_a, wrk_sol%fdn_a, wrk_sol%fup_n, wrk_sol%fdn_n, wrk_sol%tau_band)
     if (ierr /= 0) then
       err = 'Input particle radii are outside the data range.'
       return

src/radtran/clima_radtran_radiate.f90

@@ -17,7 +17,7 @@ function radiate(op, &
                    P, T_surface, T, densities, dz, &
                    pdensities, radii, &
                    rw, rz, &
-                   fup_a, fdn_a, fup_n, fdn_n) result(ierr)
+                   fup_a, fdn_a, fup_n, fdn_n, 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
@@ -52,6 +52,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) :: tau_band(:,:) !! (nz,nw) The optical depth of each layer
     integer :: ierr !! if ierr /= 0 on return, then there was an error
 
     type(Ksettings), pointer :: kset
@@ -301,7 +302,7 @@ function radiate(op, &
       enddo
       do i = 1,nz
         n = nz+1-i
-        rw%tau_band(i,l) = rz%tau_band(n) ! band optical thickness
+        tau_band(i,l) = rz%tau_band(n) ! band optical thickness
       enddo
 
     enddo

src/radtran/clima_radtran_types.f90

@@ -159,9 +159,6 @@ module function create_OpticalProperties(datadir, optype, species_names, &
     real(dp), allocatable :: w0(:,:)
     real(dp), allocatable :: qext(:,:)
     real(dp), allocatable :: gt(:,:)
-
-    !> The optical thickness in each band
-    real(dp), allocatable :: tau_band(:,:) ! (nz,nw)
 
     ! work arrays that are needed only if
     ! k_method == k_RandomOverlapResortRebin

src/radtran/clima_radtran_types_create.f90

@@ -103,9 +103,6 @@ module function create_RadiateXSWrk(op, nz) result(rw)
     allocate(rw%w0(nz,op%npart))
     allocate(rw%qext(nz,op%npart))
     allocate(rw%gt(nz,op%npart))
-
-    ! The band optical thickness
-    allocate(rw%tau_band(nz,op%nw))
 
     ! if there are k-distributions
     ! then we need to allocate some work arrays