var/mri4dvar/da_thin.f90

   1 program da_thin
   2 !----------------------------------------------------------------------
   3 ! Purpose: Thinning wrfinput by using decimation
   4 !
   5 ! Input  : wrfinput_hires -- High resolution wrfinput
   6 !
   7 ! Output : wrfinput_lores -- Thinned wrfinput
   8 !
   9 ! jliu@ucar.edu, 2011-12-15
  10 !----------------------------------------------------------------------
  11
  12   use netcdf
  13
  14   implicit none
  15
  16   integer :: i, n
  17
  18   integer :: ncidin, ncidout, varid, varid_out, status
  19   integer :: nDims, nVars, nGlobalAtts, numsAtts, nTimes
  20   integer :: dLen, attLen, xtype, dID, unlimDimID, TID
  21   integer :: divided_exactly, dimid
  22
  23   integer :: dsizes(4), start(4), stride(4)
  24
  25   integer, dimension(nf90_max_var_dims) :: vDimIDs
  26
  27   integer, dimension(:),    allocatable :: vdimsizes
  28
  29   real :: fVal
  30
  31   real,                 dimension(:,:,:,:),  allocatable :: fVar
  32   integer,              dimension(:,:,:,:),  allocatable :: iVar
  33   character (len = 19), dimension(:),        allocatable :: times
  34
  35   character (len = 14 )           :: coordinates
  36   character (len = NF90_MAX_NAME) :: vNam, dNam, attNam
  37
  38   integer             :: decimation_factor = 3
  39   integer             ::            offset = 2
  40   character (len=255) ::            filin  = "wrfinput_hires"
  41   character (len=255) ::            filout = "wrfinput_lores"
  42   character (len=255) ::               arg = ""
  43   character (len=255) ::           appname = ""
  44   character(len=8)    ::            i_char =""
  45
  46   integer iargc
  47
  48   call getarg(0, appname)
  49   n=index(appname, '/', BACK=.true.)
  50   appname = trim(appname(n+1:))
  51
  52   DO i = 1, iargc(), 2
  53     call getarg(i, arg)
  54     select case ( trim(arg) )
  55       case ("-i")
  56         call getarg(i+1, arg)
  57         filin=trim(arg)
  58       case ("-o")
  59         call getarg(i+1, arg)
  60         filout=trim(arg)
  61       case ("-thin")
  62         call getarg(i+1, arg)
  63         read(arg, '(i3)') decimation_factor
  64       case default
  65         Write(*,*) "Usage : "//trim(appname)//" [-i inputfile] [-o outputfile] [-thin decimation_factor] [-h]"
  66         Write(*,*) "  -i     Optional, input filename,    default - wrfinput_hires"
  67         Write(*,*) "  -o     Optional, output filename,   default - wrfinput_lores"
  68         Write(*,*) "  -thin  Optional, decimation factor, default - 3"
  69         Write(*,*) "  -h     Shwo this usage"
  70         call exit(0)
  71     end select
  72   END DO
  73
  74   if ( mod(decimation_factor,2) == 0 ) then
  75     Write(*,*) "\nError : decimation factor must be odd number\n"
  76     call exit(-1)
  77   endif
  78
  79   status = nf90_open(filin, NF90_NOWRITE, ncidin)
  80   if ( status /= nf90_noerr ) then
  81     Write (*,*) "File open error. Please link the input file to "//trim(filin)
  82     call exit(-1)
  83   endif
  84
  85   status = nf90_inq_dimid(ncidin, "west_east_stag", dimid)
  86   status = nf90_inquire_dimension(ncidin, dimid, len=dLen)
  87   divided_exactly = mod((dLen-1),decimation_factor)
  88
  89   status = nf90_inq_dimid(ncidin, "south_north_stag", dimid)
  90   status = nf90_inquire_dimension(ncidin, dimid, len=dLen)
  91   divided_exactly = divided_exactly + mod((dLen-1),decimation_factor)
  92
  93   if ( divided_exactly /= 0 ) then
  94     Write (*,fmt='(a,i2,a)') "Failed to thinning. Grids need to match : ( n - 1 ) mod ",decimation_factor," = 0"
  95     call exit(-1)
  96   endif
  97
  98   status = nf90_create(filout, NF90_CLOBBER, ncidout)
  99   if ( status /= nf90_noerr) call nf90_handle_err(status)
 100
 101   status = nf90_inquire(ncidin, nDims, nVars, nGlobalAtts, unlimDimID)
 102   if ( status /= nf90_noerr ) call nf90_handle_err(status)
 103
 104   write (i_char, '(i8)') decimation_factor
 105
 106   Write (*,*) "       Input file : "//trim(filin)
 107   Write (*,*) "      Output file : "//trim(filout)
 108   Write (*,*) "decimation factor : "//adjustl(i_char)
 109
 110   do i=1, nGlobalAtts
 111     status = nf90_inq_attname(ncidin, NF90_GLOBAL, i, attNam)
 112     select case (trim(attNam))
 113       case ( "WEST-EAST_GRID_DIMENSION", "SOUTH-NORTH_GRID_DIMENSION", &
 114              "WEST-EAST_PATCH_END_UNSTAG", "WEST-EAST_PATCH_END_STAG", &
 115              "SOUTH-NORTH_PATCH_END_UNSTAG", "SOUTH-NORTH_PATCH_END_STAG" )
 116         status = nf90_get_att(ncidin, NF90_GLOBAL, attNam, fVal)
 117         status = nf90_put_att(ncidout, NF90_GLOBAL, attNam, int(( fVal - 1 ) / decimation_factor + 1) )
 118       case ("DX","DY", "DT" )
 119         status = nf90_get_att(ncidin, NF90_GLOBAL, attNam, fVal)
 120         status = nf90_put_att(ncidout, NF90_GLOBAL, attNam, fVal * decimation_factor )
 121       case default
 122         status = nf90_copy_att(ncidin, NF90_GLOBAL, attNam, ncidout, NF90_GLOBAL)
 123       end select
 124   end do
 125
 126   allocate (vdimsizes(nDims), stat=status)
 127
 128   do i=1, nDims
 129
 130     status = nf90_inquire_dimension(ncidin, i, name=dNam, len = dLen)
 131
 132     vdimsizes(i) = dLen
 133     select case (trim(dNam))
 134       case ("south_north_stag", "west_east_stag")
 135         vdimsizes(i) = (dLen - 1 ) / decimation_factor + 1
 136       case ("west_east", "south_north")
 137         vdimsizes(i) = dLen / decimation_factor
 138       case ("Time")
 139         allocate(times(dLen), stat=status)
 140         vdimsizes(i) = NF90_UNLIMITED
 141         nTimes = dLen
 142         TID = i
 143     end select
 144
 145     status = nf90_def_dim(ncidout, dNam, vdimsizes(i), dID)
 146
 147   end do
 148
 149   vdimsizes(TID) = nTimes
 150
 151   do varid=1, nVars
 152     status = nf90_inquire_variable(ncidin,varid,name=vNam,xtype=xtype,ndims=nDims,dimids=vDimIDs,natts=numsAtts)
 153     status = nf90_def_var(ncidout, trim(vNam), xtype, vDimIDs(1:nDims), varid_out)
 154     if(status /= nf90_NoErr) call nf90_handle_err(status)
 155     do i=1, numsAtts
 156       status = nf90_inq_attname(ncidin, varid, i, attNam)
 157       status = nf90_copy_att(ncidin, varid, trim(attNam), ncidout, varid_out)
 158       if(status /= nf90_NoErr) call nf90_handle_err(status)
 159     end do
 160   end do
 161
 162   status = nf90_enddef(ncidout)
 163
 164   offset = (decimation_factor + 1) / 2
 165
 166   do varid=1, nVars
 167
 168     status = nf90_inquire_variable(ncidin,varid,name=vNam,xtype=xtype,ndims=nDims,dimids=vDimIDs)
 169
 170     dsizes = 1
 171     do i = 1 , nDims
 172       dsizes(i) = vdimsizes(vDimIDs(i))
 173     end do
 174
 175     status = nf90_inquire_attribute(ncidin,varid,"coordinates")
 176
 177     if ( status == nf90_noerr ) then
 178
 179       Write(*,*) "Thinning for "//trim(vNam)
 180
 181       coordinates=char(0)
 182       status = nf90_get_att(ncidin, varid, "coordinates" , coordinates)
 183       !print *, coordinates
 184
 185       stride=(/decimation_factor,decimation_factor,1,1/)
 186
 187       n = index(coordinates, char(0)) - 1
 188       if ( n < 0 ) n = len(coordinates)
 189
 190       select case (trim(coordinates(1:n)))
 191         case ("XLONG_V XLAT_V")
 192           start=(/offset,1,1,1/)
 193         case ("XLONG_U XLAT_U")
 194           start=(/1,offset,1,1/)
 195         case ("XLONG XLAT")
 196           start=(/offset,offset,1,1/)
 197         case ("XLONG XLAT XTI")
 198           start=(/offset,offset,1,1/)
 199         case default
 200           print *, "Unkown coordinates : "//coordinates
 201           call exit(-1)
 202       end select
 203
 204     else
 205
 206       stride = 1
 207       start  = 1
 208
 209       if ( trim(vNam) == 'XLONG' .or. trim(vNam) == 'XLAT' ) then
 210           stride = (/decimation_factor,decimation_factor,1,1/)
 211           start  = (/offset,offset,1,1/)
 212       endif
 213
 214     endif
 215
 216     select case (xtype)
 217       case (nf90_float)
 218         allocate(fVar( dsizes(1), dsizes(2), dsizes(3), dsizes(4) ), stat=status)
 219         status = nf90_get_var(ncidin, varid, fVar, start=start, stride=stride)
 220         if ( status == nf90_noerr ) then
 221           if ( vNam == "RDX" .or. vNam == "RDY" ) then
 222             status = nf90_put_var(ncidout, varid, fVar / decimation_factor)
 223           else
 224             status = nf90_put_var(ncidout, varid, fVar)
 225           endif
 226           if ( status /= nf90_noerr) call nf90_handle_err(status)
 227         else
 228           call nf90_handle_err(status)
 229         endif
 230         deallocate(fVar, stat=status)
 231       case (nf90_int)
 232         allocate(iVar( dsizes(1), dsizes(2), dsizes(3), dsizes(4) ), stat=status)
 233         status = nf90_get_var(ncidin, varid, iVar, start=start, stride=stride)
 234         if ( status == nf90_noerr ) then
 235           status = nf90_put_var(ncidout, varid, iVar)
 236           if ( status /= nf90_noerr) call nf90_handle_err(status)
 237         else
 238           call nf90_handle_err(status)
 239         endif
 240         deallocate(iVar, stat=status)
 241       case (nf90_char)
 242         if ( trim(vNam) == "Times") then
 243           status = nf90_get_var(ncidin, varid, times)
 244           if ( status == nf90_noerr ) then
 245             status = nf90_put_var(ncidout, varid, times)
 246             if ( status /= nf90_noerr) call nf90_handle_err(status)
 247           else
 248             call nf90_handle_err(status)
 249           endif
 250           deallocate(times, stat=status)
 251         else
 252           print *, "Unkown character variable :"//trim(vNam)
 253           call exit(-1)
 254         endif
 255       case default
 256         print *, "Unkown xtype : ", xtype
 257         call exit(-1)
 258     end select
 259   end do
 260
 261   status = nf90_close(ncidin)
 262   status = nf90_close(ncidout)
 263
 264   Write(*,*) "Completed thinning successfully"
 265
 266 contains
 267
 268   subroutine nf90_handle_err(status)
 269     integer, intent ( in) :: status
 270
 271     if(status /= nf90_noerr) then
 272       print *, trim(nf90_strerror(status))
 273       call exit(-1)
 274     end if
 275   end subroutine nf90_handle_err
 276
 277 end program da_thin