source: LMDZ.3.3/trunk/libf/dyn3d/extrapol.F @ 346

Last change on this file since 346 was 259, checked in by lmdz, 23 years ago

Nouveaux programmes pour la creation des etats initiaux et des conditions aux limites. Levan
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.