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
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	DO i = 1, kxlon
37	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	ilon = ix(k)
158	jlat = jy(k)
159	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	DO i = 1, kxlon
190	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