!
! $Header$
!
c
c
      subroutine read_reanalyse(timestep,psi
     s   ,u,v,t,q,masse,ps,mode,nlevnc)

c   mode=0 variables naturelles
c   mode=1 variabels GCM

       USE parallel
c -----------------------------------------------------------------
c   Declarations
c -----------------------------------------------------------------
      IMPLICIT NONE

c common
c ------

#include "netcdf.inc"
#include "dimensions.h"
#include "paramet.h"
#include "comgeom.h"
#include "comvert.h"
#include "guide.h"

c arguments
c ---------
      integer nlevnc
      integer timestep,mode

      real psi(iip1,jjp1)
      real u(iip1,jjp1,llm),v(iip1,jjm,llm)
      real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)
      real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)

c local
c -----
      integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
      integer ncidpl
      integer varidpl,ncidQ,varidQ,varidap,varidbp
      save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
      save ncidpl,varidap,varidbp
      save varidpl,ncidQ,varidQ

      real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
      real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
      real Qnc(iip1,jjp1,nlevnc)
      real plnc(iip1,jjp1,nlevnc)
      real apnc(nlevnc),bpnc(nlevnc)

      integer start(4),count(4),status
      integer i,j,l

      real rcode
      logical first
      save first
      INTEGER ierr

      data first/.true./

c -----------------------------------------------------------------
c   Initialisation de la lecture des fichiers
c -----------------------------------------------------------------
      if (first) then
           ncidpl=-99
           print*,'Intitialisation de read reanalsye'

c Niveaux de pression si non constants
            if (guide_modele) then
            print *,'Vous êtes entrain de lire des données sur 
     .               niveaux modèle'
            ncidpl=NCOPN('apbp.nc',NCNOWRIT,rcode)
            varidap=NCVID(ncidpl,'AP',rcode)
            varidbp=NCVID(ncidpl,'BP',rcode)
            print*,'ncidpl,varidap',ncidpl,varidap
            endif

c Vent zonal
            if (guide_u) then
               ncidu=NCOPN('u.nc',NCNOWRIT,rcode)
               varidu=NCVID(ncidu,'UWND',rcode)
               print*,'ncidu,varidu',ncidu,varidu
               if (ncidpl.eq.-99) ncidpl=ncidu
            endif

c Vent meridien
            if (guide_v) then
               ncidv=NCOPN('v.nc',NCNOWRIT,rcode)
               varidv=NCVID(ncidv,'VWND',rcode)
               print*,'ncidv,varidv',ncidv,varidv
               if (ncidpl.eq.-99) ncidpl=ncidv
            endif

c Temperature
            if (guide_T) then
               ncidt=NCOPN('T.nc',NCNOWRIT,rcode)
               varidt=NCVID(ncidt,'AIR',rcode)
               print*,'ncidt,varidt',ncidt,varidt
               if (ncidpl.eq.-99) ncidpl=ncidt
            endif

c Humidite
            if (guide_Q) then
               ncidQ=NCOPN('hur.nc',NCNOWRIT,rcode)
               varidQ=NCVID(ncidQ,'RH',rcode)
               print*,'ncidQ,varidQ',ncidQ,varidQ
               if (ncidpl.eq.-99) ncidpl=ncidQ
            endif

c Pression de surface
            if ((guide_P).OR.(guide_modele)) then
               ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)
               varidps=NCVID(ncidps,'SP',rcode)
               print*,'ncidps,varidps',ncidps,varidps
            endif

c Coordonnee verticale
            if (.not.guide_modele) then
              if (ncep) then
               print*,'Vous etes entrain de lire des donnees NCEP'
               varidpl=NCVID(ncidpl,'LEVEL',rcode)
            else
               print*,'Vous etes entrain de lire des donnees ECMWF'
               varidpl=NCVID(ncidpl,'PRESSURE',rcode)
              endif
              print*,'ncidpl,varidpl',ncidpl,varidpl
            endif
            print*,'ncidu,varidpl',ncidu,varidpl

      endif
      print*,'ok1'

c -----------------------------------------------------------------
c   lecture des champs u, v, T, ps
c -----------------------------------------------------------------

c  dimensions pour les champs scalaires et le vent zonal
c  -----------------------------------------------------

      start(1)=1
      start(2)=1
      start(3)=1
      start(4)=timestep

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

c mise a zero des tableaux 
c ------------------------
       unc(:,:,:)=0.
       vnc(:,:,:)=0.
       tnc(:,:,:)=0.
       Qnc(:,:,:)=0.
       plnc(:,:,:)=0.

c  Vent zonal
c  ----------

      if (guide_u) then
         print*,'avant la lecture de UNCEP nd de niv:',nlevnc

#ifdef NC_DOUBLE
         status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unc)
#else
         status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unc)
#endif
      print*,'WARNING!!! Correction bidon pour palier a un '
      print*,'probleme dans la creation des fichiers nc'
      call correctbid(iim,jjp1*nlevnc,unc)
      call dump2d(iip1,jjp1,unc,'UNC COUCHE 1 ')
      endif

c  Temperature
c  -----------

      print*,'ncidt=',ncidt,'varidt=',varidt,'start=',start
      print*,'count=',count
      if (guide_T) then
#ifdef NC_DOUBLE
         status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnc)
#else
         status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnc)
#endif
         call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 AAA ')
         call correctbid(iim,jjp1*nlevnc,tnc)
         call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 BBB ')
      endif

c  Humidite
c  --------

      if (guide_Q) then
#ifdef NC_DOUBLE
         status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,Qnc)
#else
         status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,Qnc)
#endif
         call correctbid(iim,jjp1*nlevnc,Qnc)
         call dump2d(iip1,jjp1,Qnc,'QNC COUCHE 1 ')
      endif

      count(2)=jjm
c  Vent meridien
c  -------------

      if (guide_v) then
#ifdef NC_DOUBLE
         status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnc)
#else
         status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnc)
#endif
         call correctbid(iim,jjm*nlevnc,vnc)
         call dump2d(iip1,jjm,vnc,'VNC COUCHE 1 ')
      endif

      start(3)=timestep
      start(4)=0
      count(2)=jjp1
      count(3)=1
      count(4)=0

c  Pression de surface
c  -------------------
      psnc(:,:)=0.
      if ((guide_P).OR.(guide_modele))  then
#ifdef NC_DOUBLE
         status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnc)
#else
         status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc)
#endif
         call dump2d(iip1,jjp1,psnc,'PSNC COUCHE 1 ')
         call correctbid(iim,jjp1,psnc)
      endif

c Calcul de la pression aux differents niveaux
c --------------------------------------------
      if (guide_modele) then
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
      status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
#else
      status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
      status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
#endif
      else
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
#else
      status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
#endif
c conversion en Pascals
      apnc=apnc*100.
      bpnc(:)=0.
      endif
      do i=1,iip1
        do j=1,jjp1
          do l=1,nlevnc
            plnc(i,j,l)=apnc(l)+bpnc(l)*psnc(i,j)
          enddo
        enddo
      enddo
      if (first) then
        print*,'Verification ordre niveaux verticaux'
        print*,'LMDZ :'
        do l=1,llm
           print*,'PL(',l,')=',(ap(l)+ap(l+1))/2.
     $                    +psi(1,jjp1)*(bp(l)+bp(l+1))/2.
        enddo
        print*,'Fichiers guidage'
        do l=1,nlevnc
           print*,'PL(',l,')=',plnc(1,1,l)
        enddo
        if (guide_u) then
          do l=1,nlevnc
            print*,'U(',l,')=',unc(1,1,l)
          enddo
        endif
        if (guide_T) then
          do l=1,nlevnc
            print*,'T(',l,')=',tnc(1,1,l)
          enddo
        endif
        print *,'inversion de l''ordre: ok_invertp=',ok_invertp
      endif

c -----------------------------------------------------------------
c  Interpollation verticale sur les niveaux modele
c -----------------------------------------------------------------
      call reanalyse2nat(nlevnc,psi,unc,vnc,tnc,Qnc,psnc,plnc,u,v,t,Q
     s    ,ps,masse,pk)

      call dump2d(iip1,jjm,v,'V COUCHE APRES ')


c -----------------------------------------------------------------
c  Passage aux variables du modele (vents covariants, temperature
c  potentielle et humidite specifique)
c -----------------------------------------------------------------
      call nat2gcm(u,v,t,Q,pk,u,v,t,Q)
      print*,'TIMESTEP ',timestep
      if(mode.ne.1) stop'mode pas egal 0'
c     call dump2d(iip1,jjm,v,'VCOV COUCHE 1 ')

c   Lignes introduites a une epoque pour un probleme oublie...
c     do l=1,llm
c        do i=1,iip1
c           v(i,1,l)=0.
c           v(i,jjm,l)=0.
c        enddo
c     enddo
      first=.false.

      return
      end


c===========================================================================
      subroutine reanalyse2nat(nlevnc,psi
     s   ,unc,vnc,tnc,qnc,psnc,plnc,u,v,t,q
     s   ,ps,masse,pk)
c===========================================================================

c -----------------------------------------------------------------
c   Inversion Nord/sud de la grille + interpollation sur les niveaux
c   verticaux du modele.
c -----------------------------------------------------------------

      implicit none

#include "dimensions.h"
#include "paramet.h"
#include "comgeom2.h"
#include "comvert.h"
#include "comconst.h"
#include "guide.h"

      integer nlevnc
      real psi(iip1,jjp1)
      real u(iip1,jjp1,llm),v(iip1,jjm,llm)
      real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)

      real plnc(iip1,jjp1,nlevnc)
      real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
      real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
      real qnc(iip1,jjp1,nlevnc)

      real zu(iip1,jjp1,llm),zv(iip1,jjm,llm)
      real zt(iip1,jjp1,llm),zq(iip1,jjp1,llm)

      real pext(iip1,jjp1,llm)
      real pbarx(iip1,jjp1,llm),pbary(iip1,jjm,llm)
      real plunc(iip1,jjp1,llm),plvnc(iip1,jjm,llm)
      real plsnc(iip1,jjp1,llm)

      REAL alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
      real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm),pks(iip1,jjp1)
      real pkf(iip1,jjp1,llm)
      real masse(iip1,jjp1,llm),pls(iip1,jjp1,llm)
      real prefkap,unskap


      integer i,j,l


c -----------------------------------------------------------------
c   calcul de la pression au milieu des couches.
c -----------------------------------------------------------------

      CALL pression( ip1jmp1, ap, bp, psi, p )
      call massdair(p,masse)
      CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)

c    ....  Calcul de pls , pression au milieu des couches ,en Pascals
      unskap=1./kappa
      prefkap =  preff  ** kappa
c     PRINT *,' Pref kappa unskap  ',preff,kappa,unskap
      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


c -----------------------------------------------------------------
c   calcul des pressions pour les grilles u et v
c -----------------------------------------------------------------

      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,jjp1+1-j,l)=pbarx(i,j,l)/aireu(i,j)
            plsnc(i,jjp1+1-j,l)=pls(i,j,l)
         enddo
      enddo
      enddo
      do l=1,llm
      do j=1,jjm
         do i=1,iip1
            plvnc(i,jjm+1-j,l)=pbary(i,j,l)/airev(i,j)
         enddo
      enddo
      enddo

c -----------------------------------------------------------------

      if (guide_P) then
      do j=1,jjp1
         do i=1,iim
            ps(i,j)=psnc(i,jjp1+1-j)
         enddo
         ps(iip1,j)=ps(1,j)
      enddo
      endif


c -----------------------------------------------------------------
      call pres2lev(unc,zu,nlevnc,llm,plnc,plunc,iip1,jjp1,ok_invertp)
      call pres2lev(vnc,zv,nlevnc,llm,plnc,plvnc,iip1,jjm,ok_invertp )
      call pres2lev(tnc,zt,nlevnc,llm,plnc,plsnc,iip1,jjp1,ok_invertp)
      call pres2lev(qnc,zq,nlevnc,llm,plnc,plsnc,iip1,jjp1,ok_invertp)

c     call dump2d(iip1,jjp1,ps,'PS    ')
c     call dump2d(iip1,jjp1,psu,'PS    ')
c     call dump2d(iip1,jjm,psv,'PS    ')
c  Inversion Nord/Sud
      do l=1,llm
         do j=1,jjp1
            do i=1,iim
               u(i,j,l)=zu(i,jjp1+1-j,l)
               t(i,j,l)=zt(i,jjp1+1-j,l)
               q(i,j,l)=zq(i,jjp1+1-j,l)
            enddo
            u(iip1,j,l)=u(1,j,l)
            t(iip1,j,l)=t(1,j,l)
            q(iip1,j,l)=q(1,j,l)
         enddo
      enddo

      do l=1,llm
         do j=1,jjm
            do i=1,iim
               v(i,j,l)=zv(i,jjm+1-j,l)
            enddo
            v(iip1,j,l)=v(1,j,l)
         enddo
      enddo

      return
      end

c===========================================================================
      subroutine nat2gcm(u,v,t,rh,pk,ucov,vcov,teta,q)
c===========================================================================

      implicit none

#include "dimensions.h"
#include "paramet.h"
#include "comgeom2.h"
#include "comconst.h"
#include "comvert.h"
#include "guide.h"

      real u(iip1,jjp1,llm),v(iip1,jjm,llm)
      real t(iip1,jjp1,llm),pk(iip1,jjp1,llm),rh(iip1,jjp1,llm)
      real ps(iip1,jjp1)

      real ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm)
      real teta(iip1,jjp1,llm),q(iip1,jjp1,llm)

      real pres(iip1,jjp1,llm),qsat(iip1,jjp1,llm)

      real unskap

      integer i,j,l


      print*,'Entree dans nat2gcm'
c    ucov(:,:,:)=0.
c    do l=1,llm
c       ucov(:,2:jjm,l)=u(:,2:jjm,l)*cu(:,2:jjm)
c    enddo
c    ucov(iip1,:,:)=ucov(1,:,:)

c    teta(:,:,:)=t(:,:,:)*cpp/pk(:,:,:)
c    teta(iip1,:,:)=teta(1,:,:)
     
c   calcul de ucov et de la temperature potentielle
      do l=1,llm
         do j=1,jjp1
            do i=1,iim
               ucov(i,j,l)=u(i,j,l)*cu(i,j)
               teta(i,j,l)=t(i,j,l)*cpp/pk(i,j,l)
            enddo
            ucov(iip1,j,l)=ucov(1,j,l)
            teta(iip1,j,l)=teta(1,j,l)
         enddo
         do i=1,iip1
            ucov(i,1,l)=0.
            ucov(i,jjp1,l)=0.
            teta(i,1,l)=teta(1,1,l)
            teta(i,jjp1,l)=teta(1,jjp1,l)
         enddo
      enddo

c   calcul de ucov
      do l=1,llm
         do j=1,jjm
            do i=1,iim
               vcov(i,j,l)=v(i,j,l)*cv(i,j)
            enddo
            vcov(iip1,j,l)=vcov(1,j,l)
         enddo
      enddo

c     call dump2d(iip1,jjp1,teta,'TETA EN BAS   ')
c     call dump2d(iip1,jjp1,teta(1,1,llm),'TETA EN HAUT   ')

c  Humidite relative -> specifique
c  -------------------------------
      if (guide_hr) then
c   FINALEMENT ON GUIDE EN HUMIDITE RELATIVE
      print*,'calcul de unskap'
      unskap   = 1./ kappa
      print*,'calcul de pres'
      pres(:,:,:)=preff*(pk(:,:,:)/cpp)**unskap
      print*,'calcul de qsat'
      call q_sat(iip1*jjp1*llm,t,pres,qsat)
      print*,'calcul de q'
c   ATTENTION : humidites relatives en %
      rh(:,:,:)=max(rh(:,:,:)*0.01,1.e-6)
      q(:,:,:)=qsat(:,:,:)*rh(:,:,:)
      print*,'calcul de q OK'
      call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1   ') 
      call dump2d(iip1,jjp1,pres,'PRESSION PREMIERE COUCHE   ')
      else
      q(:,:,:)=rh(:,:,:)
      print*,'calcul de q OK'
      call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1   ') 
      endif 

      return
      end



c===========================================================================
      subroutine correctbid(iim,nl,x)
c===========================================================================
      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)).gt.1.e10) then
               zz=0.5*(x(i-1,l)+x(i+1,l))
c              print*,'correction ',i,l,x(i,l),zz
               x(i,l)=zz
            endif
         enddo
      enddo
      return
      end