wrfv2_fire/external/io_int/diffwrf.F

   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 if(n.eq.0) n=1
 374         rmsE = sqrt ( sumE / dble( n ) )
 375         rms1 = sqrt ( sum1 / dble( n ) )
 376         rms2 = sqrt ( sum2 / dble( n ) )
 377         serr = 0.0
 378         IF ( sum2 .GT. 0.0d0 ) THEN
 379           serr = sqrt ( sumE / sum2 )
 380         ELSE
 381           IF ( sumE .GT. 0.0d0 ) serr = 1.0
 382         ENDIF
 383         perr = 0.0
 384         IF ( diff2 .GT. 0.0d0 ) THEN
 385           perr = diff1/diff2
 386         ELSE
 387           IF ( diff1 .GT. 0.0d0 ) perr = 1.0
 388         ENDIF
 389
 390         digits = 0
 391         IF ( rms1 - rms2 .EQ. 0.0d0 ) THEN
 392           digits = 15
 393         ELSE
 394           IF ( rms2 .NE. 0 ) THEN
 395             tmp1 = 1.0d0/( ( abs( rms1 - rms2 ) ) / rms2 )
 396             IF ( tmp1 .NE. 0 ) THEN
 397               digits = log10(tmp1)
 398              ENDIF
 399           ENDIF
 400         ENDIF
 401
 402         IF (IFDIFFS .NE. 0 ) THEN
 403            ! create the fort.88 and fort.98 files because regression scripts will
 404            ! look for these to see if there were differences.
 405            write(88,*)trim(varname)
 406            write(98,*)trim(varname)
 407            PRINT 77,trim(varname), IFDIFFS, ndim, rms1, rms2, digits, rmsE, perr
 408  76 FORMAT (5x,'Field ',2x,'Ndifs',4x,'Dims ',6x,'RMS (1)',12x,'RMS (2)',5x,'DIGITS',4x,'RMSE',5x,'pntwise max')
 409  77 FORMAT ( A10,1x,I9,2x,I3,1x,e18.10,1x,e18.10,1x,i3,1x,e12.4,1x,e12.4 )
 410         ENDIF
 411         deallocate(data)
 412         deallocate(data2)
 413
 414       endif
 415 1234  CONTINUE
 416       call ext_int_get_next_var (dh1, VarName, Status_next_var)
 417       call ext_int_get_next_var (dh2, VarName, Status_next_var)
 418     enddo
 419     call ext_int_get_next_time(dh1, DateStr, Status_next_time)
 420     call ext_int_get_next_time(dh2, DateStr2, Status_next_time2)
 421     IF ( DateStr .NE. DateStr2 ) THEN
 422       print*,'They differ big time.  Dates do not match'
 423       print*,'They differ big time.  Dates do not match'
 424       print*,'   ',flnm,' ',DateStr
 425       print*,'   ',flnm2,' ',DateStr2
 426       Status_next_time = 1
 427     ENDIF
 428   enddo
 429
 430 endif
 431
 432 end program readv3
 433
 434 logical function iveceq( a, b, n )
 435   implicit none
 436   integer n
 437   integer a(n), b(n)
 438   integer i
 439   iveceq = .true.
 440   do i = 1,n
 441     if ( a(i) .ne. b(i) ) iveceq = .false.
 442   enddo
 443   return
 444 end function iveceq
 445
 446
 447 ! stubs for routines called by module_wrf_error (used by this implementation of IO api)
 448 SUBROUTINE wrf_abort
 449   STOP
 450 END SUBROUTINE wrf_abort
 451
 452 SUBROUTINE get_current_time_string( time_str )
 453   CHARACTER(LEN=*), INTENT(OUT) :: time_str
 454   time_str = ''
 455 END SUBROUTINE get_current_time_string
 456
 457 SUBROUTINE get_current_grid_name( grid_str )
 458   CHARACTER(LEN=*), INTENT(OUT) :: grid_str
 459   grid_str = ''
 460 END SUBROUTINE get_current_grid_name
 461