program create_coordinates
!-----------------------------------------------------------
!
!	  to create coordinates file (e.g. 1_coordinates.nc)
!    for fine grid from coarse grid (e.g. coordinates.nc)
!
!   to make it we use a 4th order polynomial interpolation
!
!	        Created by Brice Lemaire on 12/2009.
!
!-----------------------------------------------------------
USE netcdf
USE agrif_readwrite
USE tools_brice
!------------
IMPLICIT NONE
   !
   TYPE(coordinates) :: G1, G2, G3, G4     !from agrif_types from agrif_readwrite
   TYPE(coordinates) :: fine_grid
   TYPE(mixed_coordinates) :: Gmix      !from agrif_types from agrif_readwrite
   !INTEGER :: nxG1, nyG1, nxG2, nyG2, nxG3, nyG3, nxG4, nyG4
   INTEGER :: narg,iargc
   INTEGER :: status
   CHARACTER(len=80) :: namelistname, Child_filename
   !
   !!!Read input file (namelist file)
   narg = iargc()
   !
   IF (narg == 0) THEN
	   namelistname = 'namelist.input'
   ELSE
	   CALL getarg(1,namelistname)
   ENDIF
   !
   CALL read_namelist(namelistname)      !imin,imax,jmin,jmax,parent_coordinate_file
   !
   !!!Read parent coordinate file
   IF((imin.LT.imax).AND.(jmin.LT.jmax)) THEN
      status = Read_Coordinates(TRIM(parent_coordinate_file),G1)      !from agrif_readwrite
   ELSEIF((imin.GE.imax).AND.(jmin.LT.jmax)) THEN
      status = Read_Coordinates(TRIM(parent_coordinate_file),G1)      !from agrif_readwrite
      status = Read_Coordinates(TRIM(parent_coordinate_file),G2)      !from agrif_readwrite
   ELSEIF((imin.GE.imax).AND.(jmin.GE.jmax)) THEN
	  status = Read_Coordinates(TRIM(parent_coordinate_file),G1)      !from agrif_readwrite
      status = Read_Coordinates(TRIM(parent_coordinate_file),G2)      !from agrif_readwrite
      status = Read_Coordinates(TRIM(parent_coordinate_file),G3)      !from agrif_readwrite
      status = Read_Coordinates(TRIM(parent_coordinate_file),G4)      !from agrif_readwrite
   ENDIF
   ! 
   !
   !!!Create the mixed grid with U,V,F & W values
   IF((imin.LT.imax).AND.(jmin.LT.jmax)) THEN
      CALL mixed_grid(G1,Gmix)
   ELSEIF((imin.GE.imax).AND.(jmin.LT.jmax)) THEN
      CALL mixed_grid(G1,G2,Gmix)
   ELSEIF((imin.GE.imax).AND.(jmin.GE.jmax)) THEN
      CALL mixed_grid(G1,G2,G3,G4,Gmix)      !from tools_brice
   !ELSEIF()
   ENDIF
   
   
!print*, 'G1%nav_lat= '
!DO j=jmin-1,jmax+1
!print*, G1%nav_lat(imin-1:imax+1,j)
!END DO  

!print*, 'G1%nav_lon= '
!DO j=jmin-1,jmax+1
!print*, G1%nav_lon(imin-1:imax+1,j)
!END DO 
   !
   !!!Interpolation inside the mixed grid
   IF(rho>1) THEN
      CALL interp_grid(Gmix)      !from tools_brice
   ENDIF
   !   
   !!!Define name of child coordinate file
   CALL set_child_name(parent_coordinate_file,Child_filename)      !coordinates.nc -> 1_coordinates.nc
   !
   !!!Allocation of child grid elements
   CALL agrif_grid_allocate(fine_grid,nx_fine,ny_fine)      !from agrif_types
   !
   !!!Break the mixed grid Gmix into 4 grids (G1%u, G1%v, G1%t, GI%f)
   CALL alloc_child_grid(Gmix,fine_grid)      !from tools_brice
   !
print*, 'fine_grid%glamf= '
DO j=1,ny_fine
print*, fine_grid%glamf(:,j)
END DO 
    !!!Read parent coordinate file
   status = Write_Coordinates(Child_filename,fine_grid,nx_fine,ny_fine)      !from agrif_readwrite
   !
end program create_coordinates