Changeset 2212 in lmdz_wrf for trunk/tools

Timestamp:

Nov 2, 2018, 7:01:51 PM (6 years ago)

Author:

lfita

Message:

Addding the `range_faces' with all its glory...

Location:

trunk/tools

Files:

• diag_tools.py (modified) (4 diffs)
• diagnostics.py (modified) (5 diffs)
• module_ForDiagnostics.f90 (modified) (2 diffs)
• module_ForDiagnosticsVars.f90 (modified) (8 diffs)
• variables_values.dat (modified) (2 diffs)

trunk/tools/diag_tools.py

@@ -1370 +1370 @@
     return var1, var2, vardims, varvdims
 
-def Forcompute_range_faces(lon, lat, hgt, face, dimns, dimvns):
+def Forcompute_range_faces(lon, lat, hgt, face, Nfilt, dimns, dimvns):
     """ Function to compute faces [uphill, valley, downhill] of sections of a mountain
       rage, along a given face
@@ -1377 +1377 @@
       [lat]= latitude values (assuming [y,x]) [degrees north]
       [hgt]= height values (assuming [y,x]) [m]
+      face= which face (axis along which produce slices) to use to compute the faces: WE, SN
+      Nfilt= number of grid points to use to filter the topographic values
       [dimns]= list of the name of the dimensions of [smois]
       [dimvns]= list of the name of the variables with the values of the 
@@ -1391 +1393 @@
         dy = hgt.shape[0]
 
-        facest=fdin.module_fordiagnostics.compute_range_faces(lon=lon[:].transpose(),\
-          lat=lat[:].transpose(), hgt=hgt[:].transpose(), face=face, d1=dx, d2=dy)
-        faces = facest.transpose()
+        dhgtt, peakst, valleyst, ofacest, ffacest =                                  \
+          fdin.module_fordiagnostics.compute_range_faces(lon=lon[:].transpose(),     \
+          lat=lat[:].transpose(), hgt=hgt[:].transpose(), face=face, nfilt=Nfilt,    \
+          d1=dx, d2=dy)
+        dhgt = dhgtt.transpose()
+        peaks = peakst.transpose()
+        valleys = valleyst.transpose()
+        origfaces = ofacest.transpose()
+        filtfaces = ffacest.transpose()
     else:
         print errormsg
@@ -1400 +1408 @@
         quit(-1)
 
-    return faces, vardims, varvdims
+    return dhgt, peaks, valleys, origfaces, filtfaces, vardims, varvdims
 
 ####### ###### ##### #### ### ## # END Fortran diagnostics

trunk/tools/diagnostics.py

@@ -101 +101 @@
 
 # Variables not to check
-NONcheckingvars = ['accum', 'cllmh', 'deaccum', 'face:WE', 'face:SN', 'TSrhs',       \
+NONcheckingvars = ['accum', 'cllmh', 'deaccum', 'face', 'TSrhs',                     \
   'TStd', 'TSwds', 'TSwss',                                                          \
   'WRFbils',                                                                         \
@@ -415 +415 @@
             if depvars[0] == 'accum': diagn='accum'
             for depv in depvars:
+                # Checking without extra arguments of a depending variable (':', separated)
+                if depv.find(':') != -1: depv=depv.split(':')[0]
                 if not ncobj.variables.has_key(depv) and not                         \
                   gen.searchInlist(NONcheckingvars, depv) and                        \
@@ -679 +681 @@
         ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
 
-# range_faces: LON, LAT, HGT, 'face:['WE'/'SN']'
+# range_faces: LON, LAT, HGT, 'face:['WE'/'SN']:[Nptfilt]'
     elif diagn == 'range_faces':
             
@@ -685 +687 @@
         var1 = ncobj.variables[depvars[1]][:]
         var2 = ncobj.variables[depvars[2]][:]
-        var3 = depvars[3].split(':')[1]
-
+        face = depvars[3].split(':')[1]
+        Nptfilt = depvars[3].split(':')[2]
+
+        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
+        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
         lon, lat = gen.lonlat2D(var0, var1)
         if len(var2.shape) == 3:
@@ -692 +697 @@
             print '  ' + diagn + ": shapping to 2D variable '" + depvars[2] + "' !!"
             hgt = var2[0,:,:]
+            dnamesvar.pop(0)
+            dvnamesvar.pop(0)            
         else:
             hgt = var2[:]
 
-        dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
-        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
-
-        diagout, diagoutd, diagoutvd = diag.Forcompute_range_faces(lon, lat, hgt,    \
-          var3, dnamesvar, dvnamesvar)
-
-        # Removing the nonChecking variable-dimensions from the initial list
-        varsadd = []
-        diagoutvd = list(dvnames)
-        for nonvd in NONchkvardims:
-            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
-            varsadd.append(nonvd)
-        ncvar.insert_variable(ncobj, 'range_faces', diagout, diagoutd, diagoutvd,    \
+        dhgt, peaks, valleys, origfaces, diagout, diagoutd, diagoutvd =              \
+          diag.Forcompute_range_faces(lon, lat, hgt, face, Nptfilt, dnamesvar,       \
+          dvnamesvar)
+
+        # Removing the nonChecking variable-dimensions from the initial list
+        varsadd = []
+        diagoutvd = list(dvnames)
+        for nonvd in NONchkvardims:
+            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
+            varsadd.append(nonvd)
+
+        ncvar.insert_variable(ncobj, 'hgtderivdeg', dhgt, diagoutd, diagoutvd,       \
           newnc)
+        ovar = newnc.variables['orogderivdeg']
+        if face == 'WE': axis = 'lon'
+        elif face == 'SN': axis = 'lat'
+        ncvar.set_attribute(ovar, 'deriv', axis)
+
+        ncvar.insert_variable(ncobj, 'peak', peaks, diagoutd, diagoutvd, newnc)
+        ncvar.insert_variable(ncobj, 'valley', valleys, diagoutd, diagoutvd, newnc)
+
+        ncvar.insert_variable(ncobj, 'rangefaces', diagout, diagoutd, diagoutvd,    \
+          newnc)
+        newnc.renameVariable('rangefaces', 'rangefacesfilt')
+        ovar = newnc.variables['rangefacesfilt']
+        ncvar.set_attribute(ovar, 'face', face)
+        ncvar.set_attributek(ovar, 'grid_points_filter', Nptfilt, 'I')
+
+        ncvar.insert_variable(ncobj, 'rangefaces', origfaces, diagoutd, diagoutvd,  \
+          newnc)
+        ovar = newnc.variables['rangefaces']
+        ncvar.set_attribute(ovar, 'face', face)
 
 # mrso: total soil moisture SMOIS, DZS

trunk/tools/module_ForDiagnostics.f90

@@ -949 +949 @@
   END SUBROUTINE compute_fog_FRAML50
 
-  SUBROUTINE compute_range_faces(d1, d2, lon, lat, hgt, face, faces)
+  SUBROUTINE compute_range_faces(d1, d2, lon, lat, hgt, face, Nfilt, derivhgt, peaks, valleys,        &
+    origfaces, filtfaces)
 ! Subroutine to compute faces [uphill, valleys, downhill] of a mountain range along a given face
 
     IMPLICIT NONE
 
-    INTEGER, INTENT(in)                                  :: d1, d2
+    INTEGER, INTENT(in)                                  :: d1, d2, Nfilt
     REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: lon, lat, hgt
     CHARACTER(len=*)                                     :: face
-    INTEGER, DIMENSION(d1,d2), INTENT(out)               :: faces
+    REAL(r_k), DIMENSION(d1,d2), INTENT(out)             :: derivhgt
+    INTEGER, DIMENSION(d1,d2), INTENT(out)               :: origfaces, filtfaces, peaks, valleys
  
 ! Local
     INTEGER                                              :: i, j
+    INTEGER                                              :: Npeaks, Nvalleys
+    INTEGER, DIMENSION(d1)                               :: ipeaks1, ivalleys1
+    INTEGER, DIMENSION(d2)                               :: ipeaks2, ivalleys2
 
 !!!!!!! Variables
@@ -966 +971 @@
 ! lat: latitude [degrees north]
 ! hgt: topograpical height [m]
+! face: which face (axis along which produce slices) to use to compute the faces: WE, SN
+! Nfilt: number of grid points to use to filter the topographic values
+! derivhgt: topograpic derivative along axis [m deg-1]
+! peaks: peak point
+! valleys: valley point
+! origfaces: original faces (-1, downhill; 0: valley; 1: uphill)
+! filtfaces: filtered faces (-1, downhill; 0: valley; 1: uphill)
 
     fname = 'compute_range_faces'
 
+    peaks = 0
+    valleys = 0
     IF (TRIM(face) == 'WE') THEN
       DO j=1, d2
-        CALL var_range_faces(d1, lon(:,j), lat(:,j), hgt(:,j), faces(:,j))
+        PRINT *,'Lluis:',j
+        CALL var_range_faces(d1, lon(:,j), lat(:,j), hgt(:,j), Nfilt, derivhgt(:,j), Npeaks,          &
+          ipeaks1, Nvalleys, ivalleys1, origfaces(:,j), filtfaces(:,j))
+        DO i=1, Npeaks
+          peaks(ipeaks1(i),j) = 1
+        END DO
+        DO i=1, Nvalleys
+          valleys(ivalleys1(i),j) = 1
+        END DO
       END DO
     ELSE IF (TRIM(face) == 'SN') THEN
       DO i=1, d1
-        CALL var_range_faces(d2, lon(i,:), lat(i,:), hgt(i,:), faces(i,:))
+        CALL var_range_faces(d2, lon(i,:), lat(i,:), hgt(i,:), Nfilt, derivhgt(i,:), Npeaks,          &
+          ipeaks2, Nvalleys, ivalleys2, origfaces(i,:), filtfaces(i,:))
+        DO j=1, Npeaks
+          peaks(i,ipeaks2(j)) = 1
+        END DO
+        DO j=1, Nvalleys
+          valleys(i,ivalleys2(j)) = 1
+        END DO
       END DO
     ELSE

trunk/tools/module_ForDiagnosticsVars.f90

@@ -9 +9 @@
   USE module_definitions
   USE module_generic
+  USE module_scientific
 
   IMPLICIT NONE
@@ -1579 +1580 @@
   END SUBROUTINE var_fog_FRAML50
 
-  SUBROUTINE var_range_faces(d, lon, lat, hgt, faces)
+  SUBROUTINE var_range_faces(d, lon, lat, hgt, filt, dhgt, Npeaks, ipeaks, Nvalleys, ivalleys,        &
+    faces0, faces)
 ! Subroutine to compute faces [uphill, valleys, downhill] of a monuntain range along a face
 
     IMPLICIT NONE
 
-    INTEGER, INTENT(in)                                  :: d
-    REAL, DIMENSION(d), INTENT(in)                       :: lon, lat, hgt
-    INTEGER, DIMENSION(d), INTENT(out)                   :: faces
+    INTEGER, INTENT(in)                                  :: d, filt
+    REAL(r_k), DIMENSION(d), INTENT(in)                  :: lon, lat, hgt
+    REAL(r_k), DIMENSION(d), INTENT(out)                 :: dhgt
+    INTEGER, DIMENSION(d), INTENT(out)                   :: ipeaks, ivalleys, faces0, faces
+    INTEGER, INTENT(out)                                 :: Npeaks, Nvalleys
 
 ! Local
-    INTEGER                                              :: i, Nfaces, Npeaks
+    INTEGER                                              :: i, j, k, l, m 
+    INTEGER                                              :: iface
+    INTEGER                                              :: Nfaces, Nfaces1, Npeaks1, Nvalleys1
+    INTEGER                                              :: fbeg, fend
     INTEGER, DIMENSION(1)                                :: ihmax
-    REAL                                                 :: hgtmax
-    INTEGER, DIMENSION(d)                                :: ddhgt, Ndhgt, ipeaks
-    REAL, DIMENSION(d)                                   :: dhgt, peaks
+    REAL(r_k)                                            :: hgtmax
+    INTEGER, DIMENSION(d)                                :: ddhgt, Ndhgt
+    INTEGER, DIMENSION(d)                                :: faces1, Ndhgt1, ipeaks1, ivalleys1
+    REAL(r_k), DIMENSION(d)                              :: peaks, valleys, peaks1, valleys1
+    LOGICAL                                              :: peakwithin, valleywithin
 
 !!!!!!! Variables
@@ -1599 +1608 @@
 ! lat: latitude [degrees north]
 ! hgt: topograpical height [m]
+! filt: number of grid points to use to filter the topographic values. used to define:
+!   the minimum length of a given section
+!   running mean filter
 
     fname = 'var_range_faces'
@@ -1610 +1622 @@
     PRINT *, 'height max:', hgtmax, 'location:', ihmax
 
+    ! Filtering height [NOT needed]
+    ! CALL runmean_F1D(d, hgt, filt, 'progressfill', hgtrunmean)
+
     ! range slope
     dhgt(1:d) = hgt(2:d) - hgt(1:d-1)
@@ -1615 +1630 @@
     PRINT *, 'slope:', dhgt
 
-    ! Classification
+    ! Initial classification
     Npeaks = 0
+    Nvalleys = 0
     DO i=1, d-1
       IF (dhgt(i) > 0.) THEN
-        faces(i) = 1
+        faces0(i) = 1
       ELSE IF (dhgt(i) < 0.) THEN
-        faces(i) = -1
+        faces0(i) = -1
       ELSE
-        faces(i) = 0
+        faces0(i) = 0
       END IF
       ! peaks
@@ -1631 +1647 @@
         ipeaks(Npeaks) = i+1
       END IF
+      ! valleys
+      IF (dhgt(i) < 0. .AND. dhgt(i+1) > 0.) THEN
+        Nvalleys = Nvalleys + 1
+        valleys(Nvalleys) = hgt(i+1)
+        ivalleys(Nvalleys) = i+1
+      END IF
     END DO
 
-    PRINT *, 'faces:', faces
+    PRINT *, 'faces:', faces0
     PRINT *, Npeaks, ' peaks:', peaks(1:Npeaks), ' ipeak:', ipeaks(1:Npeaks)
+    PRINT *, Nvalleys, ' valleys:', valleys(1:Nvalleys), ' ivalley:', ivalleys(1:Nvalleys)
 
     ! tendency of the slope
@@ -1641 +1664 @@
     Ndhgt(Nfaces) = 1
     DO i=2, d-1
-      IF (faces(i) /= faces(i-1)) THEN
+      IF (faces0(i) /= faces0(i-1)) THEN
         Nfaces = Nfaces + 1
         Ndhgt(Nfaces) = 1
@@ -1653 +1676 @@
     PRINT *, Nfaces, ' length faces:', Ndhgt(1:Nfaces)
 
-    PRINT *, '# lon[0] lat[1] heights[2] slope[3] faces[4] ddhgt[5]'
+    ! classification using filt. 
+    !   On that section of the hill smaller than the selected filter. Replace values by the ones from 
+    !     the previous section (except if it is already valley!) 
+    !     e.g.: filt= 3
+    !       faces0 = 0, 0, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1 -1, -1, 0, 0
+    !       faces = 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1 -1, -1, 0, 0
+    faces1 = faces0
+    Nfaces1 = 1
+    Npeaks1 = 0
+    Nvalleys1 = 0
+    Ndhgt1(1) = Ndhgt(1)
+    DO i=2, Nfaces
+      fbeg = SUM(Ndhgt(1:i-1))
+      fend = fbeg + Ndhgt(i)
+      peakwithin = .FALSE.
+      valleywithin = .FALSE.
+      DO j=1, Npeaks
+        IF (ipeaks(j) == fbeg) THEN
+          k = j
+          peakwithin = .TRUE.
+          EXIT
+        END IF
+      END DO
+      DO j=1, Nvalleys
+        IF (ivalleys(j) == fbeg) THEN
+          m = j
+          valleywithin = .TRUE.
+          EXIT
+        END IF
+      END DO
+
+      IF (Ndhgt(i) < filt) THEN
+        PRINT *, 'Lluis replacing !!!', i, ':', fbeg, fend, '<>', faces0(fbeg-1)
+        IF (faces0(fbeg) /= 0) faces1(fbeg:fend) = faces0(fbeg-1)
+        Ndhgt1(Nfaces1) = Ndhgt1(Nfaces1) + fend - fbeg + 1
+        peakwithin = .FALSE.
+        valleywithin = .FALSE.
+      ELSE
+        Nfaces1 = Nfaces1 + 1
+        Ndhgt1(Nfaces1) = Ndhgt(i)
+      END IF
+
+      IF (peakwithin) THEN
+        Npeaks1 = Npeaks1 + 1
+        peaks1(Npeaks1) = peaks(k)
+      END IF
+      IF (valleywithin) THEN
+        Nvalleys1 = Nvalleys1 + 1
+        valleys1(Nvalleys1) = valleys(m)
+      END IF
+    END DO
+
+    ! Introducing valleys between peaks lower than the highest peak in the section
+    faces = faces1
+    IF (Npeaks1 > 1) THEN
+      DO i=1,Npeaks1,2
+        fbeg = ipeaks1(i)
+        fend = ipeaks1(i+1)
+        faces(fbeg:fend) = 0
+      END DO   
+    END IF
+
+    PRINT *, '# lon[0] lat[1] heights[2] slope[3] faces0[4] faces1[5] faces1[6] ddhgt[7]'
     DO i=1, d
-      PRINT *, lon(i), lat(i), hgt(i), dhgt(i), faces(i), ddhgt(i)
+      PRINT *, lon(i), lat(i), hgt(i), dhgt(i), faces0(i), faces1(i), faces(i), ddhgt(i)
     END DO
 

trunk/tools/variables_values.dat

@@ -296 +296 @@
 HGT, orog, orography,  0., 3000., surface|altitude, m,terrain, $orog$, orog
 HGT_M, orog, orography,  0., 3000., surface|altitude, m,terrain, $orog$, orog
+orogderivdeg, orogderivdeg, orography_derivative_degrees,  0., 3000., derivative|along|degree|coordinate|of|surface|altitude, mdeg-1,seismic, $orog_{deri}^{deg}$, orog_deriv_deg
+hgtderivdeg, orogderivdeg, orography_derivative_degrees,  0., 3000., derivative|along|degree|coordinate|of|surface|altitude, mdeg-1,seismic, $orog_{deri}^{deg}$, orog_deriv_deg
+peak, peak, peak_point, 0., 1., peak|grd|point, 1, Binary, $peak$, peak
 plfc, plfc, pressure_level_free_convection, 100., 101000., pressure|of|level|free|convection, Pa, BuPu, $lfc^{p}$, lfc_p
 LFCP, plfc, pressure_level_free_convection, 100., 101000., pressure|of|level|free|convection, Pa, BuPu, $lfc^{p}$, lfc_p
@@ -539 +542 @@
 Meridional wind, va, northward_wind, -30., 30., northward|wind, ms-1, seismic, $va$, va
 LVITV, va, northward_wind, -30., 30., northward|wind, ms-1, seismic, $va$, va
+valley, valley, valley_point, 0., 1., valley|grd|point, 1, Binary, $valley$, valley
 vas, vas, northward_wind, -30., 30., northward|10m|wind, ms-1, seismic, $vas$, vas
 VAS, vas, northward_wind, -30., 30., northward|10m|wind, ms-1, seismic, $vas$, vas