! Module to interpolate values from a giving projection ! To be included in a python ! f2py -m module_ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c module_ForInterpolate.F90 module_generic.F90 >& run_f2py.log MODULE module_ForInterpolate CONTAINS SUBROUTINE CoarselonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, fraclon, fraclat, lonv, latv, per, & Nperx, Npery, ilonlat, mindiffLl) ! Function to search a given value from a coarser version of the data USE module_generic IMPLICIT NONE INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: dx, dy REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: ilon, ilat REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in) :: fraclon, fraclat REAL(r_k), INTENT(in) :: lonv, latv, per REAL(r_k), DIMENSION(2), INTENT(in) :: nxlon, nxlat INTEGER, INTENT(in) :: Nperx, Npery INTEGER, DIMENSION(2), INTENT(out) :: ilonlat REAL(r_k), INTENT(out) :: mindiffLl ! Local REAL(r_k), DIMENSION(Nperx,Npery) :: difffraclonlat REAL(r_k) :: mindifffracLl INTEGER, DIMENSION(2) :: ilonlatfrac INTEGER :: ixbeg, ixend, iybeg, iyend INTEGER :: fracx, fracy REAL(r_k) :: fraclonv, fraclatv REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: difflonlat, lon, lat CHARACTER(LEN=50) :: fname ! Variables ! ilon, ilat: original 2D matrices with the longitudes and the latitudes ! lonv, latv: longitude and latitude to find ! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat ! per: fraction of the whole domain (as percentage) ! Nper[x/y]: period (as fraction over 1) of the fractions of the original grid to use to explore ! fraclon, fraclat: longitude and latitude fractional matricies to perform the first guess fname = 'CoarselonlatFind' IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!' PRINT *,' given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )' STOP END IF IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!' PRINT *,' given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )' STOP END IF ! Initializing variables ixbeg = 0 ixend = 0 iybeg = 0 iyend = 0 fracx = int(dx*per) fracy = int(dy*per) ! PRINT *,'fraclon _______' ! PRINT *,fraclon ! PRINT *,'fraclat _______' ! PRINT *,fraclat ! Fraction point difffraclonlat = SQRT((fraclon-lonv)**2. + (fraclat-latv)**2.) mindifffracLl = MINVAL(difffraclonlat) ilonlatfrac = index2DArrayR(difffraclonlat, Nperx, Npery, mindifffracLl) ! PRINT *, 'mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac ! PRINT *, 'frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,', & ! fraclat(ilonlatfrac(1),ilonlatfrac(2)) ! PRINT *, 'values lon, lat:', lonv, latv ! Providing fraction range fraclonv = fraclon(ilonlatfrac(1),ilonlatfrac(2)) fraclatv = fraclat(ilonlatfrac(1),ilonlatfrac(2)) IF (fraclonv >= lonv .AND. fraclatv >= latv) THEN IF (ilonlatfrac(1) > 0) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE ixbeg = 0 ixend = fracx+1 END IF IF (ilonlatfrac(2) > 0) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 0 iyend = fracy+1 END IF ELSE IF (fraclonv < lonv .AND. fraclatv >= latv) THEN IF (ilonlatfrac(1) < Nperx) THEN IF (ilonlatfrac(1) /= 0) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE ixbeg = 0 ixend = fracx+1 END IF ELSE ixbeg = Nperx*fracx ixend = dx+1 END IF IF (ilonlatfrac(2) > 0) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 0 iyend = fracy+1 END IF ELSE IF (fraclonv < lonv .AND. fraclatv < latv) THEN IF (ilonlatfrac(1) < Nperx) THEN IF (ilonlatfrac(1) /= 0) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE ixbeg = 0 ixend = fracx+1 END IF ELSE ixbeg = Nperx*fracx ixend = dx+1 ENDIF IF (ilonlatfrac(2) < Npery) THEN IF (ilonlatfrac(2) /= 0) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 0 iyend = fracy+1 END IF ELSE iybeg = Npery*fracy iyend = dy+1 END IF ELSE IF (fraclonv >= lonv .AND. fraclatv < latv) THEN IF (ilonlatfrac(1) > 0) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE ixbeg = 0 ixend = fracx+1 END IF IF (ilonlatfrac(2) < Npery) THEN IF (ilonlatfrac(2) /= 0) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 0 iyend = fracy+1 END IF ELSE iybeg = Npery*fracy iyend = dy+1 END IF END IF IF (ALLOCATED(lon)) DEALLOCATE(lon) ALLOCATE(lon(ixend-ixbeg+1, iyend-iybeg+1)) IF (ALLOCATED(lat)) DEALLOCATE(lat) ALLOCATE(lat(ixend-ixbeg+1, iyend-iybeg+1)) IF (ALLOCATED(difflonlat)) DEALLOCATE(difflonlat) ALLOCATE(difflonlat(ixend-ixbeg+1, iyend-iybeg+1)) lon = ilon(ixbeg:ixend,iybeg:iyend) lat = ilat(ixbeg:ixend,iybeg:iyend) ! print *,'lon _______' ! print *,lon ! print *,'lat _______' ! print *,lat ! Find point difflonlat = SQRT((lon-lonv)**2. + (lat-latv)**2.) mindiffLl = MINVAL(difflonlat) ilonlat = index2DArrayR(difflonlat, ixend-ixbeg+1, iyend-iybeg+1, mindiffLl) ilonlat(1) = ilonlat(1) + ixbeg ilonlat(2) = ilonlat(2) + iybeg ! PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon ! PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2)) RETURN END SUBROUTINE CoarselonlatFind SUBROUTINE CoarselonlatFindExact(dx, dy, ilon, ilat, nxlon, nxlat, fracx, fracy, fraclon, fraclat, & iv, lonv, latv, per, Nperx, Npery, mindiff, ilonlat, mindiffLl) ! Function to search a given value from a coarser version of the data USE module_generic IMPLICIT NONE INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: dx, dy, iv REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: ilon, ilat INTEGER, INTENT(in) :: fracx, fracy REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in) :: fraclon, fraclat REAL(r_k), INTENT(in) :: lonv, latv, per, mindiff REAL(r_k), DIMENSION(2), INTENT(in) :: nxlon, nxlat INTEGER, INTENT(in) :: Nperx, Npery INTEGER, DIMENSION(2), INTENT(out) :: ilonlat REAL(r_k), INTENT(out) :: mindiffLl ! Local INTEGER :: i REAL(r_k), DIMENSION(Nperx,Npery) :: difffraclonlat REAL(r_k) :: mindifffracLl INTEGER, DIMENSION(2) :: ilonlatfrac INTEGER :: ixbeg, ixend, iybeg, iyend REAL(r_k) :: fraclonv, fraclatv REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: difflonlat, lon, lat CHARACTER(LEN=50) :: fname ! Variables ! ilon, ilat: original 2D matrices with the longitudes and the latitudes ! lonv, latv: longitude and latitude to find ! iv: point in the input data ! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat ! per: fraction of the whole domain (as percentage) ! Nper[x/y]: period (as fraction over 1) of the fractions of the original grid to use to explore ! frac[x/y]: Number of grid points for each fraction ! fraclon, fraclat: longitude and latitude fractional matricies to perform the first guess ! mindiff: authorized minimal distance between input and interpolated point ! ilonlat: grid point on the total lon,lat matrix ! mindiffLl: distance between input and interpolated point fname = 'CoarselonlatFindExact' IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!' PRINT *,' given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )' STOP END IF IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!' PRINT *,' given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )' STOP END IF ! Initializing variables ixbeg = 0 ixend = 0 iybeg = 0 iyend = 0 ! Fraction point difffraclonlat = SQRT((fraclon-lonv)**2. + (fraclat-latv)**2.) mindifffracLl = MINVAL(difffraclonlat) ilonlatfrac = index2DArrayR(difffraclonlat, Nperx, Npery, mindifffracLl) ! PRINT *, 'mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac ! PRINT *, 'frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,', & ! fraclat(ilonlatfrac(1),ilonlatfrac(2)) ! PRINT *, 'values lon, lat:', lonv, latv ! Providing fraction range fraclonv = fraclon(ilonlatfrac(1),ilonlatfrac(2)) fraclatv = fraclat(ilonlatfrac(1),ilonlatfrac(2)) IF (fraclonv >= lonv .AND. fraclatv >= latv) THEN PRINT *,'Lluis!',fraclonv, '>=', lonv,'&', fraclatv, '>=', latv IF (ilonlatfrac(1) > 1) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE PRINT *,'Lluis 2' ixbeg = 1 ixend = fracx+1 END IF IF (ilonlatfrac(2) > 1) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 1 iyend = fracy+1 END IF ELSE IF (fraclonv < lonv .AND. fraclatv >= latv) THEN PRINT *,'Lluis!',fraclonv, '<', lonv,'&', fraclatv, '>=', latv IF (ilonlatfrac(1) < Nperx) THEN PRINT *,'Lluis 2' IF (ilonlatfrac(1) /= 1) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE ixbeg = 1 ixend = fracx+1 END IF ELSE ixbeg = Nperx*fracx ixend = dx+1 END IF IF (ilonlatfrac(2) > 1) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 1 iyend = fracy+1 END IF ELSE IF (fraclonv < lonv .AND. fraclatv < latv) THEN PRINT *,'Lluis!',fraclonv, '<', lonv,'&', fraclatv, '<', latv IF (ilonlatfrac(1) < Nperx) THEN IF (ilonlatfrac(1) /= 1) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE ixbeg = 1 ixend = fracx+1 END IF ELSE ixbeg = Nperx*fracx ixend = dx+1 ENDIF IF (ilonlatfrac(2) < Npery) THEN IF (ilonlatfrac(2) /= 1) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 1 iyend = fracy+1 END IF ELSE iybeg = Npery*fracy iyend = dy+1 END IF ELSE IF (fraclonv >= lonv .AND. fraclatv < latv) THEN PRINT *,'Llui!',fraclonv, '>=', lonv,'&', fraclatv, '<', latv IF (ilonlatfrac(1) > 1) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE ixbeg = 1 ixend = fracx+1 END IF IF (ilonlatfrac(2) < Npery) THEN IF (ilonlatfrac(2) /= 1) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 1 iyend = fracy+1 END IF ELSE iybeg = Npery*fracy iyend = dy+1 END IF END IF IF (ALLOCATED(lon)) DEALLOCATE(lon) ALLOCATE(lon(ixend-ixbeg+1, iyend-iybeg+1)) IF (ALLOCATED(lat)) DEALLOCATE(lat) ALLOCATE(lat(ixend-ixbeg+1, iyend-iybeg+1)) IF (ALLOCATED(difflonlat)) DEALLOCATE(difflonlat) ALLOCATE(difflonlat(ixend-ixbeg+1, iyend-iybeg+1)) lon = ilon(ixbeg:ixend,iybeg:iyend) lat = ilat(ixbeg:ixend,iybeg:iyend) ! print *,'lon _______' ! print *,lon ! print *,'lat _______' ! print *,lat ! Find point difflonlat = SQRT((lon-lonv)**2. + (lat-latv)**2.) mindiffLl = MINVAL(difflonlat) IF (mindiffLl > mindiff) THEN difflonlat = SQRT((lon-lonv)**2. + (lat-latv)**2.) mindiffLl = MINVAL(difflonlat) END IF IF (mindiffLl > mindiff) THEN PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': not equivalent point closer than:',mindiff,' found!!' PRINT *,' at input point iv:', iv,' lon/lat:', lonv,', ',latv,' distance:',mindiffLl PRINT *,' Fraction values _______ (',Nperx,', ',Npery ,')' PRINT *,' fraclon' DO i=1, Nperx PRINT *,' ',fraclon(i,:) END DO PRINT *,' fraclat' DO i=1, Nperx PRINT *,' ',fraclat(i,:) END DO PRINT *,' frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,', & fraclat(ilonlatfrac(1),ilonlatfrac(2)) PRINT *,' mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac PRINT *,' Coarse values _______' PRINT *,' indices. x:', ixbeg, ', ', ixend, ' y:', iybeg, ', ', iyend PRINT *,' lon range:', '(',ilon(ixbeg,iybeg),', ',ilon(ixend,iyend),')' PRINT *,' lat range:', '(',ilat(ixbeg,iybeg),', ',ilat(ixend,iyend),')' PRINT *,' lon', UBOUND(lon) DO i=1, ixend-ixbeg+1 PRINT *,' ',lon(i,:) END DO PRINT *,' lat', UBOUND(lat) DO i=1, ixend-ixbeg+1 PRINT *,' ',lat(i,:) END DO STOP END IF ilonlat = index2DArrayR(difflonlat, ixend-ixbeg+1, iyend-iybeg+1, mindiffLl) ilonlat(1) = ilonlat(1) + ixbeg ilonlat(2) = ilonlat(2) + iybeg ! PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon ! PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2)) RETURN END SUBROUTINE CoarselonlatFindExact SUBROUTINE lonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, lonv, latv, ilonlat, mindiffLl) ! Function to search a given value from a coarser version of the data USE module_generic IMPLICIT NONE INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: dx, dy REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: ilon, ilat REAL(r_k), INTENT(in) :: lonv, latv REAL(r_k), DIMENSION(2), INTENT(in) :: nxlon, nxlat INTEGER, DIMENSION(2), INTENT(out) :: ilonlat REAL(r_k), INTENT(out) :: mindiffLl ! Local REAL(r_k), DIMENSION(dx,dy) :: difflonlat CHARACTER(LEN=50) :: fname ! Variables ! ilon, ilat: original 2D matrices with the longitudes and the latitudes ! lonv, latv: longitude and latitude to find ! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat fname = 'lonlatFind' IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!' PRINT *,' given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )' STOP END IF IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!' PRINT *,' given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )' STOP END IF ! Find point difflonlat = SQRT((ilon-lonv)**2. + (ilat-latv)**2.) mindiffLl = MINVAL(difflonlat) ilonlat = index2DArrayR(difflonlat, dx, dy, mindiffLl) ! PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon ! PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2)) RETURN END SUBROUTINE lonlatFind SUBROUTINE CoarseInterpolate(projlon, projlat, lonvs, latvs, percen, mindiff, inpt, ilonlat, & mindiffLl, dimx, dimy, Ninpts) ! Subroutine which finds the closest grid point within a projection throughout a first guest ! approche from percentages of the whole domain USE module_generic IMPLICIT NONE INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: dimx, dimy REAL(r_k), DIMENSION(dimx,dimy), INTENT(in) :: projlon, projlat INTEGER, INTENT(in) :: Ninpts REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: inpt, lonvs, latvs REAL(r_k), INTENT(in) :: mindiff, percen INTEGER, DIMENSION(Ninpts,2), INTENT(out) :: ilonlat REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: mindiffLl ! Local INTEGER :: iv,i,j INTEGER :: ierr INTEGER :: Ninpts1 REAL(r_k), DIMENSION(2) :: extremelon, extremelat REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: fractionlon, fractionlat INTEGER :: dfracdx, dfracdy, fracdx, fracdy CHARACTER(LEN=50) :: fname !!!!!!! Variables ! dimx, dimy: dimension length of the target interpolation ! proj[lon/lat]: longitudes and latitudes of the target interpolation ! Ninpts: number of points to interpolate ! [lon/lat]vs: longitudes and latitudes of the points to interpolate ! mindiff: minimal accepted distance to the target point ! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues ! inpt: whether the point has already been localized (1) or not (0) ! ilonlat: Longitude and Latitude of the input points ! mindiffLl: minimum difference between target and source longitude/latitude (in degrees) fname = 'CoarseInterpolate' Ninpts1 = Ninpts/100 extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /) extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /) PRINT *,' ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen dfracdx = INT(1./percen+1) dfracdy = INT(1./percen+1) fracdx = INT(dimx*percen) fracdy = INT(dimy*percen) PRINT *,' ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon) ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr) IF (ierr /= 0) THEN PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) //": problem allocating 'fractionlon' !!" STOP END IF IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat) ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr) IF (ierr /= 0) THEN PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) //": problem allocating 'fractionlat' !!" STOP END IF DO i=1,dfracdx DO j=1,dfracdy fractionlon(i,j) = projlon(fracdx*(i-1)+1,fracdy*(j-1)+1) fractionlat(i,j) = projlat(fracdx*(i-1)+1,fracdy*(j-1)+1) ! PRINT *,'i,j:',i,', ',j,' frac ij:',fracdx*(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),& ! ', ',fractionlat(i,j) END DO END DO ! PRINT *,' ' // TRIM(fname) // ' fractions of:' ! PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')' ! DO i=1,dfracdx ! PRINT *,fractionlon(i,:) ! END DO ! PRINT *,' lat_______' ! DO i=1,dfracdx ! PRINT *,fractionlat(i,:) ! END DO DO iv=1,Ninpts IF (inpt(iv) == 0) THEN CALL CoarselonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, fractionlon, & fractionlat, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, ilonlat(iv,:), mindiffLl(iv)) IF ((mindiffLl(iv) <= mindiff) .AND. .NOT.(ilonlat(iv,1) >= 0 .AND. ilonlat(iv,1) >= 0)) THEN PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), & ' not relocated !!' PRINT *,' mindiffl:', mindiffLl(iv), ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2) STOP END IF END IF END DO END SUBROUTINE CoarseInterpolate SUBROUTINE CoarseInterpolateExact(projlon, projlat, lonvs, latvs, percen, mindiff, ivar, newvar, & newvarin, newvarinpt, newvarindiff, dimx, dimy, Ninpts) ! Subroutine which finds the closest grid point within a projection throughout a first guest ! and then whole domain approche from percentages of the whole domain USE module_generic IMPLICIT NONE INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: dimx, dimy REAL(r_k), DIMENSION(dimx,dimy), INTENT(in) :: projlon, projlat INTEGER, INTENT(in) :: Ninpts REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: ivar, lonvs, latvs REAL(r_k), INTENT(in) :: mindiff, percen REAL(r_k), DIMENSION(dimx,dimy), INTENT(out) :: newvar INTEGER, DIMENSION(dimx,dimy), INTENT(out) :: newvarin INTEGER, DIMENSION(Ninpts), INTENT(out) :: newvarinpt REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: newvarindiff ! Local INTEGER :: iv,i,j INTEGER :: ierr INTEGER, DIMENSION(2) :: ilonlat REAL(r_k) :: mindiffLl INTEGER :: Ninpts1 REAL(r_k), DIMENSION(2) :: extremelon, extremelat REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: fractionlon, fractionlat INTEGER :: dfracdx, dfracdy, fracdx, fracdy CHARACTER(LEN=50) :: fname !!!!!!! Variables ! dimx, dimy: dimension length of the target interpolation ! proj[lon/lat]: longitudes and latitudes of the target interpolation ! Ninpts: number of points to interpolate ! [lon/lat]vs: longitudes and latitudes of the points to interpolate ! mindiff: minimal accepted distance to the target point ! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues ! ivar: values to localize in the target projection ! newvar: localisation of the [lon/lat]vs point in the target projection ! newvarin: number of point from the input data ! newvarinpt: integer value indicating if the value has been already located (0: no, 1: yes) ! newvarindiff: distance of point from the input data to the closest target point ! ncid: netCDF output file id fname = 'CoarseInterpolateExact' Ninpts1 = Ninpts/100 extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /) extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /) PRINT *,' ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen dfracdx = INT(1./percen+1) dfracdy = INT(1./percen+1) fracdx = INT(dimx*percen) fracdy = INT(dimy*percen) PRINT *,' ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon) ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr) IF (ierr /= 0) THEN PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) //": problem allocating 'fractionlon' !!" STOP END IF IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat) ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr) IF (ierr /= 0) THEN PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) //": problem allocating 'fractionlat' !!" STOP END IF DO i=1,dfracdx DO j=1,dfracdy fractionlon(i,j) = projlon(fracdx*(i-1)+1,fracdy*(j-1)+1) fractionlat(i,j) = projlat(fracdx*(i-1)+1,fracdy*(j-1)+1) ! PRINT *,'i,j:',i,', ',j,' frac ij:',fracdx*(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),& ! ', ',fractionlat(i,j) END DO END DO ! PRINT *,' ' // TRIM(fname) // ' fractions of:' ! PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')' ! DO i=1,dfracdx ! PRINT *,fractionlon(i,:) ! END DO ! PRINT *,' lat_______' ! DO i=1,dfracdx ! PRINT *,fractionlat(i,:) ! END DO DO iv=1,Ninpts IF (newvarinpt(iv) == 0) THEN CALL CoarselonlatFindExact(dimx, dimy, projlon, projlat, extremelon, extremelat, fracdx, fracdy,& fractionlon, fractionlat, iv, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, mindiff, & ilonlat, mindiffLl) IF (mindiffLl >= mindiff) THEN ! percendone(iv,Ninpts,0.5,'done:') IF (ilonlat(1) >= 0 .AND. ilonlat(1) >= 0) THEN newvar(ilonlat(1),ilonlat(2)) = ivar(iv) newvarin(ilonlat(1),ilonlat(2)) = iv newvarinpt(iv) = 1 newvarindiff(iv) = mindiffLl ! PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv), & ! ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:', & ! newvarindiff(iv), ' point:',ilonlat ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), & ' not relocated !!' PRINT *,' mindiffl:', mindiffLl, ' ilon:', ilonlat(1), ' ilat:', ilonlat(2) STOP END IF ! IF (MOD(iv,Ninpts1) == 0) newnc.sync() ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv), & ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ', & mindiff, ' !!' PRINT *,' found minimum difference:', mindiffLl STOP END IF END IF END DO END SUBROUTINE CoarseInterpolateExact SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy, & Ninpts) ! Subroutine which finds the closest grid point within a projection USE module_generic IMPLICIT NONE INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: dimx, dimy REAL(r_k), DIMENSION(dimx,dimy), INTENT(in) :: projlon, projlat INTEGER, INTENT(in) :: Ninpts REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: lonvs, latvs REAL(r_k), INTENT(in) :: mindiff INTEGER, DIMENSION(Ninpts), INTENT(inout) :: inpt REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: diffs INTEGER, DIMENSION(Ninpts,2), INTENT(out) :: ilonlat ! Local INTEGER :: iv REAL(r_k) :: mindiffLl INTEGER :: Ninpts1 REAL(r_k), DIMENSION(dimx,dimy) :: difflonlat REAL(r_k), DIMENSION(2) :: extremelon, extremelat CHARACTER(LEN=50) :: fname !!!!!!! Variables ! dimx, dimy: dimension length of the target interpolation ! proj[lon/lat]: longitudes and latitudes of the target interpolation ! Ninpts: number of points to interpolate ! [lon/lat]vs: longitudes and latitudes of the points to interpolate ! mindiff: minimal accepted distance to the target point ! inpt: whether the point has already been localized ! diffs: distance of point from the input data to the closest target point ! ilonlat: longitude and latitude of the point ! ncid: netCDF output file id fname = 'Interpolate' Ninpts1 = Ninpts/100 extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /) extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /) DO iv=1,Ninpts IF (inpt(iv) <= 0) THEN ! Not using the subroutine, not efficient! ! CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv), & ! ilonlat, mindiffLl) IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!' PRINT *,' given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )' STOP END IF IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!' PRINT *,' given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )' STOP END IF ! Find point difflonlat = SQRT((projlon-lonvs(iv))**2. + (projlat-latvs(iv))**2.) mindiffLl = MINVAL(difflonlat) ilonlat(iv,:) = index2DArrayR(difflonlat, dimx, dimy, mindiffLl) IF (mindiffLl <= mindiff) THEN ! percendone(iv,Ninpts,0.5,'done:') IF (ilonlat(iv,1) >= 0 .AND. ilonlat(iv,2) >= 0) THEN diffs(iv) = mindiffLl inpt(iv) = 1 ! PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv), & ! ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:', & ! newvarindiff(iv), ' point:',ilonlat ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), & ' not relocated !!' PRINT *,' mindiffl:', mindiffLl, ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2) STOP END IF ! IF (MOD(iv,Ninpts1) == 0) newnc.sync() ELSE ! Because doing boxes and Goode is not conitnuos, we should jump this error message PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv), & ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ', & mindiff, ' !!' PRINT *,' found minimum difference:', mindiffLl STOP END IF END IF END DO END SUBROUTINE Interpolate SUBROUTINE Interpolate1DLl(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy, & Ninpts) ! Subroutine which finds the closest grid point within a projection with 1D longitudes and latitudes USE module_generic IMPLICIT NONE INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: dimx, dimy REAL(r_k), DIMENSION(dimx), INTENT(in) :: projlon REAL(r_k), DIMENSION(dimy), INTENT(in) :: projlat INTEGER, INTENT(in) :: Ninpts REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: lonvs, latvs REAL(r_k), INTENT(in) :: mindiff INTEGER, DIMENSION(Ninpts), INTENT(inout) :: inpt REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: diffs INTEGER, DIMENSION(Ninpts,2), INTENT(out) :: ilonlat ! Local INTEGER :: iv REAL(r_k) :: mindifflo, mindiffLa, mindiffLl INTEGER :: Ninpts1 REAL(r_k), DIMENSION(dimx) :: difflon REAL(r_k), DIMENSION(dimy) :: difflat REAL(r_k), DIMENSION(2) :: extremelon, extremelat CHARACTER(LEN=50) :: fname !!!!!!! Variables ! dimx, dimy: dimension length of the target interpolation ! proj[lon/lat]: longitudes and latitudes of the target interpolation ! Ninpts: number of points to interpolate ! [lon/lat]vs: longitudes and latitudes of the points to interpolate ! mindiff: minimal accepted distance to the target point ! inpt: whether the point has already been localized ! diffs: distance of point from the input data to the closest target point ! ilonlat: longitude and latitude of the point ! ncid: netCDF output file id fname = 'Interpolate1DLl' Ninpts1 = Ninpts/100 extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /) extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /) DO iv=1,Ninpts IF (inpt(iv) <= 0) THEN ! Not using the subroutine, not efficient! ! CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv), & ! ilonlat, mindiffLl) IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!' PRINT *,' given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )' STOP END IF IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!' PRINT *,' given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )' STOP END IF ! Find point difflon = SQRT((projlon-lonvs(iv))**2.) difflat = SQRT((projlat-latvs(iv))**2.) mindifflo = MINVAL(difflon) mindiffLa = MINVAL(difflat) mindifflL = SQRT(mindifflo*mindifflo + mindiffLa*mindiffLa) ilonlat(iv,1) = index1DArrayR(difflon, dimx, mindifflo) ilonlat(iv,2) = index1DArrayR(difflat, dimy, mindiffLa) ! PRINT *,' Lluis: iv',iv,' lonvs:', lonvs(iv),' latvs:',latvs(iv) ! PRINT *,' Lluis: mindifflo:', mindifflo,' ilonlat(1):',ilonlat(iv,1) ! PRINT *,' Lluis: mindiffLa:', mindiffLa,' ilonlat(2):',ilonlat(iv,2) IF (mindiffLl <= mindiff) THEN ! percendone(iv,Ninpts,0.5,'done:') IF (ilonlat(iv,1) >= 0 .AND. ilonlat(iv,2) >= 0) THEN diffs(iv) = mindiffLl inpt(iv) = 1 ! PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv), & ! ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:', & ! newvarindiff(iv), ' point:',ilonlat ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), & ' not relocated !!' PRINT *,' mindiffl:', mindiffLl, ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2) STOP END IF ! IF (MOD(iv,Ninpts1) == 0) newnc.sync() ELSE ! Because doing boxes and Goode is not conitnuos, we should jump this error message PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv), & ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ', & mindiff, ' !!' PRINT *,' found minimum difference:', mindiffLl STOP END IF END IF END DO END SUBROUTINE Interpolate1DLl END MODULE module_ForInterpolate