MODULE guide_loc_mod

  !=======================================================================
  !   Auteur:  F.Hourdin
  !            F. Codron 01/09
  !=======================================================================

  USE getparam, ONLY: ini_getparam, fin_getparam, getpar
  USE Write_Field_loc
  USE netcdf, ONLY: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close, &
          nf90_inq_dimid, nf90_inquire_dimension, nf90_inq_dimid, &
          nf90_inquire_dimension, nf90_enddef, nf90_def_dim, nf90_put_var, nf90_noerr, nf90_close, nf90_inq_varid, &
          nf90_redef, nf90_write, nf90_unlimited, nf90_float, nf90_clobber, nf90_64bit_offset, nf90_float, &
          nf90_create, nf90_def_var, nf90_open
  USE parallel_lmdz
  USE lmdz_pres2lev, ONLY: pres2lev


  IMPLICIT NONE

  ! ---------------------------------------------
  ! Declarations des cles logiques et parametres
  ! ---------------------------------------------
  INTEGER, PRIVATE, SAVE :: iguide_read, iguide_int, iguide_sav
  INTEGER, PRIVATE, SAVE :: nlevnc, guide_plevs
  LOGICAL, PRIVATE, SAVE :: guide_u, guide_v, guide_T, guide_Q, guide_P
  LOGICAL, PRIVATE, SAVE :: guide_hr, guide_teta
  LOGICAL, PRIVATE, SAVE :: guide_BL, guide_reg, guide_add, gamma4, guide_zon
  LOGICAL, PRIVATE, SAVE :: invert_p, invert_y, ini_anal
  LOGICAL, PRIVATE, SAVE :: guide_2D, guide_sav, guide_modele
  !FC
  LOGICAL, PRIVATE, SAVE :: convert_Pa

  REAL, PRIVATE, SAVE :: tau_min_u, tau_max_u
  REAL, PRIVATE, SAVE :: tau_min_v, tau_max_v
  REAL, PRIVATE, SAVE :: tau_min_T, tau_max_T
  REAL, PRIVATE, SAVE :: tau_min_Q, tau_max_Q
  REAL, PRIVATE, SAVE :: tau_min_P, tau_max_P

  REAL, PRIVATE, SAVE :: lat_min_g, lat_max_g
  REAL, PRIVATE, SAVE :: lon_min_g, lon_max_g
  REAL, PRIVATE, SAVE :: tau_lon, tau_lat

  REAL, PRIVATE, SAVE :: plim_guide_BL

  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_u, alpha_v
  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_T, alpha_Q
  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_P, alpha_pcor

  ! ---------------------------------------------
  ! Variables de guidage
  ! ---------------------------------------------
  ! Variables des fichiers de guidage
  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: unat1, unat2
  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: vnat1, vnat2
  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: tnat1, tnat2
  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: qnat1, qnat2
  REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE, SAVE :: pnat1, pnat2
  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: psnat1, psnat2
  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: apnc, bpnc
  ! Variables aux dimensions du modele
  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: ugui1, ugui2
  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: vgui1, vgui2
  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: tgui1, tgui2
  REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE, SAVE :: qgui1, qgui2
  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: psgui1, psgui2

  INTEGER, SAVE, PRIVATE :: ijbu, ijbv, ijeu, ijev !,ijnu,ijnv
  INTEGER, SAVE, PRIVATE :: jjbu, jjbv, jjeu, jjev, jjnu, jjnv


CONTAINS
  !=======================================================================

  SUBROUTINE guide_init

    USE control_mod, ONLY: day_step
    USE serre_mod, ONLY: grossismx

    IMPLICIT NONE

    INCLUDE "dimensions.h"
    INCLUDE "paramet.h"

    INTEGER :: error, ncidpl, rid, rcod
    CHARACTER (len = 80) :: abort_message
    CHARACTER (len = 20) :: modname = 'guide_init'
    CHARACTER (len = 20) :: namedim

    ! ---------------------------------------------
    ! Lecture des parametres:
    ! ---------------------------------------------
    CALL ini_getparam("nudging_parameters_out.txt")
    ! Variables guidees
    CALL getpar('guide_u', .TRUE., guide_u, 'guidage de u')
    CALL getpar('guide_v', .TRUE., guide_v, 'guidage de v')
    CALL getpar('guide_T', .TRUE., guide_T, 'guidage de T')
    CALL getpar('guide_P', .TRUE., guide_P, 'guidage de P')
    CALL getpar('guide_Q', .TRUE., guide_Q, 'guidage de Q')
    CALL getpar('guide_hr', .TRUE., guide_hr, 'guidage de Q par H.R')
    CALL getpar('guide_teta', .FALSE., guide_teta, 'guidage de T par Teta')

    CALL getpar('guide_add', .FALSE., guide_add, 'foréage constant?')
    CALL getpar('guide_zon', .FALSE., guide_zon, 'guidage moy zonale')
    IF (guide_zon .AND. abs(grossismx - 1.) > 0.01) &
            CALL abort_gcm("guide_init", &
                    "zonal nudging requires grid regular in longitude", 1)

    !   Constantes de rappel. Unite : fraction de jour
    CALL getpar('tau_min_u', 0.02, tau_min_u, 'Cste de rappel min, u')
    CALL getpar('tau_max_u', 10., tau_max_u, 'Cste de rappel max, u')
    CALL getpar('tau_min_v', 0.02, tau_min_v, 'Cste de rappel min, v')
    CALL getpar('tau_max_v', 10., tau_max_v, 'Cste de rappel max, v')
    CALL getpar('tau_min_T', 0.02, tau_min_T, 'Cste de rappel min, T')
    CALL getpar('tau_max_T', 10., tau_max_T, 'Cste de rappel max, T')
    CALL getpar('tau_min_Q', 0.02, tau_min_Q, 'Cste de rappel min, Q')
    CALL getpar('tau_max_Q', 10., tau_max_Q, 'Cste de rappel max, Q')
    CALL getpar('tau_min_P', 0.02, tau_min_P, 'Cste de rappel min, P')
    CALL getpar('tau_max_P', 10., tau_max_P, 'Cste de rappel max, P')
    CALL getpar('gamma4', .FALSE., gamma4, 'Zone sans rappel elargie')
    CALL getpar('guide_BL', .TRUE., guide_BL, 'guidage dans C.Lim')
    CALL getpar('plim_guide_BL', 85000., plim_guide_BL, 'BL top presnivs value')

    ! Sauvegarde du forçage
    CALL getpar('guide_sav', .FALSE., guide_sav, 'sauvegarde guidage')
    CALL getpar('iguide_sav', 4, iguide_sav, 'freq. sauvegarde guidage')
    ! frequences f>0: fx/jour; f<0: tous les f jours; f=0: 1 seule fois.
    IF (iguide_sav>0) THEN
      iguide_sav = day_step / iguide_sav
    ELSE if (iguide_sav == 0) THEN
      iguide_sav = huge(0)
    ELSE
      iguide_sav = day_step * iguide_sav
    ENDIF

    ! Guidage regional seulement (sinon constant ou suivant le zoom)
    CALL getpar('guide_reg', .FALSE., guide_reg, 'guidage regional')
    CALL getpar('lat_min_g', -90., lat_min_g, 'Latitude mini guidage ')
    CALL getpar('lat_max_g', 90., lat_max_g, 'Latitude maxi guidage ')
    CALL getpar('lon_min_g', -180., lon_min_g, 'longitude mini guidage ')
    CALL getpar('lon_max_g', 180., lon_max_g, 'longitude maxi guidage ')
    CALL getpar('tau_lat', 5., tau_lat, 'raideur lat guide regional ')
    CALL getpar('tau_lon', 5., tau_lon, 'raideur lon guide regional ')

    ! Parametres pour lecture des fichiers
    CALL getpar('iguide_read', 4, iguide_read, 'freq. lecture guidage')
    CALL getpar('iguide_int', 4, iguide_int, 'freq. interpolation vert')
    IF (iguide_int==0) THEN
      iguide_int = 1
    ELSEIF (iguide_int>0) THEN
      iguide_int = day_step / iguide_int
    ELSE
      iguide_int = day_step * iguide_int
    ENDIF
    CALL getpar('guide_plevs', 0, guide_plevs, 'niveaux pression fichiers guidage')
    ! Pour compatibilite avec ancienne version avec guide_modele
    CALL getpar('guide_modele', .FALSE., guide_modele, 'niveaux pression ap+bp*psol')
    IF (guide_modele) THEN
      guide_plevs = 1
    ENDIF
    !FC
    CALL getpar('convert_Pa', .TRUE., convert_Pa, 'Convert Pressure levels in Pa')
    ! Fin raccord
    CALL getpar('ini_anal', .FALSE., ini_anal, 'Etat initial = analyse')
    CALL getpar('guide_invertp', .TRUE., invert_p, 'niveaux p inverses')
    CALL getpar('guide_inverty', .TRUE., invert_y, 'inversion N-S')
    CALL getpar('guide_2D', .FALSE., guide_2D, 'fichier guidage lat-P')

    CALL fin_getparam

    ! ---------------------------------------------
    ! Determination du nombre de niveaux verticaux
    ! des fichiers guidage
    ! ---------------------------------------------
    ncidpl = -99
    IF (guide_plevs==1) THEN
      IF (ncidpl==-99) THEN
        rcod = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
        IF (rcod/=nf90_noerr) THEN
          abort_message = ' Nudging error -> no file apbp.nc'
          CALL abort_gcm(modname, abort_message, 1)
        endif
      endif
    elseif (guide_plevs==2) THEN
      IF (ncidpl==-99) THEN
        rcod = nf90_open('P.nc', nf90_nowrite, ncidpl)
        IF (rcod/=nf90_noerr) THEN
          abort_message = ' Nudging error -> no file P.nc'
          CALL abort_gcm(modname, abort_message, 1)
        endif
      endif

    elseif (guide_u) THEN
      IF (ncidpl==-99) THEN
        rcod = nf90_open('u.nc', nf90_nowrite, ncidpl)
        IF (rcod/=nf90_noerr) THEN
          abort_message = ' Nudging error -> no file u.nc'
          CALL abort_gcm(modname, abort_message, 1)
        endif

      endif

    elseif (guide_v) THEN
      IF (ncidpl==-99) THEN
        rcod = nf90_open('v.nc', nf90_nowrite, ncidpl)
        IF (rcod/=nf90_noerr) THEN
          abort_message = ' Nudging error -> no file v.nc'
          CALL abort_gcm(modname, abort_message, 1)
        endif
      endif

    elseif (guide_T) THEN
      IF (ncidpl==-99) THEN
        rcod = nf90_open('T.nc', nf90_nowrite, ncidpl)
        IF (rcod/=nf90_noerr) THEN
          abort_message = ' Nudging error -> no file T.nc'
          CALL abort_gcm(modname, abort_message, 1)
        endif
      endif

    elseif (guide_Q) THEN
      IF (ncidpl==-99) THEN
        rcod = nf90_open('hur.nc', nf90_nowrite, ncidpl)
        IF (rcod/=nf90_noerr) THEN
          abort_message = ' Nudging error -> no file hur.nc'
          CALL abort_gcm(modname, abort_message, 1)
        endif
      endif

    endif
    error = nf90_inq_dimid(ncidpl, 'LEVEL', rid)
    IF (error/=nf90_noerr) error = nf90_inq_dimid(ncidpl, 'PRESSURE', rid)
    IF (error/=nf90_noerr) THEN
      abort_message = 'Nudging: error reading pressure levels'
      CALL abort_gcm(modname, abort_message, 1)
    ENDIF
    error = nf90_inquire_dimension(ncidpl, rid, len = nlevnc)
    WRITE(*, *)trim(modname) // ' : number of vertical levels nlevnc', nlevnc
    rcod = nf90_close(ncidpl)

    ! ---------------------------------------------
    ! Allocation des variables
    ! ---------------------------------------------
    abort_message = 'nudging allocation error'

    ALLOCATE(apnc(nlevnc), stat = error)
    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    ALLOCATE(bpnc(nlevnc), stat = error)
    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    apnc = 0.;bpnc = 0.

    ALLOCATE(alpha_pcor(llm), stat = error)
    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    ALLOCATE(alpha_u(ijb_u:ije_u), stat = error)
    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    ALLOCATE(alpha_v(ijb_v:ije_v), stat = error)
    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    ALLOCATE(alpha_T(ijb_u:ije_u), stat = error)
    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    ALLOCATE(alpha_Q(ijb_u:ije_u), stat = error)
    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    ALLOCATE(alpha_P(ijb_u:ije_u), stat = error)
    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    alpha_u = 0.;alpha_v = 0;alpha_T = 0;alpha_Q = 0;alpha_P = 0

    IF (guide_u) THEN
      ALLOCATE(unat1(iip1, jjb_u:jje_u, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(ugui1(ijb_u:ije_u, llm), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(unat2(iip1, jjb_u:jje_u, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(ugui2(ijb_u:ije_u, llm), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      unat1 = 0.;unat2 = 0.;ugui1 = 0.;ugui2 = 0.
    ENDIF

    IF (guide_T) THEN
      ALLOCATE(tnat1(iip1, jjb_u:jje_u, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(tgui1(ijb_u:ije_u, llm), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(tnat2(iip1, jjb_u:jje_u, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(tgui2(ijb_u:ije_u, llm), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      tnat1 = 0.;tnat2 = 0.;tgui1 = 0.;tgui2 = 0.
    ENDIF

    IF (guide_Q) THEN
      ALLOCATE(qnat1(iip1, jjb_u:jje_u, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(qgui1(ijb_u:ije_u, llm), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(qnat2(iip1, jjb_u:jje_u, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(qgui2(ijb_u:ije_u, llm), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      qnat1 = 0.;qnat2 = 0.;qgui1 = 0.;qgui2 = 0.
    ENDIF

    IF (guide_v) THEN
      ALLOCATE(vnat1(iip1, jjb_v:jje_v, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(vgui1(ijb_v:ije_v, llm), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(vnat2(iip1, jjb_v:jje_v, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(vgui2(ijb_v:ije_v, llm), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      vnat1 = 0.;vnat2 = 0.;vgui1 = 0.;vgui2 = 0.
    ENDIF

    IF (guide_plevs==2) THEN
      ALLOCATE(pnat1(iip1, jjb_u:jje_u, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(pnat2(iip1, jjb_u:jje_u, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      pnat1 = 0.;pnat2 = 0.;
    ENDIF

    IF (guide_P.OR.guide_plevs==1) THEN
      ALLOCATE(psnat1(iip1, jjb_u:jje_u), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(psnat2(iip1, jjb_u:jje_u), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      psnat1 = 0.;psnat2 = 0.;
    ENDIF
    IF (guide_P) THEN
      ALLOCATE(psgui2(ijb_u:ije_u), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(psgui1(ijb_u:ije_u), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      psgui1 = 0.;psgui2 = 0.
    ENDIF

    ! ---------------------------------------------
    !   Lecture du premier etat de guidage.
    ! ---------------------------------------------
    IF (guide_2D) THEN
      CALL guide_read2D(1)
    ELSE
      CALL guide_read(1)
    ENDIF
    IF (guide_v) vnat1 = vnat2
    IF (guide_u) unat1 = unat2
    IF (guide_T) tnat1 = tnat2
    IF (guide_Q) qnat1 = qnat2
    IF (guide_plevs==2) pnat1 = pnat2
    IF (guide_P.OR.guide_plevs==1) psnat1 = psnat2

  END SUBROUTINE guide_init

  !=======================================================================
  SUBROUTINE guide_main(itau, ucov, vcov, teta, q, masse, ps)
    USE exner_hyb_loc_m, ONLY: exner_hyb_loc
    USE exner_milieu_loc_m, ONLY: exner_milieu_loc
    USE parallel_lmdz
    USE control_mod
    USE write_field_loc
    USE comconst_mod, ONLY: cpp, daysec, dtvr, kappa
    USE comvert_mod, ONLY: ap, bp, preff, presnivs, pressure_exner

    IMPLICIT NONE

    INCLUDE "dimensions.h"
    INCLUDE "paramet.h"

    ! Variables entree
    INTEGER, INTENT(IN) :: itau !pas de temps
    REAL, DIMENSION (ijb_u:ije_u, llm), INTENT(INOUT) :: ucov, teta, q, masse
    REAL, DIMENSION (ijb_v:ije_v, llm), INTENT(INOUT) :: vcov
    REAL, DIMENSION (ijb_u:ije_u), INTENT(INOUT) :: ps

    ! Variables locales
    LOGICAL, SAVE :: first = .TRUE.
    !$OMP THREADPRIVATE(first)
    LOGICAL :: f_out ! sortie guidage
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :) :: f_addu ! var aux: champ de guidage
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :) :: f_addv ! var aux: champ de guidage
    ! Variables pour fonction Exner (P milieu couche)
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :, :) :: pk
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :) :: pks
    REAL :: unskap
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :) :: p ! besoin si guide_P
    ! Compteurs temps:
    INTEGER, SAVE :: step_rea, count_no_rea, itau_test ! lecture guidage
    !$OMP THREADPRIVATE(step_rea,count_no_rea,itau_test)
    REAL :: ditau, dday_step
    REAL :: tau, reste ! position entre 2 etats de guidage
    REAL, SAVE :: factt ! pas de temps en fraction de jour
    !$OMP THREADPRIVATE(factt)

    INTEGER :: i, j, l
    CHARACTER(LEN = 20) :: modname = "guide_main"

    !$OMP MASTER
    ijbu = ij_begin ; ijeu = ij_end
    jjbu = jj_begin ; jjeu = jj_end ; jjnu = jjeu - jjbu + 1
    ijbv = ij_begin ; ijev = ij_end
    jjbv = jj_begin ; jjev = jj_end ; jjnv = jjev - jjbv + 1
    IF (pole_sud) THEN
      ijeu = ij_end - iip1
      ijev = ij_end - iip1
      jjev = jj_end - 1
      jjnv = jjev - jjbv + 1
    ENDIF
    IF (pole_nord) THEN
      ijbu = ij_begin + iip1
      ijbv = ij_begin
    ENDIF
    !$OMP END MASTER
    !$OMP BARRIER

    !    PRINT *,'---> on rentre dans guide_main'
    !    CALL AllGather_Field(ucov,ip1jmp1,llm)
    !    CALL AllGather_Field(vcov,ip1jm,llm)
    !    CALL AllGather_Field(teta,ip1jmp1,llm)
    !    CALL AllGather_Field(ps,ip1jmp1,1)
    !    CALL AllGather_Field(q,ip1jmp1,llm)

    !-----------------------------------------------------------------------
    ! Initialisations au premier passage
    !-----------------------------------------------------------------------

    IF (first) THEN
      first = .FALSE.
      !$OMP MASTER
      ALLOCATE(f_addu(ijb_u:ije_u, llm))
      ALLOCATE(f_addv(ijb_v:ije_v, llm))
      ALLOCATE(pk(iip1, jjb_u:jje_u, llm))
      ALLOCATE(pks(iip1, jjb_u:jje_u))
      ALLOCATE(p(ijb_u:ije_u, llmp1))
      CALL guide_init
      !$OMP END MASTER
      !$OMP BARRIER
      itau_test = 1001
      step_rea = 1
      count_no_rea = 0
      ! Calcul des constantes de rappel
      factt = dtvr * iperiod / daysec
      !$OMP MASTER
      CALL tau2alpha(3, iip1, jjb_v, jje_v, factt, tau_min_v, tau_max_v, alpha_v)
      CALL tau2alpha(2, iip1, jjb_u, jje_u, factt, tau_min_u, tau_max_u, alpha_u)
      CALL tau2alpha(1, iip1, jjb_u, jje_u, factt, tau_min_T, tau_max_T, alpha_T)
      CALL tau2alpha(1, iip1, jjb_u, jje_u, factt, tau_min_P, tau_max_P, alpha_P)
      CALL tau2alpha(1, iip1, jjb_u, jje_u, factt, tau_min_Q, tau_max_Q, alpha_Q)
      ! correction de rappel dans couche limite
      IF (guide_BL) THEN
        alpha_pcor(:) = 1.
      else
        do l = 1, llm
          alpha_pcor(l) = (1. + tanh(((plim_guide_BL - presnivs(l)) / preff) / 0.05)) / 2.
        enddo
      endif
      !$OMP END MASTER
      !$OMP BARRIER
      ! ini_anal: etat initial egal au guidage
      IF (ini_anal) THEN
        CALL guide_interp(ps, teta)
        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        DO l = 1, llm
          IF (guide_u) ucov(ijbu:ijeu, l) = ugui2(ijbu:ijeu, l)
          IF (guide_v) vcov(ijbv:ijev, l) = ugui2(ijbv:ijev, l)
          IF (guide_T) teta(ijbu:ijeu, l) = tgui2(ijbu:ijeu, l)
          IF (guide_Q) q(ijbu:ijeu, l) = qgui2(ijbu:ijeu, l)
        ENDDO

        IF (guide_P) THEN
          !$OMP MASTER
          ps(ijbu:ijeu) = psgui2(ijbu:ijeu)
          !$OMP END MASTER
          !$OMP BARRIER
          CALL pression_loc(ijnb_u, ap, bp, ps, p)
          CALL massdair_loc(p, masse)
          !$OMP BARRIER
        ENDIF
        RETURN
      ENDIF

    ENDIF !first

    !-----------------------------------------------------------------------
    ! Lecture des fichiers de guidage ?
    !-----------------------------------------------------------------------
    IF (iguide_read/=0) THEN
      ditau = real(itau)
      dday_step = real(day_step)
      IF (iguide_read<0) THEN
        tau = ditau / dday_step / REAL(iguide_read)
      ELSE
        tau = REAL(iguide_read) * ditau / dday_step
      ENDIF
      reste = tau - AINT(tau)
      IF (reste==0.) THEN
        IF (itau_test==itau) THEN
          WRITE(*, *)trim(modname) // ' second pass in advreel at itau=', &
                  itau
          CALL abort_gcm("guide_loc_lod", "stopped", 1)
        ELSE
          !$OMP MASTER
          IF (guide_v) vnat1(:, jjbv:jjev, :) = vnat2(:, jjbv:jjev, :)
          IF (guide_u) unat1(:, jjbu:jjeu, :) = unat2(:, jjbu:jjeu, :)
          IF (guide_T) tnat1(:, jjbu:jjeu, :) = tnat2(:, jjbu:jjeu, :)
          IF (guide_Q) qnat1(:, jjbu:jjeu, :) = qnat2(:, jjbu:jjeu, :)
          IF (guide_plevs==2) pnat1(:, jjbu:jjeu, :) = pnat2(:, jjbu:jjeu, :)
          IF (guide_P.OR.guide_plevs==1) psnat1(:, jjbu:jjeu) = psnat2(:, jjbu:jjeu)
          !$OMP END MASTER
          !$OMP BARRIER
          step_rea = step_rea + 1
          itau_test = itau
          IF (is_master) THEN
            WRITE(*, *)trim(modname) // ' Reading nudging files, step ', &
                    step_rea, 'after ', count_no_rea, ' skips'
          endif
          IF (guide_2D) THEN
            !$OMP MASTER
            CALL guide_read2D(step_rea)
            !$OMP END MASTER
            !$OMP BARRIER
          ELSE
            !$OMP MASTER
            CALL guide_read(step_rea)
            !$OMP END MASTER
            !$OMP BARRIER
          ENDIF
          count_no_rea = 0
        ENDIF
      ELSE
        count_no_rea = count_no_rea + 1

      ENDIF
    ENDIF !iguide_read=0

    !-----------------------------------------------------------------------
    ! Interpolation et conversion des champs de guidage
    !-----------------------------------------------------------------------
    IF (MOD(itau, iguide_int)==0) THEN
      CALL guide_interp(ps, teta)
    ENDIF
    ! Repartition entre 2 etats de guidage
    IF (iguide_read/=0) THEN
      tau = reste
    ELSE
      tau = 1.
    ENDIF

    !    CALL WriteField_u('ucov_guide',ucov)
    !    CALL WriteField_v('vcov_guide',vcov)
    !    CALL WriteField_u('teta_guide',teta)
    !    CALL WriteField_u('masse_guide',masse)


    !-----------------------------------------------------------------------
    !   Ajout des champs de guidage
    !-----------------------------------------------------------------------
    ! Sauvegarde du guidage?
    f_out = ((MOD(itau, iguide_sav)==0).AND.guide_sav)
    IF (f_out) THEN

      !$OMP BARRIER
      CALL pression_loc(ijnb_u, ap, bp, ps, p)

      !$OMP BARRIER
      IF (pressure_exner) THEN
        CALL exner_hyb_loc(ijnb_u, ps, p, pks, pk)
      else
        CALL exner_milieu_loc(ijnb_u, ps, p, pks, pk)
      endif

      !$OMP BARRIER

      unskap = 1. / kappa
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 1, llm
        DO j = jjbu, jjeu
          DO i = 1, iip1
            p(i + (j - 1) * iip1, l) = preff * (pk(i, j, l) / cpp) ** unskap
          ENDDO
        ENDDO
      ENDDO

      CALL guide_out("SP", jjp1, llm, p(ijb_u:ije_u, 1:llm), 1.)
    ENDIF

    IF (guide_u) THEN
      IF (guide_add) THEN
        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        DO l = 1, llm
          f_addu(ijbu:ijeu, l) = (1. - tau) * ugui1(ijbu:ijeu, l) + tau * ugui2(ijbu:ijeu, l)
        ENDDO
      else
        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        DO l = 1, llm
          f_addu(ijbu:ijeu, l) = (1. - tau) * ugui1(ijbu:ijeu, l) + tau * ugui2(ijbu:ijeu, l) - ucov(ijbu:ijeu, l)
        ENDDO
      endif

      !        CALL WriteField_u('f_addu',f_addu)

      IF (guide_zon) CALL guide_zonave_u(1, llm, f_addu)
      CALL guide_addfield_u(llm, f_addu, alpha_u)
      IF (f_out) CALL guide_out("ua", jjp1, llm, (1. - tau) * ugui1(ijb_u:ije_u, :) + tau * ugui2(ijb_u:ije_u, :), factt)
      IF (f_out) CALL guide_out("u", jjp1, llm, ucov(ijb_u:ije_u, :), factt)
      IF (f_out) THEN
        ! Ehouarn: fill the gaps adequately...
        IF (ijbu>ijb_u) f_addu(ijb_u:ijbu - 1, :) = 0
        IF (ijeu<ije_u) f_addu(ijeu + 1:ije_u, :) = 0
        CALL guide_out("ucov", jjp1, llm, f_addu(ijb_u:ije_u, :) / factt, factt)
      ENDIF
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 1, llm
        ucov(ijbu:ijeu, l) = ucov(ijbu:ijeu, l) + f_addu(ijbu:ijeu, l)
      ENDDO

    endif

    IF (guide_T) THEN
      IF (guide_add) THEN
        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        DO l = 1, llm
          f_addu(ijbu:ijeu, l) = (1. - tau) * tgui1(ijbu:ijeu, l) + tau * tgui2(ijbu:ijeu, l)
        ENDDO
      else
        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        DO l = 1, llm
          f_addu(ijbu:ijeu, l) = (1. - tau) * tgui1(ijbu:ijeu, l) + tau * tgui2(ijbu:ijeu, l) - teta(ijbu:ijeu, l)
        ENDDO
      endif
      IF (guide_zon) CALL guide_zonave_u(2, llm, f_addu)
      CALL guide_addfield_u(llm, f_addu, alpha_T)
      IF (f_out) CALL guide_out("teta", jjp1, llm, f_addu(:, :) / factt, factt)
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 1, llm
        teta(ijbu:ijeu, l) = teta(ijbu:ijeu, l) + f_addu(ijbu:ijeu, l)
      ENDDO
    endif

    IF (guide_P) THEN
      IF (guide_add) THEN
        !$OMP MASTER
        f_addu(ijbu:ijeu, 1) = (1. - tau) * psgui1(ijbu:ijeu) + tau * psgui2(ijbu:ijeu)
        !$OMP END MASTER
        !$OMP BARRIER
      else
        !$OMP MASTER
        f_addu(ijbu:ijeu, 1) = (1. - tau) * psgui1(ijbu:ijeu) + tau * psgui2(ijbu:ijeu) - ps(ijbu:ijeu)
        !$OMP END MASTER
        !$OMP BARRIER
      endif
      IF (guide_zon) CALL guide_zonave_u(2, 1, f_addu(ijb_u:ije_u, 1))
      CALL guide_addfield_u(1, f_addu(ijb_u:ije_u, 1), alpha_P)
      !       IF (f_out) CALL guide_out("ps",jjp1,1,f_addu(ijb_u:ije_u,1)/factt,factt)
      !$OMP MASTER
      ps(ijbu:ijeu) = ps(ijbu:ijeu) + f_addu(ijbu:ijeu, 1)
      !$OMP END MASTER
      !$OMP BARRIER
      CALL pression_loc(ijnb_u, ap, bp, ps, p)
      CALL massdair_loc(p, masse)
      !$OMP BARRIER
    endif

    IF (guide_Q) THEN
      IF (guide_add) THEN
        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        DO l = 1, llm
          f_addu(ijbu:ijeu, l) = (1. - tau) * qgui1(ijbu:ijeu, l) + tau * qgui2(ijbu:ijeu, l)
        ENDDO
      else
        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        DO l = 1, llm
          f_addu(ijbu:ijeu, l) = (1. - tau) * qgui1(ijbu:ijeu, l) + tau * qgui2(ijbu:ijeu, l) - q(ijbu:ijeu, l)
        ENDDO
      endif
      IF (guide_zon) CALL guide_zonave_u(2, llm, f_addu)
      CALL guide_addfield_u(llm, f_addu, alpha_Q)
      IF (f_out) CALL guide_out("q", jjp1, llm, f_addu(:, :) / factt, factt)

      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 1, llm
        q(ijbu:ijeu, l) = q(ijbu:ijeu, l) + f_addu(ijbu:ijeu, l)
      ENDDO
    endif

    IF (guide_v) THEN
      IF (guide_add) THEN
        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        DO l = 1, llm
          f_addv(ijbv:ijev, l) = (1. - tau) * vgui1(ijbv:ijev, l) + tau * vgui2(ijbv:ijev, l)
        ENDDO

      else
        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        DO l = 1, llm
          f_addv(ijbv:ijev, l) = (1. - tau) * vgui1(ijbv:ijev, l) + tau * vgui2(ijbv:ijev, l) - vcov(ijbv:ijev, l)
        ENDDO

      endif

      IF (guide_zon) CALL guide_zonave_v(2, jjm, llm, f_addv(ijb_v:ije_v, :))

      CALL guide_addfield_v(llm, f_addv(ijb_v:ije_v, :), alpha_v)
      IF (f_out) CALL guide_out("v", jjm, llm, vcov(ijb_v:ije_v, :), factt)
      IF (f_out) CALL guide_out("va", jjm, llm, (1. - tau) * vgui1(ijb_v:ije_v, :) + tau * vgui2(ijb_v:ije_v, :), factt)
      IF (f_out) THEN
        ! Ehouarn: Fill in the gaps adequately
        IF (ijbv>ijb_v) f_addv(ijb_v:ijbv - 1, :) = 0
        IF (ijev<ije_v) f_addv(ijev + 1:ije_v, :) = 0
        CALL guide_out("vcov", jjm, llm, f_addv(ijb_v:ije_v, :) / factt, factt)
      ENDIF

      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 1, llm
        vcov(ijbv:ijev, l) = vcov(ijbv:ijev, l) + f_addv(ijbv:ijev, l)
      ENDDO
    endif

  END SUBROUTINE guide_main


  SUBROUTINE guide_addfield_u(vsize, field, alpha)
    ! field1=a*field1+alpha*field2

    IMPLICIT NONE
    INCLUDE "dimensions.h"
    INCLUDE "paramet.h"

    ! input variables
    INTEGER, INTENT(IN) :: vsize
    REAL, DIMENSION(ijb_u:ije_u), INTENT(IN) :: alpha
    REAL, DIMENSION(ijb_u:ije_u, vsize), INTENT(INOUT) :: field

    ! Local variables
    INTEGER :: l

    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, vsize
      field(ijbu:ijeu, l) = alpha(ijbu:ijeu) * field(ijbu:ijeu, l) * alpha_pcor(l)
    ENDDO

  END SUBROUTINE guide_addfield_u


  SUBROUTINE guide_addfield_v(vsize, field, alpha)
    ! field1=a*field1+alpha*field2

    IMPLICIT NONE
    INCLUDE "dimensions.h"
    INCLUDE "paramet.h"

    ! input variables
    INTEGER, INTENT(IN) :: vsize
    REAL, DIMENSION(ijb_v:ije_v), INTENT(IN) :: alpha
    REAL, DIMENSION(ijb_v:ije_v, vsize), INTENT(INOUT) :: field

    ! Local variables
    INTEGER :: l

    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, vsize
      field(ijbv:ijev, l) = alpha(ijbv:ijev) * field(ijbv:ijev, l) * alpha_pcor(l)
    ENDDO

  END SUBROUTINE guide_addfield_v

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

  SUBROUTINE guide_zonave_u(typ, vsize, field)

    USE comconst_mod, ONLY: pi

    IMPLICIT NONE

    INCLUDE "dimensions.h"
    INCLUDE "paramet.h"
    INCLUDE "comgeom.h"

    ! input/output variables
    INTEGER, INTENT(IN) :: typ
    INTEGER, INTENT(IN) :: vsize
    REAL, DIMENSION(ijb_u:ije_u, vsize), INTENT(INOUT) :: field

    ! Local variables
    LOGICAL, SAVE :: first = .TRUE.
    !$OMP THREADPRIVATE(first)

    INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
    !$OMP THREADPRIVATE(imin,imax)
    INTEGER :: i, j, l, ij
    REAL, DIMENSION (iip1) :: lond       ! longitude in Deg.
    REAL, DIMENSION (jjb_u:jje_u, vsize) :: fieldm     ! zon-averaged field

    IF (first) THEN
      first = .FALSE.
      !Compute domain for averaging
      lond = rlonu * 180. / pi
      imin(1) = 1;imax(1) = iip1;
      imin(2) = 1;imax(2) = iip1;
      IF (guide_reg) THEN
        DO i = 1, iim
          IF (lond(i)<lon_min_g) imin(1) = i
          IF (lond(i)<=lon_max_g) imax(1) = i
        ENDDO
        lond = rlonv * 180. / pi
        DO i = 1, iim
          IF (lond(i)<lon_min_g) imin(2) = i
          IF (lond(i)<=lon_max_g) imax(2) = i
        ENDDO
      ENDIF
    ENDIF


    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, vsize
      fieldm(:, l) = 0.
      ! Compute zonal average

      !correction bug ici
      ! ---> a verifier
      ! ym         DO j=jjbv,jjev
      DO j = jjbu, jjeu
        DO i = imin(typ), imax(typ)
          ij = (j - 1) * iip1 + i
          fieldm(j, l) = fieldm(j, l) + field(ij, l)
        ENDDO
      ENDDO
      fieldm(:, l) = fieldm(:, l) / REAL(imax(typ) - imin(typ) + 1)
      ! Compute forcing
      DO j = jjbu, jjeu
        DO i = 1, iip1
          ij = (j - 1) * iip1 + i
          field(ij, l) = fieldm(j, l)
        ENDDO
      ENDDO
    ENDDO

  END SUBROUTINE guide_zonave_u


  SUBROUTINE guide_zonave_v(typ, hsize, vsize, field)

    USE comconst_mod, ONLY: pi

    IMPLICIT NONE

    INCLUDE "dimensions.h"
    INCLUDE "paramet.h"
    INCLUDE "comgeom.h"

    ! input/output variables
    INTEGER, INTENT(IN) :: typ
    INTEGER, INTENT(IN) :: vsize
    INTEGER, INTENT(IN) :: hsize
    REAL, DIMENSION(ijb_v:ije_v, vsize), INTENT(INOUT) :: field

    ! Local variables
    LOGICAL, SAVE :: first = .TRUE.
    !$OMP THREADPRIVATE(first)
    INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
    !$OMP THREADPRIVATE(imin, imax)
    INTEGER :: i, j, l, ij
    REAL, DIMENSION (iip1) :: lond       ! longitude in Deg.
    REAL, DIMENSION (jjb_v:jjev, vsize) :: fieldm     ! zon-averaged field

    IF (first) THEN
      first = .FALSE.
      !Compute domain for averaging
      lond = rlonu * 180. / pi
      imin(1) = 1;imax(1) = iip1;
      imin(2) = 1;imax(2) = iip1;
      IF (guide_reg) THEN
        DO i = 1, iim
          IF (lond(i)<lon_min_g) imin(1) = i
          IF (lond(i)<=lon_max_g) imax(1) = i
        ENDDO
        lond = rlonv * 180. / pi
        DO i = 1, iim
          IF (lond(i)<lon_min_g) imin(2) = i
          IF (lond(i)<=lon_max_g) imax(2) = i
        ENDDO
      ENDIF
    ENDIF

    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    DO l = 1, vsize
      ! Compute zonal average
      fieldm(:, l) = 0.
      DO j = jjbv, jjev
        DO i = imin(typ), imax(typ)
          ij = (j - 1) * iip1 + i
          fieldm(j, l) = fieldm(j, l) + field(ij, l)
        ENDDO
      ENDDO
      fieldm(:, l) = fieldm(:, l) / REAL(imax(typ) - imin(typ) + 1)
      ! Compute forcing
      DO j = jjbv, jjev
        DO i = 1, iip1
          ij = (j - 1) * iip1 + i
          field(ij, l) = fieldm(j, l)
        ENDDO
      ENDDO
    ENDDO

  END SUBROUTINE guide_zonave_v

  !=======================================================================
  SUBROUTINE guide_interp(psi, teta)
    USE exner_hyb_loc_m, ONLY: exner_hyb_loc
    USE exner_milieu_loc_m, ONLY: exner_milieu_loc
    USE parallel_lmdz
    USE mod_hallo
    USE Bands
    USE comconst_mod, ONLY: cpp, kappa
    USE comvert_mod, ONLY: preff, pressure_exner, bp, ap, disvert_type
    USE lmdz_q_sat, ONLY: q_sat
    IMPLICIT NONE

    include "dimensions.h"
    include "paramet.h"
    include "comgeom2.h"

    REAL, DIMENSION (iip1, jjb_u:jje_u), INTENT(IN) :: psi ! Psol gcm
    REAL, DIMENSION (iip1, jjb_u:jje_u, llm), INTENT(IN) :: teta ! Temp. Pot. gcm

    LOGICAL, SAVE :: first = .TRUE.
    !$OMP THREADPRIVATE(first)
    ! Variables pour niveaux pression:
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :, :) :: plnc1, plnc2 !niveaux pression guidage
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :, :) :: plunc, plsnc !niveaux pression modele
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :, :) :: plvnc       !niveaux pression modele
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :, :) :: p           ! pression intercouches
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :, :) :: pls, pext   ! var intermediaire
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :, :) :: pbarx
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :, :) :: pbary
    ! Variables pour fonction Exner (P milieu couche)
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :, :) :: pk
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :) :: pks
    REAL :: unskap
    ! Pression de vapeur saturante
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :) :: qsat
    !Variables intermediaires interpolation
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :, :) :: zu1, zu2
    REAL, ALLOCATABLE, SAVE, DIMENSION (:, :, :) :: zv1, zv2

    INTEGER :: i, j, l, ij
    CHARACTER(LEN = 20), PARAMETER :: modname = "guide_interp"
    TYPE(Request), SAVE :: Req
    !$OMP THREADPRIVATE(Req)

    IF (is_master) WRITE(*, *)trim(modname) // ': interpolate nudging variables'
    ! -----------------------------------------------------------------
    ! Calcul des niveaux de pression champs guidage (pour T et Q)
    ! -----------------------------------------------------------------
    IF (first) THEN
      !$OMP MASTER
      ALLOCATE(plnc1(iip1, jjb_u:jje_u, nlevnc))
      ALLOCATE(plnc2(iip1, jjb_u:jje_u, nlevnc))
      ALLOCATE(plunc(iip1, jjb_u:jje_u, llm))
      ALLOCATE(plsnc(iip1, jjb_u:jje_u, llm))
      ALLOCATE(plvnc(iip1, jjb_v:jje_v, llm))
      ALLOCATE(p(iip1, jjb_u:jje_u, llmp1))
      ALLOCATE(pls(iip1, jjb_u:jje_u, llm))
      ALLOCATE(pext(iip1, jjb_u:jje_u, llm))
      ALLOCATE(pbarx(iip1, jjb_u:jje_u, llm))
      ALLOCATE(pbary(iip1, jjb_v:jje_v, llm))
      ALLOCATE(pk(iip1, jjb_u:jje_u, llm))
      ALLOCATE(pks (iip1, jjb_u:jje_u))
      ALLOCATE(qsat(ijb_u:ije_u, llm))
      ALLOCATE(zu1(iip1, jjb_u:jje_u, llm))
      ALLOCATE(zu2(iip1, jjb_u:jje_u, llm))
      ALLOCATE(zv1(iip1, jjb_v:jje_v, llm))
      ALLOCATE(zv2(iip1, jjb_v:jje_v, llm))
      !$OMP END MASTER
      !$OMP BARRIER
    ENDIF

    IF (guide_plevs==0) THEN
      !$OMP DO
      DO l = 1, nlevnc
        DO j = jjbu, jjeu
          DO i = 1, iip1
            plnc2(i, j, l) = apnc(l)
            plnc1(i, j, l) = apnc(l)
          ENDDO
        ENDDO
      ENDDO
    ENDIF

    IF (first) THEN
      first = .FALSE.
      !$OMP MASTER
      WRITE(*, *)trim(modname) // ' : check vertical level order'
      WRITE(*, *)trim(modname) // ' LMDZ :'
      do l = 1, llm
        WRITE(*, *)trim(modname) // ' PL(', l, ')=', (ap(l) + ap(l + 1)) / 2. &
                + psi(1, jjeu) * (bp(l) + bp(l + 1)) / 2.
      enddo
      WRITE(*, *)trim(modname) // ' nudging file :'
      SELECT CASE (guide_plevs)
      CASE (0)
        do l = 1, nlevnc
          WRITE(*, *)trim(modname) // ' PL(', l, ')=', plnc2(1, jjbu, l)
        enddo
      CASE (1)
        DO l = 1, nlevnc
          WRITE(*, *)trim(modname) // ' PL(', l, ')=', &
                  apnc(l) + bpnc(l) * psnat2(1, jjbu)
        ENDDO
      CASE (2)
        do l = 1, nlevnc
          WRITE(*, *)trim(modname) // ' PL(', l, ')=', pnat2(1, jjbu, l)
        enddo
      END SELECT
      WRITE(*, *)trim(modname) // ' invert ordering: invert_p=', invert_p
      IF (guide_u) THEN
        do l = 1, nlevnc
          WRITE(*, *)trim(modname) // ' U(', l, ')=', unat2(1, jjbu, l)
        enddo
      endif
      IF (guide_T) THEN
        do l = 1, nlevnc
          WRITE(*, *)trim(modname) // ' T(', l, ')=', tnat2(1, jjbu, l)
        enddo
      endif
      !$OMP END MASTER
    endif ! of if (first)

    IF (guide_plevs /= 1 .OR. guide_t .AND. .NOT. guide_teta &
            .OR. guide_q .AND. guide_hr) THEN
      CALL pression_loc(ijnb_u, ap, bp, psi, p)
      IF (disvert_type==1) THEN
        CALL exner_hyb_loc(ijnb_u, psi, p, pks, pk)
      else ! we assume that we are in the disvert_type==2 case
        CALL exner_milieu_loc(ijnb_u, psi, p, pks, pk)
      endif
    end if

    ! -----------------------------------------------------------------
    ! Calcul niveaux pression modele
    ! -----------------------------------------------------------------

    !    ....  Calcul de pls , pression au milieu des couches ,en Pascals
    IF (guide_plevs==1) THEN
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 1, llm
        DO j = jjbu, jjeu
          DO i = 1, iip1
            pls(i, j, l) = (ap(l) + ap(l + 1)) / 2. + psi(i, j) * (bp(l) + bp(l + 1)) / 2.
          ENDDO
        ENDDO
      ENDDO
    ELSE
      unskap = 1. / kappa
      !$OMP BARRIER
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 1, llm
        DO j = jjbu, jjeu
          DO i = 1, iip1
            pls(i, j, l) = preff * (pk(i, j, l) / cpp) ** unskap
          ENDDO
        ENDDO
      ENDDO
    ENDIF

    !   calcul des pressions pour les grilles u et v
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    do l = 1, llm
      do j = jjbu, jjeu
        do i = 1, iip1
          pext(i, j, l) = pls(i, j, l) * aire(i, j)
        enddo
      enddo
    enddo

    CALL Register_Hallo_u(pext, llm, 1, 2, 2, 1, Req)
    CALL SendRequest(Req)
    !$OMP BARRIER
    CALL WaitRequest(Req)
    !$OMP BARRIER

    CALL massbar_loc(pext, pbarx, pbary)
    !$OMP BARRIER
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    do l = 1, llm
      do j = jjbu, jjeu
        do i = 1, iip1
          plunc(i, j, l) = pbarx(i, j, l) / aireu(i, j)
          plsnc(i, j, l) = pls(i, j, l)
        enddo
      enddo
    enddo
    !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    do l = 1, llm
      do j = jjbv, jjev
        do i = 1, iip1
          plvnc(i, j, l) = pbary(i, j, l) / airev(i, j)
        enddo
      enddo
    enddo

    ! -----------------------------------------------------------------
    ! Interpolation verticale champs guidage sur niveaux modele
    ! Conversion en variables gcm (ucov, vcov...)
    ! -----------------------------------------------------------------
    IF (guide_P) THEN
      !$OMP MASTER
      do j = jjbu, jjeu
        do i = 1, iim
          ij = (j - 1) * iip1 + i
          psgui1(ij) = psnat1(i, j)
          psgui2(ij) = psnat2(i, j)
        enddo
        psgui1(iip1 * j) = psnat1(1, j)
        psgui2(iip1 * j) = psnat2(1, j)
      enddo
      !$OMP END MASTER
      !$OMP BARRIER
    endif

    IF (guide_T) THEN
      ! Calcul des nouvelles valeurs des niveaux de pression du guidage
      IF (guide_plevs==1) THEN
        !$OMP DO
        DO l = 1, nlevnc
          DO j = jjbu, jjeu
            DO i = 1, iip1
              plnc2(i, j, l) = apnc(l) + bpnc(l) * psnat2(i, j)
              plnc1(i, j, l) = apnc(l) + bpnc(l) * psnat1(i, j)
            ENDDO
          ENDDO
        ENDDO
      ELSE IF (guide_plevs==2) THEN
        !$OMP DO
        DO l = 1, nlevnc
          DO j = jjbu, jjeu
            DO i = 1, iip1
              plnc2(i, j, l) = pnat2(i, j, l)
              plnc1(i, j, l) = pnat1(i, j, l)
            ENDDO
          ENDDO
        ENDDO
      ENDIF

      ! Interpolation verticale
      !$OMP MASTER
      CALL pres2lev(tnat1(:, jjbu:jjeu, :), zu1(:, jjbu:jjeu, :), nlevnc, llm, &
              plnc1(:, jjbu:jjeu, :), plsnc(:, jjbu:jjeu, :), iip1, jjnu, invert_p)
      CALL pres2lev(tnat2(:, jjbu:jjeu, :), zu2(:, jjbu:jjeu, :), nlevnc, llm, &
              plnc2(:, jjbu:jjeu, :), plsnc(:, jjbu:jjeu, :), iip1, jjnu, invert_p)
      !$OMP END MASTER
      !$OMP BARRIER
      ! Conversion en variables GCM
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      do l = 1, llm
        do j = jjbu, jjeu
          IF (guide_teta) THEN
            do i = 1, iim
              ij = (j - 1) * iip1 + i
              tgui1(ij, l) = zu1(i, j, l)
              tgui2(ij, l) = zu2(i, j, l)
            enddo
          ELSE
            do i = 1, iim
              ij = (j - 1) * iip1 + i
              tgui1(ij, l) = zu1(i, j, l) * cpp / pk(i, j, l)
              tgui2(ij, l) = zu2(i, j, l) * cpp / pk(i, j, l)
            enddo
          ENDIF
          tgui1(j * iip1, l) = tgui1((j - 1) * iip1 + 1, l)
          tgui2(j * iip1, l) = tgui2((j - 1) * iip1 + 1, l)
        enddo
        IF (pole_nord) THEN
          do i = 1, iip1
            tgui1(i, l) = tgui1(1, l)
            tgui2(i, l) = tgui2(1, l)
          enddo
        endif
        IF (pole_sud) THEN
          do i = 1, iip1
            tgui1(ip1jm + i, l) = tgui1(ip1jm + 1, l)
            tgui2(ip1jm + i, l) = tgui2(ip1jm + 1, l)
          enddo
        endif
      enddo
    ENDIF

    IF (guide_Q) THEN
      ! Calcul des nouvelles valeurs des niveaux de pression du guidage
      IF (guide_plevs==1) THEN
        !$OMP DO
        DO l = 1, nlevnc
          DO j = jjbu, jjeu
            DO i = 1, iip1
              plnc2(i, j, l) = apnc(l) + bpnc(l) * psnat2(i, j)
              plnc1(i, j, l) = apnc(l) + bpnc(l) * psnat1(i, j)
            ENDDO
          ENDDO
        ENDDO
      ELSE IF (guide_plevs==2) THEN
        !$OMP DO
        DO l = 1, nlevnc
          DO j = jjbu, jjeu
            DO i = 1, iip1
              plnc2(i, j, l) = pnat2(i, j, l)
              plnc1(i, j, l) = pnat1(i, j, l)
            ENDDO
          ENDDO
        ENDDO
      ENDIF

      ! Interpolation verticale
      !$OMP MASTER
      CALL pres2lev(qnat1(:, jjbu:jjeu, :), zu1(:, jjbu:jjeu, :), nlevnc, llm, &
              plnc1(:, jjbu:jjeu, :), plsnc(:, jjbu:jjeu, :), iip1, jjnu, invert_p)
      CALL pres2lev(qnat2(:, jjbu:jjeu, :), zu2(:, jjbu:jjeu, :), nlevnc, llm, &
              plnc2(:, jjbu:jjeu, :), plsnc(:, jjbu:jjeu, :), iip1, jjnu, invert_p)
      !$OMP END MASTER
      !$OMP BARRIER

      ! Conversion en variables GCM
      ! On suppose qu'on a la bonne variable dans le fichier de guidage:
      ! Hum.Rel si guide_hr, Hum.Spec. sinon.
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      do l = 1, llm
        do j = jjbu, jjeu
          do i = 1, iim
            ij = (j - 1) * iip1 + i
            qgui1(ij, l) = zu1(i, j, l)
            qgui2(ij, l) = zu2(i, j, l)
          enddo
          qgui1(j * iip1, l) = qgui1((j - 1) * iip1 + 1, l)
          qgui2(j * iip1, l) = qgui2((j - 1) * iip1 + 1, l)
        enddo
        IF (pole_nord) THEN
          do i = 1, iip1
            qgui1(i, l) = qgui1(1, l)
            qgui2(i, l) = qgui2(1, l)
          enddo
        endif
        IF (pole_sud) THEN
          do i = 1, iip1
            qgui1(ip1jm + i, l) = qgui1(ip1jm + 1, l)
            qgui2(ip1jm + i, l) = qgui2(ip1jm + 1, l)
          enddo
        endif
      enddo
      IF (guide_hr) THEN
        !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        do l = 1, llm
          CALL q_sat(iip1 * jjnu, teta(:, jjbu:jjeu, l) * pk(:, jjbu:jjeu, l) / cpp, &
                  plsnc(:, jjbu:jjeu, l), qsat(ijbu:ijeu, l))
          qgui1(ijbu:ijeu, l) = qgui1(ijbu:ijeu, l) * qsat(ijbu:ijeu, l) * 0.01 !hum. rel. en %
          qgui2(ijbu:ijeu, l) = qgui2(ijbu:ijeu, l) * qsat(ijbu:ijeu, l) * 0.01
        enddo

      ENDIF
    ENDIF

    IF (guide_u) THEN
      ! Calcul des nouvelles valeurs des niveaux de pression du guidage
      IF (guide_plevs==1) THEN
        !$OMP DO
        DO l = 1, nlevnc
          DO j = jjbu, jjeu
            DO i = 1, iim
              plnc2(i, j, l) = apnc(l) + bpnc(l) * (psnat2(i, j) * aire(i, j) * alpha1p2(i, j) &
                      + psnat2(i + 1, j) * aire(i + 1, j) * alpha3p4(i + 1, j)) / aireu(i, j)
              plnc1(i, j, l) = apnc(l) + bpnc(l) * (psnat1(i, j) * aire(i, j) * alpha1p2(i, j) &
                      + psnat1(i + 1, j) * aire(i + 1, j) * alpha3p4(i + 1, j)) / aireu(i, j)
            ENDDO
            plnc2(iip1, j, l) = plnc2(1, j, l)
            plnc1(iip1, j, l) = plnc1(1, j, l)
          ENDDO
        ENDDO
      ELSE IF (guide_plevs==2) THEN
        !$OMP DO
        DO l = 1, nlevnc
          DO j = jjbu, jjeu
            DO i = 1, iim
              plnc2(i, j, l) = (pnat2(i, j, l) * aire(i, j) * alpha1p2(i, j) &
                      + pnat2(i + 1, j, l) * aire(i, j) * alpha3p4(i + 1, j)) / aireu(i, j)
              plnc1(i, j, l) = (pnat1(i, j, l) * aire(i, j) * alpha1p2(i, j) &
                      + pnat1(i + 1, j, l) * aire(i, j) * alpha3p4(i + 1, j)) / aireu(i, j)
            ENDDO
            plnc2(iip1, j, l) = plnc2(1, j, l)
            plnc1(iip1, j, l) = plnc1(1, j, l)
          ENDDO
        ENDDO
      ENDIF

      ! Interpolation verticale
      !$OMP MASTER
      CALL pres2lev(unat1(:, jjbu:jjeu, :), zu1(:, jjbu:jjeu, :), nlevnc, llm, &
              plnc1(:, jjbu:jjeu, :), plunc(:, jjbu:jjeu, :), iip1, jjnu, invert_p)
      CALL pres2lev(unat2(:, jjbu:jjeu, :), zu2(:, jjbu:jjeu, :), nlevnc, llm, &
              plnc2(:, jjbu:jjeu, :), plunc(:, jjbu:jjeu, :), iip1, jjnu, invert_p)
      !$OMP END MASTER
      !$OMP BARRIER

      ! Conversion en variables GCM
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      do l = 1, llm
        do j = jjbu, jjeu
          do i = 1, iim
            ij = (j - 1) * iip1 + i
            ugui1(ij, l) = zu1(i, j, l) * cu(i, j)
            ugui2(ij, l) = zu2(i, j, l) * cu(i, j)
          enddo
          ugui1(j * iip1, l) = ugui1((j - 1) * iip1 + 1, l)
          ugui2(j * iip1, l) = ugui2((j - 1) * iip1 + 1, l)
        enddo
        IF (pole_nord) THEN
          do i = 1, iip1
            ugui1(i, l) = 0.
            ugui2(i, l) = 0.
          enddo
        endif
        IF (pole_sud) THEN
          do i = 1, iip1
            ugui1(ip1jm + i, l) = 0.
            ugui2(ip1jm + i, l) = 0.
          enddo
        endif
      enddo
    ENDIF

    IF (guide_v) THEN
      ! Calcul des nouvelles valeurs des niveaux de pression du guidage
      IF (guide_plevs==1) THEN
        CALL Register_Hallo_u(psnat1, 1, 1, 2, 2, 1, Req)
        CALL Register_Hallo_u(psnat2, 1, 1, 2, 2, 1, Req)
        CALL SendRequest(Req)
        !$OMP BARRIER
        CALL WaitRequest(Req)
        !$OMP BARRIER
        !$OMP DO
        DO l = 1, nlevnc
          DO j = jjbv, jjev
            DO i = 1, iip1
              plnc2(i, j, l) = apnc(l) + bpnc(l) * (psnat2(i, j) * aire(i, j) * alpha2p3(i, j) &
                      + psnat2(i, j + 1) * aire(i, j + 1) * alpha1p4(i, j + 1)) / airev(i, j)
              plnc1(i, j, l) = apnc(l) + bpnc(l) * (psnat1(i, j) * aire(i, j) * alpha2p3(i, j) &
                      + psnat1(i, j + 1) * aire(i, j + 1) * alpha1p4(i, j + 1)) / airev(i, j)
            ENDDO
          ENDDO
        ENDDO
      ELSE IF (guide_plevs==2) THEN
        CALL Register_Hallo_u(pnat1, llm, 1, 2, 2, 1, Req)
        CALL Register_Hallo_u(pnat2, llm, 1, 2, 2, 1, Req)
        CALL SendRequest(Req)
        !$OMP BARRIER
        CALL WaitRequest(Req)
        !$OMP BARRIER
        !$OMP DO
        DO l = 1, nlevnc
          DO j = jjbv, jjev
            DO i = 1, iip1
              plnc2(i, j, l) = (pnat2(i, j, l) * aire(i, j) * alpha2p3(i, j) &
                      + pnat2(i, j + 1, l) * aire(i, j) * alpha1p4(i, j + 1)) / airev(i, j)
              plnc1(i, j, l) = (pnat1(i, j, l) * aire(i, j) * alpha2p3(i, j) &
                      + pnat1(i, j + 1, l) * aire(i, j) * alpha1p4(i, j + 1)) / airev(i, j)
            ENDDO
          ENDDO
        ENDDO
      ENDIF
      ! Interpolation verticale

      !$OMP MASTER
      CALL pres2lev(vnat1(:, jjbv:jjev, :), zv1(:, jjbv:jjev, :), nlevnc, llm, &
              plnc1(:, jjbv:jjev, :), plvnc(:, jjbv:jjev, :), iip1, jjnv, invert_p)
      CALL pres2lev(vnat2(:, jjbv:jjev, :), zv2(:, jjbv:jjev, :), nlevnc, llm, &
              plnc2(:, jjbv:jjev, :), plvnc(:, jjbv:jjev, :), iip1, jjnv, invert_p)
      !$OMP END MASTER
      !$OMP BARRIER
      ! Conversion en variables GCM
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      do l = 1, llm
        do j = jjbv, jjev
          do i = 1, iim
            ij = (j - 1) * iip1 + i
            vgui1(ij, l) = zv1(i, j, l) * cv(i, j)
            vgui2(ij, l) = zv2(i, j, l) * cv(i, j)
          enddo
          vgui1(j * iip1, l) = vgui1((j - 1) * iip1 + 1, l)
          vgui2(j * iip1, l) = vgui2((j - 1) * iip1 + 1, l)
        enddo
      enddo
    ENDIF

  END SUBROUTINE guide_interp

  !=======================================================================
  SUBROUTINE tau2alpha(typ, pim, jjb, jje, factt, taumin, taumax, alpha)

    ! Calcul des constantes de rappel alpha (=1/tau)

    USE comconst_mod, ONLY: pi
    USE serre_mod, ONLY: clat, clon, grossismx, grossismy

    IMPLICIT NONE

    include "dimensions.h"
    include "paramet.h"
    include "comgeom2.h"

    ! input arguments :
    INTEGER, INTENT(IN) :: typ    ! u(2),v(3), ou scalaire(1)
    INTEGER, INTENT(IN) :: pim ! dimensions en lon
    INTEGER, INTENT(IN) :: jjb, jje ! dimensions en lat
    REAL, INTENT(IN) :: factt   ! pas de temps en fraction de jour
    REAL, INTENT(IN) :: taumin, taumax
    ! output arguments:
    REAL, DIMENSION(pim, jjb:jje), INTENT(OUT) :: alpha

    !  local variables:
    LOGICAL, SAVE :: first = .TRUE.
    REAL, SAVE :: gamma, dxdy_min, dxdy_max
    REAL, DIMENSION (iip1, jjp1) :: zdx, zdy
    REAL, DIMENSION (iip1, jjp1) :: dxdys, dxdyu
    REAL, DIMENSION (iip1, jjm) :: dxdyv
    REAL dxdy_
    REAL zlat, zlon
    REAL alphamin, alphamax, xi
    INTEGER i, j, ilon, ilat
    CHARACTER(LEN = 20), parameter :: modname = "tau2alpha"

    alphamin = factt / taumax
    alphamax = factt / taumin
    IF (guide_reg.OR.guide_add) THEN
      alpha = alphamax
      !-----------------------------------------------------------------------
      ! guide_reg: alpha=alpha_min dans region, 0. sinon.
      !-----------------------------------------------------------------------
      IF (guide_reg) THEN
        do j = jjb, jje
          do i = 1, pim
            IF (typ==2) THEN
              zlat = rlatu(j) * 180. / pi
              zlon = rlonu(i) * 180. / pi
            elseif (typ==1) THEN
              zlat = rlatu(j) * 180. / pi
              zlon = rlonv(i) * 180. / pi
            elseif (typ==3) THEN
              zlat = rlatv(j) * 180. / pi
              zlon = rlonv(i) * 180. / pi
            endif
            alpha(i, j) = alphamax / 16. * &
                    (1. + tanh((zlat - lat_min_g) / tau_lat)) * &
                    (1. + tanh((lat_max_g - zlat) / tau_lat)) * &
                    (1. + tanh((zlon - lon_min_g) / tau_lon)) * &
                    (1. + tanh((lon_max_g - zlon) / tau_lon))
          enddo
        enddo
      ENDIF
    ELSE
      !-----------------------------------------------------------------------
      ! Sinon, alpha varie entre alpha_min et alpha_max suivant le zoom.
      !-----------------------------------------------------------------------
      !Calcul de l'aire des mailles
      do j = 2, jjm
        do i = 2, iip1
          zdx(i, j) = 0.5 * (cu(i - 1, j) + cu(i, j)) / cos(rlatu(j))
        enddo
        zdx(1, j) = zdx(iip1, j)
      enddo
      do j = 2, jjm
        do i = 1, iip1
          zdy(i, j) = 0.5 * (cv(i, j - 1) + cv(i, j))
        enddo
      enddo
      do i = 1, iip1
        zdx(i, 1) = zdx(i, 2)
        zdx(i, jjp1) = zdx(i, jjm)
        zdy(i, 1) = zdy(i, 2)
        zdy(i, jjp1) = zdy(i, jjm)
      enddo
      do j = 1, jjp1
        do i = 1, iip1
          dxdys(i, j) = sqrt(zdx(i, j) * zdx(i, j) + zdy(i, j) * zdy(i, j))
        enddo
      enddo
      IF (typ==2) THEN
        do j = 1, jjp1
          do i = 1, iim
            dxdyu(i, j) = 0.5 * (dxdys(i, j) + dxdys(i + 1, j))
          enddo
          dxdyu(iip1, j) = dxdyu(1, j)
        enddo
      ENDIF
      IF (typ==3) THEN
        do j = 1, jjm
          do i = 1, iip1
            dxdyv(i, j) = 0.5 * (dxdys(i, j) + dxdys(i, j + 1))
          enddo
        enddo
      ENDIF
      ! Premier appel: calcul des aires min et max et de gamma.
      IF (first) THEN
        first = .FALSE.
        ! coordonnees du centre du zoom
        CALL coordij(clon, clat, ilon, ilat)
        ! aire de la maille au centre du zoom
        dxdy_min = dxdys(ilon, ilat)
        ! dxdy maximale de la maille
        dxdy_max = 0.
        do j = 1, jjp1
          do i = 1, iip1
            dxdy_max = max(dxdy_max, dxdys(i, j))
          enddo
        enddo
        ! Calcul de gamma
        IF (abs(grossismx - 1.)<0.1.OR.abs(grossismy - 1.)<0.1) THEN
          WRITE(*, *)trim(modname) // ' ATTENTION modele peu zoome'
          WRITE(*, *)trim(modname) // ' ATTENTION on prend une constante de guidage cste'
          gamma = 0.
        else
          gamma = (dxdy_max - 2. * dxdy_min) / (dxdy_max - dxdy_min)
          WRITE(*, *)trim(modname) // ' gamma=', gamma
          IF (gamma<1.e-5) THEN
            WRITE(*, *)trim(modname) // ' gamma =', gamma, '<1e-5'
            CALL abort_gcm("guide_loc_mod", "stopped", 1)
          endif
          gamma = log(0.5) / log(gamma)
          IF (gamma4) THEN
            gamma = min(gamma, 4.)
          endif
          WRITE(*, *)trim(modname) // ' gamma=', gamma
        endif
      ENDIF !first

      do j = jjb, jje
        do i = 1, pim
          IF (typ==1) THEN
            dxdy_ = dxdys(i, j)
            zlat = rlatu(j) * 180. / pi
          elseif (typ==2) THEN
            dxdy_ = dxdyu(i, j)
            zlat = rlatu(j) * 180. / pi
          elseif (typ==3) THEN
            dxdy_ = dxdyv(i, j)
            zlat = rlatv(j) * 180. / pi
          endif
          IF (abs(grossismx - 1.)<0.1.OR.abs(grossismy - 1.)<0.1) THEN
            ! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
            alpha(i, j) = alphamin
          else
            xi = ((dxdy_max - dxdy_) / (dxdy_max - dxdy_min))**gamma
            xi = min(xi, 1.)
            IF(lat_min_g<=zlat .AND. zlat<=lat_max_g) THEN
              alpha(i, j) = xi * alphamin + (1. - xi) * alphamax
            else
              alpha(i, j) = 0.
            endif
          endif
        enddo
      enddo
    ENDIF ! guide_reg

    IF (.NOT. guide_add) alpha = 1. - exp(- alpha)

  END SUBROUTINE tau2alpha

  !=======================================================================
  SUBROUTINE guide_read(timestep)

    IMPLICIT NONE

    include "dimensions.h"
    include "paramet.h"

    INTEGER, INTENT(IN) :: timestep

    LOGICAL, SAVE :: first = .TRUE.
    ! Identification fichiers et variables NetCDF:
    INTEGER, SAVE :: ncidu, varidu, ncidv, varidv, ncidp, varidp
    INTEGER, SAVE :: ncidQ, varidQ, ncidt, varidt, ncidps, varidps
    INTEGER :: ncidpl, varidpl, varidap, varidbp, dimid, lendim
    ! Variables auxiliaires NetCDF:
    INTEGER, DIMENSION(4) :: start, count
    INTEGER :: status, rcode
    CHARACTER (len = 80) :: abort_message
    CHARACTER (len = 20) :: modname = 'guide_read'
    CHARACTER (len = 20) :: namedim
    abort_message = 'pb in guide_read'

    ! -----------------------------------------------------------------
    ! Premier appel: initialisation de la lecture des fichiers
    ! -----------------------------------------------------------------
    IF (first) THEN
      ncidpl = -99
      WRITE(*, *) trim(modname) // ': opening nudging files '
      ! Ap et Bp si Niveaux de pression hybrides
      IF (guide_plevs==1) THEN
        WRITE(*, *) trim(modname) // ' Reading nudging on model levels'
        rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no file apbp.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no AP variable in file apbp.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no BP variable in file apbp.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        WRITE(*, *) trim(modname) // ' ncidpl,varidap', ncidpl, varidap
      endif

      ! Pression si guidage sur niveaux P variables
      IF (guide_plevs==2) THEN
        rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no file P.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no PRES variable in file P.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        WRITE(*, *) trim(modname) // ' ncidp,varidp', ncidp, varidp
        IF (ncidpl==-99) ncidpl = ncidp
      endif

      ! Vent zonal
      IF (guide_u) THEN
        rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no file u.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no UWND variable in file u.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        WRITE(*, *) trim(modname) // ' ncidu,varidu', ncidu, varidu
        IF (ncidpl==-99) ncidpl = ncidu

        status = nf90_inq_dimid(ncidu, "LONU", dimid)
        status = nf90_inquire_dimension(ncidu, dimid, namedim, lendim)
        IF (lendim /= iip1) THEN
          abort_message = 'dimension LONU different from iip1 in u.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF

        status = nf90_inq_dimid(ncidu, "LATU", dimid)
        status = nf90_inquire_dimension(ncidu, dimid, namedim, lendim)
        IF (lendim /= jjp1) THEN
          abort_message = 'dimension LATU different from jjp1 in u.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF

      endif

      ! Vent meridien
      IF (guide_v) THEN
        rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no file v.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no VWND variable in file v.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        WRITE(*, *) trim(modname) // ' ncidv,varidv', ncidv, varidv
        IF (ncidpl==-99) ncidpl = ncidv

        status = nf90_inq_dimid(ncidv, "LONV", dimid)
        status = nf90_inquire_dimension(ncidv, dimid, namedim, lendim)

        IF (lendim /= iip1) THEN
          abort_message = 'dimension LONV different from iip1 in v.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF

        status = nf90_inq_dimid(ncidv, "LATV", dimid)
        status = nf90_inquire_dimension(ncidv, dimid, namedim, lendim)
        IF (lendim /= jjm) THEN
          abort_message = 'dimension LATV different from jjm in v.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF

      endif

      ! Temperature
      IF (guide_T) THEN
        rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no file T.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no AIR variable in file T.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        WRITE(*, *) trim(modname) // ' ncidT,varidT', ncidt, varidt
        IF (ncidpl==-99) ncidpl = ncidt

        status = nf90_inq_dimid(ncidt, "LONV", dimid)
        status = nf90_inquire_dimension(ncidt, dimid, namedim, lendim)
        IF (lendim /= iip1) THEN
          abort_message = 'dimension LONV different from iip1 in T.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF

        status = nf90_inq_dimid(ncidt, "LATU", dimid)
        status = nf90_inquire_dimension(ncidt, dimid, namedim, lendim)
        IF (lendim /= jjp1) THEN
          abort_message = 'dimension LATU different from jjp1 in T.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF

      endif

      ! Humidite
      IF (guide_Q) THEN
        rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no file hur.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no RH variable in file hur.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        WRITE(*, *) trim(modname) // ' ncidQ,varidQ', ncidQ, varidQ
        IF (ncidpl==-99) ncidpl = ncidQ

        status = nf90_inq_dimid(ncidQ, "LONV", dimid)
        status = nf90_inquire_dimension(ncidQ, dimid, namedim, lendim)
        IF (lendim /= iip1) THEN
          abort_message = 'dimension LONV different from iip1 in hur.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF

        status = nf90_inq_dimid(ncidQ, "LATU", dimid)
        status = nf90_inquire_dimension(ncidQ, dimid, namedim, lendim)
        IF (lendim /= jjp1) THEN
          abort_message = 'dimension LATU different from jjp1 in hur.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF

      endif
      ! Pression de surface
      IF ((guide_P).OR.(guide_plevs==1)) THEN
        rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no file ps.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidps, 'SP', varidps)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no SP variable in file ps.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        WRITE(*, *) trim(modname) // ' ncidps,varidps', ncidps, varidps
      endif
      ! Coordonnee verticale
      IF (guide_plevs==0) THEN
        rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
        IF (rcode/=0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
        WRITE(*, *) trim(modname) // ' ncidpl,varidpl', ncidpl, varidpl
      endif
      ! Coefs ap, bp pour calcul de la pression aux differents niveaux
      IF (guide_plevs==1) THEN
        status = nf90_put_var(ncidpl, varidap, apnc, [1], [nlevnc])
        status = nf90_put_var(ncidpl, varidbp, bpnc, [1], [nlevnc])
      ELSEIF (guide_plevs==0) THEN
        status = nf90_put_var(ncidpl, varidpl, apnc, [1], [nlevnc])
        !FC Pour les corrections la pression est deja en Pascals on commente la ligne ci-dessous
        IF(convert_Pa) apnc = apnc * 100.! conversion en Pascals
        bpnc(:) = 0.
      ENDIF
      first = .FALSE.
    ENDIF ! (first)

    ! -----------------------------------------------------------------
    !   lecture des champs u, v, T, Q, ps
    ! -----------------------------------------------------------------

    !  dimensions pour les champs scalaires et le vent zonal
    start(1) = 1
    start(2) = jjb_u
    start(3) = 1
    start(4) = timestep

    count(1) = iip1
    count(2) = jjnb_u
    count(3) = nlevnc
    count(4) = 1

    IF (invert_y) start(2) = jjp1 - jje_u + 1
    ! Pression
    IF (guide_plevs==2) THEN
      status = nf90_put_var(ncidp, varidp, pnat2, start, count)
      IF (invert_y) THEN
        !           PRINT*,"Invertion impossible actuellement"
        !           CALL abort_gcm(modname,abort_message,1)
        CALL invert_lat(iip1, jjnb_u, nlevnc, pnat2)
      ENDIF
    endif

    !  Vent zonal
    IF (guide_u) THEN
      status = nf90_put_var(ncidu, varidu, unat2, start, count)
      IF (invert_y) THEN
        !           PRINT*,"Invertion impossible actuellement"
        !           CALL abort_gcm(modname,abort_message,1)
        CALL invert_lat(iip1, jjnb_u, nlevnc, unat2)
      ENDIF

    endif


    !  Temperature
    IF (guide_T) THEN
      status = nf90_put_var(ncidt, varidt, tnat2, start, count)
      IF (invert_y) THEN
        !           PRINT*,"Invertion impossible actuellement"
        !           CALL abort_gcm(modname,abort_message,1)
        CALL invert_lat(iip1, jjnb_u, nlevnc, tnat2)
      ENDIF
    endif

    !  Humidite
    IF (guide_Q) THEN
      status = nf90_put_var(ncidQ, varidQ, qnat2, start, count)
      IF (invert_y) THEN
        !           PRINT*,"Invertion impossible actuellement"
        !           CALL abort_gcm(modname,abort_message,1)
        CALL invert_lat(iip1, jjnb_u, nlevnc, qnat2)
      ENDIF

    endif

    !  Vent meridien
    IF (guide_v) THEN
      start(2) = jjb_v
      count(2) = jjnb_v
      IF (invert_y) start(2) = jjm - jje_v + 1

      status = nf90_put_var(ncidv, varidv, vnat2, start, count)
      IF (invert_y) THEN
        !           PRINT*,"Invertion impossible actuellement"
        !           CALL abort_gcm(modname,abort_message,1)
        CALL invert_lat(iip1, jjnb_v, nlevnc, vnat2)
      ENDIF
    endif

    !  Pression de surface
    IF ((guide_P).OR.(guide_plevs==1))  THEN
      start(2) = jjb_u
      start(3) = timestep
      start(4) = 0
      count(2) = jjnb_u
      count(3) = 1
      count(4) = 0
      IF (invert_y) start(2) = jjp1 - jje_u + 1
      status = nf90_put_var(ncidps, varidps, psnat2, start, count)
      IF (invert_y) THEN
        !           PRINT*,"Invertion impossible actuellement"
        !           CALL abort_gcm(modname,abort_message,1)
        CALL invert_lat(iip1, jjnb_u, 1, psnat2)
      ENDIF
    endif

  END SUBROUTINE guide_read

  !=======================================================================
  SUBROUTINE guide_read2D(timestep)

    IMPLICIT NONE

    include "dimensions.h"
    include "paramet.h"

    INTEGER, INTENT(IN) :: timestep

    LOGICAL, SAVE :: first = .TRUE.
    ! Identification fichiers et variables NetCDF:
    INTEGER, SAVE :: ncidu, varidu, ncidv, varidv, ncidp, varidp
    INTEGER, SAVE :: ncidQ, varidQ, ncidt, varidt, ncidps, varidps
    INTEGER :: ncidpl, varidpl, varidap, varidbp
    ! Variables auxiliaires NetCDF:
    INTEGER, DIMENSION(4) :: start, count
    INTEGER :: status, rcode
    ! Variables for 3D extension:
    REAL, DIMENSION (jjb_u:jje_u, llm) :: zu
    REAL, DIMENSION (jjb_v:jje_v, llm) :: zv
    INTEGER :: i
    CHARACTER (len = 80) :: abort_message
    CHARACTER (len = 20) :: modname = 'guide_read2D'
    abort_message = 'pb in guide_read2D'

    ! -----------------------------------------------------------------
    ! Premier appel: initialisation de la lecture des fichiers
    ! -----------------------------------------------------------------
    IF (first) THEN
      ncidpl = -99
      WRITE(*, *)trim(modname) // ' : opening nudging files '
      ! Ap et Bp si niveaux de pression hybrides
      IF (guide_plevs==1) THEN
        WRITE(*, *)trim(modname) // ' Reading nudging on model levels'
        rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no file apbp.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no AP variable in file apbp.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no BP variable in file apbp.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        WRITE(*, *)trim(modname) // 'ncidpl,varidap', ncidpl, varidap
      endif
      ! Pression
      IF (guide_plevs==2) THEN
        rcode = nf90_open('P.nc', nf90_nowrite, ncidp)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no file P.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidp, 'PRES', varidp)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no PRES variable in file P.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        WRITE(*, *)trim(modname) // ' ncidp,varidp', ncidp, varidp
        IF (ncidpl==-99) ncidpl = ncidp
      endif
      ! Vent zonal
      IF (guide_u) THEN
        rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no file u.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no UWND variable in file u.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        WRITE(*, *)trim(modname) // ' ncidu,varidu', ncidu, varidu
        IF (ncidpl==-99) ncidpl = ncidu
      endif

      ! Vent meridien
      IF (guide_v) THEN
        rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no file v.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no VWND variable in file v.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        WRITE(*, *)trim(modname) // ' ncidv,varidv', ncidv, varidv
        IF (ncidpl==-99) ncidpl = ncidv
      endif
      ! Temperature
      IF (guide_T) THEN
        rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no file T.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no AIR variable in file T.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        WRITE(*, *)trim(modname) // ' ncidT,varidT', ncidt, varidt
        IF (ncidpl==-99) ncidpl = ncidt
      endif
      ! Humidite
      IF (guide_Q) THEN
        rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no file hur.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no RH,variable in file hur.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        WRITE(*, *)trim(modname) // ' ncidQ,varidQ', ncidQ, varidQ
        IF (ncidpl==-99) ncidpl = ncidQ
      endif
      ! Pression de surface
      IF ((guide_P).OR.(guide_plevs==1)) THEN
        rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no file ps.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        rcode = nf90_inq_varid(ncidps, 'SP', varidps)
        IF (rcode/=nf90_noerr) THEN
          abort_message = 'Nudging: error -> no SP variable in file ps.nc'
          CALL abort_gcm(modname, abort_message, 1)
        ENDIF
        WRITE(*, *)trim(modname) // ' ncidps,varidps', ncidps, varidps
      endif
      ! Coordonnee verticale
      IF (guide_plevs==0) THEN
        rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
        IF (rcode/=0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
        WRITE(*, *)trim(modname) // ' ncidpl,varidpl', ncidpl, varidpl
      endif
      ! Coefs ap, bp pour calcul de la pression aux differents niveaux
      IF (guide_plevs==1) THEN
        status = nf90_put_var(ncidpl, varidap, apnc, [1], [nlevnc])
        status = nf90_put_var(ncidpl, varidbp, bpnc, [1], [nlevnc])
      elseif (guide_plevs==0) THEN
        status = nf90_put_var(ncidpl, varidpl, apnc, [1], [nlevnc])
        apnc = apnc * 100.! conversion en Pascals
        bpnc(:) = 0.
      endif
      first = .FALSE.
    endif ! (first)

    ! -----------------------------------------------------------------
    !   lecture des champs u, v, T, Q, ps
    ! -----------------------------------------------------------------

    !  dimensions pour les champs scalaires et le vent zonal
    start(1) = 1
    start(2) = jjb_u
    start(3) = 1
    start(4) = timestep

    count(1) = 1
    count(2) = jjnb_u
    count(3) = nlevnc
    count(4) = 1

    IF (invert_y) start(2) = jjp1 - jje_u + 1
    !  Pression
    IF (guide_plevs==2) THEN
      status = nf90_put_var(ncidp, varidp, zu, start, count)
      DO i = 1, iip1
        pnat2(i, :, :) = zu(:, :)
      ENDDO

      IF (invert_y) THEN
        !           PRINT*,"Invertion impossible actuellement"
        !           CALL abort_gcm(modname,abort_message,1)
        CALL invert_lat(iip1, jjnb_u, nlevnc, pnat2)
      ENDIF
    endif
    !  Vent zonal
    IF (guide_u) THEN
      status = nf90_put_var(ncidu, varidu, zu, start, count)
      DO i = 1, iip1
        unat2(i, :, :) = zu(:, :)
      ENDDO

      IF (invert_y) THEN
        !           PRINT*,"Invertion impossible actuellement"
        !           CALL abort_gcm(modname,abort_message,1)
        CALL invert_lat(iip1, jjnb_u, nlevnc, unat2)
      ENDIF
    endif


    !  Temperature
    IF (guide_T) THEN
      status = nf90_put_var(ncidt, varidt, zu, start, count)
      DO i = 1, iip1
        tnat2(i, :, :) = zu(:, :)
      ENDDO

      IF (invert_y) THEN
        !           PRINT*,"Invertion impossible actuellement"
        !           CALL abort_gcm(modname,abort_message,1)
        CALL invert_lat(iip1, jjnb_u, nlevnc, tnat2)
      ENDIF
    endif

    !  Humidite
    IF (guide_Q) THEN
      status = nf90_put_var(ncidQ, varidQ, zu, start, count)
      DO i = 1, iip1
        qnat2(i, :, :) = zu(:, :)
      ENDDO

      IF (invert_y) THEN
        !           PRINT*,"Invertion impossible actuellement"
        !           CALL abort_gcm(modname,abort_message,1)
        CALL invert_lat(iip1, jjnb_u, nlevnc, qnat2)
      ENDIF
    endif

    !  Vent meridien
    IF (guide_v) THEN
      start(2) = jjb_v
      count(2) = jjnb_v
      IF (invert_y) start(2) = jjm - jje_v + 1
      status = nf90_put_var(ncidv, varidv, zv, start, count)
      DO i = 1, iip1
        vnat2(i, :, :) = zv(:, :)
      ENDDO

      IF (invert_y) THEN

        !           PRINT*,"Invertion impossible actuellement"
        !           CALL abort_gcm(modname,abort_message,1)
        CALL invert_lat(iip1, jjnb_v, nlevnc, vnat2)
      ENDIF
    endif

    !  Pression de surface
    IF ((guide_P).OR.(guide_plevs==1))  THEN
      start(2) = jjb_u
      start(3) = timestep
      start(4) = 0
      count(2) = jjnb_u
      count(3) = 1
      count(4) = 0
      IF (invert_y) start(2) = jjp1 - jje_u + 1
      status = nf90_put_var(ncidps, varidps, zu(:, 1), start, count)
      DO i = 1, iip1
        psnat2(i, :) = zu(:, 1)
      ENDDO

      IF (invert_y) THEN
        !           PRINT*,"Invertion impossible actuellement"
        !           CALL abort_gcm(modname,abort_message,1)
        CALL invert_lat(iip1, jjnb_u, 1, psnat2)
      ENDIF
    endif

  END SUBROUTINE guide_read2D

  !=======================================================================
  SUBROUTINE guide_out(varname, hsize, vsize, field_loc, factt)
    USE parallel_lmdz
    USE mod_hallo, ONLY: gather_field_u, gather_field_v
    USE comconst_mod, ONLY: pi
    USE comvert_mod, ONLY: presnivs
    USE netcdf95, ONLY: nf95_def_var, nf95_put_var

    IMPLICIT NONE

    INCLUDE "dimensions.h"
    INCLUDE "paramet.h"
    INCLUDE "comgeom2.h"

    ! Variables entree
    CHARACTER*(*), INTENT(IN) :: varname
    INTEGER, INTENT (IN) :: hsize, vsize
    !   REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field_loc
    REAL, DIMENSION (:, :), INTENT(IN) :: field_loc
    REAL factt

    ! Variables locales
    INTEGER, SAVE :: timestep = 0
    ! Identites fichier netcdf
    INTEGER :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
    INTEGER :: vid_lonu, vid_lonv, vid_latu, vid_latv, vid_cu, vid_cv, vid_lev
    INTEGER :: vid_au, vid_av, varid_alpha_t, varid_alpha_q
    INTEGER, DIMENSION (3) :: dim3
    INTEGER, DIMENSION (4) :: dim4, count, start
    INTEGER :: ierr, varid, l
    REAL zu(ip1jmp1), zv(ip1jm), zt(iip1, jjp1), zq(iip1, jjp1)
    REAL, ALLOCATABLE, SAVE, DIMENSION(:, :, :) :: field_glo
    CHARACTER(LEN = 20), PARAMETER :: modname = "guide_out"

    !$OMP MASTER
    ALLOCATE(field_glo(iip1, hsize, vsize))
    !$OMP END MASTER
    !$OMP BARRIER

    !    WRITE(*,*)trim(modname)//' after allocation ',hsize,vsize

    IF (hsize==jjp1) THEN
      CALL gather_field_u(field_loc, field_glo, vsize)
    ELSE IF (hsize==jjm) THEN
      CALL gather_field_v(field_loc, field_glo, vsize)
    ENDIF

    !    WRITE(*,*)trim(modname)//' after gather '
    CALL Gather_field_u(alpha_u, zu, 1)
    CALL Gather_field_u(alpha_t, zt, 1)
    CALL Gather_field_u(alpha_q, zq, 1)
    CALL Gather_field_v(alpha_v, zv, 1)

    IF (mpi_rank >  0) THEN
      !$OMP MASTER
      DEALLOCATE(field_glo)
      !$OMP END MASTER
      !$OMP BARRIER

      RETURN
    ENDIF

    !$OMP MASTER
    IF (timestep==0) THEN
      ! ----------------------------------------------
      ! initialisation fichier de sortie
      ! ----------------------------------------------
      ! Ouverture du fichier
      ierr = nf90_create("guide_ins.nc", IOR(nf90_clobber, nf90_64bit_offset), nid)
      ! Definition des dimensions
      ierr = nf90_def_dim(nid, "LONU", iip1, id_lonu)
      ierr = nf90_def_dim(nid, "LONV", iip1, id_lonv)
      ierr = nf90_def_dim(nid, "LATU", jjp1, id_latu)
      ierr = nf90_def_dim(nid, "LATV", jjm, id_latv)
      ierr = nf90_def_dim(nid, "LEVEL", llm, id_lev)
      ierr = nf90_def_dim(nid, "TIME", nf90_unlimited, id_tim)

      ! Creation des variables dimensions
      ierr = nf90_def_var(nid, "LONU", nf90_float, id_lonu, vid_lonu)
      ierr = nf90_def_var(nid, "LONV", nf90_float, id_lonv, vid_lonv)
      ierr = nf90_def_var(nid, "LATU", nf90_float, id_latu, vid_latu)
      ierr = nf90_def_var(nid, "LATV", nf90_float, id_latv, vid_latv)
      ierr = nf90_def_var(nid, "LEVEL", nf90_float, id_lev, vid_lev)
      ierr = nf90_def_var(nid, "cu", nf90_float, (/id_lonu, id_latu/), vid_cu)
      ierr = nf90_def_var(nid, "cv", nf90_float, (/id_lonv, id_latv/), vid_cv)
      ierr = nf90_def_var(nid, "au", nf90_float, (/id_lonu, id_latu/), vid_au)
      ierr = nf90_def_var(nid, "av", nf90_float, (/id_lonv, id_latv/), vid_av)
      CALL nf95_def_var(nid, "alpha_T", nf90_float, (/id_lonv, id_latu/), &
              varid_alpha_t)
      CALL nf95_def_var(nid, "alpha_q", nf90_float, (/id_lonv, id_latu/), &
              varid_alpha_q)

      ierr = nf90_enddef(nid)

      ! Enregistrement des variables dimensions
      ierr = nf90_put_var(nid, vid_lonu, rlonu * 180. / pi)
      ierr = nf90_put_var(nid, vid_lonv, rlonv * 180. / pi)
      ierr = nf90_put_var(nid, vid_latu, rlatu * 180. / pi)
      ierr = nf90_put_var(nid, vid_latv, rlatv * 180. / pi)
      ierr = nf90_put_var(nid, vid_lev, presnivs)
      ierr = nf90_put_var(nid, vid_cu, cu)
      ierr = nf90_put_var(nid, vid_cv, cv)
      ierr = nf90_put_var(nid, vid_au, zu)
      ierr = nf90_put_var(nid, vid_av, zv)
      CALL nf95_put_var(nid, varid_alpha_t, zt)
      CALL nf95_put_var(nid, varid_alpha_q, zq)
      ! --------------------------------------------------------------------
      ! Création des variables sauvegardées
      ! --------------------------------------------------------------------
      ierr = nf90_redef(nid)
      ! Pressure (GCM)
      dim4 = (/id_lonv, id_latu, id_lev, id_tim/)
      ierr = nf90_def_var(nid, "SP", nf90_float, dim4, varid)
      ! Surface pressure (guidage)
      IF (guide_P) THEN
        dim3 = (/id_lonv, id_latu, id_tim/)
        ierr = nf90_def_var(nid, "ps", nf90_float, dim3, varid)
      ENDIF
      ! Zonal wind
      IF (guide_u) THEN
        dim4 = (/id_lonu, id_latu, id_lev, id_tim/)
        ierr = nf90_def_var(nid, "u", nf90_float, dim4, varid)
        ierr = nf90_def_var(nid, "ua", nf90_float, dim4, varid)
        ierr = nf90_def_var(nid, "ucov", nf90_float, dim4, varid)
      ENDIF
      ! Merid. wind
      IF (guide_v) THEN
        dim4 = (/id_lonv, id_latv, id_lev, id_tim/)
        ierr = nf90_def_var(nid, "v", nf90_float, dim4, varid)
        ierr = nf90_def_var(nid, "va", nf90_float, dim4, varid)
        ierr = nf90_def_var(nid, "vcov", nf90_float, dim4, varid)
      ENDIF
      ! Pot. Temperature
      IF (guide_T) THEN
        dim4 = (/id_lonv, id_latu, id_lev, id_tim/)
        ierr = nf90_def_var(nid, "teta", nf90_float, dim4, varid)
      ENDIF
      ! Specific Humidity
      IF (guide_Q) THEN
        dim4 = (/id_lonv, id_latu, id_lev, id_tim/)
        ierr = nf90_def_var(nid, "q", nf90_float, dim4, varid)
      ENDIF

      ierr = nf90_enddef(nid)
      ierr = nf90_close(nid)
    ENDIF ! timestep=0

    ! --------------------------------------------------------------------
    ! Enregistrement du champ
    ! --------------------------------------------------------------------

    ierr = nf90_open("guide_ins.nc", nf90_write, nid)

    IF (varname=="SP") timestep = timestep + 1

    ierr = nf90_inq_varid(nid, varname, varid)
    SELECT CASE (varname)
    CASE ("SP", "ps")
      start = (/1, 1, 1, timestep/)
      count = (/iip1, jjp1, llm, 1/)
    CASE ("v", "va", "vcov")
      start = (/1, 1, 1, timestep/)
      count = (/iip1, jjm, llm, 1/)
    CASE DEFAULT
      start = (/1, 1, 1, timestep/)
      count = (/iip1, jjp1, llm, 1/)
    END SELECT

    !$OMP END MASTER
    !$OMP BARRIER

    SELECT CASE (varname)

    CASE("u", "ua")
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 1, llm
        field_glo(:, 2:jjm, l) = field_glo(:, 2:jjm, l) / cu(:, 2:jjm)
        field_glo(:, 1, l) = 0. ; field_glo(:, jjp1, l) = 0.
      ENDDO
    CASE("v", "va")
      !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO l = 1, llm
        field_glo(:, :, l) = field_glo(:, :, l) / cv(:, :)
      ENDDO
    END SELECT

    !    if (varname=="ua") THEN
    !    CALL dump2d(iip1,jjp1,field_glo,'ua gui1 1ere couche ')
    !    CALL dump2d(iip1,jjp1,field_glo(:,:,llm),'ua gui1 llm ')
    !    endif

    !$OMP MASTER

    ierr = nf90_put_var(nid, varid, field_glo, start, count)
    ierr = nf90_close(nid)

    DEALLOCATE(field_glo)
    !$OMP END MASTER
    !$OMP BARRIER

  END SUBROUTINE guide_out


  !===========================================================================
  SUBROUTINE correctbid(iim, nl, x)
    INTEGER iim, nl
    REAL x(iim + 1, nl)
    INTEGER i, l
    REAL zz

    do l = 1, nl
      do i = 2, iim - 1
        IF(abs(x(i, l))>1.e10) THEN
          zz = 0.5 * (x(i - 1, l) + x(i + 1, l))
          PRINT*, 'correction ', i, l, x(i, l), zz
          x(i, l) = zz
        endif
      enddo
    enddo

  END SUBROUTINE  correctbid


  !====================================================================
  ! Ascii debug output. Could be reactivated
  !====================================================================

  SUBROUTINE dump2du(var, varname)
    USE parallel_lmdz
    USE mod_hallo
    IMPLICIT NONE
    include 'dimensions.h'
    include 'paramet.h'

    CHARACTER (len = *) :: varname

    REAL, DIMENSION(ijb_u:ije_u) :: var

    REAL, DIMENSION(ip1jmp1) :: var_glob

    RETURN

    CALL barrier
    CALL Gather_field_u(var, var_glob, 1)
    CALL barrier

    IF (mpi_rank==0) THEN
      CALL dump2d(iip1, jjp1, var_glob, varname)
    endif

    CALL barrier

  END SUBROUTINE  dump2du

  !====================================================================
  ! Ascii debug output. Could be reactivated
  !====================================================================
  SUBROUTINE dumpall
    IMPLICIT NONE
    include "dimensions.h"
    include "paramet.h"
    include "comgeom.h"
    CALL barrier
    CALL dump2du(alpha_u(ijb_u:ije_u), '  alpha_u couche 1')
    CALL dump2du(unat2(:, jjbu:jjeu, nlevnc), '  unat2 couche nlevnc')
    CALL dump2du(ugui1(ijb_u:ije_u, 1) * sqrt(unscu2(ijb_u:ije_u)), '  ugui1 couche 1')

  END SUBROUTINE  dumpall

  !===========================================================================
END MODULE guide_loc_mod