source: trunk/LMDZ.COMMON/libf/evolution/iostart_pem.F90 @ 4003

MODULE iostart_pem
!-----------------------------------------------------------------------
! NAME
!     iostart_pem
!
! DESCRIPTION
!     Read and write specific netcdf files for the PEM (start and restart)
!     Provides interfaces for field and variable I/O operations
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!     Inspired by iostart from the PCM.
!
!-----------------------------------------------------------------------

    ! DECLARATION
! -----------
implicit none

! MODULE VARIABLES
! ----------------
PRIVATE
    INTEGER,SAVE :: nid_start ! NetCDF file identifier for startfi.nc file
    INTEGER,SAVE :: nid_restart ! NetCDF file identifier for restartfi.nc file

    ! restartfi.nc file dimension identifiers: (see open_restartphy())
    INTEGER,SAVE :: idim1  ! "index" dimension
    INTEGER,SAVE :: idim2  ! "physical_points" dimension
    INTEGER,SAVE :: idim3  ! "subsurface_layers" dimension
    INTEGER,SAVE :: idim4  ! "nlayer_plus_1" dimension
    INTEGER,SAVE :: idim5  ! "number_of_advected_fields" dimension
    INTEGER,SAVE :: idim6  ! "nlayer" dimension
    INTEGER,SAVE :: idim7  ! "Time" dimension
    INTEGER,SAVE :: idim8  ! "nslope" dimension
    INTEGER,SAVE :: idim9  ! "inter slope" dimension
    INTEGER,SAVE :: idim10 ! "nb_str_max" dimension
    INTEGER,SAVE :: timeindex ! current time index (for time-dependent fields)
    INTEGER,PARAMETER :: length=100 ! size of tab_cntrl array

    INTERFACE get_field
      MODULE PROCEDURE Get_field_r1,Get_field_r2,Get_field_r3
    END INTERFACE get_field

    INTERFACE get_var
      MODULE PROCEDURE get_var_r0,Get_var_r1,Get_var_r2,Get_var_r3
    END INTERFACE get_var

    INTERFACE put_field
      MODULE PROCEDURE put_field_r1,put_field_r2,put_field_r3
    END INTERFACE put_field

    INTERFACE put_var
      MODULE PROCEDURE put_var_r0,put_var_r1,put_var_r2,put_var_r3
    END INTERFACE put_var

    PUBLIC nid_start, length
    PUBLIC get_field,get_var,put_field,put_var
    PUBLIC inquire_dimension, inquire_dimension_length
    PUBLIC inquire_field, inquire_field_ndims
    PUBLIC open_startphy,close_startphy,open_restartphy,close_restartphy

CONTAINS
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!=======================================================================
  SUBROUTINE open_startphy(filename)
!-----------------------------------------------------------------------
! NAME
!     open_startphy
!
! DESCRIPTION
!     Open the startphy netCDF file for reading.

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  USE netcdf,             only: NF90_OPEN, NF90_NOERR, NF90_NOWRITE, nf90_strerror
  USE mod_phys_lmdz_para, only: is_master, bcast

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
    CHARACTER(LEN=*) :: filename

! LOCAL VARIABLES
! ---------------
    INTEGER :: ierr

! CODE
! ----
    IF (is_master) THEN
      ierr = NF90_OPEN (filename, NF90_NOWRITE, nid_start)
      IF (ierr.NE.NF90_NOERR) THEN
        write(*,*)'open_startphy: problem opening file '//trim(filename)
        write(*,*)trim(nf90_strerror(ierr))
        CALL abort_physic("open_startphy","Cannot open file",1)
      ENDIF
    ENDIF

    CALL bcast(nid_start) ! tell all procs about nid_start

  END SUBROUTINE open_startphy
!=======================================================================

!=======================================================================
  SUBROUTINE close_startphy
!-----------------------------------------------------------------------
! NAME
!     close_startphy
!
! DESCRIPTION
!     Close the startphy netCDF file.

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  USE netcdf,             only: NF90_CLOSE
  USE mod_phys_lmdz_para, only: is_master

! DECLARATION
! -----------
implicit none

! LOCAL VARIABLES
! ---------------
    INTEGER :: ierr

! CODE
! ----
    IF (is_master) THEN
        ierr = NF90_CLOSE (nid_start)
    ENDIF

  END SUBROUTINE close_startphy
!=======================================================================

!=======================================================================
  FUNCTION inquire_field(Field_name)
!-----------------------------------------------------------------------
! NAME
!     inquire_field
!
! DESCRIPTION
!     Check if a given field is present in the input file.

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  USE netcdf,             only: NF90_INQ_VARID, NF90_NOERR
  USE mod_phys_lmdz_para, only: is_master, bcast

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
    CHARACTER(LEN=*), INTENT(IN) :: Field_name

! LOCAL VARIABLES
! ---------------
    LOGICAL :: inquire_field
    INTEGER :: varid
    INTEGER :: ierr

! CODE
! ----
    IF (is_master) THEN
      ierr=NF90_INQ_VARID(nid_start,Field_name,varid)
      IF (ierr==NF90_NOERR) THEN
        Inquire_field=.TRUE.
      ELSE
        Inquire_field=.FALSE.
      ENDIF
    ENDIF

    CALL bcast(inquire_field)

  END FUNCTION inquire_field
!=======================================================================

!=======================================================================
  FUNCTION inquire_field_ndims(Field_name)
!-----------------------------------------------------------------------
! NAME
!     inquire_field_ndims
!
! DESCRIPTION
!     Give the number of dimensions of "Field_name" stored in the input file.

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  USE netcdf,             only: nf90_inq_varid, nf90_inquire_variable, &
                                NF90_NOERR, nf90_strerror
  USE mod_phys_lmdz_para, only: is_master, bcast

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
    CHARACTER(LEN=*), INTENT(IN) :: Field_name

! LOCAL VARIABLES
! ---------------
    INTEGER :: inquire_field_ndims
    INTEGER :: varid
    INTEGER :: ierr

! CODE
! ----
    IF (is_master) THEN
      ierr=nf90_inq_varid(nid_start,Field_name,varid)
      ierr=nf90_inquire_variable(nid_start,varid,&
                                  ndims=inquire_field_ndims)
      IF (ierr.NE.NF90_NOERR) THEN
        write(*,*)'inquire_field_ndims: problem obtaining ndims of '&
                  //trim(field_name)
        write(*,*)trim(nf90_strerror(ierr))
        CALL abort_physic("inquire_field_ndims","Failed to get ndims",1)
      ENDIF
    ENDIF

    CALL bcast(inquire_field_ndims)

  END FUNCTION inquire_field_ndims
!=======================================================================

!=======================================================================
  FUNCTION inquire_dimension(Field_name)
!-----------------------------------------------------------------------
! NAME
!     inquire_dimension
!
! DESCRIPTION
!     Check if a given dimension is present in the input file.

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  USE netcdf, only: nf90_inq_dimid, NF90_NOERR
  USE mod_phys_lmdz_para, only: is_master, bcast

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
    CHARACTER(LEN=*), INTENT(IN) :: Field_name

! LOCAL VARIABLES
! ---------------
    LOGICAL :: inquire_dimension
    INTEGER :: varid
    INTEGER :: ierr

! CODE
! ----
    IF (is_master) THEN
      ierr=NF90_INQ_DIMID(nid_start,Field_name,varid)
      IF (ierr==NF90_NOERR) THEN
        Inquire_dimension=.TRUE.
      ELSE
        Inquire_dimension=.FALSE.
      ENDIF
    ENDIF

    CALL bcast(inquire_dimension)

  END FUNCTION inquire_dimension
!=======================================================================

!=======================================================================
  FUNCTION inquire_dimension_length(Field_name)
!-----------------------------------------------------------------------
! NAME
!     inquire_dimension_length
!
! DESCRIPTION
!     Give the length of the "Field_name" dimension stored in the input file.

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  USE netcdf, only: nf90_inquire_dimension, nf90_inq_dimid, &
                    NF90_NOERR, nf90_strerror
  USE mod_phys_lmdz_para, only: is_master, bcast

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
    CHARACTER(LEN=*), INTENT(IN) :: Field_name

! LOCAL VARIABLES
! ---------------
    INTEGER :: inquire_dimension_length
    INTEGER :: varid
    INTEGER :: ierr

! CODE
! ----
    IF (is_master) THEN
      ierr=nf90_inq_dimid(nid_start,Field_name,varid)
      ierr=nf90_inquire_dimension(nid_start,varid,&
                                  len=inquire_dimension_length)
      IF (ierr.NE.NF90_NOERR) THEN
        write(*,*)'inquire_field_length: problem obtaining length of '&
                  //trim(field_name)
        write(*,*)trim(nf90_strerror(ierr))
        CALL abort_physic("inquire_field_ndims","Failed to get length",1)
      ENDIF
    ENDIF

    CALL bcast(inquire_dimension_length)

  END FUNCTION inquire_dimension_length
!=======================================================================

!=======================================================================
  SUBROUTINE Get_Field_r1(field_name,field,found,timeindex)
!-----------------------------------------------------------------------
! NAME
!     Get_Field_r1
!
! DESCRIPTION
!     Read a surface field (1D) from the start file.

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  use mod_grid_phy_lmdz, only: klon_glo

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN)            :: field_name
  REAL,             INTENT(INOUT)         :: field(:)
  LOGICAL,          INTENT(OUT), OPTIONAL :: found
  INTEGER,          INTENT(IN),  OPTIONAL :: timeindex

! LOCAL VARIABLES
! ---------------
    integer :: edges(4), corners(4)

! CODE
! ----
    edges(1)=klon_glo
    edges(2:4)=1
    corners(1:4)=1
    if (PRESENT(timeindex)) then
      corners(2)=timeindex
    endif

    IF (PRESENT(found)) THEN
      CALL Get_field_rgen(field_name,field,1,corners,edges,found)
    ELSE
      CALL Get_field_rgen(field_name,field,1,corners,edges)
    ENDIF

  END SUBROUTINE Get_Field_r1
!=======================================================================

!=======================================================================
  SUBROUTINE Get_Field_r2(field_name,field,found,timeindex)
!-----------------------------------------------------------------------
! NAME
!     Get_Field_r2
!
! DESCRIPTION
!     Read a "3D" horizontal-vertical field (2D) from the start file.

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  use mod_grid_phy_lmdz, only: klon_glo

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN)            :: field_name
  REAL,             INTENT(INOUT)         :: field(:,:)
  LOGICAL,          INTENT(OUT), OPTIONAL :: found
  INTEGER,          INTENT(IN),  OPTIONAL :: timeindex

! LOCAL VARIABLES
! ---------------
    integer :: edges(4), corners(4)

! CODE
! ----
    edges(1)=klon_glo
    edges(2)=size(field,2)
    edges(3:4)=1
    corners(1:4)=1
    if (PRESENT(timeindex)) then
      corners(3)=timeindex
    endif

    IF (PRESENT(found)) THEN
      CALL Get_field_rgen(field_name,field,size(field,2),&
                          corners,edges,found)
    ELSE
      CALL Get_field_rgen(field_name,field,size(field,2),&
                          corners,edges)
    ENDIF


  END SUBROUTINE Get_Field_r2
!=======================================================================

!=======================================================================
  SUBROUTINE Get_Field_r3(field_name,field,found,timeindex)
!-----------------------------------------------------------------------
! NAME
!     Get_Field_r3
!
! DESCRIPTION
!     Read a "4D" field surf/alt/?? (3D) from the start file.

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  use mod_grid_phy_lmdz, only: klon_glo

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN)            :: field_name
  REAL,             INTENT(INOUT)         :: field(:,:,:)
  LOGICAL,          INTENT(OUT), OPTIONAL :: found
  INTEGER,          INTENT(IN),  OPTIONAL :: timeindex

! LOCAL VARIABLES
! ---------------
    integer :: edges(4), corners(4)

! CODE
! ----
    edges(1)=klon_glo
    edges(2)=size(field,2)
    edges(3)=size(field,3)
    edges(4)=1
    corners(1:4)=1
    if (PRESENT(timeindex)) then
      corners(4)=timeindex
    endif

    IF (PRESENT(found)) THEN
      CALL Get_field_rgen(field_name,field,size(field,2)*size(field,3),&
                          corners,edges,found)
    ELSE
      CALL Get_field_rgen(field_name,field,size(field,2)*size(field,3),&
                          corners,edges)
    ENDIF

  END SUBROUTINE Get_Field_r3
!=======================================================================

!=======================================================================
  SUBROUTINE Get_field_rgen(field_name,field,field_size, &
                            corners,edges,found)
!-----------------------------------------------------------------------
! NAME
!     Get_field_rgen
!
! DESCRIPTION
!     Generic subroutine to read fields from start file with scatter to processes.

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  USE netcdf
  USE dimphy
  USE mod_grid_phy_lmdz
  USE mod_phys_lmdz_para

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
    CHARACTER(LEN=*), INTENT(IN)            :: field_name
    INTEGER,          INTENT(IN)            :: field_size
    REAL,             INTENT(OUT)           :: field(klon,field_size)
    INTEGER,          INTENT(IN)            :: corners(4)
    INTEGER,          INTENT(IN)            :: edges(4)
    LOGICAL,          INTENT(OUT), OPTIONAL :: found

! LOCAL VARIABLES
! ---------------
    REAL    :: field_glo(klon_glo,field_size)
    LOGICAL :: tmp_found
    INTEGER :: varid
    INTEGER :: ierr

! CODE
! ----
    IF (is_master) THEN

      ierr=NF90_INQ_VARID(nid_start,field_name,varid)

      IF (ierr==NF90_NOERR) THEN
        CALL body(field_glo)
        tmp_found=.TRUE.
      ELSE
        tmp_found=.FALSE.
      ENDIF

    ENDIF

    CALL bcast(tmp_found)

    IF (tmp_found) THEN
      CALL scatter(field_glo,field)
    ENDIF

    IF (PRESENT(found)) THEN
      found=tmp_found
    ELSE
      IF (.NOT. tmp_found) THEN
        write(*,*) 'get_field_rgen: Field <'//field_name//'> not found'
        CALL abort_physic("get_field_rgen","Failed to get field",1)
      ENDIF
    ENDIF


        CONTAINS
    !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    !=======================================================================
      SUBROUTINE body(field_glo)
    !-----------------------------------------------------------------------
    ! NAME
    !     body
    !
    ! DESCRIPTION
    !     Read the requested field from the start file into the global array.
    !
    ! AUTHORS & DATE
    !     L. Lange
    !     JB Clement, 2023-2025
    !
    ! NOTES
    !     Helper contained in Get_field_rgen.
    !-----------------------------------------------------------------------

    ! DECLARATION
    ! -----------
      ! DECLARATION
! -----------
implicit none

    ! ARGUMENTS
    ! ---------
      REAL, INTENT(OUT) :: field_glo(klon_glo*field_size) ! Buffer for global field

    ! CODE
    ! ----
          ierr=NF90_GET_VAR(nid_start,varid,field_glo,corners,edges)
          IF (ierr/=NF90_NOERR) THEN
         ! La variable exist dans le fichier mais la lecture a echouee.
         write(*,*) 'get_field_rgen: Failed reading <'//field_name//'>'

!           IF (field_name=='CLWCON' .OR. field_name=='RNEBCON' .OR. field_name=='RATQS') THEN
!              ! Essaye de lire le variable sur surface uniqument, comme fait avant
!              field_glo(:)=0.
!              ierr=NF90_GET_VAR(nid_start,varid,field_glo(1:klon_glo))
!              IF (ierr/=NF90_NOERR) THEN
!                 write(*,*) 'phyetat0: Lecture echouee aussi en 2D pour <'//field_name//'>'
!                 CALL stop_clean
!              ELSE
!                 write(*,*) 'phyetat0: La variable <'//field_name//'> lu sur surface seulement'!, selon ancien format, le reste mis a zero'
!              END IF
!           ELSE
              CALL abort_physic("get_field_rgen","Failed to read field",1)
!           ENDIF
         ENDIF

     END SUBROUTINE body

  END SUBROUTINE Get_field_rgen
!=======================================================================

!=======================================================================
  SUBROUTINE get_var_r0(var_name,var,found)
!-----------------------------------------------------------------------
! NAME
!     get_var_r0
!
! DESCRIPTION
!     Get a scalar variable from input file.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN)            :: var_name
  REAL,             INTENT(INOUT)         :: var
  LOGICAL,          INTENT(OUT), OPTIONAL :: found

! LOCAL VARIABLES
! ---------------
    REAL :: varout(1)

! CODE
! ----
    IF (PRESENT(found)) THEN
      CALL Get_var_rgen(var_name,varout,size(varout),found)
    ELSE
      CALL Get_var_rgen(var_name,varout,size(varout))
    ENDIF
    var=varout(1)

  END SUBROUTINE get_var_r0
!=======================================================================

!=======================================================================
  SUBROUTINE get_var_r1(var_name,var,found)
!-----------------------------------------------------------------------
! NAME
!     get_var_r1
!
! DESCRIPTION
!     Get a vector (1D) variable from input file.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN)            :: var_name
  REAL,             INTENT(INOUT)         :: var(:)
  LOGICAL,          INTENT(OUT), OPTIONAL :: found

! CODE
! ----
    IF (PRESENT(found)) THEN
      CALL Get_var_rgen(var_name,var,size(var),found)
    ELSE
      CALL Get_var_rgen(var_name,var,size(var))
    ENDIF

  END SUBROUTINE get_var_r1
!=======================================================================

!=======================================================================
  SUBROUTINE get_var_r2(var_name,var,found)
!-----------------------------------------------------------------------
! NAME
!     get_var_r2
!
! DESCRIPTION
!     Get a 2D field variable from input file.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN)            :: var_name
  REAL,             INTENT(OUT)           :: var(:,:)
  LOGICAL,          INTENT(OUT), OPTIONAL :: found

! CODE
! ----
    IF (PRESENT(found)) THEN
      CALL Get_var_rgen(var_name,var,size(var),found)
    ELSE
      CALL Get_var_rgen(var_name,var,size(var))
    ENDIF

  END SUBROUTINE get_var_r2
!=======================================================================

!=======================================================================
  SUBROUTINE get_var_r3(var_name,var,found)
!-----------------------------------------------------------------------
! NAME
!     get_var_r3
!
! DESCRIPTION
!     Get a 3D field variable from input file.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN)            :: var_name
  REAL,             INTENT(INOUT)         :: var(:,:,:)
  LOGICAL,          INTENT(OUT), OPTIONAL :: found

! CODE
! ----
    IF (PRESENT(found)) THEN
      CALL Get_var_rgen(var_name,var,size(var),found)
    ELSE
      CALL Get_var_rgen(var_name,var,size(var))
    ENDIF

  END SUBROUTINE get_var_r3
!=======================================================================

!=======================================================================
  SUBROUTINE Get_var_rgen(var_name,var,var_size,found)
!-----------------------------------------------------------------------
! NAME
!     Get_var_rgen
!
! DESCRIPTION
!     Generic subroutine to read variables from input file.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  USE netcdf
  USE dimphy
  USE mod_grid_phy_lmdz
  USE mod_phys_lmdz_para

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
    CHARACTER(LEN=*), INTENT(IN)            :: var_name
    INTEGER,          INTENT(IN)            :: var_size
    REAL,             INTENT(OUT)           :: var(var_size)
    LOGICAL,          INTENT(OUT), OPTIONAL :: found

! LOCAL VARIABLES
! ---------------
    LOGICAL :: tmp_found
    INTEGER :: varid
    INTEGER :: ierr

! CODE
! ----
    IF (is_mpi_root .AND. is_omp_root) THEN

      ierr=NF90_INQ_VARID(nid_start,var_name,varid)

      IF (ierr==NF90_NOERR) THEN
        ierr=NF90_GET_VAR(nid_start,varid,var)
        IF (ierr/=NF90_NOERR) THEN
          write(*,*) 'phyetat0: Failed loading <'//trim(var_name)//'>'
          CALL abort_physic("get_var_rgen","Failed to read variable",1)
        ENDIF
        tmp_found=.TRUE.
      ELSE
        tmp_found=.FALSE.
      ENDIF

    ENDIF

    CALL bcast(tmp_found)

    IF (tmp_found) THEN
      CALL bcast(var)
    ENDIF

    IF (PRESENT(found)) THEN
      found=tmp_found
    ELSE
      IF (.NOT. tmp_found) THEN
        write(*,*) 'phyetat0: Variable <'//trim(var_name)//'> not found'
        CALL abort_physic("get_var_rgen","Failed to read variable",1)
      ENDIF
    ENDIF

  END SUBROUTINE Get_var_rgen
!=======================================================================

!=======================================================================
  SUBROUTINE open_restartphy(filename)
!-----------------------------------------------------------------------
! NAME
!     open_restartphy
!
! DESCRIPTION
!     Open (or create) the restart netCDF file for writing.
!     Define dimensions and global attributes on first call.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  USE netcdf, only: NF90_CREATE, NF90_CLOBBER, NF90_64BIT_OFFSET, &
                    NF90_NOERR, nf90_strerror, &
                    NF90_PUT_ATT, NF90_GLOBAL, NF90_DEF_DIM, &
                    NF90_UNLIMITED, NF90_ENDDEF, &
                    NF90_WRITE, NF90_OPEN
  USE mod_phys_lmdz_para, only: is_master
  USE mod_grid_phy_lmdz,  only: klon_glo
  USE dimphy,             only: klev, klevp1
  USE tracer_mod,         only: nqmx
  USE soil,               only:  nsoilmx_PEM
  USE comslope_mod,       only: nslope
  use layered_deposits,   only: nb_str_max

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN) :: filename

! LOCAL VARIABLES
! ---------------
    INTEGER           :: ierr
    LOGICAL,SAVE      :: already_created=.false.

! CODE
! ----
    IF (is_master) THEN
      if (.not.already_created) then
        ! At the very first call, create the file
        ierr=NF90_CREATE(filename,IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), &
                          nid_restart)
        IF (ierr/=NF90_NOERR) THEN
          write(*,*)'open_restartphy: problem creating file '//trim(filename)
          write(*,*)trim(nf90_strerror(ierr))
          CALL abort_physic("open_restartphy","Failed creating file",1)
        ENDIF
        already_created=.true.
      else
        ! Just open the file
        ierr=NF90_OPEN(filename,NF90_WRITE,nid_restart)
        IF (ierr/=NF90_NOERR) THEN
          write(*,*)'open_restartphy: problem opening file '//trim(filename)
          write(*,*)trim(nf90_strerror(ierr))
          CALL abort_physic("open_restartphy","Failed opening file",1)
        ENDIF
        return
      endif ! of if (.not.already_created)

      ierr=NF90_PUT_ATT(nid_restart,NF90_GLOBAL,"title",&
                        "Physics start file")
      IF (ierr/=NF90_NOERR) THEN
        write(*,*)'open_restartphy: problem writing title '
        write(*,*)trim(nf90_strerror(ierr))
      ENDIF

      ierr=NF90_DEF_DIM(nid_restart,"index",length,idim1)
      IF (ierr/=NF90_NOERR) THEN
        write(*,*)'open_restartphy: problem defining index dimension '
        write(*,*)trim(nf90_strerror(ierr))
        CALL abort_physic("open_restartphy","Failed defining index",1)
      ENDIF

      ierr=NF90_DEF_DIM(nid_restart,"physical_points",klon_glo,idim2)
      IF (ierr/=NF90_NOERR) THEN
        write(*,*)'open_restartphy: problem defining physical_points dimension '
        write(*,*)trim(nf90_strerror(ierr))
        CALL abort_physic("open_restartphy","Failed defining physical_points",1)
      ENDIF

      ierr=NF90_DEF_DIM(nid_restart,"subsurface_layers",nsoilmx_PEM,idim3)
      IF (ierr/=NF90_NOERR) THEN
        write(*,*)'open_restartphy: problem defining subsurface_layers dimension '
        write(*,*)trim(nf90_strerror(ierr))
        CALL abort_physic("open_restartphy","Failed defining subsurface_layers",1)
      ENDIF

      ierr=NF90_DEF_DIM(nid_restart,"nlayer_plus_1",klevp1,idim4)
      IF (ierr/=NF90_NOERR) THEN
        write(*,*)'open_restartphy: problem defining nlayer_plus_1 dimension '
        write(*,*)trim(nf90_strerror(ierr))
        CALL abort_physic("open_restartphy","Failed defining nlayer_plus_1",1)
      ENDIF

      if (nqmx.ne.0) then
        ierr=NF90_DEF_DIM(nid_restart,"number_of_advected_fields",nqmx,idim5)
      else
        ! pretend nqmx=1 because 0 size implies unlimited dimension for NetCDF
        ierr=NF90_DEF_DIM(nid_restart,"number_of_advected_fields",1,idim5)
      endif
      IF (ierr/=NF90_NOERR) THEN
        write(*,*)'open_restartphy: problem defining number_of_advected_fields dimension '
        write(*,*)trim(nf90_strerror(ierr))
        CALL abort_physic("open_restartphy","Failed defining number_of_advected_fields",1)
      ENDIF

      ierr=NF90_DEF_DIM(nid_restart,"nlayer",klev,idim6)
      IF (ierr/=NF90_NOERR) THEN
        write(*,*)'open_restartphy: problem defining nlayer dimension '
        write(*,*)trim(nf90_strerror(ierr))
        CALL abort_physic("open_restartphy","Failed defining nlayer",1)
      ENDIF

      ierr=NF90_DEF_DIM(nid_restart,"Time",NF90_UNLIMITED,idim7)
      IF (ierr/=NF90_NOERR) THEN
        write(*,*)'open_restartphy: problem defining Time dimension '
        write(*,*)trim(nf90_strerror(ierr))
        CALL abort_physic("open_restartphy","Failed defining Time",1)
      ENDIF

      ierr=NF90_DEF_DIM(nid_restart,"nslope",nslope,idim8)
      IF (ierr/=NF90_NOERR) THEN
        write(*,*)'phyredem: problem defining nslope dimension'
        write(*,*)trim(nf90_strerror(ierr))
        CALL ABORT
      ENDIF

      ierr=NF90_DEF_DIM(nid_restart,"inter slope",nslope+1,idim9)
      IF (ierr/=NF90_NOERR) THEN
        write(*,*)'phyredem: problem defining inter slope dimension'
        write(*,*)trim(nf90_strerror(ierr))
        CALL ABORT
      ENDIF

      ierr=NF90_DEF_DIM(nid_restart,"nb_str_max",nb_str_max,idim10)
      IF (ierr/=NF90_NOERR) THEN
        write(*,*)'phyredem: problem defining nb_str_max dimension'
        write(*,*)trim(nf90_strerror(ierr))
        CALL ABORT
      ENDIF

      ierr=NF90_ENDDEF(nid_restart)
      IF (ierr/=NF90_NOERR) THEN
        write(*,*)'open_restartphy: problem ending definition mode '
        write(*,*)trim(nf90_strerror(ierr))
        CALL abort_physic("open_restartphy","Failed ending definition mode",1)
      ENDIF
    ENDIF

  END SUBROUTINE open_restartphy
!=======================================================================

!=======================================================================
  SUBROUTINE close_restartphy
!-----------------------------------------------------------------------
! NAME
!     close_restartphy
!
! DESCRIPTION
!     Close the restart netCDF file.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  USE netcdf,             only: NF90_CLOSE
  USE mod_phys_lmdz_para, only: is_master

! DECLARATION
! -----------
implicit none

! LOCAL VARIABLES
! ---------------
    INTEGER :: ierr

! CODE
! ----
    IF (is_master) THEN
      ierr = NF90_CLOSE (nid_restart)
    ENDIF

  END SUBROUTINE close_restartphy
!=======================================================================

!=======================================================================
  SUBROUTINE put_field_r1(field_name,title,field,time)
!-----------------------------------------------------------------------
! NAME
!     put_field_r1
!
! DESCRIPTION
!     Write a surface field to the restart file.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none


! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN)           :: field_name
  CHARACTER(LEN=*), INTENT(IN)           :: title
  REAL,             INTENT(IN)           :: field(:)
  REAL,             INTENT(IN), OPTIONAL :: time

! CODE
! ----
  IF (present(time)) THEN
    ! if timeindex is present, it is a time-dependent variable
    CALL put_field_rgen(field_name,title,field,1,time)
  ELSE
    CALL put_field_rgen(field_name,title,field,1)
  ENDIF

  END SUBROUTINE put_field_r1
!=======================================================================

!=======================================================================
  SUBROUTINE put_field_r2(field_name,title,field,time)
!-----------------------------------------------------------------------
! NAME
!     put_field_r2
!
! DESCRIPTION
!     Write a "3D" horizontal-vertical field (2D) to the restart file.

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN)           :: field_name
  CHARACTER(LEN=*), INTENT(IN)           :: title
  REAL,             INTENT(IN)           :: field(:,:)
  REAL,             INTENT(IN), OPTIONAL :: time

! CODE
! ----
  IF (present(time)) THEN
    ! if timeindex is present, it is a time-dependent variable
    CALL put_field_rgen(field_name,title,field,size(field,2),time)
  ELSE
    CALL put_field_rgen(field_name,title,field,size(field,2))
  ENDIF

  END SUBROUTINE put_field_r2
!=======================================================================

!=======================================================================
  SUBROUTINE put_field_r3(field_name,title,field,time)
!-----------------------------------------------------------------------
! NAME
!     put_field_r3
!
! DESCRIPTION
!     Write a "4D" field surf/alt/?? (3D) to the restart file.

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN)           :: field_name
  CHARACTER(LEN=*), INTENT(IN)           :: title
  REAL,             INTENT(IN)           :: field(:,:,:)
  REAL,             INTENT(IN), OPTIONAL :: time

! CODE
! ----
  IF (present(time)) THEN
    ! if timeindex is present, it is a time-dependent variable
    CALL put_field_rgen(field_name,title,field,size(field,2)*size(field,3),&
                        time)
  ELSE
    CALL put_field_rgen(field_name,title,field,size(field,2)*size(field,3))
  ENDIF

  END SUBROUTINE put_field_r3
!=======================================================================

!=======================================================================
  SUBROUTINE put_field_rgen(field_name,title,field,field_size,time)
!-----------------------------------------------------------------------
! NAME
!     put_field_rgen
!
! DESCRIPTION
!     Generic subroutine to write fields to restart file with gather from processes

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  USE netcdf
  USE dimphy
  USE soil,         only: nsoilmx_PEM
  USE comslope_mod, ONLY: nslope
  USE mod_grid_phy_lmdz
  USE mod_phys_lmdz_para

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*),INTENT(IN)    :: field_name
  CHARACTER(LEN=*),INTENT(IN)    :: title
  INTEGER,INTENT(IN)             :: field_size
  REAL,INTENT(IN)                :: field(klon,field_size)
  REAL,OPTIONAL,INTENT(IN)       :: time

! LOCAL VARIABLES
! ---------------
  REAL                           :: field_glo(klon_glo,field_size)
  INTEGER                        :: ierr
  INTEGER                        :: nvarid

! CODE
! ----
    CALL gather(field,field_glo)

    IF (is_master) THEN

     IF (index(field_name,"stratif") == 0) then
      IF (field_size==1) THEN
        ! input is a 1D "surface field" array
        if (.not.present(time)) then ! for a time-independent field
          ierr=NF90_REDEF(nid_restart)
#ifdef NC_DOUBLE
          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                            (/idim2/),nvarid)
#else
          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                            (/idim2/),nvarid)
#endif
          if (ierr.ne.NF90_NOERR) then
            write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
            write(*,*)trim(nf90_strerror(ierr))
          endif
          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
          ierr=NF90_ENDDEF(nid_restart)
          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo)
        else
          ! check if the variable has already been defined:
          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
          if (ierr/=NF90_NOERR) then ! variable not found, define it
            ierr=NF90_REDEF(nid_restart)
#ifdef NC_DOUBLE
            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                              (/idim2,idim7/),nvarid)
#else
            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                              (/idim2,idim7/),nvarid)
#endif
            if (ierr.ne.NF90_NOERR) then
              write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
              write(*,*)trim(nf90_strerror(ierr))
            endif
            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
            ierr=NF90_ENDDEF(nid_restart)
          endif
          ! Write the variable
          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
                            start=(/1,timeindex/))
        endif ! of if (.not.present(timeindex))

      ELSE IF (field_size==klev) THEN
        ! input is a 2D "atmospheric field" array
        if (.not.present(time)) then ! for a time-independent field
          ierr=NF90_REDEF(nid_restart)
#ifdef NC_DOUBLE
          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                            (/idim2,idim6/),nvarid)
#else
          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                            (/idim2,idim6/),nvarid)
#endif
          if (ierr.ne.NF90_NOERR) then
            write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
            write(*,*)trim(nf90_strerror(ierr))
          endif
          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
          ierr=NF90_ENDDEF(nid_restart)
          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo)
        else
          ! check if the variable has already been defined:
          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
          if (ierr/=NF90_NOERR) then ! variable not found, define it
            ierr=NF90_REDEF(nid_restart)
#ifdef NC_DOUBLE
            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                              (/idim2,idim6,idim7/),nvarid)
#else
            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                              (/idim2,idim6,idim7/),nvarid)
#endif
            if (ierr.ne.NF90_NOERR) then
              write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
              write(*,*)trim(nf90_strerror(ierr))
            endif
            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
            ierr=NF90_ENDDEF(nid_restart)
          endif
          ! Write the variable
          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
                            start=(/1,1,timeindex/))
        endif ! of if (.not.present(time))

      ELSE IF (field_size==klevp1) THEN
        ! input is a 2D "interlayer atmospheric field" array
        if (.not.present(time)) then ! for a time-independent field
          ierr=NF90_REDEF(nid_restart)
#ifdef NC_DOUBLE
          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                            (/idim2,idim4/),nvarid)
#else
          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                            (/idim2,idim4/),nvarid)
#endif
          if (ierr.ne.NF90_NOERR) then
            write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
            write(*,*)trim(nf90_strerror(ierr))
          endif
          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
          ierr = NF90_ENDDEF(nid_restart)
          ierr = NF90_PUT_VAR(nid_restart,nvarid,field_glo)
        else
          ! check if the variable has already been defined:
          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
          if (ierr/=NF90_NOERR) then ! variable not found, define it
            ierr=NF90_REDEF(nid_restart)
#ifdef NC_DOUBLE
            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                              (/idim2,idim4,idim7/),nvarid)
#else
            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                              (/idim2,idim4,idim7/),nvarid)
#endif
            if (ierr.ne.NF90_NOERR) then
              write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
              write(*,*)trim(nf90_strerror(ierr))
            endif
            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
            ierr=NF90_ENDDEF(nid_restart)
          endif
          ! Write the variable
          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
                            start=(/1,1,timeindex/))
        endif ! of if (.not.present(timeindex))

      ELSE IF (field_size==nsoilmx_PEM) THEN
        ! input is a 2D "subsurface field" array
        if (.not.present(time)) then ! for a time-independent field
          ierr = NF90_REDEF(nid_restart)
#ifdef NC_DOUBLE
          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                            (/idim2,idim3/),nvarid)
#else
          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                            (/idim2,idim3/),nvarid)
#endif
          if (ierr.ne.NF90_NOERR) then
            write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
            write(*,*)trim(nf90_strerror(ierr))
          endif
          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
          ierr = NF90_ENDDEF(nid_restart)
          ierr = NF90_PUT_VAR(nid_restart,nvarid,field_glo)
        else
          ! check if the variable has already been defined:
          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
          if (ierr/=NF90_NOERR) then ! variable not found, define it
            ierr=NF90_REDEF(nid_restart)
#ifdef NC_DOUBLE
            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                              (/idim2,idim3,idim7/),nvarid)
#else
            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                              (/idim2,idim3,idim7/),nvarid)
#endif
           if (ierr.ne.NF90_NOERR) then
              write(*,*)"put_field_rgen error: failed to define "//trim(field_name)
              write(*,*)trim(nf90_strerror(ierr))
            endif
            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
            ierr=NF90_ENDDEF(nid_restart)
          endif
          ! Write the variable
          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
                            start=(/1,1,timeindex/))

        endif ! of if (.not.present(time))


      ELSE IF (field_size==nslope) THEN
        ! input is a 2D "subsurface field" array
        if (.not.present(time)) then ! for a time-independent field
          ierr = NF90_REDEF(nid_restart)
#ifdef NC_DOUBLE
          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                            (/idim2,idim8/),nvarid)
#else
          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                            (/idim2,idim8/),nvarid)
#endif
          if (ierr.ne.NF90_NOERR) then
            write(*,*)"put_field_rgen error: failed to define"//trim(field_name)
            write(*,*)trim(nf90_strerror(ierr))
          endif
          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
          ierr = NF90_ENDDEF(nid_restart)
          ierr = NF90_PUT_VAR(nid_restart,nvarid,field_glo)
        else
          ! check if the variable has already been defined:
          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
          if (ierr/=NF90_NOERR) then ! variable not found, define it
            ierr=NF90_REDEF(nid_restart)
#ifdef NC_DOUBLE
            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                              (/idim2,idim8,idim7/),nvarid)
#else
            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                              (/idim2,idim8,idim7/),nvarid)
#endif
           if (ierr.ne.NF90_NOERR) then
              write(*,*)"put_field_rgen error: failed to define"//trim(field_name)
              write(*,*)trim(nf90_strerror(ierr))
            endif
            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
            ierr=NF90_ENDDEF(nid_restart)
          endif
          ! Write the variable
          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
                            start=(/1,1,timeindex/))

        endif ! of if (.not.present(time))

      ELSE
        write(*,*) "Error phyredem(put_field_rgen) : wrong dimension for ",trim(field_name)
        write(*,*) "  field_size =",field_size
        CALL abort_physic("put_field_rgen","wrong field dimension",1)
      ENDIF
     ELSE
        ! input is a 2D "stratification" array
        if (.not.present(time)) then ! for a time-independent field
          ierr = NF90_REDEF(nid_restart)
#ifdef NC_DOUBLE
          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                            (/idim2,idim8/),nvarid)
#else
          ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                            (/idim2,idim8/),nvarid)
#endif
          if (ierr.ne.NF90_NOERR) then
            write(*,*)"put_field_rgen error: failed to define"//trim(field_name)
            write(*,*)trim(nf90_strerror(ierr))
          endif
          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
          ierr = NF90_ENDDEF(nid_restart)
          ierr = NF90_PUT_VAR(nid_restart,nvarid,field_glo)
        else
          ! check if the variable has already been defined:
          ierr=NF90_INQ_VARID(nid_restart,field_name,nvarid)
          if (ierr/=NF90_NOERR) then ! variable not found, define it
            ierr=NF90_REDEF(nid_restart)
#ifdef NC_DOUBLE
            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                              (/idim2,idim10,idim7/),nvarid)
#else
            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                              (/idim2,idim10,idim7/),nvarid)
#endif
           if (ierr.ne.NF90_NOERR) then
              write(*,*)"put_field_rgen error: failed to define"//trim(field_name)
              write(*,*)trim(nf90_strerror(ierr))
            endif
            IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
            ierr=NF90_ENDDEF(nid_restart)
          endif
          ! Write the variable
          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
                            start=(/1,1,timeindex/))

        endif ! of if (.not.present(time))
     ENDIF

      ! Check the writting of field to file went OK
      if (ierr.ne.NF90_NOERR) then
        write(*,*) " Error phyredem(put_field_rgen) : failed writing ",trim(field_name)
        write(*,*)trim(nf90_strerror(ierr))
        CALL abort_physic("put_field_rgen","Failed writing field",1)
      endif

    ENDIF ! of IF (is_master)

  END SUBROUTINE put_field_rgen
!=======================================================================

!=======================================================================
  SUBROUTINE put_var_r0(var_name,title,var)
!-----------------------------------------------------------------------
! NAME
!     put_var_r0
!
! DESCRIPTION
!     Write a scalar variable to the restart file

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN) :: var_name
  CHARACTER(LEN=*), INTENT(IN) :: title
  REAL,             INTENT(IN) :: var

! LOCAL VARIABLES
! ---------------
     REAL :: varin(1)

! CODE
! ----
     varin(1)=var

     CALL put_var_rgen(var_name,title,varin,size(varin))

  END SUBROUTINE put_var_r0
!=======================================================================

!=======================================================================
  SUBROUTINE put_var_r1(var_name,title,var)
!-----------------------------------------------------------------------
! NAME
!    put_var_r1
!
! DESCRIPTION
!     Write a vector (1D) variable to the restart file

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN) :: var_name
  CHARACTER(LEN=*), INTENT(IN) :: title
  REAL,             INTENT(IN) :: var(:)

! CODE
! ----
  CALL put_var_rgen(var_name,title,var,size(var))

  END SUBROUTINE put_var_r1
!=======================================================================

!=======================================================================
  SUBROUTINE put_var_r2(var_name,title,var)
!-----------------------------------------------------------------------
! NAME
!     put_var_r2
!
! DESCRIPTION
!     Write a 3D field variable to the restart file

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN) :: var_name
  CHARACTER(LEN=*), INTENT(IN) :: title
  REAL,             INTENT(IN) :: var(:,:)

! CODE
! ----
  CALL put_var_rgen(var_name,title,var,size(var))

  END SUBROUTINE put_var_r2
!=======================================================================

!=======================================================================
  SUBROUTINE put_var_r3(var_name,title,var)
!-----------------------------------------------------------------------
! NAME
!     put_var_r3
!
! DESCRIPTION
!     Write a 3D field variable to the restart file
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN) :: var_name
  CHARACTER(LEN=*), INTENT(IN) :: title
  REAL,             INTENT(IN) :: var(:,:,:)

! CODE
! ----
  CALL put_var_rgen(var_name,title,var,size(var))

  END SUBROUTINE put_var_r3
!=======================================================================

!=======================================================================
  SUBROUTINE put_var_rgen(var_name,title,var,var_size)
!-----------------------------------------------------------------------
! NAME
!     put_var_rgen
!
! DESCRIPTION
!     Generic subroutine to write variables to the restart file

! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025

! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
  USE netcdf, only: NF90_REDEF, NF90_DEF_VAR, NF90_ENDDEF, NF90_PUT_VAR, &
                    NF90_FLOAT, NF90_DOUBLE, &
                    NF90_PUT_ATT, NF90_NOERR, nf90_strerror, &
                    nf90_inq_dimid, nf90_inquire_dimension, NF90_INQ_VARID
  USE soil,               only:  nsoilmx_PEM
  USE comslope_mod,       only: nslope
  USE mod_phys_lmdz_para, only: is_master

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
  CHARACTER(LEN=*), INTENT(IN) :: var_name
  CHARACTER(LEN=*), INTENT(IN) :: title
  INTEGER,          INTENT(IN) :: var_size
  REAL,             INTENT(IN) :: var(var_size)

! LOCAL VARIABLES
! ---------------
  INTEGER       :: ierr
  INTEGER       :: nvarid
  INTEGER       :: idim1d
  LOGICAL, SAVE :: firsttime = .true.

! CODE
! ----
    IF (is_master) THEN

      IF (var_name=="Time") THEN
        ! Very specific case of "Time" variable
        if (firsttime) then
          ! Create the "Time variable"
          ierr=NF90_REDEF(nid_restart)
#ifdef NC_DOUBLE
          ierr=NF90_DEF_VAR(nid_restart,var_name,NF90_DOUBLE,&
                            (/idim7/),nvarid)
#else
          ierr=NF90_DEF_VAR(nid_restart,var_name,NF90_FLOAT,&
                            (/idim7/),nvarid)
#endif
          IF (LEN_TRIM(title) > 0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
          ierr=NF90_ENDDEF(nid_restart)

          firsttime=.false.
        endif
        ! Append "time" value
        ! get current length of "Time" dimension
        ierr=nf90_inq_dimid(nid_restart,var_name,idim1d)
        ierr=nf90_inquire_dimension(nid_restart,idim1d,len=timeindex)
        timeindex=timeindex+1
        ierr=NF90_INQ_VARID(nid_restart,var_name,nvarid)
        ierr=NF90_PUT_VAR(nid_restart,nvarid,var,&
                           start=(/timeindex/))
        IF (ierr/=NF90_NOERR) THEN
          write(*,*)'put_var_rgen: problem writing Time'
          write(*,*)trim(nf90_strerror(ierr))
          CALL abort_physic("get_var_rgen","Failed to write Time",1)
        ENDIF
        return ! nothing left to do
      ELSEIF (var_size==length) THEN
        ! We know it is a "controle" kind of 1D array
        idim1d=idim1
      ELSEIF (var_size==nsoilmx_PEM) THEN
        ! We know it is an "mlayer" kind of 1D array
        idim1d=idim3
      ELSEIF (var_size==nslope+1) THEN
        ! We know it is an "inter slope" kind of 1D array
        idim1d=idim9
      ELSEIF (var_name == "nb_str_max") THEN
        ! We know it is a kind of stratification array
        idim1d = idim10
      ELSE
        write(*,*) "put_var_rgen error : wrong dimension"
        write(*,*) "  var_size =",var_size
        CALL abort_physic("get_var_rgen","Wrong variable dimension",1)

      ENDIF ! of IF (var_size==length) THEN

      ! Swich to NetCDF define mode
      ierr=NF90_REDEF (nid_restart)
      ! Define the variable
#ifdef NC_DOUBLE
      ierr=NF90_DEF_VAR(nid_restart,var_name,NF90_DOUBLE,(/idim1d/),nvarid)
#else
      ierr=NF90_DEF_VAR(nid_restart,var_name,NF90_FLOAT,(/idim1d/),nvarid)
#endif
      ! Add a "title" attribute
      IF (LEN_TRIM(title)>0) ierr=NF90_PUT_ATT(nid_restart,nvarid,"title",title)
      ! Swich out of define mode
      ierr=NF90_ENDDEF(nid_restart)
      ! Write variable to file
      ierr=NF90_PUT_VAR(nid_restart,nvarid,var)
      IF (ierr/=NF90_NOERR) THEN
        write(*,*)'put_var_rgen: problem writing '//trim(var_name)
        write(*,*)trim(nf90_strerror(ierr))
        CALL abort_physic("get_var_rgen","Failed writing variable",1)
      ENDIF
    ENDIF ! of IF (is_master)

  END SUBROUTINE put_var_rgen
!=======================================================================

END MODULE iostart_pem