Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 1356)
+++ /trunk/LMDZ.MARS/README	(revision 1357)
@@ -2147,2 +2147,5 @@
   measurements, along with dustiropacity='mcs'. Future evolutions of this option
   could include other instruments, or be used for water ice.
+
+== 20/10/2014 == JYC
+- Update and bug correction on the escape fluxes of H and H2.
Index: /trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90	(revision 1356)
+++ /trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90	(revision 1357)
@@ -3,9 +3,20 @@
 
 use tracer_mod, only: noms, mmol
-USE comcstfi_h
-
+use comgeomphy, only: airephy
 implicit none
 
+!#include "dimensions.h"
+!#include "dimphys.h"
+#include "comcstfi.h"
+!#include "callkeys.h"
+!#include "comdiurn.h"
+!#include "chimiedata.h"
+!#include "tracer.h"
+!#include "conc.h"
 #include "diffusion.h"
+
+! July 2014 JYC ADD BALISTIC Transport coupling to compute wup for H and H2
+
+
 
 !
@@ -39,5 +50,5 @@
       real*8 masseU,kBolt,RgazP,Rmars,Grav,Mmars
       real*8 rho0,D0,T0,H0,time0,dZ,time,dZraf,tdiff,Zmin,Zmax
-      real*8 FacEsc,invsgmu
+      real*8 FacEsc,invsgmu,PhiEscH,PhiEscH2
       real*8 hp(nlayer)
       real*8 pp(nlayer)
@@ -132,6 +143,6 @@
 	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
+!	if (noms(nn) .eq. 'h') i_h=n
+!	if (noms(nn) .eq. 'h2') i_h2=n
 	endif
 
@@ -151,8 +162,10 @@
 	n=n+1
 	gcmind(n)=nn
-	endif
-	enddo
-
-	print*,'gcmind',gcmind
+	if (noms(nn) .eq. 'h') i_h=n
+        if (noms(nn) .eq. 'h2') i_h2=n
+	endif
+	enddo
+
+!	print*,'gcmind',gcmind,i_h,i_h2
 
 ! find vertical index above which diffusion is computed
@@ -207,5 +220,9 @@
 	invsgmu=1d0/g/masseU
 	
-!	print*,'moldiff',i_h2,i_h,ncompdiff
+! Compute the wup(ig) for H and H2 using the balistic code from R Yelle
+
+	PhiEscH=0D0
+	PhiEscH2=0D0
+
       do ig=1,ngrid
 	pp=dble(pplay(ig,:))
@@ -393,4 +410,8 @@
 	endif
 	enddo
+
+!	print*,'wi',wi(i_h),wi(i_h2),Uthermal,Lambdaexo,mmol(gcmind(:))
+!	print*,'wi',wi
+!	stop
 
 ! Compute time step for diffusion
@@ -572,4 +593,13 @@
 	CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayer,ncompdiff)
 
+! Update total escape flux of H and H2 (if q was really qnew, but not forget we will output
+! the trend only at the end
+
+	PhiEscH=PhiEscH+wi(i_h)*Nrafk(nlraf,i_h)*airephy(ig) ! in s-1
+	PhiEscH2=PhiEscH2+wi(i_h2)*Nrafk(nlraf,i_h2)*airephy(ig) ! in s-1 (U in m/s, aire in m2, Nrafk in m-3)
+!	print*,'test',ig,wi(i_h),Nrafk(nlraf,i_h),wi(i_h2),Nrafk(nlraf,i_h2),airephy(ig),PhiEscH,PhiEscH2,i_h,i_h2
+!	stop
+
+
 	if (ig .eq. ij0) then 
 	do l=il0,nlayer
@@ -623,5 +653,5 @@
 
        enddo  ! ig loop		
-
+!	print*,'Escape flux H, H2 (s-1)',PhiEscH,PhiEscH2
 
       return
@@ -1336,5 +1366,5 @@
 	Ytot=sum(Y)
 
-	if (abs((Xtot-Ytot)/Xtot) .gt. 1d-4) then
+	if (abs((Xtot-Ytot)/Xtot) .gt. 1d-3) then
 	print*,'no conservation for mass',Xtot,Ytot,nn
 	endif
@@ -1357,5 +1387,5 @@
 !	print*,'l',l,qold(l,nn),qnew(l,nn),Mold,Mnew,DM,P(l),P(l+1)
 	ENDDO
-	IF (abs(DM/Mold) .gt. 1d-5) THEN
+	IF (abs(DM/Mold) .gt. 1d-2) THEN
 	Print*,'We dont conserve mas',nn,DM,Mold,Mnew,DM/Mold
 	ENDIF