source: trunk/WRF.COMMON/WRFV2/external/io_int/diffwrf.F @ 3567

Last change on this file since 3567 was 11, checked in by aslmd, 14 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 13.8 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 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=120) :: flnm
71  character(len=120) :: 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
133923    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.
141924    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
152print*,'Just plot ',Justplot
153
154if ( 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
235else
236  print*,'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
4131234  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
428endif
429
430end program readv3
431
432logical 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
442end function iveceq
443
444
445! stubs for routines called by module_wrf_error (used by this implementation of IO api)
446SUBROUTINE wrf_abort
447  STOP
448END SUBROUTINE wrf_abort
449
450SUBROUTINE get_current_time_string( time_str )
451  CHARACTER(LEN=*), INTENT(OUT) :: time_str
452  time_str = ''
453END SUBROUTINE get_current_time_string
454
455SUBROUTINE get_current_grid_name( grid_str )
456  CHARACTER(LEN=*), INTENT(OUT) :: grid_str
457  grid_str = ''
458END SUBROUTINE get_current_grid_name
459
Note: See TracBrowser for help on using the repository browser.