! $Id$

MODULE wxios
    USE lmdz_xios

    !Variables disponibles pendant toute l'execution du programme:
    
    INTEGER, SAVE :: g_comm
    CHARACTER(len=100), SAVE :: g_ctx_name ="LMDZ"
    TYPE(xios_context), SAVE :: g_ctx
!$OMP THREADPRIVATE(g_comm,g_ctx_name,g_ctx)
    LOGICAL, SAVE :: g_flag_xml = .FALSE.
    CHARACTER(len=100) :: g_field_name = "nofield"
!$OMP THREADPRIVATE(g_flag_xml,g_field_name)
    REAL :: missing_val_omp
    REAL :: missing_val
!$OMP THREADPRIVATE(missing_val)

#ifdef XIOS1
#error "XIOS v1 no longer supported, use XIOS v2."
#endif

    CONTAINS
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !   36day => 36d etc     !!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    SUBROUTINE reformadate(odate, ndate)
        CHARACTER(len=*), INTENT(IN) :: odate
        TYPE(xios_duration) :: ndate
        
        INTEGER :: i = 0
         !!!!!!!!!!!!!!!!!!
         ! Pour XIOS:
         !  year : y
         !  month : mo
         !  day : d
         !  hour : h
         !  minute : mi
         !  second : s
         !!!!!!!!!!!!!!!!!!

        i = INDEX(odate, "day")
        IF (i > 0) THEN
            read(odate(1:i-1),*) ndate%day
        END IF

        i = INDEX(odate, "hr")
        IF (i > 0) THEN
            read(odate(1:i-1),*) ndate%hour
        END IF

        i = INDEX(odate, "mth")
        IF (i > 0) THEN
            read(odate(1:i-1),*) ndate%month
        END IF
        
        !IF (prt_level >= 10) WRITE(lunout,*) "Xios. ", odate, " => ", ndate
    END SUBROUTINE reformadate
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !   ave(X) => average etc     !!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    CHARACTER(len=7) FUNCTION reformaop(op)
        CHARACTER(len=*), INTENT(IN) :: op
        
        INTEGER :: i = 0
        reformaop = "average"
        
        IF (op=="inst(X)") THEN
            reformaop = "instant"
        END IF
        
        IF (op=="once") THEN
            reformaop = "once"
        END IF
        
        IF (op=="t_max(X)") THEN
            reformaop = "maximum"
        END IF
        
        IF (op=="t_min(X)") THEN
            reformaop = "minimum"
        END IF
        
        !IF (prt_level >= 10) WRITE(lunout,*) "Xios. ", op, " => ", reformaop
    END FUNCTION reformaop

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Routine d'initialisation      !!!!!!!!!!!!!
    !     A lancer juste après mpi_init !!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    SUBROUTINE wxios_init(xios_ctx_name, locom, outcom, type_ocean)
        USE print_control_mod, ONLY: prt_level, lunout
        IMPLICIT NONE

      CHARACTER(len=*), INTENT(IN) :: xios_ctx_name
      INTEGER, INTENT(IN), OPTIONAL :: locom
      INTEGER, INTENT(OUT), OPTIONAL :: outcom
      CHARACTER(len=6), INTENT(IN), OPTIONAL :: type_ocean

    
        TYPE(xios_context) :: xios_ctx
        INTEGER :: xios_comm

        IF (prt_level >= 10) WRITE(lunout,*) "wxios_init: Initialization"



        IF (PRESENT(locom)) THEN
          CALL xios_initialize(xios_ctx_name, local_comm = locom, return_comm = xios_comm )
          IF (prt_level >= 10) WRITE(lunout,*) "wxios_init: ctx=",xios_ctx_name," local_comm=",locom,", return_comm=",xios_comm
        ELSE
          CALL xios_initialize(xios_ctx_name, return_comm = xios_comm )
          IF (prt_level >= 10) WRITE(lunout,*) "wxios_init: ctx=",xios_ctx_name," return_comm=",xios_comm
        END IF
        
        IF (PRESENT(outcom)) THEN
          outcom = xios_comm
          IF (prt_level >= 10) WRITE(lunout,*) "wxios_init: ctx=",xios_ctx_name," outcom=",outcom
        END IF
        
        !Enregistrement des variables globales:
        g_comm = xios_comm
        g_ctx_name = xios_ctx_name
        
!        ! Si couple alors init fait dans cpl_init
!        IF (.not. PRESENT(type_ocean)) THEN
!            CALL wxios_context_init()
!        ENDIF
         WRITE(*,*)'END of WXIOS_INIT', g_comm , g_ctx_name 

    END SUBROUTINE wxios_init

    SUBROUTINE wxios_context_init()
        USE print_control_mod, ONLY: prt_level, lunout
        USE mod_phys_lmdz_mpi_data, ONLY: COMM_LMDZ_PHY
        IMPLICIT NONE

        TYPE(xios_context) :: xios_ctx

!$OMP MASTER
        !Initialisation du contexte:
        !!CALL xios_context_initialize(g_ctx_name, g_comm)
        CALL xios_context_initialize(g_ctx_name, COMM_LMDZ_PHY)
        CALL xios_get_handle(g_ctx_name, xios_ctx)    !Récupération
        CALL xios_set_current_context(xios_ctx)            !Activation
        g_ctx = xios_ctx

        CALL wxios_add_group_init
        
        IF (prt_level >= 10) THEN
          WRITE(lunout,*) "wxios_context_init: Current context is ",trim(g_ctx_name)
          WRITE(lunout,*) "     now CALL xios_solve_inheritance()"
        ENDIF
        !Une première analyse des héritages:
        CALL xios_solve_inheritance()
!$OMP END MASTER
    END SUBROUTINE wxios_context_init


    
    SUBROUTINE wxios_add_group_init

      ! routine create by Anne Cozic (2023)
      ! This routine will create field associated to group defined without description of fields in field.xml file
      ! This routine need to be CALL before "xios_sole_inheritance" after an !$OMP MASTER directive
      
      USE infotrac_phy, ONLY: nbtr, nqtot, nqo, type_trac, tracers, niso, ntiso
      USE strings_mod,  ONLY: maxlen

      IMPLICIT NONE

      TYPE(xios_fieldgroup) :: group_handle, philev_hdl
      TYPE(xios_field) :: child
      INTEGER :: k, iq
      CHARACTER(len=12) :: nvar, name_phi    
      CHARACTER(LEN=maxlen) :: varname, dn
      CHARACTER(LEN=maxlen) :: unt
    

      ! group create for StratAER variables 
      !On ajoute les variables 3D traceurs par l interface fortran 
      CALL xios_get_handle("fields_strataer_trac_3D", group_handle)
      ! On boucle sur les traceurs pour les ajouter au groupe puis fixer les attributs
      DO iq = 1, nqtot
         IF(.NOT.(tracers(iq)%isAdvected .AND. tracers(iq)%isInPhysics)) CYCLE
         dn = 'd'//TRIM(tracers(iq)%name)//'_'

         unt = "kg kg-1"
         varname=trim(tracers(iq)%name)
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
         
         unt = "kg kg-1 s-1"
         varname=TRIM(dn)//'vdf'
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
         varname=TRIM(dn)//'the'
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
         varname=TRIM(dn)//'con'
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
         varname=TRIM(dn)//'lessi_impa'
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
         varname=TRIM(dn)//'lessi_nucl'
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
         varname=TRIM(dn)//'insc'
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
         varname=TRIM(dn)//'bcscav'
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
         varname=TRIM(dn)//'evapls'
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
         varname=TRIM(dn)//'ls'
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
         varname=TRIM(dn)//'trsp'
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
         varname=TRIM(dn)//'sscav'
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
         varname=TRIM(dn)//'sat'
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
         varname=TRIM(dn)//'uscav'
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
      END DO
      !On ajoute les variables 2D traceurs par l interface fortran 
      CALL xios_get_handle("fields_strataer_trac_2D", group_handle)
      ! On boucle sur les traceurs pour les ajouter au groupe puis fixer les attributs
      DO iq = 1, nqtot
         IF(.NOT.(tracers(iq)%isAdvected .AND. tracers(iq)%isInPhysics)) CYCLE
         
         unt = "kg m-2"
         varname='cum'//trim(tracers(iq)%name)
         
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
         
         unt = "kg m-2 s-1"
         varname='cumd'//trim(tracers(iq)%name)//'_dry'
         CALL xios_add_child(group_handle, child, varname)
         CALL xios_set_attr(child, name=varname, unit=unt)
      ENDDO
      

      ! group create for offline mass flow variables
      CALL xios_get_handle("philev_grp", philev_hdl)

      DO k=1,79
         IF (k<10) THEN
            WRITE(nvar,'(i1)') k
         ELSE IF (k<100) THEN
            WRITE(nvar,'(i2)') k
         ELSE
            WRITE(nvar,'(i3)') k
         END IF
         name_phi= "phi_lev"//TRIM(nvar)
         CALL xios_add_child(philev_hdl, child, "phi_lev"//TRIM(nvar)) 
         CALL xios_set_attr(child, name=trim(name_phi))
      ENDDO


    END SUBROUTINE wxios_add_group_init

    SUBROUTINE wxios_set_context()
        IMPLICIT NONE
        TYPE(xios_context) :: xios_ctx

       !$OMP MASTER
        CALL xios_get_handle(g_ctx_name, xios_ctx)    !Récupération
        CALL xios_set_current_context(xios_ctx)            !Activation
       !$OMP END MASTER

    END SUBROUTINE wxios_set_context

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Routine de paramétrisation !!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    SUBROUTINE wxios_set_cal(pasdetemps, calendrier, annee, mois, jour, heure, ini_an, ini_mois, ini_jour, ini_heure)
        USE print_control_mod, ONLY: prt_level, lunout
        IMPLICIT NONE

     !Paramètres:
     CHARACTER(len=*), INTENT(IN) :: calendrier
     INTEGER, INTENT(IN) :: annee, mois, jour, ini_an, ini_mois, ini_jour
     REAL, INTENT(IN) :: pasdetemps, heure, ini_heure
     
     !Variables:
     CHARACTER(len=80) :: abort_message
     CHARACTER(len=19) :: date
     INTEGER :: njour = 1
     
     !Variables pour xios:
     TYPE(xios_duration) :: mdtime
     !REAL(kind = 8) :: year=0, month=0, day=0, hour=0, minute=0, second=0
     
        mdtime%second=pasdetemps

        !Réglage du calendrier:
        SELECT CASE (calendrier)
            CASE('earth_360d')
                CALL xios_define_calendar("D360")
                IF (prt_level >= 10) WRITE(lunout,*) 'wxios_set_cal: Calendrier terrestre a 360 jours/an'
            CASE('earth_365d')
                CALL xios_define_calendar("NoLeap")
                IF (prt_level >= 10) WRITE(lunout,*) 'wxios_set_cal: Calendrier terrestre a 365 jours/an'
            CASE('gregorian')
                CALL xios_define_calendar("Gregorian")
                IF (prt_level >= 10) WRITE(lunout,*) 'wxios_set_cal: Calendrier gregorien'
            CASE DEFAULT
                abort_message = 'wxios_set_cal: Mauvais choix de calendrier'
                CALL abort_physic('Gcm:Xios',abort_message,1)
        END SELECT
        
        !Formatage de la date d'origine:
        WRITE(date, "(i4.4,'-',i2.2,'-',i2.2,' ',i2.2,':00:00')") annee, mois, jour, int(heure) 
        
        IF (prt_level >= 10) WRITE(lunout,*) "wxios_set_cal: Time origin: ", date
        CALL xios_set_time_origin(xios_date(annee,mois,jour,int(heure),0,0))

        !Formatage de la date de debut:

        WRITE(date, "(i4.4,'-',i2.2,'-',i2.2,' ',i2.2,':00:00')") ini_an, ini_mois, ini_jour, int(ini_heure)
        
        IF (prt_level >= 10) WRITE(lunout,*) "wxios_set_cal: Start date: ", date
        
        CALL xios_set_start_date(xios_date(ini_an,ini_mois,ini_jour,int(ini_heure),0,0))
        
        !Et enfin,le pas de temps:
        CALL xios_set_timestep(mdtime)
        IF (prt_level >= 10) WRITE(lunout,*) "wxios_set_cal: ts=",mdtime
    END SUBROUTINE wxios_set_cal

    SUBROUTINE wxios_set_timestep(ts)
        REAL, INTENT(IN) :: ts
        TYPE(xios_duration) :: mdtime     

        mdtime%timestep = ts

        CALL xios_set_timestep(mdtime)
    END SUBROUTINE wxios_set_timestep

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Pour initialiser un domaine !!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    SUBROUTINE wxios_domain_param(dom_id)
       USE dimphy, only: klon
       USE mod_phys_lmdz_transfert_para, ONLY: gather, bcast
       USE mod_phys_lmdz_para, only: jj_nb, jj_begin, jj_end, ii_begin, ii_end, &
                                     mpi_size, mpi_rank, klon_mpi, &
                                     is_sequential, is_south_pole_dyn
       USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, klon_glo         
       USE print_control_mod, ONLY: prt_level, lunout
       USE geometry_mod

       IMPLICIT NONE
        CHARACTER(len=*),INTENT(IN) :: dom_id ! domain identifier

        REAL   :: rlat_glo(klon_glo)
        REAL   :: rlon_glo(klon_glo)
        REAL   :: io_lat(nbp_lat)
        REAL   :: io_lon(nbp_lon)
        LOGICAL :: mask(nbp_lon,jj_nb) !Masque pour les problèmes de recouvrement MPI
        TYPE(xios_domain) :: dom
        INTEGER :: i
        LOGICAL :: boool
        


        CALL gather(latitude_deg,rlat_glo) 
        CALL bcast(rlat_glo)
        CALL gather(longitude_deg,rlon_glo)
        CALL bcast(rlon_glo)
    
  !$OMP MASTER  
        io_lat(1)=rlat_glo(1)
        io_lat(nbp_lat)=rlat_glo(klon_glo)
        IF ((nbp_lon*nbp_lat) > 1) then
          DO i=2,nbp_lat-1
            io_lat(i)=rlat_glo(2+(i-2)*nbp_lon)
          ENDDO
        ENDIF

        IF (klon_glo == 1) THEN
          io_lon(1)=rlon_glo(1)
        ELSE 
          io_lon(1:nbp_lon)=rlon_glo(2:nbp_lon+1)
        ENDIF

        
        !On récupère le handle:
        CALL xios_get_handle(dom_id, dom)
        
        !On parametrise le domaine:
        CALL xios_set_attr(dom, ni_glo=nbp_lon, ibegin=0, ni=nbp_lon, type="rectilinear")
        CALL xios_set_attr(dom, nj_glo=nbp_lat, jbegin=jj_begin-1, nj=jj_nb, data_dim=2)
        CALL xios_set_attr(dom, lonvalue_1d=io_lon(1:nbp_lon), latvalue_1d=io_lat(jj_begin:jj_end))
        CALL xios_set_domain_attr("dom_out", domain_ref=dom_id)

        !On definit un axe de latitudes pour les moyennes zonales
        IF (xios_is_valid_axis("axis_lat")) THEN
           CALL xios_set_axis_attr( "axis_lat", n_glo=nbp_lat, n=jj_nb, begin=jj_begin-1, value=io_lat(jj_begin:jj_end))
        ENDIF
        IF (xios_is_valid_axis("axis_lat_greordered")) THEN
           CALL xios_set_axis_attr( "axis_lat_greordered", n_glo=nbp_lat, n=jj_nb, begin=jj_begin-1, &
                                    value=io_lat(jj_begin:jj_end)*(-1.))
        ENDIF

        IF (.NOT.is_sequential) THEN
            mask(:,:)=.TRUE.
            if (ii_begin>1) mask(1:ii_begin-1,1) = .FALSE.
            if (ii_end<nbp_lon) mask(ii_end+1:nbp_lon,jj_nb) = .FALSE.
            ! special case for south pole
            if ((ii_end==1).and.(is_south_pole_dyn)) mask(1:nbp_lon,jj_nb)=.TRUE.
            IF (prt_level >= 10) THEN
              WRITE(lunout,*) "wxios_domain_param: mpirank=",mpi_rank," mask(:,1)=",mask(:,1)
              WRITE(lunout,*) "wxios_domain_param: mpirank=",mpi_rank," mask(:,jj_nb)=",mask(:,jj_nb)
            ENDIF
            CALL xios_set_attr(dom, mask_2d=mask)
        END IF

         CALL xios_is_defined_attr(dom,ni_glo=boool)
        !Vérification:
        IF (xios_is_valid_domain(dom_id)) THEN
            IF (prt_level >= 10) WRITE(lunout,*) "wxios_domain_param: Domain initialized: ", trim(dom_id), boool
        ELSE
            IF (prt_level >= 10) WRITE(lunout,*) "wxios_domain_param: Invalid domain: ", trim(dom_id)
        END IF
!$OMP END MASTER
        
    END SUBROUTINE wxios_domain_param
    

    SUBROUTINE wxios_domain_param_unstructured(dom_id)
        USE geometry_mod, ONLY: longitude, latitude, boundslon, boundslat,ind_cell_glo
        USE mod_grid_phy_lmdz, ONLY: nvertex, klon_glo
        USE mod_phys_lmdz_para 
        USE nrtype, ONLY: PI
        USE ioipsl_getin_p_mod, ONLY: getin_p
        IMPLICIT NONE
        CHARACTER(len=*),INTENT(IN) :: dom_id ! domain identifier
        REAL :: lon_mpi(klon_mpi)
        REAL :: lat_mpi(klon_mpi)
        REAL :: boundslon_mpi(klon_mpi,nvertex)
        REAL :: boundslat_mpi(klon_mpi,nvertex)
        INTEGER :: ind_cell_glo_mpi(klon_mpi)
        TYPE(xios_domain) :: dom
        
        LOGICAL :: remap_output

        CALL gather_omp(longitude*180/PI,lon_mpi)
        CALL gather_omp(latitude*180/PI,lat_mpi)
        CALL gather_omp(boundslon*180/PI,boundslon_mpi)
        CALL gather_omp(boundslat*180/PI,boundslat_mpi)
        CALL gather_omp(ind_cell_glo,ind_cell_glo_mpi)
        
        remap_output=.TRUE.
        CALL getin_p("remap_output",remap_output)

!$OMP MASTER
        CALL xios_get_handle(dom_id, dom)
        
        !On parametrise le domaine:
        CALL xios_set_attr(dom, ni_glo=klon_glo, ibegin=ij_begin-1, ni=ij_nb, type="unstructured")
        CALL xios_set_attr(dom, nvertex=nvertex, lonvalue_1d=lon_mpi, latvalue_1d=lat_mpi, &
                           bounds_lon_1d=TRANSPOSE(boundslon_mpi), bounds_lat_1d=TRANSPOSE(boundslat_mpi) )
        CALL xios_set_attr(dom, i_index=ind_cell_glo_mpi(:)-1)
        IF (remap_output) THEN
          CALL xios_set_domain_attr("dom_out", domain_ref="dom_regular")
          CALL xios_set_fieldgroup_attr("remap_expr", expr="@this_ref")
          CALL xios_set_fieldgroup_attr("remap_1800s", freq_op=xios_duration_convert_from_string("1800s"))
          CALL xios_set_fieldgroup_attr("remap_1h",    freq_op=xios_duration_convert_from_string("1h"))
          CALL xios_set_fieldgroup_attr("remap_3h",    freq_op=xios_duration_convert_from_string("3h"))
          CALL xios_set_fieldgroup_attr("remap_6h",    freq_op=xios_duration_convert_from_string("6h"))
          CALL xios_set_fieldgroup_attr("remap_1d",    freq_op=xios_duration_convert_from_string("1d"))
          CALL xios_set_fieldgroup_attr("remap_1mo",   freq_op=xios_duration_convert_from_string("1mo"))
        ENDIF
!$OMP END MASTER

    END SUBROUTINE wxios_domain_param_unstructured




    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Pour déclarer un axe vertical !!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    SUBROUTINE wxios_add_vaxis(axis_id, axis_size, axis_value,     &
                               positif, bnds)
        USE print_control_mod, ONLY: prt_level, lunout
        IMPLICIT NONE

        CHARACTER (len=*), INTENT(IN) :: axis_id
        INTEGER, INTENT(IN) :: axis_size
        REAL, DIMENSION(axis_size), INTENT(IN) :: axis_value
        CHARACTER (len=*), INTENT(IN), OPTIONAL :: positif
        REAL, DIMENSION(axis_size, 2), INTENT(IN), OPTIONAL :: bnds
        
!        TYPE(xios_axisgroup) :: axgroup
!        TYPE(xios_axis) :: ax
!        CHARACTER(len=50) :: axis_id 
        
!        IF (len_trim(axisgroup_id).gt.len(axis_id)) THEN
!          WRITE(lunout,*) "wxios_add_vaxis: error, size of axis_id too small!!"
!          WRITE(lunout,*) "     increase it to at least ",len_trim(axisgroup_id)
!          CALL abort_gcm("wxios_add_vaxis","len(axis_id) too small",1)
!        ENDIF
!        axis_id=trim(axisgroup_id)
        
        !On récupère le groupe d'axes qui va bien:
        !CALL xios_get_axisgroup_handle(axisgroup_id, axgroup)
        
        !On ajoute l'axe correspondant à ce fichier:
        !CALL xios_add_axis(axgroup, ax, TRIM(ADJUSTL(axis_id)))
        
        !Et on le parametrise:
        !CALL xios_set_axis_attr_hdl(ax, size=axis_size, value=axis_value)
        
        ! Ehouarn: New way to declare axis, without axis_group:
        if (PRESENT(positif) .AND. PRESENT(bnds)) then
          CALL xios_set_axis_attr(trim(axis_id),n_glo=axis_size,value=axis_value, &
                                  positive=positif, bounds=bnds)
        else if (PRESENT(positif)) then
          CALL xios_set_axis_attr(trim(axis_id),n_glo=axis_size,value=axis_value, &
                                  positive=positif)
        else if (PRESENT(bnds)) then
          CALL xios_set_axis_attr(trim(axis_id),n_glo=axis_size,value=axis_value, &
                                  bounds=bnds)
        else
          CALL xios_set_axis_attr(trim(axis_id),n_glo=axis_size,value=axis_value)
        endif

        !Vérification:
        IF (xios_is_valid_axis(TRIM(ADJUSTL(axis_id)))) THEN
            IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_vaxis: Axis created: ", TRIM(ADJUSTL(axis_id))
        ELSE
            WRITE(lunout,*) "wxios_add_vaxis: Invalid axis: ", TRIM(ADJUSTL(axis_id))
        END IF

    END SUBROUTINE wxios_add_vaxis
    
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Pour déclarer un fichier  !!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    SUBROUTINE wxios_add_file(fname, ffreq, flvl)
        USE print_control_mod, ONLY: prt_level, lunout
        IMPLICIT NONE

        CHARACTER(len=*), INTENT(IN) :: fname
        CHARACTER(len=*), INTENT(IN) :: ffreq
        INTEGER, INTENT(IN) :: flvl
        
        TYPE(xios_file) :: x_file
        TYPE(xios_filegroup) :: x_fg
        TYPE(xios_duration) :: nffreq
        
        !On regarde si le fichier n'est pas défini par XML:
        IF (.NOT.xios_is_valid_file(fname)) THEN
            !On créé le noeud:
            CALL xios_get_handle("defile", x_fg)
            CALL xios_add_child(x_fg, x_file, fname)
        
            !On reformate la fréquence:
            CALL reformadate(ffreq, nffreq)
        
            !On configure:
            CALL xios_set_attr(x_file, name="X"//fname,&
                output_freq=nffreq, output_level=flvl, enabled=.TRUE.)
       
            IF (xios_is_valid_file("X"//fname)) THEN
                IF (prt_level >= 10) THEN
                  WRITE(lunout,*) "wxios_add_file: New file: ", "X"//fname
                  WRITE(lunout,*) "wxios_add_file: output_freq=",nffreq,"; output_lvl=",flvl
                ENDIF
            ELSE
                WRITE(lunout,*) "wxios_add_file: Error, invalid file: ", "X"//trim(fname)
                WRITE(lunout,*) "wxios_add_file: output_freq=",nffreq,"; output_lvl=",flvl
            END IF
        ELSE
            IF (prt_level >= 10) THEN
              WRITE(lunout,*) "wxios_add_file: File ",trim(fname), " défined using XML."
            ENDIF
            ! Ehouarn: add an enable=.TRUE. on top of xml definitions... why???
            CALL xios_set_file_attr(fname, enabled=.TRUE.)
        END IF
    END SUBROUTINE wxios_add_file
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Pour créer un champ      !!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    SUBROUTINE wxios_add_field(fieldname, fieldgroup, fieldlongname, fieldunit)
        USE netcdf, only: nf90_fill_real

        IMPLICIT NONE
        INCLUDE 'iniprint.h'
        
        CHARACTER(len=*), INTENT(IN) :: fieldname
        TYPE(xios_fieldgroup), INTENT(IN) :: fieldgroup
        CHARACTER(len=*), INTENT(IN) :: fieldlongname
        CHARACTER(len=*), INTENT(IN) :: fieldunit
        
        TYPE(xios_field) :: field
        CHARACTER(len=10) :: newunit
        REAL(KIND=8) :: def
        
        !La valeur par défaut des champs non définis:
        def = nf90_fill_real
        
        IF (fieldunit == " ") THEN
            newunit = "-"
        ELSE
            newunit = fieldunit
        ENDIF
        
        !On ajoute le champ:
        CALL xios_add_child(fieldgroup, field, fieldname)
        !IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_field: ",fieldname,fieldgroup, fieldlongname, fieldunit
        
        !On rentre ses paramètres:
        CALL xios_set_attr(field, standard_name=fieldlongname, unit=newunit, default_value=def)
        IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_field: Field ",trim(fieldname), "cree:"
        IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_field: long_name=",trim(fieldlongname),"; unit=",trim(newunit),";  default_value=",nf90_fill_real

    END SUBROUTINE wxios_add_field
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Pour déclarer un champ      !!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    SUBROUTINE wxios_add_field_to_file(fieldname, fdim, fid, fname, fieldlongname, fieldunit, field_level, op, nam_axvert)
        USE print_control_mod, ONLY: prt_level, lunout
        IMPLICIT NONE

        CHARACTER(len=*), INTENT(IN) :: fieldname
        INTEGER, INTENT(IN)          :: fdim, fid
        CHARACTER(len=*), INTENT(IN) :: fname
        CHARACTER(len=*), INTENT(IN) :: fieldlongname
        CHARACTER(len=*), INTENT(IN) :: fieldunit
        INTEGER, INTENT(IN)          :: field_level
        CHARACTER(len=*), INTENT(IN) :: op
        
        CHARACTER(len=20) :: axis_id ! Ehouarn: dangerous...
        CHARACTER(len=20), INTENT(IN), OPTIONAL :: nam_axvert
        CHARACTER(len=100) :: operation
        TYPE(xios_file) :: f
        TYPE(xios_field) :: field
        TYPE(xios_fieldgroup) :: fieldgroup
        TYPE(xios_duration) :: freq_op

        LOGICAL :: bool=.FALSE.
        INTEGER :: lvl =0
        
        
        ! Ajout Abd pour NMC:
        IF (fid<=6) THEN
          axis_id="presnivs"
        ELSE
          axis_id="plev"
        ENDIF
 
        IF (PRESENT(nam_axvert)) THEN
           axis_id=nam_axvert
           PRINT*,'nam_axvert=',axis_id
        ENDIF
        
        !on prépare le nom de l'opération:
        operation = reformaop(op)
        
        
        !On selectionne le bon groupe de champs:
        IF (fdim==2) THEN
          CALL xios_get_handle("fields_2D", fieldgroup)
        ELSE
          CALL xios_get_handle("fields_3D", fieldgroup)
        ENDIF
        
        !On regarde si le champ à déjà été créé ou non:
        IF (xios_is_valid_field(fieldname) .AND. .NOT. g_field_name == fieldname) THEN
            !Si ce champ existe via XML (ie, dès le premier passage, ie g_field_name != fieldname) alors rien d'autre à faire
            IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_field_to_file: Field ", trim(fieldname), "exists via XML"
            g_flag_xml = .TRUE.
            g_field_name = fieldname

        ELSE IF (.NOT. g_field_name == fieldname) THEN
            !Si premier pssage et champ indéfini, alors on le créé

            IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_field_to_file: Field ", trim(fieldname), "does not exist"
            
            !On le créé:
            CALL wxios_add_field(fieldname,  fieldgroup, fieldlongname, fieldunit)
            IF (xios_is_valid_field(fieldname)) THEN
                IF (prt_level >= 10) WRITE(lunout,*) "wxios_add_field_to_file: Field ", trim(fieldname), "created"
            ENDIF

            g_flag_xml = .FALSE.
            g_field_name = fieldname

        END IF

        IF (.NOT. g_flag_xml) THEN
            !Champ existe déjà, mais pas XML, alors on l'ajoute
            !On ajoute le champ:
            CALL xios_get_handle(fname, f)
            CALL xios_add_child(f, field)
            
            
            !L'operation, sa frequence:
            freq_op%timestep=1
            CALL xios_set_attr(field, field_ref=fieldname, operation=TRIM(ADJUSTL(operation)), freq_op=freq_op, prec=4)

            
            !On rentre ses paramètres:
            CALL xios_set_attr(field, level=field_level, enabled=.TRUE.)
            
            IF (fdim==2) THEN
                !Si c'est un champ 2D:
                IF (prt_level >= 10) THEN
                  WRITE(lunout,*) "wxios_add_field_to_file: 2D Field ", trim(fieldname), " in ", "X"//trim(fname) ," configured with:"
                  WRITE(lunout,*) "wxios_add_field_to_file: op=", TRIM(ADJUSTL(operation))
                  WRITE(lunout,*) "wxios_add_field_to_file: freq_op=1ts","; lvl=",field_level
                ENDIF
            ELSE
                !Si 3D :
                !On ajoute l'axe vertical qui va bien:
                CALL xios_set_attr(field, axis_ref=TRIM(ADJUSTL(axis_id)))
                
                IF (prt_level >= 10) THEN
                  WRITE(lunout,*) "wxios_add_field_to_file: 3D Field",trim(fieldname), " in ", "X"//trim(fname), "configured with:"
                  WRITE(lunout,*) "wxios_add_field_to_file: freq_op=1ts","; lvl=",field_level
                  WRITE(lunout,*) "wxios_add_field_to_file: axis=",TRIM(ADJUSTL(axis_id))
                ENDIF
            END IF
        
        ELSE
            !Sinon on se contente de l'activer:
            CALL xios_set_field_attr(fieldname, enabled=.TRUE.)
            !NB: This will override an enable=.FALSE. set by a user in the xml file;
            !   then the only way to not output the field is by changing its
            !   output level
        ENDIF       
        
    END SUBROUTINE wxios_add_field_to_file
    
!    SUBROUTINE wxios_update_calendar(ito)
!        INTEGER, INTENT(IN) :: ito
!        CALL xios_update_calendar(ito)
!    END SUBROUTINE wxios_update_calendar

!    SUBROUTINE wxios_write_2D(fieldname, fdata)
!        CHARACTER(len=*), INTENT(IN) :: fieldname
!        REAL, DIMENSION(:,:), INTENT(IN) :: fdata

!        CALL xios_send_field(fieldname, fdata)
!    END SUBROUTINE wxios_write_2D
    
!    SUBROUTINE wxios_write_3D(fieldname, fdata)
!        CHARACTER(len=*), INTENT(IN) :: fieldname
!        REAL, DIMENSION(:,:,:), INTENT(IN) :: fdata

!        CALL xios_send_field(fieldname, fdata)
!    END SUBROUTINE wxios_write_3D
    
    SUBROUTINE wxios_closedef()
        CALL xios_close_context_definition()
!        CALL xios_update_calendar(0)
    END SUBROUTINE wxios_closedef
    
    SUBROUTINE wxios_close()
        CALL xios_context_finalize()
         CALL xios_finalize()
     END SUBROUTINE wxios_close
END MODULE wxios

