! Module to interpolate values from a giving projection and pressure interpolation ! To be included in a python ! Follow compilation instructions from Makefile ! Content ! LlInterpolateProjection: Subroutine which provides the indices for a given interpolation of a projection ! var2D_IntProj: Subroutine to interpolate a 2D variable ! var3D_IntProj: Subroutine to interpolate a 3D variable ! var4D_IntProj: Subroutine to interpolate a 4D variable ! var5D_IntProj: Subroutine to interpolate a 5D variable MODULE module_ForInterpolate USE module_generic CONTAINS SUBROUTINE CoarselonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, fraclon, fraclat, lonv, latv, per, & Nperx, Npery, ilonlat, mindiffLl) ! Function to search a given value from a coarser version of the data IMPLICIT NONE ! INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: dx, dy REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: ilon, ilat REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in) :: fraclon, fraclat REAL(r_k), INTENT(in) :: lonv, latv, per REAL(r_k), DIMENSION(2), INTENT(in) :: nxlon, nxlat INTEGER, INTENT(in) :: Nperx, Npery INTEGER, DIMENSION(2), INTENT(out) :: ilonlat REAL(r_k), INTENT(out) :: mindiffLl ! Local REAL(r_k), DIMENSION(Nperx,Npery) :: difffraclonlat REAL(r_k) :: mindifffracLl INTEGER, DIMENSION(2) :: ilonlatfrac INTEGER :: ixbeg, ixend, iybeg, iyend INTEGER :: fracx, fracy REAL(r_k) :: fraclonv, fraclatv REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: difflonlat, lon, lat CHARACTER(LEN=50) :: fname ! Variables ! ilon, ilat: original 2D matrices with the longitudes and the latitudes ! lonv, latv: longitude and latitude to find ! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat ! per: fraction of the whole domain (as percentage) ! Nper[x/y]: period (as fraction over 1) of the fractions of the original grid to use to explore ! fraclon, fraclat: longitude and latitude fractional matricies to perform the first guess fname = 'CoarselonlatFind' IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!' PRINT *,' given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )' STOP END IF IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!' PRINT *,' given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )' STOP END IF ! Initializing variables ixbeg = 0 ixend = 0 iybeg = 0 iyend = 0 fracx = int(dx*per) fracy = int(dy*per) ! PRINT *,'fraclon _______' ! PRINT *,fraclon ! PRINT *,'fraclat _______' ! PRINT *,fraclat ! Fraction point difffraclonlat = SQRT((fraclon-lonv)**2. + (fraclat-latv)**2.) mindifffracLl = MINVAL(difffraclonlat) ilonlatfrac = index2DArrayR(difffraclonlat, Nperx, Npery, mindifffracLl) ! PRINT *, 'mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac ! PRINT *, 'frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,', & ! fraclat(ilonlatfrac(1),ilonlatfrac(2)) ! PRINT *, 'values lon, lat:', lonv, latv ! Providing fraction range fraclonv = fraclon(ilonlatfrac(1),ilonlatfrac(2)) fraclatv = fraclat(ilonlatfrac(1),ilonlatfrac(2)) IF (fraclonv >= lonv .AND. fraclatv >= latv) THEN IF (ilonlatfrac(1) > 0) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE ixbeg = 0 ixend = fracx+1 END IF IF (ilonlatfrac(2) > 0) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 0 iyend = fracy+1 END IF ELSE IF (fraclonv < lonv .AND. fraclatv >= latv) THEN IF (ilonlatfrac(1) < Nperx) THEN IF (ilonlatfrac(1) /= 0) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE ixbeg = 0 ixend = fracx+1 END IF ELSE ixbeg = Nperx*fracx ixend = dx+1 END IF IF (ilonlatfrac(2) > 0) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 0 iyend = fracy+1 END IF ELSE IF (fraclonv < lonv .AND. fraclatv < latv) THEN IF (ilonlatfrac(1) < Nperx) THEN IF (ilonlatfrac(1) /= 0) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE ixbeg = 0 ixend = fracx+1 END IF ELSE ixbeg = Nperx*fracx ixend = dx+1 ENDIF IF (ilonlatfrac(2) < Npery) THEN IF (ilonlatfrac(2) /= 0) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 0 iyend = fracy+1 END IF ELSE iybeg = Npery*fracy iyend = dy+1 END IF ELSE IF (fraclonv >= lonv .AND. fraclatv < latv) THEN IF (ilonlatfrac(1) > 0) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE ixbeg = 0 ixend = fracx+1 END IF IF (ilonlatfrac(2) < Npery) THEN IF (ilonlatfrac(2) /= 0) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 0 iyend = fracy+1 END IF ELSE iybeg = Npery*fracy iyend = dy+1 END IF END IF IF (ALLOCATED(lon)) DEALLOCATE(lon) ALLOCATE(lon(ixend-ixbeg+1, iyend-iybeg+1)) IF (ALLOCATED(lat)) DEALLOCATE(lat) ALLOCATE(lat(ixend-ixbeg+1, iyend-iybeg+1)) IF (ALLOCATED(difflonlat)) DEALLOCATE(difflonlat) ALLOCATE(difflonlat(ixend-ixbeg+1, iyend-iybeg+1)) lon = ilon(ixbeg:ixend,iybeg:iyend) lat = ilat(ixbeg:ixend,iybeg:iyend) ! print *,'lon _______' ! print *,lon ! print *,'lat _______' ! print *,lat ! Find point difflonlat = SQRT((lon-lonv)**2. + (lat-latv)**2.) mindiffLl = MINVAL(difflonlat) ilonlat = index2DArrayR(difflonlat, ixend-ixbeg+1, iyend-iybeg+1, mindiffLl) ilonlat(1) = ilonlat(1) + ixbeg ilonlat(2) = ilonlat(2) + iybeg ! PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon ! PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2)) RETURN END SUBROUTINE CoarselonlatFind SUBROUTINE CoarselonlatFindExact(dx, dy, ilon, ilat, nxlon, nxlat, fracx, fracy, fraclon, fraclat, & iv, lonv, latv, per, Nperx, Npery, mindiff, ilonlat, mindiffLl) ! Function to search a given value from a coarser version of the data IMPLICIT NONE ! INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: dx, dy, iv REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: ilon, ilat INTEGER, INTENT(in) :: fracx, fracy REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in) :: fraclon, fraclat REAL(r_k), INTENT(in) :: lonv, latv, per, mindiff REAL(r_k), DIMENSION(2), INTENT(in) :: nxlon, nxlat INTEGER, INTENT(in) :: Nperx, Npery INTEGER, DIMENSION(2), INTENT(out) :: ilonlat REAL(r_k), INTENT(out) :: mindiffLl ! Local INTEGER :: i REAL(r_k), DIMENSION(Nperx,Npery) :: difffraclonlat REAL(r_k) :: mindifffracLl INTEGER, DIMENSION(2) :: ilonlatfrac INTEGER :: ixbeg, ixend, iybeg, iyend REAL(r_k) :: fraclonv, fraclatv REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: difflonlat, lon, lat CHARACTER(LEN=50) :: fname ! Variables ! ilon, ilat: original 2D matrices with the longitudes and the latitudes ! lonv, latv: longitude and latitude to find ! iv: point in the input data ! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat ! per: fraction of the whole domain (as percentage) ! Nper[x/y]: period (as fraction over 1) of the fractions of the original grid to use to explore ! frac[x/y]: Number of grid points for each fraction ! fraclon, fraclat: longitude and latitude fractional matricies to perform the first guess ! mindiff: authorized minimal distance between input and interpolated point ! ilonlat: grid point on the total lon,lat matrix ! mindiffLl: distance between input and interpolated point fname = 'CoarselonlatFindExact' IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!' PRINT *,' given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )' STOP END IF IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!' PRINT *,' given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )' STOP END IF ! Initializing variables ixbeg = 0 ixend = 0 iybeg = 0 iyend = 0 ! Fraction point difffraclonlat = SQRT((fraclon-lonv)**2. + (fraclat-latv)**2.) mindifffracLl = MINVAL(difffraclonlat) ilonlatfrac = index2DArrayR(difffraclonlat, Nperx, Npery, mindifffracLl) ! PRINT *, 'mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac ! PRINT *, 'frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,', & ! fraclat(ilonlatfrac(1),ilonlatfrac(2)) ! PRINT *, 'values lon, lat:', lonv, latv ! Providing fraction range fraclonv = fraclon(ilonlatfrac(1),ilonlatfrac(2)) fraclatv = fraclat(ilonlatfrac(1),ilonlatfrac(2)) IF (fraclonv >= lonv .AND. fraclatv >= latv) THEN PRINT *,'Lluis!',fraclonv, '>=', lonv,'&', fraclatv, '>=', latv IF (ilonlatfrac(1) > 1) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE PRINT *,'Lluis 2' ixbeg = 1 ixend = fracx+1 END IF IF (ilonlatfrac(2) > 1) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 1 iyend = fracy+1 END IF ELSE IF (fraclonv < lonv .AND. fraclatv >= latv) THEN PRINT *,'Lluis!',fraclonv, '<', lonv,'&', fraclatv, '>=', latv IF (ilonlatfrac(1) < Nperx) THEN PRINT *,'Lluis 2' IF (ilonlatfrac(1) /= 1) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE ixbeg = 1 ixend = fracx+1 END IF ELSE ixbeg = Nperx*fracx ixend = dx+1 END IF IF (ilonlatfrac(2) > 1) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 1 iyend = fracy+1 END IF ELSE IF (fraclonv < lonv .AND. fraclatv < latv) THEN PRINT *,'Lluis!',fraclonv, '<', lonv,'&', fraclatv, '<', latv IF (ilonlatfrac(1) < Nperx) THEN IF (ilonlatfrac(1) /= 1) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE ixbeg = 1 ixend = fracx+1 END IF ELSE ixbeg = Nperx*fracx ixend = dx+1 ENDIF IF (ilonlatfrac(2) < Npery) THEN IF (ilonlatfrac(2) /= 1) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 1 iyend = fracy+1 END IF ELSE iybeg = Npery*fracy iyend = dy+1 END IF ELSE IF (fraclonv >= lonv .AND. fraclatv < latv) THEN PRINT *,'Llui!',fraclonv, '>=', lonv,'&', fraclatv, '<', latv IF (ilonlatfrac(1) > 1) THEN ixbeg = (ilonlatfrac(1)-1)*fracx ixend = ilonlatfrac(1)*fracx+1 ELSE ixbeg = 1 ixend = fracx+1 END IF IF (ilonlatfrac(2) < Npery) THEN IF (ilonlatfrac(2) /= 1) THEN iybeg = (ilonlatfrac(2)-1)*fracy iyend = ilonlatfrac(2)*fracy+1 ELSE iybeg = 1 iyend = fracy+1 END IF ELSE iybeg = Npery*fracy iyend = dy+1 END IF END IF IF (ALLOCATED(lon)) DEALLOCATE(lon) ALLOCATE(lon(ixend-ixbeg+1, iyend-iybeg+1)) IF (ALLOCATED(lat)) DEALLOCATE(lat) ALLOCATE(lat(ixend-ixbeg+1, iyend-iybeg+1)) IF (ALLOCATED(difflonlat)) DEALLOCATE(difflonlat) ALLOCATE(difflonlat(ixend-ixbeg+1, iyend-iybeg+1)) lon = ilon(ixbeg:ixend,iybeg:iyend) lat = ilat(ixbeg:ixend,iybeg:iyend) ! print *,'lon _______' ! print *,lon ! print *,'lat _______' ! print *,lat ! Find point difflonlat = SQRT((lon-lonv)**2. + (lat-latv)**2.) mindiffLl = MINVAL(difflonlat) IF (mindiffLl > mindiff) THEN difflonlat = SQRT((lon-lonv)**2. + (lat-latv)**2.) mindiffLl = MINVAL(difflonlat) END IF IF (mindiffLl > mindiff) THEN PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': not equivalent point closer than:',mindiff,' found!!' PRINT *,' at input point iv:', iv,' lon/lat:', lonv,', ',latv,' distance:',mindiffLl PRINT *,' Fraction values _______ (',Nperx,', ',Npery ,')' PRINT *,' fraclon' DO i=1, Nperx PRINT *,' ',fraclon(i,:) END DO PRINT *,' fraclat' DO i=1, Nperx PRINT *,' ',fraclat(i,:) END DO PRINT *,' frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,', & fraclat(ilonlatfrac(1),ilonlatfrac(2)) PRINT *,' mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac PRINT *,' Coarse values _______' PRINT *,' indices. x:', ixbeg, ', ', ixend, ' y:', iybeg, ', ', iyend PRINT *,' lon range:', '(',ilon(ixbeg,iybeg),', ',ilon(ixend,iyend),')' PRINT *,' lat range:', '(',ilat(ixbeg,iybeg),', ',ilat(ixend,iyend),')' PRINT *,' lon', UBOUND(lon) DO i=1, ixend-ixbeg+1 PRINT *,' ',lon(i,:) END DO PRINT *,' lat', UBOUND(lat) DO i=1, ixend-ixbeg+1 PRINT *,' ',lat(i,:) END DO STOP END IF ilonlat = index2DArrayR(difflonlat, ixend-ixbeg+1, iyend-iybeg+1, mindiffLl) ilonlat(1) = ilonlat(1) + ixbeg ilonlat(2) = ilonlat(2) + iybeg ! PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon ! PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2)) RETURN END SUBROUTINE CoarselonlatFindExact SUBROUTINE lonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, lonv, latv, ilonlat, mindiffLl) ! Function to search a given value from a coarser version of the data IMPLICIT NONE ! INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: dx, dy REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: ilon, ilat REAL(r_k), INTENT(in) :: lonv, latv REAL(r_k), DIMENSION(2), INTENT(in) :: nxlon, nxlat INTEGER, DIMENSION(2), INTENT(out) :: ilonlat REAL(r_k), INTENT(out) :: mindiffLl ! Local REAL(r_k), DIMENSION(dx,dy) :: difflonlat CHARACTER(LEN=50) :: fname ! Variables ! ilon, ilat: original 2D matrices with the longitudes and the latitudes ! lonv, latv: longitude and latitude to find ! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat fname = 'lonlatFind' IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!' PRINT *,' given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )' STOP END IF IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!' PRINT *,' given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )' STOP END IF ! Find point difflonlat = SQRT((ilon-lonv)**2. + (ilat-latv)**2.) mindiffLl = MINVAL(difflonlat) ilonlat = index2DArrayR(difflonlat, dx, dy, mindiffLl) ! PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon ! PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2)) RETURN END SUBROUTINE lonlatFind SUBROUTINE CoarseInterpolate(projlon, projlat, lonvs, latvs, percen, mindiff, inpt, ilonlat, & mindiffLl, dimx, dimy, Ninpts) ! Subroutine which finds the closest grid point within a projection throughout a first guest ! approche from percentages of the whole domain IMPLICIT NONE ! INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: dimx, dimy REAL(r_k), DIMENSION(dimx,dimy), INTENT(in) :: projlon, projlat INTEGER, INTENT(in) :: Ninpts REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: inpt, lonvs, latvs REAL(r_k), INTENT(in) :: mindiff, percen INTEGER, DIMENSION(Ninpts,2), INTENT(out) :: ilonlat REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: mindiffLl ! Local INTEGER :: iv,i,j INTEGER :: ierr INTEGER :: Ninpts1 REAL(r_k), DIMENSION(2) :: extremelon, extremelat REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: fractionlon, fractionlat INTEGER :: dfracdx, dfracdy, fracdx, fracdy CHARACTER(LEN=50) :: fname !!!!!!! Variables ! dimx, dimy: dimension length of the target interpolation ! proj[lon/lat]: longitudes and latitudes of the target interpolation ! Ninpts: number of points to interpolate ! [lon/lat]vs: longitudes and latitudes of the points to interpolate ! mindiff: minimal accepted distance to the target point ! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues ! inpt: whether the point has already been localized (1) or not (0) ! ilonlat: Longitude and Latitude of the input points ! mindiffLl: minimum difference between target and source longitude/latitude (in degrees) fname = 'CoarseInterpolate' Ninpts1 = Ninpts/100 extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /) extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /) PRINT *,' ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen dfracdx = INT(1./percen+1) dfracdy = INT(1./percen+1) fracdx = INT(dimx*percen) fracdy = INT(dimy*percen) PRINT *,' ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon) ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr) IF (ierr /= 0) THEN PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) //": problem allocating 'fractionlon' !!" STOP END IF IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat) ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr) IF (ierr /= 0) THEN PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) //": problem allocating 'fractionlat' !!" STOP END IF DO i=1,dfracdx DO j=1,dfracdy fractionlon(i,j) = projlon(fracdx*(i-1)+1,fracdy*(j-1)+1) fractionlat(i,j) = projlat(fracdx*(i-1)+1,fracdy*(j-1)+1) ! PRINT *,'i,j:',i,', ',j,' frac ij:',fracdx*(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),& ! ', ',fractionlat(i,j) END DO END DO ! PRINT *,' ' // TRIM(fname) // ' fractions of:' ! PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')' ! DO i=1,dfracdx ! PRINT *,fractionlon(i,:) ! END DO ! PRINT *,' lat_______' ! DO i=1,dfracdx ! PRINT *,fractionlat(i,:) ! END DO DO iv=1,Ninpts IF (inpt(iv) == 0) THEN CALL CoarselonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, fractionlon, & fractionlat, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, ilonlat(iv,:), mindiffLl(iv)) IF ((mindiffLl(iv) <= mindiff) .AND. .NOT.(ilonlat(iv,1) >= 0 .AND. ilonlat(iv,1) >= 0)) THEN PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), & ' not relocated !!' PRINT *,' mindiffl:', mindiffLl(iv), ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2) STOP END IF END IF END DO END SUBROUTINE CoarseInterpolate SUBROUTINE CoarseInterpolateExact(projlon, projlat, lonvs, latvs, percen, mindiff, ivar, newvar, & newvarin, newvarinpt, newvarindiff, dimx, dimy, Ninpts) ! Subroutine which finds the closest grid point within a projection throughout a first guest ! and then whole domain approche from percentages of the whole domain IMPLICIT NONE ! INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: dimx, dimy REAL(r_k), DIMENSION(dimx,dimy), INTENT(in) :: projlon, projlat INTEGER, INTENT(in) :: Ninpts REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: ivar, lonvs, latvs REAL(r_k), INTENT(in) :: mindiff, percen REAL(r_k), DIMENSION(dimx,dimy), INTENT(out) :: newvar INTEGER, DIMENSION(dimx,dimy), INTENT(out) :: newvarin INTEGER, DIMENSION(Ninpts), INTENT(out) :: newvarinpt REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: newvarindiff ! Local INTEGER :: iv,i,j INTEGER :: ierr INTEGER, DIMENSION(2) :: ilonlat REAL(r_k) :: mindiffLl INTEGER :: Ninpts1 REAL(r_k), DIMENSION(2) :: extremelon, extremelat REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: fractionlon, fractionlat INTEGER :: dfracdx, dfracdy, fracdx, fracdy CHARACTER(LEN=50) :: fname !!!!!!! Variables ! dimx, dimy: dimension length of the target interpolation ! proj[lon/lat]: longitudes and latitudes of the target interpolation ! Ninpts: number of points to interpolate ! [lon/lat]vs: longitudes and latitudes of the points to interpolate ! mindiff: minimal accepted distance to the target point ! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues ! ivar: values to localize in the target projection ! newvar: localisation of the [lon/lat]vs point in the target projection ! newvarin: number of point from the input data ! newvarinpt: integer value indicating if the value has been already located (0: no, 1: yes) ! newvarindiff: distance of point from the input data to the closest target point ! ncid: netCDF output file id fname = 'CoarseInterpolateExact' Ninpts1 = Ninpts/100 extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /) extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /) PRINT *,' ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen dfracdx = INT(1./percen+1) dfracdy = INT(1./percen+1) fracdx = INT(dimx*percen) fracdy = INT(dimy*percen) PRINT *,' ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon) ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr) IF (ierr /= 0) THEN PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) //": problem allocating 'fractionlon' !!" STOP END IF IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat) ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr) IF (ierr /= 0) THEN PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) //": problem allocating 'fractionlat' !!" STOP END IF DO i=1,dfracdx DO j=1,dfracdy fractionlon(i,j) = projlon(fracdx*(i-1)+1,fracdy*(j-1)+1) fractionlat(i,j) = projlat(fracdx*(i-1)+1,fracdy*(j-1)+1) ! PRINT *,'i,j:',i,', ',j,' frac ij:',fracdx*(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),& ! ', ',fractionlat(i,j) END DO END DO ! PRINT *,' ' // TRIM(fname) // ' fractions of:' ! PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')' ! DO i=1,dfracdx ! PRINT *,fractionlon(i,:) ! END DO ! PRINT *,' lat_______' ! DO i=1,dfracdx ! PRINT *,fractionlat(i,:) ! END DO DO iv=1,Ninpts IF (newvarinpt(iv) == 0) THEN CALL CoarselonlatFindExact(dimx, dimy, projlon, projlat, extremelon, extremelat, fracdx, fracdy,& fractionlon, fractionlat, iv, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, mindiff, & ilonlat, mindiffLl) IF (mindiffLl >= mindiff) THEN ! percendone(iv,Ninpts,0.5,'done:') IF (ilonlat(1) >= 0 .AND. ilonlat(1) >= 0) THEN newvar(ilonlat(1),ilonlat(2)) = ivar(iv) newvarin(ilonlat(1),ilonlat(2)) = iv newvarinpt(iv) = 1 newvarindiff(iv) = mindiffLl ! PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv), & ! ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:', & ! newvarindiff(iv), ' point:',ilonlat ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), & ' not relocated !!' PRINT *,' mindiffl:', mindiffLl, ' ilon:', ilonlat(1), ' ilat:', ilonlat(2) STOP END IF ! IF (MOD(iv,Ninpts1) == 0) newnc.sync() ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv), & ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ', & mindiff, ' !!' PRINT *,' found minimum difference:', mindiffLl STOP END IF END IF END DO END SUBROUTINE CoarseInterpolateExact SUBROUTINE LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, & pdimx, pdimy) ! Subroutine which provides the indices for a given interpolation of a projection IMPLICIT NONE ! INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv CHARACTER(LEN=50), INTENT(in) :: intkind REAL(r_k), DIMENSION(3,16,pdimx,pdimy), INTENT(out) :: outLlw ! Local INTEGER :: i,j,iv, ix, iy REAL(r_k) :: mindiffLl, dist REAL(r_k) :: inclon, inclat, maxdiffprojlonlat, & maxdiffinlonlat REAL(r_k), DIMENSION(idimx,idimy) :: difflonlat REAL(r_k), DIMENSION(idimx,idimy) :: idifflon, idifflat REAL(r_k), DIMENSION(pdimx,pdimy) :: difflon, difflat REAL(r_k), DIMENSION(2) :: extremelon, extremelat, ipos INTEGER, DIMENSION(2) :: iLl CHARACTER(LEN=50) :: fname !!!!!!! Variables ! idimx, idimy: dimension length of the input projection ! pdimx, pdimy: dimension length of the target projection ! in[lon/lat]: longitudes and latitudes of the target projection ! proj[lon/lat]: longitudes and latitudes of the target projection ! intkind: kind of interpolation ! 'npp': nearest neighbourgh ! 'dis': weighted distance, grid-output for SW, NW, NE, SE ! outLlw: output interpolation result ! for point pi,pj; up to 16 different values of ! 1st: i-index in input projection ! 2nd: j-index in input projection ! 3rd: weight for i-index, j-index to use for ponderation (0<1.) fname = 'LlInterpolateProjection' extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /) extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /) ! Maximum distance between grid points in input projection idifflon = 0. idifflat = 0. idifflon(1:idimx-1,:) = inlonv(2:idimx,:)-inlonv(1:idimx-1,:) idifflat(:,1:idimy-1) = inlatv(:,2:idimy)-inlatv(:,1:idimy-1) maxdiffinlonlat = MAXVAL(SQRT(idifflon**2. + idifflat**2.)) ! Maximum distance between grid points in target projection difflon = 0. difflat = 0. difflon(1:pdimx-1,:) = projlon(2:pdimx,:)-projlon(1:pdimx-1,:) difflat(:,1:pdimy-1) = projlat(:,2:pdimy)-projlat(:,1:pdimy-1) maxdiffprojlonlat = MAXVAL(SQRT(difflon**2. + difflat**2.)) IF (maxdiffinlonlat > maxdiffprojlonlat) THEN PRINT *,TRIM(warnmsg) PRINT *,' ' //TRIM(fname)// '; input resolution: ', maxdiffinlonlat, ' is coarser than target:', & maxdiffprojlonlat, ' !!' END IF ! Using case outside loop to be more efficient SELECT CASE(TRIM(intkind)) CASE ('dis') inclon = inlonv(2,1) - inlonv(1,1) inclat = inlatv(1,2) - inlatv(1,1) DO i=1, pdimx DO j=1, pdimy ! Find point difflonlat = SQRT((projlon(i,j)-inlonv)**2. + (projlat(i,j)-inlatv)**2.) mindiffLl = MINVAL(difflonlat) IF ( (mindiffLl > maxdiffprojlonlat) .AND. (mindiffLl > maxdiffinlonlat)) THEN outLlw(3,:,i,j) = 0. outLlw(3,:,i,j) = -1. ELSE ! Getting the four surrounding values iLl = index2DArrayR(difflonlat, idimx, idimy, mindiffLl) IF ( (projlon(i,j) < inlonv(iLl(1),iLl(2)) .AND. inclon > 0.) .OR. & (projlon(i,j) > inlonv(iLl(1),iLl(2)) .AND. inclon < 0.) ) THEN outLlw(1,1,i,j) = MAX(iLl(1)-1,1) outLlw(1,2,i,j) = MAX(iLl(1)-1,1) outLlw(1,3,i,j) = MIN(iLl(1),idimx) outLlw(1,4,i,j) = MIN(iLl(1),idimx) ELSE outLlw(1,1,i,j) = MAX(iLl(1),1) outLlw(1,2,i,j) = MAX(iLl(1),1) outLlw(1,3,i,j) = MIN(iLl(1)+1,idimx) outLlw(1,4,i,j) = MIN(iLl(1)+1,idimx) END IF IF ( (projlat(i,j) < inlatv(iLl(2),iLl(2)) .AND. inclat > 0.) .OR. & (projlat(i,j) > inlatv(iLl(2),iLl(2)) .AND. inclat < 0.) ) THEN outLlw(2,1,i,j) = MAX(iLl(2)-1,1) outLlw(2,2,i,j) = MIN(iLl(2),idimy) outLlw(2,3,i,j) = MIN(iLl(2),idimy) outLlw(2,4,i,j) = MAX(iLl(2)-1,1) ELSE outLlw(2,1,i,j) = MAX(iLl(2),1) outLlw(2,2,i,j) = MIN(iLl(2)+1,idimy) outLlw(2,3,i,j) = MIN(iLl(2)+1,idimy) outLlw(2,4,i,j) = MAX(iLl(2),1) END IF ! Computing distances !Keep the print for future checks? ! IF (MOD(i+j,200) == 0) THEN ! PRINT *,'center point:',i,j,'=', iLl,':',projlon(i,j),projlat(i,j),' incs',inclon,' ,',inclat ! PRINT *,'SW:', outLlw(1,1,i,j), outLlw(2,1,i,j),':',inlonv(outLlw(1,1,i,j), outLlw(2,1,i,j)),& ! inlatv(outLlw(1,1,i,j), outLlw(2,1,i,j)) ! PRINT *,'NW:', outLlw(1,2,i,j), outLlw(2,2,i,j),':',inlonv(outLlw(1,2,i,j), outLlw(2,2,i,j)),& ! inlatv(outLlw(1,2,i,j), outLlw(2,2,i,j)) ! PRINT *,'NE:', outLlw(1,3,i,j), outLlw(2,3,i,j),':',inlonv(outLlw(1,3,i,j), outLlw(2,3,i,j)),& ! inlatv(outLlw(1,3,i,j), outLlw(2,3,i,j)) ! PRINT *,'SE:', outLlw(1,4,i,j), outLlw(2,4,i,j),':',inlonv(outLlw(1,4,i,j), outLlw(2,4,i,j)),& ! inlatv(outLlw(1,4,i,j), outLlw(2,4,i,j)) ! END IF DO iv=1,4 ix = INT(outLlw(1,iv,i,j)) iy = INT(outLlw(2,iv,i,j)) dist = SQRT( (projlon(i,j)-inlonv(ix,iy))**2. + (projlat(i,j)-inlatv(ix,iy))**2. ) IF ( dist /= 0.) THEN outLlw(3,iv,i,j) = 1./dist ELSE outLlw(3,iv,i,j) = 1. END IF ! IF (i+j == 2) PRINT *,'iv:',iv,'dist:',dist,'weight:',outLlw(3,iv,i,j) END DO ! Normalizing outLlw(3,:,i,j) = outLlw(3,:,i,j)/(SUM(outLlw(3,:,i,j))) ! IF (i+j == 2) PRINT *,'Normalized weights:',outLlw(3,:,i,j),':',SUM(outLlw(3,:,i,j)) END IF END DO END DO CASE ('npp') DO i=1, pdimx DO j=1, pdimy ! Find point difflonlat = SQRT((projlon(i,j)-inlonv)**2. + (projlat(i,j)-inlatv)**2.) mindiffLl = MINVAL(difflonlat) ipos = index2DArrayR(difflonlat, idimx, idimy, mindiffLl) outLlw(1,1,i,j) = REAL(ipos(1)) outLlw(2,1,i,j) = REAL(ipos(2)) ! We do not want that values larger that the maximum distance between target grid points ! PRINT *,i,j,':',mindiffLl,'maxdiffLl:',maxdiffprojlonlat IF ((mindiffLl .gt. maxdiffprojlonlat) .AND. (mindiffLl > maxdiffinlonlat)) THEN ! PRINT *,' ' // TRIM(fname) // ': reprojected minimum distance to nearest grid point:', & ! mindiffLl, ' larger than the maximum distance between grid points in target projection!!' outLlw(3,1,i,j) = -1. ELSE outLlw(3,1,i,j) = 1. END IF ix = INT(outLlw(1,1,i,j)) iy = INT(outLlw(2,1,i,j)) END DO END DO END SELECT END SUBROUTINE LlInterpolateProjection SUBROUTINE var2D_IntProj(var2Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx, & idimy, pdimx, pdimy) ! Subroutine to interpolate a 2D variable IMPLICIT NONE INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv CHARACTER(LEN=50), INTENT(in) :: intkind REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: var2Din INTEGER, DIMENSION(idimx,idimy), INTENT(in) :: mask REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(out) :: varout ! Local INTEGER :: i, j, k, iv, ix, iy REAL(r_k) :: w REAL(r_k), DIMENSION(3,16,pdimx,pdimy) :: outLlw CHARACTER(LEN=50) :: fname !!!!!!! Variables ! idimx, idimy: dimension length of the input projection ! pdimx, pdimy: dimension length of the target projection ! in[lon/lat]: longitudes and latitudes of the target projection ! proj[lon/lat]: longitudes and latitudes of the target projection ! intkind: kind of interpolation ! 'npp': nearest neighbourgh ! 'dis': weighted distance, grid-output for SW, NW, NE, SE ! outLlw: output interpolation result ! for point pi,pj; up to 16 different values of ! 1st: i-index in input projection ! 2nd: j-index in input projection ! 3rd: weight for i-index, j-index to use for ponderation (0<1.) ! var2Din: 2D variable to interpolate ! mask: mask of the intpu values (1: good, 0: none) ! varout: variable interpolated on the target projection fname = 'var2D_IntProj' CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,& pdimy) SELECT CASE (intkind) CASE('dis') DO i=1, pdimx DO j=1, pdimy IF (outLlw(3,1,i,j) == -1.) THEN varout(i,j) = fillVal64 ELSE varout(i,j) = 0. DO iv=1, 4 ix = INT(outLlw(1,iv,i,j)) iy = INT(outLlw(2,iv,i,j)) IF (mask(ix,iy) == 1) THEN w = outLlw(3,iv,i,j) varout(i,j) = varout(i,j) + w*var2Din(ix,iy) END IF END DO END IF END DO END DO CASE('npp') DO i=1, pdimx DO j=1, pdimy ix = INT(outLlw(1,1,i,j)) iy = INT(outLlw(2,1,i,j)) IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy) == 0) ) THEN varout(i,j) = fillVal64 ELSE varout(i,j) = var2Din(ix,iy)*outLlw(3,1,i,j) END IF END DO END DO END SELECT END SUBROUTINE var2D_IntProj SUBROUTINE var3D_IntProj(var3Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx, & idimy, pdimx, pdimy, d3) ! Subroutine to interpolate a 3D variable IMPLICIT NONE ! INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy, d3 REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv CHARACTER(LEN=50), INTENT(in) :: intkind REAL(r_k), DIMENSION(idimx,idimy,d3), INTENT(in) :: var3Din INTEGER, DIMENSION(idimx,idimy,d3), INTENT(in) :: mask REAL(r_k), DIMENSION(pdimx,pdimy,d3), INTENT(out) :: varout ! Local INTEGER :: i, j, k, iv, ix, iy REAL(r_k) :: w REAL(r_k), DIMENSION(3,16,pdimx,pdimy) :: outLlw CHARACTER(LEN=50) :: fname !!!!!!! Variables ! idimx, idimy: dimension length of the input projection ! pdimx, pdimy: dimension length of the target projection ! in[lon/lat]: longitudes and latitudes of the target projection ! proj[lon/lat]: longitudes and latitudes of the target projection ! intkind: kind of interpolation ! 'npp': nearest neighbourgh ! 'dis': weighted distance, grid-output for SW, NW, NE, SE ! outLlw: output interpolation result ! for point pi,pj; up to 16 different values of ! 1st: i-index in input projection ! 2nd: j-index in input projection ! 3rd: weight for i-index, j-index to use for ponderation (0<1.) ! var3Din: 3D variable to interpolate ! mask: mask of the intpu values (1: good, 0: none) ! varout: variable interpolated on the target projection fname = 'var3D_IntProj' CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,& pdimy) SELECT CASE (intkind) CASE('dis') DO i=1, pdimx DO j=1, pdimy IF (ALL(outLlw(3,:,i,j) == -1.)) THEN varout(i,j,:) = fillVal64 ELSE DO k=1, d3 varout(i,j,k) = 0. DO iv=1, 4 ix = INT(outLlw(1,iv,i,j)) iy = INT(outLlw(2,iv,i,j)) IF (mask(ix,iy,k) == 1) THEN w = outLlw(3,iv,i,j) varout(i,j,k) = varout(i,j,k) + w*var3Din(ix,iy,k) END IF END DO END DO END IF END DO END DO CASE('npp') DO i=1, pdimx DO j=1, pdimy ix = INT(outLlw(1,1,i,j)) iy = INT(outLlw(2,1,i,j)) IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy,1) == 0) ) THEN varout(i,j,:) = fillVal64 ELSE DO k=1, d3 varout(i,j,k) = var3Din(ix,iy,k)*outLlw(3,1,i,j) END DO END IF END DO END DO END SELECT END SUBROUTINE var3D_IntProj SUBROUTINE var4D_IntProj(var4Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx, & idimy, pdimx, pdimy, d3, d4) ! Subroutine to interpolate a 4D variable IMPLICIT NONE ! INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy, d3, d4 REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv CHARACTER(LEN=50), INTENT(in) :: intkind REAL(r_k), DIMENSION(idimx,idimy,d3,d4), INTENT(in) :: var4Din INTEGER, DIMENSION(idimx,idimy,d3,d4), INTENT(in) :: mask REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4), INTENT(out) :: varout ! Local INTEGER :: i, j, k, l, iv, ix, iy REAL(r_k) :: w REAL(r_k), DIMENSION(3,16,pdimx,pdimy) :: outLlw CHARACTER(LEN=50) :: fname !!!!!!! Variables ! idimx, idimy: dimension length of the input projection ! pdimx, pdimy: dimension length of the target projection ! in[lon/lat]: longitudes and latitudes of the target projection ! proj[lon/lat]: longitudes and latitudes of the target projection ! intkind: kind of interpolation ! 'npp': nearest neighbourgh ! 'dis': weighted distance, grid-output for SW, NW, NE, SE ! outLlw: output interpolation result ! for point pi,pj; up to 16 different values of ! 1st: i-index in input projection ! 2nd: j-index in input projection ! 3rd: weight for i-index, j-index to use for ponderation (0<1.) ! var4Din: 4D variable to interpolate ! mask: mask of the intpu values (1: good, 0: none) ! varout: variable interpolated on the target projection fname = 'var4D_IntProj' CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,& pdimy) SELECT CASE (intkind) CASE('dis') DO i=1, pdimx DO j=1, pdimy IF (ALL(outLlw(3,:,i,j) == -1.)) THEN varout(i,j,:,:) = fillVal64 ELSE DO k=1, d3 DO l=1, d4 varout(i,j,k,l) = 0. DO iv=1, 4 ix = INT(outLlw(1,iv,i,j)) iy = INT(outLlw(2,iv,i,j)) IF (mask(ix,iy,k,l) == 1) THEN w = outLlw(3,iv,i,j) varout(i,j,k,l) = varout(i,j,k,l) + w*var4Din(ix,iy,k,l) END IF END DO END DO END DO END IF END DO END DO CASE('npp') DO i=1, pdimx DO j=1, pdimy ix = INT(outLlw(1,1,i,j)) iy = INT(outLlw(2,1,i,j)) IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy,1,1) == 0) ) THEN varout(i,j,:,:) = fillVal64 ELSE DO k=1, d3 DO l=1, d4 varout(i,j,k,l) = var4Din(ix,iy,k,l)*outLlw(3,1,i,j) END DO END DO END IF END DO END DO END SELECT END SUBROUTINE var4D_IntProj SUBROUTINE var5D_IntProj(var5Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx, & idimy, pdimx, pdimy, d3, d4, d5) ! Subroutine to interpolate a 5D variable IMPLICIT NONE ! INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy, d3, d4, d5 REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv CHARACTER(LEN=50), INTENT(in) :: intkind REAL(r_k), DIMENSION(idimx,idimy,d3,d4,d5), INTENT(in) :: var5Din INTEGER, DIMENSION(idimx,idimy,d3,d4,d5), INTENT(in) :: mask REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4,d5), INTENT(out) :: varout ! Local INTEGER :: i, j, k, l, m, iv, ix, iy REAL(r_k) :: w REAL(r_k), DIMENSION(3,16,pdimx,pdimy) :: outLlw CHARACTER(LEN=50) :: fname !!!!!!! Variables ! idimx, idimy: dimension length of the input projection ! pdimx, pdimy: dimension length of the target projection ! in[lon/lat]: longitudes and latitudes of the target projection ! proj[lon/lat]: longitudes and latitudes of the target projection ! intkind: kind of interpolation ! 'npp': nearest neighbourgh ! 'dis': weighted distance, grid-output for SW, NW, NE, SE ! outLlw: output interpolation result ! for point pi,pj; up to 16 different values of ! 1st: i-index in input projection ! 2nd: j-index in input projection ! 3rd: weight for i-index, j-index to use for ponderation (0<1.) ! var5Din: 5D variable to interpolate ! mask: mask of the intpu values (1: good, 0: none) ! varout: variable interpolated on the target projection fname = 'var5D_IntProj' CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,& pdimy) SELECT CASE (intkind) CASE('dis') DO i=1, pdimx DO j=1, pdimy IF (ALL(outLlw(3,:,i,j) == -1.)) THEN varout(i,j,:,:,:) = fillVal64 ELSE DO k=1, d3 DO l=1, d4 DO m=1, d5 varout(i,j,k,l,m) = 0. DO iv=1, 4 ix = INT(outLlw(1,iv,i,j)) iy = INT(outLlw(2,iv,i,j)) IF (mask(ix,iy,k,l,m) == 1) THEN w = outLlw(3,iv,i,j) varout(i,j,k,l,m) = varout(i,j,k,l,m) + w*var5Din(ix,iy,k,l,m) END IF END DO END DO END DO END DO END IF END DO END DO CASE('npp') DO i=1, pdimx DO j=1, pdimy ix = INT(outLlw(1,1,i,j)) iy = INT(outLlw(2,1,i,j)) IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy,1,1,1) == 0) ) THEN varout(i,j,:,:,:) = fillVal64 ELSE DO k=1, d3 DO l=1, d4 DO m=1, d5 varout(i,j,k,l,m) = var5Din(ix,iy,k,l,m)*outLlw(3,1,i,j) END DO END DO END DO END IF END DO END DO END SELECT END SUBROUTINE var5D_IntProj SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy, & Ninpts) ! Subroutine which finds the closest grid point within a projection IMPLICIT NONE ! INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: dimx, dimy REAL(r_k), DIMENSION(dimx,dimy), INTENT(in) :: projlon, projlat INTEGER, INTENT(in) :: Ninpts REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: lonvs, latvs REAL(r_k), INTENT(in) :: mindiff INTEGER, DIMENSION(Ninpts), INTENT(inout) :: inpt REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: diffs INTEGER, DIMENSION(Ninpts,2), INTENT(out) :: ilonlat ! Local INTEGER :: iv REAL(r_k) :: mindiffLl INTEGER :: Ninpts1 REAL(r_k), DIMENSION(dimx,dimy) :: difflonlat REAL(r_k), DIMENSION(2) :: extremelon, extremelat CHARACTER(LEN=50) :: fname !!!!!!! Variables ! dimx, dimy: dimension length of the target interpolation ! proj[lon/lat]: longitudes and latitudes of the target interpolation ! Ninpts: number of points to interpolate ! [lon/lat]vs: longitudes and latitudes of the points to interpolate ! mindiff: minimal accepted distance to the target point ! inpt: whether the point has already been localized ! diffs: distance of point from the input data to the closest target point ! ilonlat: longitude and latitude of the point ! ncid: netCDF output file id fname = 'Interpolate' Ninpts1 = Ninpts/100 extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /) extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /) DO iv=1,Ninpts IF (inpt(iv) <= 0) THEN ! Not using the subroutine, not efficient! ! CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv), & ! ilonlat, mindiffLl) IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!' PRINT *,' given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )' STOP END IF IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!' PRINT *,' given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )' STOP END IF ! Find point difflonlat = SQRT((projlon-lonvs(iv))**2. + (projlat-latvs(iv))**2.) mindiffLl = MINVAL(difflonlat) ilonlat(iv,:) = index2DArrayR(difflonlat, dimx, dimy, mindiffLl) IF (mindiffLl <= mindiff) THEN ! percendone(iv,Ninpts,0.5,'done:') IF (ilonlat(iv,1) >= 0 .AND. ilonlat(iv,2) >= 0) THEN diffs(iv) = mindiffLl inpt(iv) = 1 ! PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv), & ! ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:', & ! newvarindiff(iv), ' point:',ilonlat ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), & ' not relocated !!' PRINT *,' mindiffl:', mindiffLl, ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2) STOP END IF ! IF (MOD(iv,Ninpts1) == 0) newnc.sync() ELSE ! Because doing boxes and Goode is not conitnuos, we should jump this error message PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv), & ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ', & mindiff, ' !!' PRINT *,' found minimum difference:', mindiffLl STOP END IF END IF END DO END SUBROUTINE Interpolate SUBROUTINE Interpolate1DLl(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy, & Ninpts) ! Subroutine which finds the closest grid point within a projection with 1D longitudes and latitudes IMPLICIT NONE ! INTEGER, PARAMETER :: r_k = KIND(1.d0) INTEGER, INTENT(in) :: dimx, dimy REAL(r_k), DIMENSION(dimx), INTENT(in) :: projlon REAL(r_k), DIMENSION(dimy), INTENT(in) :: projlat INTEGER, INTENT(in) :: Ninpts REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: lonvs, latvs REAL(r_k), INTENT(in) :: mindiff INTEGER, DIMENSION(Ninpts), INTENT(inout) :: inpt REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: diffs INTEGER, DIMENSION(Ninpts,2), INTENT(out) :: ilonlat ! Local INTEGER :: iv REAL(r_k) :: mindifflo, mindiffLa, mindiffLl INTEGER :: Ninpts1 REAL(r_k), DIMENSION(dimx) :: difflon REAL(r_k), DIMENSION(dimy) :: difflat REAL(r_k), DIMENSION(2) :: extremelon, extremelat CHARACTER(LEN=50) :: fname !!!!!!! Variables ! dimx, dimy: dimension length of the target interpolation ! proj[lon/lat]: longitudes and latitudes of the target interpolation ! Ninpts: number of points to interpolate ! [lon/lat]vs: longitudes and latitudes of the points to interpolate ! mindiff: minimal accepted distance to the target point ! inpt: whether the point has already been localized ! diffs: distance of point from the input data to the closest target point ! ilonlat: longitude and latitude of the point ! ncid: netCDF output file id fname = 'Interpolate1DLl' Ninpts1 = Ninpts/100 extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /) extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /) DO iv=1,Ninpts IF (inpt(iv) <= 0) THEN ! Not using the subroutine, not efficient! ! CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv), & ! ilonlat, mindiffLl) IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!' PRINT *,' given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )' STOP END IF IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN PRINT *, TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!' PRINT *,' given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )' STOP END IF ! Find point difflon = SQRT((projlon-lonvs(iv))**2.) difflat = SQRT((projlat-latvs(iv))**2.) mindifflo = MINVAL(difflon) mindiffLa = MINVAL(difflat) mindifflL = SQRT(mindifflo*mindifflo + mindiffLa*mindiffLa) ilonlat(iv,1) = index1DArrayR(difflon, dimx, mindifflo) ilonlat(iv,2) = index1DArrayR(difflat, dimy, mindiffLa) ! PRINT *,' Lluis: iv',iv,' lonvs:', lonvs(iv),' latvs:',latvs(iv) ! PRINT *,' Lluis: mindifflo:', mindifflo,' ilonlat(1):',ilonlat(iv,1) ! PRINT *,' Lluis: mindiffLa:', mindiffLa,' ilonlat(2):',ilonlat(iv,2) IF (mindiffLl <= mindiff) THEN ! percendone(iv,Ninpts,0.5,'done:') IF (ilonlat(iv,1) >= 1 .AND. ilonlat(iv,2) >= 1) THEN diffs(iv) = mindiffLl inpt(iv) = 1 ! PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv), & ! ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:', & ! newvarindiff(iv), ' point:',ilonlat ELSE PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), & ' not relocated !!' PRINT *,' mindiffl:', mindiffLl, ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2) STOP END IF ! IF (MOD(iv,Ninpts1) == 0) newnc.sync() ELSE ! Because doing boxes and Goode is not conitnuos, we should jump this error message PRINT *,TRIM(ErrWarnMsg('err')) PRINT *,' ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv), & ' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ', & mindiff, ' !!' PRINT *,' found minimum difference:', mindiffLl STOP END IF END IF END DO END SUBROUTINE Interpolate1DLl SUBROUTINE interp (data_in, pres_field, interp_levels, psfc, ter, tk, qv, LINLOG, extrapolate, & GEOPT, MISSING, data_out, ix, iy, iz, it, num_metgrid_levels) ! Interpolation subroutine from the p_interp.F90 NCAR program ! Program to read wrfout data and interpolate to pressure levels ! The program reads namelist.pinterp ! November 2007 - Cindy Bruyere ! INTEGER, INTENT(IN) :: ix, iy, iz, it INTEGER, INTENT(IN) :: num_metgrid_levels, LINLOG REAL, DIMENSION(ix,iy,iz,it), INTENT(IN) :: data_in, pres_field, tk, qv REAL, DIMENSION(ix,iy,it), INTENT(IN) :: psfc REAL, DIMENSION(ix,iy), INTENT(IN) :: ter REAL, DIMENSION(num_metgrid_levels), INTENT(IN) :: interp_levels INTEGER, INTENT(IN) :: extrapolate REAL, INTENT(IN) :: MISSING LOGICAL, INTENT(IN) :: GEOPT REAL, DIMENSION(ix,iy,num_metgrid_levels,it), & INTENT(OUT) :: data_out ! Local INTEGER :: i, j, itt, k, kk, kin REAL, DIMENSION(num_metgrid_levels) :: data_out1D REAL, DIMENSION(iz) :: data_in1D, pres_field1D REAL, DIMENSION(ix, iy, num_metgrid_levels, it) :: N REAL :: sumA, sumN, AVE_geopt !!!!!!! Variables ! data_out: interpolated field ! data_in: field to interpolate ! pres_field: pressure field [Pa] ! interp_levels: pressure levels to interpolate [hPa] ! psfc: surface pressure [Pa] ! ter: terrein height [m] ! tk: temperature [K] ! qv: mositure mizing ratio [kg/kg] ! i[x/y/z/t]: size of the matrices ! num_metgrid_levels: number of pressure values to interpolate ! LINLOG: if abs(linlog)=1 use linear interp in pressure ! if abs(linlog)=2 linear interp in ln(pressure) ! extrapolate: whether to set to missing value below/above model ground and top (0), or extrapolate (1) ! GEOPT: Wether the file is the geopotential file or not ! MISSING: Missing value N = 1.0 expon=287.04*.0065/9.81 do itt = 1, it do j = 1, iy do i = 1, ix data_in1D(:) = data_in(i,j,:,itt) pres_field1D(:) = pres_field(i,j,:,itt) CALL int1D (data_out1D, data_in1D, pres_field1D, interp_levels, iz, num_metgrid_levels, LINLOG, MISSING) data_out(i,j,:,itt) = data_out1D(:) end do end do end do ! Fill in missing values IF ( extrapolate == 0 ) RETURN !! no extrapolation - we are out of here ! First find where about 400 hPa is located kk = 0 find_kk : do k = 1, num_metgrid_levels kk = k if ( interp_levels(k) <= 40000. ) exit find_kk end do find_kk IF ( GEOPT ) THEN !! geopt is treated different below ground do itt = 1, it do k = 1, kk do j = 1, iy do i = 1, ix IF ( data_out(i,j,k,itt) == MISSING .AND. interp_levels(k) < psfc(i,j,itt) ) THEN ! We are below the first model level, but above the ground data_out(i,j,k,itt) = ((interp_levels(k) - pres_field(i,j,1,itt))*ter(i,j)*9.81 + & (psfc(i,j,itt) - interp_levels(k))*data_in(i,j,1,itt) ) / & (psfc(i,j,itt) - pres_field(i,j,1,itt)) ELSEIF ( data_out(i,j,k,itt) == MISSING ) THEN ! We are below both the ground and the lowest data level. ! First, find the model level that is closest to a "target" pressure ! level, where the "target" pressure is delta-p less that the local ! value of a horizontally smoothed surface pressure field. We use ! delta-p = 150 hPa here. A standard lapse rate temperature profile ! passing through the temperature at this model level will be used ! to define the temperature profile below ground. This is similar ! to the Benjamin and Miller (1990) method, except that for ! simplicity, they used 700 hPa everywhere for the "target" pressure. ! Code similar to what is implemented in RIP4 ptarget = (psfc(i,j,itt)*.01) - 150. dpmin=1.e4 kupper = 0 loop_kIN : do kin=iz,1,-1 kupper = kin dp=abs( (pres_field(i,j,kin,itt)*.01) - ptarget ) if (dp.gt.dpmin) exit loop_kIN dpmin=min(dpmin,dp) enddo loop_kIN pbot=max(pres_field(i,j,1,itt),psfc(i,j,itt)) zbot=min(data_in(i,j,1,itt)/9.81,ter(i,j)) tbotextrap=tk(i,j,kupper,itt)*(pbot/pres_field(i,j,kupper,itt))**expon tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt)) data_out(i,j,k,itt) = (zbot+tvbotextrap/.0065*(1.-(interp_levels(k)/pbot)**expon))*9.81 ENDIF enddo enddo enddo enddo !!! Code for filling missing data with an average - we don't want to do this !!do itt = 1, it !!loop_levels : do k = 1, num_metgrid_levels !!sumA = SUM(data_out(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING) !!sumN = SUM(N(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING) !!IF ( sumN == 0. ) CYCLE loop_levels !!AVE_geopt = sumA/sumN !!WHERE ( data_out(:,:,k,itt) == MISSING ) !!data_out(:,:,k,itt) = AVE_geopt !!END WHERE !!end do loop_levels !!end do END IF !!! All other fields and geopt at higher levels come here do itt = 1, it do j = 1, iy do i = 1, ix do k = 1, kk if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,1,itt) end do do k = kk+1, num_metgrid_levels if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,iz,itt) end do end do end do end do END SUBROUTINE interp SUBROUTINE int1D(xxout, xxin, ppin, ppout, npin, npout, LINLOG, MISSING) ! Modified from int2p - NCL code ! routine to interpolate from one set of pressure levels ! . to another set using linear or ln(p) interpolation ! ! NCL: xout = int2p (pin,xin,pout,linlog) ! This code was originally written for a specific purpose. ! . Several features were added for incorporation into NCL's ! . function suite including linear extrapolation. ! ! nomenclature: ! ! . ppin - input pressure levels. The pin can be ! . be in ascending or descending order ! . xxin - data at corresponding input pressure levels ! . npin - number of input pressure levels >= 2 ! . ppout - output pressure levels (input by user) ! . same (ascending or descending) order as pin ! . xxout - data at corresponding output pressure levels ! . npout - number of output pressure levels ! . linlog - if abs(linlog)=1 use linear interp in pressure ! . if abs(linlog)=2 linear interp in ln(pressure) ! . missing- missing data code. ! ! input types INTEGER :: npin,npout,linlog,ier real :: ppin(npin),xxin(npin),ppout(npout) real :: MISSING logical :: AVERAGE ! ! output real :: xxout(npout) INTEGER :: j1,np,nl,nin,nlmax,nplvl INTEGER :: nlsave,np1,no1,n1,n2,nlstrt real :: slope,pa,pb,pc ! automatic arrays real :: pin(npin),xin(npin),p(npin),x(npin) real :: pout(npout),xout(npout) xxout = MISSING pout = ppout p = ppin x = xxin nlmax = npin ! exact p-level matches nlstrt = 1 nlsave = 1 do np = 1,npout xout(np) = MISSING do nl = nlstrt,nlmax if (pout(np).eq.p(nl)) then xout(np) = x(nl) nlsave = nl + 1 go to 10 end if end do 10 nlstrt = nlsave end do if (LINLOG.eq.1) then do np = 1,npout do nl = 1,nlmax - 1 if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then slope = (x(nl)-x(nl+1))/ (p(nl)-p(nl+1)) xout(np) = x(nl+1) + slope* (pout(np)-p(nl+1)) end if end do end do elseif (LINLOG.eq.2) then do np = 1,npout do nl = 1,nlmax - 1 if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then pa = log(p(nl)) pb = log(pout(np)) ! special case: in case someone inadvertently enter p=0. if (p(nl+1).gt.0.d0) then pc = log(p(nl+1)) else pc = log(1.d-4) end if slope = (x(nl)-x(nl+1))/ (pa-pc) xout(np) = x(nl+1) + slope* (pb-pc) end if end do end do end if ! place results in the return array; xxout = xout END SUBROUTINE int1D FUNCTION virtual (tmp,rmix) ! This function returns virtual temperature in K, given temperature ! in K and mixing ratio in kg/kg. real :: tmp, rmix, virtual virtual=tmp*(0.622+rmix)/(0.622*(1.+rmix)) END FUNCTION virtual END MODULE module_ForInterpolate