source: LMDZ.3.3/branches/rel-LF/libf/dyn3d/extrapol.F @ 319

Last change on this file since 319 was 318, checked in by lmdzadmin, 23 years ago

Import dans la version couplée
LF

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