diff --git a/src/testprogs/ice_adjust/getdata_ice_adjust_mod.F90 b/src/testprogs/ice_adjust/getdata_ice_adjust_mod.F90
index 3637d1b4b3626b1c970540f8ac197685bae4105e..1372fcd929085a884fa12d4215b9809f555dce81 100644
--- a/src/testprogs/ice_adjust/getdata_ice_adjust_mod.F90
+++ b/src/testprogs/ice_adjust/getdata_ice_adjust_mod.F90
@@ -1,27 +1,7 @@
 MODULE GETDATA_ICE_ADJUST_MOD
 
 USE OMP_LIB
-
-INTERFACE REPLICATE
-  MODULE PROCEDURE REPLICATE3
-  MODULE PROCEDURE REPLICATE4
-END INTERFACE
-
-INTERFACE NPROMIZE
-  MODULE PROCEDURE NPROMIZE4
-  MODULE PROCEDURE NPROMIZE5
-END INTERFACE
-
-INTERFACE INTERPOLATE
-  MODULE PROCEDURE INTERPOLATE4
-  MODULE PROCEDURE INTERPOLATE5
-END INTERFACE
-
-INTERFACE SET
-  MODULE PROCEDURE SET3
-  MODULE PROCEDURE SET4
-  MODULE PROCEDURE SET5
-END INTERFACE
+USE ARRAYS_MANIP, ONLY: SETUP, XNAN, REPLICATE, NPROMIZE, INTERPOLATE, SET
 
 CONTAINS
 
@@ -30,8 +10,6 @@ SUBROUTINE GETDATA_ICE_ADJUST (NPROMA, NGPBLKS, NFLEVG, PRHODJ_B, PEXNREF_B, PRH
 & PHLI_HRI_B, PHLI_HCF_B, ZRS_B, ZZZ_B, PRS_OUT_B, PSRCS_OUT_B, PCLDFR_OUT_B, PHLC_HRC_OUT_B, PHLC_HCF_OUT_B,         &
 & PHLI_HRI_OUT_B, PHLI_HCF_OUT_B, LDVERBOSE)
 
-USE IEEE_ARITHMETIC, ONLY : IEEE_SIGNALING_NAN, IEEE_VALUE
-
 IMPLICIT NONE
 
 INTEGER, PARAMETER :: IFILE = 77
@@ -106,7 +84,8 @@ INTEGER :: NGPTOT, NPROMA, NGPBLKS, NFLEVG
 INTEGER :: IOFF, IBL
 LOGICAL :: LLEXIST
 CHARACTER(LEN=32) :: CLFILE
-REAL :: ZNAN
+
+CALL SETUP()
 
 KRR=6
 NGPTOT = NPROMA * NGPBLKS
@@ -154,9 +133,6 @@ ALLOCATE (PHLC_HCF_OUT_B  (NPROMA,1,NFLEVG,NGPBLKS))
 ALLOCATE (PHLI_HRI_OUT_B  (NPROMA,1,NFLEVG,NGPBLKS))
 ALLOCATE (PHLI_HCF_OUT_B  (NPROMA,1,NFLEVG,NGPBLKS))
 
-ZNAN = IEEE_VALUE (ZNAN, IEEE_SIGNALING_NAN)
-
-
 CALL SET (ZSIGQSAT_B    )
 CALL SET (ZICE_CLD_WGT_B)
 CALL SET (PSRCS_B       )
@@ -191,12 +167,12 @@ CALL SET (PHLI_HCF_OUT_B)
 
 ZSIGQSAT_B     = 2.0000000000000000E-002
 ZICE_CLD_WGT_B = 1.5
-PSRCS_B        = ZNAN
-PCLDFR_B       = ZNAN
-PHLI_HCF_B     = ZNAN
-PHLI_HRI_B     = ZNAN
-PHLC_HCF_B     = ZNAN
-PHLC_HRC_B     = ZNAN
+PSRCS_B        = XNAN
+PCLDFR_B       = XNAN
+PHLI_HCF_B     = XNAN
+PHLI_HRI_B     = XNAN
+PHLC_HCF_B     = XNAN
+PHLC_HRC_B     = XNAN
 
 
 IOFF = 0
@@ -339,252 +315,4 @@ CALL NPROMIZE (NPROMA, PHLI_HCF_OUT,  PHLI_HCF_OUT_B  )
 
 END SUBROUTINE 
 
-SUBROUTINE REPLICATE4 (KOFF, P)
-
-INTEGER :: KOFF
-REAL    :: P (:,:,:,:)
-
-INTEGER :: I, J
-
-DO I = KOFF+1, SIZE (P, 1)
-  J = 1 + MODULO (I - 1, KOFF)
-  P (I, :, :, :) = P (J, :, :, :)
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE REPLICATE3 (KOFF, P)
-
-INTEGER :: KOFF
-REAL    :: P (:,:,:)
-
-INTEGER :: I, J
-
-DO I = KOFF+1, SIZE (P, 1)
-  J = 1 + MODULO (I - 1, KOFF)
-  P (I, :, :) = P (J, :, :)
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE NPROMIZE4 (KPROMA, PI, PO)
-
-INTEGER :: KPROMA
-REAL, INTENT (IN)  :: PI (:,:,:,:) 
-REAL, INTENT (OUT) :: PO (:,:,:,:)
-
-INTEGER :: I, J, IGPBLK, IGPTOT, IGP, JLON, JIDIA, JFDIA
-
-IF (SIZE (PI, 4) /= 1) STOP 1
-
-IGPTOT = SIZE (PI, 1)
-IGPBLK = 1 + (IGPTOT-1) / KPROMA
-
-DO IGP = 1, IGPTOT, KPROMA
-  IBL = 1 + (IGP - 1) / KPROMA
-  JIDIA = 1
-  JFDIA = MIN (KPROMA, IGPTOT - (IBL - 1) * KPROMA)
-
-  DO JLON = JIDIA, JFDIA
-    PO (JLON, :, :, IBL) = PI (IGP + (JLON - 1), :, :, 1)
-  ENDDO
-
-  DO JLON = JFDIA+1, KPROMA
-    PO (JLON, :, :, IBL) = PO (JFDIA, :, :, IBL)
-  ENDDO
-
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE NPROMIZE5 (KPROMA, PI, PO)
-
-INTEGER :: KPROMA
-REAL, INTENT (IN)  :: PI (:,:,:,:,:) 
-REAL, INTENT (OUT) :: PO (:,:,:,:,:)
-
-INTEGER :: I, J, IGPBLK, IGPTOT, IGP, JLON, JIDIA, JFDIA
-
-IF (SIZE (PI, 5) /= 1) STOP 1
-
-IGPTOT = SIZE (PI, 1)
-IGPBLK = 1 + (IGPTOT-1) / KPROMA
-
-DO IGP = 1, IGPTOT, KPROMA
-  IBL = 1 + (IGP - 1) / KPROMA
-  JIDIA = 1
-  JFDIA = MIN (KPROMA, IGPTOT - (IBL - 1) * KPROMA)
-
-  DO JLON = JIDIA, JFDIA
-    PO (JLON, :, :, :, IBL) = PI (IGP + (JLON - 1), :, :, :, 1)
-  ENDDO
-
-  DO JLON = JFDIA+1, KPROMA
-    PO (JLON, :, :, :, IBL) = PI (JFDIA, :, :, :, IBL)
-  ENDDO
-
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE INTERPOLATE4 (KFLEVG, KOFF, P)
-
-INTEGER :: KFLEVG, KOFF
-REAL, ALLOCATABLE :: P (:,:,:,:)
-REAL :: Z (LBOUND (P, 1):UBOUND (P, 1), &
-         & LBOUND (P, 2):UBOUND (P, 2), &
-         & LBOUND (P, 3):UBOUND (P, 3), &
-         & LBOUND (P, 4):UBOUND (P, 4))
-INTEGER :: ILEV1A, ILEV1B, ILEV2, NLEV1, NLEV2
-REAL :: ZWA, ZWB, ZLEV1, ZLEV2
-
-Z = P
-
-NLEV1 = SIZE (P, 3)
-NLEV2 = KFLEVG
-
-DEALLOCATE (P)
-
-ALLOCATE (P (LBOUND (Z, 1):UBOUND (Z, 1), &
-           & LBOUND (Z, 2):UBOUND (Z, 2), &
-           & KFLEVG, &
-           & LBOUND (Z, 4):UBOUND (Z, 4)))
-
-DO ILEV2 = 1, NLEV2
-  ZLEV2 = REAL (ILEV2 - 1) / REAL (NLEV2 -1)
-  ZLEV1 = 1. + ZLEV2 * REAL (NLEV1 - 1)
-  ILEV1B = MIN (CEILING (ZLEV1), NLEV1)
-  ILEV1A = MAX (FLOOR   (ZLEV1),     1)
-
-  IF (ILEV1A == ILEV1B) THEN
-    ZWA = 1.
-    ZWB = 0.
-  ELSE
-    ZWA = REAL (ILEV1B) - ZLEV1
-    ZWB = ZLEV1 - REAL (ILEV1A)
-  ENDIF
-
-! WRITE (*, '(" ZLEV2 = ",E12.5," ZLEV1 = ",E12.5," ILEV2 = ",I4," ILEV1A = ",I4," ZWA = ",E12.5," ILEV1B = ",I4," ZWB = ",E12.5)') &
-!   & ZLEV2, ZLEV1, ILEV2, ILEV1A, ZWA, ILEV1B, ZWB
-
-  P (1:KOFF, :, ILEV2, :) = ZWA * Z (1:KOFF, :, ILEV1A, :) + ZWB * Z (1:KOFF, :, ILEV1B, :) 
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE INTERPOLATE5 (KFLEVG, KOFF, P)
-
-INTEGER :: KFLEVG, KOFF
-REAL, ALLOCATABLE :: P (:,:,:,:,:)
-REAL :: Z (LBOUND (P, 1):UBOUND (P, 1), &
-         & LBOUND (P, 2):UBOUND (P, 2), &
-         & LBOUND (P, 3):UBOUND (P, 3), &
-         & LBOUND (P, 4):UBOUND (P, 4), &
-         & LBOUND (P, 5):UBOUND (P, 5))
-INTEGER :: ILEV1A, ILEV1B, ILEV2, NLEV1, NLEV2
-REAL :: ZWA, ZWB, ZLEV1, ZLEV2
-
-Z = P
-
-NLEV1 = SIZE (P, 3)
-NLEV2 = KFLEVG
-
-DEALLOCATE (P)
-
-ALLOCATE (P (LBOUND (Z, 1):UBOUND (Z, 1), &
-           & LBOUND (Z, 2):UBOUND (Z, 2), &
-           & KFLEVG, &
-           & LBOUND (Z, 4):UBOUND (Z, 4), &
-           & LBOUND (Z, 5):UBOUND (Z, 5)))
-
-DO ILEV2 = 1, NLEV2
-  ZLEV2 = REAL (ILEV2 - 1) / REAL (NLEV2 -1)
-  ZLEV1 = 1. + ZLEV2 * REAL (NLEV1 - 1)
-  ILEV1B = MIN (CEILING (ZLEV1), NLEV1)
-  ILEV1A = MAX (FLOOR   (ZLEV1),     1)
-
-  IF (ILEV1A == ILEV1B) THEN
-    ZWA = 1.
-    ZWB = 0.
-  ELSE
-    ZWA = REAL (ILEV1B) - ZLEV1
-    ZWB = ZLEV1 - REAL (ILEV1A)
-  ENDIF
-
-! WRITE (*, '(" ZLEV2 = ",E12.5," ZLEV1 = ",E12.5," ILEV2 = ",I4," ILEV1A = ",I4," ZWA = ",E12.5," ILEV1B = ",I4," ZWB = ",E12.5)') &
-!   & ZLEV2, ZLEV1, ILEV2, ILEV1A, ZWA, ILEV1B, ZWB
-
-  P (1:KOFF, :, ILEV2, :, :) = ZWA * Z (1:KOFF, :, ILEV1A, :, :) + ZWB * Z (1:KOFF, :, ILEV1B, :, :) 
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE SET3 (P)
-
-REAL :: P (:,:,:)
-INTEGER :: IBL, IGPBLKS
-INTEGER :: NTID, ITID, JBLK1, JBLK2
-
-
-IGPBLKS = SIZE (P, 3)
-
-!$OMP PARALLEL PRIVATE (ITID, JBLK1, JBLK2, NTID)
-NTID = OMP_GET_MAX_THREADS ()
-ITID = OMP_GET_THREAD_NUM ()
-JBLK1 = 1 +  (IGPBLKS * (ITID+0)) / NTID
-JBLK2 =      (IGPBLKS * (ITID+1)) / NTID
-
-DO IBL = JBLK1, JBLK2
-  P (:,:,IBL) = ZNAN
-ENDDO
-
-!$OMP END PARALLEL
-
-END SUBROUTINE
-
-SUBROUTINE SET4 (P)
-
-REAL :: P (:,:,:,:)
-INTEGER :: IBL, IGPBLKS
-INTEGER :: NTID, ITID, JBLK1, JBLK2
-
-IGPBLKS = SIZE (P, 4)
-
-!$OMP PARALLEL PRIVATE (ITID, JBLK1, JBLK2, NTID)
-NTID = OMP_GET_MAX_THREADS ()
-ITID = OMP_GET_THREAD_NUM ()
-JBLK1 = 1 +  (IGPBLKS * (ITID+0)) / NTID
-JBLK2 =      (IGPBLKS * (ITID+1)) / NTID
-
-DO IBL = JBLK1, JBLK2
-  P (:,:,:,IBL) = ZNAN
-ENDDO
-
-!$OMP END PARALLEL
-
-END SUBROUTINE
-
-SUBROUTINE SET5 (P)
-
-REAL :: P (:,:,:,:,:)
-INTEGER :: IBL, IGPBLKS
-INTEGER :: NTID, ITID, JBLK1, JBLK2
-
-IGPBLKS = SIZE (P, 5)
-
-!$OMP PARALLEL PRIVATE (ITID, JBLK1, JBLK2, NTID)
-NTID = OMP_GET_MAX_THREADS ()
-ITID = OMP_GET_THREAD_NUM ()
-JBLK1 = 1 +  (IGPBLKS * (ITID+0)) / NTID
-JBLK2 =      (IGPBLKS * (ITID+1)) / NTID
-
-DO IBL = JBLK1, JBLK2
-  P (:,:,:,:,IBL) = ZNAN
-ENDDO
-
-!$OMP END PARALLEL
-
-END SUBROUTINE
-
-
 END  MODULE
diff --git a/src/testprogs/ice_adjust/main_ice_adjust.F90 b/src/testprogs/ice_adjust/main_ice_adjust.F90
index c488c4b41731211a88843f457b85857b600f03c7..8bf8401d9cdd7fe7340c2a69a91615407829466f 100644
--- a/src/testprogs/ice_adjust/main_ice_adjust.F90
+++ b/src/testprogs/ice_adjust/main_ice_adjust.F90
@@ -2,6 +2,7 @@ PROGRAM MAIN_ICE_ADJUST
 
 USE XRD_GETOPTIONS
 USE GETDATA_ICE_ADJUST_MOD
+USE COMPUTE_DIFF
 USE MODI_ICE_ADJUST
 USE MODD_DIMPHYEX,   ONLY: DIMPHYEX_t
 USE MODD_CST,        ONLY: CST
@@ -331,12 +332,12 @@ PRINT *, " ZTC = ", ZTC, ZTC / REAL (NPROMA*NGPBLKS*NTIME)
 IF (LLCHECK .OR. LLSTAT .OR. LLCHECKDIFF) THEN
   DO IBL = IBLOCK1, IBLOCK2
     PRINT *, " IBL = ", IBL
-    CALL DIFF ("PSRCS",    PSRCS_OUT    (:,:,:,IBL), PSRCS    (:,:,:,IBL))
-    CALL DIFF ("PCLDFR",   PCLDFR_OUT   (:,:,:,IBL), PCLDFR   (:,:,:,IBL))
-    CALL DIFF ("PHLC_HRC", PHLC_HRC_OUT (:,:,:,IBL), PHLC_HRC (:,:,:,IBL))
-    CALL DIFF ("PHLC_HCF", PHLC_HCF_OUT (:,:,:,IBL), PHLC_HCF (:,:,:,IBL))
-    CALL DIFF ("PHLI_HRI", PHLI_HRI_OUT (:,:,:,IBL), PHLI_HRI (:,:,:,IBL))
-    CALL DIFF ("PHLI_HCF", PHLI_HCF_OUT (:,:,:,IBL), PHLI_HCF (:,:,:,IBL))
+    CALL DIFF ("PSRCS",    PSRCS_OUT    (:,:,:,IBL), PSRCS    (:,:,:,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
+    CALL DIFF ("PCLDFR",   PCLDFR_OUT   (:,:,:,IBL), PCLDFR   (:,:,:,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
+    CALL DIFF ("PHLC_HRC", PHLC_HRC_OUT (:,:,:,IBL), PHLC_HRC (:,:,:,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
+    CALL DIFF ("PHLC_HCF", PHLC_HCF_OUT (:,:,:,IBL), PHLC_HCF (:,:,:,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
+    CALL DIFF ("PHLI_HRI", PHLI_HRI_OUT (:,:,:,IBL), PHLI_HRI (:,:,:,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
+    CALL DIFF ("PHLI_HCF", PHLI_HCF_OUT (:,:,:,IBL), PHLI_HCF (:,:,:,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
   ENDDO
 ENDIF
 
@@ -350,53 +351,4 @@ ENDIF
 
 STOP
 
-CONTAINS
-
-SUBROUTINE DIFF (CDNAME, PREF, POUT)
-
-CHARACTER (LEN=*) :: CDNAME
-REAL :: PREF (:,:,:)
-REAL :: POUT (:,:,:)
-
-INTEGER :: JLON, JLEV
-
-PRINT *, CDNAME
-IF (LLSTAT) THEN
-  PRINT *, MINVAL (PREF), MAXVAL (PREF), SUM (PREF) / SIZE (PREF)
-  PRINT *, MINVAL (POUT), MAXVAL (POUT), SUM (POUT) / SIZE (POUT)
-ENDIF
-
-IF (LLCHECK) THEN
-  IF (SUM (ABS (POUT) + ABS (PREF)) > 0) THEN
-  WRITE (*, '(A4)', ADVANCE='NO') ""
-  DO JLON = 1, NPROMA
-    WRITE (*, '("|",I12,A12)', ADVANCE='NO') JLON, ""
-  ENDDO
-  WRITE (*, '("|")')
-  DO JLEV = 1, KLEV
-    WRITE (*, '(I4)', ADVANCE='NO') JLEV
-    DO JLON = 1, NPROMA
-      IF (ABS (PREF (JLON, 1, JLEV)) + ABS (POUT (JLON, 1, JLEV)) == 0.) THEN
-      WRITE (*, '("|",2A12)', ADVANCE='NO') "", ""
-      ELSE
-      WRITE (*, '("|",2E12.5)', ADVANCE='NO') PREF (JLON, 1, JLEV), POUT (JLON, 1, JLEV)
-      ENDIF
-    ENDDO
-    WRITE (*, '("|")')
-  ENDDO
-  ENDIF
-ENDIF
-
-IF (LLCHECKDIFF) THEN
-  IF (SUM(ABS(POUT-PREF)) > 0.) THEN
-    PRINT*, "THERE ARE DIFF"
-    LLDIFF = .TRUE.
-  ELSE
-    PRINT*, "THERE IS NO DIFF"
-  ENDIF
-ENDIF
-
-END SUBROUTINE
-
-
 END
diff --git a/src/testprogs/rain_ice/getdata_rain_ice_mod.F90 b/src/testprogs/rain_ice/getdata_rain_ice_mod.F90
index 57c28b427d7a1b776b1ffcf68894bbf84e8b9e23..8940901fd1758dbe0374a67d5b9fdbe2511ef515 100644
--- a/src/testprogs/rain_ice/getdata_rain_ice_mod.F90
+++ b/src/testprogs/rain_ice/getdata_rain_ice_mod.F90
@@ -1,32 +1,6 @@
 MODULE GETDATA_RAIN_ICE_MOD
 
-USE OMP_LIB
-
-INTERFACE REPLICATE
-  MODULE PROCEDURE REPLICATE2
-  MODULE PROCEDURE REPLICATE3
-  MODULE PROCEDURE REPLICATE4
-  MODULE PROCEDURE REPLICATEL
-END INTERFACE
-
-INTERFACE NPROMIZE
-  MODULE PROCEDURE NPROMIZE3
-  MODULE PROCEDURE NPROMIZE4
-  MODULE PROCEDURE NPROMIZE5
-  MODULE PROCEDURE NPROMIZEL
-END INTERFACE
-
-INTERFACE INTERPOLATE
-  MODULE PROCEDURE INTERPOLATE4
-  MODULE PROCEDURE INTERPOLATE5
-  MODULE PROCEDURE INTERPOLATEL
-END INTERFACE
-
-INTERFACE SET
-  MODULE PROCEDURE SET3
-  MODULE PROCEDURE SET4
-  MODULE PROCEDURE SET5
-END INTERFACE
+USE ARRAYS_MANIP, ONLY: SETUP, XNAN, REPLICATE, NPROMIZE, INTERPOLATE, SET
 
 CONTAINS
 
@@ -35,8 +9,6 @@ SUBROUTINE GETDATA_RAIN_ICE (NPROMA, NGPBLKS, NFLEVG, LLMICRO_B, PEXNREF_B, PDZZ
   &PRS_B, PSIGS_B, PSEA_B, PTOWN_B, PCIT_OUT_B, PRS_OUT_B, ZINPRC_B, ZINPRC_OUT_B, PINPRR_B, PINPRR_OUT_B, PEVAP_B, PEVAP_OUT_B, &
   &PINPRS_B, PINPRS_OUT_B, PINPRG_B, PINPRG_OUT_B, ZINDEP_B, ZINDEP_OUT_B, ZRAINFR_B, ZRAINFR_OUT_B, PFPR_B, PFPR_OUT_B, LDVERBOSE)
 
-USE IEEE_ARITHMETIC, ONLY : IEEE_SIGNALING_NAN, IEEE_VALUE
-
 IMPLICIT NONE
 
 INTEGER, PARAMETER :: IFILE = 77
@@ -126,7 +98,8 @@ INTEGER :: NGPTOT, NPROMA, NGPBLKS, NFLEVG
 INTEGER :: IOFF, IBL
 LOGICAL :: LLEXIST
 CHARACTER(LEN=32) :: CLFILE
-REAL :: ZNAN
+
+CALL SETUP()
 
 KRR=6
 NGPTOT = NPROMA * NGPBLKS
@@ -180,9 +153,6 @@ ALLOCATE (PFPR_B          (NPROMA,1,NFLEVG,KRR,NGPBLKS))
 ALLOCATE (PFPR_OUT_B      (NPROMA,1,NFLEVG,KRR,NGPBLKS))
 
 
-ZNAN = IEEE_VALUE (ZNAN, IEEE_SIGNALING_NAN)
-
-
 !CALL SET (LLMICRO_B     )
 CALL SET (PEXNREF_B     )
 CALL SET (PDZZ_B        )
@@ -214,14 +184,14 @@ CALL SET (ZINDEP_OUT_B  )
 CALL SET (ZRAINFR_OUT_B )
 CALL SET (PFPR_OUT_B    )
 
-ZINPRC_OUT_B  =     ZNAN
-PINPRR_OUT_B  =     ZNAN
-PEVAP_OUT_B   =     ZNAN
-PINPRS_OUT_B  =     ZNAN
-PINPRG_OUT_B  =     ZNAN
-ZINDEP_OUT_B  =     ZNAN
-ZRAINFR_OUT_B =     ZNAN
-PFPR_OUT_B    =     ZNAN
+ZINPRC_OUT_B  =     XNAN
+PINPRR_OUT_B  =     XNAN
+PEVAP_OUT_B   =     XNAN
+PINPRS_OUT_B  =     XNAN
+PINPRG_OUT_B  =     XNAN
+ZINDEP_OUT_B  =     XNAN
+ZRAINFR_OUT_B =     XNAN
+PFPR_OUT_B    =     XNAN
 
 IOFF = 0
 IBL = 0
@@ -416,386 +386,4 @@ CALL NPROMIZE (NPROMA, PFPR_OUT    ,  PFPR_OUT_B      )
 
 END SUBROUTINE 
 
-SUBROUTINE REPLICATE4 (KOFF, P)
-
-INTEGER :: KOFF
-REAL    :: P (:,:,:,:)
-
-INTEGER :: I, J
-
-DO I = KOFF+1, SIZE (P, 1)
-  J = 1 + MODULO (I - 1, KOFF)
-  P (I, :, :, :) = P (J, :, :, :)
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE REPLICATE3 (KOFF, P)
-
-INTEGER :: KOFF
-REAL    :: P (:,:,:)
-
-INTEGER :: I, J
-
-DO I = KOFF+1, SIZE (P, 1)
-  J = 1 + MODULO (I - 1, KOFF)
-  P (I, :, :) = P (J, :, :)
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE REPLICATE2 (KOFF, P)
-
-INTEGER :: KOFF
-REAL    :: P (:,:)
-
-INTEGER :: I, J
-
-DO I = KOFF+1, SIZE (P, 1)
-  J = 1 + MODULO (I - 1, KOFF)
-  P (I, :) = P (J, :)
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE REPLICATEL (KOFF, L)
-
-INTEGER :: KOFF
-LOGICAL :: L (:,:,:)
-
-INTEGER :: I, J
-
-DO I = KOFF+1, SIZE (L, 1)
-  J = 1 + MODULO (I - 1, KOFF)
-  L (I, :, :) = L (J, :, :)
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE NPROMIZE3 (KPROMA, PI, PO)
-
-INTEGER :: KPROMA
-REAL, INTENT (IN)  :: PI (:,:,:) 
-REAL, INTENT (OUT) :: PO (:,:,:)
-
-INTEGER :: I, J, IGPBLK, IGPTOT, IGP, JLON, JIDIA, JFDIA
-
-IF (SIZE (PI, 3) /= 1) STOP 1
-
-IGPTOT = SIZE (PI, 1)
-IGPBLK = 1 + (IGPTOT-1) / KPROMA
-
-DO IGP = 1, IGPTOT, KPROMA
-  IBL = 1 + (IGP - 1) / KPROMA
-  JIDIA = 1
-  JFDIA = MIN (KPROMA, IGPTOT - (IBL - 1) * KPROMA)
-
-  DO JLON = JIDIA, JFDIA
-    PO (JLON, :, IBL) = PI (IGP + (JLON - 1), :, 1)
-  ENDDO
-
-  DO JLON = JFDIA+1, KPROMA
-    PO (JLON, :, IBL) = PO (JFDIA, :, IBL)
-  ENDDO
-
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE NPROMIZE4 (KPROMA, PI, PO)
-
-INTEGER :: KPROMA
-REAL, INTENT (IN)  :: PI (:,:,:,:) 
-REAL, INTENT (OUT) :: PO (:,:,:,:)
-
-INTEGER :: I, J, IGPBLK, IGPTOT, IGP, JLON, JIDIA, JFDIA
-
-IF (SIZE (PI, 4) /= 1) STOP 1
-
-IGPTOT = SIZE (PI, 1)
-IGPBLK = 1 + (IGPTOT-1) / KPROMA
-
-DO IGP = 1, IGPTOT, KPROMA
-  IBL = 1 + (IGP - 1) / KPROMA
-  JIDIA = 1
-  JFDIA = MIN (KPROMA, IGPTOT - (IBL - 1) * KPROMA)
-
-  DO JLON = JIDIA, JFDIA
-    PO (JLON, :, :, IBL) = PI (IGP + (JLON - 1), :, :, 1)
-  ENDDO
-
-  DO JLON = JFDIA+1, KPROMA
-    PO (JLON, :, :, IBL) = PO (JFDIA, :, :, IBL)
-  ENDDO
-
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE NPROMIZE5 (KPROMA, PI, PO)
-
-INTEGER :: KPROMA
-REAL, INTENT (IN)  :: PI (:,:,:,:,:) 
-REAL, INTENT (OUT) :: PO (:,:,:,:,:)
-
-INTEGER :: I, J, IGPBLK, IGPTOT, IGP, JLON, JIDIA, JFDIA
-
-IF (SIZE (PI, 5) /= 1) STOP 1
-
-IGPTOT = SIZE (PI, 1)
-IGPBLK = 1 + (IGPTOT-1) / KPROMA
-
-DO IGP = 1, IGPTOT, KPROMA
-  IBL = 1 + (IGP - 1) / KPROMA
-  JIDIA = 1
-  JFDIA = MIN (KPROMA, IGPTOT - (IBL - 1) * KPROMA)
-
-  DO JLON = JIDIA, JFDIA
-    PO (JLON, :, :, :, IBL) = PI (IGP + (JLON - 1), :, :, :, 1)
-  ENDDO
-
-  DO JLON = JFDIA+1, KPROMA
-    PO (JLON, :, :, :, IBL) = PI (JFDIA, :, :, :, IBL)
-  ENDDO
-
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE NPROMIZEL (KPROMA, LI, LO)
-
-INTEGER :: KPROMA
-LOGICAL, INTENT (IN)  :: LI (:,:,:,:) 
-LOGICAL, INTENT (OUT) :: LO (:,:,:,:)
-
-INTEGER :: I, J, IGPBLK, IGPTOT, IGP, JLON, JIDIA, JFDIA
-
-IF (SIZE (LI, 4) /= 1) STOP 1
-
-IGPTOT = SIZE (LI, 1)
-IGPBLK = 1 + (IGPTOT-1) / KPROMA
-
-DO IGP = 1, IGPTOT, KPROMA
-  IBL = 1 + (IGP - 1) / KPROMA
-  JIDIA = 1
-  JFDIA = MIN (KPROMA, IGPTOT - (IBL - 1) * KPROMA)
-
-  DO JLON = JIDIA, JFDIA
-    LO (JLON, :, :, IBL) = LI (IGP + (JLON - 1), :, :, 1)
-  ENDDO
-
-  DO JLON = JFDIA+1, KPROMA
-    LO (JLON, :, :, IBL) = LI (JFDIA, :, :, IBL)
-  ENDDO
-
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE INTERPOLATE4 (KFLEVG, KOFF, P)
-
-INTEGER :: KFLEVG, KOFF
-REAL, ALLOCATABLE :: P (:,:,:,:)
-REAL :: Z (LBOUND (P, 1):UBOUND (P, 1), &
-         & LBOUND (P, 2):UBOUND (P, 2), &
-         & LBOUND (P, 3):UBOUND (P, 3), &
-         & LBOUND (P, 4):UBOUND (P, 4))
-INTEGER :: ILEV1A, ILEV1B, ILEV2, NLEV1, NLEV2
-REAL :: ZWA, ZWB, ZLEV1, ZLEV2
-
-Z = P
-
-NLEV1 = SIZE (P, 3)
-NLEV2 = KFLEVG
-
-DEALLOCATE (P)
-
-ALLOCATE (P (LBOUND (Z, 1):UBOUND (Z, 1), &
-           & LBOUND (Z, 2):UBOUND (Z, 2), &
-           & KFLEVG, &
-           & LBOUND (Z, 4):UBOUND (Z, 4)))
-
-DO ILEV2 = 1, NLEV2
-  ZLEV2 = REAL (ILEV2 - 1) / REAL (NLEV2 -1)
-  ZLEV1 = 1. + ZLEV2 * REAL (NLEV1 - 1)
-  ILEV1B = MIN (CEILING (ZLEV1), NLEV1)
-  ILEV1A = MAX (FLOOR   (ZLEV1),     1)
-
-  IF (ILEV1A == ILEV1B) THEN
-    ZWA = 1.
-    ZWB = 0.
-  ELSE
-    ZWA = REAL (ILEV1B) - ZLEV1
-    ZWB = ZLEV1 - REAL (ILEV1A)
-  ENDIF
-
-! WRITE (*, '(" ZLEV2 = ",E12.5," ZLEV1 = ",E12.5," ILEV2 = ",I4," ILEV1A = ",I4," ZWA = ",E12.5," ILEV1B = ",I4," ZWB = ",E12.5)') &
-!   & ZLEV2, ZLEV1, ILEV2, ILEV1A, ZWA, ILEV1B, ZWB
-
-  P (1:KOFF, :, ILEV2, :) = ZWA * Z (1:KOFF, :, ILEV1A, :) + ZWB * Z (1:KOFF, :, ILEV1B, :) 
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE INTERPOLATE5 (KFLEVG, KOFF, P)
-
-INTEGER :: KFLEVG, KOFF
-REAL, ALLOCATABLE :: P (:,:,:,:,:)
-REAL :: Z (LBOUND (P, 1):UBOUND (P, 1), &
-         & LBOUND (P, 2):UBOUND (P, 2), &
-         & LBOUND (P, 3):UBOUND (P, 3), &
-         & LBOUND (P, 4):UBOUND (P, 4), &
-         & LBOUND (P, 5):UBOUND (P, 5))
-INTEGER :: ILEV1A, ILEV1B, ILEV2, NLEV1, NLEV2
-REAL :: ZWA, ZWB, ZLEV1, ZLEV2
-
-Z = P
-
-NLEV1 = SIZE (P, 3)
-NLEV2 = KFLEVG
-
-DEALLOCATE (P)
-
-ALLOCATE (P (LBOUND (Z, 1):UBOUND (Z, 1), &
-           & LBOUND (Z, 2):UBOUND (Z, 2), &
-           & KFLEVG, &
-           & LBOUND (Z, 4):UBOUND (Z, 4), &
-           & LBOUND (Z, 5):UBOUND (Z, 5)))
-
-DO ILEV2 = 1, NLEV2
-  ZLEV2 = REAL (ILEV2 - 1) / REAL (NLEV2 -1)
-  ZLEV1 = 1. + ZLEV2 * REAL (NLEV1 - 1)
-  ILEV1B = MIN (CEILING (ZLEV1), NLEV1)
-  ILEV1A = MAX (FLOOR   (ZLEV1),     1)
-
-  IF (ILEV1A == ILEV1B) THEN
-    ZWA = 1.
-    ZWB = 0.
-  ELSE
-    ZWA = REAL (ILEV1B) - ZLEV1
-    ZWB = ZLEV1 - REAL (ILEV1A)
-  ENDIF
-
-! WRITE (*, '(" ZLEV2 = ",E12.5," ZLEV1 = ",E12.5," ILEV2 = ",I4," ILEV1A = ",I4," ZWA = ",E12.5," ILEV1B = ",I4," ZWB = ",E12.5)') &
-!   & ZLEV2, ZLEV1, ILEV2, ILEV1A, ZWA, ILEV1B, ZWB
-
-  P (1:KOFF, :, ILEV2, :, :) = ZWA * Z (1:KOFF, :, ILEV1A, :, :) + ZWB * Z (1:KOFF, :, ILEV1B, :, :) 
-ENDDO
-
-END SUBROUTINE
-
-SUBROUTINE INTERPOLATEL (KFLEVG, KOFF, L)
-
-INTEGER :: KFLEVG, KOFF
-LOGICAL, ALLOCATABLE :: L (:,:,:,:)
-LOGICAL :: Z (LBOUND (L, 1):UBOUND (L, 1), &
-            & LBOUND (L, 2):UBOUND (L, 2), &
-            & LBOUND (L, 3):UBOUND (L, 3), &
-            & LBOUND (L, 4):UBOUND (L, 4))
-INTEGER :: ILEV1A, ILEV1B, ILEV2, NLEV1, NLEV2
-REAL :: ZWA, ZWB, ZLEV1, ZLEV2
-
-Z = L
-
-NLEV1 = SIZE (L, 3)
-NLEV2 = KFLEVG
-
-DEALLOCATE (L)
-
-ALLOCATE (L (LBOUND (Z, 1):UBOUND (Z, 1), &
-           & LBOUND (Z, 2):UBOUND (Z, 2), &
-           & KFLEVG, &
-           & LBOUND (Z, 4):UBOUND (Z, 4)))
-
-DO ILEV2 = 1, NLEV2
-  ZLEV2 = REAL (ILEV2 - 1) / REAL (NLEV2 -1)
-  ZLEV1 = 1. + ZLEV2 * REAL (NLEV1 - 1)
-  ILEV1B = MIN (CEILING (ZLEV1), NLEV1)
-  ILEV1A = MAX (FLOOR   (ZLEV1),     1)
-
-  IF (ILEV1A == ILEV1B) THEN
-    ZWA = 1.
-    ZWB = 0.
-  ELSE
-    ZWA = REAL (ILEV1B) - ZLEV1
-    ZWB = ZLEV1 - REAL (ILEV1A)
-  ENDIF
-
-! WRITE (*, '(" ZLEV2 = ",E12.5," ZLEV1 = ",E12.5," ILEV2 = ",I4," ILEV1A = ",I4," ZWA = ",E12.5," ILEV1B = ",I4," ZWB = ",E12.5)') &
-!   & ZLEV2, ZLEV1, ILEV2, ILEV1A, ZWA, ILEV1B, ZWB
-
-  L (1:KOFF, :, ILEV2, :) = ZWA * MERGE(1., 0., Z (1:KOFF, :, ILEV1A, :)) + ZWB * MERGE(1., 0., Z (1:KOFF, :, ILEV1B, :)) >= 0.5
-ENDDO
-
-END SUBROUTINE
-
-
-SUBROUTINE SET3 (P)
-
-REAL :: P (:,:,:)
-INTEGER :: IBL, IGPBLKS
-INTEGER :: NTID, ITID, JBLK1, JBLK2
-
-
-IGPBLKS = SIZE (P, 3)
-
-!$OMP PARALLEL PRIVATE (ITID, JBLK1, JBLK2, NTID)
-NTID = OMP_GET_MAX_THREADS ()
-ITID = OMP_GET_THREAD_NUM ()
-JBLK1 = 1 +  (IGPBLKS * (ITID+0)) / NTID
-JBLK2 =      (IGPBLKS * (ITID+1)) / NTID
-
-DO IBL = JBLK1, JBLK2
-  P (:,:,IBL) = ZNAN
-ENDDO
-
-!$OMP END PARALLEL
-
-END SUBROUTINE
-
-SUBROUTINE SET4 (P)
-
-REAL :: P (:,:,:,:)
-INTEGER :: IBL, IGPBLKS
-INTEGER :: NTID, ITID, JBLK1, JBLK2
-
-IGPBLKS = SIZE (P, 4)
-
-!$OMP PARALLEL PRIVATE (ITID, JBLK1, JBLK2, NTID)
-NTID = OMP_GET_MAX_THREADS ()
-ITID = OMP_GET_THREAD_NUM ()
-JBLK1 = 1 +  (IGPBLKS * (ITID+0)) / NTID
-JBLK2 =      (IGPBLKS * (ITID+1)) / NTID
-
-DO IBL = JBLK1, JBLK2
-  P (:,:,:,IBL) = ZNAN
-ENDDO
-
-!$OMP END PARALLEL
-
-END SUBROUTINE
-
-SUBROUTINE SET5 (P)
-
-REAL :: P (:,:,:,:,:)
-INTEGER :: IBL, IGPBLKS
-INTEGER :: NTID, ITID, JBLK1, JBLK2
-
-IGPBLKS = SIZE (P, 5)
-
-!$OMP PARALLEL PRIVATE (ITID, JBLK1, JBLK2, NTID)
-NTID = OMP_GET_MAX_THREADS ()
-ITID = OMP_GET_THREAD_NUM ()
-JBLK1 = 1 +  (IGPBLKS * (ITID+0)) / NTID
-JBLK2 =      (IGPBLKS * (ITID+1)) / NTID
-
-DO IBL = JBLK1, JBLK2
-  P (:,:,:,:,IBL) = ZNAN
-ENDDO
-
-!$OMP END PARALLEL
-
-END SUBROUTINE
-
-
 END  MODULE
diff --git a/src/testprogs/rain_ice/main_rain_ice.F90 b/src/testprogs/rain_ice/main_rain_ice.F90
index 7dc8676231d3a8a50a0356994e857f0f3a18a25a..1948e81660e4f0694aa774b9a63bc41861525a64 100644
--- a/src/testprogs/rain_ice/main_rain_ice.F90
+++ b/src/testprogs/rain_ice/main_rain_ice.F90
@@ -2,6 +2,7 @@ PROGRAM MAIN_RAIN_ICE
 
 USE XRD_GETOPTIONS
 USE GETDATA_RAIN_ICE_MOD
+USE COMPUTE_DIFF
 USE MODD_CONF
 USE MODD_DIMPHYEX,   ONLY: DIMPHYEX_t
 USE MODD_CST,        ONLY: CST_t, CST
@@ -304,20 +305,20 @@ IF (LLCHECK .OR. LLSTAT .OR. LLCHECKDIFF) THEN
     PRINT *, " IBL = ", IBL
     DO JRR=1, KRR
       WRITE (CLTEXT, '("PRS JRR=",I3.3)') JRR
-      CALL DIFF3 (CLTEXT,      PRS_OUT       (:,:,:,JRR,IBL), PRS      (:,:,:,JRR,IBL))
+      CALL DIFF3 (CLTEXT,      PRS_OUT       (:,:,:,JRR,IBL), PRS      (:,:,:,JRR,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
       IF(JRR>=2) THEN
         WRITE (CLTEXT, '("PFPR JRR=",I3.3)') JRR
-        CALL DIFF3 (CLTEXT,     PFPR_OUT      (:,:,:,JRR,IBL), PFPR     (:,:,:,JRR,IBL))
+        CALL DIFF3 (CLTEXT,     PFPR_OUT      (:,:,:,JRR,IBL), PFPR     (:,:,:,JRR,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
       ENDIF
     ENDDO
-    CALL DIFF3 ("PCIT",     PCIT_OUT      (:,:,:,IBL), PCIT     (:,:,:,IBL))
-    CALL DIFF2 ("ZINPRC",   ZINPRC_OUT    (:,:,IBL),   ZINPRC   (:,:,IBL))
-    CALL DIFF2 ("PINPRRRS", PINPRR_OUT    (:,:,IBL),   PINPRR   (:,:,IBL))
-    CALL DIFF3 ("PEVAP",    PEVAP_OUT     (:,:,:,IBL), PEVAP    (:,:,:,IBL))
-    CALL DIFF2 ("PINPRS",   PINPRS_OUT    (:,:,IBL),   PINPRS   (:,:,IBL))
-    CALL DIFF2 ("PINPRG",   PINPRG_OUT    (:,:,IBL),   PINPRG   (:,:,IBL))
-    CALL DIFF2 ("ZINDEP",   ZINDEP_OUT    (:,:,IBL),   ZINDEP   (:,:,IBL))
-    CALL DIFF3 ("ZRAINFR",  ZRAINFR_OUT   (:,:,:,IBL), ZRAINFR  (:,:,:,IBL))
+    CALL DIFF3 ("PCIT",     PCIT_OUT      (:,:,:,IBL), PCIT     (:,:,:,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
+    CALL DIFF2 ("ZINPRC",   ZINPRC_OUT    (:,:,IBL),   ZINPRC   (:,:,IBL)  , LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
+    CALL DIFF2 ("PINPRRRS", PINPRR_OUT    (:,:,IBL),   PINPRR   (:,:,IBL)  , LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
+    CALL DIFF3 ("PEVAP",    PEVAP_OUT     (:,:,:,IBL), PEVAP    (:,:,:,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
+    CALL DIFF2 ("PINPRS",   PINPRS_OUT    (:,:,IBL),   PINPRS   (:,:,IBL)  , LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
+    CALL DIFF2 ("PINPRG",   PINPRG_OUT    (:,:,IBL),   PINPRG   (:,:,IBL)  , LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
+    CALL DIFF2 ("ZINDEP",   ZINDEP_OUT    (:,:,IBL),   ZINDEP   (:,:,IBL)  , LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
+    CALL DIFF3 ("ZRAINFR",  ZRAINFR_OUT   (:,:,:,IBL), ZRAINFR  (:,:,:,IBL), LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
   ENDDO
 ENDIF
 
@@ -333,97 +334,6 @@ STOP
 
 CONTAINS
 
-SUBROUTINE DIFF3 (CDNAME, PREF, POUT)
-
-CHARACTER (LEN=*) :: CDNAME
-REAL :: PREF (:,:,:)
-REAL :: POUT (:,:,:)
-
-INTEGER :: JLON, JLEV
-
-PRINT *, CDNAME
-IF (LLSTAT) THEN
-  PRINT *, MINVAL (PREF), MAXVAL (PREF), SUM (PREF) / SIZE (PREF)
-  PRINT *, MINVAL (POUT), MAXVAL (POUT), SUM (POUT) / SIZE (POUT)
-ENDIF
-
-IF (LLCHECK) THEN
-  IF (SUM (ABS (POUT) + ABS (PREF)) > 0) THEN
-  WRITE (*, '(A4)', ADVANCE='NO') ""
-  DO JLON = 1, NPROMA
-    WRITE (*, '("|",I12,A12)', ADVANCE='NO') JLON, ""
-  ENDDO
-  WRITE (*, '("|")')
-  DO JLEV = 1, KLEV
-    WRITE (*, '(I4)', ADVANCE='NO') JLEV
-    DO JLON = 1, NPROMA
-      IF (ABS (PREF (JLON, 1, JLEV)) + ABS (POUT (JLON, 1, JLEV)) == 0.) THEN
-      WRITE (*, '("|",2A12)', ADVANCE='NO') "", ""
-      ELSE
-      WRITE (*, '("|",2E12.5)', ADVANCE='NO') PREF (JLON, 1, JLEV), POUT (JLON, 1, JLEV)
-      ENDIF
-    ENDDO
-    WRITE (*, '("|")')
-  ENDDO
-  ENDIF
-ENDIF
-
-IF (LLCHECKDIFF) THEN
-  IF (SUM(ABS(POUT-PREF)) > 0.) THEN
-    PRINT*, "THERE ARE DIFF"
-    LLDIFF = .TRUE.
-  ELSE
-    PRINT*, "THERE IS NO DIFF"
-  ENDIF
-ENDIF
-
-END SUBROUTINE
-
-SUBROUTINE DIFF2 (CDNAME, PREF, POUT)
-
-CHARACTER (LEN=*) :: CDNAME
-REAL :: PREF (:,:)
-REAL :: POUT (:,:)
-
-INTEGER :: JLON
-
-PRINT *, CDNAME
-IF (LLSTAT) THEN
-  PRINT *, MINVAL (PREF), MAXVAL (PREF), SUM (PREF) / SIZE (PREF)
-  PRINT *, MINVAL (POUT), MAXVAL (POUT), SUM (POUT) / SIZE (POUT)
-ENDIF
-
-IF (LLCHECK) THEN
-  IF (SUM (ABS (POUT) + ABS (PREF)) > 0) THEN
-  WRITE (*, '(A4)', ADVANCE='NO') ""
-  DO JLON = 1, NPROMA
-    WRITE (*, '("|",I12,A12)', ADVANCE='NO') JLON, ""
-  ENDDO
-  WRITE (*, '("|")')
-  WRITE (*, '(I4)', ADVANCE='NO') 0
-  DO JLON = 1, NPROMA
-    IF (ABS (PREF (JLON, 1)) + ABS (POUT (JLON, 1)) == 0.) THEN
-    WRITE (*, '("|",2A12)', ADVANCE='NO') "", ""
-    ELSE
-    WRITE (*, '("|",2E12.5)', ADVANCE='NO') PREF (JLON, 1), POUT (JLON, 1)
-    ENDIF
-  ENDDO
-  WRITE (*, '("|")')
-  ENDIF
-ENDIF
-
-IF (LLCHECKDIFF) THEN
-  IF (SUM(ABS(POUT-PREF)) > 0.) THEN
-    PRINT*, "THERE ARE DIFF"
-    LLDIFF = .TRUE.
-  ELSE
-    PRINT*, "THERE IS NO DIFF"
-  ENDIF
-ENDIF
-
-END SUBROUTINE
-
-
 SUBROUTINE INIT_PHYEX(KULOUT,LDWARM,CMICRO,CCSEDIM,LDCRIAUTI,&
                    PCRIAUTI,PT0CRIAUTI,PCRIAUTC)
 
diff --git a/src/testprogs/support/arrays_manip.F90 b/src/testprogs/support/arrays_manip.F90
new file mode 100644
index 0000000000000000000000000000000000000000..3f659a613ff0a7719c36c7f04d678a9993fbd6cc
--- /dev/null
+++ b/src/testprogs/support/arrays_manip.F90
@@ -0,0 +1,435 @@
+MODULE ARRAYS_MANIP
+
+USE OMP_LIB
+USE IEEE_ARITHMETIC, ONLY : IEEE_SIGNALING_NAN, IEEE_VALUE
+
+INTERFACE REPLICATE
+  MODULE PROCEDURE REPLICATE2
+  MODULE PROCEDURE REPLICATE3
+  MODULE PROCEDURE REPLICATE4
+  MODULE PROCEDURE REPLICATEL
+END INTERFACE
+
+INTERFACE NPROMIZE
+  MODULE PROCEDURE NPROMIZE3
+  MODULE PROCEDURE NPROMIZE4
+  MODULE PROCEDURE NPROMIZE5
+  MODULE PROCEDURE NPROMIZEL
+END INTERFACE                   
+
+INTERFACE INTERPOLATE
+  MODULE PROCEDURE INTERPOLATE4
+  MODULE PROCEDURE INTERPOLATE5
+  MODULE PROCEDURE INTERPOLATEL
+END INTERFACE
+
+INTERFACE SET
+  MODULE PROCEDURE SET3
+  MODULE PROCEDURE SET4
+  MODULE PROCEDURE SET5
+END INTERFACE
+
+REAL, SAVE :: XNAN
+
+CONTAINS
+
+SUBROUTINE SETUP()
+XNAN = IEEE_VALUE (1., IEEE_SIGNALING_NAN)
+END SUBROUTINE SETUP
+
+SUBROUTINE REPLICATE4 (KOFF, P)
+IMPLICIT NONE
+
+INTEGER :: KOFF
+REAL    :: P (:,:,:,:)
+
+INTEGER :: I, J
+
+DO I = KOFF+1, SIZE (P, 1)
+  J = 1 + MODULO (I - 1, KOFF)
+  P (I, :, :, :) = P (J, :, :, :)
+ENDDO
+
+END SUBROUTINE
+
+SUBROUTINE REPLICATE3 (KOFF, P)
+IMPLICIT NONE
+
+INTEGER :: KOFF
+REAL    :: P (:,:,:)
+
+INTEGER :: I, J
+
+DO I = KOFF+1, SIZE (P, 1)
+  J = 1 + MODULO (I - 1, KOFF)
+  P (I, :, :) = P (J, :, :)
+ENDDO
+
+END SUBROUTINE
+
+SUBROUTINE REPLICATE2 (KOFF, P)
+IMPLICIT NONE
+
+INTEGER :: KOFF
+REAL    :: P (:,:)
+
+INTEGER :: I, J
+
+DO I = KOFF+1, SIZE (P, 1)
+  J = 1 + MODULO (I - 1, KOFF)
+  P (I, :) = P (J, :)
+ENDDO
+
+END SUBROUTINE
+
+SUBROUTINE REPLICATEL (KOFF, L)
+IMPLICIT NONE
+
+INTEGER :: KOFF
+LOGICAL :: L (:,:,:)
+
+INTEGER :: I, J
+
+DO I = KOFF+1, SIZE (L, 1)
+  J = 1 + MODULO (I - 1, KOFF)
+  L (I, :, :) = L (J, :, :)
+ENDDO
+
+END SUBROUTINE
+
+SUBROUTINE NPROMIZE3 (KPROMA, PI, PO)
+IMPLICIT NONE
+
+INTEGER :: KPROMA
+REAL, INTENT (IN)  :: PI (:,:,:) 
+REAL, INTENT (OUT) :: PO (:,:,:)
+
+INTEGER :: I, J, IGPBLK, IGPTOT, IGP, JLON, JIDIA, JFDIA, IBL
+
+IF (SIZE (PI, 3) /= 1) STOP 1
+
+IGPTOT = SIZE (PI, 1)
+IGPBLK = 1 + (IGPTOT-1) / KPROMA
+
+DO IGP = 1, IGPTOT, KPROMA
+  IBL = 1 + (IGP - 1) / KPROMA
+  JIDIA = 1
+  JFDIA = MIN (KPROMA, IGPTOT - (IBL - 1) * KPROMA)
+
+  DO JLON = JIDIA, JFDIA
+    PO (JLON, :, IBL) = PI (IGP + (JLON - 1), :, 1)
+  ENDDO
+
+  DO JLON = JFDIA+1, KPROMA
+    PO (JLON, :, IBL) = PO (JFDIA, :, IBL)
+  ENDDO
+
+ENDDO
+
+END SUBROUTINE
+
+SUBROUTINE NPROMIZE4 (KPROMA, PI, PO)
+IMPLICIT NONE
+
+INTEGER :: KPROMA
+REAL, INTENT (IN)  :: PI (:,:,:,:) 
+REAL, INTENT (OUT) :: PO (:,:,:,:)
+
+INTEGER :: I, J, IGPBLK, IGPTOT, IGP, JLON, JIDIA, JFDIA, IBL
+
+IF (SIZE (PI, 4) /= 1) STOP 1
+
+IGPTOT = SIZE (PI, 1)
+IGPBLK = 1 + (IGPTOT-1) / KPROMA
+
+DO IGP = 1, IGPTOT, KPROMA
+  IBL = 1 + (IGP - 1) / KPROMA
+  JIDIA = 1
+  JFDIA = MIN (KPROMA, IGPTOT - (IBL - 1) * KPROMA)
+
+  DO JLON = JIDIA, JFDIA
+    PO (JLON, :, :, IBL) = PI (IGP + (JLON - 1), :, :, 1)
+  ENDDO
+
+  DO JLON = JFDIA+1, KPROMA
+    PO (JLON, :, :, IBL) = PO (JFDIA, :, :, IBL)
+  ENDDO
+
+ENDDO
+
+END SUBROUTINE
+
+SUBROUTINE NPROMIZE5 (KPROMA, PI, PO)
+IMPLICIT NONE
+
+INTEGER :: KPROMA
+REAL, INTENT (IN)  :: PI (:,:,:,:,:) 
+REAL, INTENT (OUT) :: PO (:,:,:,:,:)
+
+INTEGER :: I, J, IGPBLK, IGPTOT, IGP, JLON, JIDIA, JFDIA, IBL
+
+IF (SIZE (PI, 5) /= 1) STOP 1
+
+IGPTOT = SIZE (PI, 1)
+IGPBLK = 1 + (IGPTOT-1) / KPROMA
+
+DO IGP = 1, IGPTOT, KPROMA
+  IBL = 1 + (IGP - 1) / KPROMA
+  JIDIA = 1
+  JFDIA = MIN (KPROMA, IGPTOT - (IBL - 1) * KPROMA)
+
+  DO JLON = JIDIA, JFDIA
+    PO (JLON, :, :, :, IBL) = PI (IGP + (JLON - 1), :, :, :, 1)
+  ENDDO
+
+  DO JLON = JFDIA+1, KPROMA
+    PO (JLON, :, :, :, IBL) = PI (JFDIA, :, :, :, IBL)
+  ENDDO
+
+ENDDO
+
+END SUBROUTINE
+
+SUBROUTINE NPROMIZEL (KPROMA, LI, LO)
+IMPLICIT NONE
+
+INTEGER :: KPROMA
+LOGICAL, INTENT (IN)  :: LI (:,:,:,:) 
+LOGICAL, INTENT (OUT) :: LO (:,:,:,:)
+
+INTEGER :: I, J, IGPBLK, IGPTOT, IGP, JLON, JIDIA, JFDIA, IBL
+
+IF (SIZE (LI, 4) /= 1) STOP 1
+
+IGPTOT = SIZE (LI, 1)
+IGPBLK = 1 + (IGPTOT-1) / KPROMA
+
+DO IGP = 1, IGPTOT, KPROMA
+  IBL = 1 + (IGP - 1) / KPROMA
+  JIDIA = 1
+  JFDIA = MIN (KPROMA, IGPTOT - (IBL - 1) * KPROMA)
+
+  DO JLON = JIDIA, JFDIA
+    LO (JLON, :, :, IBL) = LI (IGP + (JLON - 1), :, :, 1)
+  ENDDO
+
+  DO JLON = JFDIA+1, KPROMA
+    LO (JLON, :, :, IBL) = LI (JFDIA, :, :, IBL)
+  ENDDO
+
+ENDDO
+
+END SUBROUTINE
+
+SUBROUTINE INTERPOLATE4 (KFLEVG, KOFF, P)
+IMPLICIT NONE
+
+INTEGER :: KFLEVG, KOFF
+REAL, ALLOCATABLE :: P (:,:,:,:)
+REAL :: Z (LBOUND (P, 1):UBOUND (P, 1), &
+         & LBOUND (P, 2):UBOUND (P, 2), &
+         & LBOUND (P, 3):UBOUND (P, 3), &
+         & LBOUND (P, 4):UBOUND (P, 4))
+INTEGER :: ILEV1A, ILEV1B, ILEV2, NLEV1, NLEV2
+REAL :: ZWA, ZWB, ZLEV1, ZLEV2
+
+Z = P
+
+NLEV1 = SIZE (P, 3)
+NLEV2 = KFLEVG
+
+DEALLOCATE (P)
+
+ALLOCATE (P (LBOUND (Z, 1):UBOUND (Z, 1), &
+           & LBOUND (Z, 2):UBOUND (Z, 2), &
+           & KFLEVG, &
+           & LBOUND (Z, 4):UBOUND (Z, 4)))
+
+DO ILEV2 = 1, NLEV2
+  ZLEV2 = REAL (ILEV2 - 1) / REAL (NLEV2 -1)
+  ZLEV1 = 1. + ZLEV2 * REAL (NLEV1 - 1)
+  ILEV1B = MIN (CEILING (ZLEV1), NLEV1)
+  ILEV1A = MAX (FLOOR   (ZLEV1),     1)
+
+  IF (ILEV1A == ILEV1B) THEN
+    ZWA = 1.
+    ZWB = 0.
+  ELSE
+    ZWA = REAL (ILEV1B) - ZLEV1
+    ZWB = ZLEV1 - REAL (ILEV1A)
+  ENDIF
+
+! WRITE (*, '(" ZLEV2 = ",E12.5," ZLEV1 = ",E12.5," ILEV2 = ",I4," ILEV1A = ",I4," ZWA = ",E12.5," ILEV1B = ",I4," ZWB = ",E12.5)') &
+!   & ZLEV2, ZLEV1, ILEV2, ILEV1A, ZWA, ILEV1B, ZWB
+
+  P (1:KOFF, :, ILEV2, :) = ZWA * Z (1:KOFF, :, ILEV1A, :) + ZWB * Z (1:KOFF, :, ILEV1B, :) 
+ENDDO
+
+END SUBROUTINE
+
+SUBROUTINE INTERPOLATE5 (KFLEVG, KOFF, P)
+IMPLICIT NONE
+
+INTEGER :: KFLEVG, KOFF
+REAL, ALLOCATABLE :: P (:,:,:,:,:)
+REAL :: Z (LBOUND (P, 1):UBOUND (P, 1), &
+         & LBOUND (P, 2):UBOUND (P, 2), &
+         & LBOUND (P, 3):UBOUND (P, 3), &
+         & LBOUND (P, 4):UBOUND (P, 4), &
+         & LBOUND (P, 5):UBOUND (P, 5))
+INTEGER :: ILEV1A, ILEV1B, ILEV2, NLEV1, NLEV2
+REAL :: ZWA, ZWB, ZLEV1, ZLEV2
+
+Z = P
+
+NLEV1 = SIZE (P, 3)
+NLEV2 = KFLEVG
+
+DEALLOCATE (P)
+
+ALLOCATE (P (LBOUND (Z, 1):UBOUND (Z, 1), &
+           & LBOUND (Z, 2):UBOUND (Z, 2), &
+           & KFLEVG, &
+           & LBOUND (Z, 4):UBOUND (Z, 4), &
+           & LBOUND (Z, 5):UBOUND (Z, 5)))
+
+DO ILEV2 = 1, NLEV2
+  ZLEV2 = REAL (ILEV2 - 1) / REAL (NLEV2 -1)
+  ZLEV1 = 1. + ZLEV2 * REAL (NLEV1 - 1)
+  ILEV1B = MIN (CEILING (ZLEV1), NLEV1)
+  ILEV1A = MAX (FLOOR   (ZLEV1),     1)
+
+  IF (ILEV1A == ILEV1B) THEN
+    ZWA = 1.
+    ZWB = 0.
+  ELSE
+    ZWA = REAL (ILEV1B) - ZLEV1
+    ZWB = ZLEV1 - REAL (ILEV1A)
+  ENDIF
+
+! WRITE (*, '(" ZLEV2 = ",E12.5," ZLEV1 = ",E12.5," ILEV2 = ",I4," ILEV1A = ",I4," ZWA = ",E12.5," ILEV1B = ",I4," ZWB = ",E12.5)') &
+!   & ZLEV2, ZLEV1, ILEV2, ILEV1A, ZWA, ILEV1B, ZWB
+
+  P (1:KOFF, :, ILEV2, :, :) = ZWA * Z (1:KOFF, :, ILEV1A, :, :) + ZWB * Z (1:KOFF, :, ILEV1B, :, :) 
+ENDDO
+
+END SUBROUTINE
+
+SUBROUTINE INTERPOLATEL (KFLEVG, KOFF, L)
+IMPLICIT NONE
+
+INTEGER :: KFLEVG, KOFF
+LOGICAL, ALLOCATABLE :: L (:,:,:,:)
+LOGICAL :: Z (LBOUND (L, 1):UBOUND (L, 1), &
+            & LBOUND (L, 2):UBOUND (L, 2), &
+            & LBOUND (L, 3):UBOUND (L, 3), &
+            & LBOUND (L, 4):UBOUND (L, 4))
+INTEGER :: ILEV1A, ILEV1B, ILEV2, NLEV1, NLEV2
+REAL :: ZWA, ZWB, ZLEV1, ZLEV2
+
+Z = L
+
+NLEV1 = SIZE (L, 3)
+NLEV2 = KFLEVG
+
+DEALLOCATE (L)
+
+ALLOCATE (L (LBOUND (Z, 1):UBOUND (Z, 1), &
+           & LBOUND (Z, 2):UBOUND (Z, 2), &
+           & KFLEVG, &
+           & LBOUND (Z, 4):UBOUND (Z, 4)))
+
+DO ILEV2 = 1, NLEV2
+  ZLEV2 = REAL (ILEV2 - 1) / REAL (NLEV2 -1)
+  ZLEV1 = 1. + ZLEV2 * REAL (NLEV1 - 1)
+  ILEV1B = MIN (CEILING (ZLEV1), NLEV1)
+  ILEV1A = MAX (FLOOR   (ZLEV1),     1)
+
+  IF (ILEV1A == ILEV1B) THEN
+    ZWA = 1.
+    ZWB = 0.
+  ELSE
+    ZWA = REAL (ILEV1B) - ZLEV1
+    ZWB = ZLEV1 - REAL (ILEV1A)
+  ENDIF
+
+! WRITE (*, '(" ZLEV2 = ",E12.5," ZLEV1 = ",E12.5," ILEV2 = ",I4," ILEV1A = ",I4," ZWA = ",E12.5," ILEV1B = ",I4," ZWB = ",E12.5)') &
+!   & ZLEV2, ZLEV1, ILEV2, ILEV1A, ZWA, ILEV1B, ZWB
+
+  L (1:KOFF, :, ILEV2, :) = ZWA * MERGE(1., 0., Z (1:KOFF, :, ILEV1A, :)) + ZWB * MERGE(1., 0., Z (1:KOFF, :, ILEV1B, :)) >= 0.5
+ENDDO
+
+END SUBROUTINE
+
+
+SUBROUTINE SET3 (P)
+IMPLICIT NONE
+
+REAL :: P (:,:,:)
+INTEGER :: IBL, IGPBLKS
+INTEGER :: NTID, ITID, JBLK1, JBLK2
+
+
+IGPBLKS = SIZE (P, 3)
+
+!$OMP PARALLEL PRIVATE (ITID, JBLK1, JBLK2, NTID)
+NTID = OMP_GET_MAX_THREADS ()
+ITID = OMP_GET_THREAD_NUM ()
+JBLK1 = 1 +  (IGPBLKS * (ITID+0)) / NTID
+JBLK2 =      (IGPBLKS * (ITID+1)) / NTID
+
+DO IBL = JBLK1, JBLK2
+  P (:,:,IBL) = XNAN
+ENDDO
+
+!$OMP END PARALLEL
+
+END SUBROUTINE
+
+SUBROUTINE SET4 (P)
+IMPLICIT NONE
+
+REAL :: P (:,:,:,:)
+INTEGER :: IBL, IGPBLKS
+INTEGER :: NTID, ITID, JBLK1, JBLK2
+
+IGPBLKS = SIZE (P, 4)
+
+!$OMP PARALLEL PRIVATE (ITID, JBLK1, JBLK2, NTID)
+NTID = OMP_GET_MAX_THREADS ()
+ITID = OMP_GET_THREAD_NUM ()
+JBLK1 = 1 +  (IGPBLKS * (ITID+0)) / NTID
+JBLK2 =      (IGPBLKS * (ITID+1)) / NTID
+
+DO IBL = JBLK1, JBLK2
+  P (:,:,:,IBL) = XNAN
+ENDDO
+
+!$OMP END PARALLEL
+
+END SUBROUTINE
+
+SUBROUTINE SET5 (P)
+IMPLICIT NONE
+
+REAL :: P (:,:,:,:,:)
+INTEGER :: IBL, IGPBLKS
+INTEGER :: NTID, ITID, JBLK1, JBLK2
+
+IGPBLKS = SIZE (P, 5)
+
+!$OMP PARALLEL PRIVATE (ITID, JBLK1, JBLK2, NTID)
+NTID = OMP_GET_MAX_THREADS ()
+ITID = OMP_GET_THREAD_NUM ()
+JBLK1 = 1 +  (IGPBLKS * (ITID+0)) / NTID
+JBLK2 =      (IGPBLKS * (ITID+1)) / NTID
+
+DO IBL = JBLK1, JBLK2
+  P (:,:,:,:,IBL) = XNAN
+ENDDO
+
+!$OMP END PARALLEL
+
+END SUBROUTINE
+
+END MODULE ARRAYS_MANIP
diff --git a/src/testprogs/support/diff.F90 b/src/testprogs/support/diff.F90
new file mode 100644
index 0000000000000000000000000000000000000000..68fade755a33ec3b94800cb67f698e328310a3a7
--- /dev/null
+++ b/src/testprogs/support/diff.F90
@@ -0,0 +1,110 @@
+MODULE COMPUTE_DIFF
+
+INTERFACE DIFF
+  MODULE PROCEDURE DIFF3
+  MODULE PROCEDURE DIFF2
+END INTERFACE DIFF
+
+CONTAINS
+
+SUBROUTINE DIFF3 (CDNAME, PREF, POUT, LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
+IMPLICIT NONE
+
+CHARACTER (LEN=*) :: CDNAME
+REAL, INTENT(IN) :: PREF (:,:,:)
+REAL, INTENT(IN) :: POUT (:,:,:)
+LOGICAL, INTENT(IN) :: LLSTAT, LLCHECK, LLCHECKDIFF
+INTEGER, INTENT(IN) :: NPROMA
+LOGICAL, INTENT(OUT) :: LLDIFF
+
+INTEGER :: JLON, JLEV, KLEV
+
+KLEV=SIZE(PREF, 3)
+
+PRINT *, CDNAME
+IF (LLSTAT) THEN
+  PRINT *, MINVAL (PREF), MAXVAL (PREF), SUM (PREF) / SIZE (PREF)
+  PRINT *, MINVAL (POUT), MAXVAL (POUT), SUM (POUT) / SIZE (POUT)
+ENDIF
+
+IF (LLCHECK) THEN
+  IF (SUM (ABS (POUT) + ABS (PREF)) > 0) THEN
+  WRITE (*, '(A4)', ADVANCE='NO') ""
+  DO JLON = 1, NPROMA
+    WRITE (*, '("|",I12,A12)', ADVANCE='NO') JLON, ""
+  ENDDO
+  WRITE (*, '("|")')
+  DO JLEV = 1, KLEV
+    WRITE (*, '(I4)', ADVANCE='NO') JLEV
+    DO JLON = 1, NPROMA
+      IF (ABS (PREF (JLON, 1, JLEV)) + ABS (POUT (JLON, 1, JLEV)) == 0.) THEN
+      WRITE (*, '("|",2A12)', ADVANCE='NO') "", ""
+      ELSE
+      WRITE (*, '("|",2E12.5)', ADVANCE='NO') PREF (JLON, 1, JLEV), POUT (JLON, 1, JLEV)
+      ENDIF
+    ENDDO
+    WRITE (*, '("|")')
+  ENDDO
+  ENDIF
+ENDIF
+
+IF (LLCHECKDIFF) THEN
+  IF (SUM(ABS(POUT-PREF)) > 0.) THEN
+    PRINT*, "THERE ARE DIFF"
+    LLDIFF = .TRUE.
+  ELSE
+    PRINT*, "THERE IS NO DIFF"
+  ENDIF
+ENDIF
+
+END SUBROUTINE
+
+SUBROUTINE DIFF2 (CDNAME, PREF, POUT, LLSTAT, LLCHECK, NPROMA, LLCHECKDIFF, LLDIFF)
+IMPLICIT NONE
+
+CHARACTER (LEN=*) :: CDNAME
+REAL, INTENT(IN) :: PREF (:,:)
+REAL, INTENT(IN) :: POUT (:,:)
+LOGICAL, INTENT(IN) :: LLSTAT, LLCHECK, LLCHECKDIFF
+INTEGER, INTENT(IN) :: NPROMA
+LOGICAL, INTENT(OUT) :: LLDIFF
+
+INTEGER :: JLON
+
+PRINT *, CDNAME
+IF (LLSTAT) THEN
+  PRINT *, MINVAL (PREF), MAXVAL (PREF), SUM (PREF) / SIZE (PREF)
+  PRINT *, MINVAL (POUT), MAXVAL (POUT), SUM (POUT) / SIZE (POUT)
+ENDIF
+
+IF (LLCHECK) THEN
+  IF (SUM (ABS (POUT) + ABS (PREF)) > 0) THEN
+  WRITE (*, '(A4)', ADVANCE='NO') ""
+  DO JLON = 1, NPROMA
+    WRITE (*, '("|",I12,A12)', ADVANCE='NO') JLON, ""
+  ENDDO
+  WRITE (*, '("|")')
+  WRITE (*, '(I4)', ADVANCE='NO') 0
+  DO JLON = 1, NPROMA
+    IF (ABS (PREF (JLON, 1)) + ABS (POUT (JLON, 1)) == 0.) THEN
+    WRITE (*, '("|",2A12)', ADVANCE='NO') "", ""
+    ELSE
+    WRITE (*, '("|",2E12.5)', ADVANCE='NO') PREF (JLON, 1), POUT (JLON, 1)
+    ENDIF
+  ENDDO
+  WRITE (*, '("|")')
+  ENDIF
+ENDIF
+
+IF (LLCHECKDIFF) THEN
+  IF (SUM(ABS(POUT-PREF)) > 0.) THEN
+    PRINT*, "THERE ARE DIFF"
+    LLDIFF = .TRUE.
+  ELSE
+    PRINT*, "THERE IS NO DIFF"
+  ENDIF
+ENDIF
+
+END SUBROUTINE
+
+END MODULE COMPUTE_DIFF