| 1 | ! |
|---|
| 2 | ! $Id: extrapol.F 1403 2010-07-01 09:02:53Z fairhead $ |
|---|
| 3 | ! |
|---|
| 4 | C |
|---|
| 5 | C |
|---|
| 6 | SUBROUTINE extrapol (pfild, kxlon, kylat, pmask, |
|---|
| 7 | . norsud, ldper, knbor, pwork) |
|---|
| 8 | IMPLICIT none |
|---|
| 9 | c |
|---|
| 10 | c OASIS routine (Adaptation: Laurent Li, le 14 mars 1997) |
|---|
| 11 | c Fill up missed values by using the neighbor points |
|---|
| 12 | c |
|---|
| 13 | INTEGER kxlon, kylat ! longitude and latitude dimensions (Input) |
|---|
| 14 | INTEGER knbor ! minimum neighbor number (Input) |
|---|
| 15 | LOGICAL norsud ! True if field is from North to South (Input) |
|---|
| 16 | LOGICAL ldper ! True if take into account the periodicity (Input) |
|---|
| 17 | REAL pmask ! mask value (Input) |
|---|
| 18 | REAL pfild(kxlon,kylat) ! field to be extrapolated (Input/Output) |
|---|
| 19 | REAL pwork(kxlon,kylat) ! working space |
|---|
| 20 | c |
|---|
| 21 | REAL zwmsk |
|---|
| 22 | INTEGER incre, idoit, i, j, k, inbor, ideb, ifin, ilon, jlat |
|---|
| 23 | INTEGER ix(9), jy(9) ! index arrays for the neighbors coordinates |
|---|
| 24 | REAL zmask(9) |
|---|
| 25 | C |
|---|
| 26 | C We search over the eight closest neighbors |
|---|
| 27 | C |
|---|
| 28 | C j+1 7 8 9 |
|---|
| 29 | C j 4 5 6 Current point 5 --> (i,j) |
|---|
| 30 | C j-1 1 2 3 |
|---|
| 31 | C i-1 i i+1 |
|---|
| 32 | c |
|---|
| 33 | c |
|---|
| 34 | IF (norsud) THEN |
|---|
| 35 | DO j = 1, kylat |
|---|
| 36 | DO i = 1, kxlon |
|---|
| 37 | pwork(i,j) = pfild(i,kylat-j+1) |
|---|
| 38 | ENDDO |
|---|
| 39 | ENDDO |
|---|
| 40 | DO j = 1, kylat |
|---|
| 41 | DO i = 1, kxlon |
|---|
| 42 | pfild(i,j) = pwork(i,j) |
|---|
| 43 | ENDDO |
|---|
| 44 | ENDDO |
|---|
| 45 | ENDIF |
|---|
| 46 | c |
|---|
| 47 | incre = 0 |
|---|
| 48 | c |
|---|
| 49 | DO j = 1, kylat |
|---|
| 50 | DO i = 1, kxlon |
|---|
| 51 | pwork(i,j) = pfild(i,j) |
|---|
| 52 | ENDDO |
|---|
| 53 | ENDDO |
|---|
| 54 | c |
|---|
| 55 | C* To avoid problems in floating point tests |
|---|
| 56 | zwmsk = pmask - 1.0 |
|---|
| 57 | c |
|---|
| 58 | 200 CONTINUE |
|---|
| 59 | incre = incre + 1 |
|---|
| 60 | DO 99999 j = 1, kylat |
|---|
| 61 | DO 99999 i = 1, kxlon |
|---|
| 62 | IF (pfild(i,j).GT. zwmsk) THEN |
|---|
| 63 | pwork(i,j) = pfild(i,j) |
|---|
| 64 | inbor = 0 |
|---|
| 65 | ideb = 1 |
|---|
| 66 | ifin = 9 |
|---|
| 67 | C |
|---|
| 68 | C* Fill up ix array |
|---|
| 69 | ix(1) = MAX (1,i-1) |
|---|
| 70 | ix(2) = i |
|---|
| 71 | ix(3) = MIN (kxlon,i+1) |
|---|
| 72 | ix(4) = MAX (1,i-1) |
|---|
| 73 | ix(5) = i |
|---|
| 74 | ix(6) = MIN (kxlon,i+1) |
|---|
| 75 | ix(7) = MAX (1,i-1) |
|---|
| 76 | ix(8) = i |
|---|
| 77 | ix(9) = MIN (kxlon,i+1) |
|---|
| 78 | C |
|---|
| 79 | C* Fill up iy array |
|---|
| 80 | jy(1) = MAX (1,j-1) |
|---|
| 81 | jy(2) = MAX (1,j-1) |
|---|
| 82 | jy(3) = MAX (1,j-1) |
|---|
| 83 | jy(4) = j |
|---|
| 84 | jy(5) = j |
|---|
| 85 | jy(6) = j |
|---|
| 86 | jy(7) = MIN (kylat,j+1) |
|---|
| 87 | jy(8) = MIN (kylat,j+1) |
|---|
| 88 | jy(9) = MIN (kylat,j+1) |
|---|
| 89 | C |
|---|
| 90 | C* Correct latitude bounds if southernmost or northernmost points |
|---|
| 91 | IF (j .EQ. 1) ideb = 4 |
|---|
| 92 | IF (j .EQ. kylat) ifin = 6 |
|---|
| 93 | C |
|---|
| 94 | C* Account for periodicity in longitude |
|---|
| 95 | C |
|---|
| 96 | IF (ldper) THEN |
|---|
| 97 | IF (i .EQ. kxlon) THEN |
|---|
| 98 | ix(3) = 1 |
|---|
| 99 | ix(6) = 1 |
|---|
| 100 | ix(9) = 1 |
|---|
| 101 | ELSE IF (i .EQ. 1) THEN |
|---|
| 102 | ix(1) = kxlon |
|---|
| 103 | ix(4) = kxlon |
|---|
| 104 | ix(7) = kxlon |
|---|
| 105 | ENDIF |
|---|
| 106 | ELSE |
|---|
| 107 | IF (i .EQ. 1) THEN |
|---|
| 108 | ix(1) = i |
|---|
| 109 | ix(2) = i + 1 |
|---|
| 110 | ix(3) = i |
|---|
| 111 | ix(4) = i + 1 |
|---|
| 112 | ix(5) = i |
|---|
| 113 | ix(6) = i + 1 |
|---|
| 114 | ENDIF |
|---|
| 115 | IF (i .EQ. kxlon) THEN |
|---|
| 116 | ix(1) = i -1 |
|---|
| 117 | ix(2) = i |
|---|
| 118 | ix(3) = i - 1 |
|---|
| 119 | ix(4) = i |
|---|
| 120 | ix(5) = i - 1 |
|---|
| 121 | ix(6) = i |
|---|
| 122 | ENDIF |
|---|
| 123 | C |
|---|
| 124 | IF (i .EQ. 1 .OR. i .EQ. kxlon) THEN |
|---|
| 125 | jy(1) = MAX (1,j-1) |
|---|
| 126 | jy(2) = MAX (1,j-1) |
|---|
| 127 | jy(3) = j |
|---|
| 128 | jy(4) = j |
|---|
| 129 | jy(5) = MIN (kylat,j+1) |
|---|
| 130 | jy(6) = MIN (kylat,j+1) |
|---|
| 131 | C |
|---|
| 132 | ideb = 1 |
|---|
| 133 | ifin = 6 |
|---|
| 134 | IF (j .EQ. 1) ideb = 3 |
|---|
| 135 | IF (j .EQ. kylat) ifin = 4 |
|---|
| 136 | ENDIF |
|---|
| 137 | ENDIF ! end for ldper test |
|---|
| 138 | C |
|---|
| 139 | C* Find unmasked neighbors |
|---|
| 140 | C |
|---|
| 141 | DO 230 k = ideb, ifin |
|---|
| 142 | zmask(k) = 0. |
|---|
| 143 | ilon = ix(k) |
|---|
| 144 | jlat = jy(k) |
|---|
| 145 | IF (pfild(ilon,jlat) .LT. zwmsk) THEN |
|---|
| 146 | zmask(k) = 1. |
|---|
| 147 | inbor = inbor + 1 |
|---|
| 148 | ENDIF |
|---|
| 149 | 230 CONTINUE |
|---|
| 150 | C |
|---|
| 151 | C* Not enough points around point P are unmasked; interpolation on P |
|---|
| 152 | C will be done in a future call to extrap. |
|---|
| 153 | C |
|---|
| 154 | IF (inbor .GE. knbor) THEN |
|---|
| 155 | pwork(i,j) = 0. |
|---|
| 156 | DO k = ideb, ifin |
|---|
| 157 | ilon = ix(k) |
|---|
| 158 | jlat = jy(k) |
|---|
| 159 | pwork(i,j) = pwork(i,j) |
|---|
| 160 | $ + pfild(ilon,jlat) * zmask(k)/ REAL(inbor) |
|---|
| 161 | ENDDO |
|---|
| 162 | ENDIF |
|---|
| 163 | C |
|---|
| 164 | ENDIF |
|---|
| 165 | 99999 CONTINUE |
|---|
| 166 | C |
|---|
| 167 | C* 3. Writing back unmasked field in pfild |
|---|
| 168 | C ------------------------------------ |
|---|
| 169 | C |
|---|
| 170 | C* pfild then contains: |
|---|
| 171 | C - Values which were not masked |
|---|
| 172 | C - Interpolated values from the inbor neighbors |
|---|
| 173 | C - Values which are not yet interpolated |
|---|
| 174 | C |
|---|
| 175 | idoit = 0 |
|---|
| 176 | DO j = 1, kylat |
|---|
| 177 | DO i = 1, kxlon |
|---|
| 178 | IF (pwork(i,j) .GT. zwmsk) idoit = idoit + 1 |
|---|
| 179 | pfild(i,j) = pwork(i,j) |
|---|
| 180 | ENDDO |
|---|
| 181 | ENDDO |
|---|
| 182 | c |
|---|
| 183 | IF (idoit .ne. 0) GOTO 200 |
|---|
| 184 | ccc PRINT*, "Number of extrapolation steps incre =", incre |
|---|
| 185 | c |
|---|
| 186 | IF (norsud) THEN |
|---|
| 187 | DO j = 1, kylat |
|---|
| 188 | DO i = 1, kxlon |
|---|
| 189 | pwork(i,j) = pfild(i,kylat-j+1) |
|---|
| 190 | ENDDO |
|---|
| 191 | ENDDO |
|---|
| 192 | DO j = 1, kylat |
|---|
| 193 | DO i = 1, kxlon |
|---|
| 194 | pfild(i,j) = pwork(i,j) |
|---|
| 195 | ENDDO |
|---|
| 196 | ENDDO |
|---|
| 197 | ENDIF |
|---|
| 198 | c |
|---|
| 199 | RETURN |
|---|
| 200 | END |
|---|
