!						********************
!						* module_pmat1.f90 *
!						* PURSER 1994/1999 *
!                                               *   FUJITA 1999    *
!						********************
!
!
!   Routines for basic algebraic operations on general matrices and vectors
!------------------------------------------------------------------------------
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!  These routines, perform basic algebraic operations on real vectors and
!  matrices. The task performed by each routine is, as far as possible,
!  encoded in each routine's name; three letters describe the
!  operation, the remainder defining the type of operand and, if needed to
!  resolve an ambiguity, the type of result.
!
!  OPERATIONS:
!   DET     evaluate log-determinant
!   DIF     differentiate
!   INT     integrate
!   INV     invert the matrix, or linear system involving the matrix operand
!   L1L     Cholesky LU decomposition, where U is just L-transpose
!   L1U     L-U decomposition of first arg, with 1's along diagonal of L and U
!   LDL     Cholesky LDU decomposition, where U is just L-transpose and D diag.
!   LDU     LDU decomposition
!   NOR     evaluate norm of operand
!   POL     polynomial (first argument) of second argument
!   POW     raise operand to some integer power
!   SWP     swap first two operands
!   TRC     evaluate trace of operand
!   U1L     back substitution with matrix decomposed into LU form, 1's on diag.
!   UDL     back substitution with matrix decomposed into LDU form
!   WRT     write out
!   ZER     set operand to zero
!
!  OPERAND TYPES:
!   B	    banded matrix
!   C	    circulant matrix
!   D	    diagonal matrix
!   H	    symmetric or hermitian matrix
!   L	    lower triangular matrix
!   M	    matrix (rectangular, in general)
!   P	    polynomial or power-series coefficient vector
!   Q	    sQuare matrix with Fortran dimension same as logical dimension
!   R	    row of a matrix
!   S	    scalar
!   T	    transpose of the matrix
!   U	    upper triangular matrix
!   V	    vector, or column of a matrix
!   X	    field of parallel X-vectors (aligned like "columns" of a matrix)
!   Y	    field of parallel Y-vectors (aligned like "rows" of a matrix)
!
!------------------------------------------------------------------------------
MODULE MODULE_pmat1
IMPLICIT NONE
INTERFACE pro333  ; MODULE PROCEDURE pro333;                 END INTERFACE
INTERFACE pro333_d; MODULE PROCEDURE dpro333;                END INTERFACE
INTERFACE cro33   ; MODULE PROCEDURE cro33;                  END INTERFACE
INTERFACE cro33_d;  MODULE PROCEDURE dcro33;                 END INTERFACE
INTERFACE norv;     MODULE PROCEDURE norv;                   END INTERFACE
INTERFACE norv_d;   MODULE PROCEDURE dnorv;                  END INTERFACE
INTERFACE norq;     MODULE PROCEDURE norq;                   END INTERFACE
INTERFACE norq_d;   MODULE PROCEDURE dnorq;                  END INTERFACE
INTERFACE swpvv;    MODULE PROCEDURE swpvv;                  END INTERFACE
INTERFACE swpvv_d;  MODULE PROCEDURE dswpvv;                 END INTERFACE
INTERFACE mulmd;    MODULE PROCEDURE mulmd;                  END INTERFACE
INTERFACE mulmd_d;  MODULE PROCEDURE dmulmd;                 END INTERFACE
INTERFACE multd;    MODULE PROCEDURE multd;                  END INTERFACE
INTERFACE multd_d;  MODULE PROCEDURE dmultd;                 END INTERFACE
INTERFACE muldm;    MODULE PROCEDURE muldm;                  END INTERFACE
INTERFACE muldm_d;  MODULE PROCEDURE dmuldm;                 END INTERFACE
INTERFACE muldt;    MODULE PROCEDURE muldt;                  END INTERFACE
INTERFACE muldt_d;  MODULE PROCEDURE dmuldt;                 END INTERFACE
INTERFACE mulpp;    MODULE PROCEDURE mulpp;                  END INTERFACE
INTERFACE mulpp_d;  MODULE PROCEDURE dmulpp;                 END INTERFACE
INTERFACE madpp;    MODULE PROCEDURE madpp;                  END INTERFACE
INTERFACE madpp_d;  MODULE PROCEDURE dmadpp;                 END INTERFACE
INTERFACE msbpp;    MODULE PROCEDURE msbpp;                  END INTERFACE
INTERFACE msbpp_d;  MODULE PROCEDURE dmsbpp;                 END INTERFACE
INTERFACE difp;     MODULE PROCEDURE difp;                   END INTERFACE
INTERFACE difp_d;   MODULE PROCEDURE ddifp;                  END INTERFACE
INTERFACE intp;     MODULE PROCEDURE intp;                   END INTERFACE
INTERFACE intp_d;   MODULE PROCEDURE dintp;                  END INTERFACE
INTERFACE invp;     MODULE PROCEDURE invp;                   END INTERFACE
INTERFACE invp_d;   MODULE PROCEDURE dinvp;                  END INTERFACE
INTERFACE prgv;     MODULE PROCEDURE prgv;                   END INTERFACE
INTERFACE prgv_d;   MODULE PROCEDURE dprgv;                  END INTERFACE
INTERFACE mulcc;    MODULE PROCEDURE mulcc;                  END INTERFACE
INTERFACE mulcc_d;  MODULE PROCEDURE dmulcc;                 END INTERFACE
INTERFACE madcc;    MODULE PROCEDURE madcc;                  END INTERFACE
INTERFACE madcc_d;  MODULE PROCEDURE dmadcc;                 END INTERFACE
INTERFACE msbcc;    MODULE PROCEDURE msbcc;                  END INTERFACE
INTERFACE msbcc_d;  MODULE PROCEDURE dmsbcc;                 END INTERFACE
INTERFACE zerl;     MODULE PROCEDURE zerl;                   END INTERFACE
INTERFACE zerl_d;   MODULE PROCEDURE dzerl;                  END INTERFACE
INTERFACE zeru;     MODULE PROCEDURE zeru;                   END INTERFACE
INTERFACE zeru_d;   MODULE PROCEDURE dzeru;                  END INTERFACE
INTERFACE ldum;     MODULE PROCEDURE ldum;                   END INTERFACE
INTERFACE ldum_d;   MODULE PROCEDURE dldum;                  END INTERFACE
INTERFACE udlmm;    MODULE PROCEDURE udlmm, udlmv;           END INTERFACE
INTERFACE udlmm_d;  MODULE PROCEDURE dudlmm,dudlmv;          END INTERFACE
INTERFACE linvan;   MODULE PROCEDURE linvan;                 END INTERFACE
INTERFACE linvan_d; MODULE PROCEDURE dlinvan;                END INTERFACE
INTERFACE copdm;    MODULE PROCEDURE copdm;                  END INTERFACE
INTERFACE copdm_d;  MODULE PROCEDURE dcopdm;                 END INTERFACE
INTERFACE condm;    MODULE PROCEDURE condm;                  END INTERFACE
INTERFACE condm_d;  MODULE PROCEDURE dcondm;                 END INTERFACE
INTERFACE copsm;    MODULE PROCEDURE copsm;                  END INTERFACE
INTERFACE copsm_d;  MODULE PROCEDURE dcopsm;                 END INTERFACE
INTERFACE consm;    MODULE PROCEDURE consm;                  END INTERFACE
INTERFACE consm_d;  MODULE PROCEDURE dconsm;                 END INTERFACE
INTERFACE addmd;    MODULE PROCEDURE addmd;                  END INTERFACE
INTERFACE addmd_d;  MODULE PROCEDURE daddmd;                 END INTERFACE
INTERFACE submd;    MODULE PROCEDURE submd;                  END INTERFACE
INTERFACE submd_d;  MODULE PROCEDURE dsubmd;                 END INTERFACE
INTERFACE addms;    MODULE PROCEDURE addms;                  END INTERFACE
INTERFACE addms_d;  MODULE PROCEDURE daddms;                 END INTERFACE
INTERFACE subms;    MODULE PROCEDURE subms;                  END INTERFACE
INTERFACE subms_d;  MODULE PROCEDURE dsubms;                 END INTERFACE
INTERFACE l1lm;     MODULE PROCEDURE l1lm;                   END INTERFACE
INTERFACE l1lm_d;   MODULE PROCEDURE dl1lm;                  END INTERFACE
INTERFACE ldlm;     MODULE PROCEDURE ldlm;                   END INTERFACE
INTERFACE ldlm_d;   MODULE PROCEDURE dldlm;                  END INTERFACE
INTERFACE invh;     MODULE PROCEDURE invh;                   END INTERFACE
INTERFACE invh_d;   MODULE PROCEDURE dinvh;                  END INTERFACE
INTERFACE invl;     MODULE PROCEDURE invl;                   END INTERFACE
INTERFACE invl_d;   MODULE PROCEDURE dinvl;                  END INTERFACE
INTERFACE linlv;    MODULE PROCEDURE linlv;                  END INTERFACE
INTERFACE linlv_d;  MODULE PROCEDURE dlinlv;                 END INTERFACE
INTERFACE linuv;    MODULE PROCEDURE linuv;                  END INTERFACE
INTERFACE linuv_d;  MODULE PROCEDURE dlinuv;                 END INTERFACE
INTERFACE powp;     MODULE PROCEDURE powp;                   END INTERFACE
INTERFACE powp_d;   MODULE PROCEDURE dpowp;                  END INTERFACE
INTERFACE polps;    MODULE PROCEDURE polps;                  END INTERFACE
INTERFACE polps_d;  MODULE PROCEDURE dpolps;                 END INTERFACE
INTERFACE polpp;    MODULE PROCEDURE polpp;                  END INTERFACE
INTERFACE polpp_d;  MODULE PROCEDURE dpolpp;                 END INTERFACE
INTERFACE trcm;     MODULE PROCEDURE trcm;                   END INTERFACE
INTERFACE trcm_d;   MODULE PROCEDURE dtrcm;                  END INTERFACE
INTERFACE inv;      MODULE PROCEDURE invmt, linmmt, linmvt;  END INTERFACE
INTERFACE inv_d;    MODULE PROCEDURE dinvmt,dlinmmt,dlinmvt; END INTERFACE

CONTAINS

FUNCTION pro333(d,e,f) RESULT(pro_res) ! TRIPLE PRODUCT OF 3 3-VECTORS
REAL(8)                :: pro_res
REAL(8),    INTENT(IN) :: d(3), e(3), f(3)
REAL(8)                :: g(3)
CALL CRO33(E,F,G)
pro_res=DOT_PRODUCT(d,g)
END FUNCTION pro333

FUNCTION dpro333(d,e,f) RESULT(pro_res) ! TRIPLE PRODUCT OF 3 3-VECTORS
REAL(8)             :: pro_res
REAL(8), INTENT(IN) :: d(3), e(3), f(3)
REAL(8)             :: g(3)
CALL CRO33_d(E,F,G)
pro_res=DOT_PRODUCT(d,g)
END FUNCTION dpro333

SUBROUTINE cro33(a,b,c) ! SPECIAL CASE OF 3-DIMENSIONS: CROSS-PRODUCT
REAL(8),    INTENT(IN) :: a(3), b(3)
REAL(8),    INTENT(OUT):: c(3)
c(1)=a(2)*b(3)-a(3)*b(2)
c(2)=a(3)*b(1)-a(1)*b(3)
c(3)=a(1)*b(2)-a(2)*b(1)
END SUBROUTINE cro33

SUBROUTINE dcro33(a,b,c) ! SPECIAL CASE OF 3-DIMENSIONS: CROSS-PRODUCT
REAL(8), INTENT(IN) :: a(3), b(3)
REAL(8), INTENT(OUT):: c(3)
c(1)=a(2)*b(3)-a(3)*b(2)
c(2)=a(3)*b(1)-a(1)*b(3)
c(3)=a(1)*b(2)-a(2)*b(1)
END SUBROUTINE dcro33


FUNCTION norv(d) RESULT(norv_res)! NORM OF VECTOR..
REAL(8)                :: norv_res
REAL(8),    INTENT(IN) :: d(:)
norv_res=SQRT(DOT_PRODUCT(D,D))
END FUNCTION norv

FUNCTION dnorv(d) ! NORM OF VECTOR..
REAL(8):: dnorv
REAL(8),    INTENT(IN) :: d(:)
dnorv=SQRT(DOT_PRODUCT(d,d))
END FUNCTION dnorv

FUNCTION norq(d) ! Norm of a matrix
REAL(8):: norq
REAL(8),INTENT(IN):: d(:,:)
INTEGER m2,i2
m2=SIZE(d,2)
norq=0.; DO i2=1,m2; norq=norq+dot_PRODUCT(d(:,i2),d(:,i2)); ENDDO
norq=SQRT(norq)
END FUNCTION norq

FUNCTION dnorq(d) ! norm of a matrix
REAL(8):: dnorq
REAL(8),INTENT(IN):: d(:,:)
INTEGER m2,i2
m2=SIZE(d,2)
dnorq=0.; DO i2=1,m2; dnorq=dnorq+dot_PRODUCT(d(:,i2),d(:,i2)); ENDDO
dnorq=SQRT(dnorq)
END FUNCTION dnorq


SUBROUTINE swpvv(d,e)
REAL(8), INTENT(INOUT) :: d(:), e(:)
REAL(8) :: t(SIZE(d))
t = d; d = e; e = t
END SUBROUTINE swpvv

SUBROUTINE dswpvv(d,e)
REAL(8), INTENT(INOUT) :: d(:), e(:)
REAL(8) :: t(SIZE(d))
t = d; d = e; e = t
END SUBROUTINE dswpvv

SUBROUTINE mulmd(a,d,b)
REAL(8), INTENT(INOUT) :: a(:,:),b(:,:) 
REAL(8), INTENT(IN)    :: d(*)
INTEGER:: m2,j
m2=SIZE(a,2)
DO j=1,m2; b(:,j)=a(:,j)*d(j); ENDDO
END SUBROUTINE mulmd

SUBROUTINE dmulmd(a,d,b)
REAL(8), INTENT(INOUT) :: a(:,:),b(:,:)
REAL(8), INTENT(IN)    :: d(*)
INTEGER:: m2,j
m2=SIZE(a,2)
DO j=1,m2; b(:,j)=a(:,j)*d(j); ENDDO
END SUBROUTINE dmulmd

SUBROUTINE multd(a,d,b)
REAL(8), INTENT(INOUT)    :: a(:,:),b(:,:) 
REAL(8), INTENT(IN)       :: d(*)
INTEGER:: m2,j
m2=SIZE(a,1)
DO j=1,m2; b(:,j) = a(j,:) * d(j); ENDDO
END SUBROUTINE multd

SUBROUTINE dmultd(a,d,b)
REAL(8), INTENT(INOUT) :: a(:,:),b(:,:) 
REAL(8), INTENT(IN)    :: d(*)
INTEGER:: m2,j
m2=SIZE(a,1)
DO j=1,m2; b(:,j) = a(j,:) * d(j); ENDDO
END SUBROUTINE dmultd

SUBROUTINE muldm(d,a,b)
REAL(8), INTENT(INOUT)    :: a(:,:),b(:,:) 
REAL(8), INTENT(IN)       :: d(*)
INTEGER                :: m1,i
m1=SIZE(a,1)
DO i=1,m1; b(i,:) = d(i)*a(i,:); ENDDO
END SUBROUTINE muldm

SUBROUTINE dmuldm(d,a,b)
REAL(8), INTENT(INOUT) :: a(:,:),b(:,:) 
REAL(8), INTENT(IN)    :: d(*)
INTEGER                :: m1,i
m1=SIZE(a,1)
DO i=1,m1; b(i,:) = d(i)*a(i,:); ENDDO
END SUBROUTINE dmuldm

SUBROUTINE muldt(d,a,b)
REAL(8), INTENT(INOUT)    :: a(:,:),b(:,:) 
REAL(8), INTENT(IN)       :: d(*)
INTEGER                :: m1,i
m1=SIZE(a,2)
DO i=1,m1; b(i,:) = d(i)*a(:,i); ENDDO
END SUBROUTINE muldt

SUBROUTINE dmuldt(d,a,b)
REAL(8), INTENT(INOUT) :: a(:,:),b(:,:) 
REAL(8), INTENT(IN)    :: d(*)
INTEGER:: m1,i
m1=SIZE(a,2)
DO i=1,m1; b(i,:) = d(i)*a(:,i); ENDDO
END SUBROUTINE dmuldt

SUBROUTINE mulpp(a,b,c) !  multiply polynomials, possibly in place
REAL(8),    INTENT(IN)    :: a(0:), b(0:)
REAL(8),    INTENT(INOUT) :: c(0:)
INTEGER                :: m,mcp, i, j
REAL(8)                   :: s, b0
m=SIZE(a)-1
mcp=mcmax(a,b,m)
c(mcp:m) = 0.0
DO j=mcp,1,-1
  s = SUM(a(j-1:0:-1)*b(0:j-1))
  c(j-1)=s
ENDDO
RETURN
ENTRY madpp(a,b,c)
m=SIZE(a)-1
mcp=mcmax(a,b,m)
DO j=mcp,1,-1
  s = SUM(a(j-1:0:-1)*b(0:j-1))
  c(j-1)=c(j-1)+s
ENDDO
RETURN
ENTRY msbpp(a,b,c)
m=SIZE(a)-1
mcp=mcmax(a,b,m)
DO j=mcp,1,-1
  s = SUM(a(j-1:0:-1)*b(0:j-1))
  c(j-1)=c(j-1)-s
ENDDO
RETURN
CONTAINS
FUNCTION mcmax(a,b,m) RESULT(mmx_res) ! This fn can be contained in mulpp().
INTEGER             :: mmx_res
INTEGER, INTENT(IN) :: m
REAL(8),    INTENT(IN) :: a(0:m), b(0:m)
INTEGER             :: ma, mb
mmx_res=0		       ! default for when ALL elements of c are zero
DO ma=m,0,-1	               ! seek last nonzero coefficient of polynomial a
  IF(a(ma) /= 0.)THEN
    DO mb=m,0,-1	       ! seek last nonzero coefficient of polynomial b
      IF(b(mb) /= 0.)THEN
        mmx_res=MIN(m,ma+mb)+1 ! hence, 1+last non-0 element of their product
        RETURN
      ENDIF
    ENDDO
    RETURN
  ENDIF
ENDDO
END FUNCTION mcmax
END SUBROUTINE mulpp

SUBROUTINE difp(a,b) ! Symbolically differentiate polynomial
REAL(8), INTENT(IN)  :: a(0:)
REAL(8), INTENT(OUT) :: b(0:)
INTEGER           :: m,mcp, i, j
REAL(8)              :: s, b0
m=SIZE(a)-1
DO i=1,m	! possibly with coincident storage for a and b
  b(i-1)=i*a(i)
ENDDO
b(m)=0.
RETURN
ENTRY intp(a,b) ! Symbolically integrate polynomial
m=SIZE(a)-1
DO i=m,1,-1	! possibly with coincident storage for a and b
  b(i)=a(i-1)/i
ENDDO
b(0)=0.
RETURN
ENTRY invp(a,b) ! Invert polynomial or power-series
m=SIZE(a)-1
b0=1./a(0)	! storage of a and b must not be the same
b(0)=b0
DO i=1,m
  s = SUM(b(i-1:0:-1)*a(1:i))
  b(i)=-b0*s
ENDDO
END SUBROUTINE difp

SUBROUTINE dmulpp(a,b,c) !  multiply polynomials, possibly in place
REAL(8), INTENT(IN)   :: a(0:), b(0:)
REAL(8), INTENT(INOUT):: c(0:)
INTEGER               :: m,mcp, i, j
REAL(8)               :: s, b0
m=SIZE(a)-1
mcp=mcmax(a,b,m)
c(mcp:m) = 0.0
DO j=mcp,1,-1
  s = SUM(a(j-1:0:-1)*b(0:j-1))
  c(j-1)=s
ENDDO
RETURN
ENTRY dmadpp(a,b,c)
m=SIZE(a)-1
mcp=mcmax(a,b,m)
DO j=mcp,1,-1
  s = SUM(a(j-1:0:-1)*b(0:j-1))
  c(j-1)=c(j-1)+s
ENDDO
RETURN
ENTRY dmsbpp(a,b,c)
m=SIZE(a)-1
mcp=mcmax(a,b,m)
DO j=mcp,1,-1
  s = SUM(a(j-1:0:-1)*b(0:j-1))
  c(j-1)=c(j-1)-s
ENDDO
RETURN
CONTAINS
FUNCTION mcmax(a,b,m) RESULT(mmx_res)
INTEGER              :: mmx_res
INTEGER,  INTENT(IN) :: m
REAL(8), INTENT(IN)  :: a(0:m), b(0:m)
INTEGER              :: ma, mb
mmx_res=0		       ! default for when all elements of c are zero
DO ma=m,0,-1	               ! seek last nonzero coefficient of polynomial a
  IF(a(ma) /= 0.d0)THEN
    DO mb=m,0,-1	       ! seek last nonzero coefficient of polynomial b
      IF(b(mb) /= 0.d0)THEN
        mmx_res=MIN(m,ma+mb)+1 ! hence, 1+last non-0 element of their product
        RETURN
      ENDIF
    ENDDO
    RETURN
  ENDIF
ENDDO
RETURN
END FUNCTION mcmax

END SUBROUTINE dmulpp

SUBROUTINE ddifp(a,b) ! Symbolically differentiate polynomial
REAL(8), INTENT(IN)   :: a(0:)
REAL(8), INTENT(INOUT):: b(0:)
INTEGER               :: m,mcp, i, j
REAL(8)               :: s, b0
m=SIZE(a)-1
DO i=1,m	 ! possibly with coincident storage for a and b
  b(i-1)=i*a(i)
ENDDO
b(m)=0.
RETURN
ENTRY dintp(a,b) ! Symbolically integrate polynomial
m=SIZE(a)-1
DO i=m,1,-1	 ! possibly with coincident storage for a and b
  b(i)=a(i-1)/i
ENDDO
b(0)=0.
RETURN
ENTRY dinvp(a,b) ! Invert polynomial or power-series
m=SIZE(a)-1
b0=1./a(0)	 ! storage of a and b must not be the same
b(0)=b0
DO i=1,m
  s = SUM(b(i-1:0:-1)*a(1:i))
  b(i)=-b0*s
ENDDO
END SUBROUTINE ddifp


SUBROUTINE prgv(d)
REAL(8), PARAMETER        :: crit=1.E-30
REAL(8), INTENT(INOUT)    :: d(:)
INTEGER                :: i,m
m=SIZE(d)
DO i=1,m; IF(ABS(d(i)) <= crit)d(i)=0.; ENDDO
END SUBROUTINE prgv

SUBROUTINE dprgv(d)
REAL(8), PARAMETER     :: crit=1.D-30
REAL(8), INTENT(INOUT) :: d(:)
INTEGER                :: i,m
m=SIZE(d)
DO i=1,m; IF(ABS(d(i)) <= crit)d(i)=0.; ENDDO
END SUBROUTINE dprgv


SUBROUTINE mulcc(a,b,c,m)  ! Multiply circulant matrices of period M
INTEGER, INTENT(IN) :: m
REAL(8), INTENT(INOUT) :: a(0:m-1), b(0:m-1), c(0:m-1)
INTEGER             :: mm, j
c(0:m-1) = 0.0
ENTRY madcc(a,b,c,m)
mm=m-1
DO j=0,mm
  c(j:m-1) = c(j:m-1) + a(0:m-j-1)*b(j)
  c(0:j-1) = c(0:j-1) + a(m-j:m-1)*b(j)
ENDDO
RETURN
ENTRY msbcc(a,b,c,m)
mm=m-1
DO j=0,mm
  c(j:m-1) = c(j:m-1) - a(0:m-j-1)*b(j)
  c(0:j-1) = c(0:j-1) - a(m-j:m-1)*b(j)
ENDDO
END SUBROUTINE mulcc

SUBROUTINE dmulcc(a,b,c,m)  ! Multiply circulant matrices of period M
INTEGER, INTENT(IN   ) :: m
REAL(8), INTENT(INOUT) :: a(0:m-1), b(0:m-1), c(0:m-1)
INTEGER                :: mm, j
c(0:m-1) = 0.0d0
ENTRY dmadcc(a,b,c,m)
mm=m-1
DO j=0,mm
  c(j:m-1) = c(j:m-1) + a(0:m-j-1)*b(j)
  c(0:j-1) = c(0:j-1) + a(m-j:m-1)*b(j)
ENDDO
RETURN
ENTRY dmsbcc(a,b,c,m)
mm=m-1
DO j=0,mm
  c(j:m-1) = c(j:m-1) - a(0:m-j-1)*b(j)
  c(0:j-1) = c(0:j-1) - a(m-j:m-1)*b(j)
ENDDO
END SUBROUTINE dmulcc

SUBROUTINE zerl(a)  ! Zero out the strictly lower triangle of elements
REAL(8),INTENT(INOUT):: a(:,:)
INTEGER           :: m,j
m=SIZE(a,1); DO j=1,m; a(j+1:m,j) = 0; ENDDO; RETURN

ENTRY zeru(a)       ! Zero out the strictly upper triangle of elements
m=SIZE(a,1); DO j=1,m; a(1:j-1,j) = 0; ENDDO
END SUBROUTINE zerl

SUBROUTINE dzerl(a) ! Zero out the strictly lower triangle of elements
REAL(8),INTENT(INOUT):: a(:,:)
INTEGER              :: m,j
m=SIZE(a,1); DO j=1,m; a(j+1:m,j) = 0; ENDDO; RETURN

ENTRY dzeru(a)      ! Zero out the strictly upper triangle of elements
m=SIZE(a,1); DO j=1,m; a(1:j-1,j) = 0; ENDDO
END SUBROUTINE dzerl

!------------------------------------------------------------------------------
!   R.J.Purser, NCEP, Washington D.C.	1996
!		    SUBROUTINE	LDUM
!  perform l-d-u decomposition of square matrix a in place with
!
!  <-> a    square matrix to be factorized
!  <-- ipiv array encoding the pivoting sequence
!  <-- d    indicator for possible sign change of determinant
!------------------------------------------------------------------------------
SUBROUTINE ldum(a,ipiv,d)
REAL(8),    INTENT(INOUT) :: a(:,:) 
REAL(8),    INTENT(OUT  ) :: d
INTEGER, INTENT(OUT  ) :: ipiv(:)
INTEGER                :: m,i, j, jp, ibig, jm
REAL(8)                   :: s(SIZE(a,1)),  aam, aa, abig,  ajj, ajji, aij
m=SIZE(a,1)
DO i=1,m
  aam=0.
  DO j=1,m
    aa=ABS(a(i,j))
    IF(aa > aam)aam=aa
  ENDDO
  IF(aam == 0.)THEN
    PRINT '(" row ",i3," of matrix in ldum vanishes")',i
    STOP
  ENDIF
  s(i)=1./aam
ENDDO
d=1.
ipiv(m)=m
DO j=1,m-1
  jp=j+1
  abig=s(j)*ABS(a(j,j))
  ibig=j
  DO i=jp,m
    aa=s(i)*ABS(a(i,j))
    IF(aa > abig)THEN
      ibig=i
      abig=aa
    ENDIF
  ENDDO
!  swap rows, recording changed sign of determinant
  ipiv(j)=ibig
  IF(ibig /= j)THEN
    d=-d
    CALL swpvv(a(j,:),a(ibig,:))
    s(ibig)=s(j)
  ENDIF
  ajj=a(j,j)
  IF(ajj == 0.)THEN
    jm=j-1
    PRINT '(" failure in ldum:"/" matrix singular, rank=",i3)',jm
    STOP
  ENDIF
  ajji=1./ajj
  DO i=jp,m
    aij=ajji*a(i,j)
    a(i,j)=aij
    a(i,jp:m) = a(i,jp:m) - aij*a(j,jp:m)
  ENDDO
ENDDO
END SUBROUTINE ldum
SUBROUTINE DLDUM(A,IPIV,D)
REAL(8), INTENT(INOUT) :: a(:,:) 
REAL(8), INTENT(OUT  ) :: d
INTEGER, INTENT(OUT  ) :: ipiv(:)
INTEGER                :: m,i, j, jp, ibig, jm
REAL(8)                :: s(SIZE(a,1)),  aam, aa, abig,  ajj, ajji, aij
m=SIZE(a,1)
DO i=1,m
  aam=0.
  DO j=1,m
    aa=ABS(a(i,j))
    IF(aa > aam)aam=aa
  ENDDO
  IF(aam == 0.d0)THEN
    PRINT '(" row ",i3," of matrix in dldum vanishes")',i
    STOP
  ENDIF
  s(i)=1./aam
ENDDO
d=1.
ipiv(m)=m
DO j=1,m-1
  jp=j+1
  abig=s(j)*ABS(a(j,j))
  ibig=j
  DO i=jp,m
    aa=s(i)*ABS(a(i,j))
    IF(aa > abig)THEN
      ibig=i
      abig=aa
    ENDIF
  ENDDO
!  swap rows, recording changed sign of determinant
  ipiv(j)=ibig
  IF(ibig /= j)THEN
    d=-d
    CALL swpvv_d(a(j,:),a(ibig,:))
    s(ibig)=s(j)
  ENDIF
  ajj=a(j,j)
  IF(ajj == 0.d0)THEN
    jm=j-1
    PRINT '(" Failure in dldum:"/" matrix singular, rank=",i3)',jm
    STOP
  ENDIF
  ajji=1./ajj
  DO i=jp,m
    aij=ajji*a(i,j)
    a(i,j)=aij
    a(i,jp:m) = a(i,jp:m) - aij*a(j,jp:m)
  ENDDO
ENDDO
END SUBROUTINE dldum


!------------------------------------------------------------------------------
!   R.J.Purser, National Meteorological Center, Washington D.C.  1993
!		    SUBROUTINE UDLMM
!  use l-u factors in A to back-substitute for mm rhs in B, using ipiv to
!  define the pivoting permutation used in the l-u decomposition.
!
!  --> A    L-D-U factorization of linear system matrux
!  <-> B    right-hand-sides on entry, corresponding matrix of solution
!	    vectors on return
!  --> IPIV array encoding the pivoting sequence
!------------------------------------------------------------------------------
SUBROUTINE udlmm(a,b,ipiv)
INTEGER, INTENT(IN)    :: ipiv(:) 
REAL(8),    INTENT(IN)    :: a(:,:) 
REAL(8),    INTENT(INOUT) :: b(:,:) 
INTEGER                :: m,mm,i, k, l
REAL(8)                   :: s,aiii
m=SIZE(a,1); mm=SIZE(b,2)
DO k=1,mm !loop over columns of b
  DO i=1,m
    l=ipiv(i)
    s=b(l,k)
    b(l,k)=b(i,k)
    s = s - SUM(b(1:i-1,k)*a(i,1:i-1))
    b(i,k)=s
  ENDDO
  b(m,k)=b(m,k)/a(m,m)
  DO i=m-1,1,-1
    aiii=1./a(i,i)
    b(i,k) = b(i,k) - SUM(b(i+1:m,k)*a(i,i+1:m))
    b(i,k)=b(i,k)*aiii
  ENDDO
ENDDO
END SUBROUTINE udlmm
SUBROUTINE dudlmm(a,b,ipiv)
INTEGER, INTENT(IN   ) :: ipiv(:) 
REAL(8), INTENT(IN   ) :: a(:,:) 
REAL(8), INTENT(INOUT) :: b(:,:) 
INTEGER                :: m,mm,i, k, l
REAL(8)                :: s,aiii
m=SIZE(a,1); mm=SIZE(b,2)
DO k=1,mm !loop over columns of b
  DO i=1,m
    l=ipiv(i)
    s=b(l,k)
    b(l,k)=b(i,k)
    s = s - SUM(b(1:i-1,k)*a(i,1:i-1))
    b(i,k)=s
  ENDDO
  b(m,k)=b(m,k)/a(m,m)
  DO i=m-1,1,-1
    aiii=1./a(i,i)
    b(i,k) = b(i,k) - SUM(b(i+1:m,k)*a(i,i+1:m))
    b(i,k)=b(i,k)*aiii
  ENDDO
ENDDO
END SUBROUTINE dudlmm

!------------------------------------------------------------------------------
!   R.J.Purser, National Meteorological Center, Washington D.C.  1993
!		    SUBROUTINE UDLMV
!  use l-u factors in A to back-substitute for mm rhs in B, using ipiv to
!  define the pivoting permutation used in the l-u decomposition.
!
!  --> A    L-D-U factorization of linear system matrux
!  <-> B    right-hand-side on entry, corresponding vector solution
!	    on return
!  --> IPIV array encoding the pivoting sequence
!------------------------------------------------------------------------------
SUBROUTINE udlmv(a,b,ipiv)
INTEGER, INTENT(IN)    :: ipiv(:) 
REAL(8),    INTENT(IN)    :: a(:,:) 
REAL(8),    INTENT(INOUT) :: b(:) 
INTEGER                :: m,i, l
REAL(8)                   :: s,aiii
m=SIZE(a,1)
DO i=1,m
   l=ipiv(i)
   s=b(l)
   b(l)=b(i)
   s = s - SUM(b(1:i-1)*a(i,1:i-1))
   b(i)=s
ENDDO
b(m)=b(m)/a(m,m)
DO i=m-1,1,-1
   aiii=1./a(i,i)
   b(i) = b(i) - SUM(b(i+1:m)*a(i,i+1:m))
   b(i)=b(i)*aiii
ENDDO
END SUBROUTINE udlmv
SUBROUTINE dudlmv(a,b,ipiv)
INTEGER,  INTENT(IN   ) :: ipiv(:) 
REAL(8),  INTENT(IN   ) :: a(:,:) 
REAL(8),  INTENT(INOUT) :: b(:) 
INTEGER                 :: m,i, l
REAL(8)                 :: s,aiii
m=SIZE(a,1)
DO i=1,m
   l=ipiv(i)
   s=b(l)
   b(l)=b(i)
   s = s - SUM(b(1:i-1)*a(i,1:i-1))
   b(i)=s
ENDDO
b(m)=b(m)/a(m,m)
DO i=m-1,1,-1
   aiii=1./a(i,i)
   b(i) = b(i) - SUM(b(i+1:m)*a(i,i+1:m))
   b(i)=b(i)*aiii
ENDDO
END SUBROUTINE dudlmv

!------------------------------------------------------------------------------
!   R.J.Purser, National Centers for Environmental Prediction, Washington D.C.
!		    SUBROUTINE LINVAN
!
!   Take square matrix W and seek row and column scalings to produce non-
!   vanishing elements of rescaled W having magnitudes as close to unity
!   as possible. The approach is make the geometric mean of the nonvanishing
!   elements of each row and of each column +1 or -1. Having rescaled the
!   matrix and the r.h.s. vector AB, compute the product P of row-vector
!   norms, then compute the determinant D and solve the linear system.
!   Rescale the solution vector (now AB) and put the conditioning indicator
!   formed by the ratio D/P into the first element of W.
!
! <-> W:    Generalized Vandermonde matrix in, conditioning indicator out.
! <-> AB:   R.h.s. vector in, solution vector of numerical coefficients out.
!------------------------------------------------------------------------------
SUBROUTINE linvan(w,ab)
INTEGER, PARAMETER  :: nit=20
REAL(8), INTENT(INOUT) :: w(:,:), ab(:)
REAL(8)                :: d1(SIZE(w,1)), d2(SIZE(w,1)), &
                       w2(SIZE(w,1),SIZE(w,1)),v(SIZE(w,1))
INTEGER             :: i, j, it, jt, ipiv(SIZE(w,1)), nc
REAL(8)                :: p, e, dw, c, d, s, d2j
REAL(8),ALLOCATABLE    :: wv(:,:) ! work variable for ab(nc) and v(nn)

nc = SIZE(w,DIM=1)
ALLOCATE(wv(nc,1))

w2=w                ! Preserve original W and AB for use
v = ab(1:nc)	    ! in later "clean-up" operation.

d1 = 1.0 	    ! Row scaling factors set to default
d2 = 1.0 	    ! Column scaling factors set to default

C=1.E-16	    ! Set initial criterion for "negligible" elements of W

! In first attempt to estimate row and column scalings, use logarithms
! to avoid the risk of under- or over-flows of the line products of W:
DO i=1,nc
  p=0.
  e=0.
  DO j=1,nc
    dw=ABS(w(i,j))
    IF(dw > c)THEN
      e=e+1.
      p=p+LOG(dw)
    ENDIF
  ENDDO
  IF(E == 0.)STOP 'W effectively singular in LINVAN'
  d1(i)=EXP(-p/e)
ENDDO
CALL muldm(d1,w2,w)

DO j=1,nc
  p=0.
  e=0.
  DO i=1,nc
    dw=ABS(w(i,j))
    IF(dw > c)THEN
      e=e+1.
      p=p+LOG(dw)
    ENDIF
  ENDDO
  IF(E == 0.)STOP 'W effectively singular in LINVAN'
  d2(j)=EXP(-p/e)
ENDDO
CALL mulmd(w,d2,w)

c=1.e-8  ! reset the criterion for "negligible" elements

! revert to iterations of the more efficient method without logarithms:
DO jt=1,2
DO it=1,nit	    !	perform nit relaxation iterations
  DO i=1,nc	    !	do rows:
    p=1.
    e=0.
    DO j=1,nc
      dw=ABS(w(i,j))
      IF(dw > c)THEN
        e=e+1.
        p=p*dw
      ENDIF
    ENDDO
    p=1./(p**(1./e))
    w(i,:) = w(i,:) * p            ! rescale this row of w..
    d1(i)=d1(i)*p			     ! ..and update d1 consistently
  ENDDO
  DO j=1,nc	    !	do columns:
    p=1.
    e=0.
    d2j=d2(j)
    DO i=1,nc
      dw=ABS(w(i,j))
      IF(dw > c)THEN
        e=e+1.
        p=p*dw
      ENDIF
    ENDDO
    p=1./(p**(1./e))
    w(:,j) = w(:,j) * p        ! rescale this column of w..
    d2(j)=d2(j)*p		       ! ..and update d2 consistently
  ENDDO
ENDDO
c=1.e-3	    ! final setting for criterion for "negligible" elements
ENDDO
ab(1:nc) = d1(1:nc) * ab(1:nc) ! rescale r.h.s vector by d1
p=1.			     ! p becomes product of row-lengths:
DO i=1,nc
   p=p*SQRT(dot_PRODUCT(w(i,:),w(i,:)))
ENDDO
CALL ldum(w,ipiv,d)
DO i=1,nc
  d=d*w(i,i)		      ! d becomes the determinant of w
ENDDO
wv(:,1) = ab ! convert shape of array
CALL udlmm(w,wv(:,1:1),ipiv)
ab = d2 * wv(:,1) ! rescale solution vector by d2
!     ab(1:nc) = d2(1:nc) * ab(1:nc) ! rescale solution vector by d2
!  note: it is very likely that round-off errors have accumulated during
!  the iterative rescaling of w. we invoke original matrix elements w2 and
!  substitute the tentative solution vector into the original (unscaled)
!  equation in order to estimate the residual components of roundoff error.

!  begin "clean-up" process. substitute solution vector in original
!  equation and leave the residual difference in v
v=v-MATMUL(w2,ab)
v = d1 * v    ! rescale the residual vector by d1
wv(:,1) = v ! convert shape of array
CALL udlmm(w,wv(:,1:1),ipiv) ! solve linear system with this rhs.
ab=ab+wv(:,1)*d2 ! add residual solution vector, 
                                      ! scaled, to ab

	  DEALLOCATE(wv)
w(1,1)=d/p  ! this ratio is an indicator of the overall conditioning
            ! when d/p is very small, treat the results with suspicion!

END SUBROUTINE linvan

!------------------------------------------------------------------------------
!   R.J.Purser, National Centers for Environmental Prediction, Washington D.C.
!   wd23jp@sun1.wwb.noaa.gov					      1996
!		    SUBROUTINE DLINVAN
!
!   Take square matrix W and seek row and column scalings to produce non-
!   vanishing elements of rescaled W having magnitudes as close to unity
!   as possible. The approach is make the geometric mean of the nonvanishing
!   elements of each row and of each column +1 or -1. Having rescaled the
!   matrix and the r.h.s. vector AB, compute the product P of row-vector
!   norms, then compute the determinant D and solve the linear system.
!   Rescale the solution vector (now AB) and put the conditioning indicator
!   formed by the ratio D/P into the first element of W.
!
! <-> W:    Generalized Vandermonde matrix in, conditioning indicator out.
! <-> AB:   R.h.s. vector in, solution vector of numerical coefficients out.
!------------------------------------------------------------------------------
SUBROUTINE dlinvan(w,ab)
INTEGER, PARAMETER     :: nit=20
REAL(8), INTENT(INOUT) :: w(:,:), ab(:)
REAL(8)                :: d1(SIZE(w,1)), d2(SIZE(w,1)), &
                          w2(SIZE(w,1),SIZE(w,1)),v(SIZE(w,1))
INTEGER                :: i, j, it, jt, ipiv(SIZE(w,1)), nc
REAL(8)                :: p, e, dw, c, d, s, d2j
REAL(8),ALLOCATABLE    :: wv(:,:) ! work variable for ab(nc) and v(nn)

nc = SIZE(w,DIM=1)
ALLOCATE(wv(nc,1))

w2=w                ! Preserve original W and AB for use
v = ab(1:nc)	    ! in later "clean-up" operation.

d1 = 1.0 	    ! Row scaling factors set to default
d2 = 1.0 	    ! Column scaling factors set to default

C=1.E-16	    ! Set initial criterion for "negligible" elements of W

! In first attempt to estimate row and column scalings, use logarithms
! to avoid the risk of under- or over-flows of the line products of W:
DO i=1,nc
  p=0.
  e=0.
  DO j=1,nc
    dw=ABS(w(i,j))
    IF(dw > c)THEN
      e=e+1.
      p=p+LOG(dw)
    ENDIF
  ENDDO
  IF(e == 0.d0)STOP 'w effectively singular in linvan'
  d1(i)=EXP(-p/e)
ENDDO
CALL muldm_d(d1,w2,w)

DO j=1,nc
  p=0.
  e=0.
  DO i=1,nc
    dw=ABS(w(i,j))
    IF(dw > c)THEN
      e=e+1.
      p=p+LOG(dw)
    ENDIF
  ENDDO
  IF(e == 0.)STOP 'w effectively singular in linvan'
  d2(j)=EXP(-p/e)
ENDDO
CALL mulmd_d(w,d2,w)

c=1.e-8  ! reset the criterion for "negligible" elements

! revert to iterations of the more efficient method without logarithms:
DO jt=1,2
DO it=1,nit	    !	perform nit relaxation iterations
  DO i=1,nc	    !	do rows:
    p=1.
    e=0.
    DO j=1,nc
      dw=ABS(w(i,j))
      IF(dw > c)THEN
        e=e+1.
        p=p*dw
      ENDIF
    ENDDO
    p=1./(p**(1./e))
    w(i,:) = w(i,:) * p            ! rescale this row of w..
    d1(i)=d1(i)*p			     ! ..and update d1 consistently
  ENDDO
  DO j=1,nc	    !	do columns:
    p=1.
    e=0.
    d2j=d2(j)
    DO i=1,nc
      dw=ABS(w(i,j))
      IF(dw > c)THEN
        e=e+1.
        p=p*dw
      ENDIF
    ENDDO
    p=1./(p**(1./e))
    w(:,j) = w(:,j) * p        ! rescale this column of w..
    d2(j)=d2(j)*p		       ! ..and update d2 consistently
  ENDDO
ENDDO
c=1.e-3	    ! final setting for criterion for "negligible" elements
ENDDO
ab(1:nc) = d1(1:nc) * ab(1:nc) ! rescale r.h.s vector by d1
p=1.			     ! p becomes product of row-lengths:
DO i=1,nc
   p=p*SQRT(dot_PRODUCT(w(i,:),w(i,:)))
ENDDO
CALL ldum_d(w,ipiv,d)
DO i=1,nc
  d=d*w(i,i)		      ! d becomes the determinant of w
ENDDO
wv(:,1) = ab ! convert shape of array
CALL udlmm_d(w,wv(:,1:1),ipiv)
ab = d2 * wv(:,1) ! rescale solution vector by d2
!     ab(1:nc) = d2(1:nc) * ab(1:nc) ! Rescale solution vector by D2
!  Note: it is very likely that round-off errors have accumulated during
!  the iterative rescaling of W. We invoke original matrix elements W2 and
!  substitute the tentative solution vector into the original (unscaled)
!  equation in order to estimate the residual components of roundoff error.

!  Begin "clean-up" process. Substitute solution vector in original
!  equation and leave the residual difference in V
v=v-MATMUL(w2,ab)
v = d1 * v    ! Rescale the residual vector by D1
wv(:,1) = v ! Convert shape of array
CALL UDLMM_d(w,wv(:,1:1),ipiv) ! Solve linear system with THIS rhs.
ab=ab+wv(:,1)*d2 ! Add residual solution vector, 
                                      ! scaled, to AB

	  DEALLOCATE(wv)
w(1,1)=d/p  ! this ratio is an indicator of the overall conditioning
            ! When D/P is very small, treat the results with suspicion!

END SUBROUTINE dlinvan

SUBROUTINE copdm(d,a)
REAL(8),DIMENSION(:),INTENT(IN)::d; REAL(8),DIMENSION(:,:),INTENT(OUT)::a; INTEGER i
                  a=0.; DO i=1,SIZE(a,1); a(i,i)= d(i); ENDDO; RETURN
ENTRY condm(d,a); a=0.; DO i=1,SIZE(a,1); a(i,i)=-d(i); ENDDO
END SUBROUTINE copdm

SUBROUTINE dcopdm(d,a)
REAL(8),DIMENSION(:),INTENT(IN)::d; REAL(8),DIMENSION(:,:),INTENT(OUT)::a
INTEGER i
                   a=0.; DO i=1,SIZE(a,1); a(i,i)= d(i); ENDDO; RETURN
ENTRY dcondm(d,a); a=0.; DO i=1,SIZE(a,1); a(i,i)=-d(i); ENDDO
END SUBROUTINE dcopdm

SUBROUTINE copsm(s,a)
REAL(8),INTENT(IN) :: s; REAL(8),DIMENSION(:,:),INTENT(OUT):: a; INTEGER i
                  a=0.; DO i=1,SIZE(a,1); a(i,i)= s; ENDDO; RETURN
ENTRY consm(s,a); a=0.; DO i=1,SIZE(a,1); a(i,i)=-s; ENDDO
END SUBROUTINE copsm

SUBROUTINE dcopsm(s,a)
REAL(8),INTENT(IN) :: s; REAL(8),DIMENSION(:,:),INTENT(OUT):: a; INTEGER i
                   a=0.; DO i=1,SIZE(a,1); a(i,i)= s; ENDDO; RETURN
ENTRY dconsm(s,a); a=0.; DO i=1,SIZE(a,1); a(i,i)=-s; ENDDO
END SUBROUTINE dcopsm

SUBROUTINE addmd(a,b,d)
REAL(8),DIMENSION(:,:),INTENT(INOUT):: a,b; REAL(8),DIMENSION(:),INTENT(IN):: d
REAL(8) s;  INTEGER i
                   b=a; DO i=1,SIZE(a,1); b(i,i)=b(i,i)+d(i); ENDDO; RETURN
ENTRY submd(a,b,d);b=a; DO i=1,SIZE(a,1); b(i,i)=b(i,i)-d(i); ENDDO; RETURN
ENTRY addms(a,b,s);b=a; DO I=1,SIZE(a,1); b(i,i)=b(i,i)+s;    ENDDO; RETURN
ENTRY SUBMS(A,B,S);b=a; DO I=1,SIZE(a,1); B(I,I)=B(I,I)-S;    ENDDO;
END SUBROUTINE addmd

SUBROUTINE daddmd(a,b,d)
REAL(8),DIMENSION(:,:),INTENT(INOUT)::A,B;REAL(8),DIMENSION(:),INTENT(IN)::D
REAL(8) s; INTEGER i
                     b=a; DO i=1,SIZE(a,1); b(i,i)=b(i,i)+d(i); ENDDO; RETURN
ENTRY DSUBMD(A,B,D); b=a; DO i=1,SIZE(a,1); b(i,i)=b(i,i)-d(i); ENDDO; RETURN
ENTRY DADDMS(A,B,S); b=a; DO i=1,SIZE(a,1); b(i,i)=b(i,i)+s;    ENDDO; RETURN
ENTRY DSUBMS(A,B,S); b=a; DO i=1,SIZE(a,1); b(i,i)=b(i,i)-s;    ENDDO;
END SUBROUTINE daddmd

SUBROUTINE l1lm(a,b) ! Cholesky, M -> L*U, U(i,j)=L(j,i)
REAL(8), INTENT(IN)    :: a(:,:)
REAL(8), INTENT(INOUT) :: b(:,:)
INTEGER             :: m,j, jm, jp, i
REAL(8)                :: s, bjji
m=SIZE(a,1)
DO j=1,m
  jm=j-1
  jp=j+1
  s = a(j,j) - SUM(b(j,1:jm)*b(j,1:jm))
  IF(S <= 0.)THEN
    PRINT '(" L1LM detects non-positivity at diagonal index",i2)',J
    STOP
  ENDIF
  b(j,j)=SQRT(s)
  bjji=1./b(j,j)
  DO i=jp,m
    s = a(i,j) - SUM(b(i,1:jm)*b(j,1:jm))
    b(i,j)=s*bjji
  ENDDO
  b(1:jm,j) = 0.0
ENDDO
END SUBROUTINE l1lm

SUBROUTINE DL1LM(A,B) ! Cholesky, M -> L*U, U(i,j)=L(j,i)
REAL(8), INTENT(IN)    :: a(:,:) 
REAL(8), INTENT(INOUT) :: b(:,:) 
INTEGER :: m,j, jm, jp, i
REAL(8) :: s, bjji
m=SIZE(a,1)
DO j=1,m
  jm=j-1
  jp=j+1
  s = a(j,j) - SUM(b(j,1:jm)*b(j,1:jm))
  IF(s <= 0.d0)THEN
    PRINT '(" L1LM detects non-positivity at diagonal index",i2)',J
    STOP
  ENDIF
  b(j,j)=SQRT(s)
  bjji=1./b(j,j)
  DO i=jp,m
    s = a(i,j) - SUM(b(i,1:jm)*b(j,1:jm))
    b(i,j)=s*bjji
  ENDDO
  b(1:jm,j) = 0.0
ENDDO
RETURN
END SUBROUTINE dl1lm

SUBROUTINE ldlm(a,b,d) ! Modified Cholesky decompose Q --> L*D*U, U(i,j)=L(j,i)
REAL(8), INTENT(IN)    :: a(:,:)
REAL(8), INTENT(INOUT) :: b(:,:)
REAL(8), INTENT(OUT)   :: d(:)
INTEGER :: m,j, jm, jp, i
REAL(8) :: s, bjji
m=SIZE(a,1)
DO j=1,m
  jm=j-1
  jp=j+1
  d(j)=a(j,j) - SUM(b(1:jm,j)*b(j,1:jm))
  
  b(j,j) = 1.
  IF(d(j) == 0.)THEN
    PRINT '(" LDLM detects singularity at diagonal index",i2)',J
    STOP
  ENDIF
  bjji=1./d(j)
  DO i=jp,m
     b(j,i)= a(i,j) - dot_PRODUCT(b(1:jm,j),b(i,1:jm))
     b(i,j)=b(j,i)*bjji
  ENDDO
ENDDO
CALL zeru(b)
RETURN
END SUBROUTINE ldlm

SUBROUTINE dldlm(a,b,d) ! Modified Cholesky  Q --> L*D*U, U(i,j)=L(j,i)
REAL(8), INTENT(IN)    :: a(:,:)
REAL(8), INTENT(INOUT) :: b(:,:)
REAL(8), INTENT(OUT)   :: d(:)
INTEGER                :: m,j, jm, jp, i
REAL(8)                :: s, bjji
m=SIZE(a,1)
DO j=1,m; jm=j-1; jp=j+1
  d(j)=a(j,j) - SUM(b(1:jm,j)*b(j,1:jm))
  b(j,j) = 1.
  IF(d(j) == 0.d0)THEN
    PRINT '(" DLDLM detects singularity at diagonal index",i2)',J
    STOP
  ENDIF
  bjji=1./d(j)
  DO i=jp,m
     b(j,i)= a(i,j) - dot_PRODUCT(b(1:jm,j),b(i,1:jm))
     b(i,j)=b(j,i)*bjji
  ENDDO
ENDDO
CALL zeru_d(b)
RETURN
END SUBROUTINE dldlm

!------------------------------------------------------------------------------
!   R.J.Purser, National Meteorological Center, Washington D.C.  1993
!		    SUBROUTINE INVH
!  Inver,t in place, a symmetric matrix
!
!  <-> A    symmetric square matrix, output as inverse of input
!
!  LIMITATION
!     This routine incorporates no pivoting - it is intended for matrices
!     that are already diagonally dominant
!------------------------------------------------------------------------------
SUBROUTINE invh(a)
REAL(8), INTENT(INOUT)      :: a(:,:) 
INTEGER                  :: m,k, kp, i, ip, j
REAL(8),DIMENSION(SIZE(a,1)):: d
m=SIZE(a,1)
!  PERFORM L.D.U DECOMPOSITION OF THE SYMMETRIC MATRIX:
CALL ldlm(a,a,d)

!  INVERT (IN PLACE) THE LOWER TRIANGULAR PART OF A, (ASSUMING UNIT
!  DIAGONAL ELEMENTS), AND INVERT THE DIAGONAL PART OF A (ASSUMING
!  ZERO OFF-DIAGONAL ELEMENTS). PUT TRANSPOSE OF LOWER, TIMES DIAGONAL,
!  INTO UPPER PART OF A.
DO k=1,m; kp=k+1
  a(k,k)=1./d(k)
  DO i=kp,m
    a(i,k) = a(i,k) + SUM(a(kp:i-1,k)*a(i,kp:i-1)) ! really??
    a(i,k)=-a(i,k)
  ENDDO
ENDDO

!  MULTIPLY: THE TRANSPOSE OF THE LOWER PART OF A (ASSUMING UNIT DIAGS),
!  TIMES THE DIAGONAL PART (ASSUMING ZERO OFF-DIAGS), TIMES THE LOWER
!  PART. THIS PRODUCT IS THE SYMMETRIC INVERSE OF THE ORIGINAL B.
DO i=2,m
  a(1:i-1,i) = a(i,1:i-1) * a(i,i) ! Really?
ENDDO
DO i=1,m
  ip=i+1
  DO j=1,i-1
    a(j,i) = a(j,i) + SUM(a(ip:ip+m-i-1,i)*a(j,ip:ip+m-i-1))
    a(i,j)=a(j,i)
  ENDDO
  a(i,i) = a(i,i) + SUM(a(ip:ip+m-i-1,i)*a(i,ip:ip+m-i-1))
ENDDO
END SUBROUTINE invh

SUBROUTINE dinvh(a)
REAL(8), INTENT(INOUT)      :: a(:,:) 
INTEGER                     :: m,k, kp, i, ip, j
REAL(8),DIMENSION(SIZE(a,1)):: d
m=SIZE(a,1)
!  PERFORM L.D.U DECOMPOSITION OF THE SYMMETRIC MATRIX:
CALL ldlm_d(a,a,d)

!  INVERT (IN PLACE) THE LOWER TRIANGULAR PART OF A, (ASSUMING UNIT
!  DIAGONAL ELEMENTS), AND INVERT THE DIAGONAL PART OF A (ASSUMING
!  ZERO OFF-DIAGONAL ELEMENTS). PUT TRANSPOSE OF LOWER, TIMES DIAGONAL,
!  INTO UPPER PART OF A.
DO k=1,m
  kp=k+1
  a(k,k)=1./d(k)
  DO i=kp,m
    a(i,k) = a(i,k) + SUM(a(kp:i-1,k)*a(i,kp:i-1)) ! really??
    a(i,k)=-a(i,k)
  ENDDO
ENDDO

!  MULTIPLY: THE TRANSPOSE OF THE LOWER PART OF A (ASSUMING UNIT DIAGS),
!  TIMES THE DIAGONAL PART (ASSUMING ZERO OFF-DIAGS), TIMES THE LOWER
!  PART. THIS PRODUCT IS THE SYMMETRIC INVERSE OF THE ORIGINAL B.
DO i=2,m
  a(1:i-1,i) = a(i,1:i-1) * a(i,i) ! really?
ENDDO
DO i=1,m
  ip=i+1
  DO j=1,i-1
    a(j,i) = a(j,i) + SUM(a(ip:ip+m-i-1,i)*a(j,ip:ip+m-i-1))
    a(i,j)=a(j,i)
  ENDDO
  a(i,i) = a(i,i) + SUM(a(ip:ip+m-i-1,i)*a(i,ip:ip+m-i-1))
ENDDO
END SUBROUTINE dinvh

!------------------------------------------------------------------------------
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE INVL
!     Invert lower triangular matrix in place if A are same
!------------------------------------------------------------------------------
SUBROUTINE invl(a)
REAL(8), INTENT(INOUT) :: a(:,:) 
INTEGER             :: m,j, i
REAL(8)                :: s
m=SIZE(a,1)
DO j=m,1,-1
  a(1:j-1,j) = 0.0
  a(j,j)=1./a(j,j)
  DO i=j+1,m
    s = SUM(a(j:i-1,j)*a(i,j:i-1))
    a(i,j)=-a(i,i)*s
  ENDDO
ENDDO
END SUBROUTINE invl

SUBROUTINE dinvl(a)
REAL(8), INTENT(INOUT) :: a(:,:) 
INTEGER                :: m,j, i
REAL(8)                :: s
m=SIZE(a,1)
DO j=m,1,-1
  a(1:j-1,j) = 0.0
  a(j,j)=1./a(j,j)
  DO i=j+1,m
    s = SUM(a(j:i-1,j)*a(i,j:i-1))
    a(i,j)=-a(i,i)*s
  ENDDO
ENDDO
END SUBROUTINE dinvl

!------------------------------------------------------------------------------
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE LINLV
!     Solve linear system involving lower triangular (LINLV) or upper
!     triangular (LINUV) matrix. u is input as right-hand-side, output
!     as the solution vector.
!------------------------------------------------------------------------------
SUBROUTINE linlv(a,u)
REAL(8), INTENT(IN)   :: a(:,:)
REAL(8), INTENT(INOUT):: u(:)
INTEGER            :: m,i, j, jp
DO i=1,SIZE(a,1);    u(i)=(u(i) - SUM(u(1:i-1)*a(i,1:i-1)))/a(i,i); ENDDO
RETURN
ENTRY linuv(a,u); m=SIZE(a,1)
DO j=m,1,-1; jp=j+1; u(j)=(u(j) - SUM(a(jp:m,j)*u(jp:m)))  /a(j,j); ENDDO
END SUBROUTINE linlv
SUBROUTINE dlinlv(a,u)
REAL(8), INTENT(IN)   :: a(:,:)
REAL(8), INTENT(INOUT):: u(:)
INTEGER :: m,i, j, jp
DO i=1,SIZE(a,1); u(i)= (u(i) - SUM(u(1:i-1)*a(i,1:i-1)))/a(i,i); ENDDO
RETURN
ENTRY dlinuv(a,u); m=SIZE(a,1)
DO j=m,1,-1; jp=j+1; u(j) = (u(j) - SUM(a(jp:m,j)*u(jp:m)))/a(j,j); ENDDO
END SUBROUTINE dlinlv


SUBROUTINE powp(a,b,n)	       ! Raise power series A to the power
INTEGER, INTENT(IN) :: n       ! of N and output as B
REAL(8),    INTENT(IN) :: a(0:)
REAL(8),    INTENT(OUT):: b(0:)
REAL(8),DIMENSION(0:SIZE(a)-1):: t; INTEGER :: k
b(0)=1.; b(1:) = 0.0; DO k=1,n; CALL mulpp(a,b,t); b=t; ENDDO
END SUBROUTINE powp
SUBROUTINE DPOWP(A,B,N)	       ! Raise power series A to the power
INTEGER,  INTENT(IN) :: n      ! of N and output as B
REAL(8), INTENT(IN) :: a(0:)
REAL(8), INTENT(OUT):: b(0:)
REAL(8),DIMENSION(0:SIZE(a)-1):: t; INTEGER :: k
B(0)=1.; b(1:) = 0.0; DO k=1,n; CALL mulpp_d(a,b,t); b=t; ENDDO
END SUBROUTINE dpowp


SUBROUTINE polps(a,s1,s2) ! Apply series A to scalar S1 to obtain S2
REAL(8),INTENT(IN) :: a(0:)
REAL(8),INTENT(IN) :: s1
REAL(8),INTENT(OUT):: s2
INTEGER m,k
m=SIZE(a)-1; s2=a(m); DO k=m-1,0,-1; s2=s2*s1+a(k); ENDDO
END SUBROUTINE polps
SUBROUTINE dpolps(a,s1,s2) ! Apply series A to scalar S1 to obtain S2
REAL(8),INTENT(IN) :: a(0:)
REAL(8),INTENT(IN) :: s1
REAL(8),INTENT(OUT):: s2
INTEGER m,k
m=SIZE(a)-1; s2=a(m); DO k=m-1,0,-1; s2=s2*s1+a(k); ENDDO
END SUBROUTINE dpolps

SUBROUTINE polpp(a,b,c) ! Apply power series A to power series B and put
                        ! the result out as power-series C.
REAL(8),INTENT(INOUT)         :: a(0:),b(0:),c(0:)
REAL(8),DIMENSION(0:SIZE(a)-1):: t
INTEGER m,k
m=SIZE(a)-1; c(0)=a(m); c(1:m) = 0.0
DO k=m-1,0,-1; CALL mulpp(b,c,t); c=t; c(0)=c(0)+a(k); ENDDO
END SUBROUTINE polpp
SUBROUTINE dpolpp(a,b,c) ! Apply power series A to power series B and put
                         ! the result out as power-series C.
REAL(8),INTENT(INOUT)         :: a(0:),b(0:),c(0:)
REAL(8),DIMENSION(0:SIZE(a)-1):: t
INTEGER m,k
m=SIZE(a)-1
c(0)=a(m); c(1:m) = 0.0
DO k=m-1,0,-1; CALL mulpp_d(b,c,t); c=t; c(0)=c(0)+a(k); ENDDO
END SUBROUTINE dpolpp

FUNCTION trcm(a) RESULT(trc_res)	    ! Trace of square matrix A
REAL(8)             :: trc_res
REAL(8), INTENT(IN) :: a(:,:)
INTEGER          :: i
trc_res=0.; DO i=1,SIZE(a,1); trc_res=trc_res+a(i,i); ENDDO
END FUNCTION trcm
FUNCTION dtrcm(a) RESULT(trc_res)	    ! Trace of square matrix A
REAL(8)             :: trc_res
REAL(8), INTENT(IN) :: a(:,:)
INTEGER              :: m,i
trc_res=0.; DO i=1,SIZE(a,1); trc_res=trc_res+a(i,i); ENDDO
END FUNCTION dtrcm

SUBROUTINE invmt(a)
REAL(8),DIMENSION(:,:),INTENT(INOUT):: a
INTEGER m,i,j,jp,l
REAL(8) d
INTEGER,DIMENSION(SIZE(a,1)):: ipiv
m=SIZE(a,1)
IF(m /= SIZE(a,2))STOP 'matrix passed to invmt is not square'
! Perform a pivoted L-D-U decomposition on matrix a:
CALL ldum(a,ipiv,d)

! Invert upper triangular portion U in place:
DO i=1,m; a(i,i)=1./a(i,i); ENDDO
DO i=1,m-1
   DO j=i+1,m; a(i,j)=-a(j,j)*DOT_PRODUCT(a(i:j-1,j),a(i,i:j-1)); ENDDO
ENDDO

! Invert lower triangular portion L in place:
DO j=1,m-1; jp=j+1
   DO i=jp,m; a(i,j)=-a(i,j)-DOT_PRODUCT(a(jp:i-1,j),a(i,jp:i-1)); ENDDO
ENDDO

!  Form the product of U**-1 and L**-1 in place
DO j=1,m-1; jp=j+1
   DO i=1,j; a(i,j)=a(i,j)+DOT_PRODUCT(a(jp:m,j),a(i,jp:m)); ENDDO
   DO i=jp,m; a(i,j)=DOT_PRODUCT(a(i:m,j),a(i,i:m));         ENDDO
ENDDO

!  Permute columns according to ipiv
DO j=m-1,1,-1; l=ipiv(j); CALL swpvv(a(:,j),a(:,l)); ENDDO
END SUBROUTINE invmt

SUBROUTINE dinvmt(a)
REAL(8),DIMENSION(:,:),INTENT(INOUT):: a
INTEGER                             :: m,i,j,jp,l
REAL(8)                             :: d
INTEGER,DIMENSION(SIZE(a,1))        :: ipiv
m=SIZE(a,1)
IF(m /= SIZE(a,2))STOP 'matrix passed to dinvmt is not square'
! Perform a pivoted L-D-U decomposition on matrix a:
CALL ldum_d(a,ipiv,d)

! Invert upper triangular portion U in place:
DO i=1,m; a(i,i)=1./a(i,i); ENDDO
DO i=1,m-1
   DO j=i+1,m; a(i,j)=-a(j,j)*DOT_PRODUCT(a(i:j-1,j),a(i,i:j-1)); ENDDO
ENDDO

! Invert lower triangular portion L in place:
DO j=1,m-1; jp=j+1
   DO i=jp,m; a(i,j)=-a(i,j)-DOT_PRODUCT(a(jp:i-1,j),a(i,jp:i-1)); ENDDO
ENDDO

!  Form the product of U**-1 and L**-1 in place
DO j=1,m-1; jp=j+1
   DO i=1,j; a(i,j)=a(i,j)+DOT_PRODUCT(a(jp:m,j),a(i,jp:m)); ENDDO
   DO i=jp,m; a(i,j)=DOT_PRODUCT(a(i:m,j),a(i,i:m));         ENDDO
ENDDO

!  Permute columns according to ipiv
DO j=m-1,1,-1; l=ipiv(j); CALL swpvv_d(a(:,j),a(:,l)); ENDDO
END SUBROUTINE dinvmt

SUBROUTINE linmmt(a,b)
REAL(8),DIMENSION(:,:),INTENT(INOUT):: a,b
INTEGER,DIMENSION(SIZE(a,1))     :: ipiv
INTEGER                          :: m
REAL(8)                             :: d
m=SIZE(a,1)
IF(m /= SIZE(a,2))STOP 'matrix passed to linmmt is not square'
IF(m /= SIZE(b,1))STOP 'matrix and vectors in linmmt have unmatched sizes'
CALL ldum(a,ipiv,d); CALL udlmm(a,b,ipiv)
END SUBROUTINE linmmt

SUBROUTINE dlinmmt(a,b)
REAL(8),DIMENSION(:,:),INTENT(INOUT):: a,b
INTEGER,DIMENSION(SIZE(a,1))        :: ipiv
INTEGER                             :: m 
REAL(8)                             :: d
m=SIZE(a,1)
IF(m /= SIZE(a,2))STOP 'matrix passed to linmmt_d is not square'
IF(m /= SIZE(b,1))STOP 'matrix and vectors in linmmt_d have unmatched sizes'
CALL ldum_d(a,ipiv,d); CALL udlmm_d(a,b,ipiv)
END SUBROUTINE dlinmmt

SUBROUTINE linmvt(a,b)
REAL(8),DIMENSION(:,:),INTENT(INOUT):: a
REAL(8),DIMENSION(:),  INTENT(INOUT):: b
INTEGER,DIMENSION(SIZE(a,1))     :: ipiv
INTEGER                          :: m
REAL(8)                             :: d
m=SIZE(a,1)
IF(m /= SIZE(a,2))STOP 'matrix passed to linmvt is not square'
IF(m /= SIZE(b))STOP 'matrix and vectors in linmvt have unmatched sizes'
CALL ldum(a,ipiv,d); CALL udlmm(a,b,ipiv)
END SUBROUTINE linmvt

SUBROUTINE dlinmvt(a,b)
REAL(8),DIMENSION(:,:),INTENT(INOUT):: a
REAL(8),DIMENSION(:),  INTENT(INOUT):: b
INTEGER,DIMENSION(SIZE(a,1))        :: ipiv
INTEGER m; REAL(8) d
m=SIZE(a,1)
IF(m /= SIZE(a,2))STOP 'matrix passed to linmvt_d is not square'
IF(m /= SIZE(b))STOP 'matrix and vectors in linmvt_d have unmatched sizes'
CALL ldum_d(a,ipiv,d); CALL udlmm_d(a,b,ipiv)
END SUBROUTINE dlinmvt

end module module_pmat1

!ENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDEND
!						********************
!						* module_pmat1.f90 *
!						* PURSER 1994/1999 *
!                                               *   FUJITA 1999    *
!						********************
!ENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDEND

!						********************
!						* module_pmat2.f90 *
!						* PURSER 1994/1999 *
!                                               *   FUJITA 1999    *
!						********************
MODULE MODULE_pmat2
USE MODULE_pmat1
IMPLICIT NONE

INTERFACE avco;   MODULE PROCEDURE avco;           END INTERFACE
INTERFACE avco_d; MODULE PROCEDURE davco;          END INTERFACE
INTERFACE dfco;   MODULE PROCEDURE dfco;           END INTERFACE
INTERFACE dfco_d; MODULE PROCEDURE ddfco;          END INTERFACE
INTERFACE dfco2;  MODULE PROCEDURE dfco2;          END INTERFACE
INTERFACE dfco2_d;MODULE PROCEDURE ddfco2;         END INTERFACE
INTERFACE clib;   MODULE PROCEDURE clib;           END INTERFACE
INTERFACE clib_d; MODULE PROCEDURE dclib;          END INTERFACE
INTERFACE cad1b;  MODULE PROCEDURE cad1b;          END INTERFACE
INTERFACE csb1b;  MODULE PROCEDURE csb1b;          END INTERFACE
INTERFACE cad2b;  MODULE PROCEDURE cad2b;          END INTERFACE
INTERFACE csb2b;  MODULE PROCEDURE csb2b;          END INTERFACE
INTERFACE copbt;  MODULE PROCEDURE copbt;          END INTERFACE
INTERFACE conbt;  MODULE PROCEDURE conbt;          END INTERFACE
INTERFACE copmb;  MODULE PROCEDURE copmb;          END INTERFACE
INTERFACE conmb;  MODULE PROCEDURE conmb;          END INTERFACE
INTERFACE copbm;  MODULE PROCEDURE copbm;          END INTERFACE
INTERFACE conbm;  MODULE PROCEDURE conbm;          END INTERFACE
INTERFACE mulbb;  MODULE PROCEDURE mulbb;          END INTERFACE
INTERFACE madbb;  MODULE PROCEDURE madbb;          END INTERFACE
INTERFACE msbbb;  MODULE PROCEDURE msbbb;          END INTERFACE
INTERFACE ldub;   MODULE PROCEDURE ldub;           END INTERFACE
INTERFACE ldub_d; MODULE PROCEDURE dldub;          END INTERFACE
INTERFACE l1ubb;  MODULE PROCEDURE l1ubb;          END INTERFACE
INTERFACE l1ubb_d;MODULE PROCEDURE dl1ubb;         END INTERFACE
INTERFACE l1ueb;  MODULE PROCEDURE l1ueb;          END INTERFACE
INTERFACE l1ueb_d;MODULE PROCEDURE dl1ueb;         END INTERFACE
INTERFACE l1lb;   MODULE PROCEDURE l1lb;           END INTERFACE
INTERFACE ldlb;   MODULE PROCEDURE ldlb;           END INTERFACE
INTERFACE ldlb_d; MODULE PROCEDURE dldlb;          END INTERFACE
INTERFACE udub;   MODULE PROCEDURE udub;           END INTERFACE
INTERFACE udub_d; MODULE PROCEDURE dudub;          END INTERFACE
INTERFACE mulbv;  MODULE PROCEDURE mulbv;          END INTERFACE
INTERFACE madbv;  MODULE PROCEDURE madbv;          END INTERFACE
INTERFACE msbbv;  MODULE PROCEDURE msbbv;          END INTERFACE
INTERFACE mulbx;  MODULE PROCEDURE mulbx;          END INTERFACE
INTERFACE madbx;  MODULE PROCEDURE madbx;          END INTERFACE
INTERFACE msbbx;  MODULE PROCEDURE msbbx;          END INTERFACE
INTERFACE mulby;  MODULE PROCEDURE mulby;          END INTERFACE
INTERFACE madby;  MODULE PROCEDURE madby;          END INTERFACE
INTERFACE msbby;  MODULE PROCEDURE msbby;          END INTERFACE
INTERFACE mulvb;  MODULE PROCEDURE mulvb;          END INTERFACE
INTERFACE madvb;  MODULE PROCEDURE madvb;          END INTERFACE
INTERFACE msbvb;  MODULE PROCEDURE msbvb;          END INTERFACE
INTERFACE mulxb;  MODULE PROCEDURE mulxb;          END INTERFACE
INTERFACE madxb;  MODULE PROCEDURE madxb;          END INTERFACE
INTERFACE msbxb;  MODULE PROCEDURE msbxb;          END INTERFACE
INTERFACE mulyb;  MODULE PROCEDURE mulyb;          END INTERFACE
INTERFACE madyb;  MODULE PROCEDURE madyb;          END INTERFACE
INTERFACE msbyb;  MODULE PROCEDURE msbyb;          END INTERFACE
INTERFACE mulbd;  MODULE PROCEDURE mulbd;          END INTERFACE
INTERFACE madbd;  MODULE PROCEDURE madbd;          END INTERFACE
INTERFACE msbbd;  MODULE PROCEDURE msbbd;          END INTERFACE
INTERFACE muldb;  MODULE PROCEDURE muldb;          END INTERFACE
INTERFACE maddb;  MODULE PROCEDURE maddb;          END INTERFACE
INTERFACE msbdb;  MODULE PROCEDURE msbdb;          END INTERFACE
INTERFACE udlbv;  MODULE PROCEDURE udlbv;          END INTERFACE
INTERFACE udlbx;  MODULE PROCEDURE udlbx;          END INTERFACE
INTERFACE udlby;  MODULE PROCEDURE udlby;          END INTERFACE
INTERFACE udlvb;  MODULE PROCEDURE udlvb;          END INTERFACE
INTERFACE udlxb;  MODULE PROCEDURE udlxb;          END INTERFACE
INTERFACE udlyb;  MODULE PROCEDURE udlyb;          END INTERFACE
INTERFACE u1lbv;  MODULE PROCEDURE u1lbv;          END INTERFACE
INTERFACE u1lbx;  MODULE PROCEDURE u1lbx;          END INTERFACE
INTERFACE u1lby;  MODULE PROCEDURE u1lby;          END INTERFACE
INTERFACE u1lvb;  MODULE PROCEDURE u1lvb;          END INTERFACE
INTERFACE u1lxb;  MODULE PROCEDURE u1lxb;          END INTERFACE
INTERFACE u1lyb;  MODULE PROCEDURE u1lyb;          END INTERFACE
INTERFACE linbv;  MODULE PROCEDURE linbv;          END INTERFACE
INTERFACE wrtb;   MODULE PROCEDURE wrtb;           END INTERFACE


CONTAINS

!=============================================================================
SUBROUTINE davco(na,nb,za,zb,z0,a,b) 
!=============================================================================
!		    SUBROUTINE DAVCO
!   R.J.Purser, National Centers for Environmental Prediction, Washington D.C.
!   jpurser@ncep.noaa.gov					      1999
!
!  Compute one row of the coefficients for the compact mid-interval
!  interpolation scheme characterized by matrix equation of the form,
!			 A.t = B.s			       (*)
!  Where s is the vector of "source" values, t the staggered "target" values.
!
! --> NA:   number of t-points operated on by this row of the A of (*)
! --> NB:   number of s-points operated on by this row of the B of (*)
! --> ZA:   coordinates of t-points used in this row of (*)
! --> ZB:   coordinates of s-points used in this row of (*)
! --> Z0:   nominal point of application of this row of (*)
! <-- A:    the NA coefficients A for this scheme
! <-- B:    the NB coefficients B for this scheme
!=============================================================================
INTEGER, INTENT(IN )          :: na,nb
REAL(8), INTENT(IN )          :: za(na),zb(nb),z0
REAL(8), INTENT(OUT)          :: a(na),b(nb)
!-----------------------------------------------------------------------------
INTEGER                       :: na1,nab,i
REAL(8),DIMENSION(na+nb,na+nb):: w
REAL(8),DIMENSION(na)         :: za0,pa
REAL(8),DIMENSION(nb)         :: zb0,pb
REAL(8),DIMENSION(na+nb)      :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0; zb0=zb-z0
pa=1.;     pb=-1.
w=0.;         ab=0.
w(1,1:na)=1.; ab(1)=1.
DO i=2,nab; w(i,1:na)=pa;    pa=pa*za0; w(i,na1:nab)=pb; pb=pb*zb0; ENDDO
CALL inv_d(w,ab)
a=ab(1:na); b=ab(na1:nab)
END SUBROUTINE davco

!=============================================================================
SUBROUTINE avco(na,nb,za,zb,z0,a,b) 
!=============================================================================
INTEGER, INTENT(IN )        :: na,nb
REAL(8),    INTENT(IN )        :: za(na),zb(nb),z0
REAL(8),    INTENT(OUT)        :: a(na),b(nb)
!-----------------------------------------------------------------------------
INTEGER                     :: na1,nab,i
REAL(8), DIMENSION(na+nb,na+nb):: w
REAL(8), DIMENSION(na)         :: za0,pa
REAL(8), DIMENSION(nb)         :: zb0,pb
REAL(8), DIMENSION(na+nb)      :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0; zb0=zb-z0
pa=1.;     pb=-1.
w=0.;         ab=0.
w(1,1:na)=1.; ab(1)=1.
DO i=2,nab; w(i,1:na)=pa;    pa=pa*za0; w(i,na1:nab)=pb; pb=pb*zb0; ENDDO
CALL inv(w,ab)
a=ab(1:na); b=ab(na1:nab)
END SUBROUTINE avco 


SUBROUTINE ddfco(na,nb,za,zb,z0,a,b) 
!=============================================================================
!   R.J.Purser, National Centers for Environmental Prediction, Washington D.C.
!   jpurser@ncep.noaa.gov					      1999
!		    SUBROUTINE DDFCO
!
!  Compute one row of the coefficients for either the compact differencing or
!  quadrature scheme characterized by matrix equation of the form,
!			 A.d = B.c			       (*)
!  In either case, d is the derivative of c.
!
! --> NA:   number of d-points operated on by this row of the A of (*)
! --> NB:   number of c-points operated on by this row of the B of (*)
! --> ZA:   coordinates of d-points used in this row of (*)
! --> ZB:   coordinates of c-points used in this row of (*)
! --> Z0:   nominal point of application of this row of (*)
! <-- A:    the A-coefficients for this scheme
! <-- B:    the B-coefficients for this scheme
!=============================================================================
INTEGER, INTENT(IN)            :: na,nb
REAL(8), INTENT(IN)            :: za(na),zb(nb),z0
REAL(8), INTENT(OUT)           :: a(na),b(nb)
!-----------------------------------------------------------------------------
INTEGER                        :: na1,nab,i
REAL(8), DIMENSION(na+nb,na+nb):: w
REAL(8), DIMENSION(na)         :: za0,pa
REAL(8), DIMENSION(nb)         :: zb0,pb
REAL(8), DIMENSION(na+nb)      :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0; zb0=zb-z0
pa=1.;     pb=-1.
w=0.;         ab=0.
w(1,1:na)=1.; ab(1)=1.
DO i=3,nab; w(i,1:na)   =pa*(i-2); pa=pa*za0; ENDDO
DO i=2,nab; w(i,na1:nab)=pb;       pb=pb*zb0; ENDDO
CALL inv_d(w,ab)
a=ab(1:na); b=ab(na1:nab)
END SUBROUTINE ddfco 

!=============================================================================
SUBROUTINE dfco(na,nb,za,zb,z0,a,b)
!=============================================================================
INTEGER, INTENT(IN )        :: na,nb
REAL(8),    INTENT(IN )        :: za(na),zb(nb),z0
REAL(8),    INTENT(OUT)        :: a(na),b(nb)
!-----------------------------------------------------------------------------
INTEGER:: na1,nab,i
REAL(8), DIMENSION(na+nb,na+nb):: w
REAL(8), DIMENSION(na)         :: za0,pa
REAL(8), DIMENSION(nb)         :: zb0,pb
REAL(8), DIMENSION(na+nb)      :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0; zb0=zb-z0
pa=1.;     pb=-1.
w=0.;         ab=0.
w(1,1:na)=1.; ab(1)=1.
DO i=3,nab; w(i,1:na)   =pa*(i-2); pa=pa*za0; ENDDO
DO i=2,nab; w(i,na1:nab)=pb;       pb=pb*zb0; ENDDO
CALL inv(w,ab)
a=ab(1:na); b=ab(na1:nab)
END SUBROUTINE dfco 

!=============================================================================
SUBROUTINE ddfco2(na,nb,za,zb,z0,a,b) 
!=============================================================================
!		    SUBROUTINE DDFCO2
!   R.J.Purser, National Centers for Environmental Prediction, Washington D.C.
!   jpurser@ncep.noaa.gov					      1999
!
!  Compute one row of the coefficients for either the compact second-
!  differencing scheme characterized by matrix equation of the form,
!			 A.d = B.c			       (*)
!  Where d is the second-derivative of c.
!
! --> NA:   number of d-points operated on by this row of the A of (*)
! --> NB:   number of c-points operated on by this row of the B of (*)
! --> ZA:   coordinates of d-points used in this row of (*)
! --> ZB:   coordinates of c-points used in this row of (*)
! --> Z0:   nominal point of application of this row of (*)
! <-- A:    the NA coefficients A for this scheme
! <-- B:    the NB coefficients B for this scheme
!=============================================================================
INTEGER, INTENT(IN )           :: na,nb
REAL(8), INTENT(IN )           :: za(na),zb(nb),z0
REAL(8), INTENT(OUT)           :: a(na),b(nb)
!-----------------------------------------------------------------------------
INTEGER                        :: na1,nab,i
REAL(8), DIMENSION(na+nb,na+nb):: w
REAL(8), DIMENSION(na)         :: za0,pa
REAL(8), DIMENSION(nb)         :: zb0,pb
REAL(8), DIMENSION(na+nb)      :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0; zb0=zb-z0
pa=1.;     pb=-1.
w=0.;         ab=0.
w(1,1:na)=1.; ab(1)=1.
DO i=4,nab; w(i,1:na)   =pa*(i-2)*(i-3); pa=pa*za0; ENDDO
DO i=2,nab; w(i,na1:nab)=pb;             pb=pb*zb0; ENDDO
CALL inv_d(w,ab)
a=ab(1:na); b=ab(na1:nab)
END SUBROUTINE ddfco2 

!=============================================================================
SUBROUTINE dfco2(na,nb,za,zb,z0,a,b) 
!=============================================================================
INTEGER, INTENT(IN )        :: na,nb
REAL(8),    INTENT(IN )        :: za(na),zb(nb),z0
REAL(8),    INTENT(OUT)        :: a(na),b(nb)
!-----------------------------------------------------------------------------
INTEGER:: na1,nab,i
REAL(8), DIMENSION(na+nb,na+nb):: w
REAL(8), DIMENSION(na)         :: za0,pa
REAL(8), DIMENSION(nb)         :: zb0,pb
REAL(8), DIMENSION(na+nb)      :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0; zb0=zb-z0
pa=1.;     pb=-1.
w=0.;         ab=0.
w(1,1:na)=1.; ab(1)=1.
DO i=4,nab; w(i,1:na)   =pa*(i-2)*(i-3); pa=pa*za0; ENDDO
DO i=2,nab; w(i,na1:nab)=pb;             pb=pb*zb0; ENDDO
CALL inv(w,ab)
a=ab(1:na); b=ab(na1:nab)
END SUBROUTINE dfco2 

!=============================================================================
SUBROUTINE clib(a,m1,m2,mah1,mah2) ! Clip the dead space of the band matrix, a
!=============================================================================
INTEGER, INTENT(IN)   :: m1, m2, mah1, mah2
REAL(8), INTENT(INOUT):: a(m1,-mah1:mah2)
INTEGER               :: j
IF(m2-m1+mah1 < 0)STOP 'In CLIB, form of band matrix implies redundant rows'
DO j=1,mah1; a(1:min(m1,j),-j)=0.; ENDDO; DO j=m2-m1+1,mah2; a(max(1,m2-j+1):m1,j)=0.; ENDDO
END SUBROUTINE clib


!=============================================================================
SUBROUTINE dclib(a,m1,m2,mah1,mah2) ! Clip dead space of the band matrix, a
!=============================================================================
INTEGER, INTENT(IN)   :: m1, m2, mah1, mah2
REAL(8), INTENT(INOUT):: a(m1,-mah1:mah2)
INTEGER               :: j
IF(m2-m1+mah1 < 0)STOP 'In CLIB_d, form of band matrix implies redundant rows'
DO j=1,mah1; a(1:min(m1,j),-j)=0.; ENDDO; DO j=m2-m1+1,mah2; a(max(1,m2-j+1):m1,j)=0.; ENDDO
END SUBROUTINE dclib

SUBROUTINE cad1b(a,m1,m2,mah1,mah2,mirror2)
!=============================================================================
! Incorporate operand symmetry near end-1 of a band matrix operator
!
! <-> A:      Input as unclipped operator, output as symmetrized and clipped.
! m1, m2:     Sizes of implied full matrix
! mah1, mah2: Left and right semi-bandwidths of A.
! mirror2:    2*location of symmetry axis relative to end-1 operand element.
!      Note: although m2 is not used here, it IS used in companion routines
!            cad2b and csb2b; it is retained in the interests of uniformity.
!=============================================================================
INTEGER,  INTENT(IN)   :: m1,m2,mah1,mah2,mirror2
REAL(8),     INTENT(INOUT):: a(0:m1-1,-mah1:mah2)
INTEGER                :: i,i2,jm,jp,jpmax
IF(mirror2+mah1 > mah2)STOP 'In cad1b, mah2 insufficient'
DO i=0,m1-1; i2=i*2; jpmax=mirror2+mah1-i2; IF(jpmax <= -mah1)EXIT
   DO jm=-mah1,mah2; jp=mirror2-jm-i2; IF(jp <= jm)EXIT
      a(i,jp)=a(i,jp)+a(i,jm) ! Reflect and add
      a(i,jm)=0.              ! zero the exterior part
   ENDDO
ENDDO
RETURN
!=============================================================================
ENTRY     csb1b(a,m1,m2,mah1,mah2,mirror2)
!=============================================================================
! Like cad1b, but for antisymmetric operand
IF(mirror2+mah1 > mah2)STOP 'In csb1b, mah2 insufficient'
DO i=0,m1-1; i2=i*2; jpmax=mirror2+mah1-i2; IF(jpmax < -mah1)EXIT
   DO jm=-mah1,mah2; jp=mirror2-jm-i2; IF(jp < jm)EXIT
      a(i,jp)=a(i,jp)-a(i,jm) ! Reflect and subtract
      a(i,jm)=0.              ! zero the exterior part
   ENDDO
ENDDO
END SUBROUTINE cad1b

!=============================================================================
SUBROUTINE cad2b(a,m1,m2,mah1,mah2,mirror2)
!=============================================================================
! Incorporate operand symmetry near end-2 of a band matrix operator
!
! <-> A:      Input as unclipped operator, output as symmetrized and clipped.
! m1, m2:     Sizes of implied full matrix
! mah1, mah2: Left and right semi-bandwidths of A.
! mirror2:    2*location of symmetry axis relative to end-2 operand element.
!=============================================================================
INTEGER,  INTENT(IN)   :: m1,m2,mah1,mah2,mirror2
REAL(8),     INTENT(INOUT):: a(1-m1:0,m1-m2-mah1:m1-m2+mah2)
INTEGER                :: i,i2,jm,jp,jmmin,nah1,nah2,mirror,j0
nah1=mah1+m2-m1; nah2=mah2+m1-m2 ! Effective 2nd-index bounds of A
IF(mirror2-nah1 > -nah2)STOP 'In cad2b, mah1 insufficient'
DO i=0,1-m1,-1; i2=i*2; jmmin=mirror2-nah2-i2; IF(jmmin >= nah2)EXIT
   DO jp=nah2,nah1,-1; jm=mirror2-jp-i2; IF(jm >= jp)EXIT
      a(i,jm)=a(i,jm)+a(i,jp) ! Reflect and add
      a(i,jp)=0.              ! zero the exterior part
   ENDDO
ENDDO
RETURN
!=============================================================================
ENTRY    csb2b(a,m1,m2,mah1,mah2,mirror2)
!=============================================================================
nah1=mah1+m2-m1; nah2=mah2+m1-m2 ! Effective 2nd-index bounds of A
IF(mirror2-nah1 > -nah2)STOP 'In csb2b, mah1 insufficient'
DO i=0,1-m1,-1; i2=i*2; jmmin=mirror2-nah2-i2; IF(jmmin > nah2)EXIT
   DO jp=nah2,nah1,-1; jm=mirror2-jp-i2; IF(jm > jp)EXIT
      a(i,jm)=a(i,jm)-a(i,jp) ! Reflect and subtract
      a(i,jp)=0.              ! zero the exterior part
   ENDDO
ENDDO
!=============================================================================
ENTRY    cex2b(a,m1,m2,mah1,mah2,mirror2)
!=============================================================================
nah1=mah1+m2-m1; nah2=mah2+m1-m2 ! Effective 2nd-index bounds of A
IF(mirror2-nah1 > -nah2)STOP 'In cex2b, mah1 insufficient'
mirror=mirror2/2
IF(mirror*2 /= mirror2)STOP 'In cex2b, mirror2 is not even'
DO i=0,1-m1,-1; i2=i*2; jmmin=mirror2-nah2-i2; IF(jmmin >= nah2)EXIT
   j0=mirror-i
   DO jp=nah2,nah1,-1; jm=mirror2-jp-i2; IF(jm >= jp)EXIT
      a(i,jm)=a(i,jm)-a(i,jp)    ! Reflect and subtract
      a(i,j0)=a(i,j0)+2.*a(i,jp) ! Apply double the coefficient to end
      a(i,jp)=0.                 ! zero the exterior part
   ENDDO
ENDDO
END SUBROUTINE cad2b

!=============================================================================
SUBROUTINE copbt(a,b,m1,m2,mah1,mah2)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE COPBT
!  Copy transpose of rectangular banded matrix A to B
!  Note: this routine expects A and B always to occupy separate storage.
!
! --> A  input matrix in banded format
! <-- B  output matrix in banded format
! --> M1 number of rows of A, columns of B
! --> M2 number of columns of A, rows of B
! --> MAH1 left-half-bandwidth of A, right-half-bandwidth of B
! --> MAH2 right-half-bandwidth of A, left-half-bandwidth of B
!=============================================================================
INTEGER,  INTENT(IN) :: m1, m2, mah1, mah2
REAL(8),     INTENT(IN) :: a(m1,-mah1:mah2)
REAL(8),     INTENT(OUT):: b(m2,-mah2:mah1)
INTEGER              :: j, i
CALL clib(b,mah2,mah1,m2,m1)
DO j=-mah1,mah2
   DO i=MAX(1,1-j),MIN(m1,m2-j); b(j+i,-j)=a(i,j); ENDDO
ENDDO
RETURN
ENTRY	 conbt(a,b,m1,m2,mah1,mah2)
CALL clib(b,mah2,mah1,m2,m1)
DO j=-mah1,mah2
   DO i=MAX(1,1-j),MIN(m1,m2-j); b(j+i,-j)=-a(i,j); ENDDO
ENDDO
END SUBROUTINE copbt

!=============================================================================
SUBROUTINE copmb(afull,aband,m1,m2,mah1,mah2)
!=============================================================================
INTEGER,                           INTENT(IN) :: m1, m2, mah1, mah2
REAL(8),     DIMENSION(m1,m2),        INTENT(IN) :: afull
REAL(8),     DIMENSION(m1,-mah1:mah2),INTENT(OUT):: aband
INTEGER                                       :: i1,i2, i, j
CALL clib(aband,m1,m2,mah1,mah2)
DO j=1,m1; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
   DO i=i1,i2; aband(i,j)= afull(i,j+i); ENDDO
ENDDO
RETURN
!=============================================================================
ENTRY      conmb(afull,aband,m1,m2,mah1,mah2)
!=============================================================================
CALL clib(aband,m1,m2,mah1,mah2)
DO j=1,m1; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
   DO i=i1,i2; aband(i,j)=-afull(i,j+i); ENDDO
ENDDO
END SUBROUTINE copmb

!=============================================================================
SUBROUTINE copbm(aband,afull,m1,m2,mah1,mah2)
!=============================================================================
INTEGER,                           INTENT(IN) :: m1, m2, mah1, mah2
REAL(8),     DIMENSION(m1,-mah1:mah2),INTENT(IN) :: aband
REAL(8),     DIMENSION(m1,m2),        INTENT(OUT):: afull
INTEGER                                       :: i1,i2, i, j
afull=0.
DO j=1,m1; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
   DO i=i1,i2; afull(i,j+i)= aband(i,j); ENDDO
ENDDO
RETURN
!=============================================================================
ENTRY      conbm(aband,afull,m1,m2,mah1,mah2)
!=============================================================================
afull=0.
DO j=1,m1; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
   DO i=i1,i2; afull(i,j+i)=-aband(i,j); ENDDO
ENDDO
END SUBROUTINE copbm
 
!=============================================================================
SUBROUTINE mulbb(a,b,c,m1,m2,mah1,mah2,mbh1,mbh2,mch1,mch2)
!=============================================================================
INTEGER,  INTENT(IN)   :: m1, m2, mah1, mah2, mbh1, mbh2, mch1, mch2
REAL(8),     INTENT(IN)   :: a(m1,-mah1:mah2), b(m2,-mbh1:mbh2)
REAL(8),     INTENT(INOUT):: c(m1,-mch1:mch2)
INTEGER                :: nch1, nch2, j, k, jpk, i1,i2
c=0.0
ENTRY      madbb(a,b,c,m1,m2,mah1,mah2,mbh1,mbh2,mch1,mch2)
nch1=mah1+mbh1; nch2=mah2+mbh2
IF(nch1 /= mch1 .OR. nch2 /= mch2)STOP 'In MULBB, dimensions inconsistent'
DO j=-mah1,mah2
   DO k=-mbh1,mbh2; jpk=j+k; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
      c(i1:i2,jpk)=c(i1:i2,jpk)+a(i1:i2,j)*b(j+i1:j+i2,k)
   ENDDO
ENDDO
END SUBROUTINE mulbb

!=============================================================================
SUBROUTINE msbbb(a,b,c,m1,m2,mah1,mah2,mbh1,mbh2,mch1,mch2)
!=============================================================================
INTEGER,  INTENT(IN) :: m1, m2, mah1, mah2, mbh1, mbh2, mch1, mch2
REAL(8),     INTENT(IN) :: a(m1,-mah1:mah2), b(m2,-mbh1:mbh2)
REAL(8),     INTENT(OUT):: c(m1,-mch1:mch2)
INTEGER              :: nch1, nch2, j, k, jpk, i1,i2
nch1=mah1+mbh1; nch2=mah2+mbh2
IF(nch1 /= mch1 .OR. nch2 /= mch2)STOP 'In MSBBB, dimensions inconsistent'
DO j=-mah1,mah2
   DO k=-mbh1,mbh2; jpk=j+k; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
      c(i1:i2,jpk)=c(i1:i2,jpk)-a(i1:i2,j)*b(j+i1:j+i2,k)
   ENDDO
ENDDO
END SUBROUTINE msbbb

!=============================================================================
SUBROUTINE LDUB(a,m,mah1,mah2)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE LDUB
!  Compute [L]*[D**-1]*[U] decomposition of asymmetric band-matrix
!
! <-> A: input as the asymmetric band matrix. On output, it contains
!     the [L]*[D**-1]*[U] factorization of the input matrix, where
!     [L] is lower triangular with unit main diagonal
!     [D] is a diagonal matrix
!     [U] is upper triangular with unit main diagonal
! --> M:    The number of rows of array A
! --> MAH1: the left half-bandwidth of fortran array A
! --> MAH2: the right half-bandwidth of fortran array A
!=============================================================================
INTEGER, INTENT(IN)   :: m,mah1, mah2 
REAL(8),    INTENT(INOUT):: a(m,-mah1:mah2) 
INTEGER               :: j, imost, jmost, jp, i
REAL(8)                  :: ajj, ajji, aij
DO j=1,m
  imost=MIN(m,j+mah1)
  jmost=MIN(m,j+mah2)
  jp=j+1
  ajj=a(j,0)
  IF(ajj == 0.)THEN
    PRINT '(" Failure in LDUB:"/" Matrix requires pivoting or is singular")'
    STOP
  ENDIF
  ajji=1./ajj
  a(j,0)=ajji
  DO i=jp,imost
    aij=ajji*a(i,j-i)
    a(i,j-i)=aij
    a(i,jp-i:jmost-i)=a(i,jp-i:jmost-i)-aij*a(j,jp-j:jmost-j)
  ENDDO
  a(j,jp-j:jmost-j)=ajji*a(j,jp-j:jmost-j)
ENDDO
END SUBROUTINE LDUB

!=============================================================================
SUBROUTINE DLDUB(a,m,mah1,mah2)
!=============================================================================
INTEGER,  INTENT(IN)   :: m,mah1, mah2 
REAL(8),  INTENT(INOUT):: a(m,-mah1:mah2) 
INTEGER                :: j, imost, jmost, jp, i
REAL(8)                :: ajj, ajji, aij
DO j=1,m
  imost=MIN(m,j+mah1)
  jmost=MIN(m,j+mah2)
  jp=j+1
  ajj=a(j,0)
  IF(ajj == 0)THEN
    PRINT '(" Fails in LDUB_d:"/" Matrix requires pivoting or is singular")'
    STOP
  ENDIF
  ajji=1./ajj
  a(j,0)=ajji
  DO i=jp,imost
    aij=ajji*a(i,j-i)
    a(i,j-i)=aij
    a(i,jp-i:jmost-i)=a(i,jp-i:jmost-i)-aij*a(j,jp-j:jmost-j)
  ENDDO
  a(j,jp-j:jmost-j)=ajji*a(j,jp-j:jmost-j)
ENDDO
END SUBROUTINE DLDUB

!=============================================================================
SUBROUTINE L1UBB(a,b,m,mah1,mah2,mbh1,mbh2)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1996
!		    SUBROUTINE L1UBB
!  Form the [L]*[D]*[U] decomposition of asymmetric band-matrix  [A] replace
!  lower triangular elements of [A] by [D**-1]*[L]*[D], the upper by [U],
!  replace matrix [B] by [D**-1]*[B].
!
! <-> A input as band matrix, output as lower and upper triangulars with 1s
!     implicitly assumed to lie on the main diagonal. The product of these
!     triangular matrices is [D**-1]*[A], where [D] is a diagonal matrix.
! <-> B in as band matrix, out as same but premultiplied by diagonal [D**-1]
! --> M    Number of rows of A and B
! --> MAH1 left half-width of fortran array A
! --> MAH2 right half-width of fortran array A
! --> MBH1 left half-width of fortran array B
! --> MBH2 right half-width of fortran array B
!=============================================================================
INTEGER, INTENT(IN) ::  m,mah1, mah2, mbh1, mbh2 
REAL(8), INTENT(INOUT) :: a(m,-mah1:mah2), b(m,-mbh1:mbh2)
INTEGER             :: j, imost, jmost, jleast, jp, i
REAL(8)                :: ajj, ajji, aij
DO j=1,m
  imost=MIN(m,j+mah1)
  jmost=MIN(m,j+mah2)
  jleast=MAX(1,j-mah1)
  jp=j+1
  ajj=a(j,0)
  IF(ajj == 0.)STOP 'failure in L1UBB'
  ajji=1./ajj
  a(j,jleast-j:jmost-j) = ajji * a(j,jleast-j:jmost-j)
  DO i=jp,imost
    aij=a(i,j-i)
    a(i,jp-i:jmost-i) = a(i,jp-i:jmost-i) - aij*a(j,jp-j:jmost-j)
  ENDDO
  a(j,0)=1.
  b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
ENDDO
END SUBROUTINE L1UBB

!=============================================================================
SUBROUTINE DL1UBB(a,b,m,mah1,mah2,mbh1,mbh2)
!=============================================================================
INTEGER                :: m,j, imost, jmost, jleast, jp, i
INTEGER,  INTENT(IN)   ::  mah1, mah2, mbh1, mbh2 
REAL(8),  INTENT(INOUT):: a(m,-mah1:mah2), b(m,-mbh1:mbh2)
REAL(8)                :: ajj, ajji, aij
DO j=1,m
  imost=MIN(m,j+mah1)
  jmost=MIN(m,j+mah2)
  jleast=MAX(1,j-mah1)
  jp=j+1
  ajj=a(j,0)
  IF(ajj == 0)STOP 'failure in DL1UBB'
  AJJI=1./AJJ
  a(j,jleast-j:jmost-j) = ajji * a(j,jleast-j:jmost-j)
  DO I=JP,IMOST
    AIJ=A(I,J-I)
    a(i,jp-i:jmost-i) = a(i,jp-i:jmost-i) - aij*a(j,jp-j:jmost-j)
  ENDDO
  A(J,0)=1.
  b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
ENDDO
END SUBROUTINE DL1UBB

!=============================================================================
SUBROUTINE l1ueb(a,b,m,mah1,mah2,mbh1,mbh2)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1998
!		    SUBROUTINE L1UEB
!  Form the [L]*[D]*[U] decomposition of asymmetric band-matrix  [A] replace
!  all but row zero of the
!  lower triangular elements of [A] by [D**-1]*[L]*[D], the upper by [U],
!  replace matrix [B] by [D**-1]*[B].
!  This is a special adaptation of L1UBB used to process quadarature weights
!  for QEDBV etc in which the initial quadrature value is provided as input
!  instead of being implicitly assumed zero (which is the case for QZDBV etc).
!
! <-> A input as band matrix, output as lower and upper triangulars with 1s
!     implicitly assumed to lie on the main diagonal. The product of these
!     triangular matrices is [D**-1]*[A], where [D] is a diagonal matrix.
! <-> B in as band matrix, out as same but premultiplied by diagonal [D**-1]
! --> M    number of rows of B, one less than the rows of A (which has "row 0")
! --> MAH1 left half-width of fortran array A
! --> MAH2 right half-width of fortran array A
! --> MBH1 left half-width of fortran array B
! --> MBH2 right half-width of fortran array B
!=============================================================================
INTEGER, INTENT(IN) :: m,mah1, mah2, mbh1, mbh2 
REAL(8), INTENT(INOUT) :: a(0:m,-mah1:mah2), b(m,-mbh1:mbh2)
INTEGER :: j, imost, jmost, jleast, jp, i
REAL(8) :: ajj, ajji, aij
DO j=1,m
  imost=MIN(m,j+mah1)
  jmost=MIN(m,j+mah2)
  jleast=MAX(0,j-mah1)
  jp=j+1
  ajj=a(j,0)
  IF(ajj == 0.)STOP 'failure in L1UEB'
  ajji=1./ajj
  a(j,jleast-j:jmost-j) = ajji * a(j,jleast-j:jmost-j)
  DO i=jp,imost
    aij=a(i,j-i)
    a(i,jp-i:jmost-i) = a(i,jp-i:jmost-i) - aij*a(j,jp-j:jmost-j)
  ENDDO
  a(j,0)=1.
  b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
ENDDO
END SUBROUTINE l1ueb

!=============================================================================
SUBROUTINE dl1ueb(a,b,m,mah1,mah2,mbh1,mbh2)
!=============================================================================
INTEGER,  INTENT(IN)   :: m,mah1, mah2, mbh1, mbh2 
REAL(8),  INTENT(INOUT):: a(0:,-mah1:), b(:,-mbh1:)
INTEGER                :: j, imost, jmost, jleast, jp, i
REAL(8)                :: ajj, ajji, aij
DO j=1,m
  imost=MIN(m,j+mah1)
  jmost=MIN(m,j+mah2)
  jleast=MAX(0,j-mah1)
  jp=j+1
  ajj=a(j,0)
  IF(ajj == 0)STOP 'failure in L1UEB_d'
  ajji=1./ajj
  a(j,jleast-j:jmost-j) = ajji * a(j,jleast-j:jmost-j)
  DO i=jp,imost
    aij=a(i,j-i)
    a(i,jp-i:jmost-i) = a(i,jp-i:jmost-i) - aij*a(j,jp-j:jmost-j)
  ENDDO
  a(j,0)=1.
  b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
ENDDO
END SUBROUTINE dl1ueb

!=============================================================================
SUBROUTINE L1LB(a,b,m,mah)	! Cholesky LU decomposition of Banded.
!=============================================================================
INTEGER,  INTENT(IN) :: m, mah
REAL(8),     INTENT(IN) :: a(m,-mah:mah)
REAL(8),     INTENT(OUT):: b(m,-mah:0)
INTEGER              :: i, j,jmi
REAL(8)                 :: s
CALL clib(b,m,m,mah,0)
DO j=1,m
   s=a(j,0)-DOT_PRODUCT(b(j,-mah:-1),b(j,-mah:-1))
   IF(s <= 0.)THEN
      PRINT '(" L1LB detects non-positivity at diagonal index",i5)',j
      STOP
   ENDIF
   s=SQRT(s); b(j,0)=s; s=1./s
   DO i=j+1,MIN(m,j+mah); jmi=j-i
      b(i,jmi)=s*(a(i,jmi)-DOT_PRODUCT(b(i,-mah:jmi-1),b(j,-mah-jmi:-1)))
   ENDDO
ENDDO
END SUBROUTINE L1LB

!=============================================================================
SUBROUTINE LDLB(a,b,d,m,mah) ! Modified Cholesky [L(D**-1)U, without sqrt]
!=============================================================================
INTEGER,  INTENT(IN) :: m, mah
REAL(8),     INTENT(IN) :: a(m,-mah:mah)
REAL(8),     INTENT(OUT):: b(m,-mah:0)
REAL(8),     INTENT(OUT):: d(m) 
INTEGER              :: i, j,k,jmi,lj,li
REAL(8)                 :: s,t
CALL clib(b,m,m,mah,0); b(:,0)=1.
DO j=1,m; lj=MAX(-mah,1-j)
   s=a(j,0)
   do k=lj,-1
      s=s-b(j,k)**2*d(k+j)
   enddo
   IF(s <= 0.)THEN
      PRINT '(" LDLB detects non-positivity at diagonal index",i5)',j
      STOP
   ENDIF
   d(j)=s; s=1./s
   DO i=j+1,MIN(m,j+mah); jmi=j-i; li=MAX(-mah,1-i); lj=li-jmi
      t=a(i,jmi)
      do k=li,jmi-1
         t=t-b(i,k)*b(j,k-jmi)*d(i+k)
      enddo
      b(i,jmi)=s*t
   ENDDO
ENDDO
d=1./d
END SUBROUTINE LDLB

!=============================================================================
SUBROUTINE DLDLB(a,b,d,m,mah) ! Modified Cholesky [L(D**-1)U, without sqrt]
!=============================================================================
INTEGER,  INTENT(IN) :: m, mah
REAL(8),  INTENT(IN) :: a(m,-mah:mah)
REAL(8),  INTENT(OUT):: b(m,-mah:0)
REAL(8),  INTENT(OUT):: d(m) 
INTEGER              :: i, j,k,jmi,lj,li
REAL(8)              :: s,t
CALL clib_d(b,m,m,mah,0); b(:,0)=1.
DO j=1,m; lj=MAX(-mah,1-j)
   s=a(j,0)
   do k=lj,-1
      s=s-b(j,k)**2*d(k+j)
   enddo
   IF(s <= 0.)THEN
      PRINT '(" DLDLB detects non-positivity at diagonal index",i5)',j
      STOP
   ENDIF
   d(j)=s; s=1./s
   DO i=j+1,MIN(m,j+mah); jmi=j-i;  
      li=MAX(-mah,1-i); 
      lj=li-jmi; 
      t=a(i,jmi)
      do k=li,jmi-1
         t=t-b(i,k)*b(j,k-jmi)*d(i+k)
      enddo
      b(i,jmi)=s*t
   ENDDO
ENDDO
d=1./d
END SUBROUTINE DLDLB

!=============================================================================
SUBROUTINE UDUB(a,b,d,m,mah) ! Modified reverse Cholesky [U(D**-1)U^t],
!=============================================================================
INTEGER,        INTENT(IN) :: m, mah
REAL(8),           INTENT(IN) :: a(m,-mah:mah)
REAL(8),           INTENT(OUT):: b(m,0:mah)
REAL(8),           INTENT(OUT):: d(m) 
REAL(8), DIMENSION(m,-mah:mah):: at
REAL(8), DIMENSION(m,-mah:0)  :: bt
REAL(8), DIMENSION(m)         :: dt
at=a(m:1:-1,mah:-mah:-1); CALL ldlb(at,bt,dt,m,mah);
b=bt(m:1:-1,0:-mah:-1); d=dt(m:1:-1)
END SUBROUTINE UDUB

!=============================================================================
SUBROUTINE DUDUB(a,b,d,m,mah) ! Modified reverse Cholesky [U(D**-1)U^t],
!=============================================================================
INTEGER,           INTENT(IN) :: m, mah
REAL(8),           INTENT(IN) :: a(m,-mah:mah)
REAL(8),           INTENT(OUT):: b(m,0:mah)
REAL(8),           INTENT(OUT):: d(m) 
REAL(8), DIMENSION(m,-mah:mah):: at
REAL(8), DIMENSION(m,-mah:0)  :: bt
REAL(8), DIMENSION(m)         :: dt
at=a(m:1:-1,mah:-mah:-1); CALL ldlb_d(at,bt,dt,m,mah);
b=bt(m:1:-1,0:-mah:-1);   d=dt(m:1:-1)
END SUBROUTINE DUDUB

!=============================================================================
SUBROUTINE mulbv(a,v1,v2, m1,m2,mah1,mah2)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE MULBV
!  MULtipication of a Banded matrix times a Vector.
!
! --> A is the matrix
! --> V1 the input vector
! <-- V2 the output vector
! --> M1 the number of rows assumed for A and for V2
! --> M2 the number of columns assumed for A and rows for V1
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
!=============================================================================
INTEGER,  INTENT(IN) :: m1, m2, mah1, mah2
REAL(8),     INTENT(IN) :: a(m1,-mah1:mah2), v1(m2)
REAL(8),     INTENT(OUT):: v2(m1)
INTEGER              :: j, i1,i2 
v2 = 0.0
!=============================================================================
ENTRY	 madbv(a,v1,v2, m1,m2,mah1,mah2)
!=============================================================================
DO j=-mah1,mah2; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
   v2(i1:i2) = v2(i1:i2) + a(i1:i2,j)*v1(j+i1:j+i2)
ENDDO
RETURN
!=============================================================================
ENTRY	 msbbv(a,v1,v2, m1,m2,mah1,mah2)
!=============================================================================
DO j=-mah1,mah2; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
   v2(i1:i2) = v2(i1:i2) - a(i1:i2,j)*v1(j+i1:j+i2)
ENDDO
END SUBROUTINE mulbv

!=============================================================================
SUBROUTINE mulbx(a,v1,v2, m1,m2,mah1,mah2,my)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE MULBX
!  MULtipication of a Banded matrix times parallel X-Vectors.
!
! --> A is the matrix
! --> V1 the array of input vectors
! <-- V2 the array of output vectors
! --> M1 the number of rows assumed for A and for V2
! --> M2 the number of columns assumed for A and rows for V1
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
! --> MY the number of parallel X-vectors
!=============================================================================
INTEGER,  INTENT(IN) :: m1, m2, mah1, mah2, my
REAL(8),     INTENT(IN) :: a(m1,-mah1:mah2), v1(m2,my)
REAL(8),     INTENT(OUT):: v2(m1,my)
INTEGER              :: i,j
v2=0.0
!=============================================================================
ENTRY	 madbx(a,v1,v2, m1,m2,mah1,mah2,my)
!=============================================================================
DO j=-mah1,mah2
   DO i=MAX(1,1-j),MIN(m1,m2-j); v2(i,:)=v2(i,:)+a(i,j)*v1(i+j,:); ENDDO
ENDDO
RETURN
!=============================================================================
ENTRY	 msbbx(a,v1,v2, m1,m2,mah1,mah2,my)
!=============================================================================
DO j=-mah1,mah2
   DO i=MAX(1,1-j),MIN(m1,m2-j); v2(i,:)=v2(i,:)-a(i,j)*v1(i+j,:); ENDDO
ENDDO
END SUBROUTINE mulbx

!=============================================================================
SUBROUTINE mulby(a,v1,v2, m1,m2,mah1,mah2,mx)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE MULBY
!  MULtipication of a Banded matrix times parallel Y-Vectors.
!
! --> A is the matrix
! --> V1 the array of input vectors
! <-- V2 the array of output vectors
! --> M1 the number of rows assumed for A and for V2
! --> M2 the number of columns assumed for A and rows for V1
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
! --> MX the length of each of the parallel Y-vectors
!=============================================================================
INTEGER,  INTENT(IN) :: m1, m2, mah1, mah2, mx
REAL(8),     INTENT(IN) :: a(m1,-mah1:mah2), v1(mx,m2)
REAL(8),     INTENT(OUT):: v2(mx,m1)
INTEGER              :: i,j
v2(1:mx,1:m1) = 0.0
ENTRY	 madby(a,v1,v2, m1,m2,mah1,mah2,mx)
DO j=-mah1,mah2
   DO i=MAX(1,1-j),MIN(m1,m2-j); v2(:,i)=v2(:,i)+a(i,j)*v1(:,i+j); ENDDO
ENDDO
RETURN
ENTRY	 msbby(a,v1,v2, m1,m2,mah1,mah2,mx)
DO j=-mah1,mah2
   DO i=MAX(1,1-j),MIN(m1,m2-j); v2(:,i)=v2(:,i)-a(i,j)*v1(:,i+j); ENDDO
ENDDO
END SUBROUTINE mulby

!=============================================================================
SUBROUTINE MULVB(v1,a,v2, m1,m2,mah1,mah2)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE MULVB
!  MULtipication of a Vector times a Banded matrix.
!
! --> V1 the input row-vector
! --> A is the matrix
! <-- V2 the output vector
! --> M1 the number of rows assumed for A and columns for V1
! --> M2 the number of columns assumed for A and for V2
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
!=============================================================================
INTEGER,  INTENT(IN) :: m1, m2, mah1, mah2
REAL(8),     INTENT(IN) :: v1(m1), a(m1,-mah1:mah2)
REAL(8),     INTENT(OUT):: v2(m2)
INTEGER              :: j, i1,i2
v2=0.0
!=============================================================================
ENTRY	 madvb(v1,a,v2, m1,m2,mah1,mah2)
!=============================================================================
DO j=-mah1,mah2; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
   v2(j+i1:j+i2)=v2(j+i1:j+i2)+v1(i1:i2)*a(i1:i2,j)
ENDDO
RETURN
!=============================================================================
ENTRY	 msbvb(v1,a,v2, m1,m2,mah1,mah2)
!=============================================================================
DO j=-mah1,mah2; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
   v2(j+i1:j+i2)=v2(j+i1:j+i2)-v1(i1:i2)*a(i1:i2,j)
ENDDO
END SUBROUTINE mulvb

!=============================================================================
SUBROUTINE mulxb(v1,a,v2, m1,m2,mah1,mah2,my)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE MULXB
!  MULtipication of X-Vectors times Banded matrix.
!
! --> V1 the array of input row-vectors
! --> A is the matrix
! <-- V2 the array of output vectors
! --> M1 the number of rows assumed for A and columns for V1
! --> M2 the number of columns assumed for A and V2
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
! --> MY the number of parallel X-vectors
!=============================================================================
INTEGER,  INTENT(IN) :: m1, m2, mah1, mah2, my
REAL(8),     INTENT(IN) :: v1(m1,my), a(m1,-mah1:mah2)
REAL(8),     INTENT(OUT):: v2(m2,my)
INTEGER              :: i,j
v2=0.0
!=============================================================================
ENTRY	 madxb(v1,a,v2, m1,m2,mah1,mah2,my)
!=============================================================================
DO j=-mah1,mah2
   DO i=MAX(1,1-j),MIN(m1,m2-j); v2(j+i,:)=v2(j+i,:)+v1(i,:)*a(i,j); ENDDO
ENDDO
RETURN
!=============================================================================
ENTRY	 msbxb(v1,a,v2, m1,m2,mah1,mah2,my)
!=============================================================================
DO j=-mah1,mah2
   DO i=MAX(1,1-j),MIN(m1,m2-j); v2(j+i,:)=v2(j+i,:)-v1(i,:)*a(i,j); ENDDO
ENDDO
END SUBROUTINE mulxb

!=============================================================================
SUBROUTINE mulyb(v1,a,v2, m1,m2,mah1,mah2,mx)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE MULYB
!  MULtipication of Y-Vectors times a Banded matrix.
!
! --> V1 the array of input row-vectors
! --> A is the matrix
! <-- V2 the array of output vectors
! --> M1 the number of rows assumed for A and columns for V1
! --> M2 the number of columns assumed for A and V2
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
! --> MX the length of each of the parallel Y-vectors
!=============================================================================
INTEGER,  INTENT(IN) :: m1, m2, mah1, mah2, mx
REAL(8),     INTENT(IN) :: v1(mx,m1), a(m1,-mah1:mah2)
REAL(8),     INTENT(OUT):: v2(mx,m2)
INTEGER              :: i,j
v2=0.0
ENTRY	 madyb(v1,a,v2, m1,m2,mah1,mah2,mx)
DO j=-mah1,mah2
   DO i=MAX(1,1-j),MIN(m1,m2-j); v2(:,j+i)=v2(:,j+i)+v1(:,i)*a(i,j); ENDDO
ENDDO
RETURN
ENTRY	 msbyb(v1,a,v2, m1,m2,mah1,mah2,mx)
DO j=-mah1,mah2
   DO i=MAX(1,1-j),MIN(m1,m2-j); v2(:,j+i)=v2(:,j+i)-v1(:,i)*a(i,j); ENDDO
ENDDO
END SUBROUTINE mulyb

!=============================================================================
SUBROUTINE mulbd(a,d,b,m1,m2,mah1,mah2)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE MULBD
! MULtipication of a Banded matrix times a Diagonal
!
! --> A is the input banded-matrix
! --> D the diagonal matrix
! <-- B the output matrix
! --> M1 the number of rows assumed for A and for B
! --> M2 number of columns assumed for A and B, number of elements of D
! --> MAH1 the left half-bandwidth of arrays A and B
! --> MAH2 the right half-bandwidth of arrays A and B
!=============================================================================
INTEGER,  INTENT(IN   ):: m1, m2, mah1, mah2
REAL(8),     INTENT(IN   ):: d(m2)
REAL(8),     INTENT(INOUT):: a(m1,-mah1:mah2),b(m1,-mah1:mah2)
INTEGER                :: j, i1,i2
CALL clib(b,m1,m2,mah1,mah2)
DO j=-mah1,mah2; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
   b(i1:i2,j)=a(i1:i2,j)*d(j+i1:j+i2)
ENDDO
RETURN
!=============================================================================
ENTRY	 madbd(a,d,b,m1,m2,mah1,mah2)
!=============================================================================
DO j=-mah1,mah2; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
   b(i1:i2,j) = b(i1:i2,j)+a(i1:i2,j)*d(j+i1:j+i2)
ENDDO
RETURN
!=============================================================================
ENTRY	 msbbd(a,d,b,m1,m2,mah1,mah2)
!=============================================================================
DO j=-mah1,mah2; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
   b(i1:i2,j) = b(i1:i2,j)-a(i1:i2,j)*d(j+i1:j+i2)
ENDDO
END SUBROUTINE mulbd

!=============================================================================
SUBROUTINE muldb(d,a,b,m1,m2,mah1,mah2)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE MULDB
!  MULtipication of a Banded matrix times a Diagonal
!
! --> D the diagonal matrix
! --> A is the input banded-matrix ! <->  if A and B are actually
! <-- B the output matrix          ! <->  equivalent arrays.
! --> M1 the number of rows assumed for A and for B
! --> M2 number of columns assumed for A and B, number of elements of D
! --> MAH1 the left half-bandwidth of arrays A and B
! --> MAH2 the right half-bandwidth of arrays A and B
!=============================================================================
INTEGER,  INTENT(IN)    :: m1, m2, mah1, mah2
REAL(8),     INTENT(IN   ) :: d(m1)
REAL(8),     INTENT(INOUT) :: a(m1,-mah1:mah2),b(m1,-mah1:mah2)
INTEGER                 :: j
CALL clib(b,m1,m2,mah1,mah2)
DO j=-mah1,mah2; b(:,j)=d(:)*a(:,j); ENDDO
END SUBROUTINE muldb

!=============================================================================
SUBROUTINE maddb(d,a,b,m1,m2,mah1,mah2)
!=============================================================================
INTEGER,  INTENT(IN)    :: m1, m2, mah1, mah2
REAL(8),     INTENT(IN   ) :: d(m1)
REAL(8),     INTENT(INOUT) :: a(m1,-mah1:mah2),b(m1,-mah1:mah2)
INTEGER                 :: j
DO j=-mah1,mah2; b(:,j)=b(:,j)+d(:)*a(:,j); ENDDO
END SUBROUTINE maddb

!=============================================================================
SUBROUTINE msbdb(d,a,b,m1,m2,mah1,mah2)
!=============================================================================
INTEGER,  INTENT(IN)    :: m1, m2, mah1, mah2
REAL(8),     INTENT(IN   ) :: d(m1) 
REAL(8),     INTENT(INOUT) :: a(m1,-mah1:mah2),b(m1,-mah1:mah2)
INTEGER                 :: j
DO j=-mah1,mah2; b(:,j)=b(:,j)-d(:)*a(:,j); ENDDO
END SUBROUTINE msbdb


!=============================================================================
SUBROUTINE udlbv(a,v, m,mah1,mah2)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE UDLBV
!  BACk-substitution step of linear inversion involving
!  Banded matrix and Vector.
!
! --> A encodes the (L)*(D**-1)*(U) factorization of the linear-system
!     matrix, as supplied by subroutine LDUB
! <-> V input as right-hand-side vector, output as solution vector
! --> M the number of rows assumed for A and for V
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
!=============================================================================
INTEGER,  INTENT(IN)   :: m, mah1, mah2
REAL(8),     INTENT(IN)   :: a(m,-mah1:mah2)
REAL(8),     INTENT(INOUT):: v(m)
INTEGER                :: i, j
REAL(8)                   :: vj
DO j=1,m
   vj=v(j)
   DO i=j+1,MIN(m,j+mah1); v(i)=v(i)-a(i,j-i)*vj; ENDDO; v(j)=a(j,0)*vj
ENDDO
DO j=m,2,-1
   vj=v(j)
   DO i=MAX(1,j-mah2),j-1; v(i)=v(i)-a(i,j-i)*vj; ENDDO
ENDDO
END SUBROUTINE udlbv

!=============================================================================
SUBROUTINE udlbx(a,v, mx,mah1,mah2,my)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE UDLBX
!  BACk-substitution step of parallel linear inversion involving
!  Banded matrix and X-Vectors.
!
! --> A encodes the (L)*(D**-1)*(U) factorization of the linear-system
!     matrix, as supplied by subroutine LDUB or, if N=NA, by LDUB
! <-> V input as right-hand-side vectors, output as solution vectors
! --> MX the number of rows assumed for A and length of
!     X-vectors stored in V
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
! --> MY number of parallel X-vectors inverted
!=============================================================================
INTEGER,  INTENT(IN)   :: mx, mah1, mah2, my
REAL(8),     INTENT(IN)   :: a(mx,-mah1:mah2)
REAL(8),     INTENT(INOUT):: v(mx,my)
INTEGER                :: jx, ix
DO jx=1,mx
   DO ix=jx+1,MIN(mx,jx+mah1); v(ix,:) = v(ix,:) - a(ix,jx-ix)*v(jx,:); ENDDO
   v(jx,:) = a(jx,0) * v(jx,:)
ENDDO
DO jx=mx,2,-1
   DO ix=MAX(1,jx-mah2),jx-1; v(ix,:) = v(ix,:) - a(ix,jx-ix)*v(jx,:); ENDDO
ENDDO
END SUBROUTINE udlbx

!=============================================================================
SUBROUTINE udlby(a,v, my,mah1,mah2,mx)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE UDLBY
!  BACk-substitution step of parallel linear inversion involving
!  Banded matrix and Y-Vectors.
!
! --> A encodes the (L)*(D**-1)*(U) factorization of the linear-system
!     matrix, as supplied by subroutine LDUB or, if N=NA, by LDUB
! <-> V input as right-hand-side vectors, output as solution vectors
! --> MY the number of rows assumed for A and length of
!     Y-vectors stored in V
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
! --> MX number of parallel Y-vectors inverted
!=============================================================================
INTEGER,  INTENT(IN)   :: my, mah1, mah2, mx
REAL(8),     INTENT(IN)   :: a(my,-mah1:mah2)
REAL(8),     INTENT(INOUT):: v(mx,my)
INTEGER                :: iy, jy
DO jy=1,my
   DO iy=jy+1,MIN(my,jy+mah1); v(:,iy) = v(:,iy)-a(iy,jy-iy)*v(:,jy); ENDDO
   v(:,jy)=a(jy,0)*v(:,jy)
ENDDO
DO jy=my,2,-1
   DO iy=MAX(1,jy-mah2),jy-1; v(:,iy)=v(:,iy)-a(iy,jy-iy)*v(:,jy); ENDDO
ENDDO
END SUBROUTINE udlby

!=============================================================================
SUBROUTINE udlvb(v,a, m,mah1,mah2)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE UDLVB
!  BACk-substitution step of linear inversion involving
!  row-Vector and Banded matrix.
!
! <-> V input as right-hand-side row-vector, output as solution vector
! --> A encodes the (L)*(D**-1)*(U) factorization of the linear-system
!     matrix, as supplied by subroutine LDUB
! --> M the number of rows assumed for A and columns for V
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
!=============================================================================
INTEGER,  INTENT(IN)   :: m, mah1, mah2
REAL(8),     INTENT(IN)   :: a(m,-mah1:mah2)
REAL(8),     INTENT(INOUT):: v(m)
INTEGER                :: i, j
REAL(8)                   :: vi
DO i=1,m
   vi=v(i)
   DO j=i+1,MIN(m,i+mah2); v(j)=v(j)-vi*a(i,j-i); ENDDO
   v(i)=vi*a(i,0)
ENDDO
DO i=m,2,-1
   vi=v(i)
   DO j=MAX(1,i-mah1),i-1; v(j)=v(j)-vi*a(i,j-i); ENDDO
ENDDO
END SUBROUTINE udlvb

!=============================================================================
SUBROUTINE udlxb(v,a, mx,mah1,mah2,my)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE UDLXB
!  BACk-substitution step of parallel linear inversion involving
!  Banded matrix and row-X-Vectors.
!
! <-> V input as right-hand-side vectors, output as solution vectors
! --> A encodes the (L)*(D**-1)*(U) factorization of the linear-system
!     matrix, as supplied by subroutine LDUB
! --> MX the number of rows assumed for A and length of
!     X-vectors stored in V
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
! --> MY number of parallel X-vectors inverted
!=============================================================================
INTEGER,  INTENT(IN)   :: mx, mah1, mah2, my
REAL(8),     INTENT(IN)   :: a(mx,-mah1:mah2)
REAL(8),     INTENT(INOUT):: v(mx,my)
INTEGER                :: ix, jx
DO ix=1,mx
   DO jx=ix+1,MIN(mx,ix+mah2); v(jx,:)=v(jx,:)-v(ix,:)*a(ix,jx-ix); ENDDO
   v(ix,:)=v(ix,:)*a(ix,0)
ENDDO
DO ix=mx,2,-1
   DO jx=MAX(1,ix-mah1),ix-1; v(jx,:)=v(jx,:)-v(ix,:)*a(ix,jx-ix); ENDDO
ENDDO
END SUBROUTINE udlxb

!=============================================================================
SUBROUTINE udlyb(v,a, my,mah1,mah2,mx)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE UDLYB
!  BACk-substitution step of parallel linear inversion involving
!  Banded matrix and row-Y-Vectors.
!
! <-> V input as right-hand-side vectors, output as solution vectors
! --> A encodes the (L)*(D**-1)*(U) factorization of the linear-system
!     matrix, as supplied by subroutine LDUB
! --> MY the number of rows assumed for A and length of
!     Y-vectors stored in V
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
! --> MX number of parallel Y-vectors inverted
!=============================================================================
INTEGER,  INTENT(IN)   :: my, mah1, mah2, mx
REAL(8),     INTENT(IN)   :: a(my,-mah1:mah2)
REAL(8),     INTENT(INOUT):: v(mx,my)
INTEGER                :: iy, jy
DO iy=1,my
   DO jy=iy+1,MIN(my,iy+mah2); v(:,jy)=v(:,jy)-v(:,iy)*a(iy,jy-iy); ENDDO
   v(:,iy)=v(:,iy)*a(iy,0)
ENDDO
DO iy=my,2,-1
   DO jy=MAX(1,iy-mah1),iy-1; v(:,jy)=v(:,jy)-v(:,iy)*a(iy,jy-iy); ENDDO
ENDDO
END SUBROUTINE udlyb

!=============================================================================
SUBROUTINE u1lbv(a,v, m,mah1,mah2)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1996
!		    SUBROUTINE U1LBV
!  BACk-substitution step ((U**-1)*(L**-1)) of linear inversion involving
!  special Banded matrix and right-Vector.
!
! --> A encodes the [L]*[U] factorization of the linear-system
!     matrix, as supplied by subroutine L1UBB
! <-> V input as right-hand-side vector, output as solution vector
! --> M the number of rows assumed for A and for V
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
!=============================================================================
INTEGER,  INTENT(IN)   :: m, mah1, mah2
REAL(8),     INTENT(IN)   :: a(m,-mah1:mah2)
REAL(8),     INTENT(INOUT):: v(m)
INTEGER                :: i, j
REAL(8)                   :: vj
DO j=1,m
   vj=v(j)
   DO i=j+1,MIN(m,j+mah1); v(i)=v(i)-a(i,j-i)*vj; ENDDO
ENDDO
DO j=m,2,-1
   vj=v(j)
   DO i=MAX(1,j-mah2),j-1; v(i)=v(i)-a(i,j-i)*vj; ENDDO
ENDDO
END SUBROUTINE u1lbv

!=============================================================================
SUBROUTINE u1lbx(a,v, mx,mah1,mah2,my)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1996
!		    SUBROUTINE U1LBX
!  Special BaCk-substitution step of parallel linear inversion involving
!  Banded matrix and X-right-Vectors.
!
! --> A encodes the [L]*[U] factorization of the linear-system
!     matrix, as supplied by subroutine L1UBB
! <-> V input as right-hand-side vectors, output as solution vectors
! --> MX the number of rows assumed for A and length of
!     X-vectors stored in V
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
! --> MY number of parallel X-vectors inverted
!=============================================================================
INTEGER,  INTENT(IN)   :: mx, mah1, mah2, my
REAL(8),     INTENT(IN)   :: a(mx,-mah1:mah2)
REAL(8),     INTENT(INOUT):: v(mx,my)
INTEGER                :: ix, jx
DO jx=1,mx
   DO ix=jx+1,MIN(mx,jx+mah1); v(ix,:)=v(ix,:)-a(ix,jx-ix)*v(jx,:); ENDDO
ENDDO
DO jx=mx,2,-1
   DO ix=MAX(1,jx-mah2),jx-1; v(ix,:)=v(ix,:)-a(ix,jx-ix)*v(jx,:); ENDDO
ENDDO
END SUBROUTINE u1lbx

!=============================================================================
SUBROUTINE u1lby(a,v, my,mah1,mah2,mx)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1996
!		    SUBROUTINE U1LBY
!  Special BaCk-substitution step of parallel linear inversion involving
!  Banded matrix and Y-right-Vectors.
!
! --> A encodes the [L]*[U] factorization of the linear-system
!     matrix, as supplied by subroutine L1UBB
! <-> V input as right-hand-side vectors, output as solution vectors
! --> MY the number of rows assumed for A and length of
!     Y-vectors stored in V
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
! --> MX number of parallel Y-vectors inverted
!=============================================================================
INTEGER,  INTENT(IN)   :: my, mah1, mah2, mx
REAL(8),     INTENT(IN)   :: a(my,-mah1:mah2)
REAL(8),     INTENT(INOUT):: v(mx,my)
INTEGER                :: iy, jy
DO jy=1,my
   DO iy=jy+1,MIN(my,jy+mah1); v(:,iy)=v(:,iy)-a(iy,jy-iy)*v(:,jy); ENDDO
ENDDO
DO jy=my,2,-1
   DO iy=MAX(1,jy-mah2),jy-1; v(:,iy)=v(:,iy)-a(iy,jy-iy)*v(:,jy); ENDDO
ENDDO
END SUBROUTINE u1lby

!=============================================================================
SUBROUTINE u1lvb(v,a, m,mah1,mah2)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1996
!		    SUBROUTINE U1LVB
!  Special BaCk-substitution step of linear inversion involving
!  left-Vector and Banded matrix.
!
! <-> V input as right-hand-side row-vector, output as solution vector
! --> A encodes the special [L]*[U] factorization of the linear-system
!     matrix, as supplied by subroutine L1UBB
! --> M the number of rows assumed for A and columns for V
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
!=============================================================================
INTEGER, INTENT(IN)   :: m, mah1, mah2
REAL(8),    INTENT(IN)   :: a(m,-mah1:mah2)
REAL(8),    INTENT(INOUT):: v(m)
INTEGER               :: i, j
REAL(8)                  :: vi
DO i=1,m
   vi=v(i)
   DO j=i+1,MIN(m,i+mah2); v(j)=v(j)-vi*a(i,j-i); ENDDO
ENDDO
DO i=m,2,-1
   vi=v(i)
   DO j=MAX(1,i-mah1),i-1; v(j)=v(j)-vi*a(i,j-i); ENDDO
ENDDO
END SUBROUTINE u1lvb

!=============================================================================
SUBROUTINE u1lxb(v,a, mx,mah1,mah2,my)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1996
!		    SUBROUTINE U1LXB
!  Special BaCk-substitution step of parallel linear inversion involving
!  Banded matrix and X-left-Vectors.
!
! <-> V input as right-hand-side vectors, output as solution vectors
! --> A encodes the special [L]*[U] factorization of the linear-system
!     matrix, as supplied by subroutine L1UBB
! --> MX the number of rows assumed for A and length of
!     X-vectors stored in V
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
! --> MY number of parallel X-vectors inverted
!=============================================================================
INTEGER,  INTENT(IN)   :: mx, mah1, mah2, my
REAL(8),     INTENT(IN)   :: a(mx,-mah1:mah2)
REAL(8),     INTENT(INOUT):: v(mx,my)
INTEGER                :: ix, jx
DO ix=1,mx
   DO jx=ix+1,MIN(mx,ix+mah2); v(jx,:)=v(jx,:)-v(ix,:)*a(ix,jx-ix); ENDDO
ENDDO
DO ix=mx,2,-1
   DO jx=MAX(1,ix-mah1),ix-1;  v(jx,:)=v(jx,:)-v(ix,:)*a(ix,jx-ix); ENDDO
ENDDO
END SUBROUTINE u1lxb

!=============================================================================
SUBROUTINE u1lyb(v,a, my,mah1,mah2,mx)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1996
!		    SUBROUTINE U1LYB
!  Special BaCk-substitution step of parallel linear inversion involving
!  special Banded matrix and Y-left-Vectors.
!
! <-> V input as right-hand-side vectors, output as solution vectors
! --> A encodes the [L]*[U] factorization of the linear-system
!     matrix, as supplied by subroutine L1UBB
! --> MY the number of rows assumed for A and length of
!     Y-vectors stored in V
! --> MAH1 the left half-bandwidth of fortran array A
! --> MAH2 the right half-bandwidth of fortran array A
! --> MX number of parallel Y-vectors inverted
!=============================================================================
INTEGER,  INTENT(IN)   :: my, mah1, mah2, mx
REAL(8),     INTENT(IN)   :: a(my,-mah1:mah2)
REAL(8),     INTENT(INOUT):: v(mx,my)
INTEGER                :: iy, jy
DO iy=1,my
   DO jy=iy+1,MIN(my,iy+mah2); v(:,jy)=v(:,jy)-v(:,iy)*a(iy,jy-iy); ENDDO
ENDDO
DO iy=my,2,-1
   DO jy=MAX(1,iy-mah1),iy-1;  v(:,jy)=v(:,jy)-v(:,iy)*a(iy,jy-iy); ENDDO
ENDDO
END SUBROUTINE u1lyb

!=============================================================================
SUBROUTINE linbv(a,v,m,mah1,mah2)
!=============================================================================
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994
!		    SUBROUTINE LINBV
!   Solve LINear system with square Banded-matrix and vector V
!
! <-> A system matrix on input, its [L]*[D**-1]*[U] factorization on exit
! <-> V vector of right-hand-sides on input, solution vector on exit
! --> M order of matrix A
! --> MAH1 left half-bandwidth of A
! --> MAH2 right half-bandwidth of A
!=============================================================================
INTEGER, INTENT(IN)    :: m, mah1, mah2
REAL(8),    INTENT(INOUT) :: a(m,-mah1:mah2), v(m)
CALL ldub(a,m,mah1,mah2)
CALL udlbv(a,v,m,mah1,mah2)
END SUBROUTINE linbv

!=============================================================================
SUBROUTINE wrtb(a,m1,m2,mah1,mah2)
!=============================================================================
INTEGER,  INTENT(IN) :: m1, m2, mah1, mah2
REAL(8),     INTENT(IN) :: a(m1,-mah1:mah2)
INTEGER              :: i1, i2, i, j1, j2, j, nj1
DO i1=1,m1,20
   i2=MIN(i1+19,m1)
   PRINT '(7x,6(i2,10x))',(j,j=-mah1,mah2)
   DO i=i1,i2
      j1=MAX(-mah1,1-i)
      j2=MIN(mah2,m2-i)
      nj1=j1+mah1
      IF(nj1==0)PRINT '(1x,i3,6(1x,e11.5))',i,(a(i,j),j=j1,j2)
      IF(nj1==1)PRINT '(1x,i3,12x,5(1x,e11.5))',i,(a(i,j),j=j1,j2)
      IF(nj1==2)PRINT '(1x,i3,24x,4(1x,e11.5))',i,(a(i,j),j=j1,j2)
      IF(nj1==3)PRINT '(1x,i3,36x,3(1x,e11.5))',i,(a(i,j),j=j1,j2)
      IF(nj1==4)PRINT '(1x,i3,48x,2(1x,e11.5))',i,(a(i,j),j=j1,j2)
      IF(nj1==5)PRINT '(1x,i3,60x,1(1x,e11.5))',i,(a(i,j),j=j1,j2)
   ENDDO
   READ(*,*)
ENDDO
END SUBROUTINE wrtb

END MODULE MODULE_pmat2


!ENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDEND
!						********************
!						* module_pmat2.f90 *
!						* PURSER 1994/1999 *
!                                               *   FUJITA 1999    *
!						********************
!ENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDEND

!						**********************
!						* module_fitcons.f90 *
!						* PURSER 1994/1999   *
!                                               *   FUJITA 1999      *
!						**********************

!============================================================================
module vkind
!============================================================================
  integer, parameter :: vp=kind(1.0d0)
!  integer, parameter :: vp=kind(1.0)
end module vkind

!============================================================================
module module_fitcons
!============================================================================
use vkind
implicit none
integer,parameter             :: noh=3,    nohm=noh-1,   nohp=noh+1,&
                                 no=noh*2, nom=no-1,     nop=no+1,   nnit=7
real(vp),parameter            :: sigc=3._vp,  sigb=2._vp
real(vp),dimension(no)        :: hunit,q,wt,dwt
real(vp),dimension(nom)       :: hunit1,hunit2,q1,wt1,dwt1
real(vp),dimension(-noh:noh)  :: qco
real(vp),dimension(-1-noh:noh):: ico,dco
real(vp)                      :: rcrit,ldsig,ldsig4
!============================================================================

contains
!============================================================================
SUBROUTINE setq(q,x,n) 
!============================================================================
!                SUBROUTINE SETQ
! Precompute the N constant denominator factors of the N-point Lagrange
! polynomial interpolation formula.
!
! <-- Q:    The N denominator constants.
! --> X:    The N abscissae.
! --> N:    The number of points involved.
!============================================================================
  use vkind
  IMPLICIT NONE
  INTEGER,          INTENT(in) :: n
  REAL(vp),DIMENSION(n),INTENT(out):: q
  REAL(vp),DIMENSION(n),INTENT(in) :: x
!-----------------------------------------------------------------------------
  INTEGER                          :: i,j
!=============================================================================
DO i=1,n
   q(i)=1.
   DO j=1,n
      IF(j /= i)q(i)=q(i)/(x(i)-x(j))
   ENDDO
ENDDO
END SUBROUTINE setq 

!============================================================================
SUBROUTINE lagw(x,xt,q,w,dw,n) 
!============================================================================
!      SUBROUTINE LAGW
! Construct the Lagrange weights and their derivatives when target abscissa
! is known and denominators Q have already been precomputed
!
! --> X:    Grid abscissae
! --> XT:   Target abscissa
! --> Q:    Q factors (denominators of the Lagrange weight formula)
! <-- W:    Lagrange weights
! <-- DW:   Derivatives, dW/dX, of Lagrange weights W
! --> N:    Number of grid points involved in the interpolation
!============================================================================
  use vkind
  IMPLICIT NONE
  INTEGER,              INTENT(in) :: n
  REAL(vp),             INTENT(in) :: xt
  REAL(vp),DIMENSION(n),INTENT(in) :: x,q
  REAL(vp),DIMENSION(n),INTENT(out):: w,dw
!-----------------------------------------------------------------------------
  REAL(vp),DIMENSION(n)            :: sdit,d,di
  INTEGER                          :: i,j
  REAL(vp)                         :: p,s,sdil,sdir
!============================================================================
p=1.       ! ...will become product of all the d(i)=xt-x(i)
DO i=1,n
   d(i)=xt-x(i)
   p=p*d(i)
ENDDO

!   test p to reveal whether any of the d(i) vanish:
IF(p == 0._vp)THEN   ! xt coincides with a grid point - use special code:
   p=1.           ! p will become the product of the nonzero d(i),
   s=0.           ! s will become the corresponding sum of q(i)/d(i)
   DO i=1,n
      IF(d(i) == 0._vp)THEN
         j=i            ! identify the grid index corresponding to present xt
         w(j)=1.        ! interpolation weighted entirely to this one.
      ELSE
         w(i)=0.
         p=p*d(i)
         dw(i)=q(i)/d(i)
         s=s+dw(i)
      ENDIF
   ENDDO
   dw(j)=-s*p
   DO i=1,n
      IF(i /= j)dw(i)=dw(i)*p
   ENDDO
ELSE             ! xt is not a grid point - use generic code:
   sdil=0.            ! will become the sum of terms to the left.
   sdir=0.            ! will become the sum of terms to the right.
   DO i=1,n
      di(i)=1./d(i)
      sdit(i)=sdil
      sdil=sdil+di(i)
      w(i)=q(i)*p*di(i)
   ENDDO
   DO i=n,1,-1
      sdit(i)=sdit(i)+sdir
      sdir=sdir+di(i)
      dw(i)=w(i)*sdit(i)
   ENDDO
ENDIF
END SUBROUTINE lagw 

!============================================================================
subroutine infit
!============================================================================
implicit none
integer :: i,l
real(vp):: divq,divd
!============================================================================
! Initialize quantities that relate to interpolations:
do i=1,no; hunit(i)=i-noh; enddo
hunit1=hunit(:nom)    ; hunit2=hunit(2:)
call setq(q,hunit,no) ; call setq(q1,hunit1,nom)
rcrit=SQRT(EPSILON(1._vp))
!------------------------------------
! Initialize coefficients for quadrature, differencing and mdpt interpolation:
divq=967680        ; divd=1024
qco(0)=862564/divq ; dco(0)=1225/divd     ; ico(0)=1225/(2*divd)
qco(1)= 57249/divq ; dco(1)=-245/(3*divd) ; ico(1)=-245/(2*divd)
qco(2)= -5058/divq ; dco(2)=  49/(5*divd) ; ico(2)=  49/(2*divd)
qco(3)=   367/divq ; dco(3)=  -5/(7*divd) ; ico(3)=  -5/(2*divd)
qco(-1:-noh:-1)  = qco(1:noh) ! complete the stencil of quadrature coeffs.
dco(-1:-nohp:-1) =-dco(0:noh) ! complete the stencil of difference coeffs
ico(-1:-nohp:-1) = ico(0:noh) ! complete the stencil of interpolation coeffs.
!------------------------------------
! Initial coefficients related to control of working grid resolution:
ldsig =log(sigc/sigb)
ldsig4=ldsig**4
end subroutine infit
end module module_fitcons

!ENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDEND
!						**********************
!						* module_fitcons.f90 *
!						* PURSER 1994/1999   *
!                                               *   FUJITA 1999      *
!						**********************
!ENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDEND

!						********************
!						* coefrf.f90       *
!						* PURSER 1994/1999 *
!                                               *   FUJITA 1999    *
!						********************

!=============================================================================
subroutine coefrf(sig,nu,n,m,bnf,lnf)
!=============================================================================
! R. J. Purser NCEP 2001
!-----------------------------------------------------------------------------
use module_pmat2
implicit none
integer,              intent(IN   ) :: n,m
real(8), dimension(n),   intent(IN   ) :: sig,nu
real(8), dimension(n),   intent(OUT  ) :: bnf
real(8), dimension(m,n), intent(OUT  ) :: lnf
!-------------------------------------------------------------------------- 
integer, parameter                  :: irmax=6
real(8), dimension(n,-m:m)             :: s
real(8), dimension(n,-m:0)             :: sl
real(8), dimension(n,-m:m,m)           :: k,l
real(8), dimension(n)                  :: eta
real(8), dimension(irmax)              :: bcofi,bcofh
integer                             :: i,i1,il,ir,ik
!--------------------------------------------------------------------------
! The coefficients bcofi are the reciprocals of the i=1 entries of TABLE 1
! of NCEP O.N. 431:
data bcofi/1., 12., 90., 560., 3150., 16632./
!=============================================================================
bcofh=.5/bcofi
do i=1,n
   eta(i)=sig(i)*sqrt(nu(i))
enddo
k=0
!-------------------------------------------------------------------------
! Set k(:, -1:1, 1) to be the K-matrix of (4.8)--(4.10) of NCEP O.N. 431: 
!--------------------------------------------------------------------------
do i=1,n-1
   k(i  , 0,1)=k(i,0,1)  +eta(i+1)/eta(i)
   k(i+1, 0,1)=k(i+1,0,1)+eta(i)/eta(i+1)
   k(i  , 1,1)=-1
   k(i+1,-1,1)=-1
enddo

!-------------------------------------------------------------------------
! Set k(:, : , ir) to be the original K-matrix raised to the power of (ir):
!--------------------------------------------------------------------------
do ir=2,m
   il=ir-1
   call mulbb(k(:,-1:1,1),k(:,-il:il,il),k(:,-ir:ir,ir),n,n,1,1,il,il,ir,ir)
enddo

!-------------------------------------------------------------------------
! Pre- and post-multiply each of the m powers of K by the diagonal matrix,
! sigma, of NCEP O.N. 431, where the elements of sigma measure the smoothing
! scale of the quasi-Gaussian filter in grid-space units.
! Also, multiply each of the resulting banded matrices by .5*b_{1,ir} for
! the appropriate index, ir, corresponding to the power by which the original
! K was raised.
!--------------------------------------------------------------------------
do ir=1,m
   call mulbd(k(:,-ir:ir,ir),sig,k(:,-ir:ir,ir),n,n,ir,ir)
   call muldb(sig,k(:,-ir:ir,ir),k(:,-ir:ir,ir),n,n,ir,ir)
   k(:,-ir:ir,ir)=k(:,-ir:ir,ir)*bcofh(ir)
enddo


s=0
s(:,0)=1.

do ir=1,m
   l(:,-ir:ir,ir)=k(:,-ir:ir,ir)
   s(:,-ir:ir)=s(:,-ir:ir)+l(:,-ir:ir,ir)
enddo
do i1=2,m
   do ir=m,i1,-1
      l(:,-ir:ir,ir)=0.
      do ik=1,ir-i1+1
         il=ir-ik
         call madbb(k(:,-ik:ik,ik),l(:,-il:il,il),l(:,-ir:ir,ir), &
              n,n,ik,ik,il,il,ir,ir)
      enddo
      l(:,-ir:ir,ir)=l(:,-ir:ir,ir)/i1
      s(:,-ir:ir)=s(:,-ir:ir)+l(:,-ir:ir,ir)
   enddo
enddo
call ldlb(s,sl,bnf,n,m)
do i1=1,m
do i=1,n
   lnf(i1,i)=sl(i,-i1)
enddo
enddo
end subroutine coefrf


!ENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDEND
!						********************
!						* coefrf.f90       *
!						* PURSER 1994/1999 *
!                                               *   FUJITA 1999    *
!						********************
!ENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDEND

!						********************
!						* hgnrf.f90        *
!						* PURSER 1994/1999 *
!                                               *   FUJITA 1999    *
!						********************

!============================================================================
subroutine ldlb1i(nol,lnf,bnf,                                              &
       ids,ide,                                                             &
       ims,ime,                                                             &
       its,ite                                                              )
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide
  INTEGER, INTENT(IN   ) :: ims,ime
  INTEGER, INTENT(IN   ) :: its,ite

  REAL(8), DIMENSION(ims:ime),                       &
           INTENT(INOUT) :: bnf
  REAL(8), DIMENSION(nol, ims:ime),                  &
           INTENT(INOUT) :: lnf
!----------------------------------------------------------------------------
  INTEGER                :: i,l,m,nola
  real(8)                   :: s
!============================================================================
do i=its,ite
   nola=min(nol,i-its)
   do l=nola,1,-1
      s=lnf(l,i)
      do m=l+1,nola
         s=s-lnf(m,i)*bnf(i-m)*lnf(m-l,i-l)
      enddo
      lnf(l,i)=s/bnf(i-l)
   enddo
   s=bnf(i)
   do l=1,nola
      s=s-lnf(l,i)**2*bnf(i-l)
   enddo
   bnf(i)=s
enddo
end subroutine ldlb1i
   
!============================================================================
subroutine ldlb2i(nol,lnf,bnf,                                              &
       ids,ide, jds,jde,                                                    &
       ims,ime, jms,jme,                                                    &
       its,ite, jts,jte                                                     )
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide, jds,jde
  INTEGER, INTENT(IN   ) :: ims,ime, jms,jme
  INTEGER, INTENT(IN   ) :: its,ite, jts,jte

  REAL(8), DIMENSION(ims:ime, jms:jme),                       &
           INTENT(INOUT) :: bnf
  REAL(8), DIMENSION(nol, ims:ime, jms:jme),                  &
           INTENT(INOUT) :: lnf
!----------------------------------------------------------------------------
  INTEGER                :: i,j,l,m,nola
  real(8)                   :: s
!============================================================================
do j=jts,jte
do i=its,ite
   nola=min(nol,i-its)
   do l=nola,1,-1
      s=lnf(l,i,j)
      do m=l+1,nola
         s=s-lnf(m,i,j)*bnf(i-m,j)*lnf(m-l,i-l,j)
      enddo
      lnf(l,i,j)=s/bnf(i-l,j)
   enddo
   s=bnf(i,j)
   do l=1,nola
      s=s-lnf(l,i,j)**2*bnf(i-l,j)
   enddo
   bnf(i,j)=s
enddo
enddo
end subroutine ldlb2i
   
!============================================================================
subroutine ldlb2j(nol,lnf,bnf,                                              &
       ids,ide, jds,jde,                                                    &
       ims,ime, jms,jme,                                                    &
       its,ite, jts,jte                                                     )
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide, jds,jde
  INTEGER, INTENT(IN   ) :: ims,ime, jms,jme
  INTEGER, INTENT(IN   ) :: its,ite, jts,jte

  REAL(8), DIMENSION(ims:ime, jms:jme),                       &
           INTENT(INOUT) :: bnf
  REAL(8), DIMENSION(nol, ims:ime, jms:jme),                  &
           INTENT(INOUT) :: lnf
!----------------------------------------------------------------------------
  INTEGER                :: i,j,l,m,nola
  real(8)                   :: s
!============================================================================
do j=jts,jte
   nola=min(nol,j-jts)
   do i=its,ite
   do l=nola,1,-1
      s=lnf(l,i,j)
      do m=l+1,nola
         s=s-lnf(m,i,j)*bnf(i,j-m)*lnf(m-l,i,j-l)
      enddo
      lnf(l,i,j)=s/bnf(i,j-l)
   enddo
   s=bnf(i,j)
   do l=1,nola
      s=s-lnf(l,i,j)**2*bnf(i,j-l)
   enddo
   bnf(i,j)=s
   enddo
enddo
end subroutine ldlb2j
   
!============================================================================
subroutine ldlb3i(nol,lnf,bnf,                                              &
       ids,ide, jds,jde, kds,kde,                                           &
       ims,ime, jms,jme, kms,kme,                                           &
       its,ite, jts,jte, kts,kte                                            )
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide, jds,jde, kds,kde 
  INTEGER, INTENT(IN   ) :: ims,ime, jms,jme, kms,kme
  INTEGER, INTENT(IN   ) :: its,ite, jts,jte, kts,kte

  REAL(8), DIMENSION(ims:ime, kms:kme, jms:jme),                       &
           INTENT(INOUT) :: bnf
  REAL(8), DIMENSION(nol, ims:ime, kms:kme, jms:jme),                  &
           INTENT(INOUT) :: lnf
!----------------------------------------------------------------------------
  INTEGER                :: i,j,k,l,m,nola
  real(8)                   :: s
!============================================================================
do j=jts,jte
do k=kts,kte
do i=its,ite
   nola=min(nol,i-its)
   do l=nola,1,-1
      s=lnf(l,i,k,j)
      do m=l+1,nola
         s=s-lnf(m,i,k,j)*bnf(i-m,k,j)*lnf(m-l,i-l,k,j)
      enddo
      lnf(l,i,k,j)=s/bnf(i-l,k,j)
   enddo
   s=bnf(i,k,j)
   do l=1,nola
      s=s-lnf(l,i,k,j)**2*bnf(i-l,k,j)
   enddo
   bnf(i,k,j)=s
enddo
enddo
enddo
end subroutine ldlb3i
   
!============================================================================
subroutine ldlb3j(nol,lnf,bnf,                                              &
       ids,ide, jds,jde, kds,kde,                                           &
       ims,ime, jms,jme, kms,kme,                                           &
       its,ite, jts,jte, kts,kte                                            )
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide, jds,jde, kds,kde 
  INTEGER, INTENT(IN   ) :: ims,ime, jms,jme, kms,kme
  INTEGER, INTENT(IN   ) :: its,ite, jts,jte, kts,kte

  REAL(8), DIMENSION(ims:ime, kms:kme, jms:jme),                       &
           INTENT(INOUT) :: bnf
  REAL(8), DIMENSION(nol, ims:ime, kms:kme, jms:jme),                  &
           INTENT(INOUT) :: lnf
!----------------------------------------------------------------------------
  INTEGER                :: i,j,k,l,m,nola
  real(8)                   :: s
!============================================================================
do j=jts,jte
   nola=min(nol,j-jts)
   do k=kts,kte
   do i=its,ite
   do l=nola,1,-1
      s=lnf(l,i,k,j)
      do m=l+1,nola
         s=s-lnf(m,i,k,j)*bnf(i,k,j-m)*lnf(m-l,i,k,j-l)
      enddo
      lnf(l,i,k,j)=s/bnf(i,k,j-l)
   enddo
   s=bnf(i,k,j)
   do l=1,nola
      s=s-lnf(l,i,k,j)**2*bnf(i,k,j-l)
   enddo
   bnf(i,k,j)=s
   enddo
   enddo
enddo
end subroutine ldlb3j
   
SUBROUTINE hbnrf1i(a,nol,lnf,bnf,                                           &
       ids,ide,                                                             &
       ims,ime,                                                             &
       its,ite                                                              )
!============================================================================
! Horizontal basic inhomogeneous recursive filter, 
! 1-dimensional, active index i
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide
  INTEGER, INTENT(IN   ) :: ims,ime
  INTEGER, INTENT(IN   ) :: its,ite

  REAL(8), DIMENSION(ims:ime),                       &
           INTENT(INOUT) :: a
  REAL(8), DIMENSION(ims:ime),                       &
           INTENT(IN   ) :: bnf
  REAL(8), DIMENSION(nol, ims:ime),                  &
           INTENT(IN   ) :: lnf
!----------------------------------------------------------------------------
  INTEGER                :: i,l,nola
!============================================================================
DO i=its+1,ite
   nola=MIN(nol,i-its)
   DO l=1,nola
      a(i)=a(i)-lnf(l,i)*a(i-l)
   ENDDO
ENDDO
DO i=its,ite
   a(i)=bnf(i)*a(i)
ENDDO
DO i=ite-1,its,-1
   nola=MIN(nol,ite-i)
   DO l=1,nola
      a(i)=a(i)-lnf(l,i+l)*a(i+l)
   ENDDO
ENDDO
END SUBROUTINE hbnrf1i

SUBROUTINE hbnrf2i(a,nol,lnf,bnf,                                           &
       ids,ide, jds,jde,                                                    &
       ims,ime, jms,jme,                                                    &
       its,ite, jts,jte                                                     )
!============================================================================
! Horizontal basic inhomogeneous recursive filter, 
! 2-dimensional, active index i
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide, jds,jde
  INTEGER, INTENT(IN   ) :: ims,ime, jms,jme
  INTEGER, INTENT(IN   ) :: its,ite, jts,jte

  REAL(8), DIMENSION(ims:ime, jms:jme),                       &
           INTENT(INOUT) :: a
  REAL(8), DIMENSION(ims:ime, jms:jme),                       &
           INTENT(IN   ) :: bnf
  REAL(8), DIMENSION(nol, ims:ime, jms:jme),                  &
           INTENT(IN   ) :: lnf
!----------------------------------------------------------------------------
  INTEGER                :: i,j,l,nola
!============================================================================
DO j=jts,jte
   DO i=its+1,ite
      nola=MIN(nol,i-its)
      DO l=1,nola
         a(i,j)=a(i,j)-lnf(l,i,j)*a(i-l,j)
      ENDDO
   ENDDO
   DO i=its,ite
      a(i,j)=bnf(i,j)*a(i,j)
   ENDDO
   DO i=ite-1,its,-1
      nola=MIN(nol,ite-i)
      DO l=1,nol
         a(i,j)=a(i,j)-lnf(l,i+l,j)*a(i+l,j)
      ENDDO
   ENDDO
ENDDO
END SUBROUTINE hbnrf2i

SUBROUTINE hbnrf2j(a,nol,lnf,bnf,                                           &
       ids,ide, jds,jde,                                                    &
       ims,ime, jms,jme,                                                    &
       its,ite, jts,jte                                                     )
!============================================================================
! Horizontal basic inhomogeneous recursive filter, 
! 2-dimensional, active index j
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide, jds,jde
  INTEGER, INTENT(IN   ) :: ims,ime, jms,jme
  INTEGER, INTENT(IN   ) :: its,ite, jts,jte

  REAL(8), DIMENSION(ims:ime, jms:jme),                       &
           INTENT(INOUT) :: a
  REAL(8), DIMENSION(ims:ime, jms:jme),                       &
           INTENT(IN   ) :: bnf
  REAL(8), DIMENSION(nol, ims:ime, jms:jme),                  &
           INTENT(IN   ) :: lnf
!----------------------------------------------------------------------------
  INTEGER                :: i,j,l,nola
!============================================================================
DO j=jts+1,jte
   nola=MIN(nol,j-jts)
   DO i=its,ite
      DO l=1,nola
         a(i,j)=a(i,j)-lnf(l,i,j)*a(i,j-l)
      ENDDO
   ENDDO
ENDDO
DO j=jts,jte
   DO i=its,ite
      a(i,j)=bnf(i,j)*a(i,j)
   ENDDO
ENDDO
DO j=jte-1,jts,-1
   nola=MIN(nol,jte-j)
   DO i=its,ite
      DO l=1,nola
         a(i,j)=a(i,j)-lnf(l,i,j+l)*a(i,j+l)
      ENDDO
   ENDDO
ENDDO
END SUBROUTINE hbnrf2j

SUBROUTINE hbnrf3i(a,nol,lnf,bnf,                                           &
       ids,ide, jds,jde, kds,kde,                                           &
       ims,ime, jms,jme, kms,kme,                                           &
       its,ite, jts,jte, kts,kte                                            )
!============================================================================
! Horizontal basic inhomogeneous recursive filter, 
! 3-dimensional, active index i
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide, jds,jde, kds,kde 
  INTEGER, INTENT(IN   ) :: ims,ime, jms,jme, kms,kme
  INTEGER, INTENT(IN   ) :: its,ite, jts,jte, kts,kte

  REAL(8), DIMENSION(ims:ime, kms:kme, jms:jme),                       &
           INTENT(INOUT) :: a
  REAL(8), DIMENSION(ims:ime, kms:kme, jms:jme),                       &
           INTENT(IN   ) :: bnf
  REAL(8), DIMENSION(nol, ims:ime, kms:kme, jms:jme),                  &
           INTENT(IN   ) :: lnf
!----------------------------------------------------------------------------
  INTEGER                :: i,j,k,l,nola
!============================================================================
DO j=jts,jte
   DO k=kts,kte
      DO i=its+1,ite
         nola=MIN(nol,i-its)
         DO l=1,nola
            a(i,k,j)=a(i,k,j)-lnf(l,i,k,j)*a(i-l,k,j)
         ENDDO
      ENDDO
      DO i=its,ite
         a(i,k,j)=bnf(i,k,j)*a(i,k,j)
      ENDDO
      DO i=ite-1,its,-1
         nola=MIN(nol,ite-i)
         DO l=1,nola
            a(i,k,j)=a(i,k,j)-lnf(l,i+l,k,j)*a(i+l,k,j)
         ENDDO
      ENDDO
   ENDDO
ENDDO
END SUBROUTINE hbnrf3i

SUBROUTINE hbnrf3j(a,nol,lnf,bnf,                                           &
       ids,ide, jds,jde, kds,kde,                                           &
       ims,ime, jms,jme, kms,kme,                                           &
       its,ite, jts,jte, kts,kte                                            )
!============================================================================
! Horizontal basic inhomogeneous recursive filter, 
! 3-dimensional, active index j
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide, jds,jde, kds,kde 
  INTEGER, INTENT(IN   ) :: ims,ime, jms,jme, kms,kme
  INTEGER, INTENT(IN   ) :: its,ite, jts,jte, kts,kte

  REAL(8), DIMENSION(ims:ime, kms:kme, jms:jme),                       &
           INTENT(INOUT) :: a
  REAL(8), DIMENSION(ims:ime, kms:kme, jms:jme),                       &
           INTENT(IN   ) :: bnf
  REAL(8), DIMENSION(nol, ims:ime, kms:kme, jms:jme),                  &
           INTENT(IN   ) :: lnf
!----------------------------------------------------------------------------
  INTEGER                :: i,j,k,l,nola
!============================================================================
DO j=jts+1,jte
   nola=MIN(nol,j-jts)
   DO k=kts,kte
      DO i=its,ite
         DO l=1,nola
            a(i,k,j)=a(i,k,j)-lnf(l,i,k,j)*a(i,k,j-l)
         ENDDO
      ENDDO
   ENDDO
ENDDO
DO j=jts,jte
   DO k=kts,kte
      DO i=its,ite
         a(i,k,j)=bnf(i,k,j)*a(i,k,j)
      ENDDO
   ENDDO
ENDDO
DO j=jte-1,jts,-1
   nola=MIN(nol,jte-j)
   DO k=kts,kte
      DO i=its,ite
         DO l=1,nola
            a(i,k,j)=a(i,k,j)-lnf(l,i,k,j+l)*a(i,k,j+l)
         ENDDO
      ENDDO
   ENDDO
ENDDO
END SUBROUTINE hbnrf3j


SUBROUTINE vbnrf1k(a,nol,lnf,bnf,                                           &
       kds,kde,                                                             &
       kms,kme,                                                             &
       kts,kte                                                              )
!============================================================================
! Vertical bounded grid inhomogeneous recursive filter, 
! 1-dimensional, active index k
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: kds,kde
  INTEGER, INTENT(IN   ) :: kms,kme
  INTEGER, INTENT(IN   ) :: kts,kte

  REAL(8), DIMENSION(kms:kme),                       &
           INTENT(INOUT) :: a
  REAL(8), DIMENSION(kms:kme),                       &
           INTENT(IN   ) :: bnf
  REAL(8), DIMENSION(nol, kms:kme),                  &
           INTENT(IN   ) :: lnf
!----------------------------------------------------------------------------
  INTEGER                :: k,l,nola
!============================================================================
DO k=kts+1,kte
   nola=MIN(nol,k-kts)
   DO l=1,nola
      a(k)=a(k)-lnf(l,k)*a(k-l)
   ENDDO
ENDDO
DO k=kts,kte
   a(k)=bnf(k)*a(k)
ENDDO
DO k=kte-1,kts,-1
   nola=MIN(nol,kte-k)
   DO l=1,nola
      a(k)=a(k)-lnf(l,k+l)*a(k+l)
   ENDDO
ENDDO
END SUBROUTINE vbnrf1k

SUBROUTINE vbnrf2k(a,nol,lnf,bnf,                                           &
       ids,ide, kds,kde,                                                    &
       ims,ime, kms,kme,                                                    &
       its,ite, kts,kte                                                     )
!============================================================================
! Vertical bounded grid inhomogeneous recursive filter, 
! 2-dimensional, active index k
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide, kds,kde 
  INTEGER, INTENT(IN   ) :: ims,ime, kms,kme
  INTEGER, INTENT(IN   ) :: its,ite, kts,kte

  REAL(8), DIMENSION(ims:ime, kms:kme),                       &
           INTENT(INOUT) :: a
  REAL(8), DIMENSION(ims:ime, kms:kme),                       &
           INTENT(IN   ) :: bnf
  REAL(8), DIMENSION(nol, ims:ime, kms:kme),                  &
           INTENT(IN   ) :: lnf
!----------------------------------------------------------------------------
  INTEGER                :: i,k,l,nola
!============================================================================
DO k=kts+1,kte
   nola=MIN(nol,k-kts)
   DO i=its,ite
      DO l=1,nola
         a(i,k)=a(i,k)-lnf(l,i,k)*a(i,k-l)
      ENDDO
   ENDDO
ENDDO
DO k=kts,kte
   DO i=its,ite
      a(i,k)=bnf(i,k)*a(i,k)
   ENDDO
ENDDO
DO k=kte-1,kts,-1
   nola=MIN(nol,kte-k)
   DO i=its,ite
      DO l=1,nola
         a(i,k)=a(i,k)-lnf(l,i,k+l)*a(i,k+l)
      ENDDO
   ENDDO
ENDDO
END SUBROUTINE vbnrf2k

SUBROUTINE vbnrf3k(a,nol,lnf,bnf,                                           &
       ids,ide, jds,jde, kds,kde,                                           &
       ims,ime, jms,jme, kms,kme,                                           &
       its,ite, jts,jte, kts,kte                                            )
!============================================================================
! Vertical bounded grid inhomogeneous recursive filter, 
! 3-dimensional, active index k
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide, jds,jde, kds,kde 
  INTEGER, INTENT(IN   ) :: ims,ime, jms,jme, kms,kme
  INTEGER, INTENT(IN   ) :: its,ite, jts,jte, kts,kte

  REAL(8), DIMENSION(ims:ime, kms:kme, jms:jme),                       &
           INTENT(INOUT) :: a
  REAL(8), DIMENSION(ims:ime, kms:kme, jms:jme),                       &
           INTENT(IN   ) :: bnf
  REAL(8), DIMENSION(nol, ims:ime, kms:kme, jms:jme),                  &
           INTENT(IN   ) :: lnf
!----------------------------------------------------------------------------
  INTEGER                :: i,j,k,l,nola
!============================================================================
DO j=jts,jte
   DO k=kts+1,kte
      nola=MIN(nol,k-kts)
      DO i=its,ite
         DO l=1,nola
            a(i,k,j)=a(i,k,j)-lnf(l,i,k,j)*a(i,k-l,j)
         ENDDO
      ENDDO
   ENDDO
   DO k=kts,kte
      DO i=its,ite
         a(i,k,j)=bnf(i,k,j)*a(i,k,j)
      ENDDO
   ENDDO
   DO k=kte-1,kts,-1
      nola=MIN(nol,kte-k)
      DO i=its,ite
         DO l=1,nola
            a(i,k,j)=a(i,k,j)-lnf(l,i,k+l,j)*a(i,k+l,j)
         ENDDO
      ENDDO
   ENDDO
ENDDO
END SUBROUTINE vbnrf3k

SUBROUTINE hbncij(a,hamp,nol,lnfi,bnfi,lnfj,bnfj,                           &
     ids,ide, jds,jde,                                                      &
     ims,ime, jms,jme,                                                      &
     its,ite, jts,jte                                                       )
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide, jds,jde
  INTEGER, INTENT(IN   ) :: ims,ime, jms,jme
  INTEGER, INTENT(IN   ) :: its,ite, jts,jte

  REAL(8), DIMENSION(ims:ime, jms:jme),                       &
           INTENT(INOUT) :: a
  REAL(8), DIMENSION(ims:ime, jms:jme),                       &
           INTENT(IN   ) :: hamp,bnfi,bnfj
  REAL(8), DIMENSION(nol, ims:ime, jms:jme),                  &
           INTENT(IN   ) :: lnfi,lnfj
!----------------------------------------------------------------------------
  INTEGER                :: i,j
!============================================================================
DO j=jts,jte
   DO i=its,ite
      a(i,j)=hamp(i,j)*a(i,j)
   ENDDO
ENDDO
!---------------
CALL hbnrf2i(a,nol,lnfi,bnfi,             &
     ids,ide, jds,jde,                    &
     ims,ime, jms,jme,                    &
     its,ite, jts,jte)
!----------
CALL hbnrf2j(a,nol,lnfj,bnfj,             &
     ids,ide, jds,jde,                    &
     ims,ime, jms,jme,                    &
     its,ite, jts,jte)
!----------
END SUBROUTINE hbncij

SUBROUTINE hbncji(a,hamp,nol,lnfi,bnfi,lnfj,bnfj,                           &
     ids,ide, jds,jde,                                                      &
     ims,ime, jms,jme,                                                      &
     its,ite, jts,jte                                                       )
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide, jds,jde
  INTEGER, INTENT(IN   ) :: ims,ime, jms,jme
  INTEGER, INTENT(IN   ) :: its,ite, jts,jte

  REAL(8), DIMENSION(ims:ime, jms:jme),                       &
           INTENT(INOUT) :: a
  REAL(8), DIMENSION(ims:ime, jms:jme),                       &
           INTENT(IN   ) :: hamp,bnfi,bnfj
  REAL(8), DIMENSION(nol, ims:ime, jms:jme),                  &
           INTENT(IN   ) :: lnfi,lnfj
!----------------------------------------------------------------------------
  INTEGER                :: i,j
!============================================================================
CALL hbnrf2j(a,nol,lnfj,bnfj,             &
     ids,ide, jds,jde,                    &
     ims,ime, jms,jme,                    &
     its,ite, jts,jte)
!----------
CALL hbnrf2i(a,nol,lnfi,bnfi,             &
     ids,ide, jds,jde,                    &
     ims,ime, jms,jme,                    &
     its,ite, jts,jte)
!---------------
DO j=jts,jte
   DO i=its,ite
      a(i,j)=hamp(i,j)*a(i,j)
   ENDDO
ENDDO
!---------------
END SUBROUTINE hbncji

SUBROUTINE hbncijk(a,hamp,nol,lnfi,bnfi,lnfj,bnfj,lnfk,bnfk,                &
       ids,ide, jds,jde, kds,kde,                                           &
       ims,ime, jms,jme, kms,kme,                                           &
       its,ite, jts,jte, kts,kte                                            )
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide, jds,jde, kds,kde 
  INTEGER, INTENT(IN   ) :: ims,ime, jms,jme, kms,kme
  INTEGER, INTENT(IN   ) :: its,ite, jts,jte, kts,kte

  REAL(8), DIMENSION(ims:ime, kms:kme, jms:jme),                       &
           INTENT(INOUT) :: a
  REAL(8), DIMENSION(ims:ime, kms:kme, jms:jme),                       &
           INTENT(IN   ) :: hamp,bnfi,bnfj,bnfk
  REAL(8), DIMENSION(nol, ims:ime, kms:kme, jms:jme),                  &
           INTENT(IN   ) :: lnfi,lnfj,lnfk
!----------------------------------------------------------------------------
  INTEGER                :: i,j,k
!============================================================================
DO j=jts,jte
   do k=kts,kte
      DO i=its,ite
         a(i,k,j)=hamp(i,k,j)*a(i,k,j)
      ENDDO
   enddo
ENDDO
!---------------
CALL hbnrf3i(a,nol,lnfi,bnfi,             &
       ids,ide, jds,jde, kds,kde,         &
       ims,ime, jms,jme, kms,kme,         &
       its,ite, jts,jte, kts,kte)
!----------
CALL hbnrf3j(a,nol,lnfj,bnfj,             &
       ids,ide, jds,jde, kds,kde,         &
       ims,ime, jms,jme, kms,kme,         &
       its,ite, jts,jte, kts,kte)
!----------
call vbnrf3k(a,nol,lnfk,bnfk,             &
       ids,ide, jds,jde, kds,kde,         &
       ims,ime, jms,jme, kms,kme,         &
       its,ite, jts,jte, kts,kte)
END SUBROUTINE hbncijk

SUBROUTINE hbnckji(a,hamp,nol,lnfi,bnfi,lnfj,bnfj,lnfk,bnfk,                &
       ids,ide, jds,jde, kds,kde,                                           &
       ims,ime, jms,jme, kms,kme,                                           &
       its,ite, jts,jte, kts,kte                                            )
!============================================================================
  IMPLICIT NONE

  INTEGER, INTENT(IN   ) :: nol
  INTEGER, INTENT(IN   ) :: ids,ide, jds,jde, kds,kde 
  INTEGER, INTENT(IN   ) :: ims,ime, jms,jme, kms,kme
  INTEGER, INTENT(IN   ) :: its,ite, jts,jte, kts,kte

  REAL(8), DIMENSION(ims:ime, kms:kme, jms:jme),                       &
           INTENT(INOUT) :: a
  REAL(8), DIMENSION(ims:ime, kms:kme, jms:jme),                       &
           INTENT(IN   ) :: hamp,bnfi,bnfj,bnfk
  REAL(8), DIMENSION(nol, ims:ime, kms:kme, jms:jme),                  &
           INTENT(IN   ) :: lnfi,lnfj,lnfk
!----------------------------------------------------------------------------
  INTEGER                :: i,j,k
!============================================================================
call vbnrf3k(a,nol,lnfk,bnfk,             &
       ids,ide, jds,jde, kds,kde,         &
       ims,ime, jms,jme, kms,kme,         &
       its,ite, jts,jte, kts,kte)
!----------
CALL hbnrf3j(a,nol,lnfj,bnfj,             &
       ids,ide, jds,jde, kds,kde,         &
       ims,ime, jms,jme, kms,kme,         &
       its,ite, jts,jte, kts,kte)
!----------
CALL hbnrf3i(a,nol,lnfi,bnfi,             &
       ids,ide, jds,jde, kds,kde,         &
       ims,ime, jms,jme, kms,kme,         &
       its,ite, jts,jte, kts,kte)
!---------------
DO j=jts,jte
   do k=kts,kte
      DO i=its,ite
         a(i,k,j)=hamp(i,k,j)*a(i,k,j)
      ENDDO
   enddo
ENDDO
!---------------
END SUBROUTINE hbnckji


!ENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDEND
!						********************
!						* hgnrf.f90        *
!						* PURSER 1994/1999 *
!                                               *   FUJITA 1999    *
!						********************
!ENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDEND

!						********************
!						* rfit.f90         *
!						* PURSER 1994/1999 *
!                                               *   FUJITA 1999    *
!						********************

!============================================================================
subroutine rfit(ng,sig,nu, ns,nw,ssig,snu,ins1,wts)  
!============================================================================
! R. J. Purser, NCEP 2001
!----------------------------------------------------------------------------
use vkind
use module_fitcons
implicit none
integer,                      intent(IN   ):: ng
real(vp),dimension(ng),       intent(IN   ):: sig,nu
integer,                      intent(OUT  ):: ns,nw
real(vp),dimension(ng),       intent(OUT  ):: ssig,snu
integer, dimension(ng),       intent(INOUT):: ins1
real(vp),dimension(no,ng),    intent(INOUT):: wts
!----------------------------------------------------------------------------
integer                                    :: i,i1,im,k,l,is,is0,is1,isn
real(vp)                                   :: t
real(vp),dimension(-nohm:ng+noh)           :: dcdg
real(vp),dimension(-noh:ng+noh)            :: cofg,cofs
real(vp),dimension(ng)                     :: dsdg,dhdg,rsnu
real(vp)                                   :: rnu
!============================================================================
nw=0
do i=1,ng
   dcdg(i)=1./sig(i)
   if(sig(i) <= sigb)then
!----------------------------------------------------------------------------
! sig(i) below threshold; cleave to original grid spacing with ds/dg and 
! dh/dg set accordingly:
!----------------------------------------------------------------------------
      dsdg(i)=1.     ;      dhdg(i)=0.
   else
!----------------------------------------------------------------------------
! sig(i) exceeds basic threshold sigb, allowing working grid with coordinate
! s to differ from original grid with coordinate g. The formula for ds/dg
! is now <1 but tends smoothly to 1 again at the threshold value, sig=sigb.
! [The function for log(ds/dg) is based on the "hyper-hyperbola":
!    y= (1+x**4)**(-4)-1, which rises very gradually from its base at x=y=0]
! Likewise, the perturbative component, dh/dg, is now < 0, but tends
! smoothly to 0 again at the threshold.
!----------------------------------------------------------------------------
      t=ldsig-sqrt(sqrt(ldsig4+(ldsig-log(sigc*dcdg(i)))**4))
      dsdg(i)=exp(t) ;      dhdg(i)=dsdg(i)*t
   endif
enddo

!----------------------------------------------------------------------------
! Apply mirror-symmetry to extrapolate beyond ends:
!----------------------------------------------------------------------------
do l=1,noh
   dcdg(1-l)=dcdg(l); dcdg(ng+l)=dcdg(ng+1-l)
enddo

!----------------------------------------------------------------------------
! Integrate dc/dg wrt g to get c(g) at each of the points of the g-grid
! which is NOT staggered relative to the boundary
!----------------------------------------------------------------------------
cofg(0)=0.
do i=1,ng
   cofg(i)=cofg(i-1)+dot_product(qco,dcdg(i-noh:i+noh))
enddo
do l=1,noh
   cofg(  -l)=-cofg(   l)
   cofg(ng+l)=-cofg(ng-l)+2*cofg(ng)
enddo

im=0
ns=0
!----------------------------------------------------------------------------
! loop over noncontiguous segments where it is numerically beneficial
! to employ a grid of relatively coarse resolution. The adoption of each
! alternative grid segment is subject to some conditions:
! 1) Each coarse-grid segment must span at least 5 points of the original grid
! 2) Each segment must shorten the tally of working grid points by at least 3.
!     Subject to the above conditions, the coarse grid is blended smoothly
! with the original grid at internal thresholds and is designed to provide
! a resolution such that the smoothing scale, sigma, never exceeds the 
! product, sigc*dg/ds, where sigc is a dimensionless parameter (e.g. sigc=3.)
! and dg/ds is the local working grid (s) spacing in units of the original
! grid (g) spacing. 
!
! Each segment found is defined by its end points in the original grid,
! i1 and im. k is the counter for segments along this line.
! ns keeps count of the number of working grid (s-grid) points found so far.
!----------------------------------------------------------------------------
cofs(0)=0.
do k=1,ng 
   do i1=im+1,ng
      if(i1< ng-3 .and. dhdg(i1) /= 0)exit
!----------------------------------------------------------------------------
! working s-grid continues to track the original g-grid; Set indices and 
! weight for the trivial "interpolation" between these coincident grids:
!----------------------------------------------------------------------------
      ns=ns+1
      ins1(i1)=-ns
      cofs(ns)=cofg(i1)
   enddo
   if(i1 > ng)exit
!----------------------------------------------------------------------------
! Having met the basic conditions for the start of a new segment in which
! the s-grid and g-grids may part company, seek the other end, im, of this
! possible segment:
!----------------------------------------------------------------------------
   do im=i1+1,ng
      if(dhdg(im) == 0)exit
   enddo
   im=im-1
   if(im < i1+4)then
!----------------------------------------------------------------------------
! Segment too short to be viable; keep s-grid and g-grids synchronized:
!----------------------------------------------------------------------------
      do i=i1,im
         ns=ns+1
         ins1(i)=-ns
         cofs(ns)=cofg(i)
      enddo
   else
!----------------------------------------------------------------------------
! Segment long enough to be potentially viable. Call jfit to determine if 
! the final condition is met, namely that the number of s-grid points 
! in this segment is smaller than the g-grid tally by at least 3. If so,
! Fit an exact integer number of s-points into this segment and compute
! the indices and weights for the associated nontrivial interpolation
! from these (and neighboring) s-points back to the g-points of the segment:
!----------------------------------------------------------------------------
      call jfit(ng,i1,im,ns,nw,cofg,dsdg,dhdg,cofs,ins1,wts)
      if(ns.lt.0) return
   endif
enddo
if(ns<no .and. nw>0)then ! <- s-grid too short; use copy of g-grid instead
   wts(:,1:nw)=0
   nw=0
   do i=1,ng
      ins1(i)=-i
      cofs(i)=cofg(i)
   enddo
   ns=ng
endif

do l=1,noh
   cofs(  -l)=-cofs(   l)
   cofs(ns+l)=-cofs(ns-l)+2*cofs(ns)
enddo
do is=1,ns
   ssig(is)=1./dot_product(dco,cofs(is-nohp:is+noh))
enddo

!----------------------------------------------------------------------------
! By applying adjoint-interpolation to the g-grid metric terms, obtain
! the corresponding metric terms for the new s-grid:
!----------------------------------------------------------------------------
call stogt(ns,ng,ins1,wts, snu,nu)

end subroutine rfit

!============================================================================
subroutine jfit(ng,ig1,igm,ns,iw,cofg,dsdg,dhdg,cofs,ins1,wts)
!============================================================================
! R. J. Purser, NCEP 2001
!----------------------------------------------------------------------------
use vkind
use module_fitcons
implicit none
integer,                         intent(IN   ):: ng,ig1,igm
integer,                         intent(INOUT):: ns,iw
real(vp),dimension(ng),          intent(IN   ):: dsdg,dhdg
real(vp),dimension(-noh:ng+noh), intent(IN   ):: cofg
real(vp),dimension(-noh:ng+noh), intent(INOUT):: cofs
integer, dimension(ng),          intent(INOUT):: ins1
real(vp),dimension(no,ng),       intent(INOUT):: wts
!----------------------------------------------------------------------------
real(vp),dimension(-noh:ng+noh) :: sofg,dsdgt
real(vp)                        :: et,estar,destar,r,dr,sm,hm
integer                         :: i,l,ie,iep,ie0,ie1,ien,ig0,is0,ism,init
!============================================================================

!----------------------------------------------------------------------------
! Form the definite integral sm, of ds/dg, within this segment:
!----------------------------------------------------------------------------
sm=sum(dsdg(ig1:igm)) 

!---------------------------------------------------------------------------
! Test whether it is worthwhile to allow s-grid to deviate from the original
! g-grid within this segment on the basis of the number of grid points that
! could potentially be eliminated (we require a saving > 3 per segment):
!---------------------------------------------------------------------------
if(sm > igm-ig1-2)then
!----------------------------------------------------------------------------
! This putative s-grid segment reduces the total number of grid points by an
! insufficient amount to justify the additional interpolations. Therefore,
! keep the original g-grid instead for this segment, and return:
!---------------------------------------------------------------------------
   do i=ig1,igm
      ns=ns+1
      ins1(i)=-ns
      cofs(ns)=cofg(i)
   enddo
   return
endif
!----------------------------------------------------------------------------
! s-grid segment achieves a worthwhile reduction of the number of points
! of the working grid. The tasks of the rest of this routine are to:
! (1) adjust the segment length in the s-metric to make it an integer;
! (2) find the s-coordinates of each g-grid points in this segment
!     and hence the nontrivial interpolation indices and weights required 
!     to go from the s-grid to the g-grid (or adjoints going the other way);
! (3) use Newton iterations to find the accurate interpolation formulae
!     that enable c(s) to be interpolated from the given c(g).
!----------------------------------------------------------------------------
ig0=ig1-1
is0=ns; ism=sm
!----------------------------------------------------------------------------
! Fractional remainder of sm, divided by the definite integral of dh/dg
! provides the adjustment factor that scales the perturbative component,
! dhdg, by exactly the amount that will make the segment integral of the 
! perturbed grid-to-grid jacobian, dsdgt, the exact integer, ism:
!----------------------------------------------------------------------------
r=(sm-ism)/sum(dhdg(ig1:igm))
do i=ig1,igm
   dsdgt(i)=dsdg(i)-r*dhdg(i)
enddo
!----------------------------------------------------------------------------
! Mirror-extrapolate adjusted ds/dg as an even-symmetry function at the 
! ends of this segment. Note that the grid on which derivatives such as
! ds/dg reside is the one staggered wrt domain boundaries and segment
! end points. The indices of this grid go from ig1 to igm inside the
! segment. (The convention for the companion grid, NOT staggered wrt 
! boundaries, is such that the two segment ends are denoted by indices,
! ig0=ig1-1 and igm.)
!----------------------------------------------------------------------------
do l=1,noh
   dsdgt(ig1-l)=dsdgt(ig0  +l)
   dsdgt(igm+l)=dsdgt(igm+1-l)
enddo
ism=is0+ism ! This integer also becomes (within round-off) the value, sofg(igm)
!----------------------------------------------------------------------------
! Set s(g) at both ends of the segment to be the appropriate integers:
!----------------------------------------------------------------------------
sofg(ig0)=is0; sofg(igm)=ism
!----------------------------------------------------------------------------
! Get s(g) inside the segment by performing a numerical quadrature of dsdgt:
!----------------------------------------------------------------------------
do i=ig1,igm
   sofg(i)=sofg(i-1)+dot_product(qco,dsdgt(i-noh:i+noh))
enddo
!----------------------------------------------------------------------------
! Mirror-extrapolate s(g) as an odd-symmetry function at segment end points.
! Note that, being an inegral, s(g) resides on the grid NOT staggered wrt
! boundaries and segment end points.
!----------------------------------------------------------------------------
do l=1,noh
   sofg(ig0-l)=2*is0-sofg(ig0+l)
   sofg(igm+l)=2*ism-sofg(igm-l)
enddo
do i=ig1,igm
   iw=iw+1 ; wts(:,iw)=0
   r=dot_product(ico,sofg(i-nohp:i+noh))+.5_vp
   ie=r            ! Take integer part...
   ins1(i)=ie-nohm ! ...hence the index of the first point in the stencil...
   r=r-ie          ! ...then the fractional part to find interpolation weights:
   call lagw(hunit1,r,q1,wt1,dwt1,nom)   ! weights for left-biased stencil
   wts(:nom,iw) =              (1-r)*wt1 !   bias weight, 1-r
   call lagw(hunit2,r,q1,wt1,dwt1,nom)   ! weights for right-biased stencil
   wts(2:   ,iw) = wts(2:   ,iw)  +r*wt1 !   bias weight, r.
!----------------------------------------------------------------------------
! Exploit the mirror symmetries to confine the weight stencil to the 
! domain interior, even though this may entail padding innermost end of
! the stencil with useless zeroes:
!----------------------------------------------------------------------------
   L=1-INS1(I)
   IF(L > 0)THEN ! FOLD LEFT OVERLAP OF L ELEMENTS BACK INSIDE: 
      WTS(1:L,IW)      =WTS(L:1:-1,IW)+WTS(L+1:L*2,IW) ! FOLD INTO 1ST L
      WTS(L+1:NO-L,IW) =WTS(L*2+1:NO,IW)               ! SHIFT THE REST LEFT
      WTS(NOP-L:NO,IW)=0 ! SET TRAILING L ELEMENTS TO ZERO
      INS1(I)=1          ! RESET INDEX OF FIRST POINT OF STENCIL
   ENDIF
   l=ins1(i)+nom-ism
   if(l > 0)then ! Fold right overlap of L elements back inside:
      wts(nop-l:no,iw)=wts(no:nop-l:-1,iw)+wts(nop-l*2:no-l,iw) ! Fold last L
      wts(l+1:no-l,iw)=wts(1:no-l*2,iw)                         ! Shift right
      wts(1:l,iw)=0      ! Set first L elements to zero
      ins1(i)=ism-nom    ! reset index of first point of stencil
   endif
enddo
ns=ism

!----------------------------------------------------------------------------
! Use Newton-Raphson iterations to locate the g-coordinates of all this
! segment's s-grid points. Then interpolate the function c to each of
! these s-grid points. (Note that, in the present context, the
! s- and g-grids are the ones NOT staggered wrt the domain boundaries.)
!----------------------------------------------------------------------------
ie=ig0
do i=is0+1,ism-1 ! Loop over s-grid target points interior to this segment
   et=i
!----------------------------------------------------------------------------
! Find the g-grid interval containing this target: 
!----------------------------------------------------------------------------
   do iep=ie+1,igm-1;  if(sofg(iep) > et)exit; enddo
   do ie=iep-1,ig1,-1; if(sofg(ie) <= et)exit; enddo

   ie1=ie-nohm;   ien=ie+noh   ! <-- Set the interpolation stencil range:

   r=(et-sofg(ie))/(sofg(ie+1)-sofg(ie)) ! Linearly estimate interval fraction

!----------------------------------------------------------------------------
! Perform Newton-Raphson iterations to refine interval fraction, r:
!----------------------------------------------------------------------------
   do init=1,nnit
      call lagw(hunit,r,q,wt,dwt,no) ! Get Lagrange weights, wt and d(wt)/dg
      estar =dot_product(wt, sofg(ie1:ien))-et ! <- Residual error, estar.
      destar=dot_product(dwt,sofg(ie1:ien))    ! <- d(estar)/dg.
      dr=-estar/destar                         ! <- Newton correction to r
      r=r+dr                                   ! <- Refined estimate, r
      if(abs(dr) <= rcrit)goto 1               ! <- Converged enough yet?
   enddo
 ! stop 'Too many Newton iterations'           ! <- It never convergenced! 
    write(*,*)' Too many Newton iterations'           ! <- It never convergenced! 
    ns=-1
    return
1  wt=wt+dr*dwt                                ! <- Final refinement to wt
   cofs(i)=dot_product(wt, cofg(ie1:ien))      ! <- Interpolate c(s)
enddo
cofs(ism)=cofg(igm)                            ! <- End value directly
end subroutine jfit

!============================================================================
subroutine stog(ns,ng,ins1,wts, as,ag) 
!============================================================================
! R. J. Purser NCEP 2001
! Forward interpolation from s-grid to g-grid
! --> ns,ng: sizes of s and g grids
! --> ins1 : array of 1st stencil indices (s-grid) for each target (g) point.
! --> wts  : interpolation weights for each target (g-grid point).
! --> as   : s-grid array of source data.
! <-- ag   : g-grid array of interpolated target data.
!============================================================================
use vkind
implicit none
integer, parameter                      :: noh=3,no=noh*2,nom=no-1
integer,                  intent(IN   ) :: ns,ng
integer, dimension(ng),   intent(IN   ) :: ins1
real(vp),dimension(no,ng),intent(IN   ) :: wts
real(vp),dimension(ns),   intent(IN   ) :: as
real(vp),dimension(ng),   intent(OUT  ) :: ag
!----------------------------------------------------------------------------
integer                                 :: i,is,iw
!============================================================================
iw=0
ag=0
do i=1,ng
   is=ins1(i)
   if(is>0)then
      iw=iw+1
      ag(i)=dot_product(wts(:,iw),as(is:is+nom))
   else
      ag(i)=as(-is)
   endif
enddo
end subroutine stog

!============================================================================
subroutine stogt(ns,ng,ins1,wts, as,ag) 
!============================================================================
! Perform the transpose of the operation defined by stog
! R. J. Purser NCEP 2001
!============================================================================
use vkind
implicit none
integer, parameter                      :: noh=3,no=noh*2,nom=no-1
integer,                  intent(IN   ) :: ns,ng
integer, dimension(ng),   intent(IN   ) :: ins1
real(vp),dimension(no,ng),intent(IN   ) :: wts
real(vp),dimension(ns),   intent(OUT  ) :: as
real(vp),dimension(ng),   intent(IN   ) :: ag
!----------------------------------------------------------------------------
integer                                 :: i,is,iw
!============================================================================
iw=0
as=0
do i=1,ng
   is=ins1(i)
   if(is>0)then
      iw=iw+1
      as(is:is+nom)=as(is:is+nom)+wts(:,iw)*ag(i)
   else
      as(-is)=as(-is)+ag(i)
   endif
enddo
end subroutine stogt

!ENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDEND
!						********************
!						* rfit.f90         *
!						* PURSER 1994/1999 *
!                                               *   FUJITA 1999    *
!						********************
!ENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDENDEND