Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/academic.h
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/academic.h	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/academic.h	(revision 5116)
@@ -3,7 +3,7 @@
 
       common/academic/tetarappel,knewt_t,kfrict,knewt_g,clat4
-      real :: tetarappel(ip1jmp1,llm)
-      real :: knewt_t(llm)
-      real :: kfrict(llm)
-      real :: knewt_g
-      real :: clat4(ip1jmp1)
+      REAL :: tetarappel(ip1jmp1,llm)
+      REAL :: knewt_t(llm)
+      REAL :: kfrict(llm)
+      REAL :: knewt_g
+      REAL :: clat4(ip1jmp1)
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/adaptdt.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/adaptdt.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/adaptdt.f90	(revision 5116)
@@ -43,6 +43,6 @@
     dtbon=dtvr/n
 
-   return
-end subroutine adaptdt
+   RETURN
+END SUBROUTINE adaptdt
 
 
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advn.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advn.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advn.f90	(revision 5116)
@@ -25,6 +25,6 @@
   !   Arguments:
   !   ----------
-  integer :: mode
-  real :: masse(ip1jmp1,llm)
+  INTEGER :: mode
+  REAL :: masse(ip1jmp1,llm)
   REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
   REAL :: q(ip1jmp1,llm)
@@ -35,21 +35,21 @@
   !
   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
+  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
+  REAL :: zzpbar,zzw
+
+  REAL :: qmin,qmax
   data qmin,qmax/0.,1./
   data temps1,temps2,temps3/0.,0.,0./
@@ -139,5 +139,5 @@
   !   Arguments:
   !   ----------
-  real :: q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm)
+  REAL :: q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm)
   !
   !  Local
@@ -146,9 +146,9 @@
   INTEGER :: ij,l
   !
-  real :: dxqu(ip1jmp1),zqu(ip1jmp1)
-  real :: zqmax(ip1jmp1),zqmin(ip1jmp1)
+  REAL :: dxqu(ip1jmp1),zqu(ip1jmp1)
+  REAL :: zqmax(ip1jmp1),zqmin(ip1jmp1)
   logical :: extremum(ip1jmp1)
 
-  integer :: mode
+  INTEGER :: mode
   save mode
   data mode/1/
@@ -156,5 +156,5 @@
   !   calcul des pentes en u:
   !   -----------------------
-  if (mode==0) then
+  if (mode==0) THEN
      do l=1,llm
         do ij=1,ip1jm
@@ -206,5 +206,5 @@
      enddo
      do ij=iip2+1,ip1jm
-        if(extremum(ij)) then
+        IF(extremum(ij)) THEN
            qg(ij,l)=q(ij,l)
            qd(ij,l)=q(ij,l)
@@ -222,5 +222,5 @@
 
      do ij=iip2+1,ip1jm
-        if(extremum(ij).and..not.extremum(ij-1)) &
+        IF(extremum(ij).and..not.extremum(ij-1)) &
               qd(ij-1,l)=q(ij,l)
      enddo
@@ -256,5 +256,5 @@
   !   Arguments:
   !   ----------
-  real :: q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm)
+  REAL :: q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm)
   !
   !  Local
@@ -263,13 +263,13 @@
   INTEGER :: ij,l
   !
-  real :: dyqv(ip1jm),zqv(ip1jm,llm)
-  real :: zqmax(ip1jm),zqmin(ip1jm)
+  REAL :: dyqv(ip1jm),zqv(ip1jm,llm)
+  REAL :: zqmax(ip1jm),zqmin(ip1jm)
   logical :: extremum(ip1jmp1)
 
-  integer :: mode
+  INTEGER :: mode
   save mode
   data mode/1/
 
-  if (mode==0) then
+  if (mode==0) THEN
      do l=1,llm
         do ij=1,ip1jmp1
@@ -315,5 +315,5 @@
 
      do ij=iip2,ip1jm
-        if(extremum(ij)) then
+        IF(extremum(ij)) THEN
            qs(ij,l)=q(ij,l)
            qn(ij,l)=q(ij,l)
@@ -352,5 +352,5 @@
   !   Arguments:
   !   ----------
-  real :: q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm)
+  REAL :: q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm)
   !
   !  Local
@@ -359,9 +359,9 @@
   INTEGER :: ij,l
   !
-  real :: dzqw(ip1jmp1,llm+1),zqw(ip1jmp1,llm+1)
-  real :: zqmax(ip1jmp1,llm),zqmin(ip1jmp1,llm)
+  REAL :: dzqw(ip1jmp1,llm+1),zqw(ip1jmp1,llm+1)
+  REAL :: zqmax(ip1jmp1,llm),zqmin(ip1jmp1,llm)
   logical :: extremum(ip1jmp1,llm)
 
-  integer :: mode
+  INTEGER :: mode
   save mode
 
@@ -371,5 +371,5 @@
   !   -----------------------
 
-  if (mode==0) then
+  if (mode==0) THEN
      do l=1,llm
         do ij=1,ip1jmp1
@@ -424,5 +424,5 @@
   do l=2,llm-1
      do ij=1,ip1jmp1
-        if(extremum(ij,l)) then
+        IF(extremum(ij,l)) THEN
            qh(ij,l)=q(ij,l)
            qb(ij,l)=q(ij,l)
@@ -435,5 +435,5 @@
   ! do l=2,llm-1
   !    do ij=1,ip1jmp1
-  !       if(extremum(ij,l)) then
+  !       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)
@@ -474,8 +474,8 @@
   !   Arguments:
   !   ----------
-  integer :: mode
-  real :: masse(ip1jmp1,llm)
-  real :: u_m( ip1jmp1,llm )
-  real :: q(ip1jmp1,llm),qd(ip1jmp1,llm),qg(ip1jmp1,llm)
+  INTEGER :: mode
+  REAL :: masse(ip1jmp1,llm)
+  REAL :: u_m( ip1jmp1,llm )
+  REAL :: q(ip1jmp1,llm),qd(ip1jmp1,llm),qg(ip1jmp1,llm)
   !
   !  Local
@@ -483,15 +483,15 @@
   !
   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
+  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
+  REAL :: prec
   save prec
 
@@ -501,5 +501,5 @@
         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
+           ! 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)
@@ -507,9 +507,9 @@
            !    qg(ij,l)=q(ij,l)
            ! endif
-           if(abs(zdq)>prec) then
+           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
+              ! 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)
@@ -530,5 +530,5 @@
    do l=1,llm
    do ij=iip2,ip1jm-1
-      if (u_m(ij,l)>=0.) then
+      if (u_m(ij,l)>=0.) THEN
          zsigp=zsigd(ij,l)
          zsigm=zsigg(ij,l)
@@ -548,14 +548,14 @@
       ladvplus(ij,l)=zu>zm
       zsig=zu/zm
-      if(zsig==0.) zsigp=0.1
-      if (mode==1) then
-         if (zsig<=zsigp) then
+      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
+         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
+         if (zsig<=zsigp) THEN
              u_mq(ij,l)=u_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
          else
@@ -565,5 +565,5 @@
          endif
       endif
-      ! if(zsig.lt.0.) then
+      ! IF(zsig.lt.0.) THEN
       !    PRINT*,'au point ij=',ij,'  l=',l,'  sig=',zsig
       !    stop
@@ -587,5 +587,5 @@
      nl(l)=0
      do ij=iip2,ip1jm
-        if(ladvplus(ij,l)) then
+        IF(ladvplus(ij,l)) THEN
            nl(l)=nl(l)+1
            u_mq(ij,l)=0.
@@ -595,5 +595,5 @@
   enddo
 
-  if(n0>1) then
+  IF(n0>1) THEN
   IF (prt_level > 9) WRITE(lunout,*) &
         'Nombre de points pour lesquels on advect plus que le' &
@@ -601,9 +601,9 @@
 
      do l=1,llm
-        if(nl(l)>0) then
+        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
+              IF(ladvplus(ij,l).and.mod(ij,iip1)/=0) THEN
                  iju=iju+1
                  indu(iju)=ij
@@ -619,5 +619,5 @@
               zu_m=u_m(ij,l)
               u_mq(ij,l)=0.
-              if(zu_m>0.) then
+              IF(zu_m>0.) THEN
                  ijq=ij
                  i=ijq-(j-1)*iip1
@@ -632,5 +632,5 @@
   !   ajout de la maille non completement advectee
          zsig=zu_m/masse(ijq,l)
-         if(zsig<=zsigd(ijq,l)) then
+         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)))
@@ -639,5 +639,5 @@
       ! goto 8888
             zz=0.5*(zsig-zsigd(ijq,l))/zsigg(ijq,l)
-            if(.not.(zz>0..and.zz<=0.5)) then
+            IF(.not.(zz>0..and.zz<=0.5)) THEN
                  WRITE(lunout,*)'probleme2 au point ij=',ij, &
                        '  l=',l
@@ -662,5 +662,5 @@
   ! 2eme MODIF SPECIFIQUE
          zsig=-zu_m/masse(ij+1,l)
-         if(zsig<=zsigg(ijq,l)) then
+         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)))
@@ -669,5 +669,5 @@
         ! goto 9999
             zz=0.5*(zsig-zsigg(ijq,l))/zsigd(ijq,l)
-            if(.not.(zz>0..and.zz<=0.5)) then
+            IF(.not.(zz>0..and.zz<=0.5)) THEN
                  WRITE(lunout,*)'probleme22 au point ij=',ij &
                        ,'  l=',l
@@ -736,7 +736,7 @@
   !   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
@@ -745,11 +745,11 @@
   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
+  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
   save prec
 
@@ -758,5 +758,5 @@
         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
+           ! 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)
@@ -764,9 +764,9 @@
            !    qs(ij,l)=q(ij,l)
            ! endif
-           if(abs(zdq)>prec) then
+           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
+              ! 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)
@@ -782,5 +782,5 @@
 
    do ij=1,ip1jm
-      if (v_m(ij,l)>=0.) then
+      if (v_m(ij,l)>=0.) THEN
          zsigp=zsign(ij+iip1)
          zsigm=zsigs(ij+iip1)
@@ -798,6 +798,6 @@
       endif
       zsig=abs(v_m(ij,l))/zm
-      if(zsig==0.) zsigp=0.1
-      if (zsig<=zsigp) then
+      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))
       else
@@ -863,7 +863,7 @@
   !   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)
 
   !
@@ -873,9 +873,9 @@
   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
+  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
 
@@ -885,5 +885,5 @@
         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
+           ! 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)
@@ -892,5 +892,5 @@
            ! endif
 
-           if(abs(zdq)>prec) then
+           IF(abs(zdq)>prec) THEN
               zsigb(ij,l)=(q(ij,l)-qh(ij,l))/zdq
               zsigh(ij,l)=1.-zsigb(ij,l)
@@ -907,5 +907,5 @@
    do l=2,llm
    do ij=1,ip1jmp1
-      if (w_m(ij,l)>=0.) then
+      if (w_m(ij,l)>=0.) THEN
          zsigp=zsigb(ij,l)
          zsigm=zsigh(ij,l)
@@ -923,6 +923,6 @@
       endif
       zsig=abs(w_m(ij,l))/zm
-      if(zsig==0.) zsigp=0.1
-      if (zsig<=zsigp) then
+      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))
       else
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advy.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advy.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advy.f90	(revision 5116)
@@ -91,5 +91,5 @@
   REAL :: cy1(llm,ntra), cyLAT(llm,ntra)
   REAL :: z1(iim), zcos(iim), zsin(iim)
-  real :: smpn,smps,s0pn,s0ps
+  REAL :: smpn,smps,s0pn,s0ps
   REAL :: SSUM
   EXTERNAL SSUM
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/com_io_dyn_mod.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/com_io_dyn_mod.F90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/com_io_dyn_mod.F90	(revision 5116)
@@ -8,24 +8,24 @@
 ! Names of various files for outputs (in the dynamics)
   ! to store instantaneous values:
-  character(len=18),parameter :: dynhist_file="dyn_hist.nc" ! on scalar grid
-  character(len=18),parameter :: dynhistv_file="dyn_histv.nc" ! on v grid
-  character(len=18),parameter :: dynhistu_file="dyn_histu.nc" ! on u grid
+  CHARACTER(LEN=18),parameter :: dynhist_file="dyn_hist.nc" ! on scalar grid
+  CHARACTER(LEN=18),parameter :: dynhistv_file="dyn_histv.nc" ! on v grid
+  CHARACTER(LEN=18),parameter :: dynhistu_file="dyn_histu.nc" ! on u grid
 
   ! to store averaged values:
-  character(len=18),parameter :: dynhistave_file="dyn_hist_ave.nc"
-  character(len=18),parameter :: dynhistvave_file="dyn_histv_ave.nc"
-  character(len=18),parameter :: dynhistuave_file="dyn_histu_ave.nc"
+  CHARACTER(LEN=18),parameter :: dynhistave_file="dyn_hist_ave.nc"
+  CHARACTER(LEN=18),parameter :: dynhistvave_file="dyn_histv_ave.nc"
+  CHARACTER(LEN=18),parameter :: dynhistuave_file="dyn_histu_ave.nc"
   
 ! Ids of various files for outputs (in the dynamics)
 
   ! instantaneous (these are set by inithist.F)
-  integer :: histid
-  integer :: histvid
-  integer :: histuid
+  INTEGER :: histid
+  INTEGER :: histvid
+  INTEGER :: histuid
   
   ! averages (these are set by initdynav.F)
-  integer :: histaveid
-  integer :: histvaveid
-  integer :: histuaveid
+  INTEGER :: histaveid
+  INTEGER :: histvaveid
+  INTEGER :: histuaveid
   
 end module com_io_dyn_mod
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/coordij.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/coordij.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/coordij.f90	(revision 5116)
@@ -22,5 +22,5 @@
   include "comgeom.h"
 
-  real :: zlon,zlat
+  REAL :: zlon,zlat
 
   zlon=lon*pi/180.
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/diagedyn.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/diagedyn.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/diagedyn.f90	(revision 5116)
@@ -138,5 +138,5 @@
   INTEGER :: ndiag     ! max number of diagnostic in parallel
   PARAMETER (ndiag=10)
-  integer :: pas(ndiag)
+  INTEGER :: pas(ndiag)
   save pas
   data pas/ndiag*0/
@@ -160,5 +160,5 @@
   !
   PRINT*,'MAIS POURQUOI DONC DIAGEDYN NE MARCHE PAS ?'
-  return
+  RETURN
   ! On ne garde les donnees que dans les colonnes i=1,iim
   DO jj = 1,jjp1
@@ -325,5 +325,5 @@
   !
   ELSE
-    write(lunout,*)'diagedyn: set to function with Earth parameters'
+    WRITE(lunout,*)'diagedyn: set to function with Earth parameters'
   ENDIF ! of if (planet_type=="earth")
   RETURN
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert.F90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert.F90	(revision 5116)
@@ -3,7 +3,7 @@
 SUBROUTINE disvert()
 
-  use ioipsl, only: getin
-  use new_unit_m, only: new_unit
-  use lmdz_assert, only: assert
+  use ioipsl, ONLY: getin
+  use new_unit_m, ONLY: new_unit
+  use lmdz_assert, ONLY: assert
   USE comvert_mod, ONLY: ap, bp, aps, bps, nivsigs, nivsig, dpres, presnivs, &
                          pseudoalt, pa, preff, scaleheight, presinter
@@ -47,7 +47,7 @@
   REAL alpha, beta, deltaz
   REAL x
-  character(len=*),parameter :: modname="disvert"
-
-  character(len=24):: vert_sampling
+  CHARACTER(LEN=*),parameter :: modname="disvert"
+
+  CHARACTER(LEN=24):: vert_sampling
   ! (allowed values are "param", "tropo", "strato" and "read")
 
@@ -62,5 +62,5 @@
   CALL getin('vert_sampling', vert_sampling)
   WRITE(lunout,*) TRIM(modname)//' vert_sampling = ' // vert_sampling
-  if (llm==39 .and. vert_sampling=="strato") then
+  if (llm==39 .and. vert_sampling=="strato") THEN
      dsigmin=0.3 ! Vieille option par défaut pour CMIP5
   else
@@ -82,5 +82,5 @@
      CLOSE(99)
      alpha=deltaz/(llm*scaleheight)
-     write(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
+     WRITE(lunout, *)trim(modname),':scaleheight, alpha, k0, k1, beta', &
                                scaleheight, alpha, k0, k1, beta
 
@@ -96,5 +96,5 @@
         dzk1=alpha*tanh(l/k0)
         dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
-        write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
+        WRITE(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
         zkm1=zk
      enddo
@@ -332,13 +332,13 @@
   ENDDO
 
-  write(lunout, *)  trim(modname),': BP '
-  write(lunout, *) bp
-  write(lunout, *)  trim(modname),': AP '
-  write(lunout, *) ap
-
-  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
-  write(lunout, *)'couches calcules pour une pression de surface =', preff
-  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de '
-  write(lunout, *) scaleheight,' km'
+  WRITE(lunout, *)  trim(modname),': BP '
+  WRITE(lunout, *) bp
+  WRITE(lunout, *)  trim(modname),': AP '
+  WRITE(lunout, *) ap
+
+  WRITE(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
+  WRITE(lunout, *)'couches calcules pour une pression de surface =', preff
+  WRITE(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de '
+  WRITE(lunout, *) scaleheight,' km'
   DO l = 1, llm
      dpres(l) = bp(l) - bp(l+1)
@@ -347,5 +347,5 @@
      presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
      pseudoalt(l) = log(preff/presnivs(l))*scaleheight
-     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
+     WRITE(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
           pseudoalt(l) &
           , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
@@ -354,9 +354,9 @@
   DO l=1, llmp1
      presinter(l)= ( ap(l)+bp(l)*preff)
-     write(lunout, *)'PRESINTER(', l, ')=', presinter(l)
+     WRITE(lunout, *)'PRESINTER(', l, ')=', presinter(l)
   ENDDO
 
-  write(lunout, *) trim(modname),': PRESNIVS '
-  write(lunout, *) presnivs
+  WRITE(lunout, *) trim(modname),': PRESNIVS '
+  WRITE(lunout, *) presnivs
 
 CONTAINS
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert_noterre.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert_noterre.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert_noterre.f90	(revision 5116)
@@ -30,20 +30,20 @@
   REAL :: snorm
   REAL :: alpha,beta,gama,delta,deltaz
-  real :: quoi,quand
+  REAL :: quoi,quand
   REAL :: zsig(llm),sig(llm+1)
   INTEGER :: np,ierr
-  integer :: ierr1,ierr2,ierr3,ierr4
+  INTEGER :: ierr1,ierr2,ierr3,ierr4
   REAL :: x
 
   REAL :: SSUM
   EXTERNAL SSUM
-  real :: newsig
+  REAL :: newsig
   REAL :: dz0,dz1,nhaut,sig1,esig,csig,zz
-  real :: tt,rr,gg, prevz
-  real :: s(llm),dsig(llm)
-
-  integer :: iz
-  real :: z, ps,p
-  character(len=*),parameter :: modname="disvert_noterre"
+  REAL :: tt,rr,gg, prevz
+  REAL :: s(llm),dsig(llm)
+
+  INTEGER :: iz
+  REAL :: z, ps,p
+  CHARACTER(LEN=*),parameter :: modname="disvert_noterre"
 
   !
@@ -55,5 +55,5 @@
   hybrid=.TRUE. ! default value for hybrid (ie: use hybrid coordinates)
   CALL getin('hybrid',hybrid)
-  write(lunout,*) trim(modname),': hybrid=',hybrid
+  WRITE(lunout,*) trim(modname),': hybrid=',hybrid
 
   ! Ouverture possible de fichiers typiquement E.T.
@@ -61,5 +61,5 @@
      open(99,file="esasig.def",status='old',form='formatted', &
            iostat=ierr2)
-     if(ierr2/=0) then
+     IF(ierr2/=0) THEN
           close(99)
           open(99,file="z2sig.def",status='old',form='formatted', &
@@ -71,6 +71,5 @@
   !   ----------------------------------------
 
-  IF(ierr2==0) then
-
+  IF(ierr2==0) THEN
      ! Lecture de esasig.def :
      ! Systeme peu souple, mais qui respecte en theorie
@@ -78,7 +77,7 @@
      ! <-> energie cinetique, d'apres la note de Frederic Hourdin...
 
-     write(lunout,*)'*****************************'
-     write(lunout,*)'WARNING reading esasig.def'
-     write(lunout,*)'*****************************'
+     WRITE(lunout,*)'*****************************'
+     WRITE(lunout,*)'WARNING reading esasig.def'
+     WRITE(lunout,*)'*****************************'
      READ(99,*) scaleheight
      READ(99,*) dz0
@@ -126,8 +125,8 @@
   !   ----------------------------------------
 
-  ELSE IF(ierr4==0) then
-     write(lunout,*)'****************************'
-     write(lunout,*)'Reading z2sig.def'
-     write(lunout,*)'****************************'
+  ELSE IF(ierr4==0) THEN
+     WRITE(lunout,*)'****************************'
+     WRITE(lunout,*)'Reading z2sig.def'
+     WRITE(lunout,*)'****************************'
 
      READ(99,*) scaleheight
@@ -146,6 +145,6 @@
   !-----------------------------------------------------------------------
   ELSE
-     write(lunout,*) 'didn t you forget something ??? '
-     write(lunout,*) 'We need file  z2sig.def ! (OR esasig.def)'
+     WRITE(lunout,*) 'didn t you forget something ??? '
+     WRITE(lunout,*) 'We need file  z2sig.def ! (OR esasig.def)'
      stop
   ENDIF
@@ -170,7 +169,7 @@
 
   if (hybrid) then  ! use hybrid coordinates
-     write(lunout,*) "*********************************"
-     write(lunout,*) "Using hybrid vertical coordinates"
-     write(lunout,*)
+     WRITE(lunout,*) "*********************************"
+     WRITE(lunout,*) "Using hybrid vertical coordinates"
+     WRITE(lunout,*)
      ! Coordonnees hybrides avec mod
      DO l = 1, llm
@@ -183,7 +182,7 @@
      ap(llmp1) = 0.
   else ! use sigma coordinates
-     write(lunout,*) "********************************"
-     write(lunout,*) "Using sigma vertical coordinates"
-     write(lunout,*)
+     WRITE(lunout,*) "********************************"
+     WRITE(lunout,*) "Using sigma vertical coordinates"
+     WRITE(lunout,*)
      ! Pour ne pas passer en coordonnees hybrides
      DO l = 1, llm
@@ -196,8 +195,8 @@
   bp(llmp1) =   0.
 
-  write(lunout,*) trim(modname),': BP '
-  write(lunout,*)  bp
-  write(lunout,*) trim(modname),': AP '
-  write(lunout,*)  ap
+  WRITE(lunout,*) trim(modname),': BP '
+  WRITE(lunout,*)  bp
+  WRITE(lunout,*) trim(modname),': AP '
+  WRITE(lunout,*)  ap
 
   ! Calcul au milieu des couches :
@@ -214,5 +213,5 @@
   ENDDO
 
-  if (hybrid) then
+  if (hybrid) THEN
      aps(llm) = aps(llm-1)**2 / aps(llm-2)
      bps(llm) = 0.5*(bp(llm) + bp(llm+1))
@@ -222,8 +221,8 @@
   end if
 
-  write(lunout,*) trim(modname),': BPs '
-  write(lunout,*)  bps
-  write(lunout,*) trim(modname),': APs'
-  write(lunout,*)  aps
+  WRITE(lunout,*) trim(modname),': BPs '
+  WRITE(lunout,*)  bps
+  WRITE(lunout,*) trim(modname),': APs'
+  WRITE(lunout,*)  aps
 
   DO l = 1, llm
@@ -232,9 +231,9 @@
   ENDDO
 
-  write(lunout,*)trim(modname),' : PRESNIVS'
-  write(lunout,*)presnivs
-  write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ', &
+  WRITE(lunout,*)trim(modname),' : PRESNIVS'
+  WRITE(lunout,*)presnivs
+  WRITE(lunout,*)'Pseudo altitude of Presnivs : (for a scale ', &
         'height of ',scaleheight,' km)'
-  write(lunout,*)pseudoalt
+  WRITE(lunout,*)pseudoalt
 
   ! --------------------------------------------------
@@ -252,5 +251,5 @@
   !   do l=2,llm
   ! approximation of scale height for Venus
-  !      if (zsig(l-1).le.55.) then
+  !      if (zsig(l-1).le.55.) THEN
   !         scaleheight = 15.5 - zsig(l-1)/55.*10.
   !      else
@@ -260,5 +259,5 @@
   !    .    log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps))
   !   END DO
-  !   write(53,'(I3,50F10.5)') iz, zsig
+  !   WRITE(53,'(I3,50F10.5)') iz, zsig
   !  END DO
   !  close(53)
@@ -296,17 +295,17 @@
 
   IMPLICIT NONE
-  real :: x1, x2, sig,pa,preff, newsig, F
-  integer :: j
+  REAL :: x1, x2, sig,pa,preff, newsig, F
+  INTEGER :: j
 
   newsig = sig
   x1=0
   x2=1
-  if (sig>=1) then
+  if (sig>=1) THEN
         newsig= sig
-  else if (sig*preff/pa>=0.25) then
+  else if (sig*preff/pa>=0.25) THEN
     DO J=1,9999  ! nombre d''iteration max
       F=((1 -pa/preff)*exp(1-1./newsig**2)+(pa/preff)*newsig)/sig
-      ! write(0,*) J, ' newsig =', newsig, ' F= ', F
-      if (F>1) then
+      ! WRITE(0,*) J, ' newsig =', newsig, ' F= ', F
+      if (F>1) THEN
           X2 = newsig
           newsig=(X1+newsig)*0.5
@@ -319,5 +318,5 @@
       IF(abs(10.*log(F))<1.E-5) goto 999
     END DO
-   else   !    if (sig*preff/pa.le.0.25) then
+   else   !    if (sig*preff/pa.le.0.25) THEN
          newsig= sig*preff/pa
    end if
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/exner_hyb_m.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/exner_hyb_m.F90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/exner_hyb_m.F90	(revision 5116)
@@ -54,15 +54,15 @@
 
     logical,save :: firstcall=.TRUE.
-    character(len=*),parameter :: modname="exner_hyb"
+    CHARACTER(LEN=*),parameter :: modname="exner_hyb"
 
     ! Sanity check
-    if (firstcall) then
+    if (firstcall) THEN
        ! sanity checks for Shallow Water case (1 vertical layer)
-       if (llm==1) then
-          if (kappa/=1) then
+       if (llm==1) THEN
+          if (kappa/=1) THEN
              CALL abort_gcm(modname, &
                   "kappa!=1 , but running in Shallow Water mode!!",42)
           endif
-          if (cpp/=r) then
+          if (cpp/=r) THEN
              CALL abort_gcm(modname, &
                   "cpp!=r , but running in Shallow Water mode!!",42)
@@ -74,6 +74,5 @@
 
     ! Specific behaviour for Shallow Water (1 vertical layer) case:
-    if (llm==1) then
-
+    if (llm==1) THEN
        ! Compute pks(:),pk(:),pkf(:)
 
@@ -83,5 +82,5 @@
        ENDDO
 
-       if (present(pkf)) then
+       if (present(pkf)) THEN
           pkf = pk
           CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
@@ -89,5 +88,5 @@
 
        ! our work is done, exit routine
-       return
+       RETURN
     endif ! of if (llm.eq.1)
 
@@ -138,5 +137,5 @@
     ENDDO
 
-    if (present(pkf)) then
+    if (present(pkf)) THEN
        !    calcul de pkf
        pkf = pk
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/exner_milieu_m.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/exner_milieu_m.F90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/exner_milieu_m.F90	(revision 5116)
@@ -51,15 +51,15 @@
 
     logical,save :: firstcall=.TRUE.
-    character(len=*),parameter :: modname="exner_milieu"
+    CHARACTER(LEN=*),parameter :: modname="exner_milieu"
 
     ! Sanity check
-    if (firstcall) then
+    if (firstcall) THEN
        ! sanity checks for Shallow Water case (1 vertical layer)
-       if (llm==1) then
-          if (kappa/=1) then
+       if (llm==1) THEN
+          if (kappa/=1) THEN
              CALL abort_gcm(modname, &
                   "kappa!=1 , but running in Shallow Water mode!!",42)
           endif
-          if (cpp/=r) then
+          if (cpp/=r) THEN
              CALL abort_gcm(modname, &
                   "cpp!=r , but running in Shallow Water mode!!",42)
@@ -71,6 +71,5 @@
 
     ! Specific behaviour for Shallow Water (1 vertical layer) case:
-    if (llm==1) then
-
+    if (llm==1) THEN
        ! Compute pks(:),pk(:),pkf(:)
 
@@ -80,5 +79,5 @@
        ENDDO
 
-       if (present(pkf)) then
+       if (present(pkf)) THEN
           pkf = pk
           CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
@@ -86,5 +85,5 @@
 
        ! our work is done, exit routine
-       return
+       RETURN
     endif ! of if (llm.eq.1)
 
@@ -117,5 +116,5 @@
     ENDDO
 
-    if (present(pkf)) then
+    if (present(pkf)) THEN
        !    calcul de pkf
        pkf = pk
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/fxhyp_m.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/fxhyp_m.F90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/fxhyp_m.F90	(revision 5116)
@@ -18,9 +18,9 @@
     ! 1., taux=0., clon=0.) est à - 180 degrés.
 
-    use lmdz_arth, only: arth
-    use invert_zoom_x_m, only: invert_zoom_x, nmax
-    use nrtype, only: pi, pi_d, twopi, twopi_d, k8
-    use principal_cshift_m, only: principal_cshift
-    use serre_mod, only: clon, grossismx, dzoomx, taux
+    use lmdz_arth, ONLY: arth
+    use invert_zoom_x_m, ONLY: invert_zoom_x, nmax
+    use nrtype, ONLY: pi, pi_d, twopi, twopi_d, k8
+    use principal_cshift_m, ONLY: principal_cshift
+    use serre_mod, ONLY: clon, grossismx, dzoomx, taux
 
     include "dimensions.h"
@@ -45,5 +45,5 @@
     print *, "Call sequence information: fxhyp"
 
-    test_iim: if (iim==1) then
+    test_iim: if (iim==1) THEN
        rlonv(1)=0.
        rlonu(1)=pi
@@ -56,5 +56,5 @@
        xprimp025(:)=1.
     else test_iim
-       test_grossismx: if (grossismx == 1.) then
+       test_grossismx: if (grossismx == 1.) THEN
           step = twopi / iim
 
@@ -193,5 +193,5 @@
              END DO
 
-             if (rlonm025(is2) < - pi) then
+             if (rlonm025(is2) < - pi) THEN
                 print *, 'Rlonm025 plus petit que - pi !'
                 STOP 1
@@ -204,5 +204,5 @@
              END DO
 
-             if (rlonm025(is2) > pi) then
+             if (rlonm025(is2) > pi) THEN
                 print *, 'Rlonm025 plus grand que pi !'
                 STOP 1
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/fyhyp_m.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/fyhyp_m.F90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/fyhyp_m.F90	(revision 5116)
@@ -16,7 +16,7 @@
     ! Il vaut mieux avoir : grossismy * dzoom < pi / 2
 
-    use lmdz_coefpoly, only: coefpoly
-    use nrtype, only: k8
-    use serre_mod, only: clat, grossismy, dzoomy, tauy
+    use lmdz_coefpoly, ONLY: coefpoly
+    use nrtype, ONLY: k8
+    use serre_mod, ONLY: clat, grossismy, dzoomy, tauy
 
     include "dimensions.h"
@@ -207,5 +207,5 @@
              yprimin = a1 + 2.*a2*yi + 3.*a3*yi2
           END DO
-          if (abs(yi-yo1) > epsilon) then
+          if (abs(yi-yo1) > epsilon) THEN
              print *, 'Pas de solution.', j, ylon2
              STOP 1
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/gr_int_dyn.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/gr_int_dyn.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/gr_int_dyn.f90	(revision 5116)
@@ -12,10 +12,10 @@
 
   INTEGER :: iim
-  integer :: ip1, jp1
+  INTEGER :: ip1, jp1
   REAL :: champin(iim, jp1)
   REAL :: champdyn(iim+1, jp1)
 
   INTEGER :: i, j
-  real :: polenord, polesud
+  REAL :: polenord, polesud
 
   !-----------------------------------------------------------------------
@@ -34,7 +34,7 @@
   do j = 1, jp1
     do i = 1, iim
-      if (j == 1) then
+      if (j == 1) THEN
         champdyn(i, j) = polenord
-      else if (j == jp1) then
+      else if (j == jp1) THEN
         champdyn(i, j) = polesud
       else
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90	(revision 5116)
@@ -226,5 +226,5 @@
   ! CLOSE netcdf file
   CALL ncclos(ncid_out,rcode_out)
-  write(*,*) "END grilles_gcm_netcdf_sub OK"
+  WRITE(*,*) "END grilles_gcm_netcdf_sub OK"
 
 END SUBROUTINE grilles_gcm_netcdf_sub
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iniconst.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iniconst.F90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iniconst.F90	(revision 5116)
@@ -21,6 +21,6 @@
   include "iniprint.h"
 
-  character(len=*),parameter :: modname="iniconst"
-  character(len=80) :: abort_message
+  CHARACTER(LEN=*),parameter :: modname="iniconst"
+  CHARACTER(LEN=80) :: abort_message
 
 
@@ -48,10 +48,10 @@
   r       = cpp * kappa
 
-  write(lunout,*) trim(modname),': R  CP  Kappa ',r,cpp,kappa
+  WRITE(lunout,*) trim(modname),': R  CP  Kappa ',r,cpp,kappa
 
   !-----------------------------------------------------------------------
 
   ! vertical discretization: default behavior depends on planet_type flag
-  if (planet_type=="earth") then
+  if (planet_type=="earth") THEN
      disvert_type=1
   else
@@ -60,17 +60,17 @@
   ! 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
+  WRITE(lunout,*) trim(modname),': disvert_type=',disvert_type
 
   pressure_exner = disvert_type == 1 ! default value
   CALL getin('pressure_exner', pressure_exner)
 
-  if (disvert_type==1) then
+  if (disvert_type==1) THEN
      ! standard case for Earth (automatic generation of levels)
      CALL disvert()
-  else if (disvert_type==2) then
+  else if (disvert_type==2) THEN
      ! standard case for planets (levels generated using z2sig.def file)
      CALL disvert_noterre
   else
-     write(abort_message,*) "Wrong value for disvert_type: ", disvert_type
+     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 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inidissip.F90	(revision 5116)
@@ -1,7 +1,6 @@
-
 ! $Id$
 
-SUBROUTINE inidissip( lstardis,nitergdiv,nitergrot,niterh  , &
-     tetagdiv,tetagrot,tetatemp, vert_prof_dissip)
+SUBROUTINE inidissip(lstardis, nitergdiv, nitergrot, niterh, &
+        tetagdiv, tetagrot, tetatemp, vert_prof_dissip)
   !=======================================================================
   !   initialisation de la dissipation horizontale
@@ -11,9 +10,10 @@
   !   -------------
 
-  USE control_mod, ONLY: dissip_period,iperiod
+  USE control_mod, ONLY: dissip_period, iperiod
   USE comconst_mod, ONLY: dissip_deltaz, dissip_factz, dissip_zref, &
-                          dtdiss, dtvr, rad
+          dtdiss, dtvr, rad
   USE comvert_mod, ONLY: preff, presnivs
   USE lmdz_filtreg, ONLY: filtreg
+  USE lmdz_libmath, ONLY: minmax
 
   IMPLICIT NONE
@@ -23,9 +23,9 @@
   include "iniprint.h"
 
-  LOGICAL,INTENT(in) :: lstardis
-  INTEGER,INTENT(in) :: nitergdiv,nitergrot,niterh
-  REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
-
-  integer, INTENT(in):: vert_prof_dissip
+  LOGICAL, INTENT(in) :: lstardis
+  INTEGER, INTENT(in) :: nitergdiv, nitergrot, niterh
+  REAL, INTENT(in) :: tetagdiv, tetagrot, tetatemp
+
+  integer, INTENT(in) :: vert_prof_dissip
   ! Vertical profile of horizontal dissipation
   ! Allowed values:
@@ -33,15 +33,15 @@
   ! 1: tanh of altitude
 
-! Local variables:
-  REAL fact,zvert(llm),zz
-  REAL zh(ip1jmp1),zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1)
-  real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1,llm)
-  REAL ullm,vllm,umin,vmin,zhmin,zhmax
+  ! Local variables:
+  REAL fact, zvert(llm), zz
+  REAL zh(ip1jmp1), zu(ip1jmp1), gx(ip1jmp1), divgra(ip1jmp1)
+  real zv(ip1jm), gy(ip1jm), deltap(ip1jmp1, llm)
+  REAL ullm, vllm, umin, vmin, zhmin, zhmax
   REAL zllm
 
-  INTEGER l,ij,idum,ii
+  INTEGER l, ij, idum, ii
   REAL tetamin
   REAL pseudoz
-  character (len=80) :: abort_message
+  character (len = 80) :: abort_message
 
   REAL ran1
@@ -53,7 +53,7 @@
   !   -----------------------------------------------------------------
 
-  crot     = -1.
-  cdivu    = -1.
-  cdivh    = -1.
+  crot = -1.
+  cdivu = -1.
+  cdivh = -1.
 
   !   calcul de la valeur propre de divgrad:
@@ -61,99 +61,99 @@
   idum = 0
   DO l = 1, llm
-     DO ij = 1, ip1jmp1
-        deltap(ij,l) = 1.
-     ENDDO
-  ENDDO
-
-  idum  = -1
-  zh(1) = RAN1(idum)-.5
-  idum  = 0
+    DO ij = 1, ip1jmp1
+      deltap(ij, l) = 1.
+    ENDDO
+  ENDDO
+
+  idum = -1
+  zh(1) = RAN1(idum) - .5
+  idum = 0
   DO ij = 2, ip1jmp1
-     zh(ij) = RAN1(idum) -.5
-  ENDDO
-
-  CALL filtreg (zh,jjp1,1,2,1,.TRUE.,1)
-
-  CALL minmax(iip1*jjp1,zh,zhmin,zhmax )
-
-  IF ( zhmin >= zhmax  )     THEN
-     write(lunout,*)'  Inidissip  zh min max  ',zhmin,zhmax
-     abort_message='probleme generateur alleatoire dans inidissip'
-     CALL abort_gcm('inidissip',abort_message,1)
+    zh(ij) = RAN1(idum) - .5
+  ENDDO
+
+  CALL filtreg (zh, jjp1, 1, 2, 1, .TRUE., 1)
+
+  CALL minmax(iip1 * jjp1, zh, zhmin, zhmax)
+
+  IF (zhmin >= zhmax)     THEN
+    WRITE(lunout, *)'  Inidissip  zh min max  ', zhmin, zhmax
+    abort_message = 'probleme generateur alleatoire dans inidissip'
+    CALL abort_gcm('inidissip', abort_message, 1)
   ENDIF
 
-  zllm = ABS( zhmax )
-  DO l = 1,50
-     IF(lstardis) THEN
-        CALL divgrad2(1,zh,deltap,niterh,divgra)
-     ELSE
-        CALL divgrad (1,zh,niterh,divgra)
-     ENDIF
-
-     zllm  = ABS(maxval(divgra))
-     zh = divgra / zllm
+  zllm = ABS(zhmax)
+  DO l = 1, 50
+    IF(lstardis) THEN
+      CALL divgrad2(1, zh, deltap, niterh, divgra)
+    ELSE
+      CALL divgrad (1, zh, niterh, divgra)
+    ENDIF
+
+    zllm = ABS(maxval(divgra))
+    zh = divgra / zllm
   ENDDO
 
   IF(lstardis) THEN
-     cdivh = 1./ zllm
+    cdivh = 1. / zllm
   ELSE
-     cdivh = zllm ** ( -1./niterh )
+    cdivh = zllm ** (-1. / niterh)
   ENDIF
 
   !   calcul des valeurs propres de gradiv (ii =1) et  nxgrarot(ii=2)
   !   -----------------------------------------------------------------
-  write(lunout,*)'inidissip: calcul des valeurs propres'
+  WRITE(lunout, *)'inidissip: calcul des valeurs propres'
 
   DO    ii = 1, 2
 
-     DO ij = 1, ip1jmp1
-        zu(ij)  = RAN1(idum) -.5
-     ENDDO
-     CALL filtreg (zu,jjp1,1,2,1,.TRUE.,1)
-     DO ij = 1, ip1jm
-        zv(ij) = RAN1(idum) -.5
-     ENDDO
-     CALL filtreg (zv,jjm,1,2,1,.FALSE.,1)
-
-     CALL minmax(iip1*jjp1,zu,umin,ullm )
-     CALL minmax(iip1*jjm, zv,vmin,vllm )
-
-     ullm = ABS ( ullm )
-     vllm = ABS ( vllm )
-
-     DO    l = 1, 50
-        IF(ii==1) THEN
-           !cccc             CALL covcont( 1,zu,zv,zu,zv )
-           IF(lstardis) THEN
-              CALL gradiv2( 1,zu,zv,nitergdiv,gx,gy )
-           ELSE
-              CALL gradiv ( 1,zu,zv,nitergdiv,gx,gy )
-           ENDIF
+    DO ij = 1, ip1jmp1
+      zu(ij) = RAN1(idum) - .5
+    ENDDO
+    CALL filtreg (zu, jjp1, 1, 2, 1, .TRUE., 1)
+    DO ij = 1, ip1jm
+      zv(ij) = RAN1(idum) - .5
+    ENDDO
+    CALL filtreg (zv, jjm, 1, 2, 1, .FALSE., 1)
+
+    CALL minmax(iip1 * jjp1, zu, umin, ullm)
+    CALL minmax(iip1 * jjm, zv, vmin, vllm)
+
+    ullm = ABS (ullm)
+    vllm = ABS (vllm)
+
+    DO    l = 1, 50
+      IF(ii==1) THEN
+        !cccc             CALL covcont( 1,zu,zv,zu,zv )
+        IF(lstardis) THEN
+          CALL gradiv2(1, zu, zv, nitergdiv, gx, gy)
         ELSE
-           IF(lstardis) THEN
-              CALL nxgraro2( 1,zu,zv,nitergrot,gx,gy )
-           ELSE
-              CALL nxgrarot( 1,zu,zv,nitergrot,gx,gy )
-           ENDIF
+          CALL gradiv (1, zu, zv, nitergdiv, gx, gy)
         ENDIF
-
-        zllm = max(abs(maxval(gx)), abs(maxval(gy)))
-        zu = gx / zllm
-        zv = gy / zllm
-     end DO
-
-     IF ( ii==1 ) THEN
+      ELSE
         IF(lstardis) THEN
-           cdivu  = 1./zllm
+          CALL nxgraro2(1, zu, zv, nitergrot, gx, gy)
         ELSE
-           cdivu  = zllm **( -1./nitergdiv )
+          CALL nxgrarot(1, zu, zv, nitergrot, gx, gy)
         ENDIF
-     ELSE
-        IF(lstardis) THEN
-           crot   = 1./ zllm
-        ELSE
-           crot   = zllm **( -1./nitergrot )
-        ENDIF
-     ENDIF
+      ENDIF
+
+      zllm = max(abs(maxval(gx)), abs(maxval(gy)))
+      zu = gx / zllm
+      zv = gy / zllm
+    end DO
+
+    IF (ii==1) THEN
+      IF(lstardis) THEN
+        cdivu = 1. / zllm
+      ELSE
+        cdivu = zllm **(-1. / nitergdiv)
+      ENDIF
+    ELSE
+      IF(lstardis) THEN
+        crot = 1. / zllm
+      ELSE
+        crot = zllm **(-1. / nitergrot)
+      ENDIF
+    ENDIF
 
   end DO
@@ -163,9 +163,9 @@
 
   !     IF(.NOT.lstardis) THEN
-  fact    = rad*24./REAL(jjm)
-  fact    = fact*fact
-  write(lunout,*)'inidissip: coef u ', fact/cdivu, 1./cdivu
-  write(lunout,*)'inidissip: coef r ', fact/crot , 1./crot
-  write(lunout,*)'inidissip: coef h ', fact/cdivh, 1./cdivh
+  fact = rad * 24. / REAL(jjm)
+  fact = fact * fact
+  WRITE(lunout, *)'inidissip: coef u ', fact / cdivu, 1. / cdivu
+  WRITE(lunout, *)'inidissip: coef r ', fact / crot, 1. / crot
+  WRITE(lunout, *)'inidissip: coef h ', fact / cdivh, 1. / cdivh
   !     ENDIF
 
@@ -174,50 +174,49 @@
   !   --------------------------------------------------
 
-  if (vert_prof_dissip == 1) then
-     do l=1,llm
-        pseudoz=8.*log(preff/presnivs(l))
-        zvert(l)=1+ &
-             (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. &
-             *(dissip_factz-1.)
-     enddo
+  if (vert_prof_dissip == 1) THEN
+    do l = 1, llm
+      pseudoz = 8. * log(preff / presnivs(l))
+      zvert(l) = 1 + &
+              (tanh((pseudoz - dissip_zref) / dissip_deltaz) + 1.) / 2. &
+                      * (dissip_factz - 1.)
+    enddo
   else
-     DO l=1,llm
-        zvert(l)=1.
-     ENDDO
-     fact=2.
-     DO l = 1, llm
-        zz      = 1. - preff/presnivs(l)
-        zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
-     ENDDO
+    DO l = 1, llm
+      zvert(l) = 1.
+    ENDDO
+    fact = 2.
+    DO l = 1, llm
+      zz = 1. - preff / presnivs(l)
+      zvert(l) = fact - (fact - 1.) / (1. + zz * zz)
+    ENDDO
   endif
 
-
-  write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale'
-
-  tetamin =  1.e+6
-
-  DO l=1,llm
-     tetaudiv(l)   = zvert(l)/tetagdiv
-     tetaurot(l)   = zvert(l)/tetagrot
-     tetah(l)      = zvert(l)/tetatemp
-
-     IF( tetamin> (1./tetaudiv(l)) ) tetamin = 1./ tetaudiv(l)
-     IF( tetamin> (1./tetaurot(l)) ) tetamin = 1./ tetaurot(l)
-     IF( tetamin> (1./   tetah(l)) ) tetamin = 1./    tetah(l)
+  WRITE(lunout, *)'inidissip: Constantes de temps de la diffusion horizontale'
+
+  tetamin = 1.e+6
+
+  DO l = 1, llm
+    tetaudiv(l) = zvert(l) / tetagdiv
+    tetaurot(l) = zvert(l) / tetagrot
+    tetah(l) = zvert(l) / tetatemp
+
+    IF(tetamin> (1. / tetaudiv(l))) tetamin = 1. / tetaudiv(l)
+    IF(tetamin> (1. / tetaurot(l))) tetamin = 1. / tetaurot(l)
+    IF(tetamin> (1. / tetah(l))) tetamin = 1. / tetah(l)
   ENDDO
 
   ! If dissip_period=0 calculate value for dissipation period, else keep value read from gcm.def
   IF (dissip_period == 0) THEN
-     dissip_period = INT( tetamin/( 2.*dtvr*iperiod) ) * iperiod
-     write(lunout,*)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ',tetamin,dtvr,iperiod,dissip_period
-     dissip_period = MAX(iperiod,dissip_period)
+    dissip_period = INT(tetamin / (2. * dtvr * iperiod)) * iperiod
+    WRITE(lunout, *)'inidissip: tetamin dtvr iperiod dissip_period(intermed) ', tetamin, dtvr, iperiod, dissip_period
+    dissip_period = MAX(iperiod, dissip_period)
   END IF
 
-  dtdiss  = dissip_period * dtvr
-  write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr
-
-  DO l = 1,llm
-     write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &
-          dtdiss*tetah(l)
+  dtdiss = dissip_period * dtvr
+  WRITE(lunout, *)'inidissip: dissip_period=', dissip_period, ' dtdiss=', dtdiss, ' dtvr=', dtvr
+
+  DO l = 1, llm
+    WRITE(lunout, *)zvert(l), dtdiss * tetaudiv(l), dtdiss * tetaurot(l), &
+            dtdiss * tetah(l)
   ENDDO
 
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inigeom.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inigeom.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inigeom.f90	(revision 5116)
@@ -16,6 +16,6 @@
   !
   !
-  use fxhyp_m, only: fxhyp
-  use fyhyp_m, only: fyhyp
+  use fxhyp_m, ONLY: fxhyp
+  use fyhyp_m, ONLY: fyhyp
   USE comconst_mod, ONLY: pi, g, omeg, rad
   USE logic_mod, ONLY: fxyhypb, ysinus
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inigrads.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inigrads.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inigrads.f90	(revision 5116)
@@ -7,15 +7,15 @@
   IMPLICIT NONE
 
-  integer :: if, im, jm, lm, i, j, l
-  real :: x(im), y(jm), z(lm), fx, fy, fz, dt
-  real :: xmin, xmax, ymin, ymax
+  INTEGER :: if, im, jm, lm, i, j, l
+  REAL :: x(im), y(jm), z(lm), fx, fy, fz, dt
+  REAL :: xmin, xmax, ymin, ymax
 
-  character(len = *), intent(in) :: file
-  character(len = *), intent(in) :: titlel
+  CHARACTER(LEN = *), intent(in) :: file
+  CHARACTER(LEN = *), intent(in) :: titlel
 
   INCLUDE "gradsdef.h"
 
   ! data unit/66,32,34,36,38,40,42,44,46,48/
-  integer :: nf
+  INTEGER :: nf
   save nf
   data nf/0/
@@ -49,6 +49,6 @@
   do i = 1, im
     xd(i, if) = x(i) * fx
-    if(xd(i, if)<xmin) iid(if) = i + 1
-    if(xd(i, if)<=xmax) ifd(if) = i
+    IF(xd(i, if)<xmin) iid(if) = i + 1
+    IF(xd(i, if)<=xmax) ifd(if) = i
   enddo
   PRINT*, 'On stoke du point ', iid(if), '  a ', ifd(if), ' en x'
@@ -59,6 +59,6 @@
   do j = 1, jm
     yd(j, if) = y(j) * fy
-    if(yd(j, if)>ymax) jid(if) = j + 1
-    if(yd(j, if)>=ymin) jfd(if) = j
+    IF(yd(j, if)>ymax) jid(if) = j + 1
+    IF(yd(j, if)>=ymin) jfd(if) = j
   enddo
   PRINT*, 'On stoke du point ', jid(if), '  a ', jfd(if), ' en y'
@@ -89,3 +89,3 @@
 
 
-end subroutine inigrads
+END SUBROUTINE inigrads
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/initfluxsto.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/initfluxsto.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/initfluxsto.f90	(revision 5116)
@@ -47,21 +47,21 @@
   !   Arguments
   !
-  character(len = *) :: infile
-  real :: tstep, t_ops, t_wrt
-  integer :: fileid, filevid, filedid
+  CHARACTER(LEN = *) :: infile
+  REAL :: tstep, t_ops, t_wrt
+  INTEGER :: fileid, filevid, filedid
 
   ! This routine needs IOIPSL to work
   !   Variables locales
   !
-  real :: nivd(1)
-  integer :: tau0
-  real :: zjulian
-  character(len = 3) :: str
-  character(len = 10) :: ctrac
-  integer :: iq
-  real :: rlong(iip1, jjp1), rlat(iip1, jjp1), rl(1, 1)
-  integer :: uhoriid, vhoriid, thoriid, zvertiid, dhoriid, dvertiid
-  integer :: ii, jj
-  integer :: zan, idayref
+  REAL :: nivd(1)
+  INTEGER :: tau0
+  REAL :: zjulian
+  CHARACTER(LEN = 3) :: str
+  CHARACTER(LEN = 10) :: ctrac
+  INTEGER :: iq
+  REAL :: rlong(iip1, jjp1), rlat(iip1, jjp1), rl(1, 1)
+  INTEGER :: uhoriid, vhoriid, thoriid, zvertiid, dhoriid, dvertiid
+  INTEGER :: ii, jj
+  INTEGER :: zan, idayref
   logical :: ok_sync
   !
@@ -211,5 +211,5 @@
   CALL histend(filevid)
   CALL histend(filedid)
-  if (ok_sync) then
+  if (ok_sync) THEN
     CALL histsync(fileid)
     CALL histsync(filevid)
@@ -217,4 +217,4 @@
   endif
 
-  return
-end subroutine initfluxsto
+  RETURN
+END SUBROUTINE initfluxsto
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inithist.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inithist.F90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inithist.F90	(revision 5116)
@@ -46,17 +46,17 @@
   !   Arguments
   !
-  integer :: day0, anne0
-  real :: tstep, t_ops, t_wrt
+  INTEGER :: day0, anne0
+  REAL :: tstep, t_ops, t_wrt
 
   ! This routine needs IOIPSL to work
   !   Variables locales
   !
-  integer :: tau0
-  real :: zjulian
-  integer :: iq
-  real :: rlong(iip1, jjp1), rlat(iip1, jjp1)
-  integer :: uhoriid, vhoriid, thoriid, zvertiid
-  integer :: ii, jj
-  integer :: zan, dayref
+  INTEGER :: tau0
+  REAL :: zjulian
+  INTEGER :: iq
+  REAL :: rlong(iip1, jjp1), rlat(iip1, jjp1)
+  INTEGER :: uhoriid, vhoriid, thoriid, zvertiid
+  INTEGER :: ii, jj
+  INTEGER :: zan, dayref
   !
   !  Initialisations
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inter_barxy_m.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inter_barxy_m.F90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inter_barxy_m.F90	(revision 5116)
@@ -1,3 +1,2 @@
-
 ! $Id$
 
@@ -15,6 +14,6 @@
   SUBROUTINE inter_barxy(dlonid, dlatid, champ, rlonimod, rlatimod, champint)
 
-    use lmdz_assert_eq, only: assert_eq
-    use lmdz_assert, only: assert
+    use lmdz_assert_eq, ONLY: assert_eq
+    use lmdz_assert, ONLY: assert
 
     include "dimensions.h"
@@ -27,13 +26,13 @@
     ! (for "aire", "apoln", "apols")
 
-    REAL, intent(in):: dlonid(:)
+    REAL, intent(in) :: dlonid(:)
     ! (longitude from input file, in rad, from -pi to pi)
 
-    REAL, intent(in):: dlatid(:), champ(:, :), rlonimod(:)
-
-    REAL, intent(in):: rlatimod(:)
+    REAL, intent(in) :: dlatid(:), champ(:, :), rlonimod(:)
+
+    REAL, intent(in) :: rlatimod(:)
     ! (latitude angle, in degrees or rad, in strictly decreasing order)
 
-    real, intent(out):: champint(:, :)
+    real, intent(out) :: champint(:, :)
     ! Si taille de la seconde dim = jjm + 1, on veut interpoler sur les
     ! jjm+1 latitudes rlatu du modele (latitudes des scalaires et de U)
@@ -55,9 +54,9 @@
 
     jnterfd = assert_eq(size(champ, 2) - 1, size(dlatid), &
-         "inter_barxy jnterfd")
+            "inter_barxy jnterfd")
     jmods = size(champint, 2)
     CALL assert(size(champ, 1) == size(dlonid), "inter_barxy size(champ, 1)")
     CALL assert((/size(rlonimod), size(champint, 1)/) == iim, &
-         "inter_barxy iim")
+            "inter_barxy iim")
     CALL assert(any(jmods == (/jjm, jjm + 1/)), 'inter_barxy jmods')
     CALL assert(size(rlatimod) == jjm, "inter_barxy size(rlatimod)")
@@ -65,16 +64,16 @@
     ! Check decreasing order for "rlatimod":
     DO i = 2, jjm
-       IF (rlatimod(i) >= rlatimod(i-1)) stop &
-            '"inter_barxy": "rlatimod" should be strictly decreasing'
+      IF (rlatimod(i) >= rlatimod(i - 1)) stop &
+              '"inter_barxy": "rlatimod" should be strictly decreasing'
     ENDDO
 
     yjmod(:jjm) = ord_coordm(rlatimod)
     IF (jmods == jjm + 1) THEN
-       IF (90. - yjmod(jjm) < 0.01) stop &
-            '"inter_barxy": with jmods = jjm + 1, yjmod(jjm) should be < 90.'
+      IF (90. - yjmod(jjm) < 0.01) stop &
+              '"inter_barxy": with jmods = jjm + 1, yjmod(jjm) should be < 90.'
     ELSE
-       ! jmods = jjm
-       IF (ABS(yjmod(jjm) - 90.) > 0.01) stop &
-            '"inter_barxy": with jmods = jjm, yjmod(jjm) should be 90.'
+      ! jmods = jjm
+      IF (ABS(yjmod(jjm) - 90.) > 0.01) stop &
+              '"inter_barxy": with jmods = jjm, yjmod(jjm) should be 90.'
     ENDIF
 
@@ -82,19 +81,19 @@
 
     DO j = 1, jnterfd + 1
-       champy(:, j) = inter_barx(dlonid, champ(:, j), rlonimod)
-    ENDDO
-
-    CALL ord_coord(dlatid, yjdat, decrois) 
+      champy(:, j) = inter_barx(dlonid, champ(:, j), rlonimod)
+    ENDDO
+
+    CALL ord_coord(dlatid, yjdat, decrois)
     IF (decrois) champy(:, :) = champy(:, jnterfd + 1:1:-1)
     DO i = 1, iim
-       champint(i, :) = inter_bary(yjdat, champy(i, :), yjmod)
+      champint(i, :) = inter_bary(yjdat, champy(i, :), yjmod)
     ENDDO
     champint(:, :) = champint(:, jmods:1:-1)
 
     IF (jmods == jjm + 1) THEN
-       ! Valeurs uniques aux poles
-       champint(:, 1) = SUM(aire(:iim,  1) * champint(:, 1)) / apoln
-       champint(:, jjm + 1) = SUM(aire(:iim, jjm + 1) &
-            * champint(:, jjm + 1)) / apols
+      ! Valeurs uniques aux poles
+      champint(:, 1) = SUM(aire(:iim, 1) * champint(:, 1)) / apoln
+      champint(:, jjm + 1) = SUM(aire(:iim, jjm + 1) &
+              * champint(:, jjm + 1)) / apols
     ENDIF
 
@@ -103,5 +102,5 @@
   !******************************
 
-  function inter_barx(dlonid, fdat, rlonimod) 
+  function inter_barx(dlonid, fdat, rlonimod)
 
     !        INTERPOLATION BARYCENTRIQUE BASEE SUR LES AIRES
@@ -117,11 +116,12 @@
     !      ( Les abscisses sont exprimees en degres)
 
-    use lmdz_assert_eq, only: assert_eq
+    USE lmdz_assert_eq, ONLY: assert_eq
+    USE lmdz_libmath, ONLY: minmax
 
     IMPLICIT NONE
 
-    REAL, intent(in):: dlonid(:)
-    real, intent(in):: fdat(:)
-    real, intent(in):: rlonimod(:)
+    REAL, intent(in) :: dlonid(:)
+    real, intent(in) :: fdat(:)
+    real, intent(in) :: rlonimod(:)
 
     real inter_barx(size(rlonimod))
@@ -130,6 +130,6 @@
 
     INTEGER idatmax, imodmax
-    REAL xxid(size(dlonid)+1), xxd(size(dlonid)+1), fdd(size(dlonid)+1)
-    REAL  fxd(size(dlonid)+1), xchan(size(dlonid)+1), fdchan(size(dlonid)+1) 
+    REAL xxid(size(dlonid) + 1), xxd(size(dlonid) + 1), fdd(size(dlonid) + 1)
+    REAL  fxd(size(dlonid) + 1), xchan(size(dlonid) + 1), fdchan(size(dlonid) + 1)
     REAL  xxim(size(rlonimod))
 
@@ -149,12 +149,12 @@
     !    A L'INTERFACE OUEST DE LA PREMIERE MAILLE DU MODELE  
     DO imod = 1, imodmax
-       xxim(imod) = rlonimod(imod)
-    ENDDO
-
-    CALL minmax( imodmax, xxim, chmin, chmax)
-    IF( chmax<6.50 )   THEN
-       DO imod = 1, imodmax
-          xxim(imod) = xxim(imod) * 180./pi
-       ENDDO
+      xxim(imod) = rlonimod(imod)
+    ENDDO
+
+    CALL minmax(imodmax, xxim, chmin, chmax)
+    IF(chmax<6.50)   THEN
+      DO imod = 1, imodmax
+        xxim(imod) = xxim(imod) * 180. / pi
+      ENDDO
     ENDIF
 
@@ -162,45 +162,45 @@
 
     DO imod = 1, imodmax
-       xxim(imod) = xxim(imod) - xim0
-    ENDDO
-
-    idatmax1 = idatmax +1
+      xxim(imod) = xxim(imod) - xim0
+    ENDDO
+
+    idatmax1 = idatmax + 1
 
     DO idat = 1, idatmax
-       xxd(idat) = dlonid(idat)
-    ENDDO
-
-    CALL minmax( idatmax, xxd, chmin, chmax)
-    IF( chmax<6.50 )  THEN
-       DO idat = 1, idatmax
-          xxd(idat) = xxd(idat) * 180./pi
-       ENDDO
+      xxd(idat) = dlonid(idat)
+    ENDDO
+
+    CALL minmax(idatmax, xxd, chmin, chmax)
+    IF(chmax<6.50)  THEN
+      DO idat = 1, idatmax
+        xxd(idat) = xxd(idat) * 180. / pi
+      ENDDO
     ENDIF
 
     DO idat = 1, idatmax
-       xxd(idat) = MOD( xxd(idat) - xim0, 360. )
-       fdd(idat) = fdat (idat)
+      xxd(idat) = MOD(xxd(idat) - xim0, 360.)
+      fdd(idat) = fdat (idat)
     ENDDO
 
     i = 2
-    DO while (xxd(i) >= xxd(i-1) .and. i < idatmax)
-       i = i + 1
-    ENDDO
-    IF (xxd(i) < xxd(i-1)) THEN
-       ichang = i
-       !  ***  reorganisation  des longitudes entre 0. et 360. degres ****
-       nid = idatmax - ichang +1
-       DO i = 1, nid
-          xchan (i) = xxd(i+ichang -1 )
-          fdchan(i) = fdd(i+ichang -1 )
-       ENDDO
-       DO i=1, ichang -1
-          xchan (i+ nid) = xxd(i)
-          fdchan(i+nid) = fdd(i) 
-       ENDDO
-       DO i =1, idatmax
-          xxd(i) = xchan(i)
-          fdd(i) = fdchan(i)
-       ENDDO
+    DO while (xxd(i) >= xxd(i - 1) .and. i < idatmax)
+      i = i + 1
+    ENDDO
+    IF (xxd(i) < xxd(i - 1)) THEN
+      ichang = i
+      !  ***  reorganisation  des longitudes entre 0. et 360. degres ****
+      nid = idatmax - ichang + 1
+      DO i = 1, nid
+        xchan (i) = xxd(i + ichang - 1)
+        fdchan(i) = fdd(i + ichang - 1)
+      ENDDO
+      DO i = 1, ichang - 1
+        xchan (i + nid) = xxd(i)
+        fdchan(i + nid) = fdd(i)
+      ENDDO
+      DO i = 1, idatmax
+        xxd(i) = xchan(i)
+        fdd(i) = fdchan(i)
+      ENDDO
     end IF
 
@@ -213,39 +213,39 @@
 
     DO idat = 1, idatmax
-       IF ( xxd( idatmax1- idat )<360.) exit
-       id1 = id1 + 1
+      IF (xxd(idatmax1 - idat)<360.) exit
+      id1 = id1 + 1
     ENDDO
 
     DO idat = 1, idatmax
-       IF (xxd(idat)>0.) exit
-       id0 = id0 + 1
+      IF (xxd(idat)>0.) exit
+      id0 = id0 + 1
     END DO
 
-    IF( id1 /= 0 ) then
-       DO idat = 1, id1
-          xxid(idat) = xxd(idatmax - id1 + idat) - 360.
-          fxd (idat) = fdd(idatmax - id1 + idat)     
-       END DO
-       DO idat = 1, idatmax - id1
-          xxid(idat + id1) = xxd(idat)
-          fxd (idat + id1) = fdd(idat)
-       END DO
+    IF(id1 /= 0) THEN
+      DO idat = 1, id1
+        xxid(idat) = xxd(idatmax - id1 + idat) - 360.
+        fxd (idat) = fdd(idatmax - id1 + idat)
+      END DO
+      DO idat = 1, idatmax - id1
+        xxid(idat + id1) = xxd(idat)
+        fxd (idat + id1) = fdd(idat)
+      END DO
     end IF
 
-    IF(id0 /= 0) then
-       DO idat = 1, idatmax - id0
-          xxid(idat) = xxd(idat + id0)
-          fxd (idat) = fdd(idat + id0)
-       END DO
-
-       DO idat = 1, id0
-          xxid (idatmax - id0 + idat) =  xxd(idat) + 360.
-          fxd  (idatmax - id0 + idat) =  fdd(idat)   
-       END DO
-    else 
-       DO idat = 1, idatmax
-          xxid(idat)  = xxd(idat)
-          fxd (idat)  = fdd(idat)
-       ENDDO
+    IF(id0 /= 0) THEN
+      DO idat = 1, idatmax - id0
+        xxid(idat) = xxd(idat + id0)
+        fxd (idat) = fdd(idat + id0)
+      END DO
+
+      DO idat = 1, id0
+        xxid (idatmax - id0 + idat) = xxd(idat) + 360.
+        fxd  (idatmax - id0 + idat) = fdd(idat)
+      END DO
+    else
+      DO idat = 1, idatmax
+        xxid(idat) = xxd(idat)
+        fxd (idat) = fdd(idat)
+      ENDDO
     end IF
     xxid(idatmax1) = xxid(1) + 360.
@@ -258,33 +258,33 @@
     ! iteration
 
-    x0   = xim0
-    dxm  = 0.
+    x0 = xim0
+    dxm = 0.
     imod = 1
     idat = 1
 
     do while (imod <= imodmax)
-       do while (xxim(imod)>xxid(idat))
-          dx   = xxid(idat) - x0
-          dxm  = dxm + dx
-          inter_barx(imod) = inter_barx(imod) + dx * fxd(idat)
-          x0   = xxid(idat)
-          idat = idat + 1
-       END DO
-       IF (xxim(imod)<xxid(idat)) THEN
-          dx   = xxim(imod) - x0
-          dxm  = dxm + dx
-          inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm
-          x0   = xxim(imod)
-          dxm  = 0.
-          imod = imod + 1
-       ELSE
-          dx   = xxim(imod) - x0
-          dxm  = dxm + dx
-          inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm
-          x0   = xxim(imod)
-          dxm  = 0.
-          imod = imod + 1
-          idat = idat + 1
-       END IF
+      do while (xxim(imod)>xxid(idat))
+        dx = xxid(idat) - x0
+        dxm = dxm + dx
+        inter_barx(imod) = inter_barx(imod) + dx * fxd(idat)
+        x0 = xxid(idat)
+        idat = idat + 1
+      END DO
+      IF (xxim(imod)<xxid(idat)) THEN
+        dx = xxim(imod) - x0
+        dxm = dxm + dx
+        inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm
+        x0 = xxim(imod)
+        dxm = 0.
+        imod = imod + 1
+      ELSE
+        dx = xxim(imod) - x0
+        dxm = dxm + dx
+        inter_barx(imod) = (inter_barx(imod) + dx * fxd(idat)) / dxm
+        x0 = xxim(imod)
+        dxm = 0.
+        imod = imod + 1
+        idat = idat + 1
+      END IF
     END DO
 
@@ -299,15 +299,15 @@
     ! L'indice 1 correspond à l'interface maille 1 -- maille 2.
 
-    use lmdz_assert, only: assert
+    use lmdz_assert, ONLY: assert
 
     IMPLICIT NONE
 
-    REAL, intent(in):: yjdat(:)
+    REAL, intent(in) :: yjdat(:)
     ! (angles, ordonnées des interfaces des mailles des données, in
     ! degrees, in increasing order)
 
-    REAL, intent(in):: fdat(:) ! champ de données
-
-    REAL, intent(in):: yjmod(:)
+    REAL, intent(in) :: fdat(:) ! champ de données
+
+    REAL, intent(in) :: yjmod(:)
     ! (ordonnées des interfaces des mailles du modèle)
     ! (in degrees, in strictly increasing order)
@@ -317,5 +317,5 @@
     ! Variables local to the procedure:
 
-    REAL y0, dy, dym 
+    REAL y0, dy, dym
     INTEGER jdat ! indice du champ de données
     integer jmod ! indice du champ du modèle
@@ -327,34 +327,34 @@
     ! Initialisation des variables
     inter_bary(:) = 0.
-    y0    = -90.
-    dym   = 0.
-    jmod  = 1
-    jdat  = 1
+    y0 = -90.
+    dym = 0.
+    jmod = 1
+    jdat = 1
 
     do while (jmod <= size(yjmod))
-       do while (yjmod(jmod) > yjdat(jdat))
-          dy         = yjdat(jdat) - y0
-          dym        = dym + dy
-          inter_bary(jmod) = inter_bary(jmod) + dy * fdat(jdat)
-          y0         = yjdat(jdat)
-          jdat       = jdat + 1
-       END DO
-       IF (yjmod(jmod) < yjdat(jdat)) THEN
-          dy         = yjmod(jmod) - y0
-          dym        = dym + dy
-          inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym
-          y0         = yjmod(jmod)
-          dym        = 0.
-          jmod       = jmod + 1
-       ELSE
-          ! {yjmod(jmod) == yjdat(jdat)}
-          dy         = yjmod(jmod) - y0
-          dym        = dym + dy
-          inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym
-          y0         = yjmod(jmod)
-          dym        = 0.
-          jmod       = jmod + 1
-          jdat       = jdat + 1
-       END IF
+      do while (yjmod(jmod) > yjdat(jdat))
+        dy = yjdat(jdat) - y0
+        dym = dym + dy
+        inter_bary(jmod) = inter_bary(jmod) + dy * fdat(jdat)
+        y0 = yjdat(jdat)
+        jdat = jdat + 1
+      END DO
+      IF (yjmod(jmod) < yjdat(jdat)) THEN
+        dy = yjmod(jmod) - y0
+        dym = dym + dy
+        inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym
+        y0 = yjmod(jmod)
+        dym = 0.
+        jmod = jmod + 1
+      ELSE
+        ! {yjmod(jmod) == yjdat(jdat)}
+        dy = yjmod(jmod) - y0
+        dym = dym + dy
+        inter_bary(jmod) = (inter_bary(jmod) + dy * fdat(jdat)) / dym
+        y0 = yjmod(jmod)
+        dym = 0.
+        jmod = jmod + 1
+        jdat = jdat + 1
+      END IF
     END DO
     ! Le test de fin suppose que l'interface 0 est commune aux deux
@@ -373,10 +373,10 @@
     ! Finally, the procedure adds 90° as the last value of the array.
 
-    use lmdz_assert_eq, only: assert_eq
-    use comconst_mod, only: pi
+    use lmdz_assert_eq, ONLY: assert_eq
+    use comconst_mod, ONLY: pi
 
     IMPLICIT NONE
 
-    REAL, intent(in):: xi(:)
+    REAL, intent(in) :: xi(:)
     ! (latitude, in degrees or radians, in increasing or decreasing order)
     ! ("xi" should contain latitudes from pole to pole.
@@ -385,6 +385,6 @@
     ! So the extreme values should not be 90° and -90°.)
 
-    REAL, intent(out):: xo(:) ! angles in degrees
-    LOGICAL, intent(out):: decrois
+    REAL, intent(out) :: xo(:) ! angles in degrees
+    LOGICAL, intent(out) :: decrois
 
     ! Variables  local to the procedure:
@@ -398,21 +398,21 @@
     decrois = xi(2) < xi(1)
     DO i = 3, nmax
-       IF (decrois .neqv. xi(i) < xi(i-1)) stop &
-            '"ord_coord":  latitudes are not monotonic'
-    ENDDO
-
-    IF (abs(xi(1)) < pi) then
-       ! "xi" contains latitudes in radians
-       xo(:nmax) = xi(:) * 180. / pi
+      IF (decrois .neqv. xi(i) < xi(i - 1)) stop &
+              '"ord_coord":  latitudes are not monotonic'
+    ENDDO
+
+    IF (abs(xi(1)) < pi) THEN
+      ! "xi" contains latitudes in radians
+      xo(:nmax) = xi(:) * 180. / pi
     else
-       ! "xi" contains latitudes in degrees
-       xo(:nmax) = xi(:)
+      ! "xi" contains latitudes in degrees
+      xo(:nmax) = xi(:)
     end IF
 
     IF (ABS(abs(xo(1)) - 90) < 0.001 .or. ABS(abs(xo(nmax)) - 90) < 0.001) THEN
-       print *, "ord_coord"
-       PRINT *, '"xi" should contain the latitudes of the boundaries of ' &
-            // 'grid cells, not the centers of grid cells.'
-       STOP
+      print *, "ord_coord"
+      PRINT *, '"xi" should contain the latitudes of the boundaries of ' &
+              // 'grid cells, not the centers of grid cells.'
+      STOP
     ENDIF
 
@@ -429,9 +429,9 @@
     ! order.
 
-    use comconst_mod, only: pi
+    use comconst_mod, ONLY: pi
 
     IMPLICIT NONE
 
-    REAL, intent(in):: xi(:) ! angle, in rad or degrees
+    REAL, intent(in) :: xi(:) ! angle, in rad or degrees
     REAL ord_coordm(size(xi)) ! angle, in degrees
 
@@ -439,9 +439,9 @@
 
     IF (xi(1) < 6.5) THEN
-       ! "xi" is in rad
-       ord_coordm(:) = xi(size(xi):1:-1) * 180. / pi
+      ! "xi" is in rad
+      ord_coordm(:) = xi(size(xi):1:-1) * 180. / pi
     else
-       ! "xi" is in degrees
-       ord_coordm(:) = xi(size(xi):1:-1)
+      ! "xi" is in degrees
+      ord_coordm(:) = xi(size(xi):1:-1)
     ENDIF
 
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/interpost.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/interpost.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/interpost.f90	(revision 5116)
@@ -12,8 +12,8 @@
 
   ! Arguments
-  real :: q(iip1,jjp1,llm)
-  real :: qppm(iim,jjp1,llm)
+  REAL :: q(iip1,jjp1,llm)
+  REAL :: qppm(iim,jjp1,llm)
   ! Local
-  integer :: l,i,j
+  INTEGER :: l,i,j
 
   ! RE-INVERSION DES NIVEAUX
@@ -41,3 +41,3 @@
    return
 
-end subroutine interpost
+END SUBROUTINE interpost
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/interpre.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/interpre.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/interpre.f90	(revision 5116)
@@ -18,24 +18,24 @@
   !---------------------------------------------------
   ! Arguments
-  real :: apppm(llm + 1), bpppm(llm + 1)
-  real :: q(iip1, jjp1, llm), qppm(iim, jjp1, llm)
+  REAL :: apppm(llm + 1), bpppm(llm + 1)
+  REAL :: q(iip1, jjp1, llm), qppm(iim, jjp1, llm)
   !---------------------------------------------------
-  real :: masse(iip1, jjp1, llm)
-  real :: massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)
-  real :: w(iip1, jjp1, llm)
-  real :: fluxwppm(iim, jjp1, llm)
-  real :: pbaru(iip1, jjp1, llm)
-  real :: pbarv(iip1, jjm, llm)
-  real :: unatppm(iim, jjp1, llm)
-  real :: vnatppm(iim, jjp1, llm)
-  real :: psppm(iim, jjp1)
+  REAL :: masse(iip1, jjp1, llm)
+  REAL :: massebx(iip1, jjp1, llm), masseby(iip1, jjm, llm)
+  REAL :: w(iip1, jjp1, llm)
+  REAL :: fluxwppm(iim, jjp1, llm)
+  REAL :: pbaru(iip1, jjp1, llm)
+  REAL :: pbarv(iip1, jjm, llm)
+  REAL :: unatppm(iim, jjp1, llm)
+  REAL :: vnatppm(iim, jjp1, llm)
+  REAL :: psppm(iim, jjp1)
   !---------------------------------------------------
   ! Local
-  real :: vnat(iip1, jjp1, llm)
-  real :: unat(iip1, jjp1, llm)
-  real :: fluxw(iip1, jjp1, llm)
-  real :: smass(iip1, jjp1)
+  REAL :: vnat(iip1, jjp1, llm)
+  REAL :: unat(iip1, jjp1, llm)
+  REAL :: fluxw(iip1, jjp1, llm)
+  REAL :: smass(iip1, jjp1)
   !----------------------------------------------------
-  integer :: l, ij, i, j
+  INTEGER :: l, ij, i, j
 
   ! CALCUL DE LA PRESSION DE SURFACE
@@ -118,6 +118,6 @@
   enddo
 
-  return
-end subroutine interpre
+  RETURN
+END SUBROUTINE interpre
 
 
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/invert_zoom_x_m.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/invert_zoom_x_m.F90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/invert_zoom_x_m.F90	(revision 5116)
@@ -9,7 +9,7 @@
   SUBROUTINE invert_zoom_x(xf, xtild, Xprimt, xlon, xprimm, xuv)
 
-    use lmdz_coefpoly, only: coefpoly
-    use nrtype, only: pi, pi_d, twopi_d, k8
-    use serre_mod, only: clon
+    use lmdz_coefpoly, ONLY: coefpoly
+    use nrtype, ONLY: pi, pi_d, twopi_d, k8
+    use serre_mod, ONLY: clon
 
     include "dimensions.h"
@@ -64,5 +64,5 @@
        end DO
 
-       if (ABS(xvrai(i) - xo1) > my_eps) then
+       if (ABS(xvrai(i) - xo1) > my_eps) THEN
           ! iter == 300
           print *, 'Pas de solution.'
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iso_verif_dyn.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iso_verif_dyn.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iso_verif_dyn.f90	(revision 5116)
@@ -5,23 +5,23 @@
 
     ! input:
-    real :: x
-    character(len=*) :: err_msg ! message d''erreur à afficher
+    REAL :: x
+    CHARACTER(LEN=*) :: err_msg ! message d''erreur à afficher
 
     ! output
-    real :: borne
+    REAL :: borne
     parameter (borne=1e19)
-    integer :: iso_verif_noNaN_nostop
+    INTEGER :: iso_verif_noNaN_nostop
 
-    if ((x>-borne).and.(x<borne)) then
+    if ((x>-borne).and.(x<borne)) THEN
             iso_verif_noNAN_nostop=0
     else
-        write(*,*) 'erreur detectee par iso_verif_nonNaN:'
-        write(*,*) err_msg
-        write(*,*) 'x=',x
+        WRITE(*,*) 'erreur detectee par iso_verif_nonNaN:'
+        WRITE(*,*) err_msg
+        WRITE(*,*) 'x=',x
         iso_verif_noNaN_nostop=1
     endif
 
-    return
-end function iso_verif_nonan_nostop
+    RETURN
+END FUNCTION iso_verif_nonan_nostop
 
   function iso_verif_egalite_nostop &
@@ -33,31 +33,31 @@
 
     ! input:
-    real :: a, b
-    character(len=*) :: err_msg ! message d''erreur à afficher
+    REAL :: a, b
+    CHARACTER(LEN=*) :: err_msg ! message d''erreur à afficher
 
     ! locals
-    real :: errmax ! erreur maximale en absolu.
-    real :: errmaxrel ! erreur maximale en relatif autorisée
+    REAL :: errmax ! erreur maximale en absolu.
+    REAL :: errmaxrel ! erreur maximale en relatif autorisée
     parameter (errmax=1e-8)
     parameter (errmaxrel=1e-3)
 
     ! output
-    integer :: iso_verif_egalite_nostop
+    INTEGER :: iso_verif_egalite_nostop
 
     iso_verif_egalite_nostop=0
 
-    if (abs(a-b)>errmax) then
+    if (abs(a-b)>errmax) THEN
       if (abs((a-b)/max(max(abs(b),abs(a)),1e-18)) &
-            >errmaxrel) then
-        write(*,*) 'erreur detectee par iso_verif_egalite:'
-        write(*,*) err_msg
-        write(*,*) 'a=',a
-        write(*,*) 'b=',b
+            >errmaxrel) THEN
+        WRITE(*,*) 'erreur detectee par iso_verif_egalite:'
+        WRITE(*,*) err_msg
+        WRITE(*,*) 'a=',a
+        WRITE(*,*) 'b=',b
         iso_verif_egalite_nostop=1
       endif
     endif
 
-    return
-end function iso_verif_egalite_nostop
+    RETURN
+END FUNCTION iso_verif_egalite_nostop
 
 
@@ -68,21 +68,21 @@
 
     ! input:
-    real :: x,q
-    integer :: iso ! 2=HDO, 1=O18
-    character(len=*) :: err_msg ! message d''erreur à afficher
+    REAL :: x,q
+    INTEGER :: iso ! 2=HDO, 1=O18
+    CHARACTER(LEN=*) :: err_msg ! message d''erreur à afficher
 
     ! locals
-    real :: qmin,deltaD
-    real :: deltaDmax,deltaDmin,tnat
+    REAL :: qmin,deltaD
+    REAL :: deltaDmax,deltaDmin,tnat
     parameter (qmin=1e-11)
     parameter (deltaDmax=200.0,deltaDmin=-999.9)
 
     ! output
-    integer :: iso_verif_aberrant_nostop
+    INTEGER :: iso_verif_aberrant_nostop
 
     iso_verif_aberrant_nostop=0
 
     ! verifier que HDO est raisonable
-     if (q>qmin) then
+     if (q>qmin) THEN
          IF(getKey('tnat', tnat, isoName(iso))) THEN
               err_msg = 'Missing isotopic parameter "tnat"'
@@ -91,16 +91,14 @@
          END IF
          deltaD=(x/q/tnat-1)*1000
-         if ((deltaD>deltaDmax).or.(deltaD<deltaDmin)) then
-              write(*,*) 'erreur detectee par iso_verif_aberrant:'
-              write(*,*) err_msg
-              write(*,*) 'q=',q
-              write(*,*) 'deltaD=',deltaD
-              write(*,*) 'iso=',iso
+         if ((deltaD>deltaDmax).or.(deltaD<deltaDmin)) THEN
+              WRITE(*,*) 'erreur detectee par iso_verif_aberrant:'
+              WRITE(*,*) err_msg
+              WRITE(*,*) 'q=',q
+              WRITE(*,*) 'deltaD=',deltaD
+              WRITE(*,*) 'iso=',iso
               iso_verif_aberrant_nostop=1
-         endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) then
-      endif !if (q(i,k,iq).gt.qmin) then
+         endif !if ((deltaD.gt.deltaDmax).or.(deltaD.lt.deltaDmin)) THEN
+      endif !if (q(i,k,iq).gt.qmin) THEN
+    RETURN
+END FUNCTION iso_verif_aberrant_nostop
 
-
-    return
-end function iso_verif_aberrant_nostop
-
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limx.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limx.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limx.f90	(revision 5116)
@@ -22,7 +22,7 @@
   !   Arguments:
   !   ----------
-  real :: pente_max
+  REAL :: pente_max
   REAL :: s0(ip1jmp1,llm),sm(ip1jmp1,llm)
-  real :: sx(ip1jmp1,llm)
+  REAL :: sx(ip1jmp1,llm)
   !
   !  Local
@@ -30,13 +30,13 @@
   !
   INTEGER :: ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
-  integer :: n0,iadvplus(ip1jmp1,llm),nl(llm)
+  INTEGER :: n0,iadvplus(ip1jmp1,llm),nl(llm)
   !
   REAL :: q(ip1jmp1,llm)
-  real :: dxq(ip1jmp1,llm)
+  REAL :: dxq(ip1jmp1,llm)
 
 
   REAL :: new_m,zm
-  real :: dxqu(ip1jmp1)
-  real :: adxqu(ip1jmp1),dxqmax(ip1jmp1)
+  REAL :: dxqu(ip1jmp1)
+  REAL :: adxqu(ip1jmp1),dxqmax(ip1jmp1)
 
   Logical :: extremum,first
@@ -44,5 +44,5 @@
 
   REAL :: SSUM
-  integer :: ismax,ismin
+  INTEGER :: ismax,ismin
   EXTERNAL  SSUM, ismin,ismax
 
@@ -84,6 +84,6 @@
 
      do ij=iip2+1,ip1jm
-        if(     dxqu(ij-1)*dxqu(ij)>0. &
-              .and. dxq(ij,l)*dxqu(ij)>0.) then
+        IF(     dxqu(ij-1)*dxqu(ij)>0. &
+              .and. dxq(ij,l)*dxqu(ij)>0.) THEN
           dxq(ij,l)= &
                 sign(min(abs(dxq(ij,l)),dxqmax(ij)),dxq(ij,l))
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limy.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limy.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limy.f90	(revision 5116)
@@ -15,4 +15,5 @@
   !   --------------------------------------------------------------------
   USE comconst_mod, ONLY: pi
+  USE lmdz_libmath, ONLY: ismax, ismin
   IMPLICIT NONE
   !
@@ -24,6 +25,6 @@
   !   Arguments:
   !   ----------
-  real :: pente_max
-  real :: s0(ip1jmp1,llm),sy(ip1jmp1,llm),sm(ip1jmp1,llm)
+  REAL :: pente_max
+  REAL :: s0(ip1jmp1,llm),sy(ip1jmp1,llm),sm(ip1jmp1,llm)
   !
   !  Local
@@ -34,6 +35,6 @@
   REAL :: q(ip1jmp1,llm)
   REAL :: airej2,airejjm,airescb(iim),airesch(iim)
-  real :: sigv,dyq(ip1jmp1),dyqv(ip1jm)
-  real :: adyqv(ip1jm),dyqmax(ip1jmp1)
+  REAL :: sigv,dyq(ip1jmp1),dyqv(ip1jm)
+  REAL :: adyqv(ip1jm),dyqmax(ip1jmp1)
   REAL :: qbyv(ip1jm,llm)
 
@@ -42,18 +43,16 @@
   save first
 
-  real :: convpn,convps,convmpn,convmps
-  real :: sinlon(iip1),sinlondlon(iip1)
-  real :: coslon(iip1),coslondlon(iip1)
+  REAL :: convpn,convps,convmpn,convmps
+  REAL :: sinlon(iip1),sinlondlon(iip1)
+  REAL :: coslon(iip1),coslondlon(iip1)
   save sinlon,coslon,sinlondlon,coslondlon
   !
   !
   REAL :: SSUM
-  integer :: ismax,ismin
-  EXTERNAL  SSUM, convflu,ismin,ismax
-  EXTERNAL filtreg
+  EXTERNAL  SSUM
 
   data first/.TRUE./
 
-  if(first) then
+  IF(first) THEN
      PRINT*,'SCHEMA AMONT NOUVEAU'
      first=.FALSE.
@@ -129,7 +128,7 @@
   !   cas ou on a un extremum au pole
 
-  ! if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
+  ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
   !    &   appn=0.
-  ! if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
+  ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
   !    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
   !    &   apps=0.
@@ -150,6 +149,6 @@
   !  enddo
 
-  if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1))<=0.) &
-        then
+  IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1))<=0.) &
+        THEN
      do ij=1,iip1
         dyqmax(ij)=0.
@@ -161,7 +160,7 @@
   endif
 
-  if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* &
+  IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)* &
         dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)<=0.) &
-        then
+        THEN
      do ij=ip1jm+1,ip1jmp1
         dyqmax(ij)=0.
@@ -176,5 +175,5 @@
 
   do ij=1,ip1jmp1 ! cf below: should it be ip1jm instead ?
-     if(dyqv(ij)*dyqv(ij-iip1)>0.) then  ! /!\ causes Warning: iteration 1056 invokes undefined behavior [-Waggressive-loop-optimizations] in 32x32x39
+     IF(dyqv(ij)*dyqv(ij-iip1)>0.) then  ! /!\ causes Warning: iteration 1056 invokes undefined behavior [-Waggressive-loop-optimizations] in 32x32x39
         dyq(ij)=sign(min(abs(dyq(ij)),dyqmax(ij)),dyq(ij))
      else
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limz.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limz.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limz.f90	(revision 5116)
@@ -22,7 +22,7 @@
   !   Arguments:
   !   ----------
-  real :: pente_max
+  REAL :: pente_max
   REAL :: s0(ip1jmp1,llm),sm(ip1jmp1,llm)
-  real :: sz(ip1jmp1,llm)
+  REAL :: sz(ip1jmp1,llm)
   !
   !  Local
@@ -30,13 +30,13 @@
   !
   INTEGER :: ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
-  integer :: n0,iadvplus(ip1jmp1,llm),nl(llm)
+  INTEGER :: n0,iadvplus(ip1jmp1,llm),nl(llm)
   !
   REAL :: q(ip1jmp1,llm)
-  real :: dzq(ip1jmp1,llm)
+  REAL :: dzq(ip1jmp1,llm)
 
 
   REAL :: new_m,zm
-  real :: dzqw(ip1jmp1)
-  real :: adzqw(ip1jmp1),dzqmax(ip1jmp1)
+  REAL :: dzqw(ip1jmp1)
+  REAL :: adzqw(ip1jmp1),dzqmax(ip1jmp1)
 
   Logical :: extremum,first
@@ -44,5 +44,5 @@
 
   REAL :: SSUM
-  integer :: ismax,ismin
+  INTEGER :: ismax,ismin
   EXTERNAL  SSUM, ismin,ismax
 
@@ -77,6 +77,6 @@
 
      do l=2,llm-1
-        if(     dzqw(l-1)*dzqw(l)>0. &
-              .and. dzq(ij,l)*dzqw(l)>0.) then
+        IF(     dzqw(l-1)*dzqw(l)>0. &
+              .and. dzq(ij,l)*dzqw(l)>0.) THEN
           dzq(ij,l)= &
                 sign(min(abs(dzq(ij,l)),dzqmax(l)),dzq(ij,l))
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/pentes_ini.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/pentes_ini.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/pentes_ini.f90	(revision 5116)
@@ -31,5 +31,5 @@
   !   Arguments:
   !   ----------
-  integer :: mode
+  INTEGER :: mode
   REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
   REAL :: q( iip1,jjp1,llm,0:3)
@@ -42,19 +42,19 @@
   REAL :: s0( iip1,jjp1,llm ),  sx( iip1,jjp1,llm )
   REAL :: sy( iip1,jjp1,llm ),  sz( iip1,jjp1,llm )
-  real :: masn,mass,zz
+  REAL :: masn,mass,zz
   INTEGER :: i,j,l,iq
 
   !  modif Fred 24 03 96
 
-  real :: sinlon(iip1),sinlondlon(iip1)
-  real :: coslon(iip1),coslondlon(iip1)
+  REAL :: sinlon(iip1),sinlondlon(iip1)
+  REAL :: coslon(iip1),coslondlon(iip1)
   save sinlon,coslon,sinlondlon,coslondlon
-  real :: dyn1,dyn2,dys1,dys2
-  real :: qpn,qps,dqzpn,dqzps
-  real :: smn,sms,s0n,s0s,sxn(iip1),sxs(iip1)
-  real :: qmin,zq,pente_max
+  REAL :: dyn1,dyn2,dys1,dys2
+  REAL :: qpn,qps,dqzpn,dqzps
+  REAL :: smn,sms,s0n,s0s,sxn(iip1),sxs(iip1)
+  REAL :: qmin,zq,pente_max
   !
   REAL :: SSUM
-  integer :: ismax,ismin,lati,latf
+  INTEGER :: ismax,ismin,lati,latf
   EXTERNAL  SSUM, ismin,ismax
   logical :: first
@@ -72,7 +72,7 @@
   limit = .TRUE.
   pente_max=2
-  ! if (mode.eq.1.or.mode.eq.3) then
-  ! if (mode.eq.1) then
-  if (mode>=1) then
+  ! if (mode.eq.1.or.mode.eq.3) THEN
+  ! if (mode.eq.1) THEN
+  if (mode>=1) THEN
     lati=2
     latf=jjm
@@ -84,5 +84,5 @@
   qmin=0.4995
   qmin=0.
-  if(first) then
+  IF(first) THEN
      PRINT*,'SCHEMA AMONT NOUVEAU'
      first=.FALSE.
@@ -181,5 +181,5 @@
    !     do i=1,iip1
    !        zq=s0(i,j,l)/sm(i,j,l)
-   !       if(zq.lt.qmin)
+   !       IF(zq.lt.qmin)
   !    ,       PRINT*,'avant advx1, s0(',i,',',j,',',l,')=',zq
    !     enddo
@@ -187,5 +187,5 @@
    ! enddo
   !CC
-   if(mode==2) then
+   IF(mode==2) THEN
       do l=1,llm
         s0s=0.
@@ -247,5 +247,5 @@
   endif
 
-  if (mode==4) then
+  if (mode==4) THEN
      do l=1,llm
         do i=1,iip1
@@ -261,5 +261,5 @@
    CALL advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
   ! CALL minmaxq(zq,1.e33,-1.e33,'avant advy     ')
-  if (mode==4) then
+  if (mode==4) THEN
      do l=1,llm
         do i=1,iip1
@@ -282,5 +282,5 @@
    CALL limz(s0,sz,sm,pente_max)
    CALL advz( limit,dtvr,w,sm,s0,sx,sy,sz )
-  if (mode==4) then
+  if (mode==4) THEN
      do l=1,llm
         do i=1,iip1
@@ -306,5 +306,5 @@
 
   ! CALL minmaxq(zq,1.e33,-1.e33,'avant advx     ')
-  if (mode==4) then
+  if (mode==4) THEN
      do l=1,llm
         do i=1,iip1
@@ -323,5 +323,5 @@
   !      do i=1,iip1
   !         zq=s0(i,j,l)/sm(i,j,l)
-  !        if(zq.lt.qmin)
+  !        IF(zq.lt.qmin)
   !    ,       PRINT*,'apres advx2, s0(',i,',',j,',',l,')=',zq
   !      enddo
@@ -346,5 +346,5 @@
   ! Traitements specifiques au pole
 
-  if(mode>=1) then
+  IF(mode>=1) THEN
   DO l=1,llm
   !   filtrages aux poles
@@ -361,5 +361,5 @@
         q( i,jjp1,llm+1-l,0)=qps
      enddo
-     if(mode==3) then
+     IF(mode==3) THEN
         dyn1=0.
         dys1=0.
@@ -382,5 +382,5 @@
         enddo
      endif
-     if(mode==1) then
+     IF(mode==1) THEN
   !   on filtre les valeurs au bord de la "grande maille pole"
         dyn1=0.
@@ -459,5 +459,5 @@
       do j=1,jjp1
        do i=1,iip1
-         if(q(i,j,l,0)<qmin) &
+         IF(q(i,j,l,0)<qmin) &
                PRINT*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0)
        enddo
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/ppm3d.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/ppm3d.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/ppm3d.f90	(revision 5116)
@@ -276,15 +276,15 @@
   !
 
-  real :: Q(IMR,JNP,NLAY,NC),PS1(IMR,JNP),PS2(IMR,JNP), &
+  REAL :: Q(IMR,JNP,NLAY,NC),PS1(IMR,JNP),PS2(IMR,JNP), &
         U(IMR,JNP,NLAY),V(IMR,JNP,NLAY),AP(NLAY+1), &
         BP(NLAY+1),W(IMR,JNP,NLAY),NDT,val(NLAY),Umax
-  integer :: IGD,IORD,JORD,KORD,NC,IMR,JNP,j1,NLAY,AE
-  integer :: IMRD2
-  real :: PT
+  INTEGER :: IGD,IORD,JORD,KORD,NC,IMR,JNP,j1,NLAY,AE
+  INTEGER :: IMRD2
+  REAL :: PT
   logical :: cross, fill, dum
   !
   ! Local dynamic arrays
   !
-  real :: CRX(IMR,JNP),CRY(IMR,JNP),xmass(IMR,JNP),ymass(IMR,JNP), &
+  REAL :: CRX(IMR,JNP),CRY(IMR,JNP),xmass(IMR,JNP),ymass(IMR,JNP), &
         fx1(IMR+1),DPI(IMR,JNP,NLAY),delp1(IMR,JNP,NLAY), &
         WK1(IMR,JNP,NLAY),PU(IMR,JNP),PV(IMR,JNP),DC2(IMR,JNP), &
@@ -294,5 +294,5 @@
   ! Local static  arrays
   !
-  real :: DTDX(Jmax), DTDX5(Jmax), acosp(Jmax), &
+  REAL :: DTDX(Jmax), DTDX5(Jmax), acosp(Jmax), &
         cosp(Jmax), cose(Jmax), DAP(kmax),DBK(Kmax)
   data NDT0, NSTEP /0, 0/
@@ -314,9 +314,9 @@
   !
   ! *********** Initialization **********************
-  if(NSTEP==1) then
-  !
-  write(6,*) '------------------------------------ '
-  write(6,*) 'NASA/GSFC Transport Core Version 4.5'
-  write(6,*) '------------------------------------ '
+  IF(NSTEP==1) THEN
+  !
+  WRITE(6,*) '------------------------------------ '
+  WRITE(6,*) 'NASA/GSFC Transport Core Version 4.5'
+  WRITE(6,*) '------------------------------------ '
   !
   WRITE(6,*) 'IMR=',IMR,' JNP=',JNP,' NLAY=',NLAY,' j1=',j1
@@ -324,21 +324,21 @@
   !
   ! controles sur les parametres
-  if(NLAY<6) then
-    write(6,*) 'NLAY must be >= 6'
+  IF(NLAY<6) THEN
+    WRITE(6,*) 'NLAY must be >= 6'
     stop
   endif
-  if (JNP<NLAY) then
-     write(6,*) 'JNP must be >= NLAY'
+  if (JNP<NLAY) THEN
+     WRITE(6,*) 'JNP must be >= NLAY'
     stop
   endif
   IMRD2=mod(IMR,2)
-  if (j1==2.and.IMRD2/=0) then
-     write(6,*) 'if j1=2 IMR must be an even integer'
+  if (j1==2.and.IMRD2/=0) THEN
+     WRITE(6,*) 'if j1=2 IMR must be an even integer'
     stop
   endif
 
   !
-  if(Jmax<JNP .or. Kmax<NLAY) then
-    write(6,*) 'Jmax or Kmax is too small'
+  IF(Jmax<JNP .or. Kmax<NLAY) THEN
+    WRITE(6,*) 'Jmax or Kmax is too small'
     stop
   endif
@@ -353,5 +353,5 @@
   DP =    PI / REAL(JMR)
   !
-  if(IGD==0) then
+  IF(IGD==0) THEN
   ! Compute analytic cosine at cell edges
         CALL cosa(cosp,cose,JNP,PI,DP)
@@ -372,19 +372,19 @@
   endif
   !
-  if(NDT0 /= NDT) then
+  IF(NDT0 /= NDT) THEN
   DT   = NDT
   NDT0 = NDT
 
-    if(Umax < 180.) then
-     write(6,*) 'Umax may be too small!'
+    IF(Umax < 180.) THEN
+     WRITE(6,*) 'Umax may be too small!'
     endif
   CR1  = abs(Umax*DT)/(DL*AE)
   MaxDT = DP*AE / abs(Umax) + 0.5
-  write(6,*)'Largest time step for max(V)=',Umax,' is ',MaxDT
-  if(MaxDT < abs(NDT)) then
-        write(6,*) 'Warning!!! NDT maybe too large!'
-  endif
-  !
-  if(CR1>=0.95) then
+  WRITE(6,*)'Largest time step for max(V)=',Umax,' is ',MaxDT
+  IF(MaxDT < abs(NDT)) THEN
+        WRITE(6,*) 'Warning!!! NDT maybe too large!'
+  endif
+  !
+  IF(CR1>=0.95) THEN
   JS0 = 0
   JN0 = 0
@@ -413,5 +413,5 @@
   DTDY5 = 0.5*DTDY
   !
-  !  write(6,*) 'J1=',J1,' J2=', J2
+  !  WRITE(6,*) 'J1=',J1,' J2=', J2
   endif
   !
@@ -429,5 +429,5 @@
 
   !
-  if(j1/=2) then
+  IF(j1/=2) THEN
   DO IC=1,NC
   DO L=1,NLAY
@@ -453,5 +453,5 @@
   do k=1,NLAY
   !
-  if(IGD==0) then
+  IF(IGD==0) THEN
   ! Convert winds on A-Grid to Courant # on C-Grid.
   CALL A2C(U(1,1,k),V(1,1,k),IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)
@@ -480,5 +480,5 @@
   do j=JS0,j1+1,-1
   do i=1,IMR
-  if(abs(CRX(i,j))>1.) then
+  IF(abs(CRX(i,j))>1.) THEN
         JS = j
         go to 2222
@@ -490,5 +490,5 @@
   do j=JN0,j2-1
   do i=1,IMR
-  if(abs(CRX(i,j))>1.) then
+  IF(abs(CRX(i,j))>1.) THEN
         JN = j
         go to 2233
@@ -498,5 +498,5 @@
 2233   continue
   !
-  if(j1/=2) then           ! Enlarged polar cap.
+  IF(j1/=2) then           ! Enlarged polar cap.
   do i=1,IMR
   DPI(i,  2,k) = 0.
@@ -586,5 +586,5 @@
   enddo
   !
-  if(j1==2) then
+  IF(j1==2) THEN
     IMH = IMR/2
   do i=1,IMH
@@ -608,5 +608,5 @@
   ! E-W advective cross term
   do j=J1,J2
-  if(J>JS  .and. J<JN) GO TO 250
+  IF(J>JS  .and. J<JN) GO TO 250
   !
   do i=1,IMR
@@ -623,5 +623,5 @@
   ru = UA(i,j) - iu
   iiu = i-iu
-  if(UA(i,j)>=0.) then
+  IF(UA(i,j)>=0.) THEN
   wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))
   else
@@ -633,5 +633,5 @@
   END DO
   !
-  if(JN/=0) then
+  IF(JN/=0) THEN
   do j=JS+1,JN-1
   !
@@ -661,7 +661,7 @@
   enddo
   !
-    if(cross) then
+    IF(cross) THEN
   ! Add cross terms in the vertical direction.
-    if(IORD >= 2) then
+    IF(IORD >= 2) THEN
             iad = 2
     else
@@ -669,5 +669,5 @@
     endif
   !
-    if(JORD >= 2) then
+    IF(JORD >= 2) THEN
             jad = 2
     else
@@ -750,5 +750,5 @@
   !
 
-  if(fill) CALL qckxyz(DQ(1,1,1,IC),DC2,IMR,JNP,NLAY,j1,j2, &
+  IF(fill) CALL qckxyz(DQ(1,1,1,IC),DC2,IMR,JNP,NLAY,j1,j2, &
         cosp,acosp,.FALSE.,IC,NSTEP)
   !
@@ -765,5 +765,5 @@
   enddo
   !
-  if(j1/=2) then
+  IF(j1/=2) THEN
   DO k=1,NLAY
   DO I=1,IMR
@@ -776,5 +776,5 @@
   END DO
   !
-  if(j1/=2) then
+  IF(j1/=2) THEN
   DO k=1,NLAY
   DO i=1,IMR
@@ -794,13 +794,13 @@
   integer,parameter :: kmax = 150
   real,parameter :: R23 = 2./3., R3 = 1./3.
-  integer :: IMR,JNP,NLAY,J1,KORD
-  real :: WZ(IMR,JNP,NLAY),P(IMR,JNP,NLAY),DC(IMR,JNP,NLAY), &
+  INTEGER :: IMR,JNP,NLAY,J1,KORD
+  REAL :: WZ(IMR,JNP,NLAY),P(IMR,JNP,NLAY),DC(IMR,JNP,NLAY), &
         wk1(IMR,*),delp(IMR,JNP,NLAY),DQ(IMR,JNP,NLAY), &
         DQDT(IMR,JNP,NLAY)
   ! Assuming JNP >= NLAY
-  real :: AR(IMR,*),AL(IMR,*),A6(IMR,*),flux(IMR,*),wk2(IMR,*), &
+  REAL :: AR(IMR,*),AL(IMR,*),A6(IMR,*),flux(IMR,*),wk2(IMR,*), &
         wz2(IMR,*)
-  integer :: JMR,IMJM,NLAYM1,LMT,K,I,J
-  real :: c0,c1,c2,tmp,qmax,qmin,a,b,fct,a1,a2,cm,cp
+  INTEGER :: JMR,IMJM,NLAYM1,LMT,K,I,J
+  REAL :: c0,c1,c2,tmp,qmax,qmin,a,b,fct,a1,a2,cm,cp
   !
   JMR = JNP - 1
@@ -871,5 +871,5 @@
   !
   ! Check if change sign
-  if(wk1(i,1)*AL(i,1)<=0.) then
+  IF(wk1(i,1)*AL(i,1)<=0.) THEN
              AL(i,1) = 0.
          flux(i,1) = 0.
@@ -887,5 +887,5 @@
   AR(i,NLAY) = wk1(i,NLAY) + fct
   AL(i,NLAY) = wk1(i,NLAY) - (fct+fct)
-  if(wk1(i,NLAY)*AR(i,NLAY)<=0.) AR(i,NLAY) = 0.
+  IF(wk1(i,NLAY)*AR(i,NLAY)<=0.) AR(i,NLAY) = 0.
   flux(i,NLAY) = AR(i,NLAY) -  wk1(i,NLAY)
   END DO
@@ -926,5 +926,5 @@
   !
   ! Interior depending on KORD
-  if(LMT<=2) &
+  IF(LMT<=2) &
         CALL lmtppm(flux(1,2),A6(1,2),AR(1,2),AL(1,2),wk1(1,2), &
         IMR*(NLAY-2),LMT)
@@ -933,5 +933,5 @@
   !
   DO i=1,IMR*NLAYM1
-  IF(wz2(i,1)>0.) then
+  IF(wz2(i,1)>0.) THEN
     CM = wz2(i,1) / wk2(i,1)
     flux(i,2) = AR(i,1)+0.5*CM*(AL(i,1)-AR(i,1)+A6(i,1)*(1.-R23*CM))
@@ -961,19 +961,19 @@
 2000   continue
   END DO
-  return
-end subroutine fzppm
+  RETURN
+END SUBROUTINE fzppm
 !
 SUBROUTINE xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ,Q,UC, &
         fx1,xmass,IORD)
   IMPLICIT NONE
-  integer :: IMR,JNP,IML,j1,j2,JN,JS,IORD
-  real :: PU,DQ,Q,UC,fx1,xmass
-  real :: dc,qtmp
-  integer :: ISAVE(IMR)
+  INTEGER :: IMR,JNP,IML,j1,j2,JN,JS,IORD
+  REAL :: PU,DQ,Q,UC,fx1,xmass
+  REAL :: dc,qtmp
+  INTEGER :: ISAVE(IMR)
   dimension UC(IMR,*),DC(-IML:IMR+IML+1),xmass(IMR,JNP) &
         ,fx1(IMR+1),DQ(IMR,JNP),qtmp(-IML:IMR+1+IML)
   dimension PU(IMR,JNP),Q(IMR,JNP)
-  integer :: jvan,j1vl,j2vl,j,i,iu,itmp,ist,imp
-  real :: rut
+  INTEGER :: jvan,j1vl,j2vl,j,i,iu,itmp,ist,imp
+  REAL :: rut
   !
   IMP = IMR + 1
@@ -990,5 +990,5 @@
   enddo
   !
-  if(j>=JN .or. j<=JS) goto 2222
+  IF(j>=JN .or. j<=JS) goto 2222
   ! ************* Eulerian **********
   !
@@ -1007,5 +1007,5 @@
   DC(0) = DC(IMR)
   !
-  if(IORD==2 .or. j<=j1vl .or. j>=j2vl) then
+  IF(IORD==2 .or. j<=j1vl .or. j>=j2vl) THEN
   DO i=1,IMR
   iu = REAL(i) - uc(i,j)
@@ -1058,10 +1058,10 @@
   !
   do i=1,IMR
-  IF(uc(i,j)>1.) then
+  IF(uc(i,j)>1.) THEN
   !DIR$ NOVECTOR
     do ist = ISAVE(i),i-1
     fx1(i) = fx1(i) + qtmp(ist)
     enddo
-  elseIF(uc(i,j)<-1.) then
+  elseIF(uc(i,j)<-1.) THEN
     do ist = i,ISAVE(i)-1
     fx1(i) = fx1(i) - qtmp(ist)
@@ -1084,26 +1084,26 @@
   !
   END DO
-  return
-end subroutine xtp
+  RETURN
+END SUBROUTINE xtp
 !
 SUBROUTINE fxppm(IMR,IML,UT,P,DC,flux,IORD)
   IMPLICIT NONE
-  integer :: IMR,IML,IORD
-  real :: UT,P,DC,flux
+  INTEGER :: IMR,IML,IORD
+  REAL :: UT,P,DC,flux
   real,parameter ::  R3 = 1./3., R23 = 2./3.
   DIMENSION UT(*),flux(*),P(-IML:IMR+IML+1),DC(-IML:IMR+IML+1)
   REAL :: AR(0:IMR),AL(0:IMR),A6(0:IMR)
-  integer :: LMT,IMP,JLVL,i
+  INTEGER :: LMT,IMP,JLVL,i
    ! logical first
    ! data first /.TRUE./
    ! SAVE LMT
-   ! if(first) then
+   ! IF(first) THEN
   !
   ! correction calcul de LMT a chaque passage pour pouvoir choisir
   ! plusieurs schemas PPM pour differents traceurs
-  !  IF (IORD.LE.0) then
-  !        if(IMR.GE.144) then
+  !  IF (IORD.LE.0) THEN
+  !        IF(IMR.GE.144) THEN
   !              LMT = 0
-  !        elseif(IMR.GE.72) then
+  !        elseif(IMR.GE.72) THEN
   !              LMT = 1
   !        else
@@ -1115,5 +1115,5 @@
   !
   LMT = IORD - 3
-   ! write(6,*) 'PPM option in E-W direction = ', LMT
+   ! WRITE(6,*) 'PPM option in E-W direction = ', LMT
    ! first = .FALSE.
    ! endif
@@ -1132,5 +1132,5 @@
   END DO
   !
-  if(LMT<=2) CALL lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT)
+  IF(LMT<=2) CALL lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT)
   !
   AL(0) = AL(IMR)
@@ -1139,5 +1139,5 @@
   !
   DO i=1,IMR
-  IF(UT(i)>0.) then
+  IF(UT(i)>0.) THEN
   flux(i) = AR(i-1) + 0.5*UT(i)*(AL(i-1) - AR(i-1) + &
         A6(i-1)*(1.-R23*UT(i)) )
@@ -1147,14 +1147,14 @@
   endif
   enddo
-  return
-end subroutine fxppm
+  RETURN
+END SUBROUTINE fxppm
 !
 SUBROUTINE xmist(IMR,IML,P,DC)
   IMPLICIT NONE
-  integer :: IMR,IML
+  INTEGER :: IMR,IML
   real,parameter :: R24 = 1./24.
-  real :: P(-IML:IMR+1+IML),DC(-IML:IMR+1+IML)
-  integer :: i
-  real :: tmp,pmax,pmin
+  REAL :: P(-IML:IMR+1+IML),DC(-IML:IMR+1+IML)
+  INTEGER :: i
+  REAL :: tmp,pmax,pmin
   !
   do i=1,IMR
@@ -1164,23 +1164,23 @@
   DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp)
   END DO
-  return
-end subroutine xmist
+  RETURN
+END SUBROUTINE xmist
 !
 SUBROUTINE ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ,P,VC,DC2 &
         ,ymass,fx,A6,AR,AL,JORD)
   IMPLICIT NONE
-  integer :: IMR,JNP,j1,j2,JORD
-  real :: acosp,RCAP,DQ,P,VC,DC2,ymass,fx,A6,AR,AL
+  INTEGER :: IMR,JNP,j1,j2,JORD
+  REAL :: acosp,RCAP,DQ,P,VC,DC2,ymass,fx,A6,AR,AL
   dimension P(IMR,JNP),VC(IMR,JNP),ymass(IMR,JNP) &
         ,DC2(IMR,JNP),DQ(IMR,JNP),acosp(JNP)
   ! Work array
   DIMENSION fx(IMR,JNP),AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)
-  integer :: JMR,len,i,jt,j
-  real :: sum1,sum2
+  INTEGER :: JMR,len,i,jt,j
+  REAL :: sum1,sum2
   !
   JMR = JNP - 1
   len = IMR*(J2-J1+2)
   !
-  if(JORD==1) then
+  IF(JORD==1) THEN
   DO i=1,len
   JT = REAL(J1) - VC(i,J1)
@@ -1191,6 +1191,5 @@
   CALL ymist(IMR,JNP,j1,P,DC2,4)
   !
-  if(JORD<=0 .or. JORD>=3) then
-
+  IF(JORD<=0 .or. JORD>=3) THEN
   CALL fyppm(VC,P,DC2,fx,IMR,JNP,j1,j2,A6,AR,AL,JORD)
 
@@ -1228,5 +1227,5 @@
   enddo
   !
-  if(j1/=2) then
+  IF(j1/=2) THEN
   do i=1,IMR
   DQ(i,  2) = sum1
@@ -1235,14 +1234,14 @@
   endif
   !
-  return
-end subroutine ytp
+  RETURN
+END SUBROUTINE ytp
 !
 subroutine  ymist(IMR,JNP,j1,P,DC,ID)
   IMPLICIT NONE
-  integer :: IMR,JNP,j1,ID
+  INTEGER :: IMR,JNP,j1,ID
   real,parameter :: R24 = 1./24.
-  real :: P(IMR,JNP),DC(IMR,JNP)
-  integer :: iimh,jmr,ijm3,imh,i
-  real :: pmax,pmin,tmp
+  REAL :: P(IMR,JNP),DC(IMR,JNP)
+  INTEGER :: iimh,jmr,ijm3,imh,i
+  REAL :: pmax,pmin,tmp
   !
   IMH = IMR / 2
@@ -1291,5 +1290,5 @@
   ENDIF
   !
-  if(j1/=2) then
+  IF(j1/=2) THEN
   do i=1,IMR
   DC(i,1) = 0.
@@ -1317,16 +1316,16 @@
   END DO
   endif
-  return
-end subroutine ymist
+  RETURN
+END SUBROUTINE ymist
 !
 SUBROUTINE fyppm(VC,P,DC,flux,IMR,JNP,j1,j2,A6,AR,AL,JORD)
   IMPLICIT NONE
-  integer :: IMR,JNP,j1,j2,JORD
+  INTEGER :: IMR,JNP,j1,j2,JORD
   real,parameter :: R3 = 1./3., R23 = 2./3.
-  real :: VC(IMR,*),flux(IMR,*),P(IMR,*),DC(IMR,*)
+  REAL :: VC(IMR,*),flux(IMR,*),P(IMR,*),DC(IMR,*)
   ! Local work arrays.
-  real :: AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)
-  integer :: LMT,i
-  integer :: IMH,JMR,j11,IMJM1,len
+  REAL :: AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)
+  INTEGER :: LMT,i
+  INTEGER :: IMH,JMR,j11,IMJM1,len
    ! logical first
    ! data first /.TRUE./
@@ -1338,9 +1337,9 @@
   IMJM1 = IMR*(J2-J1+2)
   len   = IMR*(J2-J1+3)
-   ! if(first) then
-   ! IF(JORD.LE.0) then
-   !       if(JMR.GE.90) then
+   ! IF(first) THEN
+   ! IF(JORD.LE.0) THEN
+   !       IF(JMR.GE.90) THEN
    !             LMT = 0
-   !       elseif(JMR.GE.45) then
+   !       elseif(JMR.GE.45) THEN
    !             LMT = 1
    !       else
@@ -1384,10 +1383,10 @@
   END DO
   !
-  if(LMT<=2) CALL lmtppm(DC(1,j11),A6(1,j11),AR(1,j11) &
+  IF(LMT<=2) CALL lmtppm(DC(1,j11),A6(1,j11),AR(1,j11) &
         ,AL(1,j11),P(1,j11),len,LMT)
   !
 
   DO i=1,IMJM1
-  IF(VC(i,j1)>0.) then
+  IF(VC(i,j1)>0.) THEN
   flux(i,j1) = AR(i,j11) + 0.5*VC(i,j1)*(AL(i,j11) - AR(i,j11) + &
         A6(i,j11)*(1.-R23*VC(i,j1)) )
@@ -1397,10 +1396,10 @@
   endif
   END DO
-  return
-end subroutine fyppm
+  RETURN
+END SUBROUTINE fyppm
 !
   SUBROUTINE yadv(IMR,JNP,j1,j2,p,VA,ady,wk,IAD)
     IMPLICIT NONE
-    integer :: IMR,JNP,j1,j2,IAD
+    INTEGER :: IMR,JNP,j1,j2,IAD
     REAL :: p(IMR,JNP),ady(IMR,JNP),VA(IMR,JNP)
     REAL :: WK(IMR,-1:JNP+2)
@@ -1426,24 +1425,24 @@
     wk(i+IMH,JNP+2) = p(i,JNP-2)
     enddo
-     ! write(*,*) 'toto 1'
+     ! WRITE(*,*) 'toto 1'
   ! --------------------------------
-  IF(IAD==2) then
+  IF(IAD==2) THEN
   do j=j1-1,j2+1
   do i=1,IMR
-   ! write(*,*) 'avt NINT','i=',i,'j=',j
+   ! WRITE(*,*) 'avt NINT','i=',i,'j=',j
   JP = NINT(VA(i,j))
   rv = JP - VA(i,j)
-   ! write(*,*) 'VA=',VA(i,j), 'JP1=',JP,'rv=',rv
+   ! WRITE(*,*) 'VA=',VA(i,j), 'JP1=',JP,'rv=',rv
   JP = j - JP
-   ! write(*,*) 'JP2=',JP
+   ! WRITE(*,*) 'JP2=',JP
   a1 = 0.5*(wk(i,jp+1)+wk(i,jp-1)) - wk(i,jp)
   b1 = 0.5*(wk(i,jp+1)-wk(i,jp-1))
-   ! write(*,*) 'a1=',a1,'b1=',b1
+   ! WRITE(*,*) 'a1=',a1,'b1=',b1
   ady(i,j) = wk(i,jp) + rv*(a1*rv + b1) - wk(i,j)
   enddo
   enddo
-   ! write(*,*) 'toto 2'
-  !
-  ELSEIF(IAD==1) then
+   ! WRITE(*,*) 'toto 2'
+  !
+  ELSEIF(IAD==1) THEN
     do j=j1-1,j2+1
   do i=1,imr
@@ -1454,5 +1453,5 @@
   ENDIF
   !
-    if(j1/=2) then
+    IF(j1/=2) THEN
     sum1 = 0.
     sum2 = 0.
@@ -1487,6 +1486,6 @@
     endif
   !
-    return
-end subroutine yadv
+    RETURN
+END SUBROUTINE yadv
 !
   SUBROUTINE xadv(IMR,JNP,j1,j2,p,UA,JS,JN,IML,adx,IAD)
@@ -1499,5 +1498,5 @@
     JMR = JNP-1
   do j=j1,j2
-  if(J>JS  .and. J<JN) GO TO 1309
+  IF(J>JS  .and. J<JN) GO TO 1309
   !
   do i=1,IMR
@@ -1519,10 +1518,10 @@
   adx(i,j) = qtmp(ip) + ru*(a1*ru + b1)
   enddo
-  ELSEIF(IAD==1) then
+  ELSEIF(IAD==1) THEN
   DO i=1,IMR
   iu = UA(i,j)
   ru = UA(i,j) - iu
   iiu = i-iu
-  if(UA(i,j)>=0.) then
+  IF(UA(i,j)>=0.) THEN
   adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))
   else
@@ -1560,5 +1559,5 @@
   adx(i,j) = qtmp(ip)- p(i,j) + ru*(a1*ru + b1)
   enddo
-  ELSEIF(IAD==1) then
+  ELSEIF(IAD==1) THEN
   ! 1st order
   DO i=1,IMR
@@ -1569,5 +1568,5 @@
   enddo
   !
-    if(j1/=2) then
+    IF(j1/=2) THEN
   do i=1,IMR
   adx(i,  2) = 0.
@@ -1580,6 +1579,6 @@
   adx(i,JNP) = 0.
   enddo
-    return
-end subroutine xadv
+    RETURN
+END SUBROUTINE xadv
 !
 SUBROUTINE lmtppm(DC,A6,AR,AL,P,IM,LMT)
@@ -1600,13 +1599,13 @@
   !
   real,parameter :: R12 = 1./12.
-  real :: A6(IM),AR(IM),AL(IM),P(IM),DC(IM)
-  integer :: IM,LMT
+  REAL :: A6(IM),AR(IM),AL(IM),P(IM),DC(IM)
+  INTEGER :: IM,LMT
   INTEGER :: i
   REAL :: da1,da2,a6da,fmin
   !
-  if(LMT==0) then
+  IF(LMT==0) THEN
   ! Full constraint
   do i=1,IM
-  if(DC(i)==0.) then
+  IF(DC(i)==0.) THEN
         AR(i) = p(i)
         AL(i) = p(i)
@@ -1616,8 +1615,8 @@
   da2  = da1**2
   A6DA = A6(i)*da1
-  if(A6DA < -da2) then
+  IF(A6DA < -da2) THEN
         A6(i) = 3.*(AL(i)-p(i))
         AR(i) = AL(i) - A6(i)
-  elseif(A6DA > da2) then
+  elseif(A6DA > da2) THEN
         A6(i) = 3.*(AR(i)-p(i))
         AL(i) = AR(i) - A6(i)
@@ -1625,13 +1624,13 @@
   endif
   END DO
-  elseif(LMT==1) then
+  elseif(LMT==1) THEN
   ! Semi-monotonic constraint
   do i=1,IM
-  if(abs(AR(i)-AL(i)) >= -A6(i)) go to 150
-  if(p(i)<AR(i) .and. p(i)<AL(i)) then
+  IF(abs(AR(i)-AL(i)) >= -A6(i)) go to 150
+  IF(p(i)<AR(i) .and. p(i)<AL(i)) THEN
         AR(i) = p(i)
         AL(i) = p(i)
         A6(i) = 0.
-  elseif(AR(i) > AL(i)) then
+  elseif(AR(i) > AL(i)) THEN
         A6(i) = 3.*(AL(i)-p(i))
         AR(i) = AL(i) - A6(i)
@@ -1642,14 +1641,14 @@
 150   continue
   END DO
-  elseif(LMT==2) then
+  elseif(LMT==2) THEN
   do i=1,IM
-  if(abs(AR(i)-AL(i)) >= -A6(i)) go to 250
+  IF(abs(AR(i)-AL(i)) >= -A6(i)) go to 250
   fmin = p(i) + 0.25*(AR(i)-AL(i))**2/A6(i) + A6(i)*R12
-  if(fmin>=0.) go to 250
-  if(p(i)<AR(i) .and. p(i)<AL(i)) then
+  IF(fmin>=0.) go to 250
+  IF(p(i)<AR(i) .and. p(i)<AL(i)) THEN
         AR(i) = p(i)
         AL(i) = p(i)
         A6(i) = 0.
-  elseif(AR(i) > AL(i)) then
+  elseif(AR(i) > AL(i)) THEN
         A6(i) = 3.*(AL(i)-p(i))
         AR(i) = AL(i) - A6(i)
@@ -1661,12 +1660,12 @@
   END DO
   endif
-  return
-end subroutine lmtppm
+  RETURN
+END SUBROUTINE lmtppm
 !
 SUBROUTINE A2C(U,V,IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)
   IMPLICIT NONE
-  integer :: IMR,JMR,j1,j2
-  real :: U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*),DTDY5
-  integer :: i,j
+  INTEGER :: IMR,JMR,j1,j2
+  REAL :: U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*),DTDY5
+  INTEGER :: i,j
   !
   do j=j1,j2
@@ -1683,13 +1682,13 @@
   CRY(i,2) = DTDY5*(V(i,2)+V(i,1))
   END DO
-  return
-end subroutine a2c
+  RETURN
+END SUBROUTINE a2c
 !
 SUBROUTINE cosa(cosp,cose,JNP,PI,DP)
   IMPLICIT NONE
-  integer :: JNP
-  real :: cosp(*),cose(*),PI,DP
-  integer :: JMR,j,jeq
-  real :: ph5
+  INTEGER :: JNP
+  REAL :: cosp(*),cose(*),PI,DP
+  INTEGER :: JMR,j,jeq
+  REAL :: ph5
   JMR = JNP-1
   do j=2,JNP
@@ -1699,5 +1698,5 @@
   !
   JEQ = (JNP+1) / 2
-  if(JMR == 2*(JMR/2) ) then
+  IF(JMR == 2*(JMR/2) ) THEN
   do j=JNP, JEQ+1, -1
    cose(j) =  cose(JNP+2-j)
@@ -1716,13 +1715,13 @@
   cosp(1) = 0.
   cosp(JNP) = 0.
-  return
-end subroutine cosa
+  RETURN
+END SUBROUTINE cosa
 !
 SUBROUTINE cosc(cosp,cose,JNP,PI,DP)
   IMPLICIT NONE
-  integer :: JNP
-  real :: cosp(*),cose(*),PI,DP
-  real :: phi
-  integer :: j
+  INTEGER :: JNP
+  REAL :: cosp(*),cose(*),PI,DP
+  REAL :: phi
+  INTEGER :: j
   !
   phi = -0.5*PI
@@ -1741,6 +1740,6 @@
    cosp(j) = 0.5*(cose(j)+cose(j+1))
   END DO
-  return
-end subroutine cosc
+  RETURN
+END SUBROUTINE cosc
 !
 SUBROUTINE qckxyz(Q,qtmp,IMR,JNP,NLAY,j1,j2,cosp,acosp, &
@@ -1752,5 +1751,5 @@
   logical :: cross
   INTEGER :: NLAYM1,len,ip,L,icr,ipy,ipx,i
-  real :: qup,qly,dup,sum
+  REAL :: qup,qly,dup,sum
   !
   NLAYM1 = NLAY-1
@@ -1762,12 +1761,12 @@
     icr = 1
   CALL filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
-  if(ipy==0) goto 50
+  IF(ipy==0) goto 50
   CALL filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
-  if(ipx==0) goto 50
-  !
-  if(cross) then
+  IF(ipx==0) goto 50
+  !
+  IF(cross) THEN
   CALL filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
   endif
-  if(icr==0) goto 50
+  IF(icr==0) goto 50
   !
   ! Vertical filling...
@@ -1785,11 +1784,11 @@
   !
   CALL filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
-  if(ipy==0) goto 225
+  IF(ipy==0) goto 225
   CALL filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
-  if(ipx==0) go to 225
-  if(cross) then
+  IF(ipx==0) go to 225
+  IF(cross) THEN
   CALL filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
   endif
-  if(icr==0) goto 225
+  IF(icr==0) goto 225
   !
   do i=1,len
@@ -1816,10 +1815,10 @@
   !
   CALL filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
-  if(ipy==0) goto 911
+  IF(ipy==0) goto 911
   CALL filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
-  if(ipx==0) goto 911
+  IF(ipx==0) goto 911
   !
   CALL filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
-  if(icr==0) goto 911
+  IF(icr==0) goto 911
   !
   DO  I=1,len
@@ -1841,11 +1840,11 @@
 911   continue
   !
-  if(ip>IMR) then
-  write(6,*) 'IC=',IC,' STEP=',NSTEP, &
+  IF(ip>IMR) THEN
+  WRITE(6,*) 'IC=',IC,' STEP=',NSTEP, &
         ' Vertical filling pts=',ip
   endif
   !
-  if(sum>1.e-25) then
-  write(6,*) IC,NSTEP,' Mass source from the ground=',sum
+  IF(sum>1.e-25) THEN
+  WRITE(6,*) IC,NSTEP,' Mass source from the ground=',sum
   endif
   RETURN
@@ -1854,8 +1853,8 @@
 SUBROUTINE filcr(q,IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
   IMPLICIT NONE
-  integer :: IMR,JNP,j1,j2,icr
-  real :: q(IMR,*),cosp(*),acosp(*),tiny
-  integer :: i,j
-  real :: dq,dn,d0,d1,ds,d2
+  INTEGER :: IMR,JNP,j1,j2,icr
+  REAL :: q(IMR,*),cosp(*),acosp(*),tiny
+  INTEGER :: i,j
+  REAL :: dq,dn,d0,d1,ds,d2
   icr = 0
   do j=j1+1,j2-1
@@ -1878,5 +1877,5 @@
   endif
   END DO
-  if(icr==0 .and. q(IMR,j)>=0.) goto 65
+  IF(icr==0 .and. q(IMR,j)>=0.) goto 65
   DO i=2,IMR
   IF(q(i,j)<0.) THEN
@@ -1940,5 +1939,5 @@
   !
   do i=1,IMR
-  if(q(i,j1)<0. .or. q(i,j2)<0.) then
+  IF(q(i,j1)<0. .or. q(i,j2)<0.) THEN
   icr = 1
   goto 80
@@ -1948,16 +1947,16 @@
 80   continue
   !
-  if(q(1,1)<0. .or. q(1,jnp)<0.) then
+  IF(q(1,1)<0. .or. q(1,jnp)<0.) THEN
   icr = 1
   endif
   !
-  return
-end subroutine filcr
+  RETURN
+END SUBROUTINE filcr
 !
 SUBROUTINE filns(q,IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
   IMPLICIT NONE
-  integer :: IMR,JNP,j1,j2,ipy
-  real :: q(IMR,*),cosp(*),acosp(*),tiny
-  real :: DP,CAP1,dq,dn,d0,d1,ds,d2
+  INTEGER :: IMR,JNP,j1,j2,ipy
+  REAL :: q(IMR,*),cosp(*),acosp(*),tiny
+  REAL :: DP,CAP1,dq,dn,d0,d1,ds,d2
   INTEGER :: i,j
    ! logical first
@@ -1965,5 +1964,5 @@
    ! save cap1
   !
-  !  if(first) then
+  !  IF(first) THEN
   DP = 4.*ATAN(1.)/REAL(JNP-1)
   CAP1 = IMR*(1.-COS((j1-1.5)*DP))/DP
@@ -2021,31 +2020,31 @@
   !
   ! Check Poles.
-  if(q(1,1)<0.) then
+  IF(q(1,1)<0.) THEN
   dq = q(1,1)*cap1/REAL(IMR)*acosp(j1)
   do i=1,imr
   q(i,1) = 0.
   q(i,j1) = q(i,j1) + dq
-  if(q(i,j1)<0.) ipy = 1
-  enddo
-  endif
-  !
-  if(q(1,JNP)<0.) then
+  IF(q(i,j1)<0.) ipy = 1
+  enddo
+  endif
+  !
+  IF(q(1,JNP)<0.) THEN
   dq = q(1,JNP)*cap1/REAL(IMR)*acosp(j2)
   do i=1,imr
   q(i,JNP) = 0.
   q(i,j2) = q(i,j2) + dq
-  if(q(i,j2)<0.) ipy = 1
-  enddo
-  endif
-  !
-  return
-end subroutine filns
+  IF(q(i,j2)<0.) ipy = 1
+  enddo
+  endif
+  !
+  RETURN
+END SUBROUTINE filns
 !
 SUBROUTINE filew(q,qtmp,IMR,JNP,j1,j2,ipx,tiny)
   IMPLICIT NONE
-  integer :: IMR,JNP,j1,j2,ipx
-  real :: q(IMR,*),qtmp(JNP,IMR),tiny
-  integer :: i,j
-  real :: d0,d1,d2
+  INTEGER :: IMR,JNP,j1,j2,ipx
+  REAL :: q(IMR,*),qtmp(JNP,IMR),tiny
+  INTEGER :: i,j
+  REAL :: d0,d1,d2
   !
   ipx = 0
@@ -2059,5 +2058,5 @@
   do i=2,imr-1
   do j=j1,j2
-  if(qtmp(j,i)<0.) then
+  IF(qtmp(j,i)<0.) THEN
   ipx =  1
   ! west
@@ -2077,5 +2076,5 @@
   i=1
   do j=j1,j2
-  if(qtmp(j,i)<0.) then
+  IF(qtmp(j,i)<0.) THEN
   ipx =  1
   ! west
@@ -2094,5 +2093,5 @@
   i=IMR
   do j=j1,j2
-  if(qtmp(j,i)<0.) then
+  IF(qtmp(j,i)<0.) THEN
   ipx =  1
   ! west
@@ -2110,5 +2109,5 @@
   END DO
   !
-  if(ipx/=0) then
+  IF(ipx/=0) THEN
   do j=j1,j2
   do i=1,imr
@@ -2119,17 +2118,17 @@
   !
   ! Poles.
-  if(q(1,1)<0 .or. q(1,JNP)<0.) ipx = 1
-  endif
-  return
-end subroutine filew
+  IF(q(1,1)<0 .or. q(1,JNP)<0.) ipx = 1
+  endif
+  RETURN
+END SUBROUTINE filew
 !
 SUBROUTINE zflip(q,im,km,nc)
   IMPLICIT NONE
   ! This routine flip the array q (in the vertical).
-  integer :: im,km,nc
-  real :: q(im,km,nc)
+  INTEGER :: im,km,nc
+  REAL :: q(im,km,nc)
   ! local dynamic array
-  real :: qtmp(im,km)
-  integer :: IC,k,i
+  REAL :: qtmp(im,km)
+  INTEGER :: IC,k,i
   !
   do IC = 1, nc
@@ -2145,4 +2144,4 @@
   END DO
   END DO
-  return
-end subroutine zflip
+  RETURN
+END SUBROUTINE zflip
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/prather.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/prather.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/prather.f90	(revision 5116)
@@ -33,10 +33,10 @@
   REAL :: q( iip1,jjp1,llm,0:9)
   REAL :: w( ip1jmp1,llm )
-  integer :: ordre,ilim
+  INTEGER :: ordre,ilim
 
   !   Local:
   !   ------
   LOGICAL :: limit
-  real :: zq(iip1,jjp1,llm)
+  REAL :: zq(iip1,jjp1,llm)
   REAL :: sm ( iip1,jjp1, llm )
   REAL :: s0( iip1,jjp1,llm ),  sx( iip1,jjp1,llm )
@@ -49,16 +49,16 @@
   REAL :: szz( iip1,jjp1,llm ),zz
   INTEGER :: i,j,l,indice
-  real :: sxn(iip1),sxs(iip1)
-
-  real :: sinlon(iip1),sinlondlon(iip1)
-  real :: coslon(iip1),coslondlon(iip1)
-  real :: qmin,qmax
+  REAL :: sxn(iip1),sxs(iip1)
+
+  REAL :: sinlon(iip1),sinlondlon(iip1)
+  REAL :: coslon(iip1),coslondlon(iip1)
+  REAL :: qmin,qmax
   save qmin,qmax
   save sinlon,coslon,sinlondlon,coslondlon
-  real :: dyn1,dyn2,dys1,dys2,qpn,qps,dqzpn,dqzps
-  real :: masn,mass
+  REAL :: dyn1,dyn2,dys1,dys2,qpn,qps,dqzpn,dqzps
+  REAL :: masn,mass
   !
   REAL :: SSUM
-  integer :: ismax,ismin
+  INTEGER :: ismax,ismin
   EXTERNAL  SSUM, ismin,ismax
   logical :: first
@@ -80,5 +80,5 @@
   limit = .TRUE.
 
-  if(first) then
+  IF(first) THEN
      PRINT*,'SCHEMA PRATHER'
      first=.FALSE.
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/principal_cshift_m.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/principal_cshift_m.F90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/principal_cshift_m.F90	(revision 5116)
@@ -11,6 +11,6 @@
     ! xprimm.
 
-    use nrtype, only: twopi
-    use serre_mod, only: clon
+    use nrtype, ONLY: twopi
+    use serre_mod, ONLY: clon
 
     include "dimensions.h"
@@ -22,5 +22,5 @@
     !-----------------------------------------------------
 
-    if (is2 /= 0) then
+    if (is2 /= 0) THEN
        IF (clon <= 0.) THEN
           IF (is2 /= 1) THEN
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/sortvarc.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/sortvarc.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/sortvarc.f90	(revision 5116)
@@ -70,6 +70,6 @@
   ! Ehouarn: when no initialization fields from file, resetvarc should be
        ! set to false
-   if (firstcal) then
-     if (.not.read_start) then
+   if (firstcal) THEN
+     if (.not.read_start) THEN
        resetvarc=.TRUE.
      endif
@@ -147,5 +147,5 @@
   ang   = SSUM(     llm,  angl, 1 )
 
-  IF (firstcal.and.resetvarc) then
+  IF (firstcal.and.resetvarc) THEN
      WRITE(lunout,3500) itau, rjour, heure, time
      WRITE(lunout,*) trim(modname), &
@@ -162,5 +162,5 @@
   ! compute relative changes in etot,... (except if 'reference' values
   ! are zero, which can happen when using iniacademic)
-  if (etot0/=0) then
+  if (etot0/=0) THEN
     etot= etot/etot0
   else
@@ -168,20 +168,20 @@
   endif
   rmsv= SQRT(rmsv/ptot)
-  if (ptot0/=0) then
+  if (ptot0/=0) THEN
     ptot= ptot/ptot0
   else
     ptot=1.
   endif
-  if (ztot0/=0) then
+  if (ztot0/=0) THEN
     ztot= ztot/ztot0
   else
     ztot=1.
   endif
-  if (stot0/=0) then
+  if (stot0/=0) THEN
     stot= stot/stot0
   else
     stot=1.
   endif
-  if (ang0/=0) then
+  if (ang0/=0) THEN
     ang = ang /ang0
   else
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/test_period.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/test_period.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/test_period.f90	(revision 5116)
@@ -44,5 +44,5 @@
      do ij=1,iim
       if (teta(ij,l)/=teta(1,l) &
-            .or.teta(ip1jm+ij,l)/=teta(ip1jm+1,l) ) then
+            .or.teta(ip1jm+ij,l)/=teta(ip1jm+1,l) ) THEN
       PRINT *,'STOP dans test_period car ---  TETA  ---  n est pas', &
             ' constant aux poles ! '
@@ -100,5 +100,5 @@
      do ij=1,iim
       if (p(ij,l)/=p(1,l) &
-            .or.p(ip1jm+ij,l)/=p(ip1jm+1,l) ) then
+            .or.p(ip1jm+ij,l)/=p(ip1jm+1,l) ) THEN
       PRINT *,'STOP dans test_period car ---  P     ---  n est pas', &
             ' constant aux poles ! '
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/traceurpole.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/traceurpole.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/traceurpole.f90	(revision 5116)
@@ -12,16 +12,16 @@
 
   !   Arguments
-  integer :: iq
-  real :: masse(iip1, jjp1, llm)
-  real :: q(iip1, jjp1, llm)
+  INTEGER :: iq
+  REAL :: masse(iip1, jjp1, llm)
+  REAL :: q(iip1, jjp1, llm)
 
 
   !   Locals
-  integer :: i, j, l
-  real :: sommemassen(llm)
-  real :: sommemqn(llm)
-  real :: sommemasses(llm)
-  real :: sommemqs(llm)
-  real :: qpolen(llm), qpoles(llm)
+  INTEGER :: i, j, l
+  REAL :: sommemassen(llm)
+  REAL :: sommemqn(llm)
+  REAL :: sommemasses(llm)
+  REAL :: sommemqs(llm)
+  REAL :: qpolen(llm), qpoles(llm)
 
 
@@ -56,4 +56,4 @@
   enddo
 
-  return
-end subroutine traceurpole
+  RETURN
+END SUBROUTINE traceurpole
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/ugeostr.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/ugeostr.F90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/ugeostr.F90	(revision 5116)
@@ -11,5 +11,5 @@
   ! levels are pressure levels.
 
-  use comconst_mod, only: omeg, rad
+  use comconst_mod, ONLY: omeg, rad
   
   IMPLICIT NONE
@@ -29,5 +29,5 @@
   DO j=1,jjm
 
-     if (abs(sin(rlatv(j)))<1.e-4) then
+     if (abs(sin(rlatv(j)))<1.e-4) THEN
         zlat=1.e-4
      else
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/write_grads_dyn.h
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/write_grads_dyn.h	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/write_grads_dyn.h	(revision 5116)
@@ -2,6 +2,5 @@
 ! $Header$
 
-if (callinigrads) then
-
+if (callinigrads) THEN
    string10='dyn'
    CALL inigrads(1,iip1 &
@@ -26,5 +25,5 @@
 do iq=1,nqtot
    string10='q'
-   write(string10(2:2),'(i1)') iq
+   WRITE(string10(2:2),'(i1)') iq
    CALL wrgrads(1,llm,q(:,:,iq),string10,string10)
 enddo
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writedynav.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writedynav.F90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writedynav.F90	(revision 5116)
@@ -120,5 +120,5 @@
   ! CALL histwrite(histaveid, 'phis', itau_w, phis, iip1*jjp1, ndex2d)
 
-  if (ok_sync) then
+  if (ok_sync) THEN
      CALL histsync(histaveid)
      CALL histsync(histvaveid)
Index: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writehist.f90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writehist.f90	(revision 5114)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writehist.f90	(revision 5116)
@@ -47,5 +47,5 @@
   REAL :: phis(ip1jmp1)
   REAL :: q(ip1jmp1, llm, nqtot)
-  integer :: time
+  INTEGER :: time
 
 
@@ -53,8 +53,8 @@
   !   Variables locales
   !
-  integer :: iq, ii, ll
-  integer :: ndexu(ip1jmp1 * llm), ndexv(ip1jm * llm), ndex2d(ip1jmp1)
+  INTEGER :: iq, ii, ll
+  INTEGER :: ndexu(ip1jmp1 * llm), ndexv(ip1jm * llm), ndex2d(ip1jmp1)
   logical :: ok_sync
-  integer :: itau_w
+  INTEGER :: itau_w
   REAL :: vnat(ip1jm, llm), unat(ip1jmp1, llm)
 
@@ -114,9 +114,9 @@
   !  Fin
   !
-  if (ok_sync) then
+  if (ok_sync) THEN
     CALL histsync(histid)
     CALL histsync(histvid)
     CALL histsync(histuid)
   endif
-  return
-end subroutine writehist
+  RETURN
+END SUBROUTINE writehist
