source: LMDZ.3.3/trunk/libf/dyn3d/nudge.F @ 5306

c $Header$
      subroutine nudge(itau,ucov,vcov,teta,masse,ps)

      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:  P. Le Van /L. Fairhead/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"


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

c   variables dynamiques pour les reanalyses.
      REAL ucovrea1(ip1jmp1,llm),vcovrea1(ip1jm,llm) !vts cov reas
      REAL tetarea1(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 masserea2(ip1jmp1,llm)             ! masse
      REAL psrea2(ip1jmp1)             ! ps
      real alphaT,alphaP,alphau,alphav
      real dday_step,toto,reste,itau_test
      INTEGER step_rea,count_no_rea


c
      INTEGER itau,ij,l
      integer ncidt,varidpl,nlev,status
      integer rcod,rid 
      real ditau,tau,a
      save nlev

      LOGICAL first,ncep
      integer online
      save first,online,ncep
      data first,online/.true.,1/
      data ncep/.false./

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

      save alphaT,alphau,alphav,alphaP,itau_test
      save step_rea,count_no_rea

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-----------------------------------------------------------------------

      print*,'ONLINE=',online
      if(online.eq.-1) then
          return
      endif

      if (first) then
         print*
     s   ,'1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
cnec         read(*,*) online
cnec      online=-1
cnec     print*,'Entrer les constantes de temps de rappel en jours'
cnec         read(*,*)alphaT
cnec         read(*,*)alphau
cnec         read(*,*)alphav
cnec         read(*,*)alphaP
             alphaT=0.1
             alphau=0.1
             alphav=0.1
             alphaP=1.e10
         print*,'alphaT,alphau,alphav,alphaP'
     s          ,alphaT,alphau,alphav,alphaP
         if(online.eq.-1) return
         print*,'alpha rappel pour T, u, v, P ',
     s    alphaT,alphau,alphav,alphaP
c   En fait, le coef alpha qu'on utilise est
c   x_gcm = a * x_gcm + ( 1 - a ) * x_reanalys
c   on a alors a=1-dt/tau
c   ou dt est la frequence a laquelle on corrige (ici iperiod*dtvr)
c   et tau la constante de temps lue. D'ou:
         if (online.eq.1) then
            alphaT=1.-iperiod*dtvr/(daysec*alphaT)
            alphau=1.-iperiod*dtvr/(daysec*alphau)
            alphav=1.-iperiod*dtvr/(daysec*alphav)
            alphaP=1.-iperiod*dtvr/(daysec*alphaP)
         else
            alphaT=0.
            alphau=0.
            alphav=0.
            alphaP=0.
c           physic=.false.
         endif
         print*,'alpha rappel pour T, u, v, P ',
     s    alphaT,alphau,alphav,alphaP

         itau_test=1001
         step_rea=1
         count_no_rea=0

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
         ncidt=NCOPN('T.nc',NCNOWRIT,rcod)
         if (ncep) then
          status=NF_INQ_DIMID(ncidt,'LEVEL',rid)
         else
          status=NF_INQ_DIMID(ncidt,'PRESSURE',rid)
         endif
          status=NF_INQ_DIMLEN(ncidt,rid,nlev)
         print *,'nlev', nlev 
          call ncclos(ncidt,rcod)
c   Lecture du premier etat des reanalyses.
         call read_reanalyse(1
     s   ,ucovrea2,vcovrea2,tetarea2,masserea2,ps,1,nlev)

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

         first=.false.
      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
           CALL SCOPY( ijmllm ,vcovrea2, 1, vcovrea1 , 1 )
           CALL SCOPY( ijp1llm,ucovrea2, 1, ucovrea1 , 1 )
           CALL SCOPY( ijp1llm,tetarea2,1,tetarea1   , 1 )
c          CALL SCOPY( ijp1llm,masserea2,1,masserea1   , 1 )
c          CALL SCOPY( ip1jmp1,psrea2, 1,  psrea1, 1 )

          print*,'LECTURE REANALYSES, pas ',step_rea
     s         ,'apres ',count_no_rea,' non lectures'
           step_rea=step_rea+1
           itau_test=itau
           call read_reanalyse(step_rea
     s     ,ucovrea2,vcovrea2,tetarea2,masserea2,ps,1,nlev)
        endif
      else
        count_no_rea=count_no_rea+1
      endif
 
C-----------------------------------------------------------------------
c   Corrrections
c
c    x_gcm = a * x_gcm + (1-a) * x_reanalyses
C-----------------------------------------------------------------------

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

c     print*,'alpha1',alphaT,alphaP,alphau,alphav

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

      do l=1,llm
         do ij=1,ip1jmp1
            a=(1.-tau)*ucovrea1(ij,l)+tau*ucovrea2(ij,l)
            ucov(ij,l)=alphau*ucov(ij,l)+(1-alphau)*a
            a=(1.-tau)*tetarea1(ij,l)+tau*tetarea2(ij,l)
            teta(ij,l)=alphaT*teta(ij,l)+(1-alphaT)*a
c           a=(1.-tau)*masserea1(ij,l)+tau*masserea2(ij,l)
c           masse(ij,l)=alphaP*masse(ij,l)+(1-alphaP)*a
         enddo
         do ij=1,ip1jm
            a=(1.-tau)*vcovrea1(ij,l)+tau*vcovrea2(ij,l)
            vcov(ij,l)=alphav*vcov(ij,l)+(1-alphav)*a
         enddo
      enddo

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           ')
Cmaf  on ne nudge pas sur la pression de surface
c      do ij=1,ip1jmp1
c         a=(1.-tau)*psrea1(ij)+tau*psrea2(ij)
c         ps(ij)=alphaP*ps(ij)+(1-alphaP)*a
c      enddo


      return
      end