trunk/dyn3d/extrapol.f

!
! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/extrapol.F,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $
!
C
C
      SUBROUTINE extrapol (pfild, kxlon, kylat, pmask,
     .                   norsud, ldper, knbor, pwork)
      IMPLICIT none
c
c OASIS routine (Adaptation: Laurent Li, le 14 mars 1997)
c Fill up missed values by using the neighbor points
c
      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
c
      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)
C
C  We search over the eight closest neighbors
C
C            j+1  7  8  9
C              j  4  5  6    Current point 5 --> (i,j)
C            j-1  1  2  3
C                i-1 i i+1
c
c
      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
c
      incre = 0
c
      DO j = 1, kylat
      DO i = 1, kxlon
         pwork(i,j) = pfild(i,j)
      ENDDO
      ENDDO
c
C* To avoid problems in floating point tests
      zwmsk = pmask - 1.0
c
200   CONTINUE
      incre = incre + 1
      DO 99999 j = 1, kylat
      DO 99999 i = 1, kxlon
      IF (pfild(i,j).GT. zwmsk) THEN
         pwork(i,j) = pfild(i,j)
         inbor = 0
         ideb = 1
         ifin = 9
C
C* 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)
C
C* 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)
C
C* Correct latitude bounds if southernmost or northernmost points
         IF (j .EQ. 1) ideb = 4
         IF (j .EQ. kylat) ifin = 6
C
C* Account for periodicity in longitude
C
         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
C
            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)
C
               ideb = 1
               ifin = 6
               IF (j .EQ. 1) ideb = 3
               IF (j .EQ. kylat) ifin = 4
            ENDIF
         ENDIF ! end for ldper test
C
C* Find unmasked neighbors
C
         DO 230 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
 230     CONTINUE
C
C* Not enough points around point P are unmasked; interpolation on P 
C  will be done in a future call to extrap.
C
         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)/FLOAT(inbor)
            ENDDO
         ENDIF
C
      ENDIF
99999 CONTINUE
C
C*    3. Writing back unmasked field in pfild
C        ------------------------------------
C
C* pfild then contains:
C     - Values which were not masked
C     - Interpolated values from the inbor neighbors
C     - Values which are not yet interpolated
C
      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
c
      IF (idoit .ne. 0) GOTO 200
ccc      PRINT*, "Number of extrapolation steps incre =", incre
c
      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
c
      RETURN
      END
1	guez	3	!
2			! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/extrapol.F,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $
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)/FLOAT(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