! $Id$

MODULE guide_mod

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

  USE getparam, ONLY: ini_getparam, fin_getparam, getpar
  USE lmdz_write_field
  USE netcdf, ONLY: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close, &
          nf90_inq_dimid, nf90_inquire_dimension, nf90_float, nf90_def_var, &
          nf90_create, nf90_def_dim, nf90_open, nf90_unlimited, nf90_write, nf90_enddef, nf90_redef, &
          nf90_close, nf90_inq_varid, nf90_get_var, nf90_noerr, nf90_clobber, &
          nf90_64bit_offset, nf90_inq_dimid, nf90_inquire_dimension, nf90_put_var
  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

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

  SUBROUTINE guide_init

    USE netcdf, ONLY: nf90_noerr
    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
          CALL abort_gcm(modname, &
                  ' Nudging error -> no file u.nc', 1)
        endif
      endif

    elseif (guide_v) THEN
      IF (ncidpl==-99) THEN
        rcod = nf90_open('v.nc', nf90_nowrite, ncidpl)
        IF (rcod/=nf90_noerr) THEN
          CALL abort_gcm(modname, &
                  ' Nudging error -> no file v.nc', 1)
        endif
      endif
    elseif (guide_T) THEN
      IF (ncidpl==-99) THEN
        rcod = nf90_open('T.nc', nf90_nowrite, ncidpl)
        IF (rcod/=nf90_noerr) THEN
          CALL abort_gcm(modname, &
                  ' Nudging error -> no file T.nc', 1)
        endif
      endif
    elseif (guide_Q) THEN
      IF (ncidpl==-99) THEN
        rcod = nf90_open('hur.nc', nf90_nowrite, ncidpl)
        IF (rcod/=nf90_noerr) THEN
          CALL abort_gcm(modname, &
                  ' Nudging error -> no file hur.nc', 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
      CALL abort_gcm(modname, 'Nudging: error reading pressure levels', 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(ip1jmp1), stat = error)
    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    ALLOCATE(alpha_v(ip1jm), stat = error)
    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    ALLOCATE(alpha_T(iip1, jjp1), stat = error)
    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    ALLOCATE(alpha_Q(iip1, jjp1), stat = error)
    IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
    ALLOCATE(alpha_P(ip1jmp1), 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, jjp1, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(ugui1(ip1jmp1, llm), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(unat2(iip1, jjp1, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(ugui2(ip1jmp1, 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, jjp1, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(tgui1(ip1jmp1, llm), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(tnat2(iip1, jjp1, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(tgui2(ip1jmp1, 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, jjp1, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(qgui1(ip1jmp1, llm), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(qnat2(iip1, jjp1, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(qgui2(ip1jmp1, 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, jjm, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(vgui1(ip1jm, llm), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(vnat2(iip1, jjm, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(vgui2(ip1jm, 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, jjp1, nlevnc), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(pnat2(iip1, jjp1, 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, jjp1), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(psnat2(iip1, jjp1), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      psnat1 = 0.;psnat2 = 0.;
    ENDIF
    IF (guide_P) THEN
      ALLOCATE(psgui2(ip1jmp1), stat = error)
      IF (error /= 0) CALL abort_gcm(modname, abort_message, 1)
      ALLOCATE(psgui1(ip1jmp1), 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_m, ONLY: exner_hyb
    USE exner_milieu_m, ONLY: exner_milieu
    USE control_mod, ONLY: day_step, iperiod
    USE comconst_mod, ONLY: cpp, dtvr, daysec, kappa
    USE comvert_mod, ONLY: ap, bp, preff, presnivs, pressure_exner
    USE lmdz_iniprint, ONLY: lunout, prt_level


    IMPLICIT NONE

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


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

    ! Variables locales
    LOGICAL, SAVE :: first = .TRUE.
    LOGICAL :: f_out ! sortie guidage
    REAL, DIMENSION (ip1jmp1, llm) :: f_add ! var aux: champ de guidage
    REAL :: pk(ip1jmp1, llm) ! Exner at mid-layers
    REAL :: pks(ip1jmp1) ! Exner at the surface
    REAL :: unskap ! 1./kappa
    REAL, DIMENSION (ip1jmp1, llmp1) :: p ! Pressure at inter-layers
    ! Compteurs temps:
    INTEGER, SAVE :: step_rea, count_no_rea, itau_test ! lecture guidage
    REAL :: ditau, dday_step
    REAL :: tau, reste ! position entre 2 etats de guidage
    REAL, SAVE :: factt ! pas de temps en fraction de jour

    INTEGER :: l
    CHARACTER(LEN = 20) :: modname = "guide_main"
    CHARACTER (len = 80) :: abort_message


    !-----------------------------------------------------------------------
    ! Initialisations au premier passage
    !-----------------------------------------------------------------------
    IF (first) THEN
      first = .FALSE.
      CALL guide_init
      itau_test = 1001
      step_rea = 1
      count_no_rea = 0
      ! Calcul des constantes de rappel
      factt = dtvr * iperiod / daysec
      CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
      CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
      CALL tau2alpha(1, iip1, jjp1, factt, tau_min_T, tau_max_T, alpha_T)
      CALL tau2alpha(1, iip1, jjp1, factt, tau_min_P, tau_max_P, alpha_P)
      CALL tau2alpha(1, iip1, jjp1, 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
      ! ini_anal: etat initial egal au guidage
      IF (ini_anal) THEN
        CALL guide_interp(ps, teta)
        IF (guide_u) ucov = ugui2
        IF (guide_v) vcov = ugui2
        IF (guide_T) teta = tgui2
        IF (guide_Q) q = qgui2
        IF (guide_P) THEN
          ps = psgui2
          CALL pression(ip1jmp1, ap, bp, ps, p)
          CALL massdair(p, masse)
        ENDIF

      ENDIF
      ! Verification structure guidage
      !        IF (guide_u) THEN
      !            CALL writefield('unat',unat1)
      !            CALL writefield('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/)))
      !        ENDIF
      !        IF (guide_T) THEN
      !            CALL writefield('tnat',tnat1)
      !            CALL writefield('teta',RESHAPE(teta,(/iip1,jjp1,llm/)))
      !        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(lunout, *)trim(modname) // ' second pass in advreel at itau=', &
                  itau
          abort_message = 'stopped'
          CALL abort_gcm(modname, abort_message, 1)
        ELSE
          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
          step_rea = step_rea + 1
          itau_test = itau
          WRITE(*, *)trim(modname) // ' Reading nudging files, step ', &
                  step_rea, 'after ', count_no_rea, ' skips'
          IF (guide_2D) THEN
            CALL guide_read2D(step_rea)
          ELSE
            CALL guide_read(step_rea)
          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

    !-----------------------------------------------------------------------
    !   Ajout des champs de guidage
    !-----------------------------------------------------------------------
    ! Sauvegarde du guidage?
    f_out = ((MOD(itau, iguide_sav)==0).AND.guide_sav)
    IF (f_out) THEN
      ! compute pressures at layer interfaces
      CALL pression(ip1jmp1, ap, bp, ps, p)
      IF (pressure_exner) THEN
        CALL exner_hyb(ip1jmp1, ps, p, pks, pk)
      else
        CALL exner_milieu(ip1jmp1, ps, p, pks, pk)
      endif
      unskap = 1. / kappa
      ! Now compute pressures at mid-layer
      do l = 1, llm
        p(:, l) = preff * (pk(:, l) / cpp)**unskap
      enddo
      CALL guide_out("SP", jjp1, llm, p(:, 1:llm))
    ENDIF

    IF (guide_u) THEN
      IF (guide_add) THEN
        f_add = (1. - tau) * ugui1 + tau * ugui2
      else
        f_add = (1. - tau) * ugui1 + tau * ugui2 - ucov
      endif
      IF (guide_zon) CALL guide_zonave(1, jjp1, llm, f_add)
      CALL guide_addfield(ip1jmp1, llm, f_add, alpha_u)
      IF (f_out) CALL guide_out("ua", jjp1, llm, (1. - tau) * ugui1 + tau * ugui2)
      IF (f_out) CALL guide_out("u", jjp1, llm, ucov)
      IF (f_out) CALL guide_out("ucov", jjp1, llm, f_add / factt)
      ucov = ucov + f_add
    endif

    IF (guide_T) THEN
      IF (guide_add) THEN
        f_add = (1. - tau) * tgui1 + tau * tgui2
      else
        f_add = (1. - tau) * tgui1 + tau * tgui2 - teta
      endif
      IF (guide_zon) CALL guide_zonave(2, jjp1, llm, f_add)
      CALL guide_addfield(ip1jmp1, llm, f_add, alpha_T)
      IF (f_out) CALL guide_out("teta", jjp1, llm, f_add / factt)
      teta = teta + f_add
    endif

    IF (guide_P) THEN
      IF (guide_add) THEN
        f_add(1:ip1jmp1, 1) = (1. - tau) * psgui1 + tau * psgui2
      else
        f_add(1:ip1jmp1, 1) = (1. - tau) * psgui1 + tau * psgui2 - ps
      endif
      IF (guide_zon) CALL guide_zonave(2, jjp1, 1, f_add(1:ip1jmp1, 1))
      CALL guide_addfield(ip1jmp1, 1, f_add(1:ip1jmp1, 1), alpha_P)
      !        IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt)
      ps = ps + f_add(1:ip1jmp1, 1)
      CALL pression(ip1jmp1, ap, bp, ps, p)
      CALL massdair(p, masse)
    endif

    IF (guide_Q) THEN
      IF (guide_add) THEN
        f_add = (1. - tau) * qgui1 + tau * qgui2
      else
        f_add = (1. - tau) * qgui1 + tau * qgui2 - q
      endif
      IF (guide_zon) CALL guide_zonave(2, jjp1, llm, f_add)
      CALL guide_addfield(ip1jmp1, llm, f_add, alpha_Q)
      IF (f_out) CALL guide_out("q", jjp1, llm, f_add / factt)
      q = q + f_add
    endif

    IF (guide_v) THEN
      IF (guide_add) THEN
        f_add(1:ip1jm, :) = (1. - tau) * vgui1 + tau * vgui2
      else
        f_add(1:ip1jm, :) = (1. - tau) * vgui1 + tau * vgui2 - vcov
      endif
      IF (guide_zon) CALL guide_zonave(2, jjm, llm, f_add(1:ip1jm, :))
      CALL guide_addfield(ip1jm, llm, f_add(1:ip1jm, :), alpha_v)
      IF (f_out) CALL guide_out("v", jjm, llm, vcov)
      IF (f_out) CALL guide_out("va", jjm, llm, (1. - tau) * vgui1 + tau * vgui2)
      IF (f_out) CALL guide_out("vcov", jjm, llm, f_add(1:ip1jm, :) / factt)
      vcov = vcov + f_add(1:ip1jm, :)
    endif
  END SUBROUTINE guide_main

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

    IMPLICIT NONE

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

    ! Local variables
    INTEGER :: l

    do l = 1, vsize
      field(:, l) = alpha * field(:, l) * alpha_pcor(l)
    enddo

  END SUBROUTINE guide_addfield

  !=======================================================================
  SUBROUTINE guide_zonave(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(hsize * iip1, vsize), INTENT(INOUT) :: field

    ! Local variables
    LOGICAL, SAVE :: first = .TRUE.
    INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
    INTEGER :: i, j, l, ij
    REAL, DIMENSION (iip1) :: lond       ! longitude in Deg.
    REAL, DIMENSION (hsize, 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

    fieldm = 0.
    DO l = 1, vsize
      ! Compute zonal average
      DO j = 1, hsize
        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 = 1, hsize
        DO i = 1, iip1
          ij = (j - 1) * iip1 + i
          field(ij, l) = fieldm(j, l)
        ENDDO
      ENDDO
    ENDDO

  END SUBROUTINE guide_zonave

  !=======================================================================
  SUBROUTINE guide_interp(psi, teta)

    USE exner_hyb_m, ONLY: exner_hyb
    USE exner_milieu_m, ONLY: exner_milieu
    USE comconst_mod, ONLY: kappa, cpp
    USE comvert_mod, ONLY: preff, pressure_exner, bp, ap
    USE lmdz_q_sat, ONLY: q_sat
    IMPLICIT NONE

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

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

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

    INTEGER :: i, j, l, ij
    CHARACTER(LEN = 20), PARAMETER :: modname = "guide_interp"

    WRITE(*, *)trim(modname) // ': interpolate nudging variables'
    ! -----------------------------------------------------------------
    ! Calcul des niveaux de pression champs guidage
    ! -----------------------------------------------------------------
    IF (guide_modele) THEN
      do i = 1, iip1
        do j = 1, jjp1
          do l = 1, nlevnc
            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
      do i = 1, iip1
        do j = 1, jjp1
          do l = 1, nlevnc
            plnc2(i, j, l) = apnc(l)
            plnc1(i, j, l) = apnc(l)
          enddo
        enddo
      enddo

    END IF
    IF (first) THEN
      first = .FALSE.
      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, jjp1) * (bp(l) + bp(l + 1)) / 2.
      enddo
      WRITE(*, *)trim(modname) // ' nudging file :'
      do l = 1, nlevnc
        WRITE(*, *)trim(modname) // ' PL(', l, ')=', plnc2(1, 1, l)
      enddo
      WRITE(*, *)trim(modname) // ' invert ordering: invert_p=', invert_p
      IF (guide_u) THEN
        do l = 1, nlevnc
          WRITE(*, *)trim(modname) // ' U(', l, ')=', unat2(1, 1, l)
        enddo
      endif
      IF (guide_T) THEN
        do l = 1, nlevnc
          WRITE(*, *)trim(modname) // ' T(', l, ')=', tnat2(1, 1, l)
        enddo
      endif
    endif

    ! -----------------------------------------------------------------
    ! Calcul niveaux pression modele
    ! -----------------------------------------------------------------
    CALL pression(ip1jmp1, ap, bp, psi, p)
    IF (pressure_exner) THEN
      CALL exner_hyb(ip1jmp1, psi, p, pks, pk)
    else
      CALL exner_milieu(ip1jmp1, psi, p, pks, pk)
    endif
    !    ....  Calcul de pls , pression au milieu des couches ,en Pascals
    unskap = 1. / kappa
    prefkap = preff  ** kappa
    DO l = 1, llm
      DO j = 1, jjp1
        DO i = 1, iip1
          pls(i, j, l) = preff * (pk(i, j, l) / cpp) ** unskap
        ENDDO
      ENDDO
    ENDDO

    !   calcul des pressions pour les grilles u et v
    do l = 1, llm
      do j = 1, jjp1
        do i = 1, iip1
          pext(i, j, l) = pls(i, j, l) * aire(i, j)
        enddo
      enddo
    enddo
    CALL massbar(pext, pbarx, pbary)
    do l = 1, llm
      do j = 1, jjp1
        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
    do l = 1, llm
      do j = 1, jjm
        do i = 1, iip1
          plvnc(i, j, l) = pbary(i, j, l) / airev(i, j)
        enddo
      enddo
    enddo

    ! -----------------------------------------------------------------
    ! Interpolation champs guidage sur niveaux modele (+inversion N/S)
    ! Conversion en variables gcm (ucov, vcov...)
    ! -----------------------------------------------------------------
    IF (guide_P) THEN
      do j = 1, jjp1
        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
    endif

    IF (guide_u) THEN
      CALL pres2lev(unat1, zu1, nlevnc, llm, plnc1, plunc, iip1, jjp1, invert_p)
      CALL pres2lev(unat2, zu2, nlevnc, llm, plnc2, plunc, iip1, jjp1, invert_p)
      do l = 1, llm
        do j = 1, jjp1
          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
        do i = 1, iip1
          ugui1(i, l) = 0.
          ugui1(ip1jm + i, l) = 0.
          ugui2(i, l) = 0.
          ugui2(ip1jm + i, l) = 0.
        enddo
      enddo
    ENDIF

    IF (guide_T) THEN
      CALL pres2lev(tnat1, zu1, nlevnc, llm, plnc1, plsnc, iip1, jjp1, invert_p)
      CALL pres2lev(tnat2, zu2, nlevnc, llm, plnc2, plsnc, iip1, jjp1, invert_p)
      do l = 1, llm
        do j = 1, jjp1
          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
        do i = 1, iip1
          tgui1(i, l) = tgui1(1, l)
          tgui1(ip1jm + i, l) = tgui1(ip1jm + 1, l)
          tgui2(i, l) = tgui2(1, l)
          tgui2(ip1jm + i, l) = tgui2(ip1jm + 1, l)
        enddo
      enddo
    ENDIF

    IF (guide_v) THEN

      CALL pres2lev(vnat1, zv1, nlevnc, llm, plnc1(:, :jjm, :), plvnc, iip1, jjm, invert_p)
      CALL pres2lev(vnat2, zv2, nlevnc, llm, plnc2(:, :jjm, :), plvnc, iip1, jjm, invert_p)

      do l = 1, llm
        do j = 1, jjm
          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

    IF (guide_Q) THEN
      ! On suppose qu'on a la bonne variable dans le fichier de guidage:
      ! Hum.Rel si guide_hr, Hum.Spec. sinon.
      CALL pres2lev(qnat1, zu1, nlevnc, llm, plnc1, plsnc, iip1, jjp1, invert_p)
      CALL pres2lev(qnat2, zu2, nlevnc, llm, plnc2, plsnc, iip1, jjp1, invert_p)
      do l = 1, llm
        do j = 1, jjp1
          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
        do i = 1, iip1
          qgui1(i, l) = qgui1(1, l)
          qgui1(ip1jm + i, l) = qgui1(ip1jm + 1, l)
          qgui2(i, l) = qgui2(1, l)
          qgui2(ip1jm + i, l) = qgui2(ip1jm + 1, l)
        enddo
      enddo
      IF (guide_hr) THEN
        CALL q_sat(iip1 * jjp1 * llm, teta * pk / cpp, plsnc, qsat)
        qgui1 = qgui1 * qsat * 0.01 !hum. rel. en %
        qgui2 = qgui2 * qsat * 0.01
      ENDIF
    ENDIF

  END SUBROUTINE guide_interp

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

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

    USE comconst_mod, ONLY: pi
    USE serre_mod, ONLY: clon, clat, 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, pjm ! dimensions en lat, lon
    REAL, INTENT(IN) :: factt   ! pas de temps en fraction de jour
    REAL, INTENT(IN) :: taumin, taumax
    ! output arguments:
    REAL, DIMENSION(pim, pjm), 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"
    CHARACTER (len = 80) :: abort_message

    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 = 1, pjm
          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'
            abort_message = 'stopped'
            CALL abort_gcm(modname, abort_message, 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 = 1, pjm
        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

    ! -----------------------------------------------------------------
    ! Premier appel: initialisation de la lecture des fichiers
    ! -----------------------------------------------------------------
    IF (first) THEN
      ncidpl = -99
      WRITE(*, *) trim(modname) // ': opening nudging files '
      ! Niveaux de pression si non constants
      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_modele)) 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_get_var(ncidpl, varidap, apnc, [1], [nlevnc])
        status = nf90_get_var(ncidpl, varidbp, bpnc, [1], [nlevnc])
      ELSEIF (guide_plevs==0) THEN
        status = nf90_get_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) = 1
    start(3) = 1
    start(4) = timestep

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

    ! Pression
    IF (guide_plevs==2) THEN
      status = nf90_get_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, jjp1, nlevnc, pnat2)
      ENDIF
    endif

    !  Vent zonal
    IF (guide_u) THEN
      status = nf90_get_var(ncidu, varidu, unat2, start, count)
      IF (invert_y) THEN
        CALL invert_lat(iip1, jjp1, nlevnc, unat2)
      ENDIF
    endif

    !  Temperature
    IF (guide_T) THEN
      status = nf90_get_var(ncidt, varidt, tnat2, start, count)
      IF (invert_y) THEN
        CALL invert_lat(iip1, jjp1, nlevnc, tnat2)
      ENDIF
    endif

    !  Humidite
    IF (guide_Q) THEN
      status = nf90_get_var(ncidQ, varidQ, qnat2, start, count)
      IF (invert_y) THEN
        CALL invert_lat(iip1, jjp1, nlevnc, qnat2)
      ENDIF

    endif

    !  Vent meridien
    IF (guide_v) THEN
      count(2) = jjm
      status = nf90_get_var(ncidv, varidv, vnat2, start, count)
      IF (invert_y) THEN
        CALL invert_lat(iip1, jjm, nlevnc, vnat2)
      ENDIF
    endif

    !  Pression de surface
    IF ((guide_P).OR.(guide_modele))  THEN
      start(3) = timestep
      start(4) = 0
      count(2) = jjp1
      count(3) = 1
      count(4) = 0
      status = nf90_get_var(ncidps, varidps, psnat2, start, count)
      IF (invert_y) THEN
        CALL invert_lat(iip1, jjp1, 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 (jjp1, llm) :: zu
    REAL, DIMENSION (jjm, llm) :: zv
    INTEGER :: i

    CHARACTER (len = 80) :: abort_message
    CHARACTER (len = 20) :: modname = '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_modele)) 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_get_var(ncidpl, varidap, apnc, [1], [nlevnc])
        status = nf90_get_var(ncidpl, varidbp, bpnc, [1], [nlevnc])
      elseif (guide_plevs==0) THEN
        status = nf90_get_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) = 1
    start(3) = 1
    start(4) = timestep

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

    !  Pression
    IF (guide_plevs==2) THEN
      status = nf90_get_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, jjp1, nlevnc, pnat2)
      ENDIF
    endif
    !  Vent zonal
    IF (guide_u) THEN
      status = nf90_get_var(ncidu, varidu, zu, start, count)
      DO i = 1, iip1
        unat2(i, :, :) = zu(:, :)
      ENDDO

      IF (invert_y) THEN
        CALL invert_lat(iip1, jjp1, nlevnc, unat2)
      ENDIF

    endif

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

      IF (invert_y) THEN
        CALL invert_lat(iip1, jjp1, nlevnc, tnat2)
      ENDIF

    endif

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

      IF (invert_y) THEN
        CALL invert_lat(iip1, jjp1, nlevnc, qnat2)
      ENDIF

    endif

    !  Vent meridien
    IF (guide_v) THEN
      count(2) = jjm
      status = nf90_get_var(ncidv, varidv, zv, start, count)
      DO i = 1, iip1
        vnat2(i, :, :) = zv(:, :)
      ENDDO

      IF (invert_y) THEN
        CALL invert_lat(iip1, jjm, nlevnc, vnat2)
      ENDIF

    endif

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

      IF (invert_y) THEN
        CALL invert_lat(iip1, jjp1, 1, psnat2)
      ENDIF

    endif

  END SUBROUTINE guide_read2D

  !=======================================================================
  SUBROUTINE guide_out(varname, hsize, vsize, field)

    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

    ! 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, DIMENSION (iip1, hsize, vsize) :: field2
    CHARACTER(LEN = 20), PARAMETER :: modname = "guide_out"

    WRITE(*, *)trim(modname) // ': output timestep', timestep, 'var ', varname
    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, alpha_u)
      ierr = nf90_put_var(nid, vid_av, alpha_v)
      CALL nf95_put_var(nid, varid_alpha_t, alpha_t)
      CALL nf95_put_var(nid, varid_alpha_q, alpha_q)
      ! --------------------------------------------------------------------
      ! 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

    SELECT CASE (varname)
    CASE("u", "ua")
      DO l = 1, llm ; field2(:, 2:jjm, l) = field(:, 2:jjm, l) / cu(:, 2:jjm) ;
      ENDDO
      field2(:, 1, :) = 0. ; field2(:, jjp1, :) = 0.
    CASE("v", "va")
      DO l = 1, llm ; field2(:, :, l) = field(:, :, l) / cv(:, :) ;
      ENDDO
    CASE DEFAULT
      field2 = field
    END SELECT

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

  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

  !===========================================================================
END MODULE guide_mod