!
! $Header$
!
      subroutine guide_pp(itau,ucov,vcov,teta,q,masse,ps)
      USE parallel

      IMPLICIT NONE

c      ......   Version  du 10/01/98    ..........

c             avec  coordonnees  verticales hybrides 
c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )

c=======================================================================
c
c   Auteur:  F.Hourdin
c   -------
c
c   Objet:
c   ------
c
c   GCM LMD nouvelle grille
c
c=======================================================================

c   ...  Dans inigeom , nouveaux calculs pour les elongations  cu , cv 
c        et possibilite d'appeler une fonction f(y)  a derivee tangente 
c        hyperbolique a la  place de la fonction a derivee sinusoidale.         

c   ...  Possibilite de choisir le shema de Van-leer pour l'advection de
c         q  , en faisant iadv = 10  dans   traceur  (29/04/97) .
c
c-----------------------------------------------------------------------
c   Declarations:
c   -------------

#include "dimensions.h"
#include "paramet.h"
#include "comconst.h"
#include "comdissnew.h"
#include "comvert.h"
#include "comgeom.h"
#include "logic.h"
#include "temps.h"
#include "control.h"
#include "ener.h"
#include "netcdf.inc"
#include "description.h"
#include "serre.h"
#include "tracstoke.h"
#include "guide.h"


c   variables dynamiques
      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
      REAL teta(ip1jmp1,llm)                 ! temperature potentielle 
      REAL q(ip1jmp1,llm)                 ! temperature potentielle 
      REAL ps(ip1jmp1)                       ! pression  au sol
      REAL masse(ip1jmp1,llm)                ! masse d'air

c   common passe pour des sorties
      real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)
      common/comdxdy/dxdys,dxdyu,dxdyv

c   variables dynamiques pour les reanalyses.
      REAL ucovrea1(ip1jmp1,llm),vcovrea1(ip1jm,llm) !vts cov reas
      REAL tetarea1(ip1jmp1,llm)             ! temp pot  reales
      REAL qrea1(ip1jmp1,llm)             ! temp pot  reales
      REAL masserea1(ip1jmp1,llm)             ! masse
      REAL psrea1(ip1jmp1)             ! ps
      REAL ucovrea2(ip1jmp1,llm),vcovrea2(ip1jm,llm) !vts cov reas
      REAL tetarea2(ip1jmp1,llm)             ! temp pot  reales
      REAL qrea2(ip1jmp1,llm)             ! temp pot  reales
      REAL masserea2(ip1jmp1,llm)             ! masse
      REAL psrea2(ip1jmp1)             ! ps
      real latmin

      real alpha_q(ip1jmp1)
      real alpha_T(ip1jmp1),alpha_P(ip1jmp1)
      real alpha_u(ip1jmp1),alpha_v(ip1jm)
      real dday_step,toto,reste,itau_test
      INTEGER step_rea,count_no_rea

c      real aire_min,aire_max
      integer ilon,ilat
      real factt,ztau(ip1jmp1)

      INTEGER itau,ij,l,i,j
      integer ncidpl,varidpl,nlev,status
      integer rcod,rid 
      real ditau,tau,a
      save nlev

c  TEST SUR QSAT
      real p(ip1jmp1,llmp1),pk(ip1jmp1,llm),pks(ip1jmp1)
      real pkf(ip1jmp1,llm)
      real pres(ip1jmp1,llm)
      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)

      real qsat(ip1jmp1,llm)
      real unskap
      real tnat(ip1jmp1,llm)
ccccccccccccccccc


      LOGICAL first
      save first
      data first/.true./

      save ucovrea1,vcovrea1,tetarea1,masserea1,psrea1,qrea1
      save ucovrea2,vcovrea2,tetarea2,masserea2,psrea2,qrea2

      save alpha_T,alpha_q,alpha_u,alpha_v,alpha_P,itau_test
      save step_rea,count_no_rea

      character*10 file
      integer igrads
      real dtgrads
      save igrads,dtgrads
      data igrads,dtgrads/2,100./
      integer :: ijb,ije,jjn

C-----------------------------------------------------------------------
c calcul de l'humidite saturante
C-----------------------------------------------------------------------
      ijb=ij_begin
      ije=ij_end
      jjn=jj_nb
      
      CALL pression_p( ip1jmp1, ap, bp, ps, p )
      call massdair_p(p,masse)
      CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
      tnat(ijb:ije,:)=pk(ijb:ije,:)*teta(ijb:ije,:)/cpp
      unskap   = 1./ kappa
      pres(ijb:ije,:)=preff*(pk(ijb:ije,:)/cpp)**unskap
      call q_sat(iip1*jjn*llm,tnat(ijb:ije,:),pres(ijb:ije,:),
     .            qsat(ijb:ije,:))

C-----------------------------------------------------------------------

c-----------------------------------------------------------------------
c   initialisations pour la lecture des reanalyses.
c    alpha determine la part des injections de donnees a chaque etape
c    alpha=1 signifie pas d'injection
c    alpha=0 signifie injection totale
c-----------------------------------------------------------------------

      if(online.eq.-1) then
          return
      endif

      if (first) then

         print*,'initialisation du guide '
         call conf_guide
         print*,'apres conf_guide'

         file='guide'
	 
	 if (mpi_rank==0) then
         call inigrads(igrads,iip1
     s  ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
     s  ,llm,presnivs,1.
     s  ,dtgrads,file,'dyn_zon ')
         endif
	 
         print*
     s   ,'1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'

         if(online.eq.-1) return
         if (online.eq.1) then

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Constantes de temps de rappel en jour
c  0.1 c'est en gros 2h30. 
c  1e10  est une constante infinie donc en gros pas de guidage
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c   coordonnees du centre du zoom
           call coordij(clon,clat,ilon,ilat)
c   aire de la maille au centre du zoom
           aire_min=aire(ilon+(ilat-1)*iip1)
c   aire maximale de la maille
           aire_max=0.
           do ij=1,ip1jmp1
              aire_max=max(aire_max,aire(ij))
           enddo
C  factt = pas de temps en fraction de jour
           factt=dtvr*iperiod/daysec

c     subroutine tau2alpha(type,im,jm,factt,taumin,taumax,alpha)
           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)

           call dump2d(iip1,jjp1,aire,'AIRE MAILLe ')
           call dump2d(iip1,jjp1,alpha_u,'COEFF U   ')
           call dump2d(iip1,jjp1,alpha_T,'COEFF T   ')

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c   Cas ou on force exactement par les variables analysees
         else
            alpha_T=0.
            alpha_u=0.
            alpha_v=0.
            alpha_P=0.
c           physic=.false.
         endif

         itau_test=1001
         step_rea=1
         count_no_rea=0
         ncidpl=-99
c    itau_test    montre si l'importation a deja ete faite au rang itau
c lecture d'un fichier netcdf pour determiner le nombre de niveaux
         IF (mpi_rank==0) THEN
	
         if (guide_u) then
           if (ncidpl.eq.-99) ncidpl=NCOPN('u.nc',NCNOWRIT,rcod)
         endif
c
         if (guide_v) then
           if (ncidpl.eq.-99) ncidpl=NCOPN('v.nc',NCNOWRIT,rcod)
         endif
c
         if (guide_T) then
           if (ncidpl.eq.-99) ncidpl=NCOPN('T.nc',NCNOWRIT,rcod)
         endif
c
         if (guide_Q) then
           if (ncidpl.eq.-99) ncidpl=NCOPN('hur.nc',NCNOWRIT,rcod)
         endif
c
         if (ncep) then
          status=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
         else
          status=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
         endif
          status=NF_INQ_DIMLEN(ncidpl,rid,nlev)
         print *,'nlev', nlev 
          call ncclos(ncidpl,rcod)
	 
	 ENDIF
	 
c   Lecture du premier etat des reanalyses.
         call Gather_Field(ps,ip1jmp1,1,0)

	 if (mpi_rank==0) then
	 
	 call read_reanalyse(1,ps
     s   ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)
         qrea2(:,:)=max(qrea2(:,:),0.1)
	 
	 endif
	 
	 call Broadcast_Field(ucovrea2,ip1jmp1,llm,0)
	 call Broadcast_Field(vcovrea2,ip1jm,llm,0)
	 call Broadcast_Field(tetarea2,ip1jmp1,llm,0)
	 call Broadcast_Field(qrea2,ip1jmp1,llm,0)
	 call Broadcast_Field(masserea2,ip1jmp1,llm,0)
	 call Broadcast_Field(psrea2,ip1jmp1,1,0)
	 


c-----------------------------------------------------------------------
c   Debut de l'integration temporelle:
c   ----------------------------------

      endif ! first
c
C-----------------------------------------------------------------------
C----- IMPORTATION DES VENTS,PRESSION ET TEMPERATURE REELS:
C-----------------------------------------------------------------------

      ditau=real(itau)
      DDAY_step=real(day_step)
      write(*,*)'ditau,dday_step'
      write(*,*)ditau,dday_step
      toto=4*ditau/dday_step
      reste=toto-aint(toto)
c     write(*,*)'toto,reste',toto,reste

      if (reste.eq.0.) then
        if (itau_test.eq.itau) then
          write(*,*)'deuxieme passage de advreel a itau=',itau
          stop
        else
        vcovrea1(:,:)=vcovrea2(:,:)
        ucovrea1(:,:)=ucovrea2(:,:)
        tetarea1(:,:)=tetarea2(:,:)
        qrea1(:,:)=qrea2(:,:)

          print*,'LECTURE REANALYSES, pas ',step_rea
     s         ,'apres ',count_no_rea,' non lectures'
           step_rea=step_rea+1
           itau_test=itau
	   
       call Gather_Field(ps,ip1jmp1,1,0)

       if (mpi_rank==0) then
           call read_reanalyse(step_rea,ps
     s     ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)
         qrea2(:,:)=max(qrea2(:,:),0.1)
       endif

       call Broadcast_Field(ucovrea2,ip1jmp1,llm,0)
       call Broadcast_Field(vcovrea2,ip1jm,llm,0)
       call Broadcast_Field(tetarea2,ip1jmp1,llm,0)
       call Broadcast_Field(qrea2,ip1jmp1,llm,0)
       call Broadcast_Field(masserea2,ip1jmp1,llm,0)
       call Broadcast_Field(psrea2,ip1jmp1,1,0)
       
       factt=dtvr*iperiod/daysec
       ztau(:)=factt/max(alpha_T(:),1.e-10)
      
       if (mpi_rank==0) then
         call wrgrads(igrads,1,aire   ,'aire      ','aire      ' )
         call wrgrads(igrads,1,dxdys  ,'dxdy      ','dxdy      ' )
         call wrgrads(igrads,1,alpha_u,'au        ','au        ' )
         call wrgrads(igrads,1,alpha_T,'at        ','at        ' )
         call wrgrads(igrads,1,ztau,'taut      ','taut      ' )
         call wrgrads(igrads,llm,ucov,'u         ','u         ' )
         call wrgrads(igrads,llm,ucovrea2,'ua        ','ua        ' )
         call wrgrads(igrads,llm,teta,'T         ','T         ' )
         call wrgrads(igrads,llm,tetarea2,'Ta        ','Ta        ' )
         call wrgrads(igrads,llm,qrea2,'Qa        ','Qa        ' )
         call wrgrads(igrads,llm,q,'Q         ','Q         ' )
         call wrgrads(igrads,llm,qsat,'QSAT      ','QSAT      ' )
      endif
      
        endif
      else
        count_no_rea=count_no_rea+1
      endif
 
C-----------------------------------------------------------------------
c   Guidage
c    x_gcm = a * x_gcm + (1-a) * x_reanalyses
C-----------------------------------------------------------------------

       if(ini_anal) print*,'ATTENTION !!! ON PART DU GUIDAGE'

      ditau=real(itau)
      dday_step=real(day_step)


      tau=4*ditau/dday_step
      tau=tau-aint(tau)

c  ucov
      ijb=ij_begin
      ije=ij_end
      
      if (guide_u) then
         do l=1,llm
            do ij=ijb,ije
                a=(1.-tau)*ucovrea1(ij,l)+tau*ucovrea2(ij,l)
                ucov(ij,l)=(1.-alpha_u(ij))*ucov(ij,l)+alpha_u(ij)*a
                if (first.and.ini_anal) ucov(ij,l)=a
            enddo
         enddo
      endif

c  teta
      if (guide_T) then
         do l=1,llm
            do ij=ijb,ije
                a=(1.-tau)*tetarea1(ij,l)+tau*tetarea2(ij,l)
                teta(ij,l)=(1.-alpha_T(ij))*teta(ij,l)+alpha_T(ij)*a
                if (first.and.ini_anal) teta(ij,l)=a
            enddo
         enddo
      endif

c  P
      if (guide_P) then
         do ij=ijb,ije
             a=(1.-tau)*psrea1(ij)+tau*psrea2(ij)
             ps(ij)=(1.-alpha_P(ij))*ps(ij)+alpha_P(ij)*a
             if (first.and.ini_anal) ps(ij)=a
         enddo
         CALL pression_p(ip1jmp1,ap,bp,ps,p)
         CALL massdair_p(p,masse)
      endif


c  q
      if (guide_Q) then
         do l=1,llm
            do ij=ijb,ije
                a=(1.-tau)*qrea1(ij,l)+tau*qrea2(ij,l)
c   hum relative en % -> hum specif
                a=qsat(ij,l)*a*0.01
                q(ij,l)=(1.-alpha_Q(ij))*q(ij,l)+alpha_Q(ij)*a
                if (first.and.ini_anal) q(ij,l)=a
            enddo
         enddo
      endif

c vcov
      if (pole_sud) ije=ij_end-iip1
      
      if (guide_v) then
         do l=1,llm
            do ij=ijb,ije
                a=(1.-tau)*vcovrea1(ij,l)+tau*vcovrea2(ij,l)
                vcov(ij,l)=(1.-alpha_v(ij))*vcov(ij,l)+alpha_v(ij)*a
                if (first.and.ini_anal) vcov(ij,l)=a
            enddo
         enddo
      endif

c     call dump2d(iip1,jjp1,tetarea1,'TETA REA 1     ')
c     call dump2d(iip1,jjp1,tetarea2,'TETA REA 2     ')
c     call dump2d(iip1,jjp1,teta,'TETA           ')

         first=.false.

      return
      end

c=======================================================================
      subroutine tau2alpha(type,pim,pjm,factt,taumin,taumax,alpha)
c=======================================================================

      implicit none

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

c   arguments :
      integer type
      integer pim,pjm
      real factt,taumin,taumax,dxdymin,dxdymax
      real dxdy_,alpha(pim,pjm)
      real dxdy_min,dxdy_max

c  local :
      real alphamin,alphamax,gamma,xi
      save gamma
      integer i,j,ilon,ilat

      logical first
      save first
      data first/.true./

      real cus(iip1,jjp1),cvs(iip1,jjp1)
      real cuv(iip1,jjm),cvu(iip1,jjp1)
      real zdx(iip1,jjp1),zdy(iip1,jjp1)

      real zlat
      real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)
      common/comdxdy/dxdys,dxdyu,dxdyv

      if (first) then
         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
         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
         do j=1,jjm
            do i=1,iip1
               dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
            enddo
         enddo

         call dump2d(iip1,jjp1,dxdys,'DX2DY2 SCAL  ')
         call dump2d(iip1,jjp1,dxdyu,'DX2DY2 U     ')
         call dump2d(iip1,jjp1,dxdyv,'DX2DY2 v     ')

c   coordonnees du centre du zoom
           call coordij(clon,clat,ilon,ilat)
c   aire de la maille au centre du zoom
           dxdy_min=dxdys(ilon,ilat)
c   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

         if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
             print*,'ATTENTION modele peu zoome'
             print*,'ATTENTION on prend une constante de guidage cste'
             gamma=0.
         else
            gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
            print*,'gamma=',gamma
            if (gamma.lt.1.e-5) then
              print*,'gamma =',gamma,'<1e-5'
              stop
            endif
            print*,'gamma=',gamma
            gamma=log(0.5)/log(gamma)
         endif
      endif

      alphamin=factt/taumax
      alphamax=factt/taumin

      do j=1,pjm
         do i=1,pim
            if (type.eq.1) then
               dxdy_=dxdys(i,j)
               zlat=rlatu(j)*180./pi
            elseif (type.eq.2) then
               dxdy_=dxdyu(i,j)
               zlat=rlatu(j)*180./pi
            elseif (type.eq.3) then
               dxdy_=dxdyv(i,j)
               zlat=rlatv(j)*180./pi
            endif
          if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
c  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_guide.le.zlat .and. zlat.le.lat_max_guide) then
               alpha(i,j)=xi*alphamin+(1.-xi)*alphamax
            else
               alpha(i,j)=0.
            endif
          endif
         enddo
      enddo


      return
      end