/[lmdze]/trunk/dyn3d/extrapol.f
ViewVC logotype

Diff of /trunk/dyn3d/extrapol.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/dyn3d/extrapol.f revision 76 by guez, Fri Nov 15 18:45:49 2013 UTC trunk/dyn3d/extrapol.f90 revision 81 by guez, Wed Mar 5 14:38:41 2014 UTC
# Line 1  Line 1 
1  !  
2  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/extrapol.F,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/extrapol.F,v 1.1.1.1 2004/05/19
3  !  ! 12:53:07 lmdzadmin Exp $
4  C  
5  C  
6        SUBROUTINE extrapol (pfild, kxlon, kylat, pmask,  
7       .                   norsud, ldper, knbor, pwork)  SUBROUTINE extrapol(pfild, kxlon, kylat, pmask, norsud, ldper, knbor, pwork)
8        IMPLICIT none    IMPLICIT NONE
9  c  
10  c OASIS routine (Adaptation: Laurent Li, le 14 mars 1997)    ! OASIS routine (Adaptation: Laurent Li, le 14 mars 1997)
11  c Fill up missed values by using the neighbor points    ! Fill up missed values by using the neighbor points
12  c  
13        INTEGER kxlon, kylat ! longitude and latitude dimensions (Input)    INTEGER kxlon, kylat ! longitude and latitude dimensions (Input)
14        INTEGER knbor ! minimum neighbor number (Input)    INTEGER knbor ! minimum neighbor number (Input)
15        LOGICAL norsud ! True if field is from North to South (Input)    LOGICAL norsud ! True if field is from North to South (Input)
16        LOGICAL ldper ! True if take into account the periodicity (Input)    LOGICAL ldper ! True if take into account the periodicity (Input)
17        REAL pmask ! mask value (Input)    REAL pmask ! mask value (Input)
18        REAL pfild(kxlon,kylat) ! field to be extrapolated (Input/Output)    REAL pfild(kxlon, kylat) ! field to be extrapolated (Input/Output)
19        REAL pwork(kxlon,kylat) ! working space    REAL pwork(kxlon, kylat) ! working space
20  c  
21        REAL zwmsk    REAL zwmsk
22        INTEGER incre, idoit, i, j, k, inbor, ideb, ifin, ilon, jlat    INTEGER incre, idoit, i, j, k, inbor, ideb, ifin, ilon, jlat
23        INTEGER ix(9), jy(9) ! index arrays for the neighbors coordinates    INTEGER ix(9), jy(9) ! index arrays for the neighbors coordinates
24        REAL zmask(9)    REAL zmask(9)
25  C  
26  C  We search over the eight closest neighbors    ! We search over the eight closest neighbors
27  C  
28  C            j+1  7  8  9    ! j+1  7  8  9
29  C              j  4  5  6    Current point 5 --> (i,j)    ! j  4  5  6    Current point 5 --> (i,j)
30  C            j-1  1  2  3    ! j-1  1  2  3
31  C                i-1 i i+1    ! i-1 i i+1
32  c  
33  c  
34        IF (norsud) THEN    IF (norsud) THEN
35           DO j = 1, kylat      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  
36        DO i = 1, kxlon        DO i = 1, kxlon
37           pwork(i,j) = pfild(i,j)          pwork(i, j) = pfild(i, kylat-j+1)
38        ENDDO        END DO
39        ENDDO      END DO
40  c      DO j = 1, kylat
41  C* To avoid problems in floating point tests        DO i = 1, kxlon
42        zwmsk = pmask - 1.0          pfild(i, j) = pwork(i, j)
43  c        END DO
44  200   CONTINUE      END DO
45        incre = incre + 1    END IF
46        DO 99999 j = 1, kylat  
47        DO 99999 i = 1, kxlon    incre = 0
48        IF (pfild(i,j).GT. zwmsk) THEN  
49           pwork(i,j) = pfild(i,j)    DO j = 1, kylat
50           inbor = 0      DO i = 1, kxlon
51           ideb = 1        pwork(i, j) = pfild(i, j)
52           ifin = 9      END DO
53  C    END DO
54  C* Fill up ix array  
55           ix(1) = MAX (1,i-1)    ! * To avoid problems in floating point tests
56           ix(2) = i    zwmsk = pmask - 1.0
57           ix(3) = MIN (kxlon,i+1)  
58           ix(4) = MAX (1,i-1)  200 CONTINUE
59           ix(5) = i    incre = incre + 1
60           ix(6) = MIN (kxlon,i+1)    DO j = 1, kylat
61           ix(7) = MAX (1,i-1)      DO i = 1, kxlon
62           ix(8) = i        IF (pfild(i,j)>zwmsk) THEN
63           ix(9) = MIN (kxlon,i+1)          pwork(i, j) = pfild(i, j)
64  C          inbor = 0
65  C* Fill up iy array          ideb = 1
66           jy(1) = MAX (1,j-1)          ifin = 9
67           jy(2) = MAX (1,j-1)  
68           jy(3) = MAX (1,j-1)          ! * Fill up ix array
69           jy(4) = j          ix(1) = max(1, i-1)
70           jy(5) = j          ix(2) = i
71           jy(6) = j          ix(3) = min(kxlon, i+1)
72           jy(7) = MIN (kylat,j+1)          ix(4) = max(1, i-1)
73           jy(8) = MIN (kylat,j+1)          ix(5) = i
74           jy(9) = MIN (kylat,j+1)          ix(6) = min(kxlon, i+1)
75  C          ix(7) = max(1, i-1)
76  C* Correct latitude bounds if southernmost or northernmost points          ix(8) = i
77           IF (j .EQ. 1) ideb = 4          ix(9) = min(kxlon, i+1)
78           IF (j .EQ. kylat) ifin = 6  
79  C          ! * Fill up iy array
80  C* Account for periodicity in longitude          jy(1) = max(1, j-1)
81  C          jy(2) = max(1, j-1)
82           IF (ldper) THEN          jy(3) = max(1, j-1)
83              IF (i .EQ. kxlon) THEN          jy(4) = j
84                 ix(3) = 1          jy(5) = j
85                 ix(6) = 1          jy(6) = j
86                 ix(9) = 1          jy(7) = min(kylat, j+1)
87              ELSE IF (i .EQ. 1) THEN          jy(8) = min(kylat, j+1)
88                 ix(1) = kxlon          jy(9) = min(kylat, j+1)
89                 ix(4) = kxlon  
90                 ix(7) = kxlon          ! * Correct latitude bounds if southernmost or northernmost points
91              ENDIF          IF (j==1) ideb = 4
92           ELSE          IF (j==kylat) ifin = 6
93              IF (i .EQ. 1) THEN  
94                 ix(1) = i          ! * Account for periodicity in longitude
95                 ix(2) = i + 1  
96                 ix(3) = i          IF (ldper) THEN
97                 ix(4) = i + 1            IF (i==kxlon) THEN
98                 ix(5) = i              ix(3) = 1
99                 ix(6) = i + 1              ix(6) = 1
100              ENDIF              ix(9) = 1
101              IF (i .EQ. kxlon) THEN            ELSE IF (i==1) THEN
102                 ix(1) = i -1              ix(1) = kxlon
103                 ix(2) = i              ix(4) = kxlon
104                 ix(3) = i - 1              ix(7) = kxlon
105                 ix(4) = i            END IF
106                 ix(5) = i - 1          ELSE
107                 ix(6) = i            IF (i==1) THEN
108              ENDIF              ix(1) = i
109  C              ix(2) = i + 1
110              IF (i .EQ. 1 .OR. i .EQ. kxlon) THEN              ix(3) = i
111                 jy(1) = MAX (1,j-1)              ix(4) = i + 1
112                 jy(2) = MAX (1,j-1)              ix(5) = i
113                 jy(3) = j              ix(6) = i + 1
114                 jy(4) = j            END IF
115                 jy(5) = MIN (kylat,j+1)            IF (i==kxlon) THEN
116                 jy(6) = MIN (kylat,j+1)              ix(1) = i - 1
117  C              ix(2) = i
118                 ideb = 1              ix(3) = i - 1
119                 ifin = 6              ix(4) = i
120                 IF (j .EQ. 1) ideb = 3              ix(5) = i - 1
121                 IF (j .EQ. kylat) ifin = 4              ix(6) = i
122              ENDIF            END IF
123           ENDIF ! end for ldper test  
124  C            IF (i==1 .OR. i==kxlon) THEN
125  C* Find unmasked neighbors              jy(1) = max(1, j-1)
126  C              jy(2) = max(1, j-1)
127           DO 230 k = ideb, ifin              jy(3) = j
128              zmask(k) = 0.              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)              ilon = ix(k)
158              jlat = jy(k)              jlat = jy(k)
159              IF (pfild(ilon,jlat) .LT. zwmsk) THEN              pwork(i, j) = pwork(i, j) + pfild(ilon, jlat)*zmask(k)/float( &
160                 zmask(k) = 1.                inbor)
161                 inbor = inbor + 1            END DO
162              ENDIF          END IF
163   230     CONTINUE  
164  C        END IF
165  C* Not enough points around point P are unmasked; interpolation on P      END DO
166  C  will be done in a future call to extrap.    END DO
167  C  
168           IF (inbor .GE. knbor) THEN    ! *    3. Writing back unmasked field in pfild
169              pwork(i,j) = 0.    ! ------------------------------------
170              DO k = ideb, ifin  
171                 ilon = ix(k)    ! * pfild then contains:
172                 jlat = jy(k)    ! - Values which were not masked
173                 pwork(i,j) = pwork(i,j)    ! - Interpolated values from the inbor neighbors
174       $                      + pfild(ilon,jlat) * zmask(k)/FLOAT(inbor)    ! - Values which are not yet interpolated
175              ENDDO  
176           ENDIF    idoit = 0
177  C    DO j = 1, kylat
178        ENDIF      DO i = 1, kxlon
179  99999 CONTINUE        IF (pwork(i,j)>zwmsk) idoit = idoit + 1
180  C        pfild(i, j) = pwork(i, j)
181  C*    3. Writing back unmasked field in pfild      END DO
182  C        ------------------------------------    END DO
183  C  
184  C* pfild then contains:    IF (idoit/=0) GO TO 200
185  C     - Values which were not masked    ! cc      PRINT*, "Number of extrapolation steps incre =", incre
186  C     - Interpolated values from the inbor neighbors  
187  C     - Values which are not yet interpolated    IF (norsud) THEN
188  C      DO j = 1, kylat
189        idoit = 0        DO i = 1, kxlon
190        DO j = 1, kylat          pwork(i, j) = pfild(i, kylat-j+1)
191          END DO
192        END DO
193        DO j = 1, kylat
194        DO i = 1, kxlon        DO i = 1, kxlon
195           IF (pwork(i,j) .GT. zwmsk) idoit = idoit + 1          pfild(i, j) = pwork(i, j)
196           pfild(i,j) = pwork(i,j)        END DO
197        ENDDO      END DO
198        ENDDO    END IF
199  c  
200        IF (idoit .ne. 0) GOTO 200    RETURN
201  ccc      PRINT*, "Number of extrapolation steps incre =", incre  END SUBROUTINE extrapol
 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  
Legend:
Removed from v.76  
changed lines
  Added in v.81
  ViewVC Help
Powered by ViewVC 1.1.21