-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathdirection_grid_mod.f90
105 lines (71 loc) · 2.52 KB
/
direction_grid_mod.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
module direction_grid_mod
implicit none
integer, parameter :: dir_grid_unit=70
double precision, dimension(:,:), allocatable :: direction_grid
character (len=256), parameter :: dir_grid_file = "direction_grid.bin"
logical, dimension(:,:), allocatable :: direction_mask
contains
! initialize memory
subroutine read_direction_grid_init(number_x_grid,number_y_grid)
implicit none
integer, intent(in) :: number_x_grid,number_y_grid
integer :: istat
allocate (direction_mask(number_x_grid,number_y_grid), direction_grid(number_x_grid,number_y_grid), stat=istat)
if (istat /= 0) THEN
WRITE(6,*) "allocation error in read_direction_grid_init"
stop
ENDIF
end subroutine read_direction_grid_init
! clear memory
subroutine read_direction_grid_clear()
implicit none
integer :: istat
deallocate( direction_mask, direction_grid, stat=istat)
if (istat /= 0) THEN
WRITE(6,*) "deallocation error in read_direction_grid_clear", istat
stop
ENDIF
end subroutine read_direction_grid_clear
! write out file
subroutine write_direction_grid(number_x_grid,number_y_grid)
implicit none
integer, intent(in) :: number_x_grid,number_y_grid
integer :: istat, rec_num, x_counter, y_counter
open(unit=dir_grid_unit, file=dir_grid_file, access="direct", form="unformatted", recl=8, status="replace", iostat=istat)
if (istat /= 0) THEN
WRITE(6,*) "error opening direction grid file for writing: ", trim(adjustl(dir_grid_file)), istat
stop
ENDIF
rec_num = 0
do x_counter = 1, number_x_grid
do y_counter = 1, number_y_grid
rec_num = rec_num + 1
if( .not. direction_mask(x_counter, y_counter)) THEN
write(dir_grid_unit, rec=rec_num) 0.
else
write(dir_grid_unit,rec=rec_num) direction_grid(x_counter,y_counter)
endif
end do
end do
close(unit=dir_grid_unit)
end subroutine write_direction_grid
! read in file
subroutine read_direction_grid(number_x_grid,number_y_grid)
implicit none
integer, intent(in) :: number_x_grid,number_y_grid
integer :: istat, rec_num, x_counter, y_counter
open(unit=dir_grid_unit, file=dir_grid_file, access="direct", form="unformatted", recl=8, status="old", iostat=istat)
if (istat /= 0) THEN
WRITE(6,*) "error opening direction grid file for reading: ", trim(adjustl(dir_grid_file)), istat
stop
ENDIF
rec_num = 0
do x_counter = 1, number_x_grid
do y_counter = 1, number_y_grid
rec_num = rec_num + 1
read(dir_grid_unit,rec=rec_num) direction_grid(x_counter,y_counter)
end do
end do
close(unit=dir_grid_unit)
end subroutine read_direction_grid
end module direction_grid_mod