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
