Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advn.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advn.f90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advn.f90	(revision 5118)
@@ -1,6 +1,5 @@
-
 ! $Header$
 
-SUBROUTINE advn(q,masse,w,pbaru,pbarv,pdt,mode)
+SUBROUTINE advn(q, masse, w, pbaru, pbarv, pdt, mode)
   !
   ! Auteur : F. Hourdin
@@ -15,4 +14,5 @@
   !
   !   --------------------------------------------------------------------
+  USE lmdz_iniprint, ONLY: lunout, prt_level
   IMPLICIT NONE
   !
@@ -20,5 +20,4 @@
   include "paramet.h"
   include "comgeom.h"
-  include "iniprint.h"
 
   !
@@ -26,98 +25,98 @@
   !   ----------
   INTEGER :: mode
-  REAL :: masse(ip1jmp1,llm)
-  REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
-  REAL :: q(ip1jmp1,llm)
-  REAL :: w(ip1jmp1,llm),pdt
+  REAL :: masse(ip1jmp1, llm)
+  REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
+  REAL :: q(ip1jmp1, llm)
+  REAL :: w(ip1jmp1, llm), pdt
   !
   !  Local
   !   ---------
   !
-  INTEGER :: i,ij,l,j,ii
-  INTEGER :: ijlqmin,iqmin,jqmin,lqmin
+  INTEGER :: i, ij, l, j, ii
+  INTEGER :: ijlqmin, iqmin, jqmin, lqmin
   INTEGER :: ismin
   !
-  REAL :: zm(ip1jmp1,llm),newmasse
-  REAL :: mu(ip1jmp1,llm)
-  REAL :: mv(ip1jm,llm)
-  REAL :: mw(ip1jmp1,llm+1)
-  REAL :: zq(ip1jmp1,llm),zz,qpn,qps
-  REAL :: zqg(ip1jmp1,llm),zqd(ip1jmp1,llm)
-  REAL :: zqs(ip1jmp1,llm),zqn(ip1jmp1,llm)
-  REAL :: zqh(ip1jmp1,llm),zqb(ip1jmp1,llm)
-  REAL :: temps0,temps1,temps2,temps3
-  REAL :: ztemps1,ztemps2,ssum
-  save temps1,temps2,temps3
-  REAL :: zzpbar,zzw
-
-  REAL :: qmin,qmax
-  data qmin,qmax/0.,1./
-  data temps1,temps2,temps3/0.,0.,0./
+  REAL :: zm(ip1jmp1, llm), newmasse
+  REAL :: mu(ip1jmp1, llm)
+  REAL :: mv(ip1jm, llm)
+  REAL :: mw(ip1jmp1, llm + 1)
+  REAL :: zq(ip1jmp1, llm), zz, qpn, qps
+  REAL :: zqg(ip1jmp1, llm), zqd(ip1jmp1, llm)
+  REAL :: zqs(ip1jmp1, llm), zqn(ip1jmp1, llm)
+  REAL :: zqh(ip1jmp1, llm), zqb(ip1jmp1, llm)
+  REAL :: temps0, temps1, temps2, temps3
+  REAL :: ztemps1, ztemps2, ssum
+  save temps1, temps2, temps3
+  REAL :: zzpbar, zzw
+
+  REAL :: qmin, qmax
+  data qmin, qmax/0., 1./
+  data temps1, temps2, temps3/0., 0., 0./
 
   zzpbar = 0.5 * pdt
-  zzw    = pdt
-
-  DO l=1,llm
-    DO ij = iip2,ip1jm
-        mu(ij,l)=pbaru(ij,l) * zzpbar
-     ENDDO
-     DO ij=1,ip1jm
-        mv(ij,l)=pbarv(ij,l) * zzpbar
-     ENDDO
-     DO ij=1,ip1jmp1
-        mw(ij,l)=w(ij,l) * zzw
-     ENDDO
+  zzw = pdt
+
+  DO l = 1, llm
+    DO ij = iip2, ip1jm
+      mu(ij, l) = pbaru(ij, l) * zzpbar
+    ENDDO
+    DO ij = 1, ip1jm
+      mv(ij, l) = pbarv(ij, l) * zzpbar
+    ENDDO
+    DO ij = 1, ip1jmp1
+      mw(ij, l) = w(ij, l) * zzw
+    ENDDO
   ENDDO
 
-  DO ij=1,ip1jmp1
-     mw(ij,llm+1)=0.
+  DO ij = 1, ip1jmp1
+    mw(ij, llm + 1) = 0.
   ENDDO
 
-  do l=1,llm
-     qpn=0.
-     qps=0.
-     do ij=1,iim
-        qpn=qpn+q(ij,l)*masse(ij,l)
-        qps=qps+q(ip1jm+ij,l)*masse(ip1jm+ij,l)
-     enddo
-     qpn=qpn/ssum(iim,masse(1,l),1)
-     qps=qps/ssum(iim,masse(ip1jm+1,l),1)
-     do ij=1,iip1
-        q(ij,l)=qpn
-        q(ip1jm+ij,l)=qps
-     enddo
-  enddo
-
-  do ij=1,ip1jmp1
-     mw(ij,llm+1)=0.
-  enddo
-  do l=1,llm
-     do ij=1,ip1jmp1
-        zq(ij,l)=q(ij,l)
-        zm(ij,l)=masse(ij,l)
-     enddo
+  do l = 1, llm
+    qpn = 0.
+    qps = 0.
+    do ij = 1, iim
+      qpn = qpn + q(ij, l) * masse(ij, l)
+      qps = qps + q(ip1jm + ij, l) * masse(ip1jm + ij, l)
+    enddo
+    qpn = qpn / ssum(iim, masse(1, l), 1)
+    qps = qps / ssum(iim, masse(ip1jm + 1, l), 1)
+    do ij = 1, iip1
+      q(ij, l) = qpn
+      q(ip1jm + ij, l) = qps
+    enddo
+  enddo
+
+  do ij = 1, ip1jmp1
+    mw(ij, llm + 1) = 0.
+  enddo
+  do l = 1, llm
+    do ij = 1, ip1jmp1
+      zq(ij, l) = q(ij, l)
+      zm(ij, l) = masse(ij, l)
+    enddo
   enddo
 
   ! CALL minmaxq(zq,qmin,qmax,'avant vlx     ')
-  CALL advnqx(zq,zqg,zqd)
-  CALL advnx(zq,zqg,zqd,zm,mu,mode)
-  CALL advnqy(zq,zqs,zqn)
-  CALL advny(zq,zqs,zqn,zm,mv)
-  CALL advnqz(zq,zqh,zqb)
-  CALL advnz(zq,zqh,zqb,zm,mw)
+  CALL advnqx(zq, zqg, zqd)
+  CALL advnx(zq, zqg, zqd, zm, mu, mode)
+  CALL advnqy(zq, zqs, zqn)
+  CALL advny(zq, zqs, zqn, zm, mv)
+  CALL advnqz(zq, zqh, zqb)
+  CALL advnz(zq, zqh, zqb, zm, mw)
   ! CALL vlz(zq,0.,zm,mw)
-  CALL advnqy(zq,zqs,zqn)
-  CALL advny(zq,zqs,zqn,zm,mv)
-  CALL advnqx(zq,zqg,zqd)
-  CALL advnx(zq,zqg,zqd,zm,mu,mode)
+  CALL advnqy(zq, zqs, zqn)
+  CALL advny(zq, zqs, zqn, zm, mv)
+  CALL advnqx(zq, zqg, zqd)
+  CALL advnx(zq, zqg, zqd, zm, mu, mode)
   ! CALL minmaxq(zq,qmin,qmax,'apres vlx     ')
 
-  do l=1,llm
-     do ij=1,ip1jmp1
-       q(ij,l)=zq(ij,l)
-     enddo
-     do ij=1,ip1jm+1,iip1
-        q(ij+iim,l)=q(ij,l)
-     enddo
+  do l = 1, llm
+    do ij = 1, ip1jmp1
+      q(ij, l) = zq(ij, l)
+    enddo
+    do ij = 1, ip1jm + 1, iip1
+      q(ij + iim, l) = q(ij, l)
+    enddo
   enddo
 
@@ -125,27 +124,27 @@
 END SUBROUTINE advn
 
-SUBROUTINE advnqx(q,qg,qd)
+SUBROUTINE advnqx(q, qg, qd)
   !
   ! Auteurs:   Calcul des valeurs de q aux point u.
   !
   !   --------------------------------------------------------------------
+  USE lmdz_iniprint, ONLY: lunout, prt_level
   IMPLICIT NONE
   !
   INCLUDE "dimensions.h"
   INCLUDE "paramet.h"
-  INCLUDE "iniprint.h"
   !
   !
   !   Arguments:
   !   ----------
-  REAL :: q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm)
+  REAL :: q(ip1jmp1, llm), qg(ip1jmp1, llm), qd(ip1jmp1, llm)
   !
   !  Local
   !   ---------
   !
-  INTEGER :: ij,l
-  !
-  REAL :: dxqu(ip1jmp1),zqu(ip1jmp1)
-  REAL :: zqmax(ip1jmp1),zqmin(ip1jmp1)
+  INTEGER :: ij, l
+  !
+  REAL :: dxqu(ip1jmp1), zqu(ip1jmp1)
+  REAL :: zqmax(ip1jmp1), zqmin(ip1jmp1)
   LOGICAL :: extremum(ip1jmp1)
 
@@ -157,112 +156,112 @@
   !   -----------------------
   IF (mode==0) THEN
-     do l=1,llm
-        do ij=1,ip1jm
-           qd(ij,l)=q(ij,l)
-           qg(ij,l)=q(ij,l)
-        enddo
-     enddo
+    do l = 1, llm
+      do ij = 1, ip1jm
+        qd(ij, l) = q(ij, l)
+        qg(ij, l) = q(ij, l)
+      enddo
+    enddo
   else
-  do l = 1, llm
-     do ij=iip2,ip1jm-1
-        dxqu(ij)=q(ij+1,l)-q(ij,l)
-        zqu(ij)=0.5*(q(ij+1,l)+q(ij,l))
-     enddo
-     do ij=iip1+iip1,ip1jm,iip1
-        dxqu(ij)=dxqu(ij-iim)
-        zqu(ij)=zqu(ij-iim)
-     enddo
-     do ij=iip2,ip1jm-1
-        zqu(ij)=zqu(ij)-dxqu(ij+1)/12.
-     enddo
-     do ij=iip1+iip1,ip1jm,iip1
-        zqu(ij)=zqu(ij-iim)
-     enddo
-     do ij=iip2+1,ip1jm
-        zqu(ij)=zqu(ij)+dxqu(ij-1)/12.
-     enddo
-     do ij=iip1+iip1,ip1jm,iip1
-        zqu(ij-iim)=zqu(ij)
-     enddo
-
-  !   calcul des valeurs max et min acceptees aux interfaces
-
-     do ij=iip2,ip1jm-1
-        zqmax(ij)=max(q(ij+1,l),q(ij,l))
-        zqmin(ij)=min(q(ij+1,l),q(ij,l))
-     enddo
-     do ij=iip1+iip1,ip1jm,iip1
-        zqmax(ij)=zqmax(ij-iim)
-        zqmin(ij)=zqmin(ij-iim)
-     enddo
-     do ij=iip2+1,ip1jm
-        extremum(ij)=dxqu(ij)*dxqu(ij-1)<=0.
-     enddo
-     do ij=iip1+iip1,ip1jm,iip1
-        extremum(ij-iim)=extremum(ij)
-     enddo
-     do ij=iip2,ip1jm
-        zqu(ij)=min(max(zqmin(ij),zqu(ij)),zqmax(ij))
-     enddo
-     do ij=iip2+1,ip1jm
+    do l = 1, llm
+      do ij = iip2, ip1jm - 1
+        dxqu(ij) = q(ij + 1, l) - q(ij, l)
+        zqu(ij) = 0.5 * (q(ij + 1, l) + q(ij, l))
+      enddo
+      do ij = iip1 + iip1, ip1jm, iip1
+        dxqu(ij) = dxqu(ij - iim)
+        zqu(ij) = zqu(ij - iim)
+      enddo
+      do ij = iip2, ip1jm - 1
+        zqu(ij) = zqu(ij) - dxqu(ij + 1) / 12.
+      enddo
+      do ij = iip1 + iip1, ip1jm, iip1
+        zqu(ij) = zqu(ij - iim)
+      enddo
+      do ij = iip2 + 1, ip1jm
+        zqu(ij) = zqu(ij) + dxqu(ij - 1) / 12.
+      enddo
+      do ij = iip1 + iip1, ip1jm, iip1
+        zqu(ij - iim) = zqu(ij)
+      enddo
+
+      !   calcul des valeurs max et min acceptees aux interfaces
+
+      do ij = iip2, ip1jm - 1
+        zqmax(ij) = max(q(ij + 1, l), q(ij, l))
+        zqmin(ij) = min(q(ij + 1, l), q(ij, l))
+      enddo
+      do ij = iip1 + iip1, ip1jm, iip1
+        zqmax(ij) = zqmax(ij - iim)
+        zqmin(ij) = zqmin(ij - iim)
+      enddo
+      do ij = iip2 + 1, ip1jm
+        extremum(ij) = dxqu(ij) * dxqu(ij - 1)<=0.
+      enddo
+      do ij = iip1 + iip1, ip1jm, iip1
+        extremum(ij - iim) = extremum(ij)
+      enddo
+      do ij = iip2, ip1jm
+        zqu(ij) = min(max(zqmin(ij), zqu(ij)), zqmax(ij))
+      enddo
+      do ij = iip2 + 1, ip1jm
         IF(extremum(ij)) THEN
-           qg(ij,l)=q(ij,l)
-           qd(ij,l)=q(ij,l)
+          qg(ij, l) = q(ij, l)
+          qd(ij, l) = q(ij, l)
         else
-           qd(ij,l)=zqu(ij)
-           qg(ij,l)=zqu(ij-1)
+          qd(ij, l) = zqu(ij)
+          qg(ij, l) = zqu(ij - 1)
         endif
-     enddo
-     do ij=iip1+iip1,ip1jm,iip1
-        qd(ij-iim,l)=qd(ij,l)
-        qg(ij-iim,l)=qg(ij,l)
-     enddo
-
-     goto 8888
-
-     do ij=iip2+1,ip1jm
-        IF(extremum(ij).and..not.extremum(ij-1)) &
-              qd(ij-1,l)=q(ij,l)
-     enddo
-
-     do ij=iip1+iip1,ip1jm,iip1
-        qd(ij-iim,l)=qd(ij,l)
-     enddo
-     do ij=iip2,ip1jm-1
-        IF (extremum(ij).and..not.extremum(ij+1)) &
-              qg(ij+1,l)=q(ij,l)
-     enddo
-
-     do ij=iip1+iip1,ip1jm,iip1
-        qg(ij,l)=qg(ij-iim,l)
-     enddo
-8888   continue
-  enddo
+      enddo
+      do ij = iip1 + iip1, ip1jm, iip1
+        qd(ij - iim, l) = qd(ij, l)
+        qg(ij - iim, l) = qg(ij, l)
+      enddo
+
+      goto 8888
+
+      do ij = iip2 + 1, ip1jm
+        IF(extremum(ij).and..not.extremum(ij - 1)) &
+                qd(ij - 1, l) = q(ij, l)
+      enddo
+
+      do ij = iip1 + iip1, ip1jm, iip1
+        qd(ij - iim, l) = qd(ij, l)
+      enddo
+      do ij = iip2, ip1jm - 1
+        IF (extremum(ij).and..not.extremum(ij + 1)) &
+                qg(ij + 1, l) = q(ij, l)
+      enddo
+
+      do ij = iip1 + iip1, ip1jm, iip1
+        qg(ij, l) = qg(ij - iim, l)
+      enddo
+      8888   continue
+    enddo
   ENDIF
   RETURN
 END SUBROUTINE advnqx
-SUBROUTINE advnqy(q,qs,qn)
+SUBROUTINE advnqy(q, qs, qn)
   !
   ! Auteurs:   Calcul des valeurs de q aux point v.
   !
   !   --------------------------------------------------------------------
+  USE lmdz_iniprint, ONLY: lunout, prt_level
   IMPLICIT NONE
   !
   INCLUDE "dimensions.h"
   INCLUDE "paramet.h"
-  INCLUDE "iniprint.h"
   !
   !
   !   Arguments:
   !   ----------
-  REAL :: q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm)
+  REAL :: q(ip1jmp1, llm), qs(ip1jmp1, llm), qn(ip1jmp1, llm)
   !
   !  Local
   !   ---------
   !
-  INTEGER :: ij,l
-  !
-  REAL :: dyqv(ip1jm),zqv(ip1jm,llm)
-  REAL :: zqmax(ip1jm),zqmin(ip1jm)
+  INTEGER :: ij, l
+  !
+  REAL :: dyqv(ip1jm), zqv(ip1jm, llm)
+  REAL :: zqmax(ip1jm), zqmin(ip1jm)
   LOGICAL :: extremum(ip1jmp1)
 
@@ -272,94 +271,94 @@
 
   IF (mode==0) THEN
-     do l=1,llm
-        do ij=1,ip1jmp1
-           qn(ij,l)=q(ij,l)
-           qs(ij,l)=q(ij,l)
-        enddo
-     enddo
+    do l = 1, llm
+      do ij = 1, ip1jmp1
+        qn(ij, l) = q(ij, l)
+        qs(ij, l) = q(ij, l)
+      enddo
+    enddo
   else
 
-  !   calcul des pentes en u:
-  !   -----------------------
-  do l = 1, llm
-     do ij=1,ip1jm
-        dyqv(ij)=q(ij,l)-q(ij+iip1,l)
-     enddo
-
-     do ij=iip2,ip1jm-iip1
-        zqv(ij,l)=0.5*(q(ij+iip1,l)+q(ij,l))
-        zqv(ij,l)=zqv(ij,l)+(dyqv(ij+iip1)-dyqv(ij-iip1))/12.
-     enddo
-
-     do ij=iip2,ip1jm
-        extremum(ij)=dyqv(ij)*dyqv(ij-iip1)<=0.
-     enddo
-
-  ! Pas de pentes aux poles
-     do ij=1,iip1
-        zqv(ij,l)=q(ij,l)
-        zqv(ip1jm-iip1+ij,l)=q(ip1jm+ij,l)
-        extremum(ij)=.TRUE.
-        extremum(ip1jmp1-iip1+ij)=.TRUE.
-     enddo
-
-  !   calcul des valeurs max et min acceptees aux interfaces
-     do ij=1,ip1jm
-        zqmax(ij)=max(q(ij+iip1,l),q(ij,l))
-        zqmin(ij)=min(q(ij+iip1,l),q(ij,l))
-     enddo
-
-     do ij=1,ip1jm
-        zqv(ij,l)=min(max(zqmin(ij),zqv(ij,l)),zqmax(ij))
-     enddo
-
-     do ij=iip2,ip1jm
+    !   calcul des pentes en u:
+    !   -----------------------
+    do l = 1, llm
+      do ij = 1, ip1jm
+        dyqv(ij) = q(ij, l) - q(ij + iip1, l)
+      enddo
+
+      do ij = iip2, ip1jm - iip1
+        zqv(ij, l) = 0.5 * (q(ij + iip1, l) + q(ij, l))
+        zqv(ij, l) = zqv(ij, l) + (dyqv(ij + iip1) - dyqv(ij - iip1)) / 12.
+      enddo
+
+      do ij = iip2, ip1jm
+        extremum(ij) = dyqv(ij) * dyqv(ij - iip1)<=0.
+      enddo
+
+      ! Pas de pentes aux poles
+      do ij = 1, iip1
+        zqv(ij, l) = q(ij, l)
+        zqv(ip1jm - iip1 + ij, l) = q(ip1jm + ij, l)
+        extremum(ij) = .TRUE.
+        extremum(ip1jmp1 - iip1 + ij) = .TRUE.
+      enddo
+
+      !   calcul des valeurs max et min acceptees aux interfaces
+      do ij = 1, ip1jm
+        zqmax(ij) = max(q(ij + iip1, l), q(ij, l))
+        zqmin(ij) = min(q(ij + iip1, l), q(ij, l))
+      enddo
+
+      do ij = 1, ip1jm
+        zqv(ij, l) = min(max(zqmin(ij), zqv(ij, l)), zqmax(ij))
+      enddo
+
+      do ij = iip2, ip1jm
         IF(extremum(ij)) THEN
-           qs(ij,l)=q(ij,l)
-           qn(ij,l)=q(ij,l)
-           ! if (.NOT.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)
-           ! if (.NOT.extremum(ij+iip1)) qn(ij+iip1,l)=q(ij,l)
+          qs(ij, l) = q(ij, l)
+          qn(ij, l) = q(ij, l)
+          ! if (.NOT.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)
+          ! if (.NOT.extremum(ij+iip1)) qn(ij+iip1,l)=q(ij,l)
         else
-           qs(ij,l)=zqv(ij,l)
-           qn(ij,l)=zqv(ij-iip1,l)
+          qs(ij, l) = zqv(ij, l)
+          qn(ij, l) = zqv(ij - iip1, l)
         endif
-     enddo
-
-     do ij=1,iip1
-        qs(ij,l)=q(ij,l)
-        qn(ij,l)=q(ij,l)
-        qs(ip1jm+ij,l)=q(ip1jm+ij,l)
-        qn(ip1jm+ij,l)=q(ip1jm+ij,l)
-     enddo
-
-  enddo
+      enddo
+
+      do ij = 1, iip1
+        qs(ij, l) = q(ij, l)
+        qn(ij, l) = q(ij, l)
+        qs(ip1jm + ij, l) = q(ip1jm + ij, l)
+        qn(ip1jm + ij, l) = q(ip1jm + ij, l)
+      enddo
+
+    enddo
   ENDIF
   RETURN
 END SUBROUTINE advnqy
 
-SUBROUTINE advnqz(q,qh,qb)
+SUBROUTINE advnqz(q, qh, qb)
   !
   ! Auteurs:   Calcul des valeurs de q aux point v.
   !
   !   --------------------------------------------------------------------
+  USE lmdz_iniprint, ONLY: lunout, prt_level
   IMPLICIT NONE
   !
   INCLUDE "dimensions.h"
   INCLUDE "paramet.h"
-  INCLUDE "iniprint.h"
   !
   !
   !   Arguments:
   !   ----------
-  REAL :: q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm)
+  REAL :: q(ip1jmp1, llm), qh(ip1jmp1, llm), qb(ip1jmp1, llm)
   !
   !  Local
   !   ---------
   !
-  INTEGER :: ij,l
-  !
-  REAL :: dzqw(ip1jmp1,llm+1),zqw(ip1jmp1,llm+1)
-  REAL :: zqmax(ip1jmp1,llm),zqmin(ip1jmp1,llm)
-  LOGICAL :: extremum(ip1jmp1,llm)
+  INTEGER :: ij, l
+  !
+  REAL :: dzqw(ip1jmp1, llm + 1), zqw(ip1jmp1, llm + 1)
+  REAL :: zqmax(ip1jmp1, llm), zqmin(ip1jmp1, llm)
+  LOGICAL :: extremum(ip1jmp1, llm)
 
   INTEGER :: mode
@@ -372,80 +371,80 @@
 
   IF (mode==0) THEN
-     do l=1,llm
-        do ij=1,ip1jmp1
-           qb(ij,l)=q(ij,l)
-           qh(ij,l)=q(ij,l)
-        enddo
-     enddo
+    do l = 1, llm
+      do ij = 1, ip1jmp1
+        qb(ij, l) = q(ij, l)
+        qh(ij, l) = q(ij, l)
+      enddo
+    enddo
   else
-  do l = 2, llm
-     do ij=1,ip1jmp1
-        dzqw(ij,l)=q(ij,l-1)-q(ij,l)
-        zqw(ij,l)=0.5*(q(ij,l-1)+q(ij,l))
-     enddo
-  enddo
-  do ij=1,ip1jmp1
-     dzqw(ij,1)=0.
-     dzqw(ij,llm+1)=0.
-  enddo
-  do l=2,llm
-     do ij=1,ip1jmp1
-        zqw(ij,l)=zqw(ij,l)+(dzqw(ij,l+1)-dzqw(ij,l-1))/12.
-     enddo
-  enddo
-  do l=2,llm-1
-     do ij=1,ip1jmp1
-        extremum(ij,l)=dzqw(ij,l)*dzqw(ij,l+1)<=0.
-     enddo
-  enddo
-
-  ! Pas de pentes en bas et en haut
-     do ij=1,ip1jmp1
-        zqw(ij,2)=q(ij,1)
-        zqw(ij,llm)=q(ij,llm)
-        extremum(ij,1)=.TRUE.
-        extremum(ij,llm)=.TRUE.
-     enddo
-
-  !   calcul des valeurs max et min acceptees aux interfaces
-  do l=2,llm
-     do ij=1,ip1jmp1
-        zqmax(ij,l)=max(q(ij,l-1),q(ij,l))
-        zqmin(ij,l)=min(q(ij,l-1),q(ij,l))
-     enddo
-  enddo
-
-  do l=2,llm
-     do ij=1,ip1jmp1
-        zqw(ij,l)=min(max(zqmin(ij,l),zqw(ij,l)),zqmax(ij,l))
-     enddo
-  enddo
-
-  do l=2,llm-1
-     do ij=1,ip1jmp1
-        IF(extremum(ij,l)) THEN
-           qh(ij,l)=q(ij,l)
-           qb(ij,l)=q(ij,l)
+    do l = 2, llm
+      do ij = 1, ip1jmp1
+        dzqw(ij, l) = q(ij, l - 1) - q(ij, l)
+        zqw(ij, l) = 0.5 * (q(ij, l - 1) + q(ij, l))
+      enddo
+    enddo
+    do ij = 1, ip1jmp1
+      dzqw(ij, 1) = 0.
+      dzqw(ij, llm + 1) = 0.
+    enddo
+    do l = 2, llm
+      do ij = 1, ip1jmp1
+        zqw(ij, l) = zqw(ij, l) + (dzqw(ij, l + 1) - dzqw(ij, l - 1)) / 12.
+      enddo
+    enddo
+    do l = 2, llm - 1
+      do ij = 1, ip1jmp1
+        extremum(ij, l) = dzqw(ij, l) * dzqw(ij, l + 1)<=0.
+      enddo
+    enddo
+
+    ! Pas de pentes en bas et en haut
+    do ij = 1, ip1jmp1
+      zqw(ij, 2) = q(ij, 1)
+      zqw(ij, llm) = q(ij, llm)
+      extremum(ij, 1) = .TRUE.
+      extremum(ij, llm) = .TRUE.
+    enddo
+
+    !   calcul des valeurs max et min acceptees aux interfaces
+    do l = 2, llm
+      do ij = 1, ip1jmp1
+        zqmax(ij, l) = max(q(ij, l - 1), q(ij, l))
+        zqmin(ij, l) = min(q(ij, l - 1), q(ij, l))
+      enddo
+    enddo
+
+    do l = 2, llm
+      do ij = 1, ip1jmp1
+        zqw(ij, l) = min(max(zqmin(ij, l), zqw(ij, l)), zqmax(ij, l))
+      enddo
+    enddo
+
+    do l = 2, llm - 1
+      do ij = 1, ip1jmp1
+        IF(extremum(ij, l)) THEN
+          qh(ij, l) = q(ij, l)
+          qb(ij, l) = q(ij, l)
         else
-           qh(ij,l)=zqw(ij,l+1)
-           qb(ij,l)=zqw(ij,l)
+          qh(ij, l) = zqw(ij, l + 1)
+          qb(ij, l) = zqw(ij, l)
         endif
-     enddo
-  enddo
-  ! do l=2,llm-1
-  !    do ij=1,ip1jmp1
-  !       IF(extremum(ij,l)) THEN
-  !          if (.NOT.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)
-  !          if (.NOT.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)
-  !       endif
-  !    enddo
-  ! enddo
-
-  do ij=1,ip1jmp1
-     qb(ij,1)=q(ij,1)
-     qh(ij,1)=q(ij,1)
-     qb(ij,llm)=q(ij,llm)
-     qh(ij,llm)=q(ij,llm)
-  enddo
+      enddo
+    enddo
+    ! do l=2,llm-1
+    !    do ij=1,ip1jmp1
+    !       IF(extremum(ij,l)) THEN
+    !          if (.NOT.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)
+    !          if (.NOT.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)
+    !       endif
+    !    enddo
+    ! enddo
+
+    do ij = 1, ip1jmp1
+      qb(ij, 1) = q(ij, 1)
+      qh(ij, 1) = q(ij, 1)
+      qb(ij, llm) = q(ij, llm)
+      qh(ij, llm) = q(ij, llm)
+    enddo
 
   ENDIF
@@ -454,5 +453,5 @@
 END SUBROUTINE advnqz
 
-SUBROUTINE advnx(q,qg,qd,masse,u_m,mode)
+SUBROUTINE advnx(q, qg, qd, masse, u_m, mode)
   !
   ! Auteur : F. Hourdin
@@ -465,9 +464,9 @@
   !
   !   --------------------------------------------------------------------
+  USE lmdz_iniprint, ONLY: lunout, prt_level
   IMPLICIT NONE
   !
   include "dimensions.h"
   include "paramet.h"
-  include "iniprint.h"
   !
   !
@@ -475,21 +474,21 @@
   !   ----------
   INTEGER :: mode
-  REAL :: masse(ip1jmp1,llm)
-  REAL :: u_m( ip1jmp1,llm )
-  REAL :: q(ip1jmp1,llm),qd(ip1jmp1,llm),qg(ip1jmp1,llm)
+  REAL :: masse(ip1jmp1, llm)
+  REAL :: u_m(ip1jmp1, llm)
+  REAL :: q(ip1jmp1, llm), qd(ip1jmp1, llm), qg(ip1jmp1, llm)
   !
   !  Local
   !   ---------
   !
-  INTEGER :: i,j,ij,l,indu(ip1jmp1),niju,iju,ijq
-  INTEGER :: n0,nl(llm)
-  !
-  REAL :: new_m,zu_m,zdq,zz
-  REAL :: zsigg(ip1jmp1,llm),zsigd(ip1jmp1,llm),zsig
-  REAL :: u_mq(ip1jmp1,llm)
-
-  REAL :: zm,zq,zsigm,zsigp,zqm,zqp,zu
-
-  LOGICAL :: ladvplus(ip1jmp1,llm)
+  INTEGER :: i, j, ij, l, indu(ip1jmp1), niju, iju, ijq
+  INTEGER :: n0, nl(llm)
+  !
+  REAL :: new_m, zu_m, zdq, zz
+  REAL :: zsigg(ip1jmp1, llm), zsigd(ip1jmp1, llm), zsig
+  REAL :: u_mq(ip1jmp1, llm)
+
+  REAL :: zm, zq, zsigm, zsigp, zqm, zqp, zu
+
+  LOGICAL :: ladvplus(ip1jmp1, llm)
 
   REAL :: prec
@@ -498,70 +497,70 @@
   data prec/1.e-15/
 
-  do l=1,llm
-        do ij=iip2,ip1jm
-           zdq=qd(ij,l)-qg(ij,l)
-           ! if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) THEN
-           !    PRINT*,'probleme au point ij=',ij,'  l=',l
-           !    PRINT*,qd(ij,l),q(ij,l),qg(ij,l)
-           !    qd(ij,l)=q(ij,l)
-           !    qg(ij,l)=q(ij,l)
-           ! END IF
-           IF(abs(zdq)>prec) THEN
-              zsigd(ij,l)=(q(ij,l)-qg(ij,l))/zdq
-              zsigg(ij,l)=1.-zsigd(ij,l)
-              ! IF(.NOT.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .AND.
-  !    s               zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) THEN
-              !    PRINT*,'probleme au point ij=',ij,'  l=',l
-              !    PRINT*,'sigg=',zsigg(ij,l),'  sigd=',zsigd(ij,l)
-              !    PRINT*,'q d,c,g ',qd(ij,l),q(ij,l),qg(ij,l),zdq
-              !    stop
-              ! END IF
-           else
-              zsigd(ij,l)=0.5
-              zsigg(ij,l)=0.5
-              qd(ij,l)=q(ij,l)
-              qg(ij,l)=q(ij,l)
-           endif
-        enddo
-   enddo
+  do l = 1, llm
+    do ij = iip2, ip1jm
+      zdq = qd(ij, l) - qg(ij, l)
+      ! if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) THEN
+      !    PRINT*,'probleme au point ij=',ij,'  l=',l
+      !    PRINT*,qd(ij,l),q(ij,l),qg(ij,l)
+      !    qd(ij,l)=q(ij,l)
+      !    qg(ij,l)=q(ij,l)
+      ! END IF
+      IF(abs(zdq)>prec) THEN
+        zsigd(ij, l) = (q(ij, l) - qg(ij, l)) / zdq
+        zsigg(ij, l) = 1. - zsigd(ij, l)
+        ! IF(.NOT.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .AND.
+        !    s               zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) THEN
+        !    PRINT*,'probleme au point ij=',ij,'  l=',l
+        !    PRINT*,'sigg=',zsigg(ij,l),'  sigd=',zsigd(ij,l)
+        !    PRINT*,'q d,c,g ',qd(ij,l),q(ij,l),qg(ij,l),zdq
+        !    stop
+        ! END IF
+      else
+        zsigd(ij, l) = 0.5
+        zsigg(ij, l) = 0.5
+        qd(ij, l) = q(ij, l)
+        qg(ij, l) = q(ij, l)
+      endif
+    enddo
+  enddo
 
   !   calcul de la pente maximum dans la maille en valeur absolue
 
-   do l=1,llm
-   do ij=iip2,ip1jm-1
-      IF (u_m(ij,l)>=0.) THEN
-         zsigp=zsigd(ij,l)
-         zsigm=zsigg(ij,l)
-         zqp=qd(ij,l)
-         zqm=qg(ij,l)
-         zm=masse(ij,l)
-         zq=q(ij,l)
+  do l = 1, llm
+    do ij = iip2, ip1jm - 1
+      IF (u_m(ij, l)>=0.) THEN
+        zsigp = zsigd(ij, l)
+        zsigm = zsigg(ij, l)
+        zqp = qd(ij, l)
+        zqm = qg(ij, l)
+        zm = masse(ij, l)
+        zq = q(ij, l)
       else
-         zsigm=zsigd(ij+1,l)
-         zsigp=zsigg(ij+1,l)
-         zqm=qd(ij+1,l)
-         zqp=qg(ij+1,l)
-         zm=masse(ij+1,l)
-         zq=q(ij+1,l)
-      endif
-      zu=abs(u_m(ij,l))
-      ladvplus(ij,l)=zu>zm
-      zsig=zu/zm
-      IF(zsig==0.) zsigp=0.1
+        zsigm = zsigd(ij + 1, l)
+        zsigp = zsigg(ij + 1, l)
+        zqm = qd(ij + 1, l)
+        zqp = qg(ij + 1, l)
+        zm = masse(ij + 1, l)
+        zq = q(ij + 1, l)
+      endif
+      zu = abs(u_m(ij, l))
+      ladvplus(ij, l) = zu>zm
+      zsig = zu / zm
+      IF(zsig==0.) zsigp = 0.1
       IF (mode==1) THEN
-         IF (zsig<=zsigp) THEN
-             u_mq(ij,l)=u_m(ij,l)*zqp
-         ELSE IF (mode==1) THEN
-             u_mq(ij,l)= &
-                   sign(zm,u_m(ij,l))*(zsigp*zqp+(zsig-zsigp)*zqm)
-         endif
+        IF (zsig<=zsigp) THEN
+          u_mq(ij, l) = u_m(ij, l) * zqp
+        ELSE IF (mode==1) THEN
+          u_mq(ij, l) = &
+                  sign(zm, u_m(ij, l)) * (zsigp * zqp + (zsig - zsigp) * zqm)
+        endif
       else
-         IF (zsig<=zsigp) THEN
-             u_mq(ij,l)=u_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
-         else
-            zz=0.5*(zsig-zsigp)/zsigm
-            u_mq(ij,l)=sign(zm,u_m(ij,l))*( 0.5*(zq+zqp)*zsigp &
-                  +(zsig-zsigp)*(zq+zz*(zqm-zq)) )
-         endif
+        IF (zsig<=zsigp) THEN
+          u_mq(ij, l) = u_m(ij, l) * (zqp - 0.5 * zsig / zsigp * (zqp - zq))
+        else
+          zz = 0.5 * (zsig - zsigp) / zsigm
+          u_mq(ij, l) = sign(zm, u_m(ij, l)) * (0.5 * (zq + zqp) * zsigp &
+                  + (zsig - zsigp) * (zq + zz * (zqm - zq)))
+        endif
       endif
       ! IF(zsig.lt.0.) THEN
@@ -569,12 +568,12 @@
       !    stop
       ! END IF
-  enddo
-  enddo
-
-  do l=1,llm
-   do ij=iip1+iip1,ip1jm,iip1
-      u_mq(ij,l)=u_mq(ij-iim,l)
-      ladvplus(ij,l)=ladvplus(ij-iim,l)
-   enddo
+    enddo
+  enddo
+
+  do l = 1, llm
+    do ij = iip1 + iip1, ip1jm, iip1
+      u_mq(ij, l) = u_mq(ij - iim, l)
+      ladvplus(ij, l) = ladvplus(ij - iim, l)
+    enddo
   enddo
 
@@ -583,112 +582,112 @@
   !=================================================================
   !   tris des regions a traiter
-  n0=0
-  do l=1,llm
-     nl(l)=0
-     do ij=iip2,ip1jm
-        IF(ladvplus(ij,l)) THEN
-           nl(l)=nl(l)+1
-           u_mq(ij,l)=0.
-        endif
-     enddo
-     n0=n0+nl(l)
+  n0 = 0
+  do l = 1, llm
+    nl(l) = 0
+    do ij = iip2, ip1jm
+      IF(ladvplus(ij, l)) THEN
+        nl(l) = nl(l) + 1
+        u_mq(ij, l) = 0.
+      endif
+    enddo
+    n0 = n0 + nl(l)
   enddo
 
   IF(n0>1) THEN
-  IF (prt_level > 9) WRITE(lunout,*) &
-        'Nombre de points pour lesquels on advect plus que le' &
-        ,'contenu de la maille : ',n0
-
-     do l=1,llm
-        IF(nl(l)>0) THEN
-           iju=0
-  !   indicage des mailles concernees par le traitement special
-           do ij=iip2,ip1jm
-              IF(ladvplus(ij,l).AND.mod(ij,iip1)/=0) THEN
-                 iju=iju+1
-                 indu(iju)=ij
+    IF (prt_level > 9) WRITE(lunout, *) &
+            'Nombre de points pour lesquels on advect plus que le' &
+            , 'contenu de la maille : ', n0
+
+    do l = 1, llm
+      IF(nl(l)>0) THEN
+        iju = 0
+        !   indicage des mailles concernees par le traitement special
+        do ij = iip2, ip1jm
+          IF(ladvplus(ij, l).AND.mod(ij, iip1)/=0) THEN
+            iju = iju + 1
+            indu(iju) = ij
+          endif
+        enddo
+        niju = iju
+        ! PRINT*,'niju,nl',niju,nl(l)
+
+        !  traitement des mailles
+        do iju = 1, niju
+          ij = indu(iju)
+          j = (ij - 1) / iip1 + 1
+          zu_m = u_m(ij, l)
+          u_mq(ij, l) = 0.
+          IF(zu_m>0.) THEN
+            ijq = ij
+            i = ijq - (j - 1) * iip1
+            !   accumulation pour les mailles completements advectees
+            do while(zu_m>masse(ijq, l))
+              u_mq(ij, l) = u_mq(ij, l) + q(ijq, l) * masse(ijq, l)
+              zu_m = zu_m - masse(ijq, l)
+              i = mod(i - 2 + iim, iim) + 1
+              ijq = (j - 1) * iip1 + i
+            enddo
+            !   MODIFS SPECIFIQUES DU SCHEMA
+            !   ajout de la maille non completement advectee
+            zsig = zu_m / masse(ijq, l)
+            IF(zsig<=zsigd(ijq, l)) THEN
+              u_mq(ij, l) = u_mq(ij, l) + zu_m * (qd(ijq, l) &
+                      - 0.5 * zsig / zsigd(ijq, l) * (qd(ijq, l) - q(ijq, l)))
+            else
+              ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
+              ! goto 8888
+              zz = 0.5 * (zsig - zsigd(ijq, l)) / zsigg(ijq, l)
+              IF(.NOT.(zz>0..and.zz<=0.5)) THEN
+                WRITE(lunout, *)'probleme2 au point ij=', ij, &
+                        '  l=', l
+                WRITE(lunout, *)'zz=', zz
+                stop
               endif
-           enddo
-           niju=iju
-           ! PRINT*,'niju,nl',niju,nl(l)
-
-  !  traitement des mailles
-           do iju=1,niju
-              ij=indu(iju)
-              j=(ij-1)/iip1+1
-              zu_m=u_m(ij,l)
-              u_mq(ij,l)=0.
-              IF(zu_m>0.) THEN
-                 ijq=ij
-                 i=ijq-(j-1)*iip1
-  !   accumulation pour les mailles completements advectees
-                 do while(zu_m>masse(ijq,l))
-                    u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
-                    zu_m=zu_m-masse(ijq,l)
-                    i=mod(i-2+iim,iim)+1
-                    ijq=(j-1)*iip1+i
-                 enddo
-  !   MODIFS SPECIFIQUES DU SCHEMA
-  !   ajout de la maille non completement advectee
-         zsig=zu_m/masse(ijq,l)
-         IF(zsig<=zsigd(ijq,l)) THEN
-            u_mq(ij,l)=u_mq(ij,l)+zu_m*(qd(ijq,l) &
-                  -0.5*zsig/zsigd(ijq,l)*(qd(ijq,l)-q(ijq,l)))
-         else
-            ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
-      ! goto 8888
-            zz=0.5*(zsig-zsigd(ijq,l))/zsigg(ijq,l)
-            IF(.NOT.(zz>0..and.zz<=0.5)) THEN
-                 WRITE(lunout,*)'probleme2 au point ij=',ij, &
-                       '  l=',l
-                 WRITE(lunout,*)'zz=',zz
-                 stop
+              u_mq(ij, l) = u_mq(ij, l) + masse(ijq, l) * (&
+                      0.5 * (q(ijq, l) + qd(ijq, l)) * zsigd(ijq, l) &
+                              + (zsig - zsigd(ijq, l)) * (q(ijq, l) + zz * (qg(ijq, l) - q(ijq, l))))
             endif
-            u_mq(ij,l)=u_mq(ij,l)+masse(ijq,l)*( &
-                  0.5*(q(ijq,l)+qd(ijq,l))*zsigd(ijq,l) &
-                  +(zsig-zsigd(ijq,l))*(q(ijq,l)+zz*(qg(ijq,l)-q(ijq,l))) )
-         endif
-              else
-                 ijq=ij+1
-                 i=ijq-(j-1)*iip1
-  !   accumulation pour les mailles completements advectees
-                 do while(-zu_m>masse(ijq,l))
-                    u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
-                    zu_m=zu_m+masse(ijq,l)
-                    i=mod(i,iim)+1
-                    ijq=(j-1)*iip1+i
-                 enddo
-  !   ajout de la maille non completement advectee
-  ! 2eme MODIF SPECIFIQUE
-         zsig=-zu_m/masse(ij+1,l)
-         IF(zsig<=zsigg(ijq,l)) THEN
-            u_mq(ij,l)=u_mq(ij,l)+zu_m*(qg(ijq,l) &
-                  -0.5*zsig/zsigg(ijq,l)*(qg(ijq,l)-q(ijq,l)))
-         else
-            ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
-        ! goto 9999
-            zz=0.5*(zsig-zsigg(ijq,l))/zsigd(ijq,l)
-            IF(.NOT.(zz>0..and.zz<=0.5)) THEN
-                 WRITE(lunout,*)'probleme22 au point ij=',ij &
-                       ,'  l=',l
-                 WRITE(lunout,*)'zz=',zz
-                 stop
+          else
+            ijq = ij + 1
+            i = ijq - (j - 1) * iip1
+            !   accumulation pour les mailles completements advectees
+            do while(-zu_m>masse(ijq, l))
+              u_mq(ij, l) = u_mq(ij, l) - q(ijq, l) * masse(ijq, l)
+              zu_m = zu_m + masse(ijq, l)
+              i = mod(i, iim) + 1
+              ijq = (j - 1) * iip1 + i
+            enddo
+            !   ajout de la maille non completement advectee
+            ! 2eme MODIF SPECIFIQUE
+            zsig = -zu_m / masse(ij + 1, l)
+            IF(zsig<=zsigg(ijq, l)) THEN
+              u_mq(ij, l) = u_mq(ij, l) + zu_m * (qg(ijq, l) &
+                      - 0.5 * zsig / zsigg(ijq, l) * (qg(ijq, l) - q(ijq, l)))
+            else
+              ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
+              ! goto 9999
+              zz = 0.5 * (zsig - zsigg(ijq, l)) / zsigd(ijq, l)
+              IF(.NOT.(zz>0..and.zz<=0.5)) THEN
+                WRITE(lunout, *)'probleme22 au point ij=', ij &
+                        , '  l=', l
+                WRITE(lunout, *)'zz=', zz
+                stop
+              endif
+              u_mq(ij, l) = u_mq(ij, l) - masse(ijq, l) * (&
+                      0.5 * (q(ijq, l) + qg(ijq, l)) * zsigg(ijq, l) &
+                              + (zsig - zsigg(ijq, l)) * &
+                              (q(ijq, l) + zz * (qd(ijq, l) - q(ijq, l))))
             endif
-            u_mq(ij,l)=u_mq(ij,l)-masse(ijq,l)*( &
-                  0.5*(q(ijq,l)+qg(ijq,l))*zsigg(ijq,l) &
-                  +(zsig-zsigg(ijq,l))* &
-                  (q(ijq,l)+zz*(qd(ijq,l)-q(ijq,l))) )
-         endif
-  !   fin de la modif
-              endif
-           enddo
-        endif
-     enddo
+            !   fin de la modif
+          endif
+        enddo
+      endif
+    enddo
   ENDIF  ! n0.gt.0
 
   !   bouclage en latitude
-  do l=1,llm
-    do ij=iip1+iip1,ip1jm,iip1
-       u_mq(ij,l)=u_mq(ij-iim,l)
+  do l = 1, llm
+    do ij = iip1 + iip1, ip1jm, iip1
+      u_mq(ij, l) = u_mq(ij - iim, l)
     enddo
   enddo
@@ -698,22 +697,22 @@
   !=================================================================
 
-  do l=1,llm
-     do ij=iip2+1,ip1jm
-        new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
-        q(ij,l)=(q(ij,l)*masse(ij,l)+ &
-              u_mq(ij-1,l)-u_mq(ij,l)) &
-              /new_m
-        masse(ij,l)=new_m
-     enddo
-  !   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
-     do ij=iip1+iip1,ip1jm,iip1
-        q(ij-iim,l)=q(ij,l)
-        masse(ij-iim,l)=masse(ij,l)
-     enddo
+  do l = 1, llm
+    do ij = iip2 + 1, ip1jm
+      new_m = masse(ij, l) + u_m(ij - 1, l) - u_m(ij, l)
+      q(ij, l) = (q(ij, l) * masse(ij, l) + &
+              u_mq(ij - 1, l) - u_mq(ij, l)) &
+              / new_m
+      masse(ij, l) = new_m
+    enddo
+    !   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
+    do ij = iip1 + iip1, ip1jm, iip1
+      q(ij - iim, l) = q(ij, l)
+      masse(ij - iim, l) = masse(ij, l)
+    enddo
   enddo
 
   RETURN
 END SUBROUTINE advnx
-SUBROUTINE advny(q,qs,qn,masse,v_m)
+SUBROUTINE advny(q, qs, qn, masse, v_m)
   !
   ! Auteur : F. Hourdin
@@ -726,4 +725,5 @@
   !
   !   --------------------------------------------------------------------
+  USE lmdz_iniprint, ONLY: lunout, prt_level
   IMPLICIT NONE
   !
@@ -731,23 +731,22 @@
   INCLUDE "paramet.h"
   INCLUDE "comgeom.h"
-  INCLUDE "iniprint.h"
   !
   !
   !   Arguments:
   !   ----------
-  REAL :: masse(ip1jmp1,llm)
-  REAL :: v_m( ip1jm,llm )
-  REAL :: q(ip1jmp1,llm),qn(ip1jmp1,llm),qs(ip1jmp1,llm)
+  REAL :: masse(ip1jmp1, llm)
+  REAL :: v_m(ip1jm, llm)
+  REAL :: q(ip1jmp1, llm), qn(ip1jmp1, llm), qs(ip1jmp1, llm)
   !
   !  Local
   !   ---------
   !
-  INTEGER :: ij,l
-  !
-  REAL :: new_m,zdq,zz
-  REAL :: zsigs(ip1jmp1),zsign(ip1jmp1),zsig
-  REAL :: v_mq(ip1jm,llm)
-  REAL :: convpn,convps,convmpn,convmps,massen,masses
-  REAL :: zm,zq,zsigm,zsigp,zqm,zqp
+  INTEGER :: ij, l
+  !
+  REAL :: new_m, zdq, zz
+  REAL :: zsigs(ip1jmp1), zsign(ip1jmp1), zsig
+  REAL :: v_mq(ip1jm, llm)
+  REAL :: convpn, convps, convmpn, convmps, massen, masses
+  REAL :: zm, zq, zsigm, zsigp, zqm, zqp
   REAL :: ssum
   REAL :: prec
@@ -755,91 +754,91 @@
 
   data prec/1.e-15/
-  do l=1,llm
-        do ij=1,ip1jmp1
-           zdq=qn(ij,l)-qs(ij,l)
-           ! if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) THEN
-           !    PRINT*,'probleme au point ij=',ij,'  l=',l,'  advnqx'
-           !    PRINT*,qn(ij,l),q(ij,l),qs(ij,l)
-           !    qn(ij,l)=q(ij,l)
-           !    qs(ij,l)=q(ij,l)
-           ! END IF
-           IF(abs(zdq)>prec) THEN
-              zsign(ij)=(q(ij,l)-qs(ij,l))/zdq
-              zsigs(ij)=1.-zsign(ij)
-              ! IF(.NOT.(zsign(ij).ge.0..and.zsign(ij).le.1. .AND.
-  !    s               zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) THEN
-              !    PRINT*,'probleme au point ij=',ij,'  l=',l
-              !    PRINT*,'sigs=',zsigs(ij),'  sign=',zsign(ij)
-              !    stop
-              ! END IF
-           else
-              zsign(ij)=0.5
-              zsigs(ij)=0.5
-           endif
-        enddo
-
-  !   calcul de la pente maximum dans la maille en valeur absolue
-
-   do ij=1,ip1jm
-      IF (v_m(ij,l)>=0.) THEN
-         zsigp=zsign(ij+iip1)
-         zsigm=zsigs(ij+iip1)
-         zqp=qn(ij+iip1,l)
-         zqm=qs(ij+iip1,l)
-         zm=masse(ij+iip1,l)
-         zq=q(ij+iip1,l)
+  do l = 1, llm
+    do ij = 1, ip1jmp1
+      zdq = qn(ij, l) - qs(ij, l)
+      ! if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) THEN
+      !    PRINT*,'probleme au point ij=',ij,'  l=',l,'  advnqx'
+      !    PRINT*,qn(ij,l),q(ij,l),qs(ij,l)
+      !    qn(ij,l)=q(ij,l)
+      !    qs(ij,l)=q(ij,l)
+      ! END IF
+      IF(abs(zdq)>prec) THEN
+        zsign(ij) = (q(ij, l) - qs(ij, l)) / zdq
+        zsigs(ij) = 1. - zsign(ij)
+        ! IF(.NOT.(zsign(ij).ge.0..and.zsign(ij).le.1. .AND.
+        !    s               zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) THEN
+        !    PRINT*,'probleme au point ij=',ij,'  l=',l
+        !    PRINT*,'sigs=',zsigs(ij),'  sign=',zsign(ij)
+        !    stop
+        ! END IF
       else
-         zsigm=zsign(ij)
-         zsigp=zsigs(ij)
-         zqm=qn(ij,l)
-         zqp=qs(ij,l)
-         zm=masse(ij,l)
-         zq=q(ij,l)
-      endif
-      zsig=abs(v_m(ij,l))/zm
-      IF(zsig==0.) zsigp=0.1
+        zsign(ij) = 0.5
+        zsigs(ij) = 0.5
+      endif
+    enddo
+
+    !   calcul de la pente maximum dans la maille en valeur absolue
+
+    do ij = 1, ip1jm
+      IF (v_m(ij, l)>=0.) THEN
+        zsigp = zsign(ij + iip1)
+        zsigm = zsigs(ij + iip1)
+        zqp = qn(ij + iip1, l)
+        zqm = qs(ij + iip1, l)
+        zm = masse(ij + iip1, l)
+        zq = q(ij + iip1, l)
+      else
+        zsigm = zsign(ij)
+        zsigp = zsigs(ij)
+        zqm = qn(ij, l)
+        zqp = qs(ij, l)
+        zm = masse(ij, l)
+        zq = q(ij, l)
+      endif
+      zsig = abs(v_m(ij, l)) / zm
+      IF(zsig==0.) zsigp = 0.1
       IF (zsig<=zsigp) THEN
-          v_mq(ij,l)=v_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
+        v_mq(ij, l) = v_m(ij, l) * (zqp - 0.5 * zsig / zsigp * (zqp - zq))
       else
-          zz=0.5*(zsig-zsigp)/zsigm
-          v_mq(ij,l)=sign(zm,v_m(ij,l))*( 0.5*(zq+zqp)*zsigp &
-                +(zsig-zsigp)*(zq+zz*(zqm-zq)) )
-      endif
-   enddo
-  enddo
-
-  do l=1,llm
-     do ij=iip2,ip1jm
-        new_m=masse(ij,l) &
-              +v_m(ij,l)-v_m(ij-iip1,l)
-        q(ij,l)=(q(ij,l)*masse(ij,l)+v_mq(ij,l)-v_mq(ij-iip1,l)) &
-              /new_m
-        masse(ij,l)=new_m
-     enddo
-  !.-. ancienne version
-     convpn=SSUM(iim,v_mq(1,l),1)
-     convmpn=ssum(iim,v_m(1,l),1)
-     massen=ssum(iim,masse(1,l),1)
-     new_m=massen+convmpn
-     q(1,l)=(q(1,l)*massen+convpn)/new_m
-     do ij = 1,iip1
-        q(ij,l)=q(1,l)
-        masse(ij,l)=new_m*aire(ij)/apoln
-     enddo
-
-     convps=-SSUM(iim,v_mq(ip1jm-iim,l),1)
-     convmps=-ssum(iim,v_m(ip1jm-iim,l),1)
-     masses=ssum(iim,masse(ip1jm+1,l),1)
-     new_m=masses+convmps
-     q(ip1jm+1,l)=(q(ip1jm+1,l)*masses+convps)/new_m
-     do ij = ip1jm+1,ip1jmp1
-        q(ij,l)=q(ip1jm+1,l)
-        masse(ij,l)=new_m*aire(ij)/apols
-     enddo
+        zz = 0.5 * (zsig - zsigp) / zsigm
+        v_mq(ij, l) = sign(zm, v_m(ij, l)) * (0.5 * (zq + zqp) * zsigp &
+                + (zsig - zsigp) * (zq + zz * (zqm - zq)))
+      endif
+    enddo
+  enddo
+
+  do l = 1, llm
+    do ij = iip2, ip1jm
+      new_m = masse(ij, l) &
+              + v_m(ij, l) - v_m(ij - iip1, l)
+      q(ij, l) = (q(ij, l) * masse(ij, l) + v_mq(ij, l) - v_mq(ij - iip1, l)) &
+              / new_m
+      masse(ij, l) = new_m
+    enddo
+    !.-. ancienne version
+    convpn = SSUM(iim, v_mq(1, l), 1)
+    convmpn = ssum(iim, v_m(1, l), 1)
+    massen = ssum(iim, masse(1, l), 1)
+    new_m = massen + convmpn
+    q(1, l) = (q(1, l) * massen + convpn) / new_m
+    do ij = 1, iip1
+      q(ij, l) = q(1, l)
+      masse(ij, l) = new_m * aire(ij) / apoln
+    enddo
+
+    convps = -SSUM(iim, v_mq(ip1jm - iim, l), 1)
+    convmps = -ssum(iim, v_m(ip1jm - iim, l), 1)
+    masses = ssum(iim, masse(ip1jm + 1, l), 1)
+    new_m = masses + convmps
+    q(ip1jm + 1, l) = (q(ip1jm + 1, l) * masses + convps) / new_m
+    do ij = ip1jm + 1, ip1jmp1
+      q(ij, l) = q(ip1jm + 1, l)
+      masse(ij, l) = new_m * aire(ij) / apols
+    enddo
   enddo
 
   RETURN
 END SUBROUTINE advny
-SUBROUTINE advnz(q,qh,qb,masse,w_m)
+SUBROUTINE advnz(q, qh, qb, masse, w_m)
   !
   ! Auteurs:   F.Hourdin
@@ -853,4 +852,5 @@
   !
   !   --------------------------------------------------------------------
+  USE lmdz_iniprint, ONLY: lunout, prt_level
   IMPLICIT NONE
   !
@@ -858,12 +858,11 @@
   INCLUDE "paramet.h"
   INCLUDE "comgeom.h"
-  INCLUDE "iniprint.h"
   !
   !
   !   Arguments:
   !   ----------
-  REAL :: masse(ip1jmp1,llm)
-  REAL :: w_m( ip1jmp1,llm+1)
-  REAL :: q(ip1jmp1,llm),qb(ip1jmp1,llm),qh(ip1jmp1,llm)
+  REAL :: masse(ip1jmp1, llm)
+  REAL :: w_m(ip1jmp1, llm + 1)
+  REAL :: q(ip1jmp1, llm), qb(ip1jmp1, llm), qh(ip1jmp1, llm)
 
   !
@@ -871,10 +870,10 @@
   !   ---------
   !
-  INTEGER :: ij,l
-  !
-  REAL :: new_m,zdq,zz
-  REAL :: zsigh(ip1jmp1,llm),zsigb(ip1jmp1,llm),zsig
-  REAL :: w_mq(ip1jmp1,llm+1)
-  REAL :: zm,zq,zsigm,zsigp,zqm,zqp
+  INTEGER :: ij, l
+  !
+  REAL :: new_m, zdq, zz
+  REAL :: zsigh(ip1jmp1, llm), zsigb(ip1jmp1, llm), zsig
+  REAL :: w_mq(ip1jmp1, llm + 1)
+  REAL :: zm, zq, zsigm, zsigp, zqm, zqp
   REAL :: prec
   save prec
@@ -882,68 +881,68 @@
   data prec/1.e-13/
 
-  do l=1,llm
-        do ij=1,ip1jmp1
-           zdq=qb(ij,l)-qh(ij,l)
-           ! if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) THEN
-           !    PRINT*,'probleme au point ij=',ij,'  l=',l
-           !    PRINT*,qh(ij,l),q(ij,l),qb(ij,l)
-           !    qh(ij,l)=q(ij,l)
-           !    qb(ij,l)=q(ij,l)
-           ! END IF
-
-           IF(abs(zdq)>prec) THEN
-              zsigb(ij,l)=(q(ij,l)-qh(ij,l))/zdq
-              zsigh(ij,l)=1.-zsigb(ij,l)
-              zsigb(ij,l)=min(max(zsigb(ij,l),0.),1.)
-           else
-              zsigb(ij,l)=0.5
-              zsigh(ij,l)=0.5
-           endif
-        enddo
-   enddo
-
-   ! PRINT*,'ok1'
+  do l = 1, llm
+    do ij = 1, ip1jmp1
+      zdq = qb(ij, l) - qh(ij, l)
+      ! if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) THEN
+      !    PRINT*,'probleme au point ij=',ij,'  l=',l
+      !    PRINT*,qh(ij,l),q(ij,l),qb(ij,l)
+      !    qh(ij,l)=q(ij,l)
+      !    qb(ij,l)=q(ij,l)
+      ! END IF
+
+      IF(abs(zdq)>prec) THEN
+        zsigb(ij, l) = (q(ij, l) - qh(ij, l)) / zdq
+        zsigh(ij, l) = 1. - zsigb(ij, l)
+        zsigb(ij, l) = min(max(zsigb(ij, l), 0.), 1.)
+      else
+        zsigb(ij, l) = 0.5
+        zsigh(ij, l) = 0.5
+      endif
+    enddo
+  enddo
+
+  ! PRINT*,'ok1'
   !   calcul de la pente maximum dans la maille en valeur absolue
-   do l=2,llm
-   do ij=1,ip1jmp1
-      IF (w_m(ij,l)>=0.) THEN
-         zsigp=zsigb(ij,l)
-         zsigm=zsigh(ij,l)
-         zqp=qb(ij,l)
-         zqm=qh(ij,l)
-         zm=masse(ij,l)
-         zq=q(ij,l)
+  do l = 2, llm
+    do ij = 1, ip1jmp1
+      IF (w_m(ij, l)>=0.) THEN
+        zsigp = zsigb(ij, l)
+        zsigm = zsigh(ij, l)
+        zqp = qb(ij, l)
+        zqm = qh(ij, l)
+        zm = masse(ij, l)
+        zq = q(ij, l)
       else
-         zsigm=zsigb(ij,l-1)
-         zsigp=zsigh(ij,l-1)
-         zqm=qb(ij,l-1)
-         zqp=qh(ij,l-1)
-         zm=masse(ij,l-1)
-         zq=q(ij,l-1)
-      endif
-      zsig=abs(w_m(ij,l))/zm
-      IF(zsig==0.) zsigp=0.1
+        zsigm = zsigb(ij, l - 1)
+        zsigp = zsigh(ij, l - 1)
+        zqm = qb(ij, l - 1)
+        zqp = qh(ij, l - 1)
+        zm = masse(ij, l - 1)
+        zq = q(ij, l - 1)
+      endif
+      zsig = abs(w_m(ij, l)) / zm
+      IF(zsig==0.) zsigp = 0.1
       IF (zsig<=zsigp) THEN
-          w_mq(ij,l)=w_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
+        w_mq(ij, l) = w_m(ij, l) * (zqp - 0.5 * zsig / zsigp * (zqp - zq))
       else
-          zz=0.5*(zsig-zsigp)/zsigm
-          w_mq(ij,l)=sign(zm,w_m(ij,l))*( 0.5*(zq+zqp)*zsigp &
-                +(zsig-zsigp)*(zq+zz*(zqm-zq)) )
-      endif
-  enddo
-  enddo
-
-   do ij=1,ip1jmp1
-      w_mq(ij,llm+1)=0.
-      w_mq(ij,1)=0.
-   enddo
-
-  do l=1,llm
-     do ij=1,ip1jmp1
-        new_m=masse(ij,l)+w_m(ij,l+1)-w_m(ij,l)
-        q(ij,l)=(q(ij,l)*masse(ij,l)+w_mq(ij,l+1)-w_mq(ij,l)) &
-              /new_m
-        masse(ij,l)=new_m
-     enddo
+        zz = 0.5 * (zsig - zsigp) / zsigm
+        w_mq(ij, l) = sign(zm, w_m(ij, l)) * (0.5 * (zq + zqp) * zsigp &
+                + (zsig - zsigp) * (zq + zz * (zqm - zq)))
+      endif
+    enddo
+  enddo
+
+  do ij = 1, ip1jmp1
+    w_mq(ij, llm + 1) = 0.
+    w_mq(ij, 1) = 0.
+  enddo
+
+  do l = 1, llm
+    do ij = 1, ip1jmp1
+      new_m = masse(ij, l) + w_m(ij, l + 1) - w_m(ij, l)
+      q(ij, l) = (q(ij, l) * masse(ij, l) + w_mq(ij, l + 1) - w_mq(ij, l)) &
+              / new_m
+      masse(ij, l) = new_m
+    enddo
   enddo
   ! PRINT*,'ok3'
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/diagedyn.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/diagedyn.f90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/diagedyn.f90	(revision 5118)
@@ -53,4 +53,5 @@
 
   USE control_mod, ONLY: planet_type
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -59,5 +60,4 @@
   INCLUDE "paramet.h"
   INCLUDE "comgeom.h"
-  INCLUDE "iniprint.h"
 
   ! Ehouarn: for now set these parameters to what is in Earth physics...
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert.F90	(revision 5118)
@@ -9,4 +9,5 @@
                          pseudoalt, pa, preff, scaleheight, presinter
   USE logic_mod, ONLY: ok_strato
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -14,5 +15,4 @@
   include "dimensions.h"
   include "paramet.h"
-  include "iniprint.h"
 
 !-------------------------------------------------------------------------------
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert_noterre.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert_noterre.f90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert_noterre.f90	(revision 5118)
@@ -11,4 +11,5 @@
   USE comconst_mod, ONLY: kappa
   USE logic_mod, ONLY: hybrid
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -16,5 +17,4 @@
   include "dimensions.h"
   include "paramet.h"
-  include "iniprint.h"
   !
   !=======================================================================
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/infotrac.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/infotrac.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/infotrac.F90	(revision 5118)
@@ -118,4 +118,5 @@
 #endif
     USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_INCA, CPPKEY_STRATAER
+    USE lmdz_iniprint, ONLY: lunout, prt_level
     IMPLICIT NONE
     !==============================================================================================================================
@@ -139,5 +140,4 @@
     ! Declarations:
     INCLUDE "dimensions.h"
-    INCLUDE "iniprint.h"
 
     !------------------------------------------------------------------------------------------------------------------------------
@@ -228,34 +228,34 @@
     IF(fType == 1 .AND. ANY(['inca', 'inco']==type_trac)) THEN         !=== FOUND OLD STYLE INCA "traceur.def"
       !---------------------------------------------------------------------------------------------------------------------------
-        nqo = SIZE(tracers) - nqCO2
-        CALL Init_chem_inca_trac(nqINCA)                               !--- Get nqINCA from INCA
-        nbtr = nqINCA + nqCO2                                          !--- Number of tracers passed to phytrac
-        nqtrue = nbtr + nqo                                            !--- Total number of "true" tracers
-        IF(ALL([2, 3] /= nqo)) CALL abort_gcm(modname, 'Only 2 or 3 water phases allowed ; found nqo=' // TRIM(int2str(nqo)), 1)
-        ALLOCATE(hadv(nqtrue), hadv_inca(nqINCA), conv_flg_inca(nqINCA), solsym_inca(nqINCA))
-        ALLOCATE(vadv(nqtrue), vadv_inca(nqINCA), pbl_flg_inca(nqINCA))
-        CALL init_transport(solsym_inca, conv_flg_inca, pbl_flg_inca, hadv_inca, vadv_inca)
-        ALLOCATE(ttr(nqtrue))
-        ttr(1:nqo + nqCO2) = tracers
-        ttr(1:nqo)%component = 'lmdz'
-        ttr(1 + nqo:nqCO2 + nqo)%component = 'co2i'
-        ttr(1 + nqo + nqCO2:nqtrue)%component = 'inca'
-        ttr(1 + nqo:nqtrue)%name = [('CO2     ', iq = 1, nqCO2), solsym_inca]
-        ttr(1 + nqo + nqCO2:nqtrue)%parent = tran0
-        ttr(1 + nqo + nqCO2:nqtrue)%phase = 'g'
-        lerr = getKey('hadv', had, ky = tracers(:)%keys)
-        lerr = getKey('vadv', vad, ky = tracers(:)%keys)
-        hadv(1:nqo + nqCO2) = had(:); hadv(1 + nqo + nqCO2:nqtrue) = hadv_inca
-        vadv(1:nqo + nqCO2) = vad(:); vadv(1 + nqo + nqCO2:nqtrue) = vadv_inca
-        CALL MOVE_ALLOC(FROM = ttr, TO = tracers)
-        DO iq = 1, nqtrue
-          t1 => tracers(iq)
-          CALL addKey('name', t1%name, t1%keys)
-          CALL addKey('component', t1%component, t1%keys)
-          CALL addKey('parent', t1%parent, t1%keys)
-          CALL addKey('phase', t1%phase, t1%keys)
-        END DO
-        IF(setGeneration(tracers)) CALL abort_gcm(modname, 'See above', 1) !- SET FIELDS %iGeneration, %gen0Name
-        DEALLOCATE(had, hadv_inca, vad, vadv_inca, conv_flg_inca, pbl_flg_inca, solsym_inca)
+      nqo = SIZE(tracers) - nqCO2
+      CALL Init_chem_inca_trac(nqINCA)                               !--- Get nqINCA from INCA
+      nbtr = nqINCA + nqCO2                                          !--- Number of tracers passed to phytrac
+      nqtrue = nbtr + nqo                                            !--- Total number of "true" tracers
+      IF(ALL([2, 3] /= nqo)) CALL abort_gcm(modname, 'Only 2 or 3 water phases allowed ; found nqo=' // TRIM(int2str(nqo)), 1)
+      ALLOCATE(hadv(nqtrue), hadv_inca(nqINCA), conv_flg_inca(nqINCA), solsym_inca(nqINCA))
+      ALLOCATE(vadv(nqtrue), vadv_inca(nqINCA), pbl_flg_inca(nqINCA))
+      CALL init_transport(solsym_inca, conv_flg_inca, pbl_flg_inca, hadv_inca, vadv_inca)
+      ALLOCATE(ttr(nqtrue))
+      ttr(1:nqo + nqCO2) = tracers
+      ttr(1:nqo)%component = 'lmdz'
+      ttr(1 + nqo:nqCO2 + nqo)%component = 'co2i'
+      ttr(1 + nqo + nqCO2:nqtrue)%component = 'inca'
+      ttr(1 + nqo:nqtrue)%name = [('CO2     ', iq = 1, nqCO2), solsym_inca]
+      ttr(1 + nqo + nqCO2:nqtrue)%parent = tran0
+      ttr(1 + nqo + nqCO2:nqtrue)%phase = 'g'
+      lerr = getKey('hadv', had, ky = tracers(:)%keys)
+      lerr = getKey('vadv', vad, ky = tracers(:)%keys)
+      hadv(1:nqo + nqCO2) = had(:); hadv(1 + nqo + nqCO2:nqtrue) = hadv_inca
+      vadv(1:nqo + nqCO2) = vad(:); vadv(1 + nqo + nqCO2:nqtrue) = vadv_inca
+      CALL MOVE_ALLOC(FROM = ttr, TO = tracers)
+      DO iq = 1, nqtrue
+        t1 => tracers(iq)
+        CALL addKey('name', t1%name, t1%keys)
+        CALL addKey('component', t1%component, t1%keys)
+        CALL addKey('parent', t1%parent, t1%keys)
+        CALL addKey('phase', t1%phase, t1%keys)
+      END DO
+      IF(setGeneration(tracers)) CALL abort_gcm(modname, 'See above', 1) !- SET FIELDS %iGeneration, %gen0Name
+      DEALLOCATE(had, hadv_inca, vad, vadv_inca, conv_flg_inca, pbl_flg_inca, solsym_inca)
       !---------------------------------------------------------------------------------------------------------------------------
     ELSE                                                              !=== OTHER CASES (OLD OR NEW FORMAT, NO INCA MODULE)
@@ -266,5 +266,5 @@
       nbtr = nqtrue - COUNT(delPhase(tracers(:)%gen0Name) == 'H2O' &
               .AND. tracers(:)%component == 'lmdz') !--- Number of tracers passed to phytrac
-        nqINCA = COUNT(tracers(:)%component == 'inca')
+      nqINCA = COUNT(tracers(:)%component == 'inca')
       lerr = getKey('hadv', hadv, ky = tracers(:)%keys)
       lerr = getKey('vadv', vadv, ky = tracers(:)%keys)
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iniconst.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iniconst.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iniconst.F90	(revision 5118)
@@ -1,3 +1,2 @@
-
 ! $Id$
 
@@ -7,7 +6,8 @@
   USE IOIPSL
   USE comconst_mod, ONLY: im, imp1, jm, jmp1, lllm, lllmm1, lllmp1, &
-                          unsim, pi, r, kappa, cpp, dtvr, dtphys
+          unsim, pi, r, kappa, cpp, dtvr, dtphys
   USE comvert_mod, ONLY: disvert_type, pressure_exner
-  
+  USE lmdz_iniprint, ONLY: lunout, prt_level
+
   IMPLICIT NONE
 
@@ -19,8 +19,7 @@
   include "dimensions.h"
   include "paramet.h"
-  include "iniprint.h"
 
-  CHARACTER(LEN=*),parameter :: modname="iniconst"
-  CHARACTER(LEN=80) :: abort_message
+  CHARACTER(LEN = *), parameter :: modname = "iniconst"
+  CHARACTER(LEN = 80) :: abort_message
 
 
@@ -30,23 +29,23 @@
   !   ----------------------
 
-  im      = iim
-  jm      = jjm
-  lllm    = llm
-  imp1    = iim 
-  jmp1    = jjm + 1
-  lllmm1  = llm - 1
-  lllmp1  = llm + 1
+  im = iim
+  jm = jjm
+  lllm = llm
+  imp1 = iim
+  jmp1 = jjm + 1
+  lllmm1 = llm - 1
+  lllmp1 = llm + 1
 
   !-----------------------------------------------------------------------
 
-  dtphys  = iphysiq * dtvr
-  unsim   = 1./iim
-  pi      = 2.*ASIN( 1. )
+  dtphys = iphysiq * dtvr
+  unsim = 1. / iim
+  pi = 2. * ASIN(1.)
 
   !-----------------------------------------------------------------------
 
-  r       = cpp * kappa
+  r = cpp * kappa
 
-  WRITE(lunout,*) trim(modname),': R  CP  Kappa ',r,cpp,kappa
+  WRITE(lunout, *) trim(modname), ': R  CP  Kappa ', r, cpp, kappa
 
   !-----------------------------------------------------------------------
@@ -54,11 +53,11 @@
   ! vertical discretization: default behavior depends on planet_type flag
   IF (planet_type=="earth") THEN
-     disvert_type=1
+    disvert_type = 1
   else
-     disvert_type=2
+    disvert_type = 2
   ENDIF
   ! but user can also specify using one or the other in run.def:
-  CALL getin('disvert_type',disvert_type)
-  WRITE(lunout,*) trim(modname),': disvert_type=',disvert_type
+  CALL getin('disvert_type', disvert_type)
+  WRITE(lunout, *) trim(modname), ': disvert_type=', disvert_type
 
   pressure_exner = disvert_type == 1 ! default value
@@ -66,12 +65,12 @@
 
   IF (disvert_type==1) THEN
-     ! standard case for Earth (automatic generation of levels)
-     CALL disvert()
+    ! standard case for Earth (automatic generation of levels)
+    CALL disvert()
   ELSE IF (disvert_type==2) THEN
-     ! standard case for planets (levels generated using z2sig.def file)
-     CALL disvert_noterre
+    ! standard case for planets (levels generated using z2sig.def file)
+    CALL disvert_noterre
   else
-     WRITE(abort_message,*) "Wrong value for disvert_type: ", disvert_type
-     CALL abort_gcm(modname,abort_message,0)
+    WRITE(abort_message, *) "Wrong value for disvert_type: ", disvert_type
+    CALL abort_gcm(modname, abort_message, 0)
   ENDIF
 
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inidissip.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inidissip.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inidissip.F90	(revision 5118)
@@ -17,4 +17,5 @@
   USE lmdz_libmath, ONLY: minmax
   USE lmdz_ran1, ONLY: ran1
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -22,5 +23,4 @@
   include "paramet.h"
   include "comdissipn.h"
-  include "iniprint.h"
 
   LOGICAL, INTENT(IN) :: lstardis
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/initdynav.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/initdynav.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/initdynav.F90	(revision 5118)
@@ -1,15 +1,16 @@
 ! $Id$
 
-SUBROUTINE initdynav(day0,anne0,tstep,t_ops,t_wrt)
+SUBROUTINE initdynav(day0, anne0, tstep, t_ops, t_wrt)
 
   USE IOIPSL
   USE infotrac, ONLY: nqtot
-  USE com_io_dyn_mod, ONLY: histaveid,histvaveid,histuaveid, &
-       dynhistave_file,dynhistvave_file,dynhistuave_file
+  USE com_io_dyn_mod, ONLY: histaveid, histvaveid, histuaveid, &
+          dynhistave_file, dynhistvave_file, dynhistuave_file
   USE comconst_mod, ONLY: pi
   USE comvert_mod, ONLY: presnivs
   USE temps_mod, ONLY: itau_dyn
   USE lmdz_description, ONLY: descript
-  
+  USE lmdz_iniprint, ONLY: lunout, prt_level
+
   IMPLICIT NONE
 
@@ -38,5 +39,4 @@
   include "paramet.h"
   include "comgeom.h"
-  include "iniprint.h"
 
   !   Arguments
@@ -51,7 +51,7 @@
   REAL zjulian
   INTEGER iq
-  REAL rlong(iip1,jjp1), rlat(iip1,jjp1)
+  REAL rlong(iip1, jjp1), rlat(iip1, jjp1)
   INTEGER uhoriid, vhoriid, thoriid, zvertiid
-  INTEGER ii,jj
+  INTEGER ii, jj
   INTEGER zan, dayref
 
@@ -64,5 +64,4 @@
   !  Appel a histbeg: creation du fichier netcdf et initialisations diverses
 
-
   zan = anne0
   dayref = day0
@@ -71,8 +70,8 @@
 
   do jj = 1, jjp1
-     do ii = 1, iip1
-        rlong(ii,jj) = rlonv(ii) * 180. / pi
-        rlat(ii,jj)  = rlatu(jj) * 180. / pi
-     enddo
+    do ii = 1, iip1
+      rlong(ii, jj) = rlonv(ii) * 180. / pi
+      rlat(ii, jj) = rlatu(jj) * 180. / pi
+    enddo
   enddo
 
@@ -80,7 +79,7 @@
   ! Restriction de IOIPSL: seulement 2 coordonnees dans le meme fichier
   ! Grille Scalaire
-  CALL histbeg(dynhistave_file, iip1, rlong(:,1), jjp1, rlat(1,:), &
-       1, iip1, 1, jjp1, &
-       tau0, zjulian, tstep, thoriid,histaveid)
+  CALL histbeg(dynhistave_file, iip1, rlong(:, 1), jjp1, rlat(1, :), &
+          1, iip1, 1, jjp1, &
+          tau0, zjulian, tstep, thoriid, histaveid)
 
   ! Creation du fichier histoire pour les grilles en V et U (oblige
@@ -89,33 +88,33 @@
   ! Grille V
   do jj = 1, jjm
-     do ii = 1, iip1
-        rlong(ii,jj) = rlonv(ii) * 180. / pi
-        rlat(ii,jj) = rlatv(jj) * 180. / pi
-     enddo
+    do ii = 1, iip1
+      rlong(ii, jj) = rlonv(ii) * 180. / pi
+      rlat(ii, jj) = rlatv(jj) * 180. / pi
+    enddo
   enddo
 
-  CALL histbeg(dynhistvave_file, iip1, rlong(:,1), jjm, rlat(1,:), &
-       1, iip1, 1, jjm, &
-       tau0, zjulian, tstep, vhoriid,histvaveid)
+  CALL histbeg(dynhistvave_file, iip1, rlong(:, 1), jjm, rlat(1, :), &
+          1, iip1, 1, jjm, &
+          tau0, zjulian, tstep, vhoriid, histvaveid)
   ! Grille U
   do jj = 1, jjp1
-     do ii = 1, iip1
-        rlong(ii,jj) = rlonu(ii) * 180. / pi
-        rlat(ii,jj) = rlatu(jj) * 180. / pi
-     enddo
+    do ii = 1, iip1
+      rlong(ii, jj) = rlonu(ii) * 180. / pi
+      rlat(ii, jj) = rlatu(jj) * 180. / pi
+    enddo
   enddo
 
-  CALL histbeg(dynhistuave_file, iip1, rlong(:,1),jjp1, rlat(1,:), &
-       1, iip1, 1, jjp1, &
-       tau0, zjulian, tstep, uhoriid,histuaveid)
+  CALL histbeg(dynhistuave_file, iip1, rlong(:, 1), jjp1, rlat(1, :), &
+          1, iip1, 1, jjp1, &
+          tau0, zjulian, tstep, uhoriid, histuaveid)
 
   !  Appel a histvert pour la grille verticale
 
-  CALL histvert(histaveid,'presnivs','Niveaux Pression approximatifs','mb', &
-       llm, presnivs/100., zvertiid,'down')
-  CALL histvert(histuaveid,'presnivs','Niveaux Pression approximatifs','mb', &
-       llm, presnivs/100., zvertiid,'down')
-  CALL histvert(histvaveid,'presnivs','Niveaux Pression approximatifs','mb', &
-       llm, presnivs/100., zvertiid,'down')
+  CALL histvert(histaveid, 'presnivs', 'Niveaux Pression approximatifs', 'mb', &
+          llm, presnivs / 100., zvertiid, 'down')
+  CALL histvert(histuaveid, 'presnivs', 'Niveaux Pression approximatifs', 'mb', &
+          llm, presnivs / 100., zvertiid, 'down')
+  CALL histvert(histvaveid, 'presnivs', 'Niveaux Pression approximatifs', 'mb', &
+          llm, presnivs / 100., zvertiid, 'down')
 
   !  Appels a histdef pour la definition des variables a sauvegarder
@@ -124,12 +123,12 @@
 
   CALL histdef(histuaveid, 'u', 'vent u moyen ', &
-       'm/s', iip1, jjp1, uhoriid, llm, 1, llm, zvertiid, &
-       32, 'ave(X)', t_ops, t_wrt)
+          'm/s', iip1, jjp1, uhoriid, llm, 1, llm, zvertiid, &
+          32, 'ave(X)', t_ops, t_wrt)
 
   !  Vents V
 
   CALL histdef(histvaveid, 'v', 'vent v moyen', &
-       'm/s', iip1, jjm, vhoriid, llm, 1, llm, zvertiid, &
-       32, 'ave(X)', t_ops, t_wrt)
+          'm/s', iip1, jjm, vhoriid, llm, 1, llm, zvertiid, &
+          32, 'ave(X)', t_ops, t_wrt)
 
 
@@ -137,18 +136,18 @@
 
   CALL histdef(histaveid, 'temp', 'temperature moyenne', 'K', &
-       iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
-       32, 'ave(X)', t_ops, t_wrt)
+          iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
+          32, 'ave(X)', t_ops, t_wrt)
 
   !  Temperature potentielle
 
   CALL histdef(histaveid, 'theta', 'temperature potentielle', 'K', &
-       iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
-       32, 'ave(X)', t_ops, t_wrt)
+          iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
+          32, 'ave(X)', t_ops, t_wrt)
 
   !  Geopotentiel
 
   CALL histdef(histaveid, 'phi', 'geopotentiel moyen', '-', &
-       iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
-       32, 'ave(X)', t_ops, t_wrt)
+          iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
+          32, 'ave(X)', t_ops, t_wrt)
 
   !  Traceurs
@@ -164,12 +163,12 @@
 
   CALL histdef(histaveid, 'masse', 'masse', 'kg', &
-       iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
-       32, 'ave(X)', t_ops, t_wrt)
+          iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
+          32, 'ave(X)', t_ops, t_wrt)
 
   !  Pression au sol
 
   CALL histdef(histaveid, 'ps', 'pression naturelle au sol', 'Pa', &
-       iip1, jjp1, thoriid, 1, 1, 1, -99, &
-       32, 'ave(X)', t_ops, t_wrt)
+          iip1, jjp1, thoriid, 1, 1, 1, -99, &
+          32, 'ave(X)', t_ops, t_wrt)
 
   !  Geopotentiel au sol
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/initfluxsto.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/initfluxsto.f90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/initfluxsto.f90	(revision 5118)
@@ -10,4 +10,5 @@
   USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn
   USE lmdz_description, ONLY: descript
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -43,5 +44,4 @@
   include "paramet.h"
   include "comgeom.h"
-  include "iniprint.h"
 
   !   Arguments
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inithist.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inithist.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inithist.F90	(revision 5118)
@@ -11,4 +11,5 @@
   USE temps_mod, ONLY: itau_dyn
   USE lmdz_description, ONLY: descript
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -42,5 +43,4 @@
   include "paramet.h"
   include "comgeom.h"
-  include "iniprint.h"
 
   !   Arguments
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/sortvarc.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/sortvarc.f90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/sortvarc.f90	(revision 5118)
@@ -1,16 +1,16 @@
-
 ! $Id$
 
 SUBROUTINE sortvarc &
-        (itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp,time , &
-        vcov )
+        (itau, ucov, teta, ps, masse, pk, phis, vorpot, phi, bern, dp, time, &
+        vcov)
 
   USE control_mod, ONLY: resetvarc
   USE comconst_mod, ONLY: dtvr, daysec, g, rad, omeg
   USE logic_mod, ONLY: read_start
-  USE ener_mod, ONLY: etot,ptot,ztot,stot,ang, &
-        etot0,ptot0,ztot0,stot0,ang0, &
-        rmsdpdt,rmsv
+  USE ener_mod, ONLY: etot, ptot, ztot, stot, ang, &
+          etot0, ptot0, ztot0, stot0, ang0, &
+          rmsdpdt, rmsv
   USE lmdz_filtreg, ONLY: filtreg
+  USE lmdz_iniprint, ONLY: lunout, prt_level
   IMPLICIT NONE
 
@@ -34,128 +34,127 @@
   INCLUDE "paramet.h"
   INCLUDE "comgeom.h"
-  INCLUDE "iniprint.h"
 
   !   Arguments:
   !   ----------
 
-  INTEGER,INTENT(IN) :: itau
-  REAL,INTENT(IN) :: ucov(ip1jmp1,llm)
-  REAL,INTENT(IN) :: teta(ip1jmp1,llm)
-  REAL,INTENT(IN) :: masse(ip1jmp1,llm)
-  REAL,INTENT(IN) :: vcov(ip1jm,llm)
-  REAL,INTENT(IN) :: ps(ip1jmp1)
-  REAL,INTENT(IN) :: phis(ip1jmp1)
-  REAL,INTENT(IN) :: vorpot(ip1jm,llm)
-  REAL,INTENT(IN) :: phi(ip1jmp1,llm)
-  REAL,INTENT(IN) :: bern(ip1jmp1,llm)
-  REAL,INTENT(IN) :: dp(ip1jmp1)
-  REAL,INTENT(IN) :: time
-  REAL,INTENT(IN) :: pk(ip1jmp1,llm)
+  INTEGER, INTENT(IN) :: itau
+  REAL, INTENT(IN) :: ucov(ip1jmp1, llm)
+  REAL, INTENT(IN) :: teta(ip1jmp1, llm)
+  REAL, INTENT(IN) :: masse(ip1jmp1, llm)
+  REAL, INTENT(IN) :: vcov(ip1jm, llm)
+  REAL, INTENT(IN) :: ps(ip1jmp1)
+  REAL, INTENT(IN) :: phis(ip1jmp1)
+  REAL, INTENT(IN) :: vorpot(ip1jm, llm)
+  REAL, INTENT(IN) :: phi(ip1jmp1, llm)
+  REAL, INTENT(IN) :: bern(ip1jmp1, llm)
+  REAL, INTENT(IN) :: dp(ip1jmp1)
+  REAL, INTENT(IN) :: time
+  REAL, INTENT(IN) :: pk(ip1jmp1, llm)
 
   !   Local:
   !   ------
 
-  REAL :: vor(ip1jm),bernf(ip1jmp1,llm),ztotl(llm)
-  REAL :: etotl(llm),stotl(llm),rmsvl(llm),angl(llm),ge(ip1jmp1)
-  REAL :: cosphi(ip1jm),omegcosp(ip1jm)
-  REAL :: dtvrs1j,rjour,heure,radsg,radomeg
-  REAL :: massebxy(ip1jm,llm)
+  REAL :: vor(ip1jm), bernf(ip1jmp1, llm), ztotl(llm)
+  REAL :: etotl(llm), stotl(llm), rmsvl(llm), angl(llm), ge(ip1jmp1)
+  REAL :: cosphi(ip1jm), omegcosp(ip1jm)
+  REAL :: dtvrs1j, rjour, heure, radsg, radomeg
+  REAL :: massebxy(ip1jm, llm)
   INTEGER :: l, ij, imjmp1
 
   REAL :: SSUM
-  LOGICAL,SAVE :: firstcal=.TRUE.
-  CHARACTER(LEN=*),PARAMETER :: modname="sortvarc"
+  LOGICAL, SAVE :: firstcal = .TRUE.
+  CHARACTER(LEN = *), PARAMETER :: modname = "sortvarc"
 
   !-----------------------------------------------------------------------
   ! Ehouarn: when no initialization fields from file, resetvarc should be
-       ! set to false
-   IF (firstcal) THEN
-     IF (.NOT.read_start) THEN
-       resetvarc=.TRUE.
-     endif
-   endif
-
-   dtvrs1j   = dtvr/daysec
-   rjour     = REAL( INT( itau * dtvrs1j ))
-   heure     = ( itau*dtvrs1j-rjour ) * 24.
-   imjmp1    = iim * jjp1
-   IF(ABS(heure - 24.)<=0.0001 ) heure = 0.
-  !
-   CALL massbarxy ( masse, massebxy )
+  ! set to false
+  IF (firstcal) THEN
+    IF (.NOT.read_start) THEN
+      resetvarc = .TRUE.
+    endif
+  endif
+
+  dtvrs1j = dtvr / daysec
+  rjour = REAL(INT(itau * dtvrs1j))
+  heure = (itau * dtvrs1j - rjour) * 24.
+  imjmp1 = iim * jjp1
+  IF(ABS(heure - 24.)<=0.0001) heure = 0.
+  !
+  CALL massbarxy (masse, massebxy)
 
   !   .....  Calcul  de  rmsdpdt  .....
 
-   ge(:)=dp(:)*dp(:)
-
-   rmsdpdt = SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
-  !
-   rmsdpdt = daysec* 1.e-2 * SQRT(rmsdpdt/imjmp1)
-
-   CALL SCOPY( ijp1llm,bern,1,bernf,1 )
-   CALL filtreg(bernf,jjp1,llm,-2,2,.TRUE.,1)
+  ge(:) = dp(:) * dp(:)
+
+  rmsdpdt = SSUM(ip1jmp1, ge, 1) - SSUM(jjp1, ge, iip1)
+  !
+  rmsdpdt = daysec * 1.e-2 * SQRT(rmsdpdt / imjmp1)
+
+  CALL SCOPY(ijp1llm, bern, 1, bernf, 1)
+  CALL filtreg(bernf, jjp1, llm, -2, 2, .TRUE., 1)
 
   !   .....  Calcul du moment  angulaire   .....
 
-   radsg    = rad /g
-   radomeg  = rad * omeg
-  !
-   DO ij=iip2,ip1jm
-      cosphi( ij ) = COS(rlatu((ij-1)/iip1+1))
-      omegcosp(ij) = radomeg   * cosphi(ij)
-   ENDDO
+  radsg = rad / g
+  radomeg = rad * omeg
+  !
+  DO ij = iip2, ip1jm
+    cosphi(ij) = COS(rlatu((ij - 1) / iip1 + 1))
+    omegcosp(ij) = radomeg * cosphi(ij)
+  ENDDO
 
   !  ...  Calcul  de l'energie,de l'enstrophie,de l'entropie et de rmsv  .
 
-   DO l=1,llm
-      DO ij = 1,ip1jm
-         vor(ij)=vorpot(ij,l)*vorpot(ij,l)*massebxy(ij,l)
-      ENDDO
-      ztotl(l)=(SSUM(ip1jm,vor,1)-SSUM(jjm,vor,iip1))
-
-      DO ij = 1,ip1jmp1
-         ge(ij)= masse(ij,l)*(phis(ij)+teta(ij,l)*pk(ij,l)  + &
-               bernf(ij,l)-phi(ij,l))
-      ENDDO
-      etotl(l) = SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
-
-      DO   ij   = 1, ip1jmp1
-         ge(ij) = masse(ij,l)*teta(ij,l)
-      ENDDO
-      stotl(l)= SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
-
-      DO ij=1,ip1jmp1
-         ge(ij)=masse(ij,l)*AMAX1(bernf(ij,l)-phi(ij,l),0.)
-      ENDDO
-      rmsvl(l)=2.*(SSUM(ip1jmp1,ge,1)-SSUM(jjp1,ge,iip1))
-
-      DO ij =iip2,ip1jm
-         ge(ij)=(ucov(ij,l)/cu(ij)+omegcosp(ij))*masse(ij,l) * &
-               cosphi(ij)
-      ENDDO
-      angl(l) = rad * &
-            (SSUM(ip1jm-iip1,ge(iip2),1)-SSUM(jjm-1,ge(iip2),iip1))
+  DO l = 1, llm
+    DO ij = 1, ip1jm
+      vor(ij) = vorpot(ij, l) * vorpot(ij, l) * massebxy(ij, l)
+    ENDDO
+    ztotl(l) = (SSUM(ip1jm, vor, 1) - SSUM(jjm, vor, iip1))
+
+    DO ij = 1, ip1jmp1
+      ge(ij) = masse(ij, l) * (phis(ij) + teta(ij, l) * pk(ij, l) + &
+              bernf(ij, l) - phi(ij, l))
+    ENDDO
+    etotl(l) = SSUM(ip1jmp1, ge, 1) - SSUM(jjp1, ge, iip1)
+
+    DO   ij = 1, ip1jmp1
+      ge(ij) = masse(ij, l) * teta(ij, l)
+    ENDDO
+    stotl(l) = SSUM(ip1jmp1, ge, 1) - SSUM(jjp1, ge, iip1)
+
+    DO ij = 1, ip1jmp1
+      ge(ij) = masse(ij, l) * AMAX1(bernf(ij, l) - phi(ij, l), 0.)
+    ENDDO
+    rmsvl(l) = 2. * (SSUM(ip1jmp1, ge, 1) - SSUM(jjp1, ge, iip1))
+
+    DO ij = iip2, ip1jm
+      ge(ij) = (ucov(ij, l) / cu(ij) + omegcosp(ij)) * masse(ij, l) * &
+              cosphi(ij)
+    ENDDO
+    angl(l) = rad * &
+            (SSUM(ip1jm - iip1, ge(iip2), 1) - SSUM(jjm - 1, ge(iip2), iip1))
   ENDDO
 
-      DO ij=1,ip1jmp1
-        ge(ij)= ps(ij)*aire(ij)
-      ENDDO
-  ptot  = SSUM(ip1jmp1,ge,1)-SSUM(jjp1,ge,iip1)
-  etot  = SSUM(     llm, etotl, 1 )
-  ztot  = SSUM(     llm, ztotl, 1 )
-  stot  = SSUM(     llm, stotl, 1 )
-  rmsv  = SSUM(     llm, rmsvl, 1 )
-  ang   = SSUM(     llm,  angl, 1 )
+  DO ij = 1, ip1jmp1
+    ge(ij) = ps(ij) * aire(ij)
+  ENDDO
+  ptot = SSUM(ip1jmp1, ge, 1) - SSUM(jjp1, ge, iip1)
+  etot = SSUM(llm, etotl, 1)
+  ztot = SSUM(llm, ztotl, 1)
+  stot = SSUM(llm, stotl, 1)
+  rmsv = SSUM(llm, rmsvl, 1)
+  ang = SSUM(llm, angl, 1)
 
   IF (firstcal.AND.resetvarc) THEN
-     WRITE(lunout,3500) itau, rjour, heure, time
-     WRITE(lunout,*) trim(modname), &
-           ' WARNING!!! Recomputing initial values of : '
-     WRITE(lunout,*) 'ptot,rmsdpdt,etot,ztot,stot,rmsv,ang'
-     WRITE(lunout,*) ptot,rmsdpdt,etot,ztot,stot,rmsv,ang
-     etot0 = etot
-     ptot0 = ptot
-     ztot0 = ztot
-     stot0 = stot
-     ang0  = ang
+    WRITE(lunout, 3500) itau, rjour, heure, time
+    WRITE(lunout, *) trim(modname), &
+            ' WARNING!!! Recomputing initial values of : '
+    WRITE(lunout, *) 'ptot,rmsdpdt,etot,ztot,stot,rmsv,ang'
+    WRITE(lunout, *) ptot, rmsdpdt, etot, ztot, stot, rmsv, ang
+    etot0 = etot
+    ptot0 = ptot
+    ztot0 = ztot
+    stot0 = stot
+    ang0 = ang
   END IF
 
@@ -163,41 +162,41 @@
   ! are zero, which can happen when using iniacademic)
   IF (etot0/=0) THEN
-    etot= etot/etot0
-  else
-    etot=1.
-  ENDIF
-  rmsv= SQRT(rmsv/ptot)
+    etot = etot / etot0
+  else
+    etot = 1.
+  ENDIF
+  rmsv = SQRT(rmsv / ptot)
   IF (ptot0/=0) THEN
-    ptot= ptot/ptot0
-  else
-    ptot=1.
+    ptot = ptot / ptot0
+  else
+    ptot = 1.
   ENDIF
   IF (ztot0/=0) THEN
-    ztot= ztot/ztot0
-  else
-    ztot=1.
+    ztot = ztot / ztot0
+  else
+    ztot = 1.
   ENDIF
   IF (stot0/=0) THEN
-    stot= stot/stot0
-  else
-    stot=1.
+    stot = stot / stot0
+  else
+    stot = 1.
   ENDIF
   IF (ang0/=0) THEN
-    ang = ang /ang0
-  else
-    ang=1.
+    ang = ang / ang0
+  else
+    ang = 1.
   ENDIF
 
   firstcal = .FALSE.
 
-  WRITE(lunout,3500) itau, rjour, heure, time
-  WRITE(lunout,4000) ptot,rmsdpdt,etot,ztot,stot,rmsv,ang
-
-3500   FORMAT(10("*"),4x,'pas',i7,5x,'jour',f9.0,'heure',f5.1,4x &
-             ,'date',f14.4,4x,10("*"))
-4000   FORMAT(10x,'masse',4x,'rmsdpdt',7x,'energie',2x,'enstrophie' &
-             ,2x,'entropie',3x,'rmsv',4x,'mt.ang',/,'GLOB  ' &
-             ,f10.6,e13.6,5f10.3/ &
-             )
+  WRITE(lunout, 3500) itau, rjour, heure, time
+  WRITE(lunout, 4000) ptot, rmsdpdt, etot, ztot, stot, rmsv, ang
+
+  3500   FORMAT(10("*"), 4x, 'pas', i7, 5x, 'jour', f9.0, 'heure', f5.1, 4x &
+          , 'date', f14.4, 4x, 10("*"))
+  4000   FORMAT(10x, 'masse', 4x, 'rmsdpdt', 7x, 'energie', 2x, 'enstrophie' &
+          , 2x, 'entropie', 3x, 'rmsv', 4x, 'mt.ang', /, 'GLOB  ' &
+          , f10.6, e13.6, 5f10.3/ &
+          )
 END SUBROUTINE sortvarc
 
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writedynav.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writedynav.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writedynav.F90	(revision 5118)
@@ -9,4 +9,5 @@
   USE temps_mod, ONLY: itau_dyn
   USE lmdz_description, ONLY: descript
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -33,12 +34,11 @@
   include "paramet.h"
   include "comgeom.h"
-  include "iniprint.h"
 
   !   Arguments
 
-  REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) 
-  REAL teta(ip1jmp1*llm), phi(ip1jmp1, llm), ppk(ip1jmp1*llm)     
-  REAL ps(ip1jmp1), masse(ip1jmp1, llm)                   
-  REAL phis(ip1jmp1)                  
+  REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm)
+  REAL teta(ip1jmp1 * llm), phi(ip1jmp1, llm), ppk(ip1jmp1 * llm)
+  REAL ps(ip1jmp1), masse(ip1jmp1, llm)
+  REAL phis(ip1jmp1)
   REAL q(ip1jmp1, llm, nqtot)
   INTEGER time
@@ -47,7 +47,7 @@
   !   Variables locales
 
-  INTEGER ndex2d(ip1jmp1), ndexu(ip1jmp1*llm), ndexv(ip1jm*llm)
+  INTEGER ndex2d(ip1jmp1), ndexu(ip1jmp1 * llm), ndexv(ip1jm * llm)
   INTEGER iq, ii, ll
-  REAL tm(ip1jmp1*llm)
+  REAL tm(ip1jmp1 * llm)
   REAL vnat(ip1jm, llm), unat(ip1jmp1, llm)
   LOGICAL ok_sync
@@ -74,29 +74,29 @@
   !  Vents U
 
-  CALL histwrite(histuaveid, 'u', itau_w, unat,  &
-       iip1*jjp1*llm, ndexu)
+  CALL histwrite(histuaveid, 'u', itau_w, unat, &
+          iip1 * jjp1 * llm, ndexu)
 
   !  Vents V
 
-  CALL histwrite(histvaveid, 'v', itau_w, vnat,  &
-       iip1*jjm*llm, ndexv)
+  CALL histwrite(histvaveid, 'v', itau_w, vnat, &
+          iip1 * jjm * llm, ndexv)
 
   !  Temperature potentielle moyennee
 
-  CALL histwrite(histaveid, 'theta', itau_w, teta,  &
-       iip1*jjp1*llm, ndexu)
+  CALL histwrite(histaveid, 'theta', itau_w, teta, &
+          iip1 * jjp1 * llm, ndexu)
 
   !  Temperature moyennee
 
   do ii = 1, ijp1llm
-     tm(ii) = teta(ii) * ppk(ii)/cpp
+    tm(ii) = teta(ii) * ppk(ii) / cpp
   enddo
-  CALL histwrite(histaveid, 'temp', itau_w, tm,  &
-       iip1*jjp1*llm, ndexu)
+  CALL histwrite(histaveid, 'temp', itau_w, tm, &
+          iip1 * jjp1 * llm, ndexu)
 
   !  Geopotentiel
 
-  CALL histwrite(histaveid, 'phi', itau_w, phi,  &
-       iip1*jjp1*llm, ndexu)
+  CALL histwrite(histaveid, 'phi', itau_w, phi, &
+          iip1 * jjp1 * llm, ndexu)
 
   !  Traceurs
@@ -109,10 +109,10 @@
   !  Masse
 
-  CALL histwrite(histaveid, 'masse', itau_w, masse,  &
-       iip1*jjp1*llm, ndexu)
+  CALL histwrite(histaveid, 'masse', itau_w, masse, &
+          iip1 * jjp1 * llm, ndexu)
 
   !  Pression au sol
 
-  CALL histwrite(histaveid, 'ps', itau_w, ps, iip1*jjp1, ndex2d)
+  CALL histwrite(histaveid, 'ps', itau_w, ps, iip1 * jjp1, ndex2d)
 
   ! Geopotentiel au sol
@@ -121,7 +121,7 @@
 
   IF (ok_sync) THEN
-     CALL histsync(histaveid)
-     CALL histsync(histvaveid)
-     CALL histsync(histuaveid)
+    CALL histsync(histaveid)
+    CALL histsync(histvaveid)
+    CALL histsync(histuaveid)
   ENDIF
 
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writehist.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writehist.f90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writehist.f90	(revision 5118)
@@ -8,4 +8,5 @@
   USE temps_mod, ONLY: itau_dyn
   USE lmdz_description, ONLY: descript
+  USE lmdz_iniprint, ONLY: lunout, prt_level
 
   IMPLICIT NONE
@@ -36,5 +37,4 @@
   include "paramet.h"
   include "comgeom.h"
-  include "iniprint.h"
 
   !
