        SUBROUTINE MCALC(IPR,IPD,IWF,ITD,ILV,IWE,IPWC,IWTC,WX,WY,XIX,
     1XIX2,XIY,XIY2,USX,USY,PXS,PYS,EPS1,EPS2,EPS3,DIVI,IRET,MAXG,
     2SUMSXR,SUMSYR,S2,AA,SAK,SAKR,RC,SDER,IDIRC)
	IMPLICIT REAL*8 (A-H,O-Z)
	EXTERNAL SFUNC
	CHARACTER*78 ALPHA
	DIMENSION C(10,10),V(10),FA(10),AA(10),SAK(10),RC(10,10)
	DIMENSION FX(151),SAKR(10),RX(151),RY(151),XKEEP(151)
	DIMENSION XHOLD(151),SDER(10),
     2  SPOX(10),SPOY(10),SMOX(10),SMOY(10),SMMX(10,10),
     3  SMMY(10,10),SPPX(10,10),SPPY(10,10),CINV(10,10),YN(151),
     4  CO(10,10),
     5  WX(151),WY(151),AK(10),CV(10,10)
 	PARAMETER (DM9=1.D-9)
C
C	COMMON /MORIN/ IPR,IPR1,ITYP,INX
	COMMON /IN/ N,NF,NP,MAXD,G,EPS,A(10),NOUT,NS(10),ALPHA
	COMMON /DATA/ Y(151),X(151),SX(151),SY(151)
	COMMON /EXTRA/ FAC
	COMMON /SECDER/ P1,P2,DEL,XX
	COMMON /EXTRY/ DEL1,DEL2,JJJ,KKK,XPART
	COMMON /DON/ RX,RY,XKEEP,YPART,XPORT,XHOLD
	COMMON /WTI/ WXI(151),WYI(151)
	SAVE
C
C
C	NO DEMING STEPS IF IWE = 1
	IF(IWE.EQ.1) THEN
		DEL = 1.D0
		DX = 1.D0
		GOTO 410
	ENDIF
	ISD = 0
	SDOLD = 1.D0
C
	IF(IPR.NE.0.AND.IWF.EQ.0) PRINT 180
180	FORMAT (/6X, '****** DEMING STEPS ******')
	ID = 1
680	FORMAT(/4X,19H *-*-*-*-*-*-*-*-* )
	IFIN = 0
	IF(ITD.EQ.0) THEN
        	DO 235 K = 1,NP
	   		AA(K) = 0.D0
235		CONTINUE
	ENDIF
C
200	CONTINUE
	SUMDF = 0.D0
	IF(IPR.NE.0) PRINT 220,ITD
220 	FORMAT(/5X,17HITERATION NUMBER I3/)
C
C	CLEAR C AND V
C
	DO 250 J = 1,NP
        DO 240 K = 1,NP
240      C(J,K) = 0.D0
250	V(J) = 0.D0
C
	DO 260 J = 1,NP
260	SDER(J) = 0.D0
	SB = 0.D0
	SUMSY = 0.D0
	SUMSX = 0.D0
C
	DO 380 I = 1,N
		XX = X(I)
C
C	CALCULATE DERIVATIVES W.R.T. X
C
		CALL DERXE(X,DX,I)
		DX = -DX
C
C	CALCULATE L(I)     (CALLED DL)
C	CALCULATE DERIVATIVES W.R.T. THE PARAMETERS
C
		DL = WY(I)
		DXS = DX*DX
			DL = DL + WX(I)*DXS
		DLI = 1.D0/DL
C
C	ESTIMATED PARTIALS
C
		CALL DERAE(FA)
C
		IF(NOUT.EQ.0) GOTO 300
		DO 290 J = 1,NOUT
			NSJ = NS(J)
290		FA(NSJ) = 0.D0
300		CONTINUE
C
C	CALCULATE C MATRIX
C
		DO 310 J = 1,NP
310		FA(J) = -FA(J)
		DO 320 J = 1,NP
		DO 320 K =J,NP
		C(J,K) = C(J,K) + FA(J)*FA(K)*DLI
320		CONTINUE
		IF(NOUT.EQ.0) GOTO 343
		DO 330 J = 1,NOUT
			K = NS(J)
330		C(K,K) = DLI
343		CONTINUE
C
C	CALCULATE V VECTOR
C
		YCALC = SFUNC(ISD)
		YY = Y(I)
		DF = YY - YCALC
		DFL = DF*DLI
		DO 350 K = 1,NP
			V(K) = V(K) + FA(K)*DFL
350		CONTINUE
		EXTEM = -2.D0*DFL
		DO 360 J = 1,NP
360		SDER(J) = SDER(J) + EXTEM*FA(J)
C
		SUM = 0.D0
		DO 370 J = 1,NP
370		SUM = SUM + FA(J)*AA(J)
		FLAM = (DF - SUM)*DLI
		FLAMS = FLAM*FLAM
		SUMSY = SUMSY + FLAMS*WY(I)
		SUMSX = SUMSX + DXS*FLAMS*WX(I)
380	CONTINUE
C	END OF I LOOP
C
	IF(ITD.EQ.0) THEN
		SOL = 0
	ELSE
		SOL = S2
	ENDIF
C
	S2 = SUMSY + SUMSX
	SB = S2
	SSIN = 1.D0/SB
C
C
	DO 390 J = 1,NP
390	SDER(J) = A(J)*SDER(J)*SSIN
	IF(IPR.NE.0) PRINT 400, (SDER(J),J=1,NP)
400	FORMAT(4X,'(A/S)*(DS/DA):' 1P5E11.2)
C
	IF(IWE.EQ.0) GOTO 580
C
C	%%%%%%    GENERAL METHOD PART    %%%%%%
C
410	CONTINUE
	IF(IPR.NE.0.AND.IWF.EQ.0) PRINT 420
420	FORMAT(/14X,'GENERAL METHOD STEPS:')
	SCOMP = 1.D25
	IDIRC = 0	
430	CONTINUE
	IF(IPR.NE.0) PRINT 220,IDIRC
440	CONTINUE
C
C	CALCULATE S PARTIAL DERIVATIVES
C
	IF(IPD.EQ.1) THEN
		WRITE(*,638) IDIRC,SB
		WRITE(*,*) 'LOCAL PARAMETER VALUES'
		WRITE(*,191) (A(K),K=1,NP)
	ENDIF
191	FORMAT(5X,1P5D13.4)
638	FORMAT(10X,'IT. # =',I5,5X,'S =',1PD14.5)
	JJJ = 1
	KKK = 1
	DEL1 = 0.D0
	DEL2 = 0.D0
	CALL SCALC(0,0)
	S00 = XPORT + YPART
	SUMSY = YPART
	SUMSX = XPART
	NPM1 = NP - 1
	DO 450 I = 1,N
450	FX(I) = XHOLD(I)
	DO 480 J = 1,NP
		DEL1 = DABS(FAC*A(J) + 1.D-50)
		JJJ = J
		CALL SCALC(1,0)
		SPOX(J) = XPART
		SPOY(J) = YPART
		CALL SCALC(-1,0)
		SMOX(J) = XPART
		SMOY(J) = YPART
		JP1 = J + 1
		IF(JP1.GT.NP) GOTO 473
		DO 460 K = JP1,NP
			DEL2 = DABS(FAC*A(K) + 1.D-50)
			KKK = K
			CALL SCALC(1,1)
			SPPX(J,K) = XPART
			SPPY(J,K) = YPART
			CALL SCALC(-1,-1)
			SMMX(J,K) = XPART
			SMMY(J,K) = YPART
460		CONTINUE
473		CONTINUE
480	CONTINUE
	DO 510 J = 1,NP
		DEL1 = DABS(FAC*A(J) + 1.D-50)
             V(J) = (SPOY(J) - SMOY(J) + SPOX(J) - SMOX(J))*0.5D0/DEL1
	     C(J,J) = (SPOY(J) + SMOY(J) -2.D0*SUMSY + SPOX(J) +
     1       SMOX(J) - 2.D0*SUMSX)/DEL1**2
		JP1 = J + 1
		IF(JP1.GT.NP) GOTO 500
		DO 490 K = JP1,NP
			DEL2 = DABS(FAC*A(K) + 1.D-50)
			P1 = SPOY(J) + SMOY(J) + SPOY(K) + SMOY(K) -
     1                  2.D0*SUMSY - SPPY(J,K) - SMMY(J,K)	
			P2 = SPOX(J) + SMOX(J) + SPOX(K) + SMOX(K) -
     1                  2.D0*SUMSX - SPPX(J,K) - SMMX(J,K)			
			C(J,K) = -(P1 + P2)*0.5D0/(DEL1*DEL2)
			C(K,J) = C(J,K)
490		CONTINUE
500		CONTINUE
510	CONTINUE
	IF(NOUT.EQ.0) GOTO 540
	DO 530 I = 1,NOUT
		K = NS(I)
		V(K) = 0.D0
		DO 520 II = 1,NP
			C(K,II) = 0.D0
			C(II,K) = 0.D0
520		CONTINUE
		C(K,K) = 1.D0
530	CONTINUE
540	CONTINUE
		SOL = S2
		SB = S00
		S2 = SB
		SCOMP = SB
		SSIN = 1.D0/SB
	SMDR = 0.D0
	DO 550 K = 1,NP
	        SDERO = SDER(K)
		SDER(K) = A(K)*V(K)*SSIN
		IF(SDER(K).NE.0) THEN
			SDR = (SDER(K) + SDERO)/SDER(K)
			SMDR = SMDR + SDR 
		ENDIF
C
550	CONTINUE
	IF(DABS(SMDR).LT.1.D-5) GOTO 850
	IF(IPR.NE.0) PRINT 400, (SDER(K),K=1,NP)
	DO 560 I = 1,N
560	X(I) = FX(I)
C
C
580	CONTINUE
	SDIMP = DABS((SB - SDOLD)/SDOLD)
	SDOLD = SB
	IF(IPR.NE.0) PRINT 650, (A(K),K=1,NP)
650	FORMAT(/5X,16HPARAMETER VALUES /(2X,1P4D17.9))
C
	SUMSYR = SUMSY
	SUMSXR = SUMSX
C
	STST = (S2 - SOL)/S2
	IF(IPR.NE.0) PRINT 590, SUMSXR,SUMSYR,S2
590	FORMAT(/5X,6HS(X) = 1PD17.7/ 
     1    5X,6HS(Y) = D17.7/5X,14HS(Y) + S(X) = D20.11/)
C
	IF(ILV.EQ.0.OR.IWE.EQ.1) THEN
		ALLV = 0.D0
		IDV = 1
		GOTO 663
	ELSE
	  	IF(ITD.LE.1) THEN
			ALLV = 0
			IDV = 1
		ELSEIF(ITD.EQ.2) THEN
			ALLV = 1.D-3
			IDV = 1
		ELSE
			IF(STST.GT.0) THEN
				ALLV = 1.0D1*ALLV
				ITST = 0
				IDV = 0
			ELSE
				ALLV = 0.1D0*ALLV
				IDV = 1
				IF(-STST.LT.0.1D0) ITST = ITST + 1
			ENDIF
		ENDIF
		IF(ITST.GE.1.AND.IWE.EQ.0) THEN
			ALLV = 0.D0
			IDV = 1
			ITD = MAXD
			GOTO 663
		ELSEIF(ITST.GE.2) THEN
			ALLV = 0.D0
			IDV = 1
			GOTO 663
		ENDIF
	ENDIF
663	CONTINUE
C
C	NOW FILL OUT C MATRIX
C
	DO 610 J = 1,NP
		IF(ITD.GT.0.AND.IDV.EQ.0) THEN
			A(J) = AK(J)
		ENDIF
		DO 609 K = J,NP
 			C(K,J) = C(J,K)
			IF(ITD.LE.1.OR.IDV.EQ.1) THEN
				CV(J,K) = C(J,K)
				CV(K,J) = C(J,K)
			ENDIF
			IF(IDV.EQ.0) THEN
				C(J,K) = CV(J,K)
				C(K,J) = CV(J,K)
			ENDIF
C
		   IF(J.EQ.K.AND.ILV.EQ.1) C(K,J) = C(K,J)*(1.D0 + ALLV) 
609		CONTINUE
610	CONTINUE
C
C	INVERT C MATRIX
C
	IF(NP.GT.2) GOTO 630
	IF(NP.EQ.2) GOTO 620
	CINV(1,1) = 1.D0/C(1,1)
	GOTO 640
620	DETR = C(1,1)*C(2,2) - C(1,2)**2
	DETRA = DABS(DETR)
	IF(DETRA.LE.1.D-40.OR.DETRA.GT.1.D90) THEN
		WRITE(*,864)
864    FORMAT(5X,'SINGULAR MATRIX !!: TRY DIFFERENT PARAMETER CHOICES'/)
		IRET = 1
		GOTO 518
	ENDIF
	DETRI = 1.D0/DETR
	CINV(1,1) = C(2,2)*DETRI
	CINV(2,2) = C(1,1)*DETRI
	CINV(2,1) = -C(2,1)*DETRI
	CINV(1,2) = CINV(2,1)
	GOTO 640
C
630	CONTINUE
	DO 134 I = 1,NP
	DO 135 J = 1,NP
		CO(I,J) = C(I,J)
135	CONTINUE
134	CONTINUE
	CALL DMATIN(CO,NP,CINV,ISD)
	IF(ISD.EQ.1) THEN
		WRITE(*,*) 'BAD MATRIX: TRY NEW PARAMETER CHOICES'
		ISD = 0
		PRINT 70
70	        FORMAT(/)
		IRET = 1
		GOTO 518
	ENDIF
C
640	CONTINUE
C
	DO 668 K = 1,NP
		AA(K) = 0.D0
		DO 660 J = 1,NP
			TERMCV = CINV(J,K)*V(J)*G
			AA(K) = AA(K) + TERMCV
660		CONTINUE
668	CONTINUE
	IF(IPR.NE.0) THEN
		PRINT 670, (AA(K),K=1,NP)
670		FORMAT(5X,10HINCREMENTS  /(2X,1P4D17.5))
		PRINT 680
	ENDIF
C
	IF(IWE.EQ.0) GOTO 750
C
C
C	DO GENERAL WEIGHTING SOLUTION BELOW
	IF(IWTC.EQ.5) THEN
		IPWC = IPWC + 1
		IF(IPWC.EQ.1.AND.IDIRC.EQ.0) THEN
			WRITE(*,751)
751		FORMAT(/10X,'** ITERATIVE RE-WEIGHTING SOLUTION **')
			WRITE(*,851)
851             FORMAT(/3X,'REL. CHANGE IN S',13X,'S',14X,'RWT-IT  MAIN-
     1IT')
853		        FORMAT(2X,1PD14.4,7X,1PD20.12,3X,I5,3X,I5)
		ENDIF
          IF(SDIMP.GT.3.D-11.OR.IPWC.LT.2) THEN
C
		DO 722 K = 1,NP
			A(K) = A(K) - AA(K)
			IF(A(K).EQ.0.D0) GOTO 722
722		CONTINUE
		DO 132 I = 1,N
		   IF(XIX.NE.0) THEN
			WX(I) = PXS*(DABS(X(I))**XIX2 + USX + 1.D-20)
		   ELSE
			WX(I) = PXS + USX
		   ENDIF
		   IF(XIY.NE.0) THEN
			XX = X(I)
			YN(I) = SFUNC(ISD)
 			WY(I) = PYS*(DABS(YN(I))**XIY2 + USY + 1.D-20)
		   ELSE
			WY(I) = PYS + USY 
		   ENDIF
		WXI(I) = 1.D0/WX(I)
		WYI(I) = 1.D0/WY(I)
		SX(I) = DSQRT(WX(I))
		SY(I) = DSQRT(WY(I))
132		CONTINUE
C
		IF(IPWC.GT.MAXG) THEN
			WRITE(*,157) IPWC
157			FORMAT(/5X,'!  ITERATIVE RE-WEIGHING DID NOT CON
     1VERGE IN',I3,2X,'ITERATIONS !'/)
			GOTO 850
		ENDIF
		GOTO 440
          ENDIF
		WRITE(*,853) SDIMP,SB,IDIRC,IPWC
	ENDIF
C
	DERSUM = 0.D0
	DO 700 K = 1,NP
		DERSUM = DERSUM + DABS(SDER(K))
700	CONTINUE
	IF(DERSUM.GT.EPS2) GOTO 710
	IF(SDIMP.LE.EPS3.AND.IDIRC.GT.1) GOTO 850
710	CONTINUE
	IQUIT = 0
	DO 720 K = 1,NP
		A(K) = A(K) - AA(K)
		IF(A(K).EQ.0.D0) GOTO 720
		IF(DABS(AA(K)/A(K)).GT.EPS1) IQUIT = 1
720	CONTINUE
	IF(IQUIT.EQ.0.AND.IDIRC.GT.1) GOTO 850
C
	IDIRC = IDIRC + 1
	IPWC = 0
	IF(IDIRC.LE.MAXG) GOTO 430
	IF(IPR.NE.0) PRINT 840
840	FORMAT(//5X,43HITERATION COUNT HAS REACHED MAXIMUM ALLOWED /
     1   5X, 34HAUTOMATIC BRANCH TO FINAL PRINTOUT //)
	GOTO 850
C
750	CONTINUE
C
C	HAS DEMING CONVERGENCE BEEN ATTAINED?
C
	DPARS = 0.D0
	DO 760 K = 1,NP
		DPARS = DPARS + DABS(AA(K))
		IF(A(K).EQ.0.D0) GOTO 760
		IF(DABS(AA(K)/A(K)).GT.EPS) GO TO 770
760	CONTINUE
	IF(DPARS.GT.EPS) GOTO 770
	GOTO 810
C
C	ADJUST PARAMETERS
C
770	CONTINUE
	DO 780 K = 1,NP
	   IF(ILV.EQ.0) GOTO 773
		IF(ITD.EQ.0.OR.IDV.EQ.1) THEN
			AK(K) = A(K)
		ELSE
			A(K) = AK(K)
		ENDIF
773		A(K) = A(K) - AA(K)
C
780	CONTINUE
	IF(DABS(STST).LT.EPS.AND.ITD.GT.2) GOTO 810
	SCOMP = SB
	ITD = ITD + 1
	IF(ITD.LE.MAXD) GOTO 200
C
C	DEMING CONVERGENCE REACHED
C
810	CONTINUE
	IF(ILV.EQ.1) THEN
		ALLV = 1.D-4
		ITST = 0
	ENDIF
	IF(IPR.NE.0) PRINT 820
820	FORMAT(/10X,'****** END OF DEMING ITERATIONS ******'/)
C
	IF(IPR.EQ.0.AND.IWF.EQ.0) THEN
		PRINT 825
825	FORMAT(//10X,'******  DEMING SOLUTION ******')
		PRINT 220, ITD
	CALL PRINTOUT(SUMSXR,SUMSYR,S2,A,AA,SAK,SAKR,RC,SDER,IMEND,NP)
	ENDIF
C
	IWE = 1
	CALL XCHNGE(WXI,WYI)
	DO 830 K = 1,N
830	X(K) = XHOLD(K)
	GOTO 410
C
850	CONTINUE
C
C	CALCULATE STATISTICAL MEASURES
C
	TERM = 2.D0*SB*DIVI
	DO 860 I = 1,NP
		SAK(I) = DSQRT(DABS(TERM*CINV(I,I)))
860		IF(A(I).NE.0.D0) SAKR(I) = SAK(I)/DABS(A(I))

	DO 870 I = 1,NP
	DO 870 J = 1,NP
		RC(I,J) = (CINV(I,J)*TERM)/(SAK(I)*SAK(J))
870	CONTINUE
	IF(NOUT.EQ.0) GOTO 890
	DO 880 I = 1,NOUT
		K = NS(I)
		SAK(K) = 0.D0
880		SAKR(K) = 0.D0
890	CONTINUE
C
	IF(IWTC.NE.5) WRITE(*,900)
900	FORMAT(/10X,' *** GENERAL METHOD SOLUTION ***')
	IF(IDIRC.GT.MAXG) THEN
	      WRITE(*,902) IDIRC
902	      FORMAT(8X,'!!$$$**  N0 CONVERGENCE IN',I3,' ITERATIONS  **
     1$$$!!'//
     1  5X,'### RESET PARAMETER AND DATA VALUES, THEN DO MENU 1A #5'/
     1  5X,'$$$       AND/OR:  GO TO MENU 1A #6, SET G<1, AND RERUN'/
     1  5X,'### AFTER CONVERGENCE, USE NEW PARAMETER VALUES WITH G=1'/
     1  5X,'AND RUN AGAIN AS A CHECK'/)
		IRET = 1
	ELSE
	      IF(IWTC.NE.5) WRITE(*,905) IDIRC 
905     FORMAT(5X,27H ** CONVERGENCE ACHIEVED IN,I3,14H ITERATIONS **/)
	ENDIF
C
C
518	RETURN
	END
