diff --git a/src/MNH/mode_RBK90_linearalgebra.f90 b/src/MNH/mode_RBK90_linearalgebra.f90
index 8fd23aa139e89b8c906915ccda98c36e80e102e0..de787d50d3e04f7b43d59230b4f92f7d525bf730 100644
--- a/src/MNH/mode_RBK90_linearalgebra.f90
+++ b/src/MNH/mode_RBK90_linearalgebra.f90
@@ -1,10 +1,11 @@
-!MNH_LIC Copyright 1994-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
 ! Modifications:
 !  P. Wautelet 22/02/2019: DOUBLE COMPLEX -> COMPLEX(kind(0.0d0)) to respect Fortran standard
+!  P. Wautelet 17/12/2021: convert arithmetic if to classic if (deleted from Fortran 2018 standard)
 !-----------------------------------------------------------------
 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ! 
@@ -1261,46 +1262,49 @@ END SUBROUTINE KppSolveTR
                                  
       WDOT = 0.0D0 
       IF (N .LE. 0) RETURN 
-      IF (incX .EQ. incY) IF (incX-1) 5,20,60 
+      IF (incX .EQ. incY) THEN
 !                                                                       
 !     Code for unequal or nonpositive increments.                       
 !                                                                       
-    5 IX = 1 
-      IY = 1 
-      IF (incX .LT. 0) IX = (-N+1)*incX + 1 
-      IF (incY .LT. 0) IY = (-N+1)*incY + 1 
-      DO i = 1,N 
-        WDOT = WDOT + DX(IX)*DY(IY) 
-        IX = IX + incX 
-        IY = IY + incY 
-      END DO 
-      RETURN 
+        IF ( incX < 1 ) THEN
+          IX = 1
+          IY = 1
+          IF (incX .LT. 0) IX = (-N+1)*incX + 1
+          IF (incY .LT. 0) IY = (-N+1)*incY + 1
+          DO i = 1,N
+            WDOT = WDOT + DX(IX)*DY(IY)
+            IX = IX + incX
+            IY = IY + incY
+          END DO
+        ELSE IF ( incX == 1 ) THEN
 !                                                                       
 !     Code for both increments equal to 1.                              
 !                                                                       
 !     Clean-up loop so remaining vector length is a multiple of 5.      
 !                                                                       
-   20 M = MOD(N,5) 
-      IF (M .EQ. 0) GO TO 40 
-      DO i = 1,M 
-         WDOT = WDOT + DX(i)*DY(i) 
-      END DO 
-      IF (N .LT. 5) RETURN 
-   40 MP1 = M + 1 
-      DO i = MP1,N,5 
-          WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) +  &
-                   DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4)                   
-      END DO 
-      RETURN 
+          M = MOD(N,5)
+          IF (M .EQ. 0) GO TO 40
+          DO i = 1,M
+            WDOT = WDOT + DX(i)*DY(i)
+          END DO
+          IF (N .LT. 5) RETURN
+   40     MP1 = M + 1
+          DO i = MP1,N,5
+            WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) +  &
+                     DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4)
+          END DO
+        ELSE
 !                                                                       
 !     Code for equal, positive, non-unit increments.                    
 !                                                                       
-   60 NS = N*incX 
-      DO i = 1,NS,incX 
-        WDOT = WDOT + DX(i)*DY(i) 
-      END DO 
+          NS = N*incX
+          DO i = 1,NS,incX
+            WDOT = WDOT + DX(i)*DY(i)
+          END DO
+        END IF
+      END IF
 
-      END FUNCTION WDOT                                          
+      END FUNCTION WDOT
 
 
 !--------------------------------------------------------------