! $Id$

MODULE guide_mod

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

  USE getparam, ONLY: ini_getparam, fin_getparam, getpar
  USE 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 pres2lev_mod, 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
 
    IMPLICIT NONE
  
    INCLUDE "dimensions.h"
    INCLUDE "paramet.h"
    INCLUDE "iniprint.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
  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