! THIS PROGRAM TO CONVERT THE OUTPUT FROM INTERPOLATION METHODS TO NETCDF FORMAT
	PROGRAM CONVERT
	use netcdf
	implicit none
	! For input
	integer, parameter :: SYR = 1980, EYR = 2010, NYR = EYR - SYR + 1, NMN = 12
	integer, parameter :: NLON = 35, NLAT = 63 
	integer :: NDY
	integer :: ilon, ilat, iyr, imn, idy, itmp, year
	
	character (len = 200) :: wdir, odir, filenm
	character (len = 2), parameter, dimension (NMN) :: mon = (/"01","02","03","04","05","06","07","08","09","10","11","12"/)
	character (len = 10) :: cyr, cdy, ctmp
	
	real, dimension (NLON,NLAT,15000) :: var_in
	real, dimension (NLON,NLAT) :: var_tmp
	real, parameter :: SLON = 101.5, SLAT = 8.0, INC = 0.25, NA = -99.0
	real :: lon(NLON), lat(NLAT), rtmp
	! For output
	integer :: ncid
	integer :: lon_dimid, lat_dimid, lon_varid, lat_varid
	integer :: var_varid
	integer :: NREC, irec, rec_dimid
	integer :: start(3), icount(3), dimids(3)
	
	character (len = *), parameter :: FILE_OUT = "TMP.nc"
	character (len = *), parameter :: LAT_NAME = "latitude", LON_NAME = "longitude", REC_NAME = "time"
	character (len = *), parameter :: VAR_NAME = "var"
	character (len = *), parameter :: UNITS = "units"
	character (len = *), parameter :: VAR_UNITS = "var_unit"
	character (len = *), parameter :: LAT_UNITS = "degrees_north"
	character (len = *), parameter :: LON_UNITS = "degrees_east"

	! =============== BEGIN PROGRAM =============
	wdir = "/work/users/thanhnx/R_sta_2_grid/convert_2_netcdf/"
	odir = "/work/users/thanhnx/R_sta_2_grid/convert_2_netcdf/"

	! ------ Read GRID info ------
	print*,"   Read GRID information (lon, lat) ..."
	open(3,file = trim(wdir)//"grid_point_Z.txt",status="old")
	   do ilon = 1, NLON
	      do ilat = 1 , NLAT
	         read(3,*) lon(ilon), lat(ilat), rtmp
	      enddo
	   enddo
	close(3)

	! ------ Read INPUT data ------
	print*, "   Read INPUT data ..."
	var_in = NA
	irec = 1
	do iyr = 1, NYR
	   year = SYR + iyr - 1
	   write(cyr,'(i4)') year
	   print*,"      Read data in year ...   ", year
	   do imn = 1, NMN
	      NDY = cnday(year,imn)
	      do idy = 1, NDY
	         if (idy < 10) then
	            write(cdy,'(i1)') idy
	            filenm=trim(wdir)//"input/Gsmap/"//trim(cyr)//"/"//trim(mon(imn))//"/R_day_0"//trim(cdy)//".txt"
	         else
	            write(cdy,'(i2)') idy
	            filenm=trim(wdir)//"input/Gsmap/"//trim(cyr)//"/"//trim(mon(imn))//"/R_day_"//trim(cdy)//".txt"
	         endif
	         !print*,filenm
	         open(4,file = filenm,status="old")
	            do ilon = 1, NLON
	               do ilat = 1, NLAT
	                  read(4,*) rtmp, rtmp, var_in(ilon,ilat,irec)
	               enddo
	            enddo
	         close(4)
	         irec = irec + 1
	      enddo
	   enddo
	enddo
	NREC = irec - 1

	! ------ Write OUTPUT data -----
	call check( nf90_create(FILE_OUT, nf90_clobber, ncid) )
	
	call check( nf90_def_dim(ncid, LAT_NAME, NLAT, lat_dimid) )
	call check( nf90_def_dim(ncid, LON_NAME, NLON, lon_dimid) )
	call check( nf90_def_dim(ncid, REC_NAME, NF90_UNLIMITED, rec_dimid) )
	
	call check( nf90_def_var(ncid, LAT_NAME, NF90_REAL, lat_dimid, lat_varid) )
	call check( nf90_def_var(ncid, LON_NAME, NF90_REAL, lon_dimid, lon_varid) )
	
	call check( nf90_put_att(ncid, lat_varid, UNITS, LAT_UNITS) )
	call check( nf90_put_att(ncid, lon_varid, UNITS, LON_UNITS) )
	
	dimids = (/ lon_dimid, lat_dimid, rec_dimid /)
	call check( nf90_def_var(ncid, VAR_NAME, NF90_REAL, dimids, var_varid) )
	call check( nf90_put_att(ncid, var_varid, UNITS, VAR_UNITS) )
	call check( nf90_enddef(ncid) )
	! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	call check( nf90_put_var(ncid, lat_varid, lat) )
	call check( nf90_put_var(ncid, lon_varid, lon) )
	
	icount = (/ NLON, NLAT, 1/)
	start = (/ 1, 1, 1 /)
	
	do irec = 1, NREC
	   start(3) = irec
	   var_tmp=var_in(:,:,irec)
	   call check( nf90_put_var(ncid, var_varid, var_tmp, start = start, count = icount) )
	enddo
	
	call check( nf90_close(ncid) )
	
	print *,">>>>>> SUCCESS converting file to netcdf format !!!"


	CONTAINS
	! =========================== DON'T EDIT ANY THING BELOW ======================
	integer Function cnday(year,imon)
	integer :: year, imon, check
	cnday=31
	if (imon == 2) then
	   cnday=28
	   check=year/4
	   check=check*4
	   if (check == year) then
	      cnday=29
	      check=year/100
	      check=check*100
	      if (check == year) then
	         cnday=28
	         check=year/400
	         check=check*400
	         if (check == year) then
	            cnday=29
	         endif
	      endif
	   endif
	endif
	if ((imon==4).or.(imon==6).or.(imon==9).or.(imon==11)) then
	   cnday=30
	endif
	return
	End Function

	Subroutine check(status)
	integer, intent ( in) :: status

	if (status /= nf90_noerr) then
	   print *, trim(nf90_strerror(status))
	   stop "Stopped"
	end if
	End Subroutine check
	END