--- trunk/libf/regr1_step_av.f90	2008/07/25 19:59:34	13
+++ trunk/libf/regr1_step_av.f90	2010/12/02 17:11:04	36
@@ -1,5 +1,7 @@
 module regr1_step_av_m
 
+  ! Author: Lionel GUEZ
+
   implicit none
 
   interface regr1_step_av
@@ -13,7 +15,8 @@
      ! extrapolation is allowed.
      ! The difference between the procedures is the rank of the first argument.
 
-     module procedure regr11_step_av, regr12_step_av, regr13_step_av
+     module procedure regr11_step_av, regr12_step_av, regr13_step_av, &
+          regr14_step_av
   end interface
 
   private
@@ -25,7 +28,8 @@
 
     ! "vs" has rank 1.
 
-    use numer_rec, only: assert_eq, assert, locate
+    use nr_util, only: assert_eq, assert
+    use numer_rec, only: locate
 
     real, intent(in):: vs(:) ! values of steps on the source grid
     ! (Step "is" is between "xs(is)" and "xs(is + 1)".)
@@ -82,7 +86,8 @@
 
     ! "vs" has rank 2.
 
-    use numer_rec, only: assert_eq, assert, locate
+    use nr_util, only: assert_eq, assert
+    use numer_rec, only: locate
 
     real, intent(in):: vs(:, :) ! values of steps on the source grid
     ! (Step "is" is between "xs(is)" and "xs(is + 1)".)
@@ -140,7 +145,8 @@
 
     ! "vs" has rank 3.
 
-    use numer_rec, only: assert_eq, assert, locate
+    use nr_util, only: assert_eq, assert
+    use numer_rec, only: locate
 
     real, intent(in):: vs(:, :, :) ! values of steps on the source grid
     ! (Step "is" is between "xs(is)" and "xs(is + 1)".)
@@ -193,4 +199,65 @@
 
   end function regr13_step_av
 
+  !********************************************
+
+  function regr14_step_av(vs, xs, xt) result(vt)
+
+    ! "vs" has rank 4.
+
+    use nr_util, only: assert_eq, assert
+    use numer_rec, only: locate
+
+    real, intent(in):: vs(:, :, :, :) ! values of steps on the source grid
+    ! (Step "is" is between "xs(is)" and "xs(is + 1)".)
+
+    real, intent(in):: xs(:)
+    ! (edges of steps on the source grid, in strictly increasing order)
+
+    real, intent(in):: xt(:)
+    ! (edges of cells of the target grid, in strictly increasing order)
+
+    real vt(size(xt) - 1, size(vs, 2), size(vs, 3), size(vs, 4))
+    ! (average values on the target grid)
+    ! (Cell "it" is between "xt(it)" and "xt(it + 1)".)
+
+    ! Variables local to the procedure:
+    integer is, it, ns, nt
+    real left_edge
+
+    !---------------------------------------------
+
+    ns = assert_eq(size(vs, 1), size(xs) - 1, "regr14_step_av ns")
+    nt = size(xt) - 1
+
+    ! Quick check on sort order:
+    call assert(xs(1) < xs(2), "regr14_step_av xs bad order")
+    call assert(xt(1) < xt(2), "regr14_step_av xt bad order")
+
+    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
+         "regr14_step_av extrapolation")
+
+    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
+    do it = 1, nt
+       ! 1 <= is <= ns
+       ! xs(is) <= xt(it) < xs(is + 1)
+       ! Compute "vt(it, :, :, :)":
+       left_edge = xt(it)
+       vt(it, :, :, :) = 0.
+       do while (xs(is + 1) < xt(it + 1))
+          ! 1 <= is <= ns - 1
+          vt(it, :, :, :) = vt(it, :, :, :) + (xs(is + 1) - left_edge) &
+               * vs(is, :, :, :)
+          is = is + 1
+          left_edge = xs(is)
+       end do
+       ! 1 <= is <= ns
+       vt(it, :, :, :) = (vt(it, :, :, :) + (xt(it + 1) - left_edge) &
+            * vs(is, :, :, :)) / (xt(it + 1) - xt(it))
+       if (xs(is + 1) == xt(it + 1)) is = is + 1
+       ! 1 <= is <= ns .or. it == nt
+    end do
+
+  end function regr14_step_av
+
 end module regr1_step_av_m