Changeset 764 in lmdz_wrf


Ignore:
Timestamp:
May 10, 2016, 11:04:59 AM (9 years ago)
Author:
lfita
Message:

Adding `lonlatfrac' for the 'lonlat' kind with fractions, due to the assumption that target grid-point can recieve more than one input value

Location:
trunk/tools
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/module_ForInterpolate.F90

    r742 r764  
    478478END SUBROUTINE lonlatFind
    479479
    480 SUBROUTINE CoarseInterpolate(projlon, projlat, lonvs, latvs, percen, mindiff, ivar, newvar, newvarin, &
    481   newvarinpt, newvarindiff, dimx, dimy, Ninpts)
     480SUBROUTINE CoarseInterpolate(projlon, projlat, lonvs, latvs, percen, mindiff, inpt, ilonlat,          &
     481  mindiffLl, dimx, dimy, Ninpts)
    482482! Subroutine which finds the closest grid point within a projection throughout a first guest
    483483!   approche from percentages of the whole domain
     484
     485  USE module_generic
     486
     487  IMPLICIT NONE
     488
     489  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
     490  INTEGER, INTENT(in)                                    :: dimx, dimy
     491  REAL(r_k), DIMENSION(dimx,dimy), INTENT(in)            :: projlon, projlat
     492  INTEGER, INTENT(in)                                    :: Ninpts
     493  REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: inpt, lonvs, latvs
     494  REAL(r_k), INTENT(in)                                  :: mindiff, percen
     495  INTEGER, DIMENSION(Ninpts,2), INTENT(out)              :: ilonlat
     496  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: mindiffLl
     497
     498! Local
     499  INTEGER                                                :: iv,i,j
     500  INTEGER                                                :: ierr
     501  INTEGER                                                :: Ninpts1
     502  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
     503  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: fractionlon, fractionlat
     504  INTEGER                                                :: dfracdx, dfracdy, fracdx, fracdy
     505  CHARACTER(LEN=50)                                      :: fname
     506
     507!!!!!!! Variables
     508! dimx, dimy: dimension length of the target interpolation
     509! proj[lon/lat]: longitudes and latitudes of the target interpolation
     510! Ninpts: number of points to interpolate
     511! [lon/lat]vs: longitudes and latitudes of the points to interpolate
     512! mindiff: minimal accepted distance to the target point
     513! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues
     514! inpt: whether the point has already been localized (1) or not (0)
     515! ilonlat: Longitude and Latitude of the input points
     516! mindiffLl: minimum difference between target and source longitude/latitude (in degrees)
     517
     518  fname = 'CoarseInterpolate'
     519  Ninpts1 = Ninpts/100
     520
     521  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
     522  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
     523
     524  PRINT *,'  ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen
     525
     526  dfracdx = INT(1./percen+1)
     527  dfracdy = INT(1./percen+1)
     528  fracdx = INT(dimx*percen)
     529  fracdy = INT(dimy*percen)
     530  PRINT *,'  ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy
     531
     532  IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
     533  ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr)
     534  IF (ierr /= 0) THEN
     535    PRINT *,TRIM(ErrWarnMsg('err'))
     536    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlon' !!"
     537    STOP
     538  END IF
     539  IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
     540  ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr)
     541  IF (ierr /= 0) THEN
     542    PRINT *,TRIM(ErrWarnMsg('err'))
     543    PRINT *,'  ' // TRIM(fname) //": problem allocating 'fractionlat' !!"
     544    STOP
     545  END IF
     546
     547  DO i=1,dfracdx
     548    DO j=1,dfracdy
     549      fractionlon(i,j) = projlon(fracdx*(i-1)+1,fracdy*(j-1)+1)
     550      fractionlat(i,j) = projlat(fracdx*(i-1)+1,fracdy*(j-1)+1)
     551!      PRINT *,'i,j:',i,', ',j,' frac ij:',fracdx*(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),&
     552!        ', ',fractionlat(i,j)
     553    END DO
     554  END DO
     555
     556!  PRINT *,'  ' // TRIM(fname) // ' fractions of:'
     557!  PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')'
     558!  DO i=1,dfracdx
     559!    PRINT *,fractionlon(i,:)
     560!  END DO
     561!  PRINT *,' lat_______'
     562!  DO i=1,dfracdx
     563!    PRINT *,fractionlat(i,:)
     564!  END DO
     565
     566  DO iv=1,Ninpts
     567    IF (inpt(iv) == 0) THEN
     568      CALL CoarselonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, fractionlon,        &
     569        fractionlat, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, ilonlat(iv,:), mindiffLl(iv))
     570
     571      IF ((mindiffLl(iv) <= mindiff) .AND. .NOT.(ilonlat(iv,1) >= 0 .AND. ilonlat(iv,1) >= 0)) THEN
     572        PRINT *,TRIM(ErrWarnMsg('err'))
     573        PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
     574          ' not relocated !!'
     575        PRINT *,'    mindiffl:', mindiffLl(iv), ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
     576        STOP
     577      END IF
     578
     579    END IF
     580  END DO
     581
     582END SUBROUTINE CoarseInterpolate
     583
     584SUBROUTINE CoarseInterpolateExact(projlon, projlat, lonvs, latvs, percen, mindiff, ivar, newvar,      &
     585  newvarin, newvarinpt, newvarindiff, dimx, dimy, Ninpts)
     586! Subroutine which finds the closest grid point within a projection throughout a first guest
     587!   and then whole domain approche from percentages of the whole domain
    484588
    485589  USE module_generic
    … …  
    523627! ncid: netCDF output file id
    524628
    525   fname = 'CoarseInterpolate'
     629  fname = 'CoarseInterpolateExact'
    526630  Ninpts1 = Ninpts/100
    527631
    … …  
    573677  DO iv=1,Ninpts
    574678    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
     679      CALL CoarselonlatFindExact(dimx, dimy, projlon, projlat, extremelon, extremelat, fracdx, fracdy,&
     680        fractionlon, fractionlat, iv, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, mindiff,        &
     681        ilonlat, mindiffLl)
     682
     683      IF (mindiffLl >= mindiff) THEN
    579684!        percendone(iv,Ninpts,0.5,'done:')
    580685
    … …  
    596701
    597702!        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 
    610 END SUBROUTINE CoarseInterpolate
    611 
    612 SUBROUTINE 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()
    731703      ELSE
    732704        PRINT *,TRIM(ErrWarnMsg('err'))
    … …  
    841813END SUBROUTINE Interpolate
    842814
     815SUBROUTINE Interpolate1DLl(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy, &
     816  Ninpts)
     817! Subroutine which finds the closest grid point within a projection with 1D longitudes and latitudes
     818
     819  USE module_generic
     820
     821  IMPLICIT NONE
     822
     823  INTEGER, PARAMETER                                     :: r_k = KIND(1.d0)
     824  INTEGER, INTENT(in)                                    :: dimx, dimy
     825  REAL(r_k), DIMENSION(dimx), INTENT(in)                 :: projlon
     826  REAL(r_k), DIMENSION(dimy), INTENT(in)                 :: projlat
     827  INTEGER, INTENT(in)                                    :: Ninpts
     828  REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: lonvs, latvs
     829  REAL(r_k), INTENT(in)                                  :: mindiff
     830  INTEGER, DIMENSION(Ninpts), INTENT(inout)              :: inpt
     831  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: diffs
     832  INTEGER, DIMENSION(Ninpts,2), INTENT(out)              :: ilonlat
     833
     834! Local
     835  INTEGER                                                :: iv
     836  REAL(r_k)                                              :: mindifflo, mindiffLa, mindiffLl
     837  INTEGER                                                :: Ninpts1
     838  REAL(r_k), DIMENSION(dimx)                             :: difflon
     839  REAL(r_k), DIMENSION(dimy)                             :: difflat
     840  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
     841  CHARACTER(LEN=50)                                      :: fname
     842
     843!!!!!!! Variables
     844! dimx, dimy: dimension length of the target interpolation
     845! proj[lon/lat]: longitudes and latitudes of the target interpolation
     846! Ninpts: number of points to interpolate
     847! [lon/lat]vs: longitudes and latitudes of the points to interpolate
     848! mindiff: minimal accepted distance to the target point
     849! inpt: whether the point has already been localized
     850! diffs: distance of point from the input data to the closest target point
     851! ilonlat: longitude and latitude of the point
     852! ncid: netCDF output file id
     853
     854  fname = 'Interpolate1DLl'
     855  Ninpts1 = Ninpts/100
     856
     857  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
     858  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
     859
     860  DO iv=1,Ninpts
     861    IF (inpt(iv) <= 0) THEN
     862! Not using the subroutine, not efficient!
     863!      CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv),     &
     864!        ilonlat, mindiffLl)
     865
     866      IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
     867        PRINT *, TRIM(ErrWarnMsg('err'))
     868        PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
     869        PRINT *,'    given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
     870        STOP
     871      END IF
     872      IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
     873        PRINT *, TRIM(ErrWarnMsg('err'))
     874        PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
     875        PRINT *,'    given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
     876        STOP
     877      END IF
     878
     879! Find point
     880      difflon = SQRT((projlon-lonvs(iv))**2.)
     881      difflat = SQRT((projlat-latvs(iv))**2.)
     882      mindifflo = MINVAL(difflon)
     883      mindiffLa = MINVAL(difflat)
     884      mindifflL = SQRT(mindifflo*mindifflo + mindiffLa*mindiffLa)
     885      ilonlat(iv,1) = index1DArrayR(difflon, dimx, mindifflo)
     886      ilonlat(iv,2) = index1DArrayR(difflat, dimy, mindiffLa)
     887!      PRINT *,'  Lluis: iv',iv,' lonvs:', lonvs(iv),' latvs:',latvs(iv)
     888!      PRINT *,'  Lluis: mindifflo:', mindifflo,' ilonlat(1):',ilonlat(iv,1)
     889!      PRINT *,'  Lluis: mindiffLa:', mindiffLa,' ilonlat(2):',ilonlat(iv,2)
     890
     891
     892      IF (mindiffLl <= mindiff) THEN
     893!        percendone(iv,Ninpts,0.5,'done:')
     894
     895        IF (ilonlat(iv,1) >= 0 .AND. ilonlat(iv,2) >= 0) THEN
     896          diffs(iv) = mindiffLl
     897          inpt(iv) = 1
     898!          PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv),        &
     899!            ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:',           &
     900!            newvarindiff(iv), ' point:',ilonlat   
     901        ELSE
     902          PRINT *,TRIM(ErrWarnMsg('err'))
     903          PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
     904            ' not relocated !!'
     905          PRINT *,'    mindiffl:', mindiffLl, ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
     906          STOP
     907        END IF
     908
     909!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
     910      ELSE
     911! Because doing boxes and Goode is not conitnuos, we should jump this error message
     912        PRINT *,TRIM(ErrWarnMsg('err'))
     913        PRINT *,'  ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv),   &
     914          ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ',       &
     915          mindiff, ' !!'
     916        PRINT *,'    found minimum difference:', mindiffLl
     917        STOP
     918      END IF
     919    END IF
     920  END DO
     921
     922END SUBROUTINE Interpolate1DLl
     923
    843924END MODULE module_ForInterpolate

  • trunk/tools/module_generic.F90

    r715 r764  
    33
    44  CONTAINS
     5
     6  INTEGER FUNCTION Index1DArrayR(array1D, d1, val)
     7! Function to provide the first index of a given value inside a 1D real array
     8
     9    IMPLICIT NONE
     10
     11    INTEGER, PARAMETER                                   :: r_k = KIND(1.d0)
     12    INTEGER, INTENT(in)                                  :: d1
     13    REAL(r_k), INTENT(in)                                :: val
     14    REAL(r_k), DIMENSION(d1), INTENT(in)                 :: array1D
     15
     16! Local
     17    INTEGER                                              :: i
     18    CHARACTER(LEN=50)                                    :: fname
     19
     20    fname = 'Index1DArrayR'
     21
     22    Index1DArrayR = -1
     23
     24    DO i=1,d1
     25      IF (array1d(i) == val) THEN
     26        Index1DArrayR = i
     27        EXIT
     28      END IF
     29    END DO
     30
     31  END FUNCTION Index1DArrayR
    532
    633  FUNCTION Index2DArrayR(array2D, d1, d2, val)

  • trunk/tools/nc_var_tools.py

    r763 r764  
    44674467        for ich in range(Lstr):
    44684468            chrv = stringv[ich:ich+1]
    4469             if ord(chrv) == 32: chrv = unichr(32)
    44704469            charvals[ich] = stringv[ich:ich+1]
    44714470        varo[iv,:] = charvals
    … …  
    1311313112    else:
    1311413113        print errormsg
    13115         print '  ' + fnamne + ": projection '" + project + "' not ready !!"
     13114        print '  ' + fname + ": projection '" + project + "' not ready !!"
    1311613115        quit(-1)
    1311713116
    … …  
    1363613635      foudre: f2py -m module_ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c module_ForInterpolate.F90 module_generic.F90 >& run_f2py.log
    1363713636      ciclad: f2py --f90flags="-fPIC" -L/opt/canopy-1.3.0/Canopy_64bit/System/lib/ -L/usr/lib64/ -L/opt/canopy-1.3.0/Canopy_64bit/System/lib/ -m module_ForInterpolate -c module_ForInterpolate.F90 module_generic.F90
    13638       values= [lonmame],[latname],[fillVal],[resolution],[kind],[lonlatProjfile],[ipoint],[fracd],[fracs]
     13637      values= [lonmame],[latname],[fillVal],[resolution],[kind],[lonlatProjfile],[ipoint],[fracd],[fracs],[mindiff]
    1363913638        [lonname]: name of the longitude variable
    1364013639        [latname]: name of the latitude variable
    … …  
    1364813647        [kind]: kind of target projection
    1364913648         'lonlat': standard lon/lat projection tacking number of points from [dx/dy]res at the Equator
     13649         'lonlatfrac': standard lon/lat projection tacking number of points from [dx/dy]res at the Equator with
     13650            mnultiple values at the target interpolation grid [kind]@[Ntypes]
     13651              [Ntypes]: Number of types at each grid-point
    1365013652         'lonlat_dxdyFix': projection with a fixed grid distance
    1365113653         'Goode': Goode projection
    … …  
    1365413656        [fracd]: Percentage of the fractions within perform the first guess search
    1365513657        [fracs]: Number of grid points to perform the syncronization with the file and the computed values
     13658        [mindiff]: Authorized minium distance between input and final lon,lat point
    1365613659      filen= name of the netCDF file
    1365713660      varn= name of the variable
    … …  
    1366913672
    1367013673    arguments = '[lonmame],[latname],[fillVal],[resolution],[kind],' +               \
    13671       '[lonlatProjfile],[ipoint],[fracd],[fracs]'
    13672     check_arguments(fname, values, arguments, ',')
     13674      '[lonlatProjfile],[ipoint],[fracd],[fracs],[mindiff]'
     13675    gen.check_arguments(fname, values, arguments, ',')
    1367313676
    1367413677    ofile = 'EntireGlobalMap.nc'
    … …  
    1368513688    fval0 = values.split(',')[2]
    1368613689    resolution = np.float(values.split(',')[3])
    13687     kind = values.split(',')[4]
     13690    kind0 = values.split(',')[4]
    1368813691    Projfile = values.split(',')[5]
    1368913692    ipoint = int(values.split(',')[6])
    13690     fracd = np.float(values.split(',')[7])
     13693    fracd = np.float64(values.split(',')[7])
    1369113694    fracs = int(values.split(',')[8])
     13695    mindiff = np.float64(values.split(',')[9])
    1369213696
    1369313697    if Projfile == 'None':
    … …  
    1369513699    else:
    1369613700        lonlatProjfile = Projfile
     13701
     13702    if kind0.find('@') != -1:
     13703        kind = kind0.split('@')[0]
     13704        Ntypes = np.int(kind0.split('@')[1])
     13705    else:
     13706        kind = kind0
     13707
     13708    Lkind = len(kind)
     13709    if kind[Lkind-4:Lkind] == 'frac' and kind0.find('@') == -1:
     13710        print errormsg
     13711        print '  ' + fname + ": there is a fractional kind of interpolation '" +     \
     13712          kind + "' but no number of types has been privided!!"
     13713        print "    when [kind] is provided, should be '[kind]@[Ntypes]'"
     13714        quit()
    1369713715
    1369813716    if not onc.variables.has_key(lonname):
    … …  
    1376713785        quit(-1)
    1376813786
    13769     if lonlatProjfile is None:
    13770         lonlatProjfile = kind + '.nc'
     13787    if lonlatProjfile is None or not os.path.isfile(lonlatProjfile):
     13788        if lonlatProjfile is None: lonlatProjfile = kind + '.nc'
    1377113789        print warnmsg
    1377213790        print '  ' + fname + ": no reference map !!"
    … …  
    1377413792        print '    creating it via: lonlatProj(', resolution, ',', resolution,       \
    1377513793          ', ' + kind + ', True)'
    13776         lonmap, latmap = lonlatProj(resolution, resolution, kind, True)
    13777         sub.call(['mv','lonlatProj.nc',lonlatProjfile])
     13794        if kind[Lkind-4:Lkind] == 'frac':
     13795            lonmap, latmap = lonlatProj(resolution, resolution, kind[0:Lkind-4], True)
     13796        else:
     13797            lonmap, latmap = lonlatProj(resolution, resolution, kind, True)
     13798        sub.call(['mv','lonlatProj.nc', lonlatProjfile])
    1377813799
    1377913800    oproj = NetCDFFile(lonlatProjfile, 'r')
    … …  
    1380213823        Nmaplat = olat.shape[0]
    1380313824
     13825    elif kind == 'lonlatfrac':
     13826        olon = oproj.variables['lon']
     13827        olat = oproj.variables['lat']
     13828        lonmap = olon[:]
     13829        latmap = olat[:]
     13830        projlon = olon[:]
     13831        projlat = olat[:]
     13832        omapvals = oproj.variables['points']
     13833        Ntotvals = omapvals.shape[0] * omapvals.shape[1]
     13834        Nmaplon = olon.shape[0]
     13835        Nmaplat = olat.shape[0]
     13836
    1380413837    elif kind == 'lonlat_dxdyFix':
    1380513838        oprojlon = oproj.variables['loncirc']
    … …  
    1383013863    varunits = ovar.units
    1383113864
    13832     fval = fillvalue_kind(type(ovar[0]), fval0)
    13833 
     13865    fval = gen.fillvalue_kind(type(ovar[0]), fval0)
     13866   
    1383413867# Final map creation
    1383513868##
    … …  
    1386713900
    1386813901        newvarindiff = newnc.createVariable('locindiff','f4',('inpts'))
    13869         basicvardef(newvarindiff, 'locindiff', 'distance between input point and its final location','degree')
     13902        basicvardef(newvarindiff, 'locindiff', 'distance between input point and ' + \
     13903          'its final location','degree')
    1387013904
    1387113905# map variable
    1387213906        Lmapvalsshape = len(omapvals.shape)
    1387313907        if Lmapvalsshape == 1:
    13874             newvar = newnc.createVariable(varn, nctype(vartype), ('pts'), fill_value=fval)
     13908            newvar = newnc.createVariable(varn, nctype(vartype), ('pts'),            \
     13909              fill_value=fval)
    1387513910        else:
    13876             newvar = newnc.createVariable(varn, nctype(vartype), ('lat','lon'), fill_value=fval)
     13911            if kind[Lkind-4:Lkind] == 'frac':
     13912                newdim = newnc.createDimension('Ntypes', Ntypes)
     13913                newvar = newnc.createVariable(varn, 'f4', ('Ntypes','lat','lon'),    \
     13914                  fill_value=gen.fillvalue_kind(type(1.), fval0))
     13915            else:
     13916                newvar = newnc.createVariable(varn, nctype(vartype), ('lat','lon'),  \
     13917                  fill_value=fval)
    1387713918
    1387813919        basicvardef(newvar, varn, varlongname, varunits)
    … …  
    1390913950
    1391013951# Reducing the searching points
    13911     newvarinvals = newvarinpt[:]
     13952    Nsearchpt = 499999
     13953    Ninpts=np.min([Nsearchpt, Ninpts-ipoint])
     13954    newvarinvals = newvarinpt[ipoint:ipoint+Ninpts]
    1391213955    maskpt = np.where(newvarinvals.mask == True, False, True)
    13913     points = np.arange(Ninpts)
     13956    points = np.arange(ipoint,ipoint+Ninpts,1)
    1391413957    mapoints = ma.array(points, mask=maskpt)
    1391513958    ptsf = mapoints.compressed()
     13959
     13960    Nptsf = len(ptsf)
     13961    print Ninpts,'Npoints to find:', len(ptsf), ptsf[0:10], newvarindiff[ptsf[0:10]]
     13962    print '  ' + fname + ': from:', Ninpts,'re-locating:',Nptsf,'points starting at',\
     13963      ipoint,'...'
    1391613964
    1391713965    Nptsf = len(ptsf)
    … …  
    1394313991
    1394413992    elif kind == 'lonlat':
    13945         for iv in range(Ntotvals):
    13946             difflat = np.abs(projlat - latvs[iv])
    13947             mindiffL = np.min(difflat)
    13948             ilat = index_mat(difflat, mindiffL)
    13949             difflon = np.abs(projlon - lonvs[iv])
    13950             mindiffl = difflon.min()
    13951             ilon = index_mat(difflon, mindiffl)
    13952             percendone(iv,Ntotvals,0.5,'done:')
    13953             if mindiffl > mindiff or mindiffL > mindiff:
    13954                 print errormsg
    13955                 print '  ' + fname + ': for point #', iv,'lon,lat in incomplet map:',   \
    13956                   lonvs[iv], ',', latvs[iv], 'there is not a set of lon,lat in the ' +  \
    13957                   'completed map closer than: ', mindiff, '!!'
    13958                 print '    minimum difflon:', np.min(difflon), 'difflat:', np.min(difflat)
    13959                 quit()
    13960 
    13961             if ilon >= 0 and ilat >= 0:
    13962                 newvar[ilat,ilon] = ovar[iv]         
    13963                 newnc.sync()
    13964             else:
    13965                 print errormsg
    13966                 print '  ' + fname + ': point iv:',iv,'at',lonvs[iv],',',latvs[iv],' not relocated !!'
    13967                 print '    mindiffl:',mindiffl,'mindiffL:',mindiffL,'ilon:',ilon,'ilat:',ilat
    13968                 quit(-1)
    13969 
    13970         newnc.sync()
     13993        for ir in range(ipoint,ipoint+Ninpts,fracs):
     13994            iri = ir
     13995            ire = ir + fracs
     13996            print iri,',',ire
     13997#            percendone(iri,Ninpts,0.5,'done:')
     13998            lonvss = lonvs[iri:ire+1].astype('float64')
     13999            latvss = latvs[iri:ire+1].astype('float64')
     14000            inptss = newvarinpt[iri:ire+1].astype('int32')
     14001
     14002            idiff, ilonlatv = fin.module_forinterpolate.interpolate1dll(projlon,     \
     14003              projlat, lonvss, latvss, np.float64(mindiff), inptss)
     14004            print 'Lluis shapes: newvar', newvar.shape,'ilonlatv:',ilonlatv.shape,'range:',np.min([len(idiff),fracs])
     14005            print 'Lluis shapes: ovar:',ovar.shape
     14006
     14007            for i in range(np.min([len(idiff),fracs])):
     14008#                print 'Lluis:',ilonlatv[i,1],',',ilonlatv[i,0]
     14009                newvar[ilonlatv[i,1],ilonlatv[i,0]] = ovar[iri + i]
     14010                newvarin[ilonlatv[i,1],ilonlatv[i,0]] = iri + i
     14011                newvarinpt[iri+i] = 1
     14012                newvarindiff[iri+i] = idiff[i]
     14013
     14014            newnc.sync()
     14015#            newnc.close()
     14016#            quit()
     14017
     14018    elif kind == 'lonlatfrac':
     14019        for ir in range(ipoint,ipoint+Ninpts,fracs):
     14020            iri = ir
     14021            ire = ir + fracs
     14022            print iri,',',ire
     14023#            percendone(iri,Ninpts,0.5,'done:')
     14024            lonvss = lonvs[iri:ire+1].astype('float64')
     14025            latvss = latvs[iri:ire+1].astype('float64')
     14026            inptss = newvarinpt[iri:ire+1].astype('int32')
     14027
     14028            idiff, ilonlatv = fin.module_forinterpolate.interpolate1dll(projlon,     \
     14029              projlat, lonvss, latvss, np.float64(mindiff), inptss)
     14030
     14031            for i in range(np.min([len(idiff),fracs])):
     14032                print i,'Lluis:',ovar[iri+i], ilonlatv[i,1],',',ilonlatv[i,0],newvar[ovar[iri+i],ilonlatv[i,1],ilonlatv[i,0]], type(newvar[ovar[iri+i],ilonlatv[i,1],ilonlatv[i,0]])
     14033                if newvar[ovar[iri+i],ilonlatv[i,1],ilonlatv[i,0]].mask:
     14034                    newvar[ovar[iri+i],ilonlatv[i,1],ilonlatv[i,0]].mask = False
     14035                    newvar[ovar[iri+i],ilonlatv[i,1],ilonlatv[i,0]] = 1.
     14036                else:
     14037                    newvar[ovar[iri+i],ilonlatv[i,1],ilonlatv[i,0]] =                \
     14038                      newvar[ovar[iri+i],ilonlatv[i,1],ilonlatv[i,0]] + 1
     14039                newvarin[ilonlatv[i,1],ilonlatv[i,0]] = iri + i
     14040                newvarinpt[iri+i] = 1
     14041                newvarindiff[iri+i] = idiff[i]
     14042
     14043            newnc.sync()
     14044#            newnc.close()
     14045#            quit()
    1397114046
    1397214047    elif kind == 'lonlat_dxdyFix':
Note: See TracChangeset for help on using the changeset viewer.