program aeroptical
! program for computation of aerosol opacities
! author : Antoine Bierjon, April-May 2020
! Rebranding in November 2022
!
!===============================================================================
!     PREFACE
!===============================================================================
! This program can be used to compute the opacity of any aerosol.
!
! It takes as input a GCM output file (diagfi,stats,concat) that
! contains:
!      - the Mass Mixing Ratios of the aerosols to be computed
!      - their effective radius (replaceable by a single value given by the user)
!      - the atmospheric density
! The user also inputs the wavelength to calibrate the opacities,
! the type of opacity (extinction or absorption) and the directory
! containing the ASCII files of the aerosols' optical properties.
!
! Finally, for each aerosol to be computed, the user inputs :
! <aerosol_name>,<aerosol_mmr_varname>,<aerosol_reff_varname/value>,<aerosol_rho_value>,<optpropfile_name>
!
! The program outputs a netcdf file containing the opacities [1/km] of
! the aerosols calibrated at the prescribed wavelength.
!
! NB: this program requires to compile the module aeropt_mod.F90 at the same time
!===============================================================================


use netcdf
use aeropt_mod

!===============================================================================
! Variable declarations
!===============================================================================

implicit none ! for no implicitly typed variables

! User inputs
character(len=256) :: datadir                                        ! directory containing the ASCII optical properties files
real :: wvl_val                                                      ! reference wavelength of the output opacities [m]
character(len=3) :: opatype                                          ! opacity type : extinction (ext) or absorption (abs)
character(len=3) :: compute_col                                      ! does the user want to compute the column-integrated optical depth [yes/no]
character(len=50),dimension(20) :: name_aer,mmrname_aer,reffname_aer ! aerosols' name, MMR varname, reff varname (or value)
integer,dimension(20) :: reffval_ok                                  ! to know if the user gave a variable name or a value for reff_aer
real,dimension(20) :: reffval_aer                                    ! value given by the user for reff_aer [m]
real,dimension(20) :: rho_aer                                        ! Density of the aerosols [kg.m-3]
character(len=128),dimension(20) :: optpropfile_aer                  ! name of the ASCII optical properties file
integer :: read_ok                                                   ! to know if the user input is well read
integer :: naerkind                                                  ! nb of aerosols to consider

! GCM file
character(len=128) :: gcmfile                  ! name of the netcdf GCM input file
integer :: gcmfid                              ! [netcdf] GCM input file ID number
integer :: ierr                                ! [netcdf] subroutine returned error code
character(len=256) :: error_text               ! text to display in case of error
integer :: lonlen,latlen,altlen,timelen        ! nb of grid points along longitude,latitude,altitude,time
integer :: GCM_layers                          ! number of GCM layers
integer :: layerdimout,interlayerdimout        ! "GCM_layers" and "GCM_layers+1" IDs

logical,dimension(:),allocatable :: aerok      ! to know if the needed fields are in the GCM file
integer,dimension(:),allocatable :: mmrvarid   ! stores a MMR variable ID number
integer,dimension(:),allocatable :: reffvarid  ! stores a reff variable ID number
integer :: tmpvarid                            ! temporarily stores a variable ID number
real,dimension(:,:,:,:),allocatable :: mmr     ! aerosols mass mixing ratios [kg/kg]
real,dimension(:,:,:,:),allocatable :: reff    ! aerosols effective radii [m]
real,dimension(:,:,:,:),allocatable :: rho     ! atmospheric density [kg.m-3]
integer :: iaer                                ! Loop index for aerosols
integer :: ilon,ilat,ialt,it                   ! Loop indices for coordinates

! Opacity computation
character(len=128) :: aeroptlogfile            ! name of the aeropt_mod log file
real :: Qext_val                               ! Qext after interpolations
real :: omeg_val                               ! omeg after interpolations

! Output file
character(len=128) :: outfile                  ! name of the netcdf output file
integer :: outfid                              ! [netcdf] file ID number
integer :: latdimout,londimout,altdimout,timedimout
                                               ! latdimout: stores latitude dimension ID number
                                               ! londimout: stores longitude dimension ID number
                                               ! altdimout: stores altitude dimension ID number
                                               ! timedimout: stores time dimension ID number
integer :: tmpvaridout                         ! temporarily stores a variable ID number
integer :: varshape(4)                         ! stores a variable shape (of dimensions' IDs)
character(len=16) :: tmpvarname                ! temporarily stores a variable name
character(len=128) :: tmplongname              ! temporarily stores a variable long_name attribute
real,dimension(:,:,:,:),allocatable :: opa_aer ! Aerosols opacities dtau/dz [1/km]
real,dimension(:),allocatable :: missval       ! Value to put in outfile when the reff is out of bounds
                                               ! or when there is a mising value in input file
integer :: tauvarshape(3)                      ! stores the column variable shape (of dimensions' IDs)
real,dimension(:,:,:,:),allocatable :: delta_z ! Layers' height for column computation [km]
real,dimension(:,:,:),allocatable :: tau_aer   ! Aerosols column-integrated optical depth [/]


!===============================================================================
! 0. USER INPUTS FOR NEW AEROPTICAL.DEF

!1) Name of the GCM input file
!2) Directory where the optical properties files are
!3) The wavelength at which the output opacities are calibrated (in meters)
!4) Opacity type : extinction (ext) or absorption (abs)
!5) Does the user want to compute the column-integrated optical depth? [yes/no]
!6) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
!7) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name
!...
!N) aerosol_name,aerosol_mmr(variable name),aerosol_reff(variable name or single value in m),aerosol_rho(value in kg/m3),optpropfile_name

!===============================================================================
write(*,*) "Welcome in the aerosol opacities' computation program !"
write(*,*)

! 1) Ask the GCM file name
WRITE(*,*) "Enter an input file name (GCM diagfi/stats/concat) :"
READ(*,*) gcmfile

! 2) Ask the directory containing the optical properties files
write(*,*)""
WRITE(*,*) "In which directory should we look for the optical properties' files ?"
READ(*,'(a)') datadir

! 3) Ask the wavelength to the user
write(*,*)""
WRITE(*,*) "Value of the reference wavelength for the opacities' computation (in meters) ?"
READ(*,*) wvl_val

! 4) Ask the type of opacity that has to be computed
write(*,*)""
WRITE(*,*) "Do you want to compute extinction or absorption opacities ? (ext/abs)"
READ(*,*) opatype
if ((trim(opatype).ne."ext").and.(trim(opatype).ne."abs")) then
  write(*,*) "opatype = "//opatype
  write(*,*) "Error: the opacity type must be ext or abs"
  stop
endif

! 5) Computation of the column-integrated optical depth?
write(*,*)""
WRITE(*,*) "Do you want to compute the column-integrated optical depth ? (yes/no)"
READ(*,*) compute_col
if ((trim(compute_col).ne."yes").and.(trim(compute_col).ne."no")) then
  write(*,*) "compute_col = "//compute_col
  write(*,*) "Error: please answer yes or no"
  stop
endif


! 6-N) Ask for aerosols to compute to the user
write(*,*)""
write(*,*) "For which aerosols do you want to compute the opacity?"
write(*,*) "List of aerosols (up to 20), separated by <Enter>s, with each line containing :"
write(*,*) "<aerosol_name>,<aerosol_mmr>,<aerosol_reff>,<aerosol_rho>,<optpropfile_name>"
write(*,*) "(aerosol_mmr is a variable name, aerosol_reff can be variable name or value in meter ; aerosol_rho is a value in kg/m3)"
write(*,*)""
write(*,*) "An empty line , i.e: just <Enter>, implies end of list."

reffval_ok(:) = 0 ! initialize reffval_ok to 0 for all aerosols

naerkind=0
read(*,*,iostat=read_ok) name_aer(1),mmrname_aer(1),reffname_aer(1),rho_aer(1),optpropfile_aer(1)

do while ((read_ok.eq.0).and.(name_aer(naerkind+1)/=""))
  naerkind=naerkind+1
  
  ! Check if reffname_aer is a variable name or a value (in meters)
  read(reffname_aer(naerkind),*,iostat=reffval_ok(naerkind)) reffval_aer(naerkind)
  if (reffval_ok(naerkind).eq.0) then
    write(*,*) "This value for reff of ",trim(name_aer(naerkind))," will be broadcasted to the whole grid :",reffval_aer(naerkind)
  endif

  ! Read next line
  read(*,*,iostat=read_ok) name_aer(naerkind+1),mmrname_aer(naerkind+1),reffname_aer(naerkind+1),rho_aer(naerkind+1),optpropfile_aer(naerkind+1)
enddo

if(naerkind==0) then
   write(*,*) "no aerosol... better stop now."
   stop
endif

write(*,*) "naerkind=",naerkind
write(*,*) "name_aer=",name_aer(1:naerkind)
write(*,*) "mmrname_aer=",mmrname_aer(1:naerkind)
write(*,*) "reffname_aer=",reffname_aer(1:naerkind)
write(*,*) "rho_aer=",rho_aer(1:naerkind)
write(*,*) "optpropfile_aer=",optpropfile_aer(1:naerkind)


!===============================================================================
! 1. GCM FILE & OUTPUT FILE INITIALIZATIONS
!===============================================================================

!==========================================================================
! 1.1 GCM file opening & dimensions - Output file initializations
!==========================================================================

!==========================================================================
! --> 1.1.1 Open the netcdf GCM file given by the user

ierr=nf90_open(gcmfile,nf90_nowrite,gcmfid) ! nowrite mode=the program can only read the opened file
error_text="Error: could not open file "//trim(gcmfile)
call status_check(ierr,error_text)

!==========================================================================
! --> 1.1.2 Creation of the outfile
if (opatype.eq."ext") then
  outfile=gcmfile(1:index(gcmfile, ".nc")-1)//"_OPAext.nc"
else ! opatype.eq."abs"
  outfile=gcmfile(1:index(gcmfile, ".nc")-1)//"_OPAabs.nc"
endif

 
ierr=NF90_CREATE(outfile,IOR(nf90_clobber,nf90_64bit_offset),outfid)
! NB: clobber mode=overwrite any dataset with the same file name !
! 64bit_offset enables the creation of large netcdf files with variables>2GB
error_text="Error: could not create file "//trim(outfile)
call status_check(ierr,error_text)
write(*,*)"";WRITE(*,*)"Output file is: ",trim(outfile);write(*,*)""

!==========================================================================
! --> 1.1.3 Get the dimensions and create them in the output file

! Initialize output file's lat,lon,alt,time and controle dimensions
call inidims(gcmfile,gcmfid,outfile,outfid,&
             lonlen,latlen,altlen,timelen,&
             latdimout,londimout,altdimout,timedimout,&
             GCM_layers,layerdimout,interlayerdimout)

! Length allocation for each dimension of the 4D variable delta_z
allocate(delta_z(lonlen,latlen,altlen,timelen))

! Initialize output file's aps,bps,ap,bp and phisinit variables
! + compute delta_z for the column integration
call init2(gcmfid,lonlen,latlen,altlen,timelen,GCM_layers,&
           outfid,londimout,latdimout,altdimout,timedimout,&
           layerdimout,interlayerdimout,&
           compute_col,delta_z)

!==========================================================================
! 1.2 Check GCM aerosols
!==========================================================================

! Variables length allocation
allocate(mmrvarid(naerkind))
allocate(reffvarid(naerkind))

! Initialize missing value
allocate(missval(naerkind))
missval(:)=1.e+20

! Initialize aerok to .true. for all aerosols
allocate(aerok(naerkind))
aerok(:)=.true.

! Loop on all aerosols
DO iaer = 1, naerkind
  ! MASS MIXING RATIO
  ! Check that the GCM file contains that variable
  ierr=nf90_inq_varid(gcmfid,mmrname_aer(iaer),mmrvarid(iaer))
  error_text="Failed to find variable "//trim(mmrname_aer(iaer))//&
             " in "//trim(gcmfile)//&
             " We'll skip the "//name_aer(iaer)//" aerosol."
  if (ierr.ne.nf90_noerr) then
    write(*,*)trim(error_text)
    aerok(iaer)=.false. 
  endif

  ! EFFECTIVE RADIUS
  if (aerok(iaer)) then
    IF (reffval_ok(iaer).ne.0) THEN ! if the user gave a variable's name and not a value
      ! Check that the GCM file contains that variable
      ierr=nf90_inq_varid(gcmfid,reffname_aer(iaer),reffvarid(iaer))
      error_text="Failed to find variable "//trim(reffname_aer(iaer))//&
                 " in "//trim(gcmfile)//&
                 " We'll skip the "//trim(name_aer(iaer))//" aerosol."
      if (ierr.ne.nf90_noerr) then
        write(*,*)trim(error_text)
        aerok(iaer)=.false.
      endif
    ENDIF
  endif
ENDDO !iaer

! Check if there is still at least one aerosol to compute
IF (.NOT.ANY(aerok)) THEN
  write(*,*) "No aerosol among those listed was found in the file. Better stop now..."
  stop
ENDIF

!==========================================================================
! 1.3 Get GCM atmospheric density
!==========================================================================

! Check that the GCM file contains that variable
ierr=nf90_inq_varid(gcmfid,"rho",tmpvarid)
error_text="Failed to find variable rho in "//trim(gcmfile)&
            //" We need it to compute the opacity [1/km]."
call status_check(ierr,error_text)
! Length allocation for each dimension of the 4D variable
allocate(rho(lonlen,latlen,altlen,timelen))
! Load dataset
ierr=nf90_get_var(gcmfid,tmpvarid,rho)
error_text="Failed to load atmospheric density"
call status_check(ierr,error_text)
write(*,*) "Atmospheric density rho loaded from "//trim(gcmfile)


!===============================================================================
! 2. OUTPUT FILE VARIABLES
!===============================================================================
!==========================================================================
! 2.1 Output variable's initializations and datasets loading
!==========================================================================
! Define the ID shape of the output variables
varshape(1)=londimout
varshape(2)=latdimout
varshape(3)=altdimout
varshape(4)=timedimout

tauvarshape(1)=londimout
tauvarshape(2)=latdimout
tauvarshape(3)=timedimout


! Loop on all aerosols
DO iaer = 1, naerkind
  if (aerok(iaer)) then ! check if this aerosol can be computed
    write(*,*)""
    write(*,*)"Computing ",trim(name_aer(iaer))," opacity..."
    
    ! Length allocation of the input variables
    ALLOCATE(mmr(lonlen,latlen,altlen,timelen))
    ALLOCATE(reff(lonlen,latlen,altlen,timelen))
    
    ! Datasets loading
    ! ->MMR
    ierr=NF90_GET_VAR(gcmfid,mmrvarid(iaer),mmr(:,:,:,:))
    error_text="Failed to load mass mixing ratio"
    call status_check(ierr,error_text)
    write(*,*) "Mass mixing ratio loaded from "//trim(gcmfile)
    ! Get missing value
    ierr=nf90_get_att(gcmfid,mmrvarid(iaer),"missing_value",missval(iaer))
    if (ierr.ne.nf90_noerr) then
      missval(iaer) = 1.e+20
    endif
    
    ! ->REFF
    IF (reffval_ok(iaer).ne.0) THEN ! if the user gave a variable's name and not a value
      ! Load dataset
      ierr=NF90_GET_VAR(gcmfid,reffvarid(iaer),reff(:,:,:,:))
      error_text="Failed to load effective radius"
      call status_check(ierr,error_text)
      write(*,*) "Effective radius loaded from "//trim(gcmfile)
    ELSE ! if reffval_ok(iaer)=0
      reff(:,:,:,:) = reffval_aer(iaer)
      write(*,*) "Effective radius loaded from the user input"
    ENDIF
    
  
    ! Length allocation of the output variables
    ALLOCATE(opa_aer(lonlen,latlen,altlen,timelen))
    
    if (trim(compute_col).eq."yes") then
      ALLOCATE(tau_aer(lonlen,latlen,timelen))
      tau_aer(:,:,:) = 0.
    endif ! trim(compute_col).eq."yes"
    
    
!==========================================================================
! 2.2 Reading optical properties
!==========================================================================

    write(*,*)""
    ! create the aeropt_mod log file for this aerosol
    write(aeroptlogfile,*) "aeropt_log_",trim(name_aer(iaer)),".txt"
    aeroptlogfileID = 100+iaer ! update the log file unit
    CALL create_logfile(aeroptlogfile)

    write(*,*)"Reading the optical properties..."
    ! Read the properties tables
    CALL read_optpropfile(datadir,optpropfile_aer(iaer))

!==========================================================================
! 2.3 Creation of the output aerosol opacity variable
!==========================================================================
    ! Switch to netcdf define mode
    ierr=nf90_redef(outfid)
    error_text="ERROR: Couldn't start netcdf define mode"
    call status_check(ierr,error_text)
    write(*,*)""
    
    ! Insert the definition of the variable
    tmpvarname="opa_"//trim(name_aer(iaer))
    ierr=NF90_DEF_VAR(outfid,trim(tmpvarname),nf90_float,varshape,tmpvaridout)
    error_text="ERROR: Couldn't create "//trim(tmpvarname)//" in the outfile"
    call status_check(ierr,error_text)
    write(*,*) trim(tmpvarname)//" has been created in the outfile"

    ! Write the attributes
    if (opatype.eq."ext") then
      tmplongname=trim(name_aer(iaer))//" extinction opacity at reference wavelength"
      ierr=nf90_put_att(outfid,tmpvaridout,"long_name",tmplongname)
    else ! opatype.eq."abs"
      tmplongname=trim(name_aer(iaer))//" absorption opacity at reference wavelength"
      ierr=nf90_put_att(outfid,tmpvaridout,"long_name",tmplongname)
    endif
    
    ierr=nf90_put_att(outfid,tmpvaridout,"units","opacity/km")
    ierr=nf90_put_att(outfid,tmpvaridout,"refwavelength",wvl_val)
    ierr=nf90_put_att(outfid,tmpvaridout,"missing_value",missval(iaer))
    write(*,*)"with missing value = ",missval(iaer)
    
    ! End netcdf define mode
    ierr=nf90_enddef(outfid)
    error_text="ERROR: Couldn't end netcdf define mode"
    call status_check(ierr,error_text)

!==========================================================================
! 2.4 Computation of the opacities
!==========================================================================
    IF (reffval_ok(iaer).eq.0) THEN
    ! if the user gave a reff value,
    ! do the interpolation only once
      CALL interp_wvl_reff(wvl_val,reffval_aer(iaer),missval(iaer),&
                           Qext_val,omeg_val)
      write(*,*) "Qext(reff_val)=",Qext_val," ; omeg(reff_val)=",omeg_val
    ENDIF

    do ilon=1,lonlen
     do ilat=1,latlen
      do it=1,timelen
       do ialt=1,altlen
         IF (reffval_ok(iaer).ne.0) THEN ! if the user gave a variable's name and not a value
           CALL interp_wvl_reff(wvl_val,reff(ilon,ilat,ialt,it),missval(iaer),&
                                Qext_val,omeg_val)
         ENDIF
         
         if ((Qext_val.eq.missval(iaer))) then
           ! if the user's effective radius is
           ! out of the optpropfile boundaries
           opa_aer(ilon,ilat,ialt,it)=missval(iaer)
           write(aeroptlogfileID,*) "(mmr = ",mmr(ilon,ilat,ialt,it),"kg/kg)"
           write(aeroptlogfileID,*)""
           
         else if ((mmr(ilon,ilat,ialt,it).eq.missval(iaer))) then
           ! if there is a missing value in input file
           opa_aer(ilon,ilat,ialt,it)=missval(iaer)
           
         else
           ! COMPUTATION OF THE EXTINCTION OPACITY [1/km]
           opa_aer(ilon,ilat,ialt,it)= 750.*Qext_val*mmr(ilon,ilat,ialt,it)*rho(ilon,ilat,ialt,it)&
                                               / ( rho_aer(iaer) * reff(ilon,ilat,ialt,it) )
        
           if (opatype.eq."abs") then
              ! COMPUTATION OF THE ABSORPTION OPACITY [1/km]
              opa_aer(ilon,ilat,ialt,it)= (1 - omeg_val) * opa_aer(ilon,ilat,ialt,it)
           endif ! opatype.eq."abs"

         endif
         
         if (trim(compute_col).eq."yes") then
           ! COMPUTATION OF THE COLUMN-INTEGRATED OPTICAL DEPTH [/]
           if (opa_aer(ilon,ilat,ialt,it).ne.missval(iaer)) then
             tau_aer(ilon,ilat,it)= tau_aer(ilon,ilat,it)+opa_aer(ilon,ilat,ialt,it)*delta_z(ilon,ilat,ialt,it)
           endif
         endif ! trim(compute_col).eq."yes"
             
       enddo ! ialt
       
       if (trim(compute_col).eq."yes") then
         ! COMPUTATION OF THE COLUMN-INTEGRATED OPTICAL DEPTH [/]
         if (tau_aer(ilon,ilat,it).eq.0) then
           tau_aer(ilon,ilat,it)= missval(iaer)
         endif
       endif ! trim(compute_col).eq."yes"
       
      enddo ! it
     enddo ! ilat
    enddo ! ilon

!==========================================================================
! 3.3 Writing opacity in the output file
!==========================================================================
    ierr = NF90_PUT_VAR(outfid, tmpvaridout, opa_aer(:,:,:,:))
    error_text="Error: could not write "//trim(tmpvarname)//" data in the outfile"
    call status_check(ierr,error_text)
    
!==========================================================================
! 3.4 Writing column-integrated optical depth in the output file
!==========================================================================
    if (trim(compute_col).eq."yes") then
      ! Switch to netcdf define mode
      ierr=nf90_redef(outfid)
      error_text="ERROR: Couldn't start netcdf define mode"
      call status_check(ierr,error_text)
      write(*,*)""

      ! Insert the definition of the variable
      tmpvarname="tau_"//trim(name_aer(iaer))
      ierr=NF90_DEF_VAR(outfid,trim(tmpvarname),nf90_float,tauvarshape,tmpvaridout)
      error_text="ERROR: Couldn't create "//trim(tmpvarname)//" in the outfile"
      call status_check(ierr,error_text)
      write(*,*) trim(tmpvarname)//" has been created in the outfile"

      ! Write the attributes
      if (opatype.eq."ext") then
        tmplongname=trim(name_aer(iaer))//" extinction column-integrated optical depth at ref. wavelength"
        ierr=nf90_put_att(outfid,tmpvaridout,"long_name",tmplongname)
      else ! opatype.eq."abs"
        tmplongname=trim(name_aer(iaer))//" absorption column-integrated optical depth at ref. wavelength"
        ierr=nf90_put_att(outfid,tmpvaridout,"long_name",tmplongname)
      endif

      ierr=nf90_put_att(outfid,tmpvaridout,"units","/")
      ierr=nf90_put_att(outfid,tmpvaridout,"refwavelength",wvl_val)
      ierr=nf90_put_att(outfid,tmpvaridout,"missing_value",missval(iaer))
      write(*,*)"with missing value = ",missval(iaer)

      ! End netcdf define mode
      ierr=nf90_enddef(outfid)
      error_text="ERROR: Couldn't end netcdf define mode"
      call status_check(ierr,error_text)


      ! Writing variable in output file
      ierr = NF90_PUT_VAR(outfid, tmpvaridout, tau_aer(:,:,:))
      error_text="Error: could not write "//trim(tmpvarname)//" data in the outfile"
      call status_check(ierr,error_text)
    endif ! trim(compute_col).eq."yes"
    
    
!==========================================================================
! 3.5 Prepare next loop
!==========================================================================
    DEALLOCATE(mmr)
    DEALLOCATE(reff)
    DEALLOCATE(opa_aer)
    IF (allocated(tau_aer)) DEALLOCATE(tau_aer)
    
    CALL end_aeropt_mod
    
  endif ! if aerok(iaer)

ENDDO ! iaer

!=============================================================================== 
! 4. Close the files and end the program
!===============================================================================
ierr = nf90_close(gcmfid)
error_text="Error: could not close file "//trim(gcmfile)
call status_check(ierr,error_text)

ierr = nf90_close(outfid)
error_text="Error: could not close file "//trim(outfile)
call status_check(ierr,error_text)

write(*,*)"";write(*,*)"Program aeroptical completed!"

end program aeroptical





!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! SUBROUTINES
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine status_check(ierr,error_text)

use netcdf
implicit none
!================================================================
!  Arguments:
!================================================================
integer,intent(in) :: ierr ! NetCDF status number
character(len=256),intent(in) :: error_text

if (ierr.ne.nf90_noerr) then
  write(*,*)trim(error_text)
  write(*,*)trim(nf90_strerror(ierr))
  stop
endif

end subroutine status_check


!*******************************************************************************

subroutine inidims(gcmfile,gcmfid,outfile,outfid,lonlen,latlen,altlen,timelen,&
                   latdimout,londimout,altdimout,timedimout,&
                   GCM_layers,layerdimout,interlayerdimout)

!==============================================================================
! Purpose:
! Read the dimensions of the input file
! and write them in the output file
!==============================================================================
! Remarks:
! The NetCDF files must be open
!==============================================================================
use netcdf

implicit none

!==============================================================================
! Arguments:
!==============================================================================
character(len=128),intent(in) :: gcmfile ! name of the netcdf GCM input file
integer,intent(in) :: gcmfid ! [netcdf] GCM input file ID number
character(len=128),intent(in) :: outfile ! name of the netcdf output file
integer,intent(in) :: outfid ! [netcdf] file ID number

integer,intent(out) ::  lonlen,latlen,altlen,timelen
! nb of grid points along longitude,latitude,altitude,time
integer,intent(out) :: latdimout,londimout,altdimout,timedimout
! latdimout: stores latitude dimension ID number
! londimout: stores longitude dimension ID number
! altdimout: stores altitude dimension ID number
! timedimout: stores time dimension ID number
integer,intent(out) :: GCM_layers ! number of GCM layers
integer,intent(out) :: layerdimout ! "GCM_layers" ID
integer,intent(out) :: interlayerdimout ! "GCM_layers+1" ID

!==============================================================================
! Local variables:
!==============================================================================
integer :: ierr ! [netcdf] subroutine returned error code
character(len=256) :: error_text ! text to display in case of error

integer :: tmpdimid,tmpvarid ! temporary store a dimension/variable ID number
character (len=64) :: long_name,units,positive
! long_name(): [netcdf] long_name attribute
! units(): [netcdf] units attribute
! positive(): [netcdf] positive direction attribute (for altitude)
real, dimension(:), allocatable:: lat,lon,alt,time,ctl
! lat(): array, stores latitude coordinates
! lon(): array, stores longitude coordinates
! alt(): array, stores altitude coordinates
! time(): array, stores time coordinates
! ctl(): array, stores controle variable
integer :: ctllen ! nb of elements in the controle array
integer :: tmpdimidout,tmpvaridout ! temporary stores a dimension/variable ID number

!==============================================================================
! LONGITUDE
!==============================================================================
! Get the dimension in GCM file
ierr=nf90_inq_dimid(gcmfid,"longitude",tmpdimid)
error_text="Error: Dimension <longitude> is missing in file "//trim(gcmfile)
call status_check(ierr,error_text)

ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=lonlen)
allocate(lon(lonlen))

! Create the dimension in output file
ierr=NF90_DEF_DIM(outfid,"longitude",lonlen,londimout)
error_text="Error: could not define the longitude dimension in the outfile"
call status_check(ierr,error_text)

! Get the field in GCM file
ierr=nf90_inq_varid(gcmfid,"longitude",tmpvarid)
error_text="Error: Field <longitude> is missing in file "//trim(gcmfile)
call status_check(ierr,error_text)

ierr=NF90_GET_VAR(gcmfid,tmpvarid,lon)
error_text="Failed to load longitude"
call status_check(ierr,error_text)

! Create the field in the output file
ierr=NF90_DEF_VAR(outfid,"longitude",nf90_float,(/londimout/),tmpvaridout)
error_text="Error: could not define the longitude variable in the outfile"
call status_check(ierr,error_text)

! Get the field attributes in the GCM file
ierr=nf90_get_att(gcmfid,tmpvarid,"long_name",long_name)
if (ierr.ne.nf90_noerr) then
! if no attribute "long_name", try "title"
  ierr=nf90_get_att(gcmfid,tmpvarid,"title",long_name)
endif
ierr=nf90_get_att(gcmfid,tmpvarid,"units",units)

! Put the field attributes in the output file
ierr=nf90_put_att(outfid,tmpvaridout,"long_name",long_name)
ierr=nf90_put_att(outfid,tmpvaridout,"units",units)

! Write the field values in the output file
ierr=nf90_enddef(outfid) ! end netcdf define mode
error_text="Error: could not end the define mode of the outfile"
call status_check(ierr,error_text)

ierr=NF90_PUT_VAR(outfid,tmpvaridout,lon)
error_text="Error: could not write the longitude field in the outfile"
call status_check(ierr,error_text)

!==============================================================================
! LATITUDE
!==============================================================================
! Switch to netcdf define mode
ierr=nf90_redef(outfid)
error_text="Error: could not switch to define mode in the outfile"
call status_check(ierr,error_text)

! Get the dimension in GCM file
ierr=nf90_inq_dimid(gcmfid,"latitude",tmpdimid)
error_text="Error: Dimension <latitude> is missing in file "//trim(gcmfile)
call status_check(ierr,error_text)

ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=latlen)
allocate(lat(latlen))

! Create the dimension in output file
ierr=NF90_DEF_DIM(outfid,"latitude",latlen,latdimout)
error_text="Error: could not define the latitude dimension in the outfile"
call status_check(ierr,error_text)

! Get the field in GCM file
ierr=nf90_inq_varid(gcmfid,"latitude",tmpvarid)
error_text="Error: Field <latitude> is missing in file "//trim(gcmfile)
call status_check(ierr,error_text)

ierr=NF90_GET_VAR(gcmfid,tmpvarid,lat)
error_text="Failed to load latitude"
call status_check(ierr,error_text)

! Create the field in the output file
ierr=NF90_DEF_VAR(outfid,"latitude",nf90_float,(/latdimout/),tmpvaridout)
error_text="Error: could not define the latitude variable in the outfile"
call status_check(ierr,error_text)

! Get the field attributes in the GCM file
ierr=nf90_get_att(gcmfid,tmpvarid,"long_name",long_name)
if (ierr.ne.nf90_noerr) then
! if no attribute "long_name", try "title"
  ierr=nf90_get_att(gcmfid,tmpvarid,"title",long_name)
endif
ierr=nf90_get_att(gcmfid,tmpvarid,"units",units)

! Put the field attributes in the output file
ierr=nf90_put_att(outfid,tmpvaridout,"long_name",long_name)
ierr=nf90_put_att(outfid,tmpvaridout,"units",units)

! Write the field values in the output file
ierr=nf90_enddef(outfid) ! end netcdf define mode
error_text="Error: could not end the define mode of the outfile"
call status_check(ierr,error_text)

ierr=NF90_PUT_VAR(outfid,tmpvaridout,lat)
error_text="Error: could not write the latitude field in the outfile"
call status_check(ierr,error_text)

!==============================================================================
! ALTITUDE
!==============================================================================
! Switch to netcdf define mode
ierr=nf90_redef(outfid)
error_text="Error: could not switch to define mode in the outfile"
call status_check(ierr,error_text)

! Get the dimension in GCM file
ierr=nf90_inq_dimid(gcmfid,"altitude",tmpdimid)
error_text="Error: Dimension <altitude> is missing in file "//trim(gcmfile)
call status_check(ierr,error_text)

ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=altlen)
allocate(alt(altlen))

! Create the dimension in output file
ierr=NF90_DEF_DIM(outfid,"altitude",altlen,altdimout)
error_text="Error: could not define the altitude dimension in the outfile"
call status_check(ierr,error_text)

! Get the field in GCM file
ierr=nf90_inq_varid(gcmfid,"altitude",tmpvarid)
error_text="Error: Field <altitude> is missing in file "//trim(gcmfile)
call status_check(ierr,error_text)

ierr=NF90_GET_VAR(gcmfid,tmpvarid,alt)
error_text="Failed to load altitude"
call status_check(ierr,error_text)

! Create the field in the output file
ierr=NF90_DEF_VAR(outfid,"altitude",nf90_float,(/altdimout/),tmpvaridout)
error_text="Error: could not define the altitude variable in the outfile"
call status_check(ierr,error_text)

! Get the field attributes in the GCM file
ierr=nf90_get_att(gcmfid,tmpvarid,"long_name",long_name)
if (ierr.ne.nf90_noerr) then
! if no attribute "long_name", try "title"
  ierr=nf90_get_att(gcmfid,tmpvarid,"title",long_name)
endif
ierr=nf90_get_att(gcmfid,tmpvarid,"units",units)
ierr=nf90_get_att(gcmfid,tmpvarid,"positive",positive)

! Put the field attributes in the output file
ierr=nf90_put_att(outfid,tmpvaridout,"long_name",long_name)
ierr=nf90_put_att(outfid,tmpvaridout,"units",units)
ierr=nf90_put_att(outfid,tmpvaridout,"positive",positive)

! Write the field values in the output file
ierr=nf90_enddef(outfid) ! end netcdf define mode
error_text="Error: could not end the define mode of the outfile"
call status_check(ierr,error_text)

ierr=NF90_PUT_VAR(outfid,tmpvaridout,alt)
error_text="Error: could not write the altitude field in the outfile"
call status_check(ierr,error_text)

!==============================================================================
! TIME
!==============================================================================
! Switch to netcdf define mode
ierr=nf90_redef(outfid)
error_text="Error: could not switch to define mode in the outfile"
call status_check(ierr,error_text)

! Get the dimension in GCM file
ierr=nf90_inq_dimid(gcmfid,"Time",tmpdimid)
error_text="Error: Dimension <Time> is missing in file "//trim(gcmfile)
call status_check(ierr,error_text)

ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=timelen)
allocate(time(timelen))

! Create the dimension in output file
ierr=NF90_DEF_DIM(outfid,"Time",timelen,timedimout)
error_text="Error: could not define the time dimension in the outfile"
call status_check(ierr,error_text)

! Get the field in GCM file
ierr=nf90_inq_varid(gcmfid,"Time",tmpvarid)
error_text="Error: Field <Time> is missing in file "//trim(gcmfile)
call status_check(ierr,error_text)

ierr=NF90_GET_VAR(gcmfid,tmpvarid,time)
error_text="Failed to load Time"
call status_check(ierr,error_text)

! Create the field in the output file
ierr=NF90_DEF_VAR(outfid,"Time",nf90_float,(/timedimout/),tmpvaridout)
error_text="Error: could not define the Time variable in the outfile"
call status_check(ierr,error_text)

! Get the field attributes in the GCM file
ierr=nf90_get_att(gcmfid,tmpvarid,"long_name",long_name)
if (ierr.ne.nf90_noerr) then
! if no attribute "long_name", try "title"
  ierr=nf90_get_att(gcmfid,tmpvarid,"title",long_name)
endif
ierr=nf90_get_att(gcmfid,tmpvarid,"units",units)

! Put the field attributes in the output file
ierr=nf90_put_att(outfid,tmpvaridout,"long_name",long_name)
ierr=nf90_put_att(outfid,tmpvaridout,"units",units)

! Write the field values in the output file
ierr=nf90_enddef(outfid) ! end netcdf define mode
error_text="Error: could not end the define mode of the outfile"
call status_check(ierr,error_text)

ierr=NF90_PUT_VAR(outfid,tmpvaridout,time)
error_text="Error: could not write the Time field in the outfile"
call status_check(ierr,error_text)

!==============================================================================
! CONTROLE
!==============================================================================
! Switch to netcdf define mode
ierr=nf90_redef(outfid)
error_text="Error: could not switch to define mode in the outfile"
call status_check(ierr,error_text)

! Get the dimension in GCM file
ierr=nf90_inq_dimid(gcmfid,"index",tmpdimid)
error_text="Dimension <index> is missing in file "//trim(gcmfile)&
           //". We'll skip that one."
if (ierr.ne.nf90_noerr) then
  write(*,*)trim(error_text)
  ierr=nf90_enddef(outfid) ! end netcdf define mode
  error_text="Error: could not end the define mode of the outfile"
  call status_check(ierr,error_text)
else
  ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=ctllen)
  allocate(ctl(ctllen))

  ! Create the dimension in output file
  ierr=NF90_DEF_DIM(outfid,"index",ctllen,tmpdimidout)
  error_text="Error: could not define the index dimension in the outfile"
  call status_check(ierr,error_text)

  ! Get the field in GCM file
  ierr=nf90_inq_varid(gcmfid,"controle",tmpvarid)
  error_text="Error: Field <controle> is missing in file "//trim(gcmfile)
  call status_check(ierr,error_text)

  ierr=NF90_GET_VAR(gcmfid,tmpvarid,ctl)
  error_text="Failed to load ctl"
  call status_check(ierr,error_text)

  ! Create the field in the output file
  ierr=NF90_DEF_VAR(outfid,"controle",nf90_float,(/tmpdimidout/),tmpvaridout)
  error_text="Error: could not define the controle variable in the outfile"
  call status_check(ierr,error_text)

  ! Get the field attributes in the GCM file
  ierr=nf90_get_att(gcmfid,tmpvarid,"long_name",long_name)
  if (ierr.ne.nf90_noerr) then
  ! if no attribute "long_name", try "title"
    ierr=nf90_get_att(gcmfid,tmpvarid,"title",long_name)
  endif

  ! Put the field attributes in the output file
  ierr=nf90_put_att(outfid,tmpvaridout,"long_name",long_name)

  ! Write the field values in the output file
  ierr=nf90_enddef(outfid) ! end netcdf define mode
  error_text="Error: could not end the define mode of the outfile"
  call status_check(ierr,error_text)

  ierr=NF90_PUT_VAR(outfid,tmpvaridout,ctl)
  error_text="Error: could not write the controle field in the outfile"
  call status_check(ierr,error_text)
endif


!==============================================================================
! Load size of aps() or sigma() (in case it is not altlen)
!==============================================================================
! Switch to netcdf define mode
ierr=nf90_redef(outfid)
error_text="Error: could not switch to define mode in the outfile"
call status_check(ierr,error_text)

! Default is that GCM_layers=altlen
! but for outputs of zrecast, it may be a different value
ierr=nf90_inq_dimid(gcmfid,"GCM_layers",tmpdimid)
if (ierr.ne.nf90_noerr) then
  ! didn't find a GCM_layers dimension;
  ! we look for interlayer dimension
  ierr=nf90_inq_dimid(gcmfid,"interlayer",tmpdimid)
  if (ierr.eq.nf90_noerr) then
    ! load value of GCM_layers
    ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=GCM_layers)
    GCM_layers=GCM_layers-1
  else
    ! didn't find a GCM_layers or interlayer dimension; therefore we have:
    GCM_layers=altlen
  endif 
else
  ! load value of GCM_layers
  ierr=nf90_inquire_dimension(gcmfid,tmpdimid,len=GCM_layers)
endif

! Create the dimensions in output file
ierr = NF90_DEF_DIM(outfid,"GCM_layers",GCM_layers,layerdimout)
error_text="Error: could not define the GCM_layers dimension in the outfile"
call status_check(ierr,error_text)
ierr = NF90_DEF_DIM(outfid,"GCM_interlayers",GCM_layers+1,interlayerdimout)
error_text="Error: could not define the GCM_interlayers dimension in the outfile"
call status_check(ierr,error_text)

! End netcdf define mode
ierr=nf90_enddef(outfid)
error_text="Error: could not end the define mode of the outfile"
call status_check(ierr,error_text)

end subroutine inidims


!*******************************************************************************

subroutine init2(gcmfid,lonlen,latlen,altlen,timelen,GCM_layers, &
                 outfid,londimout,latdimout,altdimout,timedimout, &
                 layerdimout,interlayerdimout,&
                 compute_col,delta_z)
!==============================================================================
! Purpose:
! Copy ap() , bp(), aps(), bps(), aire() and phisinit()
! from input file to output file
! + compute layers' height if needed
!==============================================================================
! Remarks:
! The NetCDF files must be open
!==============================================================================

use netcdf

implicit none

!==============================================================================
! Arguments:
!==============================================================================
integer, intent(in) :: gcmfid               ! NetCDF output file ID
integer, intent(in) :: lonlen               ! # of grid points along longitude
integer, intent(in) :: latlen               ! # of grid points along latitude
integer, intent(in) :: altlen               ! # of grid points along altitude
integer, intent(in) :: timelen              ! # of grid points along time
integer, intent(in) :: GCM_layers           ! # of GCM atmospheric layers
integer, intent(in) :: outfid               ! NetCDF output file ID
integer, intent(in) :: londimout            ! longitude dimension ID
integer, intent(in) :: latdimout            ! latitude dimension ID
integer, intent(in) :: altdimout            ! altitude dimension ID
integer, intent(in) :: timedimout           ! time dimension ID
integer, intent(in) :: layerdimout          ! GCM_layers dimension ID
integer, intent(in) :: interlayerdimout     ! GCM_layers+1 dimension ID
character(len=3), intent(in) :: compute_col ! does the user want to compute the column-integrated optical depth [yes/no]

real,dimension(lonlen,latlen,altlen,timelen), intent(out) :: delta_z ! altitude difference between interlayers [km]

!==============================================================================
! Local variables:
!==============================================================================
real,dimension(:),allocatable :: aps,bps    ! hybrid vertical coordinates
real,dimension(:),allocatable :: ap,bp      ! hybrid vertical coordinates
real,dimension(:),allocatable :: sigma      ! sigma levels
real,dimension(:,:),allocatable :: aire     ! mesh areas
real,dimension(:,:),allocatable :: phisinit ! Ground geopotential

integer :: ierr                  ! [netcdf] subroutine returned error code
character(len=256) :: error_text ! text to display in case of error
integer :: tmpvarid              ! temporary variable ID

logical :: area   ! is "aire" available ?
logical :: phis   ! is "phisinit" available ?
logical :: hybrid ! are "aps" and "bps" available ?
logical :: apbp   ! are "ap" and "bp" available ? + are they usable for delta_z computation ?


! for delta_z computation
real,dimension(:,:,:,:),allocatable :: zlev ! altitude of the interlayers [km]
logical :: is_zlev                          ! is zzlev available ?
integer,dimension(4) :: zlevdimids          ! dimensions of zzlev in GCM file
integer :: zlevvertlen                      ! length of the vertical dimension of zlev

real,dimension(altlen) :: alt               ! altitude coordinates
character (len=64) :: altlong_name,altunits ! altitude long_name and units attributes

real,dimension(:,:,:,:),allocatable :: plev ! pressure of the interlayers [km]
real,dimension(:,:,:),allocatable :: ps     ! surface pressure [Pa]
real,dimension(:,:,:,:),allocatable :: rho  ! atmospheric density [kg/m3]
real,parameter :: g0 =3.7257964             ! Mars gravity [m.s-2]

real,dimension(:,:,:,:),allocatable :: zlay ! altitude of the layers [km]

integer :: ilon,ilat,it,ialt                ! loop indices
real :: missval                             ! GCM variables missing value attribute

!==============================================================================
! 1. Read data from input file
!==============================================================================

! hybrid coordinate aps
!write(*,*) "aps: altlen=",altlen," GCM_layers=",GCM_layers
allocate(aps(GCM_layers),stat=ierr)
if (ierr.ne.0) then
  write(*,*) "init2: failed to allocate aps!"
  stop
endif
ierr=nf90_inq_varid(gcmfid,"aps",tmpvarid)
if (ierr.ne.nf90_noerr) then
  write(*,*) "Ooops. Failed to get aps ID. OK, will look for sigma coord."
  hybrid=.false.
else
  ierr=NF90_GET_VAR(gcmfid,tmpvarid,aps)
  hybrid=.true.
  if (ierr.ne.nf90_noerr) then
    write(*,*) "init2 Error: Failed reading aps"
    stop
  endif

  ! hybrid coordinate bps
!  write(*,*) "bps: altlen=",altlen," GCM_layers=",GCM_layers
  allocate(bps(GCM_layers),stat=ierr)
  if (ierr.ne.0) then
    write(*,*) "init2: failed to allocate bps!"
    stop
  endif
  ierr=nf90_inq_varid(gcmfid,"bps",tmpvarid)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "init2 Error: Failed to get bps ID."
    stop
  endif
  ierr=NF90_GET_VAR(gcmfid,tmpvarid,bps)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "init2 Error: Failed reading bps"
    stop
  endif
endif

! hybrid coordinate ap
allocate(ap(GCM_layers+1),stat=ierr)
if (ierr.ne.0) then
  write(*,*) "init2: failed to allocate ap!"
  stop
else
  ierr=nf90_inq_varid(gcmfid,"ap",tmpvarid)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "Ooops. Failed to get ap ID. OK."
    apbp=.false.
  else
    ierr=NF90_GET_VAR(gcmfid,tmpvarid,ap)
    apbp=.true.
    if (ierr.ne.nf90_noerr) then
      write(*,*) "Error: Failed reading ap"
      stop
    endif
  endif
endif

! hybrid coordinate bp
allocate(bp(GCM_layers+1),stat=ierr)
if (ierr.ne.0) then
  write(*,*) "init2: failed to allocate bp!"
  stop
else
  ierr=nf90_inq_varid(gcmfid,"bp",tmpvarid)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "Ooops. Failed to get bp ID. OK."
    apbp=.false.
  else
    ierr=NF90_GET_VAR(gcmfid,tmpvarid,bp)
    apbp=.true.
    if (ierr.ne.nf90_noerr) then
      write(*,*) "Error: Failed reading bp"
      stop
    endif
  endif
endif

! sigma levels (if any)
if (.not.hybrid) then
  allocate(sigma(GCM_layers),stat=ierr)
  if (ierr.ne.0) then
    write(*,*) "init2: failed to allocate sigma"
    stop
  endif
  ierr=nf90_inq_varid(gcmfid,"sigma",tmpvarid)
  ierr=NF90_GET_VAR(gcmfid,tmpvarid,sigma)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "init2 Error: Failed reading sigma"
    stop
  endif
endif ! of if (.not.hybrid)

! mesh area
allocate(aire(lonlen,latlen),stat=ierr)
if (ierr.ne.0) then
  write(*,*) "init2: failed to allocate aire!"
  stop
endif
ierr=nf90_inq_varid(gcmfid,"aire",tmpvarid)
if (ierr.ne.nf90_noerr) then
  write(*,*)"init2 warning: Failed to get aire ID."
  area = .false.
else
  ierr=NF90_GET_VAR(gcmfid,tmpvarid,aire)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "init2 Error: Failed reading aire"
    stop
  endif
  area = .true.
endif

! ground geopotential phisinit
allocate(phisinit(lonlen,latlen),stat=ierr)
if (ierr.ne.0) then
  write(*,*) "init2: failed to allocate phisinit!"
  stop
endif
ierr=nf90_inq_varid(gcmfid,"phisinit",tmpvarid)
if (ierr.ne.nf90_noerr) then
  write(*,*)"init2 warning: Failed to get phisinit ID."
  phis = .false.
else
  ierr=NF90_GET_VAR(gcmfid,tmpvarid,phisinit)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "init2 Error: Failed reading phisinit"
    stop
  endif
  phis = .true.
endif

!==============================================================================
! 2. Write
!==============================================================================

!=========================================================================
! 2.1 Hybrid coordinates ap(), bp(), aps() and bps()
!=========================================================================
if(hybrid) then 
! define aps
  ! Switch to netcdf define mode
  ierr=nf90_redef(outfid)
  ! Insert the definition of the variable
  ierr=NF90_DEF_VAR(outfid,"aps",nf90_float,(/layerdimout/),tmpvarid)
  if (ierr.ne.nf90_noerr) then
     write(*,*) "init2 Error: Failed to define the variable aps"
     stop
  endif
  ! Write the attributes
  ierr=nf90_put_att(outfid,tmpvarid,"long_name","hybrid pressure at midlayers")
  ierr=nf90_put_att(outfid,tmpvarid,"units"," ")
  ! End netcdf define mode
  ierr=nf90_enddef(outfid)

! write aps
  ierr=NF90_PUT_VAR(outfid,tmpvarid,aps)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "init2 Error: Failed to write aps"
    stop
  endif

! define bps
  ! Switch to netcdf define mode
  ierr=nf90_redef(outfid)
  ! Insert the definition of the variable
  ierr=NF90_DEF_VAR(outfid,"bps",nf90_float,(/layerdimout/),tmpvarid)
  if (ierr.ne.nf90_noerr) then
     write(*,*) "init2 Error: Failed to define the variable bps"
     stop
  endif
  ! Write the attributes
  ierr=nf90_put_att(outfid,tmpvarid,"long_name","hybrid sigma at midlayers")
  ierr=nf90_put_att(outfid,tmpvarid,"units"," ")
  ! End netcdf define mode
  ierr=nf90_enddef(outfid)
    
! write bps
  ierr=NF90_PUT_VAR(outfid,tmpvarid,bps)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "init2 Error: Failed to write bps"
    stop
  endif

  if (apbp) then
! define ap

    ! Switch to netcdf define mode
    ierr=nf90_redef(outfid)
    ! Insert the definition of the variable
    ierr=NF90_DEF_VAR(outfid,"ap",nf90_float,(/interlayerdimout/),tmpvarid)
    if (ierr.ne.nf90_noerr) then
       write(*,*) "init2 Error: Failed to define the variable ap"
       stop
    endif
    ! Write the attributes
    ierr=nf90_put_att(outfid,tmpvarid,"long_name","hybrid sigma at interlayers")
    ierr=nf90_put_att(outfid,tmpvarid,"units"," ")
    ! End netcdf define mode
    ierr=nf90_enddef(outfid)

! write ap
    ierr=NF90_PUT_VAR(outfid,tmpvarid,ap)
    if (ierr.ne.nf90_noerr) then
      write(*,*) "Error: Failed to write ap"
      stop
    endif

! define bp

    ! Switch to netcdf define mode
    ierr=nf90_redef(outfid)
    ! Insert the definition of the variable
    ierr=NF90_DEF_VAR(outfid,"bp",nf90_float,(/interlayerdimout/),tmpvarid)
    if (ierr.ne.nf90_noerr) then
       write(*,*) "init2 Error: Failed to define the variable bp"
       stop
    endif
    ! Write the attributes
    ierr=nf90_put_att(outfid,tmpvarid,"long_name","hybrid sigma at interlayers")
    ierr=nf90_put_att(outfid,tmpvarid,"units"," ")
    ! End netcdf define mode
    ierr=nf90_enddef(outfid)
  
! write bp
    ierr=NF90_PUT_VAR(outfid,tmpvarid,bp)
    if (ierr.ne.nf90_noerr) then
      write(*,*) "Error: Failed to write bp"
      stop
    endif
  endif ! of if (apbp)

else

  ! Switch to netcdf define mode
  ierr=nf90_redef(outfid)
  ! Insert the definition of the variable
  ierr=NF90_DEF_VAR(outfid,"sigma",nf90_float,(/layerdimout/),tmpvarid)
  if (ierr.ne.nf90_noerr) then
     write(*,*) "init2 Error: Failed to define the variable sigma"
     stop
  endif
  ! Write the attributes
  ierr=nf90_put_att(outfid,tmpvarid,"long_name","sigma at midlayers")
  ierr=nf90_put_att(outfid,tmpvarid,"units"," ")
  ! End netcdf define mode
  ierr=nf90_enddef(outfid)
  
! write sigma
  ierr=NF90_PUT_VAR(outfid,tmpvarid,sigma)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "init2 Error: Failed to write sigma"
    stop
  endif
endif ! of if (hybrid)

!=========================================================================
! 2.2 aire() and phisinit()
!=========================================================================

if (area) then

  ! Switch to netcdf define mode
  ierr=nf90_redef(outfid)
  ! Insert the definition of the variable
  ierr=NF90_DEF_VAR(outfid,"aire",nf90_float,(/londimout,latdimout/),tmpvarid)
  if (ierr.ne.nf90_noerr) then
     write(*,*) "init2 Error: Failed to define the variable aire"
     stop
  endif
  ! Write the attributes
  ierr=nf90_put_att(outfid,tmpvarid,"long_name","Mesh area")
  ierr=nf90_put_att(outfid,tmpvarid,"units","m2")
  ! End netcdf define mode
  ierr=nf90_enddef(outfid)
  
  ! write aire
  ierr=NF90_PUT_VAR(outfid,tmpvarid,aire)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "init2 Error: Failed to write aire"
    stop
  endif
endif ! of if (area)

if (phis) then

  ! Switch to netcdf define mode
  ierr=nf90_redef(outfid)
  ! Insert the definition of the variable
  ierr=NF90_DEF_VAR(outfid,"phisinit",nf90_float,(/londimout,latdimout/),tmpvarid)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "init2 Error: Failed to define the variable phisinit"
    stop
  endif
  ! Write the attributes
  ierr=nf90_put_att(outfid,tmpvarid,"long_name","Ground level geopotential")
  ierr=nf90_put_att(outfid,tmpvarid,"units"," ")
  ! End netcdf define mode
  ierr=nf90_enddef(outfid)

  ! write phisinit
  ierr=NF90_PUT_VAR(outfid,tmpvarid,phisinit)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "init2 Error: Failed to write phisinit"
    stop
  endif

endif ! of if (phis)


!==============================================================================
! 3. Compute delta_z for column integration
!==============================================================================
if (trim(compute_col).eq."yes") then
  write(*,*) "Computing the layers' height delta_z..."
  
!=========================================================================
! 3.1 CASE 1: the GCM file contains zzlev variable
!=========================================================================
  ! Does the field exist in GCM file?
  is_zlev=.false.
  ierr=nf90_inq_varid(gcmfid,"zzlev",tmpvarid)
  
  IF (ierr.eq.nf90_noerr) THEN
    ! zzlev is present in the GCM file
    is_zlev=.true.
    write(*,*) "zzlev was found in GCM file"
    
    ! Get the field's dimensions
    ierr=nf90_inquire_variable(gcmfid,tmpvarid,dimids=zlevdimids)
    
    ! Get the length of the vertical dimension (assuming vertical dim is #3 like other variables)
    ierr=nf90_inquire_dimension(gcmfid,zlevdimids(3),len=zlevvertlen)
    
    ! Get missing value
    ierr=nf90_get_att(gcmfid,tmpvarid,"missing_value",missval)
    if (ierr.ne.nf90_noerr) then
      missval = 1.e+20
    endif
    
    if (zlevvertlen.eq.altlen+1) then
      ! we have every interlayer altitudes, even the top one
      allocate(zlev(lonlen,latlen,altlen+1,timelen))

      ! Get zlev variable
      ierr=NF90_GET_VAR(gcmfid,tmpvarid,zlev)
      error_text="Failed to load zzlev"
      call status_check(ierr,error_text)
      
      ! Compute delta_z [km]
      ! NB : the unit for zzlev variable is "m"
      do ilon=1,lonlen
       do ilat=1,latlen
        do it=1,timelen
         do ialt=1,altlen
         
          if ((zlev(ilon,ilat,ialt+1,it).eq.missval).or.(zlev(ilon,ilat,ialt,it).eq.missval)) then
            delta_z(ilon,ilat,ialt,it) = 0
          else
            delta_z(ilon,ilat,ialt,it) = (zlev(ilon,ilat,ialt+1,it)-zlev(ilon,ilat,ialt,it)) /1000
          endif ! missval
           
         enddo ! ialt=1,altlen
        enddo ! it=1,timelen
       enddo ! ilat=1,latlen
      enddo ! ilon=1,lonlen
      
    else if (zlevvertlen.eq.altlen) then
      ! we have every interlayer altitudes, except the top one
      allocate(zlev(lonlen,latlen,altlen,timelen))

      ! Get zlev variable
      ierr=NF90_GET_VAR(gcmfid,tmpvarid,zlev)
      error_text="Failed to load zzlev"
      call status_check(ierr,error_text)
      
      ! Compute delta_z [km]
      ! NB : the unit for zrecasted altitude is "m"
      do ilon=1,lonlen
       do ilat=1,latlen
        do it=1,timelen
         do ialt=1,altlen-1
         
          if ((zlev(ilon,ilat,ialt+1,it).eq.missval).or.(zlev(ilon,ilat,ialt,it).eq.missval)) then
            delta_z(ilon,ilat,ialt,it) = 0
          else
            delta_z(ilon,ilat,ialt,it) = (zlev(ilon,ilat,ialt+1,it)-zlev(ilon,ilat,ialt,it)) /1000
          endif ! missval
           
         enddo ! ialt=1,altlen-1
          
         ! top of last layer : we assume the same height than the layer below
         delta_z(ilon,ilat,altlen,it) = delta_z(ilon,ilat,altlen-1,it)
          
        enddo ! it=1,timelen
       enddo ! ilat=1,latlen
      enddo ! ilon=1,lonlen
            
      
      write(*,*) "delta_z has been well computed from variable zzlev"
      write(*,*)""
    
    else
      ! weird case where the GCM file would have been zrecasted but not the zzlev variable?
      ! anyway, we can't use zzlev then
      write(*,*) "zzlev and altitude of GCM file don't have the same dimension."
      write(*,*) "We will try to use the altitude dimension to compute the layers' height."
      write(*,*)""
      is_zlev=.false.

    endif ! (zlevvertlen.eq.altlen+1)
      
  ENDIF ! (ierr.eq.nf90_noerr) = zzlev present in file

!=========================================================================
! 3.2 CASE 2: try to use the GCM altitude dimension
!=========================================================================
  IF (is_zlev.eq..false.) THEN
    ! Get the field in GCM file
    ierr=nf90_inq_varid(gcmfid,"altitude",tmpvarid)

    ierr=NF90_GET_VAR(gcmfid,tmpvarid,alt)
    error_text="Failed to load altitude"
    call status_check(ierr,error_text)

    ! Get the field attributes in the GCM file
    ierr=nf90_get_att(gcmfid,tmpvarid,"long_name",altlong_name)
    if (ierr.ne.nf90_noerr) then
    ! if no attribute "long_name", try "title"
      ierr=nf90_get_att(gcmfid,tmpvarid,"title",altlong_name)
    endif
    ierr=nf90_get_att(gcmfid,tmpvarid,"units",altunits)
    
    IF (trim(altlong_name).eq."pseudo-alt") THEN
      ! GCM file has not been zrecasted, so we can use ap,bp
      ! and compute delta_z from the layer mass
      ! we need ps and rho to use this method
      
      ! Check that ap,bp have the same dimension length than altlen+1
      if (GCM_layers+1.ne.altlen+1) then
        apbp=.false.
      endif  
      
      ! Load Psurf to reconstruct pressure
      ! Does the field exist in GCM file?
      ierr=nf90_inq_varid(gcmfid,"ps",tmpvarid)
      if (ierr.eq.nf90_noerr) then
        allocate(ps(lonlen,latlen,timelen))
        ! Get ps variable
        ierr=NF90_GET_VAR(gcmfid,tmpvarid,ps)
        error_text="Failed to load surface pressure ps"
        call status_check(ierr,error_text)
      else
        apbp=.false.
      endif

      ! Load rho to compute delta_z from the layer mass
      ! Does the field exist in GCM file?
      ierr=nf90_inq_varid(gcmfid,"rho",tmpvarid)
      if (ierr.eq.nf90_noerr) then
        allocate(rho(lonlen,latlen,altlen,timelen))
        ! Get temp variable
        ierr=NF90_GET_VAR(gcmfid,tmpvarid,rho)
        error_text="Failed to load atmospheric rho"
        call status_check(ierr,error_text)
      else
        apbp=.false.
      endif
        
      if (apbp) then
        allocate(plev(lonlen,latlen,altlen+1,timelen)) 
        ! Compute delta_z [km]
        do ilon=1,lonlen
         do ilat=1,latlen
          do it=1,timelen
           ! plev at first level
           plev(ilon,ilat,1,it) = ap(1)+bp(1)*ps(ilon,ilat,it)
           
           do ialt=1,altlen
              ! plev just above
              plev(ilon,ilat,ialt+1,it) = ap(ialt+1)+bp(ialt+1)*ps(ilon,ilat,it)
              ! delta_z [km] from layer mass
              delta_z(ilon,ilat,ialt,it) = (plev(ilon,ilat,ialt,it)-plev(ilon,ilat,ialt+1,it)) &
                                           / (g0*rho(ilon,ilat,ialt,it)) /1000

           enddo ! ialt=1,altlen
           
          enddo ! it=1,timelen
         enddo ! ilat=1,latlen
        enddo ! ilon=1,lonlen

        write(*,*) "delta_z has been well computed from variables ap,bp,ps,rho"
        write(*,*)""        
        
      else
        ! no ap,bp in GCM file
        ! => take the middle of the layer altitudes as the interlayer altitude
        ! NB : the unit for "pseudo-alt" altitude (not zrecasted) is "km"
        do ialt=2,altlen-1
          delta_z(:,:,ialt,:) = (alt(ialt+1)-alt(ialt-1))/2
        enddo
        delta_z(:,:,1,:) = alt(2)/2 ! bottom at 0km (ok since pseudo-alt can't be negative)
        delta_z(:,:,altlen,:) = delta_z(:,:,altlen-1,:) ! top of last layer : we assume the same height than the layer below

        write(*,*) "delta_z has been well computed from dimension pseudo-altitude"
        write(*,*)""
        
      endif ! apbp
      
    ELSE IF (trim(altunits).eq."m") THEN
      ! GCM file has been zrecasted in altitude above surface/areoid/center of planet
      ! NB : the unit for zrecasted altitude is "m"
      
      ! Take the arithmetic mean of the layer altitudes as the interlayer altitude
      do ialt=2,altlen-1
        delta_z(:,:,ialt,:) = (alt(ialt+1)-alt(ialt-1))/2 /1000
      enddo
      
      ! Bottom layer : we can't assume a bottom of the first layer at 0km
      ! since "altitude above areoid" can be negative
      ! We thus assume zlev(2)-zlev(1) = zlay(2)-zlay(1)
      delta_z(:,:,1,:) = (alt(2)-alt(1)) /1000
      
      ! Top of last layer : we assume the same height than the layer below
      delta_z(:,:,altlen,:) = delta_z(:,:,altlen-1,:)

      write(*,*) "delta_z has been well computed from dimension altitude"
      write(*,*)""
    
    ELSE IF (trim(altunits).eq."Pa") THEN
      ! GCM file has been zrecasted in pressure
      
      ! Check if "zareoid" exist in GCM file
      ierr=nf90_inq_varid(gcmfid,"zareoid",tmpvarid)
      if (ierr.eq.nf90_noerr) then
        allocate(zlay(lonlen,latlen,altlen,timelen))
        ! Get zareoid variable
        ierr=NF90_GET_VAR(gcmfid,tmpvarid,zlay)
        error_text="Failed to load zareoid"
        call status_check(ierr,error_text)
        
        ! Get missing value
        ierr=nf90_get_att(gcmfid,tmpvarid,"missing_value",missval)
        if (ierr.ne.nf90_noerr) then
          missval = 1.e+20
        endif
        
        ! Take the arithmetic mean of the layer altitudes as the interlayer altitude
        ! NB : the unit for zareoid is "m"
        ! Compute delta_z [km]
        do ilon=1,lonlen
         do ilat=1,latlen
          do it=1,timelen
          
           ! Bottom layer
           if ((zlay(ilon,ilat,2,it).eq.missval).or.(zlay(ilon,ilat,1,it).eq.missval)) then
             delta_z(ilon,ilat,1,it) = 0
           else
             ! We can't assume a bottom of the first layer at 0km
             ! since "altitude above areoid" can be negative
             ! We thus assume zlev(2)-zlev(1) = zlay(2)-zlay(1)
             delta_z(:,:,1,:) = (zlay(ilon,ilat,2,it)-zlay(ilon,ilat,1,it)) /1000
           endif ! missval
           
           ! rest of the column
           do ialt=2,altlen-1
            if ((zlay(ilon,ilat,ialt+1,it).eq.missval).or.(zlay(ilon,ilat,ialt-1,it).eq.missval)) then
              delta_z(ilon,ilat,ialt,it) = 0
            else
              delta_z(ilon,ilat,ialt,it) = (zlay(ilon,ilat,ialt+1,it)-zlay(ilon,ilat,ialt-1,it))/2 /1000
            endif ! missval
           enddo ! ialt=2,altlen-1
           
           ! Top of last layer : we assume the same height than the layer below
           delta_z(ilon,ilat,altlen,it) = delta_z(ilon,ilat,altlen-1,it)
           
          enddo ! it=1,timelen
         enddo ! ilat=1,latlen
        enddo ! ilon=1,lonlen

        write(*,*) "delta_z has been well computed from variable zareoid"
        write(*,*)""
        
      else
      
        write(*,*) "No variable zareoid present in GCM file"
        write(*,*) "Computing the layers' height with the layers' pressure coordinate"
        write(*,*) "is not (yet?) handled by the program."
        write(*,*) "The column-integrated optical depths can thus not be computed."
        write(*,*) "Better stop now..."
        write(*,*) "You may retry with putting 'no' for the column computation."
        stop
        
      endif ! (ierr.eq.nf90_noerr) = zareoid present in file
      

    ELSE
      ! abort the program
      write(*,*) "delta_z can not be computed from the GCM file's variables."
      write(*,*) "The column-integrated optical depths can thus not be computed."
      write(*,*) "Better stop now..."
      write(*,*) "You may retry with putting 'no' for the column computation."
      stop

    ENDIF ! trim(altlong_name).eq."pseudo-alt"
    
  ENDIF ! is_zlev.eq.false
  
!=========================================================================
! 3.3 Write delta_z in output file
!=========================================================================
  
  ! Switch to netcdf define mode
  ierr=nf90_redef(outfid)
  ! Insert the definition of the variable
  ierr=NF90_DEF_VAR(outfid,"delta_z",nf90_float,(/londimout,latdimout,altdimout,timedimout/),tmpvarid)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "init2 Error: Failed to define the variable delta_z"
    stop
  endif
  ! Write the attributes
  ierr=nf90_put_att(outfid,tmpvarid,"long_name","Layer height used for column optical depth computation")
  ierr=nf90_put_att(outfid,tmpvarid,"units","km")
  ! End netcdf define mode
  ierr=nf90_enddef(outfid)

  ! write delta_z
  ierr=NF90_PUT_VAR(outfid,tmpvarid,delta_z)
  if (ierr.ne.nf90_noerr) then
    write(*,*) "init2 Error: Failed to write delta_z"
    stop
  endif
  
endif ! trim(compute_col).eq."yes"

!==============================================================================
! 4. Cleanup
!==============================================================================

if (allocated(aps)) deallocate(aps)
if (allocated(bps)) deallocate(bps)
if (allocated(ap)) deallocate(ap)
if (allocated(bp)) deallocate(bp)
if (allocated(sigma)) deallocate(sigma)
if (allocated(phisinit)) deallocate(phisinit)
if (allocated(aire)) deallocate(aire)

end subroutine init2


!*******************************************************************************