Index: /LMDZ.3.3/trunk/libf/dyn3d/offlinenc.F
===================================================================
--- /LMDZ.3.3/trunk/libf/dyn3d/offlinenc.F	(revision 187)
+++ /LMDZ.3.3/trunk/libf/dyn3d/offlinenc.F	(revision 187)
@@ -0,0 +1,932 @@
+      PROGRAM offlinenc
+      USE ioipsl
+      IMPLICIT NONE
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.h"
+#include "comdissnew.h"
+#include "comvert.h"
+#include "comgeom2.h"
+#include "logic.h"
+#include "temps.h"
+#include "control.h"
+#include "ener.h"
+#include "description.h"
+#include "serre.h"
+#include "tracstoke.h"
+      integer ngridmx,nbsrf
+      PARAMETER (ngridmx=iim*(jjm-1)+2,nbsrf=4)
+c-----------------------------------------------------------------------
+
+      INTEGER*4  iday ! jour julien
+      REAL       time ! Heure de la journee en fraction d'1 jour
+      REAL dtdyn,dtav
+      INTEGER iday_step
+
+      INTEGER         longcles
+      PARAMETER     ( longcles = 20 )
+      REAL  clesphy0( longcles )
+      logical clefinit
+
+c   variables dynamiques
+      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
+      REAL q(iip1,jjp1,llm,nqmx)               ! champs advectes
+	real qmoy(iip1,jjp1,nqmx)
+      REAL pks(iip1,jjp1)                      ! exner (f pour filtre)
+      REAL phis(iip1,jjp1)                     ! geopotentiel au sol
+      REAL phi(iip1,jjp1,llm)                  ! geopotentiel
+      REAL w(iip1,jjp1,llm)                    ! vitesse verticale
+C
+C   variables physiques
+      real zmfd(ngridmx,llm),zde_d(ngridmx,llm),zen_d(ngridmx,llm)
+      REAL qfi(ngridmx,llm,nqmx)               ! champs advectes
+      real zmfu(ngridmx,llm),zde_u(ngridmx,llm),zen_u(ngridmx,llm)
+	real t(ngridmx,llm)
+      real coefkz(ngridmx,llm)
+      real frac_impa(ngridmx,llm),frac_nucl(ngridmx,llm)
+      real pphis (ngridmx)
+      REAL airefi(ngridmx)
+      REAL latfi(ngridmx),lonfi(ngridmx) ! latitude et long pour chaque point
+      real yu1(ngridmx), yv1(ngridmx)
+      real ftsol(ngridmx,nbsrf),pctsrf(ngridmx,nbsrf)
+C
+      character*10 file
+c variables dynamiques intermediaire pour le transport
+      REAL pbaru(iip1,jjp1,llm),pbarv(iip1,jjm,llm) !flux de masse
+      REAL masse(iip1,jjp1,llm)
+      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
+      REAL teta(iip1,jjp1,llm)
+      INTEGER itau,irec,nrec,ldec,itauav
+      integer recmin,recmax,recstkmin,recstkmax
+      character*1 injecte,stoke
+      character*10 nom
+
+	real tinj2,tinj1
+	integer iiinj
+      integer jour0,isplit,nsplit_dyn,nsplit_phy,nsplit
+      logical debut,lectstart,rnpb
+
+      EXTERNAL inidissip,iniconst,inifilr
+      EXTERNAL inigeom
+      EXTERNAL exner_hyb,pression
+      EXTERNAL defrun_new
+
+      REAL time_0 ! temps de debut de simulation (% journee)
+
+      character*2 str2
+      INTEGER iq,j,i,l,nq,mode
+      real paprs(ngridmx,llm+1),pplay(ngridmx,llm),tempfi(ngridmx,llm)
+      real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm)
+      real lon_stat(nqmx),lat_stat(nqmx),zlon,zlat
+      integer i_stat(nqmx),j_stat(nqmx)
+      real intensite(nqmx)
+      integer iqmin,iqmax
+      integer npasinject
+      data npasinject/1/
+      integer histid, histvid, histaveid
+      real time_step, t_wrt, t_ops
+      integer itime_step
+      character*80 dynhist_file, dynhistave_file
+
+      logical avant
+      logical debutphy,redecoupage
+      DATA debutphy/.true./
+
+      real unskap,pksurcp,smass(iip1,jjp1)
+      integer ig0,ii,jj,iii,ij
+
+      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
+      REAL ps(iip1,jjp1),pkf(ip1jmp1,llm) 
+
+	integer ntau, iinjec2,iij
+	real zttt
+c-----------------------------------------------------------------------
+c   Initialisations:
+c   ----------------
+      descript = 'Run GCM LMDZ'
+	tinj1=0.
+	tinj2=0.
+      rad = 6400000
+      omeg = 7.272205e-05
+      g = 9.8
+      kappa = 0.285716
+      daysec = 86400
+      cpp = 1004.70885
+
+Cmaf a verifier....
+      preff = 101325.
+      pa= 50000.
+C
+      print*,'PARAMETRES A ENTRER DANS offline.def'
+      print*,'Nombre effectif de traceurs'
+      read(*,*) nq
+      write(*,*) nq
+      print*,'Splitting du pas de temps dynamique'
+      read(*,*) nsplit_dyn
+      write(*,*) nsplit_dyn
+      print*,'Splitting du pas de temps physique'
+      read(*,*) nsplit_phy
+      write(*,*) nsplit_phy
+      print*,'jour0, par rapport au debut du fichier fluxmass'
+      read(*,*) jour0
+      write(*,*) jour0
+      print*,'redecoupage (T) ou standard (F)'
+      read(*,'(l1)') redecoupage
+      write(*,*) redecoupage
+      print*,'Lecture fichier start (T)?'
+      read(*,'(l1)') lectstart
+      write(*,'(l1)') lectstart
+      print*,'Traceurs 1 et 2 = Rn222 et Pb'
+      read(*,'(l1)') rnpb
+      write(*,'(l1)') rnpb
+      print*,'Avant T ou arriere F?'
+      read(*,'(l1)') avant
+      write(*,'(l1)') avant
+      print*,'Index min et max des traceurs a initaliser'
+      read(*,*) iqmin,iqmax
+      write(*,*) iqmin,iqmax
+      print*,'localisation (lon, lat, intensite)'
+
+      if(iqmax.le.0) then
+         iqmin=iqmax+1
+      endif
+
+      do iq=iqmin,iqmax
+         read(*,*) lon_stat(iq),lat_stat(iq),intensite(iq)
+         write(*,*) lon_stat(iq),lat_stat(iq),intensite(iq)
+      enddo
+
+      print*,'Records min et max pour le stokage et couche '
+     s    ,'verticale (specifique reseaux), def 0 0 1'
+      read(*,*) recstkmin,recstkmax,ldec
+      write(*,*) recstkmin,recstkmax,ldec
+
+c   redefinition du splitting
+
+      if(mod(nsplit_dyn,nsplit_phy).eq.0) then
+         nsplit=nsplit_phy
+      elseif (mod(nsplit_phy,nsplit_dyn).eq.0) then
+         nsplit=nsplit_dyn
+      else
+         stop'prendre un des deux nsplit doit diviser l autre'
+      endif
+
+      nsplit_dyn=nsplit_dyn/nsplit
+      nsplit_phy=nsplit_phy/nsplit
+
+      print*,'nsplit,nsplit_dyn,nsplit_phy'
+     s       ,nsplit,nsplit_dyn,nsplit_phy
+
+c   Lecture eventuelle 
+      time_0=0.
+      if(nq.gt.nqmx) stop'surdimensioner nqmx'
+      print*,'av initial0'
+      call initial0(nq*ijp1llm,q)
+      print*,'av lectstart'
+
+      if(lectstart) then
+         CALL dynetat0("start.nc",nq,vcov,ucov,
+     .              teta,q,masse,ps,phis, time_0)
+         print*,'Lecture du start'
+      else
+         day_ini=0
+         anne_ini=0
+      endif
+
+      day_end=day_ini+nday
+C      print*,'av defrun'
+c     open (99,file='run.def',status='old',form='formatted')
+      clefinit=.true.
+      CALL defrun_new( 99, clefinit, clesphy0 )
+c     close(99)
+C      print*,'av iniconst'
+c   lecture du jour de demarrage
+c    premiere initialisation, eventuellement bidon
+      CALL iniconst
+      CALL inigeom
+C
+      print*,'ENTREE DANS redecoupenc ou lectfluxnc0
+     s        pour irec=0'
+
+c  Lecture des entetes des fichiers de flux
+
+      if(redecoupage) then
+         call redecoupenc(0,masse,pbaru,pbarv,w,teta,phi,
+     s     nrec,avant,airefi,pphis,
+     s     t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
+     s     yu1,yv1,ftsol,pctsrf,
+     s     frac_impa,frac_nucl,phis)
+      else
+         call lectfluxnc0(0,masse,pbaru,pbarv,w,teta,phi,
+     s     nrec,avant,airefi,pphis,
+     s     t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
+     s     yu1,yv1,ftsol,pctsrf,
+     s     frac_impa,frac_nucl,phis)
+         call inigeom
+      endif
+
+c   controle sur les records
+
+      if(avant) then
+         recmin=2+(jour0-1)*iday_step
+         recmax=1+(jour0+nday-1)*iday_step
+      else
+         recmin=(jour0-nday)*iday_step+2
+         recmax=jour0*iday_step+1
+      endif
+
+      if(recmin.le.1.or.recmax.ge.nrec) then
+         print*,'Probleme de lecture '
+         print*,'recmin,recmax,nrec',recmin,recmax,nrec
+         stop
+      endif
+
+C on recalcule le nombre de pas de temps dans une journee
+c   istdyn est la frequence de stokage en "pas de temps dynamique" dans
+c   le modele on-line. dtvr est le pas de temps du meme modele en s.
+c iday-step est le nombre de lectures par jour (grand pas de temps). 
+
+      iday_step=daysec/(istdyn*dtvr)
+C
+      CALL iniconst
+c  on calcule le pas de temps 
+      dtvr    = daysec/FLOAT(iday_step)
+      dtphys=dtvr/(nsplit*nsplit_phy) ! a mettre apres iniconst!
+      print*,'Pas de temps physique ',dtphys
+C
+      call gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
+      print *,'Pas de temps dynamique', dtvr
+      print*,'Pas de temps physique ',dtphys
+
+c-----------------------------------------------------------------------
+
+C maf est ce utile ?
+      call suphec
+      CALL inifilr
+      print*,'Pas de temps physique ',dtphys
+      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
+     *                tetagdiv, tetagrot , tetatemp              )
+      print*,'Pas de temps physique ',dtphys
+
+c-----------------------------------------------------------------------
+c   temps de depart et de fin:
+c   --------------------------
+      itau = 0
+      iday = dayref
+      time = time_0
+         IF(time.GT.1.) THEN
+          time = time-1.
+          iday = iday+1
+         ENDIF
+      print*,'Pas de temps physique ',dtphys
+c=======================================================================
+c   Initialisation des fichiers de sortie
+c=======================================================================
+      print*,'Initialisation de wrgras'
+	do i=1,iip1
+	print*,'rlonu(',i,')=',rlonu(i)
+	enddo
+	do i=1,jjp1
+	print*,'rlatu(',i,')=',rlatu(i)
+	enddo
+        do i=1,iip1
+        print*,'rlonv(',i,')=',rlonv(i)
+        enddo
+        do i=1,jjm
+        print*,'rlatv(',i,')=',rlatv(i)
+        enddo
+
+      file='sol'
+      call inigrads(1,iip1
+     s  ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
+     s  ,llm,presnivs,1.
+     s  ,dtphys,file,'gcmq2 ')
+
+      file='atm'
+      print*,'Pas de temps physique ',dtphys
+      call inigrads(2,iip1
+     s  ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
+     s  ,llm,presnivs,1.
+     s  ,dtphys,file,'gcmq2 ')
+
+      file='pbu'
+      call inigrads(3,iip1
+     s  ,rlonu,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
+     s  ,llm,presnivs,1.
+     s  ,dtphys,file,'gcmq2 ')
+
+      file='pbv'
+      call inigrads(4,iip1
+     s  ,rlonv,180./pi,-180.,180.,jjm,rlatv,-90.,90.,180./pi
+     s  ,llm,presnivs,1.
+     s  ,dtphys,file,'gcmq2 ')
+
+      file='qadv'
+      call inigrads(5,iip1
+     s  ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
+     s  ,llm,presnivs,1.
+     s  ,dtphys,file,'gcmq2 ')
+
+      print*,'dtphys ap inigrads ',dtphys
+
+      print*,'Ecriture du fichier de conditions aux limites'
+      open (98,file='limit',form='unformatted')
+      write(98) float(im),float(jm),float(nq),float(nday)
+      write(98) (rlonv(i)*180./pi,i=1,iip1)
+      write(98) (rlatu(j)*180./pi,j=1,jjp1)
+      write(98) ((phis(i,j)/g,i=1,iip1),j=1,jjp1)
+      close(98)
+
+      CALL dynredem0("restart.nc",day_end,anne_ini,phis,nq)
+
+c-----------------------------------------------------------------------
+c   Debut de la boucle en temps:
+c   ----------------------------
+
+      do itau=1,nday*iday_step
+          print*,'ITAU ITAU ITAU =',itau
+         injecte='.'
+         stoke='.'
+c Gestion du temps
+         if (avant) then
+            dtdyn=dtvr/float(nsplit*nsplit_dyn)
+c            irec=1+(jour0-1)*iday_step+itau
+	     irec=(jour0-1)*iday_step+itau
+         else
+            dtdyn=-dtvr/float(nsplit*nsplit_dyn)
+            irec=1+jour0*iday_step-itau+1
+         endif
+
+         debut=itau.eq.1
+
+c-----------------------------------------------------------------------
+c   Lecture des fichiers fluxmass et physique:
+c   -----------------------------------------------------
+         print*,'CCCCC Appel a REDECOUPENC ou LECTFLUXNC
+     s           au pas itau=',itau
+      if(redecoupage) then
+         call redecoupenc(irec,masse,pbaru,pbarv,w,teta,phi,
+     s     nrec,avant,airefi,pphis,
+     s     t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
+     s     yu1,yv1,ftsol,pctsrf, 
+     s     frac_impa,frac_nucl,phis)
+      else
+         call lectfluxnc0(irec,masse,pbaru,pbarv,w,teta,phi,
+     s     nrec,avant,airefi,pphis,
+     s     t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
+     s     yu1,yv1,ftsol,pctsrf, 
+     s     frac_impa,frac_nucl,phis)
+      endif
+
+c    ...  ouverture du fichier de stockage netcdf ...
+C
+      IF (debut) then
+         dynhistave_file = 'histmoy.nc'
+         day_ini=0
+         anne_ini=0
+         t_ops =(1./48.)*daysec
+         t_wrt =(1./48.)*daysec
+         dtav=dtvr/float(nsplit)
+         mode=1
+
+         CALL initdynav(dynhistave_file,day_ini,anne_ini,dtav,
+     .              t_ops, t_wrt, nq,mode, histaveid)
+
+         pi=2.*asin(1.)
+         do iq=iqmin,iqmax
+            do l=1,llm
+               do j=1,jjp1
+                  do i=1,iip1
+                     q(i,j,l,iq)=0.
+                  enddo
+               enddo
+            enddo
+
+            zlon=lon_stat(iq)*pi/180.
+            zlat=lat_stat(iq)*pi/180.
+
+            CALL coordij(zlon,zlat,i_stat(iq),j_stat(iq))
+
+         enddo
+      ENDIF
+
+C calcul de la pression de surface.
+
+      do j=1,jjp1
+         do i=1,iip1
+            smass(i,j)=0.
+         enddo
+      enddo
+
+      do l=1,llm
+        do j=1,jjp1
+          do i=1,iip1
+             smass(i,j)=smass(i,j)+masse(i,j,l)
+          enddo
+       enddo
+      enddo
+
+      do j=1,jjp1
+         do i=1,iip1
+            ps (i,j)=smass(i,j)/aire(i,j)*g
+         end do
+      end do
+C
+C calcul de la pression et de pk (fonction d'Exner)
+C
+      CALL pression ( ip1jmp1, ap, bp, ps, p       )
+      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+C
+      do jj=1,jjp1
+         do ii=1,iip1
+           phis(ii,jj)= phi(ii,jj,1)-teta(ii,jj,1)*
+     s                             (pks(ii,jj)-pk(ii,jj,1))
+         end do
+      end do
+c=======================================================================
+c   TERMES SOURCES OU PUITS
+         do isplit=1,nsplit
+c=======================================================================
+c injection par Fred
+c	Print*,Verifier en fct de iapp_tracvl la boucle inj
+       goto 333
+c Inj pour le stokage ttes les 0.5h
+        if (itau.gt.32.and.itau.le.56)then
+        print*,'Masse de la maille inj :',
+     s          masse(i_stat(1),j_stat(1),1)
+	do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (float(nsplit)*24.*20.)
+
+        tinj1=tinj1+intensite(1)/(float(nsplit)*24.*20.)
+	enddo
+
+           print*,'Apres injection a itau ',itau
+
+            do iq=1,nq
+               write(str2,'(i2.2)') iq
+              call diagadv(q(:,:,:,iq),masse,'Traceur '//str2)
+            enddo
+
+        print*,'QT TOTALE INJECTEE TR1',tinj1
+        print*,'fin des injections'
+
+        endif
+
+c Inj pour le stokage ttes les 1h
+
+        if (itau.gt.16.and.itau.le.28)then
+        print*,'Masse de la maille inj :',
+     s          masse(i_stat(1),j_stat(1),1)
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (float(nsplit)*12.*20.)
+
+        tinj1=tinj1+intensite(1)/(float(nsplit)*12.*20.)
+
+        enddo
+
+           print*,'Apres injection a itau ',itau
+
+            do iq=1,nq
+               write(str2,'(i2.2)') iq
+              call diagadv(q(:,:,:,iq),masse,'Traceur '//str2)
+            enddo
+
+        print*,'QT TOTALE INJECTEE TR1',tinj1
+        print*,'fin des injections'
+
+        endif
+
+c Inj pour le stokage ttes les 1.5h
+
+        if (itau.eq.11.and.isplit.gt.2)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+        endif
+
+        if (itau.ge.12.and.itau.le.18)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+        endif
+
+        if (itau.eq.19.and.isplit.le.2)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+
+        endif
+
+c Inj pour le stokage ttes les 4h
+
+        if (itau.ge.5.and.itau.le.7)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+        endif
+
+c Inj pour le stokage ttes les 3h
+
+        if (itau.eq.6.and.isplit.gt.2)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+        endif
+
+        if (itau.ge.7.and.itau.le.9)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+        endif
+
+        if (itau.eq.10.and.isplit.le.2)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+        endif
+
+c  Inj pour le stokage ttes les 6h
+
+        if (itau.eq.3.and.isplit.gt.8)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+        endif
+
+        if (itau.eq.4)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+        endif
+
+        if (itau.eq.5.and.isplit.le.8)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+        endif
+
+c Inj sur les 4 1ere couches pour le stokage ttes les 3h
+
+        if (itau.eq.6.and.isplit.gt.2)then
+	do l=1,4
+        q(i_stat(1),j_stat(1),l,1)=
+     s      q(i_stat(1),j_stat(1),l,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),l)/
+     s   (24.*4.)
+	enddo
+        endif
+
+        if (itau.ge.7.and.itau.le.9)then
+	do l=1,4
+        q(i_stat(1),j_stat(1),l,1)=
+     s      q(i_stat(1),j_stat(1),l,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),l)/
+     s   (24.*4.)
+	enddo
+        endif
+
+        if (itau.eq.10.and.isplit.le.2)then
+        do l=1,4
+        q(i_stat(1),j_stat(1),l,1)=
+     s      q(i_stat(1),j_stat(1),l,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),l)/
+     s   (24.*4.)
+        enddo
+        endif
+
+c Inj pour le stokage ttes les 12h
+
+	if (itau.eq.2.and.isplit.gt.8)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+        endif
+
+        if (itau.eq.3.and.isplit.le.8)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+        endif
+
+c Inj pour le stokage ttes les 24h
+
+        if (itau.eq.1.and.isplit.gt.32)then
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/24.
+        endif
+
+        if (itau.eq.2.and.isplit.le.8)then
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/24.
+        endif     
+
+c Inj pour le stokage ttes les 3h
+
+        if (itau.eq.6.and.isplit.gt.2)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+        endif
+
+        if (itau.ge.7.and.itau.le.9)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+        endif
+
+        if (itau.eq.10.and.isplit.le.2)then
+        do iiinj=1,20
+        q(i_stat(1),j_stat(1),1,1)=
+     s      q(i_stat(1),j_stat(1),1,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),1)/
+     s   (24.*20.)
+        enddo
+        endif
+
+333     continue
+
+        if (itau.eq.6.and.isplit.gt.2)then
+
+        do l=1,2
+        q(i_stat(1),j_stat(1),l,1)=
+     s      q(i_stat(1),j_stat(1),l,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),l)/
+     s   (24.*2.)
+        enddo
+        endif
+
+        if (itau.ge.7.and.itau.le.9)then
+        do l=1,2
+        q(i_stat(1),j_stat(1),l,1)=
+     s      q(i_stat(1),j_stat(1),l,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),l)/
+     s   (24.*2.)
+        enddo
+        endif
+
+        if (itau.eq.10.and.isplit.le.2)then
+        do l=1,2
+        q(i_stat(1),j_stat(1),l,1)=
+     s      q(i_stat(1),j_stat(1),l,1)
+     s  +intensite(1)/masse(i_stat(1),j_stat(1),l)/
+     s   (24.*2.)
+        enddo
+        endif
+
+c=======================================================================
+c   FIN DE LA PARTIE TERMES SOURCES
+c=======================================================================
+c-----------------------------------------------------------------------
+c   Advection
+c   ----------
+	print*,'ENTREE DANS VLSPLT'
+c=======================================================================
+c         do isplit=1,nsplit
+c=======================================================================
+CC MAF dans le cas du offline on ne fait pas la difference masse et massem
+            do iii=1,nsplit_dyn
+		  nom='pbaru'
+                  call wrgrads(3,llm,pbaru(:,:,1),nom,nom)
+                  nom='pbarv'
+                  call wrgrads(4,llm,pbarv(:,:,1),nom,nom)
+                  nom='masse'
+                  call wrgrads(3,llm,masse(:,:,1),nom,nom)
+                  nom='w'
+                  call wrgrads(3,llm,w(:,:,1),nom,nom)
+           print*,'Entree dans vlsplt a itau et au temps',itau,time
+	   print*,'Pas advect :',dtdyn
+           print*,'Avant vlsplt'
+            do iq=1,nq
+               write(str2,'(i2.2)') iq
+              call diagadv(q(:,:,:,iq),masse,'Traceur '//str2)
+            enddo
+               do iq=1,nq
+ctest abder   CALL vlsplt(q(1,1,1,iq),2.,masse,w,pbaru,pbarv,dtdyn)
+              CALL vlsplt(q(:,:,:,iq),2.,masse,w,pbaru,pbarv,dtdyn)
+               enddo
+
+               do iq=iqmin,iqmax
+                  write(str2,'(i2.2)') iq
+                  nom='q'//str2
+                  call wrgrads(5,1,q(:,:,1,iq),nom,nom)
+               enddo
+            print*,'Apres vlsplt'
+            do iq=1,nq
+               write(str2,'(i2.2)') iq
+              call diagadv(q(:,:,:,iq),masse,'Traceur '//str2)
+            enddo
+
+c   mise a jour du champ de masse tenant compte de l'advection sur
+c   le pas de temps considere.
+
+            enddo
+
+c-----------------------------------------------------------------------
+c   Physique  (lessivage dans phytrac maintenant) MAF
+c   ---------------------------------------------------------
+        if(physic) then
+         if(debut) then
+         latfi(1)=rlatu(1)*180./pi
+         lonfi(1)=0.
+         DO j=2,jjm
+            DO i=1,iim
+               latfi((j-2)*iim+1+i)= rlatu(j)*180./pi
+               lonfi((j-2)*iim+1+i)= rlonv(i)*180./pi
+            ENDDO
+         ENDDO
+         latfi(ngridmx)= rlatu(jjp1)*180./pi
+         lonfi(ngridmx)= 0.
+         endif
+c
+C Calcul de paprs et pplay et de temp
+c   -----------------------------------------------------------------
+c     .... paprs  definis aux (llm +1) interfaces des couches  ....
+c     .... pplay  definis aux (  llm )    milieux des couches  ....
+c   -----------------------------------------------------------------
+C
+      DO l = 1, llmp1
+        paprs( 1,l ) = p(1,1,l)
+        ig0 = 2
+          DO j = 2, jjm
+             DO i =1, iim
+              paprs( ig0,l ) = p(i,j,l)
+              ig0 = ig0 +1
+             ENDDO
+          ENDDO
+        paprs( ngridmx,l ) = p(1,jjp1,l)
+      ENDDO
+c
+      unskap   = 1./ kappa
+      DO l=1,llm
+         pksurcp     =  pk(1,1,l) / cpp
+         pplay(1,l)  =  preff * pksurcp ** unskap
+         tempfi(1,l)   =  teta(1,1,l) *  pksurcp
+         ig0         = 2
+         DO j = 2, jjm
+            DO i = 1, iim
+              pksurcp        = pk(i,j,l) / cpp
+              pplay(ig0,l)   = preff * pksurcp ** unskap
+              tempfi(ig0,l)    = teta(i,j,l)  * pksurcp
+              ig0            = ig0 + 1
+            ENDDO
+         ENDDO
+         pksurcp       = pk(1,jjp1,l) / cpp
+         pplay(ig0,l)  = preff * pksurcp ** unskap
+         tempfi (ig0,l)  = teta(1,jjp1,l)  * pksurcp
+      ENDDO
+
+c On passe le traceur et le geopot de surf  sur la grille physique
+             call gr_dyn_fi(llm*nq,iip1,jjp1,ngridmx,q,qfi)
+c             call gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,pphis)
+             do iii=1,nsplit_phy
+C
+      print*,'dtphys avant phytrac ',dtphys
+             call phytrac(rnpb,
+     I                   ecritphy, debutphy,
+     I                   nq,
+     I                   ngridmx,llm,dtphys,
+     I                   t,paprs,pplay,
+     I                   zmfu, zmfd, zen_u, zde_u, zen_d, zde_d,
+     I                   coefkz,yu1,yv1,ftsol,pctsrf,latfi,
+     I                   frac_impa, frac_nucl,
+     I                   lonfi,presnivs,airefi,pphis,
+     O                   qfi)
+C
+             debutphy= .FALSE.
+C
+             enddo ! nsplit_phy
+C
+c   On passe le traceur sur la grille dynamique.
+             call gr_fi_dyn(llm*nq,ngridmx,iip1,jjp1,qfi,q)
+
+             endif
+
+             itauav=(itau-1)*nsplit+isplit
+c             itauav=itau*nsplit+isplit
+             CALL writedynav(histaveid, nq,mode, itauav,vcov ,
+     ,                   ucov,teta,pk,phi,q,masse,ps,phis)
+c            qmoy(:,:,:)=qmoy(:,:,:)+q(:,:,1,:)
+         do iq=1,nq
+            do j=1,jjp1
+               do i=1,iip1
+                  qmoy(i,j,iq)=qmoy(i,j,iq)+q(i,j,1,iq)
+               enddo
+            enddo
+         enddo
+         enddo ! isplit=1,nsplit
+
+         call histsync
+
+c-----------------------------------------------------------------------
+c    Calcul de ucov et vcov pour les sorties
+c   -------------------------------------
+
+         call massbar(masse, massebx, masseby )
+         do l=1,llm
+            do j=1,jjp1
+               do i=1,iip1
+               ucov(i,j,l)=pbaru(i,j,l)/massebx(i,j,l)*cu(i,j)*cu(i,j)
+     s         /istdyn
+               enddo
+            enddo
+
+            do j=1,jjm
+               do i=1,iip1
+               vcov(i,j,l)=pbarv(i,j,l)/masseby(i,j,l)*cv(i,j)*cv(i,j)
+     s         /istdyn
+               enddo
+            enddo
+         enddo
+         iday= dayref+itau/iday_step
+         time=
+     s   float(itau-(iday-dayref)*iday_step)/iday_step+time_0
+         IF(time.GT.1.) THEN
+            time = time-1.
+            iday = iday+1
+         ENDIF
+C
+c=======================================================================
+               do iq=iqmin,iqmax
+                  write(str2,'(i2.2)') iq
+                  nom='q'//str2
+                  call wrgrads(1,1,qmoy(:,:,iq),nom,nom)
+                  call wrgrads(2,1,q(:,:,ldec,iq),nom,nom)
+               enddo
+
+         do iq=1,nq
+            do j=1,jjp1
+               do i=1,iip1
+                  qmoy(i,j,iq)=0.
+               enddo
+            enddo
+         enddo
+
+         print*,'Record=',irec,'   stoke=',stoke,
+     s      '      injecte=',injecte
+
+      enddo ! itau
+
+      CALL dynredem1("restart.nc",0.0,
+     .          vcov,ucov,teta,q,nq,masse,ps)
+
+c=======================================================================
+C fermeture des fichiers netcdf
+          CALL histclo
+C
+ 300  FORMAT('1',/,15x,'run du pas',i7,2x,'au pas',i7,2x,
+     s 'c"est a dire du jour',i7,3x,'au jour',i7,/,/)
+C----------------------------------------------------------------------
+
+      stop
+      end
