c interpolate_matrix subroutine c This program interpolates the matrix grid given the proper parameters c created 9/96 subroutine interpolate_matrix(grid,grerror,nx,ny,xres,yres, 1 xleft,xright,ytop,ybot) parameter (ngx=1000, ngy=700) integer nx,ny,choice real xres,yres,xleft,xright,ytop,ybot real grid(ngy,ngx),grerror(ngy,ngx),degrees(ngx),position(ngx) call prompt_interp(choice) if (choice .eq. 1) then c assign distance and lat or long to each column in the grid call interp_sdis(position,degrees,grid,nx,xres) c interpolate distance grid to desired lat or lon grid call interp_data(position,degrees,grid,grerror,nx,ny, 1 xleft,xright,xres) c print *,nx,ny,xleft,xright,xres end if return end c prompt_interp subroutine c This program asks the user whether they actually want to interpolate the c matrix or not. created 9/96 subroutine prompt_interp(choice) c c <- choice c character*1, yesorno integer choice print *," " print *,"Do you want to interpolate from kilometers to degrees? (y/n)" read (5,*)yesorno if (yesorno .eq. "y") then choice = 1 else choice = 2 end if return end c interpolate_sdis subroutine c this program finds the latitude or longitude of each grid column c using the stations.dat file and header from gridded file c created 9/96 c subroutine interp_sdis(grdis,grdeg,grid,nx,xres) c c -> grid input gridded data file c -> nx input number of points in x c -> xres resolution in x c <- grdis distance for each grid column c <- grdeg latitude or longitude of each grid column c parameter (ngx=1000, ngy=700) character*40, statn integer nx,choice,iost,fnum,snum real xres,fdeg,sdeg,grid(ngy,ngx),grdeg(ngx) real fdis,flat,flon,sdis,slat,slon,junk,grdis(ngx) c prompt user for input file print *,"Enter the name of the stations file" read *,statn iost = 2 open(iost,file=statn,status='old',form='formatted',iostat=icd) if (icd .ne. 0) then write(6,*) 'error opening stations file iostat= ',icd stop 99 end if c Ask user choice of latitude or longitude choice = 1 print *,"Interpolate to longitude (1) or latitude (2)?" print *,"(default = longitude)" read *,choice read(iost,*)fnum,fdis,flat,flon,junk read(iost,*)snum,sdis,slat,slon,junk if ( flon .lt. 0) flon = flon + 360.0 if ( slon .lt. 0) slon = slon + 360.0 c print *,flon,slon,choice if (choice .eq. 2) then fdeg = flat else fdeg = flon endif c c set the first lat or lon coordinate grdeg(1) = fdeg grdis(1) = 0 c c c start main loop do j = 2,nx c find the current position grdis(j) = (j-1)*xres c while the current pointer is in between the left and right pointers c do while (fdis .lt. grdis(j)) if (choice .eq. 2) then fdeg = flat sdeg = slat else fdeg = flon sdeg = slon endif c print *,j,fnum,snum,fdeg,sdeg,fdis,sdis,grdis(j) c if (sdis .gt. grdis(j)) then grdeg(j) = fdeg + ( sdeg-fdeg)*(grdis(j)-fdis)/(sdis-fdis) go to 20 else c save current station to old station fnum = snum fdis = sdis flat = slat flon = slon endif read(iost,*,end=100) snum,sdis,slat,slon,junk if ( slon .lt. 0) slon = slon + 360.0 enddo c end-of-while c 20 continue c print *,j,grdeg(j),grdis(j) enddo c end-of-main-loop 100 continue if (choice .eq. 2) then grdeg(nx) = slat else grdeg(nx) = slon endif c c close(iost) return end c interpolate data by the given parameters subroutine interp_data(grdis,grdeg,grid,grerror,nx,ny, 1 xleft,xright,xres) c c <-> grid input grid, changed on output c <-> grerror input error grid, changed on output c <-> nx input columns, changed on output c -> ny input rows c <-> xleft input left distance, changed to degrees on output c <-> xright input right distance, changed to degrees on output c --> grdis input distance for each column c --> grdeg input degrees for each column parameter (ngx=1000, ngy=700) integer nx,ny,i,j,lctr,nctr real grdis(ngx),grdeg(ngx),newgrid(ngy,ngx),grid(ngy,ngx), 1 newerr(ngy,ngx),grerror(ngy,ngx),newlon(ngx) real xleft,xright,resolutn,start_l,current_l,end_l c resolutn = 1. start_l = grdeg(1) end_l = grdeg(nx) c c ask user about other parameters print *,"Enter the resolution for output in decimal degrees" print *,"(default = 1 degree)" read *,resolutn c print *,'Start and end positions are: ',grdeg(1),grdeg(nx) print *,' ' print *,"Enter the starting latitude or longitude for output" print *,"(default = Nearest degree to starting point)" read *,start_l c print *,"Enter the ending latitude or longitude for output" print *,"(default = Existing ending point)" read *,end_l c c check to see if the resolution is valid if ((resolutn .le. 0) .or. (resolutn .gt. 1 (grdeg(nx)-grdeg(1)))) 2 resolutn = 1.0 c c check to see if the starting point is valid if ((start_l .lt. grdeg(1)) .or. 1 (start_l .gt. grdeg(nx))) 2 start_l = int(grdeg(1)) + 1 c c check to see if the ending point is valid if ((end_l .gt. grdeg(nx)) .or. 1 (end_l .lt. grdeg(1))) 2 end_l = grdeg(nx) c c check to see if going in the right direction in the matrix if ( (end_l - start_l) * (grdeg(nx) - grdeg(1) ) . 1 lt. 0.) then print *,'Going in wrong direction - flip' sn = start_l start_l = end_l end_l = sn endif print *,'Starting and ending degrees: ',start_l, end_l c c initialize counters lctr = 2 nctr = 1 current_l = start_l c c while the lat/lon counter has not reached the end of the spectrum do while ((lctr .le. nx) .and. (current_l .le. end_l)) c if the current lat/lon pointer is in between a greater and lesser degree if (grdeg(lctr) .ge. current_l) then fac = (current_l-grdeg(lctr-1))/ 1 (grdeg(lctr)-grdeg(lctr-1)) do i = 1, ny newgrid(i,nctr) = grid(i,lctr-1) + 1 ( grid(i,lctr)-grid(i,lctr-1) ) * fac newerr(i,nctr) = grerror(i,lctr-1) + 1 ( grerror(i,lctr)-grerror(i,lctr-1) ) * fac enddo newlon(nctr) = current_l nctr = nctr + 1 current_l = current_l + resolutn else lctr = lctr + 1 end if enddo c end-of-while c c replace old grid with new grid do i = 1,ny do j = 1,nctr - 1 grid(i,j) = newgrid(i,j) grerror(i,j) = newerr(i,j) enddo enddo c c assign new values to the column length and resolution nx = nctr - 1 xres = resolutn c find the new starting and ending positions xleft = start_l xright = end_l return end