Index: LMDZ6/trunk/libf/dyn3d_common/advx.F
===================================================================
--- LMDZ6/trunk/libf/dyn3d_common/advx.F	(revision 5079)
+++ LMDZ6/trunk/libf/dyn3d_common/advx.F	(revision 5084)
@@ -95,7 +95,7 @@
 
 C  -------------------------------------
-      DO j = 1,jjp1 
+      DO 300 j = 1,jjp1 
          NUM(j) = 1
-      END DO
+  300 CONTINUE
       sqi = 0.
       sqf = 0.
@@ -121,12 +121,10 @@
 C  ugri est en kg/s
 
-      DO l = 1,llm
-         DO j = 1,jjm+1
-            DO i = 1,iip1
+      DO 500 l = 1,llm
+         DO 500 j = 1,jjm+1
+            DO 500 i = 1,iip1  
 C            ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
              ugri (i,j,llm+1-l) = pbaru (i,j,l)
-      END DO
-      END DO
-      END DO
+  500 CONTINUE
 
 
@@ -139,6 +137,6 @@
 C  boucle principale sur les niveaux et les latitudes
 C
-      DO L=1,NIV
-      DO K=lati,latf
+      DO 1 L=1,NIV
+      DO 1 K=lati,latf
 C
 C  initialisation
@@ -146,7 +144,7 @@
 C  program assumes periodic boundaries in X
 C
-      DO I=2,LON
+      DO 10 I=2,LON
          SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
-      END DO
+ 10   CONTINUE
       SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
 C
@@ -156,26 +154,26 @@
       LONK=LON/NUMK
 C
-      IF(NUMK>1) THEN
-C
-      DO I=1,LON
+      IF(NUMK.GT.1) THEN
+C
+      DO 111 I=1,LON
          TM(I)=0.
-      END DO
-      DO JV=1,NTRA
-      DO I=1,LON
+ 111  CONTINUE
+      DO 112 JV=1,NTRA
+      DO 1120 I=1,LON
          T0(I,JV)=0.
          TX(I,JV)=0.
          TY(I,JV)=0.
          TZ(I,JV)=0.
-      END DO
-      END DO
-C
-      DO I2=1,NUMK
-C
-         DO I=1,LONK
+ 1120 CONTINUE
+ 112  CONTINUE
+C
+      DO 11 I2=1,NUMK
+C
+         DO 113 I=1,LONK
             I3=(I-1)*NUMK+I2
             TM(I)=TM(I)+SM(I3,K,L)
             ALF(I)=SM(I3,K,L)/TM(I)
             ALF1(I)=1.-ALF(I)
-      END DO
+ 113     CONTINUE
 C
          DO  JV=1,NTRA
@@ -192,25 +190,25 @@
          ENDDO
 C
-      END DO
+ 11   CONTINUE
 C
       ELSE
 C
-      DO I=1,LON
+      DO 115 I=1,LON
          TM(I)=SM(I,K,L)
-      END DO
-      DO JV=1,NTRA
-      DO I=1,LON
+ 115  CONTINUE
+      DO 116 JV=1,NTRA
+      DO 1160 I=1,LON
          T0(I,JV)=S0(I,K,L,JV)
          TX(I,JV)=sx(I,K,L,JV)
          TY(I,JV)=sy(I,K,L,JV)
          TZ(I,JV)=sz(I,K,L,JV)
-      END DO
-      END DO
-C
-      ENDIF
-C
-      DO I=1,LONK
+ 1160 CONTINUE
+ 116  CONTINUE
+C
+      ENDIF
+C
+      DO 117 I=1,LONK
          UEXT(I)=UGRI(I*NUMK,K,L)
-      END DO
+ 117  CONTINUE
 C
 C  place limits on appropriate moments before transport
@@ -219,9 +217,9 @@
       IF(.NOT.LIMIT) GO TO 13
 C
-      DO JV=1,NTRA
-      DO I=1,LONK
+      DO 12 JV=1,NTRA
+      DO 120 I=1,LONK
         TX(I,JV)=SIGN(AMIN1(AMAX1(T0(I,JV),0.),ABS(TX(I,JV))),TX(I,JV))
-      END DO
-      END DO
+ 120  CONTINUE
+ 12   CONTINUE
 C
  13   CONTINUE
@@ -233,14 +231,14 @@
 C  flux from IP to I if U(I).lt.0
 C
-      DO I=1,LONK-1
-         IF(UEXT(I)<0.) THEN
+      DO 140 I=1,LONK-1
+         IF(UEXT(I).LT.0.) THEN
            FM(I)=-UEXT(I)*DTX
            ALF(I)=FM(I)/TM(I+1)
            TM(I+1)=TM(I+1)-FM(I)
          ENDIF
-      END DO
+ 140  CONTINUE
 C
       I=LONK
-      IF(UEXT(I)<0.) THEN
+      IF(UEXT(I).LT.0.) THEN
         FM(I)=-UEXT(I)*DTX
         ALF(I)=FM(I)/TM(1)
@@ -250,22 +248,22 @@
 C  flux from I to IP if U(I).gt.0
 C
-      DO I=1,LONK
-         IF(UEXT(I)>=0.) THEN
+      DO 141 I=1,LONK
+         IF(UEXT(I).GE.0.) THEN
            FM(I)=UEXT(I)*DTX
            ALF(I)=FM(I)/TM(I)
            TM(I)=TM(I)-FM(I)
          ENDIF
-      END DO
-C
-      DO I=1,LONK
+ 141  CONTINUE
+C
+      DO 142 I=1,LONK
          ALFQ(I)=ALF(I)*ALF(I)
          ALF1(I)=1.-ALF(I)
          ALF1Q(I)=ALF1(I)*ALF1(I)
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LONK-1
-C
-         IF(UEXT(I)<0.) THEN
+ 142  CONTINUE
+C
+      DO 150 JV=1,NTRA
+      DO 1500 I=1,LONK-1
+C
+         IF(UEXT(I).LT.0.) THEN
 C
            F0(I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*TX(I+1,JV) )
@@ -281,11 +279,11 @@
          ENDIF
 C
-      END DO
-      END DO
+ 1500 CONTINUE
+ 150  CONTINUE
 C
       I=LONK
-      IF(UEXT(I)<0.) THEN
-C
-        DO JV=1,NTRA
+      IF(UEXT(I).LT.0.) THEN
+C
+        DO 151 JV=1,NTRA
 C
            F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*TX(1,JV) )
@@ -299,12 +297,12 @@
            TZ(1,JV)=TZ(1,JV)-FZ(I,JV)
 C
-      END DO
-C
-      ENDIF
-C
-      DO JV=1,NTRA
-      DO I=1,LONK
-C
-         IF(UEXT(I)>=0.) THEN
+ 151    CONTINUE
+C
+      ENDIF
+C
+      DO 152 JV=1,NTRA
+      DO 1520 I=1,LONK
+C
+         IF(UEXT(I).GE.0.) THEN
 C
            F0(I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*TX(I,JV) )
@@ -320,37 +318,37 @@
          ENDIF
 C
-      END DO
-      END DO
+ 1520 CONTINUE
+ 152  CONTINUE
 C
 C  puts the temporary moments Fi into appropriate neighboring boxes
 C
-      DO I=1,LONK
-         IF(UEXT(I)<0.) THEN
+      DO 160 I=1,LONK
+         IF(UEXT(I).LT.0.) THEN
            TM(I)=TM(I)+FM(I)
            ALF(I)=FM(I)/TM(I)
          ENDIF
-      END DO
-C
-      DO I=1,LONK-1
-         IF(UEXT(I)>=0.) THEN
+ 160  CONTINUE
+C
+      DO 161 I=1,LONK-1
+         IF(UEXT(I).GE.0.) THEN
            TM(I+1)=TM(I+1)+FM(I)
            ALF(I)=FM(I)/TM(I+1)
          ENDIF
-      END DO
+ 161  CONTINUE
 C
       I=LONK
-      IF(UEXT(I)>=0.) THEN
+      IF(UEXT(I).GE.0.) THEN
         TM(1)=TM(1)+FM(I)
         ALF(I)=FM(I)/TM(1)
       ENDIF
 C
-      DO I=1,LONK
+      DO 162 I=1,LONK
          ALF1(I)=1.-ALF(I)
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LONK
-C
-         IF(UEXT(I)<0.) THEN
+ 162  CONTINUE
+C
+      DO 170 JV=1,NTRA
+      DO 1700 I=1,LONK
+C
+         IF(UEXT(I).LT.0.) THEN
 C
            TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*F0(I,JV)
@@ -362,11 +360,11 @@
          ENDIF
 C
-      END DO
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LONK-1
-C
-         IF(UEXT(I)>=0.) THEN
+ 1700 CONTINUE
+ 170  CONTINUE
+C
+      DO 171 JV=1,NTRA
+      DO 1710 I=1,LONK-1
+C
+         IF(UEXT(I).GE.0.) THEN
 C
            TEMPTM=ALF(I)*T0(I+1,JV)-ALF1(I)*F0(I,JV)
@@ -378,10 +376,10 @@
          ENDIF
 C
-      END DO
-      END DO
+ 1710 CONTINUE
+ 171  CONTINUE
 C
       I=LONK
-      IF(UEXT(I)>=0.) THEN
-        DO JV=1,NTRA
+      IF(UEXT(I).GE.0.) THEN
+        DO 172 JV=1,NTRA
            TEMPTM=ALF(I)*T0(1,JV)-ALF1(I)*F0(I,JV)
            T0(1,JV)=T0(1,JV)+F0(I,JV)
@@ -389,14 +387,14 @@
            TY(1,JV)=TY(1,JV)+FY(I,JV)
            TZ(1,JV)=TZ(1,JV)+FZ(I,JV)
-      END DO
+ 172    CONTINUE
       ENDIF
 C
 C  retour aux mailles d'origine (passage des Tij aux Sij)
 C
-      IF(NUMK>1) THEN
-C
-      DO I2=1,NUMK
-C
-         DO I=1,LONK
+      IF(NUMK.GT.1) THEN
+C
+      DO 180 I2=1,NUMK
+C
+         DO 180 I=1,LONK
 C
             I3=I2+(I-1)*NUMK
@@ -409,6 +407,5 @@
             ALF1Q(I)=ALF1(I)*ALF1(I)
 C
-      END DO
-      END DO
+ 180     CONTINUE
 C
          DO  JV=1,NTRA
@@ -434,20 +431,19 @@
       ELSE
 C
-      DO I=1,LON
+      DO 190 I=1,LON
          SM(I,K,L)=TM(I)
-      END DO
-      DO JV=1,NTRA
-      DO I=1,LON
+ 190  CONTINUE
+      DO 191 JV=1,NTRA
+      DO 1910 I=1,LON
          S0(I,K,L,JV)=T0(I,JV)
          sx(I,K,L,JV)=TX(I,JV)
          sy(I,K,L,JV)=TY(I,JV)
          sz(I,K,L,JV)=TZ(I,JV)
-      END DO
-      END DO
-C
-      ENDIF
-C
-      END DO
-      END DO
+ 1910 CONTINUE
+ 191  CONTINUE
+C
+      ENDIF
+C
+ 1    CONTINUE
 C
 C ----------- AA Test en fin de ADVX ------ Controle des S*
Index: LMDZ6/trunk/libf/dyn3d_common/advxp.F
===================================================================
--- LMDZ6/trunk/libf/dyn3d_common/advxp.F	(revision 5079)
+++ LMDZ6/trunk/libf/dyn3d_common/advxp.F	(revision 5084)
@@ -126,7 +126,7 @@
 c test
 c  -------------------------------------
-        DO j =1,jjp1
+        DO 300 j =1,jjp1
          NUM(j) =1 
-      END DO
+ 300  CONTINUE
 c       DO l=1,llm
 c      NUM(2,l)=6
@@ -150,11 +150,9 @@
 C  ugri est en kg/s
 
-       DO l = 1,llm
-       DO j = 1,jjp1
-       DO i = 1,iip1
+       DO 500 l = 1,llm
+       DO 500 j = 1,jjp1
+       DO 500 i = 1,iip1
        ugri (i,j,llm+1-l) =pbaru (i,j,l) 
-      END DO
-      END DO
-      END DO
+ 500   CONTINUE
 
 C  ---------------------------------------------------------
@@ -163,6 +161,6 @@
 C  boucle principale sur les niveaux et les latitudes
 C     
-      DO L=1,NIV
-      DO K=lati,latf
+      DO 1 L=1,NIV
+      DO 1 K=lati,latf
 
 C
@@ -171,7 +169,7 @@
 C  program assumes periodic boundaries in X
 C
-      DO I=2,LON
+      DO 10 I=2,LON
          SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
-      END DO
+ 10   CONTINUE
       SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
 C
@@ -181,11 +179,11 @@
       LONK=LON/NUMK
 C
-      IF(NUMK>1) THEN
-C
-      DO I=1,LON
+      IF(NUMK.GT.1) THEN
+C
+      DO 111 I=1,LON
          TM(I)=0.
-      END DO
-      DO JV=1,NTRA
-      DO I=1,LON
+ 111  CONTINUE
+      DO 112 JV=1,NTRA
+      DO 1120 I=1,LON
          T0 (I,JV)=0.
          TX (I,JV)=0.
@@ -198,10 +196,10 @@
          TYZ(I,JV)=0.
          TZZ(I,JV)=0.
-      END DO
-      END DO
-C
-      DO I2=1,NUMK
-C
-         DO I=1,LONK
+ 1120 CONTINUE
+ 112  CONTINUE
+C
+      DO 11 I2=1,NUMK
+C
+         DO 113 I=1,LONK
             I3=(I-1)*NUMK+I2
             TM(I)=TM(I)+SM(I3,K,L)
@@ -212,8 +210,8 @@
             ALF2(I)=ALF1(I)-ALF(I)
             ALF3(I)=ALF(I)*ALF1(I)
-      END DO
-C
-         DO JV=1,NTRA
-         DO I=1,LONK
+ 113     CONTINUE
+C
+         DO 114 JV=1,NTRA
+         DO 1140 I=1,LONK
             I3=(I-1)*NUMK+I2
             TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*S0(I3,K,L,JV)
@@ -231,16 +229,16 @@
             TYZ(I,JV)=TYZ(I,JV)+SYZ(I3,K,L,JV)
             TZZ(I,JV)=TZZ(I,JV)+SZZ(I3,K,L,JV)
-      END DO
-      END DO
-C
-      END DO
+ 1140    CONTINUE
+ 114     CONTINUE
+C
+ 11   CONTINUE
 C
       ELSE
 C
-      DO I=1,LON
+      DO 115 I=1,LON
          TM(I)=SM(I,K,L)
-      END DO
-      DO JV=1,NTRA
-      DO I=1,LON
+ 115  CONTINUE
+      DO 116 JV=1,NTRA
+      DO 1160 I=1,LON
          T0 (I,JV)=S0 (I,K,L,JV)
          TX (I,JV)=SSX (I,K,L,JV)
@@ -253,12 +251,12 @@
          TYZ(I,JV)=SYZ(I,K,L,JV)
          TZZ(I,JV)=SZZ(I,K,L,JV)
-      END DO
-      END DO
+ 1160 CONTINUE
+ 116  CONTINUE
 C
       ENDIF
 C
-      DO I=1,LONK
+      DO 117 I=1,LONK
          UEXT(I)=UGRI(I*NUMK,K,L)
-      END DO
+ 117  CONTINUE
 C
 C  place limits on appropriate moments before transport
@@ -267,7 +265,7 @@
       IF(.NOT.LIMIT) GO TO 13
 C
-      DO JV=1,NTRA
-      DO I=1,LONK
-        IF(T0(I,JV)>0.) THEN
+      DO 12 JV=1,NTRA
+      DO 120 I=1,LONK
+        IF(T0(I,JV).GT.0.) THEN
           SLPMAX=T0(I,JV)
           S1MAX=1.5*SLPMAX
@@ -285,6 +283,6 @@
           TXZ(I,JV)=0.
         ENDIF
-      END DO
-      END DO
+ 120  CONTINUE
+ 12   CONTINUE
 C
  13   CONTINUE
@@ -296,14 +294,14 @@
 C  flux from IP to I if U(I).lt.0
 C
-      DO I=1,LONK-1
-         IF(UEXT(I)<0.) THEN
+      DO 140 I=1,LONK-1
+         IF(UEXT(I).LT.0.) THEN
            FM(I)=-UEXT(I)*DTX
            ALF(I)=FM(I)/TM(I+1)
            TM(I+1)=TM(I+1)-FM(I)
          ENDIF
-      END DO
+ 140  CONTINUE
 C
       I=LONK
-      IF(UEXT(I)<0.) THEN
+      IF(UEXT(I).LT.0.) THEN
         FM(I)=-UEXT(I)*DTX
         ALF(I)=FM(I)/TM(1)
@@ -313,13 +311,13 @@
 C  flux from I to IP if U(I).gt.0
 C
-      DO I=1,LONK
-         IF(UEXT(I)>=0.) THEN
+      DO 141 I=1,LONK
+         IF(UEXT(I).GE.0.) THEN
            FM(I)=UEXT(I)*DTX
            ALF(I)=FM(I)/TM(I)
            TM(I)=TM(I)-FM(I)
          ENDIF
-      END DO
-C
-      DO I=1,LONK
+ 141  CONTINUE
+C
+      DO 142 I=1,LONK
          ALFQ(I)=ALF(I)*ALF(I)
          ALF1(I)=1.-ALF(I)
@@ -328,10 +326,10 @@
          ALF3(I)=ALF(I)*ALFQ(I)
          ALF4(I)=ALF1(I)*ALF1Q(I)
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LONK-1
-C
-         IF(UEXT(I)<0.) THEN
+ 142  CONTINUE
+C
+      DO 150 JV=1,NTRA
+      DO 1500 I=1,LONK-1
+C
+         IF(UEXT(I).LT.0.) THEN
 C
            F0 (I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*
@@ -360,11 +358,11 @@
          ENDIF
 C
-      END DO
-      END DO
+ 1500 CONTINUE
+ 150  CONTINUE
 C
       I=LONK
-      IF(UEXT(I)<0.) THEN
-C
-        DO JV=1,NTRA
+      IF(UEXT(I).LT.0.) THEN
+C
+        DO 151 JV=1,NTRA
 C
            F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*
@@ -391,12 +389,12 @@
            TXZ(1,JV)=ALF1Q(I)*TXZ(1,JV)
 C
-      END DO
+ 151    CONTINUE
 C
       ENDIF
 C
-      DO JV=1,NTRA
-      DO I=1,LONK
-C
-         IF(UEXT(I)>=0.) THEN
+      DO 152 JV=1,NTRA
+      DO 1520 I=1,LONK
+C
+         IF(UEXT(I).GE.0.) THEN
 C
            F0 (I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*
@@ -425,30 +423,30 @@
          ENDIF
 C
-      END DO
-      END DO
+ 1520 CONTINUE
+ 152  CONTINUE
 C
 C  puts the temporary moments Fi into appropriate neighboring boxes
 C
-      DO I=1,LONK
-         IF(UEXT(I)<0.) THEN
+      DO 160 I=1,LONK
+         IF(UEXT(I).LT.0.) THEN
            TM(I)=TM(I)+FM(I)
            ALF(I)=FM(I)/TM(I)
          ENDIF
-      END DO
-C
-      DO I=1,LONK-1
-         IF(UEXT(I)>=0.) THEN
+ 160  CONTINUE
+C
+      DO 161 I=1,LONK-1
+         IF(UEXT(I).GE.0.) THEN
            TM(I+1)=TM(I+1)+FM(I)
            ALF(I)=FM(I)/TM(I+1)
          ENDIF
-      END DO
+ 161  CONTINUE
 C
       I=LONK
-      IF(UEXT(I)>=0.) THEN
+      IF(UEXT(I).GE.0.) THEN
         TM(1)=TM(1)+FM(I)
         ALF(I)=FM(I)/TM(1)
       ENDIF
 C
-      DO I=1,LONK
+      DO 162 I=1,LONK
          ALF1(I)=1.-ALF(I)
          ALFQ(I)=ALF(I)*ALF(I)
@@ -456,10 +454,10 @@
          ALF2(I)=ALF1(I)-ALF(I)
          ALF3(I)=ALF(I)*ALF1(I)
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LONK
-C
-         IF(UEXT(I)<0.) THEN
+ 162  CONTINUE
+C
+      DO 170 JV=1,NTRA
+      DO 1700 I=1,LONK
+C
+         IF(UEXT(I).LT.0.) THEN
 C
            TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*F0(I,JV)
@@ -480,11 +478,11 @@
          ENDIF
 C
-      END DO
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LONK-1
-C
-         IF(UEXT(I)>=0.) THEN
+ 1700 CONTINUE
+ 170  CONTINUE
+C
+      DO 171 JV=1,NTRA
+      DO 1710 I=1,LONK-1
+C
+         IF(UEXT(I).GE.0.) THEN
 C
            TEMPTM=ALF(I)*T0(I+1,JV)-ALF1(I)*F0(I,JV)
@@ -505,10 +503,10 @@
          ENDIF
 C
-      END DO
-      END DO
+ 1710 CONTINUE
+ 171  CONTINUE
 C
       I=LONK
-      IF(UEXT(I)>=0.) THEN
-        DO JV=1,NTRA
+      IF(UEXT(I).GE.0.) THEN
+        DO 172 JV=1,NTRA
            TEMPTM=ALF(I)*T0(1,JV)-ALF1(I)*F0(I,JV)
            T0 (1,JV)=T0(1,JV)+F0(I,JV)
@@ -525,14 +523,14 @@
            TYZ(1,JV)=TYZ(1,JV)+FYZ(I,JV)
            TZZ(1,JV)=TZZ(1,JV)+FZZ(I,JV)
-      END DO
+ 172    CONTINUE
       ENDIF
 C
 C  retour aux mailles d'origine (passage des Tij aux Sij)
 C
-      IF(NUMK>1) THEN
-C
-      DO I2=1,NUMK
-C
-         DO I=1,LONK
+      IF(NUMK.GT.1) THEN
+C
+      DO 18 I2=1,NUMK
+C
+         DO 180 I=1,LONK
 C
             I3=I2+(I-1)*NUMK
@@ -548,8 +546,8 @@
             ALF4(I)=ALF1(I)*ALF1Q(I)
 C
-      END DO
-C
-         DO JV=1,NTRA
-         DO I=1,LONK
+ 180     CONTINUE
+C
+         DO 181 JV=1,NTRA
+         DO 181 I=1,LONK
 C
             I3=I2+(I-1)*NUMK
@@ -579,16 +577,15 @@
             TXZ(I,JV)=ALF1Q(I)*TXZ(I,JV)
 C
-      END DO
-      END DO
-C
-      END DO
+ 181     CONTINUE
+C
+ 18   CONTINUE
 C
       ELSE
 C
-      DO I=1,LON
+      DO 190 I=1,LON
          SM(I,K,L)=TM(I)
-      END DO
-      DO JV=1,NTRA
-      DO I=1,LON
+ 190  CONTINUE
+      DO 191 JV=1,NTRA
+      DO 1910 I=1,LON
          S0 (I,K,L,JV)=T0 (I,JV)
          SSX (I,K,L,JV)=TX (I,JV)
@@ -601,11 +598,10 @@
          SYZ(I,K,L,JV)=TYZ(I,JV)
          SZZ(I,K,L,JV)=TZZ(I,JV)
-      END DO
-      END DO
+ 1910 CONTINUE
+ 191  CONTINUE
 C
       ENDIF
 C
-      END DO
-      END DO
+ 1    CONTINUE
 C
 C ----------- AA Test en fin de ADVX ------ Controle des S*
Index: LMDZ6/trunk/libf/dyn3d_common/advy.F
===================================================================
--- LMDZ6/trunk/libf/dyn3d_common/advy.F	(revision 5079)
+++ LMDZ6/trunk/libf/dyn3d_common/advy.F	(revision 5084)
@@ -121,5 +121,5 @@
       enddo
 
-      DO L=1,NIV
+      DO 1 L=1,NIV
 C
 C  place limits on appropriate moments before transport
@@ -128,12 +128,11 @@
       IF(.NOT.LIMIT) GO TO 11
 C
-      DO JV=1,NTRA
-      DO K=1,LAT
-      DO I=1,LON
+      DO 10 JV=1,NTRA
+      DO 10 K=1,LAT
+      DO 100 I=1,LON
          sy(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.),
      +                           ABS(sy(I,K,L,JV))),sy(I,K,L,JV))
-      END DO
-      END DO
-      END DO
+ 100  CONTINUE
+ 10   CONTINUE
 C
  11   CONTINUE
@@ -142,11 +141,11 @@
 C
       SM0=0.
-      DO JV=1,NTRA
+      DO 20 JV=1,NTRA
          S00(JV)=0.
-      END DO
-C
-      DO I=1,LON
-C
-         IF(VGRI(I,0,L)<=0.) THEN
+ 20   CONTINUE
+C
+      DO 21 I=1,LON
+C
+         IF(VGRI(I,0,L).LE.0.) THEN
            FM(I,0)=-VGRI(I,0,L)*DTY
            ALF(I,0)=FM(I,0)/SM(I,1,L)
@@ -159,10 +158,10 @@
          ALF1Q(I,0)=ALF1(I,0)*ALF1(I,0)
 C
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LON
-C
-         IF(VGRI(I,0,L)<=0.) THEN
+ 21   CONTINUE
+C
+      DO 22 JV=1,NTRA
+      DO 220 I=1,LON
+C
+         IF(VGRI(I,0,L).LE.0.) THEN
 C
            F0(I,0,JV)=ALF(I,0)*
@@ -177,27 +176,27 @@
          ENDIF
 C
-      END DO
-      END DO
-C
-      DO I=1,LON
-         IF(VGRI(I,0,L)>0.) THEN
+ 220  CONTINUE
+ 22   CONTINUE
+C
+      DO 23 I=1,LON
+         IF(VGRI(I,0,L).GT.0.) THEN
            FM(I,0)=VGRI(I,0,L)*DTY
            ALF(I,0)=FM(I,0)/SM0
          ENDIF
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LON
-         IF(VGRI(I,0,L)>0.) THEN
+ 23   CONTINUE
+C
+      DO 24 JV=1,NTRA
+      DO 240 I=1,LON
+         IF(VGRI(I,0,L).GT.0.) THEN
            F0(I,0,JV)=ALF(I,0)*S00(JV)
          ENDIF
-      END DO
-      END DO
+ 240  CONTINUE
+ 24   CONTINUE
 C
 C  puts the temporary moments Fi into appropriate neighboring boxes
 C
-      DO I=1,LON
-C
-         IF(VGRI(I,0,L)>0.) THEN
+      DO 25 I=1,LON
+C
+         IF(VGRI(I,0,L).GT.0.) THEN
            SM(I,1,L)=SM(I,1,L)+FM(I,0)
            ALF(I,0)=FM(I,0)/SM(I,1,L)
@@ -206,10 +205,10 @@
          ALF1(I,0)=1.-ALF(I,0)
 C
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LON
-C
-         IF(VGRI(I,0,L)>0.) THEN
+ 25   CONTINUE
+C
+      DO 26 JV=1,NTRA
+      DO 260 I=1,LON
+C
+         IF(VGRI(I,0,L).GT.0.) THEN
 C
          TEMPTM=ALF(I,0)*S0(I,1,L,JV)-ALF1(I,0)*F0(I,0,JV)
@@ -219,6 +218,6 @@
          ENDIF
 C
-      END DO
-      END DO
+ 260  CONTINUE
+ 26   CONTINUE
 C
 C  calculate flux and moments between adjacent boxes
@@ -228,9 +227,9 @@
 C  flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
 C
-      DO K=1,LAT-1
+      DO 30 K=1,LAT-1
       KP=K+1
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)<0.) THEN
+      DO 300 I=1,LON
+C
+         IF(VGRI(I,K,L).LT.0.) THEN
            FM(I,K)=-VGRI(I,K,L)*DTY
            ALF(I,K)=FM(I,K)/SM(I,KP,L)
@@ -246,13 +245,13 @@
          ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
 C
-      END DO
-      END DO
-C
-      DO JV=1,NTRA
-      DO K=1,LAT-1
+ 300  CONTINUE
+ 30   CONTINUE
+C
+      DO 31 JV=1,NTRA
+      DO 31 K=1,LAT-1
       KP=K+1
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)<0.) THEN
+      DO 310 I=1,LON
+C
+         IF(VGRI(I,K,L).LT.0.) THEN
 C
            F0(I,K,JV)=ALF (I,K)*
@@ -282,15 +281,14 @@
          ENDIF
 C
-      END DO
-      END DO
-      END DO
+ 310  CONTINUE
+ 31   CONTINUE
 C
 C  puts the temporary moments Fi into appropriate neighboring boxes
 C
-      DO K=1,LAT-1
+      DO 32 K=1,LAT-1
       KP=K+1
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)<0.) THEN
+      DO 320 I=1,LON
+C
+         IF(VGRI(I,K,L).LT.0.) THEN
            SM(I,K,L)=SM(I,K,L)+FM(I,K)
            ALF(I,K)=FM(I,K)/SM(I,K,L)
@@ -302,13 +300,13 @@
          ALF1(I,K)=1.-ALF(I,K)
 C
-      END DO
-      END DO
-C
-      DO JV=1,NTRA
-      DO K=1,LAT-1
+ 320  CONTINUE
+ 32   CONTINUE
+C
+      DO 33 JV=1,NTRA
+      DO 33 K=1,LAT-1
       KP=K+1
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)<0.) THEN
+      DO 330 I=1,LON
+C
+         IF(VGRI(I,K,L).LT.0.) THEN
 C
          TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
@@ -330,7 +328,6 @@
          ENDIF
 C
-      END DO
-      END DO
-      END DO
+ 330  CONTINUE
+ 33   CONTINUE
 C
 C  traitement special pour le pole Sud (idem pole Nord)
@@ -339,11 +336,11 @@
 C
       SM0=0.
-      DO JV=1,NTRA
+      DO 40 JV=1,NTRA
          S00(JV)=0.
-      END DO
-C
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)>=0.) THEN
+ 40   CONTINUE
+C
+      DO 41 I=1,LON
+C
+         IF(VGRI(I,K,L).GE.0.) THEN
            FM(I,K)=VGRI(I,K,L)*DTY
            ALF(I,K)=FM(I,K)/SM(I,K,L)
@@ -356,10 +353,10 @@
          ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
 C
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)>=0.) THEN
+ 41   CONTINUE
+C
+      DO 42 JV=1,NTRA
+      DO 420 I=1,LON
+C
+         IF(VGRI(I,K,L).GE.0.) THEN
            F0 (I,K,JV)=ALF(I,K)*
      +                ( S0(I,K,L,JV)+ALF1(I,K)*sy(I,K,L,JV) )
@@ -372,27 +369,27 @@
          ENDIF
 C
-      END DO
-      END DO
-C
-      DO I=1,LON
-         IF(VGRI(I,K,L)<0.) THEN
+ 420  CONTINUE
+ 42   CONTINUE
+C
+      DO 43 I=1,LON
+         IF(VGRI(I,K,L).LT.0.) THEN
            FM(I,K)=-VGRI(I,K,L)*DTY
            ALF(I,K)=FM(I,K)/SM0
          ENDIF
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LON
-         IF(VGRI(I,K,L)<0.) THEN
+ 43   CONTINUE
+C
+      DO 44 JV=1,NTRA
+      DO 440 I=1,LON
+         IF(VGRI(I,K,L).LT.0.) THEN
            F0(I,K,JV)=ALF(I,K)*S00(JV)
          ENDIF
-      END DO
-      END DO
+ 440  CONTINUE
+ 44   CONTINUE
 C
 C  puts the temporary moments Fi into appropriate neighboring boxes
 C
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)<0.) THEN
+      DO 45 I=1,LON
+C
+         IF(VGRI(I,K,L).LT.0.) THEN
            SM(I,K,L)=SM(I,K,L)+FM(I,K)
            ALF(I,K)=FM(I,K)/SM(I,K,L)
@@ -401,10 +398,10 @@
          ALF1(I,K)=1.-ALF(I,K)
 C
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)<0.) THEN
+ 45   CONTINUE
+C
+      DO 46 JV=1,NTRA
+      DO 460 I=1,LON
+C
+         IF(VGRI(I,K,L).LT.0.) THEN
 C
          TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
@@ -414,8 +411,8 @@
          ENDIF
 C
-      END DO
-      END DO
-C
-      END DO
+ 460  CONTINUE
+ 46   CONTINUE
+C
+ 1    CONTINUE
 C
       RETURN
Index: LMDZ6/trunk/libf/dyn3d_common/advyp.F
===================================================================
--- LMDZ6/trunk/libf/dyn3d_common/advyp.F	(revision 5079)
+++ LMDZ6/trunk/libf/dyn3d_common/advyp.F	(revision 5084)
@@ -153,11 +153,9 @@
 C-AA 20/10/94  le signe -1 est necessaire car indexation opposee
 
-      DO l = 1,llm
-         DO j = 1,jjm
-            DO i = 1,iip1
+      DO 500 l = 1,llm
+         DO 500 j = 1,jjm
+            DO 500 i = 1,iip1  
             vgri (i,j,llm+1-l)=-1.*pbarv (i,j,l)
-      END DO
-      END DO
-      END DO
+  500 CONTINUE
 
 CAA Initialisation de flux fictifs aux bords sup. des boites pol.
@@ -173,5 +171,5 @@
 C  boucle sur les niveaux
 C
-      DO L=1,NIV
+      DO 1 L=1,NIV
 C
 C  place limits on appropriate moments before transport
@@ -180,8 +178,8 @@
       IF(.NOT.LIMIT) GO TO 11
 C
-      DO JV=1,NTRA
-      DO K=1,LAT
-      DO I=1,LON
-         IF(S0(I,K,L,JV)>0.) THEN
+      DO 10 JV=1,NTRA
+      DO 10 K=1,LAT
+      DO 100 I=1,LON
+         IF(S0(I,K,L,JV).GT.0.) THEN
            SLPMAX=AMAX1(S0(I,K,L,JV),0.)
            S1MAX=1.5*SLPMAX
@@ -199,7 +197,6 @@
            SYZ(I,K,L,JV)=0.
          ENDIF
-      END DO
-      END DO
-      END DO
+ 100  CONTINUE
+ 10   CONTINUE
 C
  11   CONTINUE
@@ -208,11 +205,11 @@
 C
       SM0=0.
-      DO JV=1,NTRA
+      DO 20 JV=1,NTRA
          S00(JV)=0.
-      END DO
-C
-      DO I=1,LON
-C
-         IF(VGRI(I,0,L)<=0.) THEN
+ 20   CONTINUE
+C
+      DO 21 I=1,LON
+C
+         IF(VGRI(I,0,L).LE.0.) THEN
            FM(I,0)=-VGRI(I,0,L)*DTY
            ALF(I,0)=FM(I,0)/SM(I,1,L)
@@ -228,11 +225,11 @@
          ALF4(I,0)=ALF1(I,0)*ALF1Q(I,0)
 C
-      END DO
+ 21   CONTINUE
 c     print*,'ADVYP 21'
 C
-      DO JV=1,NTRA
-      DO I=1,LON
-C
-         IF(VGRI(I,0,L)<=0.) THEN
+      DO 22 JV=1,NTRA
+      DO 220 I=1,LON
+C
+         IF(VGRI(I,0,L).LE.0.) THEN
 C
            F0(I,0,JV)=ALF(I,0)* ( S0(I,1,L,JV)-ALF1(I,0)*
@@ -256,28 +253,28 @@
          ENDIF
 C
-      END DO
-      END DO
-C
-      DO I=1,LON
-         IF(VGRI(I,0,L)>0.) THEN
+ 220  CONTINUE
+ 22   CONTINUE
+C
+      DO 23 I=1,LON
+         IF(VGRI(I,0,L).GT.0.) THEN
            FM(I,0)=VGRI(I,0,L)*DTY
            ALF(I,0)=FM(I,0)/SM0
          ENDIF
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LON
-         IF(VGRI(I,0,L)>0.) THEN
+ 23   CONTINUE
+C
+      DO 24 JV=1,NTRA
+      DO 240 I=1,LON
+         IF(VGRI(I,0,L).GT.0.) THEN
            F0(I,0,JV)=ALF(I,0)*S00(JV)
          ENDIF
-      END DO
-      END DO
+ 240  CONTINUE
+ 24   CONTINUE
 C
 C  puts the temporary moments Fi into appropriate neighboring boxes
 C
 c     print*,'av ADVYP 25'
-      DO I=1,LON
-C
-         IF(VGRI(I,0,L)>0.) THEN
+      DO 25 I=1,LON
+C
+         IF(VGRI(I,0,L).GT.0.) THEN
            SM(I,1,L)=SM(I,1,L)+FM(I,0)
            ALF(I,0)=FM(I,0)/SM(I,1,L)
@@ -290,11 +287,11 @@
          ALF3(I,0)=ALF1(I,0)*ALF(I,0)
 C
-      END DO
+ 25   CONTINUE
 c     print*,'av ADVYP 25'
 C
-      DO JV=1,NTRA
-      DO I=1,LON
-C
-         IF(VGRI(I,0,L)>0.) THEN
+      DO 26 JV=1,NTRA
+      DO 260 I=1,LON
+C
+         IF(VGRI(I,0,L).GT.0.) THEN
 C
          TEMPTM=ALF(I,0)*S0(I,1,L,JV)-ALF1(I,0)*F0(I,0,JV)
@@ -308,6 +305,6 @@
          ENDIF
 C
-      END DO
-      END DO
+ 260  CONTINUE
+ 26   CONTINUE
 C
 C  calculate flux and moments between adjacent boxes
@@ -318,9 +315,9 @@
 C
 c     print*,'av ADVYP 30'
-      DO K=1,LAT-1
+      DO 30 K=1,LAT-1
       KP=K+1
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)<0.) THEN
+      DO 300 I=1,LON
+C
+         IF(VGRI(I,K,L).LT.0.) THEN
            FM(I,K)=-VGRI(I,K,L)*DTY
            ALF(I,K)=FM(I,K)/SM(I,KP,L)
@@ -339,14 +336,14 @@
          ALF4(I,K)=ALF1(I,K)*ALF1Q(I,K)
 C
-      END DO
-      END DO
+ 300  CONTINUE
+ 30   CONTINUE
 c     print*,'ap ADVYP 30'
 C
-      DO JV=1,NTRA
-      DO K=1,LAT-1
+      DO 31 JV=1,NTRA
+      DO 31 K=1,LAT-1
       KP=K+1
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)<0.) THEN
+      DO 310 I=1,LON
+C
+         IF(VGRI(I,K,L).LT.0.) THEN
 C
            F0 (I,K,JV)=ALF (I,K)* ( S0(I,KP,L,JV)-ALF1(I,K)*
@@ -406,16 +403,15 @@
          ENDIF
 C
-      END DO
-      END DO
-      END DO
+ 310  CONTINUE
+ 31   CONTINUE
 c     print*,'ap ADVYP 31'
 C
 C  puts the temporary moments Fi into appropriate neighboring boxes
 C
-      DO K=1,LAT-1
+      DO 32 K=1,LAT-1
       KP=K+1
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)<0.) THEN
+      DO 320 I=1,LON
+C
+         IF(VGRI(I,K,L).LT.0.) THEN
            SM(I,K,L)=SM(I,K,L)+FM(I,K)
            ALF(I,K)=FM(I,K)/SM(I,K,L)
@@ -431,14 +427,14 @@
          ALF3(I,K)=ALF1(I,K)*ALF(I,K)
 C
-      END DO
-      END DO
+ 320  CONTINUE
+ 32   CONTINUE
 c     print*,'ap ADVYP 32'
 C
-      DO JV=1,NTRA
-      DO K=1,LAT-1
+      DO 33 JV=1,NTRA
+      DO 33 K=1,LAT-1
       KP=K+1
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)<0.) THEN
+      DO 330 I=1,LON
+C
+         IF(VGRI(I,K,L).LT.0.) THEN
 C
          TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
@@ -478,7 +474,6 @@
          ENDIF
 C
-      END DO
-      END DO
-      END DO
+ 330  CONTINUE
+ 33   CONTINUE
 c     print*,'ap ADVYP 33'
 C
@@ -488,11 +483,11 @@
 C
       SM0=0.
-      DO JV=1,NTRA
+      DO 40 JV=1,NTRA
          S00(JV)=0.
-      END DO
-C
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)>=0.) THEN
+ 40   CONTINUE
+C
+      DO 41 I=1,LON
+C
+         IF(VGRI(I,K,L).GE.0.) THEN
            FM(I,K)=VGRI(I,K,L)*DTY
            ALF(I,K)=FM(I,K)/SM(I,K,L)
@@ -508,11 +503,11 @@
          ALF4(I,K)=ALF1(I,K)*ALF1Q(I,K)
 C
-      END DO
+ 41   CONTINUE
 c     print*,'ap ADVYP 41'
 C
-      DO JV=1,NTRA
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)>=0.) THEN
+      DO 42 JV=1,NTRA
+      DO 420 I=1,LON
+C
+         IF(VGRI(I,K,L).GE.0.) THEN
            F0 (I,K,JV)=ALF(I,K)* ( S0(I,K,L,JV)+ALF1(I,K)*
      +             ( SY(I,K,L,JV)+ALF2(I,K)*SYY(I,K,L,JV) ) )
@@ -532,29 +527,29 @@
          ENDIF
 C
-      END DO
-      END DO
+ 420  CONTINUE
+ 42   CONTINUE
 c     print*,'ap ADVYP 42'
 C
-      DO I=1,LON
-         IF(VGRI(I,K,L)<0.) THEN
+      DO 43 I=1,LON
+         IF(VGRI(I,K,L).LT.0.) THEN
            FM(I,K)=-VGRI(I,K,L)*DTY
            ALF(I,K)=FM(I,K)/SM0
          ENDIF
-      END DO
+ 43   CONTINUE
 c     print*,'ap ADVYP 43'
 C
-      DO JV=1,NTRA
-      DO I=1,LON
-         IF(VGRI(I,K,L)<0.) THEN
+      DO 44 JV=1,NTRA
+      DO 440 I=1,LON
+         IF(VGRI(I,K,L).LT.0.) THEN
            F0(I,K,JV)=ALF(I,K)*S00(JV)
          ENDIF
-      END DO
-      END DO
+ 440  CONTINUE
+ 44   CONTINUE
 C
 C  puts the temporary moments Fi into appropriate neighboring boxes
 C
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)<0.) THEN
+      DO 45 I=1,LON
+C
+         IF(VGRI(I,K,L).LT.0.) THEN
            SM(I,K,L)=SM(I,K,L)+FM(I,K)
            ALF(I,K)=FM(I,K)/SM(I,K,L)
@@ -567,11 +562,11 @@
          ALF3(I,K)=ALF1(I,K)*ALF(I,K)
 C
-      END DO
+ 45   CONTINUE
 c     print*,'ap ADVYP 45'
 C
-      DO JV=1,NTRA
-      DO I=1,LON
-C
-         IF(VGRI(I,K,L)<0.) THEN
+      DO 46 JV=1,NTRA
+      DO 460 I=1,LON
+C
+         IF(VGRI(I,K,L).LT.0.) THEN
 C
          TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
@@ -585,9 +580,9 @@
          ENDIF
 C
-      END DO
-      END DO
+ 460  CONTINUE
+ 46   CONTINUE
 c     print*,'ap ADVYP 46'
 C
-      END DO
+ 1    CONTINUE
 
 c--------------------------------------------------
Index: LMDZ6/trunk/libf/dyn3d_common/advz.F
===================================================================
--- LMDZ6/trunk/libf/dyn3d_common/advz.F	(revision 5079)
+++ LMDZ6/trunk/libf/dyn3d_common/advz.F	(revision 5084)
@@ -117,7 +117,7 @@
 C  Conversion du flux de masse en kg.s-1
 
-      DO l = 1,llm
-         DO j = 1,jjp1
-            DO i = 1,iip1  
+      DO 500 l = 1,llm
+         DO 500 j = 1,jjp1
+            DO 500 i = 1,iip1  
 c            wgri (i,j,llm+1-l) =  w (i,j,l) / g 
                wgri (i,j,llm+1-l) =  w (i,j,l) 
@@ -125,7 +125,5 @@
 c             wgri (i,j,l) = 0.1               !    w (i,j,l) 
 c             wgri (i,j,llm) = 0.              ! a detruire ult.
-      END DO
-      END DO
-      END DO
+  500 CONTINUE
          DO  j = 1,jjp1
             DO i = 1,iip1  
@@ -139,5 +137,5 @@
 C  boucle sur les latitudes
 C
-      DO K=1,LAT
+      DO 1 K=1,LAT
 C
 C  place limits on appropriate moments before transport
@@ -146,12 +144,11 @@
       IF(.NOT.LIMIT) GO TO 101
 C
-      DO JV=1,NTRA
-      DO L=1,NIV
-         DO I=1,LON
+      DO 10 JV=1,NTRA
+      DO 10 L=1,NIV
+         DO 100 I=1,LON
             sz(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.),
      +                              ABS(sz(I,K,L,JV))),sz(I,K,L,JV))
-      END DO
-      END DO
-      END DO
+ 100     CONTINUE
+ 10   CONTINUE
 C
  101  CONTINUE
@@ -165,10 +162,10 @@
 C  2- reajusts moments remaining in the box
 C
-      DO L=1,NIV-1
+      DO 11 L=1,NIV-1
       LP=L+1
 C
-      DO I=1,LON
-C
-         IF(WGRI(I,K,L)<0.) THEN
+      DO 110 I=1,LON
+C
+         IF(WGRI(I,K,L).LT.0.) THEN
            FM(I,L)=-WGRI(I,K,L)*DTZ
            ALF(I)=FM(I,L)/SM(I,K,LP)
@@ -184,10 +181,10 @@
          ALF1Q(I)=ALF1(I)*ALF1(I)
 C
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LON
-C
-         IF(WGRI(I,K,L)<0.) THEN
+ 110  CONTINUE
+C
+      DO 111 JV=1,NTRA
+      DO 1110 I=1,LON
+C
+         IF(WGRI(I,K,L).LT.0.) THEN
 C
            F0(I,L,JV)=ALF (I)*( S0(I,K,LP,JV)-ALF1(I)*sz(I,K,LP,JV) )
@@ -215,17 +212,17 @@
          ENDIF
 C
-      END DO
-      END DO
-C
-      END DO
+ 1110 CONTINUE
+ 111  CONTINUE
+C
+ 11   CONTINUE
 C
 C  puts the temporary moments Fi into appropriate neighboring boxes
 C
-      DO L=1,NIV-1
+      DO 12 L=1,NIV-1
       LP=L+1
 C
-      DO I=1,LON
-C
-         IF(WGRI(I,K,L)<0.) THEN
+      DO 120 I=1,LON
+C
+         IF(WGRI(I,K,L).LT.0.) THEN
            SM(I,K,L)=SM(I,K,L)+FM(I,L)
            ALF(I)=FM(I,L)/SM(I,K,L)
@@ -239,10 +236,10 @@
          ALF1Q(I)=ALF1(I)*ALF1(I)
 C
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LON
-C
-         IF(WGRI(I,K,L)<0.) THEN
+ 120  CONTINUE
+C
+      DO 121 JV=1,NTRA
+      DO 1210 I=1,LON
+C
+         IF(WGRI(I,K,L).LT.0.) THEN
 C
            TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)
@@ -263,12 +260,12 @@
          ENDIF
 C
-      END DO
-      END DO
-C
-      END DO
+ 1210 CONTINUE
+ 121  CONTINUE
+C
+ 12   CONTINUE
 C
 C  fin de la boucle principale sur les latitudes
 C
-      END DO
+ 1    CONTINUE
 C
 C-------------------------------------------------------------
Index: LMDZ6/trunk/libf/dyn3d_common/advzp.F
===================================================================
--- LMDZ6/trunk/libf/dyn3d_common/advzp.F	(revision 5079)
+++ LMDZ6/trunk/libf/dyn3d_common/advzp.F	(revision 5084)
@@ -135,11 +135,9 @@
 C  Conversion des flux de masses en kg
 
-      DO l = 1,llm
-         DO j = 1,jjp1
-            DO i = 1,iip1  
+      DO 500 l = 1,llm
+         DO 500 j = 1,jjp1
+            DO 500 i = 1,iip1  
             wgri (i,j,llm+1-l) = w (i,j,l)  
-      END DO
-      END DO
-      END DO
+  500 CONTINUE
       do j=1,jjp1
          do i=1,iip1
@@ -156,5 +154,5 @@
 C  boucle sur les latitudes
 C
-      DO K=1,LAT
+      DO 1 K=1,LAT
 C
 C  place limits on appropriate moments before transport
@@ -163,8 +161,8 @@
       IF(.NOT.LIMIT) GO TO 101
 C
-      DO JV=1,NTRA
-      DO L=1,NIV
-         DO I=1,LON
-            IF(S0(I,K,L,JV)>0.) THEN
+      DO 10 JV=1,NTRA
+      DO 10 L=1,NIV
+         DO 100 I=1,LON
+            IF(S0(I,K,L,JV).GT.0.) THEN
               SLPMAX=S0(I,K,L,JV)
               S1MAX =1.5*SLPMAX
@@ -182,7 +180,6 @@
               SYZ(I,K,L,JV)=0.
             ENDIF
-      END DO
-      END DO
-      END DO
+ 100     CONTINUE
+ 10   CONTINUE
 C
  101  CONTINUE
@@ -196,10 +193,10 @@
 C  2- reajusts moments remaining in the box
 C
-      DO L=1,NIV-1
+      DO 11 L=1,NIV-1
       LP=L+1
 C
-      DO I=1,LON
-C
-         IF(WGRI(I,K,L)<0.) THEN
+      DO 110 I=1,LON
+C
+         IF(WGRI(I,K,L).LT.0.) THEN
            FM(I,L)=-WGRI(I,K,L)*DTZ
            ALF(I)=FM(I,L)/SM(I,K,LP)
@@ -218,10 +215,10 @@
          ALF4 (I)=ALF1(I)*ALF1Q(I)
 C
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LON
-C
-         IF(WGRI(I,K,L)<0.) THEN
+ 110  CONTINUE
+C
+      DO 111 JV=1,NTRA
+      DO 1110 I=1,LON
+C
+         IF(WGRI(I,K,L).LT.0.) THEN
 C
            F0 (I,L,JV)=ALF (I)* ( S0(I,K,LP,JV)-ALF1(I)*
@@ -276,17 +273,17 @@
          ENDIF
 C
-      END DO
-      END DO
-C
-      END DO
+ 1110 CONTINUE
+ 111  CONTINUE
+C
+ 11   CONTINUE
 C
 C  puts the temporary moments Fi into appropriate neighboring boxes
 C
-      DO L=1,NIV-1
+      DO 12 L=1,NIV-1
       LP=L+1
 C
-      DO I=1,LON
-C
-         IF(WGRI(I,K,L)<0.) THEN
+      DO 120 I=1,LON
+C
+         IF(WGRI(I,K,L).LT.0.) THEN
            SM(I,K,L)=SM(I,K,L)+FM(I,L)
            ALF(I)=FM(I,L)/SM(I,K,L)
@@ -302,10 +299,10 @@
          ALF3(I)=ALF1(I)-ALF(I)
 C
-      END DO
-C
-      DO JV=1,NTRA
-      DO I=1,LON
-C
-         IF(WGRI(I,K,L)<0.) THEN
+ 120  CONTINUE
+C
+      DO 121 JV=1,NTRA
+      DO 1210 I=1,LON
+C
+         IF(WGRI(I,K,L).LT.0.) THEN
 C
            TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)
@@ -345,12 +342,12 @@
          ENDIF
 C
-      END DO
-      END DO
-C
-      END DO
+ 1210 CONTINUE
+ 121  CONTINUE
+C
+ 12   CONTINUE
 C
 C  fin de la boucle principale sur les latitudes
 C
-      END DO
+ 1    CONTINUE
 C
       DO l = 1,llm
Index: LMDZ6/trunk/libf/dyn3d_common/comdissip.h
===================================================================
--- LMDZ6/trunk/libf/dyn3d_common/comdissip.h	(revision 5079)
+++ LMDZ6/trunk/libf/dyn3d_common/comdissip.h	(revision 5084)
@@ -6,5 +6,5 @@
 
       COMMON/comdissip/                                                 &
-     &    coefdis,tetavel,tetatemp,gamdissip,niterdis
+     &    niterdis,coefdis,tetavel,tetatemp,gamdissip
 
 
Index: LMDZ6/trunk/libf/dyn3d_common/extrapol.F
===================================================================
--- LMDZ6/trunk/libf/dyn3d_common/extrapol.F	(revision 5079)
+++ LMDZ6/trunk/libf/dyn3d_common/extrapol.F	(revision 5084)
@@ -58,7 +58,7 @@
 200   CONTINUE
       incre = incre + 1
-      DO j = 1, kylat
-      DO i = 1, kxlon
-      IF (pfild(i,j)> zwmsk) THEN
+      DO 99999 j = 1, kylat
+      DO 99999 i = 1, kxlon
+      IF (pfild(i,j).GT. zwmsk) THEN
          pwork(i,j) = pfild(i,j)
          inbor = 0
@@ -89,15 +89,15 @@
 C
 C* Correct latitude bounds if southernmost or northernmost points
-         IF (j == 1) ideb = 4
-         IF (j == kylat) ifin = 6
+         IF (j .EQ. 1) ideb = 4
+         IF (j .EQ. kylat) ifin = 6
 C
 C* Account for periodicity in longitude
 C
          IF (ldper) THEN 
-            IF (i == kxlon) THEN
+            IF (i .EQ. kxlon) THEN
                ix(3) = 1
                ix(6) = 1
                ix(9) = 1
-            ELSE IF (i == 1) THEN
+            ELSE IF (i .EQ. 1) THEN
                ix(1) = kxlon
                ix(4) = kxlon
@@ -105,5 +105,5 @@
             ENDIF
          ELSE
-            IF (i == 1) THEN
+            IF (i .EQ. 1) THEN
                ix(1) = i
                ix(2) = i + 1
@@ -113,5 +113,5 @@
                ix(6) = i + 1
             ENDIF 
-            IF (i == kxlon) THEN
+            IF (i .EQ. kxlon) THEN
                ix(1) = i -1
                ix(2) = i
@@ -122,5 +122,5 @@
             ENDIF
 C
-            IF (i == 1 .OR. i == kxlon) THEN 
+            IF (i .EQ. 1 .OR. i .EQ. kxlon) THEN 
                jy(1) = MAX (1,j-1)
                jy(2) = MAX (1,j-1)
@@ -132,6 +132,6 @@
                ideb = 1
                ifin = 6
-               IF (j == 1) ideb = 3
-               IF (j == kylat) ifin = 4
+               IF (j .EQ. 1) ideb = 3
+               IF (j .EQ. kylat) ifin = 4
             ENDIF
          ENDIF ! end for ldper test
@@ -139,18 +139,18 @@
 C* Find unmasked neighbors
 C
-         DO k = ideb, ifin
+         DO 230 k = ideb, ifin
             zmask(k) = 0.
             ilon = ix(k)
             jlat = jy(k)
-            IF (pfild(ilon,jlat) < zwmsk) THEN
+            IF (pfild(ilon,jlat) .LT. zwmsk) THEN
                zmask(k) = 1.
                inbor = inbor + 1
             ENDIF
-      END DO
+ 230     CONTINUE
 C
 C* Not enough points around point P are unmasked; interpolation on P 
 C  will be done in a future call to extrap.
 C
-         IF (inbor >= knbor) THEN
+         IF (inbor .GE. knbor) THEN
             pwork(i,j) = 0.
             DO k = ideb, ifin
@@ -163,6 +163,5 @@
 C
       ENDIF
-      END DO
-      END DO
+99999 CONTINUE
 C
 C*    3. Writing back unmasked field in pfild
@@ -177,10 +176,10 @@
       DO j = 1, kylat
       DO i = 1, kxlon
-         IF (pwork(i,j) > zwmsk) idoit = idoit + 1
+         IF (pwork(i,j) .GT. zwmsk) idoit = idoit + 1
          pfild(i,j) = pwork(i,j)
       ENDDO
       ENDDO
 c
-      IF (idoit /= 0) GOTO 200
+      IF (idoit .ne. 0) GOTO 200
 ccc      PRINT*, "Number of extrapolation steps incre =", incre
 c
Index: LMDZ6/trunk/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90
===================================================================
--- LMDZ6/trunk/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90	(revision 5079)
+++ LMDZ6/trunk/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90	(revision 5084)
@@ -14,6 +14,5 @@
   USE comconst_mod, ONLY: cpp, kappa, g, omeg, daysec, rad, pi
   USE comvert_mod, ONLY: presnivs, preff, pa
-  USE lmdz_netcdf, ONLY: nf90_def_var, nf90_int, nf90_float, nf90_put_var, nf_enddef, &
-      nf_put_att_text,nf_def_dim,nf_64bit_offset,nf_clobber,nf_create
+  use netcdf, only: nf90_def_var, nf90_int, nf90_float, nf90_put_var
   
   IMPLICIT NONE
@@ -22,4 +21,5 @@
   INCLUDE "paramet.h"
   INCLUDE "comgeom.h"
+  INCLUDE "netcdf.inc"
 
 !======================== 
@@ -232,8 +232,8 @@
 
 SUBROUTINE handle_err(status)
-  USE lmdz_netcdf, ONLY: nf_strerror
+  INCLUDE "netcdf.inc"
 
   INTEGER status
-  IF (status/=nf_noerr) THEN
+  IF (status.NE.nf_noerr) THEN
      PRINT *,NF_STRERROR(status)
      CALL abort_gcm('grilles_gcm_netcdf','netcdf error',1)
Index: LMDZ6/trunk/libf/dyn3d_common/ppm3d.F
===================================================================
--- LMDZ6/trunk/libf/dyn3d_common/ppm3d.F	(revision 5079)
+++ LMDZ6/trunk/libf/dyn3d_common/ppm3d.F	(revision 5084)
@@ -68,5 +68,5 @@
       implicit none
 
-c     rajout de d�clarations
+c     rajout de déclarations
 c      integer Jmax,kmax,ndt0,nstep,k,j,i,ic,l,js,jn,imh,iad,jad,krd
 c      integer iu,iiu,j2,jmr,js0,jt
@@ -315,5 +315,5 @@
 C
 C *********** Initialization **********************
-      if(NSTEP==1) then
+      if(NSTEP.eq.1) then
 c
       write(6,*) '------------------------------------ '
@@ -325,14 +325,14 @@
 C
 C controles sur les parametres
-      if(NLAY<6) then
+      if(NLAY.LT.6) then
         write(6,*) 'NLAY must be >= 6'
         stop
       endif
-      if (JNP<NLAY) then
+      if (JNP.LT.NLAY) then
          write(6,*) 'JNP must be >= NLAY'
         stop
       endif
       IMRD2=mod(IMR,2)
-      if (j1==2.and.IMRD2/=0) then
+      if (j1.eq.2.and.IMRD2.NE.0) then
          write(6,*) 'if j1=2 IMR must be an even integer'
         stop
@@ -340,5 +340,5 @@
 
 C
-      if(Jmax<JNP .or. Kmax<NLAY) then
+      if(Jmax.lt.JNP .or. Kmax.lt.NLAY) then
         write(6,*) 'Jmax or Kmax is too small'
         stop
@@ -354,5 +354,5 @@
       DP =    PI / REAL(JMR)
 C
-      if(IGD==0) then
+      if(IGD.eq.0) then
 C Compute analytic cosine at cell edges
             call cosa(cosp,cose,JNP,PI,DP)
@@ -362,7 +362,6 @@
       endif
 C
-      do J=2,JMR
-      acosp(j) = 1. / cosp(j)
-      END DO
+      do 15 J=2,JMR
+15    acosp(j) = 1. / cosp(j)
 C
 C Inverse of the Scaled polar cap area.
@@ -373,9 +372,9 @@
       endif
 C
-      if(NDT0 /= NDT) then
+      if(NDT0 .ne. NDT) then
       DT   = NDT
       NDT0 = NDT
 
-	if(Umax < 180.) then
+	if(Umax .lt. 180.) then
          write(6,*) 'Umax may be too small!'
 	endif
@@ -383,9 +382,9 @@
       MaxDT = DP*AE / abs(Umax) + 0.5
       write(6,*)'Largest time step for max(V)=',Umax,' is ',MaxDT
-      if(MaxDT < abs(NDT)) then
+      if(MaxDT .lt. abs(NDT)) then
             write(6,*) 'Warning!!! NDT maybe too large!'
       endif
 C
-      if(CR1>=0.95) then
+      if(CR1.ge.0.95) then
       JS0 = 0
       JN0 = 0
@@ -430,47 +429,37 @@
           
 C
-      if(j1/=2) then
-      DO IC=1,NC
-      DO L=1,NLAY
-      DO I=1,IMR
+      if(j1.ne.2) then
+      DO 40 IC=1,NC
+      DO 40 L=1,NLAY
+      DO 40 I=1,IMR
       Q(I,  2,L,IC) = Q(I,  1,L,IC)
-      Q(I,JMR,L,IC) = Q(I,JNP,L,IC)
-      END DO
-      END DO
-      END DO
+40    Q(I,JMR,L,IC) = Q(I,JNP,L,IC)
       endif
 C
 C Compute "tracer density"
-      DO IC=1,NC
-      DO k=1,NLAY
-      DO j=1,JNP
-      DO i=1,IMR
-      DQ(i,j,k,IC) = Q(i,j,k,IC)*delp1(i,j,k)
-      END DO
-      END DO
-      END DO
-      END DO
-C
-      do k=1,NLAY
-C
-      if(IGD==0) then
+      DO 550 IC=1,NC
+      DO 44 k=1,NLAY
+      DO 44 j=1,JNP
+      DO 44 i=1,IMR
+44    DQ(i,j,k,IC) = Q(i,j,k,IC)*delp1(i,j,k)
+550	continue
+C
+      do 1500 k=1,NLAY
+C
+      if(IGD.eq.0) then
 C 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)
       else
 C Convert winds on C-grid to Courant #
-      do j=j1,j2
-      do i=2,IMR
-      CRX(i,J) = dtdx(j)*U(i-1,j,k)
-      END DO
-      END DO
+      do 45 j=j1,j2
+      do 45 i=2,IMR
+45    CRX(i,J) = dtdx(j)*U(i-1,j,k)
    
 C
-      do j=j1,j2
-      CRX(1,J) = dtdx(j)*U(IMR,j,k)
-      END DO
-C
-      do i=1,IMR*JMR
-      CRY(i,2) = DTDY*V(i,1,k)
-      END DO
+      do 50 j=j1,j2
+50    CRX(1,J) = dtdx(j)*U(IMR,j,k)
+C
+      do 55 i=1,IMR*JMR
+55    CRY(i,2) = DTDY*V(i,1,k)
       endif
 C     
@@ -481,5 +470,5 @@
       do j=JS0,j1+1,-1
       do i=1,IMR
-      if(abs(CRX(i,j))>1.) then
+      if(abs(CRX(i,j)).GT.1.) then
             JS = j
             go to 2222
@@ -491,5 +480,5 @@
       do j=JN0,j2-1
       do i=1,IMR
-      if(abs(CRX(i,j))>1.) then
+      if(abs(CRX(i,j)).GT.1.) then
             JN = j
             go to 2233
@@ -499,5 +488,5 @@
 2233  continue
 C
-      if(j1/=2) then           ! Enlarged polar cap.
+      if(j1.ne.2) then           ! Enlarged polar cap.
       do i=1,IMR
       DPI(i,  2,k) = 0.
@@ -516,9 +505,7 @@
       enddo
 C
-      do j=j1,j2
-      DO i=1,IMR
-      DPI(i,j,k) = (ymass(i,j) - ymass(i,j+1)) * acosp(j)
-      END DO
-      END DO
+      do 95 j=j1,j2
+      DO 95 i=1,IMR
+95    DPI(i,j,k) = (ymass(i,j) - ymass(i,j+1)) * acosp(j)
 C
 C Poles
@@ -549,19 +536,14 @@
       enddo
 C
-      do j=j1,j2
-      DO i=1,IMR
-      xmass(i,j) = PU(i,j)*CRX(i,j)
-      END DO
-      END DO
-C
-      DO j=j1,j2
-      DO i=1,IMR-1
-      DPI(i,j,k) = DPI(i,j,k) + xmass(i,j) - xmass(i+1,j)
-      END DO
-      END DO
-C
-      DO j=j1,j2
-      DPI(IMR,j,k) = DPI(IMR,j,k) + xmass(IMR,j) - xmass(1,j)
-      END DO
+      do 110 j=j1,j2
+      DO 110 i=1,IMR
+110   xmass(i,j) = PU(i,j)*CRX(i,j)
+C
+      DO 120 j=j1,j2
+      DO 120 i=1,IMR-1
+120   DPI(i,j,k) = DPI(i,j,k) + xmass(i,j) - xmass(i+1,j)
+C
+      DO 130 j=j1,j2
+130   DPI(IMR,j,k) = DPI(IMR,j,k) + xmass(IMR,j) - xmass(1,j)
 C
       DO j=j1,j2
@@ -587,5 +569,5 @@
       enddo
 C
-      if(j1==2) then
+      if(j1.eq.2) then
 	IMH = IMR/2
       do i=1,IMH
@@ -600,5 +582,5 @@
 C
 C ****6***0*********0*********0*********0*********0*********0**********72
-      do IC=1,NC
+      do 1000 IC=1,NC
 C
       do i=1,IMJM
@@ -608,6 +590,6 @@
 C
 C E-W advective cross term
-      do j=J1,J2
-      if(J>JS  .and. J<JN) GO TO 250
+      do 250 j=J1,J2
+      if(J.GT.JS  .and. J.LT.JN) GO TO 250
 C
       do i=1,IMR
@@ -620,9 +602,9 @@
       enddo
 C
-      DO i=1,IMR
+      DO 230 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).GE.0.) then
       wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))
       else
@@ -630,9 +612,8 @@
       endif
       wk1(i,j,1) = wk1(i,j,1) - qtmp(i)
-      END DO
+230   continue
 250   continue
-      END DO
-C
-      if(JN/=0) then
+C
+      if(JN.ne.0) then
       do j=JS+1,JN-1
 C
@@ -664,5 +645,5 @@
 	if(cross) then
 C Add cross terms in the vertical direction.
-	if(IORD >= 2) then
+	if(IORD .GE. 2) then
 		iad = 2
 	else
@@ -670,5 +651,5 @@
 	endif
 C
-	if(JORD >= 2) then
+	if(JORD .GE. 2) then
 		jad = 2
 	else
@@ -690,6 +671,6 @@
      &  DC2,ymass,WK1(1,1,3),wk1(1,1,4),WK1(1,1,5),WK1(1,1,6),JORD)
 C
-      END DO
-      END DO
+1000  continue
+1500  continue
 C
 C ******* Compute vertical mass flux (same unit as PS) ***********
@@ -697,20 +678,16 @@
 C 1st step: compute total column mass CONVERGENCE.
 C
-      do j=1,JNP
-      do i=1,IMR
-      CRY(i,j) = DPI(i,j,1)
-      END DO
-      END DO
-C
-      do k=2,NLAY
-      do j=1,JNP
-      do i=1,IMR
+      do 320 j=1,JNP
+      do 320 i=1,IMR
+320   CRY(i,j) = DPI(i,j,1)
+C
+      do 330 k=2,NLAY
+      do 330 j=1,JNP
+      do 330 i=1,IMR
       CRY(i,j)  = CRY(i,j) + DPI(i,j,k)
-      END DO
-      END DO
-      END DO
-C
-      do j=1,JNP
-      do i=1,IMR
+330   continue
+C
+      do 360 j=1,JNP
+      do 360 i=1,IMR
 C
 C 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.
@@ -723,25 +700,20 @@
       W(i,j,1) = DPI(i,j,1) - DBK(1)*CRY(i,j)
       W(i,j,NLAY) = 0.
-      END DO
-      END DO
-C
-      do k=2,NLAY-1
-      do j=1,JNP
-      do i=1,IMR
+360   continue
+C
+      do 370 k=2,NLAY-1
+      do 370 j=1,JNP
+      do 370 i=1,IMR
       W(i,j,k) = W(i,j,k-1) + DPI(i,j,k) - DBK(k)*CRY(i,j)
-      END DO
-      END DO
-      END DO
-C
-      DO k=1,NLAY
-      DO j=1,JNP
-      DO i=1,IMR
+370   continue
+C
+      DO 380 k=1,NLAY
+      DO 380 j=1,JNP
+      DO 380 i=1,IMR
       delp2(i,j,k) = DAP(k) + DBK(k)*PS2(i,j)
-      END DO
-      END DO
-      END DO
+380   continue
 C
 	KRD = max(3, KORD)
-      do IC=1,NC
+      do 4000 IC=1,NC
 C
 C****6***0*********0*********0*********0*********0*********0**********72
@@ -766,22 +738,20 @@
       enddo
 C     
-      if(j1/=2) then
-      DO k=1,NLAY
-      DO I=1,IMR
-c     j=1 c'est le p�le Sud, j=JNP c'est le p�le Nord
+      if(j1.ne.2) then
+      DO 400 k=1,NLAY
+      DO 400 I=1,IMR
+c     j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord
       Q(I,  2,k,IC) = Q(I,  1,k,IC)
       Q(I,JMR,k,IC) = Q(I,JNP,k,IC)
-      END DO
-      END DO
-      endif
-      END DO
-C
-      if(j1/=2) then
-      DO k=1,NLAY
-      DO i=1,IMR
+400   CONTINUE
+      endif
+4000  continue
+C
+      if(j1.ne.2) then
+      DO 5000 k=1,NLAY
+      DO 5000 i=1,IMR
       W(i,  2,k) = W(i,  1,k)
       W(i,JMR,k) = W(i,JNP,k)
-      END DO
-      END DO
+5000  continue
       endif
 C
@@ -815,12 +785,11 @@
 C ****6***0*********0*********0*********0*********0*********0**********72
 C
-      do k=1,NLAYM1
-      do i=1,IMJM
+      do 1000 k=1,NLAYM1
+      do 1000 i=1,IMJM
       DQDT(i,1,k) = P(i,1,k+1) - P(i,1,k)
-      END DO
-      END DO
-C
-      DO k=2,NLAYM1
-      DO I=1,IMJM
+1000  continue
+C
+      DO 1220 k=2,NLAYM1
+      DO 1220 I=1,IMJM    
        c0 =  delp(i,1,k) / (delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1))
        c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k))    
@@ -830,6 +799,5 @@
       Qmin = P(i,1,k) - min(P(i,1,k-1),P(i,1,k),P(i,1,k+1))
       DC(i,1,k) = sign(min(abs(tmp),Qmax,Qmin), tmp)   
-      END DO
-      END DO
+1220  CONTINUE
      
 C     
@@ -838,6 +806,6 @@
 C ****6***0*********0*********0*********0*********0*********0**********72
 C
-      DO j=1,JNP
-      if((j==2 .or. j==JMR) .and. j1/=2) goto 2000
+      DO 2000 j=1,JNP
+      if((j.eq.2 .or. j.eq.JMR) .and. j1.ne.2) goto 2000
 C
       DO k=1,NLAY
@@ -860,5 +828,5 @@
 C
 C First guess top edge value
-      DO i=1,IMR
+      DO 10 i=1,IMR
 C three-cell PPM
 C Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp
@@ -872,5 +840,5 @@
 C
 C Check if change sign
-      if(wk1(i,1)*AL(i,1)<=0.) then
+      if(wk1(i,1)*AL(i,1).le.0.) then
 		 AL(i,1) = 0.
              flux(i,1) = 0.
@@ -878,8 +846,8 @@
              flux(i,1) =  wk1(i,1) - AL(i,1)
 	endif
-      END DO
+10    continue
 C
 C Bottom
-      DO i=1,IMR
+      DO 15 i=1,IMR
 C 2-cell PPM with zero gradient right at the surface
 C
@@ -888,7 +856,7 @@
       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).le.0.) AR(i,NLAY) = 0.
       flux(i,NLAY) = AR(i,NLAY) -  wk1(i,NLAY)
-      END DO
+15    continue
      
 C
@@ -897,6 +865,6 @@
 C****6***0*********0*********0*********0*********0*********0**********72
 C
-      DO k=3,NLAYM1
-      DO i=1,IMR
+      DO 14 k=3,NLAYM1
+      DO 12 i=1,IMR
       c1 =  DQDT(i,j,k-1)*wk2(i,k-1) / (wk2(i,k-1)+wk2(i,k))
       c2 =  2. / (wk2(i,k-2)+wk2(i,k-1)+wk2(i,k)+wk2(i,k+1))
@@ -907,16 +875,16 @@
      &          wk2(i,k-1)*A1*flux(i,k)  )
 C      print *,'AL1',i,k, AL(i,k)
-      END DO
-      END DO
-C
-      do i=1,IMR*NLAYM1
+12    CONTINUE
+14    continue
+C
+      do 20 i=1,IMR*NLAYM1
       AR(i,1) = AL(i,2)
 C      print *,'AR1',i,AR(i,1)
-      END DO
-C
-      do i=1,IMR*NLAY
+20    continue
+C
+      do 30 i=1,IMR*NLAY
       A6(i,1) = 3.*(wk1(i,1)+wk1(i,1) - (AL(i,1)+AR(i,1)))
 C      print *,'A61',i,A6(i,1)
-      END DO
+30    continue
 C
 C****6***0*********0*********0*********0*********0*********0**********72
@@ -927,5 +895,5 @@
 C
 C Interior depending on KORD
-      if(LMT<=2)
+      if(LMT.LE.2)
      &  call lmtppm(flux(1,2),A6(1,2),AR(1,2),AL(1,2),wk1(1,2),
      &              IMR*(NLAY-2),LMT)
@@ -933,6 +901,6 @@
 C****6***0*********0*********0*********0*********0*********0**********72
 C
-      DO i=1,IMR*NLAYM1
-      IF(wz2(i,1)>0.) then
+      DO 140 i=1,IMR*NLAYM1
+      IF(wz2(i,1).GT.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))
@@ -944,22 +912,19 @@
 C        print *,'test2',i, AL(i,2),AR(i,2),A6(i,2),R23
       endif
-      END DO
-C
-      DO i=1,IMR*NLAYM1
+140   continue
+C
+      DO 250 i=1,IMR*NLAYM1
       flux(i,2) = wz2(i,1) * flux(i,2)
-      END DO
-C
-      do i=1,IMR
+250   continue
+C
+      do 350 i=1,IMR
       DQ(i,j,   1) = DQ(i,j,   1) - flux(i,   2)
       DQ(i,j,NLAY) = DQ(i,j,NLAY) + flux(i,NLAY)
-      END DO
-C
-      do k=2,NLAYM1
-      do i=1,IMR
-      DQ(i,j,k) = DQ(i,j,k) + flux(i,k) - flux(i,k+1)
-      END DO
-      END DO
+350   continue
+C
+      do 360 k=2,NLAYM1
+      do 360 i=1,IMR
+360   DQ(i,j,k) = DQ(i,j,k) + flux(i,k) - flux(i,k+1)
 2000  continue
-      END DO
       return
       end
@@ -985,5 +950,5 @@
       j2vl = j2-jvan
 C
-      do j=j1,j2
+      do 1310 j=j1,j2
 C
       do i=1,IMR
@@ -991,5 +956,5 @@
       enddo
 C
-      if(j>=JN .or. j<=JS) goto 2222
+      if(j.ge.JN .or. j.le.JS) goto 2222
 C ************* Eulerian **********
 C
@@ -999,18 +964,16 @@
       qtmp(IMP+1) = q(2,J)
 C
-      IF(IORD==1 .or. j==j1 .or. j==j2) THEN
-      DO i=1,IMR
+      IF(IORD.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN
+      DO 1406 i=1,IMR
       iu = REAL(i) - uc(i,j)
-      fx1(i) = qtmp(iu)
-      END DO
+1406  fx1(i) = qtmp(iu)
       ELSE
       call xmist(IMR,IML,Qtmp,DC)
       DC(0) = DC(IMR)
 C
-      if(IORD==2 .or. j<=j1vl .or. j>=j2vl) then
-      DO i=1,IMR
+      if(IORD.eq.2 .or. j.le.j1vl .or. j.ge.j2vl) then
+      DO 1408 i=1,IMR
       iu = REAL(i) - uc(i,j)
-      fx1(i) = qtmp(iu) + DC(iu)*(sign(1.,uc(i,j))-uc(i,j))
-      END DO
+1408  fx1(i) = qtmp(iu) + DC(iu)*(sign(1.,uc(i,j))-uc(i,j))
       else
       call fxppm(IMR,IML,UC(1,j),Qtmp,DC,fx1,IORD)
@@ -1019,7 +982,6 @@
       ENDIF
 C
-      DO i=1,IMR
-      fx1(i) = fx1(i)*xmass(i,j)
-      END DO
+      DO 1506 i=1,IMR
+1506  fx1(i) = fx1(i)*xmass(i,j)
 C
       goto 1309
@@ -1034,11 +996,10 @@
       enddo
 C
-      IF(IORD==1 .or. j==j1 .or. j==j2) THEN
-      DO i=1,IMR
+      IF(IORD.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN
+      DO 1306 i=1,IMR
       itmp = INT(uc(i,j))
       ISAVE(i) = i - itmp
       iu = i - uc(i,j)
-      fx1(i) = (uc(i,j) - itmp)*qtmp(iu)
-      END DO
+1306  fx1(i) = (uc(i,j) - itmp)*qtmp(iu)
       ELSE
       call xmist(IMR,IML,Qtmp,DC)
@@ -1049,20 +1010,19 @@
       enddo
 C
-      DO i=1,IMR
+      DO 1307 i=1,IMR
       itmp = INT(uc(i,j))
       rut  = uc(i,j) - itmp
       ISAVE(i) = i - itmp
       iu = i - uc(i,j)
-      fx1(i) = rut*(qtmp(iu) + DC(iu)*(sign(1.,rut) - rut))
-      END DO
+1307  fx1(i) = rut*(qtmp(iu) + DC(iu)*(sign(1.,rut) - rut))
       ENDIF
 C
-      do i=1,IMR
-      IF(uc(i,j)>1.) then
+      do 1308 i=1,IMR
+      IF(uc(i,j).GT.1.) then
 CDIR$ NOVECTOR
         do ist = ISAVE(i),i-1
         fx1(i) = fx1(i) + qtmp(ist)
         enddo
-      elseIF(uc(i,j)<-1.) then
+      elseIF(uc(i,j).LT.-1.) then
         do ist = i,ISAVE(i)-1
         fx1(i) = fx1(i) - qtmp(ist)
@@ -1070,5 +1030,5 @@
 CDIR$ VECTOR
       endif
-      END DO
+1308  continue
       do i=1,IMR
       fx1(i) = PU(i,j)*fx1(i)
@@ -1078,11 +1038,10 @@
 C
 1309  fx1(IMP) = fx1(1)
-      DO i=1,IMR
-      DQ(i,j) =  DQ(i,j) + fx1(i)-fx1(i+1)
-      END DO
+      DO 1215 i=1,IMR
+1215  DQ(i,j) =  DQ(i,j) + fx1(i)-fx1(i+1)
 C
 C ***************************************
 C
-      END DO
+1310  continue
       return
       end
@@ -1120,18 +1079,15 @@
 C      endif
 C
-      DO i=1,IMR
-      AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3
-      END DO
-C
-      do i=1,IMR-1
-      AR(i) = AL(i+1)
-      END DO
+      DO 10 i=1,IMR
+10    AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3
+C
+      do 20 i=1,IMR-1
+20    AR(i) = AL(i+1)
       AR(IMR) = AL(1)
 C
-      do i=1,IMR
-      A6(i) = 3.*(p(i)+p(i)  - (AL(i)+AR(i)))
-      END DO
-C
-      if(LMT<=2) call lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT)
+      do 30 i=1,IMR
+30    A6(i) = 3.*(p(i)+p(i)  - (AL(i)+AR(i)))
+C
+      if(LMT.LE.2) call lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT)
 C
       AL(0) = AL(IMR)
@@ -1140,5 +1096,5 @@
 C
       DO i=1,IMR
-      IF(UT(i)>0.) then
+      IF(UT(i).GT.0.) then
       flux(i) = AR(i-1) + 0.5*UT(i)*(AL(i-1) - AR(i-1) +
      &                 A6(i-1)*(1.-R23*UT(i)) )
@@ -1159,10 +1115,9 @@
       real :: tmp,pmax,pmin
 C
-      do i=1,IMR
+      do 10  i=1,IMR
       tmp = R24*(8.*(p(i+1) - p(i-1)) + p(i-2) - p(i+2))
       Pmax = max(P(i-1), p(i), p(i+1)) - p(i)
       Pmin = p(i) - min(P(i-1), p(i), p(i+1))
-      DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp)
-      END DO
+10    DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp)
       return
       end
@@ -1183,34 +1138,29 @@
       len = IMR*(J2-J1+2)
 C
-      if(JORD==1) then
-      DO i=1,len
+      if(JORD.eq.1) then
+      DO 1000 i=1,len
       JT = REAL(J1) - VC(i,J1)
-      fx(i,j1) = p(i,JT)
-      END DO
+1000  fx(i,j1) = p(i,JT)
       else
    
       call ymist(IMR,JNP,j1,P,DC2,4)
 C
-      if(JORD<=0 .or. JORD>=3) then
+      if(JORD.LE.0 .or. JORD.GE.3) then
    
       call fyppm(VC,P,DC2,fx,IMR,JNP,j1,j2,A6,AR,AL,JORD)
     
       else
-      DO i=1,len
+      DO 1200 i=1,len
       JT = REAL(J1) - VC(i,J1)
-      fx(i,j1) = p(i,JT) + (sign(1.,VC(i,j1))-VC(i,j1))*DC2(i,JT)
-      END DO
-      endif
-      endif
-C
-      DO i=1,len
-      fx(i,j1) = fx(i,j1)*ymass(i,j1)
-      END DO
-C
-      DO j=j1,j2
-      DO i=1,IMR
-      DQ(i,j) = DQ(i,j) + (fx(i,j) - fx(i,j+1)) * acosp(j)
-      END DO
-      END DO
+1200  fx(i,j1) = p(i,JT) + (sign(1.,VC(i,j1))-VC(i,j1))*DC2(i,JT)
+      endif
+      endif
+C
+      DO 1300 i=1,len
+1300  fx(i,j1) = fx(i,j1)*ymass(i,j1)
+C
+      DO 1400 j=j1,j2
+      DO 1400 i=1,IMR
+1400  DQ(i,j) = DQ(i,j) + (fx(i,j) - fx(i,j+1)) * acosp(j)
 C
 C Poles
@@ -1229,5 +1179,5 @@
       enddo
 C
-      if(j1/=2) then
+      if(j1.ne.2) then
       do i=1,IMR
       DQ(i,  2) = sum1
@@ -1251,13 +1201,13 @@
       IJM3 = IMR*(JMR-3)
 C
-      IF(ID==2) THEN
-      do i=1,IMR*(JMR-1)
+      IF(ID.EQ.2) THEN
+      do 10 i=1,IMR*(JMR-1)
       tmp = 0.25*(p(i,3) - p(i,1))
       Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2)
       Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3))
       DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp)
-      END DO
+10    CONTINUE
       ELSE
-      do i=1,IMH
+      do 12 i=1,IMH
 C J=2
       tmp = (8.*(p(i,3) - p(i,1)) + p(i+IMH,2) - p(i,4))*R24
@@ -1270,6 +1220,6 @@
       Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP))
       DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp)
-      END DO
-      do i=IMH+1,IMR
+12    CONTINUE
+      do 14 i=IMH+1,IMR
 C J=2
       tmp = (8.*(p(i,3) - p(i,1)) + p(i-IMH,2) - p(i,4))*R24
@@ -1282,15 +1232,15 @@
       Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP))
       DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp)
-      END DO
-C
-      do i=1,IJM3
+14    CONTINUE
+C
+      do 15 i=1,IJM3
       tmp = (8.*(p(i,4) - p(i,2)) + p(i,1) - p(i,5))*R24
       Pmax = max(p(i,2),p(i,3),p(i,4)) - p(i,3)
       Pmin = p(i,3) - min(p(i,2),p(i,3),p(i,4))
       DC(i,3) = sign(min(abs(tmp),Pmin,Pmax),tmp)
-      END DO
+15    CONTINUE
       ENDIF
 C
-      if(j1/=2) then
+      if(j1.ne.2) then
       do i=1,IMR
       DC(i,1) = 0.
@@ -1300,5 +1250,5 @@
 C Determine slopes in polar caps for scalars!
 C
-      do i=1,IMH
+      do 13 i=1,IMH
 C South
       tmp = 0.25*(p(i,2) - p(i+imh,2))
@@ -1311,10 +1261,10 @@
       Pmin = p(i,JNP) - min(p(i+imh,JMR),p(i,jnp), p(i,JMR))
       DC(i,JNP) = sign(min(abs(tmp),Pmax,pmin),tmp)
-      END DO
-C
-      do i=imh+1,IMR
+13    continue
+C
+      do 25 i=imh+1,IMR
       DC(i,  1) =  - DC(i-imh,  1)
       DC(i,JNP) =  - DC(i-imh,JNP)
-      END DO
+25    continue
       endif
       return
@@ -1358,8 +1308,8 @@
       LMT = JORD - 3      
 C
-      DO i=1,IMR*JMR
+      DO 10 i=1,IMR*JMR        
       AL(i,2) = 0.5*(p(i,1)+p(i,2)) + (DC(i,1) - DC(i,2))*R3
       AR(i,1) = AL(i,2)
-      END DO
+10    CONTINUE
 C
 CPoles:
@@ -1381,14 +1331,13 @@
       
            
-      do i=1,len
-      A6(i,j11) = 3.*(p(i,j11)+p(i,j11)  - (AL(i,j11)+AR(i,j11)))
-      END DO
-C
-      if(LMT<=2) call lmtppm(DC(1,j11),A6(1,j11),AR(1,j11)
+      do 30 i=1,len
+30    A6(i,j11) = 3.*(p(i,j11)+p(i,j11)  - (AL(i,j11)+AR(i,j11)))
+C
+      if(LMT.le.2) call lmtppm(DC(1,j11),A6(1,j11),AR(1,j11)
      &                       ,AL(1,j11),P(1,j11),len,LMT)
 C
      
-      DO i=1,IMJM1
-      IF(VC(i,j1)>0.) then
+      DO 140 i=1,IMJM1
+      IF(VC(i,j1).GT.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,5 +1346,5 @@
      &                        A6(i,j1)*(1.+R23*VC(i,j1)))
       endif
-      END DO
+140   continue
       return
       end
@@ -1429,5 +1378,5 @@
 c        write(*,*) 'toto 1' 
 C --------------------------------
-      IF(IAD==2) then
+      IF(IAD.eq.2) then
       do j=j1-1,j2+1
       do i=1,IMR
@@ -1446,5 +1395,5 @@
 c      write(*,*) 'toto 2'
 C
-      ELSEIF(IAD==1) then
+      ELSEIF(IAD.eq.1) then
 	do j=j1-1,j2+1
       do i=1,imr
@@ -1455,5 +1404,5 @@
       ENDIF
 C
-	if(j1/=2) then
+	if(j1.ne.2) then
 	sum1 = 0.
 	sum2 = 0.
@@ -1499,6 +1448,6 @@
 C
 	JMR = JNP-1
-      do j=j1,j2
-      if(J>JS  .and. J<JN) GO TO 1309
+      do 1309 j=j1,j2
+      if(J.GT.JS  .and. J.LT.JN) GO TO 1309
 C
       do i=1,IMR
@@ -1511,5 +1460,5 @@
       enddo
 C
-      IF(IAD==2) THEN
+      IF(IAD.eq.2) THEN
       DO i=1,IMR
       IP = NINT(UA(i,j))
@@ -1520,10 +1469,10 @@
       adx(i,j) = qtmp(ip) + ru*(a1*ru + b1)
       enddo
-      ELSEIF(IAD==1) then
+      ELSEIF(IAD.eq.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).GE.0.) then
       adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))
       else
@@ -1537,5 +1486,4 @@
       enddo
 1309  continue
-      END DO
 C
 C Eulerian upwind
@@ -1550,5 +1498,5 @@
       qtmp(IMR+1) = p(1,J)
 C
-      IF(IAD==2) THEN
+      IF(IAD.eq.2) THEN
       qtmp(-1)     = p(IMR-1,J)
       qtmp(IMR+2) = p(2,J)
@@ -1561,5 +1509,5 @@
       adx(i,j) = qtmp(ip)- p(i,j) + ru*(a1*ru + b1)
       enddo
-      ELSEIF(IAD==1) then
+      ELSEIF(IAD.eq.1) then
 C 1st order
       DO i=1,IMR
@@ -1570,5 +1518,5 @@
       enddo
 C
-	if(j1/=2) then
+	if(j1.ne.2) then
       do i=1,IMR
       adx(i,  2) = 0.
@@ -1606,8 +1554,8 @@
       REAL da1,da2,a6da,fmin
 C
-      if(LMT==0) then
+      if(LMT.eq.0) then
 C Full constraint
-      do i=1,IM
-      if(DC(i)==0.) then
+      do 100 i=1,IM
+      if(DC(i).eq.0.) then
             AR(i) = p(i)
             AL(i) = p(i)
@@ -1617,22 +1565,22 @@
       da2  = da1**2
       A6DA = A6(i)*da1
-      if(A6DA < -da2) then
+      if(A6DA .lt. -da2) then
             A6(i) = 3.*(AL(i)-p(i))
             AR(i) = AL(i) - A6(i)
-      elseif(A6DA > da2) then
+      elseif(A6DA .gt. da2) then
             A6(i) = 3.*(AR(i)-p(i))
             AL(i) = AR(i) - A6(i)
       endif
       endif
-      END DO
-      elseif(LMT==1) then
+100   continue
+      elseif(LMT.eq.1) then
 C 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
+      do 150 i=1,IM
+      if(abs(AR(i)-AL(i)) .GE. -A6(i)) go to 150
+      if(p(i).lt.AR(i) .and. p(i).lt.AL(i)) then
             AR(i) = p(i)
             AL(i) = p(i)
             A6(i) = 0.
-      elseif(AR(i) > AL(i)) then
+      elseif(AR(i) .gt. AL(i)) then
             A6(i) = 3.*(AL(i)-p(i))
             AR(i) = AL(i) - A6(i)
@@ -1642,15 +1590,14 @@
       endif
 150   continue
-      END DO
-      elseif(LMT==2) then
-      do i=1,IM
-      if(abs(AR(i)-AL(i)) >= -A6(i)) go to 250
+      elseif(LMT.eq.2) then
+      do 250 i=1,IM
+      if(abs(AR(i)-AL(i)) .GE. -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.ge.0.) go to 250
+      if(p(i).lt.AR(i) .and. p(i).lt.AL(i)) then
             AR(i) = p(i)
             AL(i) = p(i)
             A6(i) = 0.
-      elseif(AR(i) > AL(i)) then
+      elseif(AR(i) .gt. AL(i)) then
             A6(i) = 3.*(AL(i)-p(i))
             AR(i) = AL(i) - A6(i)
@@ -1660,5 +1607,4 @@
       endif
 250   continue
-      END DO
       endif
       return
@@ -1671,17 +1617,13 @@
       integer i,j
 C
-      do j=j1,j2
-      do i=2,IMR
-      CRX(i,J) = dtdx5(j)*(U(i,j)+U(i-1,j))
-      END DO
-      END DO
-C
-      do j=j1,j2
-      CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j))
-      END DO
-C
-      do i=1,IMR*JMR
-      CRY(i,2) = DTDY5*(V(i,2)+V(i,1))
-      END DO
+      do 35 j=j1,j2
+      do 35 i=2,IMR
+35    CRX(i,J) = dtdx5(j)*(U(i,j)+U(i-1,j))
+C
+      do 45 j=j1,j2
+45    CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j))
+C
+      do 55 i=1,IMR*JMR
+55    CRY(i,2) = DTDY5*(V(i,2)+V(i,1))
       return
       end
@@ -1694,11 +1636,10 @@
       real ph5
       JMR = JNP-1
-      do j=2,JNP
+      do 55 j=2,JNP
         ph5  =  -0.5*PI + (REAL(J-1)-0.5)*DP
-        cose(j) = cos(ph5)
-      END DO
+55      cose(j) = cos(ph5)
 C
       JEQ = (JNP+1) / 2
-      if(JMR == 2*(JMR/2) ) then
+      if(JMR .eq. 2*(JMR/2) ) then
       do j=JNP, JEQ+1, -1
        cose(j) =  cose(JNP+2-j)
@@ -1712,7 +1653,6 @@
       endif
 C
-      do j=2,JMR
-      cosp(j) = 0.5*(cose(j)+cose(j+1))
-      END DO
+      do 66 j=2,JMR
+66    cosp(j) = 0.5*(cose(j)+cose(j+1))
       cosp(1) = 0.
       cosp(JNP) = 0.
@@ -1728,18 +1668,17 @@
 C
       phi = -0.5*PI
-      do j=2,JNP-1
+      do 55 j=2,JNP-1
       phi  =  phi + DP
-      cosp(j) = cos(phi)
-      END DO
+55    cosp(j) = cos(phi)
         cosp(  1) = 0.
         cosp(JNP) = 0.
 C
-      do j=2,JNP
+      do 66 j=2,JNP
         cose(j) = 0.5*(cosp(j)+cosp(j-1))
-      END DO
-C
-      do j=2,JNP-1
+66    CONTINUE
+C
+      do 77 j=2,JNP-1
        cosp(j) = 0.5*(cose(j)+cose(j+1))
-      END DO
+77    CONTINUE
       return
       end
@@ -1763,16 +1702,16 @@
 	icr = 1
       call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
-      if(ipy==0) goto 50
+      if(ipy.eq.0) goto 50
       call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
-      if(ipx==0) goto 50
+      if(ipx.eq.0) goto 50
 C
       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.eq.0) goto 50
 C
 C Vertical filling...
       do i=1,len
-      IF( Q(i,j1,1)<0.) THEN
+      IF( Q(i,j1,1).LT.0.) THEN
       ip = ip + 1
           Q(i,j1,2) = Q(i,j1,2) + Q(i,j1,1)
@@ -1782,18 +1721,18 @@
 C
 50    continue
-      DO L = 2,NLAYM1
+      DO 225 L = 2,NLAYM1
       icr = 1
 C
       call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
-      if(ipy==0) goto 225
+      if(ipy.eq.0) goto 225
       call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
-      if(ipx==0) go to 225
+      if(ipx.eq.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.eq.0) goto 225
 C
       do i=1,len
-      IF( Q(I,j1,L)<0.) THEN
+      IF( Q(I,j1,L).LT.0.) THEN
 C
       ip = ip + 1
@@ -1810,5 +1749,4 @@
       ENDDO
 225   CONTINUE
-      END DO
 C
 C BOTTOM LAYER
@@ -1817,13 +1755,13 @@
 C
       call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
-      if(ipy==0) goto 911
+      if(ipy.eq.0) goto 911
       call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
-      if(ipx==0) goto 911
+      if(ipx.eq.0) goto 911
 C
       call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
-      if(icr==0) goto 911
+      if(icr.eq.0) goto 911
 C
       DO  I=1,len
-      IF( Q(I,j1,L)<0.) THEN
+      IF( Q(I,j1,L).LT.0.) THEN
       ip = ip + 1
 c
@@ -1842,10 +1780,10 @@
 911   continue
 C
-      if(ip>IMR) then
+      if(ip.gt.IMR) then
       write(6,*) 'IC=',IC,' STEP=',NSTEP,
      &           ' Vertical filling pts=',ip
       endif
 C
-      if(sum>1.e-25) then
+      if(sum.gt.1.e-25) then
       write(6,*) IC,NSTEP,' Mass source from the ground=',sum
       endif
@@ -1860,7 +1798,7 @@
       real :: dq,dn,d0,d1,ds,d2
       icr = 0
-      do j=j1+1,j2-1
-      DO i=1,IMR-1
-      IF(q(i,j)<0.) THEN
+      do 65 j=j1+1,j2-1
+      DO 50 i=1,IMR-1
+      IF(q(i,j).LT.0.) THEN
       icr =  1
       dq  = - q(i,j)*cosp(j)
@@ -1878,8 +1816,8 @@
       q(i,j) = (d2 - dq)*acosp(j) + tiny
       endif
-      END DO
-      if(icr==0 .and. q(IMR,j)>=0.) goto 65
-      DO i=2,IMR
-      IF(q(i,j)<0.) THEN
+50    continue
+      if(icr.eq.0 .and. q(IMR,j).ge.0.) goto 65
+      DO 55 i=2,IMR
+      IF(q(i,j).LT.0.) THEN
       icr =  1
       dq  = - q(i,j)*cosp(j)
@@ -1897,9 +1835,9 @@
       q(i,j) = (d2 - dq)*acosp(j) + tiny
       endif
-      END DO
+55    continue
 C *****************************************
 C i=1
       i=1
-      IF(q(i,j)<0.) THEN
+      IF(q(i,j).LT.0.) THEN
       icr =  1
       dq  = - q(i,j)*cosp(j)
@@ -1920,5 +1858,5 @@
 C i=IMR
       i=IMR
-      IF(q(i,j)<0.) THEN
+      IF(q(i,j).LT.0.) THEN
       icr =  1
       dq  = - q(i,j)*cosp(j)
@@ -1938,8 +1876,7 @@
 C *****************************************
 65    continue
-      END DO
 C
       do i=1,IMR
-      if(q(i,j1)<0. .or. q(i,j2)<0.) then
+      if(q(i,j1).lt.0. .or. q(i,j2).lt.0.) then
       icr = 1
       goto 80
@@ -1949,5 +1886,5 @@
 80    continue
 C
-      if(q(1,1)<0. .or. q(1,jnp)<0.) then
+      if(q(1,1).lt.0. .or. q(1,jnp).lt.0.) then
       icr = 1
       endif
@@ -1973,7 +1910,7 @@
 C
       ipy = 0
-      do j=j1+1,j2-1
-      DO i=1,IMR
-      IF(q(i,j)<0.) THEN
+      do 55 j=j1+1,j2-1
+      DO 55 i=1,IMR
+      IF(q(i,j).LT.0.) THEN
       ipy =  1
       dq  = - q(i,j)*cosp(j)
@@ -1991,9 +1928,8 @@
       q(i,j) = (d2 - dq)*acosp(j) + tiny
       endif
-      END DO
-      END DO
+55    continue
 C
       do i=1,imr
-      IF(q(i,j1)<0.) THEN
+      IF(q(i,j1).LT.0.) THEN
       ipy =  1
       dq  = - q(i,j1)*cosp(j1)
@@ -2009,5 +1945,5 @@
       j = j2
       do i=1,imr
-      IF(q(i,j)<0.) THEN
+      IF(q(i,j).LT.0.) THEN
       ipy =  1
       dq  = - q(i,j)*cosp(j)
@@ -2022,19 +1958,19 @@
 C
 C Check Poles.
-      if(q(1,1)<0.) then
+      if(q(1,1).lt.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
-C
-      if(q(1,JNP)<0.) then
+      if(q(i,j1).lt.0.) ipy = 1
+      enddo
+      endif
+C
+      if(q(1,JNP).lt.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
+      if(q(i,j2).lt.0.) ipy = 1
       enddo
       endif
@@ -2052,13 +1988,11 @@
       ipx = 0
 C Copy & swap direction for vectorization.
-      do i=1,imr
-      do j=j1,j2
-      qtmp(j,i) = q(i,j)
-      END DO
-      END DO
-C
-      do i=2,imr-1
-      do j=j1,j2
-      if(qtmp(j,i)<0.) then
+      do 25 i=1,imr
+      do 25 j=j1,j2
+25    qtmp(j,i) = q(i,j)
+C
+      do 55 i=2,imr-1
+      do 55 j=j1,j2
+      if(qtmp(j,i).lt.0.) then
       ipx =  1
 c west
@@ -2073,10 +2007,9 @@
       qtmp(j,i) = qtmp(j,i) + d2 + tiny
       endif
-      END DO
-      END DO
+55    continue
 c
       i=1
-      do j=j1,j2
-      if(qtmp(j,i)<0.) then
+      do 65 j=j1,j2
+      if(qtmp(j,i).lt.0.) then
       ipx =  1
 c west
@@ -2092,8 +2025,8 @@
       qtmp(j,i) = qtmp(j,i) + d2 + tiny
       endif
-      END DO
+65    continue
       i=IMR
-      do j=j1,j2
-      if(qtmp(j,i)<0.) then
+      do 75 j=j1,j2
+      if(qtmp(j,i).lt.0.) then
       ipx =  1
 c west
@@ -2109,16 +2042,14 @@
       qtmp(j,i) = qtmp(j,i) + d2 + tiny
       endif
-      END DO
-C
-      if(ipx/=0) then
-      do j=j1,j2
-      do i=1,imr
-      q(i,j) = qtmp(j,i)
-      END DO
-      END DO
+75    continue
+C
+      if(ipx.ne.0) then
+      do 85 j=j1,j2
+      do 85 i=1,imr
+85    q(i,j) = qtmp(j,i)
       else
 C
 C Poles.
-      if(q(1,1)<0 .or. q(1,JNP)<0.) ipx = 1
+      if(q(1,1).lt.0. or. q(1,JNP).lt.0.) ipx = 1
       endif
       return
@@ -2134,16 +2065,14 @@
       integer IC,k,i
 C
-      do IC = 1, nc
-C
-      do k=1,km
-      do i=1,im
+      do 4000 IC = 1, nc
+C
+      do 1000 k=1,km
+      do 1000 i=1,im
       qtmp(i,k) = q(i,km+1-k,IC)
-      END DO
-      END DO
-C
-      do i=1,im*km
-      q(i,1,IC) = qtmp(i,1)
-      END DO
-      END DO
+1000  continue
+C
+      do 2000 i=1,im*km
+2000  q(i,1,IC) = qtmp(i,1)
+4000  continue
       return
       end
