-
Notifications
You must be signed in to change notification settings - Fork 22
/
Copy pathReadMatrix.F90
295 lines (196 loc) · 7.2 KB
/
ReadMatrix.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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
#include "ppdef.h"
!--------------------------------------------------------
!----- subroutine readMatrix ----------------------------
!--------------------------------------------------------
! reads the data and mask files and creates 'matrix',
! containing data and missing data, no land values
!
!--------------------------------------------------------
subroutine ReadMatrix(sstfile,matrix,m,n,valex,maskfile,norma,imax,jmax,first,fileMean,fileStd)
use ufileformat
use initfile
implicit none
character(len=200),intent(in) :: sstfile(:),maskfile(:)
real,pointer :: matrix(:,:)
integer,intent(in) :: norma
integer,intent(out) :: m,n
real,intent(out) :: valex
integer, pointer :: imax(:),jmax(:),first(:)
real,pointer :: fileMean(:),fileStd(:)
integer, allocatable :: land(:),sea(:),nlines(:)
real,pointer :: file(:,:,:),mask(:,:)
real, parameter :: valexc = 9999.
integer :: i,j,s,t,k,d,q,nbvars,testN
real :: Y,maskvalex,testmean
real,pointer :: valexN
character(len=200) :: initfilename
interface
subroutine stat(file,mask,mean,stddev)
real,pointer :: file(:,:,:),mask(:,:)
real,intent(out) :: mean,stddev
end subroutine stat
end interface
interface
subroutine norm(file,mean,stddev)
real,pointer :: file(:,:,:)
real,intent(in) :: mean,stddev
end subroutine norm
end interface
nbvars=size(sstfile)
allocate(imax(nbvars),jmax(nbvars),land(nbvars),sea(nbvars),nlines(nbvars),first(nbvars),fileMean(nbvars),fileStd(nbvars))
call getarg(1,initfilename)
write(stdout,*)'initfilename: ',initfilename
!if( presentInitValue(initfilename,'mask')) then
!if(maskfile = 'nomaskfile') then
! write(stdout,*)'no mask for reading '
!endif
do q=1,nbvars
call uload(sstfile(q),file, valex)
if( presentInitValue(initfilename,'mask')) then
call uload(maskfile(q),mask,maskvalex)
else
allocate(mask(size(file,1),size(file,2)))
mask = mask+1
endif
where (mask.eq.0) mask = valexc
!where (spread(mask,3,size(file,3)).eq.valexc) file = valex
do t=1,size(file,3)
do j=1,size(file,2)
do i=1,size(file,1)
if (mask(i,j).eq.valexc) file(i,j,t) = valexc
enddo
enddo
enddo
! ------------------------------------
! check if all input matrices
! have the same size as their masks
! ------------------------------------
if(size(file,1).ne.size(mask,1).or.size(file,2).ne.size(mask,2)) then
write(stdout,*) 'File', trim(sstfile(q)), 'does not have the same size as its mask,', trim(maskfile(q)), 'please check.'
STOP
end if
imax(q)=size(mask,1)
jmax(q)=size(mask,2)
! -----------------------------------------
! temporal dimension, same for all files!!
! -----------------------------------------
nlines(q)=size(file,3)
! -----------------------------------------
! counts land points
! -----------------------------------------
land(q)=count(mask.eq.valexc)
! -----------------------------------------
! calculates first dimension of the matrix
! (spatial points)
! -----------------------------------------
sea(q)=(imax(q)*jmax(q))-land(q)
if(q.eq.1) then
first(q)=1
else
first(q)=first(q-1)+sea(q-1)
end if
deallocate(file,mask)
end do
if(any(nlines.ne.nlines(1))) then
write(stdout,*) 'The temporal sizes of both files provided are not the same'
write(stdout,*) 'I have found the following temporal dimensions,', nlines(:)
write(stdout,*) 'Please bear in mind that this routine is constructed for the'
write(stdout,*) 'case when the temporal size of all files being reconstructed is the same.'
STOP
end if
write(stdout,*)
write(stdout,*) '********************************************************************'
write(stdout,*) 'Now some statistics about your data:'
write(stdout,*)
do q=1,nbvars
write(stdout,400)' Number of mask land points: ',land(q)
write(stdout,300)' Dimension of file ',q,': ','',imax(q),'x',jmax(q),'x',nlines(q)
write(stdout,*)
enddo
m=sum(sea(:))
allocate(matrix(m,nlines(1)))
do q=1,nbvars
call uload(sstfile(q),file, valex)
! where (file.eq.valex) file = valexc
do t=1,size(file,3)
do j=1,size(file,2)
do i=1,size(file,1)
if (file(i,j,t).eq.valex) file(i,j,t)=valexc
enddo
enddo
enddo
if( presentInitValue(initfilename,'mask')) then
call uload(maskfile(q),mask,maskvalex)
else
allocate(mask(size(file,1),size(file,2)))
mask = mask+1
endif
! where (mask.eq.0) mask = valexc
do j=1,size(mask,2)
do i=1,size(mask,1)
if (mask(i,j).eq.0) mask(i,j)=valexc
enddo
enddo
! -----------------------------------
! normalisation of the input matrices
! if norma = 1
! ------------------------------------
if(norma.eq.1) then
call stat(file,mask,fileMean(q),fileStd(q))
call norm(file,fileMean(q),fileStd(q))
!fileNorm=file
else
call stat(file,mask,fileMean(q),fileStd(q))
! fileNorm = file
! deallocate(file)
end if
!where(fileNorm.eq.9999) fileNorm=valexc
do t=1,size(file,3)
do j=1,size(file,2)
do i=1,size(file,1)
if (file(i,j,t).eq.9999) file(i,j,t)=valexc
enddo
enddo
enddo
! where (spread(mask,3,nlines(q)).eq.valexc) file = valexc
! testN = count(file.ne.valexc)
! testmean = sum(file,file.ne.valexc)/testN
! ---------------------------------------
! reshape input matrices into one matrix
! size M*N (M spatial dimension,
! N temporal dimension)
! ---------------------------------------
do t=1,nlines(q)
s = first(q)
do i=1,imax(q)
do j=1,jmax(q)
if (mask(i,j).eq.valexc) file(i,j,t) = valexc
if(mask(i,j).eq.valexc.and.file(i,j,t).ne.valexc) then
write(stdout,*) 'WARNING ',i,j
end if
if(mask(i,j).ne.valexc) then
matrix(s,t)=file(i,j,t)
s = s+1
endif
enddo
enddo
enddo
deallocate(file,mask)
call flush(stdout,istat)
! -----------------------------------
! write mean and standard deviation
! of input matrices
! ------------------------------------
write(stdout,*)
write(stdout,500)' Mean: ',fileMean(q)
write(stdout,500)' Standard deviation: ',fileStd(q)
write(stdout,*)
enddo
valex=valexc
n=nlines(1)
deallocate(land,sea,nlines)
200 format(a20,f6.2)
300 format(a26,i2,a2,a15,i6,a3,i6,a3,i6)
400 format(a30,i39)
500 format(a30,f39.4)
end subroutine ReadMatrix