subroutine datint(fuldates, nful, hstart, ntimes, interval, out_format, prefix) ! ! !*****************************************************************************! ! ! ! interpolate missing data in time ! out_format: requested output format ! ! !*****************************************************************************! use gridinfo use storage_module use module_debug use misc_definitions_module implicit none integer :: nful integer :: interval character(len=*), dimension(nful) :: fuldates character(len=*) :: hstart integer :: ntimes character(len=24) :: hdate = "0000-00-00_00:00:00.0000" character(len=24) :: hdate_output, jdate character(len=9) :: field character(len=25) :: units character(len=46) :: desc character(LEN=3) :: out_format character(LEN=MAX_FILENAME_LEN) :: prefix real :: xfcst real :: level real, allocatable, dimension(:,:) :: scr2d, bfr2d integer :: iful, intervala, intervalb, ifv real :: awt integer :: itime ! DATELEN: length of date strings to use for our output file names. integer :: datelen ! Decide the length of date strings to use for output file names. ! DATELEN is 13 for hours, 16 for minutes, and 19 for seconds. if (mod(interval,3600) == 0) then datelen = 13 else if (mod(interval, 60) == 0) then datelen = 16 else datelen = 19 end if call mprintf(.true.,STDOUT,"Subroutine DATINT: Interpolating 3-d files to fill in any missing data...") call mprintf(.true.,LOGFILE,"Subroutine DATINT: Interpolating 3-d files to fill in any missing data...") TIMELOOP : do itime = 1, ntimes call geth_newdate(hdate(1:19), hstart(1:19), (itime-1)*interval) call mprintf(.true.,STDOUT,"Looking for data at time %s",s1=hdate(1:datelen)) call mprintf(.true.,LOGFILE,"Looking for data at time %s",s1=hdate(1:datelen)) do iful = 1, nful if (fuldates(iful).eq.hdate) then call mprintf(.true.,STDOUT,"Found file: %s:%s", & s1=trim(prefix),s2=hdate(1:datelen)) call mprintf(.true.,LOGFILE,"Found file: %s:%s", & s1=trim(prefix),s2=hdate(1:datelen)) cycle TIMELOOP else if ((fuldates(iful).lt.hdate) .and. & (fuldates(iful+1).gt.hdate) )then call mprintf(.true.,STDOUT,"Found surrounding files: %s: %s and %s: %s", & s1=trim(prefix),s2=fuldates(iful)(1:datelen), & s3=trim(prefix),s4=fuldates(iful+1)(1:datelen)) call mprintf(.true.,LOGFILE,"Found surrounding files: %s: %s and %s: %s", & s1=trim(prefix),s2=fuldates(iful)(1:datelen), & s3=trim(prefix),s4=fuldates(iful+1)(1:datelen)) call mprintf(.true.,STDOUT,"Interpolating to create file: %s: %s", & s1=trim(prefix),s2=hdate(1:datelen)) call mprintf(.true.,LOGFILE,"Interpolating to create file: %s: %s", & s1=trim(prefix),s2=hdate(1:datelen)) call geth_idts(hdate(1:19), fuldates(iful)(1:19), intervalA) call mprintf(.true.,STDOUT,"A Time Difference = %f",f1=float(intervalA) / 3600.) call mprintf(.true.,LOGFILE,"A Time Difference = %f",f1=float(intervalA) / 3600.) call geth_idts(fuldates(iful+1)(1:19), hdate(1:19), intervalB) call mprintf(.true.,STDOUT,"B Time Difference = %f",f1=float(intervalB) / 3600.) call mprintf(.true.,LOGFILE,"B Time Difference = %f",f1=float(intervalB) / 3600.) AWT = 1. - (float(intervalA)/float(intervalA+intervalB)) open(10, file=trim(prefix)//':'//fuldates(iful)(1:datelen), form='unformatted', & status='old') call clear_storage READLOOP1 : do read(10, end=44) ifv if ( ifv .eq. 5) then ! WPS read (10) jdate, xfcst, map%source, field, units, desc, level, & map%nx, map%ny, map%igrid select case (map%igrid) case (0, 4) read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%r_earth case (3) read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, & map%lov, map%truelat1, map%truelat2, map%r_earth case (5) read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, & map%lov, map%truelat1, map%r_earth case default call mprintf(.true.,ERROR, & "Unrecognized map%%igrid: %i in DATINT 1",i1=map%igrid) end select read (10) map%grid_wind else if ( ifv .eq. 4 ) then ! SI read (10) jdate, xfcst, map%source, field, units, desc, level, & map%nx, map%ny, map%igrid select case (map%igrid) case (0, 4) read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx case (3) read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, & map%lov, map%truelat1, map%truelat2 case (5) read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, & map%lov, map%truelat1 case default call mprintf(.true.,ERROR, & "Unrecognized map%%igrid: %i in DATINT 2",i1=map%igrid) end select else if ( ifv .eq. 3 ) then ! MM5 read(10) jdate, xfcst, field, units, desc, level,& map%nx, map%ny, map%igrid select case (map%igrid) case (3) ! lamcon read (10) map%lat1, map%lon1, map%dx, map%dy, map%lov, & map%truelat1, map%truelat2 case (5) ! Polar Stereographic read (10) map%lat1, map%lon1, map%dx, map%dy, map%lov, & map%truelat1 case (0, 4) ! lat/lon read (10) map%lat1, map%lon1, map%dy, map%dx case (1) ! Mercator read (10) map%lat1, map%lon1, map%dy, map%dx, map%truelat1 case default call mprintf(.true.,ERROR, & "Unrecognized map%%igrid: %i in DATINT 3",i1=map%igrid) end select else call mprintf(.true.,ERROR, & "Unknown out_format: %i in DATINT ",i1=ifv) endif allocate(scr2d(map%nx, map%ny)) read (10) scr2d call put_storage(nint(level), field, scr2d, map%nx, map%ny) deallocate(scr2d) enddo READLOOP1 44 close(10) open(10, file=trim(prefix)//':'//fuldates(iful+1)(1:datelen), status='old', & form = 'unformatted') open(11, file=trim(prefix)//':'//hdate(1:datelen), status='new', form='unformatted') READLOOP2 : do read (10,END=45) ifv if ( ifv .eq. 5) then ! WPS read (10) jdate, xfcst, map%source, field, units, desc, level, & map%nx, map%ny, map%igrid select case (map%igrid) case (0, 4) read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%r_earth case (3) read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, & map%lov, map%truelat1, map%truelat2, map%r_earth case (5) read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, & map%lov, map%truelat1, map%r_earth case default call mprintf(.true.,ERROR, & "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid) end select read (10) map%grid_wind else if ( ifv .eq. 4 ) then ! SI read (10) jdate, xfcst, map%source, field, units, desc, level, & map%nx, map%ny, map%igrid select case (map%igrid) case (0, 4) read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx case (1) read(10) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%truelat1 case (3) read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, & map%lov, map%truelat1, map%truelat2 case (5) read (10) map%startloc, map%lat1, map%lon1, map%dx, map%dy, & map%lov, map%truelat1 case default call mprintf(.true.,ERROR, & "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid) end select else if ( ifv .eq. 3 ) then ! MM5 read(10) jdate, xfcst, field, units, desc, level,& map%nx, map%ny, map%igrid select case (map%igrid) case (3) ! lamcon read (10) map%lat1, map%lon1, map%dx, map%dy, map%lov, & map%truelat1, map%truelat2 case (5) ! Polar Stereographic read (10) map%lat1, map%lon1, map%dx, map%dy, map%lov, & map%truelat1 case (0, 4) ! lat/lon read (10) map%lat1, map%lon1, map%dy, map%dx case (1) ! Mercator read (10) map%lat1, map%lon1, map%dy, map%dx, map%truelat1 case default call mprintf(.true.,ERROR, & "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid) end select else call mprintf(.true.,ERROR, & "Unknown out_format: %i in DATINT ",i1=ifv) endif allocate(scr2d(map%nx, map%ny)) read (10) scr2d if (is_there(nint(level), field)) then allocate(bfr2d(map%nx,map%ny)) call get_storage(nint(level), field, bfr2d, map%nx, map%ny) scr2d = bfr2d * (AWT) + scr2d * (1.-AWT) hdate_output = hdate if (out_format(1:2) .eq. 'SI') then write(11) ifv write(11) hdate_output, xfcst, map%source, field, units, desc, & level, map%nx, map%ny, map%igrid if (map%igrid == 0 .or. map%igrid == 4) then write(11) map%startloc, map%lat1, map%lon1, map%dy, map%dx elseif (map%igrid == 1) then write(11) map%startloc, map%lat1, map%lon1, map%dy, map%dx, map%truelat1 elseif (map%igrid == 3) then write (11) map%startloc, map%lat1, map%lon1, map%dx, map%dy, & map%lov, map%truelat1, map%truelat2 elseif (map%igrid == 5) then write (11) map%startloc, map%lat1, map%lon1, map%dx, map%dy, & map%lov, map%truelat1 else call mprintf(.true.,ERROR, & "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid) endif else if (out_format(1:2) .eq. 'WP') then write(11) ifv write(11) hdate_output, xfcst, map%source, field, units, desc, & level, map%nx, map%ny, map%igrid if (map%igrid == 0 .or. map%igrid == 4) then write(11) map%startloc, map%lat1, map%lon1, map%dy, map%dx, & map%r_earth elseif (map%igrid == 1) then write(11) map%startloc, map%lat1, map%lon1, map%dy, map%dx, & map%truelat1, map%r_earth elseif (map%igrid == 3) then write (11) map%startloc, map%lat1, map%lon1, map%dx, map%dy, & map%lov, map%truelat1, map%truelat2, map%r_earth elseif (map%igrid == 5) then write (11) map%startloc, map%lat1, map%lon1, map%dx, map%dy, & map%lov, map%truelat1, map%r_earth else call mprintf(.true.,ERROR, & "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid) endif write(11) map%grid_wind else if (out_format(1:2) .eq. 'MM') then write (11) ifv write (11) hdate_output, xfcst, field, units, Desc, level,& map%nx, map%ny, map%igrid if (map%igrid .eq. 3) then ! lamcon write (11) map%lat1, map%lon1, map%dx, map%dy, map%lov, & map%truelat1, map%truelat2 elseif (map%igrid .eq. 5) then ! Polar Stereographic write (11) map%lat1, map%lon1, map%dx, map%dy, map%lov, & map%truelat1 elseif (map%igrid .eq. 0 .or. map%igrid .eq. 4)then ! lat/lon write (11) map%lat1, map%lon1, map%dy, map%dx elseif (map%igrid.eq.1)then ! Mercator write (11) map%lat1, map%lon1, map%dy, map%dx, map%truelat1 else call mprintf(.true.,ERROR, & "Unrecognized map%%igrid: %i in DATINT ",i1=map%igrid) endif endif write(11) scr2d else call mprintf(.true.,ERROR, & "hdate = %s , fuldates = %s %s, Field = %s",s1=hdate,s2=fuldates(iful),s3=fuldates(iful+1),s4=field) endif deallocate(scr2d, bfr2d) enddo READLOOP2 45 close(10) close(11) cycle TIMELOOP endif enddo call mprintf(.true.,ERROR, & "Data not found: %s",s1=hdate) enddo TIMELOOP call mprintf(.true.,STDOUT, & "End Subroutine DATINT.") call mprintf(.true.,LOGFILE, & "End Subroutine DATINT.") end subroutine datint