读取全国160个站资料用于GrADS绘站点图（源程序）

! ######################################################################################
! Description: This programe is designed for reading the temperature &
!    rainfall data from 160 national meteorological stations. The output
!    is in GrADS station format including the CTL file for STNMAP.
!
! Version: 1.0
! Author: cwystar
! Date: 2005-06-16 For Ver 1.0
!
! History: 2005-06-16 For Ver 1.0
! #######################################################################################

program main
implicit none

! #######################################################################################
! MODIFIED PARAMETER --------------------------------------------------------------------
integer, parameter :: dataFirstYer = 1951, dataLastYer = 2001 ! 数据起止年份

character ( len = * ), parameter :: dataInfo = "D:/data/160station/全国160个站资料说明.txt" ! 输入的台站信息文件
character ( len = * ), parameter :: dataPath = "D:/data/160station/temperature" ! 输入的数据存放路径，末尾带/
character ( len = 1 ), parameter :: dataType = "T"  ! R -> precipation T -> temperature ! 温度或降水资料

character ( len = * ), parameter :: outputPath = "E:/20050615/readStaData"  ! 输出文件路径，以/结尾
character ( len = * ), parameter :: outputFileName = "t160"  ! 不含扩展名的输出文件名，用于CTL、MAP文件等

! END OF MODIFIED PARAMETER -------------------------------------------------------------
! #########################################################################

integer :: i, yer, mon
integer :: stdTemp, nlev, nflag
real :: tim
real, dimension( 160 ) :: lat, lon
real, dimension( dataFirstYer : dataLastYer, 12, 160 ) :: val
character( len = 8 ), dimension( 160 ) :: stdid
character( len = 100 ) :: temp
character( len = 4 ), dimension( 12 ) :: fileName = (/ "1601", "1602", "1603", "1604", "1605", &
 & "1606", "1607", "1608", "1609", "1610", "1611", "1612" /)

! 读入台站信息 -------------------
open ( 1, file = dataInfo, status = "old" )
do i = 1, 21
 read ( 1, "(a)" ) temp
end do

do i = 22, 181
 read( 1, "(a)" ) temp
 if ( i <= 9 ) then
  temp( 1:12 ) = ""
 else if ( i >= 10 .and. i <= 99 ) then
  temp( 1:13 ) = ""
 else if ( i >= 100 ) then
  temp( 1:14 ) = ""
 end if
 
 read( temp, * ) stdTemp, lat( i - 21 ), lon( i - 21 )
 write( stdid( i - 21 ), "(i8)" ) stdTemp
 stdid( i - 21 ) = trim( adjustl( stdid( i - 21 ) ) )
end do
close( 1 )

! 读入各台站所有时次的数据 -------------------
do mon = 1, 12
 open ( 2, file = dataPath // dataType // fileName( mon ) // ".txt", status = "old" )
 read( 2, * ) ( ( val( yer, mon, i ), i = 1, 160 ), yer = dataFirstYer, dataLastYer )
 close( 2 )
end do

! 输出台站数据 -------------------
tim = 0.
nflag = 1
open ( 3, file = outputFileName // ".dat", status = "unknown", form = "unformatted" )
do yer = dataFirstYer, dataLastYer
 do mon = 1, 12
  nlev = 1
  do i = 1, 160
   write( 3 ) stdid( i ), lat( i ), lon( i ), tim, nlev, nflag
   write( 3 ) val( yer, mon, i )
   if ( i == 160 ) then
    nlev = 0
    write( 3 ) stdid( i ), lat( i ), lon( i ), tim, nlev, nflag
   end if
  end do

 end do
end do
close( 3 )

! 输出相应的CTL文件 -------------------
open ( 4, file = outputPath // outputFileName // ".ctl", status = "unknown" )
write( 4, "(a)" ) "dset " // outputPath // outputFileName // ".dat"
write( 4, "(a)" ) "dtype station"
write( 4, "(a)" ) "stnmap " // outputPath // outputFileName // ".map"
write( 4, "(a)" ) "options sequential"
write( 4, "(a)" ) "undef -999.0"

if ( dataType == "R" ) then
 write( 4, "(a)" ) "title 160 station rainfall"
else if ( dataType == "T" ) then
 write( 4, "(a)" ) "title 160 station temperature"
end if

write( temp, "(i)" ) ( dataLastYer - dataFirstYer + 1 ) * 12
write( 4, "(a)" ) "tdef " // trim( adjustl( temp ) ) // " linear jan1951 1mo"
write( 4, "(a)" ) "vars 1"

if ( dataType == "R" ) then
 write( 4, "(a)" ) "r 0 99 station rainfall"
else if ( dataType == "T" ) then
 write( 4, "(a)" ) "t 0 99 station temperature"
end if

write( 4, "(a)" ) "endvars"
close( 4 )

write( *,* ) "--- The data and relatived ctl file have been outputed successfully. ---"
write( *,* ) "--- Please run STNMAP next. ---"
write( *,* )

end program main