MODULE regr_conserv_m

  USE assert_eq_m,   ONLY: assert_eq
  USE assert_m,      ONLY: assert
  USE interpolation, ONLY: locate

  IMPLICIT NONE

! Purpose: Each procedure regrids a piecewise linear function (not necessarily
!          continuous) by averaging it. This is a conservative regridding.
!          The regridding operation is done along dimension "ix" of the input
!          field "vs". Input are positions of cell edges.
! Remarks:
!   * The target grid should be included in the source grid:
!     no extrapolation is allowed.
!   * Field on target grid "vt" has same rank, slopes and averages as "vs".
!   * If optional argument "slope" is not given, 0 is assumed for slopes.
!     Then the regridding is first order instead of second.
!   * "vs" and "vt" have the same dimensions except Nr. ix (ns for vs, nt for vt)

! Argument                         Type         Description
!-------------------------------------------------------------------------------
!  INTEGER, INTENT(IN)  :: ix      Scalar       dimension regridded <=rank(vs)
!  REAL,    INTENT(IN)  :: vs(*)   Rank>=1      averages in source grid cells
!  REAL,    INTENT(IN)  :: xs(:)   Vector(ns+1) edges of source grid, asc. order
!  REAL,    INTENT(IN)  :: xt(:)   Vector(nt+1) edges of target grid, asc. order
!  REAL,    INTENT(OUT) :: vt(*)   Rank>=1      regridded field
! [REAL,    INTENT(IN)  :: slope]  Rank>=1      slopes in source grid cells

  INTERFACE regr_conserv
    ! The procedures differ only from the rank of the input/output fields.
    MODULE PROCEDURE regr1_conserv, regr2_conserv, regr3_conserv, &
                     regr4_conserv, regr5_conserv
  END INTERFACE

  PRIVATE
  PUBLIC :: regr_conserv

CONTAINS

!-------------------------------------------------------------------------------
!
SUBROUTINE regr1_conserv(ix, vs, xs, xt, vt, slope)
!
!-------------------------------------------------------------------------------
! Arguments:
  INTEGER,        INTENT(IN)  :: ix
  REAL,           INTENT(IN)  :: vs(:)
  REAL,           INTENT(IN)  :: xs(:), xt(:)
  REAL,           INTENT(OUT) :: vt(:)
  REAL, OPTIONAL, INTENT(IN)  :: slope(:)
!-------------------------------------------------------------------------------
! Local variables:
  INTEGER :: is, it, ns, nt
  REAL    :: co, idt
  LOGICAL :: lslope
!-------------------------------------------------------------------------------
  lslope=PRESENT(slope)
  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
  is=locate(xs,xt(1))                    !--- 1<= is <= ns (no extrapolation)
  vt(:)=0.
  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
      CALL mean_lin(xt(it),xt(it+1))
    ELSE
      CALL mean_lin(xt(it),xs(is+1)); is=is+1
      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
      CALL mean_lin(xs(is),xs(is+1)); is=is+1
      END DO
      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
    END IF
    IF(xs(is+1)==xt(it+1)) is=is+1       !--- 1<=is<=ns .OR. it==nt
  END DO

CONTAINS

!-------------------------------------------------------------------------------
SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
!-------------------------------------------------------------------------------
! Arguments:
  REAL, INTENT(IN)    :: a, b
!-------------------------------------------------------------------------------
  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
    vt(it) = vt(it)+idt*(b-a)*(vs(is)+co*slope(is)/2.)
  ELSE
    vt(it) = vt(it)+idt*(b-a)* vs(is)
  END IF

END SUBROUTINE mean_lin
!-------------------------------------------------------------------------------

END SUBROUTINE regr1_conserv
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE regr2_conserv(ix, vs, xs, xt, vt, slope)
!
!-------------------------------------------------------------------------------
! Arguments:
  INTEGER,        INTENT(IN)  :: ix
  REAL,           INTENT(IN)  :: vs(:,:)
  REAL,           INTENT(IN)  :: xs(:), xt(:)
  REAL,           INTENT(OUT) :: vt(:,:)
  REAL, OPTIONAL, INTENT(IN)  :: slope(:,:)
!-------------------------------------------------------------------------------
! Local variables:
  INTEGER :: is, it, ns, nt
  REAL    :: co, idt
  LOGICAL :: lslope
!-------------------------------------------------------------------------------
  lslope=PRESENT(slope)
  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
  is=locate(xs,xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
  vt(:,:)=0.
  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
      CALL mean_lin(xt(it),xt(it+1))
    ELSE
      CALL mean_lin(xt(it),xs(is+1)); is=is+1
      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
      CALL mean_lin(xs(is),xs(is+1)); is=is+1
      END DO
      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
    END IF
    IF(xs(is+1)==xt(it+1)) is=is+1       !--- 1<=is<=ns .OR. it==nt
  END DO

CONTAINS

!-------------------------------------------------------------------------------
SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
!-------------------------------------------------------------------------------
! Arguments:
  REAL, INTENT(IN)    :: a, b
!-------------------------------------------------------------------------------
  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
    IF(ix==1) vt(it,:) = vt(it,:)+idt*(b-a)*(vs(is,:)+co*slope(is,:)/2.)
    IF(ix==2) vt(:,it) = vt(:,it)+idt*(b-a)*(vs(:,is)+co*slope(:,is)/2.)
  ELSE
    IF(ix==1) vt(it,:) = vt(it,:)+idt*(b-a)* vs(is,:)
    IF(ix==2) vt(:,it) = vt(:,it)+idt*(b-a)* vs(:,is)
  END IF

END SUBROUTINE mean_lin
!-------------------------------------------------------------------------------

END SUBROUTINE regr2_conserv
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE regr3_conserv(ix, vs, xs, xt, vt, slope)
!
!-------------------------------------------------------------------------------
! Arguments:
  INTEGER,        INTENT(IN)  :: ix
  REAL,           INTENT(IN)  :: vs(:,:,:)
  REAL,           INTENT(IN)  :: xs(:), xt(:)
  REAL,           INTENT(OUT) :: vt(:,:,:)
  REAL, OPTIONAL, INTENT(IN)  :: slope(:,:,:)
!-------------------------------------------------------------------------------
! Local variables:
  INTEGER :: is, it, ns, nt
  REAL    :: co, idt
  LOGICAL :: lslope
!-------------------------------------------------------------------------------
  lslope=PRESENT(slope)
  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
  is=locate(xs,xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
  vt(:,:,:)=0.
  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
      CALL mean_lin(xt(it),xt(it+1))
    ELSE
      CALL mean_lin(xt(it),xs(is+1)); is=is+1
      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
      CALL mean_lin(xs(is),xs(is+1)); is=is+1
      END DO
      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
    END IF
    IF(xs(is+1)==xt(it+1)) is=is+1       !--- 1<=is<=ns .OR. it==nt
  END DO

CONTAINS

!-------------------------------------------------------------------------------
SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
!-------------------------------------------------------------------------------
! Arguments:
  REAL, INTENT(IN)    :: a, b
!-------------------------------------------------------------------------------
  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
    IF(ix==1) vt(it,:,:) = vt(it,:,:)+idt*(b-a)*(vs(is,:,:)+co*slope(is,:,:)/2.)
    IF(ix==2) vt(:,it,:) = vt(:,it,:)+idt*(b-a)*(vs(:,is,:)+co*slope(:,is,:)/2.)
    IF(ix==3) vt(:,:,it) = vt(:,:,it)+idt*(b-a)*(vs(:,:,is)+co*slope(:,:,is)/2.)
  ELSE
    IF(ix==1) vt(it,:,:) = vt(it,:,:)+idt*(b-a)* vs(is,:,:)
    IF(ix==2) vt(:,it,:) = vt(:,it,:)+idt*(b-a)* vs(:,is,:)
    IF(ix==3) vt(:,:,it) = vt(:,:,it)+idt*(b-a)* vs(:,:,is)
  END IF

END SUBROUTINE mean_lin
!-------------------------------------------------------------------------------

END SUBROUTINE regr3_conserv
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE regr4_conserv(ix, vs, xs, xt, vt, slope)
!
!-------------------------------------------------------------------------------
! Arguments:
  INTEGER,        INTENT(IN)  :: ix
  REAL,           INTENT(IN)  :: vs(:,:,:,:)
  REAL,           INTENT(IN)  :: xs(:), xt(:)
  REAL,           INTENT(OUT) :: vt(:,:,:,:)
  REAL, OPTIONAL, INTENT(IN)  :: slope(:,:,:,:)
!-------------------------------------------------------------------------------
! Local variables:
  INTEGER :: is, it, ns, nt
  REAL    :: co, idt
  LOGICAL :: lslope
!-------------------------------------------------------------------------------
  lslope=PRESENT(slope)
  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
  is=locate(xs,xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
  vt(:,:,:,:)=0.
  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
      CALL mean_lin(xt(it),xt(it+1))
    ELSE
      CALL mean_lin(xt(it),xs(is+1)); is=is+1
      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
      CALL mean_lin(xs(is),xs(is+1)); is=is+1
      END DO
      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
    END IF
    IF(xs(is+1)==xt(it+1)) is=is+1       !--- 1<=is<=ns .OR. it==nt
  END DO

CONTAINS

!-------------------------------------------------------------------------------
SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
!-------------------------------------------------------------------------------
! Arguments:
  REAL, INTENT(IN) :: a, b
!-------------------------------------------------------------------------------
  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
    IF(ix==1) vt(it,:,:,:) = vt(it,:,:,:)+idt*(b-a)*(vs(is,:,:,:)+co*slope(is,:,:,:)/2.)
    IF(ix==2) vt(:,it,:,:) = vt(:,it,:,:)+idt*(b-a)*(vs(:,is,:,:)+co*slope(:,is,:,:)/2.)
    IF(ix==3) vt(:,:,it,:) = vt(:,:,it,:)+idt*(b-a)*(vs(:,:,is,:)+co*slope(:,:,is,:)/2.)
    IF(ix==4) vt(:,:,:,it) = vt(:,:,:,it)+idt*(b-a)*(vs(:,:,:,is)+co*slope(:,:,:,is)/2.)
  ELSE
    IF(ix==1) vt(it,:,:,:) = vt(it,:,:,:)+idt*(b-a)* vs(is,:,:,:)
    IF(ix==2) vt(:,it,:,:) = vt(:,it,:,:)+idt*(b-a)* vs(:,is,:,:)
    IF(ix==3) vt(:,:,it,:) = vt(:,:,it,:)+idt*(b-a)* vs(:,:,is,:)
    IF(ix==4) vt(:,:,:,it) = vt(:,:,:,it)+idt*(b-a)* vs(:,:,:,is)
  END IF

END SUBROUTINE mean_lin
!-------------------------------------------------------------------------------

END SUBROUTINE regr4_conserv
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE regr5_conserv(ix, vs, xs, xt, vt, slope)
!
!-------------------------------------------------------------------------------
! Arguments:
  INTEGER,        INTENT(IN)  :: ix
  REAL,           INTENT(IN)  :: vs(:,:,:,:,:)
  REAL,           INTENT(IN)  :: xs(:), xt(:)
  REAL,           INTENT(OUT) :: vt(:,:,:,:,:)
  REAL, OPTIONAL, INTENT(IN)  :: slope(:,:,:,:,:)
!-------------------------------------------------------------------------------
! Local variables:
  INTEGER :: is, it, ns, nt
  REAL    :: co, idt
  LOGICAL :: lslope
!-------------------------------------------------------------------------------
  lslope=PRESENT(slope)
  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
  is=locate(xs,xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
  vt(:,:,:,:,:)=0.
  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
      CALL mean_lin(xt(it),xt(it+1))
    ELSE
      CALL mean_lin(xt(it),xs(is+1)); is=is+1
      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
      CALL mean_lin(xs(is),xs(is+1)); is=is+1
      END DO
      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
    END IF
    IF(xs(is+1)==xt(it+1)) is=is+1     !--- 1<=is<=ns .OR. it==nt
  END DO

CONTAINS

!-------------------------------------------------------------------------------
SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
!-------------------------------------------------------------------------------
! Arguments:
  REAL, INTENT(IN) :: a, b
!-------------------------------------------------------------------------------
  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
    IF(ix==1) vt(it,:,:,:,:) = vt(it,:,:,:,:)+idt*(b-a)*(vs(is,:,:,:,:)+co*slope(is,:,:,:,:)/2.)
    IF(ix==2) vt(:,it,:,:,:) = vt(:,it,:,:,:)+idt*(b-a)*(vs(:,is,:,:,:)+co*slope(:,is,:,:,:)/2.)
    IF(ix==3) vt(:,:,it,:,:) = vt(:,:,it,:,:)+idt*(b-a)*(vs(:,:,is,:,:)+co*slope(:,:,is,:,:)/2.)
    IF(ix==4) vt(:,:,:,it,:) = vt(:,:,:,it,:)+idt*(b-a)*(vs(:,:,:,is,:)+co*slope(:,:,:,is,:)/2.)
    IF(ix==5) vt(:,:,:,:,it) = vt(:,:,:,:,it)+idt*(b-a)*(vs(:,:,:,:,is)+co*slope(:,:,:,:,is)/2.)
  ELSE
    IF(ix==1) vt(it,:,:,:,:) = vt(it,:,:,:,:)+idt*(b-a)* vs(is,:,:,:,:)
    IF(ix==2) vt(:,it,:,:,:) = vt(:,it,:,:,:)+idt*(b-a)* vs(:,is,:,:,:)
    IF(ix==3) vt(:,:,it,:,:) = vt(:,:,it,:,:)+idt*(b-a)* vs(:,:,is,:,:)
    IF(ix==4) vt(:,:,:,it,:) = vt(:,:,:,it,:)+idt*(b-a)* vs(:,:,:,is,:)
    IF(ix==5) vt(:,:,:,:,it) = vt(:,:,:,:,it)+idt*(b-a)* vs(:,:,:,:,is)
  END IF

END SUBROUTINE mean_lin
!-------------------------------------------------------------------------------

END SUBROUTINE regr5_conserv
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE check_size(ix,svs,svt,xs,xt,ns,nt)
!
!-------------------------------------------------------------------------------
! Arguments:
  INTEGER, INTENT(IN)  :: ix, svs(:), svt(:)
  REAL,    INTENT(IN)  :: xs(:), xt(:)
  INTEGER, INTENT(OUT) :: ns,    nt
!-------------------------------------------------------------------------------
! Local variables:
  INTEGER :: rk, is, ii, n
  CHARACTER(LEN=80) :: sub, msg
!-------------------------------------------------------------------------------
  WRITE(sub,'(a,2i0,a)')"regr",SIZE(svs),ix,"_conserv"
  rk=assert_eq(SIZE(svs),SIZE(svt),TRIM(sub)//': inconsistent ranks')
  WRITE(msg,'(a,2(i0,a))')TRIM(sub)//": ix (",ix,") exceeds field rank (",rk,")"
  CALL assert(ix>=1.AND.ix<=rk,msg)
  ii=0
  DO is=1,rk; IF(is==ix) CYCLE
    WRITE(msg,'(a,i0)')TRIM(sub)//" n",is
    n=assert_eq(svs(is),svt(is),msg)
    ii=ii+1
  END DO
  ns=assert_eq(svs(ix),SIZE(xs)-1,TRIM(sub)//" ns")
  nt=assert_eq(svt(ix),SIZE(xt)-1,TRIM(sub)//" nt")

  !--- Quick check on sort order:
  CALL assert(xs(1)<xs(2), TRIM(sub)//": xs bad order")
  CALL assert(xt(1)<xt(2), TRIM(sub)//": xt bad order")
  CALL assert(xs(1)<=xt(1).AND.xt(nt+1)<=xs(ns+1), TRIM(sub)//": extrapolation")

END SUBROUTINE check_size
!
!-------------------------------------------------------------------------------


END MODULE regr_conserv_m