subroutine llintp(f1,slon1,slat1,dlon1,dlat1,im1,jm1,i1,j1, + f2,slon2,slat2,dlon2,dlat2,im2,jm2,i2,j2, + izp,ierr) c c This routine linearly interpolates f1 on an evenly spaced c lat,lon grid to obtain f2 on a different lat,lon grid. c The subroutine arguments are defined as follows: c c f1(im1,jm1): Contains the function to be interpolated. The c first index represents longitude (increasing c eastward) and the second index represents c latitude (increasing northward). c c f2(im2,jm2): The interpolated function with indices defined c the same as for f1. c c slon1 The first longitude of f1 (deg E positive) c slat1 The first latitude of f1 (deg N positive) c dlon1 The longitude increment of f1 (deg) c dlat1 The latitude increment of f1 (deg) c slon2 The first longitude of f2 (deg E positive) c slat2 The first latitude of f2 (deg N positive) c dlon2 The longitude increment of f2 (deg) c dlat2 The latitude increment of f2 (deg) c c im1,jm1 The dimensions of f1 c im2,jm2 The dimensions of f2 c c i1,j1 The number of longitude,latitude points of f1 c to use in the interpolation c i2,j2 The number of longitude,latitude points of f2 c to interpolate c c izp Zonal Periodic flag: c =1 when f1 is periodic in the zonal direction c (normally used only when f1 spans 360 deg c of longitue) c =0 if not periodic c c ierr Error flag: =0 for normal return c =1 if f2 domain exceeds f1 domain c (non fatal) c =2 if indices i1,j1 or i2,j2 exceed c dimension of f1 or f2 (fatal) c =3 if dlon or dlat .le. 0.0 (fatal) c c Written: Jack Dostalek 2004 c Modified: Jack Dostalek 1/30/2006 c Correction applied to prevent a potential problem when c rlat is close to either pole c c dimension f1(im1,jm1),f2(im2,jm2) c c Initialize error flag ierr=0 c c Check indices and dlat,dlon if (i1 .gt. im1 .or. j1 .gt. jm1 .or. + i2 .gt. im2 .or. j2 .gt. jm2) then ierr=2 return endif c if (dlat1 .le. 0.0 .or. dlon1 .le. 0.0 .or. + dlat2 .le. 0.0 .or. dlon2 .le. 0.0) then ierr=3 return endif c c Specify needed constants pi = 3.14159265 dtr = pi/180.0 erad = 6371.0 adtr = erad*dtr c c Calculate min and max long,lat of f1 domain rlonn1 = slon1 rlonx1 = slon1 + dlon1*float(i1-1) rlatn1 = slat1 rlatx1 = slat1 + dlat1*float(j1-1) c c Start loop for f2 points do 10 j=1,j2 do 10 i=1,i2 rlon = slon2 + dlon2*float(i-1) rlat = slat2 + dlat2*float(j-1) c c Check if current f2 point is outside of f1 domain. c If yes, move the f2 point to the nearest point in the c f1 domain and set error flag. c if (izp .ne. 1) then c Adjust f2 longitude for case without zonal periodicity if (rlon .gt. rlonx1) then rlon = rlonx1 ierr=1 endif c if (rlon .lt. rlonn1) then rlon = rlonn1 ierr=1 endif else c Zonal periodic case if (rlon .ge. 360.0) rlon=rlon-360.0 endif c if (rlat .gt. rlatx1) then rlat = rlatx1 ierr=1 endif c if (rlat .lt. rlatn1) then rlat = rlatn1 ierr=1 endif c c Find the indices of the f1 point closest to, c but with lon,lat less than the current f2 point. i00 = 1 + ifix( (rlon-rlonn1)/dlon1 ) j00 = 1 + ifix( (rlat-rlatn1)/dlat1 ) if (i00 .lt. 1) i00= 1 if (izp .ne. 1) then if (i00 .gt. i1-1) i00=i1-1 endif if (j00 .lt. 1) j00= 1 if (j00 .gt. j1-1) j00=j1-1 c c Define the four f1 values to be used in the linear c interpolation. f00 = f1(i00 ,j00 ) f01 = f1(i00 ,j00+1) c if (izp .eq. 1 .and. i00 .eq. i1) then f10 = f1( 1,j00 ) f11 = f1( 1,j00+1) else f10 = f1(i00+1,j00 ) f11 = f1(i00+1,j00+1) endif c c Calculate the lon,lat of the point i00,j00 rlon00 = rlonn1 + dlon1*float(i00-1) rlat00 = rlatn1 + dlat1*float(j00-1) c In order to correct a potential problem when rlat is close to c either pole, introduce variables arg0 arg1, and arg. In the c subsequent calculations of dx0,dx1,and x, I'll use the fact c that cosine is an even function to avoid if statements. c The problem occurred for situations when the spacing between c points of grid 1 was greater than the spacing between points c of grid 2. When rlat was between the pole and the latitude of c the first row equatorward of the pole in grid 1, dx0 would be c measured at the pole, so it is 0, or sometimes very small. The c variable x was not zero in these cases so the second term in c f2 became very large. When rlat was at the pole, the problem c did not occur because even though dx0 was still small, x=0, so c the second term of f2 did not contribute. This was found when c the NCEP/NCAR reanalysis (2.5 x 2.5) was interpolated to a global c 2.0 x 2.0 grid. argmax=89.5 arg0=abs(rlat00) arg1=abs(rlat00+dlat1) arg=abs(rlat) arg0=min(arg0,argmax) arg1=min(arg1,argmax) arg=min(arg,argmax) c c Calculate the x,y distances between the four f1 points c where x,y = 0,0 at i00,j00 dx0 = dlon1*adtr*cos( dtr*(arg0) ) dx1 = dlon1*adtr*cos( dtr*(arg1) ) dy = dlat1*adtr c c Calculate the x,y coordinates of the current f2 point x = adtr*(rlon-rlon00)*cos(dtr*arg) y = adtr*(rlat-rlat00) c c Calculate the coefficients for the linear interpolation a = f00 b = (f10-f00)/dx0 c = (f01-f00)/dy d = (f11 - f00 - b*dx1 - c*dy)/(dx1*dy) c c Perform interpolation and then go to the next f2 point f2(i,j) = a + b*x + c*y + d*x*y 10 continue c return end