C  @(#)fitcon.for	19.1 (ESO-DMD) 02/25/03 13:17:46 
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
      SUBROUTINE FITCON(ALGOR,BNDTYP,ISTAT)                                    
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++        
C                                                                               
C.MODULE                                                                        
C       FIT                                                                     
C                                                                               
C.NAME                                                                          
C       FIT_CON                                                                 
C                                                                               
C.PURPOSE                                                                       
C   makes the interface between the MIDAS data structures and the NAG           
C   routines for non linear LSQ problems with constraints                       
C                                                                               
C.KEYWORDS                                                                      
C  Constrained Non-linear Least Square, Constrained Non-linear Fitting,         
C  Constrained Non-linear Regression, FIT/TABLE, FIT/IMAGE,                     
C  Corrected Gauss-Newton Method, Quasi-Newton Method, Modified                 
C  Gauss-Newton.                                                                
C                                                                               
C.DESCRIPTION                                                                   
C       Recovers qualifiers and parameters from MIDAS keywords;                 
C       Recovers bounds;                                                        
C       Give virtual memory for calculations in NAG routines;                   
C       Call the NAG routines for non-linear least squares problems;            
C       Display information messages;                                           
C       Compute errors on parameters.                                           
C ****
C CHECK CALLING SEQUENCE OF E04KAF,E04KBF,E04KCF,E04KDF
C     LAST TWO ARGUMENTS...
C ****
C                                                                               
C.RESTRICTIONS                                                                  
C       ...                                                                     
C                                                                               
C.LANGUAGE                                                                      
C       FORTRAN                                                                 
C                                                                               
C.CALLING SEQUENCE                                                              
C       CALL FITCON(ALGOR,BNDTYP,ISTAT)                                        
C                                                                               
C.INPUT PARAMETERS                                                              
C       ALGOR    (*)        CHARACTER   algorithm to be used.                   
C                           (CGNND,QN,MGN)                                      
C       BNDTYP   (*)        CHARACTER   Bounds types                            
C                           (P,G,I)                                             
C       ISTAT               INTEGER     exit status                             
C                                                                               
C.MODIFIED PARAMETERS                                                           
C       None                                                                    
C                                                                               
C.OUTPUT PARAMETERS                                                             
C       None                                                                    
C                                                                               
C.FILES                                                                         
C       FIT_NAG.INC/NOLIST                                                      
C                                                                               
C.MODULES CALLED                                                                
C       NAG library     (E04___)                                                
C       MINMSQ          (message display)                                       
C                                                                               
C.AUTHOR                                                                        
C       Ph. DEFERT,      Feb 1986                                               
C        M. Peron        Jul 1992                                             
C.MODIFICATIONS                                                                 
C                                                                               
C                                                                               
C.BIBLIOGRAPHY                                                                  
C       NAG Mark II Vol 3 (E04)                                                 
C       Algorithms for the Solution of non-linear least squares problem         
C          SIAM J. on Num. An.,15,977-992,1978                                  
C.VERSION
C 021122	last modif
C                                                                               
C-----------------------------------------------------------------------        
C 
C      IMPLICIT NONE                                                            
C     .. Scalar Arguments ..                                                    
      CHARACTER ALGOR* (*),BNDTYP* (*)                                          
      INTEGER ISTAT ,KUN,KNUL,TID,TID1
C     ..                                                                        
C     .. Scalars in Common ..                                                   
      INTEGER NRCOL,ISTAR                                                       
      CHARACTER WGTTYP*1                                                        
C     ..                                                                        
C     .. Arrays in Common ..                                                    
      INTEGER ICOL(10)                 
      INTEGER MADRID(1)                 
C     ..                                                                        
C     .. Local Scalars ..                                                       
      INTEGER I,IDAT,IFAIL1,IFAIL3,J,K,LIW,LW,NC,NROW,NRFEVA,           
     +        NRVM,NS,STATUS
      INTEGER IBOUND,LHSL,IMIN,IMAX,NRBND,ROUTIN                     
      INTEGER NAC,NAR,STAT,NAC1,NAR1,NBYT
      INTEGER*8 PTRW,PTRIW,PTRDLT,PTRHSL,PTRHSD,PTRSTA,PTRGRD

      DOUBLE PRECISION FSUMSQ
      DOUBLE PRECISION HMIN,HMAX

      LOGICAL VALID,ISEL,NOMPAR,VCMATPR                                         
      CHARACTER WGTKEY*10,LINE*78,BNDTAB*60                           

C     ..                                                                        
C     .. Local Arrays ..                                                        
      REAL XVAL(10)                                                             
      LOGICAL NULL(10)                                                          
C     ..                                                                        
C     ..                                                                        
C     .. External Files ..                                                      
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      INCLUDE 'MID_INCLUDE:FITNAGI.INC'                                 
      COMMON/VMR/MADRID
      COMMON /LSQFUN/ICOL,NRCOL,ISTAR,WGTTYP 
      INCLUDE 'MID_INCLUDE:FITNAGC.INC'                                 
      EXTERNAL E04LBS,E04JBQ,FUNCT3,FUNCT4,MINMON

      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
      DATA WGTKEY(1:10)/'FITCHAR   '/                                           

C                                                                               
C  Executable Statements                                                        
C                                                                               
C  Initialize parameters and compute number of fixed parameters                 
C                                                                               
      VCMATPR = METPAR(1) .LT. 0.                                               
      METPAR(1) = ABS(METPAR(1))                                                
      WRITE(6,*)METPAR(1)
      NRPFIX = 0                                                                
      DO 10 J = 1,NRPARAM                                                       
          IF (METPAR(2).GE.0) PARAM(J) = PARINI(J)                              
          IF (FIXPAR(J).GE.0) NRPFIX = NRPFIX + 1                               
   10 CONTINUE                                                                  
      METPAR(2) = ABS(METPAR(2))                                                
C                                                                               
C  Get the type of weighting from FITCHAR 13:13                                 
C                                                                               
      CALL STKRDC(WGTKEY,1,FZIWGT,1,K,WGTTYP,
     +    KUN,KNUL,ISTAT)
      CALL UPCAS(WGTTYP,WGTTYP)
      IF (WGTTYP(1:1).NE.'C' .AND. WGTTYP(1:1).NE.'W' .AND.                     
     +    WGTTYP(1:1).NE.'I' .AND. WGTTYP(1:1).NE.'S') CALL STETER(12,          
     +    'FIT/TABLE: Weight type mismatch')                                    
C                                                                               
      IF (FILTYP(1:1).EQ.'T') THEN                                              
C                                                                               
C  Get the selected columns for indep and dep variables                         
C                                                                               
          IF (WGTCOL.EQ.0) THEN                                                 
              ISTAR  = 2                                                        
              NRCOL  = NRIND + 1                                                
              ICOL(1) = 0                                                       
                                                                                
          ELSE                                                                  
              ISTAR  = 1                                                        
              NRCOL  = NRIND + 2                                                
              ICOL(1) = WGTCOL                                                  
          END IF                                                                
                                                                                
          ICOL(2) = DEPCOL                                                      
          NULL(1) = .FALSE.                                                     
          DO 20 I = 1,NRIND                                                     
              ICOL(I+2) = INDCOL(I)                                             
   20     CONTINUE                                                              
C                                                                               
C  Get number of selected rows in the table (nr of residuals)                   
C                                                                               
          NRDATA = 0                                                            
          CALL TBIGET(FZIDEN,NC,NROW,NS,NAC,NAR,STAT)                                 
          DO 40 IDAT = 1,NROW                                                   
              CALL TBSGET(FZIDEN,IDAT,ISEL,ISTAT)                              
              IF ( .NOT. ISEL) GO TO 40                                         
              CALL TBRRDR(FZIDEN,IDAT,NRCOL,ICOL(ISTAR),                  
     +                    XVAL(ISTAR),NULL,ISTAT)                              
              VALID  = ISEL                                                     
              DO 30 K = 1,NRCOL                                                 
                  VALID  = VALID .AND. ( .NOT. NULL(K))                         
   30         CONTINUE                                                          
              IF (VALID) NRDATA = NRDATA + 1                                    
   40     CONTINUE                                                              
                                                                                
      ELSE                                                                      
C                                                                               
C  Get the number of residuals for FIT/IMAG                                     
C                                                                               
          NRDATA = NRPIX(1)*NRPIX(2)*NRPIX(3)                                   
          NROW   = NRDATA                                                       
      END IF                                                                    
C                                                                               
C     Get qualifiers of the FIT command : PRINT,MAX. NR. Of F. EVAL,            
C     PRECISION,ETA and STEP. NOMPAR detects if method parameters               
C     have been given or are all to be defaulted.                               
C                                                                               
      NOMPAR = METPAR(1) .LT. 0.1                                               
      IF (ABS(ABS(METPAR(3))-999999.).LT.1.E-6) THEN                            
          METPAR(3) = 0.0                                                       
                                                                                
      ELSE                                                                      
          NOMPAR = .FALSE.                                                      
      END IF                                                                    
                                                                                
      IF (ABS(ABS(METPAR(2))-999999.).LT.1.E-6) THEN                            
          IF (ALGOR(1:5).EQ.'CGNND') THEN                                       
              METPAR(2) = 40*NRPARAM* (NRPARAM+5)                               
                                                                                
          ELSE IF (ALGOR(1:2).EQ.'QN') THEN                                     
              METPAR(2) = 100*NRPARAM                                           
                                                                                
          ELSE IF (ALGOR(1:3).EQ.'MGN') THEN                                    
              METPAR(2) = 50*NRPARAM                                            
          END IF                                                                
                                                                                
      ELSE                                                                      
          NOMPAR = .FALSE.                                                      
      END IF                                                                    
                                                                                
      IF (ABS(ABS(METPAR(4))-999999.).LT.1.E-6) THEN                            
          IF (NRPARAM.EQ.1) THEN                                                
              METPAR(4) = 0.                                                    
                                                                                
          ELSE                                                                  
              IF (ALGOR(1:5).EQ.'CGNND') THEN                                   
                  METPAR(4) = 0.5                                               
                                                                                
              ELSE IF (ALGOR(1:2).EQ.'QN') THEN                                 
                  METPAR(4) = 0.9                                               
                                                                                
              ELSE IF (ALGOR(1:3).EQ.'MGN') THEN                                
                  IF (NRPARAM.LT.10) THEN                                       
                      METPAR(4) = 0.5                                           
                                                                                
                  ELSE IF (NRPARAM.LE.20) THEN                                  
                      METPAR(4) = 0.1                                           
                                                                                
                  ELSE                                                          
                      METPAR(4) = 0.01                                          
                  END IF                                                        
                                                                                
              END IF                                                            
                                                                                
          END IF                                                                
                                                                                
      ELSE                                                                      
          NOMPAR = .FALSE.                                                      
      END IF                                                                    
                                                                                
      IF (ABS(ABS(METPAR(5))-999999.).LT.1.E-6) THEN                            
          METPAR(5) = 1.E05                                                     
                                                                                
      ELSE                                                                      
          NOMPAR = .FALSE.                                                      
      END IF                                                                    
C                                                                               
C   Get the bounds and store them in PARLOW and PARUPP                          
C                                                                               
      IF (BNDTYP(1:1).EQ.'G') THEN                                              
          IBOUND = 3                                                            
          CALL STKRDC('BNDTAB',1,1,60,K,BNDTAB,KUN,KNUL,ISTAT)                  
          IF (BNDTAB(1:1).EQ.'+') THEN
            CALL STTPUT('Bounds table not defined',ISTAT)
            CALL STSEPI
          ENDIF 
          CALL TBTOPN(BNDTAB,F_IO_MODE,TID,ISTAT)                             
          CALL TBIGET(TID,NC,NRBND,NS,NAC1,NAR1,ISTAT)                           
          IF (NC.LT.2) THEN                                                     
              CALL STETER(12,'FIT/... : Error in bounds table')                 
                                                                                
          ELSE IF (NC.GT.2) THEN                                                
              CALL STTPUT('*** WARN-FIT : Bounds table more than'//             
     +                    ' 2 columns, take first ones ***',ISTAT)
          END IF                                                                
                                                                                
          IF (NRBND.NE.1) THEN                                                  
              CALL STTPUT('*** WARN-FIT : Bounds table more than'//             
     +                    ' 1 row, global bounds taken as first one ***'        
     +                    ,ISTAT)                                              
          END IF                                                                
                                                                                
          CALL TBSGET(TID,1,ISEL,ISTAT)                                     
          ICOL(1) = 1                                                           
          ICOL(2) = 2                                                           
          CALL TBRRDR(TID,1,2,ICOL,                  
     +                    XVAL,NULL,ISTAT)                        
          VALID  = ISEL .AND. ((.NOT.NULL(1)) .OR. (.NOT.NULL(2)))              
          IF ( .NOT. VALID) CALL STETER(8,                                      
     +                           ' FIT/TABLE : Global bounds not founds'        
     +                                  )                                       
          IF ( .NOT. NULL(1)) THEN                                              
              PARLOW(1) = XVAL(1)                                               
                                                                                
          ELSE                                                                  
              PARLOW(1) = -1.D+6                                                
          END IF                                                                
                                                                                
          IF ( .NOT. NULL(2)) THEN                                              
              PARUPP(1) = XVAL(2)                                               
                                                                                
          ELSE                                                                  
              PARUPP(1) = 1.D+6                                                 
          END IF                                                                
                                                                                
          PARLOW(1) = XVAL(1)                                                   
          PARUPP(1) = XVAL(2)                                                   
          CALL STTPUT('   Lower Bound     Upper Bound',ISTAT)
          WRITE (LINE,9000) PARLOW(1),PARUPP(1)                                 
          CALL STTPUT(LINE,ISTAT)                                              
          CALL STTPUT(' ',ISTAT)                                               
                                                                                
      ELSE IF (BNDTYP(1:1).EQ.'I') THEN                                         
          IBOUND = 0                                                            
          CALL STKRDC('BNDTAB',1,1,60,K,BNDTAB,KUN,KNUL,ISTAT)   
          IF (BNDTAB(1:1).EQ.'+') THEN
            CALL STTPUT('Bounds table not defined',ISTAT)
            CALL STSEPI
          ENDIF 
          CALL TBTOPN(BNDTAB,F_IO_MODE,TID1,ISTAT)                          
          CALL TBIGET(TID1,NC,NRBND,NS,NAC1,NAR1,ISTAT)
          IF (NC.LT.2) THEN                                                     
              CALL STETER(12,'FIT/... : Error in bounds table')                 
                                                                                
          ELSE IF (NC.GT.2) THEN                                                
              CALL STTPUT('*** WARN-FIT : Bounds table more than'//             
     +                    ' 2 columns, take first ones ***',STATUS)             
          END IF                                                                
                                                                                
          IF (NRBND.GT.NRPARAM) THEN                                            
              CALL STTPUT('*** WARN-FIT : Bounds table has more '//             
     +                    'rows than parameters, take first ones ***',          
     +                    STATUS)                                               
          END IF                                                                
                                                                                
          DO 50 IDAT = 1,MIN0(NRBND,NRPARAM)                                    
              CALL TBSGET(TID1,IDAT,ISEL,STATUS)                              
              IF (ISEL) THEN                                                    
                  ICOL(1) = 1                                                   
                  ICOL(2) = 2                                                   
                  CALL TBRRDR(TID1,IDAT,2,ICOL,                  
     +                    XVAL,NULL,ISTAT)                        
               IF ( .NOT. NULL(1)) THEN                                      
                      PARLOW(IDAT) = XVAL(1)                                    
                                                                                
                  ELSE                                                          
                      PARLOW(IDAT) = -1.D+6                                     
                  END IF                                                        
                                                                                
                  IF ( .NOT. NULL(2)) THEN                                      
                      PARUPP(IDAT) = XVAL(2)                                    
                                                                                
                  ELSE                                                          
                      PARUPP(IDAT) = 1.D+6                                      
                  END IF                                                        
                                                                                
              ELSE                                                              
                  PARLOW(IDAT) = -1.D+6                                         
                  PARUPP(IDAT) = 1.D+6                                          
              END IF                                                            
                                                                                
   50     CONTINUE                                                              
          IF (NRBND.LT.NRPARAM) THEN                                            
              DO 60 IDAT = NRBND + 1,NRPARAM                                    
                  PARLOW(IDAT) = -1.D+6                                         
                  PARUPP(IDAT) = 1.D+6                                          
   60         CONTINUE                                                          
          END IF                                                                
                                                                                
          CALL STTPUT(                                                          
     +             '   Initial values     Lower Bounds     Upper Bounds'        
     +                ,STATUS)                                                  
          DO 70 K = 1,NRPARAM                                                   
              WRITE (LINE,9010) PARAM(K),PARLOW(K),PARUPP(K)                    
              CALL STTPUT(LINE,STATUS)                                          
   70     CONTINUE                                                              
          CALL STTPUT(' ',STATUS)                                               
                                                                                
      ELSE IF (BNDTYP(1:1).EQ.'P') THEN                                         
          CALL STTPUT(' Parameters must be positive',STATUS)                    
          IBOUND = 2                                                            
      END IF                                                                    
C                                                                               
C  Call to the NAG routines after having displayed the parameters and           
C  get the necessary V.M. space.                                                
C                                                                               
      IFAIL1 = 1                                                                
      XVAL(1) = 1.                                                              
      IF (ALGOR(1:5).EQ.'CGNND') THEN                                           
          CALL STTPUT(' ',STATUS)                                               
          CALL STTPUT(                                                          
     +          ' Constrained Non-linear Least Squares Fitting : Corr. '        
     +                //'Gauss-Newton M. (no der)',STATUS)                      
          CALL STTPUT(                                                          
     +          ' _____________________________________________________'        
     +                //'________________________',STATUS)                      
          CALL STTPUT(' ',STATUS)                                               
          CALL STTPUT(' ',STATUS)                                               
          CALL STTPUT(                                                          
     +  ' Print   Max. Fct. eval.  Prec. on Param.  Lin. Search  Domain'        
     +                //'  Weight',STATUS)                                      
          WRITE (LINE,9020) NINT(METPAR(1)),NINT(METPAR(2)),METPAR(3),          
     +      METPAR(4),METPAR(5),WGTTYP                                          
          CALL STTPUT(LINE,STATUS)                                              
          CALL STTPUT(' ',STATUS)                                               
                                                                                
          IF (NOMPAR) THEN                                                      
              ROUTIN = 1                                                        
              IF (NRPARAM.EQ.1) THEN                                            
                  LW     = 13                                                   
                                                                                
              ELSE                                                              
                  LW     = 12*NRPARAM + NRPARAM* (NRPARAM-1)/2                  
              END IF                                                            
                                                                                
              LIW    = NRPARAM + 2                                              
                                                                                
              NRVM   = 2*LW + LIW                                               
              NBYT = 4*NRVM
              CALL TDMGET(NBYT,PTRW,STATUS)
              PTRIW  = PTRW + 2*LW                                              
                                                                                
              CALL E04JAF(NRPARAM,IBOUND,PARLOW,PARUPP,PARAM,FSUMSQ,            
     +                    MADRID(PTRIW),LIW,MADRID(PTRW),LW,IFAIL1)             
                                                                                
          ELSE                                                                  
              ROUTIN = 2                                                        
              LW     = 9*NRPARAM                                                
              LHSL   = MAX0(NRPARAM* (NRPARAM-1)/2,1)                           
              LIW    = 10                                                       
                                                                                
              NRVM   = 2* (LW+LHSL+4*NRPARAM) + LIW                             

C             CALL STIPUT('WORK ',10,9,1,1,NRVM,1.,XVAL,' ',' ',PTRW,            
C    +                    STATUS)                                               
              NBYT = 4*NRVM
              CALL TDMGET(NBYT,PTRW,STATUS)
              PTRDLT = PTRW + 2*LW                                              
              PTRHSL = PTRDLT + 2*NRPARAM                                       
              PTRHSD = PTRHSL + 2*LHSL                                          
              PTRSTA = PTRHSD + 2*NRPARAM                                       
              PTRGRD = PTRSTA + NRPARAM                                       
              PTRIW  = PTRGRD + 2*NRPARAM                                       
                                                                                
              IFAIL3 = 1                                                        
              CALL E04HBF(NRPARAM,FUNCT3,PARAM,NRFEVA,MADRID(PTRDLT),          
     +                    MADRID(PTRHSL),LHSL,MADRID(PTRHSD),FSUMSQ,            
     +                    MADRID(PTRGRD),MADRID(PTRIW),LIW,MADRID(PTRW),        
     +                    LW,IFAIL3)                                            
                                                                                
                                                                                
              CALL MNMX(MADRID(PTRHSD),NRPARAM,HMIN,HMAX,IMIN,IMAX)             
              IF (HMAX.LT.1.D+4*HMIN) THEN                                      
                  WRITE (LINE,9030) HMIN,HMAX                                   
                  CALL STETER(21,LINE)                                          
              END IF                                                            
                                                                                
              CALL E04JBF(NRPARAM,FUNCT3,MINMON,NINT(METPAR(1)),.TRUE.,         
     +                    1,E04JBQ,NINT(METPAR(2)),DBLE(METPAR(4)),             
     +                    DBLE(METPAR(3)),DBLE(METPAR(5)),0.D0,                 
     +                    MADRID(PTRDLT),IBOUND,PARLOW,PARUPP,PARAM,            
     +                    MADRID(PTRHSL),LHSL,MADRID(PTRHSD),                   
     +                    MADRID(PTRSTA),FSUMSQ,MADRID(PTRGRD),                 
     +                    MADRID(PTRIW),LIW,MADRID(PTRW),LW,IFAIL1)             
                                                                                
          END IF                                                                
                                                                                
      ELSE IF (ALGOR(1:2).EQ.'QN') THEN                                         
          CALL STTPUT(' ',STATUS)                                               
          CALL STTPUT(                                                          
     +   ' Constrained Non-linear Least Squares Fitting : Quasi-Newton '        
     +                //'Method',STATUS)                                        
          CALL STTPUT(                                                          
     +   ' ____________________________________________________________'        
     +                //'______',STATUS)                                        
          CALL STTPUT(' ',STATUS)                                               
          CALL STTPUT(' ',STATUS)                                               
          CALL STTPUT(                                                          
     +  ' Print   Max. Fct. eval.  Prec. on Param.  Lin. Search  Domain'        
     +                //'  Weight',STATUS)                                      
          WRITE (LINE,9020) NINT(METPAR(1)),NINT(METPAR(2)),METPAR(3),          
     +      METPAR(4),METPAR(5),WGTTYP                                          
          CALL STTPUT(LINE,STATUS)                                              
          CALL STTPUT(' ',STATUS)                                               
                                                                                
          IF (NOMPAR) THEN                                                      
              ROUTIN = 3                                                        
              IF (NRPARAM.EQ.1) THEN                                            
                  LW     = 11                                                   
                                                                                
              ELSE                                                              
                  LW     = 10*NRPARAM + NRPARAM* (NRPARAM-1)/2                  
              END IF                                                            
                                                                                
              LIW    = NRPARAM                                              
              NRVM   = LW*2 + LIW+ 2*NRPARAM                                               
C              CALL STIPUT('WORK ',10,9,1,1,NRVM,1.,XVAL,' ',' ',PTRW,            
C     +                    STATUS)                                               
              NBYT=NRVM*4
              CALL TDMGET(NBYT,PTRW,STATUS)
              PTRIW  = PTRW + 2*LW                                              
              PTRGRD = PTRIW + 2*NRPARAM
              CALL E04KAF(NRPARAM,IBOUND,PARLOW,PARUPP,PARAM,FSUMSQ,            
     +                    MADRID(PTRGRD),MADRID(PTRIW),LIW,MADRID(PTRW),        
     +                    LW,IFAIL1)                                            
                                                                                
          ELSE                                                                  
              ROUTIN = 4                                                        
              LW     = 9*NRPARAM                                                
              LIW    = 10                                                       
              LHSL   = MAX(NRPARAM* (NRPARAM-1)/2,1)                            
                                                                                
              NRVM   = (LW+3*NRPARAM+LHSL+1)*2 + LIW + NRPARAM                  
              NBYT = NRVM*4
              CALL TDMGET(NBYT,PTRW,STATUS)
C              CALL STIPUT('WORK ',10,9,1,1,NRVM,1.,XVAL,' ',' ',PTRW,            
C     +                    STATUS)                                               
              PTRHSL = PTRW + 3*NRPARAM                                         
              PTRHSD = PTRHSL + 2*LHSL                                          
              PTRSTA = PTRHSD + 2*NRPARAM                                       
              PTRGRD = PTRSTA + NRPARAM                                       
              PTRIW  = PTRGRD + 2*NRPARAM                                       
                                                                                
                                                                                
              IFAIL3 = 1                                                        
              CALL E04HCF(NRPARAM,FUNCT4,PARINI,FSUMSQ,MADRID(PTRGRD),          
     +                    MADRID(PTRIW),LIW,MADRID(PTRW),LW,IFAIL3)             
                                                                                
              IF (IFAIL3.EQ.2) THEN                                             
                  CALL MINMSG(ALGOR(1:5),ROUTIN,10)                             
                  GO TO 100                                                     
                                                                                
              END IF                                                            
                                                                                
              CALL SETZER(MADRID(PTRHSL),LHSL)                                 
              CALL SETONE(MADRID(PTRHSD),NRPARAM)                               
                                                                                
              CALL E04KBF(NRPARAM,FUNCT4,MINMON,NINT(METPAR(1)),.TRUE.,         
     +                    1,E04LBS,NINT(METPAR(2)),DBLE(METPAR(4)),             
     +                    DBLE(METPAR(3)),DBLE(METPAR(5)),0.D0,IBOUND,          
     +                    PARLOW,PARUPP,PARAM,MADRID(PTRHSL),LHSL,              
     +                    MADRID(PTRHSD),MADRID(PTRSTA),FSUMSQ,                 
     +                    MADRID(PTRGRD),MADRID(PTRIW),LIW,MADRID(PTRW),        
     +                    LW,IFAIL1)                                            
                                                                                
          END IF                                                                
                                                                                
      ELSE IF (ALGOR(1:3).EQ.'MGN') THEN                                        
          CALL STTPUT(' ',STATUS)                                               
          CALL STTPUT(                                                          
     +       ' Constrained Non-linear Least Squares Fitting : Modified '        
     +                //'Gauss-Newton Method',STATUS)                           
          CALL STTPUT(                                                          
     +       ' ________________________________________________________'        
     +                //'___________________',STATUS)                           
          CALL STTPUT(' ',STATUS)                                               
          CALL STTPUT(' ',STATUS)                                               
          CALL STTPUT(                                                          
     +  ' Print   Max. Fct. eval.  Prec. on Param.  Lin. Search  Domain'        
     +                //'  Weight',STATUS)                                      
          WRITE (LINE,9020) NINT(METPAR(1)),NINT(METPAR(2)),METPAR(3),          
     +      METPAR(4),METPAR(5),WGTTYP                                          
          CALL STTPUT(LINE,STATUS)                                              
          CALL STTPUT(' ',STATUS)                                               
                                                                                
          IF (NOMPAR) THEN                                                      
              ROUTIN = 5                                                        
              IF (NRPARAM.EQ.1) THEN                                            
                  LW     = 11                                                   
                                                                                
              ELSE                                                              
                  LW     = 10*NRPARAM + NRPARAM* (NRPARAM-1)/2                  
              END IF                                                            
                                                                                
              LIW    = NRPARAM + 4                                              
              NRVM   = LW*2 + LIW                                               
C             CALL STIPUT('WORK ',10,9,1,1,NRVM,1.,XVAL,' ',' ',PTRW,            
C     +                    STATUS)                  
              NBYT = NRVM*4
              CALL TDMGET(NBYT,PTRW,STATUS)                             
              PTRIW  = PTRW + 2*LW                                              
              PTRGRD = PTRIW + 2 + NRPARAM
              CALL E04KCF(NRPARAM,IBOUND,PARLOW,PARUPP,PARAM,FSUMSQ,            
     +                    MADRID(PTRGRD),MADRID(PTRIW),LIW,MADRID(PTRW),        
     +                    LW,IFAIL1)                                            
                                                                                
          ELSE                                                                  
              ROUTIN = 6                                                        
              IF (NRPARAM.EQ.1) THEN                                            
                  LW     = 8                                                    
                                                                                
              ELSE                                                              
                  LW     = 7*NRPARAM + NRPARAM* (NRPARAM-1)/2                   
              END IF                                                            
                                                                                
              LIW    = 10                                                       
              LHSL   = MAX(NRPARAM* (NRPARAM-1)/2,1)                            
                                                                                
              NRVM   = (LW+3*NRPARAM+LHSL+1)*2 + LIW + NRPARAM                  
              NBYT = NRVM*4
              CALL TDMGET(NBYT,PTRW,STATUS)
C              CALL STIPUT('WORK ',10,9,1,1,NRVM,1.,XVAL,' ',' ',PTRW,            
C     +                    STATUS)                                               
              PTRHSL = PTRW + 2*LW                                              
              PTRHSD = PTRHSL + 2*LHSL                                          
              PTRSTA = PTRHSD + 2*NRPARAM                                       
              PTRGRD = PTRSTA + NRPARAM                                       
              PTRIW  = PTRGRD + 2*NRPARAM                                       
                                                                                
                                                                                
              IFAIL3 = 1                                                        
              CALL E04HCF(NRPARAM,FUNCT4,PARINI,FSUMSQ,MADRID(PTRGRD),          
     +                    MADRID(PTRIW),LIW,MADRID(PTRW),LW,IFAIL3)             
                                                                                
              IF (IFAIL3.EQ.2) THEN                                             
                  CALL MINMSG(ALGOR(1:5),ROUTIN,10)                             
                  GO TO 100                                                     
                                                                                
              END IF                                                            
                                                                                
              CALL E04KDF(NRPARAM,FUNCT4,MINMON,NINT(METPAR(1)),                
     +                    NINT(METPAR(2)),DBLE(METPAR(4)),                      
     +                    DBLE(METPAR(3)),0.D0,DBLE(METPAR(5)),IBOUND,          
     +                    PARLOW,PARUPP,PARAM,MADRID(PTRHSL),LHSL,              
     +                    MADRID(PTRHSD),MADRID(PTRSTA),FSUMSQ,                 
     +                    MADRID(PTRGRD),MADRID(PTRIW),LIW,MADRID(PTRW),        
     +                    LW,IFAIL1)                                            
                                                                                
          END IF                                                                
                                                                                
      END IF                                                                    
                                                                                
      CALL MINMSG(ALGOR(1:5),ROUTIN,IFAIL1)                                     
      CALL STTPUT(' ',STATUS)                                                   
                                                                                
      IF (IFAIL1.NE.1 .AND. IFAIL1.NE.4 .AND. IFAIL1.NE.9) THEN                 
          IF (NOMPAR) THEN                                                      
              IF (ALGOR(1:5).EQ.'CGNND') THEN                                   
                  PTRHSD = PTRW + 8* (6*NRPARAM+2*NRDATA+NRDATA*NRPARAM+        
     +                     MAX0(1, (NRPARAM* (NRPARAM-1)/2)))                   
                  PTRHSD = PTRHSD + 8*NRPARAM                                   
                                                                                
              ELSE                                                              
                  PTRHSD = PTRW + 8* (7*NRPARAM+2*NRDATA+NRDATA*NRPARAM+        
     +                     NRPARAM* (NRPARAM+1)/2+                              
     +                     MAX0(1, (NRPARAM* (NRPARAM-1)/2)))                   
                  PTRHSD = PTRHSD + 8*NRPARAM                                   
                                                                                
              END IF                                                            
                                                                                
          END IF                                                                
                                                                                
                                                                                
      ELSE                                                                      
          CALL STTPUT(                                                          
     +        ' *** Unable to evaluate the errors on the parameters ***'        
     +                ,STATUS)                                                  
      END IF                                                                    
                                                                                
  100 CONTINUE                                                                  
C                                                                               
C  Pass the final reduced chi squares to the COMMON block                       
C                                                                               
      NRITER = NRFEVA                                                          
      FINCHI = FSUMSQ/DBLE(NRDATA-NRPARAM+NRPFIX)                             
C                                                                               
C  Exit                                                                         
C                                                                               
      RETURN                                                                    
                                                                                
 9000 FORMAT (5X,G12.5,10X,G12.5)                                               
 9010 FORMAT (5X,G12.5,10X,G12.5,10X,G12.5)                                     
 9020 FORMAT (I6,8X,I7,7X,1PE9.1,8X,0PF4.2,6X,0PF8.0,5X,A1)                     
 9030 FORMAT ('FIT/TABLE : Bad Scaling of the problem ','HSD from ',            
     +       1PD9.2,' to ',1PD9.2)                                              
 9040 FORMAT ('*** Detected only ',I3,' linear independant parameters ',        
     +       'on the total of ',I3,' ***')                                      
      END                                                                       
