! ------------------------------------------------------------------------
! Program that creates a 3D file with information about orography height.
! Useful to depict orography in vertical sections.
!
! File needed: gridded data with pressure surface information.
!
! If this information is not in a normal binary file (e.g., netCDF, GRIB or other
! kind of data file) you may use GrADS itself to convert it. You may do this with the
! fwrite option for set gxout command, like this:
!
!    set gxout fwrite
!    set fwrite output.bin
!    set lat ...     |- if necessary.
!    set lon ...     |
!    display pressure_surface
!    disable fwrite
!
! For this you must have a file with same levels and horizontal resolution of the
! data file that you want to obtain vertical profiles.
! Below is a fragment of GrADS script that plots the orography:
!
!    set lat -30
!    set lon 280 330
!    set lev 1000 300
!    ....
!    <plot your fields here>
!    ....
!    set clab off
!    set cint 5
!    set cthick 20
!    set cmax 0
!    set ccolor 1
!    d orog(t=1)
!
! As you probably open other files along with orography you have to put
! the appropriate number for variable orog. The t=1 prevents problems
! related with difference in times of the files.
!
! Mateus da Silva Teixeira.
! April / 2006
! ----------------------------------------------------------------------------
program orography
implicit none

! Required information for CTL file
integer, parameter :: xx=25,yy=13,zz=8 ! file dimensions
integer :: levs(8)=(/1000,925,850,700,600,500,400,300/) ! vertical levels
real :: rx=2.5,ry=2.5,ilat=-40.,ilon=270    ! horizontal resolution and lat,lon initial

! Variables
integer :: lev=1,i,j,k
real :: psfc(xx,yy),oro(xx,yy)

! Opening files
open(28,file='orog.ctl',status='unknown')
open(29,file='orog.bin',status='unknown',access='direct',form='unformatted',recl=xx*yy)
open(30,file='psfc.bin',status='old',access='direct',form='unformatted',recl=xx*yy)

read(30,rec=1) psfc ! Reading surface pressure data

do k=1,8
   do j=1,13
      do i=1,25
         oro(i,j)=-levs(k)+psfc(i,j)   ! Calculation of orography elevation
      end do
   end do
   write(29,rec=lev) oro
   lev=lev+1
end do

! Creating CTL file
write(28,'(A14)') 'DSET ^orog.bin'
write(28,'(A15)') 'TITLE Orography'
write(28,'(A15)') 'UNDEF -9.99e+33'
write(28,'(A4,1X,I2,1X,A6,1X,F5.1,1X,F3.1)') 'XDEF',xx,'LINEAR',ilon,rx
write(28,'(A4,1X,I2,1X,A6,1X,F5.1,1X,F3.1)') 'YDEF',yy,'LINEAR',ilat,ry
write(28,'(A4,1X,I2,1X,A6,1X,8I5)') 'ZDEF',zz,'LEVELS',(levs(i),i=1,zz)
write(28,'(A31)') 'TDEF 1  LINEAR 16Z06APR2006 1DY'
write(28,'(A6)') 'VARS 1'
write(28,'(A4,1X,I2,1X,A12)') 'orog',zz,'99 Orografia'
write(28,'(A7)') 'ENDVARS'

end program orography

