c****************************************************************************
c Version: 0.3
c Last changed: Mario Monteiro, 2005-11
c These set of subroutines uses the data from the following inputs:
c    GONG, FGONG, OSC
c to writte the AMDL (binary) and FAMDL (ascii) output files.
c****************************************************************************
	subroutine gong_amdl (itype,afileout,xin,ndimin,nparin,
     *                        xdata1in,xdata2in,xdata3in,aheaderin)
c
	implicit double precision (b-h,o-y)
	implicit integer (i-n)
	parameter (npt=5000,ncol=40,ndat=200)
	character afileout*80
	dimension ndimin(ndat),xin(ncol,npt),nparin(ndat)
	dimension xdata1in(ndat),xdata2in(ndat),xdata3in(ndat)
	dimension ndimout(ndat),xout(ncol,npt),xdataout(ndat)
c
	common /datanum/c4pi,cl10
	common /dataphy/cg,csig,cmu
c
c------- Translation of the OSC input into the AMDL output:
c
c------- Dimensions:
	ndimout(1)=0
	ndimout(2)=ndimin(3)
	ndimout(3)=5
c
c------- Global parameters:
	xdataout(1)=xdata1in(23)
	xdataout(2)=xdata1in(24)
	xdataout(3)=xdata2in(6)
	xdataout(4)=xdata2in(10)
	xdataout(6)=-xdata2in(14)*xdataout(2)**2/xdataout(4)
	xdataout(7)=0.0d0
	xdataout(8)=0.0d0
c
c------- Mesh points of the model:
	fac1=cg*xdataout(1)/xdataout(2)
	fac2=c4pi*xdataout(2)**3/xdataout(1)
	if (xin(2,1).lt.xin(2,ndimout(2))) then
	   icentre=1
	   xdataout(5)=-xdata2in(13)*xdataout(2)**2/xdataout(3)/
     *                 xin(14,icentre)
c
	   xout(1,1)=0.0d0
	   xout(2,1)=(fac2/3.0d0)*xdata2in(10)
	   xout(3,1)=0.0d0
	   xout(4,1)=xin(14,icentre)
	   xout(5,1)=xin(18,icentre)
	   xout(6,1)=3.0d0
	   do 40 i=2,ndimout(2)
	      rx=exp(cl10*xin(2,i))*1.0d11/xdataout(2)
	      qm=exp(cl10*xin(1,i))
	      xout(1,i)=rx
	      xout(2,i)=qm/rx**3
	      xout(3,i)=fac1*(qm/rx)*xin(10,i)/(xin(14,i)*xin(9,i))
	      xout(4,i)=xin(14,i)
	      xout(5,i)=xin(18,i)
	      xout(6,i)=fac2*xin(10,i)/xout(2,i)
 40	   continue
	else
	   icentre=ndimout(2)
	   xdataout(5)=-xdata2in(13)*xdataout(2)**2/xdataout(3)/
     *                 xin(14,icentre)
c
	   xout(1,1)=0.0d0
	   xout(2,1)=(fac2/3.0d0)*xdata2in(10)
	   xout(3,1)=0.0d0
	   xout(4,1)=xin(14,icentre)
	   xout(5,1)=xin(18,icentre)
	   xout(6,1)=3.0d0
	   do 50 i=2,ndimout(2)
	      rx=exp(cl10*xin(2,icentre-i+1))*1.0d11/xdataout(2)
	      qm=exp(cl10*xin(1,icentre-i+1))
	      xout(1,i)=rx
	      xout(2,i)=qm/rx**3
	      xout(3,i)=fac1*(qm/rx)*xin(10,icentre-i+1)/
     *                  (xin(14,icentre-i+1)*xin(9,icentre-i+1))
	      xout(4,i)=xin(14,icentre-i+1)
	      xout(5,i)=xin(18,icentre-i+1)
	      xout(6,i)=fac2*xin(10,icentre-i+1)/xout(2,i)
 50	continue
	endif
c
	if (itype.eq.14) then
	   call write_amdl (afileout,xout,ndimout,xdataout)   
	else if (itype.eq.15) then
	   call write_famdl (afileout,xout,ndimout,xdataout)
	else
	   stop
	endif
c
	return
	end
c***************************************************************************
	subroutine fgong_amdl (itype,afileout,xin,ndimin,xdatain,aheaderin)
c
	implicit double precision (b-h,o-y)
	implicit integer (i-n)
	parameter (npt=5000,ncol=40,ndat=200)
	character afileout*80
	dimension ndimin(ndat),xin(ncol,npt),xdatain(ndat)
	dimension ndimout(ndat),xout(ncol,npt),xdataout(ndat)
c
	common /datanum/c4pi,cl10
	common /dataphy/cg,csig,cmu
c
c------- Translation of the OSC input into the AMDL output:
c
c------- Dimensions:
	ndimout(1)=0
	ndimout(2)=ndimin(1)
	ndimout(3)=5
c
c------- Global parameters:
	xdataout(1)=xdatain(1)
	xdataout(2)=xdatain(2)
	xdataout(6)=-xdatain(12)
	xdataout(7)=0.0d0
	xdataout(8)=0.0d0
c
c------- Mesh points of the model:
	fac1=cg*xdataout(1)/xdataout(2)
	fac2=c4pi*xdataout(2)**3/xdataout(1)
	if (xin(2,1).lt.xin(2,ndimout(2))) then
	   icentre=1
	   xdataout(3)=xin(4,icentre)
	   xdataout(4)=xin(5,icentre)
	   xdataout(5)=-xdatain(11)/xin(10,icentre)
c
	   xout(1,1)=0.0d0
	   xout(2,1)=(fac2/3.0d0)*xin(5,icentre)
	   xout(3,1)=0.0d0
	   xout(4,1)=xin(10,icentre)
	   xout(5,1)=xin(15,icentre)
	   xout(6,1)=3.0d0
	   do 40 i=2,ndimout(2)
	      rx=xin(1,i)/xdataout(2)
	      qm=exp(xin(2,i))
	      xout(1,i)=rx
	      xout(2,i)=qm/rx**3
	      xout(3,i)=fac1*(qm/rx)*xin(5,i)/(xin(10,i)*xin(4,i))
	      xout(4,i)=xin(10,i)
	      xout(5,i)=xin(15,i)
	      xout(6,i)=fac2*xin(5,i)/xout(2,i)
 40	   continue
	else
	   icentre=ndimout(2)
	   xdataout(3)=xin(4,icentre)
	   xdataout(4)=xin(5,icentre)
	   xdataout(5)=-xdatain(11)/xin(10,icentre)
c
	   xout(1,1)=0.0d0
	   xout(2,1)=(fac2/3.0d0)*xin(5,icentre)
	   xout(3,1)=0.0d0
	   xout(4,1)=xin(10,icentre)
	   xout(5,1)=xin(15,icentre)
	   xout(6,1)=3.0d0
	   do 50 i=2,ndimout(2)
	      rx=xin(1,icentre-i+1)/xdataout(2)
	      qm=exp(xin(2,icentre-i+1))
	      xout(1,i)=rx
	      xout(2,i)=qm/rx**3
	      xout(3,i)=fac1*(qm/rx)*xin(5,icentre-i+1)/
     *                  (xin(10,icentre-i+1)*xin(4,icentre-i+1))
	      xout(4,i)=xin(10,icentre-i+1)
	      xout(5,i)=xin(15,icentre-i+1)
	      xout(6,i)=fac2*xin(5,icentre-i+1)/xout(2,i)
 50	continue
	endif
c
	if (itype.eq.24) then
	   call write_amdl (afileout,xout,ndimout,xdataout)   
	else if (itype.eq.25) then
	   call write_famdl (afileout,xout,ndimout,xdataout)
	else
	   stop
	endif
c
	return
	end
c***************************************************************************
	subroutine osc_amdl  (itype,afileout,xin,ndimin,xdatain)
c
	implicit double precision (b-h,o-y)
	implicit integer (i-n)
	parameter (npt=5000,ncol=40,ndat=200)
	character afileout*80
	dimension ndimin(ndat),xin(ncol,npt),xdatain(ndat)
	dimension ndimout(ndat),xout(ncol,npt),xdataout(ndat)
c
	common /datanum/c4pi,cl10
	common /dataphy/cg,csig,cmu
c
c------- Translation of the OSC input into the AMDL output:
c
c------- Dimensions:
	ndimout(1)=0
	ndimout(2)=ndimin(1)
	ndimout(3)=5
c
c------- Global parameters:
	xdataout(1)=xdatain(1)
	xdataout(2)=xdatain(2)
	xdataout(7)=0.0d0
	xdataout(8)=0.0d0
c
c------- Mesh points of the model:
	fac1=cg*xdataout(1)/xdataout(2)
	fac2=c4pi*xdataout(2)**3/xdataout(1)
	if (xin(1,1).lt.xin(1,ndimout(2))) then
	   icentre=1
	   xdataout(3)=xin(4,icentre)
	   xdataout(4)=xin(5,icentre)
	   xdataout(5)=-xdatain(9)/xin(10,icentre)
	   xdataout(6)=-xdatain(10)
c
	   xout(1,1)=0.0d0
	   xout(2,1)=(fac2/3.0d0)*xin(5,icentre)
	   xout(3,1)=0.0d0
	   xout(4,1)=xin(10,icentre)
	   xout(5,1)=xin(15,icentre)
	   xout(6,1)=3.0d0
	   do 40 i=2,ndimout(2)
	      qm=exp(xin(2,i))
	      xout(1,i)=xin(1,i)/xdataout(2)
	      xout(2,i)=qm/(xout(1,i))**3
	      xout(3,i)=fac1*qm*xin(5,i)/(xin(10,i)*xin(4,i)*xout(1,i))
	      xout(4,i)=xin(10,i)
	      xout(5,i)=xin(15,i)
	      xout(6,i)=fac2*xin(5,i)/xout(2,i)
 40	   continue
	else
	   icentre=ndimout(2)
	   xdataout(3)=xin(4,icentre)
	   xdataout(4)=xin(5,icentre)
	   xdataout(5)=-xdatain(9)/xin(10,icentre)
	   xdataout(6)=-xdatain(10)
c
	   xout(1,1)=0.0d0
	   xout(2,1)=(fac2/3.0d0)*xin(5,icentre)
	   xout(3,1)=0.0d0
	   xout(4,1)=xin(10,icentre)
	   xout(5,1)=xin(15,icentre)
	   xout(6,1)=3.0d0
	   do 50 i=2,ndimout(2)
	      qm=exp(xin(2,icentre-i+1))
	      xout(1,i)=xin(1,icentre-i+1)/xdataout(2)
	      xout(2,i)=qm/(xout(1,i))**3
	      xout(3,i)=fac1*qm*xin(5,icentre-i+1)/(xin(10,icentre-i+1)*
     *                  xin(4,icentre-i+1)*xout(1,i))
	      xout(4,i)=xin(10,icentre-i+1)
	      xout(5,i)=xin(15,icentre-i+1)
	      xout(6,i)=fac2*xin(5,icentre-i+1)/xout(2,i)
 50	continue
	endif

c
	if (itype.eq.34) then
	   call write_amdl (afileout,xout,ndimout,xdataout)   
	else if (itype.eq.35) then
	   call write_famdl (afileout,xout,ndimout,xdataout)
	else
	   stop
	endif
c
	return
	end
c***************************************************************************
	subroutine srox_amdl  (itype,afileout,xin,ndimin,xdatain)
c
	implicit double precision (b-h,o-y)
	implicit integer (i-n)
	parameter (npt=5000,ncol=40,ndat=200)
	character afileout*80
	dimension ndimin(ndat),xin(ncol,npt),xdatain(ndat)
	dimension ndimout(ndat),xout(ncol,npt),xdataout(ndat)
c
	common /datanum/c4pi,cl10
	common /dataphy/cg,csig,cmu
c
c------- Translation of the STAROX input into the AMDL output:
c
c------- Dimensions:
	ndimout(1)=0
	ndimout(2)=ndimin(1)
	ndimout(3)=5
c
c------- Global parameters:
	xdataout(1)=xdatain(3)
	xdataout(2)=xdatain(2)
	xdataout(5)=-xdatain(5)
	xdataout(6)=-xdatain(4)
	xdataout(7)=0.0d0
	xdataout(8)=0.0d0
c
c------- Mesh points of the model:
	fac1=xdatain(1)*xdatain(3)/xdatain(2)
c	fac1=cg*xdataout(1)/xdataout(2)
	fac2=c4pi*xdataout(2)**3/xdataout(1)
	if (xin(1,1).lt.xin(1,ndimout(2))) then
	   icentre=1
	   xdataout(3)=xin(4,icentre)
	   xdataout(4)=xin(5,icentre)
c
	   xout(1,1)=0.0d0
	   xout(2,1)=(fac2/3.0d0)*xin(4,icentre)
	   xout(3,1)=0.0d0
	   xout(4,1)=xin(5,icentre)
	   xout(5,1)=0.0d0
	   xout(6,1)=3.0d0
	   do 40 i=2,ndimout(2)
	      xout(1,i)=xin(1,i)
	      xout(2,i)=xin(2,i)/(xin(1,i))**3
	      dlpdlr=-fac1*xin(2,i)*xin(4,i)/(xin(1,i)*xin(3,i))
	      xout(3,i)=-dlpdlr/xin(5,i)
c	      xout(3,i)=fac1*xin(2,i)*xin(4,i)/(xin(1,i)*xin(3,i)*xin(5,i))
	      xout(4,i)=xin(5,i)
	      xout(5,i)=xin(6,i)*dlpdlr
	      xout(6,i)=fac2*xin(4,i)/xout(2,i)
 40	   continue
	else
	   icentre=ndimout(2)
	   xdataout(3)=xin(4,icentre)
	   xdataout(4)=xin(5,icentre)
c
	   xout(1,1)=0.0d0
	   xout(2,1)=(fac2/3.0d0)*xin(4,icentre)
	   xout(3,1)=0.0d0
	   xout(4,1)=xin(5,icentre)
	   xout(5,1)=0.0d0
	   xout(6,1)=3.0d0
	   do 50 i=2,ndimout(2)
	      xout(1,i)=xin(1,icentre-i+1)
	      xout(2,i)=xin(2,icentre-i+1)/(xin(1,icentre-i+1))**3
	      dlpdlr=-fac1*xin(2,icentre-i+1)*xin(4,icentre-i+1)/
     *               (xin(1,icentre-i+1)*xin(3,icentre-i+1))
	      xout(3,i)=-dlpdlr/xin(5,icentre-i+1)
	      xout(4,i)=xin(5,icentre-i+1)
	      xout(5,i)=xin(6,icentre-i+1)*dlpdlr
	      xout(6,i)=fac2*xin(4,icentre-i+1)/xout(2,i)
 50	continue
	endif

c
	if (itype.eq.64) then
	   call write_amdl (afileout,xout,ndimout,xdataout)   
	else if (itype.eq.65) then
	   call write_famdl (afileout,xout,ndimout,xdataout)
	else
	   stop
	endif
c
	return
	end
c***************************************************************************
	subroutine write_amdl  (afileout,xout,ndimout,xdataout)
c
c Subroutine to write out the AMDL file (binary).
c
	implicit double precision (b-h,o-y)
	implicit integer (i-n)
	parameter (npt=5000,ncol=40,ndat=200)
	character afileout*80
	dimension ndimout(ndat),xout(ncol,npt),xdataout(ndat)
c
	write (*,1000) ndimout(2),afileout
 1000	format ('      Writing AMDL model with ',i4,
     *          ' mesh points, in file: ',a20,/)
c
	open (2,file=afileout,form='unformatted',status='unknown')
c
	write (2) (ndimout(jn),jn=1,2),(xdataout(jd),jd=1,8),
     *            ((xout(j,i),j=1,ndimout(3)+1),i=1,ndimout(2))
c
	close (2)
c
	return
	end
c**************************************************************************
	subroutine write_famdl  (afileout,xout,ndimout,xdataout)
c
c Subroutine to write out the FAMDL file (ascii).
c
	implicit double precision (b-h,o-y)
	implicit integer (i-n)
	parameter (npt=5000,ncol=40,ndat=200)
	character afileout*80
	dimension ndimout(ndat),xout(ncol,npt),xdataout(ndat)
c
	write (*,1000) ndimout(2),afileout
 1000	format ('      Writing FAMDL model with ',i4,
     *          ' mesh points, in file: ',a20,/)
c
	open (2,file=afileout,form='formatted',status='unknown')
c
	write (2,1010) (ndimout(j),j=1,3)
 1010	format (3i10)
c
	write (2,1020) (xdataout(k),k=1,8),
     *                 ((xout(j,i),j=1,ndimout(3)+1),i=1,ndimout(2))
 50	continue
 1020	format (1p4e20.13)
c
	close (2)
c
	return
	end
c**************************************************************************