program extractcolumnfor1D

! program to extract a column from a 3D output file (like diagfi.nc or stats.nc)
! and generate a start file for the 1D model.
! usage
! extractcolumnfor1D -i # -j # -l # -v list infile
! list == comma separated list of variables names
use netcdf

implicit none

! input netcdf file:
character(len=256) :: infile
! output files "profile_*"
character(len=256) :: outfile

! NetCDF stuff
integer :: status ! NetCDF return code
integer :: inid ! NetCDF file IDs
integer :: lon_dimid, lat_dimid, alt_dimid, time_dimid
integer :: varid ! to store the ID of a variable
integer :: datashape(4) ! to store dimension IDs of a given dataset
integer :: corner(4),edges(4) ! to read data with a time axis

! array sizes of dataset
integer :: lonlen
integer :: latlen
integer :: altlen
integer :: timelen

integer :: i,j
integer :: lonindex,latindex,timeindex ! indexes for extraction

integer :: nbvar ! number of variables
character(len=90),allocatable :: varnames(:)
character(len=90) :: surfvarname ! for surface variable
character(len=256) :: varlist ! comma separated list of variables
integer :: comma_index(100) ! index of commas in varlist

integer :: nbarg ! # of program arguments
character(len=256) :: arg,arg2 ! to store a program argument
integer,parameter :: expected_nbarg=9 ! expected number of arguments

!variables to read and write:
real, allocatable :: variable(:)
real :: surface ! surface value 

! 0. Preliminary stuff, check arguments
! get number of arguments this program was called with
nbarg=command_argument_count()

if (nbarg.ne.expected_nbarg) then
  write(*,*) "Wrong number of arguments!"
  write(*,*) "Usage:"
  write(*,*) "extractcolumnfor1D.e -i # -j # -l # -v list_of_variables infile"
  stop
else
  ! parse the argument list 
  do i=1,expected_nbarg-2,2
    call get_command_argument(i,arg) ! get argument # i
    select case (trim(arg))
      case ("-i")
        call get_command_argument(i+1,arg2)
        read(arg2,*) lonindex
      case ("-j")
        call get_command_argument(i+1,arg2)
        read(arg2,*) latindex
      case ("-l")
        call get_command_argument(i+1,arg2)
        read(arg2,*) timeindex
      case ("-v")
        call get_command_argument(i+1,arg2)
        varlist=adjustl(trim(arg2))
      case default
        write(*,*) "invalid option:", trim(arg)
    end select
  enddo
  call get_command_argument(expected_nbarg,arg) ! last argument (input file)
  infile=trim(arg)
endif

! 1.1 open netcdf input file
status=nf90_open(infile,NF90_NOWRITE,inid)
if (status.ne.nf90_noerr) then
  write(*,*)"Failed to open datafile ",trim(infile)
  write(*,*)trim(nf90_strerror(status))
 stop
endif
write(*,*) "Reading input file: ",trim(infile)

! 1.2. get dimensions from file
!longitude
status=nf90_inq_dimid(inid,"longitude",lon_dimid)
if (status.ne.nf90_noerr) then
  write(*,*)"Failed to find longitude dimension"
  write(*,*)trim(nf90_strerror(status))
 stop
endif
status=nf90_inquire_dimension(inid,lon_dimid,len=lonlen)
if (status.ne.nf90_noerr) then
  write(*,*)"Failed to find longitude length"
  write(*,*)trim(nf90_strerror(status))
endif

! latitude
status=nf90_inq_dimid(inid,"latitude",lat_dimid)
if (status.ne.nf90_noerr) then
  write(*,*)"Failed to find latitude dimension"
  write(*,*)trim(nf90_strerror(status))
 stop
endif
status=nf90_inquire_dimension(inid,lat_dimid,len=latlen)
if (status.ne.nf90_noerr) then
  write(*,*)"Failed to find latitude length"
  write(*,*)trim(nf90_strerror(status))
endif

!altitude
status=nf90_inq_dimid(inid,"altitude",alt_dimid)
if (status.ne.nf90_noerr) then
  write(*,*)"Failed to find altitude dimension"
  write(*,*)trim(nf90_strerror(status))
 stop
endif
status=nf90_inquire_dimension(inid,alt_dimid,len=altlen)
if (status.ne.nf90_noerr) then
  write(*,*)"Failed to find altitude length"
  write(*,*)trim(nf90_strerror(status))
endif
allocate(variable(altlen))

!time
status=nf90_inq_dimid(inid,"Time",time_dimid)
if (status.ne.nf90_noerr) then
  write(*,*)"Failed to find Time dimension"
  write(*,*)trim(nf90_strerror(status))
 stop
endif
status=nf90_inquire_dimension(inid,time_dimid,len=timelen)
if (status.ne.nf90_noerr) then
  write(*,*)"Failed to find Time length"
  write(*,*)trim(nf90_strerror(status))
endif

! 1.3 scan varnames(:) to extract number of variablas and their names
write(*,*) "list of variables to process: ",trim(varlist)
nbvar=1
do i=1,len(varlist)
  if(varlist(i:i)==",") then
    comma_index(nbvar)=i
    ! write(*,*)"nbvar=",nbvar," comma_index(nbvar)=",comma_index(nbvar)
    nbvar=nbvar+1
  endif
enddo
allocate(varnames(nbvar))

if (nbvar>1) then
 read(varlist(1:comma_index(1)-1),'(a)') varnames(1)
 ! write(*,*)"varnames(1):",trim(varnames(1))
 do i=2,nbvar-1
  read(varlist(comma_index(i-1)+1:comma_index(i)-1),'(a)') varnames(i)
  ! write(*,*)"varnames(i):",trim(varnames(i))
 enddo
 read(varlist(comma_index(nbvar-1)+1:),'(a)') varnames(nbvar)
 ! write(*,*)"varnames(nbvar):",trim(varnames(nbvar))
else
 varnames(1)=adjustl(trim(varlist))
 ! write(*,*)"varnames(1):",trim(varnames(1))
endif ! of if (nbvar>1)

! 2. Load data and write it to file

do i=1,nbvar
  status=nf90_inq_varid(inid,trim(varnames(i)),varid)
  if (status.ne.nf90_noerr) then
    write(*,*)"Failed to find variable ",trim(varnames(i))
  else
    status=nf90_get_var(inid,varid,variable(1:altlen),&
                        start=[lonindex,latindex,1,timeindex],&
                        count=[1,1,altlen,1])
    
    ! 2.2 try to find a corresponding "surface" value
    if(trim(varnames(i))=="temp") then
      surfvarname="tsurf" ! surface temperature
    else
      surfvarname=trim(varnames(i))//"_s"
    endif
    status=nf90_inq_varid(inid,surfvarname,varid)
    if (status.ne.nf90_noerr) then
      write(*,*)"Failed to find surface variable ",trim(surfvarname)
      write(*,*)" set its value to zero"
      surface=0.
    else
      status=nf90_get_var(inid,varid,surface,&
                          start=[lonindex,latindex,timeindex])
    endif
    
    ! 2.3 write to file
    outfile="profile_"//trim(varnames(i))
    open(42,file=trim(outfile))
    ! add a first line for surface value
    write(42,*) surface
    do j=1,altlen
      write(42,*) variable(j)
    enddo
    close(42)
    write(*,*) "wrote file profile_"//trim(varnames(i))
  endif

enddo ! of do i=1,nbvar

end program extractcolumnfor1D