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

Last change on this file since 762 was 742, checked in by lfita, 9 years ago

Final seem to work version II
File size: 31.7 KB
RevLine 
[715]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
[735]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
[742]218  INTEGER                                                :: i
[735]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
[733]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
[730]480SUBROUTINE CoarseInterpolate(projlon, projlat, lonvs, latvs, percen, mindiff, ivar, newvar, newvarin, &
481  newvarinpt, newvarindiff, dimx, dimy, Ninpts)
[715]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)               :: ivar, lonvs, latvs
[733]494  REAL(r_k), INTENT(in)                                  :: mindiff, percen
[715]495  REAL(r_k), DIMENSION(dimx,dimy), INTENT(out)           :: newvar
496  INTEGER, DIMENSION(dimx,dimy), INTENT(out)             :: newvarin
497  INTEGER, DIMENSION(Ninpts), INTENT(out)                :: newvarinpt
498  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: newvarindiff
[728]499
[715]500! Local
[728]501  INTEGER                                                :: iv,i,j
502  INTEGER                                                :: ierr
[715]503  INTEGER, DIMENSION(2)                                  :: ilonlat
504  REAL(r_k)                                              :: mindiffLl
505  INTEGER                                                :: Ninpts1
506  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
507  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: fractionlon, fractionlat
[728]508  INTEGER                                                :: dfracdx, dfracdy, fracdx, fracdy
[715]509  CHARACTER(LEN=50)                                      :: fname
510
511!!!!!!! Variables
512! dimx, dimy: dimension length of the target interpolation
513! proj[lon/lat]: longitudes and latitudes of the target interpolation
514! Ninpts: number of points to interpolate
515! [lon/lat]vs: longitudes and latitudes of the points to interpolate
516! mindiff: minimal accepted distance to the target point
517! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues
518! ivar: values to localize in the target projection
519! newvar: localisation of the [lon/lat]vs point in the target projection
520! newvarin: number of point from the input data
521! newvarinpt: integer value indicating if the value has been already located (0: no, 1: yes)
522! newvarindiff: distance of point from the input data to the closest target point
523! ncid: netCDF output file id
524
525  fname = 'CoarseInterpolate'
526  Ninpts1 = Ninpts/100
527
528  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
529  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
530
[729]531  PRINT *,'  ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen
[728]532
[729]533  dfracdx = INT(1./percen+1)
534  dfracdy = INT(1./percen+1)
[715]535  fracdx = INT(dimx*percen)
536  fracdy = INT(dimy*percen)
[729]537  PRINT *,'  ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy
[715]538
539  IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
[728]540  ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr)
541  IF (ierr /= 0) THEN
542    PRINT *,TRIM(ErrWarnMsg('err'))
543    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlon' !!"
544    STOP
545  END IF
[715]546  IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
[728]547  ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr)
548  IF (ierr /= 0) THEN
549    PRINT *,TRIM(ErrWarnMsg('err'))
550    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlat' !!"
551    STOP
552  END IF
[715]553
[728]554  DO i=1,dfracdx
555    DO j=1,dfracdy
[729]556      fractionlon(i,j) = projlon(fracdx*(i-1)+1,fracdy*(j-1)+1)
557      fractionlat(i,j) = projlat(fracdx*(i-1)+1,fracdy*(j-1)+1)
558!      PRINT *,'i,j:',i,', ',j,' frac ij:',fracdx*(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),&
559!        ', ',fractionlat(i,j)
[728]560    END DO
561  END DO
[715]562
[733]563!  PRINT *,'  ' // TRIM(fname) // ' fractions of:'
564!  PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')'
565!  DO i=1,dfracdx
566!    PRINT *,fractionlon(i,:)
567!  END DO
568!  PRINT *,' lat_______'
569!  DO i=1,dfracdx
570!    PRINT *,fractionlat(i,:)
571!  END DO
[728]572
[715]573  DO iv=1,Ninpts
574    IF (newvarinpt(iv) == 0) THEN
575      CALL CoarselonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, fractionlon,        &
[728]576        fractionlat, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, ilonlat, mindiffLl)
[715]577
578      IF (mindiffLl <= mindiff) THEN
579!        percendone(iv,Ninpts,0.5,'done:')
580
581        IF (ilonlat(1) >= 0 .AND. ilonlat(1) >= 0) THEN
[730]582          newvar(ilonlat(1),ilonlat(2)) = ivar(iv)
583          newvarin(ilonlat(1),ilonlat(2)) = iv
584          newvarinpt(iv) = 1
585          newvarindiff(iv) = mindiffLl
[716]586!          PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv),        &
587!            ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:',           &
588!            newvarindiff(iv), ' point:',ilonlat   
[715]589        ELSE
590          PRINT *,TRIM(ErrWarnMsg('err'))
591          PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
592            ' not relocated !!'
593          PRINT *,'    mindiffl:', mindiffLl, ' ilon:', ilonlat(1), ' ilat:', ilonlat(2)
594          STOP
595        END IF
596
597!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
[716]598!      ELSE
599! Because doing boxes and Goode is not conitnuos, we should jump this error message
600!        PRINT *,TRIM(ErrWarnMsg('err'))
601!        PRINT *,'  ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv),   &
602!          ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ',       &
603!          mindiff, ' !!'
604!        PRINT *,'    found minimum difference:', mindiffLl
605!        STOP
[715]606      END IF
607    END IF
608  END DO
609
610END SUBROUTINE CoarseInterpolate
611
[735]612SUBROUTINE CoarseInterpolateExact(projlon, projlat, lonvs, latvs, percen, mindiff, ivar, newvar,      &
613  newvarin, newvarinpt, newvarindiff, dimx, dimy, Ninpts)
614! Subroutine which finds the closest grid point within a projection throughout a first guest
615!   and then whole domain approche from percentages of the whole domain
616
617  USE module_generic
618
619  IMPLICIT NONE
620
621  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
622  INTEGER, INTENT(in)                                    :: dimx, dimy
623  REAL(r_k), DIMENSION(dimx,dimy), INTENT(in)            :: projlon, projlat
624  INTEGER, INTENT(in)                                    :: Ninpts
625  REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: ivar, lonvs, latvs
626  REAL(r_k), INTENT(in)                                  :: mindiff, percen
627  REAL(r_k), DIMENSION(dimx,dimy), INTENT(out)           :: newvar
628  INTEGER, DIMENSION(dimx,dimy), INTENT(out)             :: newvarin
629  INTEGER, DIMENSION(Ninpts), INTENT(out)                :: newvarinpt
630  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: newvarindiff
631
632! Local
633  INTEGER                                                :: iv,i,j
634  INTEGER                                                :: ierr
635  INTEGER, DIMENSION(2)                                  :: ilonlat
636  REAL(r_k)                                              :: mindiffLl
637  INTEGER                                                :: Ninpts1
638  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
639  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: fractionlon, fractionlat
640  INTEGER                                                :: dfracdx, dfracdy, fracdx, fracdy
641  CHARACTER(LEN=50)                                      :: fname
642
643!!!!!!! Variables
644! dimx, dimy: dimension length of the target interpolation
645! proj[lon/lat]: longitudes and latitudes of the target interpolation
646! Ninpts: number of points to interpolate
647! [lon/lat]vs: longitudes and latitudes of the points to interpolate
648! mindiff: minimal accepted distance to the target point
649! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues
650! ivar: values to localize in the target projection
651! newvar: localisation of the [lon/lat]vs point in the target projection
652! newvarin: number of point from the input data
653! newvarinpt: integer value indicating if the value has been already located (0: no, 1: yes)
654! newvarindiff: distance of point from the input data to the closest target point
655! ncid: netCDF output file id
656
657  fname = 'CoarseInterpolateExact'
658  Ninpts1 = Ninpts/100
659
660  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
661  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
662
663  PRINT *,'  ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen
664
665  dfracdx = INT(1./percen+1)
666  dfracdy = INT(1./percen+1)
667  fracdx = INT(dimx*percen)
668  fracdy = INT(dimy*percen)
669  PRINT *,'  ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy
670
671  IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
672  ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr)
673  IF (ierr /= 0) THEN
674    PRINT *,TRIM(ErrWarnMsg('err'))
675    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlon' !!"
676    STOP
677  END IF
678  IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
679  ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr)
680  IF (ierr /= 0) THEN
681    PRINT *,TRIM(ErrWarnMsg('err'))
682    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlat' !!"
683    STOP
684  END IF
685
686  DO i=1,dfracdx
687    DO j=1,dfracdy
688      fractionlon(i,j) = projlon(fracdx*(i-1)+1,fracdy*(j-1)+1)
689      fractionlat(i,j) = projlat(fracdx*(i-1)+1,fracdy*(j-1)+1)
690!      PRINT *,'i,j:',i,', ',j,' frac ij:',fracdx*(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),&
691!        ', ',fractionlat(i,j)
692    END DO
693  END DO
694
695!  PRINT *,'  ' // TRIM(fname) // ' fractions of:'
696!  PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')'
697!  DO i=1,dfracdx
698!    PRINT *,fractionlon(i,:)
699!  END DO
700!  PRINT *,' lat_______'
701!  DO i=1,dfracdx
702!    PRINT *,fractionlat(i,:)
703!  END DO
704
705  DO iv=1,Ninpts
706    IF (newvarinpt(iv) == 0) THEN
707      CALL CoarselonlatFindExact(dimx, dimy, projlon, projlat, extremelon, extremelat, fracdx, fracdy,&
708        fractionlon, fractionlat, iv, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, mindiff,        &
709        ilonlat, mindiffLl)
710
711      IF (mindiffLl >= mindiff) THEN
712!        percendone(iv,Ninpts,0.5,'done:')
713
714        IF (ilonlat(1) >= 0 .AND. ilonlat(1) >= 0) THEN
715          newvar(ilonlat(1),ilonlat(2)) = ivar(iv)
716          newvarin(ilonlat(1),ilonlat(2)) = iv
717          newvarinpt(iv) = 1
718          newvarindiff(iv) = mindiffLl
719!          PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv),        &
720!            ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:',           &
721!            newvarindiff(iv), ' point:',ilonlat   
722        ELSE
723          PRINT *,TRIM(ErrWarnMsg('err'))
724          PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
725            ' not relocated !!'
726          PRINT *,'    mindiffl:', mindiffLl, ' ilon:', ilonlat(1), ' ilat:', ilonlat(2)
727          STOP
728        END IF
729
730!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
731      ELSE
732        PRINT *,TRIM(ErrWarnMsg('err'))
733        PRINT *,'  ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv),   &
734          ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ',       &
735          mindiff, ' !!'
736        PRINT *,'    found minimum difference:', mindiffLl
737        STOP
738      END IF
739    END IF
740  END DO
741
742END SUBROUTINE CoarseInterpolateExact
743
[742]744SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy,     &
745  Ninpts)
[733]746! Subroutine which finds the closest grid point within a projection
747
748  USE module_generic
749
750  IMPLICIT NONE
751
752  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
753  INTEGER, INTENT(in)                                    :: dimx, dimy
754  REAL(r_k), DIMENSION(dimx,dimy), INTENT(in)            :: projlon, projlat
755  INTEGER, INTENT(in)                                    :: Ninpts
[742]756  REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: lonvs, latvs
[733]757  REAL(r_k), INTENT(in)                                  :: mindiff
[742]758  INTEGER, DIMENSION(Ninpts), INTENT(inout)              :: inpt
759  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: diffs
[740]760  INTEGER, DIMENSION(Ninpts,2), INTENT(out)              :: ilonlat
[733]761
762! Local
[742]763  INTEGER                                                :: iv
[733]764  REAL(r_k)                                              :: mindiffLl
765  INTEGER                                                :: Ninpts1
[740]766  REAL(r_k), DIMENSION(dimx,dimy)                        :: difflonlat
[733]767  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
768  CHARACTER(LEN=50)                                      :: fname
769
770!!!!!!! Variables
771! dimx, dimy: dimension length of the target interpolation
772! proj[lon/lat]: longitudes and latitudes of the target interpolation
773! Ninpts: number of points to interpolate
774! [lon/lat]vs: longitudes and latitudes of the points to interpolate
775! mindiff: minimal accepted distance to the target point
[742]776! inpt: whether the point has already been localized
777! diffs: distance of point from the input data to the closest target point
[740]778! ilonlat: longitude and latitude of the point
[733]779! ncid: netCDF output file id
780
781  fname = 'Interpolate'
782  Ninpts1 = Ninpts/100
783
784  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
785  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
786
787  DO iv=1,Ninpts
[742]788    IF (inpt(iv) <= 0) THEN
[733]789! Not using the subroutine, not efficient!
790!      CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv),     &
791!        ilonlat, mindiffLl)
792
793      IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
794        PRINT *, TRIM(ErrWarnMsg('err'))
795        PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
796        PRINT *,'    given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
797        STOP
798      END IF
799      IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
800        PRINT *, TRIM(ErrWarnMsg('err'))
801        PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
802        PRINT *,'    given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
803        STOP
804      END IF
805
806! Find point
807      difflonlat = SQRT((projlon-lonvs(iv))**2. + (projlat-latvs(iv))**2.)
808      mindiffLl = MINVAL(difflonlat)
[740]809      ilonlat(iv,:) = index2DArrayR(difflonlat, dimx, dimy, mindiffLl)
[733]810
811      IF (mindiffLl <= mindiff) THEN
812!        percendone(iv,Ninpts,0.5,'done:')
813
[740]814        IF (ilonlat(iv,1) >= 0 .AND. ilonlat(iv,2) >= 0) THEN
[742]815          diffs(iv) = mindiffLl
816          inpt(iv) = 1
[733]817!          PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv),        &
818!            ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:',           &
819!            newvarindiff(iv), ' point:',ilonlat   
820        ELSE
821          PRINT *,TRIM(ErrWarnMsg('err'))
822          PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
823            ' not relocated !!'
[740]824          PRINT *,'    mindiffl:', mindiffLl, ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
[733]825          STOP
826        END IF
827
828!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
[735]829      ELSE
[733]830! Because doing boxes and Goode is not conitnuos, we should jump this error message
[735]831        PRINT *,TRIM(ErrWarnMsg('err'))
832        PRINT *,'  ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv),   &
833          ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ',       &
834          mindiff, ' !!'
835        PRINT *,'    found minimum difference:', mindiffLl
836        STOP
[733]837      END IF
838    END IF
839  END DO
840
841END SUBROUTINE Interpolate
842
[715]843END MODULE module_ForInterpolate
Note: See TracBrowser for help on using the repository browser.