      SUBROUTINE LTPBVP(INTL,IORD,NPTS,FPRM,COFX,FRHS,YSOL,WORK,PRTB,   
     1IERR)                                                             
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C                                                                       
C                                                                       
C     PURPOSE . . .                                                     
C     LTPBVP SOLVES FOR EITHER (1) THE SECOND ORDER FINITE DIFFERENCE   
C     APPROXIMATION OR (2) A FOURTH ORDER APPROXIMATION TO THE SOLUTION 
C     OF THE LINEAR (SINGULAR OR NON-SINGULAR) TWO POINT BOUNDARY       
C     VALUE PROBLEM                                                     
C                                                                       
C          A(X)*YXX(X) + B(X)*YX(X) + C(X)*YF(X) = F(X)                 
C                                                                       
C     ON THE REGION A.LE.X.LE.B.  THE POSSIBLE BOUNDARY CONDTIIONS ARE  
C     (1) PERIODIC, YF(X+B-A) = YF(X) FOR ALL X (OR)                    
C     (2) MIXED, ALPHA(A)*YX(A)+BETA(A)*YF(A) = BDY(A)                  
C              ALPHA(B)*YX(B)+BETA(B)*YF(B) = BDY(B).                   
C     ACCESS CARDS . . .                                                
C     *FORTRAN,S=PLIB,SN=LTPBVP,P=XLIB                                  
C     USAGE . . .                                                       
C     CALL LTPBVP(INTL,IORD,NPTS,FPRM,COFX,FRHS,YSOL,WORK,PRTB)         
C     ARGUMENTS . . .                                                   
C     ON INPUT . . .                                                    
C     INTL                                                              
C          = 0 ON INITIAL ENTRY TO LTPBVP                               
C          = 1 IF LTPBVP HAS BEEN CALLED BEFORE AND ONLY THE VALUES IN  
C              FRHS (SEE BELOW) HAVE CHANGED.                           
C     IORD                                                              
C          = 2 IF A SECOND ORDER APPROXIMATION IS SOUGHT.               
C          = 4 IF A FOURTH ORDER APPROXIMATION IS SOUGHT.               
C     NPTS                                                              
C          THE NUMBER OF EQUALLY SPACED POINTS IN THE INTERVAL (A,B).   
C          THESE ARE GIVEN BY X(I) = A+(I-1)*DLX, I=1,...,NPTS WHERE    
C          DLX = (B-A)/(NPTS-1).  NPTS MUST BE AT LEAST 6.              
C     FPRM                                                              
C          A FLOATING POINT VECTOR OF LENGTH 8 USED TO EFFICIENTLY INPUT
C          PARAMETERS.  FPRM IS SET INTERNALLY AND DEFINED AS FOLLOWS . 
C          A=FPRM(1),B=FPRM(2)                                          
C               THE RANGE OF THE INDEPENDENT VARIABLE X.  A MUST BE     
C               LESS THAN B.                                            
C          BDA=FPRM(3),BDB=FPRM(4),ALPHAA=FPRM(5),BETAA=FPRM(6),        
C          ALPHAB=FPRM(7),BETAB=FPRM(8)                                 
C               SCALARS USED TO INPUT MIXED OR PERIODIC BOUNDARY        
C               CONDITIONS AT X=A,X=B (SEE PURPOSE SECTION ABOVE).      
C               HERE BDA=BDY(A),BDB=BDY(B),ALPHAA=ALPHA(A),             
C               BETAA=BETA(A),ALPHAB=ALPHA(B),BETAB=BETA(B).            
C               PERIODIC BOUNDARY CONDITIONS ARE FLAGGED BY INPUTING    
C               FPRM(5)=FPRM(6)=FPRM(7)=FPRM(8)=0.0 (FLOATING POINT     
C               ZERO).  IN THIS CASE FPRM(3),FPRM(4) ARE DUMMY          
C               PARAMETERS.                                             
C     COFX                                                              
C          A SUBPROGRAM WITH PARAMETERS X,AF,BF,CF WHICH EVALUATES      
C          A(X),B(X),C(X) (SEE FORM OF EQUATION IN SECTION PURPOSE).    
C          COFX MUST BE DECLARED EXTERNAL IN THE CALLING ROUTINE.       
C          NOTE . . .                                                   
C          THE COEFFICIENTS PROVIDED MAY LEAD TO A TRIDIAGONAL MATRIX   
C          WHICH IS NOT DIAGONALLY DOMINANT (PARTICULARLY IF BF         
C          IS LARGE RELATIVE TO AF FOR SOME X) AND PRODUCE ERRONEOUS    
C          RESULTS.  THIS CANNOT HAPPEN IN THE LIMIT AS DLX APPROACHES  
C          ZERO AND HENCE MAYBE REMEDIED BY INCREASING NPTS.            
C     FRHS                                                              
C          A VECTOR OF LENGTH NPTS THAT SPECIFIES THE RIGHT HAND        
C          SIDE OF THE EQUATION AT EACH X(I) , I=1,...,NPTS             
C          (I.E., FRHS(I)=F(X) FOR X=X(I)).                             
C     WORK                                                              
C          A VECTOR THAT MUST BE PROVIDED FOR WORK SPACE.  WORK         
C          SHOULD BE AT LEAST 8*NPTS WORDS LONG IN THE CALLING          
C          ROUTINE.                                                     
C     ON OUTPUT . . .                                                   
C     YSOL                                                              
C          A VECTOR OF LENGTH NPTS THAT CONTAINS THE APPROXIMATE        
C          SOLUTION TO THE LINEAR TWO POINT BOUNDARY VALUE PROBLEM.     
C          YSOL(I) IS THE APPROXIMATION TO YF(X) AT X=X(I) FOR          
C          I=1,...,NPTS.  THE APPROXIMATION HAS ERROR O(DLX**IORD).     
C          YSOL MAY BE EQUIVALENCED WITH FRHS IF IORD=2.                
C     WORK                                                              
C          CONTAINS INTERMEDIATE VALUES THAT MUST NOT BE DESTROYED IF   
C          LTPBVP IS CALLED AGAIN WITH INTL=1.                          
C     PRTB                                                              
C          ORDINARILY PRTB=0.0 IS RETURNED.  HOWEVER IF PERIODIC OR     
C          PURE DERIVATIVE (ALPHAA=FPRM(5)=0.0, ALPHAB=FPRM(7)=0.0)     
C          BOUNDARY CONDITIONS ARE SPECIFIED AND IF C(X)=CF=0.0         
C          FOR ALL X THEN THE LINEAR TWO POINT BOUNDARY VALUE           
C          PROBLEM IS SINGULAR.  IN THIS CASE A SOLUTION TO THE         
C          DISCRETIZED TRI-DIAGONAL SYSTEM MAY NOT EXITT.  PTRB         
C          IS A CONSTANT CALCULATED AND SUBTRACTED FROM EACH COMPONENT  
C          OF FRHS WHICH INSURES THAT A SOLUTION DOES EXIST.  LTPBVP    
C          THEN SOLVES THE RESULTING CONSISTENT MATRIX EQUATION.        
C          THIS YIELDS A WEIGHTED MINIMAL LEAST SQUARES SOLUTION        
C          TO THE ORIGINAL PROBLEM.  THIS MAY OCCURR FOR IORD=2 OR 4.   
C     IERR                                                              
C          AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS        
C          = 0 IF NO ERROR                                              
C          = 1 IF A.GE.B                                                
C          = 2 IF NPTS.LT.6                                             
C          = 3 IF IORD.NE.2,4                                           
C          = 4 IF A(X)=0 FOR SOME X=X(I),1.LE.I.LE.NPTS.                
C          = 5 IF INTL.NE.0,1                                           
C          = 6 IF THERE IS AMBIGUITY FLAGGING PERIODIC B.C.             
C              FOR EXAMPLE HOW IS ALPHAA=BETAA=0.0,ALPHAB.NE.0.         
C              TO BE INTERPRETED.                                       
C     ENTRY POINTS                                                      
C          LTPBVP,TWPT,CHKPR,DEFERR,DXDIF,ORTHO,MINSO,CHKSN,TRISLV      
C     SPECIAL CONDITIONS . . . NONE                                     
C     COMMON BLOCKS . . . LWTWP                                         
C     I/O . . . NONE                                                    
C     PRECISION . . . SINGLE                                            
C     REQUIRED ULIB ROUTINES . . . NONE                                 
C     SPECIALIST . . . JOHN C. ADAMS, NCAR, BOULDER, COLO. 80303        
C     LANGUAGE . . . FORTRAN                                            
C     HISTORY . . .                                                     
C          DEVELOPED AS A ONE DIMENSIONAL VERSION OF THE                
C          ULIB ROUTINE SEPELI DURING 1976-77 BY JOHN ADAMS.            
C     ALGORITHIM . . .                                                  
C          LTPBVP AUTOMATICALLY DISCRETIZES THE LINEAR TWO POINT        
C          BOUNDARY VALUE PROBLEM INPUT.  THE RESULTING TRIDIAGONAL     
C          SYSTEM IS SOLVED DIRECTLY (WITHOUT PIVOTING).  THE FOURTH    
C          ORDER APPROXIMATION IS OBTAINED USING DEFERRED CORRECTIONS.  
C          SINGULAR TWO POINT ODES ARE DETECTED AND THE RESULTING       
C          TRIDIAGONAL MATRIX EQUATION IS MADE CONSISTENT BY            
C          ORTHOGONALIZING THE RIGHT HAND SIDE WITH RESPECT TO THE      
C          NULL SPACE OF THE TRANSPOSE OF THE COEFFICIENT MATRIX.       
C          DIRECTLY SOLVING THE RESULTANT CONSISTENT SYSTEM YIELDS      
C          A WEIGHTED MINIMAL LEAST SQUARES SOLUTION TO THE ORIGINAL    
C          EQUATION (SEE PARAMETER PRTB).                               
C     SPACE REQUIRED . . . 2423(BASE8) = 1235 DECIMAL WORDS             
C     PORTABILITY . . . THERE ARE NO MACHINE DEPENDENT CONSTANTS        
C     REQUIRED RESIDENT ROUTINES . . . NONE                             
C                                                                       
C                                                                       
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
      DIMENSION YSOL(NPTS),FRHS(NPTS),WORK(1),FPRM(8)                   
      COMMON/ LWTWP/A,B,BDA,BDB,ALPHAA,BETAA,ALPHAB,BETAB,MIT,M,MS,K,IS,
     1  INT ,IORDER,IXA,IXB,DLX,DLXX,TDLX3,DLX4                         

c Define the subroutine COFX to be external for passing as argument below.
      EXTERNAL COFX

C SET PARAMETERS INTERNALLY                                         
      INT = INTL                                                        
      IORDER = IORD                                                     
      K = NPTS                                                          
      M=K-1                                                             
      A=FPRM(1)                                                         
      B=FPRM(2)                                                         
      BDA=FPRM(3)                                                       
      BDB=FPRM(4)                                                       
      ALPHAA=FPRM(5)                                                    
      BETAA=FPRM(6)                                                     
      ALPHAB=FPRM(7)                                                    
      BETAB=FPRM(8)                                                     
C     CHECK INPUT PARAMETERS ON AN INITIAL (INTL=0) CALL                
C                                                                       
      IERR=5                                                            
      IF(INT.NE.0 .AND. INT.NE.1) RETURN                                
      IERR=0                                                            
      IF(INTL.EQ.0) CALL CHKPR(COFX,IERR)                               
      IF(IERR.NE.0) RETURN                                              
C                                                                       
C     SET WORK SPACE INDICES                                            
C                                                                       
      I1=1                                                              
      I2=I1+K                                                           
      I3=I2+K                                                           
      I4=I3+K                                                           
      I5=I4+K                                                           
C     COMPUTE APPROXIMATION                                             
      CALL TWPT(COFX,FRHS,YSOL,WORK(I5),PRTB,IERR,WORK(I1),WORK(I2),    
     1WORK(I3),WORK(I4))                                                
      RETURN                                                            
      END                                                               


      SUBROUTINE TWPT(COFX,FRHS,YSOL,WORK,PRTB,IERR,AM,BM,CM,ZM)        
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C                                                                       
C     THIS SUBROUTINE DISCRETIZES THE TWO PT. B.V. PROBLEM INCORPORATING
C     BOUNDARY CONDITIONS APPROPRIATELY.  SECOND OR FOURTH ORDER        
C     APPROXIMATONS ARE COMPUTED FOR SINGULAR OR NON-SINGULAR           
C     DIFFERENTIAL OPERATORS                                            
C                                                                       
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
      DIMENSION FRHS(1),YSOL(1),WORK(1),AM(1),BM(1),CM(1),ZM(1)         
      COMMON/ LWTWP/A,B,BDA,BDB,ALPHAA,BETAA,ALPHAB,BETAB,MIT,M,MS,K,IS,
     1  INT ,IORDER,IXA,IXB,DLX,DLXX,TDLX3,DLX4                         
      LOGICAL SINGLR                                                    

	external COFX

      IF(INT.EQ.1) GO TO 130                                            
      I1=1                                                              
      I2=I1+K                                                           
      I3=I2+K                                                           
      I4=I3+K                                                           
C                                                                       
C     SET DLX AND SIZE OF TRI-DIAGONAL MATRIX                           
      MIT=K-1                                                           
      IF(IXA.EQ.1 .AND. IXB.EQ.1) MIT=K-2                               
      IF(IXA.EQ.2 .AND. IXB.EQ.2) MIT=K                                 
      DLX=(B-A)/FLOAT(M)                                                
      DLXX=DLX*DLX                                                      
      TDLX3=(DLXX+DLXX)*DLX                                             
      DLX4=DLXX*DLXX                                                    
C                                                                       
C     SET SUBSCRIPT LIMITS                                              
      IS=1                                                              
      IF(IXA.EQ.1) IS=2                                                 
      MS=MIT+IS-1                                                       
C                                                                       
C     CONSTRUCT MATRIX                                                  
      DO 111 I=1,MIT                                                    
      X=A  +FLOAT(IS+I-2)*DLX                                           
      CALL COFX(X,AI,BI,CI)                                             
      AXI=(AI/DLX-0.5*BI)/DLX                                           
      BXI=-2.0*AI/DLXX+CI                                               
      CXI=(AI/DLX+0.5*BI)/DLX                                           
      AM(I)=AXI                                                         
      BM(I)=BXI                                                         
      CM(I)=CXI                                                         
  111 CONTINUE                                                          
C                                                                       
C     ADJUST AT X=A, X=B UNLESS PERIODIC                                
      IF(IXA.EQ.0) GO TO 117                                            
      AX1=AM(1)                                                         
      CXM=CM(MIT)                                                       
      AM(1)=0.0                                                         
      CM(MIT)=0.0                                                       
      IF(IXA.NE.2) GO TO 118                                            
C     MIXED AT X=A                                                      
      BM(1)=BM(1)+2.0*DLX*BETAA*AX1/ALPHAA                              
      CM(1)=CM(1)+AX1                                                   
  118 IF(IXB.NE.2) GO TO 117                                            
C     MIXED AT X=B                                                      
      AM(MIT)=AM(MIT)+CXM                                               
      BM(MIT)=BM(MIT)-2.0*DLX*BETAB*CXM/ALPHAB                          
  117 CONTINUE                                                          

C     CHECK IF OPERATOR IS SINGULAR                                     
      CALL CHKSN(COFX,SINGLR)                                           
C                                                                       
      IF(.NOT.SINGLR) GO TO 140                                         
C     COMPUTE POSITIVE VECTOR IN NULL SPACE OF TRANSPOSE                
C                                                                       
      KSING=0                                                           
      CALL TRISLV(AM,BM,CM,WORK(I1),WORK(I2),WORK(I3),WORK(I4),ZM,MIT,  
     1KSING,SINGLR,INT)                                                 
      INT=1                                                             
  140 CONTINUE                                                          
  130 CONTINUE                                                          
C                                                                       
C     SET  SPECIFIED (DIRICHLET) BOUNDARIES                             
      IF(IXA.EQ.1) YSOL(1) = BDA/BETAA                                  
      IF(IXB.EQ.1) YSOL(K) = BDB/BETAB                                  
C                                                                       
C     SET RHS VALUES IN YSOL ON INTERIOR AND NON-SPECIFIED BNDYS        
C                                                                       
      DO 102 I=IS,MS                                                    
  102 YSOL(I)=FRHS(I)                                                   
C                                                                       
C     ADJUST R.H.S. IN YSOL                                             
      IF(IXA.EQ.2) YSOL(IS)=YSOL(IS)+2.0*DLX*AX1*BDA/ALPHAA             
      IF(IXA.EQ.1) YSOL(IS) = YSOL(IS)-AX1*YSOL(1)                      
      IF(IXB.EQ.2) YSOL(MS) = YSOL(MS)-2.0*DLX*CXM*BDB/ALPHAB           
      IF(IXB.EQ.1) YSOL(MS) = YSOL(MS)-CXM*YSOL(K)                      
C                                                                       
      IORD=IORDER                                                       
C     SAVE ADJUSTED CORNERS IN FRHS IF IORD=4                           
      IF(IORD.EQ.2) GO TO 135                                           
      FRHS(IS)=YSOL(IS)                                                 
      FRHS(MS)=YSOL(MS)                                                 
  135 CONTINUE                                                          
      PRTB=0.0                                                          
C                                                                       
  136 CONTINUE                                                          
C     MAKE R.H.S. CONSISTENT IF SINGULAR                                
      IF(SINGLR) CALL ORTHO(YSOL,ZM,PRTB)                               
C                                                                       
C     SOLVE TRIDIAGONAL SYSTEM                                          
      KSING=1                                                           
      CALL TRISLV(AM,BM,CM,WORK(I1),WORK(I2),WORK(I3),WORK(I4),YSOL(IS),
     1MIT,KSING,SINGLR,INT)                                             
      INT=1                                                             
C     SET PERIODIC B.C. IF NECESSARY                                    
      IF(IXA.EQ.0) YSOL(K) = YSOL(1)                                    
C                                                                       
C     MINIMIZE SOL W.R.T. WEIGHTED LEAST SQUARES NORM IF SINGULAR       
      IF(SINGLR) CALL MINSO(YSOL,ZM,PRBB)                               
C     RETURN IF 2ND ORDER APPROX IS FLAGGED.  OTHERWISE                 
C     PROCEED WITH DEFERRED CORRECTIONS FOR FOURTH ORDER APPROX.        
      IF(IORD.EQ.2) RETURN                                              
      IORD=2                                                            
C                                                                       
C     COMPUTE NEW R.H.S FOR DEFERRED CORRECTION FOURTH ORDER SOL        
      CALL DEFERR(COFX,YSOL,FRHS)                                       

      GO TO 136                                                         
      END                                                               
      SUBROUTINE CHKPR(COFX,IERROR)                                     
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C                                                                       
C     THIS SUBROUTINE CHECKS INPUT PARAMETERS FOR ERRORS AND SETS       
C     INTERNAL FLAGS FOR BOUNDARY CONDITIONS IN IXA,IXB                 
C     IXA=IXB=0 MEANS PERIODIC B.C.                                     
C     IXA=1,IXB=1 MEANS B.C. AT X=A,X=B IS SPECIFIED                    
C     IXA=2,IXB=2 MEANS B.C. AT X=A,X=B IS FULLY MIXED (ALPHAA.NE.0.    
C     , ALPHAB.NE.0.)                                                   
C                                                                       
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
      LOGICAL LA,LB                                                     
      COMMON/ LWTWP/A,B,BDA,BDB,ALPHAA,BETAA,ALPHAB,BETAB,MIT,M,MS,K,IS,
     1  INT ,IORDER,IXA,IXB,DLX,DLXX,TDLX3,DLX4                         
      IERROR=1                                                          
      IF(A.GE.B) RETURN                                                 
      IERROR=2                                                          
      IF(K.LT.6) RETURN                                                 
      IERROR=3                                                          
      IF(IORDER.NE.2 .AND. IORDER.NE.4) RETURN                          
C     CHECK THAT THE COEFFICIENT OF YXX IS NOT ZERO ON THE INTERIOR     
      IERROR=4                                                          
      DLX=(B-A)/FLOAT(M)                                                
      DO 1 I=2,M                                                        
      X=A+FLOAT(I-1)*DLX                                                
      CALL COFX(X,AI,BI,CI)                                             
      IF(AI.EQ.0.0) RETURN                                              
    1 CONTINUE                                                          
C     CHECK AND FLAG (VIA IXA,IXB) BOUNDARY CONDITIONS                  
      LA=(ALPHAA.EQ.0.0 .AND. BETAA.EQ.0.0)                             
      LB=(ALPHAB.EQ.0.0 .AND. BETAB.EQ.0.0)                             
      IERROR=6                                                          
      IF((.NOT.LA .AND. LB) .OR. (LA .AND. .NOT.LB)) RETURN             
      IERROR=0                                                          
      IXA=0                                                             
      IXB=0                                                             
      IF(LA .AND. LB) RETURN                                            
      LA=(ALPHAA.NE.0.0)                                                
      IXA=1                                                             
      IF (LA) IXA=2                                                     
      LB = (ALPHAB.NE.0.0)                                              
      IXB=1                                                             
      IF(LB) IXB=2                                                      
      RETURN                                                            
      END                                                               

      SUBROUTINE DEFERR(COFX,YSOL,FRHS)                                 
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C                                                                       
C     THIS SUBROUTINE APPROXIMATES THE TRUNCATION ERROR                 
C     GIVEN BY DLX**2*TX WHERE TX=AF(X)*YXXXX/12.0+BF(X)*YXXX/6.0       
C     ON THE INTERIOR AND AT THE BOUNDARIES IF                          
C     PERIODIC. TX=AF(X)/3.0*(YXXXX/4.0+YXXX/(3.0*DLX)) IF MIXED AT X=A 
C     AND TX=AF(B)/3.0*(YXXXX/4.0-YXXX/(3.*DLX))IF MIXED AT X=B         
C     TX=0. AT SPECIFIED BOUNDARIES.                                    
C     THE SECOND ORDER APPROXIMATION IN YSOL IS USED                    
C     TO ESTIMATE THE THIRD, FOURTH DERIVATIVES YXXX,YXXXX              
C     VIA FINITE DIFFERENCING. THIS SUBROUTINE RETURNS                  
C     YSOL=FRHS+DLX**2*TX                                               
C                                                                       
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
      COMMON/ LWTWP/A,B,BDA,BDB,ALPHAA,BETAA,ALPHAB,BETAB,MIT,M,MS,K,IS,
     1  INT ,IORDER,IXA,IXB,DLX,DLXX,TDLX3,DLX4                         
      DIMENSION YSOL(1),FRHS(1)                                         

	external COFX

C     ESTIMATE TRUNCATION ERROR AND ADD TO FRHS                         
C     THE LOOP TO DO THIS IS 'UNROLLED ' TO OPTIMIZE                    
C     X.LE.A+DLX                                                        
      DO 103 I=IS,2                                                     
      X=A+FLOAT(I-1)*DLX                                                
      CALL COFX(X,AI,BI,CI)                                             
      CALL DXDIF(YSOL,I,YXXX,YXXXX)                                     
C     PERIODIC OR INTERIOR POINT TRUNCATION ERROR                       
      IF(IXA.EQ.0 .OR. I.EQ.2) TX=(AI*YXXXX*0.5+BI*YXXX)/6.0            
C     NON-PERIODIC BOUNDARY POINT TRUNCATION ERROR                      
      IF(IXA.NE.0 .AND. I.EQ.1) TX=AI/3.0*(YXXXX/4.0+YXXX/DLX)          
      FRHS(I)=FRHS(I)+DLXX*TX                                           
  103 CONTINUE                                                          
C     A+DLX.LT.X.LT.B-DLX                                               
      MM1=M-1                                                           
      DO 105 I=3,MM1                                                    
      X=A+FLOAT(I-1)*DLX                                                
      CALL COFX(X,AI,BI,CI)                                             
      YXXX=(-YSOL(I-2)+2.*(YSOL(I-1)-YSOL(I+1))+YSOL(I+2))/TDLX3        
      YXXXX=(YSOL(I-2)-4.*(YSOL(I-1)+YSOL(I+1))+6.*YSOL(I)+YSOL(I+2))/  
     1DLX4                                                              
      TX=(AI*YXXXX*0.5+BI*YXXX)/6.0                                     
      FRHS(I)=FRHS(I)+DLXX*TX                                           
  105 CONTINUE                                                          
C     B-DLX.LE.X                                                        
      DO 104 I=M,MS                                                     
      X=A+FLOAT(I-1)*DLX                                                
      CALL COFX(X,AI,BI,CI)                                             
      CALL DXDIF(YSOL,I,YXXX,YXXXX)                                     
      IF(I.LT.K) TX=(AI*YXXXX*0.5+BI*YXXX)/6.0                          
C     RIGHT BOUNDARY TRUNCATION ERROR                                   
      IF(I.EQ.K) TX=AI/3.0*(YXXXX/4.0-YXXX/DLX)                         
      FRHS(I)=FRHS(I)+DLXX*TX                                           
  104 CONTINUE                                                          
C     RESET RHS IN YSOL                                                 
      DO 106 I=IS,MS                                                    
      YSOL(I)=FRHS(I)                                                   
  106 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE DXDIF(Y,I,YXXX,YXXXX)                                  
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C                                                                       
C     ESTIMATE PARTIAL DERIVATIVES AT BOUNDARIES AND ONE POINT IN       
C     FROM BOUNDARIES USING NON-SYMMETRIC DIFFERENCE FORMULA FOR        
C     NON-PERIODIC BOUNDARY CONDITION AND SYMMETRIC FORMULA             
C     WITH ADJUSTED INDICES FOR PERIODIC B.C.                           
C                                                                       
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
      DIMENSION Y(1)                                                    
      COMMON/ LWTWP/A,B,BDA,BDB,ALPHAA,BETAA,ALPHAB,BETAB,MIT,M,MS,K,IS,
     1  INT ,IORDER,IXA,IXB,DLX,DLXX,TDLX3,DLX4                         
      IF(I.EQ.1) GO TO 101                                              
      IF(I.EQ.2) GO TO 103                                              
      IF(I.EQ.K-1) GO TO 106                                            
      IF(I.EQ.K) GO TO 108                                              
C     AT X=A                                                            
  101 IF(IXA.EQ.0) GO TO 102                                            
      YXXX=(-5.*Y(1)+18.*Y(2)-24.*Y(3)+14.*Y(4)-3.*Y(5))/TDLX3          
      YXXXX=(3.*Y(1)-14.*Y(2)+26.*Y(3)-24.*Y(4)+11.*Y(5)-2.*Y(6))/DLX4  
      RETURN                                                            
  102 YXXX=(-Y(K-2)+2.*Y(K-1)-2.*Y(2)+Y(3))/TDLX3                       
      YXXXX=(Y(K-2)-4.*Y(K-1)+6.*Y(1)-4.*Y(2)+Y(3))/DLX4                
      RETURN                                                            
C     AT X=A+DLX                                                        
  103 IF(IXA.EQ.0) GO TO 104                                            
      YXXX=(-3.*Y(1)+10.*Y(2)-12.*Y(3)+6.*Y(4)-Y(5))/TDLX3              
      YXXXX=(2.*Y(1)-9.*Y(2)+16.*Y(3)-14.*Y(4)+6.*Y(5)-Y(6))/DLX4       
      RETURN                                                            
  104 YXXX=(-Y(K-1)+2.*Y(1)-2.*Y(3)+Y(4))/TDLX3                         
      YXXXX=(Y(K-1)-4.*Y(1)+6.*Y(2)-4.*Y(3)+Y(4))/DLX4                  
      RETURN                                                            
C     AT X=B-DLX                                                        
  106 IF(IXA.EQ.0) GO TO 107                                            
      YXXX=(Y(K-4)-6.*Y(K-3)+12.*Y(K-2)-10.*Y(K-1)+3.*Y(K))/TDLX3       
      YXXXX=(-Y(K-5)+6.*Y(K-4)-14.*Y(K-3)+16.*Y(K-2)-9.*Y(K-1)+2.*Y(K))/
     1DLX4                                                              
      RETURN                                                            
  107 YXXX=(-Y(K-3)+2.*Y(K-2)-2.*Y(1)+Y(2))/TDLX3                       
      YXXXX=(Y(K-3)-4.*Y(K-2)+6.*Y(K-1)-4.*Y(1)+Y(2))/DLX4              
      RETURN                                                            
C     AT X=B                                                            
  108 YXXX=  (3.*Y(K-4)-14.*Y(K-3)+24.*Y(K-2)-18.*Y(K-1)+5.*Y(K))/TDLX3 
      YXXXX=(-2.*Y(K-5)+11.*Y(K-4)-24.*Y(K-3)+26.*Y(K-2)-14.*Y(K-1)+3.*Y
     1(K))/DLX4                                                         
      RETURN                                                            
      END                                                               
      SUBROUTINE ORTHO(YSOL,ZM,PERTRB)                                  
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C                                                                       
C     THIS SUBROUTINE ADJUSTS THE RHS SO EQUATIONS ARE CONSISTENT       
C                                                                       
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
      DIMENSION YSOL(1),ZM(1)                                           
      COMMON/ LWTWP/A,B,BDA,BDB,ALPHAA,BETAA,ALPHAB,BETAB,MIT,M,MS,K,IS,
     1  INT ,IORDER,IXA,IXB,DLX,DLXX,TDLX3,DLX4                         
      ISTRT=IS                                                          
      IFNL=MS                                                           
C                                                                       
C     COMPUTE WEIGHTED INNER PRODUCT                                    
      YTE=0.0                                                           
      ETE=0.0                                                           
      DO 103 I=IS,MS                                                    
      II=I-IS+1                                                         
      ETE=ETE+ZM(II)                                                    
      YTE=YTE+ YSOL(I)*ZM(II)                                           
  103 CONTINUE                                                          
C                                                                       
C     SET PERTURBATION PARAMETER AND SUBTRACT OFF                       
C                                                                       
      PERTRB=YTE/ETE                                                    
      DO 105 I=ISTRT,IFNL                                               
      YSOL(I)=YSOL(I)-PERTRB                                            
  105 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE MINSO(YSOL,ZM,PERTRB)                                  
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C                                                                       
C     MINIMIZE FINAL SOLUTION WITH RESPECT TO WEIGHTED LEAST SQUARES    
C     NORM DEFINED BY ZM                                                
C                                                                       
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
      DIMENSION YSOL(1),ZM(1)                                           
      COMMON/ LWTWP/A,B,BDA,BDB,ALPHAA,BETAA,ALPHAB,BETAB,MIT,M,MS,K,IS,
     1  INT,IORDER,IXA,IXB,DLX,DLXX,TDLX3,DLX4                          
      YTE=0.0                                                           
      ETE=0.0                                                           
      DO 103 I=IS,MS                                                    
      II=I-IS+1                                                         
      ETE=ETE+ZM(II)                                                    
       YTE=YTE+YSOL(I)*ZM(II)                                           
  103 CONTINUE                                                          
      PERTRB=YTE/ETE                                                    
      DO 105 I=1,K                                                      
      YSOL(I)=YSOL(I)-PERTRB                                            
  105 CONTINUE                                                          
      RETURN                                                            
      END                                                               

      SUBROUTINE CHKSN(COFX,SINGLR)                                     
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C                                                                       
C     THIS SUBROUTINE CHECKS IF THE LINEAR TWO POINT B.V.               
C     PROBLEM IS SINGULAR.                                              
C                                                                       
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
      COMMON/ LWTWP/A,B,BDA,BDB,ALPHAA,BETAA,ALPHAB,BETAB,MIT,M,MS,K,IS,
     1  INT ,IORDER,IXA,IXB,DLX,DLXX,TDLX3,DLX4                         

	external COFX
      LOGICAL SINGLR                                                    

      SINGLR=.FALSE.                                                    

C     RETURN IF SCALR MULTIPLYING Y IS NONZERO AT X=A OR X=B            
      IF(BETAA.NE.0.0 .OR. BETAB.NE.0.0) RETURN                         
C     RETURN IF THE COEFFICIENT MULTIPLYING Y IN THE DIFFERENTIAL       
C     OPERATOR IS NONZERO                                               
      DO 103 I=IS,MS                                                    
      X=A+FLOAT(I-1)*DLX                                                
      CALL COFX(X,AI,BI,CI)                                             
      IF(CI.NE.0.0) RETURN                                              
  103 CONTINUE                                                          
C     DIFFERENTIAL OPERATOR IS SINGULAR IF THIS POINT IS REACHED        
      SINGLR=.TRUE.                                                     
      RETURN                                                            
      END                                                               
      SUBROUTINE TRISLV(A,B,C,D,E,U,V,R,N,KSING,SINGLR,INTL)            
C     ******************************************************************
C                                                                       
C     A SUBROUTINE TO SOLVE A TRI-DIAGONAL PERIODIC SYSTEM OF EQUATIONS.
C     OR COMPUTE A POSITIVE EIGENVECTOR IN THE NULL SPACE OF THE TRANSPO
C     THIS SUBROUTINE MAY BE MODIFIED SO THAT ONLY D NEED BE COMPUTED   
C     IN THE LU DECOMPOSITION TO SAVE SPACE.  THE VALUES IN E,U,V       
C     CAN BE COMPUTED RECURSIVELY AS NEEDED.   TO SAVE COMPUTATION      
C     THIS IS NOT DONE HERE.                                            
C     NO PIVOTING IS DONE IS THIS ROUTINE.  IT IS POSSIBLE THAT THE     
C     DISCRETIZATION OF THE TWO POINT B.V. PROBLEM WILL YIELD           
C     A TRIDIAGONAL SYSTEM WHICH IS NOT DIAGONALLY DOMINANT OR          
C     WHICH HAS ZERO DIAGONAL ELEMENT (E.G., CONSIDER YXX+2Y=F          
C     WHERE DLX=1).  IN THIS LATTER CASE A DIVISION BY ZERO WILL        
C     CAUSE FLOATING POINT OVERFLOW.                                    
C                                                                       
C     ******************************************************************
      DIMENSION A(N),B(N),C(N),D(N),E(N),U(N),V(N),R(N)                 
      LOGICAL SINGLR                                                    
C     BYPASS LU DECOMPOSITION IF INTL=1                                 
      IF(INTL.EQ.1) GO TO 10                                            
      NM1=N-1                                                           
      D(1)=B(1)                                                         
      U(1)=A(1)                                                         
      V(1)=C(N)/D(1)                                                    
      SUM=V(1)*U(1)                                                     
      DO 1 I=2,NM1                                                      
      E(I)=A(I)/D(I-1)                                                  
      D(I)=B(I)-E(I)*C(I-1)                                             
      U(I)=-E(I)*U(I-1)                                                 
      V(I)=-V(I-1)*C(I-1)/D(I)                                          
      SUM=SUM+V(I)*U(I)                                                 
    1 CONTINUE                                                          
C     READJUST SUM                                                      
      SUM=SUM-V(NM1)*U(NM1)                                             
      U(NM1)=C(NM1)+U(NM1)                                              
      E(N)=A(N)/D(NM1)+V(NM1)                                           
      D(N)=B(N)-(SUM+E(N)*U(NM1))                                       
   10 CONTINUE                                                          
      IF(KSING.EQ.0) GO TO 20                                           
C     FORWARD SWEEP                                                     
      SUM=V(1)*R(1)                                                     
      DO 2 I=2,NM1                                                      
      R(I)=R(I)-E(I)*R(I-1)                                             
      SUM=SUM+V(I)*R(I)                                                 
    2 CONTINUE                                                          
      SUM=SUM+(E(N)-V(NM1))*R(NM1)                                      
      R(N)=R(N)-SUM                                                     
C     BACKWARD SWEEP                                                    
C     START WITH DIVIDE ONLY IF ODE IS NOT DIAGNOSED AS SINGULAR        
      IF (.NOT. SINGLR) R(N) = R(N)/D(N)                                
      R(NM1)=(R(NM1)-U(NM1)*R(N))/D(NM1)                                
      DO 3 IB=2,NM1                                                     
      I=N-IB                                                            
      R(I)=(R(I)-C(I)*R(I+1)-U(I)*R(N))/D(I)                            
    3 CONTINUE                                                          
      RETURN                                                            
   20 CONTINUE                                                          
C     COMPUTE POSITIVE EIGENVECTOR                                      
      R(N)=+1.0                                                         
      R(NM1)=-E(N)*R(N)                                                 
      DO 4 IB=2,NM1                                                     
      I=N-IB                                                            
      R(I)=-(V(I)+E(I+1)*R(I+1))                                        
    4 CONTINUE                                                          
      RETURN                                                            
      END                                                               
