Changeset 2215 in lmdz_wrf for trunk/tools

Timestamp:

Nov 6, 2018, 8:14:36 PM (6 years ago)

Author:

lfita

Message:

Final working version of `range_faces'

Location:

trunk/tools

Files:

• diag_tools.py (modified) (4 diffs)
• diagnostics.inf (modified) (1 diff)
• diagnostics.py (modified) (10 diffs)
• module_ForDiagnostics.f90 (modified) (4 diffs)
• module_ForDiagnosticsVars.f90 (modified) (10 diffs)
• variables_values.dat (modified) (3 diffs)

trunk/tools/diag_tools.py

@@ -1370 +1370 @@
     return var1, var2, vardims, varvdims
 
-def Forcompute_range_faces(lon, lat, hgt, face, Nfilt, newrng, dimns, dimvns):
+def Forcompute_range_faces(lon, lat, hgt, ds, face, dsfilt, dsnewrng, hvalleyrng,    \
+  dimns, dimvns):
     """ Function to compute faces [uphill, valley, downhill] of sections of a mountain
       rage, along a given face
@@ -1377 +1378 @@
       [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
-      newrng= number of grid points as distance to start a new mountain range
+      [ds]= distance between grid points (assuming [y,x]) [m]
+      face= which face (axis along which produce slices) to use to compute the 
+        faces: WE, SN
+      dsfilt= distance to filter orography smaller scale of it [m]
+      dsnewrng= distance to start a new mountain range [m]
+      hvalleyrng: maximum height of a valley to mark change of range [m]
       [dimns]= list of the name of the dimensions of [smois]
       [dimvns]= list of the name of the variables with the values of the 
@@ -1385 +1389 @@
     """
     fname = 'Forcompute_range_faces'
+
+    print fname + 'Lluis:', dsfilt, dsnewrng, hvalleyrng
 
     vardims = dimns[:]
@@ -1397 +1403 @@
           ptrnghgtmaxt =                                                             \
           fdin.module_fordiagnostics.compute_range_faces(lon=lon[:].transpose(),     \
-          lat=lat[:].transpose(), hgt=hgt[:].transpose(), face=face, nfilt=Nfilt,    \
-          newrange=newrng, d1=dx, d2=dy)
+          lat=lat[:].transpose(), hgt=hgt[:].transpose(), dist=ds[:].transpose(),    \
+          face=face, dsfilt=dsfilt, dsnewrange=dsnewrng, hvalrng=hvalleyrng, d1=dx,  \
+          d2=dy)
+        print 'Finished!'
         hgtmax = hgtmaxt.transpose()
         pthgtmax = pthgtmaxt.transpose()

trunk/tools/diagnostics.inf

@@ -34 +34 @@
 psl, WRFpsl_ecmwf, PSFC@HGT@WRFt@WRFp@ZNU@ZNW
 psl, WRFpsl_ptarget, WRFp@PSFC@WRFt@HGT@QVAPOR
+rangefaces, range_faces, XLONG@XLAT@HGT@WRFdxdy@face:WE:10.:30.:200.
 rvor, rvor, U@V
 rvors, WRFrvors, U10@V10@WRFpos

trunk/tools/diagnostics.py

@@ -104 +104 @@
   'TStd', 'TSwds', 'TSwss',                                                          \
   'WRFbils',                                                                         \
-  'WRFclivi', 'WRFclwvi', 'WRFdens', 'WRFgeop',                                      \
+  'WRFclivi', 'WRFclwvi', 'WRFdens', 'WRFdx', 'WRFdxdy', 'WRFdy', 'WRFgeop',         \
   'WRFp', 'WRFtd',                                                                   \
   'WRFpos', 'WRFprc', 'WRFprls', 'WRFrh', 'LMDZrh', 'LMDZrhs',                       \
   'WRFrhs', 'WRFrvors',                                                              \
-  'WRFt', 'WRFtime', 'WRFua', 'WRFva', 'WRFwds', 'WRFwss', 'WRFheight', 'WRFz']
+  'WRFt', 'WRFtime', 'WRFua', 'WRFva', 'WRFwds', 'WRFwss', 'WRFheight', 'WRFz',      \
+  'WRFzg']
 
 # diagnostics not to check their dependeny
@@ -165 +166 @@
 WRFtime_compute = False
 WRFz_compute = False
+WRFdxdy_compute = False
 
 # File creation
@@ -189 +191 @@
     if dnv == 'WRFtime': WRFtime_compute = True
     if dnv == 'WRFz':WRFz_compute = True
+    if dnv == 'WRFdxdy':WRFdxdy_compute = True
 
 # diagnostics to compute
@@ -217 +220 @@
         if gen.searchInlist(depvars, 'WRFtime'): WRFtime_compute = True
         if gen.searchInlist(depvars, 'WRFz'): WRFz_compute = True
+        if gen.searchInlist(depvars, 'WRFdxdy'): WRFdxdy_compute = True
 
 # Dictionary with the new computed variables to be able to add them
@@ -396 +400 @@
     # Attributes of the variable
     Vvals = gen.variables_values('WRFz')
-    dictcompvars['WRFz'] = {'name': Vvals[0], 'standard_name': Vvals[1],          \
+    dictcompvars['WRFz'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
+      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
+
+if WRFdxdy_compute:
+    print '  ' + main + ': Retrieving dxdy: real distance between grid points ' +    \
+      'from WRF as dx=(XLONG(i+1)-XLONG(i))*DX/MAPFAC_M, dy=(XLAT(j+1)-XLAT(i))*DY/'+\
+      'MAPFAC_M, ds=sqrt(dx**2+dy**2)'
+    dimv = ncobj.variables['XLONG'].shape
+    WRFlon = ncobj.variables['XLONG'][0,:,:]
+    WRFlat = ncobj.variables['XLAT'][0,:,:]
+    WRFmapfac_m = ncobj.variables['MAPFAC_M'][0,:,:]
+    DX = ncobj.DX
+    DY = ncobj.DY
+
+    dimx = dimv[2]
+    dimy = dimv[1]
+
+    WRFdx = np.zeros((dimy,dimx), dtype=np.float)
+    WRFdy = np.zeros((dimy,dimx), dtype=np.float)
+
+    WRFdx[:,0:dimx-1]=(WRFlon[:,1:dimx]-WRFlon[:,0:dimx-1])*DX/WRFmapfac_m[:,0:dimx-1]
+    WRFdy[0:dimy-1,:]=(WRFlat[1:dimy,:]-WRFlat[0:dimy-1,:])*DY/WRFmapfac_m[0:dimy-1,:]
+    WRFds = np.sqrt(WRFdx**2 + WRFdy**2)
+
+    # Attributes of the variable
+    Vvals = gen.variables_values('WRFdx')
+    dictcompvars['WRFdx'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
+      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
+    Vvals = gen.variables_values('WRFdy')
+    dictcompvars['WRFdy'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
+      'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
+    Vvals = gen.variables_values('WRFds')
+    dictcompvars['WRFds'] = {'name': Vvals[0], 'standard_name': Vvals[1],             \
       'long_name': Vvals[4].replace('|',' '), 'units': Vvals[5]}
 
@@ -681 +717 @@
         ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
 
-# range_faces: LON, LAT, HGT, 'face:['WE'/'SN']:[Nptfilt]:[Nptnewrange]'
+# range_faces: LON, LAT, HGT, WRFdxdy, 'face:['WE'/'SN']:[dsfilt]:[dsnewrange]:[hvalleyrange]'
     elif diagn == 'range_faces':
             
@@ -687 +723 @@
         var1 = ncobj.variables[depvars[1]][:]
         var2 = ncobj.variables[depvars[2]][:]
-        face = depvars[3].split(':')[1]
-        Nptfilt = int(depvars[3].split(':')[2])
-        Nptnewrange = int(depvars[3].split(':')[3])
+        face = depvars[4].split(':')[1]
+        dsfilt = np.float(depvars[4].split(':')[2])
+        dsnewrange = np.float(depvars[4].split(':')[3])
+        hvalleyrange = np.float(depvars[4].split(':')[4])
 
         dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
@@ -705 +742 @@
         orogmax, ptorogmax, dhgt, peaks, valleys, origfaces, diagout, diagoutd,      \
           diagoutvd, rngorogmax, ptrngorogmax = diag.Forcompute_range_faces(lon, lat,\
-          hgt, face, Nptfilt, Nptnewrange, 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)
+          hgt, WRFds, face, dsfilt, dsnewrange, hvalleyrange, 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, 'dx', WRFdx, diagoutd, diagoutvd, newnc)
+        ncvar.insert_variable(ncobj, 'dy', WRFdy, diagoutd, diagoutvd, newnc)
+        ncvar.insert_variable(ncobj, 'ds', WRFds, diagoutd, diagoutvd, newnc)
 
         # adding variables to output file
@@ -748 +789 @@
         ncvar.set_attribute(ovar, 'deriv', axis)
 
-        ncvar.insert_variable(ncobj, 'hgtderivdeg', dhgt, diagoutd, diagoutvd,       \
-          newnc)
-        ovar = newnc.variables['orogderivdeg']
+        ncvar.insert_variable(ncobj, 'orogderiv', dhgt, diagoutd, diagoutvd, newnc)
+        ovar = newnc.variables['orogderiv']
         ncvar.set_attribute(ovar, 'deriv', axis)
 
@@ -761 +801 @@
         ovar = newnc.variables['rangefacesfilt']
         ncvar.set_attribute(ovar, 'face', face)
-        ncvar.set_attributek(ovar, 'grid_points_filter', Nptfilt, 'I')
+        ncvar.set_attributek(ovar, 'dist_filter', dsfilt, 'F')
 
         ncvar.insert_variable(ncobj, 'rangefaces', origfaces, diagoutd, diagoutvd,  \
-          newnc)
+          newnc, fill=gen.fillValueI)
         ovar = newnc.variables['rangefaces']
         ncvar.set_attribute(ovar, 'face', face)
+        ncvar.set_attributek(ovar, 'dist_newrange', dsnewrange, 'F')
+        ncvar.set_attributek(ovar, 'h_valley_newrange', hvalleyrange, 'F')
 
 # 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, Nfilt, newrange, hgtmax, pthgtmax,      &
-    derivhgt, peaks, valleys, origfaces, filtfaces, rangeshgtmax, ptrangeshgtmax)
+  SUBROUTINE compute_range_faces(d1, d2, lon, lat, hgt, dist, face, dsfilt, dsnewrange, hvalrng,      &
+    hgtmax, pthgtmax, derivhgt, peaks, valleys, origfaces, filtfaces, rangeshgtmax, ptrangeshgtmax)
 ! Subroutine to compute faces [uphill, valleys, downhill] of a mountain range along a given face
 
     IMPLICIT NONE
 
-    INTEGER, INTENT(in)                                  :: d1, d2, Nfilt, newrange
-    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: lon, lat, hgt
+    INTEGER, INTENT(in)                                  :: d1, d2
+    REAL(r_k), INTENT(in)                                :: dsfilt, dsnewrange, hvalrng
+    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: lon, lat, hgt, dist
     CHARACTER(len=*)                                     :: face
     REAL(r_k), DIMENSION(d1,d2), INTENT(out)             :: derivhgt, hgtmax, rangeshgtmax
@@ -977 +978 @@
 ! 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
-! newrange: number of grid points as distance to start a new mountain range
+! dsfilt: distance to filter orography smaller scale of it [m]
+! dsnewrange: distance to start a new mountain range [m]
+! hvalrng: maximum height of a valley to mark change of range [m]
 ! hgtmax: maximum height of the face [m]
 ! pthgtmax: grid point of the maximum height [1]
@@ -994 +996 @@
     valleys = 0
     pthgtmax = 0
+    rangeshgtmax = fillVal64
     IF (TRIM(face) == 'WE') THEN
       DO j=1, d2
-        CALL var_range_faces(d1, lon(:,j), lat(:,j), hgt(:,j), Nfilt, newrange, hgtmax1, pthgtmax1,   &
-          derivhgt(:,j), Npeaks, ipeaks1, Nvalleys, ivalleys1, origfaces(:,j), filtfaces(:,j),        &
-          Nranges, ranges1, rangeshgtmax1, irangeshgtmax1)
+        PRINT *,'Lluis:', j-1, '***'
+        CALL var_range_faces(d1, lon(:,j), lat(:,j), hgt(:,j), dist(:,j), dsfilt, dsnewrange, hvalrng,&
+          hgtmax1, pthgtmax1, derivhgt(:,j), Npeaks, ipeaks1, Nvalleys, ivalleys1, origfaces(:,j),    &
+          filtfaces(:,j), Nranges, ranges1, rangeshgtmax1, irangeshgtmax1)
         hgtmax(:,j) = hgtmax1
         pthgtmax(pthgtmax1,j) = 1
@@ -1014 +1018 @@
     ELSE IF (TRIM(face) == 'SN') THEN
       DO i=1, d1
-        CALL var_range_faces(d2, lon(i,:), lat(i,:), hgt(i,:), Nfilt, newrange, hgtmax1, pthgtmax1,   &
-          derivhgt(i,:), Npeaks, ipeaks2, Nvalleys, ivalleys2, origfaces(i,:), filtfaces(i,:),        &
-          Nranges, ranges2, rangeshgtmax2, irangeshgtmax2)
+        CALL var_range_faces(d2, lon(i,:), lat(i,:), hgt(i,:), dist(i,:), dsfilt, dsnewrange, hvalrng,&
+          hgtmax1, pthgtmax1, derivhgt(i,:), Npeaks, ipeaks2, Nvalleys, ivalleys2, origfaces(i,:),    &
+          filtfaces(i,:), Nranges, ranges2, rangeshgtmax2, irangeshgtmax2)
         hgtmax(i,:) = hgtmax1
         pthgtmax(i,pthgtmax1) = 1

trunk/tools/module_ForDiagnosticsVars.f90

@@ -1580 +1580 @@
   END SUBROUTINE var_fog_FRAML50
 
-  SUBROUTINE var_range_faces(d, lon, lat, hgt, filt, newrange, hgtmax, ihgtmax, dhgt, Npeaks, ipeaks, &
-    Nvalleys, ivalleys, faces0, faces, Nranges, ranges, rangeshgtmax, irangeshgtmax)
+  SUBROUTINE var_range_faces(d, lon, lat, hgt, dist, filt, newrange, hvalleyrange, hgtmax, ihgtmax,   &
+    dhgt, Npeaks, ipeaks, Nvalleys, ivalleys, faces0, faces, Nranges, ranges, rangeshgtmax,           &
+    irangeshgtmax)
 ! Subroutine to compute faces [uphill, valleys, downhill] of a monuntain range along a face
 
     IMPLICIT NONE
 
-    INTEGER, INTENT(in)                                  :: d, filt, newrange
-    REAL(r_k), DIMENSION(d), INTENT(in)                  :: lon, lat, hgt
+    INTEGER, INTENT(in)                                  :: d
+    REAL(r_k), INTENT(in)                                :: filt, newrange, hvalleyrange
+    REAL(r_k), DIMENSION(d), INTENT(in)                  :: lon, lat, hgt, dist
     INTEGER, INTENT(out)                                 :: ihgtmax, Npeaks, Nvalleys, Nranges
     INTEGER, DIMENSION(d), INTENT(out)                   :: ipeaks, ivalleys, faces0, faces,          &
@@ -1600 +1602 @@
     INTEGER                                              :: Nfaces, Nfaces1, Npeaks1, Nvalleys1
     INTEGER                                              :: fbeg, fend
+    REAL(r_k)                                            :: dd
     INTEGER, DIMENSION(1)                                :: ihmax
     INTEGER, DIMENSION(d)                                :: ddhgt, Ndhgt
@@ -1606 +1609 @@
     REAL(r_k), DIMENSION(d)                              :: sortedpeaks
     INTEGER, DIMENSION(d)                                :: isortedpeaks
-    LOGICAL                                              :: rangewithin, peakwithin, valleywithin
+    LOGICAL                                              :: firstvalley, rangewithin, peakwithin,      &
+      valleywithin
 
 !!!!!!! Variables
@@ -1612 +1616 @@
 ! lat: latitude [degrees north]
 ! hgt: topograpical height [m]
-! filt: number of grid points to use to filter the topographic values. used to define:
+! dist: distance between grid points [m]
+! filt: distance to use to filter the topographic values [m]. used to define:
 !   the minimum length of a given section 
-!   running mean filter
-! newrange: number of grid points as distance to start a new mountain range
+! newrange: distance to start a new mountain range [m]
+! hvalleyrange: maximum height of a valley to mark change of range [m]
 
     fname = 'var_range_faces'
@@ -1632 +1637 @@
 
     ! range slope
-    dLl = SQRT((lon(2:d)-lon(1:d-1))**2+(lat(2:d)-lat(1:d-1))**2)
     dhgt(1:d) = hgt(2:d) - hgt(1:d-1)
-    dhgt = dhgt/dLl
+    dhgt = dhgt/dist(1:d)
 
     !PRINT *, 'slope:', dhgt
@@ -1643 +1647 @@
     DO i=1, d-1
       IF (dhgt(i) > 0.) THEN
-        faces0(i) = 1
+        IF (hgt(i) < hvalleyrange) THEN
+          faces0(i) = fillValI
+        ELSE
+          faces0(i) = 1
+        END IF
       ELSE IF (dhgt(i) < 0.) THEN
-        faces0(i) = -1
+        IF (hgt(i) < hvalleyrange) THEN
+          faces0(i) = fillValI
+        ELSE
+          faces0(i) = -1
+        END IF
       ELSE
-        faces0(i) = 0
+        IF (hgt(i) == zeroRK .AND. hgt(i+1) == zeroRK) THEN
+          faces0(i) = fillValI
+        ELSE
+          faces0(i) = 0
+        END IF
       END IF
       ! peaks
@@ -1686 +1702 @@
     ! Defining quantitiy of ranges within the face 
     !   sorting peaks within the face and defining ranges as maximum peaks distanced > newrage
-    sortedpeaks = peaks
-    CALL SortR_K(sortedpeaks, Npeaks)
-    isortedpeaks = 0
-    DO i=1, Npeaks
-      DO j=1, Npeaks
-        IF (peaks(j) == sortedpeaks(i)) isortedpeaks(i) = ipeaks(j)
+
+    !sortedpeaks = peaks
+    !CALL SortR_K(sortedpeaks, Npeaks)
+    !isortedpeaks = 0
+    !DO i=1, Npeaks
+    !  DO j=1, Npeaks
+    !    IF (peaks(j) == sortedpeaks(i)) isortedpeaks(i) = ipeaks(j)
+    !  END DO
+    !END DO
+ 
+    ranges = -1
+    Nranges = 0
+    l = 0
+    DO i=1,d
+      IF (hgt(i) >= hvalleyrange) THEN
+        IF (Nranges == 0) THEN
+          Nranges = 1
+          ranges(1,Nranges) = i
+        ELSE
+          IF (ranges(2,Nranges) /= -1) THEN
+            Nranges = Nranges + 1
+            ranges(1,Nranges) = i
+          END IF
+        END IF
+      ELSE
+        IF ((ranges(1,Nranges) /= -1) .AND. (ranges(2,Nranges) == -1)) ranges(2,Nranges) = i-1
+      END IF
+    END DO
+    IF (ranges(2,Nranges) == -1) ranges(2,Nranges) = d
+
+    ! Getting characteristics of ranges
+    irangeshgtmax = -1
+    rangeshgtmax = -1000.
+    k = 1
+    DO i=1, Nranges
+      DO j=ranges(1,i), ranges(2,i)
+       IF (hgt(j) > rangeshgtmax(i)) THEN
+         irangeshgtmax(i) = j
+         rangeshgtmax(i) = hgt(j)
+        END IF
       END DO
     END DO
- 
-    ranges = -1
-    Nranges = 1
-    ranges(1,Nranges) = ipeaks(1)
-    ranges(2,Nranges) = ipeaks(1)
-    irangeshgtmax(Nranges) = ipeaks(1)
-    rangeshgtmax(Nranges) = peaks(1)
-    IF (Npeaks >= 2) THEN
-      DO j=2, Npeaks
-        rangewithin = .FALSE.
-        IF (ranges(1,Nranges) - isortedpeaks(j) < newrange) THEN
-          ranges(2,Nranges) = ipeaks(j)
-          IF (peaks(j) > rangeshgtmax(Nranges)) THEN
-            rangeshgtmax(Nranges) = peaks(j)
-            irangeshgtmax(Nranges) = ipeaks(j)
-          END IF
-        ELSE
-          Nranges = Nranges + 1
-          ranges(1,Nranges) = ipeaks(j)
-          ranges(2,Nranges) = ipeaks(j)
-          rangeshgtmax(Nranges) = peaks(j)
-          irangeshgtmax(Nranges) = ipeaks(j)
-        END IF
-      END DO
-    END IF
     !PRINT *, Nranges, ' _______'
-    !DO i=, Nranges
-    !   PRINT *,i,':', ranges(:,i),' max:', rangeshgtmax(i), ' ,', irangeshgtmax(i)
+    !DO i=1, Nranges
+    !  PRINT *,i,':', ranges(:,i),' max:', rangeshgtmax(i), ' ,', irangeshgtmax(i)
     !END DO
 
@@ -1743 +1768 @@
     END IF
 
+    ! Using defined ranges for the last refinement. 
+    !   Uphills after range maximum peak as valley: if dhgt[i] > 0. and i > rangemax --> face[i] = 0
+    !   Downhills before range maximum peak as valley: if dhgt[i] < 0. and i < rangemax --> face[i] = 0
+    !   face values before first range, removed: if dhgt[i] != 0. and i < range[1,1] --> face[i] = fillvalueI
+    !   face values after last range, removed: if dhgt[i] != 0. and i < range[2,Nranges] --> face[i] = fillvalueI
+    DO i=1,Nranges
+      DO j=ranges(1,i), ranges(2,i)
+        IF (dhgt(j) > 0. .AND. j > irangeshgtmax(i)) faces0(j) = 0
+        IF (dhgt(j) < 0. .AND. j < irangeshgtmax(i)) faces0(j) = 0
+      END DO
+    END DO
+    IF (dhgt(j) /= 0. .AND. j < ranges(1,1)) faces0(j) = fillValI
+    IF (dhgt(j) /= 0. .AND. j > ranges(2,Nranges)) faces0(j) = fillValI
+
     ! classification using filt. 
     !   On that section of the hill smaller than the selected filter. Replace values by the ones from 
@@ -1749 +1788 @@
     !       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
+    !!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
+
+    !!  dd = SUM(dist(fbeg:fend))
+    !!  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
@@ -1805 +1845 @@
     END IF
 
+
+    ! Removing above sea points
+    DO i=1, d-1
+      IF (hgt(i) == zeroRK .AND. hgt(i+1) == zeroRK) dhgt(i) = fillVal64
+    END DO
+
     !PRINT *, '# lon[0] lat[1] heights[2] slope[3] faces0[4] faces1[5] faces1[6] ddhgt[7]'
     !DO i=1, d

trunk/tools/variables_values.dat

@@ -123 +123 @@
 dqwak, dqwak, wake_water_vapor_tendency, -3.e-7, 3.e-7, wake|water|vapor|tendency, kg/kg/s, seismic, $\partial_{t}{qv}^{wak}$, dtqvwak
 LDQWAK, dqwak, wake_water_vapor_tendency, -3.e-7, 3.e-7, wake|water|vapor|tendency, kg/kg/s, seismic, $\partial_{t}{qv}^{wak}$, dtqvwak
+ds, ds, horizontal_grid_distance, 0., 200., total|distance|between|grid|points, m, Reds, $ds$, ds
+WRFds, ds, horizontal_grid_distance, 0., 200., total|distance|between|grid|points, m, Reds, $ds$, ds
 dta, dta, tendency_air_temperature, -3.e-3, 3.e-3, tendency|of|air|temperature, K/s, seismic, $\partial_{t}{ta}$, dtta
 tnt, dta, tendency_air_temperature, -3.e-3, 3.e-3, tendency|of|air|temperature, K/s, seismic, $\partial_{t}{ta}$, dtta
@@ -173 +175 @@
 dvvdf, dvvdf, vertical_difussion_northward_wind_tendency, -3.e-3, 3.e-3, vertical|difussion|northward|wind|tendency, ms-2, seismic, $\partial_{t}{va}^{zdiff}$, dtvazdiff
 LDVVDF, dvvdf, vertical_difussion_northward_wind_tendency, -3.e-3, 3.e-3, vertical|difussion|northward|wind|tendency, ms-2, seismic, $\partial_{t}{va}^{zdiff}$, dtvazdiff
+dx, dx, eastward_horizontal_grid_distance, 0., 200., eastward|distance|between|grid|points, m, Reds, $dx$, dx
+WRFdx, dx, eastward_horizontal_grid_distance, 0., 200., eastward|distance|between|grid|points, m, Reds, $dx$, dx
+dy, dy, northward_horizontal_grid_distance, 0., 200., northward|distance|between|grid|points, m, Blues, $dy$, dy
+WRFdy, dy, northward_horizontal_grid_distance, 0., 200., northward|distance|between|grid|points, m, Blues, $dy$, dy
 etau, etau, etau, 0., 1, eta|values|on|half|(mass)|levels, -, reds, $\eta_{hm}$, eta_hm
 ZNU, etau, etau, 0., 1, eta|values|on|half|(mass)|levels, -, reds, $\eta_{hm}$, eta_hm
@@ -296 +302 @@
 HGT, orog, orography,  0., 3000., surface|altitude, m,terrain, $orog$, orog
 HGT_M, orog, orography,  0., 3000., surface|altitude, m,terrain, $orog$, orog
+orogderiv, orogderiv, orography_derivative_degrees,  0., 3000., derivative|along|degree|coordinate|of|surface|altitude, m m-1,seismic, $orog_{deri}$, orog_deriv
+hgtderiv, orogderiv, orography_derivative_degrees,  0., 3000., derivative|along|degree|coordinate|of|surface|altitude, m m-1,seismic, $orog_{deri}$, orog_deriv
 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