저번 포스팅과 이어집니다
저번 포스팅에서 포트란에 바람 방향을 입력시켜주기 위한 각도 변환 과정을 알아봤습니다
오늘은 직접 코드로 계산을 해보겠습니다
준비물)
오늘 코딩은 미리 준비물이 필요합니다
1시간 간격으로 24시간 동안 관측된 하루 동안의 바람을 가지고 계산을 해볼 것이기 때문에
시간대 별 풍속과 풍향 자료가 필요합니다
0:00 1.1 290
1:00 0.5 20
2:00 0.5 200
3:00 2.9 50
4:00 3.3 50
5:00 1.6 340
6:00 1.5 320
7:00 2.0 360
8:00 2.1 50
9:00 2.3 70
10:00 2.1 90
11:00 1.6 70
12:00 2.3 90
13:00 2.3 90
14:00 1.6 90
15:00 1.7 90
16:00 1.7 50
17:00 0.5 110
18:00 0.3 0
19:00 0.2 0
20:00 0.9 160
21:00 1.8 140
22:00 0.7 200
23:00 1.1 180저희가 사용할 데이터 세트입니다
순서대로 관측 시간, 풍속, 풍향을 나타냅니다
관측시간은 시작시간을 기준으로 나타내기 때문에 00:00 은 00:00부터 00:59까지를 말합니다
이렇게 데이터 세트를 만들면 모두 끝났습니다
코드)
다음은 코드 입니다
REAL, DIMENSION(24) :: ws, wd
REAL, DIMENSION(24) :: u, v, cpwd
REAL, parameter :: pi = 4. * atan(1.)
REAL :: ap1, ap2, ap3, ap4
CHARACTER(2), parameter :: ug=char(int(Z'C2'))//char(int(Z'B0'))
open(1, file = 'time.dat')
open(2, file = 'time2.dat')
write(2, 60) 'u component', 'v compnent', 'compose component'
60 format(a, 3x, a, 3x, a)
do i = 1, 24
read(1, 10) ws(i), wd(i)
10 format(6x, f3.1, 1x, f3.0)
u(i) = ws(i) * cos((270. - wd(i)) * (pi/180))
v(i) = ws(i) * sin((270. - wd(i)) * (pi/180))
cpwd(i) = atan2(v(i), u(i))
ap1 = ap1 + u(i)
ap2 = ap2 + v(i)
ap3 = ap3 + cpwd(i)
write(2, 20) u(i), v(i), cpwd(i)
20 format(f10.4, 2x, f10.4, 5x, f10.4)
end do
rp1 = ap1/24.
rp2 = ap2/24.
rp3 = (270. - ((atan2(rp2, rp1)) * (180/pi))) - 360.
write(2, 30) "compose u direction componet :", rp1
30 format(a,f10.4, 1x, 'm/s')
write(2, 40) "compose v direction component :", rp2
40 format(a, f10.4, 1x, 'm/s')
write(2, 50) "compose wind speed :", sqrt(rp1**2 + rp2**2)
50 format(a,f10.4, 1x, 'm/s')
write(2, 70) 'compose wind direction', rp3, ug
70 format(a, f8.2, 1x, a)
end이번 코드는 길기 때문에 설명을 잘 따라오셔야 합니다
세부 분석
1) 변수 선언과 DIMENSION 선언
REAL, DIMENSION(24) :: ws, wd
REAL, DIMENSION(24) :: u, v, cpwd
REAL, parameter :: pi = 4. * atan(1.)
REAL :: ap1, ap2, ap3, ap4
CHARACTER(2), parameter :: ug=char(int(Z'C2'))//char(int(Z'B0'))첫 번째 부분입니다
저희가 작성한 데이터를 읽어와야 하는데 저런 자료들을 차원, 즉 DIMENSION을 가지고 있습니다
1,1 크기가 아니면 포트란은 거의 대부분 DIMENSION을 가진다고 생각하시면 됩니다
24시간 동안이니까 크기는 24로 잡아준 모습입니다
만약 2일동안으로 만들어서 하겠다 하시면 DIMENSION(48) 로 하시면 됩니다
다른 변수들은 새로 변수 만들고 pi를 선언하는 부분이라 큰 차이는 없습니다
마지막 줄은 각도 표시를 위해 만든 변수입니다
저 부분은 굳이 안하셔도 되지만 혹시 각도 표시가 필요하신 분은 그대로 적어서 쓰시면 됩니다
2) 파일 읽어오기
open(1, file = 'time.dat')
open(2, file = 'time2.dat')저희가 만든 파일을 읽어와야겠죠??
open 메소드를 통해서 1번에 데이터를 읽어오고 2번에 결과를 작성해주겠습니다
3) 라인 이름 만들어주기
write(2, 60) 'u component', 'v compnent', 'compose component'
60 format(a, 3x, a, 3x, a)이 부분은 제가 결과가 출력될 때 이 줄이 어떤 값이다라는 걸 알기 쉽게 하려고 추가한 부분입니다
필요 없다고 생각하시는 분들은 넘어가셔도 좋습니다
4) 반복문으로 파일들 읽어오기
do i = 1, 24
read(1, 10) ws(i), wd(i)
10 format(6x, f3.1, 1x, f3.0)
u(i) = ws(i) * cos((270. - wd(i)) * (pi/180))
v(i) = ws(i) * sin((270. - wd(i)) * (pi/180))
cpwd(i) = atan2(v(i), u(i))다음은 반복문을 사용해서 데이터를 읽어오는 부분입니다
i라는 변수에 1부터 24까지 먼저 할당을 시켜줍니다
read 메소드를 통해서 1번에서 ws(i)와 wd(i)를 읽어오겠다 라는 뜻입니다
format문에서도 중요한게 올라가서 데이터 세트를 보시면 저희가 필요한 부분은 풍속과 풍향만이지 시간은 필요가 없습니다
그래서 읽어오는 형식을 지정할 때 6x로 6칸을 띄고 읽어와야 풍속과 풍향만 읽어올 수 있습니다
그 다음은 1부터 24까지의 대입하면서 계산하는 과정입니다
0시 자료를 예로 들어서 설명해보겠습니다
i = 1이 0시이기 때문에 0:00 1.1 290 여기서 1.1이 ws(1), 290이 wd(1)이 됩니다
이 값들이 24까지 순서대로 들어가면서 u(i), v(i), cpwd(i)를 계산합니다
여기서 중요한 부분은 cpwd를 계산할 때 꼭 그냥 atan가 아니라 atan2를 사용해야합니다
atan과 atan2의 차이점은 atan2는 값의 범위를 고려해서 계산을 합니다
혹시 0이나 음수가 나올 때는 전부 다르게 계산을 하는 것이죠
바람 방향을 벡터로 나타내 때 음수가 나오기 때문에 atan2를 사용합니다
혹시 atan2를 더 알아보고 싶으신 분들은 위에 링크로 가시면 됩니다
그리고 반드시 atan2를 사용할 때는 (x, y)가 아니라 (y, x)로 작성을 해야합니다
저희가 계산할 공식이 atan(v/u)이기 때문에 atan2(v, u)로 작성을 해주면 됩니다
4) 출력할 Format 설정해주기
ap1 = ap1 + u(i)
ap2 = ap2 + v(i)
ap3 = ap3 + cpwd(i)
write(2, 20) u(i), v(i), cpwd(i)
20 format(f10.4, 2x, f10.4, 5x, f10.4)
end do
rp1 = ap1/24.
rp2 = ap2/24.
rp3 = (270. - ((atan2(rp2, rp1)) * (180/pi))) - 360.
write(2, 30) "compose u direction componet :", rp1
30 format(a,f10.4, 1x, 'm/s')
write(2, 40) "compose v direction component :", rp2
40 format(a, f10.4, 1x, 'm/s')
write(2, 50) "compose wind speed :", sqrt(rp1**2 + rp2**2)
50 format(a,f10.4, 1x, 'm/s')
write(2, 70) 'compose wind direction', rp3, ug
70 format(a, f8.2, 1x, a)
end복잡해 보이지만 정말 별거 없는 부분입니다
사진 설명을 입력하세요.
이 공식을 사용하기 위해 값을 계산하는 것 뿐입니다
ap1, ap2, ap3라는 변수에 u(i), v(i), cpwd(i)를 하나씩 더해서 누적합을 구합니다
나중에 평균을 구할 때 필요한 부분입니다
그리고 이제 write문을 통해 출력을 해줍니다
저는 위에 라인 이름들을 표시해주려고 했기 때문에 띄어쓰기를 넣어서 맞춰주려고 했습니다
마지막에 출력된 데이터를 보여드리면서 다시 한 번 말씀드리겠습니다
rp1, rp2에 u성분과 v성분을 24로 나눠서 평균을 구해줬습니다
rp3은 합성 풍향을 구하는 부분인데 저 부분은 쉽게 되지 않아서 이렇게 저렇게 하다가 만든 코드라
추후 수정해서 다시 올리겠습니다
저렇게 해도 알맞은 값이 나오니 참고만 하시면 될 것 같습니다
다른 부분은 이제 전부 format문이라 읽어보시면 될 것 같구여
sqrt(rp1**2 + rp2**2) 이 부분만 설명드리면 될 것 같습니다
sqrt는 square root의 약자로 말 그대로 루트입니다
합성 풍속 공식이 루트로 계산이 되게 때문에 사용했습니다
ug는 아까 맨 위에서 선언했던 파라미터 중에 각도 표시입니다
결과물)
제 코드로 출력한 결과물입니다
u component v compnent compose component
1.0337 -0.3762 -0.3491
-0.1710 -0.4698 -1.9199
0.1710 0.4698 1.2217
-2.2215 -1.8641 -2.4435
-2.5279 -2.1212 -2.4435
0.5472 -1.5035 -1.2217
0.9642 -1.1491 -0.8727
-0.0000 -2.0000 -1.5708
-1.6087 -1.3499 -2.4435
-2.1613 -0.7866 -2.7925
-2.1000 -0.0000 -3.1416
-1.5035 -0.5472 -2.7925
-2.3000 -0.0000 -3.1416
-2.3000 -0.0000 -3.1416
-1.6000 -0.0000 -3.1416
-1.7000 -0.0000 -3.1416
-1.3023 -1.0927 -2.4435
-0.4698 0.1710 2.7925
0.0000 -0.3000 -1.5708
0.0000 -0.2000 -1.5708
-0.3078 0.8457 1.9199
-1.1570 1.3789 2.2689
0.2394 0.6578 1.2217
-0.0000 1.1000 1.5708
compose u direction componet : -0.8531 m/s
compose v direction component : -0.3807 m/s
compose wind speed : 0.9342 m/s
compose wind direction 65.95 °
저는 약간 가시성이 좋은 출력물을 선호하는 타입이라 이렇게 만들었지만
혼자서 해보시는 분들은 이렇게까지 안하셔도 좋습니다
계산 결과를 만들어냈다는 것부터 잘하신 겁니다
익숙해지시면 저처럼 보기 좋게 만들어보시는 것도 좋습니다
오늘도 긴 글 읽어주셔서 감사합니다!
'Fortran' 카테고리의 다른 글
| Fortran으로 netCDF 파일 읽고 평균 구하기(gfortran 컴파일 에러) (0) | 2023.08.24 |
|---|---|
| 작업7 : Fortran으로 파장 별 방출에너지 구하기 (0) | 2022.10.26 |
| 작업5 : Fortran으로 합성 풍향 풍속 구하기1 (0) | 2022.10.26 |
| 작업4 : Fortran으로 시간방정식 계산하기 (0) | 2022.10.26 |
| 작업3: Fortran으로 태양의 적위 계산하기 (0) | 2022.10.26 |