Updates to solve error

Author: Cansu Uluseker

data.example/constants_test-scenario.nml

@@ -54,7 +54,7 @@
     river_attachment_efficiency = 0.01
     estuary_attachment_efficiency = 0.01
     estuary_tidal_m2 = 2.25
-    estuary_ tidal_s2 = 0.62
+    estuary_tidal_s2 = 0.62
     estuary_mean_depth_expa = 19.4
     estuary_mean_depth_expb = 1.46e-05
     estuary_width_expa = 12548
@@ -82,4 +82,4 @@
     urban_gardens = 150
     urban_parks = 250
     vertical_distribution = 50, 35, 15
-/
+/

fpm.toml

@@ -17,4 +17,5 @@ source-dir = "src"
 main = "main.f90"
 
 [build]
-link = ["netcdf", "netcdff"]
+link = ["netcdf", "netcdff"]
+

src/BedSediment/BedSedimentModule.f90

@@ -105,98 +105,113 @@ subroutine getMTCMatrix1(me, djdep, djres)
     !! Initialised `BedSediment` object, including all layers and included `FineSediment`
     !! objects
     function createBedSediment1(me, x, y, w) result(r)
-    class(BedSediment) :: me                                    !! Self-reference
-    integer :: x                                                !! x index of the containing water body
-    integer :: y                                                !! y index of the containing water body
-    integer :: w                                                !! w index of the containing water body
-    type(Result) :: r                                           !! Returned `Result` object
-    type(BedSedimentLayer), allocatable :: bsl1                 ! LOCAL object of type BedSedimentLayer, for implementation of polymorphism
-    integer :: L                                                ! LOCAL loop counter
-    integer :: allst                                            ! LOCAL array allocation status
-    character(len=256) :: tr                                    ! LOCAL error trace
-    character(len=16), parameter :: ms = "Allocation error"     ! LOCAL allocation error message
-    type(ErrorInstance) :: err(1)
+        class(BedSediment) :: me
+        integer            :: x, y, w
+        type(Result)       :: r
+        type(BedSedimentLayer), allocatable :: bsl1
+        integer            :: L, allst
+        character(len=256) :: tr
+        character(len=16), parameter :: ms = "Allocation error"
+        type(ErrorInstance) :: err(1)
+        integer :: nx, ny
+        logical :: inbounds
 
-    me%name = trim(ref('BedSediment', x, y, w))
-    me%nSizeClasses = C%nSizeClassesSpm                         ! set number of size classes from global value
-    me%nfComp = C%nFracCompsSpm                                 ! set number of compositional fractions from global value
-    tr = trim(me%name) // "%createBedSediment1"                 ! procedure name as trace
+        me%name = trim(ref('BedSediment', x, y, w))
+        ! >>> FIX: set indices so later DATASET(x,y,...) lookups are valid
+        me%x = x
+        me%y = y
+        ! <<<
 
-    ! Initialise Contaminant mass pools matrix
-    allocate(me%m_contaminant(C%nSedimentLayers + 3), stat=allst)
-    if (allst /= 0) then
-        err(1) = ErrorInstance(code=1, message=ms, trace=[tr])
-        call r%addError(err(1))
-        call LOGR%toFile(errors=err)
-        return
-    end if
-    do L = 1, C%nSedimentLayers + 3
-        ! Pass DATASET%nc and provide required arguments from DATASET
-        r = me%m_contaminant(L)%create_from_data( &
-            data=DATASET%nc, &
-            compartment='sediment', &
-            contaminantDensity=DATASET%contaminantDensity, &
-            soilAttachmentEfficiency=DATASET%soilConstantAttachmentEfficiency, &
-            riverAttachmentEfficiency=DATASET%riverAttachmentEfficiency, &
-            estuaryAttachmentEfficiency=DATASET%estuaryAttachmentEfficiency, &
-            k_diss_pristine=DATASET%contaminant_k_diss_pristine, &
-            k_diss_transformed=DATASET%contaminant_k_diss_transformed, &
-            k_transform_pristine=DATASET%contaminant_k_transform_pristine, &
-            waterTemperature=DATASET%waterTemperature(1) &  ! Use first timestep or adjust as needed
-        )
-        if (r%hasCriticalError()) then
-            call LOGR%toFile(errors=r%getErrors())
+        me%nSizeClasses = C%nSizeClassesSpm
+        me%nfComp       = C%nFracCompsSpm
+        tr = trim(me%name) // "%createBedSediment1"
+
+        ! Initialise Contaminant mass pools matrix
+        allocate(me%m_contaminant(C%nSedimentLayers + 3), stat=allst)
+        if (allst /= 0) then
+            err(1) = ErrorInstance(code=1, message=ms, trace=[tr])
+            call r%addError(err(1))
+            call LOGR%toFile(errors=err)
             return
         end if
-        if (L > 2 .and. allocated(DATASET%initialContaminantConcsSediment)) then
-            me%m_contaminant(L)%c = DATASET%initialContaminantConcsSediment(me%x, me%y, :, :, :)
-            if (allocated(DATASET%initialDissolvedConcsSediment)) then
-                me%m_contaminant(L)%m_dissolved = DATASET%initialDissolvedConcsSediment(me%x, me%y)
+
+        do L = 1, C%nSedimentLayers + 3
+            r = me%m_contaminant(L)%create_from_data( &
+                compartment='sediment', &  ! keep as used before; 'soil' also acceptable if your API expects it
+                contaminantDensity=DATASET%contaminantDensity, &
+                soilAttachmentEfficiency=DATASET%soilConstantAttachmentEfficiency, &
+                riverAttachmentEfficiency=DATASET%riverAttachmentEfficiency, &
+                estuaryAttachmentEfficiency=DATASET%estuaryAttachmentEfficiency, &
+                k_diss_pristine=DATASET%contaminant_k_diss_pristine, &
+                k_diss_transformed=DATASET%contaminant_k_diss_transformed, &
+                k_transform_pristine=DATASET%contaminant_k_transform_pristine, &
+                waterTemperature=real(DATASET%waterTemperature(1), dp) )
+            if (r%hasCriticalError()) then
+                call LOGR%toFile(errors=r%getErrors())
+                return
             end if
-        end if
-    end do
 
-    allocate(me%colBedSedimentLayers(C%nSedimentLayers), stat=allst)
-    if (allst /= 0) then
-        err(1) = ErrorInstance(code=1, message=ms, trace=[tr])
-        call r%addError(err(1))
-        call LOGR%toFile(errors=err)
-        return
-    end if
-    me%n_delta_sed = C%nSedimentLayers + 3
-    allocate(me%delta_sed(C%nSedimentLayers + 3, C%nSedimentLayers + 3, me%nSizeClasses), stat=allst)
-    if (allst /= 0) then
-        err(1) = ErrorInstance(code=1, message=ms, trace=[tr])
-        call r%addError(err(1))
-        call LOGR%toFile(errors=err)
-        return
-    end if
-    me%delta_sed = 0.0_dp
-    allocate(me%delta_sed_csr(C%nSizeClassesSpm), stat=allst)
-    if (allst /= 0) then
-        err(1) = ErrorInstance(code=1, message=ms, trace=[tr])
-        call r%addError(err(1))
-        call LOGR%toFile(errors=err)
-        return
-    end if
+            ! Seed initial concentrations for actual sediment layers (L>2) if provided
+            if (L > 2 .and. allocated(DATASET%initialContaminantConcsSediment)) then
+                nx = size(DATASET%initialContaminantConcsSediment, 1)
+                ny = size(DATASET%initialContaminantConcsSediment, 2)
+                inbounds = (me%x>=1 .and. me%y>=1 .and. me%x<=nx .and. me%y<=ny)
+                if (inbounds) then
+                    me%m_contaminant(L)%c = DATASET%initialContaminantConcsSediment(me%x, me%y, :, :, :)
+                    if (allocated(DATASET%initialDissolvedConcsSediment)) then
+                        me%m_contaminant(L)%m_dissolved = DATASET%initialDissolvedConcsSediment(me%x, me%y)
+                    end if
+                else
+                    me%m_contaminant(L)%c          = 0.0_dp
+                    me%m_contaminant(L)%m_dissolved = 0.0_dp
+                end if
+            end if
+        end do
 
-    do L = 1, C%nSedimentLayers
-        allocate(bsl1)
-        call r%addErrors(.errors. bsl1%create(L))
-        allocate(me%colBedSedimentLayers(L)%item, source=bsl1, stat=allst)
-        deallocate(bsl1)
+        allocate(me%colBedSedimentLayers(C%nSedimentLayers), stat=allst)
         if (allst /= 0) then
             err(1) = ErrorInstance(code=1, message=ms, trace=[tr])
             call r%addError(err(1))
             call LOGR%toFile(errors=err)
             return
         end if
-        if (r%hasCriticalError()) then
-            call r%addToTrace(tr)
+
+        me%n_delta_sed = C%nSedimentLayers + 3
+        allocate(me%delta_sed(C%nSedimentLayers + 3, C%nSedimentLayers + 3, me%nSizeClasses), stat=allst)
+        if (allst /= 0) then
+            err(1) = ErrorInstance(code=1, message=ms, trace=[tr])
+            call r%addError(err(1))
+            call LOGR%toFile(errors=err)
+            return
+        end if
+        me%delta_sed = 0.0_dp
+
+        allocate(me%delta_sed_csr(C%nSizeClassesSpm), stat=allst)
+        if (allst /= 0) then
+            err(1) = ErrorInstance(code=1, message=ms, trace=[tr])
+            call r%addError(err(1))
+            call LOGR%toFile(errors=err)
             return
         end if
-    end do
-end function
+
+        do L = 1, C%nSedimentLayers
+            allocate(bsl1)
+            call r%addErrors(.errors. bsl1%create(L))
+            allocate(me%colBedSedimentLayers(L)%item, source=bsl1, stat=allst)
+            deallocate(bsl1)
+            if (allst /= 0) then
+                err(1) = ErrorInstance(code=1, message=ms, trace=[tr])
+                call r%addError(err(1))
+                call LOGR%toFile(errors=err)
+                return
+            end if
+            if (r%hasCriticalError()) then
+                call r%addToTrace(tr)
+                return
+            end if
+        end do
+    end function
+
 
     !> **Function purpose**                                         <br>
     !! Deallocate all allocatable variables and call destroy methods for all

src/ConstantsDefaultsModule.f90

@@ -2,28 +2,31 @@
 module ConstantsDefaultsModule
     implicit none
 
+    integer, private, parameter :: dp = selected_real_kind(15, 307)
+
     ! Defaults for constants
-    real, parameter :: defaultSoilAttachmentEfficiency = 0.0
-    real, parameter :: defaultRiverAttachmentEfficiency = 0.0  ! Attachment efficiency for NM to SPM in rivers
-    real, parameter :: defaultSoilDarcyVelocity = 9e-6      ! [m/s] Tufenkji et al, 2004: https://doi.org/10.1021/es034049r 
-    real, parameter :: default_k_diss_pristine = 0.0        ! Dissolution rate for pristine contaminant [s-1]
-    real, parameter :: default_k_diss_transformed = 0.0     ! Dissolution rate for transformed contaminant [s-1]
-    real, parameter :: default_k_transform_pristine = 0.0   ! Transformation rate for pristine contaminant [s-1]
-    real, parameter :: default_rho_contaminant = 1000.0     ! Density of contaminant [kg/m3]
-    real, parameter :: defaultShearRate = 10.0              ! Arvidsson et al, 2009: https://doi.org/10.1080/10807039.2011.538639
-    real, parameter :: defaultMinWaterTemperature = 4.0     ! Thames River
-    real, parameter :: defaultMaxWaterTemperature = 21.0    ! Thames River
-    integer, parameter :: defaultMinWaterTemperatureDayOfYear = 32  ! Thames River
-    real, parameter :: defaultSedimentTransport_a = 2.0e-9
-    real, parameter :: defaultSedimentTransport_b = 0.0
-    real, parameter :: defaultSedimentTransport_c = 0.2
-    real, parameter :: defaultSedimentEnrichment_k = 1.0
-    real, parameter :: defaultSedimentEnrichment_a = 0.0
-    real, parameter :: defaultSlope = 0.0005
-    real, parameter :: defaultDepositionAlpha = 38.1        ! Zhiyao et al, 2008: https://doi.org/10.1016/S1674-2370(15)30017-X
-    real, parameter :: defaultDepositionBeta = 0.93         ! Zhiyao et al, 2008: https://doi.org/10.1016/S1674-2370(15)30017-X
-    real, parameter :: defaultBankErosionAlpha = 1.0e-9     ! [kg/m5] Loosely based on Lazar et al, 2010: https://doi.org/10.1016/j.scitotenv.2010.02.030
-    real, parameter :: defaultBankErosionBeta = 1.0         ! [-] Loosely based on Lazar et al, 2010: https://doi.org/10.1016/j.scitotenv.2010.02.030
+    real(dp), parameter :: defaultSoilAttachmentEfficiency = 0.0_dp
+    real(dp), parameter :: defaultRiverAttachmentEfficiency = 0.0_dp  ! Attachment efficiency for NM to SPM in rivers
+    real(dp), parameter :: defaultEstuaryAttachmentEfficiency = 0.0_dp
+    real(dp), parameter :: defaultSoilDarcyVelocity = 9.0e-6_dp      ! [m/s] Tufenkji et al, 2004: https://doi.org/10.1021/es034049r 
+    real(dp), parameter :: default_k_diss_pristine = 0.0_dp        ! Dissolution rate for pristine contaminant [s-1]
+    real(dp), parameter :: default_k_diss_transformed = 0.0_dp     ! Dissolution rate for transformed contaminant [s-1]
+    real(dp), parameter :: default_k_transform_pristine = 0.0_dp   ! Transformation rate for pristine contaminant [s-1]
+    real(dp), parameter :: default_rho_contaminant = 1000.0_dp     ! Density of contaminant [kg/m3]
+    real(dp), parameter :: defaultShearRate = 10.0_dp              ! Arvidsson et al, 2009: https://doi.org/10.1080/10807039.2011.538639
+    real(dp), parameter :: defaultMinWaterTemperature = 4.0_dp     ! Thames River
+    real(dp), parameter :: defaultMaxWaterTemperature = 21.0_dp    ! Thames River
+    integer, parameter  :: defaultMinWaterTemperatureDayOfYear = 32  ! Thames River
+    real(dp), parameter :: defaultSedimentTransport_a = 2.0e-9_dp
+    real(dp), parameter :: defaultSedimentTransport_b = 0.0_dp
+    real(dp), parameter :: defaultSedimentTransport_c = 0.2_dp
+    real(dp), parameter :: defaultSedimentEnrichment_k = 1.0_dp
+    real(dp), parameter :: defaultSedimentEnrichment_a = 0.0_dp
+    real(dp), parameter :: defaultSlope = 0.0005_dp
+    real(dp), parameter :: defaultDepositionAlpha = 38.1_dp        ! Zhiyao et al, 2008: https://doi.org/10.1016/S1674-2370(15)30017-X
+    real(dp), parameter :: defaultDepositionBeta = 0.93_dp         ! Zhiyao et al, 2008: https://doi.org/10.1016/S1674-2370(15)30017-X
+    real(dp), parameter :: defaultBankErosionAlpha = 1.0e-9_dp     ! [kg/m5] Loosely based on Lazar et al, 2010: https://doi.org/10.1016/j.scitotenv.2010.02.030
+    real(dp), parameter :: defaultBankErosionBeta = 1.0_dp         ! [-] Loosely based on Lazar et al, 2010: https://doi.org/10.1016/j.scitotenv.2010.02.030
     ! Default for contaminant form distribution (pristine, transformed, dissolved)
-    real, parameter :: defaultContaminantFormDistribution(3) = [1.0, 0.0, 0.0]  ! All in pristine form
+    real(dp), parameter :: defaultContaminantFormDistribution(3) = [1.0_dp, 0.0_dp, 0.0_dp]  ! All in pristine form
 end module