source: lmdz_wrf/trunk/tools/module_ForInterpolate.F90 @ 843

Last change on this file since 843 was 768, checked in by lfita, 9 years ago

Getting there

  • Starting from the last localized point in 'lonlatfrac'
File size: 35.1 KB
Line 
1! Module to interpolate values from a giving projection
2! To be included in a python
3! f2py -m module_ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c module_ForInterpolate.F90 module_generic.F90 >& run_f2py.log
4MODULE module_ForInterpolate
5
6  CONTAINS
7
8SUBROUTINE CoarselonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, fraclon, fraclat, lonv, latv, per,      &
9  Nperx, Npery, ilonlat, mindiffLl)
10! Function to search a given value from a coarser version of the data
11
12  USE module_generic
13
14  IMPLICIT NONE
15
16  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
17  INTEGER, INTENT(in)                                    :: dx, dy
18  REAL(r_k), DIMENSION(dx,dy), INTENT(in)                :: ilon, ilat
19  REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in)          :: fraclon, fraclat
20  REAL(r_k), INTENT(in)                                  :: lonv, latv, per
21  REAL(r_k), DIMENSION(2), INTENT(in)                    :: nxlon, nxlat
22  INTEGER, INTENT(in)                                    :: Nperx, Npery
23  INTEGER, DIMENSION(2), INTENT(out)                     :: ilonlat
24  REAL(r_k), INTENT(out)                                 :: mindiffLl
25! Local
26  REAL(r_k), DIMENSION(Nperx,Npery)                      :: difffraclonlat
27  REAL(r_k)                                              :: mindifffracLl
28  INTEGER, DIMENSION(2)                                  :: ilonlatfrac
29  INTEGER                                                :: ixbeg, ixend, iybeg, iyend
30  INTEGER                                                :: fracx, fracy
31  REAL(r_k)                                              :: fraclonv, fraclatv
32  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: difflonlat, lon, lat
33  CHARACTER(LEN=50)                                      :: fname
34 
35! Variables
36! ilon, ilat: original 2D matrices with the longitudes and the latitudes
37! lonv, latv: longitude and latitude to find
38! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat
39! per: fraction of the whole domain (as percentage)
40! Nper[x/y]: period (as fraction over 1) of the fractions of the original grid to use to explore
41! fraclon, fraclat: longitude and latitude fractional matricies to perform the first guess
42
43  fname = 'CoarselonlatFind'
44
45  IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN
46    PRINT *, TRIM(ErrWarnMsg('err'))
47    PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
48    PRINT *,'    given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )'
49    STOP
50  END IF
51  IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN
52    PRINT *, TRIM(ErrWarnMsg('err'))
53    PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
54    PRINT *,'    given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )'
55    STOP
56  END IF
57
58! Initializing variables
59  ixbeg = 0
60  ixend = 0
61  iybeg = 0
62  iyend = 0
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
197END SUBROUTINE CoarselonlatFind
198
199SUBROUTINE CoarselonlatFindExact(dx, dy, ilon, ilat, nxlon, nxlat, fracx, fracy, fraclon, fraclat,    &
200  iv, lonv, latv, per, Nperx, Npery, mindiff, ilonlat, mindiffLl)
201! Function to search a given value from a coarser version of the data
202
203  USE module_generic
204
205  IMPLICIT NONE
206
207  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
208  INTEGER, INTENT(in)                                    :: dx, dy, iv
209  REAL(r_k), DIMENSION(dx,dy), INTENT(in)                :: ilon, ilat
210  INTEGER, INTENT(in)                                    :: fracx, fracy
211  REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in)          :: fraclon, fraclat
212  REAL(r_k), INTENT(in)                                  :: lonv, latv, per, mindiff
213  REAL(r_k), DIMENSION(2), INTENT(in)                    :: nxlon, nxlat
214  INTEGER, INTENT(in)                                    :: Nperx, Npery
215  INTEGER, DIMENSION(2), INTENT(out)                     :: ilonlat
216  REAL(r_k), INTENT(out)                                 :: mindiffLl
217! Local
218  INTEGER                                                :: i
219  REAL(r_k), DIMENSION(Nperx,Npery)                      :: difffraclonlat
220  REAL(r_k)                                              :: mindifffracLl
221  INTEGER, DIMENSION(2)                                  :: ilonlatfrac
222  INTEGER                                                :: ixbeg, ixend, iybeg, iyend
223  REAL(r_k)                                              :: fraclonv, fraclatv
224  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: difflonlat, lon, lat
225  CHARACTER(LEN=50)                                      :: fname
226 
227! Variables
228! ilon, ilat: original 2D matrices with the longitudes and the latitudes
229! lonv, latv: longitude and latitude to find
230! iv: point in the input data
231! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat
232! per: fraction of the whole domain (as percentage)
233! Nper[x/y]: period (as fraction over 1) of the fractions of the original grid to use to explore
234! frac[x/y]: Number of grid points for each fraction
235! fraclon, fraclat: longitude and latitude fractional matricies to perform the first guess
236! mindiff: authorized minimal distance between input and interpolated point
237! ilonlat: grid point on the total lon,lat matrix
238! mindiffLl: distance between input and interpolated point
239
240  fname = 'CoarselonlatFindExact'
241
242  IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN
243    PRINT *, TRIM(ErrWarnMsg('err'))
244    PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
245    PRINT *,'    given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )'
246    STOP
247  END IF
248  IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN
249    PRINT *, TRIM(ErrWarnMsg('err'))
250    PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
251    PRINT *,'    given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )'
252    STOP
253  END IF
254
255! Initializing variables
256  ixbeg = 0
257  ixend = 0
258  iybeg = 0
259  iyend = 0
260
261! Fraction point
262  difffraclonlat = SQRT((fraclon-lonv)**2. + (fraclat-latv)**2.)
263  mindifffracLl = MINVAL(difffraclonlat)
264  ilonlatfrac = index2DArrayR(difffraclonlat, Nperx, Npery, mindifffracLl)
265
266!  PRINT *, 'mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac
267!  PRINT *, 'frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,',                            &
268!    fraclat(ilonlatfrac(1),ilonlatfrac(2))
269!  PRINT *, 'values lon, lat:', lonv, latv
270     
271! Providing fraction range
272  fraclonv = fraclon(ilonlatfrac(1),ilonlatfrac(2))
273  fraclatv = fraclat(ilonlatfrac(1),ilonlatfrac(2))
274
275  IF (fraclonv >= lonv .AND. fraclatv >= latv) THEN
276    PRINT *,'Lluis!',fraclonv, '>=', lonv,'&', fraclatv, '>=', latv
277    IF (ilonlatfrac(1) > 1) THEN
278      ixbeg = (ilonlatfrac(1)-1)*fracx
279      ixend = ilonlatfrac(1)*fracx+1
280    ELSE
281      PRINT *,'Lluis 2'
282      ixbeg = 1
283      ixend = fracx+1
284    END IF
285    IF (ilonlatfrac(2) > 1) THEN
286      iybeg = (ilonlatfrac(2)-1)*fracy
287      iyend = ilonlatfrac(2)*fracy+1
288    ELSE
289      iybeg = 1
290      iyend = fracy+1
291    END IF
292  ELSE IF (fraclonv < lonv .AND. fraclatv >= latv) THEN
293    PRINT *,'Lluis!',fraclonv, '<', lonv,'&', fraclatv, '>=', latv
294    IF (ilonlatfrac(1) < Nperx) THEN
295      PRINT *,'Lluis 2'
296      IF (ilonlatfrac(1) /= 1) THEN
297        ixbeg = (ilonlatfrac(1)-1)*fracx
298        ixend = ilonlatfrac(1)*fracx+1
299      ELSE
300        ixbeg = 1
301        ixend = fracx+1
302      END IF
303    ELSE
304      ixbeg = Nperx*fracx
305      ixend = dx+1
306    END IF
307    IF (ilonlatfrac(2) > 1) THEN
308      iybeg = (ilonlatfrac(2)-1)*fracy
309      iyend = ilonlatfrac(2)*fracy+1
310    ELSE
311      iybeg = 1
312      iyend = fracy+1
313    END IF   
314  ELSE IF (fraclonv < lonv .AND. fraclatv < latv) THEN
315    PRINT *,'Lluis!',fraclonv, '<', lonv,'&', fraclatv, '<', latv
316    IF (ilonlatfrac(1) < Nperx) THEN
317      IF (ilonlatfrac(1) /= 1) THEN
318        ixbeg = (ilonlatfrac(1)-1)*fracx
319        ixend = ilonlatfrac(1)*fracx+1
320      ELSE
321        ixbeg = 1
322        ixend = fracx+1
323      END IF
324    ELSE
325      ixbeg = Nperx*fracx
326      ixend = dx+1
327    ENDIF
328    IF (ilonlatfrac(2) < Npery) THEN
329      IF (ilonlatfrac(2) /= 1) THEN
330        iybeg = (ilonlatfrac(2)-1)*fracy
331        iyend = ilonlatfrac(2)*fracy+1
332      ELSE
333        iybeg = 1
334        iyend = fracy+1
335      END IF
336    ELSE
337      iybeg = Npery*fracy
338      iyend = dy+1
339    END IF
340  ELSE IF (fraclonv >= lonv .AND. fraclatv < latv) THEN
341    PRINT *,'Llui!',fraclonv, '>=', lonv,'&', fraclatv, '<', latv
342    IF (ilonlatfrac(1) > 1) THEN
343      ixbeg = (ilonlatfrac(1)-1)*fracx
344      ixend = ilonlatfrac(1)*fracx+1
345    ELSE
346      ixbeg = 1
347      ixend = fracx+1
348    END IF
349    IF (ilonlatfrac(2) < Npery) THEN
350      IF (ilonlatfrac(2) /= 1) THEN
351        iybeg = (ilonlatfrac(2)-1)*fracy
352        iyend = ilonlatfrac(2)*fracy+1
353      ELSE
354        iybeg = 1
355        iyend = fracy+1
356      END IF
357    ELSE
358      iybeg = Npery*fracy
359      iyend = dy+1
360    END IF
361  END IF
362
363  IF (ALLOCATED(lon)) DEALLOCATE(lon)
364  ALLOCATE(lon(ixend-ixbeg+1, iyend-iybeg+1))
365  IF (ALLOCATED(lat)) DEALLOCATE(lat)
366  ALLOCATE(lat(ixend-ixbeg+1, iyend-iybeg+1))
367  IF (ALLOCATED(difflonlat)) DEALLOCATE(difflonlat)
368  ALLOCATE(difflonlat(ixend-ixbeg+1, iyend-iybeg+1))
369
370  lon = ilon(ixbeg:ixend,iybeg:iyend)
371  lat = ilat(ixbeg:ixend,iybeg:iyend)
372
373!  print *,'lon _______'
374!  print *,lon
375!  print *,'lat _______'
376!  print *,lat
377
378! Find point
379  difflonlat = SQRT((lon-lonv)**2. + (lat-latv)**2.)
380  mindiffLl = MINVAL(difflonlat)
381
382  IF (mindiffLl > mindiff) THEN
383    difflonlat = SQRT((lon-lonv)**2. + (lat-latv)**2.)
384    mindiffLl = MINVAL(difflonlat)
385  END IF
386
387  IF (mindiffLl > mindiff) THEN
388    PRINT *,TRIM(ErrWarnMsg('err'))
389    PRINT *,'  ' // TRIM(fname) // ': not equivalent point closer than:',mindiff,' found!!'
390    PRINT *,'    at input point iv:', iv,' lon/lat:', lonv,', ',latv,' distance:',mindiffLl
391    PRINT *,'    Fraction values _______ (',Nperx,', ',Npery ,')'
392    PRINT *,'      fraclon'
393    DO i=1, Nperx
394      PRINT *,'      ',fraclon(i,:)
395    END DO
396    PRINT *,'      fraclat'
397    DO i=1, Nperx
398      PRINT *,'      ',fraclat(i,:)
399    END DO
400    PRINT *,'      frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,',                     &
401     fraclat(ilonlatfrac(1),ilonlatfrac(2))
402    PRINT *,'      mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac
403    PRINT *,'    Coarse values _______'
404    PRINT *,'      indices. x:', ixbeg, ', ', ixend, ' y:', iybeg, ', ', iyend
405    PRINT *,'      lon range:', '(',ilon(ixbeg,iybeg),', ',ilon(ixend,iyend),')'
406    PRINT *,'      lat range:', '(',ilat(ixbeg,iybeg),', ',ilat(ixend,iyend),')'
407    PRINT *,'      lon', UBOUND(lon)
408    DO i=1, ixend-ixbeg+1
409      PRINT *,'      ',lon(i,:)
410    END DO
411    PRINT *,'      lat', UBOUND(lat)
412    DO i=1, ixend-ixbeg+1
413      PRINT *,'      ',lat(i,:)
414    END DO
415    STOP
416  END IF
417
418  ilonlat = index2DArrayR(difflonlat, ixend-ixbeg+1, iyend-iybeg+1, mindiffLl)
419
420  ilonlat(1) = ilonlat(1) + ixbeg
421  ilonlat(2) = ilonlat(2) + iybeg
422
423!  PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon
424!  PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2))
425
426  RETURN
427
428END SUBROUTINE CoarselonlatFindExact
429
430SUBROUTINE lonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, lonv, latv, ilonlat, mindiffLl)
431! Function to search a given value from a coarser version of the data
432
433  USE module_generic
434
435  IMPLICIT NONE
436
437  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
438  INTEGER, INTENT(in)                                    :: dx, dy
439  REAL(r_k), DIMENSION(dx,dy), INTENT(in)                :: ilon, ilat
440  REAL(r_k), INTENT(in)                                  :: lonv, latv
441  REAL(r_k), DIMENSION(2), INTENT(in)                    :: nxlon, nxlat
442  INTEGER, DIMENSION(2), INTENT(out)                     :: ilonlat
443  REAL(r_k), INTENT(out)                                 :: mindiffLl
444! Local
445  REAL(r_k), DIMENSION(dx,dy)                            :: difflonlat
446  CHARACTER(LEN=50)                                      :: fname
447 
448! Variables
449! ilon, ilat: original 2D matrices with the longitudes and the latitudes
450! lonv, latv: longitude and latitude to find
451! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat
452
453  fname = 'lonlatFind'
454
455  IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN
456    PRINT *, TRIM(ErrWarnMsg('err'))
457    PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
458    PRINT *,'    given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )'
459    STOP
460  END IF
461  IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN
462    PRINT *, TRIM(ErrWarnMsg('err'))
463    PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
464    PRINT *,'    given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )'
465    STOP
466  END IF
467
468! Find point
469  difflonlat = SQRT((ilon-lonv)**2. + (ilat-latv)**2.)
470  mindiffLl = MINVAL(difflonlat)
471  ilonlat = index2DArrayR(difflonlat, dx, dy, mindiffLl)
472
473!  PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon
474!  PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2))
475
476  RETURN
477
478END SUBROUTINE lonlatFind
479
480SUBROUTINE CoarseInterpolate(projlon, projlat, lonvs, latvs, percen, mindiff, inpt, ilonlat,          &
481  mindiffLl, dimx, dimy, Ninpts)
482! Subroutine which finds the closest grid point within a projection throughout a first guest
483!   approche from percentages of the whole domain
484
485  USE module_generic
486
487  IMPLICIT NONE
488
489  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
490  INTEGER, INTENT(in)                                    :: dimx, dimy
491  REAL(r_k), DIMENSION(dimx,dimy), INTENT(in)            :: projlon, projlat
492  INTEGER, INTENT(in)                                    :: Ninpts
493  REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: inpt, lonvs, latvs
494  REAL(r_k), INTENT(in)                                  :: mindiff, percen
495  INTEGER, DIMENSION(Ninpts,2), INTENT(out)              :: ilonlat
496  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: mindiffLl
497
498! Local
499  INTEGER                                                :: iv,i,j
500  INTEGER                                                :: ierr
501  INTEGER                                                :: Ninpts1
502  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
503  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: fractionlon, fractionlat
504  INTEGER                                                :: dfracdx, dfracdy, fracdx, fracdy
505  CHARACTER(LEN=50)                                      :: fname
506
507!!!!!!! Variables
508! dimx, dimy: dimension length of the target interpolation
509! proj[lon/lat]: longitudes and latitudes of the target interpolation
510! Ninpts: number of points to interpolate
511! [lon/lat]vs: longitudes and latitudes of the points to interpolate
512! mindiff: minimal accepted distance to the target point
513! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues
514! inpt: whether the point has already been localized (1) or not (0)
515! ilonlat: Longitude and Latitude of the input points
516! mindiffLl: minimum difference between target and source longitude/latitude (in degrees)
517
518  fname = 'CoarseInterpolate'
519  Ninpts1 = Ninpts/100
520
521  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
522  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
523
524  PRINT *,'  ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen
525
526  dfracdx = INT(1./percen+1)
527  dfracdy = INT(1./percen+1)
528  fracdx = INT(dimx*percen)
529  fracdy = INT(dimy*percen)
530  PRINT *,'  ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy
531
532  IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
533  ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr)
534  IF (ierr /= 0) THEN
535    PRINT *,TRIM(ErrWarnMsg('err'))
536    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlon' !!"
537    STOP
538  END IF
539  IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
540  ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr)
541  IF (ierr /= 0) THEN
542    PRINT *,TRIM(ErrWarnMsg('err'))
543    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlat' !!"
544    STOP
545  END IF
546
547  DO i=1,dfracdx
548    DO j=1,dfracdy
549      fractionlon(i,j) = projlon(fracdx*(i-1)+1,fracdy*(j-1)+1)
550      fractionlat(i,j) = projlat(fracdx*(i-1)+1,fracdy*(j-1)+1)
551!      PRINT *,'i,j:',i,', ',j,' frac ij:',fracdx*(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),&
552!        ', ',fractionlat(i,j)
553    END DO
554  END DO
555
556!  PRINT *,'  ' // TRIM(fname) // ' fractions of:'
557!  PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')'
558!  DO i=1,dfracdx
559!    PRINT *,fractionlon(i,:)
560!  END DO
561!  PRINT *,' lat_______'
562!  DO i=1,dfracdx
563!    PRINT *,fractionlat(i,:)
564!  END DO
565
566  DO iv=1,Ninpts
567    IF (inpt(iv) == 0) THEN
568      CALL CoarselonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, fractionlon,        &
569        fractionlat, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, ilonlat(iv,:), mindiffLl(iv))
570
571      IF ((mindiffLl(iv) <= mindiff) .AND. .NOT.(ilonlat(iv,1) >= 0 .AND. ilonlat(iv,1) >= 0)) THEN
572        PRINT *,TRIM(ErrWarnMsg('err'))
573        PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
574          ' not relocated !!'
575        PRINT *,'    mindiffl:', mindiffLl(iv), ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
576        STOP
577      END IF
578
579    END IF
580  END DO
581
582END SUBROUTINE CoarseInterpolate
583
584SUBROUTINE CoarseInterpolateExact(projlon, projlat, lonvs, latvs, percen, mindiff, ivar, newvar,      &
585  newvarin, newvarinpt, newvarindiff, dimx, dimy, Ninpts)
586! Subroutine which finds the closest grid point within a projection throughout a first guest
587!   and then whole domain approche from percentages of the whole domain
588
589  USE module_generic
590
591  IMPLICIT NONE
592
593  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
594  INTEGER, INTENT(in)                                    :: dimx, dimy
595  REAL(r_k), DIMENSION(dimx,dimy), INTENT(in)            :: projlon, projlat
596  INTEGER, INTENT(in)                                    :: Ninpts
597  REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: ivar, lonvs, latvs
598  REAL(r_k), INTENT(in)                                  :: mindiff, percen
599  REAL(r_k), DIMENSION(dimx,dimy), INTENT(out)           :: newvar
600  INTEGER, DIMENSION(dimx,dimy), INTENT(out)             :: newvarin
601  INTEGER, DIMENSION(Ninpts), INTENT(out)                :: newvarinpt
602  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: newvarindiff
603
604! Local
605  INTEGER                                                :: iv,i,j
606  INTEGER                                                :: ierr
607  INTEGER, DIMENSION(2)                                  :: ilonlat
608  REAL(r_k)                                              :: mindiffLl
609  INTEGER                                                :: Ninpts1
610  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
611  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: fractionlon, fractionlat
612  INTEGER                                                :: dfracdx, dfracdy, fracdx, fracdy
613  CHARACTER(LEN=50)                                      :: fname
614
615!!!!!!! Variables
616! dimx, dimy: dimension length of the target interpolation
617! proj[lon/lat]: longitudes and latitudes of the target interpolation
618! Ninpts: number of points to interpolate
619! [lon/lat]vs: longitudes and latitudes of the points to interpolate
620! mindiff: minimal accepted distance to the target point
621! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues
622! ivar: values to localize in the target projection
623! newvar: localisation of the [lon/lat]vs point in the target projection
624! newvarin: number of point from the input data
625! newvarinpt: integer value indicating if the value has been already located (0: no, 1: yes)
626! newvarindiff: distance of point from the input data to the closest target point
627! ncid: netCDF output file id
628
629  fname = 'CoarseInterpolateExact'
630  Ninpts1 = Ninpts/100
631
632  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
633  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
634
635  PRINT *,'  ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen
636
637  dfracdx = INT(1./percen+1)
638  dfracdy = INT(1./percen+1)
639  fracdx = INT(dimx*percen)
640  fracdy = INT(dimy*percen)
641  PRINT *,'  ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy
642
643  IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
644  ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr)
645  IF (ierr /= 0) THEN
646    PRINT *,TRIM(ErrWarnMsg('err'))
647    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlon' !!"
648    STOP
649  END IF
650  IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
651  ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr)
652  IF (ierr /= 0) THEN
653    PRINT *,TRIM(ErrWarnMsg('err'))
654    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlat' !!"
655    STOP
656  END IF
657
658  DO i=1,dfracdx
659    DO j=1,dfracdy
660      fractionlon(i,j) = projlon(fracdx*(i-1)+1,fracdy*(j-1)+1)
661      fractionlat(i,j) = projlat(fracdx*(i-1)+1,fracdy*(j-1)+1)
662!      PRINT *,'i,j:',i,', ',j,' frac ij:',fracdx*(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),&
663!        ', ',fractionlat(i,j)
664    END DO
665  END DO
666
667!  PRINT *,'  ' // TRIM(fname) // ' fractions of:'
668!  PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')'
669!  DO i=1,dfracdx
670!    PRINT *,fractionlon(i,:)
671!  END DO
672!  PRINT *,' lat_______'
673!  DO i=1,dfracdx
674!    PRINT *,fractionlat(i,:)
675!  END DO
676
677  DO iv=1,Ninpts
678    IF (newvarinpt(iv) == 0) THEN
679      CALL CoarselonlatFindExact(dimx, dimy, projlon, projlat, extremelon, extremelat, fracdx, fracdy,&
680        fractionlon, fractionlat, iv, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, mindiff,        &
681        ilonlat, mindiffLl)
682
683      IF (mindiffLl >= mindiff) THEN
684!        percendone(iv,Ninpts,0.5,'done:')
685
686        IF (ilonlat(1) >= 0 .AND. ilonlat(1) >= 0) THEN
687          newvar(ilonlat(1),ilonlat(2)) = ivar(iv)
688          newvarin(ilonlat(1),ilonlat(2)) = iv
689          newvarinpt(iv) = 1
690          newvarindiff(iv) = mindiffLl
691!          PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv),        &
692!            ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:',           &
693!            newvarindiff(iv), ' point:',ilonlat   
694        ELSE
695          PRINT *,TRIM(ErrWarnMsg('err'))
696          PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
697            ' not relocated !!'
698          PRINT *,'    mindiffl:', mindiffLl, ' ilon:', ilonlat(1), ' ilat:', ilonlat(2)
699          STOP
700        END IF
701
702!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
703      ELSE
704        PRINT *,TRIM(ErrWarnMsg('err'))
705        PRINT *,'  ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv),   &
706          ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ',       &
707          mindiff, ' !!'
708        PRINT *,'    found minimum difference:', mindiffLl
709        STOP
710      END IF
711    END IF
712  END DO
713
714END SUBROUTINE CoarseInterpolateExact
715
716SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy,     &
717  Ninpts)
718! Subroutine which finds the closest grid point within a projection
719
720  USE module_generic
721
722  IMPLICIT NONE
723
724  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
725  INTEGER, INTENT(in)                                    :: dimx, dimy
726  REAL(r_k), DIMENSION(dimx,dimy), INTENT(in)            :: projlon, projlat
727  INTEGER, INTENT(in)                                    :: Ninpts
728  REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: lonvs, latvs
729  REAL(r_k), INTENT(in)                                  :: mindiff
730  INTEGER, DIMENSION(Ninpts), INTENT(inout)              :: inpt
731  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: diffs
732  INTEGER, DIMENSION(Ninpts,2), INTENT(out)              :: ilonlat
733
734! Local
735  INTEGER                                                :: iv
736  REAL(r_k)                                              :: mindiffLl
737  INTEGER                                                :: Ninpts1
738  REAL(r_k), DIMENSION(dimx,dimy)                        :: difflonlat
739  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
740  CHARACTER(LEN=50)                                      :: fname
741
742!!!!!!! Variables
743! dimx, dimy: dimension length of the target interpolation
744! proj[lon/lat]: longitudes and latitudes of the target interpolation
745! Ninpts: number of points to interpolate
746! [lon/lat]vs: longitudes and latitudes of the points to interpolate
747! mindiff: minimal accepted distance to the target point
748! inpt: whether the point has already been localized
749! diffs: distance of point from the input data to the closest target point
750! ilonlat: longitude and latitude of the point
751! ncid: netCDF output file id
752
753  fname = 'Interpolate'
754  Ninpts1 = Ninpts/100
755
756  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
757  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
758
759  DO iv=1,Ninpts
760    IF (inpt(iv) <= 0) THEN
761! Not using the subroutine, not efficient!
762!      CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv),     &
763!        ilonlat, mindiffLl)
764
765      IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
766        PRINT *, TRIM(ErrWarnMsg('err'))
767        PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
768        PRINT *,'    given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
769        STOP
770      END IF
771      IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
772        PRINT *, TRIM(ErrWarnMsg('err'))
773        PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
774        PRINT *,'    given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
775        STOP
776      END IF
777
778! Find point
779      difflonlat = SQRT((projlon-lonvs(iv))**2. + (projlat-latvs(iv))**2.)
780      mindiffLl = MINVAL(difflonlat)
781      ilonlat(iv,:) = index2DArrayR(difflonlat, dimx, dimy, mindiffLl)
782
783      IF (mindiffLl <= mindiff) THEN
784!        percendone(iv,Ninpts,0.5,'done:')
785
786        IF (ilonlat(iv,1) >= 0 .AND. ilonlat(iv,2) >= 0) THEN
787          diffs(iv) = mindiffLl
788          inpt(iv) = 1
789!          PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv),        &
790!            ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:',           &
791!            newvarindiff(iv), ' point:',ilonlat   
792        ELSE
793          PRINT *,TRIM(ErrWarnMsg('err'))
794          PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
795            ' not relocated !!'
796          PRINT *,'    mindiffl:', mindiffLl, ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
797          STOP
798        END IF
799
800!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
801      ELSE
802! Because doing boxes and Goode is not conitnuos, we should jump this error message
803        PRINT *,TRIM(ErrWarnMsg('err'))
804        PRINT *,'  ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv),   &
805          ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ',       &
806          mindiff, ' !!'
807        PRINT *,'    found minimum difference:', mindiffLl
808        STOP
809      END IF
810    END IF
811  END DO
812
813END SUBROUTINE Interpolate
814
815SUBROUTINE Interpolate1DLl(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy, &
816  Ninpts)
817! Subroutine which finds the closest grid point within a projection with 1D longitudes and latitudes
818
819  USE module_generic
820
821  IMPLICIT NONE
822
823  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
824  INTEGER, INTENT(in)                                    :: dimx, dimy
825  REAL(r_k), DIMENSION(dimx), INTENT(in)                 :: projlon
826  REAL(r_k), DIMENSION(dimy), INTENT(in)                 :: projlat
827  INTEGER, INTENT(in)                                    :: Ninpts
828  REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: lonvs, latvs
829  REAL(r_k), INTENT(in)                                  :: mindiff
830  INTEGER, DIMENSION(Ninpts), INTENT(inout)              :: inpt
831  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: diffs
832  INTEGER, DIMENSION(Ninpts,2), INTENT(out)              :: ilonlat
833
834! Local
835  INTEGER                                                :: iv
836  REAL(r_k)                                              :: mindifflo, mindiffLa, mindiffLl
837  INTEGER                                                :: Ninpts1
838  REAL(r_k), DIMENSION(dimx)                             :: difflon
839  REAL(r_k), DIMENSION(dimy)                             :: difflat
840  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
841  CHARACTER(LEN=50)                                      :: fname
842
843!!!!!!! Variables
844! dimx, dimy: dimension length of the target interpolation
845! proj[lon/lat]: longitudes and latitudes of the target interpolation
846! Ninpts: number of points to interpolate
847! [lon/lat]vs: longitudes and latitudes of the points to interpolate
848! mindiff: minimal accepted distance to the target point
849! inpt: whether the point has already been localized
850! diffs: distance of point from the input data to the closest target point
851! ilonlat: longitude and latitude of the point
852! ncid: netCDF output file id
853
854  fname = 'Interpolate1DLl'
855  Ninpts1 = Ninpts/100
856
857  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
858  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
859
860  DO iv=1,Ninpts
861    IF (inpt(iv) <= 0) THEN
862! Not using the subroutine, not efficient!
863!      CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv),     &
864!        ilonlat, mindiffLl)
865
866      IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
867        PRINT *, TRIM(ErrWarnMsg('err'))
868        PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
869        PRINT *,'    given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
870        STOP
871      END IF
872      IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
873        PRINT *, TRIM(ErrWarnMsg('err'))
874        PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
875        PRINT *,'    given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
876        STOP
877      END IF
878
879! Find point
880      difflon = SQRT((projlon-lonvs(iv))**2.)
881      difflat = SQRT((projlat-latvs(iv))**2.)
882      mindifflo = MINVAL(difflon)
883      mindiffLa = MINVAL(difflat)
884      mindifflL = SQRT(mindifflo*mindifflo + mindiffLa*mindiffLa)
885      ilonlat(iv,1) = index1DArrayR(difflon, dimx, mindifflo)
886      ilonlat(iv,2) = index1DArrayR(difflat, dimy, mindiffLa)
887!      PRINT *,'  Lluis: iv',iv,' lonvs:', lonvs(iv),' latvs:',latvs(iv)
888!      PRINT *,'  Lluis: mindifflo:', mindifflo,' ilonlat(1):',ilonlat(iv,1)
889!      PRINT *,'  Lluis: mindiffLa:', mindiffLa,' ilonlat(2):',ilonlat(iv,2)
890
891
892      IF (mindiffLl <= mindiff) THEN
893!        percendone(iv,Ninpts,0.5,'done:')
894
895        IF (ilonlat(iv,1) >= 1 .AND. ilonlat(iv,2) >= 1) THEN
896          diffs(iv) = mindiffLl
897          inpt(iv) = 1
898!          PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv),        &
899!            ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:',           &
900!            newvarindiff(iv), ' point:',ilonlat   
901        ELSE
902          PRINT *,TRIM(ErrWarnMsg('err'))
903          PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
904            ' not relocated !!'
905          PRINT *,'    mindiffl:', mindiffLl, ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
906          STOP
907        END IF
908
909!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
910      ELSE
911! Because doing boxes and Goode is not conitnuos, we should jump this error message
912        PRINT *,TRIM(ErrWarnMsg('err'))
913        PRINT *,'  ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv),   &
914          ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ',       &
915          mindiff, ' !!'
916        PRINT *,'    found minimum difference:', mindiffLl
917        STOP
918      END IF
919    END IF
920  END DO
921
922END SUBROUTINE Interpolate1DLl
923
924END MODULE module_ForInterpolate
Note: See TracBrowser for help on using the repository browser.