CurvatureMicroMet Subroutine

SUBROUTINE CurvatureMicroMet &
!
(dem, curve_len_scale, curvature)

IMPLICIT NONE

!arguments with intent in
TYPE(grid_real),    INTENT(in) :: dem
REAL (KIND = float), INTENT(in) :: curve_len_scale

!arguments with intent out
TYPE (grid_real), INTENT (out) :: curvature

!local declarations:
 REAL (KIND = float) :: deltaxy
 REAL (KIND = float) :: z, zw, ze, zs, zn, zsw, zne, znw, zse
 REAL (KIND = float) :: curve_max
 INTEGER (KIND = short) :: inc
 INTEGER (KIND = short) :: i, j


!-----------------------------end of declarations------------------------------
!allocate grid
CALL NewGrid (curvature, dem)


!compute curvature
DO i = 1, dem % idim
  DO j = 1, dem % jdim
   IF (dem % mat (i,j) /= dem % nodata) THEN
   
     !average grid increment
     deltaxy = CellArea (dem,i,j) ** 0.5

     !Convert the length scale to an appropriate grid increment.
     inc = max(1,nint(curve_len_scale/deltaxy))
     
     !srt z
     z = dem % mat(i,j)
     
     !set zn
     IF ( .NOT. IsOutOfGrid (i-inc,j,dem) ) THEN
        IF ( dem % mat (i-inc,j) /= dem % nodata) THEN
           zn = dem % mat (i-inc,j)
        ELSE
           zn = dem % mat (i,j)
        END IF  
      ELSE 
           zn = dem % mat (i,j)
      END IF
      
      !set zne
      IF ( .NOT. IsOutOfGrid (i-inc,j+inc,dem) ) THEN
         IF ( dem % mat (i-inc,j+inc) /= dem % nodata) THEN
            zne = dem % mat (i-inc,j+inc)
         ELSE
            zne = dem % mat (i,j)
         END IF  
       ELSE 
            zne = dem % mat (i,j)
       END IF
       
       !set ze
       IF ( .NOT. IsOutOfGrid (i,j+inc,dem) ) THEN
          IF ( dem % mat (i,j+inc) /= dem % nodata) THEN
             ze = dem % mat (i,j+inc)
          ELSE
             ze = dem % mat (i,j)
          END IF  
        ELSE 
             ze = dem % mat (i,j)
        END IF
        
        !set zse
        IF ( .NOT. IsOutOfGrid (i+inc,j+inc,dem) ) THEN
           IF ( dem % mat (i+inc,j+inc) /= dem % nodata) THEN
              zse = dem % mat (i+inc,j+inc)
           ELSE
              zse = dem % mat (i,j)
           END IF  
         ELSE 
              zse = dem % mat (i,j)
         END IF
         
         !set zs
         IF ( .NOT. IsOutOfGrid (i+inc,j,dem) ) THEN
            IF ( dem % mat (i+inc,j) /= dem % nodata) THEN
               zs = dem % mat (i+inc,j)
            ELSE
               zs = dem % mat (i,j)
            END IF  
          ELSE 
               zs = dem % mat (i,j)
          END IF
          
          !set zsw
          IF ( .NOT. IsOutOfGrid (i+inc,j-inc,dem) ) THEN
             IF ( dem % mat (i+inc,j-inc) /= dem % nodata) THEN
                zsw = dem % mat (i+inc,j-inc)
             ELSE
                zsw = dem % mat (i,j)
             END IF  
           ELSE 
                zsw = dem % mat (i,j)
           END IF
           
           !set zw
           IF ( .NOT. IsOutOfGrid (i,j-inc,dem) ) THEN
              IF ( dem % mat (i,j-inc) /= dem % nodata) THEN
                 zw = dem % mat (i,j-inc)
              ELSE
                 zw = dem % mat (i,j)
              END IF  
            ELSE 
                 zw = dem % mat (i,j)
            END IF
            
            !set znw
            IF ( .NOT. IsOutOfGrid (i-inc,j-inc,dem) ) THEN
               IF ( dem % mat (i-inc,j-inc) /= dem % nodata) THEN
                  znw = dem % mat (i-inc,j-inc)
               ELSE
                  znw = dem % mat (i,j)
               END IF  
             ELSE 
                  znw = dem % mat (i,j)
             END IF
             
             
             !curvature
              curvature % mat(i,j) = (4.0 * z - znw - zse - zsw - zne ) / &
                   (sqrt(2.0) * 16.0 * real(inc) * deltaxy) + &
                   (4.0 * z - zs - zn - ze - zw) / &
                   (16.0 * real(inc) * deltaxy)
      
   END IF
  END DO
END DO


! Scale the curvature such that the max abs(curvature) has a value
!   of abs(0.5).  Include a 1 mm curvature in curve_max to prevent
!   divisions by zero in flat terrain where the curvature is zero.

curve_max = 0.0 + 0.001

DO i = 1, dem % idim
  DO j = 1, dem % jdim
    IF (dem % mat (i,j) /= dem % nodata) THEN 
       curve_max = MAX(curve_max,ABS(curvature%mat(i,j)))
    END IF
  END DO
END DO

DO i = 1, dem % idim
  DO j = 1, dem % jdim
    IF (dem % mat (i,j) /= dem % nodata) THEN 
       curvature%mat(i,j) = curvature%mat(i,j) / (2.0 * curve_max)
    END IF
  END DO
END DO
       

RETURN

END SUBROUTINE CurvatureMicroMet