util/src/int2nc.F

   1 program int2nc
   2
   3 !  use netcdf
   4    use module_debug
   5    use misc_definitions_module
   6    use read_met_module
   7
   8    implicit none
   9
  10   include "netcdf.inc"
  11
  12    integer, parameter :: NDIMS = 2
  13    integer :: istatus, dim, i, varid, ablevel, proj, nproj
  14    real :: fcst, slat, slon, dlat, dlon, nlat, dxn, dyn
  15    real :: xloncen, tlat1, tlat2, radius, si, sj
  16    logical :: windrot
  17    character (len=132) field, name, cablevel, date, source, units, desc, flnm, nfile
  18    real, allocatable, dimension(:,:) :: data
  19    integer :: ncid
  20    integer, dimension(20) :: dval
  21    integer :: dimids(NDIMS)
  22    integer :: tmp_dims(NDIMS)
  23
  24    type (met_data)             :: fg_data
  25
  26    character (len=*),dimension(10),parameter :: dname = (/"i1","j1","i2","j2","i3","j3","i4","j4","i5","j5"/)
  27    character (len=*),parameter :: DATEV = "date"
  28    character (len=*),parameter :: FCSTV = "forecast"
  29    character (len=*),parameter :: SOURCEV = "map_source"
  30    character (len=*),parameter :: LEVELV = "level"
  31    character (len=*),parameter :: FIELDV = "field"
  32    character (len=*),parameter :: UNITSV = "units"
  33    character (len=*),parameter :: DESCV = "description"
  34    character (len=*),parameter :: NX = "nx"
  35    character (len=*),parameter :: NY = "ny"
  36    character (len=*),parameter :: IPROJ = "projection"
  37    character (len=*),parameter :: STARTI = "starti"
  38    character (len=*),parameter :: STARTJ = "startj"
  39    character (len=*),parameter :: STARTLAT = "startlat"
  40    character (len=*),parameter :: STARTLON = "startlon"
  41    character (len=*),parameter :: DELTALAT = "deltalat"
  42    character (len=*),parameter :: DELTALON = "deltalon"
  43    character (len=*),parameter :: NLATS = "nlats"
  44    character (len=*),parameter :: DX = "dx"
  45    character (len=*),parameter :: DY = "dy"
  46    character (len=*),parameter :: XLONC = "xlonc"
  47    character (len=*),parameter :: TRUELAT1 = "truelat1"
  48    character (len=*),parameter :: TRUELAT2 = "truelat2"
  49    character (len=*),parameter :: EARTH_RADIUS = "earth_radius"
  50    character (len=*),parameter :: IS_WIND_GRID_REL = "is_wind_grid_rel"
  51    character (len=*),parameter :: FILLVALUE = "_FillValue"
  52
  53    dval = 0
  54
  55    !  Get the input file name from the command line.
  56    call getarg (1,flnm)
  57
  58    if (flnm(1:1) == ' ') then
  59       print *,'USAGE: int2nc.exe <filename>'
  60       print *,'       where <filename> is the name of an intermediate-format file'
  61       stop
  62    end if
  63    nfile = trim(adjustl(flnm))//".nc"
  64
  65    call set_debug_level(WARN)
  66
  67    call read_met_init(trim(flnm), .true., '0000-00-00_00', istatus)
  68    call check(nf_create(trim(nfile),nf_clobber,ncid))
  69
  70    print*, 'OPENING FILE: ',trim(adjustl(flnm))
  71
  72    if(istatus == 0) then
  73      call read_next_met_field(fg_data,istatus)
  74      do while (istatus == 0)
  75        tmp_dims(1) = fg_data%nx
  76        tmp_dims(2) = fg_data%ny
  77          do dim = 1,2
  78            i = 1
  79            CHECKDIMS : DO
  80              if(dval(i) == 0) then
  81                dval(i) = tmp_dims(dim)
  82                call check(nf_def_dim(ncid,dname(i),dval(i),i))
  83                EXIT CHECKDIMS
  84              else
  85                if (dval(i) == tmp_dims(dim) ) then
  86                  EXIT CHECKDIMS
  87                end if
  88              end if
  89              i = i+1
  90              CYCLE CHECKDIMS
  91            END DO CHECKDIMS
  92          end do
  93        call read_next_met_field(fg_data,istatus)
  94      end do
  95    else
  96      print *, 'File = ',trim(flnm)
  97      print *, 'Problem with input file, I can''t open it'
  98      stop
  99    end if
 100
 101
 102    call read_met_close()
 103    call check(nf_close(ncid))
 104
 105    call read_met_init(trim(flnm), .true., '0000-00-00_00', istatus)
 106    call check(nf_open(trim(nfile),nf_write,ncid))
 107
 108    if (istatus == 0) then
 109       call read_next_met_field(fg_data, istatus)
 110       do while (istatus == 0)
 111
 112          date = trim(adjustl(fg_data%hdate))
 113          fcst = fg_data%xfcst
 114          source = fg_data%map_source
 115          field = fg_data%field
 116          units = fg_data%units
 117          desc = fg_data%desc
 118          ablevel = fg_data%xlvl
 119          write(cablevel,'(I6)') ablevel
 120          name = trim(adjustl(field))//"__0"//trim(adjustl(cablevel))
 121          print *,"Reading Field, Level: ",trim(adjustl(field)),", ",trim(adjustl(cablevel))
 122
 123          nproj = fg_data%iproj
 124          proj = nproj
 125          if(nproj == 1) proj = 3
 126          if(nproj == 3) proj = 1
 127          if(proj == 0) then         ! Cylindrical Equidistand
 128            si = fg_data%starti
 129            sj = fg_data%startj
 130            slat = fg_data%startlat
 131            slon = fg_data%startlon
 132            dlat = fg_data%deltalat
 133            dlon = fg_data%deltalon
 134            radius = fg_data%earth_radius
 135          else if(proj == 1) then    ! Mercator
 136            si = fg_data%starti
 137            sj = fg_data%startj
 138            slat = fg_data%startlat
 139            slon = fg_data%startlon
 140            dxn = fg_data%dx
 141            dyn = fg_data%dy
 142            tlat1 = fg_data%truelat1
 143            radius = fg_data%earth_radius
 144          else if(proj == 3) then    ! Lambert Conformal
 145            si = fg_data%starti
 146            sj = fg_data%startj
 147            slat = fg_data%startlat
 148            slon = fg_data%startlon
 149            dxn = fg_data%dx
 150            dyn = fg_data%dy
 151            xloncen = fg_data%xlonc
 152            tlat1 = fg_data%truelat1
 153            tlat2 = fg_data%truelat2
 154            radius = fg_data%earth_radius
 155          else if(proj == 4) then    ! Gaussian
 156            si = fg_data%starti
 157            sj = fg_data%startj
 158            slat = fg_data%startlat
 159            slon = fg_data%startlon
 160            nlat = fg_data%deltalat
 161            dlon = fg_data%deltalon
 162            radius = fg_data%earth_radius
 163          else if(proj == 5) then    ! Polar Stereographic
 164            si = fg_data%starti
 165            sj = fg_data%startj
 166            slat = fg_data%startlat
 167            slon = fg_data%startlon
 168            dxn = fg_data%dx
 169            dyn = fg_data%dy
 170            xloncen = fg_data%xlonc
 171            tlat1 = fg_data%truelat1
 172            radius = fg_data%earth_radius
 173          end if
 174          windrot = fg_data%is_wind_grid_rel
 175
 176          if(allocated(data)) deallocate(data)
 177          allocate(data(fg_data%nx,fg_data%ny))
 178          data = fg_data%slab
 179
 180          tmp_dims(1) = fg_data%nx
 181          tmp_dims(2) = fg_data%ny
 182          do dim = 1,2
 183            i = 1
 184            CHECKDIMS2 : DO
 185              if (dval(i) == tmp_dims(dim) ) then
 186                dimids(dim) = i
 187                EXIT CHECKDIMS2
 188              end if
 189              i = i+1
 190              CYCLE CHECKDIMS2
 191            END DO CHECKDIMS2
 192          end do
 193
 194          call check(nf_redef(ncid))
 195          call check(nf_def_var(ncid,name,NF_REAL,NDIMS,dimids,varid))
 196          call check(nf_put_att_text(ncid,varid,DATEV,132,date))
 197          call check(nf_put_att_real(ncid,varid,FCSTV,nf_float,1,fcst))
 198          call check(nf_put_att_text(ncid,varid,SOURCEV,132,source))
 199          call check(nf_put_att_text(ncid,varid,FIELDV,132,field))
 200          call check(nf_put_att_text(ncid,varid,UNITSV,132,units))
 201          call check(nf_put_att_text(ncid,varid,DESCV,132,desc))
 202          call check(nf_put_att_int (ncid,varid,LEVELV,nf_int,1,ablevel))
 203          call check(nf_put_att_int (ncid,varid,NX,nf_int,1,fg_data%nx))
 204          call check(nf_put_att_int (ncid,varid,NY,nf_int,1,fg_data%ny))
 205          call check(nf_put_att_int (ncid,varid,IPROJ,nf_int,1,proj))
 206          call check(nf_put_att_real(ncid,varid,FILLVALUE,nf_real,1,-1e30))
 207
 208
 209          if(proj == 0) then
 210            call check(nf_put_att_real(ncid,varid,STARTI,nf_float,1,si))
 211            call check(nf_put_att_real(ncid,varid,STARTJ,nf_float,1,sj))
 212            call check(nf_put_att_real(ncid,varid,STARTLAT,nf_float,1,slat))
 213            call check(nf_put_att_real(ncid,varid,STARTLON,nf_float,1,slon))
 214            call check(nf_put_att_real(ncid,varid,DELTALAT,nf_float,1,dlat))
 215            call check(nf_put_att_real(ncid,varid,DELTALON,nf_float,1,dlon))
 216            call check(nf_put_att_real(ncid,varid,EARTH_RADIUS,nf_float,1,radius))
 217          else if(proj == 1) then
 218            call check(nf_put_att_real(ncid,varid,STARTI,nf_float,1,si))
 219            call check(nf_put_att_real(ncid,varid,STARTJ,nf_float,1,sj))
 220            call check(nf_put_att_real(ncid,varid,STARTLAT,nf_float,1,slat))
 221            call check(nf_put_att_real(ncid,varid,STARTLON,nf_float,1,slon))
 222            call check(nf_put_att_real(ncid,varid,DX,nf_float,1,dxn))
 223            call check(nf_put_att_real(ncid,varid,DY,nf_float,1,dyn))
 224            call check(nf_put_att_real(ncid,varid,TRUELAT1,nf_float,1,tlat1))
 225            call check(nf_put_att_real(ncid,varid,EARTH_RADIUS,nf_float,1,radius))
 226          else if(proj == 3) then
 227            call check(nf_put_att_real(ncid,varid,STARTI,nf_float,1,si))
 228            call check(nf_put_att_real(ncid,varid,STARTJ,nf_float,1,sj))
 229            call check(nf_put_att_real(ncid,varid,STARTLAT,nf_float,1,slat))
 230            call check(nf_put_att_real(ncid,varid,STARTLON,nf_float,1,slon))
 231            call check(nf_put_att_real(ncid,varid,DX,nf_float,1,dxn))
 232            call check(nf_put_att_real(ncid,varid,DY,nf_float,1,dyn))
 233            call check(nf_put_att_real(ncid,varid,XLONC,nf_float,1,xloncen))
 234            call check(nf_put_att_real(ncid,varid,TRUELAT1,nf_float,1,tlat1))
 235            call check(nf_put_att_real(ncid,varid,TRUELAT2,nf_float,1,tlat2))
 236            call check(nf_put_att_real(ncid,varid,EARTH_RADIUS,nf_float,1,radius))
 237          else if(proj == 4) then
 238            call check(nf_put_att_real(ncid,varid,STARTI,nf_float,1,si))
 239            call check(nf_put_att_real(ncid,varid,STARTJ,nf_float,1,sj))
 240            call check(nf_put_att_real(ncid,varid,STARTLAT,nf_float,1,slat))
 241            call check(nf_put_att_real(ncid,varid,STARTLON,nf_float,1,slon))
 242            call check(nf_put_att_real(ncid,varid,NLATS,nf_float,1,nlat))
 243            call check(nf_put_att_real(ncid,varid,DELTALON,nf_float,1,dlon))
 244            call check(nf_put_att_real(ncid,varid,EARTH_RADIUS,nf_float,1,radius))
 245          else if(proj == 5) then
 246            call check(nf_put_att_real(ncid,varid,STARTI,nf_float,1,si))
 247            call check(nf_put_att_real(ncid,varid,STARTJ,nf_float,1,sj))
 248            call check(nf_put_att_real(ncid,varid,STARTLAT,nf_float,1,slat))
 249            call check(nf_put_att_real(ncid,varid,STARTLON,nf_float,1,slon))
 250            call check(nf_put_att_real(ncid,varid,DX,nf_float,1,dxn))
 251            call check(nf_put_att_real(ncid,varid,DY,nf_float,1,dyn))
 252            call check(nf_put_att_real(ncid,varid,XLONC,nf_float,1,xloncen))
 253            call check(nf_put_att_real(ncid,varid,TRUELAT1,nf_float,1,tlat1))
 254            call check(nf_put_att_real(ncid,varid,EARTH_RADIUS,nf_float,1,radius))
 255          end if
 256
 257          call check(nf_enddef(ncid))
 258          call check(nf_put_var_real(ncid,varid,data))
 259
 260          call read_next_met_field(fg_data,istatus)
 261
 262       end do
 263
 264       call read_met_close()
 265
 266    end if
 267
 268    call check(nf_close(ncid))
 269
 270    print *,'SUCCESSFUL COMPLETION OF PROGRAM INT2NC, ',trim(nfile),' WRITTEN.'
 271
 272 contains
 273   subroutine check(status)
 274     integer, intent ( in) :: status
 275
 276     if(status /= nf_noerr) then
 277       print *, trim(nf_strerror(status))
 278       stop "Stopped"
 279     end if
 280   end subroutine check
 281 end program int2nc