Changeset 735 in lmdz_wrf


Ignore:
Timestamp:
Apr 28, 2016, 2:43:28 PM (9 years ago)
Author:
lfita
Message:

Seems to working version (from scratch) as `Partialmap_EntiremapForExact' directly without coarse matrices

Location:
trunk/tools
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/module_ForInterpolate.F90

    r733 r735  
    197197END SUBROUTINE CoarselonlatFind
    198198
     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
    199430SUBROUTINE lonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, lonv, latv, ilonlat, mindiffLl)
    200431! Function to search a given value from a coarser version of the data
    … …  
    345576        fractionlat, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, ilonlat, mindiffLl)
    346577
    347       PRINT *,'  Lluis: iv',iv,', ', mindiffLl,'<= ',mindiff,' ilonlat:',ilonlat
    348578      IF (mindiffLl <= mindiff) THEN
    349579!        percendone(iv,Ninpts,0.5,'done:')
    … …  
    380610END SUBROUTINE CoarseInterpolate
    381611
    382 SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, ivar, newvar, newvarin, newvarinpt,   &
    383   newvarindiff, dimx, dimy, Ninpts)
    384 ! Subroutine which finds the closest grid point within a projection
     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
    385616
    386617  USE module_generic
    … …  
    393624  INTEGER, INTENT(in)                                    :: Ninpts
    394625  REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: ivar, lonvs, latvs
    395   REAL(r_k), INTENT(in)                                  :: mindiff
     626  REAL(r_k), INTENT(in)                                  :: mindiff, percen
    396627  REAL(r_k), DIMENSION(dimx,dimy), INTENT(out)           :: newvar
    397628  INTEGER, DIMENSION(dimx,dimy), INTENT(out)             :: newvarin
    … …  
    407638  INTEGER                                                :: Ninpts1
    408639  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
     640  REAL(r_k), ALLOCATABLE, DIMENSION(:,:)                 :: fractionlon, fractionlat
     641  INTEGER                                                :: dfracdx, dfracdy, fracdx, fracdy
    409642  CHARACTER(LEN=50)                                      :: fname
    410643
    … …  
    415648! [lon/lat]vs: longitudes and latitudes of the points to interpolate
    416649! 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
    417651! ivar: values to localize in the target projection
    418652! newvar: localisation of the [lon/lat]vs point in the target projection
    … …  
    422656! ncid: netCDF output file id
    423657
    424   fname = 'Interpolate'
     658  fname = 'CoarseInterpolateExact'
    425659  Ninpts1 = Ninpts/100
    426660
    … …  
    428662  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
    429663
     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
    430706  DO iv=1,Ninpts
    431707    IF (newvarinpt(iv) == 0) THEN
    432 ! Not using the subroutine, not efficient!
    433 !      CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv),     &
    434 !        ilonlat, mindiffLl)
    435 
    436       IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
    437         PRINT *, TRIM(ErrWarnMsg('err'))
    438         PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
    439         PRINT *,'    given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
    440         STOP
    441       END IF
    442       IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
    443         PRINT *, TRIM(ErrWarnMsg('err'))
    444         PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
    445         PRINT *,'    given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
    446         STOP
    447       END IF
    448 
    449 ! Find point
    450       difflonlat = SQRT((projlon-lonvs(iv))**2. + (projlat-latvs(iv))**2.)
    451       mindiffLl = MINVAL(difflonlat)
    452       ilonlat = index2DArrayR(difflonlat, dimx, dimy, mindiffLl)
    453 
    454       IF (mindiffLl <= mindiff) 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
    455713!        percendone(iv,Ninpts,0.5,'done:')
    456714
    … …  
    472730
    473731!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
    474 !      ELSE
     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, newvar, newvarin, newvarinpt,   &
     746  newvarindiff, 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  REAL(r_k), DIMENSION(dimx,dimy), INTENT(out)           :: newvar
     760  INTEGER, DIMENSION(dimx,dimy), INTENT(out)             :: newvarin
     761  INTEGER, DIMENSION(Ninpts), INTENT(out)                :: newvarinpt
     762  REAL(r_k), DIMENSION(Ninpts), INTENT(out)              :: newvarindiff
     763
     764! Local
     765  INTEGER                                                :: iv,i,j
     766  INTEGER                                                :: ierr
     767  INTEGER, DIMENSION(2)                                  :: ilonlat
     768  REAL(r_k), DIMENSION(dimx,dimy)                        :: difflonlat
     769  REAL(r_k)                                              :: mindiffLl
     770  INTEGER                                                :: Ninpts1
     771  REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat
     772  CHARACTER(LEN=50)                                      :: fname
     773
     774!!!!!!! Variables
     775! dimx, dimy: dimension length of the target interpolation
     776! proj[lon/lat]: longitudes and latitudes of the target interpolation
     777! Ninpts: number of points to interpolate
     778! [lon/lat]vs: longitudes and latitudes of the points to interpolate
     779! mindiff: minimal accepted distance to the target point
     780! ivar: values to localize in the target projection
     781! newvar: localisation of the [lon/lat]vs point in the target projection
     782! newvarin: number of point from the input data
     783! newvarinpt: integer value indicating if the value has been already located (0: no, 1: yes)
     784! newvarindiff: distance of point from the input data to the closest target point
     785! ncid: netCDF output file id
     786
     787  fname = 'Interpolate'
     788  Ninpts1 = Ninpts/100
     789
     790  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
     791  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
     792
     793  DO iv=1,Ninpts
     794    IF (newvarinpt(iv) == 0) THEN
     795! Not using the subroutine, not efficient!
     796!      CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv),     &
     797!        ilonlat, mindiffLl)
     798
     799      IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
     800        PRINT *, TRIM(ErrWarnMsg('err'))
     801        PRINT *,'  ' // TRIM(fname) // ': longitude outside data range!!'
     802        PRINT *,'    given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
     803        STOP
     804      END IF
     805      IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
     806        PRINT *, TRIM(ErrWarnMsg('err'))
     807        PRINT *,'  ' // TRIM(fname) // ': latitude outside data range!!'
     808        PRINT *,'    given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
     809        STOP
     810      END IF
     811
     812! Find point
     813      difflonlat = SQRT((projlon-lonvs(iv))**2. + (projlat-latvs(iv))**2.)
     814      mindiffLl = MINVAL(difflonlat)
     815      ilonlat = index2DArrayR(difflonlat, dimx, dimy, mindiffLl)
     816
     817      IF (mindiffLl <= mindiff) THEN
     818!        percendone(iv,Ninpts,0.5,'done:')
     819
     820        IF (ilonlat(1) >= 0 .AND. ilonlat(1) >= 0) THEN
     821          newvar(ilonlat(1),ilonlat(2)) = ivar(iv)
     822          newvarin(ilonlat(1),ilonlat(2)) = iv
     823          newvarinpt(iv) = 1
     824          newvarindiff(iv) = mindiffLl
     825!          PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv),        &
     826!            ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:',           &
     827!            newvarindiff(iv), ' point:',ilonlat   
     828        ELSE
     829          PRINT *,TRIM(ErrWarnMsg('err'))
     830          PRINT *,'  ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv),        &
     831            ' not relocated !!'
     832          PRINT *,'    mindiffl:', mindiffLl, ' ilon:', ilonlat(1), ' ilat:', ilonlat(2)
     833          STOP
     834        END IF
     835
     836!        IF (MOD(iv,Ninpts1) == 0) newnc.sync()
     837      ELSE
    475838! Because doing boxes and Goode is not conitnuos, we should jump this error message
    476 !        PRINT *,TRIM(ErrWarnMsg('err'))
    477 !        PRINT *,'  ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv),   &
    478 !          ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ',       &
    479 !          mindiff, ' !!'
    480 !        PRINT *,'    found minimum difference:', mindiffLl
    481 !        STOP
    482       END IF
     839        PRINT *,TRIM(ErrWarnMsg('err'))
     840        PRINT *,'  ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv),   &
     841          ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ',       &
     842          mindiff, ' !!'
     843        PRINT *,'    found minimum difference:', mindiffLl
     844        STOP
     845      END IF
     846    ELSE
     847      PRINT *,TRIM(ErrWarnMsg('err'))
     848      PRINT *,'  ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv),     &
     849        ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ',         &
     850        mindiff, ' !!'
     851      PRINT *,'    found minimum difference:', mindiffLl
     852      STOP
    483853    END IF
    484854  END DO

  • trunk/tools/nc_var.py

    r716 r735  
    3535  'maskvar',                                                                         \
    3636  'ncreplace', 'ncstepdiff', 'netcdf_concatenation', 'netcdf_fold_concatenation',    \
    37   'Partialmap_Entiremap', 'Partialmap_EntiremapFor', 'remapnn',                      \
     37  'Partialmap_Entiremap', 'Partialmap_EntiremapFor', 'Partialmap_EntiremapForExact', \
     38  'remapnn',                                                                         \
    3839  'seasmean', 'sellonlatbox', 'sellonlatlevbox', 'selvar', 'setvar_asciivalues',     \
    3940  'sorttimesmat', 'spacemean', 'SpatialWeightedMean', 'statcompare_files', 'submns', \
    … …  
    212213elif oper == 'Partialmap_EntiremapFor':
    213214    ncvar.Partialmap_EntiremapFor(opts.values, opts.ncfile, opts.varname)
     215elif oper == 'Partialmap_EntiremapForExact':
     216    ncvar.Partialmap_EntiremapForExact(opts.values, opts.ncfile, opts.varname)
    214217elif oper == 'seasmean':
    215218    ncvar.seasmean(timename, opts.ncfile, opts.varname)

  • trunk/tools/nc_var_tools.py

    r734 r735  
    1984019840    return
    1984119841
     19842def Partialmap_EntiremapForExact(values, filen, varn):
     19843    """ Function to transform from a partial global map (e.g.: only land points) to an entire one using Fortran code with exact location
     19844      Coincidence of points is done throughout a first guess from fractions of the total domain of search
     19845      Using fortran codes: module_ForInterpolate.F90, module_generic.F90
     19846      foudre: f2py -m module_ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c module_ForInterpolate.F90 module_generic.F90 >& run_f2py.log
     19847      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
     19848      values= [lonmame],[latname],[fillVal],[resolution],[kind],[lonlatProjfile],[ipoint],[fracd],[fracs],[mindiff]
     19849        [lonname]: name of the longitude variable
     19850        [latname]: name of the latitude variable
     19851        [fillVal]: value for '_FillValue', 'std' for the standard value
     19852           Float = 1.e20
     19853           Character = '-'
     19854           Integer = -99999
     19855           Float64 = 1.e20
     19856           Integer32 = -99999
     19857        [resolution]: resolution of the map
     19858        [kind]: kind of target projection
     19859         'lonlat': standard lon/lat projection tacking number of points from [dx/dy]res at the Equator
     19860         'lonlat_dxdyFix': projection with a fixed grid distance
     19861         'Goode': Goode projection
     19862        [lonlatProjfile]: file with the lon,lat of the desired projection. 'None' to be computed and written on fly
     19863        [ipoint]: initial point to use for the interpolation (0, the first)
     19864        [fracd]: Percentage of the fractions within perform the first guess search
     19865        [fracs]: Number of grid points to perform the syncronization with the file and the computed values
     19866        [mindiff]: Authorized minium distance between input and final lon,lat point
     19867      filen= name of the netCDF file
     19868      varn= name of the variable
     19869    """
     19870    import module_ForInterpolate as fin
     19871    import numpy.ma as ma
     19872    import subprocess as sub
     19873
     19874    fname = 'Partialmap_EntiremapFortran'
     19875
     19876    if values == 'h':
     19877        print fname + '_____________________________________________________________'
     19878        print Partialmap_Entiremap.__doc__
     19879        quit()
     19880
     19881    arguments = '[lonmame],[latname],[fillVal],[resolution],[kind],' +               \
     19882      '[lonlatProjfile],[ipoint],[fracd],[fracs],[mindiff]'
     19883    check_arguments(fname, values, arguments, ',')
     19884
     19885    ofile = 'EntireGlobalMap.nc'
     19886
     19887    onc = NetCDFFile(filen, 'r')
     19888    if not onc.variables.has_key(varn):
     19889        print errormsg
     19890        print '  ' + fname + ": file '" + filen + "' does not have " +       \
     19891          "variable '" + varn + "' !!"
     19892        quit(-1)
     19893
     19894    lonname = values.split(',')[0]
     19895    latname = values.split(',')[1]
     19896    fval0 = values.split(',')[2]
     19897    resolution = np.float(values.split(',')[3])
     19898    kind = values.split(',')[4]
     19899    Projfile = values.split(',')[5]
     19900    ipoint = int(values.split(',')[6])
     19901    fracd = np.float64(values.split(',')[7])
     19902    fracs = int(values.split(',')[8])
     19903    mindiff = np.float64(values.split(',')[9])
     19904
     19905    if Projfile == 'None':
     19906        lonlatProjfile = None
     19907    else:
     19908        lonlatProjfile = Projfile
     19909
     19910    if not onc.variables.has_key(lonname):
     19911        print errormsg
     19912        print '  ' + fname + ": file '" + filen + "' does not have " +       \
     19913          "longitude '" + lonname + "' !!"
     19914        quit(-1)
     19915    if not onc.variables.has_key(latname):
     19916        print errormsg
     19917        print '  ' + fname + ": file '" + filen + "' does not have " +       \
     19918          "latitude '" + latname + "' !!"
     19919        quit(-1)
     19920
     19921    olon = onc.variables[lonname]
     19922    olat = onc.variables[latname]
     19923
     19924    lonvs = olon[:]
     19925    latvs = olat[:]
     19926
     19927    Ninpts = len(lonvs)
     19928
     19929    nlat = np.min(latvs)
     19930    nlon = np.min(lonvs)
     19931
     19932    minlat = np.min(np.abs(latvs))
     19933    minlon = np.min(np.abs(lonvs))
     19934
     19935    print '  ' + fname + ': Closest latitude to Equator:', minlat
     19936    print '  ' + fname + ': Closest longitude to Greenwich Meridian:', minlon
     19937    print '  ' + fname + ': Minium distance for point coincidence:', mindiff
     19938
     19939# Existing longitudes along/closest to the Equator
     19940    eqlons = []
     19941    for i in range(len(latvs)):
     19942        if latvs[i] == minlat: eqlons.append(lonvs[i])
     19943
     19944# Existing latitudes along/closest to the Greenwich Meridian
     19945    grlats = []
     19946    for i in range(len(lonvs)):
     19947        if lonvs[i] == minlon: grlats.append(latvs[i])
     19948
     19949    sorteqlons = np.sort(eqlons)
     19950    sortgrlats = np.sort(grlats)
     19951 
     19952    Neqlons = len(sorteqlons)
     19953    Ngrlats = len(sortgrlats)
     19954
     19955    if Neqlons > 1:
     19956        diffeqlons = sorteqlons[1:Neqlons-1] - sorteqlons[0:Neqlons-2]     
     19957        mindeqlon = np.min(diffeqlons)
     19958        print 'N sorteqlons:',Neqlons,'min deqlon:', mindeqlon
     19959    else:
     19960        mindeqlon = None
     19961
     19962    if Ngrlats > 1:
     19963        mindgrlat = np.min(diffgrlats)
     19964        diffgrlats = sortgrlats[1:Ngrlats-1] - sortgrlats[0:Ngrlats-2]
     19965        print 'N sortgrlats:',Ngrlats,'min dgrlat:', mindgrlat
     19966        latmap = np.range(0,360.,360./mindeqlon)
     19967    else:
     19968        mindgrlat = None
     19969
     19970# Fixing in case it has not worked
     19971    if mindeqlon is not None and mindgrlat is None: mindgrlat = mindeqlon
     19972    if mindeqlon is None and mindgrlat is not None: mindeqlon = mindgrlat
     19973    if mindeqlon is None and mindgrlat is None:
     19974        print errormsg
     19975        print fname + ': Not enough values along the Equator and Greenwich!!'
     19976        quit(-1)
     19977
     19978    if lonlatProjfile is None:
     19979        lonlatProjfile = kind + '.nc'
     19980        print warnmsg
     19981        print '  ' + fname + ": no reference map !!"
     19982        print "    creation of '" + lonlatProjfile + "'"
     19983        print '    creating it via: lonlatProj(', resolution, ',', resolution,       \
     19984          ', ' + kind + ', True)'
     19985        lonmap, latmap = lonlatProj(resolution, resolution, kind, True)
     19986        sub.call(['mv','lonlatProj.nc',lonlatProjfile])
     19987
     19988    oproj = NetCDFFile(lonlatProjfile, 'r')
     19989    if kind == 'Goode':
     19990        olon = oproj.variables['lon']
     19991        olat = oproj.variables['lat']
     19992        lonmap = olon[:]
     19993        latmap = olat[:]
     19994        projlon = olon[:].astype('float64')
     19995        projlat = olat[:].astype('float64')
     19996        omapvals = oproj.variables['ptid']
     19997        Ntotvals = omapvals.shape[0] * omapvals.shape[1]
     19998        Nmaplon = olon.shape[1]
     19999        Nmaplat = olat.shape[0]
     20000
     20001    elif kind == 'lonlat':
     20002        olon = oproj.variables['lon']
     20003        olat = oproj.variables['lat']
     20004        lonmap = olon[:]
     20005        latmap = olat[:]
     20006        projlon = olon[:]
     20007        projlat = olat[:]
     20008        omapvals = oproj.variables['points']
     20009        Ntotvals = omapvals.shape[0] * omapvals.shape[1]
     20010        Nmaplon = olon.shape[0]
     20011        Nmaplat = olat.shape[0]
     20012
     20013    elif kind == 'lonlat_dxdyFix':
     20014        oprojlon = oproj.variables['loncirc']
     20015        oprojlat = oproj.variables['latcirc']
     20016        oNprojlon = oproj.variables['Nloncirc']
     20017        projlat = oprojlat[:]
     20018        Nprojlon = oNprojlon[:]
     20019        Ntotvals = len(oproj.dimensions['Npts'])
     20020        olon = oproj.variables['lon']
     20021        olat = oproj.variables['lat']
     20022        lonmap = olon[:]
     20023        latmap = olat[:]
     20024        omapvals = oproj.variables['points']
     20025
     20026#    Ntotvals = len(lonmap)
     20027        Nmaplon = Ntotvals
     20028        Nmaplat = Ntotvals
     20029    else:
     20030        print errormsg
     20031        print '  ' + fname + ": projection kind '" + kind + "' not ready!!"
     20032        quit(-1)
     20033
     20034# Matrix map creation
     20035##
     20036    ovar = onc.variables[varn]
     20037    vartype = type(ovar[0])
     20038    varlongname = ovar.long_name
     20039    varunits = ovar.units
     20040
     20041    fval = fillvalue_kind(type(ovar[0]), fval0)
     20042
     20043# Final map creation
     20044##
     20045    if not os.path.isfile(ofile):
     20046        print '  ' + fname + "File '" + ofile + "' does not exist !!"
     20047        print '    cration of output file'
     20048        newnc = NetCDFFile(ofile, 'w')
     20049# dimensions
     20050        newdim = newnc.createDimension('lon',Nmaplon)
     20051        newdim = newnc.createDimension('lat',Nmaplat)
     20052        newdim = newnc.createDimension('inpts', Ninpts)
     20053        newdim = newnc.createDimension('pts', Ntotvals)
     20054
     20055# Variables
     20056        if kind == 'Goode':
     20057            newvar = newnc.createVariable('lon','f8',('lat', 'lon'))
     20058        else:
     20059            newvar = newnc.createVariable('lon','f8',('lon'))
     20060        basicvardef(newvar, 'lon', 'Longitudes','degrees_East')
     20061        newvar[:] = lonmap
     20062        newvar.setncattr('axis', 'X')
     20063        newvar.setncattr('_CoordinateAxisType', 'Lon')   
     20064
     20065        if kind == 'Goode':
     20066            newvar = newnc.createVariable('lat','f8',('lat','lon'))
     20067        else:
     20068            newvar = newnc.createVariable('lat','f8',('lat'))
     20069        basicvardef(newvar, 'lat', 'Latitudes','degrees_North')
     20070        newvar[:] = latmap
     20071        newvar.setncattr('axis', 'Y')
     20072        newvar.setncattr('_CoordinateAxisType', 'Lat')   
     20073
     20074        newvarinpt = newnc.createVariable('locinpt','i',('inpts'))
     20075        basicvardef(newvarinpt, 'locinpt', 'input point located: 0: no, 1: yes','-')
     20076
     20077        newvarindiff = newnc.createVariable('locindiff','f4',('inpts'))
     20078        basicvardef(newvarindiff, 'locindiff', 'distance between input point and its final location','degree')
     20079
     20080# map variable
     20081        Lmapvalsshape = len(omapvals.shape)
     20082        if Lmapvalsshape == 1:
     20083            newvar = newnc.createVariable(varn, nctype(vartype), ('pts'), fill_value=fval)
     20084        else:
     20085            newvar = newnc.createVariable(varn, nctype(vartype), ('lat','lon'), fill_value=fval)
     20086
     20087        basicvardef(newvar, varn, varlongname, varunits)
     20088        newvar.setncattr('coordinates', 'lon lat')
     20089
     20090        if Lmapvalsshape == 1:
     20091            newvarin = newnc.createVariable('inpts', 'i4', ('pts'))   
     20092        else:
     20093            newvarin = newnc.createVariable('inpts', 'i4', ('lat','lon'))
     20094
     20095        basicvardef(newvarin, 'inpts', 'Equivalent point from the input source', '-')
     20096        newvar.setncattr('coordinates', 'lon lat')
     20097
     20098    else:
     20099        print '  ' + fname + "File '" + ofile + "' exists !!"
     20100        print '    reading values from file'
     20101        newnc = NetCDFFile(ofile, 'a')
     20102        newvar = newnc.variables[varn]
     20103        newvarin = newnc.variables['inpts']
     20104        newvarinpt = newnc.variables['locinpt']
     20105        newvarindiff = newnc.variables['locindiff']
     20106
     20107    amsk = np.arange(3)
     20108    amsk = ma.masked_equal(amsk, 0)
     20109
     20110#    fraclon = projlon[::Nmaplon*0.1]
     20111#    fraclat = projlat[::Nmaplat*0.1]
     20112
     20113#    print 'fraclon________________', fraclon.shape
     20114#    print fraclon
     20115
     20116#    print 'fraclat________________', fraclat.shape
     20117#    print fraclat
     20118
     20119# Reducing the searching points
     20120    newvarinvals = newvarinpt[:]
     20121    maskpt = np.where(newvarinvals.mask == True, False, True)
     20122    points = np.arange(Ninpts)
     20123    mapoints = ma.array(points, mask=maskpt)
     20124    ptsf = mapoints.compressed()
     20125
     20126    Nptsf = len(ptsf)
     20127    print Ninpts,'Npoints to find:', len(ptsf), ptsf[0:10], newvarindiff[ptsf[0:10]]
     20128# Error at 150024, 150025, 151709, 153421
     20129    print '  ' + fname + ': from:', Ninpts,'re-locating:',Nptsf,'points...'
     20130    if kind == 'Goode':
     20131#        newvar,newvarin,newvarinpt,newvarindiff =                                    \
     20132#          fin.module_forinterpolate.coarseinterpolate(projlon, projlat, lonvs,       \
     20133#          latvs, percen, mindiff, ivar, dimx, dimy, ninpts)
     20134        for ir in range(ptsf[0],Ninpts,fracs):
     20135            iri = ir
     20136            ire = ir + fracs + 1
     20137            print iri,',',ire
     20138#            percendone(iri,Ninpts,0.5,'done:')
     20139            pts = np.arange(iri,ire,1, dtype='int32')
     20140            lonvss = lonvs[iri:ire].astype('float64')
     20141            latvss = latvs[iri:ire].astype('float64')
     20142            ovars = ovar[iri:ire].astype('float64')
     20143
     20144#            newvar,newvarin,newvarinpt[pts],newvarindiff[pts] =                      \
     20145#              fin.module_forinterpolate.coarseinterpolateexact(projlon, projlat,     \
     20146#              lonvss, latvss, fracd, mindiff, ovars)
     20147# Slow way
     20148            newvar,newvarin,newvarinpt[pts],newvarindiff[pts] =                      \
     20149              fin.module_forinterpolate.interpolate(projlon, projlat,                \
     20150              lonvss, latvss, np.float64(mindiff), ovars)
     20151            newnc.sync()
     20152
     20153    elif kind == 'lonlat':
     20154        for iv in range(Ntotvals):
     20155            difflat = np.abs(projlat - latvs[iv])
     20156            mindiffL = np.min(difflat)
     20157            ilat = index_mat(difflat, mindiffL)
     20158            difflon = np.abs(projlon - lonvs[iv])
     20159            mindiffl = difflon.min()
     20160            ilon = index_mat(difflon, mindiffl)
     20161            percendone(iv,Ntotvals,0.5,'done:')
     20162            if mindiffl > mindiff or mindiffL > mindiff:
     20163                print errormsg
     20164                print '  ' + fname + ': for point #', iv,'lon,lat in incomplet map:',   \
     20165                  lonvs[iv], ',', latvs[iv], 'there is not a set of lon,lat in the ' +  \
     20166                  'completed map closer than: ', mindiff, '!!'
     20167                print '    minimum difflon:', np.min(difflon), 'difflat:', np.min(difflat)
     20168                quit()
     20169
     20170            if ilon >= 0 and ilat >= 0:
     20171                newvar[ilat,ilon] = ovar[iv]         
     20172                newnc.sync()
     20173            else:
     20174                print errormsg
     20175                print '  ' + fname + ': point iv:',iv,'at',lonvs[iv],',',latvs[iv],' not relocated !!'
     20176                print '    mindiffl:',mindiffl,'mindiffL:',mindiffL,'ilon:',ilon,'ilat:',ilat
     20177                quit(-1)
     20178
     20179        newnc.sync()
     20180
     20181    elif kind == 'lonlat_dxdyFix':
     20182
     20183        for iv in range(Ntotvals):
     20184#        for iv in range(15):
     20185#            print 'Lluis:',iv,'lon lat:',lonvs[iv],',',latvs[iv]
     20186            difflat = np.abs(projlat - latvs[iv])
     20187            mindiffL = np.min(difflat)
     20188            ilat = index_mat(difflat, mindiffL)
     20189#            print '  Lluis mindiffL:',mindiffL,'<> ilat:',ilat,'lat_ilat:',projlat[ilat],'Nprojlon:',Nprojlon[ilat],'shape oporjlon:',oprojlon[ilat,:].shape, type(oprojlon[ilat,:])
     20190            loncirc = np.zeros((Nprojlon[ilat]), dtype=np.float)
     20191            loncirc[:] = np.asarray(oprojlon[ilat,0:int(Nprojlon[ilat])].flatten())
     20192            difflon = np.abs(loncirc - lonvs[iv])
     20193            mindiffl = difflon.min()
     20194            ilon = index_mat(difflon, mindiffl)
     20195#            print '  difflon:',type(difflon),'shape difflon',difflon.shape,'shape loncirc:', loncirc.shape
     20196#            print '  Lluis mindiffl:',mindiffl,'<> ilon:',ilon,'oprojlon:', loncirc[ilon]
     20197
     20198            percendone(iv,Ntotvals,0.5,'done:')
     20199            if mindiffl > mindiff or mindiffL > mindiff:
     20200                print errormsg
     20201                print '  ' + fname + ': for point #', iv,'lon,lat in incomplet map:',   \
     20202                  lonvs[iv], ',', latvs[iv], 'there is not a set of lon,lat in the ' +  \
     20203                  'completed map closer than: ', mindiff, '!!'
     20204                print '    minimum difflon:', np.min(difflon), 'difflat:', np.min(difflat)
     20205                quit()
     20206
     20207            if ilon >= 0 and ilat >= 0:
     20208                if Lmapvalsshape ==1:
     20209                    newvar[iv] = ovar[iv]           
     20210                newnc.sync()
     20211            else:
     20212                print errormsg
     20213                print '  ' + fname + ': point iv:',iv,'at',lonvs[iv],',',latvs[iv],' not relocated !!'
     20214                print '    mindiffl:',mindiffl,'mindiffL:',mindiffL,'ilon:',ilon,'ilat:',ilat
     20215                quit(-1)
     20216
     20217        newnc.sync()
     20218# Global attributes
     20219##
     20220    newnc.setncattr('script',  fname)
     20221    newnc.setncattr('version',  '1.0')
     20222    newnc.setncattr('author',  'L. Fita')
     20223    newnc.setncattr('institution',  'Laboratoire de Meteorology Dynamique')
     20224    newnc.setncattr('university',  'Pierre et Marie Curie')
     20225    newnc.setncattr('country',  'France')
     20226    for attrs in onc.ncattrs():
     20227        attrv = onc.getncattr(attrs)
     20228        attr = set_attribute(newnc, attrs, attrv)
     20229
     20230    newnc.sync()
     20231    onc.close()
     20232    newnc.close()
     20233
     20234    print fname + ": Successfull written of file: '" + ofile + "' !!"
     20235
     20236    return
     20237
    1984220238#Partialmap_Entiremap('longitude,latitude,std,5000.,Goode,Goode_5km.nc', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/#carteveg5km.nc', 'vegetation_map')
    1984320239#Partialmap_Entiremap('longitude,latitude,std,5000.,lonlat_dxdyFix', '/home/lluis/etudes/DYNAMICO/ORCHIDEE/interpolation/data/carteveg5km.nc', 'vegetation_map')
Note: See TracChangeset for help on using the changeset viewer.