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

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

Final version of the 'EntireGlobalMap?'?
File size: 32.5 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,j
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, ivar, newvar, newvarin, &
481  newvarinpt, newvarindiff, 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)               :: ivar, lonvs, latvs
494  REAL(r_k), INTENT(in)                                  :: mindiff, percen
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
499
500! Local
501  INTEGER                                                :: iv,i,j
502  INTEGER                                                :: ierr
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
508  INTEGER                                                :: dfracdx, dfracdy, fracdx, fracdy
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
531  PRINT *,'  ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen
532
533  dfracdx = INT(1./percen+1)
534  dfracdy = INT(1./percen+1)
535  fracdx = INT(dimx*percen)
536  fracdy = INT(dimy*percen)
537  PRINT *,'  ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy
538
539  IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
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
546  IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
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
553
554  DO i=1,dfracdx
555    DO j=1,dfracdy
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)
560    END DO
561  END DO
562
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
572
573  DO iv=1,Ninpts
574    IF (newvarinpt(iv) == 0) THEN
575      CALL CoarselonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, fractionlon,        &
576        fractionlat, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, ilonlat, mindiffLl)
577
578      IF (mindiffLl <= mindiff) THEN
579!        percendone(iv,Ninpts,0.5,'done:')
580
581        IF (ilonlat(1) >= 0 .AND. ilonlat(1) >= 0) THEN
582          newvar(ilonlat(1),ilonlat(2)) = ivar(iv)
583          newvarin(ilonlat(1),ilonlat(2)) = iv
584          newvarinpt(iv) = 1
585          newvarindiff(iv) = mindiffLl
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   
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()
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
606      END IF
607    END IF
608  END DO
609
610END SUBROUTINE CoarseInterpolate
611
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), DIMENSION(dimx,dimy)                        :: difflonlat
637  REAL(r_k)                                              :: mindiffLl
638  INTEGER                                                :: Ninpts1
639  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
640  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: fractionlon, fractionlat
641  INTEGER                                                :: dfracdx, dfracdy, fracdx, fracdy
642  CHARACTER(LEN=50)                                      :: fname
643
644!!!!!!! Variables
645! dimx, dimy: dimension length of the target interpolation
646! proj[lon/lat]: longitudes and latitudes of the target interpolation
647! Ninpts: number of points to interpolate
648! [lon/lat]vs: longitudes and latitudes of the points to interpolate
649! mindiff: minimal accepted distance to the target point
650! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues
651! ivar: values to localize in the target projection
652! newvar: localisation of the [lon/lat]vs point in the target projection
653! newvarin: number of point from the input data
654! newvarinpt: integer value indicating if the value has been already located (0: no, 1: yes)
655! newvarindiff: distance of point from the input data to the closest target point
656! ncid: netCDF output file id
657
658  fname = 'CoarseInterpolateExact'
659  Ninpts1 = Ninpts/100
660
661  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
662  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
663
664  PRINT *,'  ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen
665
666  dfracdx = INT(1./percen+1)
667  dfracdy = INT(1./percen+1)
668  fracdx = INT(dimx*percen)
669  fracdy = INT(dimy*percen)
670  PRINT *,'  ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy
671
672  IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
673  ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr)
674  IF (ierr /= 0) THEN
675    PRINT *,TRIM(ErrWarnMsg('err'))
676    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlon' !!"
677    STOP
678  END IF
679  IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
680  ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr)
681  IF (ierr /= 0) THEN
682    PRINT *,TRIM(ErrWarnMsg('err'))
683    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlat' !!"
684    STOP
685  END IF
686
687  DO i=1,dfracdx
688    DO j=1,dfracdy
689      fractionlon(i,j) = projlon(fracdx*(i-1)+1,fracdy*(j-1)+1)
690      fractionlat(i,j) = projlat(fracdx*(i-1)+1,fracdy*(j-1)+1)
691!      PRINT *,'i,j:',i,', ',j,' frac ij:',fracdx*(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),&
692!        ', ',fractionlat(i,j)
693    END DO
694  END DO
695
696!  PRINT *,'  ' // TRIM(fname) // ' fractions of:'
697!  PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')'
698!  DO i=1,dfracdx
699!    PRINT *,fractionlon(i,:)
700!  END DO
701!  PRINT *,' lat_______'
702!  DO i=1,dfracdx
703!    PRINT *,fractionlat(i,:)
704!  END DO
705
706  DO iv=1,Ninpts
707    IF (newvarinpt(iv) == 0) THEN
708      CALL CoarselonlatFindExact(dimx, dimy, projlon, projlat, extremelon, extremelat, fracdx, fracdy,&
709        fractionlon, fractionlat, iv, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, mindiff,        &
710        ilonlat, mindiffLl)
711
712      IF (mindiffLl >= mindiff) THEN
713!        percendone(iv,Ninpts,0.5,'done:')
714
715        IF (ilonlat(1) >= 0 .AND. ilonlat(1) >= 0) THEN
716          newvar(ilonlat(1),ilonlat(2)) = ivar(iv)
717          newvarin(ilonlat(1),ilonlat(2)) = iv
718          newvarinpt(iv) = 1
719          newvarindiff(iv) = mindiffLl
720!          PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv),        &
721!            ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:',           &
722!            newvarindiff(iv), ' point:',ilonlat   
723        ELSE
724          PRINT *,TRIM(ErrWarnMsg('err'))
725          PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
726            ' not relocated !!'
727          PRINT *,'    mindiffl:', mindiffLl, ' ilon:', ilonlat(1), ' ilat:', ilonlat(2)
728          STOP
729        END IF
730
731!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
732      ELSE
733        PRINT *,TRIM(ErrWarnMsg('err'))
734        PRINT *,'  ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv),   &
735          ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ',       &
736          mindiff, ' !!'
737        PRINT *,'    found minimum difference:', mindiffLl
738        STOP
739      END IF
740    END IF
741  END DO
742
743END SUBROUTINE CoarseInterpolateExact
744
745SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, ivar, newvarinpt,                     &
746  newvarindiff, ilonlat, dimx, dimy, Ninpts)
747! Subroutine which finds the closest grid point within a projection
748
749  USE module_generic
750
751  IMPLICIT NONE
752
753  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
754  INTEGER, INTENT(in)                                    :: dimx, dimy
755  REAL(r_k), DIMENSION(dimx,dimy), INTENT(in)            :: projlon, projlat
756  INTEGER, INTENT(in)                                    :: Ninpts
757  REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: ivar, lonvs, latvs
758  REAL(r_k), INTENT(in)                                  :: mindiff
759  INTEGER, DIMENSION(Ninpts), INTENT(out)                :: newvarinpt
760  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: newvarindiff
761  INTEGER, DIMENSION(Ninpts,2), INTENT(out)              :: ilonlat
762
763! Local
764  INTEGER                                                :: iv,i,j
765  INTEGER                                                :: ierr
766  REAL(r_k)                                              :: mindiffLl
767  INTEGER                                                :: Ninpts1
768  REAL(r_k), DIMENSION(dimx,dimy)                        :: difflonlat
769  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
770  CHARACTER(LEN=50)                                      :: fname
771
772!!!!!!! Variables
773! dimx, dimy: dimension length of the target interpolation
774! proj[lon/lat]: longitudes and latitudes of the target interpolation
775! Ninpts: number of points to interpolate
776! [lon/lat]vs: longitudes and latitudes of the points to interpolate
777! mindiff: minimal accepted distance to the target point
778! ivar: values to localize in the target projection
779! newvar: localisation of the [lon/lat]vs point in the target projection
780! newvarin: number of point from the input data
781! newvarinpt: integer value indicating if the value has been already located (0: no, 1: yes)
782! newvarindiff: distance of point from the input data to the closest target point
783! ilonlat: longitude and latitude of the point
784! ncid: netCDF output file id
785
786  fname = 'Interpolate'
787  Ninpts1 = Ninpts/100
788
789  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
790  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
791
792  DO iv=1,Ninpts
793    IF (newvarinpt(iv) == 0) THEN
794! Not using the subroutine, not efficient!
795!      CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv),     &
796!        ilonlat, mindiffLl)
797
798      IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
799        PRINT *, TRIM(ErrWarnMsg('err'))
800        PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
801        PRINT *,'    given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
802        STOP
803      END IF
804      IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
805        PRINT *, TRIM(ErrWarnMsg('err'))
806        PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
807        PRINT *,'    given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
808        STOP
809      END IF
810
811! Find point
812      difflonlat = SQRT((projlon-lonvs(iv))**2. + (projlat-latvs(iv))**2.)
813      mindiffLl = MINVAL(difflonlat)
814      ilonlat(iv,:) = index2DArrayR(difflonlat, dimx, dimy, mindiffLl)
815
816      IF (mindiffLl <= mindiff) THEN
817!        percendone(iv,Ninpts,0.5,'done:')
818
819        IF (ilonlat(iv,1) >= 0 .AND. ilonlat(iv,2) >= 0) THEN
820          newvarinpt(iv) = 1
821          newvarindiff(iv) = mindiffLl
822!          PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv),        &
823!            ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:',           &
824!            newvarindiff(iv), ' point:',ilonlat   
825        ELSE
826          PRINT *,TRIM(ErrWarnMsg('err'))
827          PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
828            ' not relocated !!'
829          PRINT *,'    mindiffl:', mindiffLl, ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
830          STOP
831        END IF
832
833!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
834      ELSE
835! Because doing boxes and Goode is not conitnuos, we should jump this error message
836        PRINT *,TRIM(ErrWarnMsg('err'))
837        PRINT *,'  ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv),   &
838          ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ',       &
839          mindiff, ' !!'
840        PRINT *,'    found minimum difference:', mindiffLl
841        STOP
842      END IF
843    ELSE
844      PRINT *,TRIM(ErrWarnMsg('err'))
845      PRINT *,'  ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv),     &
846        ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ',         &
847        mindiff, ' !!'
848      PRINT *,'    found minimum difference:', mindiffLl
849      STOP
850    END IF
851  END DO
852
853END SUBROUTINE Interpolate
854
855END MODULE module_ForInterpolate
Note: See TracBrowser for help on using the repository browser.