var/mri4dvar/da_bilin.f90

   1 program da_bilin
   2
   3 !----------------------------------------------------------------------
   4 ! Purpose: Regridding increment from low-resolution to high-resolution
   5 !                     by using bilinear interpolation
   6 !
   7 ! Input     : fg                   -- low resolution first guess file
   8 !             wrfvar_output        -- low resolution analysis file
   9 !             wrfinput_hires       -- high resolution first guess file
  10 !
  11 ! Output    : wrfvar_output_hires  -- regridded high resolution analysis
  12 !
  13 ! Increment           = an_lores  - fg_lores
  14 ! wrfvar_output_hires = Increment + wrfinput_hires
  15 !
  16 ! In order to keep the domain size, it needs to match ( n - 1 )*ns + 1
  17 !
  18 ! where n  is the grid number in x or y
  19 !       ns is the refinement ratio between two resulotions
  20 !
  21 ! Compile:
  22 !
  23 !   pgf90 -o da_bilin.exe -I$NETCDF/include -L$NETCDF/lib -lnetcdf da_bilin.f90
  24 !
  25 ! Usage:
  26 !
  27 !   da_bilin.exe  [-h] [-fg_lores filename] [-an_lores filename]
  28 !                 [-fg_hires filename] [-ns n ] [-o outputfile]
  29 !
  30 !     -fg_lores    Optional, low resulotion first guess file,                 default - fg"
  31 !     -an_lores    Optional, low resulotion analysis file comes from wrfvar,  default - wrfvar_output"
  32 !     -fg_hires    Optional, high resultion first guess file,                 default - wrfinput_hires"
  33 !     -ns          Optional, the refinement ratio between two resulotions,    default - 3"
  34 !     -o           Optional, output high resulotion analysis file,            default - wrfvar_output_hires"
  35 !     -h           Show this help"
  36 !
  37 ! jliu@ucar.edu , 2011-12-15
  38 !----------------------------------------------------------------------
  39
  40   use netcdf
  41
  42   implicit none
  43
  44   !These variables' incremental will be regridded by default
  45   character (len=6), dimension(1:19) :: vNam
  46
  47   integer :: i, j, k, n, status
  48   integer :: nLat, nLon, oLat, oLon
  49   integer :: sLat, eLat, sLon, eLon
  50   integer :: rLat, rLon
  51
  52   integer :: ncidfg, ncidan, ncidout
  53   integer :: varid, nDims, dLen, varid_fg, varid_an, dimid
  54   integer :: regridsize, domainsize_out
  55
  56   real, dimension(:,:,:,:), allocatable :: fg, an, increment, var_out
  57   real, dimension(:,:),     allocatable :: iVar, oVar
  58
  59   integer, dimension(nf90_max_var_dims) :: vDimIDs
  60   integer, dimension(4)                 :: vdimsizes
  61
  62   character (len = 19), dimension(:),        allocatable :: times
  63   character (len = 255)  :: appname   = ""
  64   character (len = 255)  :: arg       = ""
  65   character (len = 255)  :: fg_lores  = "fg"
  66   character (len = 255)  :: an_lores  = "wrfvar_output"
  67   character (len = 255)  :: fg_hires  = "wrfinput_hires"
  68   character (len = 255)  :: f_out     = "wrfvar_output_hires"
  69   character (len = 255)  :: errmsg    = ""
  70   character (len = 8)    :: i_char    = ""
  71
  72   integer                :: ns        = 3
  73   !integer                :: cloud_cv_options = 0
  74   !integer                :: cv_w      = 0
  75
  76   LOGICAL :: file_exists
  77
  78   integer iargc
  79
  80   !These variables' incremental will be regridded by default
  81
  82   vNam(1)="U"
  83   vNam(2)="V"
  84   vNam(3)="T"
  85   vNam(4)="QVAPOR"
  86   vNam(5)="PH"
  87   vNam(6)="P"
  88   vNam(7)="MU"
  89   vNam(8)="U10"
  90   vNam(9)="V10"
  91   vNam(10)="T2"
  92   vNam(11)="Q2"
  93   vNam(12)="PSFC"
  94   vNam(13)="TH2"
  95
  96   vNam(14)="QCLOUD"
  97   vNam(15)="QRAIN"
  98   vNam(16)="QICE"
  99   vNam(17)="QSNOW"
 100   vNam(18)="QGRAUP"
 101   vNam(19)="W"
 102
 103   call getarg(0, appname)
 104   n=index(appname, '/', BACK=.true.)
 105   appname = trim(appname(n+1:))
 106
 107   DO i = 1, iargc(), 2
 108     arg=""
 109     call getarg(i, arg)
 110     select case ( trim(arg) )
 111       case ("-fg_lores")
 112         call getarg(i+1, arg)
 113         fg_lores=trim(arg)
 114       case ("-an_lores")
 115         call getarg(i+1, arg)
 116         an_lores=trim(arg)
 117       case ("-fg_hires")
 118         call getarg(i+1, arg)
 119         fg_hires=trim(arg)
 120       case ("-ns")
 121         call getarg(i+1, arg)
 122         read(arg, '(i3)') ns
 123       case ("-o")
 124         call getarg(i+1, arg)
 125         f_out=trim(arg)
 126       !case ("-cloud_cv_options")
 127       !  call getarg(i+1, arg)
 128       !  read(arg, '(i3)') cloud_cv_options
 129       !case ("-cv_w")
 130       !  call getarg(i+1, arg)
 131       !  read(arg, '(i3)') cv_w
 132       case default
 133         call show_usage()
 134         call exit(0)
 135     end select
 136   END DO
 137
 138   write (i_char, '(i8)') ns
 139
 140   inquire(FILE=trim(fg_hires), EXIST=file_exists)
 141
 142   if ( .not. file_exists ) then
 143     Write(*,*) "\nError: "//trim(fg_hires)//" not exists\n"
 144     call show_usage()
 145     call exit(-1)
 146   endif
 147
 148   call system("cp "//fg_hires//" "//f_out)
 149
 150   status = nf90_open(fg_lores, NF90_NOWRITE, ncidfg)
 151   errmsg = trim(fg_lores)
 152   if ( status /= nf90_noerr ) call nf90_handle_err(status, errmsg)
 153
 154   status = nf90_open(an_lores, NF90_NOWRITE, ncidan)
 155   errmsg = trim(an_lores)
 156   if ( status /= nf90_noerr ) call nf90_handle_err(status, errmsg)
 157
 158   status = nf90_open(f_out, NF90_WRITE, ncidout)
 159   errmsg= trim(f_out)
 160   if ( status /= nf90_noerr ) call nf90_handle_err(status, errmsg)
 161
 162   status = nf90_inq_dimid(ncidout, "west_east_stag", dimid)
 163   status = nf90_inquire_dimension(ncidout, dimid, len=dLen)
 164   domainsize_out = dLen
 165
 166   status = nf90_inq_dimid(ncidout, "south_north_stag", dimid)
 167   status = nf90_inquire_dimension(ncidout, dimid, len=dLen)
 168   domainsize_out = domainsize_out * dLen
 169
 170   status = nf90_inq_dimid(ncidfg, "west_east_stag", dimid)
 171   status = nf90_inquire_dimension(ncidfg, dimid, len=dLen)
 172   regridsize = (dLen-1)*ns+1
 173
 174   status = nf90_inq_dimid(ncidfg, "south_north_stag", dimid)
 175   status = nf90_inquire_dimension(ncidfg, dimid, len=dLen)
 176   regridsize = regridsize * ( (dLen-1)*ns+1 )
 177
 178   if ( regridsize /= domainsize_out ) then
 179     Write(*,'(a,i2,a)') "Error : It needs to match m = (n-1)*",ns, &
 180                         "+1 where n is coarse grid number in x or y, "// &
 181                         "m is fine grid number in x or y."
 182     call exit(-1)
 183   end if
 184
 185   write (i_char, '(i8)') ns
 186
 187   Write(*,*) " Input :"
 188   Write(*,*) "   Low resolution first guess  : "//trim(fg_lores)
 189   Write(*,*) "   Low resolution analysis     : "//trim(an_lores)
 190   Write(*,*) "   High resolution first guess : "//trim(fg_hires)
 191   Write(*,*) "   ns                          : "//adjustl(i_char)
 192   Write(*,*) "Output :"
 193   Write(*,*) "   High resolution analysis    : "//trim(f_out)
 194
 195   errmsg = ""
 196
 197   n = ubound(vNam,1)
 198   do i=1,n
 199
 200     Write (*,*) "Regridding increment for "//trim(vNam(i))
 201
 202     status = nf90_inq_varid(ncidout, trim(vNam(i)), varid)
 203     status = nf90_inquire_variable(ncidout, varid, ndims=nDims,dimids=vDimIDs)
 204
 205     vdimsizes = 1
 206     do j=1, nDims
 207       status = nf90_inquire_dimension(ncidout, vDimIDs(j), len = dLen )
 208       vdimsizes(j) = dLen
 209     end do
 210
 211     allocate(var_out(vdimsizes(1), vdimsizes(2), vdimsizes(3), vdimsizes(4)), stat=status)
 212
 213     status = nf90_get_var(ncidout, varid, var_out)
 214
 215     status = nf90_inq_varid(ncidfg, trim(vNam(i)), varid_fg)
 216     status = nf90_inq_varid(ncidan, trim(vNam(i)), varid_an)
 217
 218     status = nf90_inquire_variable(ncidfg, varid_fg, ndims=nDims,dimids=vDimIDs)
 219
 220     vdimsizes = 1
 221     do j=1, nDims
 222       status = nf90_inquire_dimension(ncidfg, vDimIDs(j), len = dLen )
 223       vdimsizes(j) = dLen
 224     end do
 225
 226     allocate(fg(vdimsizes(1), vdimsizes(2), vdimsizes(3), vdimsizes(4)), stat=status)
 227     allocate(an(vdimsizes(1), vdimsizes(2), vdimsizes(3), vdimsizes(4)), stat=status)
 228     allocate(increment(vdimsizes(1), vdimsizes(2), vdimsizes(3), vdimsizes(4)), stat=status)
 229
 230     status = nf90_get_var(ncidfg, varid_fg, fg)
 231     status = nf90_get_var(ncidan, varid_an, an)
 232
 233     increment = an - fg
 234
 235     nLon = vdimsizes(1)
 236     nLat = vdimsizes(2)
 237
 238     if ( trim(vNam(i) ) == "U" ) then
 239       rLat = nLat * ns
 240       rLon = (nLon-1) * ns + 1
 241       nLat = nLat + 2
 242     else
 243       rLon = nLon * ns
 244       rLat = (nLat-1) * ns + 1
 245       nLon = nLon + 2
 246       if ( trim(vNam(i)) /= "V" ) then
 247         rLat = nLat * ns
 248         nLat =  nLat + 2
 249       endif
 250     endif
 251
 252     oLon = ( nLon - 1 ) * ns + 1
 253     oLat = ( nLat - 1 ) * ns + 1
 254
 255     elat = (oLat - rLat) / 2
 256     slat =  oLat - rLat - elat + 1
 257
 258     elon = (oLon - rLon) / 2
 259     slon =  oLon - rLon - elon + 1
 260
 261     allocate(iVar(nLon, nLat), stat=status)
 262     allocate(oVar(oLon, oLat), stat=status)
 263
 264     do j=1, vdimsizes(4)
 265       do k=1, vdimsizes(3)
 266
 267         iVar = 0
 268         oVar = 0
 269
 270         select case ( trim(vNam(i)) )
 271           case ("U")
 272             iVar(:,2:nlat-1) = increment(:,:,k,j)
 273             iVar(:,1)        = iVar(:,2)
 274             iVar(:,nlat)     = iVar(:,nlat-1)
 275           case ("V")
 276             iVar(2:nlon-1,:) = increment(:,:,k,j)
 277             iVar(1,:)        = iVar(2,:)
 278             iVar(nlon,:)     = iVar(nlon-1,:)
 279           case default
 280             iVar(2:nlon-1,2:nlat-1) = increment(:,:,k,j)
 281             iVar(1,:)               = iVar(2,:)
 282             iVar(nlon,:)            = iVar(nlon-1,:)
 283             iVar(:,1)               = iVar(:,2)
 284             iVar(:,nlat)            = iVar(:,nlat-1)
 285         end select
 286
 287         call bilin(iVar, nLon, nLat, ns, oVar, oLon, oLat)
 288
 289         select case ( trim(vNam(i)) )
 290           case ("U")
 291             var_out(:,:,k,j) = var_out(:,:,k,j) + oVar(:,slat:olat-elat)
 292           case ("V")
 293             var_out(:,:,k,j) = var_out(:,:,k,j) + oVar(slon:olon-elon,:)
 294           case default
 295             var_out(:,:,k,j) = var_out(:,:,k,j) + oVar(slon:olon-elon,slat:olat-elat)
 296         end select
 297
 298       end do
 299     end do
 300
 301     status = nf90_put_var(ncidout, varid, var_out)
 302
 303     deallocate(var_out, stat=status)
 304     deallocate(iVar, stat=status)
 305     deallocate(oVar, stat=status)
 306     deallocate(fg, stat=status)
 307     deallocate(an, stat=status)
 308     deallocate(increment, stat=status)
 309
 310   end do
 311
 312   status = nf90_close(ncidfg)
 313   status = nf90_close(ncidan)
 314   status = nf90_close(ncidout)
 315
 316   Write(*,*) "Regridding increment completed successfully"
 317
 318 contains
 319   subroutine show_usage()
 320      Write(*,*) 'Usage :'//trim(appname)// &
 321      '[-h] [-fg_lores filename] [-an_lores filename] [-fg_hires filename] [-ns n ] [-o outputfile]'
 322      Write(*,*) "  -fg_lores    Optional, low resulotion first guess file,                 default - fg"
 323      Write(*,*) "  -an_lores    Optional, low resulotion analysis file comes from wrfvar,  default - wrfvar_output"
 324      Write(*,*) "  -fg_hires    Optional, high resultion first guess file,                 default - wrfinput_hires"
 325      Write(*,*) "  -ns          Optional, the refinement ratio between two resulotions,    default - 3"
 326      Write(*,*) "  -o           Optional, output high resulotion analysis file,            default - wrfvar_output_hires"
 327      Write(*,*) "  -h           Show this help"
 328   end subroutine show_usage
 329
 330   subroutine nf90_handle_err(status, errmsg)
 331     integer,          intent(in) :: status
 332     character(len=*), intent(in) :: errmsg
 333
 334     if(status /= nf90_noerr) then
 335       print *, trim(nf90_strerror(status))//" : "//trim(errmsg)
 336       Stop
 337     end if
 338   end subroutine nf90_handle_err
 339
 340   subroutine bilin(old,xi,yi,ns,new,xo,yo)
 341
 342     implicit none
 343
 344     integer,                 intent(in) :: xi,yi,xo,yo
 345     real, dimension(xi,yi),  intent(in) :: old
 346     integer,                 intent(in) :: ns
 347     real, dimension(xo,yo),  intent(out):: new
 348
 349     real                                :: im(1:ns+1,2)
 350     integer:: i,j,jm1,im1,ix1,ix2,iy1,iy2
 351
 352     forall(i=1:ns+1) im(i,2) = real(i-1)/ns
 353       im(:,1) = 1 - im(:,2)
 354
 355     do j=2,yi
 356       jm1 = j - 1
 357       iy2 = jm1 * ns + 1
 358       iy1 = iy2 - ns
 359       do i=2,xi
 360         im1 = i - 1
 361         ix2 = im1 * ns + 1
 362         ix1 = ix2 - ns
 363         new(ix1:ix2,iy1:iy2) = matmul(im,matmul(old(im1:i,jm1:j),transpose(im)))
 364       end do
 365     end do
 366
 367   end subroutine bilin
 368
 369 end program da_bilin