Changeset 733 in lmdz_wrf for trunk

Timestamp:

Apr 27, 2016, 11:28:38 AM (9 years ago)

Author:

lfita

Message:

Adding new modules on 'module_ForInterpolate':

• Interpolate: finds the closest grid point within a projection
• lonlatFind: Function to search a given value from a coarser version of the data

Location:

trunk/tools

Files:

• module_ForInterpolate.F90 (modified) (5 diffs)
• nc_var_tools.py (modified) (1 diff)

trunk/tools/module_ForInterpolate.F90

@@ -197 +197 @@
 END SUBROUTINE CoarselonlatFind
 
+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, ivar, newvar, newvarin, &
   newvarinpt, newvarindiff, dimx, dimy, Ninpts)
@@ -211 +261 @@
   INTEGER, INTENT(in)                                    :: Ninpts
   REAL(r_k), DIMENSION(Ninpts), INTENT(in)               :: ivar, lonvs, latvs
-  REAL(r_k)                                              :: mindiff, percen
+  REAL(r_k), INTENT(in)                                  :: mindiff, percen
   REAL(r_k), DIMENSION(dimx,dimy), INTENT(out)           :: newvar
   INTEGER, DIMENSION(dimx,dimy), INTENT(out)             :: newvarin
@@ -280 +330 @@
   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
+!  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
@@ -295 +345 @@
         fractionlat, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, ilonlat, mindiffLl)
 
+      PRINT *,'  Lluis: iv',iv,', ', mindiffLl,'<= ',mindiff,' ilonlat:',ilonlat
       IF (mindiffLl <= mindiff) THEN 
 !        percendone(iv,Ninpts,0.5,'done:')
@@ -329 +380 @@
 END SUBROUTINE CoarseInterpolate
 
+SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, ivar, newvar, newvarin, newvarinpt,   &
+  newvarindiff, 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)               :: ivar, lonvs, latvs
+  REAL(r_k), INTENT(in)                                  :: mindiff
+  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), DIMENSION(dimx,dimy)                        :: difflonlat
+  REAL(r_k)                                              :: mindiffLl
+  INTEGER                                                :: Ninpts1
+  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
+! 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 = 'Interpolate'
+  Ninpts1 = Ninpts/100
+
+  extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
+  extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
+
+  DO iv=1,Ninpts
+    IF (newvarinpt(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 = index2DArrayR(difflonlat, dimx, dimy, 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 Interpolate
+
 END MODULE module_ForInterpolate

trunk/tools/nc_var_tools.py

@@ -19746 +19746 @@
             ovars = ovar[iri:ire].astype('float64')
 
+#            newvar,newvarin,newvarinpt[pts],newvarindiff[pts] =                      \
+#              fin.module_forinterpolate.coarseinterpolate(projlon, projlat,          \
+#              lonvss, latvss, np.float64(fracd), np.float64(mindiff), ovars)
             newvar,newvarin,newvarinpt[pts],newvarindiff[pts] =                      \
-              fin.module_forinterpolate.coarseinterpolate(projlon, projlat,          \
-              lonvss, latvss, np.float64(fracd), np.float64(mindiff), ovars)
+              fin.module_forinterpolate.interpolate(projlon, projlat,                \
+              lonvss, latvss, np.float64(mindiff), ovars)
             newnc.sync()