external/io_netcdfpar/diffwrf.F90

   1 module read_util_module
   2
   3 contains
   4
   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
  47
  48
  49
  50  program readv3
  51   use wrf_data_ncpar
  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
  64
  65   integer :: flag, flag2
  66   integer :: iunit, iunit2
  67
  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
  86
  87   character (len=80), dimension(3)  ::  dimnames
  88   character (len=80) :: SysDepInfo
  89
  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
  95
  96   real, allocatable, dimension(:,:,:,:) :: data,data2
  97   real, allocatable, dimension(:,:)     :: xlat,xlong
  98
  99   integer :: ierr, ierr2, ier, ier2, Status, Status_next_time, Status_next_time2, Status_next_var, Status_next_var_2
 100   integer :: nargs
 101
 102   logical :: newtime = .TRUE.
 103   logical :: justplot, efound
 104
 105   logical, external :: iveceq
 106
 107   levlim = -1
 108
 109   call ext_ncdpar_ioinit(SysDepInfo,Status)
 110   call set_wrf_debug_level ( 1 )
 111
 112   nargs = command_argument_count()
 113
 114   Justplot = .false.
 115   searchcoords = .false.
 116 ! get arguments
 117   if ( nargs .ge. 2 ) then
 118     call get_command_argument(number=1, value=flnm)
 119     call get_command_argument(number=2, value=flnm2)
 120     IF ( flnm2(1:4) .EQ. '-lat' ) THEN
 121 print*,'reading ',TRIM(flnm2(5:))
 122       read(flnm2(5:),*)searchlat
 123       call get_command_argument(number=3, value=flnm2)
 124       IF ( flnm2(1:5) .EQ. '-long' ) THEN
 125 print*,'reading ',TRIM(flnm2(6:))
 126         read(flnm2(6:),*)searchlong
 127       ELSE
 128         write(*,*)'missing -long argument (no spaces after -lat or -long, either)'
 129         STOP
 130       ENDIF
 131       nargs = 0
 132       Justplot = .true.
 133       searchcoords = .true.
 134       call ext_ncdpar_open_for_read( trim(flnm), 0, 0, "", dh1, Status)
 135       goto 924
 136     ENDIF
 137     ierr = 0
 138     call ext_ncdpar_open_for_read( trim(flnm), 0, 0, "", dh1, Status)
 139     if ( Status /= 0 ) then
 140       print*,'error opening ',flnm, ' Status = ', Status ; stop
 141     endif
 142     call ext_ncdpar_open_for_read( trim(flnm2), 0, 0, "", dh2, Status)
 143     if ( Status /= 0 ) go to 923
 144     goto 924
 145 923    continue
 146
 147 ! bounce here if second name is not openable -- this would mean that
 148 ! it is a field name instead.
 149
 150     print*,'could not open ',flnm2
 151     name = flnm2
 152     Justplot = .true.
 153 924    continue
 154   if ( nargs .eq. 3 ) then
 155     call get_command_argument(number=3, value=arg3)
 156     read(arg3,*)levlim
 157     print*,'LEVLIM = ',LEVLIM
 158   endif
 159   else
 160      print*,'Usage: command file1 file2'
 161      stop
 162   endif
 163
 164 print*,'Just plot ',Justplot
 165
 166 if ( Justplot ) then
 167   print*, 'flnm = ', trim(flnm)
 168   first = .TRUE.
 169
 170   call ext_ncdpar_get_next_time(dh1, DateStr, Status_next_time)
 171
 172   ntimes = 0
 173   DO WHILE ( Status_next_time .eq. 0 )
 174     write(*,*)'Next Time ',TRIM(Datestr)
 175     ntimes = ntimes + 1
 176     call ext_ncdpar_get_next_var (dh1, VarName, Status_next_var)
 177     DO WHILE ( Status_next_var .eq. 0 )
 178 !    write(*,*)'Next Var |',TRIM(VarName),'|'
 179
 180       start_index = 1
 181       end_index = 1
 182       call ext_ncdpar_get_var_info (dh1,VarName,ndim,ordering,staggering,start_index,end_index, WrfType, ierr )
 183       if(WrfType /= WRF_REAL .AND. WrfType /= WRF_DOUBLE) then
 184         call ext_ncdpar_get_next_var (dh1, VarName, Status_next_var)
 185         cycle
 186       endif
 187       IF ( .NOT. searchcoords ) THEN
 188         write(*,'(A9,1x,I1,3(1x,I5),1x,A,1x,A)')&
 189                  VarName, ndim, end_index(1), end_index(2), end_index(3), &
 190                  trim(ordering), trim(DateStr)
 191       ENDIF
 192
 193       if ( VarName .eq. name .OR. TRIM(VarName) .EQ. 'XLAT' .OR. TRIM(VarName) .EQ. 'XLONG' ) then
 194         write(*,*)'Writing fort.88 file for ', trim(name)
 195
 196         allocate(data(end_index(1), end_index(2), end_index(3), 1))
 197
 198         if ( ndim .eq. 3 ) then
 199           ord = 'XYZ'
 200         else if ( ndim .eq. 2 ) then
 201           ord = 'XY'
 202         else if ( ndim .eq. 1 ) then
 203           ord = 'Z'
 204         else if ( ndim .eq. 0 ) then
 205           ord = '0'
 206         endif
 207
 208         call ext_ncdpar_read_field(dh1,DateStr,TRIM(VarName),data,WRF_REAL,0,0,0,ord, &
 209                             staggering, dimnames ,                      &
 210                             start_index,end_index,                      & !dom
 211                             start_index,end_index,                      & !mem
 212                             start_index,end_index,                      & !pat
 213                             ierr)
 214
 215         if ( ierr/=0 ) then
 216              write(*,*)'error reading data record'
 217              write(*,*)'  ndim = ', ndim
 218              write(*,*)'  end_index(1) ',end_index(1)
 219              write(*,*)'  end_index(2) ',end_index(2)
 220              write(*,*)'  end_index(3) ',end_index(3)
 221         endif
 222
 223 write(*,*)'name: ',TRIM(VarName)
 224         IF ( TRIM(VarName) .EQ. 'XLAT' .AND. .NOT. ALLOCATED(xlat)) THEN
 225 write(*,*)'allocating xlat'
 226            ALLOCATE(xlat(end_index(1), end_index(2)))
 227            xlat = data(:,:,1,1)
 228         ENDIF
 229         IF ( TRIM(VarName) .EQ. 'XLONG' .AND. .NOT. ALLOCATED(xlong)) THEN
 230 write(*,*)'allocating xlong'
 231            ALLOCATE(xlong(end_index(1), end_index(2)))
 232            xlong = data(:,:,1,1)
 233         ENDIF
 234
 235
 236         if ( VarName .eq. name ) then
 237 #if 0
 238 ! uncomment this to have the code give i-slices
 239         do i = 1, end_index(1)
 240           if ( levlim .eq. -1 .or. i .eq. levlim ) then
 241             write(88,*)end_index(2),end_index(3),' ',trim(name),' ',k,' time ',TRIM(Datestr)
 242             do k = start_index(3), end_index(3)
 243             do j = 1, end_index(2)
 244                 write(88,*) data(i,j,k,1)
 245               enddo
 246             enddo
 247           endif
 248         enddo
 249 #else
 250 ! give k-slices
 251         do k = start_index(3), end_index(3)
 252           if ( levlim .eq. -1 .or. k .eq. levlim ) then
 253             write(88,*)end_index(1),end_index(2),' ',trim(name),' ',k,' time ',TRIM(Datestr)
 254             do j = 1, end_index(2)
 255               do i = 1, end_index(1)
 256                 write(88,*) data(i,j,k,1)
 257               enddo
 258             enddo
 259           endif
 260         enddo
 261 #endif
 262         endif
 263
 264         deallocate(data)
 265       endif
 266       call ext_ncdpar_get_next_var (dh1, VarName, Status_next_var)
 267       IF ( ntimes .EQ. 1 .AND. ALLOCATED(xlong) .AND. ALLOCATED(xlat) .AND. first ) THEN
 268         first = .FALSE.
 269         icenter = 1
 270         jcenter = 1
 271         ntries = 0
 272         prev_icenter = 0
 273         prev_jcenter = 0
 274         DO WHILE ( ntries .LT. 10 .AND. (icenter .NE. prev_icenter .OR. jcenter .NE. prev_jcenter ))
 275           prev_icenter = icenter
 276           prev_jcenter = jcenter
 277           DO j = start_index(2), end_index(2)-1
 278             IF ( xlat(icenter,j) .LE. searchlat .AND. searchlat .LT. xlat(icenter,j+1) ) THEN
 279               jcenter = j
 280 !write(*,*)'xlat ',ntries,icenter,jcenter,xlat(icenter,j),searchlat
 281               exit
 282             ENDIF
 283           ENDDO
 284           DO i = start_index(1), end_index(1)-1
 285             IF ( xlong(i,jcenter) .LE. searchlong .AND. searchlong .LT. xlong(i+1,jcenter)) THEN
 286               icenter = i
 287 !write(*,*)'xlon ',ntries,icenter,jcenter,xlong(i,jcenter),searchlong
 288               exit
 289             ENDIF
 290           ENDDO
 291           ntries = ntries + 1
 292         ENDDO
 293         write(*,*)'Lon ',searchlong,' Lat ',searchlat,' : ',icenter,jcenter
 294         write(*,*)'Coordinates at that point ',xlong(icenter,jcenter),xlat(icenter,jcenter)
 295         write(*,*)'Coordinates at next point ',xlong(icenter+1,jcenter+1),xlat(icenter+1,jcenter+1)
 296         write(*,*)'Ntries : ',ntries
 297         if ( ntries .GE. 10 ) write(*,*)'max tries exceeded. Probably did not find'
 298       ENDIF
 299     enddo
 300     call ext_ncdpar_get_next_time(dh1, DateStr, Status_next_time)
 301   enddo
 302 else
 303   write (6,FMT='(4A)') 'Diffing ',trim(flnm),' ',trim(flnm2)
 304
 305   call ext_ncdpar_get_next_time(dh1, DateStr, Status_next_time)
 306   call ext_ncdpar_get_next_time(dh2, DateStr2, Status_next_time2)
 307
 308   IF ( DateStr .NE. DateStr2 ) THEN
 309     print*,'They differ big time.  Dates do not match'
 310     print*,'   ',flnm,' ',DateStr
 311     print*,'   ',flnm2,' ',DateStr2
 312     Status_next_time = 1
 313   ENDIF
 314
 315   DO WHILE ( Status_next_time .eq. 0 .AND. Status_next_time2 .eq. 0 )
 316     write(*,*)'Next Time ',TRIM(Datestr)
 317     print 76
 318     call ext_ncdpar_get_next_var (dh1, VarName, Status_next_var)
 319     DO WHILE ( Status_next_var .eq. 0 )
 320 !    write(*,*)'Next Var |',TRIM(VarName),'|'
 321
 322       start_index = 1
 323       end_index = 1
 324       start_index2 = 1
 325       end_index2 = 1
 326
 327       call ext_ncdpar_get_var_info (dh1,VarName,ndim,ordering,staggering,start_index,end_index, WrfType, ierr )
 328       call ext_ncdpar_get_var_info (dh2,VarName,ndim2,ordering2,staggering2,start_index2,end_index2, WrfType2, ierr )
 329       IF ( ierr /= 0 ) THEN
 330         write(*,*)'Big difference: ',VarName,' not found in ',flnm2
 331         GOTO 1234
 332       ENDIF
 333       IF ( ndim /= ndim2 ) THEN
 334         write(*,*)'Big difference: Number of dimensions for ',Varname,' differs in ',flnm2,'(',ndim,') /= (',ndim2
 335         GOTO 1234
 336       ENDIF
 337       IF ( WrfType /= WrfType2 ) THEN
 338         write(*,*)'Big difference: The types do not match'
 339         GOTO 1234
 340       ENDIF
 341       if( WrfType == WRF_REAL) then
 342         DO i = 1, ndim
 343           IF ( end_index(i) /= end_index2(i) ) THEN
 344             write(*,*)'Big difference: dim ',i,' lengths differ for ',Varname,' differ in ',flnm2
 345             GOTO 1234
 346           ENDIF
 347         ENDDO
 348         DO i = ndim+1,3
 349           start_index(i) = 1
 350           end_index(i) = 1
 351           start_index2(i) = 1
 352           end_index2(i) = 1
 353         ENDDO
 354
 355 !        write(*,'(A9,1x,I1,3(1x,I3),1x,A,1x,A)')&
 356 !                 VarName, ndim, end_index(1), end_index(2), end_index(3), &
 357 !                 trim(ordering), trim(DateStr)
 358
 359         allocate(data (end_index(1), end_index(2), end_index(3), 1))
 360         allocate(data2(end_index(1), end_index(2), end_index(3), 1))
 361
 362         if ( ndim .eq. 3 ) then
 363           ord = 'XYZ'
 364         else if ( ndim .eq. 2 ) then
 365           ord = 'XY'
 366         else if ( ndim .eq. 1 ) then
 367           ord = 'Z'
 368         else if ( ndim .eq. 0 ) then
 369           ord = '0'
 370         endif
 371
 372         call ext_ncdpar_read_field(dh1,DateStr,TRIM(VarName),data,WRF_REAL,0,0,0,ord,&
 373                             staggering, dimnames ,                      &
 374                             start_index,end_index,                      & !dom
 375                             start_index,end_index,                      & !mem
 376                             start_index,end_index,                      & !pat
 377                             ierr)
 378
 379         IF ( ierr /= 0 ) THEN
 380           write(*,*)'Error reading ',Varname,' from ',flnm
 381           write(*,*)'  ndim = ', ndim
 382           write(*,*)'  end_index(1) ',end_index(1)
 383           write(*,*)'  end_index(2) ',end_index(2)
 384           write(*,*)'  end_index(3) ',end_index(3)
 385         ENDIF
 386         call ext_ncdpar_read_field(dh2,DateStr,TRIM(VarName),data2,WRF_REAL,0,0,0,ord,&
 387                             staggering, dimnames ,                      &
 388                             start_index,end_index,                      & !dom
 389                             start_index,end_index,                      & !mem
 390                             start_index,end_index,                      & !pat
 391                             ierr)
 392         IF ( ierr /= 0 ) THEN
 393           write(*,*)'Error reading ',Varname,' from ',flnm2
 394           write(*,*)'  ndim = ', ndim
 395           write(*,*)'  end_index(1) ',end_index(1)
 396           write(*,*)'  end_index(2) ',end_index(2)
 397           write(*,*)'  end_index(3) ',end_index(3)
 398         ENDIF
 399
 400         IFDIFFS=0
 401         sumE = 0.0
 402         sum1 = 0.0
 403         sum2 = 0.0
 404         diff1 = 0.0
 405         diff2 = 0.0
 406         n = 0
 407         DO K = 1,end_index(3)-start_index(3)+1
 408          IF (LEVLIM.EQ.-1.OR.K.EQ.LEVLIM.OR.NDIM.eq.2) THEN
 409           cross = 0
 410           IKDIFFS = 0
 411           do i = 1, end_index(1)-cross
 412             do j = 1, end_index(2)-cross
 413               a = data(I,J,K,1)
 414               b = data2(I,J,K,1)
 415               ! borrowed from  Thomas Oppe's comp program
 416               sumE = sumE + ( a - b ) * ( a - b )
 417               sum1 = sum1 + a * a
 418               sum2 = sum2 + b * b
 419               diff1 = max ( diff1 , abs ( a - b ) )
 420               diff2 = max ( diff2 , abs ( b ) )
 421               n = n + 1
 422               IF (a .ne. b) then
 423                 IKDIFFS = IKDIFFS + 1
 424                 IFDIFFS = IFDIFFS + 1
 425               ENDIF
 426             ENDDO
 427           ENDDO
 428          ENDIF
 429         enddo
 430         rmsE = sqrt ( sumE / dble( n ) )
 431         rms1 = sqrt ( sum1 / dble( n ) )
 432         rms2 = sqrt ( sum2 / dble( n ) )
 433         serr = 0.0
 434         IF ( sum2 .GT. 0.0d0 ) THEN
 435           serr = sqrt ( sumE / sum2 )
 436         ELSE
 437           IF ( sumE .GT. 0.0d0 ) serr = 1.0
 438         ENDIF
 439         perr = 0.0
 440         IF ( diff2 .GT. 0.0d0 ) THEN
 441           perr = diff1/diff2
 442         ELSE
 443           IF ( diff1 .GT. 0.0d0 ) perr = 1.0
 444         ENDIF
 445
 446         IF ( rms1 - rms2 .EQ. 0.0d0 ) THEN
 447           digits = 15
 448         ELSE
 449           IF ( rms2 .NE. 0 ) THEN
 450             tmp1 = 1.0d0/( ( abs( rms1 - rms2 ) ) / rms2 )
 451             IF ( tmp1 .NE. 0 ) THEN
 452               digits = log10(tmp1)
 453             ENDIF
 454           ENDIF
 455         ENDIF
 456
 457         IF (IFDIFFS .NE. 0 ) THEN
 458            ! create the fort.88 and fort.98 files because regression scripts will
 459            ! look for these to see if there were differences.
 460            write(88,*)trim(varname)
 461            write(98,*)trim(varname)
 462            PRINT 77,trim(varname), IFDIFFS, ndim, rms1, rms2, digits, rmsE, perr
 463  76 FORMAT (5x,'Field ',2x,'Ndifs',4x,'Dims ',6x,'RMS (1)',12x,'RMS (2)',5x,'DIGITS',4x,'RMSE',5x,'pntwise max')
 464  77 FORMAT ( A10,1x,I9,2x,I3,1x,e18.10,1x,e18.10,1x,i3,1x,e12.4,1x,e12.4 )
 465         ENDIF
 466         deallocate(data)
 467         deallocate(data2)
 468
 469       endif
 470  1234 CONTINUE
 471       call ext_ncdpar_get_next_var (dh1, VarName, Status_next_var)
 472     enddo
 473     call ext_ncdpar_get_next_time(dh1, DateStr, Status_next_time)
 474     call ext_ncdpar_get_next_time(dh2, DateStr2, Status_next_time2)
 475     IF ( DateStr .NE. DateStr2 ) THEN
 476       print*,'They differ big time.  Dates do not match'
 477       print*,'They differ big time.  Dates do not match'
 478       print*,'   ',flnm,' ',DateStr
 479       print*,'   ',flnm2,' ',DateStr2
 480       Status_next_time = 1
 481     ENDIF
 482   enddo
 483
 484 endif
 485
 486 end program readv3
 487
 488 logical function wrf_dm_on_monitor()
 489   wrf_dm_on_monitor=.true.
 490 end function wrf_dm_on_monitor
 491
 492 logical function iveceq( a, b, n )
 493   implicit none
 494   integer n
 495   integer a(n), b(n)
 496   integer i
 497   iveceq = .true.
 498   do i = 1,n
 499     if ( a(i) .ne. b(i) ) iveceq = .false.
 500   enddo
 501   return
 502 end function iveceq
 503
 504 ! stubs for routines called by module_wrf_error (used by netcdf implementation of IO api)
 505 SUBROUTINE wrf_abort
 506   STOP
 507 END SUBROUTINE wrf_abort
 508
 509 SUBROUTINE get_current_time_string( time_str )
 510   CHARACTER(LEN=*), INTENT(OUT) :: time_str
 511   time_str = ''
 512 END SUBROUTINE get_current_time_string
 513
 514 SUBROUTINE get_current_grid_name( grid_str )
 515   CHARACTER(LEN=*), INTENT(OUT) :: grid_str
 516   grid_str = ''
 517 END SUBROUTINE get_current_grid_name
 518