Changeset 2213 in lmdz_wrf

Timestamp:

Nov 5, 2018, 3:57:05 PM (6 years ago)

Author:

lfita

Message:

Improving/Enhancing? the range_faces

Location:

trunk/tools

Files:

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

trunk/tools/diag_tools.py

@@ -1370 +1370 @@
     return var1, var2, vardims, varvdims
 
-def Forcompute_range_faces(lon, lat, hgt, face, Nfilt, dimns, dimvns):
+def Forcompute_range_faces(lon, lat, hgt, face, Nfilt, newrng, dimns, dimvns):
     """ Function to compute faces [uphill, valley, downhill] of sections of a mountain
       rage, along a given face
@@ -1379 +1379 @@
       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
       [dimns]= list of the name of the dimensions of [smois]
       [dimvns]= list of the name of the variables with the values of the 
@@ -1393 +1394 @@
         dy = hgt.shape[0]
 
-        dhgtt, peakst, valleyst, ofacest, ffacest =                                  \
+        hgtmaxt, pthgtmaxt, 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)
+          newrange=newrng, d1=dx, d2=dy)
+        hgtmax = hgtmaxt.transpose()
+        pthgtmax = pthgtmaxt.transpose()
         dhgt = dhgtt.transpose()
         peaks = peakst.transpose()
@@ -1408 +1411 @@
         quit(-1)
 
-    return dhgt, peaks, valleys, origfaces, filtfaces, vardims, varvdims
+    return hgtmax, pthgtmax, dhgt, peaks, valleys, origfaces, filtfaces, vardims,    \
+      varvdims
 
 ####### ###### ##### #### ### ## # END Fortran diagnostics

trunk/tools/diagnostics.py

@@ -681 +681 @@
         ncvar.insert_variable(ncobj, 'hurs', diagout, diagoutd, diagoutvd, newnc)
 
-# range_faces: LON, LAT, HGT, 'face:['WE'/'SN']:[Nptfilt]'
+# range_faces: LON, LAT, HGT, 'face:['WE'/'SN']:[Nptfilt]:[Nptnewrange]'
     elif diagn == 'range_faces':
             
@@ -688 +688 @@
         var2 = ncobj.variables[depvars[2]][:]
         face = depvars[3].split(':')[1]
-        Nptfilt = depvars[3].split(':')[2]
+        Nptfilt = int(depvars[3].split(':')[2])
+        Nptnewrange = int(depvars[3].split(':')[3])
 
         dnamesvar = list(ncobj.variables[depvars[0]].dimensions)
@@ -702 +703 @@
             hgt = var2[:]
 
-        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)
+        orogmax, ptorogmax, dhgt, peaks, valleys, origfaces, diagout, diagoutd,      \
+          diagoutvd = 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)
+
+        # adding variables to output file
+        if face == 'WE': axis = 'lon'
+        elif face == 'SN': axis = 'lat'
+
+        ncvar.insert_variable(ncobj, 'orogmax', orogmax, diagoutd, diagoutvd,        \
+          newnc)
+        ovar = newnc.variables['orogmax']
+        ncvar.set_attribute(ovar, 'deriv', axis)
+
+        ncvar.insert_variable(ncobj, 'ptorogmax', ptorogmax, diagoutd, diagoutvd,    \
+          newnc)
+        ovar = newnc.variables['ptorogmax']
+        ncvar.set_attribute(ovar, 'deriv', axis)
 
         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)
 

trunk/tools/module_ForDiagnostics.f90

@@ -949 +949 @@
   END SUBROUTINE compute_fog_FRAML50
 
-  SUBROUTINE compute_range_faces(d1, d2, lon, lat, hgt, face, Nfilt, derivhgt, peaks, valleys,        &
-    origfaces, filtfaces)
+  SUBROUTINE compute_range_faces(d1, d2, lon, lat, hgt, face, Nfilt, newrange, hgtmax, pthgtmax,      &
+    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, Nfilt
+    INTEGER, INTENT(in)                                  :: d1, d2, Nfilt, newrange
     REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: lon, lat, hgt
     CHARACTER(len=*)                                     :: face
-    REAL(r_k), DIMENSION(d1,d2), INTENT(out)             :: derivhgt
-    INTEGER, DIMENSION(d1,d2), INTENT(out)               :: origfaces, filtfaces, peaks, valleys
+    REAL(r_k), DIMENSION(d1,d2), INTENT(out)             :: derivhgt, hgtmax
+    INTEGER, DIMENSION(d1,d2), INTENT(out)               :: pthgtmax, origfaces, filtfaces, peaks,    &
+      valleys
  
 ! Local
     INTEGER                                              :: i, j
-    INTEGER                                              :: Npeaks, Nvalleys
+    INTEGER                                              :: pthgtmax1, Npeaks, Nvalleys
+    REAL(r_k)                                            :: hgtmax1
     INTEGER, DIMENSION(d1)                               :: ipeaks1, ivalleys1
     INTEGER, DIMENSION(d2)                               :: ipeaks2, ivalleys2
@@ -973 +975 @@
 ! 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
+! hgtmax: maximum height of the face [m]
+! pthgtmax: grid point of the maximum height [1]
 ! derivhgt: topograpic derivative along axis [m deg-1]
 ! peaks: peak point
@@ -983 +988 @@
     peaks = 0
     valleys = 0
+    pthgtmax = 0
     IF (TRIM(face) == 'WE') THEN
       DO j=1, d2
-        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))
+        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))
+        hgtmax(:,j) = hgtmax1
+        pthgtmax(pthgtmax1,j) = 1
         DO i=1, Npeaks
           peaks(ipeaks1(i),j) = 1
@@ -997 +1004 @@
     ELSE IF (TRIM(face) == 'SN') THEN
       DO i=1, d1
-        CALL var_range_faces(d2, lon(i,:), lat(i,:), hgt(i,:), Nfilt, derivhgt(i,:), Npeaks,          &
-          ipeaks2, Nvalleys, ivalleys2, origfaces(i,:), filtfaces(i,:))
+        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,:))
+        hgtmax(i,:) = hgtmax1
+        pthgtmax(i,pthgtmax1) = 1
         DO j=1, Npeaks
           peaks(i,ipeaks2(j)) = 1

trunk/tools/module_ForDiagnosticsVars.f90

@@ -1580 +1580 @@
   END SUBROUTINE var_fog_FRAML50
 
-  SUBROUTINE var_range_faces(d, lon, lat, hgt, filt, dhgt, Npeaks, ipeaks, Nvalleys, ivalleys,        &
-    faces0, faces)
+  SUBROUTINE var_range_faces(d, lon, lat, hgt, filt, newrange, hgtmax, ihgtmax, 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, filt
+    INTEGER, INTENT(in)                                  :: d, filt, newrange
     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
+    INTEGER, INTENT(out)                                 :: ihgtmax
+    REAL(r_k), INTENT(out)                               :: hgtmax
+    REAL(r_k), DIMENSION(d), INTENT(out)                 :: dhgt, rangeshgtmax
+    INTEGER, DIMENSION(d), INTENT(out)                   :: ipeaks, ivalleys, faces0, faces,          &
+      irangeshgtmax
+    INTEGER, INTENT(out)                                 :: Npeaks, Nvalleys, Nranges
 
 ! Local
-    INTEGER                                              :: i, j, k, l, m 
-    INTEGER                                              :: iface
+    INTEGER                                              :: i, j, j1, k, l, m 
+    INTEGER                                              :: iface, Nranges
     INTEGER                                              :: Nfaces, Nfaces1, Npeaks1, Nvalleys1
     INTEGER                                              :: fbeg, fend
     INTEGER, DIMENSION(1)                                :: ihmax
-    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
+    REAL(r_k), DIMENSION(d)                              :: dLl, peaks, valleys, peaks1, valleys1
     LOGICAL                                              :: peakwithin, valleywithin
 
@@ -1609 +1611 @@
 ! 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
+!   the minimum length of a given section 
 !   running mean filter
+! newrange: number of grid points as distance to start a new mountain range
 
     fname = 'var_range_faces'
 
-    PRINT *, 'heights:', hgt
+    !PRINT *, 'heights:', hgt
 
     ! Looking for the maximum height
     hgtmax = MAXVAL(hgt)
     ihmax = MAXLOC(hgt)
-
-    PRINT *, 'height max:', hgtmax, 'location:', ihmax
+    ihgtmax = ihmax(1)
+
+    !PRINT *, 'height max:', hgtmax, 'location:', ihmax
 
     ! Filtering height [NOT needed]
@@ -1626 +1630 @@
 
     ! 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)
-
-    PRINT *, 'slope:', dhgt
+    dhgt = dhgt/dLl
+
+    !PRINT *, 'slope:', dhgt
 
     ! Initial classification
@@ -1655 +1661 @@
     END DO
 
-    PRINT *, 'faces:', faces0
-    PRINT *, Npeaks, ' peaks:', peaks(1:Npeaks), ' ipeak:', ipeaks(1:Npeaks)
-    PRINT *, Nvalleys, ' valleys:', valleys(1:Nvalleys), ' ivalley:', ivalleys(1:Nvalleys)
+    !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
@@ -1674 +1680 @@
     END DO
 
-    PRINT *, Nfaces, ' length faces:', Ndhgt(1:Nfaces)
+    !PRINT *, Nfaces, ' length faces:', Ndhgt(1:Nfaces)
+
+    ! Defining quantitiy of ranges within the face
+    Nranges = 1
+    IF (Npeaks > 2) THEN
+      DO i=1, Npeaks-1
+        IF (ipeaks(j+1)-ipeaks(j) < newrange)) THEN
+         AQUI
+        END IF
+      END DO
+    END IF
+
+    ! Defining valleys as that consecutive grid points below surroungind peaks from the same range
+    IF (Npeaks > 1) THEN
+      j = 1
+      j1 = 2
+      DO i=2, d
+        IF (i >= ipeaks(j) .AND. i < ipeaks(j+1) .AND. (ipeaks(j+1)-ipeaks(j) < newrange)) THEN
+          IF (hgt(i) < peaks(j) .AND. hgt(i) < peaks(j+1)) faces0(i) = 0
+          IF (i == ipeaks(j) .AND. hgt(i) < hgtmax) faces0(i) = 0
+        ELSE IF (i == ipeaks(j+1)) THEN
+          j = j1
+          j1 = j1 + 1
+          IF (j1 > Npeaks) EXIT
+        END IF
+      END DO
+    END IF
 
     ! classification using filt. 
@@ -1738 +1770 @@
     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), faces0(i), faces1(i), faces(i), ddhgt(i)
-    END DO
+    !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), faces0(i), faces1(i), faces(i), ddhgt(i)
+    !END DO
 
     RETURN

trunk/tools/variables_values.dat

@@ -387 +387 @@
 slp, psl, air_pressure_at_sea_level, 85000., 104000., sea|level|pressure, Pa, Greens, $psl$, psl
 WRFmslp, psl, air_pressure_at_sea_level, 85000., 104000., mean|sea|level|pressure, Pa, Greens, $psl$, psl
+ptorog, ptorog, orography_point,  0., 3000., surface|altitude|point, m, terrain, $ptorog$, ptorog
 qth, qth, thermal_plume_total_water_content, 0., 25., total|water|cotent|in|thermal|plume, mm, YlOrRd, $prw^{th}_{plume}$, prwth_plume
 q_th, qth, thermal_plume_total_water_content, 0., 25., total|water|cotent|in|thermal|plume, mm, YlOrRd, $prw^{th}_{plume}$, prwth_plume