Merge pull request #22 from wirc-sjsu/develop-w21
[WRF-SFIRE.git] / external / io_netcdf / diffwrf.F90
blob4c72cfed7f65f2784c399eb808ae6b15044eaa50
1 module read_util_module
3 contains
5    subroutine arguments(v2file, lmore)
6      implicit none
7      character(len=*) :: v2file
8      character(len=120) :: harg
9      logical :: lmore
10    
11      integer :: ierr, i, numarg
12    
13      numarg = command_argument_count()
14    
15      i = 1
16      lmore = .false.
17    
18      do while ( i < numarg) 
19         call get_command_argument(number=i, value=harg)
20         print*, 'harg = ', trim(harg)
21    
22         if (harg == "-v") then
23            i = i + 1
24            lmore = .true.
25         elseif (harg == "-h") then
26            call help
27         endif
28    
29      enddo
30    
31      call get_command_argument(number=i, value=harg)
32      v2file = harg
33    end subroutine arguments
34    
35    subroutine help
36      implicit none
37      character(len=120) :: cmd
38      call get_command_argument(number=0, value=cmd)
39    
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.",/)')
44      stop
45    end subroutine help
46 end module read_util_module
50  program readv3
51   use wrf_data
52   use read_util_module
53   implicit none
54 #include "wrf_status_codes.h"
55 #include "netcdf.inc"
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
63   integer dh1, dh2
65   integer :: flag, flag2
66   integer :: iunit, iunit2
68   integer :: i,j,k
69   integer :: levlim
70   integer :: cross
71   integer :: ndim, ndim2
72   integer :: WrfType, WrfType2
73   real :: time, time2
74   real*8 :: a, b
75   real*8 :: sumE, sum1, sum2, diff1, diff2, serr, perr, rmse, rms1, rms2, tmp1, tmp2
76   integer digits,d1, d2
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
101   integer :: nargs
103   logical :: newtime = .TRUE.
104   logical :: justplot, efound
106   logical, external :: iveceq
108   levlim = -1
110   call ext_ncd_ioinit(SysDepInfo,Status)
111   call set_wrf_debug_level ( 1 )
113   nargs = command_argument_count()
115   Justplot = .false.
116   searchcoords = .false.
117 ! get arguments
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
128       ELSE
129         write(*,*)'missing -long argument (no spaces after -lat or -long, either)'
130         STOP
131       ENDIF
132       nargs = 0
133       Justplot = .true.
134       searchcoords = .true.
135       call ext_ncd_open_for_read( trim(flnm), 0, 0, "", dh1, Status)
136       goto 924
137     ENDIF
138     ierr = 0
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 
142     endif
143     call ext_ncd_open_for_read( trim(flnm2), 0, 0, "", dh2, Status)
144     if ( Status /= 0 ) go to 923
145     goto 924
146 923    continue
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
152     name = flnm2
153     Justplot = .true.
154 924    continue
155   if ( nargs .eq. 3 ) then
156     call get_command_argument(number=3, value=arg3)
157     read(arg3,*)levlim
158     print*,'LEVLIM = ',LEVLIM
159   endif
160   else
161      print*,'Usage: command file1 file2'
162      stop
163   endif
165 print*,'Just plot ',Justplot
167 if ( Justplot ) then
168   print*, 'flnm = ', trim(flnm)
169   first = .TRUE.
171   call ext_ncd_get_next_time(dh1, DateStr, Status_next_time)
173   ntimes = 0
174   DO WHILE ( Status_next_time .eq. 0 )
175     write(*,*)'Next Time ',TRIM(Datestr)
176     ntimes = ntimes + 1
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),'|'
181       start_index = 1
182       end_index = 1
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) 
186         cycle 
187       endif 
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)
192       ENDIF
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
200           ord = 'XYZ'
201         else if ( ndim .eq. 2 ) then
202           ord = 'XY'
203         else if ( ndim .eq. 1 ) then
204           ord = 'Z'
205         else if ( ndim .eq. 0 ) then
206           ord = '0'
207         endif
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
214                             ierr)
216         if ( ierr/=0 ) then
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)
222         endif
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)))
228            xlat = data(:,:,1,1)
229         ENDIF
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)
234         ENDIF
237         if ( VarName .eq. name ) then
238 #if 0
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)
246               enddo
247             enddo
248           endif
249         enddo
250 #else
251 ! give k-slices 
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)
258               enddo
259             enddo
260           endif
261         enddo
262 #endif
263         endif
265         deallocate(data)
266       endif
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
269         first = .FALSE.
270         icenter = 1 
271         jcenter = 1
272         ntries = 0
273         prev_icenter = 0 
274         prev_jcenter = 0  
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
280               jcenter = j
281 !write(*,*)'xlat ',ntries,icenter,jcenter,xlat(icenter,j),searchlat
282               exit
283             ENDIF
284           ENDDO
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
287               icenter = i
288 !write(*,*)'xlon ',ntries,icenter,jcenter,xlong(i,jcenter),searchlong
289               exit
290             ENDIF
291           ENDDO
292           ntries = ntries + 1
293         ENDDO
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'
299       ENDIF
300     enddo
301     call ext_ncd_get_next_time(dh1, DateStr, Status_next_time)
302   enddo
303 else
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
313     Status_next_time = 1
314   ENDIF
316   DO WHILE ( Status_next_time .eq. 0 .AND. Status_next_time2 .eq. 0 )
317     write(*,*)'Next Time ',TRIM(Datestr)
318     print 76
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),'|'
323       start_index = 1
324       end_index = 1
325       start_index2 = 1
326       end_index2 = 1
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
332         GOTO 1234
333       ENDIF
334       IF ( ndim /= ndim2 ) THEN
335         write(*,*)'Big difference: Number of dimensions for ',Varname,' differs in ',flnm2,'(',ndim,') /= (',ndim2
336         GOTO 1234
337       ENDIF
338       IF ( WrfType /= WrfType2 ) THEN
339         write(*,*)'Big difference: The types do not match'
340         GOTO 1234
341       ENDIF
342       if( WrfType == WRF_REAL .OR. WrfType == WRF_INTEGER ) then
343         DO i = 1, ndim
344           IF ( end_index(i) /= end_index2(i) ) THEN
345             write(*,*)'Big difference: dim ',i,' lengths differ for ',Varname,' differ in ',flnm2
346             GOTO 1234
347           ENDIF
348         ENDDO
349         DO i = ndim+1,3
350           start_index(i) = 1
351           end_index(i) = 1
352           start_index2(i) = 1
353           end_index2(i) = 1
354         ENDDO
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))
366         endif
368         if ( ndim .eq. 3 ) then
369           ord = 'XYZ'
370         else if ( ndim .eq. 2 ) then
371           ord = 'XY'
372         else if ( ndim .eq. 1 ) then
373           ord = 'Z'
374         else if ( ndim .eq. 0 ) then
375           ord = '0'
376         endif
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
384                                ierr)
385    
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)
392            ENDIF
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
398                                ierr)
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)
405            ENDIF
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
412                                ierr)
413    
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)
420            ENDIF
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
426                                ierr)
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)
433            ENDIF
434         endif
436         IFDIFFS=0
437         sumE = 0.0
438         sum1 = 0.0
439         sum2 = 0.0
440         diff1 = 0.0
441         diff2 = 0.0
442         n = 0 
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
445           cross = 0 
446           IKDIFFS = 0
448           if      (WrfType == WRF_REAL ) then
449              do i = 1, end_index(1)-cross
450                do j = 1, end_index(2)-cross
451                  a = data(I,J,K,1)
452                  b = data2(I,J,K,1)
453                  ! borrowed from  Thomas Oppe's comp program
454                  sumE = sumE + ( a - b ) * ( a - b )
455                  sum1 = sum1 + a * a
456                  sum2 = sum2 + b * b
457                  diff1 = max ( diff1 , abs ( a - b ) )
458                  diff2 = max ( diff2 , abs ( b ) )
459                  n = n + 1
460                  IF (a .ne. b) then
461                    IKDIFFS = IKDIFFS + 1
462                    IFDIFFS = IFDIFFS + 1
463                  ENDIF
464                ENDDO
465              ENDDO
466           else if (WrfType == WRF_INTEGER ) then
467              do i = 1, end_index(1)-cross
468                do j = 1, end_index(2)-cross
469                  a = idata(I,J,K,1)
470                  b = idata2(I,J,K,1)
471                  ! borrowed from  Thomas Oppe's comp program
472                  sumE = sumE + ( a - b ) * ( a - b )
473                  sum1 = sum1 + a * a
474                  sum2 = sum2 + b * b
475                  diff1 = max ( diff1 , abs ( a - b ) )
476                  diff2 = max ( diff2 , abs ( b ) )
477                  n = n + 1
478                  IF (a .ne. b) then
479                    IKDIFFS = IKDIFFS + 1
480                    IFDIFFS = IFDIFFS + 1
481                  ENDIF
482                ENDDO
483              ENDDO
484           endif
486          ENDIF
487         enddo
488         rmsE = sqrt ( sumE / dble( n ) )
489         rms1 = sqrt ( sum1 / dble( n ) )
490         rms2 = sqrt ( sum2 / dble( n ) )
491         serr = 0.0
492         IF ( sum2 .GT. 0.0d0 ) THEN
493           serr = sqrt ( sumE / sum2 )
494         ELSE
495           IF ( sumE .GT. 0.0d0 ) serr = 1.0
496         ENDIF
497         perr = 0.0
498         IF ( diff2 .GT. 0.0d0 ) THEN
499           perr = diff1/diff2
500         ELSE
501           IF ( diff1 .GT. 0.0d0 ) perr = 1.0
502         ENDIF
504         IF ( rms1 - rms2 .EQ. 0.0d0 ) THEN
505           digits = 15
506         ELSE
507           IF ( rms2 .NE. 0 ) THEN
508             tmp1 = 1.0d0/( ( abs( rms1 - rms2 ) ) / rms2 )
509             IF ( tmp1 .NE. 0 ) THEN
510               digits = log10(tmp1)
511             ENDIF
512           ENDIF
513         ENDIF
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 )
523         ENDIF
524         if      (WrfType == WRF_REAL    ) then
525           deallocate( data )
526           deallocate( data2)
527         else if (WrfType == WRF_INTEGER ) then
528           deallocate(idata )
529           deallocate(idata2)
530         endif
532       endif
533  1234 CONTINUE
534       call ext_ncd_get_next_var (dh1, VarName, Status_next_var)
535     enddo
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
543       Status_next_time = 1
544     ENDIF
545   enddo
547 endif
549 end program readv3
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 )
556   implicit none
557   integer n
558   integer a(n), b(n)
559   integer i
560   iveceq = .true.
561   do i = 1,n
562     if ( a(i) .ne. b(i) ) iveceq = .false.
563   enddo
564   return
565 end function iveceq
567 ! stubs for routines called by module_wrf_error (used by netcdf implementation of IO api)
568 SUBROUTINE wrf_abort
569   STOP
570 END SUBROUTINE wrf_abort
572 SUBROUTINE get_current_time_string( time_str )
573   CHARACTER(LEN=*), INTENT(OUT) :: time_str
574   time_str = ''
575 END SUBROUTINE get_current_time_string
577 SUBROUTINE get_current_grid_name( grid_str )
578   CHARACTER(LEN=*), INTENT(OUT) :: grid_str
579   grid_str = ''
580 END SUBROUTINE get_current_grid_name