Index: trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90	(revision 1036)
+++ trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90	(revision 1047)
@@ -1,16 +1,17 @@
-subroutine moldiff_red(pplay,pplev,pt,pdt,pq,pdq,ptimestep,zzlay,pdteuv,pdtconduc,pdqdiff)
-
-use tracer_mod, only: nqmx, noms, mmol
+subroutine moldiff_red(ngrid,nlayer,nq,pplay,pplev,pt,pdt,pq,pdq,&
+                       ptimestep,zzlay,pdteuv,pdtconduc,pdqdiff)
+
+use tracer_mod, only: noms, mmol
 
 implicit none
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
-#include "callkeys.h"
-#include "comdiurn.h"
-#include "chimiedata.h"
+!#include "callkeys.h"
+!#include "comdiurn.h"
+!#include "chimiedata.h"
 !#include "tracer.h"
-#include "conc.h"
+!#include "conc.h"
 #include "diffusion.h"
 
@@ -19,15 +20,18 @@
 ! Input/Output
 !
+      integer,intent(in) :: ngrid ! number of atmospheric columns
+      integer,intent(in) :: nlayer ! number of atmospheric layers
+      integer,intent(in) :: nq ! number of advected tracers
       real ptimestep
-      real pplay(ngridmx,nlayermx)
-      real zzlay(ngridmx,nlayermx)
-      real pplev(ngridmx,nlayermx+1)
-      real pq(ngridmx,nlayermx,nqmx)
-      real pdq(ngridmx,nlayermx,nqmx)
-      real pt(ngridmx,nlayermx)
-      real pdt(ngridmx,nlayermx)
-      real pdteuv(ngridmx,nlayermx)
-      real pdtconduc(ngridmx,nlayermx)
-      real pdqdiff(ngridmx,nlayermx,nqmx)
+      real pplay(ngrid,nlayer)
+      real zzlay(ngrid,nlayer)
+      real pplev(ngrid,nlayer+1)
+      real pq(ngrid,nlayer,nq)
+      real pdq(ngrid,nlayer,nq)
+      real pt(ngrid,nlayer)
+      real pdt(ngrid,nlayer)
+      real pdteuv(ngrid,nlayer)
+      real pdtconduc(ngrid,nlayer)
+      real pdqdiff(ngrid,nlayer,nq)
 !
 ! Local
@@ -37,5 +41,5 @@
 !      real hco2(ncompdiff),ho
 
-      integer,dimension(nqmx) :: indic_diff
+      integer,dimension(nq) :: indic_diff
       integer ig,iq,nz,l,k,n,nn,p,ij0
       integer istep,il,gcn,ntime,nlraf
@@ -44,14 +48,14 @@
       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)
-      real*8 pint(nlayermx)
-      real*8 tt(nlayermx),tnew(nlayermx),tint(nlayermx)
-      real*8 zz(nlayermx)
+      real*8 hp(nlayer)
+      real*8 pp(nlayer)
+      real*8 pint(nlayer)
+      real*8 tt(nlayer),tnew(nlayer),tint(nlayer)
+      real*8 zz(nlayer)
       real*8,dimension(:,:),allocatable,save :: qq,qnew,qint,FacMass
       real*8,dimension(:,:),allocatable,save :: rhoK,rhokinit
-      real*8 rhoT(nlayermx)
-      real*8 dmmeandz(nlayermx)
-      real*8 massemoy(nlayermx)
+      real*8 rhoT(nlayer)
+      real*8 dmmeandz(nlayer)
+      real*8 massemoy(nlayer)
       real*8,dimension(:),allocatable :: Praf,Traf,Rraf,Mraf,Nraf,Draf,Hraf,Wraf
       real*8,dimension(:),allocatable :: Zraf,Tdiffraf
@@ -130,7 +134,7 @@
 
 	ncompdiff=0
-	indic_diff(1:nqmx)=0
+	indic_diff(1:nq)=0
 	
-	do nn=1,nqmx
+	do nn=1,nq
 	do n=1,14
 	if (ListeDiff(n) .eq. noms(nn)) then
@@ -151,5 +155,5 @@
 ! Store gcm indexes in gcmind
 	n=0
-	do nn=1,nqmx
+	do nn=1,nq
 	if (indic_diff(nn) .eq. 1) then
 	n=n+1
@@ -162,5 +166,5 @@
 ! find vertical index above which diffusion is computed
 
-	do l=1,nlayermx
+	do l=1,nlayer
 	if (pplay(1,l) .gt. Pdiff) then
 	il0=l
@@ -176,10 +180,10 @@
 
 ! 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(qq(nlayer,ncompdiff))
+	allocate(qnew(nlayer,ncompdiff))
+	allocate(qint(nlayer,ncompdiff))
+	allocate(FacMass(nlayer,ncompdiff))
+	allocate(rhok(nlayer,ncompdiff))
+	allocate(rhokinit(nlayer,ncompdiff))
 
 	allocate(wi(ncompdiff))
@@ -212,5 +216,5 @@
 	
 !	print*,'moldiff',i_h2,i_h,ncompdiff
-      do ig=1,ngridmx
+      do ig=1,ngrid
 	pp=dble(pplay(ig,:))
 
@@ -218,6 +222,6 @@
 
 !	CALL TMNEW(pt(ig,:),pdt(ig,:),pdtconduc(ig,:),pdteuv(ig,:)      & 
-!     &  ,tt,ptimestep,nlayermx,ig)
-        do l=1,nlayermx
+!     &  ,tt,ptimestep,nlayer,ig)
+        do l=1,nlayer
           tt(l)=pt(ig,l)*1D0+(pdt(ig,l)*dble(ptimestep)+                &
                               pdtconduc(ig,l)*dble(ptimestep)+          &
@@ -228,37 +232,37 @@
               pdt(ig,l),pdtconduc(ig,l),pdteuv(ig,l),dble(ptimestep)
           endif
-        enddo ! of do l=1,nlayermx
+        enddo ! of do l=1,nlayer
 
 ! Update the mass mixing ratios modified by other processes
 
-!	CALL QMNEW(pq(ig,:,:),pdq(ig,:,:),qq,ptimestep,nlayermx,        &
+!	CALL QMNEW(pq(ig,:,:),pdq(ig,:,:),qq,ptimestep,nlayer,        &
 !     &  ncompdiff,gcmind,ig)
         do iq=1,ncompdiff
-         do l=1,nlayermx
+         do l=1,nlayer
           qq(l,iq)=pq(ig,l,gcmind(iq))*1D0+(                            &
                            pdq(ig,l,gcmind(iq))*dble(ptimestep))
           qq(l,iq)=max(qq(l,iq),1d-30)
-         enddo ! of do l=1,nlayermx
+         enddo ! of do l=1,nlayer
         enddo ! of do iq=1,ncompdiff
 
 ! Compute the Pressure scale height
 
-	CALL HSCALE(pp,hp,nlayermx)
+	CALL HSCALE(pp,hp,nlayer)
 
 ! Compute the atmospheric mass (in Dalton)
 
-	CALL MMOY(massemoy,mmol,qq,gcmind,nlayermx,ncompdiff)
+	CALL MMOY(massemoy,mmol,qq,gcmind,nlayer,ncompdiff)
 
 ! Compute the vertical gradient of atmospheric mass
 
-	CALL DMMOY(massemoy,hp,dmmeandz,nlayermx)
+	CALL DMMOY(massemoy,hp,dmmeandz,nlayer)
 
 ! Compute the altitude of each layer
 
-	CALL ZVERT(pp,tt,massemoy,zz,nlayermx,ig)
+	CALL ZVERT(pp,tt,massemoy,zz,nlayer,ig)
 
 ! Compute the total mass density (kg/m3)
 	
-	CALL RHOTOT(pp,tt,massemoy,qq,RHOT,RHOK,nlayermx,ncompdiff)
+	CALL RHOTOT(pp,tt,massemoy,qq,RHOT,RHOK,nlayer,ncompdiff)
 	RHOKINIT=RHOK
 
@@ -271,5 +275,5 @@
 	Mtot1(1:ncompdiff)=0d0
 
-	do l=il0,nlayermx
+	do l=il0,nlayer
 	do nn=1,ncompdiff
 	Mtot1(nn)=Mtot1(nn)+1d0/g*qq(l,nn)*                             &
@@ -279,5 +283,5 @@
 
 	Zmin=zz(il0)
-	Zmax=zz(nlayermx)
+	Zmax=zz(nlayer)
 
 
@@ -286,5 +290,5 @@
 	if (Zmax .gt. 4000000.) then
 	Print*,'Zmax too high',ig,zmax,zmin
-	do l=1,nlayermx
+	do l=1,nlayer
 	print*,'old',zz(l),pt(ig,l),pdteuv(ig,l),pdq(ig,l,:)
 	print*,'l',l,rhot(l),tt(l),pp(l),massemoy(l),qq(l,:)
@@ -321,5 +325,5 @@
 	CALL UPPER_RESOL(pp,tt,zz,massemoy,RHOT,RHOK,                   &
      &  qq,mmol,gcmind,Praf,Traf,Qraf,Mraf,Zraf,                        &
-     &  Nraf,Nrafk,Rraf,Rrafk,il0,nlraf,ncompdiff,nlayermx,ig)
+     &  Nraf,Nrafk,Rraf,Rrafk,il0,nlraf,ncompdiff,nlayer,ig)
 
 	Prafold=Praf
@@ -330,9 +334,9 @@
 
 	CALL GCMGRID_P(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qint,tt,tint  &
-     &    ,pp,mmol,gcmind,nlraf,ncompdiff,nlayermx,ig)
+     &    ,pp,mmol,gcmind,nlraf,ncompdiff,nlayer,ig)
 
 ! We compute the mass correction factor of each specie at each pressure level
 
-	CALL CORRMASS(qq,qint,FacMass,nlayermx,ncompdiff)
+	CALL CORRMASS(qq,qint,FacMass,nlayer,ncompdiff)
 
 ! Altitude step
@@ -367,5 +371,5 @@
 !	enddo
 
-!	do l=1,nlayermx
+!	do l=1,nlayer
 !	print*,'l',l,zz(l),pp(l),tt(l),sum(qq(l,:)),massemoy(l)
 !	enddo
@@ -374,5 +378,5 @@
 ! No change below il0
 	
-	do l=1,nlayermx
+	do l=1,nlayer
          qnew(l,:)=qq(l,:) ! No effet below il0
        enddo
@@ -527,5 +531,5 @@
 	print*,'Mraf',Mraf
 	stop
-!	pdqdiff(1:ngridmx,1:nlayermx,1:nqmx)=0.
+!	pdqdiff(1:ngrid,1:nlayer,1:nq)=0.
 !	return
 !	Rrafk(l,nn)=1D-30*Rraf(l)
@@ -572,10 +576,10 @@
 
 	CALL GCMGRID_P2(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qnew,tt,tnew,&
-     &   pp,mmol,gcmind,nlraf,ncompdiff,nlayermx,FacMass,ig)
-
-	CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayermx,ncompdiff)
+     &   pp,mmol,gcmind,nlraf,ncompdiff,nlayer,FacMass,ig)
+
+	CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayer,ncompdiff)
 
 	if (ig .eq. ij0) then 
-	do l=il0,nlayermx
+	do l=il0,nlayer
         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,:)),&
@@ -590,5 +594,5 @@
 	Mtot2(1:ncompdiff)=0d0
 
-        do l=il0,nlayermx
+        do l=il0,nlayer
         do nn=1,ncompdiff
         Mtot2(nn)=Mtot2(nn)+1d0/g*qnew(l,nn)*                           &
@@ -600,10 +604,10 @@
 
 !	do nn=1,ncompdiff
-!	CALL CheckMass2(qq,qnew,pplev(ig,:),il0,nlayermx,nn,ncompdiff)
+!	CALL CheckMass2(qq,qnew,pplev(ig,:),il0,nlayer,nn,ncompdiff)
 !	enddo
 
 ! Compute the diffusion trends du to diffusion
 
-	do l=1,nlayermx
+	do l=1,nlayer
 	do nn=1,ncompdiff
 	pdqdiff(ig,l,gcmind(nn))=(qnew(l,nn)-qq(l,nn))/ptimestep
