source: LMDZ5/trunk/libf/dyn3dpar/extrapol.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.0 KB
Line 
1!
2! $Id: extrapol.F 1907 2013-11-26 13:10:46Z lguez $
3!
4C
5C
6      SUBROUTINE extrapol (pfild, kxlon, kylat, pmask,
7     .                   norsud, ldper, knbor, pwork)
8      IMPLICIT none
9c
10c OASIS routine (Adaptation: Laurent Li, le 14 mars 1997)
11c Fill up missed values by using the neighbor points
12c
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
20c
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)
25C
26C  We search over the eight closest neighbors
27C
28C            j+1  7  8  9
29C              j  4  5  6    Current point 5 --> (i,j)
30C            j-1  1  2  3
31C                i-1 i i+1
32c
33c
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
46c
47      incre = 0
48c
49      DO j = 1, kylat
50      DO i = 1, kxlon
51         pwork(i,j) = pfild(i,j)
52      ENDDO
53      ENDDO
54c
55C* To avoid problems in floating point tests
56      zwmsk = pmask - 1.0
57c
58200   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
67C
68C* 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)
78C
79C* 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)
89C
90C* Correct latitude bounds if southernmost or northernmost points
91         IF (j .EQ. 1) ideb = 4
92         IF (j .EQ. kylat) ifin = 6
93C
94C* Account for periodicity in longitude
95C
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
123C
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)
131C
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
138C
139C* Find unmasked neighbors
140C
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
150C
151C* Not enough points around point P are unmasked; interpolation on P
152C  will be done in a future call to extrap.
153C
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)/ REAL(inbor)
161            ENDDO
162         ENDIF
163C
164      ENDIF
16599999 CONTINUE
166C
167C*    3. Writing back unmasked field in pfild
168C        ------------------------------------
169C
170C* pfild then contains:
171C     - Values which were not masked
172C     - Interpolated values from the inbor neighbors
173C     - Values which are not yet interpolated
174C
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
182c
183      IF (idoit .ne. 0) GOTO 200
184ccc      PRINT*, "Number of extrapolation steps incre =", incre
185c
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
198c
199      RETURN
200      END
Note: See TracBrowser for help on using the repository browser.