module read_met_module use constants_module use module_debug use misc_definitions_module use met_data_module ! State variables? integer :: input_unit character (len=MAX_FILENAME_LEN) :: filename contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: read_met_init ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_met_init(fg_source, source_is_constant, datestr, istatus) implicit none ! Arguments integer, intent(out) :: istatus logical, intent(in) :: source_is_constant character (len=*), intent(in) :: fg_source character (len=*), intent(in) :: datestr ! Local variables integer :: io_status logical :: is_used istatus = 0 ! 1) BUILD FILENAME BASED ON TIME filename = ' ' if (.not. source_is_constant) then write(filename, '(a)') trim(fg_source)//':'//trim(datestr) else write(filename, '(a)') trim(fg_source) end if ! 2) OPEN FILE do input_unit=10,100 inquire(unit=input_unit, opened=is_used) if (.not. is_used) exit end do call mprintf((input_unit > 100),ERROR,'In read_met_init(), couldn''t find an available Fortran unit.') open(unit=input_unit, file=trim(filename), status='old', form='unformatted', iostat=io_status) if (io_status > 0) istatus = 1 return end subroutine read_met_init !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: read_next_met_field ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_next_met_field(fg_data, istatus) implicit none ! Arguments type (met_data), intent(inout) :: fg_data integer, intent(out) :: istatus ! Local variables character (len=8) :: startloc istatus = 1 ! 1) READ FORMAT VERSION read(unit=input_unit,err=1001,end=1001) fg_data % version ! PREGRID if (fg_data % version == 3) then read(unit=input_unit) fg_data % hdate, & fg_data % xfcst, & fg_data % field, & fg_data % units, & fg_data % desc, & fg_data % xlvl, & fg_data % nx, & fg_data % ny, & fg_data % iproj fg_data % map_source = ' ' if (fg_data % field == 'HGT ') fg_data % field = 'GHT ' fg_data % starti = 1.0 fg_data % startj = 1.0 ! Cylindrical equidistant if (fg_data % iproj == 0) then fg_data % iproj = PROJ_LATLON read(unit=input_unit,err=1001,end=1001) fg_data % startlat, & fg_data % startlon, & fg_data % deltalat, & fg_data % deltalon ! Mercator else if (fg_data % iproj == 1) then fg_data % iproj = PROJ_MERC read(unit=input_unit,err=1001,end=1001) fg_data % startlat, & fg_data % startlon, & fg_data % dx, & fg_data % dy, & fg_data % truelat1 ! Lambert conformal else if (fg_data % iproj == 3) then fg_data % iproj = PROJ_LC read(unit=input_unit,err=1001,end=1001) fg_data % startlat, & fg_data % startlon, & fg_data % dx, & fg_data % dy, & fg_data % xlonc, & fg_data % truelat1, & fg_data % truelat2 ! Polar stereographic else if (fg_data % iproj == 5) then fg_data % iproj = PROJ_PS read(unit=input_unit,err=1001,end=1001) fg_data % startlat, & fg_data % startlon, & fg_data % dx, & fg_data % dy, & fg_data % xlonc, & fg_data % truelat1 ! ????????? else call mprintf(.true.,ERROR,'Unrecognized projection code %i when reading from %s', & i1=fg_data % iproj, s1=filename) end if fg_data % earth_radius = EARTH_RADIUS_M / 1000. fg_data % is_wind_grid_rel = .true. allocate(fg_data % slab(fg_data % nx, fg_data % ny)) read(unit=input_unit,err=1001,end=1001) fg_data % slab istatus = 0 ! GRIB_PREP else if (fg_data % version == 4) then read(unit=input_unit) fg_data % hdate, & fg_data % xfcst, & fg_data % map_source, & fg_data % field, & fg_data % units, & fg_data % desc, & fg_data % xlvl, & fg_data % nx, & fg_data % ny, & fg_data % iproj if (fg_data % field == 'HGT ') fg_data % field = 'GHT ' ! Cylindrical equidistant if (fg_data % iproj == 0) then fg_data % iproj = PROJ_LATLON read(unit=input_unit,err=1001,end=1001) startloc, & fg_data % startlat, & fg_data % startlon, & fg_data % deltalat, & fg_data % deltalon ! Mercator else if (fg_data % iproj == 1) then fg_data % iproj = PROJ_MERC read(unit=input_unit,err=1001,end=1001) startloc, & fg_data % startlat, & fg_data % startlon, & fg_data % dx, & fg_data % dy, & fg_data % truelat1 ! Lambert conformal else if (fg_data % iproj == 3) then fg_data % iproj = PROJ_LC read(unit=input_unit,err=1001,end=1001) startloc, & fg_data % startlat, & fg_data % startlon, & fg_data % dx, & fg_data % dy, & fg_data % xlonc, & fg_data % truelat1, & fg_data % truelat2 ! Polar stereographic else if (fg_data % iproj == 5) then fg_data % iproj = PROJ_PS read(unit=input_unit,err=1001,end=1001) startloc, & fg_data % startlat, & fg_data % startlon, & fg_data % dx, & fg_data % dy, & fg_data % xlonc, & fg_data % truelat1 ! ????????? else call mprintf(.true.,ERROR,'Unrecognized projection code %i when reading from %s', & i1=fg_data % iproj, s1=filename) end if if (startloc == 'CENTER ') then fg_data % starti = real(fg_data % nx)/2. fg_data % startj = real(fg_data % ny)/2. else if (startloc == 'SWCORNER') then fg_data % starti = 1.0 fg_data % startj = 1.0 end if fg_data % earth_radius = EARTH_RADIUS_M / 1000. fg_data % is_wind_grid_rel = .true. allocate(fg_data % slab(fg_data % nx, fg_data % ny)) read(unit=input_unit,err=1001,end=1001) fg_data % slab istatus = 0 ! WPS else if (fg_data % version == 5) then read(unit=input_unit) fg_data % hdate, & fg_data % xfcst, & fg_data % map_source, & fg_data % field, & fg_data % units, & fg_data % desc, & fg_data % xlvl, & fg_data % nx, & fg_data % ny, & fg_data % iproj if (fg_data % field == 'HGT ') fg_data % field = 'GHT ' ! Cylindrical equidistant if (fg_data % iproj == 0) then fg_data % iproj = PROJ_LATLON read(unit=input_unit,err=1001,end=1001) startloc, & fg_data % startlat, & fg_data % startlon, & fg_data % deltalat, & fg_data % deltalon, & fg_data % earth_radius ! Mercator else if (fg_data % iproj == 1) then fg_data % iproj = PROJ_MERC read(unit=input_unit,err=1001,end=1001) startloc, & fg_data % startlat, & fg_data % startlon, & fg_data % dx, & fg_data % dy, & fg_data % truelat1, & fg_data % earth_radius ! Lambert conformal else if (fg_data % iproj == 3) then fg_data % iproj = PROJ_LC read(unit=input_unit,err=1001,end=1001) startloc, & fg_data % startlat, & fg_data % startlon, & fg_data % dx, & fg_data % dy, & fg_data % xlonc, & fg_data % truelat1, & fg_data % truelat2, & fg_data % earth_radius ! Gaussian else if (fg_data % iproj == 4) then fg_data % iproj = PROJ_GAUSS read(unit=input_unit,err=1001,end=1001) startloc, & fg_data % startlat, & fg_data % startlon, & fg_data % deltalat, & fg_data % deltalon, & fg_data % earth_radius ! Polar stereographic else if (fg_data % iproj == 5) then fg_data % iproj = PROJ_PS read(unit=input_unit,err=1001,end=1001) startloc, & fg_data % startlat, & fg_data % startlon, & fg_data % dx, & fg_data % dy, & fg_data % xlonc, & fg_data % truelat1, & fg_data % earth_radius ! ????????? else call mprintf(.true.,ERROR,'Unrecognized projection code %i when reading from %s', & i1=fg_data % iproj, s1=filename) end if if (startloc == 'CENTER ') then fg_data % starti = real(fg_data % nx)/2. fg_data % startj = real(fg_data % ny)/2. else if (startloc == 'SWCORNER') then fg_data % starti = 1.0 fg_data % startj = 1.0 end if read(unit=input_unit,err=1001,end=1001) fg_data % is_wind_grid_rel allocate(fg_data % slab(fg_data % nx, fg_data % ny)) read(unit=input_unit,err=1001,end=1001) fg_data % slab istatus = 0 else call mprintf(.true.,ERROR,'Didn''t recognize format version of data in %s.\n'// & 'Found version %i but expected either 3, 4, or 5. This could be an endian problem.', & s1=filename, i1=fg_data % version) end if return 1001 return end subroutine read_next_met_field !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: read_met_close ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine read_met_close() implicit none close(unit=input_unit) filename = 'UNINITIALIZED_FILENAME' end subroutine read_met_close end module read_met_module