util/src/elev_angle.F

   1 program elev_angle
   2
   3    implicit none
   4
   5    integer, external :: iargc
   6
   7    integer :: istatus
   8    integer :: i, j, n_bins, we_dim, sn_dim
   9    real, pointer, dimension(:,:) :: hgt
  10    real, allocatable, dimension(:,:,:) :: topo_angle
  11    character (len=1024) :: filename
  12
  13    if (iargc() /= 1) then
  14       write(6,*) ' '
  15       write(6,*) 'Usage: elev_angle.exe <geogrid output file>'
  16       write(6,*) ' '
  17       stop
  18    end if
  19
  20    call getarg(1,filename)
  21
  22    !
  23    ! Read in topography field from geogrid output file
  24    !
  25    call read_topo_field(filename, hgt, we_dim, sn_dim, istatus)
  26    if (istatus /= 0) stop
  27
  28    write(6,*) 'Read HGT_M field dimensioned ',we_dim,sn_dim
  29    do j=1,sn_dim,10
  30    do i=1,we_dim,10
  31       write(6,'(a6,i3,a1,i3,a2,f13.5)') 'HGT_M(',i,',',j,')=',hgt(i,j)
  32    end do
  33    end do
  34
  35    n_bins = 180
  36    allocate(topo_angle(we_dim, sn_dim, n_bins))
  37
  38    !
  39    ! Compute elevation angles for each azimuth angle bin
  40    !
  41    topo_angle = 10.0
  42
  43    !
  44    ! Write elevation angle data to intermediate file
  45    !
  46    call write_elev_angles(topo_angle, we_dim, sn_dim, n_bins, istatus)
  47
  48
  49    deallocate(topo_angle)
  50    deallocate(hgt)
  51
  52    stop
  53
  54 end program elev_angle
  55
  56
  57 subroutine read_topo_field(filename, hgt, we_dim, sn_dim, istatus)
  58
  59    use netcdf
  60
  61    implicit none
  62
  63    ! Arguments
  64    character (len=*), intent(in) :: filename
  65    real, pointer, dimension(:,:) :: hgt
  66    integer, intent(out)          :: we_dim, sn_dim
  67    integer, intent(out)          :: istatus
  68
  69    ! Local variables
  70    integer :: ncid, topo_varid, we_dimid, sn_dimid
  71    character (len=NF90_MAX_NAME) :: toponame, we_name, sn_name
  72
  73
  74    istatus = nf90_open(trim(filename), 0, ncid)
  75    if (istatus /= NF90_NOERR) then
  76       write(6,*) ' '
  77       write(6,*) 'Error: Could not open file '//trim(filename)
  78       write(6,*) ' '
  79       istatus = 1
  80       return
  81    end if
  82
  83    sn_name = 'south_north'
  84    istatus = nf90_inq_dimid(ncid, sn_name, sn_dimid)
  85    if (istatus /= NF90_NOERR) then
  86       write(6,*) ' '
  87       write(6,*) 'Error: Could not get ID of dimension south_north'
  88       write(6,*) ' '
  89       istatus = 1
  90       return
  91    end if
  92
  93    istatus = nf90_inquire_dimension(ncid, sn_dimid, sn_name, sn_dim)
  94    if (istatus /= NF90_NOERR) then
  95       write(6,*) ' '
  96       write(6,*) 'Error: Could not get south_north dimension'
  97       write(6,*) ' '
  98       istatus = 1
  99       return
 100    end if
 101
 102    we_name = 'west_east'
 103    istatus = nf90_inq_dimid(ncid, we_name, we_dimid)
 104    if (istatus /= NF90_NOERR) then
 105       write(6,*) ' '
 106       write(6,*) 'Error: Could not get ID of dimension west_east'
 107       write(6,*) ' '
 108       istatus = 1
 109       return
 110    end if
 111
 112    istatus = nf90_inquire_dimension(ncid, we_dimid, we_name, we_dim)
 113    if (istatus /= NF90_NOERR) then
 114       write(6,*) ' '
 115       write(6,*) 'Error: Could not get west_east dimension'
 116       write(6,*) ' '
 117       istatus = 1
 118       return
 119    end if
 120
 121    toponame = 'HGT_M'
 122    istatus = nf90_inq_varid(ncid, toponame, topo_varid)
 123    if (istatus /= NF90_NOERR) then
 124       write(6,*) ' '
 125       write(6,*) 'Error: Could not get ID of variable HGT_M'
 126       write(6,*) ' '
 127       istatus = 1
 128       return
 129    end if
 130
 131    allocate(hgt(we_dim,sn_dim))
 132
 133    istatus = nf90_get_var(ncid, topo_varid, hgt)
 134    if (istatus /= NF90_NOERR) then
 135       write(6,*) ' '
 136       write(6,*) 'Error: Could not read HGT_M field'
 137       write(6,*) ' '
 138       istatus = 1
 139       deallocate(hgt)
 140       return
 141    end if
 142
 143    istatus = nf90_close(ncid)
 144
 145    istatus = 0
 146
 147 end subroutine read_topo_field
 148
 149
 150 subroutine write_elev_angles(field, we_dim, sn_dim, kdim, istatus)
 151
 152    use write_met_module
 153
 154    implicit none
 155
 156    ! Arguments
 157    real, dimension(we_dim, sn_dim, kdim), intent(in) :: field
 158    integer, intent(in)                               :: we_dim, sn_dim, kdim
 159    integer, intent(out)                              :: istatus
 160
 161    ! Local variables
 162    integer :: i
 163    type (met_data) :: met_angle
 164
 165    call write_met_init('ELEVANGLES', .true., '0000-00-00_00:00:00', istatus)
 166
 167    if (istatus /= 0) then
 168       write(6,*) ' '
 169       write(6,*) 'Error opening output file ELEVANGLES'
 170       write(6,*) ' '
 171       istatus = 1
 172       return
 173    end if
 174
 175    met_angle%version = 5
 176    met_angle%iproj = PROJ_MERC
 177    met_angle%field = 'TOPO_ELEV'
 178    met_angle%units = 'degrees'
 179    met_angle%desc  = 'Topography elevation'
 180 !   met_angle%iproj = ...
 181 !   met_angle%truelat1 = ...
 182 !   met_angle%blah = ...
 183 !   met_angle%blah2 = ...
 184
 185    allocate(met_angle%slab(we_dim, sn_dim))
 186
 187    do i=1,kdim
 188
 189       met_angle%xlvl = real(i)
 190       met_angle%slab(:,:) = field(:,:,i)
 191       call write_next_met_field(met_angle, istatus)
 192
 193       if (istatus /= 0) then
 194          write(6,*) 'Error writing data to output file ELEVANGLES'
 195          istatus = 1
 196          deallocate(met_angle%slab)
 197          return
 198       end if
 199
 200    end do
 201
 202    call write_met_close()
 203
 204    istatus = 0
 205
 206 end subroutine write_elev_angles