Changeset 1180 in lmdz_wrf for trunk/tools

Timestamp:

Oct 12, 2016, 1:08:30 AM (9 years ago)

Author:

lfita

Message:

Working version with `npp' of 'reproject'

Location:

trunk/tools

Files:

• module_ForInterpolate.F90 (modified) (12 diffs)
• nc_var_tools.py (modified) (6 diffs)

trunk/tools/module_ForInterpolate.F90

@@ -739 +739 @@
   REAL(r_k)                                              :: mindiffLl, dist
   REAL(r_k), DIMENSION(idimx,idimy)                      :: difflonlat
+  REAL(r_k), DIMENSION(idimx-1,idimy)                    :: difflon
+  REAL(r_k), DIMENSION(idimx,idimy-1)                    :: difflat
   REAL(r_k), DIMENSION(2)                                :: extremelon, extremelat, ipos
   INTEGER, DIMENSION(2)                                  :: iLl
@@ -760 +762 @@
   extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
   extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
-
+ 
   ! Using case outside loop to be more efficient
   SELECT CASE(TRIM(intkind))
@@ -815 +817 @@
           mindiffLl = MINVAL(difflonlat)
           ipos = index2DArrayR(difflonlat, idimx, idimy, mindiffLl)
-          outLlw(1:2,1,i,j) = ipos
+          outLlw(1,1,i,j) = REAL(ipos(1))
+          outLlw(2,1,i,j) = REAL(ipos(2))
           outLlw(3,1,i,j) = 1.
-          ix = outLlw(1,1,i,j)
-          iy = outLlw(2,1,i,j)
-          PRINT *,i,j,projlon(i,j),projlat(i,j),':',mindiffLl,'!',outLlw(1:2,1,i,j),'|',inlonv(ix,iy), &
-            inlatv(ix,iy), '<>', MINVAL(inlonv), MINVAL(inlatv)
-          IF (i > 1) THEN
-            PRINT *,inlonv(1,1),inlatv(1,1),':',index2DArrayR(inlonv, idimx, idimy, MINVAL(inlonv))
-            STOP
-          END IF
+          ix = INT(outLlw(1,1,i,j))
+          iy = INT(outLlw(2,1,i,j))
         END DO
       END DO
@@ -876 +873 @@
     CASE('dis')
       DO i=1, pdimx
-        DO j=1, pdimx
+        DO j=1, pdimy
           varout(i,j) = 0.
           DO iv=1, 4
@@ -888 +885 @@
     CASE('npp')
       DO i=1, pdimx
-        DO j=1, pdimx
+        DO j=1, pdimy
           ix = INT(outLlw(1,1,i,j))
           iy = INT(outLlw(2,1,i,j))
@@ -943 +940 @@
     CASE('dis')
       DO i=1, pdimx
-        DO j=1, pdimx
+        DO j=1, pdimy
           DO k=1, d3
             varout(i,j,k) = 0.
@@ -957 +954 @@
     CASE('npp')
       DO i=1, pdimx
-        DO j=1, pdimx
+        DO j=1, pdimy
+          ix = INT(outLlw(1,1,i,j))
+          iy = INT(outLlw(2,1,i,j))
           DO k=1, d3
-            ix = INT(outLlw(1,1,i,j))
-            iy = INT(outLlw(2,1,i,j))
             varout(i,j,k) = var3Din(ix,iy,k)
           END DO
@@ -982 +979 @@
   REAL(r_k), DIMENSION(idimx,idimy), INTENT(in)          :: inlonv, inlatv
   CHARACTER(LEN=50), INTENT(in)                          :: intkind
-  REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4), INTENT(in)    :: var4Din
+  REAL(r_k), DIMENSION(idimx,idimy,d3,d4), INTENT(in)    :: var4Din
   REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4), INTENT(out)   :: varout
 
@@ -1014 +1011 @@
     CASE('dis')
       DO i=1, pdimx
-        DO j=1, pdimx
+        DO j=1, pdimy
           DO k=1, d3
             DO l=1, d4
@@ -1030 +1027 @@
     CASE('npp')
       DO i=1, pdimx
-        DO j=1, pdimx
+        DO j=1, pdimy
           ix = INT(outLlw(1,1,i,j))
           iy = INT(outLlw(2,1,i,j))
@@ -1090 +1087 @@
     CASE('dis')
       DO i=1, pdimx
-        DO j=1, pdimx
+        DO j=1, pdimy
           DO k=1, d3
             DO l=1, d4
@@ -1108 +1105 @@
     CASE('npp')
       DO i=1, pdimx
-        DO j=1, pdimx
+        DO j=1, pdimy
           ix = INT(outLlw(1,1,i,j))
           iy = INT(outLlw(2,1,i,j))

trunk/tools/nc_var_tools.py

@@ -15860 +15860 @@
     olatvals0 = oprojlat[:]
 
+#    ilonvals0 = np.arange(0.,360.,360./oinlon.shape[0])
+#    ilatvals0 = np.arange(-90.,90.,180./oinlat.shape[0])
+#    olonvals0 = np.arange(0.5,359.5,359./oprojlon.shape[0])
+#    olatvals0 = np.arange(-89.5,89.5,179./oprojlat.shape[0])
+
     # Getting 2D longitudes and latitudes
     ilonvals, ilatvals = gen.lonlat2D(ilonvals0, ilatvals0)
@@ -15905 +15910 @@
     rminlon = np.min(rlonvals)
     rmaxlon = np.max(rlonvals)
+    iminlat = np.min(ilatvals)
+    imaxlat = np.max(ilatvals)
+    rminlat = np.min(rlatvals)
+    rmaxlat = np.max(rlatvals)
 
     if np.abs(imaxlon - rmaxlon) > 100.: 
@@ -15913 +15922 @@
         print '    shifting input longitudes...'
         if rmaxlon > 180.:
-            ilonvals = np.where(ilonvals < 0., ilonvals + 360., ilonvals)
+            ilonvals = np.where(ilonvals < 0., 360. + ilonvals, ilonvals)
         else:
             ilonvals = np.where(ilonvals > 180., ilonvals - 360., ilonvals)
+
+    flipy = False
+    if np.abs(ilatvals[0,0] - rlatvals[0,0]) > 90.: 
+        print warnmsg
+        print '  ' + fname + ': reshaping latitudes!'
+        print '    input minimum lat:', ilatvals[0,0]
+        print '    target minimum lat:', rlatvals[0,0]
+        print '    flip input latitudes...'
+        ilatvals = ilatvals[::-1,:]
+        flipy = True
 
     idx = ilonvals.shape[1]
@@ -15973 +15992 @@
         print '  ' + fname + ": re-projecting '" + varn + "' ..."
         ovarin = inf.variables[varn]
-        varin = ovarin[:]
+        if flipy:
+            varin = ovarin[...,::-1,:]
+        else:
+            varin = ovarin[:]
 
         Nvarind = len(varin.shape)
@@ -15986 +16008 @@
 
         if Nvarind == 2:
-            varout = fin.module_forinterpolate.var2d_intproj(var2din=varint,         \
+            varoutt = fin.module_forinterpolate.var2d_intproj(var2din=varint,        \
               inlonv=ilont, inlatv=ilatt, projlon=rlont, projlat=rlatt, intkind=kind,\
               idimx=idx, idimy=idy, pdimx=rdx, pdimy=rdy)
         elif Nvarind == 3:
-            varout = fin.module_forinterpolate.var3d_intproj(var3din=varint,         \
+            varoutt = fin.module_forinterpolate.var3d_intproj(var3din=varint,        \
               inlonv=ilont, inlatv=ilatt, projlon=rlont, projlat=rlatt, intkind=kind,\
               idimx=idx, idimy=idy, pdimx=rdx, pdimy=rdy, d3=varin.shape[0])
         elif Nvarind == 4:
-            varout = fin.module_forinterpolate.var4d_intproj(var4din=varint,         \
+            varoutt = fin.module_forinterpolate.var4d_intproj(var4din=varint,        \
               inlonv=ilont, inlatv=ilatt, projlon=rlont, projlat=rlatt, intkind=kind,\
               idimx=idx, idimy=idy, pdimx=rdx, pdimy=rdy, d3=varin.shape[0],         \
               d4=varin.shape[1])
         elif Nvarind == 5:
-            varout = fin.module_forinterpolate.var5d_intproj(var5din=varint,         \
+            varoutt = fin.module_forinterpolate.var5d_intproj(var5din=varint,        \
               inlonv=ilont, inlatv=ilatt, projlon=rlont, projlat=rlatt, intkind=kind,\
               idimx=idx, idimy=idy, pdimx=rdx, pdimy=rdy, d3=varin.shape[0],         \
@@ -16007 +16029 @@
             print '  ' + fname + ': rank:', Nvarind, 'for input data not ready !!'
             quit(-1)
-    
+
+        # Recovering 'python' C-like space
+        varout = varoutt.transpose()
+
         # Variable's dimensions and their related variables
         vdims = list(ovarin.dimensions)