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

Annotation of /trunk/dyn3d/extrapol.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 81 - (hide annotations)
Wed Mar 5 14:38:41 2014 UTC (10 years, 3 months ago) by guez
Original Path: trunk/dyn3d/extrapol.f90
File size: 4802 byte(s) 
 Converted to free source form files which were still in fixed source
form. The conversion was done using the polish mode of the NAG Fortran
Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).
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
  ViewVC Help
Powered by ViewVC 1.1.21