1 module read_util_module
5 subroutine arguments(v2file, lmore)
7 character(len=*) :: v2file
8 character(len=120) :: harg
11 integer :: ierr, i, numarg
13 numarg = command_argument_count()
18 do while ( i < numarg)
19 call get_command_argument(number=i, value=harg)
20 print*, 'harg = ', trim(harg)
22 if (harg == "-v") then
25 elseif (harg == "-h") then
31 call get_command_argument(number=i, value=harg)
33 end subroutine arguments
37 character(len=120) :: cmd
38 call get_command_argument(number=0, value=cmd)
40 write(*,'(/,"Usage: ", A, " [-v] v2file ")') trim(cmd)
41 write(*,'(8x, "-v : Print extra info")')
42 write(*,'(8x, "v3file : MM5v3 file name to read.")')
43 write(*,'(8x, "-h : print this help message and exit.",/)')
46 end module read_util_module
54 #include "wrf_status_codes.h"
56 character(len=255) :: flnm
57 character(len=255) :: flnm2
58 character(len=120) :: arg3
59 character(len=19) :: DateStr
60 character(len=19) :: DateStr2
61 character(len=31) :: VarName
62 character(len=31) :: VarName2
65 integer :: flag, flag2
66 integer :: iunit, iunit2
71 integer :: ndim, ndim2
72 integer :: WrfType, WrfType2
75 real*8 :: sumE, sum1, sum2, diff1, diff2, serr, perr, rmse, rms1, rms2, tmp1, tmp2
77 integer, dimension(4) :: start_index, end_index, start_index2, end_index2
78 integer , Dimension(3) :: MemS,MemE,PatS,PatE
79 character (len= 4) :: staggering, staggering2
80 character (len= 3) :: ordering, ordering2, ord
81 character (len=24) :: start_date, start_date2
82 character (len=24) :: current_date, current_date2
83 character (len=31) :: name, name2, tmpname
84 character (len=25) :: units, units2
85 character (len=46) :: description, description2
87 character (len=80), dimension(3) :: dimnames
88 character (len=80) :: SysDepInfo
90 logical :: first, searchcoords
91 integer :: l, n, ntimes
92 integer :: ikdiffs, ifdiffs
93 integer :: icenter, prev_icenter, jcenter, prev_jcenter,ntries
94 real :: searchlat, searchlong
96 real, allocatable, dimension(:,:,:,:) :: data, data2
97 integer, allocatable, dimension(:,:,:,:) :: idata,idata2
98 real, allocatable, dimension(:,:) :: xlat,xlong
100 integer :: ierr, ierr2, ier, ier2, Status, Status_next_time, Status_next_time2, Status_next_var, Status_next_var_2
103 logical :: newtime = .TRUE.
104 logical :: justplot, efound
106 logical, external :: iveceq
110 call ext_ncd_ioinit(SysDepInfo,Status)
111 call set_wrf_debug_level ( 1 )
113 nargs = command_argument_count()
116 searchcoords = .false.
118 if ( nargs .ge. 2 ) then
119 call get_command_argument(number=1, value=flnm)
120 call get_command_argument(number=2, value=flnm2)
121 IF ( flnm2(1:4) .EQ. '-lat' ) THEN
122 print*,'reading ',TRIM(flnm2(5:))
123 read(flnm2(5:),*)searchlat
124 call get_command_argument(number=3, value=flnm2)
125 IF ( flnm2(1:5) .EQ. '-long' ) THEN
126 print*,'reading ',TRIM(flnm2(6:))
127 read(flnm2(6:),*)searchlong
129 write(*,*)'missing -long argument (no spaces after -lat or -long, either)'
134 searchcoords = .true.
135 call ext_ncd_open_for_read( trim(flnm), 0, 0, "", dh1, Status)
139 call ext_ncd_open_for_read( trim(flnm), 0, 0, "", dh1, Status)
140 if ( Status /= 0 ) then
141 print*,'error opening ',flnm, ' Status = ', Status ; stop
143 call ext_ncd_open_for_read( trim(flnm2), 0, 0, "", dh2, Status)
144 if ( Status /= 0 ) go to 923
148 ! bounce here if second name is not openable -- this would mean that
149 ! it is a field name instead.
151 print*,'could not open ',flnm2
155 if ( nargs .eq. 3 ) then
156 call get_command_argument(number=3, value=arg3)
158 print*,'LEVLIM = ',LEVLIM
161 print*,'Usage: command file1 file2'
165 print*,'Just plot ',Justplot
168 print*, 'flnm = ', trim(flnm)
171 call ext_ncd_get_next_time(dh1, DateStr, Status_next_time)
174 DO WHILE ( Status_next_time .eq. 0 )
175 write(*,*)'Next Time ',TRIM(Datestr)
177 call ext_ncd_get_next_var (dh1, VarName, Status_next_var)
178 DO WHILE ( Status_next_var .eq. 0 )
179 ! write(*,*)'Next Var |',TRIM(VarName),'|'
183 call ext_ncd_get_var_info (dh1,VarName,ndim,ordering,staggering,start_index,end_index, WrfType, ierr )
184 if(WrfType /= WRF_REAL .AND. WrfType /= WRF_DOUBLE .AND. WrfType /= WRF_INTEGER) then
185 call ext_ncd_get_next_var (dh1, VarName, Status_next_var)
188 IF ( .NOT. searchcoords ) THEN
189 write(*,'(A9,1x,I1,3(1x,I5),1x,A,1x,A)')&
190 VarName, ndim, end_index(1), end_index(2), end_index(3), &
191 trim(ordering), trim(DateStr)
194 if ( VarName .eq. name .OR. TRIM(VarName) .EQ. 'XLAT' .OR. TRIM(VarName) .EQ. 'XLONG' ) then
195 write(*,*)'Writing fort.88 file for ', trim(name)
197 allocate(data(end_index(1), end_index(2), end_index(3), 1))
199 if ( ndim .eq. 3 ) then
201 else if ( ndim .eq. 2 ) then
203 else if ( ndim .eq. 1 ) then
205 else if ( ndim .eq. 0 ) then
209 call ext_ncd_read_field(dh1,DateStr,TRIM(VarName),data,WRF_REAL,0,0,0,ord, &
210 staggering, dimnames , &
211 start_index,end_index, & !dom
212 start_index,end_index, & !mem
213 start_index,end_index, & !pat
217 write(*,*)'error reading data record'
218 write(*,*)' ndim = ', ndim
219 write(*,*)' end_index(1) ',end_index(1)
220 write(*,*)' end_index(2) ',end_index(2)
221 write(*,*)' end_index(3) ',end_index(3)
224 write(*,*)'name: ',TRIM(VarName)
225 IF ( TRIM(VarName) .EQ. 'XLAT' .AND. .NOT. ALLOCATED(xlat)) THEN
226 write(*,*)'allocating xlat'
227 ALLOCATE(xlat(end_index(1), end_index(2)))
230 IF ( TRIM(VarName) .EQ. 'XLONG' .AND. .NOT. ALLOCATED(xlong)) THEN
231 write(*,*)'allocating xlong'
232 ALLOCATE(xlong(end_index(1), end_index(2)))
233 xlong = data(:,:,1,1)
237 if ( VarName .eq. name ) then
239 ! uncomment this to have the code give i-slices
240 do i = 1, end_index(1)
241 if ( levlim .eq. -1 .or. i .eq. levlim ) then
242 write(88,*)end_index(2),end_index(3),' ',trim(name),' ',k,' time ',TRIM(Datestr)
243 do k = start_index(3), end_index(3)
244 do j = 1, end_index(2)
245 write(88,*) data(i,j,k,1)
252 do k = start_index(3), end_index(3)
253 if ( levlim .eq. -1 .or. k .eq. levlim ) then
254 write(88,*)end_index(1),end_index(2),' ',trim(name),' ',k,' time ',TRIM(Datestr)
255 do j = 1, end_index(2)
256 do i = 1, end_index(1)
257 write(88,*) data(i,j,k,1)
267 call ext_ncd_get_next_var (dh1, VarName, Status_next_var)
268 IF ( ntimes .EQ. 1 .AND. ALLOCATED(xlong) .AND. ALLOCATED(xlat) .AND. first ) THEN
275 DO WHILE ( ntries .LT. 10 .AND. (icenter .NE. prev_icenter .OR. jcenter .NE. prev_jcenter ))
276 prev_icenter = icenter
277 prev_jcenter = jcenter
278 DO j = start_index(2), end_index(2)-1
279 IF ( xlat(icenter,j) .LE. searchlat .AND. searchlat .LT. xlat(icenter,j+1) ) THEN
281 !write(*,*)'xlat ',ntries,icenter,jcenter,xlat(icenter,j),searchlat
285 DO i = start_index(1), end_index(1)-1
286 IF ( xlong(i,jcenter) .LE. searchlong .AND. searchlong .LT. xlong(i+1,jcenter)) THEN
288 !write(*,*)'xlon ',ntries,icenter,jcenter,xlong(i,jcenter),searchlong
294 write(*,*)'Lon ',searchlong,' Lat ',searchlat,' : ',icenter,jcenter
295 write(*,*)'Coordinates at that point ',xlong(icenter,jcenter),xlat(icenter,jcenter)
296 write(*,*)'Coordinates at next point ',xlong(icenter+1,jcenter+1),xlat(icenter+1,jcenter+1)
297 write(*,*)'Ntries : ',ntries
298 if ( ntries .GE. 10 ) write(*,*)'max tries exceeded. Probably did not find'
301 call ext_ncd_get_next_time(dh1, DateStr, Status_next_time)
304 write (6,FMT='(4A)') 'Diffing ',trim(flnm),' ',trim(flnm2)
306 call ext_ncd_get_next_time(dh1, DateStr, Status_next_time)
307 call ext_ncd_get_next_time(dh2, DateStr2, Status_next_time2)
309 IF ( DateStr .NE. DateStr2 ) THEN
310 print*,'They differ big time. Dates do not match'
311 print*,' ',flnm,' ',DateStr
312 print*,' ',flnm2,' ',DateStr2
316 DO WHILE ( Status_next_time .eq. 0 .AND. Status_next_time2 .eq. 0 )
317 write(*,*)'Next Time ',TRIM(Datestr)
319 call ext_ncd_get_next_var (dh1, VarName, Status_next_var)
320 DO WHILE ( Status_next_var .eq. 0 )
321 ! write(*,*)'Next Var |',TRIM(VarName),'|'
328 call ext_ncd_get_var_info (dh1,VarName,ndim,ordering,staggering,start_index,end_index, WrfType, ierr )
329 call ext_ncd_get_var_info (dh2,VarName,ndim2,ordering2,staggering2,start_index2,end_index2, WrfType2, ierr )
330 IF ( ierr /= 0 ) THEN
331 write(*,*)'Big difference: ',VarName,' not found in ',flnm2
334 IF ( ndim /= ndim2 ) THEN
335 write(*,*)'Big difference: Number of dimensions for ',Varname,' differs in ',flnm2,'(',ndim,') /= (',ndim2
338 IF ( WrfType /= WrfType2 ) THEN
339 write(*,*)'Big difference: The types do not match'
342 if( WrfType == WRF_REAL .OR. WrfType == WRF_INTEGER ) then
344 IF ( end_index(i) /= end_index2(i) ) THEN
345 write(*,*)'Big difference: dim ',i,' lengths differ for ',Varname,' differ in ',flnm2
356 ! write(*,'(A9,1x,I1,3(1x,I3),1x,A,1x,A)')&
357 ! VarName, ndim, end_index(1), end_index(2), end_index(3), &
358 ! trim(ordering), trim(DateStr)
360 if (WrfType == WRF_REAL ) then
361 allocate( data (end_index(1), end_index(2), end_index(3), 1))
362 allocate( data2(end_index(1), end_index(2), end_index(3), 1))
363 else if (WrfType == WRF_INTEGER ) then
364 allocate(idata (end_index(1), end_index(2), end_index(3), 1))
365 allocate(idata2(end_index(1), end_index(2), end_index(3), 1))
368 if ( ndim .eq. 3 ) then
370 else if ( ndim .eq. 2 ) then
372 else if ( ndim .eq. 1 ) then
374 else if ( ndim .eq. 0 ) then
378 if (WrfType == WRF_REAL ) then
379 call ext_ncd_read_field(dh1,DateStr,TRIM(VarName),data,WRF_REAL,0,0,0,ord,&
380 staggering, dimnames , &
381 start_index,end_index, & !dom
382 start_index,end_index, & !mem
383 start_index,end_index, & !pat
386 IF ( ierr /= 0 ) THEN
387 write(*,*)'Error reading ',Varname,' from ',flnm
388 write(*,*)' ndim = ', ndim
389 write(*,*)' end_index(1) ',end_index(1)
390 write(*,*)' end_index(2) ',end_index(2)
391 write(*,*)' end_index(3) ',end_index(3)
393 call ext_ncd_read_field(dh2,DateStr,TRIM(VarName),data2,WRF_REAL,0,0,0,ord,&
394 staggering, dimnames , &
395 start_index,end_index, & !dom
396 start_index,end_index, & !mem
397 start_index,end_index, & !pat
399 IF ( ierr /= 0 ) THEN
400 write(*,*)'Error reading ',Varname,' from ',flnm2
401 write(*,*)' ndim = ', ndim
402 write(*,*)' end_index(1) ',end_index(1)
403 write(*,*)' end_index(2) ',end_index(2)
404 write(*,*)' end_index(3) ',end_index(3)
406 else if (WrfType == WRF_INTEGER ) then
407 call ext_ncd_read_field(dh1,DateStr,TRIM(VarName),idata,WRF_INTEGER,0,0,0,ord,&
408 staggering, dimnames , &
409 start_index,end_index, & !dom
410 start_index,end_index, & !mem
411 start_index,end_index, & !pat
414 IF ( ierr /= 0 ) THEN
415 write(*,*)'Error reading ',Varname,' from ',flnm
416 write(*,*)' ndim = ', ndim
417 write(*,*)' end_index(1) ',end_index(1)
418 write(*,*)' end_index(2) ',end_index(2)
419 write(*,*)' end_index(3) ',end_index(3)
421 call ext_ncd_read_field(dh2,DateStr,TRIM(VarName),idata2,WRF_INTEGER,0,0,0,ord,&
422 staggering, dimnames , &
423 start_index,end_index, & !dom
424 start_index,end_index, & !mem
425 start_index,end_index, & !pat
427 IF ( ierr /= 0 ) THEN
428 write(*,*)'Error reading ',Varname,' from ',flnm2
429 write(*,*)' ndim = ', ndim
430 write(*,*)' end_index(1) ',end_index(1)
431 write(*,*)' end_index(2) ',end_index(2)
432 write(*,*)' end_index(3) ',end_index(3)
443 DO K = 1,end_index(3)-start_index(3)+1
444 IF (LEVLIM.EQ.-1.OR.K.EQ.LEVLIM.OR.NDIM.eq.2) THEN
448 if (WrfType == WRF_REAL ) then
449 do i = 1, end_index(1)-cross
450 do j = 1, end_index(2)-cross
453 ! borrowed from Thomas Oppe's comp program
454 sumE = sumE + ( a - b ) * ( a - b )
457 diff1 = max ( diff1 , abs ( a - b ) )
458 diff2 = max ( diff2 , abs ( b ) )
461 IKDIFFS = IKDIFFS + 1
462 IFDIFFS = IFDIFFS + 1
466 else if (WrfType == WRF_INTEGER ) then
467 do i = 1, end_index(1)-cross
468 do j = 1, end_index(2)-cross
471 ! borrowed from Thomas Oppe's comp program
472 sumE = sumE + ( a - b ) * ( a - b )
475 diff1 = max ( diff1 , abs ( a - b ) )
476 diff2 = max ( diff2 , abs ( b ) )
479 IKDIFFS = IKDIFFS + 1
480 IFDIFFS = IFDIFFS + 1
488 rmsE = sqrt ( sumE / dble( n ) )
489 rms1 = sqrt ( sum1 / dble( n ) )
490 rms2 = sqrt ( sum2 / dble( n ) )
492 IF ( sum2 .GT. 0.0d0 ) THEN
493 serr = sqrt ( sumE / sum2 )
495 IF ( sumE .GT. 0.0d0 ) serr = 1.0
498 IF ( diff2 .GT. 0.0d0 ) THEN
501 IF ( diff1 .GT. 0.0d0 ) perr = 1.0
504 IF ( rms1 - rms2 .EQ. 0.0d0 ) THEN
507 IF ( rms2 .NE. 0 ) THEN
508 tmp1 = 1.0d0/( ( abs( rms1 - rms2 ) ) / rms2 )
509 IF ( tmp1 .NE. 0 ) THEN
515 IF (IFDIFFS .NE. 0 ) THEN
516 ! create the fort.88 and fort.98 files because regression scripts will
517 ! look for these to see if there were differences.
518 write(88,*)trim(varname)
519 write(98,*)trim(varname)
520 PRINT 77,trim(varname), IFDIFFS, ndim, rms1, rms2, digits, rmsE, perr
521 76 FORMAT (5x,'Field ',2x,'Ndifs',4x,'Dims ',6x,'RMS (1)',12x,'RMS (2)',5x,'DIGITS',4x,'RMSE',5x,'pntwise max')
522 77 FORMAT ( A10,1x,I9,2x,I3,1x,e18.10,1x,e18.10,1x,i3,1x,e12.4,1x,e12.4 )
524 if (WrfType == WRF_REAL ) then
527 else if (WrfType == WRF_INTEGER ) then
534 call ext_ncd_get_next_var (dh1, VarName, Status_next_var)
536 call ext_ncd_get_next_time(dh1, DateStr, Status_next_time)
537 call ext_ncd_get_next_time(dh2, DateStr2, Status_next_time2)
538 IF ( DateStr .NE. DateStr2 ) THEN
539 print*,'They differ big time. Dates do not match'
540 print*,'They differ big time. Dates do not match'
541 print*,' ',flnm,' ',DateStr
542 print*,' ',flnm2,' ',DateStr2
551 logical function wrf_dm_on_monitor()
552 wrf_dm_on_monitor=.true.
553 end function wrf_dm_on_monitor
555 logical function iveceq( a, b, n )
562 if ( a(i) .ne. b(i) ) iveceq = .false.
567 ! stubs for routines called by module_wrf_error (used by netcdf implementation of IO api)
570 END SUBROUTINE wrf_abort
572 SUBROUTINE get_current_time_string( time_str )
573 CHARACTER(LEN=*), INTENT(OUT) :: time_str
575 END SUBROUTINE get_current_time_string
577 SUBROUTINE get_current_grid_name( grid_str )
578 CHARACTER(LEN=*), INTENT(OUT) :: grid_str
580 END SUBROUTINE get_current_grid_name