source: trunk/WRF.COMMON/WRFV3/external/io_netcdf/diffwrf.F90

Last change on this file was 2759, checked in by aslmd, 2 years ago

adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 13.7 KB
Line 
1module read_util_module
2
3#ifdef crayx1
4#define iargc ipxfargc
5#endif
6
7contains
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
62end 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
138923    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.
146924    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
157print*,'Just plot ',Justplot
158
159if ( 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
244else
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
426endif
427
428end program readv3
429
430logical 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
440end function iveceq
441
442! stubs for routines called by module_wrf_error (used by netcdf implementation of IO api)
443SUBROUTINE wrf_abort
444  STOP
445END SUBROUTINE wrf_abort
446
447SUBROUTINE get_current_time_string( time_str )
448  CHARACTER(LEN=*), INTENT(OUT) :: time_str
449  time_str = ''
450END SUBROUTINE get_current_time_string
451
452SUBROUTINE get_current_grid_name( grid_str )
453  CHARACTER(LEN=*), INTENT(OUT) :: grid_str
454  grid_str = ''
455END SUBROUTINE get_current_grid_name
456
Note: See TracBrowser for help on using the repository browser.