text: program da_thin
  !----------------------------------------------------------------------
  ! Purpose: Thinning wrfinput by using decimation
  !
  ! Input  : wrfinput_hires -- High resolution wrfinput
  !
  ! Output : wrfinput_lores -- Thinned wrfinput
  !
  ! jliu@ucar.edu, 2011-12-15
  !----------------------------------------------------------------------
  
    use netcdf
  
    implicit none
  
    integer :: i, n
  
    integer :: ncidin, ncidout, varid, varid_out, status
    integer :: nDims, nVars, nGlobalAtts, numsAtts, nTimes
    integer :: dLen, attLen, xtype, dID, unlimDimID, TID 
    integer :: divided_exactly, dimid
  
    integer :: dsizes(4), start(4), stride(4)
  
    integer, dimension(nf90_max_var_dims) :: vDimIDs
  
    integer, dimension(:),    allocatable :: vdimsizes
  
    real :: fVal
  
    real,                 dimension(:,:,:,:),  allocatable :: fVar
    integer,              dimension(:,:,:,:),  allocatable :: iVar
    character (len = 19), dimension(:),        allocatable :: times
  
    character (len = 14 )           :: coordinates
    character (len = NF90_MAX_NAME) :: vNam, dNam, attNam
  
    integer             :: decimation_factor = 3
    integer             ::            offset = 2
    character (len=255) ::            filin  = "wrfinput_hires"
    character (len=255) ::            filout = "wrfinput_lores"
    character (len=255) ::               arg = ""
    character (len=255) ::           appname = ""
    character(len=8)    ::            i_char =""
  
    integer iargc
  
    call getarg(0, appname)
    n=index(appname, '/', BACK=.true.)
    appname = trim(appname(n+1:))
      
    DO i = 1, iargc(), 2
      call getarg(i, arg)
      select case ( trim(arg) )
        case ("-i")
          call getarg(i+1, arg)
          filin=trim(arg)
        case ("-o")
          call getarg(i+1, arg)
          filout=trim(arg)
        case ("-thin")
          call getarg(i+1, arg)
          read(arg, '(i3)') decimation_factor
        case default 
          Write(*,*) "Usage : "//trim(appname)//" [-i inputfile] [-o outputfile] [-thin decimation_factor] [-h]"
          Write(*,*) "  -i     Optional, input filename,    default - wrfinput_hires"
          Write(*,*) "  -o     Optional, output filename,   default - wrfinput_lores"
          Write(*,*) "  -thin  Optional, decimation factor, default - 3"
          Write(*,*) "  -h     Shwo this usage"
          call exit(0)
      end select
    END DO
  
    if ( mod(decimation_factor,2) == 0 ) then
      Write(*,*) "\nError : decimation factor must be odd number\n"
      call exit(-1)
    endif
  
    status = nf90_open(filin, NF90_NOWRITE, ncidin)
    if ( status /= nf90_noerr ) then
      Write (*,*) "File open error. Please link the input file to "//trim(filin)
      call exit(-1)
    endif
  
    status = nf90_inq_dimid(ncidin, "west_east_stag", dimid)
    status = nf90_inquire_dimension(ncidin, dimid, len=dLen)
    divided_exactly = mod((dLen-1),decimation_factor)
  
    status = nf90_inq_dimid(ncidin, "south_north_stag", dimid)
    status = nf90_inquire_dimension(ncidin, dimid, len=dLen)
    divided_exactly = divided_exactly + mod((dLen-1),decimation_factor)
  
    if ( divided_exactly /= 0 ) then
      Write (*,fmt='(a,i2,a)') "Failed to thinning. Grids need to match : ( n - 1 ) mod ",decimation_factor," = 0"
      call exit(-1)
    endif
  
    status = nf90_create(filout, NF90_CLOBBER, ncidout)
    if ( status /= nf90_noerr) call nf90_handle_err(status)
  
    status = nf90_inquire(ncidin, nDims, nVars, nGlobalAtts, unlimDimID)
    if ( status /= nf90_noerr ) call nf90_handle_err(status)
  
    write (i_char, '(i8)') decimation_factor
  
    Write (*,*) "       Input file : "//trim(filin)
    Write (*,*) "      Output file : "//trim(filout)
    Write (*,*) "decimation factor : "//adjustl(i_char)
  
    do i=1, nGlobalAtts
      status = nf90_inq_attname(ncidin, NF90_GLOBAL, i, attNam)
      select case (trim(attNam))
        case ( "WEST-EAST_GRID_DIMENSION", "SOUTH-NORTH_GRID_DIMENSION", &
               "WEST-EAST_PATCH_END_UNSTAG", "WEST-EAST_PATCH_END_STAG", &
               "SOUTH-NORTH_PATCH_END_UNSTAG", "SOUTH-NORTH_PATCH_END_STAG" )
          status = nf90_get_att(ncidin, NF90_GLOBAL, attNam, fVal)
          status = nf90_put_att(ncidout, NF90_GLOBAL, attNam, int(( fVal - 1 ) / decimation_factor + 1) )
        case ("DX","DY", "DT" )
          status = nf90_get_att(ncidin, NF90_GLOBAL, attNam, fVal)
          status = nf90_put_att(ncidout, NF90_GLOBAL, attNam, fVal * decimation_factor )
        case default
          status = nf90_copy_att(ncidin, NF90_GLOBAL, attNam, ncidout, NF90_GLOBAL)
        end select
    end do
  
    allocate (vdimsizes(nDims), stat=status)
  
    do i=1, nDims
  
      status = nf90_inquire_dimension(ncidin, i, name=dNam, len = dLen)
  
      vdimsizes(i) = dLen
      select case (trim(dNam))
        case ("south_north_stag", "west_east_stag")
          vdimsizes(i) = (dLen - 1 ) / decimation_factor + 1
        case ("west_east", "south_north")
          vdimsizes(i) = dLen / decimation_factor
        case ("Time")
          allocate(times(dLen), stat=status)
          vdimsizes(i) = NF90_UNLIMITED
          nTimes = dLen
          TID = i
      end select
  
      status = nf90_def_dim(ncidout, dNam, vdimsizes(i), dID)
  
    end do
  
    vdimsizes(TID) = nTimes
  
    do varid=1, nVars
      status = nf90_inquire_variable(ncidin,varid,name=vNam,xtype=xtype,ndims=nDims,dimids=vDimIDs,natts=numsAtts)
      status = nf90_def_var(ncidout, trim(vNam), xtype, vDimIDs(1:nDims), varid_out)
      if(status /= nf90_NoErr) call nf90_handle_err(status)
      do i=1, numsAtts
        status = nf90_inq_attname(ncidin, varid, i, attNam)
        status = nf90_copy_att(ncidin, varid, trim(attNam), ncidout, varid_out)
        if(status /= nf90_NoErr) call nf90_handle_err(status)
      end do
    end do
  
    status = nf90_enddef(ncidout)
  
    offset = (decimation_factor + 1) / 2
  
    do varid=1, nVars
  
      status = nf90_inquire_variable(ncidin,varid,name=vNam,xtype=xtype,ndims=nDims,dimids=vDimIDs)
  
      dsizes = 1
      do i = 1 , nDims
        dsizes(i) = vdimsizes(vDimIDs(i))
      end do
  
      status = nf90_inquire_attribute(ncidin,varid,"coordinates")
  
      if ( status == nf90_noerr ) then
  
        Write(*,*) "Thinning for "//trim(vNam)
  
        coordinates=char(0)
        status = nf90_get_att(ncidin, varid, "coordinates" , coordinates)
        !print *, coordinates
  
        stride=(/decimation_factor,decimation_factor,1,1/)
  
        n = index(coordinates, char(0)) - 1
        if ( n < 0 ) n = len(coordinates)
  
        select case (trim(coordinates(1:n)))
          case ("XLONG_V XLAT_V")
            start=(/offset,1,1,1/)
          case ("XLONG_U XLAT_U")
            start=(/1,offset,1,1/)
          case ("XLONG XLAT")
            start=(/offset,offset,1,1/)
          case ("XLONG XLAT XTI")
            start=(/offset,offset,1,1/)
          case default
            print *, "Unkown coordinates : "//coordinates
            call exit(-1)
        end select
  
      else
  
        stride = 1
        start  = 1 
        
        if ( trim(vNam) == 'XLONG' .or. trim(vNam) == 'XLAT' ) then
            stride = (/decimation_factor,decimation_factor,1,1/)
            start  = (/offset,offset,1,1/)
        endif 
  
      endif
  
      select case (xtype)
        case (nf90_float)
          allocate(fVar( dsizes(1), dsizes(2), dsizes(3), dsizes(4) ), stat=status)
          status = nf90_get_var(ncidin, varid, fVar, start=start, stride=stride)
          if ( status == nf90_noerr ) then
            if ( vNam == "RDX" .or. vNam == "RDY" ) then
              status = nf90_put_var(ncidout, varid, fVar / decimation_factor)
            else
              status = nf90_put_var(ncidout, varid, fVar)
            endif
            if ( status /= nf90_noerr) call nf90_handle_err(status)
          else
            call nf90_handle_err(status)
          endif
          deallocate(fVar, stat=status)
        case (nf90_int)
          allocate(iVar( dsizes(1), dsizes(2), dsizes(3), dsizes(4) ), stat=status)
          status = nf90_get_var(ncidin, varid, iVar, start=start, stride=stride)
          if ( status == nf90_noerr ) then
            status = nf90_put_var(ncidout, varid, iVar)
            if ( status /= nf90_noerr) call nf90_handle_err(status)
          else
            call nf90_handle_err(status)
          endif
          deallocate(iVar, stat=status)
        case (nf90_char)
          if ( trim(vNam) == "Times") then
            status = nf90_get_var(ncidin, varid, times)
            if ( status == nf90_noerr ) then
              status = nf90_put_var(ncidout, varid, times)
              if ( status /= nf90_noerr) call nf90_handle_err(status)
            else
              call nf90_handle_err(status)
            endif
            deallocate(times, stat=status)
          else
            print *, "Unkown character variable :"//trim(vNam)
            call exit(-1)  
          endif
        case default
          print *, "Unkown xtype : ", xtype
          call exit(-1)
      end select
    end do
  
    status = nf90_close(ncidin)
    status = nf90_close(ncidout)
  
    Write(*,*) "Completed thinning successfully"
  
  contains
  
    subroutine nf90_handle_err(status)
      integer, intent ( in) :: status
  
      if(status /= nf90_noerr) then
        print *, trim(nf90_strerror(status))
        call exit(-1) 
      end if
    end subroutine nf90_handle_err
  
  end program da_thin
  