-
QQ截图20170915175358.png (39.82 KB, 下载次数: 506)
Program weather
Implicit None
Integer l, i, j, m, daysinyear, k, mm, e, z
Integer, Allocatable :: site1(:), site2(:), year(:), month(:), day(:)
Real, Allocatable :: lat(:), lon(:), h(:), tmax(:), tmin(:), lonx(:), laty(:)
Integer :: ave, ier
Real :: h2
Character (45) :: filename = 'SURF_CLI_CHN_MUL_DAY-TEM-12001-123456.txt'
Character (8) :: outfile = '12345678'
Character (Len=20) :: line
!=====================================
Open (5, File='D:\气象数据New\经纬度高程\经纬度高程.txt')
e = 2421
Allocate (site2(e), lonx(e), laty(e))
Do z = 1, e
Read (5, *) site2(z), lonx(z), laty(z), h2
write(*,*) site2(z), lonx(z), laty(z), h2
End Do
Deallocate (site2, lonx, laty)
Close (5)
Do l = 2010, 2015
Do j = 1, 12
Write (filename(32:37), '(i6)') l*100 + j
Open (Unit=11, File='D:\气象数据\气温\'//filename, Action='READ')
m = 0
Do
Read (11, *, Iostat=ier) line
If (ier==-1) Exit ! end-of-file
m = m + 1
End Do
Allocate (site1(m), lat(m), lon(m), h(m), year(m), month(m), day(m), tmax(m), tmin(m))
Rewind (11)
mm = 0
Do i = 1, m
mm = mm + 1
Read (11, *) site1(mm), lat(mm), lon(mm), h(mm), year(mm), month(mm), day(mm), ave, tmax(mm), tmin(mm)
lat(mm) = lat(mm)/100
lon(mm) = lon(mm)/100
tmax(mm) = tmax(mm)*0.1
tmin(mm) = tmin(mm)*0.1
If (h(mm)>=100000) Then
h(mm) = (h(mm)-100000)
End If
h(mm) = h(mm)*0.1
If (tmax(mm)==3270.0 .Or. tmax(mm)==3276.6 .Or. tmin(mm)==3270.0 .Or. tmin(mm)==3276.6) mm = mm - 1
End Do
Close (Unit=11)
Do k = 1, daysinyear(l, j)
Write (outfile(1:8), '(i8)') l*10000 + j*100 + k
Open (Unit=15, File='D:\气象数据\tmax\'//outfile//'.txt')
Open (Unit=16, File='D:\气象数据\tmin\'//outfile//'.txt')
Do i = 1, mm
If (l==year(i) .And. j==month(i) .And. k==day(i)) Then
Do While (site1(i)==site2(z))
Write (15, '(i5,2x,f6.2,2x,f5.2,2x,f12.4,2x,f12.4,2x,f6.2)') &
&site1(i), lon(i), lat(i), lonx(z), laty(z), tmax(i) + (h(i)/100)*0.6
! "(i5,2x,f5.2,2x,f6.2,2x,f5.0,2x,i5,2x,i2,2x,i2,2x,f6.2)"
Write (16, '(i5,2x,f6.2,2x,f5.2,2x,f12.4,2x,f12.4,2x,f6.2)') &
&site1(i), lon(i), lat(i), lonx(z), laty(z), tmin(i) + (h(i)/100)*0.6
Print *, site1(i), year(i), month(i), day(i)
End Do
End If
End Do
Close (Unit=15)
End Do
Deallocate (site1, lat, lon, h, year, month, day, tmax, tmin)
End Do
End Do
End Program weather
!=======================================
Integer Function daysinyear(year, month)
Integer, Intent (In) :: year, month
Integer :: daysinmonth(12) = [ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 ]
If (((mod(year,4)==0) .And. (mod(year,100)/=0)) .Or. (mod(year,400)==0)) Then
daysinmonth(2) = 29
Else
daysinmonth(2) = 28
End If
daysinyear = daysinmonth(month)
End Function daysinyear
QQ截图20170915175358.png (39.82 KB, 下载次数: 506)
chiangtp 发表于 2017-9-15 21:31
line24: CLOSE(UNIT=5)?
不要用 UNIT= 0, 5, 6 對應到 diskfile (尤其是又有 READ/WRITE(*,...)時, compile ...
wxy 发表于 2017-9-16 00:15
额 我换成其他10,13等还是同样错误
wxy 发表于 2017-9-16 00:15
额 我换成其他10,13等还是同样错误，我在运行时，到第24行运行都没问题，是16-24行这部分循环运行结束时报 ...
chiangtp 发表于 2017-9-16 00:44
抱歉,
delete line 23: Deallocate (site2, lonx, laty)
END PROGRAM之前, 方可Deallocate (site2, lonx, ...
REWIND(UNIT=11)
mm = 0
DO i = 1, m
mm = mm + 1
READ(11,*) site1(mm), lat(mm), lon(mm), h(mm), year(mm), month(mm), day(mm), ave, tmax(mm), tmin(mm)
IF( tmax(mm)==3270.0 .OR. tmax(mm)==3276.6 .OR. tmin(mm)==3270.0 .OR. tmin(mm)==3276.6 ) THEN
mm = mm - 1
CYCLE
END IF
IF( h(mm) >= 100000.0 ) h(mm) = h(mm) - 100000.0
h (mm) = h (mm)*0.1
tmax(mm) = tmax(mm)*0.1
tmin(mm) = tmin(mm)*0.1
lon (mm) = lon (mm)/100.0
lat (mm) = lat (mm)/100.0
END DO
CLOSE(UNIT=11)
DO k = 1, daysinyear(l,j)
WRITE(outfile(1:8),'(i8)') l*10000 + j*100 + k
OPEN( UNIT=15, FILE='...' )
OPEN( UNIT=16, FILE='...' )
100 FORMAT(i5,2x,f6.2,2x,f5.2,2x,f12.4,2x,f12.4,2x,f6.2)
DO i = 1, mm
IF( l/=year(i) .OR. j/=month(i) .OR. k==day(i) ) CYCLE
DO z = 1, e
IF( site1(i) == site2(z) ) THEN
WRITE(15,100) site1(i), lon(i), lat(i), lonx(z), laty(z), tmax(i)+(h(i)/100)*0.6
WRITE(16,100) site1(i), lon(i), lat(i), lonx(z), laty(z), tmin(i)+(h(i)/100)*0.6
WRITE(*,*) site1(i), year(i), month(i), day(i)
END IF
END DO
END DO
CLOSE(UNIT=15)
END DO
chiangtp 发表于 2017-9-16 01:05
[Fortran] 纯文本查看 复制代码REWIND(UNIT=11)
mm = 0
DO i = 1, m
[mw_shl_code=fortran,true]
CLOSE(UNIT=15)
CLOSE(UNIT=16)
END DO
e = 2421
ALLOCATE( site2(e), lonx(e), laty(e) )
OPEN( UNIT=15, FILE='...' )
DO z = 1, e
READ(15,*) site2(z), lonx(z), laty(z), h2
WRITE(*,*) site2(z), lonx(z), laty(z), h2
END DO
CLOSE( UNIT=15 )
!-----------------
Do l = 2010, 2015
Do j = 1, 12
WRITE(filename(32:37), '(i6)') l*100 + j
OPEN( UNIT=11, FILE='...', ACTION='READ' )
m = 0
DO
READ(11,*,IOSTAT=ier) line
IF( ier == -1 ) EXIT ! end-of-file
m = m + 1
END DO
ALLOCATE( site1(m), lat(m), lon(m), h(m), year(m), month(m), day(m), tmax(m), tmin(m) )
REWIND( UNIT=11 )
mm = 0
DO i = 1, m
mm = mm + 1
READ(11,*) site1(mm), lat(mm), lon(mm), h(mm), year(mm), month(mm), day(mm), ave, tmax(mm), tmin(mm)
IF( tmax(mm)==3270.0 .OR. tmax(mm)==3276.6 .OR. tmin(mm)==3270.0 .OR. tmin(mm)==3276.6 ) THEN
mm = mm - 1
CYCLE
END IF
IF( h(mm) >= 100000.0 ) h(mm) = h(mm) - 100000.0
h (mm) = h (mm)*0.1
tmax(mm) = tmax(mm)*0.1
tmin(mm) = tmin(mm)*0.1
lon (mm) = lon (mm)/100.0
lat (mm) = lat (mm)/100.0
END DO
CLOSE( UNIT=11 )
!-----
DO k = 1, daysinyear(l,j)
WRITE(outfile(1:8),'(i8)') l*10000 + j*100 + k
OPEN( UNIT=15, FILE='...' )
OPEN( UNIT=16, FILE='...' )
100 FORMAT(i5,2x,f6.2,2x,f5.2,2x,f12.4,2x,f12.4,2x,f6.2)
DO i = 1, mm
IF( l/=year(i) .OR. j/=month(i) .OR. k==day(i) ) CYCLE
DO z = 1, e
IF( site1(i) == site2(z) ) THEN
WRITE(15,100) site1(i), lon(i), lat(i), lonx(z), laty(z), tmax(i)+(h(i)/100)*0.6
WRITE(16,100) site1(i), lon(i), lat(i), lonx(z), laty(z), tmin(i)+(h(i)/100)*0.6
WRITE(*,*) site1(i), year(i), month(i), day(i)
END IF
END DO
END DO
CLOSE( UNIT=15 )
CLOSE( UNIT=16 )
END DO
DEALLOCATE( site1, lat, lon, h, year, month, day, tmax, tmin )
END DO
END DO
DEALLOCATE( site2, lonx, laty )
END PROGRAM weather
chiangtp 发表于 2017-9-16 01:18
[Fortran] 纯文本查看 复制代码e = 2421
ALLOCATE( site2(e), lonx(e), laty(e) )
1. 還是有點小錯: "k==day(i)" ---> "k/=day(i)"
2. 改善Loop的效率: "IF(l/=year(i) .OR. j/=month(i))"往前提 (not function of "iday" loop)
3. 重要的subscripts: "l, j, k" ---> "iyear, imonth, iday", 可讀性高
4. 數理/Fortran的文化: "z"不宜為subscript "iz"可也
5. for REAL: "100" ---> "100.0"
[mw_shl_code=fortran,true] INTEGER :: iyear, imonth, iday, nn
!-----------------
nn = 2421
ALLOCATE( site2(nn), lonx(nn), laty(nn) )
OPEN( UNIT=11, FILE='...' )
DO i = 1, nn
READ(11,*) site2(i), lonx(i), laty(i)
END DO
CLOSE( UNIT=11 )
!-----------------
DO iyear = 2010, 2015
DO imonth = 1, 12
WRITE(filename(32:37), '(i6)') iyear*100 + imonth
OPEN( UNIT=11, FILE='...', ACTION='READ' )
m = 0
DO
READ(11,*,IOSTAT=ier) line
IF( ier == -1 ) EXIT ! end-of-file
m = m + 1
END DO
ALLOCATE( site1(m), lat(m), lon(m), h(m), year(m), month(m), day(m), tmax(m), tmin(m) )
REWIND( UNIT=11 )
mm = 0
DO i = 1, m
mm = mm + 1
READ(11,*) site1(mm), lat(mm), lon(mm), h(mm), year(mm), month(mm), day(mm), ave, tmax(mm), tmin(mm)
IF( year(mm)/=iyear .OR. month(mm)/=imonth .OR. &
tmax(mm)==3270.0 .OR. tmax(mm)==3276.6 .OR. tmin(mm)==3270.0 .OR. tmin(mm)==3276.6 ) THEN
mm = mm - 1
CYCLE
END IF
IF( h(mm) >= 100000.0 ) h(mm) = h(mm) - 100000.0
h (mm) = h (mm)*0.1
tmax(mm) = tmax(mm)*0.1
tmin(mm) = tmin(mm)*0.1
lon (mm) = lon (mm)/100.0
lat (mm) = lat (mm)/100.0
END DO
CLOSE( UNIT=11 )
!-----
IF( mm == 0 ) CYCLE
!-----
DO iday = 1, daysinyear(iyear,imonth)
WRITE(outfile(1:8),'(i8)') iyear*10000 + imonth*100 + iday
OPEN( UNIT=15, FILE='...' )
OPEN( UNIT=16, FILE='...' )
100 FORMAT(i5,2x,f6.2,2x,f5.2,2x,f12.4,2x,f12.4,2x,f6.2)
DO i = 1, mm
IF( iday /= day(i) ) CYCLE
DO j = 1, nn
IF( site1(i) /= site2(j) ) CYCLE
WRITE(15,100) site1(i), lon(i), lat(i), lonx(j), laty(j), tmax(i)+(h(i)/100.0)*0.6
WRITE(16,100) site1(i), lon(i), lat(i), lonx(j), laty(j), tmin(i)+(h(i)/100.0)*0.6
WRITE(*,*) site1(i), year(i), month(i), day(i)
END DO
END DO
CLOSE( UNIT=15 )
CLOSE( UNIT=16 )
END DO
DEALLOCATE( site1, lat, lon, h, year, month, day, tmax, tmin )
END DO
END DO
DEALLOCATE( site2, lonx, laty )
END PROGRAM weather
chiangtp 发表于 2017-9-16 11:38
1. 還是有點小錯: "k==day(i)" ---> "k/=day(i)"
2. 改善Loop的效率: "IF(l/=year(i) .OR. j/=month(i) ...
wxy 发表于 2017-9-16 13:11
不好意思刚看到 非常感谢 ，我还是不太明白这个报错是什么意思，为甚末之前将deallocate放在下一个do循环 ...
DO i = 1, m
IF( a /= b ) CYCLE
... ! do something for (a==b)
END DO
or,
DO i = 1, m
IF( a == b ) THEN
... ! do something for (a==b)
END IF
END DO
REWIND( UNIT=11 )
mm = 0
DO i = 1, m
mm = mm + 1
READ(11,*) site1(mm), lat(mm), lon(mm), h(mm), year(mm), month(mm), day(mm), ave, tmax(mm), tmin(mm)
IF( year(mm)/=iyear .OR. month(mm)/=imonth .OR. &
tmax(mm)==3270.0 .OR. tmax(mm)==3276.6 .OR. tmin(mm)==3270.0 .OR. tmin(mm)==3276.6 ) mm = mm - 1
END DO
CLOSE( UNIT=11 )
WHERE( h(1:mm) >= 100000.0 ) h(1:mm) = h(1:mm) - 100000.0
h (1:mm) = h (1:mm)*0.1
tmax(1:mm) = tmax(1:mm)*0.1
tmin(1:mm) = tmin(1:mm)*0.1
lon (1:mm) = lon (1:mm)/100.0
lat (1:mm) = lat (1:mm)/100.0
chiangtp 发表于 2017-9-16 13:57
(1) 为甚末之前将deallocate放在下一个do循环之前就会报错，放在end program 之前就没事
DEALLOCATE(site ...
| 欢迎光临 Fortran Coder (http://bbs.fcode.cn/) | Powered by Discuz! X3.2 |
