PROGRAM rearrange_startphy

USE netcdf

IMPLICIT NONE

INTEGER :: ncid
INTEGER :: dimids(4)
INTEGER :: varid
INTEGER :: ierr
INTEGER :: nvar
INTEGER :: ndim
CHARACTER(LEN=100) :: varname
INTEGER :: varname_dimids(4)

INTEGER :: i,j,k,nb_cells
INTEGER,ALLOCATABLE :: cell_index(:)
REAL,ALLOCATABLE :: ref_lon(:)
REAL,ALLOCATABLE :: ref_lat(:)
REAL,ALLOCATABLE :: lon(:)
REAL,ALLOCATABLE :: lat(:)
REAL,ALLOCATABLE :: ref_field(:),field(:)
REAL :: lati,loni
REAL :: diff_lat,diff_lon
CHARACTER(LEN=*),PARAMETER :: ref_file="startphy_icosa_ref.nc"
CHARACTER(LEN=*),PARAMETER :: file="startphy_icosa.nc"

! load coordinates from files
  ierr=NF90_OPEN(ref_file, NF90_NOWRITE, ncid)
  ierr=NF90_INQ_DIMID(ncid,"points_physiques",dimids(1))
  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(1), len=nb_cells)
  write(*,*) "nb_cells=",nb_cells,"dimids(1)=",dimids(1)
  allocate(ref_lat(nb_cells),ref_lon(nb_cells))
  
  ierr=NF90_INQ_VARID(ncid,"latitude",varid)
  ierr=NF90_GET_VAR(ncid,varid,ref_lat)
  if (ierr /= nf90_noerr) then
    write(*,*) "cannot load latitude from ",trim(ref_file)
  else
    write(*,*) "ref_lat(1:5)=",ref_lat(1:5)
  endif
  ierr=NF90_INQ_VARID(ncid,"longitude",varid)
  ierr=NF90_GET_VAR(ncid,varid,ref_lon)
  if (ierr /= nf90_noerr) then
    write(*,*) "cannot load longitude from ",trim(ref_file)
  else
    write(*,*) "ref_lon(1:5)=",ref_lon(1:5)
  endif
  ierr=NF90_CLOSE(ncid)
  
  ierr=NF90_OPEN(file, NF90_WRITE, ncid)
  allocate(lat(nb_cells),lon(nb_cells))
  ierr=NF90_INQ_VARID(ncid,"lat",varid)
  ierr=NF90_GET_VAR(ncid,varid,lat)
  if (ierr /= nf90_noerr) then
    write(*,*) "cannot load lat from ",trim(file)
  else
    write(*,*) "lat(1:5)=",lat(1:5)
  endif
  ierr=NF90_INQ_VARID(ncid,"lon",varid)
  ierr=NF90_GET_VAR(ncid,varid,lon)
  if (ierr /= nf90_noerr) then
    write(*,*) "cannot load lat from ",trim(file)
  else
    write(*,*) "lon(1:5)=",lon(1:5)
  endif
  
  ! find correspondances between lon/ref_lon & lat/ref_lat
  allocate(cell_index(nb_cells))
  cell_index(:)=0
  do i=1,nb_cells
    lati=lat(i) ; loni=lon(i)
    do j=1,nb_cells
     if (lati/=0) then ! use relative difference, if possible
       diff_lat=abs((lati-ref_lat(j))/lati)
     else
       diff_lat=abs(lati-ref_lat(j))
     endif
     if (loni/=0) then ! use relative difference, if possible
       diff_lon=abs((loni-ref_lon(j))/loni)
     else
       diff_lon=abs(loni-ref_lon(j))
     endif
     
     if ((diff_lat <= 1.e-6).and.(diff_lon <= 1.e-6)) then
      cell_index(i)=j
      write(*,*)"j=",j," lati=",lati," ref_lat(j)=",ref_lat(j)
      write(*,*)"j=",j," loni=",loni," ref_lon(j)=",ref_lon(j)
     endif
    enddo ! of do j=1,nb_cells
    write(*,*) ">> i=",i,"cell_index(i)=",cell_index(i)

    ! sanity check:
    if (cell_index(i)==0) then
      write(*,*) "Error, could not find lon-lat match for i=",i
      write(*,*) "       lat(i)=",lat(i)," lon(i)=",lon(i)
      stop
    endif
  enddo ! of do i=1,nb_cells
  
!  do i=1,nb_cells
!    write(*,*) "i=",i," cell_index(i)=",cell_index(i)
!    write(*,*) "  lat(i)=",lat(i),"ref_lat(cell_index(i))=",ref_lat(cell_index(i))
!  enddo
  
  ! load, rearrange and write variables
  ierr=NF90_INQ_DIMID(ncid,"cell",dimids(1))
  ierr=NF90_INQUIRE(ncid,nVariables=nvar)
  write(*,*) "nvar=",nvar
  allocate(ref_field(nb_cells),field(nb_cells))
  ! loop over variables:
  do k=1,nvar
    ierr=NF90_INQUIRE_VARIABLE(ncid,k,name=varname,ndims=ndim,dimids=varname_dimids)
    if (ierr /= nf90_noerr) then
      write(*,*) "error for variable k=",k
      write(*,*) nf90_strerror(ierr)
    endif
    if ((ndim==1).and.(varname_dimids(1)==dimids(1))) then
      write(*,*) "processing ",trim(varname)
      ! load field
      ierr=NF90_INQ_VARID(ncid,varname,varid)
      ierr=NF90_GET_VAR(ncid,varid,ref_field)
      ! rearrange field
      do i=1,nb_cells
        field(cell_index(i))=ref_field(i)
      enddo
      ! write field
      ierr=NF90_PUT_VAR(ncid,varid,field)
    endif
  enddo
  
END PROGRAM rearrange_startphy