program genCol

! ----------------------------------------------------------------------------
! This program can re-compute the column density (kg m-2)
! of any tracer from its mmr, vmr or concentration
! input : diagfi.nc  / concat.nc / stats.nc kind of files
! output : input_col.nc including lon,lat,alt,time and all
! specified tracer_col
! author: F. Forget, M. Alhossani, J.Mauxion
! ----------------------------------------------------------------------------

implicit none

include "netcdf.inc" ! NetCDF definitions

!!! I/O variables
character (len=128) :: infile ! input file name (diagfi.nc or stats.nc format)
character (len=128) :: infile2 ! second input file (may be needed for 'aire' and 'cv')
integer infid ! NetCDF input file ID (of diagfi.nc or stats.nc format)
integer infid2 ! NetCDF input file which contains 'phisini' dataset (diagfi.nc)
character (len=100) :: filename
integer :: ierr ! ierr: [netcdf] subroutine returned error code
integer :: nout,latdimout,londimout,altdimout,timedimout
! nout: [netcdf] output file ID
! latdimout: [netcdf] output latitude (dimension) ID
! londimout: [netcdf] output longitude (dimension) ID
! altdimout: [netcdf] output altitude (dimension) ID
! timedimout: [netcdf] output time (dimension) ID

!!! Temporary variables
character (len=64) :: tmpvarname ! temporary store a variable name
character (len=64) :: tmpunit ! temporary store a unit
integer :: tmpvarid ! varid: temporary store a variable ID #
integer :: tmpdimid ! temporarily store a dimension ID
integer tmpndims ! temporarily store # of dimensions of a variable
real :: tmpmm ! Temporarily store the molar mass of current specie

!!! Auto-detection and tracer selection
integer nbvarinfile ! # of variables in input file
integer nbtrinfile ! # of 4D mmr, vmr or concentration tracers in input file
integer nbvar ! # of variables to process
character (len=64), dimension(:), allocatable :: var
! var(): name(s) of variable(s) that will be concatenated

!!! Grid-related variables
real, dimension(:), allocatable:: lat,lon,alt,time
! lat(): array, stores latitude coordinates
! lon(): array, stores longitude coordinates
! alt(): array, stores altitude coordinates
! time(): array, stores time coordinates
integer lonlength ! # of grid points along longitude
integer latlength ! # of grid points along latitude
integer altlength ! # of grid point along altitude (of input datasets)
integer timelength ! # of points along time
real,dimension(:),allocatable :: ap,bp ! hybrid vertical coordinates of levels
real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates of layerS

!!! Loops-related variables
integer :: ilat,ilon,ialt,it,i
! ilat: index for loops on latitude
! ilon: index for loops on longitude 
! ialt: index for loops on altitude
! it: index for loops on time
! i: index for loops

!!! Constants
real :: g ! Martian gravitational acceleration
real :: nav ! Avogadro number
real :: rgp ! Perfect gas constant
real :: Rmean ! Mean atmospheric molar mass

!!! Booleans to handle various setup configuration
logical :: isco, ismmr, isvmr ! is the tracer co, mmr, vmr ?
logical :: ismmcst ! true if mm assumed constant, false otherwise
logical :: foundmm ! to store whether we find molar mass in input file
logical :: foundpress ! to store whether we find pressure in input file
logical :: foundrho ! to store whether we find density in input file

!!! Incrementing
real :: incMass ! mass of air in each layer
real :: tracerMass ! cumulated column density
real :: cstmm ! Constant astmospheric molar mass

!!! Variable to read from input file
real,dimension(:,:,:),allocatable :: ps ! GCM surface pressure
real,dimension(:,:,:,:),allocatable :: tracer ! vmr,mmr or concentration
real,dimension(:,:,:,:),allocatable :: temp ! atmospheric temperature
real,dimension(:,:,:,:),allocatable :: rho ! atmospheric density
real,dimension(:,:,:,:),allocatable :: press ! atmospheric pressure
real,dimension(:,:,:,:),allocatable :: locmm ! local atmospheric molar mass

!!! New variable to write in output file
integer :: tracvarout ! var id for column density
real,dimension(:,:,:),allocatable :: tracerCol ! total tracer column

!===============================================================================
! 1.1 Input file
!===============================================================================

write(*,*) ""
write(*,*) "Enter input file name:"

read(*,'(a128)') infile
write(*,*) ""

! open input file

ierr = NF_OPEN(infile,NF_NOWRITE,infid)
if (ierr.ne.NF_NOERR) then
   write(*,*) 'ERROR: Pb opening input file'
   stop ""
endif

!===============================================================================
! 1.2 Get # and names of tracers in input file
!===============================================================================

! Get number of var in input file
ierr=NF_INQ_NVARS(infid,nbvarinfile)
if (ierr.ne.NF_NOERR) then
   write(*,*) 'ERROR: Failed getting number of variables from file'
   write(*,*) NF_STRERROR(ierr)
   stop ""
endif

! Find which var is a tracer
write(*,*)" The following tracers have been found:"
nbtrinfile=0
do i=1,nbvarinfile
   ! Get name of variable # i
   ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
   ! Get dimension and unit
   ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
   ierr=NF_GET_ATT_TEXT(infid,i,'units',tmpunit)
   ismmr = (tmpunit.eq.'kg/kg').or.(tmpunit.eq.'kg.kg-1')
   isvmr = (tmpunit.eq.'mol/mol').or.(tmpunit.eq.'mol.mol-1')
   isco = (tmpunit.eq.'cm-3').or.(tmpunit.eq.'m-3')
   if ((tmpndims.eq.4).and.(ismmr.or.isvmr.or.isco)) then
      nbtrinfile=nbtrinfile+1
      write(*,*) trim(tmpvarname), " ", trim(tmpunit)
   endif
enddo

! Allocate array for kept tracers
allocate(var(nbtrinfile),stat=ierr)
if (ierr.ne.0) then
  write(*,*) "Error: failed memory allocation of var(nbtrinfile)"
  stop ""
endif

! Ask for tracers to compute
write(*,*) ""
write(*,*) "Which variable do you want to keep?"
write(*,*) "all or list of <variables> (separated by <Return>s)"
write(*,*) "(an empty line , i.e: just <Return>, implies end of list)"
nbvar=0
read(*,'(a64)') tmpvarname
do while ((tmpvarname.ne.' ').and.(trim(tmpvarname).ne.'all'))
   ! check if tmpvarname is valid
   ierr=NF_INQ_VARID(infid,tmpvarname,tmpvarid)
   if (ierr.eq.NF_NOERR) then ! valid name
      ! also check that it is indeed 4D tracer with consistent unit
      ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims)
      ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'units',tmpunit)
      ismmr = (tmpunit.eq.'kg/kg').or.(tmpunit.eq.'kg.kg-1')
      isvmr = (tmpunit.eq.'mol/mol').or.(tmpunit.eq.'mol.mol-1')
      isco = (tmpunit.eq.'cm-3').or.(tmpunit.eq.'m-3')
      if ((tmpndims.eq.4).and.(ismmr.or.isvmr.or.isco)) then
         nbvar=nbvar+1
         var(nbvar)=tmpvarname
      else ! invalid dimension or unit
         write(*,*) 'Error: ',trim(tmpvarname),' has wrong dimension and/or unit'
         write(*,*) '       (we''ll skip that one)'
      endif
   else ! invalid name
      write(*,*) 'Error: ',trim(tmpvarname),' is not a valid name'
      write(*,*) '       (we''ll skip that one)'
   endif
   read(*,'(a64)') tmpvarname
enddo

! handle "all" case
if (tmpvarname.eq.'all') then
   nbvar=0
   do i=1,nbvarinfile
      ! look for 4D tracers with good unit
      ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
      ierr=NF_GET_ATT_TEXT(infid,i,'units',tmpunit)
      ismmr = (tmpunit.eq.'kg/kg').or.(tmpunit.eq.'kg.kg-1')
      isvmr = (tmpunit.eq.'mol/mol').or.(tmpunit.eq.'mol.mol-1')
      isco = (tmpunit.eq.'cm-3').or.(tmpunit.eq.'m-3')
      if ((tmpndims.eq.4).and.(ismmr.or.isvmr.or.isco)) then
         nbvar=nbvar+1
         ! get the corresponding name
         ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
         var(nbvar)=tmpvarname
      endif
   enddo
endif

! Check that there is at least 1 variable to process
if (nbvar.eq.0) then
   write(*,*) 'No variables to process !?'
   write(*,*) 'Might as well stop here'
   stop ""
else
   write(*,*) ""
   write(*,*) 'OK, the following variables will be processed:'
   do i=1,nbvar
      write(*,*) var(i)
   enddo
endif


!===============================================================================
! 1.3 Get grids in lon,lat,alt,time,
!     as well as hybrid coordinates aps() and bps()
!===============================================================================

! 1.3.1 longitude, latitude, altitude and time

! latitude
ierr=NF_INQ_DIMID(infid,"latitude",tmpdimid)
if (ierr.ne.NF_NOERR) then
   stop "Error: Failed to get latitude dimension ID"
else
   ierr=NF_INQ_VARID(infid,"latitude",tmpvarid)
   if (ierr.ne.NF_NOERR) then
      stop "Error: Failed to get latitude ID"
   else
      ierr=NF_INQ_DIMLEN(infid,tmpdimid,latlength)
      if (ierr.ne.NF_NOERR) then
         stop "Error: Failed to get latitude length"
      else
         allocate(lat(latlength))
         ierr=NF_GET_VAR_REAL(infid,tmpvarid,lat)
         if (ierr.ne.NF_NOERR) then
            write(*,*) "Error: Failed reading latitude"
            stop ""
         endif
      endif
   endif
endif

! longitude
ierr=NF_INQ_DIMID(infid,"longitude",tmpdimid)
if (ierr.ne.NF_NOERR) then
   stop "Error: Failed to get longitude dimension ID"
else
   ierr=NF_INQ_VARID(infid,"longitude",tmpvarid)
   if (ierr.ne.NF_NOERR) then
      stop "Error: Failed to get longitude ID"
   else
      ierr=NF_INQ_DIMLEN(infid,tmpdimid,lonlength)
      if (ierr.ne.NF_NOERR) then
         stop "Error: Failed to get longitude length"
      else
         allocate(lon(lonlength))
         ierr=NF_GET_VAR_REAL(infid,tmpvarid,lon)
         if (ierr.ne.NF_NOERR) then
            write(*,*) "Error: Failed reading longitude"
            stop ""      
         endif
      endif
   endif
endif

! time
ierr=NF_INQ_DIMID(infid,"Time",tmpdimid)
if (ierr.ne.NF_NOERR) then
   stop "Error: Failed to get Time dimension ID"
else
   ierr=NF_INQ_VARID(infid,"Time",tmpvarid)
   if (ierr.ne.NF_NOERR) then
      stop "Error: Failed to get Time ID"
   else
      ierr=NF_INQ_DIMLEN(infid,tmpdimid,timelength)
      if (ierr.ne.NF_NOERR) then
         stop "Error: Failed to get Time length"
      else
         allocate(time(timelength))
         ierr=NF_GET_VAR_REAL(infid,tmpvarid,time)
         if (ierr.ne.NF_NOERR) then
            write(*,*) "Error: Failed reading Time"
            stop ""
         endif
      endif
   endif
endif

! altitude
ierr=NF_INQ_DIMID(infid,"altitude",tmpdimid)
if (ierr.ne.NF_NOERR) then
   stop "Error: Failed to get altitude dimension ID"
else
   ierr=NF_INQ_VARID(infid,"altitude",tmpvarid)
   if (ierr.ne.NF_NOERR) then
      stop "Error: Failed to get altitude length"
   else
      ierr=NF_INQ_DIMLEN(infid,tmpdimid,altlength)
      if (ierr.ne.NF_NOERR) then
         stop "Error: Failed to get altitude length"
      else   
         allocate(alt(altlength))
         ierr=NF_GET_VAR_REAL(infid,tmpvarid,alt)
         if (ierr.ne.NF_NOERR) then
            write(*,*) "Error: Failed reading Altitude"
            stop ""
         endif
      endif
   endif
endif

! 1.3.2 Get hybrid coordinates

! Look for hybrid coordinates ap and bp
allocate(ap(altlength+1))
allocate(bp(altlength+1))
ierr=NF_INQ_VARID(infid,"ap",tmpvarid)
if (ierr.ne.NF_NOERR) then
   ! If ap not in input file, it is in a diagfi file
   write(*,*) "Failed to get ap ID from input file ",trim(infile)
   infile2="diagfi.nc"
   write(*,*) "Trying file ",trim(infile2)
   ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
   if (ierr.ne.NF_NOERR) then
      ! If diagfi not found, try diagfi1
      infile2="diagfi1.nc"
      write(*,*) "Trying file ",trim(infile2)
      ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
      if (ierr.ne.NF_NOERR) then
         write(*,*) "Error: Could not open these files. I'll stop."
         stop ""
      endif
   endif
   ! Get ap and bp (assumed to be in the same file)
   ! ap
   ierr=NF_INQ_VARID(infid2,"ap",tmpvarid)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed to get ap ID"
      stop ""
   endif
   ierr=NF_GET_VAR_REAL(infid2,tmpvarid,ap)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed reading ap"
      stop ""
   endif
   ! bp
   ierr=NF_INQ_VARID(infid2,"bp",tmpvarid)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed to get bp ID"
      stop ""
   endif
   ierr=NF_GET_VAR_REAL(infid2,tmpvarid,bp)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed reading bp"
      stop ""
   endif
   ! Close file
   write(*,*) 'OK, got ap and bp'
   ierr=NF_CLOSE(infid2)
else
   ! Get ap and bp (assumed to be in the same file)
   ! ap (already got ID)
   ierr=NF_GET_VAR_REAL(infid,tmpvarid,ap)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed reading ap"
      stop ""
   endif
   ! bp
   ierr=NF_INQ_VARID(infid,"bp",tmpvarid)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed to get bp ID"
      stop ""
   endif
   ierr=NF_GET_VAR_REAL(infid,tmpvarid,bp)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed reading bp"
      stop ""
   endif
   write(*,*) 'OK, got ap and bp'
endif

! Look for hybrid coordinates aps and bps
allocate(aps(altlength))
allocate(bps(altlength))
ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
if (ierr.ne.NF_NOERR) then
   ! If ap not in input file, it is in a diagfi file
   write(*,*) "Failed to get ap ID from input file ",trim(infile)
   infile2="diagfi.nc"
   write(*,*) "Trying file ",trim(infile2)
   ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
   if (ierr.ne.NF_NOERR) then
      ! If diagfi not found, try diagfi1
      infile2="diagfi1.nc"
      write(*,*) "Trying file ",trim(infile2)
      ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
      if (ierr.ne.NF_NOERR) then
         write(*,*) "Error: Could not open these files. I'll stop."
         stop ""
      endif
   endif
   ! Get aps and bps (assumed to be in the same file)
   ! aps
   ierr=NF_INQ_VARID(infid2,"aps",tmpvarid)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed to get aps ID"
      stop ""
   endif
   ierr=NF_GET_VAR_REAL(infid2,tmpvarid,aps)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed reading aps"
      stop ""
   endif
   ! bps
   ierr=NF_INQ_VARID(infid2,"bps",tmpvarid)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed to get bps ID"
      stop ""
   endif
   ierr=NF_GET_VAR_REAL(infid2,tmpvarid,bps)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed reading bps"
      stop ""
   endif
   ! Close file
   write(*,*) 'OK, got aps and bps'
   ierr=NF_CLOSE(infid2)
else
   ! Get aps and bps (assumed to be in the same file)
   ! aps (already got ID)
   ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed reading aps"
      stop ""
   endif
   ! bps
   ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed to get bps ID"
      stop ""
   endif
   ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed reading bps"
      stop ""
   endif
   write(*,*) 'OK, got aps and bps'
endif

!===============================================================================
! 1.4 Read surface pressure for vertical integration
!===============================================================================

ierr=NF_INQ_VARID(infid,"ps",tmpvarid)
if (ierr.ne.NF_NOERR) then
   write(*,*) "Error: Failed to get ps ID"
   stop ""
else
   allocate(ps(lonlength,latlength,timelength))
   ierr=NF_GET_VAR_REAL(infid,tmpvarid,ps)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed reading ps"
      stop ""
   endif
endif

!===============================================================================
! 1.5 Read atmospheric density to convert concentration to mass
!===============================================================================

foundrho = .false.
ierr=NF_INQ_VARID(infid,"rho",tmpvarid)
if (ierr.ne.NF_NOERR) then
   write(*,*) "Failed to get rho ID"
else
   allocate(rho(lonlength,latlength,altlength,timelength))
   ierr=NF_GET_VAR_REAL(infid,tmpvarid,rho)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Failed to read rho"
   else
      foundrho = .true.
   endif
endif

!===============================================================================
! 1.6 Read/compute or assume MMEAN cst
!     (for vmr and concentration conversion)
!===============================================================================

rgp = 8.3143 ! universal gas constant in J mol-1 K-1
Rmean = 191 ! mean gas constant in J kg-1 K-1

ismmcst = .false.
foundmm = .false.
foundpress = .true.

! 1.5.1 Look for mmean
ierr=NF_INQ_VARID(infid,"mmean",tmpvarid)
if (ierr.ne.NF_NOERR) then ! can't find mean id
   write(*,*) "Failed to get mmean ID."
   write(*,*) "I'll try to compute it from rho, temp and pressure."
   ! allocate locmm for later computation
   allocate(locmm(lonlength,latlength,altlength,timelength))
   ! Look for rho
   if (foundrho) then
      write(*,*) "Previously found rho, looking for temp and pressure"
      ! Look for temp
      ierr=NF_INQ_VARID(infid,"temp",tmpvarid)
      if (ierr.ne.NF_NOERR) then
         write(*,*) "Failed to get temp ID."
         write(*,*) "Can't compute mmean, I'll assume it is constant."
         cstmm = rgp/Rmean
         ismmcst = .true.
      else
         allocate(temp(lonlength,latlength,altlength,timelength))
         ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp)
         if (ierr.ne.NF_NOERR) then
            write(*,*) "Failed to read temp."
            write(*,*) "Can't compute mmean, I'll assume it is constant."
            cstmm = rgp/Rmean
            ismmcst = .true.
         else
            ! Look for pressure
            ierr=NF_INQ_VARID(infid,"pressure",tmpvarid)
            if (ierr.ne.NF_NOERR) then
               write(*,*) "Failed to get pressure ID."
               write(*,*) "I'll compute it from aps, bps and ps"
               ! allocate press for later computation
               allocate(press(lonlength,latlength,altlength,timelength))
               foundpress=.false.
            else
               allocate(press(lonlength,latlength,altlength,timelength))
               ierr=NF_GET_VAR_REAL(infid,tmpvarid,press)
               if (ierr.ne.NF_NOERR) then
                  write(*,*) "Failed to read pressure."
                  write(*,*) "I'll compute it from aps, bps and ps."
                  foundpress=.false.
               endif
            endif ! end look for pressure
         endif
      endif ! end look for temp
   else ! not found rho
      write(*,*) "Failed to read rho."
      write(*,*) "Can't compute mmean, I'll assume it is constant."
      cstmm = rgp/Rmean
      ismmcst = .true.
   endif ! end if found rho
else ! if found mmean id, try to read it
   allocate(locmm(lonlength,latlength,altlength,timelength))
   ierr=NF_GET_VAR_REAL(infid,tmpvarid,locmm)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Failed to read mmean ID."
      write(*,*) "I'll try to compute it from rho, temp and pressure."
      ! Look for rho
      if (foundrho) then
         write(*,*) "Previously found rho, looking for temp and pressure"
         ! Look for temp
         ierr=NF_INQ_VARID(infid,"temp",tmpvarid)
         if (ierr.ne.NF_NOERR) then
            write(*,*) "Failed to get temp ID."
            write(*,*) "Can't compute mmean, I'll assume it is constant."
            cstmm = rgp/Rmean
            ismmcst = .true.
         else
            allocate(temp(lonlength,latlength,altlength,timelength))
            ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp)
            if (ierr.ne.NF_NOERR) then
               write(*,*) "Failed to read temp."
               write(*,*) "Can't compute mmean, I'll assume it is constant."
               cstmm = rgp/Rmean
               ismmcst = .true.
            else
               ! Look for pressure
               ierr=NF_INQ_VARID(infid,"pressure",tmpvarid)
               if (ierr.ne.NF_NOERR) then
                  write(*,*) "Failed to get pressure ID."
                  write(*,*) "I'll compute it from aps, bps and ps"
                  ! allocate press for later computation
                  allocate(press(lonlength,latlength,altlength,timelength))
                  foundpress=.false.
               else
                  allocate(press(lonlength,latlength,altlength,timelength))
                  ierr=NF_GET_VAR_REAL(infid,tmpvarid,press)
                  if (ierr.ne.NF_NOERR) then
                     write(*,*) "Failed to read pressure."
                     write(*,*) "I'll compute it from aps, bps and ps."
                     foundpress=.false.
                  endif
               endif ! end look for pressure
            endif
         endif ! end look for temp
      else ! not found rho
         write(*,*) "Failed to read rho."
         write(*,*) "Can't compute mmean, I'll assume it is constant."
         cstmm = rgp/Rmean
         ismmcst = .true.
      endif ! end look for rho
   else
      foundmm=.true.
   endif
endif ! end look for mmean

! 1.5.2 Compute press and mm if required
if (.not.ismmcst) then
   if (.not.foundmm) then
      if (.not.foundpress) then
         ! Compute press
         do it=1,timelength
            do ialt=1,altlength
               do ilat=1,latlength
                  do ilon=1,lonlength
                     press(ilon,ilat,ialt,it)=&
                        aps(ialt)+bps(ialt)*ps(ilon,ialt,it)
                  enddo
               enddo
            enddo
         enddo
         ! end compute press
      endif ! end if not found press

      ! Compute local atmospheric molecular mass
      do it=1,timelength
         do ialt=1,altlength
            do ilat=1,latlength
               do ilon=1,lonlength
                  locmm(ilon,ilat,ialt,it)=&
                     rgp*rho(ilon,ilat,ialt,it)*&
                     temp(ilon,ilat,ialt,it)/&
                     press(ilon,ilat,ialt,it)
               enddo
            enddo
         enddo
      enddo
      ! end compute local atmospheric molecular mass
   endif ! end if not found mm
endif ! end if mm not cst

!==============================================================================
! 1.6. Create and initialize output file
!==============================================================================

filename=infile(1:len_trim(infile)-3)//"_col.nc"
write(*,*)'The output file has name:'
write(*,*) filename

! Initialize output file's lat,lon,alt and time dimensions
call initiate(filename,lat,lon,alt,time,nout,&
     latdimout,londimout,altdimout,timedimout)
! ! ! Initialize output file's aps,bps variables
! call init2(infid,lonlength,latlength,altlength,&
!      nout,londimout,latdimout,altdimout, interlayerdimout)

!==============================================================================
! 2. PROCESS DATA
!==============================================================================

!This is where you can work on data 

! Constants
g = 3.72
nav = 6.023d23

!=============================================================
allocate(tracer(lonlength,latlength,altlength,timelength))
allocate(tracerCol(lonlength,latlength,timelength))

do i=1,nbvar
   ! Read tracer
   ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed reading ID for ", trim(var(i))
      stop ""
   endif
   ierr=NF_GET_VAR_REAL(infid,tmpvarid,tracer)
   if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Failed reading ", trim(var(i))
      stop ""
   endif
   ! check tracer unit
   ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'units',tmpunit)
   ismmr = (tmpunit.eq.'kg/kg').or.(tmpunit.eq.'kg.kg-1')
   isvmr = (tmpunit.eq.'mol/mol').or.(tmpunit.eq.'mol.mol-1')
   isco = (tmpunit.eq.'cm-3').or.(tmpunit.eq.'m-3')
   ! Compute column in mmr case
   if (ismmr) then
      do ilat=1,latlength
         do ilon=1,lonlength
            do it=1,timelength
               tracerMass = 0.0
               do ialt=1,altlength
                  incMass = ((ap(ialt) - ap(ialt + 1)) + (bp(ialt) - bp(ialt + 1))*ps(ilon,ilat,it))/g
                  tracerMass = tracerMass + incMass*tracer(ilon,ilat,ialt,it)
               enddo
               tracerCol(ilon,ilat,it) = tracerMass
            enddo
         enddo
      enddo
   ! Compute column in vmr case
   else if (isvmr) then
      write(*,*) trim(var(i)), " is a vmr tracer"
      ! Get current specie's molecular mass
      call getmm(var(i),tmpmm,foundmm)
      if (.not.foundmm) then
         write(*,*) "Cannot compute column from vmr variable ", trim(var(i))
         write(*,*) "Check unit, name, or add it to getmm. I'll stop here"
         stop ""
      else
         write(*,*) "Getting molecular mass", tmpmm, " kg/mol"
         do ilat=1,latlength
            do ilon=1,lonlength
               do it=1,timelength
                  tracerMass = 0.0
                  do ialt=1,altlength
                     incMass = ((ap(ialt) - ap(ialt + 1)) + (bp(ialt) - bp(ialt + 1))*ps(ilon,ilat,it))/g
                     ! Handle cst mm case
                     if (ismmcst) then
                        tracerMass = tracerMass + incMass*tracer(ilon,ilat,ialt,it)*&
                           tmpmm/cstmm
                     ! Handle local mm case
                     else
                        ! mmean is stored in g/mol so we multiply by 1.e3
                        tracerMass = tracerMass + incMass*tracer(ilon,ilat,ialt,it)*&
                           tmpmm/locmm(ilon,ilat,ialt,it)*1.e3 
                     endif
                  enddo
                  tracerCol(ilon,ilat,it) = tracerMass
               enddo
            enddo
         enddo
      endif
   ! Compute column in concentration case
   else if(isco) then   
      write(*,*) trim(var(i)), " is a num tracer"
      write(*,*) "Checking rho is available..."
      if (.not.foundrho) then
         write(*,*) "Error: rho not available. Can't proceed"
         stop ""
      endif
      ! Get current specie's molecular mass
      call getmm(var(i),tmpmm,foundmm)
      if (.not.foundmm) then
         write(*,*) "Cannot compute column from concentration variable ", trim(var(i))
         write(*,*) "Check unit, name, or add it to getmm. I'll stop here"
         stop ""
      else
         write(*,*) "Getting molecular mass", tmpmm, " kg/mol"
         do ilat=1,latlength
            do ilon=1,lonlength
               do it=1,timelength
                  tracerMass = 0.0
                  do ialt=1,altlength
                     incMass = ((ap(ialt) - ap(ialt + 1)) + (bp(ialt) - bp(ialt + 1))*ps(ilon,ilat,it))/g
                     tracerMass = tracerMass + incMass/rho(ilon,ilat,ialt,it)*tracer(ilon,ilat,ialt,it)*tmpmm/nav
                  enddo
                  ! Is it in S.I. ?
                  if (tmpunit.eq.'m-3') then
                     tracerCol(ilon,ilat,it) = tracerMass
                  else
                     tracerCol(ilon,ilat,it) = tracerMass*1.d6
                  endif
               enddo
            enddo
         enddo
      endif
   ! If none of previous case matches, skip
   else
      write(*,*) trim(var(i)), " is neither mmr, vmr or concentration"
      write(*,*) "Will skip"
   endif ! end of ismmr, isvmr, isco

   ! Define variable to write
   tmpvarname = trim(var(i))//"_col"
   call def_var(nout,tmpvarname,tmpvarname,"kg/m2",3,&
               (/londimout,latdimout,timedimout/),tracvarout,ierr)
   if (ierr.ne.NF_NOERR) then
      write(*,*) 'Error, could not define variable ', trim(var(i))
      stop ""
   endif
   ! Write in the output file
   ierr=NF_PUT_VAR_REAL(nout,tracvarout,tracerCol)
   if (ierr.ne.NF_NOERR) then
      write(*,*)'Error, Failed to write variable ', trim(var(i))
      stop ""
   endif
enddo ! end of do nbvar


!----------------------------------
! Close input file
ierr=NF_CLOSE(infid)

! Close output file
ierr=NF_CLOSE(nout)

contains

!******************************************************************************
Subroutine initiate (filename,lat,lon,alt,time,&
                     nout,latdimout,londimout,altdimout,timedimout)
!==============================================================================
! Purpose:
! Create and initialize a data file (NetCDF format)
! Include grid and time (lon,lat,alt,time)
!==============================================================================
! Remarks:
! The NetCDF file (created in this subroutine) remains open
!==============================================================================

implicit none

include "netcdf.inc" ! NetCDF definitions

!==============================================================================
! Arguments:
!==============================================================================
character (len=*), intent(in):: filename
! filename(): the file's name
real, dimension(:), intent(in):: lat
! lat(): latitude
real, dimension(:), intent(in):: lon
! lon(): longitude
real, dimension(:), intent(in):: alt
! alt(): altitude
real, dimension(:), intent(in):: time
! time(): Time
integer, intent(out):: nout
! nout: [netcdf] file ID
integer, intent(out):: latdimout
! latdimout: [netcdf] lat() (i.e.: latitude)  ID
integer, intent(out):: londimout
! londimout: [netcdf] lon()  ID
integer, intent(out):: altdimout
! altdimout: [netcdf] alt()  ID
integer, intent(out):: timedimout
! timedimout: [netcdf] "Time"  ID

!==============================================================================
! Local variables:
!==============================================================================
integer :: nvarid,ierr
! nvarid: [netcdf] ID of a variable
! ierr: [netcdf]  return error code (from called subroutines)

!==============================================================================
! 1. Create (and open) output file
!==============================================================================
write(*,*) "creating "//trim(adjustl(filename))//'...'
ierr = NF_CREATE(filename,NF_CLOBBER,nout)
! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
if (ierr.NE.NF_NOERR) then
   WRITE(*,*)'ERROR: Impossible to create the file.'
   stop ""
endif

!==============================================================================
! 2. Define/write "dimensions" and get their IDs
!==============================================================================

ierr = NF_DEF_DIM(nout, "latitude", size(lat), latdimout)
ierr = NF_DEF_DIM(nout, "longitude", size(lon), londimout)
ierr = NF_DEF_DIM(nout, "altitude", size(alt), altdimout)
ierr = NF_DEF_DIM(nout, "Time", size(time), timedimout)

! End netcdf define mode
ierr = NF_ENDDEF(nout)

!==============================================================================
! 3. Write "Time" (attributes)
!==============================================================================

call def_var(nout,"Time","Time","days since 0000-00-0 00:00:00",1,&
             (/timedimout/),nvarid,ierr)

#ifdef NC_DOUBLE
ierr = NF_PUT_VAR_DOUBLE(nout,nvarid,time)
#else
ierr = NF_PUT_VAR_REAL(nout,nvarid,time)
#endif

!==============================================================================
! 4. Write "latitude" (data and attributes)
!==============================================================================

call def_var(nout,"latitude","latitude","degrees_north",1,&
             (/latdimout/),nvarid,ierr)

#ifdef NC_DOUBLE
ierr = NF_PUT_VAR_DOUBLE(nout,nvarid,lat)
#else
ierr = NF_PUT_VAR_REAL(nout,nvarid,lat)
#endif

!==============================================================================
! 4. Write "longitude" (data and attributes)
!==============================================================================

call def_var(nout,"longitude","East longitude","degrees_east",1,&
             (/londimout/),nvarid,ierr)

#ifdef NC_DOUBLE
ierr = NF_PUT_VAR_DOUBLE(nout,nvarid,lon)
#else
ierr = NF_PUT_VAR_REAL(nout,nvarid,lon)
#endif

!==============================================================================
! 4. Write "altitude" (data and attributes)
!==============================================================================

call def_var(nout,"altitude","Altitude","km",1,&
             (/altdimout/),nvarid,ierr)

#ifdef NC_DOUBLE
ierr = NF_PUT_VAR_DOUBLE(nout,nvarid,alt)
#else
ierr = NF_PUT_VAR_REAL(nout,nvarid,alt)
#endif 

end Subroutine initiate

subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
!==============================================================================
! Purpose: Write a variable (i.e: add a variable to a dataset)
! called "name"; along with its attributes "title", "units"...
! to a file (following the NetCDF format)
!==============================================================================
! Remarks:
! The NetCDF file must be open
!==============================================================================

implicit none

include "netcdf.inc" ! NetCDF definitions

!==============================================================================
! Arguments:
!==============================================================================
integer, intent(in) :: nid
! nid: [netcdf] file ID #
character (len=*), intent(in) :: name
! name(): [netcdf] variable's name
character (len=*), intent(in) :: title
! title(): [netcdf] variable's "title" attribute
character (len=*), intent(in) :: units
! unit(): [netcdf] variable's "units" attribute
integer, intent(in) :: nbdim
! nbdim: number of dimensions of the variable
integer, dimension(nbdim), intent(in) :: dim
! dim(nbdim): [netcdf] dimension(s) ID(s)
integer, intent(out) :: nvarid
! nvarid: [netcdf] ID # of the variable
integer, intent(out) :: ierr
! ierr: [netcdf] subroutines returned error code

! Switch to netcdf define mode
ierr=NF_REDEF(nid)

! Insert the definition of the variable
#ifdef NC_DOUBLE
ierr=NF_DEF_VAR(nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
#else
ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
#endif

! Write the attributes
ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",len_trim(adjustl(title)),adjustl(title))
ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))

! End netcdf define mode
ierr=NF_ENDDEF(nid)

end subroutine def_var

!******************************************************************************
subroutine getmm(invarname, tmpmm, foundmm)
!==============================================================================
! Purpose:
! Get the molar mass tmpmm coresponding to a given varname
!==============================================================================
! Remarks:
! - varname is assumed to have the format *_name_* or *_name
!   (based on the usual format vmr_..., num_...)
! - Feel free to add any missing specie
!==============================================================================

implicit none

!==============================================================================
! Arguments:
!==============================================================================
character(len=*), intent(in) :: invarname ! name of the input variable
real, intent(out) :: tmpmm ! molar mass to output
logical, intent(out) :: foundmm ! did we find a correspondance ?

!==============================================================================
! Local variables
!==============================================================================
integer :: first, last, pos, length, count
! first: beginning of specie name, last: end of specie name
! pos: cursor position in string varname. length: varname length
! count: number of encountered separator "_" (2 at most)
character(len=len_trim(invarname)) :: varname ! name of the specie

!==============================================================================
! Find specie name
!==============================================================================

length = len_trim(invarname)
first = 1
last = length
pos=1
count=0
! Cursor goes through varname
do while (pos < length)
   if ((invarname(pos:pos).eq.'_').and.(count.eq.0)) then
      ! First "_" locates the beginning of the specie's name
      count = count + 1
      first = pos + 1
   else if ((invarname(pos:pos).eq.'_').and.(count.eq.1)) then
      ! Second "_" locates the end (default is length)
      count = count + 1
      last = pos - 1
   endif
   pos = pos + 1
enddo
varname = invarname(first:last)

!==============================================================================
! Check name among all possibilities
!==============================================================================

foundmm=.false.

if (varname.eq."co2") then
   tmpmm=44.
   foundmm=.true.
endif
if (varname.eq."co") then
   tmpmm=28.
   foundmm=.true.
endif
if (varname.eq."o") then
   tmpmm=16.
   foundmm=.true.
endif
if (varname.eq."o1d") then
   tmpmm=16.
   foundmm=.true.
endif
if (varname.eq."o2") then
   tmpmm=32.
   foundmm=.true.
endif
if (varname.eq."o3") then
   tmpmm=48.
   foundmm=.true.
endif
if (varname.eq."h") then
   tmpmm=1.
   foundmm=.true.
endif
if (varname.eq."h2") then
   tmpmm=2.
   foundmm=.true.
endif
if (varname.eq."oh") then
   tmpmm=17.
   foundmm=.true.
endif
if (varname.eq."ho2") then
   tmpmm=33.
   foundmm=.true.
endif
if (varname.eq."h2o2") then
   tmpmm=34.
   foundmm=.true.
endif
if (varname.eq."n2") then
   tmpmm=28.
   foundmm=.true.
endif
if (varname.eq."ch4") then
   tmpmm=16.
   foundmm=.true.
endif
if (varname.eq."ar") then
   tmpmm=40.
   foundmm=.true.
endif
if (varname.eq."n") then
   tmpmm=14.
   foundmm=.true.
endif
if (varname.eq."no") then
   tmpmm=30.
   foundmm=.true.
endif
if (varname.eq."no2") then
   tmpmm=46.
   foundmm=.true.
endif
if (varname.eq."n2d") then
   tmpmm=28.
   foundmm=.true.
endif
if (varname.eq."he") then
   tmpmm=4.
   foundmm=.true.
endif
if (varname.eq."co2plus") then
   tmpmm=44.
   foundmm=.true.
endif
if (varname.eq."oplus") then
   tmpmm=16.
   foundmm=.true.
endif
if (varname.eq."o2plus") then
   tmpmm=32.
   foundmm=.true.
endif
if (varname.eq."coplus") then
   tmpmm=28.
   foundmm=.true.
endif
if (varname.eq."cplus") then
   tmpmm=12.
   foundmm=.true.
endif
if (varname.eq."nplus") then
   tmpmm=14.
   foundmm=.true.
endif
if (varname.eq."noplus") then
   tmpmm=30.
   foundmm=.true.
endif
if (varname.eq."n2plus") then
   tmpmm=28.
   foundmm=.true.
endif
if (varname.eq."hplus") then
   tmpmm=1.
   foundmm=.true.
endif
if (varname.eq."hco2plus") then
   tmpmm=45.
   foundmm=.true.
endif
if (varname.eq."hcoplus") then
   tmpmm=29.
   foundmm=.true.
endif
if (varname.eq."h2plus") then
   tmpmm=18.
   foundmm=.true.
endif
if (varname.eq."h3oplus") then
   tmpmm=19.
   foundmm=.true.
endif
if (varname.eq."ohplus") then
   tmpmm=17.
   foundmm=.true.
endif
if (varname.eq."elec") then
   tmpmm=1./1822.89
   foundmm=.true.
endif
if (varname.eq."h2ovap") then ! tmp
   tmpmm=18.
   foundmm=.true.
endif
if (varname.eq."h2o") then ! tmp
   tmpmm=18.
   foundmm=.true.
endif
if (varname.eq."hdovap") then
   tmpmm=19.
   foundmm=.true.
endif
if (varname.eq."od") then
   tmpmm=18.
   foundmm=.true.
endif
if (varname.eq."d") then
   tmpmm=2.
   foundmm=.true.
endif
if (varname.eq."hd") then
   tmpmm=3.
   foundmm=.true.
endif
if (varname.eq."do2") then
   tmpmm=34.
   foundmm=.true.
endif
if (varname.eq."hdo2") then
   tmpmm=35.
   foundmm=.true.
endif
if (varname.eq."co2ice") then
   tmpmm=44.
   foundmm=.true.
endif
if (varname.eq."h2oice") then
   tmpmm=18.
   foundmm=.true.
endif
if (varname.eq."hdoice") then
   tmpmm=19.
   foundmm=.true.
endif
if (varname.eq."Ar_N2") then
   tmpmm=30.
   foundmm=.true.
endif        

if (foundmm) then
   ! Back in kg
   tmpmm = tmpmm/1.e3
else
   write(*,*) "Found no molar mass corresponding to variable ", trim(varname)
endif

end subroutine getmm

end program