Context Navigation

interpolate.F90 @ 1292

Last change on this file since 1292 was 1184, checked in by lfita, 8 years ago

Adding:

Use of generic fillval64' and r_k' definition from 'module_generic.F90'
Working version of the npp' and dis' interpolation in 3D using mask also variable
`add_varattrs': Function to add variables attributes from a given netCDF object to other variables in another one

File size: 10.9 KB

Line
1	PROGRAM interpolate
2	! Program to interpolate values from a giving projection
3	! To be included in a python
4	! f2py -m ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c interpolate.F90 module_generic.F90
5
6	IMPLICIT NONE
7
8	CHARACTER(LEN=50) :: main, ErrWarnMsg
9
10	main='interpolate'
11
12	END PROGRAM interpolate
13
14	SUBROUTINE CoarselonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, fraclon, fraclat, lonv, latv, per, &
15	Nperx, Npery, ilonlat, mindiffLl)
16	! Function to search a given value from a coarser version of the data
17
18	USE module_generic
19
20	IMPLICIT NONE
21
22	! INTEGER, PARAMETER :: r_k = KIND(1.d0)
23	INTEGER, INTENT(in) :: dx, dy
24	REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: ilon, ilat
25	REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in) :: fraclon, fraclat
26	REAL(r_k), INTENT(in) :: lonv, latv, per
27	REAL(r_k), DIMENSION(2), INTENT(in) :: nxlon, nxlat
28	INTEGER, INTENT(in) :: Nperx, Npery
29	INTEGER, DIMENSION(2), INTENT(out) :: ilonlat
30	REAL(r_k), INTENT(out) :: mindiffLl
31	! Local
32	REAL(r_k), DIMENSION(Nperx,Npery) :: difffraclonlat
33	REAL(r_k) :: mindifffracLl
34	INTEGER, DIMENSION(2) :: ilonlatfrac
35	INTEGER :: ixbeg, ixend, iybeg, iyend
36	INTEGER :: fracx, fracy
37	REAL(r_k) :: fraclonv, fraclatv
38	REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: difflonlat, lon, lat
39	CHARACTER(LEN=50) :: fname
40
41	! Variables
42	! ilon, ilat: original 2D matrices with the longitudes and the latitudes
43	! lonv, latv: longitude and latitude to find
44	! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat
45	! per: fraction of the whole domain (as percentage)
46	! Nper[x/y]: period (as fraction over 1) of the fractions of the original grid to use to explore
47	! fraclon, fraclat: longitude and latitude fractional matricies to perform the first guess
48
49	fname = 'CoarselonlatFind'
50
51	IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN
52	PRINT *, TRIM(ErrWarnMsg('err'))
53	PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!'
54	PRINT *,' given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )'
55	STOP
56	END IF
57	IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN
58	PRINT *, TRIM(ErrWarnMsg('err'))
59	PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!'
60	PRINT *,' given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )'
61	STOP
62	END IF
63
64	fracx = int(dx*per)
65	fracy = int(dy*per)
66
67	! PRINT *,'fraclon _______'
68	! PRINT *,fraclon
69
70	! PRINT *,'fraclat _______'
71	! PRINT *,fraclat
72
73	! Fraction point
74	difffraclonlat = SQRT((fraclon-lonv)2. + (fraclat-latv)2.)
75	mindifffracLl = MINVAL(difffraclonlat)
76	ilonlatfrac = index2DArrayR(difffraclonlat, Nperx, Npery, mindifffracLl)
77
78	! PRINT *, 'mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac
79	! PRINT *, 'frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,', &
80	! fraclat(ilonlatfrac(1),ilonlatfrac(2))
81	! PRINT *, 'values lon, lat:', lonv, latv
82
83	! Providing fraction range
84	fraclonv = fraclon(ilonlatfrac(1),ilonlatfrac(2))
85	fraclatv = fraclat(ilonlatfrac(1),ilonlatfrac(2))
86
87	IF (fraclonv >= lonv .AND. fraclatv >= latv) THEN
88	IF (ilonlatfrac(1) > 0) THEN
89	ixbeg = (ilonlatfrac(1)-1)*fracx
90	ixend = ilonlatfrac(1)*fracx+1
91	ELSE
92	ixbeg = 0
93	ixend = fracx+1
94	END IF
95	IF (ilonlatfrac(2) > 0) THEN
96	iybeg = (ilonlatfrac(2)-1)*fracy
97	iyend = ilonlatfrac(2)*fracy+1
98	ELSE
99	iybeg = 0
100	iyend = fracy+1
101	END IF
102	ELSE IF (fraclonv < lonv .AND. fraclatv >= latv) THEN
103	IF (ilonlatfrac(1) < Nperx) THEN
104	IF (ilonlatfrac(1) /= 0) THEN
105	ixbeg = (ilonlatfrac(1)-1)*fracx
106	ixend = ilonlatfrac(1)*fracx+1
107	ELSE
108	ixbeg = 0
109	ixend = fracx+1
110	END IF
111	ELSE
112	ixbeg = Nperx*fracx
113	ixend = dx+1
114	END IF
115	IF (ilonlatfrac(2) > 0) THEN
116	iybeg = (ilonlatfrac(2)-1)*fracy
117	iyend = ilonlatfrac(2)*fracy+1
118	ELSE
119	iybeg = 0
120	iyend = fracy+1
121	END IF
122	ELSE IF (fraclonv < lonv .AND. fraclatv < latv) THEN
123	IF (ilonlatfrac(1) < Nperx) THEN
124	IF (ilonlatfrac(1) /= 0) THEN
125	ixbeg = (ilonlatfrac(1)-1)*fracx
126	ixend = ilonlatfrac(1)*fracx+1
127	ELSE
128	ixbeg = 0
129	ixend = fracx+1
130	END IF
131	ELSE
132	ixbeg = Nperx*fracx
133	ixend = dx+1
134	ENDIF
135	IF (ilonlatfrac(2) < Npery) THEN
136	IF (ilonlatfrac(2) /= 0) THEN
137	iybeg = (ilonlatfrac(2)-1)*fracy
138	iyend = ilonlatfrac(2)*fracy+1
139	ELSE
140	iybeg = 0
141	iyend = fracy+1
142	END IF
143	ELSE
144	iybeg = Npery*fracy
145	iyend = dy+1
146	END IF
147	ELSE IF (fraclonv >= lonv .AND. fraclatv < latv) THEN
148	IF (ilonlatfrac(1) > 0) THEN
149	ixbeg = (ilonlatfrac(1)-1)*fracx
150	ixend = ilonlatfrac(1)*fracx+1
151	ELSE
152	ixbeg = 0
153	ixend = fracx+1
154	END IF
155	IF (ilonlatfrac(2) < Npery) THEN
156	IF (ilonlatfrac(2) /= 0) THEN
157	iybeg = (ilonlatfrac(2)-1)*fracy
158	iyend = ilonlatfrac(2)*fracy+1
159	ELSE
160	iybeg = 0
161	iyend = fracy+1
162	END IF
163	ELSE
164	iybeg = Npery*fracy
165	iyend = dy+1
166	END IF
167	END IF
168
169	IF (ALLOCATED(lon)) DEALLOCATE(lon)
170	ALLOCATE(lon(ixend-ixbeg+1, iyend-iybeg+1))
171	IF (ALLOCATED(lat)) DEALLOCATE(lat)
172	ALLOCATE(lat(ixend-ixbeg+1, iyend-iybeg+1))
173	IF (ALLOCATED(difflonlat)) DEALLOCATE(difflonlat)
174	ALLOCATE(difflonlat(ixend-ixbeg+1, iyend-iybeg+1))
175
176	lon = ilon(ixbeg:ixend,iybeg:iyend)
177	lat = ilat(ixbeg:ixend,iybeg:iyend)
178
179	! print *,'lon _______'
180	! print *,lon
181	! print *,'lat _______'
182	! print *,lat
183
184	! Find point
185	difflonlat = SQRT((lon-lonv)2. + (lat-latv)2.)
186	mindiffLl = MINVAL(difflonlat)
187	ilonlat = index2DArrayR(difflonlat, ixend-ixbeg+1, iyend-iybeg+1, mindiffLl)
188
189	ilonlat(1) = ilonlat(1) + ixbeg
190	ilonlat(2) = ilonlat(2) + iybeg
191
192	! PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon
193	! PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2))
194
195	RETURN
196
197	END SUBROUTINE CoarselonlatFind
198
199	SUBROUTINE CoarseInterpolate(dimx, dimy, projlon, projlat, Ninpts, lonvs, latvs, percen, mindiff, &
200	ivar, newvar, newvarin, newvarinpt, newvarindiff, ncid)
201	! Subroutine which finds the closest grid point within a projection throughout a first guest
202	! approche from percentages of the whole domain
203
204	USE module_generic
205
206	IMPLICIT NONE
207
208	! INTEGER, PARAMETER :: r_k = KIND(1.d0)
209	INTEGER, INTENT(in) :: dimx, dimy
210	REAL(r_k), DIMENSION(dimx,dimy), INTENT(in) :: projlon, projlat
211	INTEGER, INTENT(in) :: Ninpts
212	REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: ivar, lonvs, latvs
213	REAL(r_k) :: mindiff, percen
214	INTEGER, INTENT(in) :: ncid
215	REAL(r_k), DIMENSION(dimx,dimy), INTENT(inout) :: newvar
216	INTEGER, DIMENSION(dimx,dimy), INTENT(inout) :: newvarin
217	INTEGER, DIMENSION(Ninpts), INTENT(inout) :: newvarinpt
218	REAL(r_k), DIMENSION(Ninpts), INTENT(inout) :: newvarindiff
219	! Local
220	INTEGER :: iv
221	INTEGER, DIMENSION(2) :: ilonlat
222	REAL(r_k) :: mindiffLl
223	INTEGER :: Ninpts1
224	REAL(r_k), DIMENSION(2) :: extremelon, extremelat
225	REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: fractionlon, fractionlat
226	INTEGER :: fracdx, fracdy
227	CHARACTER(LEN=50) :: fname, errormsg
228
229	!!!!!!! Variables
230	! dimx, dimy: dimension length of the target interpolation
231	! proj[lon/lat]: longitudes and latitudes of the target interpolation
232	! Ninpts: number of points to interpolate
233	! [lon/lat]vs: longitudes and latitudes of the points to interpolate
234	! mindiff: minimal accepted distance to the target point
235	! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues
236	! ivar: values to localize in the target projection
237	! newvar: localisation of the [lon/lat]vs point in the target projection
238	! newvarin: number of point from the input data
239	! newvarinpt: integer value indicating if the value has been already located (0: no, 1: yes)
240	! newvarindiff: distance of point from the input data to the closest target point
241	! ncid: netCDF output file id
242
243	fname = 'CoarseInterpolate'
244	Ninpts1 = Ninpts/100
245
246	extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
247	extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
248
249	fracdx = INT(dimx*percen)
250	fracdy = INT(dimy*percen)
251
252	IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
253	ALLOCATE(fractionlon(fracdx, fracdy))
254	IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
255	ALLOCATE(fractionlat(fracdx, fracdy))
256
257	fractionlon = projlon(::fracdx,::fracdy)
258	fractionlat = projlat(::fracdx,::fracdy)
259
260	DO iv=1,Ninpts
261	IF (newvarinpt(iv) == 0) THEN
262	CALL CoarselonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, fractionlon, &
263	fractionlat, lonvs(iv), latvs(iv), percen, fracdx, fracdy, ilonlat, mindiffLl)
264
265	IF (mindiffLl <= mindiff) THEN
266	! percendone(iv,Ninpts,0.5,'done:')
267
268	IF (ilonlat(1) >= 0 .AND. ilonlat(1) >= 0) THEN
269	newvar(ilonlat(1),ilonlat(2)) = ivar(iv)
270	newvarin(ilonlat(1),ilonlat(2)) = iv
271	newvarinpt(iv) = 1
272	newvarindiff(iv) = mindiffLl
273	PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv), &
274	' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:', &
275	newvarindiff(iv), ' point:',ilonlat
276	ELSE
277	PRINT *,TRIM(ErrWarnMsg('err'))
278	PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), &
279	' not relocated !!'
280	PRINT *,' mindiffl:', mindiffLl, ' ilon:', ilonlat(1), ' ilat:', ilonlat(2)
281	STOP
282	END IF
283
284	! IF (MOD(iv,Ninpts1) == 0) newnc.sync()
285	ELSE
286	PRINT *,TRIM(ErrWarnMsg('err'))
287	PRINT *,' ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv), &
288	' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ', &
289	mindiff, ' !!'
290	PRINT *,' found minimum difference:', mindiffLl
291	STOP
292	END IF
293	END IF
294	END DO
295
296	END SUBROUTINE CoarseInterpolate

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format