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

! 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 CoarseInterpolate(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 
!   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)                                              :: 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 = '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 (newvarinpt(iv) == 0) THEN
      CALL CoarselonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, fractionlon,        &
        fractionlat, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, 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
! 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 CoarseInterpolate

END MODULE module_ForInterpolate