Index: trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90	(revision 552)
+++ trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90	(revision 645)
@@ -33,12 +33,13 @@
       
 
-      real hco2(ncompdiff),ho
-
-      integer ig,nz,l,k,n,nn,iq,kn,p,il0,kl,igen,igen2,kl2,ngen,iter,il1,ij0
+!      real hco2(ncompdiff),ho
+
+      integer,dimension(nqmx) :: indic_diff
+      integer ig,nz,l,k,n,nn,p,ij0
       integer istep,il,gcn,ntime,nlraf
       real*8 masse
       real*8 masseU,kBolt,RgazP,Rmars,Grav,Mmars
-      real*8 rho0,D0,T0,H0,time0,dZ,time,dZraf,tdiff,Hup,Zmin,Zmax
-      real*8 hh,dcoef,dcoef1,ptfac, ntot, dens,mgcn,FacEsc
+      real*8 rho0,D0,T0,H0,time0,dZ,time,dZraf,tdiff,Zmin,Zmax
+      real*8 FacEsc,invsgmu
       real*8 hp(nlayermx)
       real*8 pp(nlayermx)
@@ -46,9 +47,7 @@
       real*8 tt(nlayermx),tnew(nlayermx),tint(nlayermx)
       real*8 zz(nlayermx)
-      real*8 qq(nlayermx,ncompdiff)
-      real*8 qnew(nlayermx,ncompdiff)
-      real*8 qint(nlayermx,ncompdiff),FacMass(nlayermx,ncompdiff)
-      real*8 rhoK(nlayermx,ncompdiff),rhoT(nlayermx)
-      real*8 rhoKinit(nlayermx,ncompdiff)
+      real*8,dimension(:,:),allocatable :: qq,qnew,qint,FacMass
+      real*8,dimension(:,:),allocatable :: rhoK,rhokinit
+      real*8 rhoT(nlayermx)
       real*8 dmmeandz(nlayermx)
       real*8 massemoy(nlayermx)
@@ -61,184 +60,116 @@
       real*8,dimension(:),allocatable :: Atri,Btri,Ctri,Dtri,Xtri,Tad,Dad,Zad,rhoad
       real*8,dimension(:),allocatable :: alpha,beta,gama,delta,eps
-      real*8 wi(ncompdiff),Wad(ncompdiff),Uthermal(ncompdiff)
-      real*8 Lambdaexo(ncompdiff)
-      real*8 MToT1(ncompdiff),Mtot2(ncompdiff)
-      real*8 MRaf1(ncompdiff),Mraf2(ncompdiff)
+      real*8,dimension(:),allocatable ::  wi,Wad,Uthermal,Lambdaexo,Hspecie
+      real*8,dimension(:),allocatable :: Mtot1,Mtot2,Mraf1,Mraf2
+      character(len=20),dimension(14) :: ListeDiff
 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 !     tracer numbering in the molecular diffusion
 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-!  Atomic oxygen must always be the LAST species of the list as
-! it is the dominant species at high altitudes. 
-      integer,parameter :: i_o   = 1
-      integer,parameter :: i_n2   = 2
-      integer,parameter :: i_co   = 3
-      integer,parameter :: i_ar  = 4
-      integer,parameter :: i_h2   = 5
-      integer,parameter :: i_h    = 6
-      integer,parameter :: i_o2   = 7
-      integer,parameter :: i_oh  = 8
-      integer,parameter :: i_ho2  = 9
-      integer,parameter :: i_h2o = 10
-      integer,parameter :: i_h2o2  = 11
-      integer,parameter :: i_o1d   = 12
-!      integer,parameter :: i_o3   = 13
-      integer,parameter :: i_n    = 13
-      integer,parameter :: i_no   = 14
-      integer,parameter :: i_no2  = 15
-!      integer,parameter :: i_n2d  = 17
-!      integer,parameter :: i_oplus = 18
-      integer,parameter :: i_co2    = 16
-!      integer,parameter :: i_oplus = 17
-!      integer,parameter :: i_hplus = 18
+!  We need the index of escaping species 
+
+      integer,save :: i_h2  
+      integer,save :: i_h
+! vertical index limit for the molecular diffusion
+      integer,save :: il0  
 
 ! Tracer indexes in the GCM:
-      integer,save :: g_co2=0
-      integer,save :: g_co=0
-      integer,save :: g_o=0
-      integer,save :: g_o1d=0
-      integer,save :: g_o2=0
-      integer,save :: g_o3=0
-      integer,save :: g_h=0
-      integer,save :: g_h2=0
-      integer,save :: g_oh=0
-      integer,save :: g_ho2=0
-      integer,save :: g_h2o2=0
-      integer,save :: g_n2=0
-      integer,save :: g_ar=0
-      integer,save :: g_h2o=0
-      integer,save :: g_n=0
-      integer,save :: g_no=0
-      integer,save :: g_no2=0
-      integer,save :: g_n2d=0
+!      integer,save :: g_co2=0
+!      integer,save :: g_co=0
+!      integer,save :: g_o=0
+!      integer,save :: g_o1d=0
+!      integer,save :: g_o2=0
+!      integer,save :: g_o3=0
+!      integer,save :: g_h=0
+!      integer,save :: g_h2=0
+!      integer,save :: g_oh=0
+!      integer,save :: g_ho2=0
+!      integer,save :: g_h2o2=0
+!      integer,save :: g_n2=0
+!      integer,save :: g_ar=0
+!      integer,save :: g_h2o=0
+!      integer,save :: g_n=0
+!      integer,save :: g_no=0
+!      integer,save :: g_no2=0
+!      integer,save :: g_n2d=0
 !      integer,save :: g_oplus=0
 !      integer,save :: g_hplus=0
 
-      integer,save :: gcmind(ncompdiff) ! array of GCM indexes
+      integer,save :: ncompdiff
+      integer,dimension(:),allocatable,save :: gcmind ! array of GCM indexes
       integer ierr
 
       logical,save :: firstcall=.true.
-      real abfac(ncompdiff)
-      real,save :: dij(ncompdiff,ncompdiff)
+!      real abfac(ncompdiff)
+      real,dimension(:,:),allocatable,save :: dij
       real,save :: step
 
 ! Initializations at first call
       if (firstcall) then
-        call moldiffcoeff_red(dij)
+
+! Liste des especes qui sont diffuses si elles sont presentes
+!	ListeDiff=['co2','o','n2','ar','co','h2','h','d2','d','hd','o2','oh','ho2','h2o_vap','h2o2','o1d','n','no','no2']
+	ListeDiff(1)='co2'
+	ListeDiff(2)='o'
+	ListeDiff(3)='n2'
+	ListeDiff(4)='ar'
+	ListeDiff(5)='co'
+	ListeDiff(6)='h2'
+	ListeDiff(7)='h'
+	ListeDiff(8)='d2'
+	ListeDiff(9)='hd'
+	ListeDiff(10)='d'
+	ListeDiff(11)='o2'
+	ListeDiff(12)='h2o_vap'
+	ListeDiff(13)='o3'
+	ListeDiff(14)='n'
+	i_h=1000
+	i_h2=1000
+! On regarde quelles especes diffusables sont presentes
+
+	ncompdiff=0
+	indic_diff(1:nqmx)=0
+	
+	do nn=1,nqmx
+	do n=1,14
+	if (ListeDiff(n) .eq. noms(nn)) then
+	indic_diff(nn)=1
+	if (noms(nn) .eq. 'h') i_h=n
+	if (noms(nn) .eq. 'h2') i_h2=n
+	endif
+
+	enddo
+	if (indic_diff(nn) .eq. 1) then
+	print*,'specie', noms(nn), 'diffused in moldiff_red'
+	ncompdiff=ncompdiff+1
+	endif
+	enddo
+	print*,'number of diffused species:',ncompdiff
+	allocate(gcmind(ncompdiff))
+
+! Store gcm indexes in gcmind
+	n=0
+	do nn=1,nqmx
+	if (indic_diff(nn) .eq. 1) then
+	n=n+1
+	gcmind(n)=nn
+	endif
+	enddo
+
+	print*,'gcmind',gcmind
+
+! find vertical index above which diffusion is computed
+
+	do l=1,nlayermx
+	if (pplay(1,l) .gt. Pdiff) then
+	il0=l
+	endif
+	enddo
+	il0=il0+1
+	print*,'vertical index for diffusion',il0,pplay(1,il0)
+
+	allocate(dij(ncompdiff,ncompdiff))
+
+        call moldiffcoeff_red(dij,indic_diff,gcmind,ncompdiff)
         print*,'MOLDIFF  EXO'
-
-        ! identify the indexes of the tracers we'll need
-        g_co2=igcm_co2
-        if (g_co2.eq.0) then
-          write(*,*) "moldiff: Error; no CO2 tracer !!!"
-          stop
-        endif
-        g_n2=igcm_n2
-        if (g_n2.eq.0) then
-          write(*,*) "moldiff: Error; no N2 tracer !!!"
-          stop
-        endif
-        g_ar=igcm_ar
-        if (g_ar.eq.0) then
-          write(*,*) "moldiff: Error; no Ar tracer !!!"
-          stop
-        endif
-        g_h2=igcm_h2
-        if (g_h2.eq.0) then
-          write(*,*) "moldiff: Error; no H2 tracer !!!"
-          stop
-        endif
-        g_h=igcm_h
-        if (g_h.eq.0) then
-          write(*,*) "moldiff: Error; no H tracer !!!"
-          stop
-        endif
-        g_co=igcm_co
-        if (g_co.eq.0) then
-          write(*,*) "moldiff: Error; no CO tracer !!!"
-          stop
-        endif
-        g_o2=igcm_o2
-        if (g_o2.eq.0) then
-          write(*,*) "moldiff: Error; no O2 tracer !!!"
-          stop
-        endif
-        g_oh=igcm_oh
-        if (g_oh.eq.0) then
-          write(*,*) "moldiff: Error; no OH tracer !!!"
-          stop
-        endif
-        g_ho2=igcm_ho2
-        if (g_ho2.eq.0) then
-          write(*,*) "moldiff: Error; no HO2 tracer !!!"
-          stop
-        endif
-        g_h2o=igcm_h2o_vap
-        if (g_h2o.eq.0) then
-          write(*,*) "moldiff: Error; no H2O tracer !!!"
-          stop
-        endif
-        g_h2o2=igcm_h2o2
-        if (g_h2o2.eq.0) then
-          write(*,*) "moldiff: Error; no H2O2 tracer !!!"
-          stop
-        endif
-        g_o1d=igcm_o1d
-        if (g_o1d.eq.0) then
-          write(*,*) "moldiff: Error; no O(1D) tracer !!!"
-          stop
-        endif
-        g_o3=igcm_o3
-        if (g_o3.eq.0) then
-          write(*,*) "moldiff: Error; no O3 tracer !!!"
-          stop
-        endif
-        g_n=igcm_n
-        if (g_n.eq.0) then
-          write(*,*) "moldiff: Error; no N tracer !!!"
-          stop
-        endif
-        g_no=igcm_no
-        if (g_no.eq.0) then
-          write(*,*) "moldiff: Error; no NO tracer !!!"
-          stop
-        endif
-        g_no2=igcm_no2
-        if (g_no2.eq.0) then
-          write(*,*) "moldiff: Error; no NO2 tracer !!!"
-          stop
-        endif
-        g_n2d=igcm_n2d
-        if (g_n2d.eq.0) then
-          write(*,*) "moldiff: Error; no N2(D) tracer !!!"
-          stop
-        endif
-        g_o=igcm_o
-        if (g_o.eq.0) then
-          write(*,*) "moldiff: Error; no O tracer !!!"
-          stop
-        endif 
-
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-!   fill array to relate local indexes to gcm indexes
-!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-        gcmind(i_co2)= g_co2
-        gcmind(i_n2)= g_n2
-        gcmind(i_ar)= g_ar
-        gcmind(i_h2)= g_h2
-        gcmind(i_h)= g_h
-        gcmind(i_co)= g_co
-        gcmind(i_o2)= g_o2
-        gcmind(i_oh)= g_oh
-        gcmind(i_ho2)= g_ho2
-        gcmind(i_h2o)= g_h2o
-        gcmind(i_h2o2)= g_h2o2
-        gcmind(i_o1d)= g_o1d
-!        gcmind(i_o3)=g_o3
-        gcmind(i_n)=g_n
-        gcmind(i_no)=g_no
-        gcmind(i_no2)=g_no2
-!        gcmind(i_n2d)=g_n2d
-        gcmind(i_o)=g_o
 
         firstcall= .false.
@@ -257,9 +188,27 @@
 	Grav=6.67d-11
 	Mmars=6.4d23
-      il0=idifflow ! Diffusion limited to regions above il0 JYC  
-	ij0=302 ! For test
+	ij0=6000 ! For test
+	invsgmu=1d0/g/masseU
 	
+! allocatation des tableaux dependants du nombre d especes diffusees
+	allocate(qq(nlayermx,ncompdiff))
+	allocate(qnew(nlayermx,ncompdiff))
+	allocate(qint(nlayermx,ncompdiff))
+	allocate(FacMass(nlayermx,ncompdiff))
+	allocate(rhok(nlayermx,ncompdiff))
+	allocate(rhokinit(nlayermx,ncompdiff))
+
+	allocate(wi(ncompdiff))
+	allocate(wad(ncompdiff))
+	allocate(uthermal(ncompdiff))
+	allocate(lambdaexo(ncompdiff))
+	allocate(Hspecie(ncompdiff))
+	allocate(Mtot1(ncompdiff))
+	allocate(Mtot2(ncompdiff))
+	allocate(Mraf1(ncompdiff))
+	allocate(Mraf2(ncompdiff))
+
+!	print*,'moldiff',i_h2,i_h,ncompdiff
       do ig=1,ngridmx
-
 	pp=dble(pplay(ig,:))
 
@@ -310,5 +259,5 @@
 	enddo
 
-	Zmin=zz(idifflow)
+	Zmin=zz(il0)
 	Zmax=zz(nlayermx)
 
@@ -348,5 +297,5 @@
         allocate(alpha(nlraf),beta(nlraf),gama(nlraf),delta(nlraf),eps(nlraf))
 
-! before beginning, I use a better vertical resolution above il1, 
+! before beginning, I use a better vertical resolution above il0, 
 ! altitude grid reinterpolation 
 ! The diffusion is solved on an altitude grid with constant step dz
@@ -404,5 +353,5 @@
 !	enddo
 
-! The diffusion is computed above il0 (defined in diffusion.h)
+! The diffusion is computed above il0 computed from Pdiff in diffusion.h
 ! No change below il0
 	
@@ -435,4 +384,7 @@
 ! Loop on species
 
+	T0=Traf(nlraf)
+	rho0=1d0
+
         do nn=1,ncompdiff
 	masse=dble(mmol(gcmind(nn)))
@@ -443,4 +395,5 @@
 	CALL HSCALEREAL(nn,Nrafk,Dzraf,Hraf,nlraf,ncompdiff)
 	Hrafmol(:,nn)=Hraf(:)
+!	Hspecie(nn)=kbolt*T0/masse*invsgmu
 ! Computation of the diffusion vertical velocity of the specie
 	CALL VELVERT(nn,Traf,Hraf,Draf,Dzraf,masse,Wraf,nlraf)
@@ -451,14 +404,23 @@
 	enddo
 ! We use a lower time step
-	Tdiff=minval(Tdiffrafmol)/10d0
-
+	Tdiff=minval(Tdiffrafmol)
+	Tdiff=minval(Tdiffrafmol(nlraf,:))*Mraf(nlraf)
 ! Some problems when H is dominant 
 ! The time step is chosen function of atmospheric mass at higher level 
 ! It is not very nice
 
-	if (tdiff .lt. tdiffmin) tdiff=tdiffmin*Mraf(nlraf)
+!	if (ig .eq. ij0) then 
+!	print*,'test',ig,tdiff,tdiffmin,minloc(Tdiffrafmol),minloc(Tdiffrafmol(nlraf,:))
+!	endif
+	if (tdiff .lt. tdiffmin*Mraf(nlraf)) tdiff=tdiffmin*Mraf(nlraf)
 
 ! Number of time step
 	ntime=anint(dble(ptimestep)/tdiff)
+!	print*,'ptime',ig,ptimestep,tdiff,ntime,tdiffmin,Mraf(nlraf)
+! Adimensionned temperature
+
+	do l=1,nlraf
+	Tad(l)=Traf(l)/T0
+	enddo
 
 	do istep=1,ntime
@@ -469,7 +431,5 @@
 ! Parameters to adimension the problem 
 
-	rho0=1d0
-	T0=Traf(nlraf)
-	H0=kbolt*T0/dble(mmol(gcmind(nn)))/g/masseU
+	H0=kbolt*T0/dble(mmol(gcmind(nn)))*invsgmu
 	D0=Draf(nlraf)
 	Time0=H0*H0/D0
@@ -479,11 +439,14 @@
 
 	do l=1,nlraf
-	Tad(l)=Traf(l)/T0
 	Dad(l)=Draf(l)/D0
 	Zad(l)=Zraf(l)/H0
-	RhoAd(l)=Rrafk(l,nn)/rho0
 	enddo
 	Wad(nn)=wi(nn)*Time0/H0
 	DZ=Zad(2)-Zad(1)
+	FacEsc=exp(-wad(nn)*DZ/Dad(nlraf-1))
+
+	do l=1,nlraf
+	RhoAd(l)=Rrafk(l,nn)/rho0
+	enddo
 
 ! Compute intermediary coefficients
@@ -494,6 +457,4 @@
 ! Compute escape factor (exp(-ueff*dz/D(nz)))
 
-	FacEsc=exp(-wad(nn)*DZ/Dad(nlraf-1))
-
 ! Compute matrix coefficients
 
@@ -509,8 +470,8 @@
 !	Xtri=rhoAd
 
-!	if (nn .eq. 6) then
+!	if (ig .eq. ij0 .and. (nn .eq. 1 .or. nn .eq. 5 .or. nn .eq. 6 .or. nn .eq. 16)) then
 !	do l=1,nlraf
 !	if (Xtri(l) .lt. 0.) then
-!	print*,'l',l,rhoAd(l)*rho0,Xtri(l)*rho0,nn
+!	print*,'l',l,rhoAd(l)*rho0,Xtri(l)*rho0,nn,Tad(l),Zad(l),Dad(l)
 !	stop
 !	endif
@@ -533,15 +494,19 @@
 	print*,'pb moldiff',istep,ig,l,nn,Rrafk(l,nn),tdiff,& 
      &  rho0*Rhoad(l),Zmin,Zmax,ntime
-!	print*,'Atri',Atri
-!	print*,'Btri',Btri
-!	print*,'Ctri',Ctri
-!	print*,'Dtri',Dtri
-!	print*,'Xtri',Xtri
-!	print*,'alpha',alpha
-!	print*,'beta',beta
-!	print*,'gama',gama
-!	print*,'delta',delta
-!	print*,'eps',eps
-!	print*,'Dad',Dad
+	print*,'Atri',Atri
+	print*,'Btri',Btri
+	print*,'Ctri',Ctri
+	print*,'Dtri',Dtri
+	print*,'Xtri',Xtri
+	print*,'alpha',alpha
+	print*,'beta',beta
+	print*,'gama',gama
+	print*,'delta',delta
+	print*,'eps',eps
+	print*,'Dad',Dad
+	print*,'Temp',Traf
+	print*,'alt',Zraf
+	print*,'Mraf',Mraf
+	stop
 !	pdqdiff(1:ngridmx,1:nlayermx,1:nqmx)=0.
 !	return
@@ -593,9 +558,13 @@
 	CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayermx,ncompdiff)
 
-!	if (ig .eq. ij0) then 
-!	do l=il0,nlayermx
-!        write(*,'(i2,1x,19(e12.4,1x))') l,zz(l),tt(l),RHOK(l,16)/sum(RHOK(l,:)),RHOKINIT(l,16)/sum(RHOKINIT(l,:))
-!        enddo
-!	endif
+	if (ig .eq. ij0) then 
+	do l=il0,nlayermx
+        write(*,'(i2,1x,19(e12.4,1x))') l,zz(l),tt(l),RHOK(l,1)/sum(RHOK(l,:)),RHOKINIT(l,1)/sum(RHOKINIT(l,:)),&
+     &  RHOK(l,5)/sum(RHOK(l,:)),RHOKINIT(l,5)/sum(RHOKINIT(l,:)),&
+     &  RHOK(l,3)/sum(RHOK(l,:)),RHOKINIT(l,3)/sum(RHOKINIT(l,:)),&
+     &  RHOK(l,6)/sum(RHOK(l,:)),RHOKINIT(l,6)/sum(RHOKINIT(l,:)),&
+     &  RHOK(l,7)/sum(RHOK(l,:)),RHOKINIT(l,7)/sum(RHOKINIT(l,:))
+        enddo
+	endif
 
 ! Compute total mass of each specie on the GCM vertical grid
@@ -640,4 +609,12 @@
 
        enddo  ! ig loop		
+
+	deallocate(qq,qnew,qint)
+	deallocate(FacMass)
+	deallocate(rhok,rhokinit)
+	deallocate(wi,wad,uthermal,lambdaexo)
+	deallocate(Hspecie)
+	deallocate(Mtot1,Mtot2,Mraf1,Mraf2)
+
       return
       end
@@ -1242,40 +1219,60 @@
 	INTEGER :: nl,l
 	REAL*8,DIMENSION(nl) :: T,D
-	REAL*8 :: DZ
+	REAL*8 :: DZ,DZinv
 	REAL*8,DIMENSION(nl) :: alpha,beta,gama,delta,eps
 
 ! Compute alpha,beta and delta
 ! lower altitude values
-
-	alpha(1)=1D0/T(1)*(-3D0*T(1)+4D0*T(2)-T(3))/(2D0*DZ)
-	delta(1)=(-3D0*D(1)+4D0*D(2)-D(3))/(2D0*DZ)
-	beta(1)=1D0/T(1)
-
-	do l=2,nl-1
-	alpha(l)=1D0/T(l)*(T(l+1)-T(l-1))/(2D0*DZ)
-	delta(l)=(D(l+1)-D(l-1))/(2D0*DZ)
+	DZinv=1d0/(2D0*DZ)
+
+	beta(1)=1d0/T(1)
+	alpha(1)=beta(1)*(-3D0*T(1)+4D0*T(2)-T(3))*Dzinv
+	delta(1)=(-3D0*D(1)+4D0*D(2)-D(3))*Dzinv
+
+	beta(2)=1d0/T(2)
+	alpha(2)=beta(2)*(T(3)-T(1))*Dzinv
+	delta(2)=(D(3)-D(1))*Dzinv
+
+!	do l=2,nl-1
+!	beta(l)=1D0/T(l)
+!	alpha(l)=beta(l)*(T(l+1)-T(l-1))*Dzinv
+!	delta(l)=(D(l+1)-D(l-1))*Dzinv
+!	end do   
+
+	do l=3,nl-1
 	beta(l)=1D0/T(l)
-	end do   
+	alpha(l)=beta(l)*(T(l+1)-T(l-1))*Dzinv
+	delta(l)=(D(l+1)-D(l-1))*Dzinv
+	gama(l-1)=(beta(l)-beta(l-2))*Dzinv
+	eps(l-1)=(alpha(l)-alpha(l-2))*Dzinv
+	enddo
 
 ! Upper altitude values
 
-	alpha(nl)=1D0/T(nl)*(3D0*T(nl)-4D0*T(nl-1)+T(nl-2))/(2D0*DZ)     
-	delta(nl)=(3D0*D(nl)-4D0*D(nl-1)+D(nl-2))/(2D0*DZ)
 	beta(nl)=1D0/T(nl)
+	alpha(nl)=beta(nl)*(3D0*T(nl)-4D0*T(nl-1)+T(nl-2))*Dzinv     
+	delta(nl)=(3D0*D(nl)-4D0*D(nl-1)+D(nl-2))*Dzinv
 
 ! Compute the gama and eps coefficients
 ! Lower altitude values
 
-	gama(1)=(-3D0*beta(1)+4D0*beta(2)-beta(3))/(2d0*dz)
-	eps(1)=(-3D0*alpha(1)+4D0*alpha(2)-alpha(3))/(2d0*dz)
-
-	do l=2,nl-1
-	gama(l)=(beta(l+1)-beta(l-1))/(2d0*Dz)
-	eps(l)=(alpha(l+1)-alpha(l-1))/(2d0*Dz)
-	end do
-
-	gama(nl)=(3D0*beta(nl)-4D0*beta(nl-1)+beta(nl-2))/(2D0*DZ)
-	eps(nl)=(3D0*alpha(nl)-4D0*alpha(nl-1)+alpha(nl-2))/(2D0*DZ)
-
+	gama(1)=(-3D0*beta(1)+4D0*beta(2)-beta(3))*Dzinv
+	eps(1)=(-3D0*alpha(1)+4D0*alpha(2)-alpha(3))*Dzinv
+
+	gama(nl-1)=(beta(nl)-beta(nl-2))*Dzinv
+	eps(nl-1)=(alpha(nl)-alpha(nl-2))*Dzinv
+
+!	do l=2,nl-1
+!	gama(l)=(beta(l+1)-beta(l-1))*Dzinv
+!	eps(l)=(alpha(l+1)-alpha(l-1))*Dzinv
+!	end do
+
+	gama(nl)=(3D0*beta(nl)-4D0*beta(nl-1)+beta(nl-2))*Dzinv
+	eps(nl)=(3D0*alpha(nl)-4D0*alpha(nl-1)+alpha(nl-2))*Dzinv
+
+!	do l=1,nl
+!	print*,'test diffparam',alpha(l),beta(l),delta(l),gama(l),eps(l)
+!	enddo
+!	stop
 
 	END
@@ -1541,2 +1538,3 @@
 
         END
+        
