!$$$ program documentation block
! . . . .
! program: diag2grad create GrADS files
! prgmmr: tahara org: np20 date: 2003-01-01
!
! abstract: This program create a GrADS station data file and its
! control file a GSI radiance diagnostic file.
!
! program history log:
! 2010-12-22 treadon - add this doc block
!
! contains
! adjust_lon - adjust longitude to degrees west
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
program diag2grads
! Include modules
use Utilities_Time
use read_diag
use grads_hdr
use write_station
implicit none
! Declare and set parameters
integer,parameter :: ft_namelist = 1 ! (in) namelist
integer,parameter :: max_n_chan = 900
! real(4),parameter :: Criterion_SolarZenAng = 96._8 ! 96 degrees is better
real(4),parameter :: Criterion_SolarZenAng = 90._8
! Declare namelist with defaults
logical :: retrieval = .false.
integer :: idsat = 14 ! HIRS14
integer :: n_chan = 19
integer :: ft_diag_top = 10 ! diagnostic files
integer :: ft_diag_end = 10
integer :: datetime(5) = (/ 2000, 1, 1, 0, 0 /)
integer :: tint = 6
real(4) :: North = 90. ! Statistics area
real(4) :: South = -90. ! Statistics area
real(4) :: West = 0. ! Statistics area
real(4) :: East = 360. ! Statistics area
integer :: qc = 0 ! 0:all 1:only qc pass -1:only rejected
integer :: solar = 0 ! 0:all 1:night 2:day
character(len=5) :: cqc(-1:1) = (/ 'ALL ', 'PASS ', 'RJCT ' /)
character(len=5) :: clandsea(0:2) = (/ 'ALL ', 'LAND ', 'SEA ' /)
character(len=5) :: csolar(0:2) = (/ 'ALL ', 'NIGHT', 'DAY ' /)
character(len=GRADS_MAXLEN_COMMENT) :: comment = ''
character(len=GRADS_MAXLEN_FILENAME) :: filename = ''
character(len=GRADS_MAXLEN_FILENAME+4) :: dsetname = ''
namelist / PARM / comment, ft_diag_top, ft_diag_end, filename, &
idsat, n_chan, datetime, tint, South, North, West, East, &
qc, solar, retrieval
! Declare variables for reading / processing satellite data
type(diag_header_fix_list ) :: header_fix
type(diag_header_chan_list),allocatable :: header_chan(:)
type(diag_data_name_list ) :: data_name
type(diag_data_fix_list ) :: data_fix
type(diag_data_chan_list ),allocatable :: data_chan(:)
type(diag_data_extra_list ),allocatable :: data_extra(:,:)
real(4),allocatable :: nuchan(:)
real(4),allocatable :: data_sfc(:)
real(4),allocatable :: data_lvl(:,:)
character(len=GRADS_LEN_ID) :: id
! Other variables
character(8):: stid
integer :: ft_diag
integer :: lsflag, iceflag, scan, mype
integer :: iflag, iexist, nelem_sfc,nelem_lvl
integer :: ttlmin_file, ttlmin_assign, ttlmin_assign0
integer :: time_work(5)
integer :: iobs_total
integer :: i, k, n
integer :: npred_radiag, n_bcor
integer :: lunstn,lunctl
integer :: nlev,nflag,itime_terminator
real(4):: rlat,rlon,rtim
npred_radiag = 5
lunstn = 21
lunctl = 22
! Read namelist
read(ft_namelist,PARM)
! Adjust longitude
call adjust_lon( East )
call adjust_lon( West )
if( East == West )then
West = 0.
East = 360.
endif
! Print namelist
write(6,*)' '
write(6,*)'*** NAMELIST ***************************************************'
write(6,*) 'Title : ', trim(comment)
write(6,*) 'FT # of diag files : ', ft_diag_top, ft_diag_end
write(6,*) 'GrADS data file name : ', trim(filename)
write(6,*) 'SAT ID : ', idsat
write(6,*) '# of channels : ', n_chan
write(6,*) 'Time of the 1st file : ', datetime
write(6,*) 'Interval of files(hr): ', tint
write(6,*) 'Area(S,N,W,E) : ', real( (/South, North, West, East/), 4 )
write(6,*) '0:All 1:used -1:rjct -2:All w/o any qc : ', qc
write(6,*) '0:All 1:night 2:day : ', solar
! Check number of channels
if (n_chan > max_n_chan) then
write(6,*)'### ERROR: n_chan=',n_chan,' > max_n_chan=',max_n_chan
stop 99
endif
! Open file to GrADS station data file
dsetname = trim(filename) // '.dat'
open(lunstn,file=dsetname,form='unformatted')
write(6,*)'Open unit lunstn=',lunstn,' to GrADS station data file ',trim(dsetname)
! Allocate channel list array
allocate( nuchan(n_chan) )
! Set GrADS data id with sat id
write(id,'(i8)') idsat
! Compute total minutes for the first diagnostic file
call Convert_DateTime_TotalMin( ttlmin_assign0, datetime )
! Big loop over diagnostic files
iobs_total = 0
do ft_diag = ft_diag_top, ft_diag_end
itime_terminator=0
! Read header
call read_radiag_header( ft_diag, npred_radiag, retrieval, header_fix, header_chan, data_name, iflag )
if (iflag/=0) then
write(6,*)'***ERROR*** problem reading header, iflag=',iflag
stop 99
endif
! Check for satellite and sensor
! if( header_fix%isat /= idsat )cycle
! if( header_fix%isat /= mod(idsat,100) )cycle !### TEMPORAL ###
! Check number of channels
if( header_fix%nchan /= n_chan )then
write(6,*) '### ERROR: UNEXPECTED NUMBER OF CHANNELS IS READ'
write(6,*) 'THE NUMBER IS', header_fix%nchan
stop 99
endif
! Check time
time_work(1) = header_fix%idate / 1000000
time_work(2) = mod( header_fix%idate, 1000000 ) / 10000
time_work(3) = mod( header_fix%idate, 10000 ) / 100
time_work(4) = mod( header_fix%idate, 100 )
time_work(5) = 0
call Convert_DateTime_TotalMin( ttlmin_file, time_work )
ttlmin_assign = ttlmin_assign0 + (ft_diag - ft_diag_top) * tint * 60
if( ttlmin_file /= ttlmin_assign )then
write(6,*) '### ERROR: TIMES ARE INCONSISTENT'
write(6,*) 'FILE :', header_fix%idate
write(6,*) 'ASSIGNED:', datetime, '+', (ft_diag - ft_diag_top) * tint, 'hours'
stop 99
endif
! Extract channel list
do n=1, n_chan
nuchan(n) = real( header_chan(n)%nuchan, 4 )
enddo
! Set number of level (channel) and surface (fix) data
nelem_lvl=header_fix%idiag
nelem_sfc=ireal_radiag-2
if (allocated(data_sfc)) deallocate(data_sfc)
allocate(data_sfc(nelem_sfc))
write(6,*)'Allocate rank-1 data_sfc as nelem_sfc=',nelem_sfc
if (allocated(data_lvl)) deallocate(data_lvl)
allocate(data_lvl(nelem_lvl,n_chan))
write(6,*)'Allocate rank-2 data_lvl as nelem_lvl=',nelem_lvl,' by n_chan=',n_chan
! Print header
write(6,*)' '
write(6,*) '*** HEADER OF DIAGNOSTIC FILE ***'
write(6,*) 'FIX PART:', header_fix
write(6,*) 'n_chan=',n_chan,' n_bcor=',n_bcor
do n=1, n_chan
write(6,*) 'CH', n, ':', header_chan(n)
enddo
! Loop over data in current diagnostic file
write(6,*)' '
write(6,*) '*** DATA OF DIAGNOSTIC FILE ***'
do
! Read data record
call read_radiag_data( ft_diag, header_fix, retrieval, data_fix, data_chan, data_extra, iflag )
if( iflag /= 0 )then
write(6,*)'***ERROR*** problem reading data iflag=',iflag,' iexist=',iexist
exit ! EOF or error
endif
! Check area
if( data_fix%lat > North .or. data_fix%lat < South )cycle
call adjust_lon( data_fix%lon )
if( West < East )then
if( data_fix%lon < West .or. East < data_fix%lon )cycle
else if( West > East )then
if( East < data_fix%lon .and. data_fix%lon < West )cycle
endif
! Check solar
if( ( solar == 1 .and. data_fix%solzen_ang < Criterion_SolarZenAng ) .or. &
( solar == 2 .and. data_fix%solzen_ang >= Criterion_SolarZenAng ) )then
cycle
endif
! Initialize surface data to missing
do i=1,nelem_sfc
data_sfc(i)=GRADS_MISSING
end do
! Load surface data
rlat = data_fix%lat
rlon = data_fix%lon
rtim = 0.
data_sfc(1) = data_fix%zsges
data_sfc(2) = data_fix%obstime
data_sfc(3) = data_fix%senscn_pos
data_sfc(4) = data_fix%satzen_ang
data_sfc(5) = data_fix%satazm_ang
data_sfc(6) = data_fix%solzen_ang
data_sfc(7) = data_fix%solazm_ang
data_sfc(8) = data_fix%sungln_ang
data_sfc(9) = data_fix%water_frac
data_sfc(10) = data_fix%land_frac
data_sfc(11) = data_fix%ice_frac
data_sfc(12) = data_fix%snow_frac
data_sfc(13) = data_fix%water_temp
data_sfc(14) = data_fix%land_temp
data_sfc(15) = data_fix%ice_temp
data_sfc(16) = data_fix%snow_temp
data_sfc(17) = data_fix%soil_temp
data_sfc(18) = data_fix%soil_mois
data_sfc(19) = data_fix%land_type
data_sfc(20) = data_fix%veg_frac
data_sfc(21) = data_fix%snow_depth
data_sfc(22) = data_fix%sfc_wndspd
data_sfc(23) = data_fix%qcdiag1
data_sfc(24) = data_fix%qcdiag2
! Initialize the channel (level) data
do n=1,n_chan
do i=1,nelem_lvl
data_lvl(i,n) = GRADS_MISSING
end do
end do
! Loop over channels
iexist = 0
do n = 1, n_chan
! Check Tb
if( qc /= -2 .and. ( abs(data_chan(n)%omgnbc) > 200. .or. &
data_chan(n)%tbobs < 50. .or. data_chan(n)%tbobs > 500. ) )cycle
! Check uss flag
if( ( qc == -1 .and. data_chan(n)%errinv >= 1.e-6 ) .or. &
( qc == 1 .and. data_chan(n)%errinv < 1.e-6 ) )cycle
! There is data to be written
iexist = 1
! Load arrays with level (channel) dependent data
data_lvl(1,n) = data_chan(n)%tbobs
data_lvl(2,n) = data_chan(n)%omgbc
data_lvl(3,n) = data_chan(n)%omgnbc
data_lvl(4,n) = data_chan(n)%errinv
data_lvl(5,n) = data_chan(n)%qcmark
data_lvl(6,n) = data_chan(n)%emiss
data_lvl(7,n) = data_chan(n)%tlap
if (header_fix%iversion < iversion_radiag) then
data_lvl( 8,n)=data_chan(n)%bifix(1)
data_lvl( 9,n)=data_chan(n)%bilap
data_lvl(10,n)=data_chan(n)%bilap2
data_lvl(11,n)=data_chan(n)%bicons
data_lvl(12,n)=data_chan(n)%biang
data_lvl(13,n)=data_chan(n)%biclw
if (retrieval) data_lvl(13,n)=data_chan(n)%bisst
else
data_lvl(8,n) = data_chan(n)%bicons
data_lvl(9,n) = data_chan(n)%biang
data_lvl(10,n) = data_chan(n)%biclw
data_lvl(11,n) = data_chan(n)%bilap2
data_lvl(12,n) = data_chan(n)%bilap
do i=1,header_fix%angord+1
data_lvl(12+i,n) = data_chan(n)%bifix(i)
end do
data_lvl(12+header_fix%angord+2,n)=data_chan(n)%bisst
endif
! End loop over channels
enddo
! Check if data should be written
if( iexist == 0 )cycle
iobs_total = iobs_total + 1
! Write station data
call write_station_data(lunstn,&
id, & ! id
rlat, & ! lat
rlon, & ! lon
rtim, & ! t
n_chan, & ! nlev(# of chan)
nelem_lvl, &
nelem_sfc, &
nuchan, & ! ch # assigned as hight
data_lvl, &
data_sfc )
! Periodically print data
if( mod(iobs_total,2000) == 0 )then
write(6,*) 'DATA #:', iobs_total
write(6,*) 'FIX PART:', data_fix
do n=1, n_chan, n_chan/2
write(6,*) 'CH SEQ #:', n, ':', data_lvl(:,n)
enddo
endif
! End loop over data records
enddo
! Write a record to indicate the end of time group
if (itime_terminator==0) then
write(6,*)'Write GrADS station data terminator for time ft_diag=',ft_diag
stid='time_end'
rlat=0.0
rlon=0.0
rtim=0.0
nlev=0
nflag=0
write(lunstn) stid,rlat,rlon,rtim,nlev,nflag
itime_terminator=1
endif
! End loop over times
enddo
! Close GrADS station file
if (itime_terminator==0) then
write(6,*)'Write GrADS station data terminator'
stid='end'
rlat=0.0
rlon=0.0
rtim=0.0
nlev=0
nflag=0
write(lunstn) stid,rlat,rlon,rtim,nlev,nflag
endif
write(6,*)'Close GrADS station data file attached to lunstn=',lunstn
close(lunstn)
! Create GrADS control file
dsetname = trim(filename) // '.ctl'
open(lunctl,file=dsetname,form='formatted')
write(6,*)'Open unit lunctl=',lunctl,' to GrADS control file ',dsetname
call write_station_ctl(lunctl,comment, &
filename, &
datetime, &
tint, &
ft_diag_end - ft_diag_top + 1, &
n_chan, &
nelem_lvl, &
nelem_sfc, &
data_name)
close(lunctl)
! Check total number of data
if( iobs_total <= 0 )then
write(6,*) '### NO DATA ###'
stop
endif
! End of programk !--- end of program
write(6,*)' '
write(6,*) '=== NORMAL END ==='
contains
!-------------------------------------------------------------
! Adjust longitude between [0:360)
!-------------------------------------------------------------
subroutine adjust_lon( lon )
real(4),intent(inout) :: lon
if( lon < 0. )then
lon = lon + 360.
else if( lon >= 360. )then
lon = lon - 360.
endif
end subroutine adjust_lon
end program diag2grads