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