source: lmdz_wrf/trunk/tools/interpolate.f90 @ 2024

Last change on this file since 2024 was 1654, checked in by lfita, 7 years ago

Fixing:

  • calcualtion of mean weighted-center of the polygon at `polygon_properties'
  • Adding compilation at single resolution
File size: 10.8 KB
RevLine 
[1140]1PROGRAM 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
12END PROGRAM interpolate
13
14SUBROUTINE 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
[1608]18  USE module_definitions
[1140]19  USE module_generic
20
21  IMPLICIT NONE
22
[1184]23!  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
[1140]24  INTEGER, INTENT(in)                                    :: dx, dy
25  REAL(r_k), DIMENSION(dx,dy), INTENT(in)                :: ilon, ilat
26  REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in)          :: fraclon, fraclat
27  REAL(r_k), INTENT(in)                                  :: lonv, latv, per
28  REAL(r_k), DIMENSION(2), INTENT(in)                    :: nxlon, nxlat
29  INTEGER, INTENT(in)                                    :: Nperx, Npery
30  INTEGER, DIMENSION(2), INTENT(out)                     :: ilonlat
31  REAL(r_k), INTENT(out)                                 :: mindiffLl
32! Local
33  REAL(r_k), DIMENSION(Nperx,Npery)                      :: difffraclonlat
34  REAL(r_k)                                              :: mindifffracLl
35  INTEGER, DIMENSION(2)                                  :: ilonlatfrac
36  INTEGER                                                :: ixbeg, ixend, iybeg, iyend
37  INTEGER                                                :: fracx, fracy
38  REAL(r_k)                                              :: fraclonv, fraclatv
39  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: difflonlat, lon, lat
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)
[1654]76!  ilonlatfrac = index2DArrayR_K(difffraclonlat, Nperx, Npery, mindifffracLl)
[1140]77  ilonlatfrac = index2DArrayR(difffraclonlat, Nperx, Npery, mindifffracLl)
78
79!  PRINT *, 'mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac
80!  PRINT *, 'frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,',                            &
81!    fraclat(ilonlatfrac(1),ilonlatfrac(2))
82!  PRINT *, 'values lon, lat:', lonv, latv
83     
84! Providing fraction range
85  fraclonv = fraclon(ilonlatfrac(1),ilonlatfrac(2))
86  fraclatv = fraclat(ilonlatfrac(1),ilonlatfrac(2))
87
88  IF (fraclonv >= lonv .AND. fraclatv >= latv) THEN
89    IF (ilonlatfrac(1) > 0) THEN
90      ixbeg = (ilonlatfrac(1)-1)*fracx
91      ixend = ilonlatfrac(1)*fracx+1
92    ELSE
93      ixbeg = 0
94      ixend = fracx+1
95    END IF
96    IF (ilonlatfrac(2) > 0) THEN
97      iybeg = (ilonlatfrac(2)-1)*fracy
98      iyend = ilonlatfrac(2)*fracy+1
99    ELSE
100      iybeg = 0
101      iyend = fracy+1
102    END IF
103  ELSE IF (fraclonv < lonv .AND. fraclatv >= latv) THEN
104    IF (ilonlatfrac(1) < Nperx) THEN
105      IF (ilonlatfrac(1) /= 0) THEN
106        ixbeg = (ilonlatfrac(1)-1)*fracx
107        ixend = ilonlatfrac(1)*fracx+1
108      ELSE
109        ixbeg = 0
110        ixend = fracx+1
111      END IF
112    ELSE
113      ixbeg = Nperx*fracx
114      ixend = dx+1
115    END IF
116    IF (ilonlatfrac(2) > 0) THEN
117      iybeg = (ilonlatfrac(2)-1)*fracy
118      iyend = ilonlatfrac(2)*fracy+1
119    ELSE
120      iybeg = 0
121      iyend = fracy+1
122    END IF   
123  ELSE IF (fraclonv < lonv .AND. fraclatv < latv) THEN
124    IF (ilonlatfrac(1) < Nperx) THEN
125      IF (ilonlatfrac(1) /= 0) THEN
126        ixbeg = (ilonlatfrac(1)-1)*fracx
127        ixend = ilonlatfrac(1)*fracx+1
128      ELSE
129        ixbeg = 0
130        ixend = fracx+1
131      END IF
132    ELSE
133      ixbeg = Nperx*fracx
134      ixend = dx+1
135    ENDIF
136    IF (ilonlatfrac(2) < Npery) THEN
137      IF (ilonlatfrac(2) /= 0) THEN
138        iybeg = (ilonlatfrac(2)-1)*fracy
139        iyend = ilonlatfrac(2)*fracy+1
140      ELSE
141        iybeg = 0
142        iyend = fracy+1
143      END IF
144    ELSE
145      iybeg = Npery*fracy
146      iyend = dy+1
147    END IF
148  ELSE IF (fraclonv >= lonv .AND. fraclatv < latv) THEN
149    IF (ilonlatfrac(1) > 0) THEN
150      ixbeg = (ilonlatfrac(1)-1)*fracx
151      ixend = ilonlatfrac(1)*fracx+1
152    ELSE
153      ixbeg = 0
154      ixend = fracx+1
155    END IF
156    IF (ilonlatfrac(2) < Npery) THEN
157      IF (ilonlatfrac(2) /= 0) THEN
158        iybeg = (ilonlatfrac(2)-1)*fracy
159        iyend = ilonlatfrac(2)*fracy+1
160      ELSE
161        iybeg = 0
162        iyend = fracy+1
163      END IF
164    ELSE
165      iybeg = Npery*fracy
166      iyend = dy+1
167    END IF
168  END IF
169
170  IF (ALLOCATED(lon)) DEALLOCATE(lon)
171  ALLOCATE(lon(ixend-ixbeg+1, iyend-iybeg+1))
172  IF (ALLOCATED(lat)) DEALLOCATE(lat)
173  ALLOCATE(lat(ixend-ixbeg+1, iyend-iybeg+1))
174  IF (ALLOCATED(difflonlat)) DEALLOCATE(difflonlat)
175  ALLOCATE(difflonlat(ixend-ixbeg+1, iyend-iybeg+1))
176
177  lon = ilon(ixbeg:ixend,iybeg:iyend)
178  lat = ilat(ixbeg:ixend,iybeg:iyend)
179
180!  print *,'lon _______'
181!  print *,lon
182!  print *,'lat _______'
183!  print *,lat
184
185! Find point
186  difflonlat = SQRT((lon-lonv)**2. + (lat-latv)**2.)
187  mindiffLl = MINVAL(difflonlat)
188  ilonlat = index2DArrayR(difflonlat, ixend-ixbeg+1, iyend-iybeg+1, mindiffLl)
189
190  ilonlat(1) = ilonlat(1) + ixbeg
191  ilonlat(2) = ilonlat(2) + iybeg
192
193!  PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon
194!  PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2))
195
196  RETURN
197
198END SUBROUTINE CoarselonlatFind
199
200SUBROUTINE CoarseInterpolate(dimx, dimy, projlon, projlat, Ninpts, lonvs, latvs, percen, mindiff,     &
201  ivar, newvar, newvarin, newvarinpt, newvarindiff, ncid)
202! Subroutine which finds the closest grid point within a projection throughout a first guest
203!   approche from percentages of the whole domain
204
205  USE module_generic
206
207  IMPLICIT NONE
208
[1184]209!  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
[1140]210  INTEGER, INTENT(in)                                    :: dimx, dimy
211  REAL(r_k), DIMENSION(dimx,dimy), INTENT(in)            :: projlon, projlat
212  INTEGER, INTENT(in)                                    :: Ninpts
213  REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: ivar, lonvs, latvs
214  REAL(r_k)                                              :: mindiff, percen
215  INTEGER, INTENT(in)                                    :: ncid
216  REAL(r_k), DIMENSION(dimx,dimy), INTENT(inout)         :: newvar
217  INTEGER, DIMENSION(dimx,dimy), INTENT(inout)           :: newvarin
218  INTEGER, DIMENSION(Ninpts), INTENT(inout)              :: newvarinpt
219  REAL(r_k), DIMENSION(Ninpts), INTENT(inout)            :: newvarindiff
220! Local
221  INTEGER                                                :: iv
222  INTEGER, DIMENSION(2)                                  :: ilonlat
223  REAL(r_k)                                              :: mindiffLl
224  INTEGER                                                :: Ninpts1
225  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
226  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: fractionlon, fractionlat
227  INTEGER                                                :: fracdx, fracdy
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
296END SUBROUTINE CoarseInterpolate
Note: See TracBrowser for help on using the repository browser.