! 
! $Header$
!
c
c
#define IO_DEBUG

!#undef CPP_IOIPSL

      SUBROUTINE leapfrog_p(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0,
     &                    time_0)

       USE misc_mod
       USE parallel
       USE times
       USE mod_hallo
       USE Bands
       USE Write_Field
       USE Write_Field_p
       USE vampir

      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
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 pour l'advection de
c        q  , en modifiant iadv dans traceur.def  (10/02) .
c
c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
c      Pour Van-Leer iadv=10 
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 "description.h"
#include "serre.h"
#include "com_io_dyn.h"
#include "iniprint.h"

c#include "tracstoke.h"

#include "academic.h"
#include "clesphys.h"
#include "advtrac.h"
      
      include 'mpif.h'
      integer nq

      INTEGER         longcles
      PARAMETER     ( longcles = 20 )
      REAL  clesphy0( longcles )

      real zqmin,zqmax
      INTEGER nbetatmoy, nbetatdem,nbetat

c   variables dynamiques
      REAL :: vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
      REAL :: teta(ip1jmp1,llm)                 ! temperature potentielle 
      REAL :: q(ip1jmp1,llm,nqmx)               ! champs advectes
      REAL :: ps(ip1jmp1)                       ! pression  au sol
      REAL,SAVE :: p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
      REAL,SAVE :: pks(ip1jmp1)                      ! exner au  sol
      REAL,SAVE :: pk(ip1jmp1,llm)                   ! exner au milieu des couches
      REAL,SAVE :: pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
      REAL :: masse(ip1jmp1,llm)                ! masse d'air
      REAL :: phis(ip1jmp1)                     ! geopotentiel au sol
      REAL,SAVE :: phi(ip1jmp1,llm)                  ! geopotentiel
      REAL,SAVE :: w(ip1jmp1,llm)                    ! vitesse verticale

c variables dynamiques intermediaire pour le transport
      REAL,SAVE :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse

c   variables dynamiques au pas -1
      REAL,SAVE :: vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
      REAL,SAVE :: tetam1(ip1jmp1,llm),psm1(ip1jmp1)
      REAL,SAVE :: massem1(ip1jmp1,llm)

c   tendances dynamiques
      REAL,SAVE :: dv(ip1jm,llm),du(ip1jmp1,llm)
      REAL,SAVE :: dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)

c   tendances de la dissipation
      REAL,SAVE :: dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
      REAL,SAVE :: dtetadis(ip1jmp1,llm)

c   tendances physiques
      REAL,SAVE :: dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
      REAL,SAVE :: dtetafi(ip1jmp1,llm)
      REAL,SAVE :: dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)

c   variables pour le fichier histoire
      REAL dtav      ! intervalle de temps elementaire

      REAL tppn(iim),tpps(iim),tpn,tps
c
      INTEGER itau,itaufinp1,iav
      INTEGER*4  iday ! jour julien
      REAL       time ! Heure de la journee en fraction d'1 jour

      REAL  SSUM
      REAL time_0 
      REAL,SAVE :: finvmaold(ip1jmp1,llm)

cym      LOGICAL  lafin
      LOGICAL :: lafin
      INTEGER ij,iq,l
      INTEGER ik

      real time_step, t_wrt, t_ops

      REAL rdayvrai,rdaym_ini
      LOGICAL first,callinigrads

      data callinigrads/.true./
      character*10 string10

      REAL,SAVE :: alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
#ifdef INCA
      REAL,SAVE :: flxw(ip1jmp1,llm)
#endif

c+jld variables test conservation energie
      REAL,SAVE :: ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
C     Tendance de la temp. potentiel d (theta)/ d t due a la 
C     tansformation d'energie cinetique en energie thermique
C     cree par la dissipation
      REAL,SAVE :: dtetaecdt(ip1jmp1,llm)
      REAL,SAVE :: vcont(ip1jm,llm),ucont(ip1jmp1,llm)
      REAL,SAVE :: vnat(ip1jm,llm),unat(ip1jmp1,llm)
      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
      CHARACTER*15 ztit
!      INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
!      SAVE      ip_ebil_dyn
!      DATA      ip_ebil_dyn/0/
c-jld 

      character*80 dynhist_file, dynhistave_file
      character*20 modname
      character*80 abort_message

C Calendrier
      LOGICAL true_calendar
      PARAMETER (true_calendar = .false.)

      logical,PARAMETER :: dissip_conservative=.TRUE.
 
      INTEGER testita
      PARAMETER (testita = 9)
      
c declaration liees au parallelisme
      INTEGER :: ierr
      LOGICAL :: FirstCaldyn
      LOGICAL :: FirstPhysic
      INTEGER :: ijb,ije,j,i
      type(Request) :: TestRequest
      type(Request) :: Request_Dissip
      type(Request) :: Request_physic
      REAL,SAVE :: dvfi_tmp(iip1,llm),dufi_tmp(iip1,llm)
      REAL,SAVE :: dtetafi_tmp(iip1,llm),dqfi_tmp(iip1,llm,nqmx)
      REAL,SAVE :: dpfi_tmp(iip1)

      INTEGER :: true_itau
      LOGICAL :: verbose=.true.
      INTEGER :: iapptrac
      INTEGER :: AdjustCount
      INTEGER :: var_time
      LOGICAL :: ok_start_timer=.FALSE.

c$OMP MASTER
      ItCount=0
c$OMP END MASTER      
      true_itau=0
      FirstCaldyn=.TRUE.
      FirstPhysic=.TRUE.
      iapptrac=0
      AdjustCount = 0
      lafin=.false.
      
      itaufin   = nday*day_step
      itaufinp1 = itaufin +1


      itau = 0
      iday = day_ini+itau/day_step
      time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
         IF(time.GT.1.) THEN
          time = time-1.
          iday = iday+1
         ENDIF

c$OMP MASTER
      OMP_CHUNK=5
c$OMP END MASTER	  
c-----------------------------------------------------------------------
c   On initialise la pression et la fonction d'Exner :
c   --------------------------------------------------

c$OMP MASTER
      dq=0.
      CALL pression ( ip1jmp1, ap, bp, ps, p       )
      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
c$OMP END MASTER
c-----------------------------------------------------------------------
c   Debut de l'integration temporelle:
c   ----------------------------------
c et du parallelisme !!

   1  CONTINUE

c$OMP MASTER
      call MPI_BARRIER(COMM_LMDZ,ierr)
c$OMP END MASTER
c$OMP BARRIER

#ifdef CPP_IOIPSL
c$OMP MASTER
      if (ok_guide.and.(itaufin-itau-1)*dtvr.gt.21600) then
        call guide_pp(itau,ucov,vcov,teta,q,masse,ps)
      else
        IF(prt_level>9)WRITE(*,*)'attention on ne guide pas les ',
     .    '6 dernieres heures'
      endif
c$OMP END MASTER
#endif
c
c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
c       CALL  test_period ( ucov,vcov,teta,q,p,phis )
c       PRINT *,' ----   Test_period apres continue   OK ! -----', itau
c     ENDIF 
c
cym      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
cym      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
cym      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
cym      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
cym      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )

       if (FirstCaldyn) then
c$OMP MASTER
         ucovm1=ucov
         vcovm1=vcov
         tetam1= teta
         massem1= masse
         psm1= ps
         
         finvmaold = masse
         CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
c$OMP END MASTER
c$OMP BARRIER
       else
          
         ijb=ij_begin
         ije=ij_end

c$OMP MASTER           
         psm1     (ijb:ije) = ps    (ijb:ije)
c$OMP END MASTER

c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
	 DO l=1,llm      
           ije=ij_end
	   ucovm1   (ijb:ije,l) = ucov  (ijb:ije,l)
           tetam1   (ijb:ije,l) = teta  (ijb:ije,l)
           massem1  (ijb:ije,l) = masse (ijb:ije,l)
           finvmaold(ijb:ije,l)=masse(ijb:ije,l)
	         
           if (pole_sud) ije=ij_end-iip1
           vcovm1(ijb:ije,l) = vcov  (ijb:ije,l)
       

         ENDDO
c$OMP ENDDO	  


          CALL filtreg_p ( finvmaold ,jj_begin,jj_end,jjp1, 
     .                    llm, -2,2, .TRUE., 1 )

       endif
       
      forward = .TRUE.
      leapf   = .FALSE.
      dt      =  dtvr

c   ...    P.Le Van .26/04/94  ....

cym      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
cym      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )

cym  ne sert a rien
cym      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)

   2  CONTINUE

c$OMP MASTER
      ItCount=ItCount+1
      if (MOD(ItCount,1)==1) then
        debug=.true.
      else
        debug=.false.
      endif
c$OMP END MASTER
c-----------------------------------------------------------------------

c   date:
c   -----


c   gestion des appels de la physique et des dissipations:
c   ------------------------------------------------------
c
c   ...    P.Le Van  ( 6/02/95 )  ....

      apphys = .FALSE.
      statcl = .FALSE.
      conser = .FALSE.
      apdiss = .FALSE.
c      idissip=1
      IF( purmats ) THEN
         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
         IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward 
     s          .and. iflag_phys.NE.0                 ) apphys = .TRUE.
      ELSE
         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.NE.0) apphys=.TRUE.
      END IF

cym    ---> Pour le moment      
cym      apphys = .FALSE.
      statcl = .FALSE.
      conser = .FALSE.
      
      if (firstCaldyn) then
c$OMP MASTER
          call SetDistrib(jj_Nb_Caldyn)
c$OMP END MASTER
c$OMP BARRIER
          firstCaldyn=.FALSE.
cym          call InitTime
c$OMP MASTER
          call Init_timer
c$OMP END MASTER
      endif

c$OMP MASTER      
      IF (ok_start_timer) THEN
        CALL InitTime
	ok_start_timer=.FALSE.
      ENDIF      
c$OMP END MASTER      
     
      if (Adjust) then
c$OMP MASTER 
        AdjustCount=AdjustCount+1
        if (iapptrac==iapp_tracvl .and. (forward. OR . leapf)
     &         .and. itau/iphysiq>2 .and. Adjustcount>30) then
           AdjustCount=0
           call allgather_timer_average

        if (Verbose) then
        
	print *,'*********************************'
	print *,'******    TIMER CALDYN     ******'
	do i=0,mpi_size-1
	  print *,'proc',i,' :   Nb Bandes  :',jj_nb_caldyn(i),
     &            '  : temps moyen :',
     &             timer_average(jj_nb_caldyn(i),timer_caldyn,i),
     &            '+-',timer_delta(jj_nb_caldyn(i),timer_caldyn,i)
        enddo
      
        print *,'*********************************'
	print *,'******    TIMER VANLEER    ******'
	do i=0,mpi_size-1
	  print *,'proc',i,' :   Nb Bandes  :',jj_nb_vanleer(i),
     &            '  : temps moyen :',
     &             timer_average(jj_nb_vanleer(i),timer_vanleer,i),
     &            '+-',timer_delta(jj_nb_vanleer(i),timer_vanleer,i)
        enddo
      
        print *,'*********************************'
	print *,'******    TIMER DISSIP    ******'
	do i=0,mpi_size-1
	  print *,'proc',i,' :   Nb Bandes  :',jj_nb_dissip(i),
     &            '  : temps moyen :',
     &             timer_average(jj_nb_dissip(i),timer_dissip,i),
     &             '+-',timer_delta(jj_nb_dissip(i),timer_dissip,i)
        enddo
        
	if (mpi_rank==0) call WriteBands
	
       endif
       
	 call AdjustBands_caldyn
         if (mpi_rank==0) call WriteBands
	 
	 call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(ucovm1,ucovm1,ip1jmp1,llm,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(vcovm1,vcovm1,ip1jm,llm,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(tetam1,tetam1,ip1jmp1,llm,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(massem1,massem1,ip1jmp1,llm,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(ps,ps,ip1jmp1,1,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(psm1,psm1,ip1jmp1,1,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(pkf,pkf,ip1jmp1,llm,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(pks,pks,ip1jmp1,1,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
     &                                jj_Nb_caldyn,0,0,TestRequest)
         call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
     &                                jj_Nb_caldyn,0,0,TestRequest)
 
	do j=1,nqmx
         call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
     &                                jj_nb_caldyn,0,0,TestRequest)
        enddo

         call SetDistrib(jj_nb_caldyn)
         call SendRequest(TestRequest)
         call WaitRequest(TestRequest)
         
	call AdjustBands_dissip
	call AdjustBands_physic

      endif
c$OMP END MASTER  
      endif       
     
      
      
c-----------------------------------------------------------------------
c   calcul des tendances dynamiques:
c   --------------------------------
c$OMP BARRIER
c$OMP MASTER
       call VTb(VThallo)
       call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,TestRequest)
       call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,TestRequest)
       call Register_Hallo(teta,ip1jmp1,llm,1,1,1,1,TestRequest)
       call Register_Hallo(ps,ip1jmp1,1,1,2,2,1,TestRequest)
       call Register_Hallo(pkf,ip1jmp1,llm,1,1,1,1,TestRequest)
       call Register_Hallo(pk,ip1jmp1,llm,1,1,1,1,TestRequest)
       call Register_Hallo(pks,ip1jmp1,1,1,1,1,1,TestRequest)
       call Register_Hallo(p,ip1jmp1,llmp1,1,1,1,1,TestRequest)
       
c       do j=1,nqmx
c         call Register_Hallo(q(1,1,j),ip1jmp1,llm,1,1,1,1,
c     *                       TestRequest)
c        enddo

       call SendRequest(TestRequest)
       call WaitRequest(TestRequest)
       call VTe(VThallo)

      
      if (debug) then	
        call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
        call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
        call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
        call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
        call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
	call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
        call WriteField_p('pks',reshape(pks,(/iip1,jmp1/)))
        call WriteField_p('pkf',reshape(pkf,(/iip1,jmp1,llm/)))
        call WriteField_p('phis',reshape(phis,(/iip1,jmp1/)))
        do j=1,nqmx
          call WriteField_p('q'//trim(int2str(j)),
     .                reshape(q(:,:,j),(/iip1,jmp1,llm/)))
        enddo        
      endif
c$OMP END MASTER 
c$OMP BARRIER
      
      True_itau=True_itau+1

c$OMP MASTER
      print*,"Iteration No",True_itau


      call start_timer(timer_caldyn)

      CALL geopot_p  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )

      
      call VTb(VTcaldyn)
c$OMP END MASTER
c$OMP BARRIER
      var_time=time+iday-day_ini
cc$OMP PARALLEL DEFAULT(SHARED)
cc$OMP+         PRIVATE(rdaym_ini,rdayvrai,ijb,ije,
cc$OMP+                 tppn,tpn,tpps,tps)

cc$OMP+         SHARED(itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
cc$OMP+                phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, 
cc$OMP+                time, iday,day_ini,forward,leapf, iapptrac,
cc$OMP+                q,dq,p,VTcaldyn,offline,dtvr,itau)      

      CALL caldyn_p 
     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )

ccc$OMP END PARALLEL      
c$OMP MASTER
      call VTe(VTcaldyn)
c$OMP END MASTER      
c      call WriteField_p('du',reshape(du,(/iip1,jmp1,llm/)))
c      call WriteField_p('dv',reshape(dv,(/iip1,jjm,llm/)))
c      call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
c      call WriteField_p('dp',reshape(dp,(/iip1,jmp1/)))
c      call WriteField_p('w',reshape(w,(/iip1,jmp1,llm/)))
c      call WriteField_p('pbaru',reshape(pbaru,(/iip1,jmp1,llm/)))
c      call WriteField_p('pbarv',reshape(pbarv,(/iip1,jjm,llm/)))

c-----------------------------------------------------------------------
c   calcul des tendances advection des traceurs (dont l'humidite)
c   -------------------------------------------------------------

      IF( forward. OR . leapf )  THEN
cc$OMP PARALLEL DEFAULT(SHARED) 
c
#ifdef INCA
             CALL caladvtrac_p(q,pbaru,pbarv,
     *                      p, masse, dq,  teta,
     .             flxw,
     .             pk, 
     .             iapptrac)
#else
             CALL caladvtrac_p(q,pbaru,pbarv,
     *                      p, masse, dq,  teta,
     .             pk,iapptrac)
#endif



c      do j=1,nqmx
c        call WriteField_p('q'//trim(int2str(j)),
c     .                reshape(q(:,:,j),(/iip1,jmp1,llm/)))
c        call WriteField_p('dq'//trim(int2str(j)),
c     .                reshape(dq(:,:,j),(/iip1,jmp1,llm/)))
c      enddo
       IF (offline) THEN
Cmaf stokage du flux de masse pour traceurs OFF-LINE
#undef CPP_IOIPSL
#ifdef CPP_IOIPSL
           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
     .   dtvr, itau)
#endif


         ENDIF
c
      ENDIF
cc$OMP END PARALLEL

c-----------------------------------------------------------------------
c   integrations dynamique et traceurs:
c   ----------------------------------

c$OMP MASTER 
       call VTb(VTintegre)
c$OMP END MASTER
c      call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
c      call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
c      call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
c      call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
c      call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
c      call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
c      call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
c      call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
cc$OMP PARALLEL DEFAULT(SHARED)
       CALL integrd_p ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
     $              finvmaold                                    )


c      call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
c      call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
c      call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
c      call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
c      call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
c      call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
c      call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
c      call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))

c      call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
c$OMP MASTER 
       call VTe(VTintegre)
c$OMP END MASTER
c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
c
c-----------------------------------------------------------------------
c   calcul des tendances physiques:
c   -------------------------------
c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
c
       IF( purmats )  THEN
          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
       ELSE
          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
       ENDIF

cc$OMP END PARALLEL

c
c
       IF( apphys )  THEN
c
c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
c
cc$OMP PARALLEL DEFAULT(SHARED)
cc$OMP+         PRIVATE(rdaym_ini,rdayvrai,ijb,ije)

c$OMP MASTER
         call suspend_timer(timer_caldyn)
         print*,'Entree dans la physique : Iteration No ',true_itau
c$OMP END MASTER

c$OMP BARRIER
         CALL pression_p (  ip1jmp1, ap, bp, ps,  p      )

         CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
           rdaym_ini  = itau * dtvr / daysec
           rdayvrai   = rdaym_ini  + day_ini


c rajout debug
c       lafin = .true.


c   Inbterface avec les routines de phylmd (phymars ... )
c   -----------------------------------------------------

#ifdef CPP_PHYS
c+jld

c  Diagnostique de conservation de l'energie : initialisation
      IF (ip_ebil_dyn.ge.1 ) THEN 
          ztit='bil dyn'
          CALL diagedyn(ztit,2,1,1,dtphys
     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
      ENDIF 
c-jld
c$OMP BARRIER
c$OMP MASTER
        call VTb(VThallo)
        call SetTag(Request_physic,800)
	
        call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
     *                               jj_Nb_physic,2,2,Request_physic)
        
        call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
     *                               jj_Nb_physic,2,2,Request_physic)
        
        call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
     *                               jj_Nb_physic,2,2,Request_physic)
        
        call Register_SwapFieldHallo(p,p,ip1jmp1,llmp1,
     *                               jj_Nb_physic,2,2,Request_physic)
        
        call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
     *                               jj_Nb_physic,2,2,Request_physic)
        
        call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
     *                               jj_Nb_physic,2,2,Request_physic)
        
        call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
     *                               jj_Nb_physic,2,2,Request_physic)
        
        call Register_SwapFieldHallo(w,w,ip1jmp1,llm,
     *                               jj_Nb_physic,2,2,Request_physic)
        
c	call SetDistrib(jj_nb_vanleer)
        do j=1,nqmx
 
          call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
     *                               jj_Nb_physic,2,2,Request_physic)
        enddo
#ifdef INCA
        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
     *                               jj_Nb_physic,2,2,Request_physic)
#endif
        call SetDistrib(jj_nb_Physic)
	
        call SendRequest(Request_Physic)
        call WaitRequest(Request_Physic)       
	
	call VTe(VThallo)
	
	call VTb(VTphysiq)
c$OMP END MASTER
c$OMP BARRIER

cc$OMP MASTER	
c      call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
c      call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
c      call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
c      call WriteField_p('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
c      call WriteField_p('pkfi',reshape(pk,(/iip1,jmp1,llm/)))
cc$OMP END MASTER
cc$OMP BARRIER
	
        CALL calfis_p( nq, lafin ,rdayvrai,time  ,
     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
     $               du,dv,dteta,dq,w,
#ifdef INCA
     $               flxw,
#endif
     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
        ijb=ij_begin
        ije=ij_end  
        if ( .not. pole_nord) then
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
          DO l=1,llm
          dufi_tmp(1:iip1,l)   = dufi(ijb:ijb+iim,l) 
          dvfi_tmp(1:iip1,l)   = dvfi(ijb:ijb+iim,l)  
          dtetafi_tmp(1:iip1,l)= dtetafi(ijb:ijb+iim,l)  
          dqfi_tmp(1:iip1,l,:) = dqfi(ijb:ijb+iim,l,:)  
          ENDDO
c$OMP END DO NOWAIT

c$OMP MASTER
          dpfi_tmp(1:iip1)     = dpfi(ijb:ijb+iim)  
c$OMP END MASTER
	endif

c$OMP BARRIER
c$OMP MASTER
	call SetDistrib(jj_nb_Physic_bis)

        call VTb(VThallo)
 
        call Register_Hallo(dufi,ip1jmp1,llm,
     *                      1,0,0,1,Request_physic)
        
        call Register_Hallo(dvfi,ip1jm,llm,
     *                      1,0,0,1,Request_physic)
        
        call Register_Hallo(dtetafi,ip1jmp1,llm,
     *                      1,0,0,1,Request_physic)

        call Register_Hallo(dpfi,ip1jmp1,1,
     *                      1,0,0,1,Request_physic)

        do j=1,nqmx
          call Register_Hallo(dqfi(1,1,j),ip1jmp1,llm,
     *                        1,0,0,1,Request_physic)
        enddo
        
        call SendRequest(Request_Physic)
        call WaitRequest(Request_Physic)
             
        call VTe(VThallo)
 
        call SetDistrib(jj_nb_Physic)
c$OMP END MASTER
c$OMP BARRIER	
	        ijb=ij_begin
	if (.not. pole_nord) then
	
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
          DO l=1,llm
            dufi(ijb:ijb+iim,l) = dufi(ijb:ijb+iim,l)+dufi_tmp(1:iip1,l)
            dvfi(ijb:ijb+iim,l) = dvfi(ijb:ijb+iim,l)+dvfi_tmp(1:iip1,l) 
            dtetafi(ijb:ijb+iim,l) = dtetafi(ijb:ijb+iim,l)
     &                              +dtetafi_tmp(1:iip1,l)
            dqfi(ijb:ijb+iim,l,:) = dqfi(ijb:ijb+iim,l,:)
     &                              + dqfi_tmp(1:iip1,l,:)
          ENDDO
c$OMP END DO NOWAIT

c$OMP MASTER
          dpfi(ijb:ijb+iim)   = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
c$OMP END MASTER
	  
        endif
c$OMP BARRIER
cc$OMP MASTER	
c      call WriteField_p('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
c      call WriteField_p('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
c      call WriteField_p('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
c      call WriteField_p('dpfi',reshape(dpfi,(/iip1,jmp1/)))
cc$OMP END MASTER
c      
c      do j=1,nqmx
c        call WriteField_p('dqfi'//trim(int2str(j)),
c     .                reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
c      enddo

c      ajout des tendances physiques:
c      ------------------------------
          CALL addfi_p( nqmx, dtphys, leapf, forward   ,
     $                  ucov, vcov, teta , q   ,ps ,
     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )

c$OMP BARRIER
c$OMP MASTER
        call VTe(VTphysiq)

        call VTb(VThallo)

        call SetTag(Request_physic,800)
        call Register_SwapField(ucov,ucov,ip1jmp1,llm,
     *                               jj_Nb_caldyn,Request_physic)
        
        call Register_SwapField(vcov,vcov,ip1jm,llm,
     *                               jj_Nb_caldyn,Request_physic)
        
        call Register_SwapField(teta,teta,ip1jmp1,llm,
     *                               jj_Nb_caldyn,Request_physic)
        
        call Register_SwapField(p,p,ip1jmp1,llmp1,
     *                               jj_Nb_caldyn,Request_physic)
        
        call Register_SwapField(pk,pk,ip1jmp1,llm,
     *                               jj_Nb_caldyn,Request_physic)
        
        call Register_SwapField(phis,phis,ip1jmp1,1,
     *                               jj_Nb_caldyn,Request_physic)
        
        call Register_SwapField(phi,phi,ip1jmp1,llm,
     *                               jj_Nb_caldyn,Request_physic)
        
        call Register_SwapField(w,w,ip1jmp1,llm,
     *                               jj_Nb_caldyn,Request_physic)

        do j=1,nqmx
        
          call Register_SwapField(q(1,1,j),q(1,1,j),ip1jmp1,llm,
     *                               jj_Nb_caldyn,Request_physic)
        
        enddo

        call SendRequest(Request_Physic)
        call WaitRequest(Request_Physic)     

        call VTe(VThallo)

        call SetDistrib(jj_Nb_caldyn)
c$OMP END MASTER
c$OMP BARRIER
c
c  Diagnostique de conservation de l'energie : difference
      IF (ip_ebil_dyn.ge.1 ) THEN 
          ztit='bil phys'
          CALL diagedyn(ztit,2,1,1,dtphys
     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
      ENDIF 

cc$OMP MASTER      
c      if (debug) then
c       call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
c       call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
c       call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
c      endif
cc$OMP END MASTER

#else

c   Calcul academique de la physique = Rappel Newtonien + fritcion 
c   --------------------------------------------------------------
cym       teta(:,:)=teta(:,:)
cym     s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
       ijb=ij_begin
       ije=ij_end
       teta(ijb:ije,:)=teta(ijb:ije,:)
     s  -iphysiq*dtvr*(teta(ijb:ije,:)-tetarappel(ijb:ije,:))/taurappel

       call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
       call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Physic)
       call SendRequest(Request_Physic)
       call WaitRequest(Request_Physic)     

       call friction_p(ucov,vcov,iphysiq*dtvr)

#endif

c-jld
c$OMP MASTER
         call resume_timer(timer_caldyn)
	 if (FirstPhysic) then
	   ok_start_timer=.TRUE.
	   FirstPhysic=.false.
	 endif
c$OMP END MASTER
       ENDIF

c$OMP BARRIER
        CALL pression_p ( ip1jmp1, ap, bp, ps, p                  )
        CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
c$OMP BARRIER

cc$OMP END PARALLEL

c-----------------------------------------------------------------------
c   dissipation horizontale et verticale  des petites echelles:
c   ----------------------------------------------------------

      IF(apdiss) THEN
cc$OMP  PARALLEL DEFAULT(SHARED) 
cc$OMP+          PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
c$OMP MASTER
        call suspend_timer(timer_caldyn)
	
        print*,'Entree dans la dissipation : Iteration No ',true_itau
c   calcul de l'energie cinetique avant dissipation
        print *,'Passage dans la dissipation'

        call VTb(VThallo)
c$OMP END MASTER

c$OMP BARRIER
c$OMP MASTER
        call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
     *                          jj_Nb_dissip,1,1,Request_dissip)

        call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
     *                          jj_Nb_dissip,1,1,Request_dissip)

        call Register_SwapField(teta,teta,ip1jmp1,llm,
     *                          jj_Nb_dissip,Request_dissip)

        call Register_SwapField(p,p,ip1jmp1,llmp1,
     *                          jj_Nb_dissip,Request_dissip)

        call Register_SwapField(pk,pk,ip1jmp1,llm,
     *                          jj_Nb_dissip,Request_dissip)

        call SendRequest(Request_dissip)       
        call WaitRequest(Request_dissip)       
        call SetDistrib(jj_Nb_dissip)
        
	call VTe(VThallo)

        call VTb(VTdissipation)
        
	call start_timer(timer_dissip)
c$OMP END MASTER
c$OMP BARRIER

        call covcont_p(llm,ucov,vcov,ucont,vcont)
        call enercin_p(vcov,ucov,vcont,ucont,ecin0)

c   dissipation

        CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
         
	ijb=ij_begin
	ije=ij_end
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)	
	DO l=1,llm
	  ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
        ENDDO
c$OMP END DO NOWAIT	
	if (pole_sud) ije=ije-iip1
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)	
	DO l=1,llm
	  vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
        ENDDO
c$OMP END DO NOWAIT	

c       teta=teta+dtetadis


c------------------------------------------------------------------------
        if (dissip_conservative) then
C       On rajoute la tendance due a la transform. Ec -> E therm. cree
C       lors de la dissipation
c$OMP BARRIER
c$OMP MASTER
            call suspend_timer(timer_dissip)
	    call VTb(VThallo)

            call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Request_Dissip)
            call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Dissip)
            call SendRequest(Request_Dissip)
            call WaitRequest(Request_Dissip)
            call VTe(VThallo)
            call resume_timer(timer_dissip)
c$OMP END MASTER
c$OMP BARRIER	    
            call covcont_p(llm,ucov,vcov,ucont,vcont)
            call enercin_p(vcov,ucov,vcont,ucont,ecin)
            
	    ijb=ij_begin
	    ije=ij_end
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)	    
	    do l=1,llm
	      do ij=ijb,ije
	        dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
		dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
              enddo
	    enddo
c$OMP END DO NOWAIT	    
       endif

       ijb=ij_begin
       ije=ij_end
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)	    
	 do l=1,llm
	   do ij=ijb,ije
              teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
	   enddo
	 enddo
c$OMP END DO NOWAIT	 
c------------------------------------------------------------------------


c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
c

        ijb=ij_begin
        ije=ij_end
         
        if (pole_nord) then
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
          DO l  =  1, llm
            DO ij =  1,iim
             tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
            ENDDO
             tpn  = SSUM(iim,tppn,1)/apoln

            DO ij = 1, iip1
             teta(  ij    ,l) = tpn
            ENDDO
          ENDDO
c$OMP END DO NOWAIT

c$OMP MASTER               
          DO ij =  1,iim
            tppn(ij)  = aire(  ij    ) * ps (  ij    )
          ENDDO
            tpn  = SSUM(iim,tppn,1)/apoln
  
          DO ij = 1, iip1
            ps(  ij    ) = tpn
          ENDDO
c$OMP END MASTER
        endif
        
        if (pole_sud) then
c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
          DO l  =  1, llm
            DO ij =  1,iim
             tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
            ENDDO
             tps  = SSUM(iim,tpps,1)/apols

            DO ij = 1, iip1
             teta(ij+ip1jm,l) = tps
            ENDDO
          ENDDO
c$OMP END DO NOWAIT

c$OMP MASTER               
          DO ij =  1,iim
            tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
          ENDDO
            tps  = SSUM(iim,tpps,1)/apols
  
          DO ij = 1, iip1
            ps(ij+ip1jm) = tps
          ENDDO
c$OMP END MASTER
        endif


c$OMP BARRIER
c$OMP MASTER
        call VTe(VTdissipation)

        call stop_timer(timer_dissip)
	
        call VTb(VThallo)

        call Register_SwapField(ucov,ucov,ip1jmp1,llm,
     *                          jj_Nb_caldyn,Request_dissip)

        call Register_SwapField(vcov,vcov,ip1jm,llm,
     *                          jj_Nb_caldyn,Request_dissip)

        call Register_SwapField(teta,teta,ip1jmp1,llm,
     *                          jj_Nb_caldyn,Request_dissip)

        call Register_SwapField(p,p,ip1jmp1,llmp1,
     *                          jj_Nb_caldyn,Request_dissip)

        call Register_SwapField(pk,pk,ip1jmp1,llm,
     *                          jj_Nb_caldyn,Request_dissip)

        call SendRequest(Request_dissip)       
        call WaitRequest(Request_dissip)       
        call SetDistrib(jj_Nb_caldyn)
        call VTe(VThallo)
        call resume_timer(timer_caldyn)
	print *,'fin dissipation'
c$OMP END MASTER
      END IF

cc$OMP END PARALLEL

c ajout debug
c              IF( lafin ) then  
c                abort_message = 'Simulation finished'
c                call abort_gcm(modname,abort_message,0)
c              ENDIF
        
c   ********************************************************************
c   ********************************************************************
c   .... fin de l'integration dynamique  et physique pour le pas itau ..
c   ********************************************************************
c   ********************************************************************

c   preparation du pas d'integration suivant  ......
cym      call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
cym      call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
c$OMP MASTER      
      call stop_timer(timer_caldyn)
c$OMP END MASTER
      IF (itau==itaumax) then
c$OMP MASTER
            call allgather_timer_average

      if (mpi_rank==0) then
        
        print *,'*********************************'
	print *,'******    TIMER CALDYN     ******'
	do i=0,mpi_size-1
	  print *,'proc',i,' :   Nb Bandes  :',jj_nb_caldyn(i),
     &            '  : temps moyen :',
     &             timer_average(jj_nb_caldyn(i),timer_caldyn,i)
        enddo
      
        print *,'*********************************'
	print *,'******    TIMER VANLEER    ******'
	do i=0,mpi_size-1
	  print *,'proc',i,' :   Nb Bandes  :',jj_nb_vanleer(i),
     &            '  : temps moyen :',
     &             timer_average(jj_nb_vanleer(i),timer_vanleer,i)
        enddo
      
        print *,'*********************************'
	print *,'******    TIMER DISSIP    ******'
	do i=0,mpi_size-1
	  print *,'proc',i,' :   Nb Bandes  :',jj_nb_dissip(i),
     &            '  : temps moyen :',
     &             timer_average(jj_nb_dissip(i),timer_dissip,i)
        enddo
	
	print *,'*********************************'
	print *,'******    TIMER PHYSIC    ******'
	do i=0,mpi_size-1
	  print *,'proc',i,' :   Nb Bandes  :',jj_nb_physic(i),
     &            '  : temps moyen :',
     &             timer_average(jj_nb_physic(i),timer_physic,i)
        enddo
	
      endif  
      
      print *,'Taille du Buffer MPI (REAL*8)',MaxBufferSize
      print *,'Taille du Buffer MPI utilise (REAL*8)',MaxBufferSize_Used
      print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
      print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()

        call finalize_parallel
c$OMP END MASTER
	RETURN
      ENDIF
      
      IF ( .NOT.purmats ) THEN
c       ........................................................
c       ..............  schema matsuno + leapfrog  ..............
c       ........................................................

            IF(forward. OR. leapf) THEN
              itau= itau + 1
              iday= day_ini+itau/day_step
              time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
                IF(time.GT.1.) THEN
                  time = time-1.
                  iday = iday+1
                ENDIF
            ENDIF


            IF( itau. EQ. itaufinp1 ) then  

c$OMP MASTER
              call finalize_parallel
c$OMP END MASTER
              abort_message = 'Simulation finished'
              call abort_gcm(modname,abort_message,0)
            ENDIF
c-----------------------------------------------------------------------
c   ecriture du fichier histoire moyenne:
c   -------------------------------------

            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
c$OMP BARRIER
c$OMP MASTER
               IF(itau.EQ.itaufin) THEN
                  iav=1
               ELSE
                  iav=0
               ENDIF
#ifdef CPP_IOIPSL
             call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
              call SendRequest(TestRequest)
              call WaitRequest(TestRequest)

              CALL writedynav_p(histaveid, nqmx, itau,vcov ,
     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
c               call bilan_dyn_p(2,dtvr*iperiod,dtvr*day_step*periodav,
c     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
#endif
c$OMP END MASTER
            ENDIF

c-----------------------------------------------------------------------
c   ecriture de la bande histoire:
c   ------------------------------

c      IF( MOD(itau,iecri         ).EQ.0) THEN

           IF( MOD(itau,iecri*day_step).EQ.0) THEN
c$OMP BARRIER
c$OMP MASTER
               nbetat = nbetatdem
        CALL geopot_p  ( ip1jmp1, teta  , pk , pks,  phis  , phi )
       
cym	unat=0.
	
	ijb=ij_begin
	ije=ij_end
	
	if (pole_nord) then
	  ijb=ij_begin+iip1
	  unat(1:iip1,:)=0.
	endif
	
	if (pole_sud) then 
	  ije=ij_end-iip1
	  unat(ij_end-iip1+1:ij_end,:)=0.
	endif
	    
        do l=1,llm
           unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
        enddo

	ijb=ij_begin
	ije=ij_end
	if (pole_sud) ije=ij_end-iip1
	
	do l=1,llm
           vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
        enddo
	
#ifdef CPP_IOIPSL
 
        CALL writehist_p(histid,histvid, nqmx,itau,vcov, 
     s                       ucov,teta,phi,q,masse,ps,phis)
c#else
c       call Gather_Field(unat,ip1jmp1,llm,0)
c       call Gather_Field(vnat,ip1jm,llm,0)
c       call Gather_Field(teta,ip1jmp1,llm,0)
c       call Gather_Field(ps,ip1jmp1,1,0)
c       do iq=1,nqmx
c        call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
c       enddo
c      
c       if (mpi_rank==0) then
c#include "write_grads_dyn.h"
c       endif
#endif
c$OMP END MASTER
           ENDIF

            IF(itau.EQ.itaufin) THEN

c$OMP BARRIER
c$OMP MASTER

c#ifdef CPP_IOIPSL

       CALL dynredem1_p("restart.nc",0.0,
     ,                     vcov,ucov,teta,q,nqmx,masse,ps)
c#endif

              CLOSE(99)
c$OMP END MASTER
            ENDIF

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

            IF( MOD(itau,iperiod).EQ.0 )    THEN
                    GO TO 1
            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN

                   IF( forward )  THEN
c      fin du pas forward et debut du pas backward

                      forward = .FALSE.
                        leapf = .FALSE.
                           GO TO 2

                   ELSE
c      fin du pas backward et debut du premier pas leapfrog

                        leapf =  .TRUE.
                        dt  =  2.*dtvr
                        GO TO 2
                   END IF
            ELSE

c      ......   pas leapfrog  .....

                 leapf = .TRUE.
                 dt  = 2.*dtvr
                 GO TO 2
            END IF

      ELSE

c       ........................................................
c       ..............       schema  matsuno        ...............
c       ........................................................
            IF( forward )  THEN

             itau =  itau + 1
             iday = day_ini+itau/day_step
             time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0

                  IF(time.GT.1.) THEN
                   time = time-1.
                   iday = iday+1
                  ENDIF

               forward =  .FALSE.
               IF( itau. EQ. itaufinp1 ) then  
c$OMP MASTER
	         call finalize_parallel
c$OMP END MASTER
                 abort_message = 'Simulation finished'
                 call abort_gcm(modname,abort_message,0)
               ENDIF
               GO TO 2

            ELSE

            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
               IF(itau.EQ.itaufin) THEN
                  iav=1
               ELSE
                  iav=0
               ENDIF
#ifdef CPP_IOIPSL
c$OMP BARRIER
c$OMP MASTER

              call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
              call SendRequest(TestRequest)
              call WaitRequest(TestRequest)

              CALL writedynav_p(histaveid, nqmx, itau,vcov ,
     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
c               call bilan_dyn_p (2,dtvr*iperiod,dtvr*day_step*periodav,
c     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
c$OMP END MASTER
#endif
            ENDIF

c               IF(MOD(itau,iecri         ).EQ.0) THEN
              IF(MOD(itau,iecri*day_step).EQ.0) THEN
c$OMP BARRIER
c$OMP MASTER
                  nbetat = nbetatdem
       CALL geopot_p( ip1jmp1, teta  , pk , pks,  phis  , phi   )

cym        unat=0.
	ijb=ij_begin
	ije=ij_end
	
	if (pole_nord) then
	  ijb=ij_begin+iip1
	  unat(1:iip1,:)=0.
	endif
	
	if (pole_sud) then 
	  ije=ij_end-iip1
	  unat(ij_end-iip1+1:ij_end,:)=0.
	endif
	    
        do l=1,llm
           unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
        enddo

	ijb=ij_begin
	ije=ij_end
	if (pole_sud) ije=ij_end-iip1
	
	do l=1,llm
           vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
        enddo

#ifdef CPP_IOIPSL

       CALL writehist_p( histid, histvid, nqmx, itau,vcov , 
     ,                           ucov,teta,phi,q,masse,ps,phis)
c#else
c      call Gather_Field(unat,ip1jmp1,llm,0)
c      call Gather_Field(vnat,ip1jm,llm,0)
c      call Gather_Field(teta,ip1jmp1,llm,0)
c      call Gather_Field(ps,ip1jmp1,1,0)
c      do iq=1,nqmx
c        call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
c      enddo
c      
c      if (mpi_rank==0) then
c#include "write_grads_dyn.h"
c      endif
#endif

c$OMP END MASTER
               ENDIF

                 IF(itau.EQ.itaufin) THEN
c$OMP MASTER
                   CALL dynredem1_p("restart.nc",0.0,
     .                               vcov,ucov,teta,q,nqmx,masse,ps)
c$OMP END MASTER
                 ENDIF
                 forward = .TRUE.
                 GO TO  1

            ENDIF

      END IF
c$OMP MASTER
        call finalize_parallel
c$OMP END MASTER
	RETURN
      END
