SUBROUTINE MERS(WW,WEND,YY,N,ACCC,ACCC0,H)
C
C
C   **********************************************************************
C   SUBROUTINE  "MERS" CONTAINS  A GENERIC NUMERICAL INTEGRATOR (4TH ORDER 
C   RUNGE-KUTTA-MERSON, WITH STEP-SIZE CONTROL) USED TO INTEGRATE A SET OF 
C   ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM
C
C        D YI/D X = F_I(X,Y1,Y2,....YN)
C
C   WHERE X IS THE INDEPENDENT VARIABLE  AND YI ARE N DEPENDENT VARIABLES.
C   "WW" DENOTES THE INITIAL AND "WEND" THE FINAL VALUE OF THE INDEPENDENT 
C   VARIABLE.  ARRAY "YY" STORES  THE VALUES OF  THE  DEPENDENT  VARIABLES  
C   ("N" SPECIFIES THEIR NUMBER).  "ACCC" SPECIFIES  THE  MAXIMUM  ALLOWED 
C   RELATIVE  ERROR  FOR ANY  DEPENDENT VARIABLE, "H" THE INTEGRATION STEP 
C   SIZE; (IF  "ACCC<ACCC0"  THE STEP SIZE  IS NOT INCREASED;  USEFUL AT A
C   BOUNDARY).
C   **********************************************************************
C
C
	IMPLICIT REAL*8 (A-H,O-Z)		
	DIMENSION YY(N),F(N),YZ(N),A(N),B(N),U(N)
	DATA RZERO/1.E-30/
C
C  DEFINITION OF LOCAL VARIABLES
C
	DO 100 K=1,N
 100	  U(K)=YY(K)
	ZZ=WW
 10	CONTINUE
C
C  ----------------------------------------
C  RUNGE-KUTTA-MERSON INTEGRATION ESTIMATES
C  ----------------------------------------
C
C  K1 = H*F1(WW,YY)
C
	DO 101 K=1,N
 101	  YZ(K)=U(K)
C
	CALL EXT(ZZ,U,F,N)
C
C  K2 = H*F2(WW+H/3,YY+K1/3)
C
	H3=1./3.*H
	ZZ=ZZ+H3
	DO 102 K=1,N
	  A(K)=H3*F(K)
 102	  U(K)=A(K)+YZ(K)
C
	CALL EXT(ZZ,U,F,N)
C
C  K3 = H*F3(WW+H/3,YY+K1/6+K2/6)
C
	DO 103 K=1,N
	  A(K)=0.5*A(K)
 103	  U(K)=0.5*H3*F(K)+A(K)+YZ(K)
C
	CALL EXT(ZZ,U,F,N)
C
C  K4 = H*F4(WW+H/2,YY+K1/8+3*K3/8)
C
	ZZ=ZZ+0.5*H3
	DO 104 K=1,N
	  B(K)=4.5*H3*F(K)
 104	  U(K)=0.25*B(K)+0.75*A(K)+YZ(K)
C
	CALL EXT(ZZ,U,F,N)
C
C  K5 = H*F5(WW+H,YY+K1/2-3*K3/2+2*K4)
C
	ZZ=ZZ+0.5*H
	DO 105 K=1,N
	  A(K)=2.*H3*F(K)+A(K)
 105	  U(K)=3.*A(K)-B(K)+YZ(K)
C
	CALL EXT(ZZ,U,F,N)
C 20
C
C  ---------------
C  ERROR ESTIMATES
C  ---------------
C
	DO 200 K=1,N
C
C  ESTIMATE OF THE LOCAL ERROR = B(K)/5
C
	  B(K)=-0.5*H3*F(K)-B(K)+2.*A(K)
C
C  VALUES OF INDEPENDENT VARIABLES AT WW+H
C
	  U(K)=U(K)-B(K)
C
C  LOCAL RELATIVE ERROR LARGER THAN ACCC?
C
	  A(K)=ABS(5.*ACCC*U(K))
	  B(K)=ABS(B(K))
	  IF(ABS(U(K)).LE.RZERO) GOTO 200
	  IF(B(K).GT.A(K)) GOTO 40
 200	  CONTINUE
C
C  STEP SIZE UNNECESSARILY SMALL?
C
	DO 201 K=1,N
	  IF(B(K).GT.0.03125*A(K)) GOTO 50
 201	  CONTINUE
C 30
C
C  DOUBLING OF THE STEP SIZE
C  ACCC IS USED TO AVOID INCREMENT OF H NEAR THE OUTER FITTING POINT
C
	IF(ACCC.LT.ACCC0) GOTO 50
	H=H+H
	GOTO 50
C
 40   	CONTINUE
C
C  HALVING OF STEP SIZE (UNLESS H SMALLER THAN 0.00001)
C
	IF(H.LE.0.00001) GOTO 50
	H=0.5*H
	DO 400 K=1,N
 400	  U(K)=YZ(K)
	ZZ=WW
C
C  REDO THE LAST INTEGRATION WITH SMALLER STEP SIZE
C  
	GOTO 10
C
 50	CONTINUE
C
C  PASS LOCAL VARIABLES TO GLOBAL VARIABLES
C
	WEND=ZZ
	DO 500 K=1,N
 500	  YY(K)=U(K)
	RETURN
	END