diffwrf_netcdf.F
References to this file elsewhere.
1 module read_util_module
2
3 #ifdef crayx1
4 #define iargc ipxfargc
5 #endif
6
7 contains
8
9 #ifdef crayx1
10 subroutine getarg(i, harg)
11 implicit none
12 character(len=*) :: harg
13 integer :: ierr, ilen, i
14
15 call pxfgetarg(i, harg, ilen, ierr)
16 return
17 end subroutine getarg
18 #endif
19
20 subroutine arguments(v2file, lmore)
21 implicit none
22 character(len=*) :: v2file
23 character(len=120) :: harg
24 logical :: lmore
25
26 integer :: ierr, i, numarg
27 integer, external :: iargc
28
29 numarg = iargc()
30
31 i = 1
32 lmore = .false.
33
34 do while ( i < numarg)
35 call getarg(i, harg)
36 print*, 'harg = ', trim(harg)
37
38 if (harg == "-v") then
39 i = i + 1
40 lmore = .true.
41 elseif (harg == "-h") then
42 call help
43 endif
44
45 enddo
46
47 call getarg(i,harg)
48 v2file = harg
49 end subroutine arguments
50
51 subroutine help
52 implicit none
53 character(len=120) :: cmd
54 call getarg(0, cmd)
55
56 write(*,'(/,"Usage: ", A, " [-v] v2file ")') trim(cmd)
57 write(*,'(8x, "-v : Print extra info")')
58 write(*,'(8x, "v3file : MM5v3 file name to read.")')
59 write(*,'(8x, "-h : print this help message and exit.",/)')
60 stop
61 end subroutine help
62 end module read_util_module
63
64
65
66 program readv3
67 use wrf_data
68 use read_util_module
69 implicit none
70 #include "wrf_status_codes.h"
71 #include "netcdf.inc"
72 character(len=255) :: flnm
73 character(len=255) :: flnm2
74 character(len=120) :: arg3
75 character(len=19) :: DateStr
76 character(len=19) :: DateStr2
77 character(len=31) :: VarName
78 character(len=31) :: VarName2
79 integer dh1, dh2
80
81 integer :: flag, flag2
82 integer :: iunit, iunit2
83
84 integer :: i,j,k
85 integer :: levlim
86 integer :: cross
87 integer :: ndim, ndim2
88 integer :: WrfType, WrfType2
89 real :: time, time2
90 real*8 :: a, b
91 real*8 :: sumE, sum1, sum2, diff1, diff2, serr, perr, rmse, rms1, rms2, tmp1, tmp2
92 integer digits,d1, d2
93 integer, dimension(4) :: start_index, end_index, start_index2, end_index2
94 integer , Dimension(3) :: MemS,MemE,PatS,PatE
95 character (len= 4) :: staggering, staggering2
96 character (len= 3) :: ordering, ordering2, ord
97 character (len=24) :: start_date, start_date2
98 character (len=24) :: current_date, current_date2
99 character (len=31) :: name, name2, tmpname
100 character (len=25) :: units, units2
101 character (len=46) :: description, description2
102
103 character (len=80), dimension(3) :: dimnames
104 character (len=80) :: SysDepInfo
105
106 integer :: l, n
107 integer :: ikdiffs, ifdiffs
108
109 real, allocatable, dimension(:,:,:,:) :: data,data2
110
111 integer :: ierr, ierr2, ier, ier2, Status, Status_next_time, Status_next_time2, Status_next_var, Status_next_var_2
112
113 logical :: newtime = .TRUE.
114 logical :: justplot, efound
115
116 integer, external :: iargc
117 logical, external :: iveceq
118
119 levlim = -1
120
121 call ext_ncd_ioinit(SysDepInfo,Status)
122 call set_wrf_debug_level ( 1 )
123
124
125 Justplot = .false.
126 ! get arguments
127 if ( iargc() .ge. 2 ) then
128 call getarg(1,flnm)
129 call getarg(2,flnm2)
130 ierr = 0
131 call ext_ncd_open_for_read( trim(flnm), 0, 0, "", dh1, Status)
132 if ( Status /= 0 ) then
133 print*,'error opening ',flnm, ' Status = ', Status ; stop
134 endif
135 call ext_ncd_open_for_read( trim(flnm2), 0, 0, "", dh2, Status)
136 if ( Status /= 0 ) go to 923
137 goto 924
138 923 continue
139
140 ! bounce here if second name is not openable -- this would mean that
141 ! it is a field name instead.
142
143 print*,'could not open ',flnm2
144 name = flnm2
145 Justplot = .true.
146 924 continue
147 if ( iargc() .eq. 3 ) then
148 call getarg(3,arg3)
149 read(arg3,*)levlim
150 print*,'LEVLIM = ',LEVLIM
151 endif
152 else
153 print*,'Usage: command file1 file2'
154 stop
155 endif
156
157 print*,'Just plot ',Justplot
158
159 if ( Justplot ) then
160 print*, 'flnm = ', trim(flnm)
161
162 call ext_ncd_get_next_time(dh1, DateStr, Status_next_time)
163
164 DO WHILE ( Status_next_time .eq. 0 )
165 write(*,*)'Next Time ',TRIM(Datestr)
166 call ext_ncd_get_next_var (dh1, VarName, Status_next_var)
167 DO WHILE ( Status_next_var .eq. 0 )
168 ! write(*,*)'Next Var |',TRIM(VarName),'|'
169
170 start_index = 1
171 end_index = 1
172 call ext_ncd_get_var_info (dh1,VarName,ndim,ordering,staggering,start_index,end_index, WrfType, ierr )
173 if(WrfType /= WRF_REAL .AND. WrfType /= WRF_DOUBLE) then
174 call ext_ncd_get_next_var (dh1, VarName, Status_next_var)
175 cycle
176 endif
177 write(*,'(A9,1x,I1,3(1x,I5),1x,A,1x,A)')&
178 VarName, ndim, end_index(1), end_index(2), end_index(3), &
179 trim(ordering), trim(DateStr)
180
181 if ( VarName .eq. name ) then
182 write(*,*)'Writing fort.88 file for ', trim(name)
183
184 allocate(data(end_index(1), end_index(2), end_index(3), 1))
185
186 if ( ndim .eq. 3 ) then
187 ord = 'XYZ'
188 else if ( ndim .eq. 2 ) then
189 ord = 'XY'
190 else if ( ndim .eq. 1 ) then
191 ord = 'Z'
192 else if ( ndim .eq. 0 ) then
193 ord = '0'
194 endif
195
196 call ext_ncd_read_field(dh1,DateStr,TRIM(name),data,WRF_REAL,0,0,0,ord, &
197 staggering, dimnames , &
198 start_index,end_index, & !dom
199 start_index,end_index, & !mem
200 start_index,end_index, & !pat
201 ierr)
202
203 if ( ierr/=0 ) then
204 write(*,*)'error reading data record'
205 write(*,*)' ndim = ', ndim
206 write(*,*)' end_index(1) ',end_index(1)
207 write(*,*)' end_index(2) ',end_index(2)
208 write(*,*)' end_index(3) ',end_index(3)
209 endif
210
211
212 #if 0
213 ! uncomment this to have the code give i-slices
214 do i = 1, end_index(1)
215 if ( levlim .eq. -1 .or. i .eq. levlim ) then
216 write(88,*)end_index(2),end_index(3),' ',trim(name),' ',k,' time ',TRIM(Datestr)
217 do k = start_index(3), end_index(3)
218 do j = 1, end_index(2)
219 write(88,*) data(i,j,k,1)
220 enddo
221 enddo
222 endif
223 enddo
224 #else
225 ! give k-slices
226 do k = start_index(3), end_index(3)
227 if ( levlim .eq. -1 .or. k .eq. levlim ) then
228 write(88,*)end_index(1),end_index(2),' ',trim(name),' ',k,' time ',TRIM(Datestr)
229 do j = 1, end_index(2)
230 do i = 1, end_index(1)
231 write(88,*) data(i,j,k,1)
232 enddo
233 enddo
234 endif
235 enddo
236 #endif
237
238 deallocate(data)
239 endif
240 call ext_ncd_get_next_var (dh1, VarName, Status_next_var)
241 enddo
242 call ext_ncd_get_next_time(dh1, DateStr, Status_next_time)
243 enddo
244 else
245 write (6,FMT='(4A)') 'Diffing ',trim(flnm),' ',trim(flnm2)
246
247 call ext_ncd_get_next_time(dh1, DateStr, Status_next_time)
248 call ext_ncd_get_next_time(dh2, DateStr2, Status_next_time2)
249
250 IF ( DateStr .NE. DateStr2 ) THEN
251 print*,'They differ big time. Dates do not match'
252 print*,' ',flnm,' ',DateStr
253 print*,' ',flnm2,' ',DateStr2
254 Status_next_time = 1
255 ENDIF
256
257 DO WHILE ( Status_next_time .eq. 0 .AND. Status_next_time2 .eq. 0 )
258 write(*,*)'Next Time ',TRIM(Datestr)
259 print 76
260 call ext_ncd_get_next_var (dh1, VarName, Status_next_var)
261 DO WHILE ( Status_next_var .eq. 0 )
262 ! write(*,*)'Next Var |',TRIM(VarName),'|'
263
264 start_index = 1
265 end_index = 1
266 start_index2 = 1
267 end_index2 = 1
268
269 call ext_ncd_get_var_info (dh1,VarName,ndim,ordering,staggering,start_index,end_index, WrfType, ierr )
270 call ext_ncd_get_var_info (dh2,VarName,ndim2,ordering2,staggering2,start_index2,end_index2, WrfType2, ierr )
271 IF ( ierr /= 0 ) THEN
272 write(*,*)'Big difference: ',VarName,' not found in ',flnm2
273 GOTO 1234
274 ENDIF
275 IF ( ndim /= ndim2 ) THEN
276 write(*,*)'Big difference: Number of dimensions for ',Varname,' differs in ',flnm2,'(',ndim,') /= (',ndim2
277 GOTO 1234
278 ENDIF
279 IF ( WrfType /= WrfType2 ) THEN
280 write(*,*)'Big difference: The types do not match'
281 GOTO 1234
282 ENDIF
283 if( WrfType == WRF_REAL) then
284 DO i = 1, ndim
285 IF ( end_index(i) /= end_index2(i) ) THEN
286 write(*,*)'Big difference: dim ',i,' lengths differ for ',Varname,' differ in ',flnm2
287 GOTO 1234
288 ENDIF
289 ENDDO
290 DO i = ndim+1,3
291 start_index(i) = 1
292 end_index(i) = 1
293 start_index2(i) = 1
294 end_index2(i) = 1
295 ENDDO
296
297 ! write(*,'(A9,1x,I1,3(1x,I3),1x,A,1x,A)')&
298 ! VarName, ndim, end_index(1), end_index(2), end_index(3), &
299 ! trim(ordering), trim(DateStr)
300
301 allocate(data (end_index(1), end_index(2), end_index(3), 1))
302 allocate(data2(end_index(1), end_index(2), end_index(3), 1))
303
304 if ( ndim .eq. 3 ) then
305 ord = 'XYZ'
306 else if ( ndim .eq. 2 ) then
307 ord = 'XY'
308 else if ( ndim .eq. 1 ) then
309 ord = 'Z'
310 else if ( ndim .eq. 0 ) then
311 ord = '0'
312 endif
313
314 call ext_ncd_read_field(dh1,DateStr,TRIM(VarName),data,WRF_REAL,0,0,0,ord,&
315 staggering, dimnames , &
316 start_index,end_index, & !dom
317 start_index,end_index, & !mem
318 start_index,end_index, & !pat
319 ierr)
320
321 IF ( ierr /= 0 ) THEN
322 write(*,*)'Error reading ',Varname,' from ',flnm
323 write(*,*)' ndim = ', ndim
324 write(*,*)' end_index(1) ',end_index(1)
325 write(*,*)' end_index(2) ',end_index(2)
326 write(*,*)' end_index(3) ',end_index(3)
327 ENDIF
328 call ext_ncd_read_field(dh2,DateStr,TRIM(VarName),data2,WRF_REAL,0,0,0,ord,&
329 staggering, dimnames , &
330 start_index,end_index, & !dom
331 start_index,end_index, & !mem
332 start_index,end_index, & !pat
333 ierr)
334 IF ( ierr /= 0 ) THEN
335 write(*,*)'Error reading ',Varname,' from ',flnm2
336 write(*,*)' ndim = ', ndim
337 write(*,*)' end_index(1) ',end_index(1)
338 write(*,*)' end_index(2) ',end_index(2)
339 write(*,*)' end_index(3) ',end_index(3)
340 ENDIF
341
342 IFDIFFS=0
343 sumE = 0.0
344 sum1 = 0.0
345 sum2 = 0.0
346 diff1 = 0.0
347 diff2 = 0.0
348 n = 0
349 DO K = 1,end_index(3)-start_index(3)+1
350 IF (LEVLIM.EQ.-1.OR.K.EQ.LEVLIM.OR.NDIM.eq.2) THEN
351 cross = 0
352 IKDIFFS = 0
353 do i = 1, end_index(1)-cross
354 do j = 1, end_index(2)-cross
355 a = data(I,J,K,1)
356 b = data2(I,J,K,1)
357 ! borrowed from Thomas Oppe's comp program
358 sumE = sumE + ( a - b ) * ( a - b )
359 sum1 = sum1 + a * a
360 sum2 = sum2 + b * b
361 diff1 = max ( diff1 , abs ( a - b ) )
362 diff2 = max ( diff2 , abs ( b ) )
363 n = n + 1
364 IF (a .ne. b) then
365 IKDIFFS = IKDIFFS + 1
366 IFDIFFS = IFDIFFS + 1
367 ENDIF
368 ENDDO
369 ENDDO
370 ENDIF
371 enddo
372 rmsE = sqrt ( sumE / dble( n ) )
373 rms1 = sqrt ( sum1 / dble( n ) )
374 rms2 = sqrt ( sum2 / dble( n ) )
375 serr = 0.0
376 IF ( sum2 .GT. 0.0d0 ) THEN
377 serr = sqrt ( sumE / sum2 )
378 ELSE
379 IF ( sumE .GT. 0.0d0 ) serr = 1.0
380 ENDIF
381 perr = 0.0
382 IF ( diff2 .GT. 0.0d0 ) THEN
383 perr = diff1/diff2
384 ELSE
385 IF ( diff1 .GT. 0.0d0 ) perr = 1.0
386 ENDIF
387
388 IF ( rms1 - rms2 .EQ. 0.0d0 ) THEN
389 digits = 15
390 ELSE
391 IF ( rms2 .NE. 0 ) THEN
392 tmp1 = 1.0d0/( ( abs( rms1 - rms2 ) ) / rms2 )
393 IF ( tmp1 .NE. 0 ) THEN
394 digits = log10(tmp1)
395 ENDIF
396 ENDIF
397 ENDIF
398
399 IF (IFDIFFS .NE. 0 ) THEN
400 ! create the fort.88 and fort.98 files because regression scripts will
401 ! look for these to see if there were differences.
402 write(88,*)trim(varname)
403 write(98,*)trim(varname)
404 PRINT 77,trim(varname), IFDIFFS, ndim, rms1, rms2, digits, rmsE, perr
405 76 FORMAT (5x,'Field ',2x,'Ndifs',4x,'Dims ',6x,'RMS (1)',12x,'RMS (2)',5x,'DIGITS',4x,'RMSE',5x,'pntwise max')
406 77 FORMAT ( A10,1x,I9,2x,I3,1x,e18.10,1x,e18.10,1x,i3,1x,e12.4,1x,e12.4 )
407 ENDIF
408 deallocate(data)
409 deallocate(data2)
410
411 endif
412 1234 CONTINUE
413 call ext_ncd_get_next_var (dh1, VarName, Status_next_var)
414 enddo
415 call ext_ncd_get_next_time(dh1, DateStr, Status_next_time)
416 call ext_ncd_get_next_time(dh2, DateStr2, Status_next_time2)
417 IF ( DateStr .NE. DateStr2 ) THEN
418 print*,'They differ big time. Dates do not match'
419 print*,'They differ big time. Dates do not match'
420 print*,' ',flnm,' ',DateStr
421 print*,' ',flnm2,' ',DateStr2
422 Status_next_time = 1
423 ENDIF
424 enddo
425
426 endif
427
428 end program readv3
429
430 logical function iveceq( a, b, n )
431 implicit none
432 integer n
433 integer a(n), b(n)
434 integer i
435 iveceq = .true.
436 do i = 1,n
437 if ( a(i) .ne. b(i) ) iveceq = .false.
438 enddo
439 return
440 end function iveceq
441
442 ! stubs for routines called by module_wrf_error (used by netcdf implementation of IO api)
443 SUBROUTINE wrf_abort
444 STOP
445 END SUBROUTINE wrf_abort
446
447 SUBROUTINE get_current_time_string( time_str )
448 CHARACTER(LEN=*), INTENT(OUT) :: time_str
449 time_str = ''
450 END SUBROUTINE get_current_time_string
451
452 SUBROUTINE get_current_grid_name( grid_str )
453 CHARACTER(LEN=*), INTENT(OUT) :: grid_str
454 grid_str = ''
455 END SUBROUTINE get_current_grid_name
456