forked from metno/NORTRIP_multiroad
-
Notifications
You must be signed in to change notification settings - Fork 0
/
read_netcdf4_terrain_file.f90
83 lines (60 loc) · 3.28 KB
/
read_netcdf4_terrain_file.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
!read_netcdf_terrain_file.f90
subroutine read_netcdf4_terrain_file(filename_nc_sub,ncols_sub,nrows_sub,cellsize_sub)
use NORTRIP_multiroad_index_definitions
!Update to netcdf 4 and 64 bit in this version 2 of read_netcdf_terrain_file
use netcdf
implicit none
!include 'netcdf.inc'
integer status_nc !Error message
integer id_nc
integer dim_id_nc(2)
character(256) dimname_temp
integer var_id_nc
character(*) filename_nc_sub
character(256) temp_str
integer i,j,ii,jj
integer ncols_sub,nrows_sub
real cellsize_sub
!Open the netcdf file for reading
write(unit_logfile,'(2A)') ' Opening netcdf terrain file: ',trim(filename_nc_sub)
status_nc = NF90_OPEN (filename_nc_sub, NF90_NOWRITE, id_nc)
if (status_nc .NE. NF90_NOERR) write(unit_logfile,'(A,I)') 'ERROR opening netcdf file: ',status_nc
!Find out the x,y and time dimmensions of the file by looking at pressure variable
status_nc = NF90_INQ_DIMID (id_nc,dim_name_nc(x_index),dim_id_nc(x_index))
status_nc = NF90_INQUIRE_DIMENSION (id_nc,dim_id_nc(x_index),dimname_temp,dim_length_nc(x_index))
status_nc = NF90_INQ_DIMID (id_nc,dim_name_nc(y_index),dim_id_nc(y_index))
status_nc = NF90_INQUIRE_DIMENSION (id_nc,dim_id_nc(y_index),dimname_temp,dim_length_nc(y_index))
ncols_sub=dim_length_nc(x_index)
nrows_sub=dim_length_nc(y_index)
dim_start_nc=1
if (.not.allocated(array)) allocate(array(ncols_sub,nrows_sub))
if (.not.allocated(x_array)) allocate(x_array(ncols_sub))
if (.not.allocated(y_array)) allocate(y_array(nrows_sub))
allocate (var1d_nc(num_dims_nc,max(ncols_sub,nrows_sub))) !x and y and time maximum dimmensions
!Read the x, y and time values
do i=1,num_dims_terrain_nc
status_nc = NF90_INQ_VARID (id_nc, trim(dim_name_terrain_nc(i)), var_id_nc)
status_nc = NF90_GET_VAR (id_nc, var_id_nc, var1d_nc(i,:), start=(/dim_start_nc(i)/), count=(/dim_length_nc(i)/))
write(unit_logfile,'(3A,2f12.2)') ' ',trim(dim_name_terrain_nc(i)),' (min, max in km): ' &
,minval(var1d_nc(i,1:dim_length_nc(i))/1000.),maxval(var1d_nc(i,1:dim_length_nc(i))/1000.)
enddo
!x_array=var1d_nc(x_index,1:nrows_sub)
!y_array=var1d_nc(y_index,1:ncols_sub)
x_array=var1d_nc(x_index,1:ncols_sub)
y_array=var1d_nc(y_index,1:nrows_sub)
!Assumes x and y are the same
cellsize_sub=x_array(2)-x_array(1)
i=terrain_index
status_nc = NF90_INQ_VARID (id_nc, trim(var_name_terrain_nc(i)), var_id_nc)
status_nc = NF90_GET_VAR (id_nc, var_id_nc, array(:,:),start=(/dim_start_nc(1:2)/), count=(/dim_length_nc(1:2)/))
!write(*,*) array
write(unit_logfile,'(3A,2f12.4)') ' ',trim(var_name_terrain_nc(i)),' (min, max): ' &
,minval(array(:,:)) &
,maxval(array(:,:))
status_nc = NF90_CLOSE (id_nc)
array=max(0.,array)
write(unit_logfile,'(2a10,1a12)')'ncols','nrows','cellsize'
write(unit_logfile,'(2i10,1f12.1)')ncols_sub,nrows_sub,cellsize_sub
write(unit_logfile,'(a,f6.1,a,f6.1,a,f6.1)') ' Min array val: ',minval(array),' Max array val: ',maxval(array),' Mean array val: ',sum(array)/ncols_sub/nrows_sub
deallocate (var1d_nc)
end subroutine read_netcdf4_terrain_file