PROGRAM GILL
C
	PARAMETER(NY=43,NX=21)
	REAL*8 S0(NY),GY(NY),GY2(NY),DUM(NX),YG(NY),XG(NX)
	COMPLEX*16 S(NY),VS(NY),US(NY),W(NY)
	COMPLEX*16 AI,FAC,SUM
	REAL*8 H, HB,HM, LEQ, K, Y,YM,YM2,DY
C
C   ***  Assign Unit 11 to the Input Parameter File   ***
C
	OPEN(UNIT=11,FILE='gill_params.txt',STATUS='OLD')
C
C           ***  Read in PARAMETERS  ***
C
	READ(11,*)
	READ(11,*)
	READ(11,*)CHI
	READ(11,*)B
	READ(11,*)YMAX
	READ(11,*)ALPHA
	READ(11,*)K
C
	SA=-100.0
	AI=(0.0,1.0)
	ANY2=FLOAT((NY+1)/2)
	NY2=(NY+1)/2
c	WL=2.*3.14159/K
	WL=2.*3.14159
	AL=1.5/SQRT(B)
C
	DO I=1,NY
	 Y=2.*AL*(FLOAT(I-NY2))/ANY2
	 YM=Y-AL/ANY2
	 YM2=Y+AL/ANY2
	 CALL GFUNC(SA,B,Y,YMAX,G)
	 S0(I)=G
	 CALL GFUNC(SA,B,YM,YMAX,G)
	 GY(I)=G
	 CALL GFUNC(SA,B,YM2,YMAX,G)
	 GY2(I)=G
	END DO	
	US(1)=0.0
	S(NY2)=0.0
	SUM=0.0
	DO I=NY2+1,NY
	 Y=2.*AL*(FLOAT(I-NY2))/ANY2
	 YM=Y-AL/ANY2
	 DY=2.*AL/ANY2
	 SUM=SUM+(YM**(1.+AI*K/ALPHA))*GY(I)*EXP(-(0.5*CHI
     1  /ALPHA)*YM*YM)*DY
	 S(I)=SUM*(Y**(-AI*K/ALPHA))*EXP(0.5*CHI*Y*Y/ALPHA)
	 VS(I)=-AI*K*S(I)/Y
	 US(I)=(S(I)-S(I-1))/(DY*YM)
	END DO
	SUM=0.0
	DO I=NY2-1,1,-1
	 Y=2.*AL*(FLOAT(I-NY2))/ANY2
	 YM=Y+AL/ANY2
	 DY=-2.*AL/ANY2
	 SUM=SUM+(YM**(1.+AI*K/ALPHA))*GY2(I)*EXP(-(0.5*CHI/
     1   ALPHA)*YM*YM)*DY
	 S(I)=SUM*(Y**(-AI*K/ALPHA))*EXP(0.5*CHI*Y*Y/ALPHA)
	 VS(I)=-AI*K*S(I)/Y
	 US(I+1)=(S(I)-S(I+1))/(DY*YM)
	END DO
	US(NY2+2)=2.*US(NY2+3)-US(NY2+4)
	US(NY2+1)=US(NY2+2)
	US(NY2-1)=2.*US(NY2-2)-US(NY2-3)
	US(NY2)=US(NY2-1)
	DO I=1,NY-1
	 W(I)=S0(I)-CHI*S(I)+0.5*ALPHA*(US(I)+US(I+1))
	END DO
	W(NY)=S0(NY)-CHI*S(NY)+ALPHA*US(NY)
	US(1)=0.0
	VS(NY2)=0.5*(VS(NY2+1)+VS(NY2-1))
c	W(NY2)=1.5*W(NY2+1)-0.5*W(NY2+3)
C
	OPEN(UNIT=10,FILE='u.out')
	OPEN(UNIT=11,FILE='v.out')	
	OPEN(UNIT=12,FILE='w.out')
	OPEN(UNIT=13,FILE='s.out')
	OPEN(UNIT=14,FILE='s0.out')
	OPEN(UNIT=15,FILE='x.out')
	OPEN(UNIT=16,FILE='y.out')
	OPEN(UNIT=17,FILE='params.out')
	OPEN(UNIT=18,FILE='params2.out')
c
	WRITE(17,250) 'CHI= ',CHI,'  B=',B,'  K=',K,'  ALPHA=',ALPHA,
     1  '   (u,v,s)'
	WRITE(18,250) 'CHI= ',CHI,'  B=',B,'  K=',K,'  ALPHA=',ALPHA,
     1  '   (s_0,w)'
C
	DO I=1,NY
	 Y=2.*AL*(FLOAT(I-NY2))/ANY2
	 YG(I)=Y
	 DO J=1,NX
	  X=WL*(0.25+FLOAT(J-1)/(FLOAT(NX-1)))
	  XG(J)=(X-0.25*WL)
	  DUM(J)=0.1*REAL(EXP(AI*K*X)*US(I))
	 END DO
	 WRITE(10,100)(DUM(J),J=1,NX)
	 DO J=1,NX
	  X=WL*(0.25+FLOAT(J-1)/(FLOAT(NX-1)))
	  DUM(J)=0.2*REAL(EXP(AI*K*X)*VS(I))
	 END DO
	 WRITE(11,100)(DUM(J),J=1,NX)
	 DO J=1,NX
	  X=WL*(0.25+FLOAT(J-1)/(FLOAT(NX-1)))
	  DUM(J)=REAL(EXP(AI*K*X)*W(I))
	 END DO
	 WRITE(12,100)(DUM(J),J=1,NX)
	 DO J=1,NX
	  X=WL*(0.25+FLOAT(J-1)/(FLOAT(NX-1)))
	  DUM(J)=REAL(EXP(AI*K*X)*S(I))
	 END DO
	 WRITE(13,100)(DUM(J),J=1,NX)
	 DO J=1,NX
	  X=WL*(0.25+FLOAT(J-1)/(FLOAT(NX-1)))
	  DUM(J)=REAL(EXP(AI*K*X)*S0(I))
	 END DO
	 WRITE(14,100)(DUM(J),J=1,NX)
	END DO
	WRITE(16,200)(YG(I),I=1,NY)
	WRITE(15,200)(XG(I),I=1,NX)
  100	FORMAT(2X,101(F7.2,1X))
  200	FORMAT(2X,101(F6.3,1X))
  250 FORMAT(A5,F4.1,A4,F4.1,A4,F4.1,A8,F4.1,A10)
C
	STOP
	END
C
	SUBROUTINE GFUNC(SA,B,Y,YMAX,G)
	REAL*8 Y
c	G=SA*(1.+2.*SQRT(B)*Y)*EXP(-B*Y*Y)
C	G=SA*EXP(-B*Y*Y)
	G=SA*EXP(-B*(Y-YMAX)**2)
	RETURN
	END