! $Id: extrapol.f90 5159 2024-08-02 19:58:25Z dcugnet $



SUBROUTINE extrapol(pfild, kxlon, kylat, pmask, &
        norsud, ldper, knbor, pwork)
  IMPLICIT NONE

  ! OASIS routine (Adaptation: Laurent Li, le 14 mars 1997)
  ! Fill up missed values by using the neighbor points

  INTEGER :: kxlon, kylat ! longitude and latitude dimensions (Input)
  INTEGER :: knbor ! minimum neighbor number (Input)
  LOGICAL :: norsud ! True if field is from North to South (Input)
  LOGICAL :: ldper ! True if take into account the periodicity (Input)
  REAL :: pmask ! mask value (Input)
  REAL :: pfild(kxlon,kylat) ! field to be extrapolated (Input/Output)
  REAL :: pwork(kxlon,kylat) ! working space

  REAL :: zwmsk
  INTEGER :: incre, idoit, i, j, k, inbor, ideb, ifin, ilon, jlat
  INTEGER :: ix(9), jy(9) ! index arrays for the neighbors coordinates
  REAL :: zmask(9)

  !  We search over the eight closest neighbors

  !        j+1  7  8  9
  !          j  4  5  6    Current point 5 --> (i,j)
  !        j-1  1  2  3
  !            i-1 i i+1


  IF (norsud) THEN
     DO j = 1, kylat
     DO i = 1, kxlon
        pwork(i,j) = pfild(i,kylat-j+1)
     ENDDO
     ENDDO
     DO j = 1, kylat
     DO i = 1, kxlon
        pfild(i,j) = pwork(i,j)
     ENDDO
     ENDDO
  ENDIF

  incre = 0

  DO j = 1, kylat
  DO i = 1, kxlon
     pwork(i,j) = pfild(i,j)
  ENDDO
  ENDDO

  !* To avoid problems in floating point tests
  zwmsk = pmask - 1.0

200   CONTINUE
  incre = incre + 1
  DO j = 1, kylat
  DO i = 1, kxlon
  IF (pfild(i,j)> zwmsk) THEN
     pwork(i,j) = pfild(i,j)
     inbor = 0
     ideb = 1
     ifin = 9

  !* Fill up ix array
     ix(1) = MAX (1,i-1)
     ix(2) = i
     ix(3) = MIN (kxlon,i+1)
     ix(4) = MAX (1,i-1)
     ix(5) = i
     ix(6) = MIN (kxlon,i+1)
     ix(7) = MAX (1,i-1)
     ix(8) = i
     ix(9) = MIN (kxlon,i+1)

  !* Fill up iy array
     jy(1) = MAX (1,j-1)
     jy(2) = MAX (1,j-1)
     jy(3) = MAX (1,j-1)
     jy(4) = j
     jy(5) = j
     jy(6) = j
     jy(7) = MIN (kylat,j+1)
     jy(8) = MIN (kylat,j+1)
     jy(9) = MIN (kylat,j+1)

  !* Correct latitude bounds if southernmost or northernmost points
     IF (j == 1) ideb = 4
     IF (j == kylat) ifin = 6

  !* Account for periodicity in longitude

     IF (ldper) THEN
        IF (i == kxlon) THEN
           ix(3) = 1
           ix(6) = 1
           ix(9) = 1
        ELSE IF (i == 1) THEN
           ix(1) = kxlon
           ix(4) = kxlon
           ix(7) = kxlon
        ENDIF
     ELSE
        IF (i == 1) THEN
           ix(1) = i
           ix(2) = i + 1
           ix(3) = i
           ix(4) = i + 1
           ix(5) = i
           ix(6) = i + 1
        ENDIF
        IF (i == kxlon) THEN
           ix(1) = i -1
           ix(2) = i
           ix(3) = i - 1
           ix(4) = i
           ix(5) = i - 1
           ix(6) = i
        ENDIF

        IF (i == 1 .OR. i == kxlon) THEN
           jy(1) = MAX (1,j-1)
           jy(2) = MAX (1,j-1)
           jy(3) = j
           jy(4) = j
           jy(5) = MIN (kylat,j+1)
           jy(6) = MIN (kylat,j+1)

           ideb = 1
           ifin = 6
           IF (j == 1) ideb = 3
           IF (j == kylat) ifin = 4
        ENDIF
     ENDIF ! end for ldper test

  !* Find unmasked neighbors

     DO k = ideb, ifin
        zmask(k) = 0.
        ilon = ix(k)
        jlat = jy(k)
        IF (pfild(ilon,jlat) < zwmsk) THEN
           zmask(k) = 1.
           inbor = inbor + 1
        ENDIF
  END DO

  !* Not enough points around point P are unmasked; interpolation on P
  !  will be done in a future CALL to extrap.

     IF (inbor >= knbor) THEN
        pwork(i,j) = 0.
        DO k = ideb, ifin
           ilon = ix(k)
           jlat = jy(k)
           pwork(i,j) = pwork(i,j) &
                 + pfild(ilon,jlat) * zmask(k)/ REAL(inbor)
        ENDDO
     ENDIF

  ENDIF
  END DO
  END DO

  !*    3. Writing back unmasked field in pfild
  !    ------------------------------------

  !* pfild then contains:
  ! - Values which were not masked
  ! - Interpolated values from the inbor neighbors
  ! - Values which are not yet interpolated

  idoit = 0
  DO j = 1, kylat
  DO i = 1, kxlon
     IF (pwork(i,j) > zwmsk) idoit = idoit + 1
     pfild(i,j) = pwork(i,j)
  ENDDO
  ENDDO

  IF (idoit /= 0) GOTO 200
  !cc      PRINT*, "Number of extrapolation steps incre =", incre

  IF (norsud) THEN
     DO j = 1, kylat
     DO i = 1, kxlon
        pwork(i,j) = pfild(i,kylat-j+1)
     ENDDO
     ENDDO
     DO j = 1, kylat
     DO i = 1, kxlon
        pfild(i,j) = pwork(i,j)
     ENDDO
     ENDDO
  ENDIF

  RETURN
END SUBROUTINE extrapol