source: trunk/LMDZ.VENUS/Tools/zrecast.F90 @ 937

program zrecast

! SL 12/2009:
! Adaptation of the martian zrecast for use with Titan and Venus outputs.

! 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 gcm_layers_dimid ! NetCDF dimension ID for # of layers in GCM
integer sigma_varid,aps_varid,bps_varid,ls_varid
integer za_varid,p_varid ! above areoid and pressure data IDs



integer ps_varid,phis_varid ! surface pressure and surface geopotential data ID
integer,dimension(4) :: datashape ! shape of 4D datasets
integer,dimension(3) :: surfdatashape ! shape of 3D (surface+time) datasets
integer,dimension(2) :: surfdatashape2d ! shape of 2D (surface) datasets

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 :: ls ! solar longitude
real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
real,dimension(:,:,:),allocatable :: ps ! GCM surface pressure
real,dimension(:,:,:),allocatable :: ls2d ! solar longitude (3D field in input - Titan...)
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 (m)
real,dimension(:,:,:,:),allocatable :: zs_gcm ! GCM above surface heights (m)
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  ztype=3: above local surface
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
real,dimension(:),allocatable :: zsurface ! Above surface heights for output
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)
logical :: have_hyb   ! Flag: true if aps and bps hybrid levels are known 
                      !  (false if sigma levels or press (3D field) are used)
logical :: auto_vert_levels ! Flag: true if the positions of vertical levels
                            ! has to be computed; false if these are given
                            ! by the user

include "planet.h"

!===============================================================================
! 1. Input parameters
!===============================================================================

write(*,*) ""
write(*,*) "You are working on the atmosphere of ",planet

!===============================================================================
! 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
    ! also check that it is indeed a 4D variable
    ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims)
    if (tmpndims.eq.4) then
      nbvar=nbvar+1
      var(nbvar)=tmpvarname    
    else
      write(*,*) 'Error: ',trim(tmpvarname),' is not a 4D variable'
      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 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
  write(*,*) "Could not get latitude dimension ID"
  write(*,*) "  looking for lat dimension instead... "
  ierr=NF_INQ_DIMID(infid,"lat",tmpdimid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed to get lat dimension ID"
  else
    ierr=NF_INQ_VARID(infid,"lat",tmpvarid)
    if (ierr.ne.NF_NOERR) then
      stop "Error: Failed to get lat ID"
    else
      ierr=NF_INQ_DIMLEN(infid,tmpdimid,latlength)
      if (ierr.ne.NF_NOERR) then
        stop "Error: Failed to get lat length"
      else
        allocate(lat(latlength))
        ierr=NF_GET_VAR_REAL(infid,tmpvarid,lat)
        if (ierr.ne.NF_NOERR) then
          stop "Error: Failed reading lat"
        endif
      endif
    endif
  endif
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
  write(*,*) "Could not get longitude dimension ID"
  write(*,*) "  looking for lon dimension instead... "
  ierr=NF_INQ_DIMID(infid,"lon",tmpdimid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed to get lon dimension ID"
  else
    ierr=NF_INQ_VARID(infid,"lon",tmpvarid)
    if (ierr.ne.NF_NOERR) then
      stop "Error: Failed to get lon ID"
    else
      ierr=NF_INQ_DIMLEN(infid,tmpdimid,lonlength)
      if (ierr.ne.NF_NOERR) then
        stop "Error: Failed to get lon 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
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
  write(*,*) "Could not get Time dimension ID"
  write(*,*) "  looking for time dimension instead... "
  ierr=NF_INQ_DIMID(infid,"time",tmpdimid)
 if (ierr.ne.NF_NOERR) then
  write(*,*) "Could not get time dimension ID"
  write(*,*) "  looking for time_counter dimension instead... "
  ierr=NF_INQ_DIMID(infid,"time_counter",tmpdimid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed to get time_counter dimension ID"
  else
    ierr=NF_INQ_VARID(infid,"time_counter",tmpvarid)
    if (ierr.ne.NF_NOERR) then
      stop "Error: Failed to get time_counter ID"
    else
      ierr=NF_INQ_DIMLEN(infid,tmpdimid,timelength)
      if (ierr.ne.NF_NOERR) then
        stop "Error: Failed to get time_counter length"
      else
        allocate(time(timelength))
        ierr=NF_GET_VAR_REAL(infid,tmpvarid,time)
        if (ierr.ne.NF_NOERR) then
          stop "Error: Failed reading time_counter"
        endif
      endif
    endif
  endif
 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
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
  write(*,*) "Could not get altitude dimension ID"
  write(*,*) "  looking for sigma dimension instead... "
  ierr=NF_INQ_DIMID(infid,"sigma",tmpdimid)
 if (ierr.ne.NF_NOERR) then
  write(*,*) "Could not get sigma dimension ID"
  write(*,*) "  looking for presnivs dimension instead... "
  ierr=NF_INQ_DIMID(infid,"presnivs",tmpdimid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed to get presnivs dimension ID"
  else
    ierr=NF_INQ_DIMLEN(infid,tmpdimid,altlength)
    if (ierr.ne.NF_NOERR) then
      stop "Error: Failed to get presnivs length"
    endif
  endif
 else
    ierr=NF_INQ_DIMLEN(infid,tmpdimid,altlength)
    if (ierr.ne.NF_NOERR) then
      stop "Error: Failed to get sigma length"
    endif
 endif
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
    have_hyb=.false.
    write(*,*) "Could not find aps coordinates... will look for 3D pres field"
  else
    have_hyb=.true.
    allocate(aps(altlength))
    ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
    if (ierr.ne.NF_NOERR) then
      stop "Error: Failed reading aps"
    endif
  endif
 if (have_hyb) then
  ! 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
endif !of if (.not.have_sigma)

if ((.not.have_sigma).and.(.not.have_hyb)) then
  ! direct press field
  ierr=NF_INQ_VARID(infid,"pres",tmpvarid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed to get pres ID"
  else
    allocate(press(lonlength,latlength,altlength,timelength))
    ierr=NF_GET_VAR_REAL(infid,tmpvarid,press)
    if (ierr.ne.NF_NOERR) then
      stop "Error: Failed reading pres"
    endif
  endif
endif !of if ((.not.have_sigma).and.(.not.have_hyb))

! 1.3.3 ground geopotential phisinit

allocate(phisinit(lonlength,latlength),stat=ierr)
if (ierr.ne.0) then
  write(*,*) "Failed allocation of phisinit(lonlength,latlength) !!!"
  write(*,*) "lonlength=",lonlength," latlength=",latlength
endif
! look for 'phisinit' in current file
ierr=NF_INQ_VARID(infid,"phis",tmpvarid) 
if (ierr.ne.NF_NOERR) then
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
else
  write(*,*) 'OK, got phis for phisinit'
  ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed reading phis"
  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.lt.1).or.(ztype.gt.3))
  write(*,*) "(1: pressure, 2: above areoid altitude 3: above local surface)"
  read(*,*)ztype
enddo

text="dummy" ! dummy initialization
auto_vert_levels=.true. ! dummy initialization to get rid of compiler warning
do while ((trim(text).ne."yes").and.(trim(text).ne."no"))
  write(*,*) ""
  write(*,*) "Automatic generation of vertical levels distribution? (yes/no)"
  write(*,*) "(yes: you only provide min, max and number of levels)"
  write(*,*) "(no: you provide values for each level)"
  read(*,'(a64)') text
  if (trim(text).eq."yes") then
    auto_vert_levels=.true.
  else
    auto_vert_levels=.false.
  endif
enddo

if (auto_vert_levels) then
  ! ask for # of points and end values for pressure or above areoid cases
  write(*,*) ""
  if (ztype.le.2) then
    write(*,*) "Enter min and max of vertical coordinate (Pa or m)"
    write(*,*) " (in that order and on the same line)"
    if (ztype.eq.1) then ! pressure coordinate
      read(*,*) pmin,pmax
    else ! above areoid coordinate
      read(*,*) zamin,zamax
    endif
  endif

  ! Build corresponding vertical coordinates
  if (ztype.eq.1) then ! pressure coordinate
    write(*,*) "Number of levels along vertical coordinate?"
    read(*,*) nblev
    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 if (ztype.eq.2) then ! above areoid heights
    write(*,*) "Number of levels along vertical coordinate?"
    read(*,*) nblev
    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
  else ! above local surface
    ! set nblev to # of levels from input data files
    nblev=altlength
    allocate(zsurface(nblev))
    ! build specific above local surface altitudes
    if (have_sigma.or.have_hyb) then
      call build_zs(nblev,have_sigma,sigma,aps,bps,zsurface)
    else
      stop "Case not dealt with yet..."
    endif
  endif
else ! auto_vert_levels=.false. ; user provides values
  ! ask for # of points along the vertical
  write(*,*) ""
  write(*,*) "Number of levels along vertical coordinate?"
  read(*,*) nblev
  if (ztype.eq.1) then ! pressure coordinate
    allocate(plevel(nblev))
    write(*,*) "Enter Pressure (Pa) of levels, ordered"
    write(*,*) " from max (near-surface) to min (top of atmosphere),"
    write(*,*) " (one value per line)"
    do i=1,nblev
      read(*,*) plevel(i)
    enddo
  else if (ztype.eq.2) then ! above areoid heights
    allocate(zareoid(nblev))
    write(*,*) "Enter altitude (m) above areoid of levels, ordered"
    write(*,*) " from min to max (one value per line)"
    do i=1,nblev
      read(*,*) zareoid(i)
    enddo
  else ! above local surface
    allocate(zsurface(nblev))
    write(*,*) "Enter altitude (m) above surface of levels, ordered"
    write(*,*) " from min to max (one value per line)"
    do i=1,nblev
      read(*,*) zsurface(i)
    enddo
  endif
endif ! of if (auto_vert_levels)

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

!===============================================================================
! 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
  write(*,*) "Could not get ps ID"
  write(*,*) "  looking for psol instead... "
 ierr=NF_INQ_VARID(infid,"psol",tmpvarid)
 if (ierr.ne.NF_NOERR) then
  stop "Error: Failed to get psol 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
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.1 If Titan, solar longitude ls2d
!===============================================================================
if (planet.eq."Titan") then

ierr=NF_INQ_VARID(infid,"ls",tmpvarid)
if (ierr.ne.NF_NOERR) then
  write(*,*) "Could not get ls ID"
    stop "Error: Failed reading solar longitude"
else
  allocate(ls2d(lonlength,latlength,timelength))
  ierr=NF_GET_VAR_REAL(infid,tmpvarid,ls2d)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Failed reading solar longitude"
  endif
  allocate(ls(timelength))
  ls = ls2d(1,1,:)
endif

endif

!===============================================================================
! 2.1.2 Atmospheric pressure
!===============================================================================
if (have_sigma.or.have_hyb) then
! we compute the 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

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 (or above surface) altitudes of GCM nodes
!===============================================================================

if (have_rho) then
  if (ztype.le.2) then ! above areoid altitudes (also needed for ztype=1)
    allocate(za_gcm(lonlength,latlength,altlength,timelength))
    call build_gcm_za(lonlength,latlength,altlength,timelength, &
                      phisinit,ps,press,temp,rho,za_gcm)
  else ! above local surface altitudes
    allocate(zs_gcm(lonlength,latlength,altlength,timelength))
    call build_gcm_zs(lonlength,latlength,altlength,timelength, &
                      phisinit,ps,press,temp,rho,zs_gcm)
  endif
else
  write(*,*)"Warning: Using constant R to integrate hydrostatic equation"
  if (ztype.le.2) then ! above areoid altitudes (also needed for ztype=1)
    allocate(za_gcm(lonlength,latlength,altlength,timelength))
    call crude_gcm_za(lonlength,latlength,altlength,timelength, &
                      phisinit,ps,press,temp,za_gcm)
  else ! above local surface altitudes
    allocate(zs_gcm(lonlength,latlength,altlength,timelength))
    call crude_gcm_zs(lonlength,latlength,altlength,timelength, &
                      phisinit,ps,press,temp,zs_gcm)
  endif
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
if (ztype.eq.1) then
  ierr=NF_DEF_DIM(outfid,"presnivs",nblev,alt_dimid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Could not define presnivs dimension"
  endif
else
  ierr=NF_DEF_DIM(outfid,"altitude",nblev,alt_dimid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Could not define altitude dimension"
  endif
endif

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

! GCM layers (for sigma or aps and bps)
ierr=NF_DEF_DIM(outfid,"GCM_layers",altlength,gcm_layers_dimid)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not define GCM_layers dimension"
endif


!===============================================================================
! 3.3. Define variables and their attributes
!===============================================================================

! 3.3.1 Define 1D variables

! 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
if (ztype.eq.1) then
  ierr=NF_DEF_VAR(outfid,"presnivs",NF_REAL,1,datashape(3),alt_varid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Could not define presnivs variable"
  endif
else
  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
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 presnivs"
  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 presnivs"
  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 presnivs"
  endif
else if (ztype.eq.2) then ! 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
else ! above surface vertical coordinate
  text='Altitude above local surface'
  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 

! GCM_layers
!ierr=NF_DEF_VAR(outfid,"gcm_layers",NF_REAL,1,,sigma_varid)

! sigma levels or hybrid coordinates
if (have_sigma) then
  ierr=NF_DEF_VAR(outfid,"sigma",NF_REAL,1,gcm_layers_dimid,sigma_varid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Could not define sigma variable"
  endif
elseif (have_hyb) then ! hybrid coordinates
  ierr=NF_DEF_VAR(outfid,"aps",NF_REAL,1,gcm_layers_dimid,aps_varid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Could not define aps variable"
  endif
  ierr=NF_DEF_VAR(outfid,"bps",NF_REAL,1,gcm_layers_dimid,bps_varid)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Could not define bps variable"
  endif
endif

! sigma levels (or hybrid coordinates) attributes
if (have_sigma) then
  text="sigma levels"
  ierr=NF_PUT_ATT_TEXT(outfid,sigma_varid,'long_name',len_trim(text),text)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing long_name for sigma"
  endif
elseif (have_hyb) then  ! hybrid coordinates
  text="hybrid pressure at midlayers"
  ierr=NF_PUT_ATT_TEXT(outfid,aps_varid,'long_name',len_trim(text),text)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing long_name for aps"
  endif
  text="hybrid sigma at midlayers"
  ierr=NF_PUT_ATT_TEXT(outfid,bps_varid,'long_name',len_trim(text),text)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Problem writing long_name for bps"
  endif
endif ! of if (have_sigma)

! 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

! Solar longitude for Titan
if (planet.eq."Titan") then

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

! Ls attributes
text='Ls'
ierr=NF_PUT_ATT_TEXT(outfid,ls_varid,'long_name',len_trim(text),text)
if (ierr.ne.NF_NOERR) then
  stop "Error: Problem writing long_name for Ls"
endif
text='Degrees'
ierr=NF_PUT_ATT_TEXT(outfid,ls_varid,'units',len_trim(text),text)
if (ierr.ne.NF_NOERR) then
  stop "Error: Problem writing units for Ls"
endif

endif ! ls for Titan

! 3.3.2 Define 2D (surface only) and 3D variables (ie: surface+time variables)

! Surface geopotential
surfdatashape2d(1)=lon_dimid
surfdatashape2d(2)=lat_dimid
ierr=NF_DEF_VAR(outfid,"phis",NF_REAL,2,surfdatashape2d,phis_varid)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not define phis variable"
endif

! Surface geopotential attributes
text='Surface geopotential'
ierr=NF_PUT_ATT_TEXT(outfid,phis_varid,'long_name',len_trim(text),text)
if (ierr.ne.NF_NOERR) then
  stop "Error: Problem writing long_name for surface geopotential"
endif
text='m2 s-2'
ierr=NF_PUT_ATT_TEXT(outfid,phis_varid,'units',len_trim(text),text)
if (ierr.ne.NF_NOERR) then
  stop "Error: Problem writing units for surface geopotential"
endif

! Surface pressure
surfdatashape(1)=lon_dimid
surfdatashape(2)=lat_dimid
surfdatashape(3)=time_dimid
ierr=NF_DEF_VAR(outfid,"ps",NF_REAL,3,surfdatashape,ps_varid)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not define ps variable"
endif

! Surface pressure attributes
text='Surface pressure'
ierr=NF_PUT_ATT_TEXT(outfid,ps_varid,'long_name',len_trim(text),text)
if (ierr.ne.NF_NOERR) then
  stop "Error: Problem writing long_name for surface pressure"
endif
text='Pa'
ierr=NF_PUT_ATT_TEXT(outfid,ps_varid,'units',len_trim(text),text)
if (ierr.ne.NF_NOERR) then
  stop "Error: Problem writing units for surface pressure"
endif

! 3.3.3 Define 4D variables

! 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 or above local surface 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

! add zs_gcm
if (ztype.eq.3) then
endif

! variables requested by user
allocate(var_id(nbvar))
do i=1,nbvar
  write(*,*) ""
  write(*,*) "Creating ",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
    else ! no 'title' attribute
      ! try to find a "Physics_diagnostic" attribute (UK GCM outputs)
      text='   '
      ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'Physics_diagnostic',text)
      if (ierr.eq.NF_NOERR) then ! found 'Physics_diagnostic' attribute
        write(*,*) "Found Physics_diagnostic ",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 Physics_diagnostic attribute:",trim(text)
          stop
        endif
      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
    if (miss_val.ne.miss_val_def) then
       miss_val=miss_val_def
    write(*,*) "set to default:",miss_val_def
    endif
    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 (and 1D variables)
!===============================================================================
! 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 pressure data to output file"
  endif
else if (ztype.eq.2) then ! 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
else ! above local surface
  ierr=NF_PUT_VAR_REAL(outfid,alt_varid,zsurface)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Could not write altitude data to output file"
  endif
endif

! Write sigma levels (or hybrid coordinates)
if (have_sigma) then
  ierr=NF_PUT_VAR_REAL(outfid,sigma_varid,sigma)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Could not write sigma data to output file"
  endif
elseif (have_hyb) then ! hybrid coordinates
  ierr=NF_PUT_VAR_REAL(outfid,aps_varid,aps)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Could not write aps data to output file"
  endif
  ierr=NF_PUT_VAR_REAL(outfid,bps_varid,bps)
  if (ierr.ne.NF_NOERR) then
    stop "Error: Could not write bps 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

! Solar longitude for Titan
if (planet.eq."Titan") then

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

endif ! ls for Titan

!===============================================================================
! 3.5 Write 2D and 3D variables
!===============================================================================

! Write surface geopotential

ierr=NF_PUT_VAR_REAL(outfid,phis_varid,phisinit)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not write phis data to output file"
endif

! Write surface pressure

ierr=NF_PUT_VAR_REAL(outfid,ps_varid,ps)
if (ierr.ne.NF_NOERR) then
  stop "Error: Could not write ps data to output file"
endif

!===============================================================================
! 4. Interpolate and write 4D 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 z_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 z_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 If output is in above local surface altitude
if (ztype.eq.3) then
  do i=1,nbvar ! loop on 4D variable to process
    write(*,*) " "
    write(*,*) "Processing "//trim(var(i))
    ! 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 zs_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
                        miss_val,zs_gcm,indata,phisinit,press,temp,rho,&
                        j,zsurface,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
  write(*,*) ""
  write(*,*) "Processing pressure"
  call zs_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
                      miss_val,zs_gcm,press,phisinit,press,temp,rho,&
                      1,zsurface,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.3)

! 4.4 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_zs(lonlength,latlength,altlength,timelength, &
                         phis,ps,press,temp,rho,zs_gcm)
!==============================================================================
! Purpose: Integrate hydrostatic equation in order to build the "above local
!          surface 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) :: zs_gcm
! zs_gcm(:,:,:,:) are the above local surface altitudes 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

include "planet.h"

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)
      zs_gcm(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))
        ! NB: zareoid=zsurface+phis/g0
        gz=g0*(a0**2)/ &
           (a0+zs_gcm(iloop,jloop,kloop-1,tloop)+phis(iloop,jloop)/g0)**2
        
        ! compute zs_gcm(iloop,jloop,kloop,tloop)
        zs_gcm(iloop,jloop,kloop,tloop)=zs_gcm(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_zs

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

!#include"build_gcm_alt.F90"
subroutine build_gcm_za(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

include "planet.h"

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_za

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

subroutine build_zs(altlength,have_sigma,sigma,aps,bps,zsurface)
!==============================================================================
! Build generic above surface altitudes, using either sigma levels
! or hybrid vertical coordinates.
! In order to do so, we need to use scale heights that vary with altitude.
! The scale heights distribution is then set to vary from a min. to a max.
! following a tanh() law over an imposed transition range (see below).
! Similarily, a reference mean surface pressure (610 Pa) is used to
! compute generic above surface altitudes.
!
! MARS ONLY
!==============================================================================
implicit none
!==============================================================================
! Arguments:
!==============================================================================
integer,intent(in) :: altlength ! # of points along altitude
logical,intent(in) :: have_sigma ! true if sigma(:) are known
! have_sigma is false if aps() and bps(:) are given instead
real,dimension(altlength),intent(in) :: sigma ! sigma levels
real,dimension(altlength),intent(in) :: aps ! hybrid pressure levels
real,dimension(altlength),intent(in) :: bps ! hybrid sigma levels
real,dimension(altlength),intent(out) :: zsurface ! altitudes (m)

!===============================================================================
! Local variables:
!===============================================================================
real,dimension(:),allocatable :: H ! scale heights
real,parameter :: H_low=9650   ! scale height at low altitudes
real,parameter :: H_high=15000 ! scale height at high altitudes
real,parameter :: trans_window=10 ! # of levels over which H(:) goes
                                  ! from H_low to H_high
real,parameter :: lev_trans=32+trans_window/2
! lev_trans: level at which H(lev_trans)=(H_low+H_high)/2
! N.B.: keep lev_trans and lev_trans as reals to avoid truncation issues

real,parameter :: P_ref=610 ! reference surface pressure used to build zsurface
integer i 

! 1. Build scale heights
allocate(H(altlength))
do i=1,altlength
  H(i)=H_low+(H_high-H_low)*0.5*(1.0+tanh(6.*(i-lev_trans)/trans_window))
enddo

! 2. Compute zsurface(:)
if (have_sigma) then ! use sigma levels
  do i=1,altlength
    zsurface(i)=-H(i)*log(sigma(i)*P_ref)
  enddo
else ! use hybrid coordinates
  do i=1,altlength
    zsurface(i)=-H(i)*log((aps(i)/P_ref)+bps(i))
  enddo
endif

! Cleanup
deallocate(H)

end subroutine build_zs

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

subroutine crude_gcm_zs(lonlength,latlength,altlength,timelength, &
                         phis,ps,press,temp,zs_gcm)
!==============================================================================
! Purpose: Integrate hydrostatic equation in order to build the "above local
!          surface 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) :: zs_gcm
! zs_gcm(:,:,:,:) are the above local surface altitudes of GCM levels

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

include "planet.h"

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)
      zs_gcm(iloop,jloop,1,tloop)=-log(sigma(1))* &
                                              R0*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))
        ! NB: zareoid=zsurface+phis/g0
        gz=g0*(a0**2)/ &
           (a0+zs_gcm(iloop,jloop,kloop-1,tloop)+phis(iloop,jloop)/g0)**2
        
        ! compute zs_gcm(iloop,jloop,kloop,tloop)
        zs_gcm(iloop,jloop,kloop,tloop)=zs_gcm(iloop,jloop,kloop-1,tloop)- &
                log(sigma(kloop)/sigma(kloop-1))*R0*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_zs

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

!#include"crude_gcm_alt.F90"
subroutine crude_gcm_za(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,dimension(:,:,:,:),allocatable :: zlocal ! local above surface altitude
real :: Tmean ! "mean" value of temperature for a given level
integer iloop,jloop,kloop,tloop

include "planet.h"

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))* &
                                              R0*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))*R0*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_za

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

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

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

!#include"za_coord_interp.F90"
subroutine z_coord_interp(lonlen,latlen,altlen,tlen,newaltlen, &
                          missing,z_gcm,gcmdata,flag,z_new,newdata)
!==============================================================================
! Purpose:
! Recast a 4D (spatio-temporal) GCM dataset 'gcmdata', for which corresponding
! grid values of altitude are known 'z_gcm', onto a new altitude grid 'z_new'.
! "Altitudes" can be above areoid or above local surface altitudes, as long as
! both 'z_gcm' and 'znew' refer to altitudes of a same type.
! Note: If altitudes in 'znew' fall outside of the range of altitudes
!       in 'z_gcm' then corresponding 'newdata' value is set to 'missing'.
!       'z_gcm' and 'znew' altitudes must be monotone increasing sequences.
!==============================================================================
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) :: z_gcm
!z_gcm(:,:,:,) GCM grid points altitude (above areoid or above surface)
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) :: z_new
! z_new(:) altitudes (above areoid or surface) 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 z_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
        za(kloop)=z_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 z_new(kloop) is within range
      if ((z_new(kloop).ge.z_gcm(iloop,jloop,1,tloop)).and. &
          (z_new(kloop).le.z_gcm(iloop,jloop,altlen,tloop))) then
        ! z_new(kloop) is within range
        x=z_new(kloop)
        call interpolf(x,y,missing,za,q,altlen)
        newdata(iloop,jloop,kloop,tloop)=y
      else ! z_new(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 
        za(kloop)=z_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 z_new(kloop) is within range
      if ((z_new(kloop).ge.z_gcm(iloop,jloop,1,tloop)).and. &
          (z_new(kloop).le.z_gcm(iloop,jloop,altlen,tloop))) then
        ! z_new(kloop) is within range
        x=z_new(kloop)
        call interpolf(x,y,missing,za,logq,altlen)
        newdata(iloop,jloop,kloop,tloop)=exp(y)
      else  ! z_new(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 z_coord_interp

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

subroutine zs_coord_interp(lonlen,latlen,altlen,tlen,newaltlen, &
                          missing,z_gcm,gcmdata,phis,press,temp,rho, &
                          flag,z_new,newdata)
!==============================================================================
! Purpose:
! Recast a 4D (spatio-temporal) GCM dataset 'gcmdata', for which corresponding
! grid values of altitude are known 'z_gcm', onto a new altitude grid 'z_new'.
! "Altitudes" must be above local surface altitudes.
! Notes: 
!       'z_gcm' and 'znew' altitudes must be monotone increasing sequences.
!       If altitudes in 'znew(i)' fall below those in 'z_gcm(:,:,1,:)', then
!       'newdata(:,:,i,:)' is set to constant value 'gcmdata(:,:,1,:)' if
!       flag=0, and extrapolated (exponentially) if flag=1.
!       If altitudes in 'znew(i)' are above those in 'z_gcm(:,:,altlen,:)', then
!       'newdata(:,:,i,:)' is set to constant value 'gcmdata(:,:,altlen,:)' 
!       if flag=0, and extrapolated (exponentially) if flag=1.
!==============================================================================
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) :: z_gcm
!z_gcm(:,:,:,) GCM grid points altitude (above areoid or above surface)
real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: gcmdata ! GCM dataset
real,dimension(lonlen,latlen),intent(in) :: phis
! phis(:,:) is the ground geopotential
real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: press
! press(:,:,:,:) is atmospheric pressure on GCM levels
real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: temp
! temp(:,:,:,:) is atmospheric temperature on GCM levels
real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: rho
! rho(:,:,:,:) is atmospheric density on GCM levels
integer,intent(in) :: flag ! flag (==1 for 'log' interpolation)
! flag==0 (standard linear interpolation)
real,dimension(newaltlen),intent(in) :: z_new
! z_new(:) altitudes (above areoid or surface) 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 :: z,q
real,dimension(:),allocatable :: logq
! z(:): used to store z_gcm at a given location
! q(:): used to store field values along the vertical direction
real,dimension(:,:,:),allocatable :: Rbottom,Rtop
! values of R (gas constant) below and above GCM layers
real :: x,y ! temporary variables
integer :: iloop,jloop,kloop,tloop
include "planet.h"
real :: gz 
! gz: gravitational acceleration at a given (above areoid) altitude

allocate(z(altlen))
allocate(q(altlen))
allocate(logq(altlen))

! 1. Build Rbottom and Rtop (only necessary if flag=1)
if (flag.eq.1) then
  allocate(Rbottom(lonlen,latlen,tlen))
  allocate(Rtop(lonlen,latlen,tlen))
  do iloop=1,lonlen
    do jloop=1,latlen
      do tloop=1,tlen
        Rbottom(iloop,jloop,tloop)=press(iloop,jloop,1,tloop)/ &
                                            (rho(iloop,jloop,1,tloop)* &
                                             temp(iloop,jloop,1,tloop))
        Rtop(iloop,jloop,tloop)=press(iloop,jloop,altlen,tloop)/ &
                                            (rho(iloop,jloop,altlen,tloop)* &
                                             temp(iloop,jloop,altlen,tloop))
      enddo
    enddo
  enddo
endif

! 2. Interpolation
if (flag.eq.0) then
 do tloop=1,tlen
  do jloop=1,latlen
    do iloop=1,lonlen
      ! preliminary stuff
      do kloop=1,altlen
        ! extract the vertical coordinates
        z(kloop)=z_gcm(iloop,jloop,kloop,tloop)
        ! store values along altitude
        q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
      enddo !kloop
      
      ! Interpolation
      do kloop=1,newaltlen
        if (z_new(kloop).lt.z_gcm(iloop,jloop,1,tloop)) then
          ! z_new(kloop) is below lowest GCM level
          newdata(iloop,jloop,kloop,tloop)=gcmdata(iloop,jloop,1,tloop)
        else if (z_new(kloop).gt.z_gcm(iloop,jloop,altlen,tloop)) then
          ! z_new(kloop) is above highest GCM level
          newdata(iloop,jloop,kloop,tloop)=gcmdata(iloop,jloop,altlen,tloop)
        else ! z_new(kloop) is within range
          x=z_new(kloop)
          call interpolf(x,y,missing,z,q,altlen)
          newdata(iloop,jloop,kloop,tloop)=y
        endif
      enddo !kloop
    enddo !iloop
  enddo !jloop
 enddo !tloop 
else ! case when flag=1 (i.e.: rho or pressure)
 do tloop=1,tlen
  do jloop=1,latlen
    do iloop=1,lonlen
      ! preliminary stuff
      do kloop=1,altlen
        ! extract the vertical coordinates 
        z(kloop)=z_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
        if (z_new(kloop).lt.z_gcm(iloop,jloop,1,tloop)) then
          ! z_new(kloop) is below lowest GCM level
          newdata(iloop,jloop,kloop,tloop)=gcmdata(iloop,jloop,1,tloop)* &
                   exp(((z_gcm(iloop,jloop,1,tloop)-z_new(kloop))*g0)/ &
                       (Rbottom(iloop,jloop,tloop)* &
                        temp(iloop,jloop,1,tloop)))
        else if (z_new(kloop).gt.z_gcm(iloop,jloop,altlen,tloop)) then
          ! z_new(kloop) is above highest GCM level
          ! NB: zareoid=zsurface+phis/g0
          gz=g0*a0*a0/(a0+z_new(kloop)+phis(iloop,jloop)/g0)**2
          newdata(iloop,jloop,kloop,tloop)=gcmdata(iloop,jloop,altlen,tloop)* &
                   exp(((z_gcm(iloop,jloop,altlen,tloop)-z_new(kloop))*gz)/ &
                       (Rtop(iloop,jloop,tloop)* &
                        temp(iloop,jloop,altlen,tloop)))
        else ! z_new(kloop) is within range
          x=z_new(kloop)
          call interpolf(x,y,missing,z,logq,altlen)
          newdata(iloop,jloop,kloop,tloop)=exp(y)
        endif
      enddo !kloop
    enddo !iloop
  enddo !jloop
 enddo !tloop 
endif

! Cleanup
deallocate(z)
deallocate(q)
deallocate(logq)
if (flag.eq.1) then
  deallocate(Rbottom)
  deallocate(Rtop)
endif

end subroutine zs_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