WindMicromet Subroutine

SUBROUTINE WindMicromet &
!
(speed, dir, slope, curvature, slope_az, slopewt, curvewt, grid, winddir)

USE Units, ONLY: &
!Imported parameters
degToRad, radToDeg, pi

IMPLICIT NONE

!Arguments with intent (in):
TYPE (ObservationalNetwork), INTENT(IN) :: speed !windspeed observations at stations 
TYPE (ObservationalNetwork), INTENT(IN) :: dir !wind direction observations at stations 
TYPE (grid_real), INTENT(IN) :: slope !terrain slope grid
TYPE (grid_real), INTENT(IN) :: curvature !terrain curvature grid
TYPE (grid_real), INTENT(IN) :: slope_az !terrain slope azimuth
REAL (KIND = float),  INTENT(IN) :: slopewt !slope weighting factor
REAL (KIND = float),  INTENT(IN) :: curvewt !curvature weighting factor

!Arguments with intent (inout):
TYPE (grid_real), INTENT (INOUT) :: grid !wind speed grid [m/s]
TYPE (grid_real), OPTIONAL, INTENT (INOUT) :: winddir !wind direction [deg]


!Local declarations
TYPE (ObservationalNetwork) :: temp_speed, temp_dir, active_speed, active_dir
TYPE (ObservationalNetwork) :: uwind, vwind
INTEGER :: count_speed, count_dir, i, j
TYPE (grid_real) :: uwind_grid !zonal wind [m/s]
TYPE (grid_real) :: vwind_grid !meridional wind [m/s]
TYPE (grid_real) :: winddir_grid !wind direction [rad]
TYPE (grid_real) :: wind_slope !slope in the direction of the wind [rad]
REAL (KIND = float) :: wslope_max ![rad]
REAL (KIND = float) :: windwt !wind weighting factor
REAL (KIND = float) :: thetad !diverting factor
REAL (KIND = float),  PARAMETER :: windspd_min = 0.00001 ![m/s]


!----------------end of declarations-------------------------------------------

!check for available measurements. Both speed and direction must exist
!first check nodata. If only speed or direction are missing, set both to missing 
! assume dir and speed sations are in memory in the same order
CALL CopyNetwork (speed,temp_speed)
CALL CopyNetwork (dir,temp_dir)

DO i = 1, temp_speed %countObs
  IF (temp_speed % obs (i) % value == temp_speed % nodata) THEN
      !set dir to nodata
      temp_dir % obs (i) % value = temp_dir % nodata
  ELSE IF (temp_dir % obs (i) % value == temp_dir % nodata) THEN
      !set speed to nodata
      temp_speed % obs (i) % value = temp_speed % nodata
  END IF
END DO

!now remove nodata
CALL ActualObservations (temp_speed, count_speed, active_speed)
CALL ActualObservations (temp_dir, count_dir, active_dir)

!initialize zonal and meridional component network
CALL CopyNetwork (active_speed, uwind)
CALL CopyNetwork (active_speed, vwind)

!initialize grid
CALL NewGrid (uwind_grid, grid)
CALL NewGrid (vwind_grid, grid)
CALL NewGrid (winddir_grid, grid)
CALL NewGrid (wind_slope, grid)

!compute u and v at stations
DO i = 1, active_speed % countObs
  uwind % obs (i) % value = - active_speed % obs (i) % value * &
                          SIN ( active_dir % obs (i) % value * degToRad)
  vwind % obs (i) % value = - active_speed % obs (i) % value * &
                          COS ( active_dir % obs (i) % value * degToRad)
END DO

!interpolate u and v grid 
CALL IDW (network=uwind, grid=uwind_grid, method=1, power=2.0, near=5)
!CALL BarnesOI (network=uwind, grid=uwind_grid, gammapar = 0.6)
     
CALL IDW (network=vwind, grid=vwind_grid, method=1, power=2.0, near=5)
!CALL BarnesOI (network=vwind, grid=vwind_grid, gammapar = 0.6)
   
!compute wind direction and windspeed grid 
DO i = 1, grid % idim
  DO j = 1, grid % jdim
      IF (grid % mat (i,j) /= grid % nodata) THEN
           
        ! Some compilers do not allow both u and v to be 0.0 in
        !   the atan2 computation.
        IF (ABS(uwind_grid % mat (i,j)) < 1e-10) &
              uwind_grid % mat (i,j) = 1e-10
           
        winddir_grid %mat (i,j) = &
            ATAN2(uwind_grid % mat (i,j),vwind_grid % mat (i,j))
  
        IF (winddir_grid %mat (i,j) >= pi) THEN
          winddir_grid % mat (i,j) = winddir_grid % mat (i,j) - pi
        ELSE
          winddir_grid % mat (i,j) = winddir_grid % mat (i,j) + pi
        END IF
           
        !windspeed
        grid % mat (i,j) = &
            SQRT(uwind_grid % mat (i,j)**2 + vwind_grid % mat (i,j)**2)
            
      END IF
  END DO
END DO

! Compute the slope in the direction of the wind.
DO i = 1, grid % idim
  DO j = 1, grid % jdim
      IF (grid % mat (i,j) /= grid % nodata) THEN
     
      wind_slope % mat (i,j) = slope % mat (i,j) * &
            COS ( winddir_grid % mat (i,j) - slope_az % mat (i,j) )
         
      END IF
  END DO
END DO

! Scale the wind slope such that the max abs(wind slope) has a value
! of abs(0.5).  Include a 1 mm slope in slope_max to prevent
! divisions by zero in flat terrain where the slope is zero.
wslope_max = 0.0 + 0.001
DO i = 1, grid % idim
  DO j = 1, grid % jdim
    IF (grid % mat (i,j) /= grid % nodata) THEN
        wslope_max = MAX (wslope_max,ABS(wind_slope % mat (i,j)))
    END IF
  END DO
END DO
      
DO i = 1, grid % idim
  DO j = 1, grid % jdim
    IF (grid % mat (i,j) /= grid % nodata) THEN
        wind_slope % mat (i,j) = wind_slope % mat (i,j) / (2.0 * wslope_max)
    END IF
  END DO
END DO

! Calculate the wind speed and direction adjustments.  The
! curvature and wind_slope values range between -0.5 and +0.5.
! Valid slopewt and curvewt values are between 0 and 1, with
! values of 0.5 giving approximately equal weight to slope and
! curvature.  I suggest that slopewt and curvewt be set such
! that slopewt + curvewt = 1.0.  This will limit the total
! wind weight to between 0.5 and 1.5 (but this is not required).
DO i = 1, grid % idim
  DO j = 1, grid % jdim
    IF (grid % mat (i,j) /= grid % nodata) THEN
            
      ! Compute the wind weighting factor.
      windwt = 1.0 + slopewt * wind_slope % mat (i,j) + &
      curvewt * curvature % mat (i,j)
      !if (windwt>0) write(*,*) windwt,  wind_slope % mat (i,j), curvature % mat (i,j)
      !pause
            
      ! Generate the terrain-modified wind speed.
      grid % mat (i,j) = windwt * grid % mat (i,j)
      
      !compute diverting factor
      thetad = - 0.5 * wind_slope % mat (i,j) * &
          SIN ( 2. * (slope_az % mat (i,j) - winddir_grid % mat (i,j)) )
      
      IF (PRESENT (winddir)) THEN
      
         winddir % mat (i,j) = ( winddir_grid % mat (i,j) + thetad ) * radToDeg
         
      END IF   
              
    END IF
  END DO
END DO

!deallocate memory
CALL DestroyNetwork (temp_speed)
CALL DestroyNetwork (temp_dir)
CALL DestroyNetwork (active_speed)
CALL DestroyNetwork (active_dir)
CALL DestroyNetwork (uwind)
CALL DestroyNetwork (vwind)

CALL GridDestroy (uwind_grid)
CALL GridDestroy (vwind_grid)
CALL GridDestroy (winddir_grid)
CALL GridDestroy (wind_slope)

RETURN
END SUBROUTINE WindMicromet