metgrid/src/read_met_module.F

   1 module read_met_module
   2
   3    use constants_module
   4    use module_debug
   5    use misc_definitions_module
   6    use met_data_module
   7
   8    ! State variables?
   9    integer :: input_unit
  10    character (len=MAX_FILENAME_LEN) :: filename
  11
  12    contains
  13
  14    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  15    ! Name: read_met_init
  16    !
  17    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  18    subroutine read_met_init(fg_source, source_is_constant, datestr, istatus)
  19
  20       implicit none
  21
  22       ! Arguments
  23       integer, intent(out) :: istatus
  24       logical, intent(in) :: source_is_constant
  25       character (len=*), intent(in) :: fg_source
  26       character (len=*), intent(in) :: datestr
  27
  28       ! Local variables
  29       integer :: io_status
  30       logical :: is_used
  31
  32       istatus = 0
  33
  34       !  1) BUILD FILENAME BASED ON TIME
  35       filename = ' '
  36       if (.not. source_is_constant) then
  37          write(filename, '(a)') trim(fg_source)//':'//trim(datestr)
  38       else
  39          write(filename, '(a)') trim(fg_source)
  40       end if
  41
  42       !  2) OPEN FILE
  43       do input_unit=10,100
  44          inquire(unit=input_unit, opened=is_used)
  45          if (.not. is_used) exit
  46       end do
  47       call mprintf((input_unit > 100),ERROR,'In read_met_init(), couldn''t find an available Fortran unit.')
  48       open(unit=input_unit, file=trim(filename), status='old', form='unformatted', iostat=io_status)
  49
  50       if (io_status > 0) istatus = 1
  51
  52       return
  53
  54
  55    end subroutine read_met_init
  56
  57
  58    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  59    ! Name: read_next_met_field
  60    !
  61    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  62    subroutine read_next_met_field(fg_data, istatus)
  63
  64       implicit none
  65
  66       ! Arguments
  67       type (met_data), intent(inout) :: fg_data
  68       integer, intent(out) :: istatus
  69
  70       ! Local variables
  71       character (len=8) :: startloc
  72
  73       istatus = 1
  74
  75       !  1) READ FORMAT VERSION
  76       read(unit=input_unit,err=1001,end=1001) fg_data % version
  77
  78       ! PREGRID
  79       if (fg_data % version == 3) then
  80
  81          read(unit=input_unit) fg_data % hdate, &
  82                                fg_data % xfcst, &
  83                                fg_data % field, &
  84                                fg_data % units, &
  85                                fg_data % desc,  &
  86                                fg_data % xlvl,  &
  87                                fg_data % nx,    &
  88                                fg_data % ny,    &
  89                                fg_data % iproj
  90
  91          fg_data % map_source = ' '
  92
  93          if (fg_data % field == 'HGT      ') fg_data % field = 'GHT      '
  94
  95          fg_data % starti = 1.0
  96          fg_data % startj = 1.0
  97
  98          ! Cylindrical equidistant
  99          if (fg_data % iproj == 0) then
 100             fg_data % iproj = PROJ_LATLON
 101             read(unit=input_unit,err=1001,end=1001) fg_data % startlat, &
 102                                                     fg_data % startlon, &
 103                                                     fg_data % deltalat, &
 104                                                     fg_data % deltalon
 105
 106          ! Mercator
 107          else if (fg_data % iproj == 1) then
 108             fg_data % iproj = PROJ_MERC
 109             read(unit=input_unit,err=1001,end=1001) fg_data % startlat, &
 110                                                     fg_data % startlon, &
 111                                                     fg_data % dx,       &
 112                                                     fg_data % dy,       &
 113                                                     fg_data % truelat1
 114
 115          ! Lambert conformal
 116          else if (fg_data % iproj == 3) then
 117             fg_data % iproj = PROJ_LC
 118             read(unit=input_unit,err=1001,end=1001) fg_data % startlat, &
 119                                                     fg_data % startlon, &
 120                                                     fg_data % dx,       &
 121                                                     fg_data % dy,       &
 122                                                     fg_data % xlonc,    &
 123                                                     fg_data % truelat1, &
 124                                                     fg_data % truelat2
 125
 126          ! Polar stereographic
 127          else if (fg_data % iproj == 5) then
 128             fg_data % iproj = PROJ_PS
 129             read(unit=input_unit,err=1001,end=1001) fg_data % startlat, &
 130                                                     fg_data % startlon, &
 131                                                     fg_data % dx,       &
 132                                                     fg_data % dy,       &
 133                                                     fg_data % xlonc,    &
 134                                                     fg_data % truelat1
 135
 136          ! ?????????
 137          else
 138             call mprintf(.true.,ERROR,'Unrecognized projection code %i when reading from %s', &
 139                          i1=fg_data % iproj, s1=filename)
 140
 141          end if
 142
 143          fg_data % earth_radius = EARTH_RADIUS_M / 1000.
 144
 145 #if (defined _GEOGRID) || (defined _METGRID)
 146          fg_data % dx = fg_data % dx * 1000.
 147          fg_data % dy = fg_data % dy * 1000.
 148
 149          if (fg_data % xlonc    > 180.) fg_data % xlonc    = fg_data%xlonc    - 360.
 150
 151          if (fg_data % startlon > 180.) fg_data % startlon = fg_data%startlon - 360.
 152
 153          if (fg_data % startlat < -90.) fg_data % startlat = -90.
 154          if (fg_data % startlat >  90.) fg_data % startlat = 90.
 155 #endif
 156
 157          fg_data % is_wind_grid_rel = .true.
 158
 159          allocate(fg_data % slab(fg_data % nx, fg_data % ny))
 160          read(unit=input_unit,err=1001,end=1001) fg_data % slab
 161
 162          istatus = 0
 163
 164       ! GRIB_PREP
 165       else if (fg_data % version == 4) then
 166
 167          read(unit=input_unit) fg_data % hdate,      &
 168                                fg_data % xfcst,      &
 169                                fg_data % map_source, &
 170                                fg_data % field,      &
 171                                fg_data % units,      &
 172                                fg_data % desc,       &
 173                                fg_data % xlvl,       &
 174                                fg_data % nx,         &
 175                                fg_data % ny,         &
 176                                fg_data % iproj
 177
 178          if (fg_data % field == 'HGT      ') fg_data % field = 'GHT      '
 179
 180          ! Cylindrical equidistant
 181          if (fg_data % iproj == 0) then
 182             fg_data % iproj = PROJ_LATLON
 183             read(unit=input_unit,err=1001,end=1001) startloc, &
 184                                                     fg_data % startlat, &
 185                                                     fg_data % startlon, &
 186                                                     fg_data % deltalat, &
 187                                                     fg_data % deltalon
 188
 189          ! Mercator
 190          else if (fg_data % iproj == 1) then
 191             fg_data % iproj = PROJ_MERC
 192             read(unit=input_unit,err=1001,end=1001) startloc, &
 193                                                     fg_data % startlat, &
 194                                                     fg_data % startlon, &
 195                                                     fg_data % dx,       &
 196                                                     fg_data % dy,       &
 197                                                     fg_data % truelat1
 198
 199          ! Lambert conformal
 200          else if (fg_data % iproj == 3) then
 201             fg_data % iproj = PROJ_LC
 202             read(unit=input_unit,err=1001,end=1001) startloc, &
 203                                                     fg_data % startlat, &
 204                                                     fg_data % startlon, &
 205                                                     fg_data % dx,       &
 206                                                     fg_data % dy,       &
 207                                                     fg_data % xlonc,    &
 208                                                     fg_data % truelat1, &
 209                                                     fg_data % truelat2
 210
 211          ! Polar stereographic
 212          else if (fg_data % iproj == 5) then
 213             fg_data % iproj = PROJ_PS
 214             read(unit=input_unit,err=1001,end=1001) startloc, &
 215                                                     fg_data % startlat, &
 216                                                     fg_data % startlon, &
 217                                                     fg_data % dx,       &
 218                                                     fg_data % dy,       &
 219                                                     fg_data % xlonc,    &
 220                                                     fg_data % truelat1
 221
 222          ! ?????????
 223          else
 224             call mprintf(.true.,ERROR,'Unrecognized projection code %i when reading from %s', &
 225                          i1=fg_data % iproj, s1=filename)
 226
 227          end if
 228
 229          if (startloc == 'CENTER  ') then
 230             fg_data % starti = real(fg_data % nx)/2.
 231             fg_data % startj = real(fg_data % ny)/2.
 232          else if (startloc == 'SWCORNER') then
 233             fg_data % starti = 1.0
 234             fg_data % startj = 1.0
 235          end if
 236
 237          fg_data % earth_radius = EARTH_RADIUS_M / 1000.
 238
 239 #if (defined _GEOGRID) || (defined _METGRID)
 240          fg_data % dx = fg_data % dx * 1000.
 241          fg_data % dy = fg_data % dy * 1000.
 242
 243          if (fg_data % xlonc    > 180.) fg_data % xlonc    = fg_data % xlonc    - 360.
 244
 245          if (fg_data % startlon > 180.) fg_data % startlon = fg_data % startlon - 360.
 246
 247          if (fg_data % startlat < -90.) fg_data % startlat = -90.
 248          if (fg_data % startlat >  90.) fg_data % startlat = 90.
 249 #endif
 250
 251          fg_data % is_wind_grid_rel = .true.
 252
 253          allocate(fg_data % slab(fg_data % nx, fg_data % ny))
 254          read(unit=input_unit,err=1001,end=1001) fg_data % slab
 255
 256          istatus = 0
 257
 258       ! WPS
 259       else if (fg_data % version == 5) then
 260
 261          read(unit=input_unit) fg_data % hdate,      &
 262                                fg_data % xfcst,      &
 263                                fg_data % map_source, &
 264                                fg_data % field,      &
 265                                fg_data % units,      &
 266                                fg_data % desc,       &
 267                                fg_data % xlvl,       &
 268                                fg_data % nx,         &
 269                                fg_data % ny,         &
 270                                fg_data % iproj
 271
 272          if (fg_data % field == 'HGT      ') fg_data % field = 'GHT      '
 273
 274          ! Cylindrical equidistant
 275          if (fg_data % iproj == 0) then
 276             fg_data % iproj = PROJ_LATLON
 277             read(unit=input_unit,err=1001,end=1001) startloc, &
 278                                                     fg_data % startlat, &
 279                                                     fg_data % startlon, &
 280                                                     fg_data % deltalat, &
 281                                                     fg_data % deltalon, &
 282                                                     fg_data % earth_radius
 283
 284          ! Mercator
 285          else if (fg_data % iproj == 1) then
 286             fg_data % iproj = PROJ_MERC
 287             read(unit=input_unit,err=1001,end=1001) startloc, &
 288                                                     fg_data % startlat, &
 289                                                     fg_data % startlon, &
 290                                                     fg_data % dx,       &
 291                                                     fg_data % dy,       &
 292                                                     fg_data % truelat1, &
 293                                                     fg_data % earth_radius
 294
 295          ! Lambert conformal
 296          else if (fg_data % iproj == 3) then
 297             fg_data % iproj = PROJ_LC
 298             read(unit=input_unit,err=1001,end=1001) startloc, &
 299                                                     fg_data % startlat, &
 300                                                     fg_data % startlon, &
 301                                                     fg_data % dx,       &
 302                                                     fg_data % dy,       &
 303                                                     fg_data % xlonc,    &
 304                                                     fg_data % truelat1, &
 305                                                     fg_data % truelat2, &
 306                                                     fg_data % earth_radius
 307
 308          ! Gaussian
 309          else if (fg_data % iproj == 4) then
 310             fg_data % iproj = PROJ_GAUSS
 311             read(unit=input_unit,err=1001,end=1001) startloc, &
 312                                                     fg_data % startlat, &
 313                                                     fg_data % startlon, &
 314                                                     fg_data % deltalat, &
 315                                                     fg_data % deltalon, &
 316                                                     fg_data % earth_radius
 317
 318          ! Polar stereographic
 319          else if (fg_data % iproj == 5) then
 320             fg_data % iproj = PROJ_PS
 321             read(unit=input_unit,err=1001,end=1001) startloc, &
 322                                                     fg_data % startlat, &
 323                                                     fg_data % startlon, &
 324                                                     fg_data % dx,       &
 325                                                     fg_data % dy,       &
 326                                                     fg_data % xlonc,    &
 327                                                     fg_data % truelat1, &
 328                                                     fg_data % earth_radius
 329
 330          ! ?????????
 331          else
 332             call mprintf(.true.,ERROR,'Unrecognized projection code %i when reading from %s', &
 333                          i1=fg_data % iproj, s1=filename)
 334
 335          end if
 336
 337          if (startloc == 'CENTER  ') then
 338             fg_data % starti = real(fg_data % nx)/2.
 339             fg_data % startj = real(fg_data % ny)/2.
 340          else if (startloc == 'SWCORNER') then
 341             fg_data % starti = 1.0
 342             fg_data % startj = 1.0
 343          end if
 344
 345 #if (defined _GEOGRID) || (defined _METGRID)
 346          fg_data % dx = fg_data % dx * 1000.
 347          fg_data % dy = fg_data % dy * 1000.
 348
 349          if (fg_data % xlonc    > 180.) fg_data % xlonc    = fg_data % xlonc    - 360.
 350
 351          if (fg_data % startlon > 180.) fg_data % startlon = fg_data % startlon - 360.
 352
 353          if (fg_data % startlat < -90.) fg_data % startlat = -90.
 354          if (fg_data % startlat >  90.) fg_data % startlat =  90.
 355 #endif
 356
 357          read(unit=input_unit,err=1001,end=1001) fg_data % is_wind_grid_rel
 358
 359          allocate(fg_data % slab(fg_data % nx, fg_data % ny))
 360          read(unit=input_unit,err=1001,end=1001) fg_data % slab
 361
 362          istatus = 0
 363
 364       else
 365          call mprintf(.true.,ERROR,'Didn''t recognize format version of data in %s.\n'// &
 366                                    'Found version %i but expected either 3, 4, or 5. This could be an endian problem.', &
 367                                    s1=filename, i1=fg_data % version)
 368       end if
 369
 370       return
 371
 372    1001 return
 373
 374    end subroutine read_next_met_field
 375
 376
 377    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 378    ! Name: read_met_close
 379    !
 380    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 381    subroutine read_met_close()
 382
 383       implicit none
 384
 385       close(unit=input_unit)
 386       filename = 'UNINITIALIZED_FILENAME'
 387
 388    end subroutine read_met_close
 389
 390 end module read_met_module