Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 709)
+++ trunk/LMDZ.MARS/README	(revision 710)
@@ -1728,2 +1728,9 @@
 == 22/06/12 == EM
 - Enforce that imposed TES ice albedo can never be greater than 0.9
+
+== 28/06/12 == JYC+EM
+- Decreased time step for molecular diffusion (diffusion.h) for better stability
+- Added test in moldiff_red.F90 to enforce that the system can't be stuck in an 
+  infinite loop.
+- The creation and output of coefficients in file 'coeffs.dat' in moldiff_red.F
+  is now controled by an internal flag (by default set to false).
Index: trunk/LMDZ.MARS/libf/aeronomars/diffusion.h
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/diffusion.h	(revision 709)
+++ trunk/LMDZ.MARS/libf/aeronomars/diffusion.h	(revision 710)
@@ -10,5 +10,5 @@
 
       parameter (Pdiff=15.)      ! pressure below which diffusion is computed
-      parameter (tdiffmin=10d0)
+      parameter (tdiffmin=5d0)
       parameter (dzres=2d0)	 ! grid resolution (km) for diffusion
       
Index: trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90	(revision 709)
+++ trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90	(revision 710)
@@ -97,5 +97,5 @@
       integer,save :: ncompdiff
       integer,dimension(:),allocatable,save :: gcmind ! array of GCM indexes
-      integer ierr
+      integer ierr,compteur
 
       logical,save :: firstcall=.true.
@@ -312,5 +312,5 @@
 
 	CALL GCMGRID_P(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qint,tt,tint  &
-     &    ,pp,mmol,gcmind,nlraf,ncompdiff,nlayermx)
+     &    ,pp,mmol,gcmind,nlraf,ncompdiff,nlayermx,ig)
 
 ! We compute the mass correction factor of each specie at each pressure level
@@ -554,5 +554,5 @@
 
 	CALL GCMGRID_P2(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qnew,tt,tnew,&
-     &   pp,mmol,gcmind,nlraf,ncompdiff,nlayermx,FacMass)
+     &   pp,mmol,gcmind,nlraf,ncompdiff,nlayermx,FacMass,ig)
 
 	CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayermx,ncompdiff)
@@ -561,7 +561,7 @@
 	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,2)/sum(RHOK(l,:)),RHOKINIT(l,2)/sum(RHOKINIT(l,:)),&
+     &  RHOK(l,6)/sum(RHOK(l,:)),RHOKINIT(l,6)/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
@@ -1352,8 +1352,8 @@
 
 	SUBROUTINE GCMGRID_P(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew,             &
-     &    pp,M,gc,nl,nq,nlx)
+     &    pp,M,gc,nl,nq,nlx,ig)
 	IMPLICIT NONE
 #include "dimensions.h"
-	INTEGER :: nl,nq,nlx,il,nn,iP
+	INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
 	INTEGER,DIMENSION(1) :: indP
 	INTEGER,DIMENSION(nq) :: gc 
@@ -1412,5 +1412,7 @@
 	Pnew2=Pnew
 
+	compteur=0
 	do while (pnew2 .ge. pp(il))	
+	compteur=compteur+1
 	do nn=1,nq
 	Hi=Konst*T(nl)/dble(M(gc(nn)))/g
@@ -1421,4 +1423,9 @@
 	Znew=Znew2
 	Znew2=Znew2+Dz
+	if (compteur .ge. 100000) then
+	print*,'error moldiff_red infinite loop'
+	print*,ig,il,pp(il),tt(nl),pnew2,qnew(il,:),Znew2
+	stop
+	endif
 !	print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
 	enddo
@@ -1446,8 +1453,8 @@
 
 	SUBROUTINE GCMGRID_P2(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew             &
-     &   ,pp,M,gc,nl,nq,nlx,facM)
+     &   ,pp,M,gc,nl,nq,nlx,facM,ig)
         IMPLICIT NONE
 #include "dimensions.h"
-        INTEGER :: nl,nq,nlx,il,nn,iP
+        INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
         INTEGER,DIMENSION(1) :: indP
         INTEGER,DIMENSION(nq) :: gc
@@ -1505,5 +1512,7 @@
         Pnew2=Pnew
 
+	compteur=0
         do while (pnew2 .ge. pp(il))    
+	compteur=compteur+1
         do nn=1,nq
         Hi=Konst*T(nl)/dble(M(gc(nn)))/g
@@ -1514,4 +1523,10 @@
         Znew=Znew2
         Znew2=Znew2+Dz
+	if (compteur .ge. 100000) then
+	print*,'pb moldiff_red infinite loop'
+	print*,ig,nl,T(nl),pnew2,qnew(il,:),Znew2
+	stop
+	endif
+
 !       print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
         enddo
Index: trunk/LMDZ.MARS/libf/aeronomars/moldiffcoeff_red.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/moldiffcoeff_red.F	(revision 709)
+++ trunk/LMDZ.MARS/libf/aeronomars/moldiffcoeff_red.F	(revision 710)
@@ -87,4 +87,6 @@
       real dnh
       logical,save :: firstcall=.true.
+      logical,parameter :: outputcoeffs=.false. ! to output 'coeffs.dat' file,
+                                                ! set outputcoeffs=.true.
 
 ! Initializations at first call (and some sanity checks)
@@ -253,8 +255,8 @@
 	if (noms(n) .eq. 'o') g_o=n
 	enddo
-	print*,'gh2',g_h2,g_h,g_o
-
-       print*,'moldiffcoef: COEFF CALC'
-       open(56,file='coeffs.dat',status='unknown')
+	print*,'moldiffcoeff_red: gh2',g_h2,g_h,g_o
+
+       print*,'moldiffcoeff_red: COEFF CALC'
+
       do n=1,ncompdiff2
 	dijref=0.
@@ -283,10 +285,14 @@
       enddo 
 
-      do n=1,ncompdiff2
+      if (outputcoeffs) then
+       ! output coefficients in 'coeffs.dat' file
+       open(56,file='coeffs.dat',status='unknown')
+       do n=1,ncompdiff2
         do nn=n,ncompdiff2
           write(56,*) n,nn,dij(n,nn)	!*1.e5/1.381e-23/(273**1.75)
         enddo
-      enddo
-      close(56)
+       enddo
+       close(56)
+      endif