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 $


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)
      END DO
    END DO
    DO j = 1, kylat
      DO i = 1, kxlon
        pfild(i, j) = pwork(i, j)
      END DO
    END DO
  END IF

  incre = 0

  DO j = 1, kylat
    DO i = 1, kxlon
      pwork(i, j) = pfild(i, j)
    END DO
  END DO

  ! * 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)>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==1) ideb = 4
        IF (j==kylat) ifin = 6

        ! * Account for periodicity in longitude

        IF (ldper) THEN
          IF (i==kxlon) THEN
            ix(3) = 1
            ix(6) = 1
            ix(9) = 1
          ELSE IF (i==1) THEN
            ix(1) = kxlon
            ix(4) = kxlon
            ix(7) = kxlon
          END IF
        ELSE
          IF (i==1) THEN
            ix(1) = i
            ix(2) = i + 1
            ix(3) = i
            ix(4) = i + 1
            ix(5) = i
            ix(6) = i + 1
          END IF
          IF (i==kxlon) THEN
            ix(1) = i - 1
            ix(2) = i
            ix(3) = i - 1
            ix(4) = i
            ix(5) = i - 1
            ix(6) = i
          END IF

          IF (i==1 .OR. i==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==1) ideb = 3
            IF (j==kylat) ifin = 4
          END IF
        END IF ! end for ldper test

        ! * Find unmasked neighbors

        DO k = ideb, ifin
          zmask(k) = 0.
          ilon = ix(k)
          jlat = jy(k)
          IF (pfild(ilon,jlat)<zwmsk) THEN
            zmask(k) = 1.
            inbor = inbor + 1
          END IF
        END DO

        ! * Not enough points around point P are unmasked; interpolation on P
        ! will be done in a future call to extrap.

        IF (inbor>=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)
          END DO
        END IF

      END IF
    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)>zwmsk) idoit = idoit + 1
      pfild(i, j) = pwork(i, j)
    END DO
  END DO

  IF (idoit/=0) GO TO 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)
      END DO
    END DO
    DO j = 1, kylat
      DO i = 1, kxlon
        pfild(i, j) = pwork(i, j)
      END DO
    END DO
  END IF

  RETURN
END SUBROUTINE extrapol
1	guez	81
2			! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/extrapol.F,v 1.1.1.1 2004/05/19
3			! 12:53:07 lmdzadmin Exp $
4
5
6
7			SUBROUTINE extrapol(pfild, kxlon, kylat, pmask, norsud, ldper, knbor, pwork)
8			IMPLICIT NONE
9
10			! OASIS routine (Adaptation: Laurent Li, le 14 mars 1997)
11			! Fill up missed values by using the neighbor points
12
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
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
26			! We search over the eight closest neighbors
27
28			! j+1 7 8 9
29			! j 4 5 6 Current point 5 --> (i,j)
30			! j-1 1 2 3
31			! i-1 i i+1
32
33
34			IF (norsud) THEN
35			DO j = 1, kylat
36	guez	3	DO i = 1, kxlon
37	guez	81	pwork(i, j) = pfild(i, kylat-j+1)
38			END DO
39			END DO
40			DO j = 1, kylat
41			DO i = 1, kxlon
42			pfild(i, j) = pwork(i, j)
43			END DO
44			END DO
45			END IF
46
47			incre = 0
48
49			DO j = 1, kylat
50			DO i = 1, kxlon
51			pwork(i, j) = pfild(i, j)
52			END DO
53			END DO
54
55			! * To avoid problems in floating point tests
56			zwmsk = pmask - 1.0
57
58			200 CONTINUE
59			incre = incre + 1
60			DO j = 1, kylat
61			DO i = 1, kxlon
62			IF (pfild(i,j)>zwmsk) THEN
63			pwork(i, j) = pfild(i, j)
64			inbor = 0
65			ideb = 1
66			ifin = 9
67
68			! * 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
79			! * 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
90			! * Correct latitude bounds if southernmost or northernmost points
91			IF (j==1) ideb = 4
92			IF (j==kylat) ifin = 6
93
94			! * Account for periodicity in longitude
95
96			IF (ldper) THEN
97			IF (i==kxlon) THEN
98			ix(3) = 1
99			ix(6) = 1
100			ix(9) = 1
101			ELSE IF (i==1) THEN
102			ix(1) = kxlon
103			ix(4) = kxlon
104			ix(7) = kxlon
105			END IF
106			ELSE
107			IF (i==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			END IF
115			IF (i==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			END IF
123
124			IF (i==1 .OR. i==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
132			ideb = 1
133			ifin = 6
134			IF (j==1) ideb = 3
135			IF (j==kylat) ifin = 4
136			END IF
137			END IF ! end for ldper test
138
139			! * Find unmasked neighbors
140
141			DO k = ideb, ifin
142			zmask(k) = 0.
143			ilon = ix(k)
144			jlat = jy(k)
145			IF (pfild(ilon,jlat)<zwmsk) THEN
146			zmask(k) = 1.
147			inbor = inbor + 1
148			END IF
149			END DO
150
151			! * Not enough points around point P are unmasked; interpolation on P
152			! will be done in a future call to extrap.
153
154			IF (inbor>=knbor) THEN
155			pwork(i, j) = 0.
156			DO k = ideb, ifin
157	guez	3	ilon = ix(k)
158			jlat = jy(k)
159	guez	81	pwork(i, j) = pwork(i, j) + pfild(ilon, jlat)*zmask(k)/float( &
160			inbor)
161			END DO
162			END IF
163
164			END IF
165			END DO
166			END DO
167
168			! * 3. Writing back unmasked field in pfild
169			! ------------------------------------
170
171			! * pfild then contains:
172			! - Values which were not masked
173			! - Interpolated values from the inbor neighbors
174			! - Values which are not yet interpolated
175
176			idoit = 0
177			DO j = 1, kylat
178			DO i = 1, kxlon
179			IF (pwork(i,j)>zwmsk) idoit = idoit + 1
180			pfild(i, j) = pwork(i, j)
181			END DO
182			END DO
183
184			IF (idoit/=0) GO TO 200
185			! cc PRINT*, "Number of extrapolation steps incre =", incre
186
187			IF (norsud) THEN
188			DO j = 1, kylat
189	guez	3	DO i = 1, kxlon
190	guez	81	pwork(i, j) = pfild(i, kylat-j+1)
191			END DO
192			END DO
193			DO j = 1, kylat
194			DO i = 1, kxlon
195			pfild(i, j) = pwork(i, j)
196			END DO
197			END DO
198			END IF
199
200			RETURN
201			END SUBROUTINE extrapol