! ! $Id: extrapol.f90 5246 2024-10-21 12:58:45Z abarral $ ! ! ! 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).GT. 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 .EQ. 1) ideb = 4 IF (j .EQ. kylat) ifin = 6 ! !* Account for periodicity in longitude ! IF (ldper) THEN IF (i .EQ. kxlon) THEN ix(3) = 1 ix(6) = 1 ix(9) = 1 ELSE IF (i .EQ. 1) THEN ix(1) = kxlon ix(4) = kxlon ix(7) = kxlon ENDIF ELSE IF (i .EQ. 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 .EQ. 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 .EQ. 1 .OR. i .EQ. 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 .EQ. 1) ideb = 3 IF (j .EQ. 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) .LT. 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 .GE. 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) .GT. zwmsk) idoit = idoit + 1 pfild(i,j) = pwork(i,j) ENDDO ENDDO ! IF (idoit .ne. 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