program zrecast

! This program reads 4D (lon-lat-alt-time) fields from GCM output files
! (ie: diagfi.nc time series or concat.nc or stats.nc files) and, by
! integrating the hydrostatic equation, recasts data along the vertical
! direction.
! The vertical coordinate can be either "above areoid altitudes" or
! "pressure". Some interpolation along the vertical direction is also
! done, following instructions given by user.
! For "above areoid altitudes" output, Atmospheric pressure is added to
! output dataset; for "pressure coordinate" outputs, the above areoid
! altitude of pressure is added to output dataset.
!
! Minimal requirements and dependencies:
! The dataset must include the following data:
! - surface pressure
! - atmospheric temperature
! - hybrid coordinates aps() and bps(), or sigma levels() (see section 1.3.2)
! - ground geopotential (in input file; if not found, it is sought
!   in a 'diagfi.nc' file. If not found there, it is then sought in
!   a 'phisinit.nc' file  (see section 1.3.3 of program)
!
! - When integration the hydrostatic equation, we assume that R, the molecular
!   Gas Constant, may not be constant, so it is computed as
!      R=P/(rho*T) (P=Pressure, rho=density, T=temperature)
!   If 'rho' is not available, then we use a constant R (see section 2.2)
!
! WARNING: Asking for many points along the vertical direction quickly
!          leads to HUGE output files.
!
! EM 09/2006
! EM 10/2006 : Modified program so that it can now process 'stats.nc'
!              files obtained from British GCM (ie: vertical coordinate
!              given as sigma levels and geopotential read from file
!              'phisinit.nc')
! EM 01/2006 : Corrected a bug in vertical (log) interpolation for pressure
!              and density

implicit none

include "netcdf.inc" ! NetCDF definitions

character (len=128) :: infile ! input file name (diagfi.nc or stats.nc format)
character (len=128) :: infile2 ! second input file (may be needed for 'phisini')
character (len=128) :: outfile ! output file name

character (len=64) :: text ! to store some text
character (len=64) :: tmpvarname ! temporarily store a variable name
integer tmpvarid ! temporarily store a variable ID
integer tmpdimid ! temporarily store a dimension ID
integer tmpndims ! temporarily store # of dimensions of a variable
integer infid ! NetCDF input file ID (of diagfi.nc or stats.nc format)
integer infid2 ! NetCDF input file which contains 'phisini' dataset (diagfi.nc)
integer nbvarinfile ! # of variables in input file
integer nbattr ! # of attributes of a given variable in input file
integer nbvar4dinfile ! # of 4D (lon,lat,alt,time) variables in input file
integer outfid ! NetCDF output file ID
integer lon_dimid,lat_dimid,alt_dimid,time_dimid ! NetCDF dimension IDs
integer lon_varid,lat_varid,alt_varid,time_varid
integer za_varid,p_varid ! above areoid and pressure data IDs
integer,dimension(4) :: datashape

real :: miss_val=-9.99e+33 ! special "missing value" to specify missing data
real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
character (len=64), dimension(:), allocatable :: var
! var(): names of variables that will be processed
integer nbvar ! # of variables to process
integer,dimension(:),allocatable :: var_id ! IDs of variables var() (in outfile)
real,dimension(:),allocatable :: lon ! longitude
integer lonlength ! # of grid points along longitude
real,dimension(:),allocatable :: lat ! latitude
integer latlength ! # of grid points along latitude
integer altlength ! # of grid point along altitude (of input datasets)
real,dimension(:),allocatable :: time ! time
integer timelength ! # of points along time
real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
real,dimension(:),allocatable :: sigma ! sigma levels
real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
real,dimension(:,:,:),allocatable :: ps ! GCM surface pressure
real,dimension(:,:,:,:),allocatable :: press ! GCM atmospheric pressure 
real,dimension(:,:,:,:),allocatable :: temp ! GCM atmospheric temperature
real,dimension(:,:,:,:),allocatable :: rho ! GCM atmospheric density
real,dimension(:,:,:,:),allocatable :: za_gcm ! GCM above areoid levels
real,dimension(:,:,:,:),allocatable :: indata ! to store a GCM 4D dataset
real,dimension(:,:,:,:),allocatable :: outdata ! to store a 4D dataset
integer ierr ! NetCDF routines return code
integer i,j,ilon,ilat,ilev,itim ! for loops

integer ztype ! Flag for vertical coordinate of output
! ztype=1: pressure  ztype=2: above areoid
integer nblev ! # of levels (along vertical coordinate) for output data
real pmin,pmax ! min and max values for output pressure coordinate
real,dimension(:),allocatable :: plevel ! Pressure levels for output data
real zamin,zamax ! min and max values for output above areoid coordinate
real,dimension(:),allocatable :: zareoid ! Above areoid heights for output data
logical :: have_rho ! Flag: true if density 'rho' is available
logical :: have_sigma ! Flag: true if sigma levels are known (false if hybrid
                      !       coordinates are used)
!===============================================================================
! 1. Input parameters
!===============================================================================

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

write(*,*) ""
write(*,*) " Program valid for diagfi.nc, concatnc.nc and stats.nc files"
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 file ',trim(infile)
   stop ""
endif

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

ierr=NF_INQ_NVARS(infid,nbvarinfile)
if (ierr.ne.NF_NOERR) then
  write(*,*) 'ERROR: Failed geting number of variables from file'
  stop
endif

write(*,*)" The following variables have been found:"
nbvar4dinfile=0
do i=1,nbvarinfile
  ! get name of variable # i
  ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
  ! check if it is a 4D variable
  ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
  if (tmpndims.eq.4) then
    nbvar4dinfile=nbvar4dinfile+1
    write(*,*) trim(tmpvarname)
  endif
enddo

allocate(var(nbvar4dinfile))

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
    nbvar=nbvar+1
    var(nbvar)=tmpvarname
  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 variables
    ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
    if (tmpndims.eq.4) 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() (or sigma levels sigma())
!     and phisinit() from input file
!===============================================================================

! 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
        stop "Error: Failed reading latitude"
      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
        stop "Error: Failed reading longitude"
      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
        stop "Error: Failed reading Time"
      endif
    endif
  endif
endif

! altlength
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_DIMLEN(infid,tmpdimid,altlength)
  if (ierr.ne.NF_NOERR) then
      stop "Error: Failed to get altitude length"
  endif
endif

! 1.3.2 Get hybrid coordinates (or sigma levels)

! start by looking for sigma levels
ierr=NF_INQ_VARID(infid,"sigma",tmpvarid)
if (ierr.ne.NF_NOERR) then
  have_sigma=.false.
  write(*,*) "Could not find sigma levels... will look for hybrid coordinates"
else
  have_sigma=.true.
  allocate(sigma(altlength))
  ierr=NF_GET_VAR_REAL(infid,tmpvarid,sigma)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed reading sigma"
  endif
endif

! if no sigma levels, look for hybrid coordinates
if (.not.have_sigma) then
  ! hybrid coordinate aps
  ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed to get aps ID"
  else
    allocate(aps(altlength))
    ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
    if (ierr.ne.NF_NOERR) then
      stop "Error: Failed reading aps"
    endif
  endif

  ! hybrid coordinate bps
  ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed to get bps ID"
  else
    allocate(bps(altlength))
    ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
    if (ierr.ne.NF_NOERR) then
      stop "Error: Failed reading bps"
    endif
  endif
endif !of if (.not.have_sigma)

! 1.3.3 ground geopotential phisinit

allocate(phisinit(lonlength,latlength))
! look for 'phisinit' in current file
ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid) 
if (ierr.ne.NF_NOERR) then
  write(*,*) "Warning: Failed to get phisinit ID from file ",trim(infile)
  infile2="diagfi.nc"
  write(*,*) "         Trying file ",trim(infile2)
  ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
  if (ierr.ne.NF_NOERR) then
    write(*,*) "Problem: Could not find/open that file"
    infile2="phisinit.nc"
    write(*,*) "         Trying file ",trim(infile2)
    ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
    if (ierr.ne.NF_NOERR) then
      write(*,*) "Error: Could not open that file either"
      stop "Might as well stop here"
    endif
  endif

  ! Get ID for phisinit
  ierr=NF_INQ_VARID(infid2,"phisinit",tmpvarid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed to get phisinit ID"
  endif
  ! Get physinit
  ierr=NF_GET_VAR_REAL(infid2,tmpvarid,phisinit)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed reading phisinit"
  endif
  ! Close file
  write(*,*) 'OK, got phisinit'
  ierr=NF_CLOSE(infid2)
else
  ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed reading phisinit"
  endif
endif

!===============================================================================
! 1.4 Choose and build the new vertical coordinate
!===============================================================================

write(*,*) ""
write(*,*) "Which vertical coordinate should the output be in?"
ztype=0
do while ((ztype.ne.1).and.(ztype.ne.2))
  write(*,*) "(1: pressure, 2: above areoid altitude)"
  read(*,*)ztype
enddo

! ask for # of points and end values
write(*,*) ""
write(*,*) "Enter min and max of vertical coordinate (Pa or m)"
if (ztype.eq.1) then ! pressure coordinate
  read(*,*) pmin,pmax
else  ! above areoid coordinate
  read(*,*) zamin,zamax
endif
write(*,*) "Number of levels along vertical coordinate?"
read(*,*) nblev

! Build corresponding vertical coordinates
if (ztype.eq.1) then ! pressure coordinate
  allocate(plevel(nblev))
  if (nblev.eq.1) then ! in case only one level is asked for
    plevel(nblev)=pmin
  else
    do i=1,nblev
      !    build exponentially spread layers
      plevel(i)=exp(log(pmax)+(log(pmin)-log(pmax))* &
                    ((real(i)-1.0)/(real(nblev)-1.0)))
    enddo
  endif
else ! above areoid heights
  allocate(zareoid(nblev))
  if (nblev.eq.1) then ! in case only one level is asked for
    zareoid(nblev)=zamin
  else
    do i=1,nblev
      zareoid(i)=zamin+(real(i)-1.0)*((zamax-zamin)/(real(nblev)-1.0))
    enddo
  endif
endif

!===============================================================================
! 1.5 Get output file name
!===============================================================================
write(*,*) ""
write(*,*) "Enter output file name"
read(*,*) outfile
! Francois:
!if (ztype.eq.1) then ! pressure coordinate
!  outfile=infile(1:len_trim(infile)-3)//"_P.nc"
!else  ! above areoid coordinate
!  outfile=infile(1:len_trim(infile)-3)//"_A.nc"
!endif


!===============================================================================
! 2.1 Build/store GCM fields which will be used later
!===============================================================================

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

!===============================================================================
! 2.1.2 Atmospheric pressure
!===============================================================================
allocate(press(lonlength,latlength,altlength,timelength))

if (have_sigma) then ! sigma coordinate
  do itim=1,timelength
    do ilev=1,altlength
      do ilat=1,latlength
        do ilon=1,lonlength
          press(ilon,ilat,ilev,itim)=sigma(ilev)*ps(ilon,ilat,itim)
        enddo
      enddo
    enddo
  enddo
else ! hybrid coordinates
  do itim=1,timelength
    do ilev=1,altlength
      do ilat=1,latlength
        do ilon=1,lonlength
          press(ilon,ilat,ilev,itim)=aps(ilev)+bps(ilev)*ps(ilon,ilat,itim)
        enddo
      enddo
    enddo
  enddo
endif

!===============================================================================
! 2.1.3 Atmospheric temperature
!===============================================================================
allocate(temp(lonlength,latlength,altlength,timelength))

ierr=NF_INQ_VARID(infid,"temp",tmpvarid)
if (ierr.ne.NF_NOERR) then
  ! stop "Error: Failed to get temp ID"
  ! try "t" for temperature
  ierr=NF_INQ_VARID(infid,"t",tmpvarid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed to get t ID"
  else
    ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp)
    if (ierr.ne.NF_NOERR) then
      stop "Error: Failed reading atmospheric temperature"
    endif
  endif
else
  ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed reading atmospheric temperature"
  endif
endif

!===============================================================================
! 2.1.4 Atmospheric density
!===============================================================================

ierr=NF_INQ_VARID(infid,"rho",tmpvarid)
if (ierr.ne.NF_NOERR) then
  write(*,*) "Warning: Failed to get rho ID"
  have_rho=.false.
else
  have_rho=.true.
  allocate(rho(lonlength,latlength,altlength,timelength))
  ierr=NF_GET_VAR_REAL(infid,tmpvarid,rho)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed reading atmospheric density"
  endif
endif

!===============================================================================
! 2.2 Build GCM Above areoid altitudes of GCM nodes
!===============================================================================
allocate(za_gcm(lonlength,latlength,altlength,timelength))

if (have_rho) then
  call build_gcm_alt(lonlength,latlength,altlength,timelength, &
                     phisinit,ps,press,temp,rho,za_gcm)
else
  write(*,*)"Warning: Using constant R to integrate hydrostatic equation"
  call crude_gcm_alt(lonlength,latlength,altlength,timelength, &
                     phisinit,ps,press,temp,za_gcm)
endif

!===============================================================================
! 3. Create output file and initialize definitions of variables and dimensions
!===============================================================================

!===============================================================================
! 3.1. Output file
!===============================================================================

! Create output file
ierr=NF_CREATE(outfile,NF_CLOBBER,outfid)
if (ierr.ne.NF_NOERR) then
  write(*,*)"Error: could not create file ",outfile
  stop
endif

!===============================================================================
! 3.2. Define dimensions
!===============================================================================
! longitude
ierr=NF_DEF_DIM(outfid,"longitude",lonlength,lon_dimid)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not define longitude dimension"
endif

! latitude
ierr=NF_DEF_DIM(outfid,"latitude",latlength,lat_dimid)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not define latitude dimension"
endif

! altitude
ierr=NF_DEF_DIM(outfid,"altitude",nblev,alt_dimid)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not define altitude dimension"
endif

! time
ierr=NF_DEF_DIM(outfid,"Time",timelength,time_dimid)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not define latitude dimension"
endif

!===============================================================================
! 3.3. Define variables and their attributes
!===============================================================================
! longitude 
datashape(1)=lon_dimid
ierr=NF_DEF_VAR(outfid,"longitude",NF_REAL,1,datashape(1),lon_varid)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not define longitude variable"
endif

! longitude attributes
text='east longitude'
ierr=NF_PUT_ATT_TEXT(outfid,lon_varid,'long_name',len_trim(text),text)
if (ierr.ne.NF_NOERR) then
  stop "Error: Problem writing long_name for longitude"
endif
text='degrees_east'
ierr=NF_PUT_ATT_TEXT(outfid,lon_varid,'units',len_trim(text),text)
if (ierr.ne.NF_NOERR) then
  stop "Error: Problem writing units for longitude"
endif

! latitude
datashape(2)=lat_dimid
ierr=NF_DEF_VAR(outfid,"latitude",NF_REAL,1,datashape(2),lat_varid)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not define latitude variable"
endif

! latitude attributes
text='north latitude'
ierr=NF_PUT_ATT_TEXT(outfid,lat_varid,'long_name',len_trim(text),text)
if (ierr.ne.NF_NOERR) then
  stop "Error: Problem writing long_name for latitude"
endif
text='degrees_north'
ierr=NF_PUT_ATT_TEXT(outfid,lat_varid,'units',len_trim(text),text)
if (ierr.ne.NF_NOERR) then
  stop "Error: Problem writing units for latitude"
endif

! altitude
datashape(3)=alt_dimid
ierr=NF_DEF_VAR(outfid,"altitude",NF_REAL,1,datashape(3),alt_varid)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not define altitude variable"
endif

!altitude attributes
if (ztype.eq.1) then ! pressure vertical coordinate
  text='Pressure levels'
  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'long_name',len_trim(text),text)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing long_name for altitude"
  endif
  text='Pa'
  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'units',len_trim(text),text)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing units for altitude"
  endif
  text='down'
  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'positive',len_trim(text),text)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing positive for altitude"
  endif
else ! above areoid vertical coordinate
  text='Altitude above areoid'
  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'long_name',len_trim(text),text)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing long_name for altitude"
  endif
  text='m'
  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'units',len_trim(text),text)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing units for altitude"
  endif
  text='up'
  ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'positive',len_trim(text),text)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing positive for altitude"
  endif
endif

! time
datashape(4)=time_dimid
ierr=NF_DEF_VAR(outfid,"Time",NF_REAL,1,datashape(4),time_varid)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not define Time variable"
endif

! time attributes
text='Time'
ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'long_name',len_trim(text),text)
if (ierr.ne.NF_NOERR) then
  stop "Error: Problem writing long_name for Time"
endif
text='days since 0000-01-1 00:00:00'
ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'units',len_trim(text),text)
if (ierr.ne.NF_NOERR) then
  stop "Error: Problem writing units for Time"
endif

! add pressure or zareoid
if (ztype.eq.1) then ! pressure vertical coordinate
  ! zareoid dataset
  ierr=NF_DEF_VAR(outfid,"zareoid",NF_REAL,4,datashape,za_varid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Could not define zareoid variable"
  endif
  ! zareoid attributes
  text='altitude above areoid'
  ierr=NF_PUT_ATT_TEXT(outfid,za_varid,'long_name',len_trim(text),text)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing long_name for zareoid"
  endif
  text='m'
  ierr=NF_PUT_ATT_TEXT(outfid,za_varid,'units',len_trim(text),text)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing units for zareoid"
  endif
  ! zareoid missing value
  ierr=NF_PUT_ATT_REAL(outfid,za_varid,'missing_value',NF_REAL,1,miss_val)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing missing_value for zareoid"
  endif
else ! above areoid vertical coordinate
  ! pressure dataset
  ierr=NF_DEF_VAR(outfid,"pressure",NF_REAL,4,datashape,p_varid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Could not define pressure variable"
  endif
  ! pressure attributes
  text='Atmospheric pressure'
  ierr=NF_PUT_ATT_TEXT(outfid,p_varid,'long_name',len_trim(text),text)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing long_name for pressure"
  endif
  text='Pa'
  ierr=NF_PUT_ATT_TEXT(outfid,p_varid,'units',len_trim(text),text)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing units for pressure"
  endif
  ! pressure missing value
  ierr=NF_PUT_ATT_REAL(outfid,p_varid,'missing_value',NF_REAL,1,miss_val)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing missing_value for pressure"
  endif
endif

! Handle 4D variables
allocate(var_id(nbvar))
do i=1,nbvar
  write(*,*) ""
  write(*,*) "Processing ",trim(var(i))
  ! define the variable
  ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,4,datashape,var_id(i))
  if (ierr.ne.NF_NOERR) then
    write(*,*) 'Error, could not define variable ',trim(var(i))
    stop ""
  endif

  ! Get the (input file) ID for the variable and
  ! the # of attributes associated to that variable
  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
  if (ierr.ne.NF_NOERR) then
    write(*,*) 'Error, failed to get ID for variable ',trim(var(i))
    stop ""
  endif
  ierr=NF_INQ_VARNATTS(infid,tmpvarid,nbattr)
  if (ierr.ne.NF_NOERR) then
    write(*,*) 'Error, could not get number of attributes for variable ',&
               trim(var(i))
    stop ""
  endif
  ! inititialize j == number of attributes written to output
  j=0

  ! look for a "long_name" attribute
  text='   '
  ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'long_name',text)
  if (ierr.ne.NF_NOERR) then ! no long_name attribute
    ! try to find an equivalent 'title' attribute
    text='   '
    ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'title',text)
    if (ierr.eq.NF_NOERR) then ! found 'title' attribute
      write(*,*) "Found title ",trim(text)
      j=j+1
      ! write it as a 'long_name' attribute
      ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'long_name',len_trim(text),text)
      if (ierr.ne.NF_NOERR) then
        write(*,*) "Error failed to copy title attribute:",trim(text)
      stop ""
      endif
    endif
  else ! found long_name; write it to outfile
    write(*,*) "Found long_name ",trim(text)
    j=j+1
    ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'long_name',len_trim(text),text)
    if (ierr.ne.NF_NOERR) then
      write(*,*) "Error failed to copy long_name attribute:",trim(text)
      stop""
    endif
  endif
  
  ! look for a "units" attribute
  text='   '
  ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'units',text)
  if (ierr.eq.NF_NOERR) then ! found 'units' attribute
    write(*,*) "Found units ",trim(text)
    j=j+1
    ! write it to output
    ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'units',len_trim(text),text)
    if (ierr.ne.NF_NOERR) then
      write(*,*) "Error failed to copy units attribute:",trim(text)
      stop""
    endif
  endif
  
  ! look for a "missing_value" attribute
  ierr=NF_GET_ATT_REAL(infid,tmpvarid,"missing_value",miss_val)
  if (ierr.eq.NF_NOERR) then ! found 'missing_value' attribute
    write(*,*) "Found missing_value ",miss_val
    j=j+1
  else ! no 'missing_value' attribute, set miss_val to default
    miss_val=miss_val_def
  endif
  
  ! write the missing_value attribute to output
  ierr=NF_PUT_ATT_REAL(outfid,var_id(i),'missing_value',NF_REAL,1,miss_val)
  if (ierr.ne.NF_NOERR) then
    stop "Error, failed to write missing_value attribute"
  endif

  ! warn if some attributes were missed
  if (j.ne.nbattr) then
    write(*,*)'Warning, it seems some attributes of variable ',trim(var(i))
    write(*,*)"were not transfered to the new file"
    write(*,*)'nbattr:',nbattr,' j:',j
  endif
    
enddo ! of do=1,nbvar


!===============================================================================
! 3.4. Write dimensions
!===============================================================================
! Switch out of NetCDF define mode
ierr=NF_ENDDEF(outfid)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not switch out of define mode"
endif

! Write longitude
ierr=NF_PUT_VAR_REAL(outfid,lon_varid,lon)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not write longitude data to output file"
endif

! Write latitude
ierr=NF_PUT_VAR_REAL(outfid,lat_varid,lat)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not write latitude data to output file"
endif

! Write altitude
if (ztype.eq.1) then ! pressure vertical coordinate
  ierr=NF_PUT_VAR_REAL(outfid,alt_varid,plevel)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Could not write altitude data to output file"
  endif
else ! above areoid altitude
  ierr=NF_PUT_VAR_REAL(outfid,alt_varid,zareoid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Could not write altitude data to output file"
  endif
endif

! Write time
ierr=NF_PUT_VAR_REAL(outfid,time_varid,time)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not write Time data to output file"
endif

!===============================================================================
! 4. Interpolate and write variables
!===============================================================================

! 4.0 Allocations
!indata() to store input (on GCM grid) data
allocate(indata(lonlength,latlength,altlength,timelength))
! outdata() to store output (on new vertical grid) data
allocate(outdata(lonlength,latlength,nblev,timelength))

! 4.1 If output is in pressure coordinate
if (ztype.eq.1) then
  do i=1,nbvar ! loop on 4D variable to process
  ! identify and read a dataset
  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
  if (ierr.ne.NF_NOERR) then
    write(*,*) 'Error, failed to get ID for variable ',var(i)
    stop
  endif
  ierr=NF_GET_VAR_REAL(infid,tmpvarid,indata)
  if (ierr.ne.NF_NOERR) then
    write(*,*) 'Error, failed to load variable ',var(i)
    stop
  endif
  
  ! interpolate dataset onto new grid
  call p_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
                      miss_val,ps,press,indata,plevel,outdata)
  
  ! write new dataset to output file
  ierr=NF_PUT_VAR_REAL(outfid,var_id(i),outdata)
  if (ierr.ne.NF_NOERR) then
    write(*,*)'Error, Failed to write variable ',var(i)
    stop
  endif
  
  enddo ! of do i=1,nbvar
  
  ! interpolate zareoid onto new grid
  call p_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
                      miss_val,ps,press,za_gcm,plevel,outdata)
  ! write result to output file
  ierr=NF_PUT_VAR_REAL(outfid,za_varid,outdata)
  if (ierr.ne.NF_NOERR) then
    stop "Error, Failed to write zareoid to output file"
  endif
endif ! of if (ztype.eq.1)

! 4.2 If output is in above areoid altitude
if (ztype.eq.2) then
  do i=1,nbvar ! loop on 4D variable to process
  ! identify and read a dataset
  ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
  if (ierr.ne.NF_NOERR) then
    write(*,*) 'Error, failed to get ID for variable ',var(i)
    stop
  endif
  ierr=NF_GET_VAR_REAL(infid,tmpvarid,indata)
  if (ierr.ne.NF_NOERR) then
    write(*,*) 'Error, failed to load variable ',var(i)
    stop
  endif
  
  ! interpolate dataset onto new grid
  ! check if variable is "rho" (to set flag for interpolation below)
  if (var(i).eq.'rho') then
    j=1
  else
    j=0
  endif
  
  call za_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
                      miss_val,za_gcm,indata,j,zareoid,outdata)

  ! write new dataset to output file
  ierr=NF_PUT_VAR_REAL(outfid,var_id(i),outdata)
  if (ierr.ne.NF_NOERR) then
    write(*,*)'Error, Failed to write variable ',var(i)
    stop
  endif

  enddo ! of do i=1,nbvar
  
  ! interpolate pressure onto new grid
  call za_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
                      miss_val,za_gcm,press,1,zareoid,outdata)

  ! write result to output file
  ierr=NF_PUT_VAR_REAL(outfid,p_varid,outdata)
  if (ierr.ne.NF_NOERR) then
    stop "Error, Failed to write pressure to output file"
  endif
endif ! of if (ztype.eq.2)

! 4.3 Close output file
ierr=NF_CLOSE(outfid)
if (ierr.ne.NF_NOERR) then
  write(*,*) 'Error, failed to close output file ',outfile
endif

end program

!===============================================================================

subroutine build_gcm_alt(lonlength,latlength,altlength,timelength, &
                         phis,ps,press,temp,rho,za_gcm)
!==============================================================================
! Purpose: Integrate hydrostatic equation in order to build the "above areoid
!          altitudes" corresponding to GCM atmospheric levels
!==============================================================================
implicit none
!==============================================================================
! Arguments:
!==============================================================================
integer,intent(in) :: lonlength ! # of points along longitude
integer,intent(in) :: latlength ! # of points along latitude
integer,intent(in) :: altlength ! # of points along altitude
integer,intent(in) :: timelength ! # of stored time steps
real,dimension(lonlength,latlength),intent(in) :: phis
! phis(:,:) is the ground geopotential
real,dimension(lonlength,latlength,timelength),intent(in) :: ps
! ps(:,:) is the surface pressure
real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: press
! press(:,:,:,:) is atmospheric pressure
real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: temp
! temp(:,:,:,:) is atmospheric temperature
real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: rho
! rho(:,:,:,:) is atmospheric density

real,dimension(lonlength,latlength,altlength,timelength),intent(out) :: za_gcm
! za_gcm(:,:,:,:) are the above aroid heights of GCM levels

!===============================================================================
! Local variables:
!===============================================================================
real,dimension(:),allocatable :: sigma ! GCM sigma levels
real,dimension(:,:,:,:),allocatable :: R ! molecular gas constant
real,dimension(:,:,:,:),allocatable :: zlocal ! local above surface altitude
real :: Rmean ! mean value of R for a given level
real :: Tmean ! "mean" value of temperature for a given level
integer iloop,jloop,kloop,tloop

! Parameters needed to integrate hydrostatic equation:
real :: g0 ! exact mean gravity at radius=3396.km (Lemoine et al. 2001)
data    g0 /3.7257964/
real :: a0  
data    a0 /3396.E3/  ! radius=3396.km
real :: gz 
! gz: gravitational acceleration at a given (zareoid) altitude

!===============================================================================
! 1. Various initialisations
!===============================================================================
allocate(sigma(altlength))
allocate(R(lonlength,latlength,altlength,timelength))
allocate(zlocal(lonlength,latlength,altlength,timelength))

!==============================================================================
! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point 
!   --later used to integrate the hydrostatic equation--
!==============================================================================
do tloop=1,timelength
  do kloop=1,altlength
    do jloop=1,latlength
      do iloop=1,lonlength
        R(iloop,jloop,kloop,tloop)=press(iloop,jloop,kloop,tloop)/ &
                                            (rho(iloop,jloop,kloop,tloop)* &
                                             temp(iloop,jloop,kloop,tloop))
      enddo
    enddo
  enddo
enddo

!===============================================================================
! 3. Integrate hydrostatic equation to compute zlocal and za_gcm
!===============================================================================
do tloop=1,timelength
  do jloop=1,latlength
    do iloop=1,lonlength
      ! handle case of first altitude level
      sigma(1)=press(iloop,jloop,1,tloop)/ps(iloop,jloop,tloop)
      zlocal(iloop,jloop,1,tloop)=-log(sigma(1))*R(iloop,jloop,1,tloop)* &
                                                 temp(iloop,jloop,1,tloop)/g0
      za_gcm(iloop,jloop,1,tloop)=zlocal(iloop,jloop,1,tloop)+ &
                                  phis(iloop,jloop)/g0
      do kloop=2,altlength
        ! compute sigma level of layer
        sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
        
        ! compute "mean" temperature of the layer
        if(temp(iloop,jloop,kloop,tloop).eq.  &
           temp(iloop,jloop,kloop-1,tloop)) then
          Tmean=temp(iloop,jloop,kloop,tloop)
        else
          Tmean=(temp(iloop,jloop,kloop,tloop)-      &
                 temp(iloop,jloop,kloop-1,tloop))/   &
                 log(temp(iloop,jloop,kloop,tloop)/  &
                     temp(iloop,jloop,kloop-1,tloop))
        endif
        
        ! compute mean value of R of the layer
        Rmean=0.5*(R(iloop,jloop,kloop,tloop)+R(iloop,jloop,kloop-1,tloop))
        
        ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
        gz=g0*(a0**2)/(a0+za_gcm(iloop,jloop,kloop-1,tloop))**2
        
        ! compute zlocal(kloop)
        zlocal(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop-1,tloop)- &
                log(sigma(kloop)/sigma(kloop-1))*Rmean*Tmean/gz
        
        ! compute za_gcm(kloop)
        za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
                                        phis(iloop,jloop)/g0
      enddo ! kloop
    enddo ! iloop
  enddo ! jloop
enddo ! tloop

! Cleanup
deallocate(sigma)
deallocate(R)
deallocate(zlocal)

end subroutine build_gcm_alt

!===============================================================================

subroutine crude_gcm_alt(lonlength,latlength,altlength,timelength, &
                         phis,ps,press,temp,za_gcm)
!==============================================================================
! Purpose: Integrate hydrostatic equation in order to build the "above areoid
!          altitudes" corresponding to GCM atmospheric levels
!==============================================================================
implicit none
!==============================================================================
! Arguments:
!==============================================================================
integer,intent(in) :: lonlength ! # of points along longitude
integer,intent(in) :: latlength ! # of points along latitude
integer,intent(in) :: altlength ! # of points along altitude
integer,intent(in) :: timelength ! # of stored time steps
real,dimension(lonlength,latlength),intent(in) :: phis
! phis(:,:) is the ground geopotential
real,dimension(lonlength,latlength,timelength),intent(in) :: ps
! ps(:,:) is the surface pressure
real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: press
! press(:,:,:,:) is atmospheric pressure
real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: temp
! temp(:,:,:,:) is atmospheric temperature

real,dimension(lonlength,latlength,altlength,timelength),intent(out) :: za_gcm
! za_gcm(:,:,:,:) are the above aroid heights of GCM levels

!===============================================================================
! Local variables:
!===============================================================================
real,dimension(:),allocatable :: sigma ! GCM sigma levels
real,parameter :: R=190 ! molecular gas constant
real,dimension(:,:,:,:),allocatable :: zlocal ! local above surface altitude
real :: Tmean ! "mean" value of temperature for a given level
integer iloop,jloop,kloop,tloop

! Parameters needed to integrate hydrostatic equation:
real :: g0 ! exact mean gravity at radius=3396.km (Lemoine et al. 2001)
data    g0 /3.7257964/
real :: a0  
data    a0 /3396.E3/  ! radius=3396.km
real :: gz 
! gz: gravitational acceleration at a given (zareoid) altitude

!===============================================================================
! 1. Various initialisations
!===============================================================================
allocate(sigma(altlength))
allocate(zlocal(lonlength,latlength,altlength,timelength))

!===============================================================================
! 2. Integrate hydrostatic equation to compute zlocal and za_gcm
!===============================================================================
do tloop=1,timelength
  do jloop=1,latlength
    do iloop=1,lonlength
      ! handle case of first altitude level
      sigma(1)=press(iloop,jloop,1,tloop)/ps(iloop,jloop,tloop)
      zlocal(iloop,jloop,1,tloop)=-log(sigma(1))* &
                                               R*temp(iloop,jloop,1,tloop)/g0
      za_gcm(iloop,jloop,1,tloop)=zlocal(iloop,jloop,1,tloop)+ &
                                  phis(iloop,jloop)/g0
      do kloop=2,altlength
        ! compute sigma level of layer
        sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
        
        ! compute "mean" temperature of the layer
        if(temp(iloop,jloop,kloop,tloop).eq.  &
           temp(iloop,jloop,kloop-1,tloop)) then
          Tmean=temp(iloop,jloop,kloop,tloop)
        else
          Tmean=(temp(iloop,jloop,kloop,tloop)-      &
                 temp(iloop,jloop,kloop-1,tloop))/   &
                 log(temp(iloop,jloop,kloop,tloop)/  &
                     temp(iloop,jloop,kloop-1,tloop))
        endif
                
        ! compute gravitational acceleration (at altitude zaeroid(kloop-1))
        gz=g0*(a0**2)/(a0+za_gcm(iloop,jloop,kloop-1,tloop))**2
        
        ! compute zlocal(kloop)
        zlocal(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop-1,tloop)- &
                log(sigma(kloop)/sigma(kloop-1))*R*Tmean/gz
        
        ! compute za_gcm(kloop)
        za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
                                        phis(iloop,jloop)/g0
      enddo ! kloop
    enddo ! iloop
  enddo ! jloop
enddo ! tloop

! Cleanup
deallocate(sigma)
deallocate(zlocal)

end subroutine crude_gcm_alt

!===============================================================================

subroutine p_coord_interp(lonlen,latlen,altlen,tlen,newaltlen, &
                          missing,ps,press,gcmdata,plevels,newdata)
!==============================================================================
! Purpose:
! Recast a 4D (spatio-temporal) GCM dataset, which has a vertical coordinate
! in pseudo-altitude, into a dataset which has a vertical coordinate at given
! pressure levels
!==============================================================================
implicit none
!==============================================================================
! Arguments:
!==============================================================================
integer,intent(in) :: lonlen ! # of points along longitude (GCM dataset)
integer,intent(in) :: latlen ! # of points along latitude (GCM dataset)
integer,intent(in) :: altlen ! # of points along altitude (GCM dataset)
integer,intent(in) :: tlen   ! # of stored time steps (GCM dataset)
integer,intent(in) :: newaltlen ! # of points along altitude
real,intent(in) :: missing ! missing value
real,dimension(lonlen,latlen,tlen),intent(in) :: ps ! GCM surface pressure
real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: press ! GCM pressure
real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: gcmdata ! GCM dataset
real,dimension(newaltlen),intent(in) :: plevels
! plevels(:) pressure levels at which gcmdata is to be interpolated

real,dimension(lonlen,latlen,newaltlen,tlen),intent(out) :: newdata
! newdata(:,:,:,:) gcmdata recasted along vertical coordinate
!==============================================================================
! Local variables:
!==============================================================================
real,dimension(:),allocatable :: lnp
! lnp(:): used to store log(pressure) values
real,dimension(:),allocatable :: q
! q(:): used to store values along vertical direction
real :: x,y
! x,y: temporary variables
integer :: iloop,jloop,kloop,tloop

allocate(lnp(altlen))
allocate(q(altlen))

do tloop=1,tlen
  do jloop=1,latlen
    do iloop=1,lonlen
      ! build lnp() and q() along vertical coordinate
      do kloop=1,altlen
        lnp(kloop)=-log(press(iloop,jloop,kloop,tloop))
        q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
      enddo
      
      ! Interpolation  
      do kloop=1,newaltlen
        ! compute, by interpolation, value at pressure level plevels(kloop)
        if ((plevels(kloop).le.ps(iloop,jloop,tloop)).and. &
            (plevels(kloop).ge.press(iloop,jloop,altlen,tloop))) then
          x=-log(plevels(kloop))
          call interpolf(x,y,missing,lnp,q,altlen)
          newdata(iloop,jloop,kloop,tloop) = y
        else ! if plevels(kloop) is out of range,
          ! assign a "missing_value" at this node
          newdata(iloop,jloop,kloop,tloop) = missing
        endif
      enddo
      
    enddo !iloop
  enddo !jloop
enddo !tloop

! Cleanup
deallocate(lnp)
deallocate(q)

end subroutine p_coord_interp

!===============================================================================

subroutine za_coord_interp(lonlen,latlen,altlen,tlen,newaltlen, &
                          missing,za_gcm,gcmdata,flag,zareoid,newdata)
!==============================================================================
! Purpose:
! Recast a 4D (spatio-temporal) GCM dataset, for which corresponding grid
! values of above areoid altitude are known, along chosen above areoid
! altitude
!==============================================================================
implicit none
!==============================================================================
! Arguments:
!==============================================================================
integer,intent(in) :: lonlen ! # of points along longitude (GCM dataset)
integer,intent(in) :: latlen ! # of points along latitude (GCM dataset)
integer,intent(in) :: altlen ! # of points along altitude (GCM dataset)
integer,intent(in) :: tlen   ! # of stored time steps (GCM dataset)
integer,intent(in) :: newaltlen ! # of points along altitude
real,intent(in) :: missing ! missing value
real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: za_gcm
!za_gcm(:,:,:,) above areoid altitude of GCM grid points
real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: gcmdata ! GCM dataset
integer,intent(in) :: flag ! flag (==1 for 'log' interpolation)
! flag==0 (standard linear interpolation)
real,dimension(newaltlen),intent(in) :: zareoid
! zareoid(:) above areoid altitudes at which data must be recast

real,dimension(lonlen,latlen,newaltlen,tlen),intent(out) :: newdata
! newdata(:,:,:,:) gcmdata recasted along vertical coordinate

!==============================================================================
! Local variables:
!==============================================================================
real,dimension(:),allocatable :: za,q
real,dimension(:),allocatable :: logq
! za(:): used to store za_gcm at a given location
! q(:): used to store field values along the vertical direction
real :: x,y ! temporary variables
integer :: iloop,jloop,kloop,tloop

allocate(za(altlen))
allocate(q(altlen))
allocate(logq(altlen))

if (flag.eq.0) then
 do tloop=1,tlen
  do jloop=1,latlen
    do iloop=1,lonlen
      do kloop=1,altlen
        ! extract the vertical coordinates zaij(), or build p_za()
        za(kloop)=za_gcm(iloop,jloop,kloop,tloop)
        ! store values along altitude
        q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
      enddo !kloop
      ! Interpolation
      do kloop=1,newaltlen
      ! Check if zareoid(kloop) is within range
      if ((zareoid(kloop).ge.za_gcm(iloop,jloop,1,tloop)).and. &
          (zareoid(kloop).le.za_gcm(iloop,jloop,altlen,tloop))) then
        ! zareoid(kloop) is within range
        x=zareoid(kloop)
        call interpolf(x,y,missing,za,q,altlen)
        newdata(iloop,jloop,kloop,tloop)=y
      else ! zareoid(kloop) is out of range
        newdata(iloop,jloop,kloop,tloop)=missing
      endif 
      enddo !kloop
    enddo !iloop
  enddo !jloop
 enddo !tloop 
else ! case when flag=1 (i.e.: rho)
 do tloop=1,tlen
  do jloop=1,latlen
    do iloop=1,lonlen
      do kloop=1,altlen
        ! extract the vertical coordinates zaij(), or build p_za()
        za(kloop)=za_gcm(iloop,jloop,kloop,tloop)
        ! extract the vertical coordinates zaij(), or build p_za()
        za(kloop)=za_gcm(iloop,jloop,kloop,tloop)
        ! store log values along altitude
        logq(kloop)=log(gcmdata(iloop,jloop,kloop,tloop))
      enddo !kloop

      ! Interpolation
      do kloop=1,newaltlen
      ! Check if zareoid(kloop) is within range
      if ((zareoid(kloop).ge.za_gcm(iloop,jloop,1,tloop)).and. &
          (zareoid(kloop).le.za_gcm(iloop,jloop,altlen,tloop))) then
        ! zareoid(kloop) is within range
        x=zareoid(kloop)
        call interpolf(x,y,missing,za,logq,altlen)
        newdata(iloop,jloop,kloop,tloop)=exp(y)
      else  ! zareoid(kloop) is out of range
        newdata(iloop,jloop,kloop,tloop)=missing
      endif
      enddo !kloop
    enddo !iloop
  enddo !jloop
 enddo !tloop 
endif

! Cleanup
deallocate(za)
deallocate(q)
deallocate(logq)

end subroutine za_coord_interp

!===============================================================================

subroutine interpolf(x,y,missing,xd,yd,nd)
!==============================================================================
! Purpose:
! Yield y=f(x), where xd() end yd() are arrays of known values,
! using linear interpolation
! If x is not included in the interval spaned by xd(), then y is set
! to a default value 'missing'
! Note:
! Array xd() should contain ordered (either increasing or decreasing) abscissas
!==============================================================================
implicit none
!==============================================================================
! Arguments:
!==============================================================================
real,intent(in) :: x ! abscissa at which interpolation is to be done
real,intent(in) :: missing ! missing value (if no interpolation is performed)
integer :: nd ! size of arrays
real,dimension(nd),intent(in) :: xd ! array of known absissas
real,dimension(nd),intent(in) :: yd ! array of correponding values

real,intent(out) :: y ! interpolated value
!==============================================================================
! Local variables:
!==============================================================================
integer :: i

! default: set y to 'missing'
y=missing

   do i=1,nd-1
     if (((x.ge.xd(i)).and.(x.le.xd(i+1))).or.&
          ((x.le.xd(i)).and.(x.ge.xd(i+1)))) then
        y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i))
        exit
     endif
   enddo


end subroutine interpolf