Merge pull request #15 from NCAR/develop-interpolate-error

add detail to interpolator errors

src/cross_section.F90

@@ -291,7 +291,9 @@ subroutine process_file( this, config, grid_warehouse, parameters )
         work_parameter(:) =                                                   &
             interpolator%interpolate( x_target = wavelength_grid%edge_,       &
                                       x_source = data_lambda,                 &
-                                      y_source = data_parameter )
+                                      y_source = data_parameter,              &
+                                      requested_by =                          &
+                                      "cross section wavelength grid" )
         parameters%array( i_low:i_high, i_param ) =                           &
                                   parameters%array( i_low:i_high, i_param ) + &
                                   work_parameter( i_low:i_high )

src/cross_sections/hno3-oh_no2.F90

@@ -133,7 +133,9 @@ function constructor( config, grid_warehouse, profile_warehouse )           &
           this%cross_section_parms( fileNdx )%array( :, parmNdx ) =         &
               interpolator%interpolate( x_target = lambdaGrid%edge_,        &
                                         x_source = data_lambda,             &
-                                        y_source = data_parameter )
+                                        y_source = data_parameter,          &
+                                        requested_by =                      &
+                               "HNO3 OH NO2 cross section wavelength grid" )
         enddo
       else
         this%cross_section_parms( fileNdx )%array = netcdf_obj%parameters

src/cross_sections/no2_tint.F90

@@ -170,7 +170,9 @@ function constructor( config, grid_warehouse, profile_warehouse )           &
           Xsection%array( :, parmNdx ) =                                    &
               interpolator%interpolate( x_target = lambdaGrid%edge_,        &
                                         x_source = data_lambda,             &
-                                        y_source = data_parameter )
+                                        y_source = data_parameter,          &
+                                        requested_by =                      &
+                "NO2 temperature integrated cross section wavelength grid" )
         enddo
       else
         Xsection%array = netcdf_obj%parameters

src/cross_sections/o3_tint.F90

@@ -195,7 +195,9 @@ function constructor( config, grid_warehouse, profile_warehouse )           &
           Xsection%array( :, parmNdx ) =                                      &
                 interpolator%interpolate( x_target = lambdaGrid%edge_,        &
                                           x_source = data_lambda,             &
-                                          y_source = data_parameter )
+                                          y_source = data_parameter,          &
+                                          requested_by =                      &
+                   "O3 temperature integrated cross section wavelength grid" )
         enddo
         deallocate( interpolator )
       else

src/cross_sections/rono2.F90

@@ -143,7 +143,9 @@ function constructor( config, grid_warehouse, profile_warehouse )           &
           this%cross_section_parms( fileNdx )%array( :, parmNdx ) =           &
               interpolator%interpolate( x_target = lambdaGrid%edge_,          &
                                         x_source = data_lambda,               &
-                                        y_source = data_parameter )
+                                        y_source = data_parameter,            &
+                                        requested_by =                        &
+                                        "RONO2 cross section wavelength grid" )
         enddo
       else
         this%cross_section_parms( fileNdx )%array = netcdf_obj%parameters

src/cross_sections/temperature_based.F90

@@ -266,7 +266,9 @@ function calculate( this, grid_warehouse, profile_warehouse, at_mid_point ) &
       cross_section( i_height, : ) =                                          &
           this%interpolator_%interpolate( x_target = wavelengths%edge_,       &
                                           x_source = this%raw_wavelengths_,   &
-                                          y_source = raw_data )
+                                          y_source = raw_data,                &
+                                          requested_by =                      &
+                           "temperature based cross section wavelength grid" )
     end do
     deallocate( temperatures )
     deallocate( wavelengths  )

src/cross_sections/tint.F90

@@ -169,7 +169,9 @@ function constructor( config, grid_warehouse, profile_warehouse )           &
           Xsection%array( :, parmNdx ) =                                    &
               interpolator%interpolate( x_target = lambdaGrid%edge_,        &
                                         x_source = data_lambda,             &
-                                        y_source = data_parameter )
+                                        y_source = data_parameter,          &
+                                        requested_by =                      &
+                    "temperature integrated cross section wavelength grid" )
         enddo
       else
         Xsection%array = netcdf_obj%parameters

src/grids/equal_delta.F90

@@ -62,15 +62,16 @@ function constructor( config ) result ( this )
     call config%get( 'units', this%units_, Iam )
 
     this%ncells_ = ( Upper_val - Lower_val ) / Delta_val
-    if( mod( ( Upper_val - Lower_val ), Delta_val ) /= 0._dk ) then
+    if( Upper_val - ( Lower_val + this%ncells_ * Delta_val )                  &
+        > spacing( Lower_val + this%ncells_ * Delta_val ) ) then
       this%ncells_ = this%ncells_ + 1
     endif
     allocate( this%mid_( this%ncells_ ) )
     allocate( this%delta_( this%ncells_ ) )
     allocate( this%edge_( this%ncells_ + 1 ) )
     do n = 1, this%ncells_ + 1
       this%edge_( n ) =                                                       &
-        min( real( ( n - 1 ), kind = dk ) * Delta_val + Lower_val, Upper_val )
+        min( Lower_val + ( n - 1 ) * Delta_val, Upper_val )
     enddo
     this%mid_(:) = .5_dk * &
       ( this%edge_( 1 : this%ncells_ ) + this%edge_( 2 : this%ncells_ + 1 ) )

src/interpolate.F90

@@ -62,18 +62,20 @@ module tuvx_interpolate
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-  function interpolate( this, x_target, x_source, y_source ) result( y_target )
+  function interpolate( this, x_target, x_source, y_source, requested_by )    &
+      result( y_target )
     ! Interpolates data in y_source on the x_source axis, to y_target on the
     ! x_target axis
 
     use musica_constants, only : dk => musica_dk
 
     import interpolator_t
 
-    class(interpolator_t), intent(in) :: this        ! Interpolator
-    real(dk),              intent(in) :: x_target(:) ! Target axis
-    real(dk),              intent(in) :: x_source(:) ! Source axis
-    real(dk),              intent(in) :: y_source(:) ! Source data
+    class(interpolator_t), intent(in) :: this         ! Interpolator
+    real(dk),              intent(in) :: x_target(:)  ! Target axis
+    real(dk),              intent(in) :: x_source(:)  ! Source axis
+    real(dk),              intent(in) :: y_source(:)  ! Source data
+    character(len=*),      intent(in) :: requested_by ! Calling function
 
     real(dk), allocatable :: y_target(:) ! Target data
 
@@ -132,8 +134,8 @@ end function constructor
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-  function interpolate_linear( this, x_target, x_source, y_source )           &
-      result( y_target )
+  function interpolate_linear( this, x_target, x_source, y_source,            &
+      requested_by ) result( y_target )
     !  Standard linear interpolation
     !
     !  Map input data given on single, discrete points, onto a discrete target
@@ -151,10 +153,11 @@ function interpolate_linear( this, x_target, x_source, y_source )           &
     !  If the input data does not encompass the target grid, use ADDPNT to
     !  expand the input array.
 
-    class(interpolator_linear_t), intent(in) :: this        ! Interpolator
-    real(dk),                     intent(in) :: x_target(:) ! Target axis
-    real(dk),                     intent(in) :: x_source(:) ! Source axis
-    real(dk),                     intent(in) :: y_source(:) ! Source data
+    class(interpolator_linear_t), intent(in) :: this         ! Interpolator
+    real(dk),                     intent(in) :: x_target(:)  ! Target axis
+    real(dk),                     intent(in) :: x_source(:)  ! Source axis
+    real(dk),                     intent(in) :: y_source(:)  ! Source data
+    character(len=*),             intent(in) :: requested_by ! Calling function
 
     real(dk), allocatable :: y_target(:)
 
@@ -191,8 +194,8 @@ end function interpolate_linear
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-  function interpolate_conserving( this, x_target, x_source, y_source )       &
-      result( y_target )
+  function interpolate_conserving( this, x_target, x_source, y_source,        &
+      requested_by ) result( y_target )
     !  Area-conserving linear interpolation
     !
     !  Map input data given on single, discrete points onto a set of target
@@ -216,10 +219,11 @@ function interpolate_conserving( this, x_target, x_source, y_source )       &
 
     use musica_assert, only : die_msg
 
-    class(interpolator_conserving_t), intent(in) :: this        ! Interpolator
-    real(dk),                         intent(in) :: x_target(:) ! Target axis
-    real(dk),                         intent(in) :: x_source(:) ! Source axis
-    real(dk),                         intent(in) :: y_source(:) ! Source data
+    class(interpolator_conserving_t), intent(in) :: this         ! Interpolator
+    real(dk),                         intent(in) :: x_target(:)  ! Target axis
+    real(dk),                         intent(in) :: x_source(:)  ! Source axis
+    real(dk),                         intent(in) :: y_source(:)  ! Source data
+    character(len=*),                 intent(in) :: requested_by ! Calling function
 
     real(dk), allocatable   :: y_target(:)
 
@@ -236,11 +240,15 @@ function interpolate_conserving( this, x_target, x_source, y_source )       &
 
     ! test for correct ordering of data, by increasing value of x
     if( any( x_source( 1 : n - 1 ) >= x_source( 2 : n ) ) ) then
-      call die_msg( 996313430,'src grid must be monotonically increasing' )
+      call die_msg( 996313430,                                                &
+                    "Source grid must be monotonically increasing for '"//    &
+                    requested_by//"' interpolation" )
     endif
 
     if( any( x_target( 1 : ng - 1 ) >= x_target( 2 : ng ) ) ) then
-      call die_msg( 208631776,'target grid must be monotonically increasing' )
+      call die_msg( 208631776,                                                &
+                    "Target grid must be monotonically increasing for '"//    &
+                    requested_by//"' interpolation" )
     endif
 
     ! check for xg-values outside the x-range
@@ -311,8 +319,8 @@ end function interpolate_conserving
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-  function interpolate_fractional_source( this, x_target, x_source, y_source )&
-      result( y_target )
+  function interpolate_fractional_source( this, x_target, x_source, y_source, &
+      requested_by ) result( y_target )
     !  Interpolation based on fractional overlap of grid sections relative
     !  to the source grid section width.
     !
@@ -340,10 +348,11 @@ function interpolate_fractional_source( this, x_target, x_source, y_source )&
 
     use musica_assert, only : die_msg
 
-    class(interpolator_fractional_source_t), intent(in) :: this        ! Interpolator
-    real(dk),                                intent(in) :: x_target(:) ! Target axis
-    real(dk),                                intent(in) :: x_source(:) ! Source axis
-    real(dk),                                intent(in) :: y_source(:) ! Source data
+    class(interpolator_fractional_source_t), intent(in) :: this         ! Interpolator
+    real(dk),                                intent(in) :: x_target(:)  ! Target axis
+    real(dk),                                intent(in) :: x_source(:)  ! Source axis
+    real(dk),                                intent(in) :: y_source(:)  ! Source data
+    character(len=*),                        intent(in) :: requested_by ! Calling function
 
     real(dk), allocatable   :: y_target(:)
 
@@ -413,8 +422,8 @@ end function interpolate_fractional_source
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-  function interpolate_fractional_target( this, x_target, x_source, y_source )&
-      result( y_target )
+  function interpolate_fractional_target( this, x_target, x_source, y_source, &
+      requested_by ) result( y_target )
     !  Interpolation based on fractional overlap of grid sections relative
     !  to the target grid section width.
     !
@@ -445,6 +454,7 @@ function interpolate_fractional_target( this, x_target, x_source, y_source )&
     real(dk),                                intent(in) :: x_target(:) ! Target axis
     real(dk),                                intent(in) :: x_source(:) ! Source axis
     real(dk),                                intent(in) :: y_source(:) ! Source data
+    character(len=*),                        intent(in) :: requested_by ! Calling function
 
     real(dk), allocatable :: y_target(:)
 

src/profiles/air.F90

@@ -140,7 +140,8 @@ function constructor( config, grid_warehouse ) result ( this )
     nData = size( zdata )
     zdata( nData ) = zdata( nData ) + .001_dk
     airlog = log( profile )
-    this%edge_val_ = theInterpolator%interpolate( zGrid%edge_, zdata, airlog )
+    this%edge_val_ = theInterpolator%interpolate( zGrid%edge_, zdata, airlog, &
+                                     this%handle_%val_//" profile height grid" )
     this%edge_val_ = exp( this%edge_val_ )
 
     this%mid_val_ = .5_dk * ( this%edge_val_( 1 : this%ncells_ ) +            &

src/profiles/extraterrestrial_flux.F90

@@ -218,7 +218,8 @@ function constructor( config, grid_warehouse ) result ( this )
 
       ! interpolate from source to model wavelength grid
       interpolatedEtfl = theInterpolator%interpolate( &
-        lambdaGrid%edge_, inputGrid, inputData )
+        lambdaGrid%edge_, inputGrid, inputData,                               &
+        this%handle_%val_//" profile height grid" )
       if( .not. allocated( this%mid_val_ ) ) then
         allocate( this%mid_val_,mold=interpolatedEtfl )
         this%mid_val_ = 0.0_dk

src/profiles/from_csv_file.F90

@@ -145,7 +145,8 @@ function constructor( config, grid_warehouse ) result( this )
           // " not a valid selection" )
     end select
 
-    this%edge_val_ = theInterpolator%interpolate( grid%edge_, zdata, profile )
+    this%edge_val_ = theInterpolator%interpolate( grid%edge_, zdata, profile, &
+                                     this%handle_%val_//" profile height grid" )
 
     this%mid_val_ = .5_dk * ( this%edge_val_( 1 : this%ncells_ ) +            &
                               this%edge_val_( 2 : this%ncells_ + 1 ) )

src/profiles/o2.F90

@@ -142,7 +142,8 @@ function constructor( config, grid_warehouse ) result( this )
     nData = size( zdata )
     zdata( nData ) = zdata( nData ) + .001_dk
     airlog = log( profile )
-    this%edge_val_ = theInterpolator%interpolate( zGrid%edge_, zdata, airlog )
+    this%edge_val_ = theInterpolator%interpolate( zGrid%edge_, zdata, airlog, &
+                                     this%handle_%val_//" profile height grid" )
     this%edge_val_ = exp( this%edge_val_ )
     this%edge_val_ = o2Vmr * this%edge_val_
 

src/profiles/o3.F90

@@ -158,7 +158,8 @@ function constructor( config, grid_warehouse ) result( this )
                                   // " not a valid selection" )
     end select
 
-    this%edge_val_ = theInterpolator%interpolate( zGrid%edge_, zdata, profile )
+    this%edge_val_ = theInterpolator%interpolate( zGrid%edge_, zdata, profile,&
+                                     this%handle_%val_//" profile height grid" )
 
     this%mid_val_ = .5_dk * ( this%edge_val_( 1 : this%ncells_ ) +            &
                               this%edge_val_( 2 : this%ncells_ + 1 ) )

src/quantum_yield.F90

@@ -210,7 +210,9 @@ subroutine base_constructor( this, config, grid_warehouse,                  &
             this%quantum_yield_parms( fileNdx )%array( :, parmNdx ) =         &
                 interpolator%interpolate( x_target = lambdaGrid%edge_,        &
                                           x_source = data_lambda,             &
-                                          y_source = data_parameter )
+                                          y_source = data_parameter,          &
+                                          requested_by =                      &
+                                          "quantum yield wavelength grid" )
           enddo
         else
           this%quantum_yield_parms( fileNdx )%array = netcdf_obj%parameters

src/quantum_yields/no2_tint.F90

@@ -156,7 +156,9 @@ function constructor( config, grid_warehouse, profile_warehouse )           &
           this%parameters( fileNdx )%array( :, parmNdx ) =                    &
                 interpolator%interpolate( x_target = lambdaGrid%edge_,        &
                                           x_source = data_lambda,             &
-                                          y_source = data_parameter )
+                                          y_source = data_parameter,          &
+                                          requested_by =                      &
+                "NO2 temperature interpolated quantum yield wavelength grid" )
         enddo
       else
         this%parameters( fileNdx )%array = netcdf_obj%parameters

src/quantum_yields/tint.F90

@@ -171,8 +171,10 @@ function constructor( config, grid_warehouse, profile_warehouse )           &
             this%parameters( fileNdx )%array( :, parmNdx ) =                  &
                 interpolator%interpolate( x_target = lambdaGrid%edge_,        &
                                           x_source = data_lambda,             &
-                                          y_source = data_parameter )
-          enddo
+                                          y_source = data_parameter,          &
+                                          requested_by =                      &
+                    "temperature interpolated quantum yield wavelength grid" )
+          enddo 
         else
           this%parameters( fileNdx )%array = netcdf_obj%parameters
         endif