Index: /LMDZ6/branches/Amaury_dev/libf/dyn3d/leapfrog.F
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3d/leapfrog.F	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3d/leapfrog.F	(revision 5081)
@@ -301,20 +301,20 @@
       IF( purmats ) THEN
       ! Purely Matsuno time stepping
-         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
-         IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward ) 
+         IF( MOD(itau,iconser) ==0.AND.  forward    ) conser = .TRUE.
+         IF( MOD(itau,dissip_period )==0.AND..NOT.forward )
      s        apdiss = .TRUE.
-         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward 
+         IF( MOD(itau,iphysiq )==0.AND..NOT.forward
      s          .and. physic                        ) apphys = .TRUE.
       ELSE
       ! Leapfrog/Matsuno time stepping 
-         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
-         IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
+         IF( MOD(itau   ,iconser) == 0              ) conser = .TRUE.
+         IF( MOD(itau+1,dissip_period)==0 .AND. .NOT. forward )
      s        apdiss = .TRUE.
-         IF( MOD(itau+1,iphysiq).EQ.0.AND.physic       ) apphys=.TRUE.
+         IF( MOD(itau+1,iphysiq)==0.AND.physic       ) apphys=.TRUE.
       END IF
 
 ! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
 !          supress dissipation step
-      if (llm.eq.1) then
+      if (llm==1) then
         apdiss=.false.
       endif
@@ -387,5 +387,5 @@
 c
        IF( purmats )  THEN
-          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
+          IF( itau==itaufin.AND..NOT.forward ) lafin = .TRUE.
        ELSE
           IF( itau+1. EQ. itaufin )              lafin = .TRUE.
@@ -418,5 +418,5 @@
      &          (itau+1)/day_step
 
-           IF (planet_type .eq."generic") THEN
+           IF (planet_type =="generic") THEN
               ! AS: we make jD_cur to be pday
               jD_cur = int(day_ini + itau/day_step)
@@ -441,8 +441,8 @@
 
 c  Diagnostique de conservation de l'energie : initialisation
-         IF (ip_ebil_dyn.ge.1 ) THEN 
+         IF (ip_ebil_dyn>=1 ) THEN
           ztit='bil dyn'
 ! Ehouarn: be careful, diagedyn is Earth-specific!
-           IF (planet_type.eq."earth") THEN
+           IF (planet_type=="earth") THEN
             CALL diagedyn(ztit,2,1,1,dtphys
      &    , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
@@ -488,7 +488,7 @@
 c
 c  Diagnostique de conservation de l'energie : difference
-         IF (ip_ebil_dyn.ge.1 ) THEN 
+         IF (ip_ebil_dyn>=1 ) THEN
           ztit='bil phys'
-          IF (planet_type.eq."earth") THEN
+          IF (planet_type=="earth") THEN
            CALL diagedyn(ztit,2,1,1,dtphys
      &     , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
@@ -498,5 +498,5 @@
        ENDIF ! of IF( apphys )
 
-      IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
+      IF(iflag_phys==2) THEN ! "Newtonian" case
 !   Academic case : Simple friction and Newtonan relaxation 
 !   -------------------------------------------------------
@@ -508,5 +508,5 @@
         ENDDO ! of DO l=1,llm 
         
-        if (planet_type.eq."giant") then
+        if (planet_type=="giant") then
           ! add an intrinsic heat flux at the base of the atmosphere
           teta(:,1)=teta(:,1)+dtvr*aire(:)*ihf/cpp/masse(:,1)
@@ -663,6 +663,6 @@
 c   -------------------------------------
 
-            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
-               IF(itau.EQ.itaufin) THEN
+            IF(MOD(itau,iperiod)==0 .OR. itau==itaufin) THEN
+               IF(itau==itaufin) THEN
                   iav=1
                ELSE
@@ -694,5 +694,5 @@
 c   ------------------------------
 
-            IF( MOD(itau,iecri).EQ.0) THEN
+            IF( MOD(itau,iecri)==0) THEN
              ! Ehouarn: output only during LF or Backward Matsuno
              if (leapf.or.(.not.leapf.and.(.not.forward))) then
@@ -721,5 +721,5 @@
             ENDIF ! of IF(MOD(itau,iecri).EQ.0)
 
-            IF(itau.EQ.itaufin) THEN
+            IF(itau==itaufin) THEN
 
 
@@ -743,5 +743,5 @@
 c   ------------------------------------
 
-            IF( MOD(itau,iperiod).EQ.0 )    THEN
+            IF( MOD(itau,iperiod)==0 )    THEN
                     GO TO 1
             ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limy.F
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limy.F	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3d_common/limy.F	(revision 5081)
@@ -150,5 +150,5 @@
 c      enddo
 
-      if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
+      if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1))<=0.)
      &   then
          do ij=1,iip1
@@ -162,5 +162,5 @@
 
       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.)
+     & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)<=0.)
      &then
          do ij=ip1jm+1,ip1jmp1
@@ -175,6 +175,6 @@
 c   calcul des pentes limitees
 
-      do ij=1,ip1jmp1
-         if(dyqv(ij)*dyqv(ij-iip1).gt.0.) then
+      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
             dyq(ij)=sign(min(abs(dyq(ij)),dyqmax(ij)),dyq(ij))
          else
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3d_common/massbarxy.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3d_common/massbarxy.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3d_common/massbarxy.F90	(revision 5081)
@@ -21,7 +21,7 @@
     DO ij=1,ip1jm-1
       massebxy(ij,l)=masse(ij     ,l)*alpha2(ij     ) + &
-     +               masse(ij+1   ,l)*alpha3(ij+1   ) + &
-     +               masse(ij+iip1,l)*alpha1(ij+iip1) + &
-     +               masse(ij+iip2,l)*alpha4(ij+iip2)
+                     masse(ij+1   ,l)*alpha3(ij+1   ) + &
+                     masse(ij+iip1,l)*alpha1(ij+iip1) + &
+                     masse(ij+iip2,l)*alpha4(ij+iip2)
     END DO
     DO ij=iip1,ip1jm,iip1; massebxy(ij,l)=massebxy(ij-iim,l); END DO
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/guide_loc_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/guide_loc_mod.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/guide_loc_mod.F90	(revision 5081)
@@ -1598,8 +1598,8 @@
     if (first) then
          ncidpl=-99
-         write(*,*),trim(modname)//': opening nudging files '
+         write(*,*) trim(modname)//': opening nudging files '
 ! Ap et Bp si Niveaux de pression hybrides
          if (guide_plevs==1) then
-             write(*,*),trim(modname)//' Reading nudging on model levels'
+             write(*,*) trim(modname)//' Reading nudging on model levels'
              rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
              IF (rcode/=NF_NOERR) THEN
@@ -1617,5 +1617,5 @@
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
-             write(*,*),trim(modname)//' ncidpl,varidap',ncidpl,varidap
+             write(*,*) trim(modname)//' ncidpl,varidap',ncidpl,varidap
          endif
          
@@ -1632,5 +1632,5 @@
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
-             write(*,*),trim(modname)//' ncidp,varidp',ncidp,varidp
+             write(*,*) trim(modname)//' ncidp,varidp',ncidp,varidp
              if (ncidpl==-99) ncidpl=ncidp
          endif
@@ -1648,5 +1648,5 @@
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
-             write(*,*),trim(modname)//' ncidu,varidu',ncidu,varidu
+             write(*,*) trim(modname)//' ncidu,varidu',ncidu,varidu
              if (ncidpl==-99) ncidpl=ncidu
 
@@ -1680,5 +1680,5 @@
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
-             write(*,*),trim(modname)//' ncidv,varidv',ncidv,varidv
+             write(*,*) trim(modname)//' ncidv,varidv',ncidv,varidv
              if (ncidpl==-99) ncidpl=ncidv
              
@@ -1713,5 +1713,5 @@
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
-             write(*,*),trim(modname)//' ncidT,varidT',ncidt,varidt
+             write(*,*) trim(modname)//' ncidT,varidT',ncidt,varidt
              if (ncidpl==-99) ncidpl=ncidt
 
@@ -1744,5 +1744,5 @@
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
-             write(*,*),trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
+             write(*,*) trim(modname)//' ncidQ,varidQ',ncidQ,varidQ
              if (ncidpl==-99) ncidpl=ncidQ
 
@@ -1776,5 +1776,5 @@
               CALL abort_gcm(modname,abort_message,1)
              ENDIF
-             write(*,*),trim(modname)//' ncidps,varidps',ncidps,varidps
+             write(*,*) trim(modname)//' ncidps,varidps',ncidps,varidps
          endif
 ! Coordonnee verticale
@@ -1782,5 +1782,5 @@
               rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
               IF (rcode/=0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
-              write(*,*),trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
+              write(*,*) trim(modname)//' ncidpl,varidpl',ncidpl,varidpl
          endif
 ! Coefs ap, bp pour calcul de la pression aux differents niveaux
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/massbarxy_loc.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/massbarxy_loc.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/massbarxy_loc.F90	(revision 5081)
@@ -27,7 +27,7 @@
     DO ij=ijb,ije-1
       massebxy(ij,l)=masse(ij     ,l)*alpha2(ij     ) + &
-     +               masse(ij+1   ,l)*alpha3(ij+1   ) + &
-     +               masse(ij+iip1,l)*alpha1(ij+iip1) + &
-     +               masse(ij+iip2,l)*alpha4(ij+iip2)
+                     masse(ij+1   ,l)*alpha3(ij+1   ) + &
+                     masse(ij+iip1,l)*alpha1(ij+iip1) + &
+                     masse(ij+iip2,l)*alpha4(ij+iip2)
     END DO
     DO ij=ijb+iip1-1,ije+iip1-1,iip1; massebxy(ij,l)=massebxy(ij-iim,l); END DO
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/vlspltqs_loc.F
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/vlspltqs_loc.F	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/vlspltqs_loc.F	(revision 5081)
@@ -62,5 +62,5 @@
       if (pole_sud.and.ije==ip1jmp1)  ije=ije-iip1
       
-      IF (pente_max.gt.-1.e-5) THEN
+      IF (pente_max>-1.e-5) THEN
 c     IF (pente_max.gt.10) THEN
 
@@ -104,5 +104,5 @@
      ,         cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
 #else
-               IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
+               IF(dxqu(ij-1)*dxqu(ij)>0) THEN
                   dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
                ELSE
@@ -135,5 +135,5 @@
                zz(ij)=dxqu(ij-1)*dxqu(ij)
                zz(ij)=zz(ij)+zz(ij)
-               IF(zz(ij).gt.0) THEN
+               IF(zz(ij)>0) THEN
                   dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
                ELSE
@@ -205,5 +205,5 @@
       DO l=1,llm
        DO ij=ijb,ije-1
-          IF (u_m(ij,l).gt.0.) THEN
+          IF (u_m(ij,l)>0.) THEN
              zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
              u_mq(ij,l)=u_m(ij,l)*
@@ -225,5 +225,5 @@
       DO l=1,llm
          DO ij=ijb,ije-1
-            IF(zdum(ij,l).lt.0) THEN
+            IF(zdum(ij,l)<0) THEN
                iadvplus(ij,l)=1
                u_mq(ij,l)=0.
@@ -269,9 +269,9 @@
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
          DO l=1,llm
-            IF(nl(l).gt.0) THEN
+            IF(nl(l)>0) THEN
                iju=0
 c   indicage des mailles concernees par le traitement special
                DO ij=ijb,ije
-                  IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
+                  IF(iadvplus(ij,l)==1.and.mod(ij,iip1)/=0) THEN
                      iju=iju+1
                      indu(iju)=ij
@@ -287,9 +287,9 @@
                   zu_m=u_m(ij,l)
                   u_mq(ij,l)=0.
-                  IF(zu_m.gt.0.) THEN
+                  IF(zu_m>0.) THEN
                      ijq=ij
                      i=ijq-(j-1)*iip1
 c   accumulation pour les mailles completements advectees
-                     do while(zu_m.gt.masse(ijq,l,iq))
+                     do while(zu_m>masse(ijq,l,iq))
                         u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
      &                     *masse(ijq,l,iq) 
@@ -305,5 +305,5 @@
                      i=ijq-(j-1)*iip1
 c   accumulation pour les mailles completements advectees
-                     do while(-zu_m.gt.masse(ijq,l,iq))
+                     do while(-zu_m>masse(ijq,l,iq))
                         u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
      &                   *masse(ijq,l,iq)
@@ -345,5 +345,5 @@
             !MVals: veiller a ce qu'on n'ait pas de denominateur nul
             masse(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
-            if (q(ij,l,iq).gt.min_qParent) then ! modif 13 nov 2020
+            if (q(ij,l,iq)>min_qParent) then ! modif 13 nov 2020
               Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
             else
@@ -478,5 +478,5 @@
       ij=3525
       l=3
-      if ((ij.ge.ijb).and.(ij.le.ije)) then
+      if ((ij>=ijb).and.(ij<=ije)) then
         !write(*,*) 'vlyqs 480: ij,l,iq,ijb,q(ij,l,:)=',
 !     &             ij,l,iq,ijb,q(ij,l,:)
@@ -576,5 +576,5 @@
         fn=1.
         DO ij=1,iim
-          IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
+          IF(pente_max*adyqv(ij)<abs(dyq(ij,l))) THEN
             fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
           ENDIF
@@ -608,5 +608,5 @@
         fs=1.
         DO ij=1,iim
-        IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
+        IF(pente_max*adyqv(ij+ip1jm-iip1)<abs(dyq(ij+ip1jm,l))) THEN
          fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
         ENDIF
@@ -694,5 +694,5 @@
 
       DO ij=ijb,ije
-         IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
+         IF(dyqv(ij)*dyqv(ij-iip1)>0.) THEN
             dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
          ELSE
@@ -712,5 +712,5 @@
       DO l=1,llm
        DO ij=ijb,ije
-         IF( masse_adv_v(ij,l).GT.0. ) THEN
+         IF( masse_adv_v(ij,l)>0. ) THEN
            qbyv(ij,l,iq)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq )  +
      ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)
@@ -757,5 +757,5 @@
             !MVals: veiller a ce qu'on n'ait pas de denominateur nul
             !write(lunout,*) 'ij,l,q(ij,l,iq)=',ij,l,q(ij,l,iq)
-            if (q(ij,l,iq).gt.min_qParent) then ! modif 13 nov 2020
+            if (q(ij,l,iq)>min_qParent) then ! modif 13 nov 2020
               Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
             else
@@ -806,5 +806,5 @@
          IF (pole_sud) THEN
          
-           convps  = -SSUM(iim,qbyv(ip1jm-iim,l,iq),iq,1)/apols
+           convps  = -SSUM(iim,qbyv(ip1jm-iim,l,iq),iq,1)/apols  ! /!\ called with 4 args ???
            convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
            DO ij = ip1jm+1,ip1jmp1
Index: /LMDZ6/branches/Amaury_dev/libf/misc/cray.F
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/misc/cray.F	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/misc/cray.F	(revision 5081)
@@ -15,9 +15,9 @@
       iy=1
       ix=1
-      do 10 i=1,n
+      DO i=1,n
          sy(iy)=sx(ix)
          ix=ix+incx
          iy=iy+incy
-10    continue
+      END DO
 c
       return
@@ -33,8 +33,8 @@
       ssum=0.
       ix=1
-      do 10 i=1,n
+      do i=1,n
          ssum=ssum+sx(ix)
          ix=ix+incx
-10    continue
+      end do
 c
       return
Index: /LMDZ6/branches/Amaury_dev/libf/phy_common/mod_phys_lmdz_mpi_transfert.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phy_common/mod_phys_lmdz_mpi_transfert.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phy_common/mod_phys_lmdz_mpi_transfert.F90	(revision 5081)
@@ -65,5 +65,5 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-!! -- Les chaine de charactère -- !!
+!! -- Les chaine de charact�re -- !!
 
   SUBROUTINE bcast_mpi_c(var1)
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/cond_evap_tstep_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/cond_evap_tstep_mod.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/cond_evap_tstep_mod.F90	(revision 5081)
@@ -227,5 +227,5 @@
 !     H2SO4 mass fraction in aerosol
       WH2=R2SO4*1.0E-2
-      IF(WH2.EQ.0.0) RETURN
+      IF(WH2==0.0) RETURN
 !                               ACTIVITY COEFFICIENT(SEE GIAUQUE,1951)
 !                               AYERS ET AL (1980)
@@ -324,18 +324,18 @@
         ff(IK,:)=0.0 
         DO k=1, nbtr_bin 
-          IF (k.LE.(nbtr_bin-1)) THEN 
-            IF (Vbin_wet(k).LE.Vnew.AND.Vnew.LT.Vbin_wet(k+1)) THEN
+          IF (k<=(nbtr_bin-1)) THEN
+            IF (Vbin_wet(k)<=Vnew.AND.Vnew<Vbin_wet(k+1)) THEN
               ff(IK,k)= Vbin_wet(k)/Vnew*(Vbin_wet(k+1)-Vnew)/(Vbin_wet(k+1)-Vbin_wet(k))
             ENDIF
           ENDIF
-          IF (k.EQ.1.AND.Vnew.LE.Vbin_wet(k)) THEN
+          IF (k==1.AND.Vnew<=Vbin_wet(k)) THEN
             ff(IK,k)= 1.
           ENDIF
-          IF (k.GT.1) THEN 
-            IF (Vbin_wet(k-1).LT.Vnew.AND.Vnew.LT.Vbin_wet(k)) THEN
+          IF (k>1) THEN
+            IF (Vbin_wet(k-1)<Vnew.AND.Vnew<Vbin_wet(k)) THEN
               ff(IK,k)= 1.-ff(IK,k-1)
             ENDIF
           ENDIF
-          IF (k.EQ.nbtr_bin.AND.Vnew.GE.Vbin_wet(k)) THEN
+          IF (k==nbtr_bin.AND.Vnew>=Vbin_wet(k)) THEN
             ff(IK,k)= 1.
           ENDIF
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/YOECUMF.h
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/YOECUMF.h	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/YOECUMF.h	(revision 5081)
@@ -12,7 +12,8 @@
 !
       COMMON /YOECUMF/                                                  &
-     &                 LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV,              &
      &                 ENTRPEN,ENTRSCV,ENTRMID,ENTRDD,CMFCTOP,          &
-     &                 CMFCMAX,CMFCMIN,CMFDEPS,RHCDD,CPRCON
+     &                 CMFCMAX,CMFCMIN,CMFDEPS,RHCDD,CPRCON,            &
+     &                 LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV
+
 
       LOGICAL          LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/aaam_bud.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/aaam_bud.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/aaam_bud.F90	(revision 5081)
@@ -127,5 +127,5 @@
   hadley = 1.E18
   hadday = 1.E18*24.*3600.
-  IF(klon_glo.EQ.1) THEN
+  IF(klon_glo==1) THEN
     dlat = xpi
   ELSE
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/acama_gwd_rando_m.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/acama_gwd_rando_m.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/acama_gwd_rando_m.F90	(revision 5081)
@@ -358,5 +358,5 @@
          CMAX*2.*(MOD(TT(II, JW+4*(JJ-1)+JJ)**2, 1.)-0.5)*SQRT(3.)/SQRT(NA*1.)
                 END DO
-                IF (CPHA.LT.0.)  THEN
+                IF (CPHA<0.)  THEN
                    CPHA = -1.*CPHA
                    ZP(JW,II) = ZP(JW,II) + RPI
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/add_phys_tend_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/add_phys_tend_mod.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/add_phys_tend_mod.F90	(revision 5081)
@@ -226,5 +226,5 @@
   ENDDO
 
-  IF (fl_ebil .GT. 0) THEN
+  IF (fl_ebil > 0) THEN
     ! ------------------------------------------------
     ! Compute vertical sum for each atmospheric column
@@ -283,8 +283,8 @@
 !=====================================================================================
 
-  IF (jbad .GT. 0) THEN
+  IF (jbad > 0) THEN
      DO j = 1, jbad
         i=jadrs(j)
-        IF (prt_level.ge.debug_level) THEN
+        IF (prt_level>=debug_level) THEN
           print*,'PLANTAGE POUR LE POINT i lon lat =',&
                  i,longitude_deg(i),latitude_deg(i),text
@@ -301,9 +301,9 @@
 ! Impression, warning et correction en cas de probleme moins important
 !=====================================================================================
-  IF (jqbad .GT. 0) THEN
+  IF (jqbad > 0) THEN
       done(:) = .false.                         !jyg
       DO j = 1, jqbad
         i=jqadrs(j)
-          if(prt_level.ge.debug_level) THEN
+          if(prt_level>=debug_level) THEN
            print*,'WARNING  : EAU POUR LE POINT i lon lat =',&
                   i,longitude_deg(i),latitude_deg(i),text
@@ -325,5 +325,5 @@
                 zqp_int = zqp_int + zqp(k)     *(paprs(i,k)-paprs(i,k+1))/Rg
               ENDDO
-              IF (prt_level.ge.debug_level) THEN
+              IF (prt_level>=debug_level) THEN
                print*,' cas q_seri<1.e-15 i k zq_int zqp_int zq_int/zqp_int :', &
                                     i, kqadrs(j), zq_int, zqp_int, zq_int/zqp_int
@@ -340,7 +340,7 @@
             DO k = 1, klev
               zq=q_seri(i,k)+zdq(i,k)
-              IF (zq.lt.1.e-15) THEN
-                 IF (q_seri(i,k).lt.1.e-15) THEN
-                   IF (prt_level.ge.debug_level) THEN
+              IF (zq<1.e-15) THEN
+                 IF (q_seri(i,k)<1.e-15) THEN
+                   IF (prt_level>=debug_level) THEN
                     print*,' cas q_seri<1.e-15 i k q_seri zq zdq :',i,k,q_seri(i,k),zq,zdq(i,k)
                    ENDIF
@@ -383,9 +383,9 @@
      ENDDO
   ENDDO
-  IF (jbad .GT. 0) THEN
+  IF (jbad > 0) THEN
       DO j = 1, jbad
          i=jadrs(j)
          k=kadrs(j)
-         if(prt_level.ge.debug_level) THEN
+         if(prt_level>=debug_level) THEN
           print*,'PLANTAGE2 POUR LE POINT i itap lon lat txt jbad zdt t',&
                  i,itap,longitude_deg(i),latitude_deg(i),text,jbad, &
@@ -401,9 +401,9 @@
   ENDIF 
 !
-  IF (jqbad .GT. 0) THEN
+  IF (jqbad > 0) THEN
       DO j = 1, jqbad
          i=jqadrs(j)
          k=kqadrs(j)
-         IF (prt_level.ge.debug_level) THEN
+         IF (prt_level>=debug_level) THEN
           print*,'WARNING  : EAU2 POUR LE POINT i itap lon lat txt jqbad zdq q zdql ql',&
                  i,itap,longitude_deg(i),latitude_deg(i),text,jqbad,&
@@ -441,5 +441,5 @@
 !======================================================================
 
-  IF (fl_ebil .GT. 0) THEn
+  IF (fl_ebil > 0) THEn
   
     ! ------------------------------------------------
@@ -584,5 +584,5 @@
   ENDDO
 
-  IF (fl_ebil .GT. 0) THEN
+  IF (fl_ebil > 0) THEN
     ! ------------------------------------------------
     ! Compute vertical sum for each atmospheric column
@@ -613,5 +613,5 @@
 !======================================================================
 
-  IF (fl_ebil .GT. 0) THEN
+  IF (fl_ebil > 0) THEN
   
     ! ------------------------------------------------
@@ -761,5 +761,5 @@
 
 !!print *,'prt_level:',prt_level,' fl_ebil:',fl_ebil,' fl_cor_ebil:',fl_cor_ebil
-IF ((fl_ebil .GT. 0) .AND. (klon .EQ. 1)) THEN
+IF ((fl_ebil > 0) .AND. (klon == 1)) THEN
 
   bilq_bnd = 0.
@@ -793,10 +793,10 @@
   bilh_error = d_h_col(1) - bilh_bnd
 ! are the errors too large?
-  IF (abs(bilq_error) .GT. bilq_seuil) bilq_ok=1
-  IF (abs(bilh_error) .GT. bilh_seuil) bilh_ok=1
+  IF (abs(bilq_error) > bilq_seuil) bilq_ok=1
+  IF (abs(bilh_error) > bilh_seuil) bilh_ok=1
 !
 ! Print diagnostics
 ! =================
-  IF ( (bilq_ok .eq. 0).AND.(bilh_ok .eq. 0) ) THEN
+  IF ( (bilq_ok == 0).AND.(bilh_ok == 0) ) THEN
     status="enerbil-OK"
   ELSE
@@ -804,12 +804,12 @@
   ENDIF
 
-  IF (prt_level .GE. 3) THEN
+  IF (prt_level >= 3) THEN
     write(*,9010) text,status," itap:",itap,"enerbilERROR: Q", bilq_error,"  H", bilh_error
 9010  format (1x,A8,2x,A12,A6,I4,A18,E15.6,A5,E15.6)
   ENDIF
-  IF (prt_level .GE. 3) THEN
+  IF (prt_level >= 3) THEN
     write(*,9000) text,"enerbil: Q,H,KE budget", d_qt_col(1), d_h_col(1),d_ek_col(1)
   ENDIF
-  IF (prt_level .GE. 5) THEN
+  IF (prt_level >= 5) THEN
     write(*,9000) text,"enerbil at boundaries: Q, H",bilq_bnd, bilh_bnd
     write(*,9000) text,"enerbil: water budget",d_qt_col(1),d_qw_col(1),d_ql_col(1),d_qs_col(1), d_qbs_col(1)
@@ -819,15 +819,15 @@
   specific_diag: SELECT CASE (text)
   CASE("vdf") specific_diag
-    IF (prt_level .GE. 5) THEN
+    IF (prt_level >= 5) THEN
       write(*,9000) text,"enerbil: d_h, bilh, sens,t_seri", d_h_col(1), bilh_bnd, sens(1), t_seri(1,1)
       write(*,9000) text,"enerbil: d_h_col_vdf, f_h, diff",d_h_col_vdf, f_h_bnd, bilh_bnd-sens(1)
     ENDIF
   CASE("lsc") specific_diag
-    IF (prt_level .GE. 5) THEN
+    IF (prt_level >= 5) THEN
       write(*,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_lsc(1), rlvtt * rain_lsc(1), -(rcw-rcpd)*t_seri(1,1) * rain_lsc(1)
       write(*,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_lsc(1), rlstt * snow_lsc(1), -(rcs-rcpd)*t_seri(1,1) * snow_lsc(1)
     ENDIF
   CASE("convection") specific_diag
-    IF (prt_level .GE. 5) THEN
+    IF (prt_level >= 5) THEN
       write(*,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_con(1), rlvtt * rain_con(1), -(rcw-rcpd)*t_seri(1,1) * rain_con(1)
       write(*,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_con(1), rlstt * snow_con(1), -(rcs-rcpd)*t_seri(1,1) * snow_con(1)
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/add_wake_tend.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/add_wake_tend.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/add_wake_tend.F90	(revision 5081)
@@ -1,5 +1,5 @@
 SUBROUTINE add_wake_tend(zddeltat, zddeltaq, zds, zdas, zddensw, zddensaw, zoccur, text, abortphy)
 !===================================================================
-! Ajoute les tendances liées aux diverses parametrisations physiques aux
+! Ajoute les tendances li�es aux diverses parametrisations physiques aux
 ! variables d'etat des poches froides.
 !===================================================================
@@ -43,5 +43,5 @@
          DO l = 1, klev
            DO i = 1, klon
-             IF (zoccur(i) .GE. 1) THEN
+             IF (zoccur(i) >= 1) THEN
                wake_deltat(i, l) = wake_deltat(i, l) + zddeltat(i,l)
                wake_deltaq(i, l) = wake_deltaq(i, l) + zddeltaq(i,l)
@@ -53,5 +53,5 @@
          END DO
          DO i = 1, klon
-           IF (zoccur(i) .GE. 1) THEN
+           IF (zoccur(i) >= 1) THEN
              wake_s(i)     = wake_s(i)    + zds(i)
              awake_s(i)    = awake_s(i)    + zdas(i)
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/alpale.h
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/alpale.h	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/alpale.h	(revision 5081)
@@ -14,5 +14,5 @@
 
       common/calpale1/iflag_trig_bl,iflag_clos_bl,tau_trig_shallow,tau_trig_deep,iflag_strig
-      common/calpale2/s_trig,iflag_coupl,iflag_clos,iflag_wake,alp_bl_k,h_trig
+      common/calpale2/alp_bl_k,s_trig,h_trig,iflag_coupl,iflag_clos,iflag_wake
 
 !$OMP THREADPRIVATE(/calpale1/,/calpale2/)
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/alpale_wk.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/alpale_wk.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/alpale_wk.F90	(revision 5081)
@@ -58,5 +58,5 @@
 !!  print *,'alpale_wk: sigmaw(1), wdens(1) ', sigmaw(1), wdens(1)
   DO i = 1,klon
-    IF (zoccur(i) .GE. 1) THEN
+    IF (zoccur(i) >= 1) THEN
       wkrad(i) = sqrt(sigmaw(i)/(rpi*wdens(i)))
     ELSE
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/clesphys.h
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/clesphys.h	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/clesphys.h	(revision 5081)
@@ -133,4 +133,5 @@
      &     , qsol0,albsno0,evap0                                        &
      &     , co2_ppm0                                                   &
+     &     , tau_thermals                                               &
 !FC
      &     , Cd_frein,zrel_oro_t,zpmm_orodr_t,zpmm_orolf_t,zstd_orodr_t &
@@ -163,5 +164,5 @@
      &     , ok_daily_climoz, ok_all_xml, ok_lwoff                      &
      &     , iflag_phytrac, ok_new_lscp, ok_bs, ok_rad_bs               &
-     &     ,  iflag_thermals,nsplit_thermals, tau_thermals              &
+     &     ,  iflag_thermals,nsplit_thermals              &
      &     , iflag_physiq, ok_3Deffect, ok_water_mass_fixer
        save /clesphys/
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/climb_hq_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/climb_hq_mod.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/climb_hq_mod.F90	(revision 5081)
@@ -7,5 +7,4 @@
 
   IMPLICIT NONE
-  SAVE 
   PRIVATE
   PUBLIC :: climb_hq_down, climb_hq_up, d_h_col_vdf, f_h_bnd
@@ -222,5 +221,5 @@
 !!! jyg le 07/02/2012
 !!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
-       IF (mod(iflag_pbl_split,10) .ge.1) THEN
+       IF (mod(iflag_pbl_split,10) >=1) THEN
 !!! nrlmd le 02/05/2011
     DO k= 1, klev
@@ -231,5 +230,5 @@
         Dcoef_Q_out(i,k) = Dcoef_Q(i,k)
         Kcoef_hq_out(i,k) = Kcoefhq(i,k)
-          IF (k.eq.1) THEN
+          IF (k==1) THEN
             gama_h_out(i,k)  = 0.
             gama_q_out(i,k)  = 0.
@@ -379,5 +378,5 @@
 !!! jyg le 07/02/2012
 !!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
-       IF (mod(iflag_pbl_split,10) .ge.1) THEN
+       IF (mod(iflag_pbl_split,10) >=1) THEN
 !!! nrlmd le 02/05/2011
     DO i = 1, knon
@@ -394,5 +393,5 @@
         Dcoef_Q(i,k)=Dcoef_Q_in(i,k)
         Kcoefhq(i,k)=Kcoef_hq_in(i,k)
-          IF (k.gt.1) THEN
+          IF (k>1) THEN
             gamah(i,k)=gama_h_in(i,k)
             gamaq(i,k)=gama_q_in(i,k)
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/conf_phys_m.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/conf_phys_m.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/conf_phys_m.F90	(revision 5081)
@@ -1352,5 +1352,5 @@
     !Config Help =
     !
-    ok_ice_sursat_omp = 0
+    ok_ice_sursat_omp = .FALSE.
     CALL getin('ok_ice_sursat',ok_ice_sursat_omp)
 
@@ -2574,16 +2574,16 @@
 
     !--test on radiative scheme 
-    IF (iflag_rrtm .EQ. 0) THEN 
-      IF (NSW.NE.2) THEN 
+    IF (iflag_rrtm == 0) THEN
+      IF (NSW/=2) THEN
         WRITE(lunout,*) ' ERROR iflag_rrtm=0 and NSW<>2 not possible'
         CALL abort_physic('conf_phys','choice NSW not valid',1)
       ENDIF
-    ELSE IF (iflag_rrtm .EQ. 1) THEN
-      IF (NSW.NE.2.AND.NSW.NE.4.AND.NSW.NE.6) THEN
+    ELSE IF (iflag_rrtm == 1) THEN
+      IF (NSW/=2.AND.NSW/=4.AND.NSW/=6) THEN
         WRITE(lunout,*) ' ERROR iflag_rrtm=1 and NSW<>2,4,6 not possible'
         CALL abort_physic('conf_phys','choice NSW not valid',1)
       ENDIF
-   ELSE IF (iflag_rrtm .EQ. 2) THEN
-      IF (NSW.NE.2.AND.NSW.NE.4.AND.NSW.NE.6) THEN
+   ELSE IF (iflag_rrtm == 2) THEN
+      IF (NSW/=2.AND.NSW/=4.AND.NSW/=6) THEN
         WRITE(lunout,*) ' ERROR iflag_rrtm=1 and NSW<>2,4,6 not possible'
         CALL abort_physic('conf_phys','choice NSW not valid',1)
@@ -2610,5 +2610,5 @@
 
     !--test on ocean surface albedo
-    IF (iflag_albedo.LT.0.OR.iflag_albedo.GT.2) THEN
+    IF (iflag_albedo<0.OR.iflag_albedo>2) THEN
        WRITE(lunout,*) ' ERROR iflag_albedo<>0,1'
        CALL abort_physic('conf_phys','choice iflag_albedo not valid',1)
@@ -2617,5 +2617,5 @@
     ! Flag_aerosol cannot be set to zero if aerosol direct effect (ade) or aerosol indirect effect (aie) are activated
     IF (ok_ade .OR. ok_aie) THEN 
-       IF ( flag_aerosol .EQ. 0 ) THEN
+       IF ( flag_aerosol == 0 ) THEN
           CALL abort_physic('conf_phys','flag_aerosol=0 not compatible avec ok_ade ou ok_aie=.TRUE.',1)
        ENDIF
@@ -2623,15 +2623,15 @@
 
     ! Flag_aerosol cannot be set to zero if we are in coupled mode for aerosol
-    IF (aerosol_couple .AND. flag_aerosol .EQ. 0 ) THEN
+    IF (aerosol_couple .AND. flag_aerosol == 0 ) THEN
        CALL abort_physic('conf_phys', 'flag_aerosol cannot be to zero if aerosol_couple=y ', 1)
     ENDIF
 
     ! Read_climoz needs to be set zero if we are in couple mode for chemistry 
-    IF (chemistry_couple .AND. read_climoz .ne. 0) THEN 
+    IF (chemistry_couple .AND. read_climoz /= 0) THEN
        CALL abort_physic('conf_phys', 'read_climoz need to be to zero if chemistry_couple=y ', 1) 
     ENDIF
 
     ! flag_aerosol need to be different to zero if ok_cdnc is activated
-    IF (ok_cdnc .AND. flag_aerosol .EQ. 0) THEN
+    IF (ok_cdnc .AND. flag_aerosol == 0) THEN
        CALL abort_physic('conf_phys', 'flag_aerosol cannot be to zero if ok_cdnc is activated ', 1)
     ENDIF
@@ -2643,27 +2643,27 @@
 
     ! flag_aerosol=7 => MACv2SP climatology 
-    IF (flag_aerosol.EQ.7.AND. iflag_rrtm.NE.1) THEN
+    IF (flag_aerosol==7.AND. iflag_rrtm/=1) THEN
        CALL abort_physic('conf_phys', 'flag_aerosol=7 (MACv2SP) can only be activated with RRTM',1)
     ENDIF
-    IF (flag_aerosol.EQ.7.AND. NSW.NE.6) THEN
+    IF (flag_aerosol==7.AND. NSW/=6) THEN
        CALL abort_physic('conf_phys', 'flag_aerosol=7 (MACv2SP) can only be activated with NSW=6',1)
     ENDIF
 
     ! BC internal mixture is only possible with RRTM & NSW=6 & flag_aerosol=6 or aerosol_couple
-    IF (flag_bc_internal_mixture .AND. NSW.NE.6) THEN 
+    IF (flag_bc_internal_mixture .AND. NSW/=6) THEN
        CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with NSW=6',1)
     ENDIF
-    IF (flag_bc_internal_mixture .AND. iflag_rrtm.NE.1) THEN 
+    IF (flag_bc_internal_mixture .AND. iflag_rrtm/=1) THEN
        CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with RRTM',1)
     ENDIF
-    IF (flag_bc_internal_mixture .AND. flag_aerosol.NE.6) THEN 
+    IF (flag_bc_internal_mixture .AND. flag_aerosol/=6) THEN
        CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with flag_aerosol=6',1)
     ENDIF
 
     ! test sur flag_volc_surfstrat
-    IF (flag_volc_surfstrat.LT.0.OR.flag_volc_surfstrat.GT.2) THEN
+    IF (flag_volc_surfstrat<0.OR.flag_volc_surfstrat>2) THEN
        CALL abort_physic('conf_phys', 'flag_volc_surfstrat can only be 0 1 or 2',1)
     ENDIF 
-    IF ((.NOT.ok_volcan.OR..NOT.ok_ade.OR..NOT.ok_aie).AND.flag_volc_surfstrat.GT.0) THEN
+    IF ((.NOT.ok_volcan.OR..NOT.ok_ade.OR..NOT.ok_aie).AND.flag_volc_surfstrat>0) THEN
        CALL abort_physic('conf_phys', 'ok_ade, ok_aie, ok_volcan need to be activated if flag_volc_surfstrat is 1 or 2',1)
     ENDIF 
@@ -2679,5 +2679,5 @@
 
     ! Test on chemistry cycle 
-    IF ((type_trac .ne. "inca" .AND. type_trac .ne. "inco") .AND. ( dms_cycle_cpl .OR. n2o_cycle_cpl)  ) THEN 
+    IF ((type_trac /= "inca" .AND. type_trac /= "inco") .AND. ( dms_cycle_cpl .OR. n2o_cycle_cpl)  ) THEN
        CALL abort_physic('conf_phys', 'dms_cycle_cpl or n2o_cycle_cpl has to be TRUE only with INCA coupling model',1)
     ENDIF
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/m_mrgrnk.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/m_mrgrnk.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/m_mrgrnk.F90	(revision 5081)
@@ -33,6 +33,4 @@
          IRNGT (1) = 1
          Return
-      Case Default
-         Continue
       End Select
 !
@@ -233,6 +231,4 @@
          IRNGT (1) = 1
          Return
-      Case Default
-         Continue
       End Select
 !
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/optics_lib.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/optics_lib.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/optics_lib.F90	(revision 5081)
@@ -38,15 +38,15 @@
   
 ! ----- INPUTS -----
-  real*8, intent(in) :: freq,tk
+  real(kind=8), intent(in) :: freq,tk
   
 ! ----- OUTPUTS -----
-  real*8, intent(out) :: n_r, n_i
+  real(kind=8), intent(out) :: n_r, n_i
 
 ! ----- INTERNAL -----    
-  real*8 ld,es,ei,a,ls,sg,tm1,cos1,sin1
-  real*8 e_r,e_i
-  real*8 pi
-  real*8 tc
-  complex*16 e_comp, sq
+  real(kind=8) ld,es,ei,a,ls,sg,tm1,cos1,sin1
+  real(kind=8) e_r,e_i
+  real(kind=8) pi
+  real(kind=8) tc
+  complex(kind=8) e_comp, sq
 
   tc = tk - 273.15
@@ -102,17 +102,17 @@
 
 ! ----- INPUTS -----
-  real*8, intent(in) :: freq, t
+  real(kind=8), intent(in) :: freq, t
   
 ! ----- OUTPUTS -----  
-  real*8, intent(out) :: n_r,n_i
+  real(kind=8), intent(out) :: n_r,n_i
 
 ! Parameters:
-  integer*2 :: i,lt1,lt2,nwl,nwlt
+  integer(kind=2) :: i,lt1,lt2,nwl,nwlt
   parameter(nwl=468,nwlt=62)
 
-  real*8 :: alam,cutice,pi,t1,t2,wlmax,wlmin, &
+  real(kind=8) :: alam,cutice,pi,t1,t2,wlmax,wlmin, &
             x,x1,x2,y,y1,y2,ylo,yhi,tk
 
-  real*8 :: &
+  real(kind=8) :: &
        tabim(nwl),tabimt(nwlt,4),tabre(nwl),tabret(nwlt,4),temref(4), &
        wl(nwl),wlt(nwlt)
@@ -540,10 +540,10 @@
     if(tk < temref(4)) tk=temref(4)
     do 11 i=2,4
-      if(tk.ge.temref(i)) go to 12
+      if(tk>=temref(i)) go to 12
     11 continue
     12 lt1=i
     lt2=i-1
     do 13 i=2,nwlt
-      if(alam.le.wlt(i)) go to 14
+      if(alam<=wlt(i)) go to 14
     13 continue
     14 x1=log(wlt(i-1))
@@ -586,64 +586,64 @@
       Subroutine MieInt(Dx, SCm, Inp, Dqv, Dqxt, Dqsc, Dbsc, Dg, Xs1, Xs2, DPh, Error)
 
-      Integer * 2  Imaxx
+      Integer (kind=2)  Imaxx
       Parameter (Imaxx = 12000)
-      Real * 4     RIMax          ! largest real part of refractive index
+      Real (kind=4)     RIMax          ! largest real part of refractive index
       Parameter (RIMax = 2.5)
-      Real * 4     IRIMax         ! largest imaginary part of refractive index
+      Real (kind=4)     IRIMax         ! largest imaginary part of refractive index
       Parameter (IRIMax = -2)
-      Integer * 2  Itermax
+      Integer (kind=2)  Itermax
       Parameter (Itermax = 12000 * 2.5)
                                 ! must be large enough to cope with the
                                 ! largest possible nmx = x * abs(scm) + 15
                                 ! or nmx =  Dx + 4.05*Dx**(1./3.) + 2.0
-      Integer * 2  Imaxnp
+      Integer (kind=2)  Imaxnp
       Parameter (Imaxnp = 10000)  ! Change this as required
 !     INPUT
-      Real * 8     Dx
-      Complex * 16  SCm
-      Integer * 4  Inp
-      Real * 8     Dqv(Inp)
+      Real (kind=8)     Dx
+      Complex (kind=8)  SCm
+      Integer (kind=4)  Inp
+      Real (kind=8)     Dqv(Inp)
 !     OUTPUT
-      Complex * 16  Xs1(InP)
-      Complex * 16  Xs2(InP)
-      Real * 8     Dqxt
-      Real * 8     Dqsc
-      Real * 8     Dg
-      Real * 8     Dbsc
-      Real * 8     DPh(InP)
-      Integer * 4  Error
+      Complex (kind=8)  Xs1(InP)
+      Complex (kind=8)  Xs2(InP)
+      Real (kind=8)     Dqxt
+      Real (kind=8)     Dqsc
+      Real (kind=8)     Dg
+      Real (kind=8)     Dbsc
+      Real (kind=8)     DPh(InP)
+      Integer (kind=4)  Error
 !     LOCAL
-      Integer * 2  I
-      Integer * 2  NStop
-      Integer * 2  NmX
-      Integer * 4  N    ! N*N > 32767 ie N > 181
-      Integer * 4  Inp2
-      Real * 8     Chi,Chi0,Chi1
-      Real * 8     APsi,APsi0,APsi1
-      Real * 8     Pi0(Imaxnp)
-      Real * 8     Pi1(Imaxnp)
-      Real * 8     Taun(Imaxnp)
-      Real * 8     Psi,Psi0,Psi1
-      Complex * 8  Ir
-      Complex * 16 Cm
-      Complex * 16 A,ANM1,APB
-      Complex * 16 B,BNM1,AMB
-      Complex * 16 D(Itermax)
-      Complex * 16 Sp(Imaxnp)
-      Complex * 16 Sm(Imaxnp)
-      Complex * 16 Xi,Xi0,Xi1
-      Complex * 16 Y
+      Integer (kind=2)  I
+      Integer (kind=2)  NStop
+      Integer (kind=2)  NmX
+      Integer (kind=4)  N    ! N*N > 32767 ie N > 181
+      Integer (kind=4)  Inp2
+      Real (kind=8)     Chi,Chi0,Chi1
+      Real (kind=8)     APsi,APsi0,APsi1
+      Real (kind=8)     Pi0(Imaxnp)
+      Real (kind=8)     Pi1(Imaxnp)
+      Real (kind=8)     Taun(Imaxnp)
+      Real (kind=8)     Psi,Psi0,Psi1
+      Complex (kind=4)  Ir
+      Complex (kind=8) Cm
+      Complex (kind=8) A,ANM1,APB
+      Complex (kind=8) B,BNM1,AMB
+      Complex (kind=8) D(Itermax)
+      Complex (kind=8) Sp(Imaxnp)
+      Complex (kind=8) Sm(Imaxnp)
+      Complex (kind=8) Xi,Xi0,Xi1
+      Complex (kind=8) Y
 !     ACCELERATOR VARIABLES
-      Integer * 2  Tnp1
-      Integer * 2  Tnm1
-      Real * 8     Dn
-      Real * 8     Rnx
-      Real * 8     S(Imaxnp)
-      Real * 8     T(Imaxnp)
-      Real * 8     Turbo
-      Real * 8     A2
-      Complex * 16 A1
+      Integer (kind=2)  Tnp1
+      Integer (kind=2)  Tnm1
+      Real (kind=8)     Dn
+      Real (kind=8)     Rnx
+      Real (kind=8)     S(Imaxnp)
+      Real (kind=8)     T(Imaxnp)
+      Real (kind=8)     Turbo
+      Real (kind=8)     A2
+      Complex (kind=8) A1
       
-      If ((Dx.Gt.Imaxx) .Or. (InP.Gt.ImaxNP)) Then
+      If ((Dx>Imaxx) .Or. (InP>ImaxNP)) Then
         Error = 1
         Return
@@ -652,11 +652,11 @@
       Ir = 1 / Cm
       Y =  Dx * Cm
-      If (Dx.Lt.0.02) Then
+      If (Dx<0.02) Then
          NStop = 2
       Else
-         If (Dx.Le.8.0) Then
+         If (Dx<=8.0) Then
             NStop = Dx + 4.00*Dx**(1./3.) + 2.0
          Else
-            If (Dx.Lt. 4200.0) Then
+            If (Dx< 4200.0) Then
                NStop = Dx + 4.05*Dx**(1./3.) + 2.0
             Else
@@ -666,5 +666,5 @@
       End If
       NmX = Max(Real(NStop),Real(Abs(Y))) + 15.
-      If (Nmx .gt. Itermax) then
+      If (Nmx > Itermax) then
           Error = 1
           Return
@@ -709,5 +709,5 @@
          Dqxt = Tnp1 *      Dble(A + B)          + Dqxt
          Dqsc = Tnp1 * (A*Conjg(A) + B*Conjg(B)) + Dqsc
-         If (N.Gt.1) then
+         If (N>1) then
          Dg = Dg + (dN*dN - 1) * Dble(ANM1*Conjg(A) + BNM1 * Conjg(B)) / dN + TNM1 * Dble(ANM1*Conjg(BNM1)) / (dN*dN - dN)
          End If
@@ -717,5 +717,5 @@
          AMB = A2 * (A - B)
          Do I = 1,Inp2
-            If (I.GT.Inp) Then
+            If (I>Inp) Then
                S(I) = -Pi1(I)
             Else
@@ -736,5 +736,5 @@
          Xi1 = Dcmplx(APsi1,Chi1)
       End Do
-      If (Dg .GT.0) Dg = 2 * Dg / Dqsc
+      If (Dg >0) Dg = 2 * Dg / Dqsc
       Dqsc =  2 * Dqsc / Dx**2
       Dqxt =  2 * Dqxt / Dx**2
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/cosp2/mrgrnk.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/cosp2/mrgrnk.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/cosp2/mrgrnk.F90	(revision 5081)
@@ -68,6 +68,4 @@
        IRNGT (1) = 1
        Return
-    Case Default
-       Continue
     End Select
     !
@@ -268,6 +266,4 @@
        IRNGT (1) = 1
        Return
-    Case Default
-       Continue
     End Select
     !
@@ -467,6 +463,4 @@
        IRNGT (1) = 1
        Return
-    Case Default
-       Continue
     End Select
     !
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_optics.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_optics.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_optics.F90	(revision 5081)
@@ -72,8 +72,8 @@
     varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
     do j=1,dim2
-       where(flag(:,j,:) .eq. 1)
+       where(flag(:,j,:) == 1)
           varOUT(:,j,:) = varIN2
        endwhere
-       where(flag(:,j,:) .eq. 2)
+       where(flag(:,j,:) == 2)
           varOUT(:,j,:) = varIN1
        endwhere
@@ -96,8 +96,8 @@
     
     varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
-   where(flag(:,:,:) .eq. 1)
+   where(flag(:,:,:) == 1)
        varOUT(:,:,:) = varIN2
     endwhere
-    where(flag(:,:,:) .eq. 2)
+    where(flag(:,:,:) == 2)
        varOUT(:,:,:) = varIN1
     endwhere
@@ -295,9 +295,9 @@
 
     ! Which LIDAR frequency are we using?
-    if (lidar_freq .eq. 355) then
+    if (lidar_freq == 355) then
        Cmol   = Cmol_355nm
        rdiffm = rdiffm_355nm
     endif
-    if (lidar_freq .eq. 532) then
+    if (lidar_freq == 532) then
        Cmol   = Cmol_532nm
        rdiffm = rdiffm_532nm
@@ -336,9 +336,9 @@
     
     ! LS and CONV Ice water coefficients
-    if (ice_type .eq. 0) then
+    if (ice_type == 0) then
        polpart(INDX_LSICE,1:5) = polpartLSICE0
        polpart(INDX_CVICE,1:5) = polpartCVICE0
     endif
-    if (ice_type .eq. 1) then
+    if (ice_type == 1) then
        polpart(INDX_LSICE,1:5) = polpartLSICE1
        polpart(INDX_CVICE,1:5) = polpartCVICE1
@@ -393,5 +393,5 @@
     ! Polynomials kp_lidar derived from Mie theory
     do i = 1, npart
-       where (rad_part(1:npoints,1:nlev,i) .gt. 0.0)
+       where (rad_part(1:npoints,1:nlev,i) > 0.0)
           kp_part(1:npoints,1:nlev,i) = &
                polpart(i,1)*(rad_part(1:npoints,1:nlev,i)*1e6)**4 &
@@ -426,5 +426,5 @@
        ! Alpha of particles in each subcolumn:
        do i = 1, npart
-          where (rad_part(1:npoints,1:nlev,i) .gt. 0.0)
+          where (rad_part(1:npoints,1:nlev,i) > 0.0)
              alpha_part(1:npoints,1:nlev,i) = 3._wp/4._wp * Qscat &
                   * rhoair(1:npoints,1:nlev) * qpart(1:npoints,1:nlev,i) &
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/mrgrnk.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/mrgrnk.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/mrgrnk.F90	(revision 5081)
@@ -68,6 +68,4 @@
        IRNGT (1) = 1
        Return
-    Case Default
-       Continue
     End Select
     !
@@ -268,6 +266,4 @@
        IRNGT (1) = 1
        Return
-    Case Default
-       Continue
     End Select
     !
@@ -467,6 +463,4 @@
        IRNGT (1) = 1
        Return
-    Case Default
-       Continue
     End Select
     !
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/optics_lib.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/optics_lib.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/optics_lib.F90	(revision 5081)
@@ -558,10 +558,10 @@
        if(tk < temref(4)) tk=temref(4)
        do i=2,4
-          if(tk.ge.temref(i)) go to 12
+          if(tk>=temref(i)) go to 12
        enddo
 12     lt1 = i
        lt2 = i-1
        do i=2,nwlt
-          if(alam.le.wlt(i)) go to 14
+          if(alam<=wlt(i)) go to 14
        enddo
 14     x1  = log(wlt(i-1))
@@ -652,5 +652,5 @@
     Complex(wp) :: A1
     
-    If ((Dx.Gt.Imaxx) .Or. (InP.Gt.ImaxNP)) Then
+    If ((Dx>Imaxx) .Or. (InP>ImaxNP)) Then
        Error = 1
        Return
@@ -659,11 +659,11 @@
     Ir = 1 / Cm
     Y =  Dx * Cm
-    If (Dx.Lt.0.02) Then
+    If (Dx<0.02) Then
        NStop = 2
     Else
-       If (Dx.Le.8.0) Then
+       If (Dx<=8.0) Then
           NStop = Dx + 4.00*Dx**(1./3.) + 2.0
        Else
-          If (Dx.Lt. 4200.0) Then
+          If (Dx< 4200.0) Then
              NStop = Dx + 4.05*Dx**(1./3.) + 2.0
           Else
@@ -673,5 +673,5 @@
     End If
     NmX = Max(Real(NStop),Real(Abs(Y))) + 15.
-    If (Nmx .gt. Itermax) then
+    If (Nmx > Itermax) then
        Error = 1
        Return
@@ -726,5 +726,5 @@
 !ds       Dqxt = Tnp1 *      Dble(A + B)          + Dqxt
        Dqsc = Tnp1 * (A*Conjg(A) + B*Conjg(B)) + Dqsc
-       If (N.Gt.1) then
+       If (N>1) then
           Dg = Dg + (dN*dN - 1) * (ANM1*Conjg(A) + BNM1 * Conjg(B)) / dN + TNM1 *(ANM1*Conjg(BNM1)) / (dN*dN - dN)
 !ds          Dg = Dg + (dN*dN - 1) * Dble(ANM1*Conjg(A) + BNM1 * Conjg(B)) / dN + TNM1 * Dble(ANM1*Conjg(BNM1)) / (dN*dN - dN)
@@ -735,5 +735,5 @@
        AMB = A2 * (A - B)
        Do I = 1,Inp2
-          If (I.GT.Inp) Then
+          If (I>Inp) Then
              S(I) = -Pi1(I)
           Else
@@ -756,5 +756,5 @@
     End Do
 
-    If (Dg .GT.0) Dg = 2 * Dg / Dqsc
+    If (Dg >0) Dg = 2 * Dg / Dqsc
     Dqsc =  2 * Dqsc / Dx**2
     Dqxt =  2 * Dqxt / Dx**2
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam.F90	(revision 5081)
@@ -179,5 +179,5 @@
           
           ! Attenuation due to gaseous absorption between radar and volume
-          if ((rcfg%use_gas_abs == 1) .or. (rcfg%use_gas_abs == 2 .and. pr .eq. 1)) then
+          if ((rcfg%use_gas_abs == 1) .or. (rcfg%use_gas_abs == 2 .and. pr == 1)) then
              if (d_gate==1) then
                 if (k>1) then
@@ -272,5 +272,5 @@
 
     ! Which platforms to create diagnostics for?
-    if (platform .eq. 'cloudsat') lcloudsat=.true.
+    if (platform == 'cloudsat') lcloudsat=.true.
 
     ! Create Cloudsat diagnostics.
@@ -289,5 +289,5 @@
              enddo
           enddo
-          where(cfad_ze .ne. R_UNDEF) cfad_ze = cfad_ze/Ncolumns
+          where(cfad_ze /= R_UNDEF) cfad_ze = cfad_ze/Ncolumns
 
           ! Compute cloudsat near-surface precipitation diagnostics
@@ -306,5 +306,5 @@
              enddo
           enddo
-          where(cfad_ze .ne. R_UNDEF) cfad_ze = cfad_ze/Ncolumns
+          where(cfad_ze /= R_UNDEF) cfad_ze = cfad_ze/Ncolumns
        endif
     endif
@@ -402,6 +402,6 @@
        do pr=1,Ncolumns
           ! 1) Compute the PIA in all profiles containing hydrometeors
-          if ( (Ze_non_out(i,pr,cloudsat_preclvl).gt.-100) .and. (Ze_out(i,pr,cloudsat_preclvl).gt.-100) ) then
-             if ( (Ze_non_out(i,pr,cloudsat_preclvl).lt.100) .and. (Ze_out(i,pr,cloudsat_preclvl).lt.100) ) then
+          if ( (Ze_non_out(i,pr,cloudsat_preclvl)>-100) .and. (Ze_out(i,pr,cloudsat_preclvl)>-100) ) then
+             if ( (Ze_non_out(i,pr,cloudsat_preclvl)<100) .and. (Ze_out(i,pr,cloudsat_preclvl)<100) ) then
                 cloudsat_precip_pia(i,pr) = Ze_non_out(i,pr,cloudsat_preclvl) - Ze_out(i,pr,cloudsat_preclvl)
              endif
@@ -412,13 +412,13 @@
           ! 2a) Oceanic points.
           ! ################################################################################
-          if (land(i) .eq. 0) then
+          if (land(i) == 0) then
 !             print*, 'aaa i, pr, fracPrecipIce(i,pr) : ', i, pr, fracPrecipIce(i,pr) !Artem
              ! Snow
-             if(fracPrecipIce(i,pr).gt.0.9) then
-                if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(2)) then
+             if(fracPrecipIce(i,pr)>0.9) then
+                if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(2)) then
                    cloudsat_pflag(i,pr) = pClass_Snow2                   ! TSL: Snow certain
                 endif
-                if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(4).and. &
-                     Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(2)) then
+                if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(4).and. &
+                     Ze_non_out(i,pr,cloudsat_preclvl)<=Zenonbinval(2)) then
                    cloudsat_pflag(i,pr) = pClass_Snow1                   ! TSL: Snow possible
                 endif
@@ -426,10 +426,10 @@
              
              ! Mixed
-             if(fracPrecipIce(i,pr).gt.0.1.and.fracPrecipIce(i,pr).le.0.9) then
-                if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(2)) then
+             if(fracPrecipIce(i,pr)>0.1.and.fracPrecipIce(i,pr)<=0.9) then
+                if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(2)) then
                    cloudsat_pflag(i,pr) = pClass_Mixed2                  ! TSL: Mixed certain
                 endif
-                if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(4).and. &
-                     Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(2)) then
+                if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(4).and. &
+                     Ze_non_out(i,pr,cloudsat_preclvl)<=Zenonbinval(2)) then
                    cloudsat_pflag(i,pr) = pClass_Mixed1                  ! TSL: Mixed possible
                 endif
@@ -437,17 +437,17 @@
              
              ! Rain
-             if(fracPrecipIce(i,pr).le.0.1) then
-                if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(1)) then
+             if(fracPrecipIce(i,pr)<=0.1) then
+                if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(1)) then
                    cloudsat_pflag(i,pr) = pClass_Rain3                   ! TSL: Rain certain
                 endif
-                if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(3).and. &
-                     Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(1)) then
+                if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(3).and. &
+                     Ze_non_out(i,pr,cloudsat_preclvl)<=Zenonbinval(1)) then
                    cloudsat_pflag(i,pr) = pClass_Rain2                   ! TSL: Rain probable
                 endif
-                if(Ze_non_out(i,pr,cloudsat_preclvl).gt.Zenonbinval(4).and. &
-                     Ze_non_out(i,pr,cloudsat_preclvl).le.Zenonbinval(3)) then
+                if(Ze_non_out(i,pr,cloudsat_preclvl)>Zenonbinval(4).and. &
+                     Ze_non_out(i,pr,cloudsat_preclvl)<=Zenonbinval(3)) then
                    cloudsat_pflag(i,pr) = pClass_Rain1                   ! TSL: Rain possible
                 endif
-                if(cloudsat_precip_pia(i,pr).gt.40) then
+                if(cloudsat_precip_pia(i,pr)>40) then
                    cloudsat_pflag(i,pr) = pClass_Rain4                   ! TSL: Heavy Rain
                 endif
@@ -455,5 +455,5 @@
              
              ! No precipitation
-             if(Ze_non_out(i,pr,cloudsat_preclvl).le.-15) then
+             if(Ze_non_out(i,pr,cloudsat_preclvl)<=-15) then
                 cloudsat_pflag(i,pr) = pClass_noPrecip                   ! TSL: Not Raining
              endif
@@ -463,15 +463,15 @@
           ! 2b) Land points.
           ! ################################################################################
-          if (land(i) .eq. 1) then
+          if (land(i) == 1) then
              ! Find Zmax, the maximum reflectivity value in the attenuated profile (Ze_out);
              Zmax=maxval(Ze_out(i,pr,:))
 
              ! Snow (T<273)
-             if(t2m(i) .lt. 273._wp) then
-                if(Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(5)) then
+             if(t2m(i) < 273._wp) then
+                if(Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(5)) then
                    cloudsat_pflag(i,pr) = pClass_Snow2                      ! JEK: Snow certain
                 endif
-                if(Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6) .and. &
-                     Ze_out(i,pr,cloudsat_preclvl).le.Zbinvallnd(5)) then
+                if(Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(6) .and. &
+                     Ze_out(i,pr,cloudsat_preclvl)<=Zbinvallnd(5)) then
                    cloudsat_pflag(i,pr) = pClass_Snow1                      ! JEK: Snow possible
                 endif
@@ -479,12 +479,12 @@
              
              ! Mized phase (273<T<275)
-             if(t2m(i) .ge. 273._wp .and. t2m(i) .le. 275._wp) then
-                if ((Zmax .gt. Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr).gt.30) .or. &
-                     (Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(4))) then
+             if(t2m(i) >= 273._wp .and. t2m(i) <= 275._wp) then
+                if ((Zmax > Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr)>30) .or. &
+                     (Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(4))) then
                    cloudsat_pflag(i,pr) = pClass_Mixed2                     ! JEK: Mixed certain
                 endif
-                if ((Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6)  .and. &
-                     Ze_out(i,pr,cloudsat_preclvl) .le. Zbinvallnd(4)) .and. &
-                     (Zmax .gt. Zbinvallnd(5)) ) then
+                if ((Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(6)  .and. &
+                     Ze_out(i,pr,cloudsat_preclvl) <= Zbinvallnd(4)) .and. &
+                     (Zmax > Zbinvallnd(5)) ) then
                    cloudsat_pflag(i,pr) = pClass_Mixed1                     ! JEK: Mixed possible
                 endif
@@ -492,18 +492,18 @@
 
              ! Rain (T>275)
-             if(t2m(i) .gt. 275) then
-                if ((Zmax .gt. Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr).gt.30) .or. &
-                     (Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(2))) then
+             if(t2m(i) > 275) then
+                if ((Zmax > Zbinvallnd(1) .and. cloudsat_precip_pia(i,pr)>30) .or. &
+                     (Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(2))) then
                    cloudsat_pflag(i,pr) = pClass_Rain3                      ! JEK: Rain certain
                 endif
-                if((Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6)) .and. &
-                     (Zmax .gt. Zbinvallnd(3))) then
+                if((Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(6)) .and. &
+                     (Zmax > Zbinvallnd(3))) then
                    cloudsat_pflag(i,pr) = pClass_Rain2                      ! JEK: Rain probable
                 endif
-                if((Ze_out(i,pr,cloudsat_preclvl) .gt. Zbinvallnd(6)) .and. &
-                     (Zmax.lt.Zbinvallnd(3))) then
+                if((Ze_out(i,pr,cloudsat_preclvl) > Zbinvallnd(6)) .and. &
+                     (Zmax<Zbinvallnd(3))) then
                    cloudsat_pflag(i,pr) = pClass_Rain1                      ! JEK: Rain possible
                 endif
-                if(cloudsat_precip_pia(i,pr).gt.40) then
+                if(cloudsat_precip_pia(i,pr)>40) then
                    cloudsat_pflag(i,pr) = pClass_Rain4                      ! JEK: Heavy Rain
                 endif
@@ -511,5 +511,5 @@
              
              ! No precipitation
-             if(Ze_out(i,pr,cloudsat_preclvl).le.-15) then
+             if(Ze_out(i,pr,cloudsat_preclvl)<=-15) then
                 cloudsat_pflag(i,pr) =  pClass_noPrecip                     ! JEK: Not Precipitating
              endif         
@@ -526,6 +526,6 @@
        ! Gridmean precipitation fraction for each precipitation type
        do k=1,nCloudsatPrecipClass
-          if (any(cloudsat_pflag(i,:) .eq. k-1)) then
-             cloudsat_precip_cover(i,k) = count(cloudsat_pflag(i,:) .eq. k-1)
+          if (any(cloudsat_pflag(i,:) == k-1)) then
+             cloudsat_precip_cover(i,k) = count(cloudsat_pflag(i,:) == k-1)
           endif
        enddo
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam_optics.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam_optics.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/quickbeam_optics.F90	(revision 5081)
@@ -172,5 +172,5 @@
                 
                 ! Compute effective radius from number concentration and distribution parameters
-                if (Re_internal .eq. 0) then
+                if (Re_internal == 0) then
                    call calc_Re(hm_matrix(pr,k,tp),Np_matrix(pr,k,tp),rho_a, &
                         sd%dtype(tp),sd%apm(tp),sd%bpm(tp),sd%rho(tp),sd%p1(tp),sd%p2(tp),sd%p3(tp),Re)
@@ -187,5 +187,5 @@
                 ! Index into particle size dimension of scaling tables 
                 iRe_type=1
-                if(Re.gt.0) then
+                if(Re>0) then
                    ! Determine index in to scale LUT
                    ! Distance between Re points (defined by "base" and "step") for
@@ -197,10 +197,10 @@
                    base = rcfg%base_list(n+1)
                    iRe_type=Re/step
-                   if (iRe_type.lt.1) iRe_type=1
+                   if (iRe_type<1) iRe_type=1
                    Re=step*(iRe_type+0.5_wp)    ! set value of Re to closest value allowed in LUT.
                    iRe_type=iRe_type+base-int(n*Re_BIN_LENGTH/step)
                    
                    ! Make sure iRe_type is within bounds
-                   if (iRe_type.ge.nRe_types) then
+                   if (iRe_type>=nRe_types) then
                       !write(*,*) 'Warning: size of Re exceed value permitted ', &
                       !            'in Look-Up Table (LUT).  Will calculate. '
@@ -405,5 +405,5 @@
     ! Exponential is same as modified gamma with vu =1
     ! if Np is specified then we will just treat as modified gamma
-    if(dtype .eq. 2 .and. Np .gt. 0) then
+    if(dtype == 2 .and. Np > 0) then
        local_dtype = 1
        local_p3    = 1
@@ -441,9 +441,9 @@
        endif
        
-       if( Np.eq.0 .and. p2+1 > 1E-8) then     ! use default value for MEAN diameter as first default  
+       if( Np==0 .and. p2+1 > 1E-8) then     ! use default value for MEAN diameter as first default
           dm = p2             ! by definition, should have units of microns
           D0 = gamma(vu)/gamma(vu+1)*dm
        else   ! use value of Np
-          if(Np.eq.0) then
+          if(Np==0) then
              if( abs(p1+1) > 1E-8 ) then  !   use default number concentration   
                 local_Np = p1 ! total number concentration / pa --- units kg^-1
@@ -525,5 +525,5 @@
        
        ! get rg ... 
-       if( Np.eq.0 .and. (abs(p2+1) > 1E-8) ) then ! use default value of rg
+       if( Np==0 .and. (abs(p2+1) > 1E-8) ) then ! use default value of rg
           rg = p2     
        else
@@ -826,5 +826,5 @@
           log_sigma_g = p3
           tmp2 = (bpm*log_sigma_g)*(bpm*log_sigma_g)
-          if(Re.le.0) then 
+          if(Re<=0) then
              rg = p2
           else
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/scops.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/scops.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/scops.F90	(revision 5081)
@@ -75,5 +75,5 @@
 
     ! Test for valid input overlap assumption
-    if (overlap .ne. 1 .and. overlap .ne. 2 .and. overlap .ne. 3) then
+    if (overlap /= 1 .and. overlap /= 2 .and. overlap /= 3) then
        overlap=default_overlap
        call errorMessage('ERROR(scops): Invalid overlap assumption provided. Using default overlap assumption (max/ran)')
@@ -92,5 +92,5 @@
     tca(1:npoints,1:nlev) = cc(1:npoints,1:nlev)
     
-    if (ncolprint.ne.0) then
+    if (ncolprint/=0) then
        write (6,'(a)') 'frac_out_pp_rev:'
        do j=1,npoints,1000
@@ -102,5 +102,5 @@
        write (6,'(I3)') ncol
     endif
-    if (ncolprint.ne.0) then
+    if (ncolprint/=0) then
        write (6,'(a)') 'last_frac_pp:'
        do j=1,npoints,1000
@@ -122,7 +122,7 @@
        
        ! Initialise threshold
-       IF (ilev.eq.1) then
+       IF (ilev==1) then
           ! If max overlap 
-          IF (overlap.eq.1) then
+          IF (overlap==1) then
              ! Select pixels spread evenly across the gridbox
              threshold(1:npoints,1:ncol)=boxpos(1:npoints,1:ncol)
@@ -137,5 +137,5 @@
              enddo
           ENDIF
-          IF (ncolprint.ne.0) then
+          IF (ncolprint/=0) then
              write (6,'(a)') 'threshold_nsf2:'
              do j=1,npoints,1000
@@ -147,5 +147,5 @@
        ENDIF
        
-       IF (ncolprint.ne.0) then
+       IF (ncolprint/=0) then
           write (6,'(a)') 'ilev:'
           write (6,'(I2)') ilev
@@ -157,5 +157,5 @@
           !maxocc(1:npoints,ibox) = merge(1,0, conv(1:npoints,ilev) .gt. boxpos(1:npoints,ibox))
           do j=1,npoints
-             if (boxpos(j,ibox).le.conv(j,ilev)) then
+             if (boxpos(j,ibox)<=conv(j,ilev)) then
                 maxocc(j,ibox) = 1
              else
@@ -165,5 +165,5 @@
           
           ! Max overlap
-          if (overlap.eq.1) then 
+          if (overlap==1) then
              threshold_min(1:npoints,ibox) = conv(1:npoints,ilev)
              maxosc(1:npoints,ibox)        = 1               
@@ -171,10 +171,10 @@
           
           ! Random overlap
-          if (overlap.eq.2) then 
+          if (overlap==2) then
              threshold_min(1:npoints,ibox) = conv(1:npoints,ilev)
              maxosc(1:npoints,ibox)        = 0
           endif
           ! Max/Random overlap
-          if (overlap.eq.3) then 
+          if (overlap==3) then
              ! DS2014 START: The bounds on tca are not valid when ilev=1.
              !threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)))
@@ -182,14 +182,14 @@
              !     min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)) .and. &
              !     (threshold(1:npoints,ibox).gt.conv(1:npoints,ilev)))
-             if (ilev .ne. 1) then
+             if (ilev /= 1) then
                 threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)))
-                maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) .lt. &
+                maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) < &
                      min(tca(1:npoints,ilev-1),tca(1:npoints,ilev)) .and. &
-                     (threshold(1:npoints,ibox).gt.conv(1:npoints,ilev)))
+                     (threshold(1:npoints,ibox)>conv(1:npoints,ilev)))
              else
                 threshold_min(1:npoints,ibox) = max(conv(1:npoints,ilev),min(0._wp,tca(1:npoints,ilev)))
-                maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) .lt. &
+                maxosc(1:npoints,ibox) = merge(1,0,threshold(1:npoints,ibox) < &
                      min(0._wp,tca(1:npoints,ilev)) .and. &
-                     (threshold(1:npoints,ibox).gt.conv(1:npoints,ilev)))
+                     (threshold(1:npoints,ibox)>conv(1:npoints,ilev)))
              endif
           endif
@@ -205,13 +205,13 @@
           
           ! Fill frac_out with 1's where tca is greater than the threshold
-          frac_out(1:npoints,ibox,ilev) = merge(1,0,tca(1:npoints,ilev).gt.threshold(1:npoints,ibox))
+          frac_out(1:npoints,ibox,ilev) = merge(1,0,tca(1:npoints,ilev)>threshold(1:npoints,ibox))
           
           ! Code to partition boxes into startiform and convective parts goes here
-          where(threshold(1:npoints,ibox).le.conv(1:npoints,ilev) .and. conv(1:npoints,ilev).gt.0.) frac_out(1:npoints,ibox,ilev)=2
+          where(threshold(1:npoints,ibox)<=conv(1:npoints,ilev) .and. conv(1:npoints,ilev)>0.) frac_out(1:npoints,ibox,ilev)=2
        ENDDO ! ibox
        
        
        ! Set last_frac to tca at this level, so as to be tca from last level next time round
-       if (ncolprint.ne.0) then
+       if (ncolprint/=0) then
           do j=1,npoints ,1000
              write(6,'(a10)') 'j='
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/cv30param.h
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/cv30param.h	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/cv30param.h	(revision 5081)
@@ -20,11 +20,11 @@
       real betad
 
-      COMMON /cv30param/  noff, minorig, nl, nlp, nlm &
-                      ,  sigd, spfac &
+      COMMON /cv30param/  sigd, spfac &
 !IM cf. FH : pour compatibilite avec conema3 TEMPORAIRE  :                ,pbcrit, ptcrit, epmax
                       ,pbcrit, ptcrit &
                       ,omtrain &
                       ,dtovsh, dpbase, dttrig &
-                      ,dtcrit, tau, beta, alpha, delta, betad
+                      ,dtcrit, tau, beta, alpha, delta, betad &
+                      ,noff, minorig, nl, nlp, nlm
 
 !$OMP THREADPRIVATE(/cv30param/)
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/cv3param.h
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/cv3param.h	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/cv3param.h	(revision 5081)
@@ -38,6 +38,7 @@
                       ,delta, betad  &
                       ,ejectliq, ejectice &
+                      ,flag_wb &
                       ,flag_epKEorig &
-                      ,flag_wb, cv_flag_feed &
+                      ,cv_flag_feed &
                       ,noff, minorig, nl, nlp, nlm
 !$OMP THREADPRIVATE(/cv3param/)
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/cva_driver.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/cva_driver.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/cva_driver.F90	(revision 5081)
@@ -858,5 +858,5 @@
 !   is assumed useless.
 !
-  compress = ncum .lt. len*comp_threshold
+  compress = ncum < len*comp_threshold
 !
   IF (.not. compress) THEN
@@ -1100,5 +1100,5 @@
         Print *, 'cva_driver after cv3_unsat:mp , water, ice, evap, fondue '
         DO k = 1,nd
-        write (6, '(i4,5(1x,e13.6))'), &
+        write (6, '(i4,5(1x,e13.6))') &
           k, mp(igout,k), water(igout,k), ice(igout,k), &
            evap(igout,k), fondue(igout,k)
@@ -1106,5 +1106,5 @@
         Print *, 'cva_driver after cv3_unsat: wdtrainA, wdtrainS, wdtrainM '     !!jygprl
         DO k = 1,nd
-        write (6, '(i4,3(1x,e13.6))'), &
+        write (6, '(i4,3(1x,e13.6))') &
            k, wdtrainA(igout,k), wdtrainS(igout,k), wdtrainM(igout,k)            !!jygprl
         ENDDO
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/cvparam.h
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/cvparam.h	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/cvparam.h	(revision 5081)
@@ -21,9 +21,9 @@
       real delta
 
-      COMMON /cvparam/ noff, minorig, nl, nlp, nlm &
-                      ,elcrit, tlcrit &
+      COMMON /cvparam/ elcrit, tlcrit &
                       ,entp, sigs, sigd &
                       ,omtrain, omtsnow, coeffr, coeffs &
-                      ,dtmax, cu, betad, alpha, damp, delta
+                      ,dtmax, cu, betad, alpha, damp, delta &
+                      ,noff, minorig, nl, nlp, nlm
 
 !$OMP THREADPRIVATE(/cvparam/)
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/flux_arp.h
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/flux_arp.h	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/flux_arp.h	(revision 5081)
@@ -15,5 +15,5 @@
       real :: tg
 
-      common /flux_arp/fsens,flat,ust,tg,ok_flux_surf,ok_prescr_ust,ok_prescr_beta,betaevap,ok_forc_tsurf
+      common /flux_arp/fsens,flat,betaevap,ust,tg,ok_flux_surf,ok_prescr_ust,ok_prescr_beta,ok_forc_tsurf
 
 !$OMP THREADPRIVATE(/flux_arp/)
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_ratqs_ini.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_ratqs_ini.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_ratqs_ini.F90	(revision 5081)
@@ -3,20 +3,18 @@
 IMPLICIT NONE
 
-save
+INTEGER :: lunout
 
-integer :: lunout
+INTEGER, PROTECTED :: nbsrf,is_lic,is_ter
+REAL, PROTECTED :: RG,RV,RD,RCPD,RLSTT,RLVTT,RTT
+REAL, PROTECTED :: a_ratqs_cv
+REAL, PROTECTED :: tau_var
+REAL, PROTECTED :: fac_tau
+REAL, PROTECTED :: tau_cumul
+REAL, PROTECTED :: a_ratqs_wake
+INTEGER, PROTECTED :: dqimpl
 
-INTEGER, SAVE, PROTECTED :: nbsrf,is_lic,is_ter
-REAL, SAVE, PROTECTED :: RG,RV,RD,RCPD,RLSTT,RLVTT,RTT
-REAL, SAVE, PROTECTED :: a_ratqs_cv
-REAL, SAVE, PROTECTED :: tau_var
-REAL, SAVE, PROTECTED :: fac_tau
-REAL, SAVE, PROTECTED :: tau_cumul
-REAL, SAVE, PROTECTED :: a_ratqs_wake
-INTEGER, SAVE, PROTECTED :: dqimpl
-
-real, allocatable, SAVE :: povariance(:,:)
+REAL, ALLOCATABLE :: povariance(:,:)
 !$OMP THREADPRIVATE(povariance)
-real, allocatable, SAVE :: var_conv(:,:)
+REAL, ALLOCATABLE :: var_conv(:,:)
 !$OMP THREADPRIVATE(var_conv)
 
@@ -78,5 +76,5 @@
 !--------------------------------------------------------
 
-if (klon.eq.1) then
+if (klon==1) then
     do k=1,klev
       do i=1,klon 
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_ini.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_ini.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_ini.F90	(revision 5081)
@@ -2,7 +2,4 @@
 
 IMPLICIT NONE
-
-save
-
 
 integer, protected :: dvdq=1,dqimpl=-1,prt_level=0,lunout
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/rrtm/rrtm_taumol1.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/rrtm/rrtm_taumol1.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/rrtm/rrtm_taumol1.F90	(revision 5081)
@@ -190,5 +190,5 @@
 IF (LHOOK) CALL DR_HOOK('RRTM_TAUMOL1',0,ZHOOK_HANDLE)
 !--ajout OB 
-IF (K_LAYTROP.GT.100) THEN 
+IF (K_LAYTROP>100) THEN
 PRINT *,'ATTENTION KLAY_TROP > 100 PROBLEME ARRAY DANS RRTM ON ARRETE'
 STOP
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/tlift.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/tlift.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/tlift.F90	(revision 5081)
@@ -213,6 +213,4 @@
         tvp(i) = tpk(i)*(1.+qsat_new/eps-rr(nk))
         ! jyg2
-      ELSE
-        CONTINUE
       END IF
 
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/yamada_ini_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/yamada_ini_mod.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/yamada_ini_mod.F90	(revision 5081)
@@ -5,6 +5,4 @@
 
 implicit none
-
-save
 
   LOGICAL :: new_yamada4
Index: /LMDZ6/branches/Amaury_dev/libf/phylmdiso/add_phys_tend_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmdiso/add_phys_tend_mod.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmdiso/add_phys_tend_mod.F90	(revision 5081)
@@ -297,5 +297,5 @@
   END DO
 
-  if (fl_ebil .GT. 0) then
+  if (fl_ebil > 0) then
     ! ------------------------------------------------
     ! Compute vertical sum for each atmospheric column
@@ -408,8 +408,8 @@
 !=====================================================================================
 
-IF (jbad .GT. 0) THEN
+IF (jbad > 0) THEN
       DO j = 1, jbad
          i=jadrs(j)
-         if(prt_level.ge.debug_level) THEN
+         if(prt_level>=debug_level) THEN
           print*,'PLANTAGE POUR LE POINT i lon lat =',&
                  i,longitude_deg(i),latitude_deg(i),text
@@ -426,9 +426,9 @@
 ! Impression, warning et correction en cas de probleme moins important
 !=====================================================================================
-IF (jqbad .GT. 0) THEN
+IF (jqbad > 0) THEN
       done(:) = .false.                         !jyg
       DO j = 1, jqbad
         i=jqadrs(j)
-          if(prt_level.ge.debug_level) THEN
+          if(prt_level>=debug_level) THEN
            print*,'WARNING  : EAU POUR LE POINT i lon lat =',&
                   i,longitude_deg(i),latitude_deg(i),text
@@ -465,5 +465,5 @@
               enddo
 #endif
-              if(prt_level.ge.debug_level) THEN
+              if(prt_level>=debug_level) THEN
                print*,' cas q_seri<1.e-15 i k zq_int zqp_int zq_int/zqp_int :', &
                                     i, kqadrs(j), zq_int, zqp_int, zq_int/zqp_int
@@ -487,7 +487,7 @@
             DO k = 1, klev
               zq=q_seri(i,k)+zdq(i,k)
-              if (zq.lt.1.e-15) then
-                 if (q_seri(i,k).lt.1.e-15) then
-                  if(prt_level.ge.debug_level) THEN
+              if (zq<1.e-15) then
+                 if (q_seri(i,k)<1.e-15) then
+                  if(prt_level>=debug_level) THEN
                    print*,' cas q_seri<1.e-15 i k q_seri zq zdq :',i,k,q_seri(i,k),zq,zdq(i,k)
                   endif
@@ -555,9 +555,9 @@
          ENDDO
       ENDDO
-IF (jbad .GT. 0) THEN
+IF (jbad > 0) THEN
       DO j = 1, jbad
          i=jadrs(j)
          k=kadrs(j)
-         if(prt_level.ge.debug_level) THEN
+         if(prt_level>=debug_level) THEN
           print*,'PLANTAGE2 POUR LE POINT i itap lon lat txt jbad zdt t',&
                  i,itap,longitude_deg(i),latitude_deg(i),text,jbad, &
@@ -573,9 +573,9 @@
 ENDIF 
 !
-IF (jqbad .GT. 0) THEN
+IF (jqbad > 0) THEN
       DO j = 1, jqbad
          i=jqadrs(j)
          k=kqadrs(j)
-         if(prt_level.ge.debug_level) THEN
+         if(prt_level>=debug_level) THEN
           print*,'WARNING  : EAU2 POUR LE POINT i itap lon lat txt jqbad zdq q zdql ql',&
                  i,itap,longitude_deg(i),latitude_deg(i),text,jqbad,&
@@ -613,5 +613,5 @@
 !======================================================================
 
-  if (fl_ebil .GT. 0) then
+  if (fl_ebil > 0) then
   
     ! ------------------------------------------------
@@ -764,5 +764,5 @@
   END DO
 
-  if (fl_ebil .GT. 0) then
+  if (fl_ebil > 0) then
     ! ------------------------------------------------
     ! Compute vertical sum for each atmospheric column
@@ -795,5 +795,5 @@
 !======================================================================
 
-  if (fl_ebil .GT. 0) then
+  if (fl_ebil > 0) then
   
     ! ------------------------------------------------
@@ -944,5 +944,5 @@
 
 !!print *,'prt_level:',prt_level,' fl_ebil:',fl_ebil,' fl_cor_ebil:',fl_cor_ebil
-if ( (fl_ebil .GT. 0) .and. (klon .EQ. 1)) then
+if ( (fl_ebil > 0) .and. (klon == 1)) then
 
   bilq_bnd = 0.
@@ -976,10 +976,10 @@
   bilh_error = d_h_col(1) - bilh_bnd
 ! are the errors too large?
-  if ( abs(bilq_error) .gt. bilq_seuil) bilq_ok=1
-  if ( abs(bilh_error) .gt. bilh_seuil) bilh_ok=1
+  if ( abs(bilq_error) > bilq_seuil) bilq_ok=1
+  if ( abs(bilh_error) > bilh_seuil) bilh_ok=1
 !
 ! Print diagnostics
 ! =================
-  if ( (bilq_ok .eq. 0).and.(bilh_ok .eq. 0) ) then
+  if ( (bilq_ok == 0).and.(bilh_ok == 0) ) then
     status="enerbil-OK"
   else
@@ -987,12 +987,12 @@
   end if
 
-  if ( prt_level .GE. 3) then
+  if ( prt_level >= 3) then
     write(*,9010) text,status," itap:",itap,"enerbilERROR: Q", bilq_error,"  H", bilh_error
 9010  format (1x,A8,2x,A12,A6,I4,A18,E15.6,A5,E15.6)
   end if
-  if ( prt_level .GE. 3) then
+  if ( prt_level >= 3) then
     write(*,9000) text,"enerbil: Q,H,KE budget", d_qt_col(1), d_h_col(1),d_ek_col(1)
   end if
-  if ( prt_level .GE. 5) then
+  if ( prt_level >= 5) then
     write(*,9000) text,"enerbil at boundaries: Q, H",bilq_bnd, bilh_bnd
     write(*,9000) text,"enerbil: water budget",d_qt_col(1),d_qw_col(1),d_ql_col(1),d_qs_col(1), d_qbs_col(1)
@@ -1002,15 +1002,15 @@
   specific_diag: SELECT CASE (text)
   CASE("vdf") specific_diag
-    if ( prt_level .GE. 5) then
+    if ( prt_level >= 5) then
       write(*,9000) text,"enerbil: d_h, bilh, sens,t_seri", d_h_col(1), bilh_bnd, sens(1), t_seri(1,1)
       write(*,9000) text,"enerbil: d_h_col_vdf, f_h, diff",d_h_col_vdf, f_h_bnd, bilh_bnd-sens(1)
     end if
   CASE("lsc") specific_diag
-    if ( prt_level .GE. 5) then
+    if ( prt_level >= 5) then
       write(*,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_lsc(1), rlvtt * rain_lsc(1), -(rcw-rcpd)*t_seri(1,1) * rain_lsc(1)
       write(*,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_lsc(1), rlstt * snow_lsc(1), -(rcs-rcpd)*t_seri(1,1) * snow_lsc(1)
     end if
   CASE("convection") specific_diag
-    if ( prt_level .GE. 5) then
+    if ( prt_level >= 5) then
       write(*,9000) text,"enerbil: rain, bil_lat, bil_sens", rain_con(1), rlvtt * rain_con(1), -(rcw-rcpd)*t_seri(1,1) * rain_con(1)
       write(*,9000) text,"enerbil: snow, bil_lat, bil_sens", snow_con(1), rlstt * snow_con(1), -(rcs-rcpd)*t_seri(1,1) * snow_con(1)
Index: /LMDZ6/branches/Amaury_dev/libf/phylmdiso/physiq_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmdiso/physiq_mod.F90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmdiso/physiq_mod.F90	(revision 5081)
@@ -6175,6 +6175,6 @@
     itaprad = itaprad + 1
 
-    IF (iflag_radia.eq.0) THEN
-       IF (prt_level.ge.9) THEN
+    IF (iflag_radia==0) THEN
+       IF (prt_level>=9) THEN
           PRINT *,'--------------------------------------------------'
           PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
@@ -6258,5 +6258,5 @@
     ! a l'echelle sous-maille:
     !
-    IF (prt_level .GE.10) THEN
+    IF (prt_level >=10) THEN
        print *,' call orography ? ', ok_orodr
     ENDIF
@@ -6274,5 +6274,5 @@
           ! zpmm_orodr_t and zstd_orodr_t are activation thresholds set by F. Lott to
           ! earn computation time but they are not physical.
-          IF (((zpic(i)-zmea(i)).GT.zpmm_orodr_t).AND.(zstd(i).GT.zstd_orodr_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
+          IF (((zpic(i)-zmea(i))>zpmm_orodr_t).AND.(zstd(i)>zstd_orodr_t).AND.(zrel_oro(i)<=zrel_oro_t)) THEN
              itest(i)=1
              igwd=igwd+1
@@ -6331,5 +6331,5 @@
           !    such as ice sheets (work by V. Wiener)
           zrel_oro(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
-          IF (((zpic(i)-zmea(i)).GT.zpmm_orolf_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
+          IF (((zpic(i)-zmea(i))>zpmm_orolf_t).AND.(zrel_oro(i)<=zrel_oro_t)) THEN
              itest(i)=1
              igwd=igwd+1
@@ -6555,5 +6555,5 @@
        CALL getin_p('addtkeoro',addtkeoro) 
       
-       IF (prt_level.ge.5) &
+       IF (prt_level>=5) &
             print*,'addtkeoro', addtkeoro
             
@@ -6576,5 +6576,5 @@
 !
 
-       IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
+       IF (addtkeoro > 0 .AND. ok_orodr ) THEN
 !      -------------------------------------------
 
@@ -6583,10 +6583,10 @@
 
 
-  IF (addtkeoro .EQ. 1 ) THEN
+  IF (addtkeoro == 1 ) THEN
 
             duadd(:,:)=alphatkeoro*d_u_oro(:,:)
             dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
 
-  ELSE IF (addtkeoro .EQ. 2) THEN
+  ELSE IF (addtkeoro == 2) THEN
 
      IF (smallscales_tkeoro) THEN
@@ -6597,5 +6597,5 @@
 ! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
 ! Mais attention, cela ne va pas dans le sens de la conservation de l'energie! 
-          IF ((zstd(i).GT.1.0) .AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
+          IF ((zstd(i)>1.0) .AND.(zrel_oro(i)<=zrel_oro_t)) THEN
              itest(i)=1
              igwd=igwd+1
@@ -6609,5 +6609,5 @@
        DO i=1,klon
           itest(i)=0
-        IF (((zpic(i)-zmea(i)).GT.zpmm_orodr_t).AND.(zstd(i).GT.zstd_orodr_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
+        IF (((zpic(i)-zmea(i))>zpmm_orodr_t).AND.(zstd(i)>zstd_orodr_t).AND.(zrel_oro(i)<=zrel_oro_t)) THEN
              itest(i)=1
              igwd=igwd+1
@@ -6831,5 +6831,5 @@
   IF (ok_airs) then
 
-  IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
+  IF (itap==1.or.MOD(itap,NINT(freq_airs/phys_tstep))==0) THEN
      write(*,*) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
      CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
@@ -6938,5 +6938,5 @@
     IF (offline) THEN
 
-       IF (prt_level.ge.9) &
+       IF (prt_level>=9) &
             print*,'Attention on met a 0 les thermiques pour phystoke'
        CALL phystokenc ( &
@@ -6959,5 +6959,5 @@
     !
     !IM global posePB BEG
-    IF(1.EQ.0) THEN
+    IF(1==0) THEN
        !
        CALL transp_lay (paprs,zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, &
@@ -7118,5 +7118,5 @@
     ! Convertir les incrementations en tendances
     !
-    IF (prt_level .GE.10) THEN
+    IF (prt_level >=10) THEN
        print *,'Convertir les incrementations en tendances '
     ENDIF
@@ -7141,9 +7141,9 @@
           ENDIF
           !--ice_sursat: nqo=4, on ajoute rneb
-          IF (nqo.ge.4 .and. ok_ice_sursat) THEN
+          IF (nqo>=4 .and. ok_ice_sursat) THEN
              d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
           ENDIF
 
-           IF (nqo.ge.4 .and. ok_bs) THEN
+           IF (nqo>=4 .and. ok_bs) THEN
              d_qx(i,k,ibs) = ( qbs_seri(i,k) - qx(i,k,ibs) ) / phys_tstep
           ENDIF
@@ -7216,5 +7216,5 @@
     !==========================================================================
 
-    IF (prt_level.ge.1) THEN
+    IF (prt_level>=1) THEN
        write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
        write(lunout,*) &
Index: /LMDZ6/branches/Amaury_dev/tools/netcdf95/nf95_abort.f90
===================================================================
--- /LMDZ6/branches/Amaury_dev/tools/netcdf95/nf95_abort.f90	(revision 5080)
+++ /LMDZ6/branches/Amaury_dev/tools/netcdf95/nf95_abort.f90	(revision 5081)
@@ -7,7 +7,7 @@
   subroutine nf95_abort(message, ncerr, ncid, varid)
 
-    use, intrinsic:: iso_fortran_env
+    use, intrinsic:: iso_fortran_env, ONLY: error_unit
 
-    ! Libraries:
+    ! Libraries:mod_phys_lmdz_mpi_transfert
     use lmdz_netcdf, only: nf90_strerror
 
