Index: /trunk/LMDZ.VENUS/deftank/run.def
===================================================================
--- /trunk/LMDZ.VENUS/deftank/run.def	(revision 1300)
+++ /trunk/LMDZ.VENUS/deftank/run.def	(revision 1301)
@@ -6,7 +6,7 @@
 ## Calendrier
 #calend=venus
-## Jour de l'etat initial ( = 350  si 20 Decembre ,par expl. ,comme ici )
+## Jour de l etat initial ( = 350  si 20 Decembre ,par expl. ,comme ici )
 dayref=1
-##  Annee de l'etat  initial (   avec  4  chiffres   )
+##  Annee de l etat  initial (   avec  4  chiffres   )
 anneeref=1111
 ## Remise a zero de la date initiale
@@ -14,10 +14,15 @@
 ## Reinit des variables de controle
 resetvarc=n
-## Nombre de jours d'integration
+## Nombre de jours d integration
 nday=2
-## Si on veut moins d'un jour :
+## Si on veut moins d un jour :
 less1day=n
 ## si less1day T, fraction de jour du run :
 fractday=0.01
+## si less1day T et que les runs sont enchaines, 
+## indiquer ici la fraction de jour de depart de l etat initial precedent
+## Il faut mettre dans ce cas raz_date=1
+## Defaut=0.
+#starttime=0. 
 ## periode de sortie des variables de controle (en pas)
 iconser=500
@@ -25,5 +30,5 @@
 ## sorties instantanees dans la dynamique (fichiers dyn_hist.nc and co.)
 ok_dyn_ins=n
-## periode d'ecriture des sorties instantanees dans la dynamique
+## periode d ecriture des sorties instantanees dans la dynamique
 ## (en pas dynamiques)
 iecri=500 
Index: /trunk/LMDZ.VENUS/libf/phyvenus/ajsec.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/ajsec.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/ajsec.F	(revision 1301)
@@ -147,5 +147,5 @@
               ENDIF
 c ----------------------------
-c TEST --- PAS DE MELANGE DE U
+c TEST --- PAS DE MELANGE DE U ni Q
 c             zalpha=0.
 c ----------------------------
Index: /trunk/LMDZ.VENUS/libf/phyvenus/clmain.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/clmain.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/clmain.F	(revision 1301)
@@ -368,15 +368,4 @@
          ycoefm(1:knon,2:klev)=ykmm(1:knon,2:klev)
          ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
-c --
-c VENUS: on met le coefficient K =0 au-dessus de 40 km (2.5e5 Pa)
-         do i=1,knon
-         do k=2,klev
-           if (ypaprs(i,k).lt.2.5e5) then
-             ycoefm(i,k)=1.e-7
-             ycoefh(i,k)=1.e-7
-           endif
-         enddo
-         enddo
-c --
 
 c-------------------------------------------------
Index: /trunk/LMDZ.VENUS/libf/phyvenus/comcstVE.h
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/comcstVE.h	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/comcstVE.h	(revision 1301)
@@ -2,15 +2,16 @@
 ! INCLUDE comcstVE.h
 
-        integer nnuve,nbztopve,nbpsve,nbmat
-        parameter (nnuve=68)   ! fichiers Vincent et Bullock
-        parameter (nbpsve=16)   ! number of psurf in Vincent's matrixes
-        parameter (nbztopve=4) ! number of ztop  in Vincent's matrixes
-        parameter (nbmat=nbztopve*nbpsve) ! number of matrixes in Vincent's file
+        integer nnuve,nbmat
+        parameter (nnuve=68)  ! fichiers Vincent et Bullock
+        parameter (nbmat=210) ! Max number of matrixes in Vincent's file
 
-      common/comcstVE/al,bl,ztopve,psurfve
+      common/comcstVE/al,bl,nlatve,indexve,nbpsve,nbszave,               &
+     & psurfve,szave
 
-      real   al(nnuve),bl(nnuve)      ! for Planck luminances calculations
-      real   ztopve(nbmat)   ! cloud top altitude in matrixes (km)
-      real   psurfve(nbmat)            ! surface pressure in matrixes (Pa)
+      real   al(nnuve),bl(nnuve)     ! for Planck luminances calculations
+! Structure of the ksi matrixes
+      integer nlatve,indexve(5),nbpsve(5),nbszave(5)
+      real   psurfve(16,5)           ! surface pressure in matrixes (Pa)
+      real   szave(8,5)              ! converted in mu0
 
 
Index: /trunk/LMDZ.VENUS/libf/phyvenus/drag_noro.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/drag_noro.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/drag_noro.F	(revision 1301)
@@ -30,5 +30,5 @@
 c paprs---input-R-Pressure in semi layers    (Pa)
 c pplay---input-R-Pressure model-layers      (Pa)
-c pgeop---input-R-Geopotential model layers      (m)
+c pgeop---input-R-Geopotential model layers (reference to ground)
 c pn2-----input-R-Brunt-Vaisala freq.^2 at 1/2 layers
 c t-------input-R-temperature (K)
@@ -66,5 +66,5 @@
 c ================
 c
-c zgeom-----R: Altitude of layer above ground
+c zgeom-----R: Altitude (m) of layer above ground (from top to bottom)
 c pt, pu, pv --R: t u v from top to bottom
 c pdtdt, pdudt, pdvdt --R: t u v tendencies (from top to bottom) 
@@ -137,5 +137,5 @@
       DO k = klev, 1, -1
       DO i = 1, klon
-         zgeom(i,k) = pgeop(i,klev-k+1)
+         zgeom(i,k) = pgeop(i,klev-k+1)/RG
 	 zn2(i,k)   = pn2(i,klev-k+1)
       ENDDO
Index: /trunk/LMDZ.VENUS/libf/phyvenus/grid_noro.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/grid_noro.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/grid_noro.F	(revision 1301)
@@ -93,7 +93,11 @@
       REAL zpic(imar+1,jmar),zval(imar+1,jmar)
       real num_tot(2200,1100),num_lan(2200,1100)
-c
+
       REAL a(2200),b(2200),c(1100),d(1100)
-c
+
+c pas defini puisque pas de physique dans newstart...
+      RPI=2.*ASIN(1.)
+      RA=6051300.
+
       print *,' parametres de l orographie a l echelle sous maille' 
 
Index: /trunk/LMDZ.VENUS/libf/phyvenus/ini_histmth.h
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/ini_histmth.h	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/ini_histmth.h	(revision 1301)
@@ -151,15 +151,15 @@
      .                32, "ave(X)", zsto,zout)
 c
-         CALL histdef(nid_mth, "fluxvdf", "PBL net flux","W/m2",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert,
-     .                32, "ave(X)", zsto,zout)
-c
-         CALL histdef(nid_mth, "fluxdyn", "Dyn. net flux","W/m2",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert,
-     .                32, "ave(X)", zsto,zout)
-c
-         CALL histdef(nid_mth, "fluxajs", "Dry adj. net flux","W/m2",
-     .                iim,jj_nb,nhori, klev,1,klev,nvert,
-     .                32, "ave(X)", zsto,zout)
+c        CALL histdef(nid_mth, "fluxvdf", "PBL net flux","W/m2",
+c    .                iim,jj_nb,nhori, klev,1,klev,nvert,
+c    .                32, "ave(X)", zsto,zout)
+c
+c        CALL histdef(nid_mth, "fluxdyn", "Dyn. net flux","W/m2",
+c    .                iim,jj_nb,nhori, klev,1,klev,nvert,
+c    .                32, "ave(X)", zsto,zout)
+c
+c        CALL histdef(nid_mth, "fluxajs", "Dry adj. net flux","W/m2",
+c    .                iim,jj_nb,nhori, klev,1,klev,nvert,
+c    .                32, "ave(X)", zsto,zout)
 c
 c        CALL histdef(nid_mth, "fluxec", "Cin. net flux","W/m2",
Index: /trunk/LMDZ.VENUS/libf/phyvenus/load_ksi.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/load_ksi.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/load_ksi.F	(revision 1301)
@@ -23,5 +23,6 @@
 C     MODIFICATIONS.
 C     --------------
-C        version multimatrice (topographie, sommet nuages): 20/12/2006
+C        
+c   New ksi matrix: possibility of different cloud model fct of lat   05/2014
 C     ------------------------------------------------------------------
 C
@@ -30,17 +31,16 @@
 c inputs
       real   psurf(klon)           ! Surface pressure
-      real   ztop(klon)            ! Altitude of the top of cloud deck (km)
 c outputs
       real   ksive(0:klev+1,0:klev+1,nnuve,nbmat)  ! ksi matrixes in Vincent's file
 
 c local variables
-      integer i,j,ig,band,pve,nlve
-      integer mat,Nb,m,Nmat,nl_init,mat0
-      parameter(nl_init=8)
+      integer i,j,isza,ips,band,pve,sve,nlve
+      integer lat,Nb,m,mat
       character*9 tmp1
       character*100 file
+      CHARACTER*2 str2
       real   lambda(nnuve)            ! wavelenght in table (mu->m, middle of interval)
       real   lambdamin(nnuve),lambdamax(nnuve) ! in microns
-      real   dlambda                  ! cm-1  
+      real   dlambda                  ! in microns 
 
       nlve = klev
@@ -53,28 +53,32 @@
       open(10,file=file)
       
-      do i=1,nl_init-1
-         read(10,*)
-      enddo
-      read(10,*) (tmp1,i=1,4),Nmat
+      read(10,*)
+      read(10,*) nlatve
+      read(10,*)
 
-      if (nbmat.ne.Nmat) then
-         write(*,*) 'This is subroutine load_ksi'
-         print*,'Probleme de dimension entre ksi.txt et le param nbmat'
-         print*,'Nb matrices = ',nbmat,Nmat
-         stop
-      endif
-
-      do mat=1,nbmat
+      write(*,*) 'This is subroutine load_ksi'
+      write(*,*) 'Nb of lat bands:',nlatve
+      
+      do lat=1,nlatve
+        read(10,*) !line for lat range
+        read(10,*) indexve(lat)
+        read(10,*) nbpsve(lat)
         read(10,*)
-        read(10,*) 
+        read(10,*) nbszave(lat)
+        read(10,*)
+	
+	do isza=1,nbszave(lat)
+	  do ips=1,nbpsve(lat)
+	  
+        read(10,*) (tmp1,j=1,3),mat    !line for matrix number
         read(10,*) (tmp1,j=1,2),pve
-        psurfve(mat) = pve*1.e5  ! pve en bar, psurfve en Pa
-        read(10,*) (tmp1,j=1,7),ztopve(mat)
-        ztopve(mat) = ztopve(mat)*1.e-3 ! passage en km
+        psurfve(ips,lat) = pve*1.e5    ! pve in bar, psurfve in Pa
+        read(10,*) (tmp1,j=1,3),sve 
+	szave(isza,lat) = cos(sve*RPI/180.) ! conversion in mu0
         read(10,*)
         read(10,*) m,Nb
         if (m.ne.nlve) then
          write(*,*) 'This is subroutine load_ksi'
-         print*,'Probleme de dimension entre ksi.txt et le param nlve'
+         print*,'Dimension problem between ksi.txt and nlve'
          print*,'N levels = ',m,nlve
          stop
@@ -82,40 +86,53 @@
         if (Nb.ne.nnuve) then
          write(*,*) 'This is subroutine load_ksi'
-         print*,'Probleme de dimension entre ksi.txt et le param nnuve'
+         print*,'Dimension problem between ksi.txt and nnuve'
          print*,'N freq = ',Nb,nnuve
          stop
         endif
 c     Now reading ksi matrix index "mat"
+        write(str2,'(i2.2)') m+2
         do band=1,Nb
          read(10,*) lambdamin(band),lambdamax(band)
          do i=0,m+1
-            read(10,'(100e17.9)') (ksive(i,j,band,mat),j=0,m+1) ! sr/µm/cm¯¹
+            read(10,'('//str2//'e17.9)') (ksive(i,j,band,mat),j=0,m+1) ! no unit
          enddo                  ! i
         enddo                     ! band
 c       print*,"Matrice ",mat," lue"
-c       print*,"   psurf=",psurfve(mat)," bars, Ztop=",ztopve(mat)," km"
-      enddo  ! mat
+c       print*,"   psurf=",psurfve(ips,lat)," bars"
+        if (mat+1.gt.nbmat) then
+         write(*,*) 'This is subroutine load_ksi'
+         print*,'Dimension problem between ksi.txt and nbmat'
+         print*,'(max number of matrixes)'
+         print*,'nbmat (in comcstVE.h) should be raised'
+         stop
+        endif
+
+	  enddo    ! ips
+	enddo      ! isza
+      enddo        ! lat
+      
+      write(*,*) 'Total nb of matrixes:',mat
       
       close(10)
 
-c longueur d'onde centrale et largeur de chaque bande
+c central wavelength and wavelength width
       do band=1,nnuve
-         lambda(band)=(lambdamin(band)+lambdamax(band))/2.*1.e-6   ! en m
-         dlambda     =(1./lambdamin(band)-1./lambdamax(band))*1.e4 ! en cm-1
+         lambda(band)=(lambdamin(band)+lambdamax(band))/2.*1.e-6   ! in m
+         dlambda     =(lambdamax(band)-lambdamin(band))            ! in microns
 c        print*,band,lambdamin(band),dlambda,lambdamax(band)
 
-c changement de convention (signe) pour ksi,
-c et prise en compte de la largeur de bande (en cm-1):
+c sign convention for ksi,
+c and taking into account the wavelength width (in microns):
          do mat=1,nbmat
          do i=0,nlve+1
            do j=0,nlve+1
-              ksive(i,j,band,mat) = -ksive(i,j,band,mat)*dlambda
+              ksive(i,j,band,mat) = +ksive(i,j,band,mat)*dlambda    ! in µm
            enddo
          enddo
          enddo
-c calcul des coeff al et bl pour luminance Planck 
+c computing coeff al and bl for Planck luminance 
          al(band) = 2.*RHPLA*RCLUM*RCLUM/(lambda(band))**5.
-c cette luminance doit etre en W/m²/sr/µm pour correspondre au calcul
-c des ksi. Ici, elle est en W/m²/sr/m donc il faut mettre un facteur 1.e-6
+c in W/m²/m
+c We need W/m²/µm :
      .                * 1.e-6
          bl(band) = RHPLA*RCLUM/(RKBOL*lambda(band))
Index: /trunk/LMDZ.VENUS/libf/phyvenus/lw_venus_ve.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/lw_venus_ve.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/lw_venus_ve.F	(revision 1301)
@@ -42,5 +42,5 @@
 c output
 
-      REAL   PCOOL(klev) ! LONGWAVE COOLING (K/VENUSDAY) within each layer
+      REAL   PCOOL(klev)  ! LONGWAVE COOLING (K/s) within each layer
       REAL   PTOPLW       ! LONGWAVE FLUX AT T.O.A. (net, + vers le haut)
       REAL   PSOLLW       ! LONGWAVE FLUX AT SURFACE (net, + vers le haut)
@@ -51,7 +51,4 @@
 C* LOCAL VARIABLES:
 C
-      real    dureejour
-      parameter (dureejour=10.087e6)
-      
       integer i,j,p
       real   zlnet(klev+1)    ! net thermal flux (W/m**2)
@@ -127,8 +124,4 @@
 c     print*,PCOOL
 
-      do j=1,klev
-        PCOOL(j) = PCOOL(j)*dureejour ! K/Venusday
-      enddo
-
       return
       end
Index: /trunk/LMDZ.VENUS/libf/phyvenus/newstart.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/newstart.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/newstart.F	(revision 1301)
@@ -69,5 +69,6 @@
       REAL tsurf(ngridmx),tsoil(ngridmx,nsoilmx)
       REAL albe(ngridmx),radsol(ngridmx),sollw(ngridmx)
-      real solsw(ngridmx),dlw(ngridmx)
+      real solsw(ngridmx),fder(ngridmx)
+      real sollwdown(ngridmx),dlw(ngridmx)
       REAL zmea(ngridmx), zstd(ngridmx)
       REAL zsig(ngridmx), zgam(ngridmx), zthe(ngridmx)
@@ -84,5 +85,6 @@
       real tsurfS(iip1,jjp1),tsoilS(iip1,jjp1,nsoilmx)
       real albeS(ip1jmp1),radsolS(ip1jmp1),sollwS(ip1jmp1)
-      real solswS(ip1jmp1),dlwS(ip1jmp1)
+      real solswS(ip1jmp1),fderS(ip1jmp1)
+      real sollwdownS(ip1jmp1),dlwS(ip1jmp1)
       real zmeaS(ip1jmp1),zstdS(ip1jmp1),zsigS(ip1jmp1)
       real zgamS(ip1jmp1),ztheS(ip1jmp1),zpicS(ip1jmp1)
@@ -116,5 +118,6 @@
       real, dimension(:,:),   allocatable :: albeold,radsolold
       real, dimension(:,:),   allocatable :: sollwold,solswold
-      real, dimension(:,:),   allocatable :: dlwold
+      real, dimension(:,:),   allocatable :: fderold
+      real, dimension(:,:),   allocatable :: dlwold,sollwdownold
       real, dimension(:,:),   allocatable :: zmeaold,zstdold,zsigold
       real, dimension(:,:),   allocatable :: zgamold,ztheold,zpicold
@@ -135,5 +138,5 @@
       integer, dimension(4) :: start,counter
       REAL phisinverse(iip1,jjp1)  ! geopotentiel au sol avant inversion
-      logical topoflag,albedoflag,razvitu
+      logical topoflag,albedoflag,razvitu,razvitv 
       real    albedo
       
@@ -318,5 +321,6 @@
       allocate(albeold(imold+1,jmold+1),radsolold(imold+1,jmold+1))
       allocate(sollwold(imold+1,jmold+1),solswold(imold+1,jmold+1))
-      allocate(dlwold(imold+1,jmold+1))
+      allocate(fderold(imold+1,jmold+1))
+      allocate(sollwdownold(imold+1,jmold+1),dlwold(imold+1,jmold+1))
       allocate(zmeaold(imold+1,jmold+1),zstdold(imold+1,jmold+1))
       allocate(zsigold(imold+1,jmold+1),zgamold(imold+1,jmold+1))
@@ -616,4 +620,19 @@
       ENDIF
 c
+      ierr = NF_INQ_VARID (nid, "fder", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <fder> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,fderold)
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,fderold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <fder>"
+         CALL abort
+      ENDIF
+c
       ierr = NF_INQ_VARID (nid, "dlw", nvarid)
       IF (ierr .NE. NF_NOERR) THEN
@@ -628,4 +647,19 @@
       IF (ierr .NE. NF_NOERR) THEN
          PRINT*, "lect_start_archive: Lecture echouee pour <dlw>"
+         CALL abort
+      ENDIF
+c
+      ierr = NF_INQ_VARID (nid, "sollwdown", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Le champ <sollwdown> est absent"
+         CALL abort
+      ENDIF
+#ifdef NC_DOUBLE
+      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,counter,sollwdownold)
+#else
+      ierr = NF_GET_VARA_REAL(nid, nvarid,start,counter,sollwdownold)
+#endif
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: Lecture echouee pour <sollwdown>"
          CALL abort
       ENDIF
@@ -844,4 +878,6 @@
        topoflag = . FALSE .
        CALL getin('topoflag',topoflag)
+
+        print*,zmeaold(2,1:10)
 
        IF ( topoflag ) then
@@ -952,7 +988,15 @@
       call gr_dyn_fi (1,iip1,jjp1,ngridmx,solswS,solsw)
 
+      call interp_horiz (fderold,fderS,imold,jmold,iim,jjm,1,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iip1,jjp1,ngridmx,fderS,fder)
+
       call interp_horiz (dlwold,dlwS,imold,jmold,iim,jjm,1,
      &                   rlonuold,rlatvold,rlonu,rlatv)
       call gr_dyn_fi (1,iip1,jjp1,ngridmx,dlwS,dlw)
+
+      call interp_horiz (sollwdownold,sollwdownS,imold,jmold,iim,jjm,1,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iip1,jjp1,ngridmx,sollwdownS,sollwdown)
 
       print*,"Nouvelles var physiques OK"
@@ -1044,4 +1088,6 @@
        razvitu = . FALSE .
        CALL getin('razvitu',razvitu)
+       razvitv = . FALSE .
+       CALL getin('razvitv',razvitv)
 
 c calcul des champ de vent; passage en vent covariant
@@ -1053,5 +1099,5 @@
       call interp_horiz(var,us,imold,jmold,iim,jjm,llm,
      &                   rlonuold,rlatvold,rlonu,rlatv)
-      write (*,*) 'us ', us (1,2,1)   ! INFO
+      write (*,*) 'us ', us (1,2,1)    ! INFO
 
       call interp_vert
@@ -1078,4 +1124,9 @@
       write (*,*) 'ucov ', ucov (1,2,1)  ! INFO
 c     write(48,*) 'ucov',ucov
+! Reseting v=0
+      if (razvitv) then
+           vnat(:,:,:) = 0.
+      endif
+      write (*,*) 'vnat ', vnat (1,2,1)    ! INFO
       do l=1,llm
         do j = 1, jjm
@@ -1114,6 +1165,6 @@
 
       call writerestartphy('restartphy.nc',tab_cntrl_fi,ngridmx,llm,
-     .           rlat,rlon,tsurf,tsoil,albe,solsw, sollw,dlw,
-     .           radsol,
+     .           rlat,rlon,tsurf,tsoil,albe,solsw, sollw,fder,dlw,
+     .           sollwdown,radsol,
      .           zmea,zstd,zsig,zgam,zthe,zpic,zval,
      .           t_fi)
Index: /trunk/LMDZ.VENUS/libf/phyvenus/phyetat0.F90
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/phyetat0.F90	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/phyetat0.F90	(revision 1301)
@@ -231,12 +231,11 @@
       PRINT*,'Rayonnement net au sol radsol:', xmin, xmax
 
-! Lecture de l'orographie sous-maille si ok_orodr:
-
-      if(ok_orodr) then
-     
+! Lecture de l'orographie sous-maille:
+
       CALL get_field("ZMEA",zmea,found)
       IF (.not.found) THEN
          PRINT*, 'phyetat0: Le champ <ZMEA> est absent'
-         CALL abort
+         PRINT*, 'mis a zero'
+         zmea=0.
       ENDIF
       xmin = 1.0E+20
@@ -251,5 +250,6 @@
       IF (.not.found) THEN
          PRINT*, 'phyetat0: Le champ <ZSTD> est absent'
-         CALL abort
+         PRINT*, 'mis a zero'
+         zstd=0.
       ENDIF
       xmin = 1.0E+20
@@ -264,5 +264,6 @@
       IF (.not.found) THEN
          PRINT*, 'phyetat0: Le champ <ZSIG> est absent'
-         CALL abort
+         PRINT*, 'mis a zero'
+         zsig=0.
       ENDIF
       xmin = 1.0E+20
@@ -277,5 +278,6 @@
       IF (.not.found) THEN
          PRINT*, 'phyetat0: Le champ <ZGAM> est absent'
-         CALL abort
+         PRINT*, 'mis a zero'
+         zgam=0.
       ENDIF
       xmin = 1.0E+20
@@ -290,5 +292,6 @@
       IF (.not.found) THEN
          PRINT*, 'phyetat0: Le champ <ZTHE> est absent'
-         CALL abort
+         PRINT*, 'mis a zero'
+         zthe=0.
       ENDIF
       xmin = 1.0E+20
@@ -303,5 +306,6 @@
       IF (.not.found) THEN
          PRINT*, 'phyetat0: Le champ <ZPIC> est absent'
-         CALL abort
+         PRINT*, 'mis a zero'
+         zpic=0.
       ENDIF
       xmin = 1.0E+20
@@ -316,5 +320,6 @@
       IF (.not.found) THEN
          PRINT*, 'phyetat0: Le champ <ZVAL> est absent'
-         CALL abort
+         PRINT*, 'mis a zero'
+         zval=0.
       ENDIF
       xmin = 1.0E+20
@@ -325,15 +330,4 @@
       ENDDO
       PRINT*,'OROGRAPHIE SOUS-MAILLE zval:', xmin, xmax
-
-      else 
-         zmea = 0.
-         zstd = 0.
-         zsig = 0.
-         zgam = 0.
-         zthe = 0.
-         zpic = 0.
-         zval = 0.
-
-      endif   ! fin test sur ok_orodr
 
 ! Lecture de TANCIEN:
Index: /trunk/LMDZ.VENUS/libf/phyvenus/physiq.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/physiq.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/physiq.F	(revision 1301)
@@ -7,5 +7,5 @@
      .            paprs,pplay,ppk,pphi,pphis,presnivs,
      .            u,v,t,qx,
-     .            omega,
+     .            flxmw,
      .            d_u, d_v, d_t, d_qx, d_ps)
 
@@ -46,5 +46,5 @@
 c qx------input-R-mass mixing ratio traceurs (kg/kg) 
 c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
-c omega---input-R-vitesse verticale en Pa/s
+c flxmw---input-R-flux de masse vertical en kg/s
 c
 c d_u-----output-R-tendance physique de "u" (m/s/s)
@@ -125,5 +125,5 @@
       REAL d_t_dyn(klon,klev)
 
-      REAL omega(klon,klev)
+      REAL flxmw(klon,klev)
 
       REAL d_u(klon,klev)
@@ -146,7 +146,8 @@
 c Variables propres a la physique
 c
-      REAL,save,allocatable :: rlev(:,:) ! altitude a chaque niveau (interface inferieure de la couche) 
       INTEGER,save :: itap        ! compteur pour la physique
       REAL delp(klon,klev)        ! epaisseur d'une couche
+      REAL omega(klon,klev)       ! vitesse verticale en Pa/s
+
       
       INTEGER igwd,idx(klon),itest(klon)
@@ -242,5 +243,6 @@
 c
       REAL zphi(klon,klev)
-c
+      REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
+
 c Variables du changement
 c
@@ -362,5 +364,4 @@
 c========================
       IF (debut) THEN
-         allocate(rlev(klon,klevp1))
          allocate(source(klon,nqmax))
 
@@ -621,4 +622,12 @@
       ENDIF
 c====================================================================
+
+c Calcule de vitesse verticale a partir de flux de masse verticale
+      DO k = 1, klev
+       DO i = 1, klon
+        omega(i,k) = RG*flxmw(i,k) / airephy(i)
+       END DO
+      END DO
+
 c
 c Ajouter le geopotentiel du sol:
@@ -629,4 +638,27 @@
       ENDDO
       ENDDO
+
+c   calcul du geopotentiel aux niveaux intercouches
+c   ponderation des altitudes au niveau des couches en dp/p
+
+      DO k=1,klev
+         DO i=1,klon
+            zzlay(i,k)=zphi(i,k)/RG
+         ENDDO
+      ENDDO
+      DO i=1,klon
+         zzlev(i,1)=pphis(i)/RG
+      ENDDO
+      DO k=2,klev
+         DO i=1,klon
+            z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k))
+            z2=(paprs(i,k)  +pplay(i,k))/(paprs(i,k)  -pplay(i,k))
+            zzlev(i,k)=(z1*zzlay(i,k-1)+z2*zzlay(i,k))/(z1+z2)
+         ENDDO
+      ENDDO
+      DO i=1,klon
+         zzlev(i,klev+1)=zzlay(i,klev)+(zzlay(i,klev)-zzlev(i,klev))
+      ENDDO
+
 c====================================================================
 c
@@ -729,8 +761,4 @@
       fder = dlw
 
-c     print*,"radsol avant clmain=",radsol(klon/2)
-c     print*,"solsw avant clmain=",solsw(klon/2)
-c     print*,"sollw avant clmain=",sollw(klon/2)
-
 ! ADAPTATION GCM POUR CP(T)
 
@@ -749,8 +777,4 @@
      s            ycoefh,yu1,yv1) 
 
-c     print*,"radsol apres clmain=",radsol(klon/2)
-c     print*,"solsw apres clmain=",solsw(klon/2)
-c     print*,"sollw apres clmain=",sollw(klon/2)
-
 CXXX Incrementation des flux
       DO i = 1, klon
@@ -770,10 +794,4 @@
       ENDDO
       ENDDO
-c
-c        print*,"d_t_vdf1=",d_t_vdf(1,:)*dtime
-c        print*,"d_t_vdf2=",d_t_vdf(klon/2,:)*dtime
-c        print*,"d_t_vdf3=",d_t_vdf(klon,:)*dtime
-c        print*,"d_u_vdf=",d_u_vdf(klon/2,:)*dtime
-c        print*,"d_v_vdf=",d_v_vdf(klon/2,:)*dtime
 
 C TRACEURS
@@ -810,9 +828,7 @@
 c Incrementer la temperature du sol
 c
-c     print*,'Tsol avant clmain:',ftsol(1)
       DO i = 1, klon
          ftsol(i) = ftsol(i) + d_ts(i)
       ENDDO
-c     print*,'Tsol apres clmain:',ftsol(1)
 
 c Calculer la derive du flux infrarouge
@@ -878,14 +894,9 @@
          v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
          d_v_ajs(:,:)= d_v_ajs(:,:)/dtime          ! (m/s)/s
-      if (iflag_trac.eq.1) then
+
+         if (iflag_trac.eq.1) then
            tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
            d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime  ! /s
-      endif
-
-c        print*,"d_t_ajs1=",d_t_ajs(1,:)*dtime
-c        print*,"d_t_ajs2=",d_t_ajs(klon/2,:)*dtime
-c        print*,"d_t_ajs3=",d_t_ajs(klon,:)*dtime
-c        print*,"d_u_ajs=",d_u_ajs(klon/2,:)*dtime
-c        print*,"d_v_ajs=",d_v_ajs(klon/2,:)*dtime
+         endif
 
       endif
@@ -919,11 +930,6 @@
 c     PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
             
-c     print*,"radsol avant radlwsw=",radsol(klon/2)
-c     print*,"solsw avant radlwsw=",solsw(klon/2)
-c     print*,"sollw avant radlwsw=",sollw(klon/2)
-c     print*,"avant radlwsw"
-
       CALL radlwsw
-     e            (dist, rmu0, fract, 
+     e            (dist, rmu0, fract, zzlev,
      e             paprs, pplay,ftsol, t_seri,
      s             heat,cool,radsol,
@@ -932,19 +938,11 @@
      s             lwnet, swnet)
 
-c     print*,"apres radlwsw"
-c     print*,"radsol apres radlwsw=",radsol(klon/2)
-c     print*,"solsw apres radlwsw=",solsw(klon/2)
-c     print*,"sollw apres radlwsw=",sollw(klon/2)
       itaprad = 0
       DO k = 1, klev
       DO i = 1, klon
-         dtrad(i,k) = heat(i,k)-cool(i,k)
-         dtrad(i,k) = dtrad(i,k)/RDAY  !K/s
-      ENDDO
-      ENDDO
-c        print*,"dtrad1=",dtrad(1,:)*dtime
-c        print*,"dtrad2=",dtrad(klon/2,:)*dtime
-c        print*,"dtrad3=",dtrad(klon,:)*dtime
-      
+         dtrad(i,k) = heat(i,k)-cool(i,k)  !K/s
+      ENDDO
+      ENDDO
+
       ENDIF
       itaprad = itaprad + 1
Index: /trunk/LMDZ.VENUS/libf/phyvenus/phytrac_emiss.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/phytrac_emiss.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/phytrac_emiss.F	(revision 1301)
@@ -77,5 +77,6 @@
       integer,parameter :: nbsrc=2,nblat=5,nblon=4
       integer,parameter :: Nemiss=1   ! duree emission (Ed)
-      real,save :: source_volcan(nbsrc)
+      integer,save :: Nemiss(nbsrc)      ! duration emission (Ed)
+      real,save :: source_volcan(nbsrc)  ! flux emission (kg/s)
       real,save :: lat_volcan(nblat),lon_volcan(nblon)
       real,save :: area_emiss(nblat,nblon)
@@ -131,4 +132,7 @@
          source_volcan(1) = 1.
          source_volcan(2) = 1000.
+c duration in Ed
+         Nemiss(1) = 1
+         Nemiss(2) = 10
 c localisation volcan
          lat_volcan(1) =  70.
@@ -173,10 +177,12 @@
            do ilat  = 1,nblat
             do ilon  = 1,nblon
-             it=min(iemiss*ilat*ilon,nqtot)
+             it=(iemiss-1)*nblat*nblon+(ilat-1)*nblon+ilon
+             it=min(it,nqtot)
+             deltatr(i,1,it) = 0.
 
              if (i .eq. ig_volcan(ilat,ilon)) then
 
 c source appliquee pendant Nemiss Ed
-               if (timesimu .lt. 86400.*Nemiss) then
+               if (timesimu .lt. 86400.*Nemiss(iemiss)) then
 
 c source en kg/kg/s
@@ -185,8 +191,5 @@
            tr_seri(i,1,it) = tr_seri(i,1,it) + deltatr(i,1,it)*pdtphys
 
-               else
-           deltatr(i,1,it) = 0.
                end if  ! duree emission
-
              end if ! i localisation
             end do
Index: /trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.F	(revision 1301)
@@ -2,5 +2,5 @@
 ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/radlwsw.F,v 1.2 2004/10/27 10:14:46 lmdzadmin Exp $
 !
-      SUBROUTINE radlwsw(dist, rmu0, fract, 
+      SUBROUTINE radlwsw(dist, rmu0, fract, zzlev,
      .                  paprs, pplay,tsol, t,
      .                  heat,cool,radsol,
@@ -17,10 +17,11 @@
 c fract----input-R- duree d'ensoleillement normalisee
 c solaire--input-R- constante solaire (W/m**2) (dans clesphys.h)
+c zzlev----input-R- altitude a inter-couche (m)
 c paprs----input-R- pression a inter-couche (Pa)
 c pplay----input-R- pression au milieu de couche (Pa)
 c tsol-----input-R- temperature du sol (en K)
 c t--------input-R- temperature (K)
-c heat-----output-R- echauffement atmospherique (visible) (K/jour)
-c cool-----output-R- refroidissement dans l'IR (K/jour)
+c heat-----output-R- echauffement atmospherique (visible) (K/s)
+c cool-----output-R- refroidissement dans l'IR (K/s)
 c radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas)
 c topsw----output-R- flux solaire net au sommet de l'atm. (+ vers le bas)
@@ -36,4 +37,5 @@
 c   S. Lebonnois    20/12/2006
 c   corrections     13/07/2007
+c   New ksi matrix: possibility of different cloud model fct of lat   05/2014
 
 c======================================================================
@@ -46,7 +48,11 @@
 #include "clesphys.h" 
 #include "comcstVE.h"
-c
+
+!===========
+! Arguments
+!===========
       real rmu0(klon), fract(klon), dist
-c
+
+      REAL zzlev(klon,klev+1)
       real paprs(klon,klev+1), pplay(klon,klev)
       real tsol(klon)
@@ -57,13 +63,16 @@
       real sollwdown(klon)
       REAL swnet(klon,klev+1),lwnet(klon,klev+1)
-c
+
+!===========
+! Local
+!===========
       INTEGER k, kk, i, j, band
-c
+
       REAL   PPB(klev+1)
-c
+
       REAL   zfract, zrmu0
-c
+
       REAL   zheat(klev), zcool(klev)
-      real  temp(klev)
+      real   temp(klev),znivs(klev+1)
       REAL   ZFSNET(klev+1),ZFLNET(klev+1)
       REAL   ztopsw, ztoplw
@@ -73,14 +82,12 @@
 cIM END
       real,save,allocatable :: ksive(:,:,:,:) ! ksi matrixes in Vincent's file
-      real,save,allocatable :: ztop(:) ! in km
 
       real    psi(0:klev+1,0:klev+1)
       real    deltapsi(0:klev+1,0:klev+1)
-      real    latdeg
       real    pt0(0:klev+1)
       real    bplck(0:klev+1,nnuve)    ! Planck luminances in table layers
-      real    y(0:klev,nnuve)          ! intermediaire Planck
-      real    zdblay(0:klev+1,nnuve)   ! gradient en temperature de planck
-      integer mat,mat0
+      real    y(0:klev,nnuve)          ! temporary variable for Planck
+      real    zdblay(0:klev+1,nnuve)   ! temperature gradient of planck function
+      integer mat0,lat,ips,isza,ips0,isza0
       real    factp,factz,ksi
 
@@ -98,24 +105,4 @@
         allocate(ksive(0:klev+1,0:klev+1,nnuve,nbmat))
         call load_ksi(ksive)
-c ---------- ztop --------------
-        allocate(ztop(klon))
-        DO i = 1, klon
-             ztop(i) = 70.
-        ENDDO !i
-c ztop: d'apres fit à figure 16 du papier Zavosa et al (tmp) traitant des
-c       donnees Venera
-c       DO i = 1, klon
-c         latdeg = abs(rlatd(i))
-c         if (latdeg.lt.15) then
-c            ztop(i) = 70.
-c         elseif (latdeg.lt.50) then
-c            ztop(i) = 63.95+6*cos((latdeg-15)*RPI/2./50.)
-c         else 
-c            ztop(i) = min(63.95+6*cos((latdeg-15)*RPI/2./50.),
-c    .                     63.95-5.9*sin((latdeg-60)*RPI/2/30))
-c         endif 
-c       print*,'lat(',i,')=',latdeg,'  ztop=',ztop(i)
-c       ENDDO !i
-c ---------- ztop --------------
 
       endif ! firstcall
@@ -191,25 +178,60 @@
 c --------------------------
        mat0 = 0
-       do mat=1,nbmat-nbztopve
-         if (  (psurfve(mat).ge.paprs(j,1))
-     .    .and.(psurfve(mat+nbztopve).lt.paprs(j,1))
-     .    .and.(ztopve(mat).lt.ztop(j))
-     .    .and.(ztopve(mat+1).ge.ztop(j)) ) then
-              mat0  = mat
-c             print*,'ig=',j,'  mat0=',mat
-              factp = (paprs(j,1)           -psurfve(mat))
-     .               /(psurfve(mat+nbztopve)-psurfve(mat))
-              factz = (ztop(j)      -ztopve(mat))
-     .               /(ztopve(mat+1)-ztopve(mat))
+       if (nlatve.eq.1) then   ! clouds are taken as uniform
+         lat=1
+       else
+         if (abs(rlatd(j)).le.50.) then
+	   lat=1
+	 elseif (abs(rlatd(j)).le.60.) then
+	   lat=2
+	 elseif (abs(rlatd(j)).le.70.) then
+	   lat=3
+	 elseif (abs(rlatd(j)).le.80.) then
+	   lat=4
+	 else
+	   lat=5
+	 endif
+       endif
+       
+       ips0=0
+       do ips=1,nbpsve(lat)-1
+         if (  (psurfve(ips,lat).ge.paprs(j,1))
+     .    .and.(psurfve(ips+1,lat).lt.paprs(j,1)) ) then
+              ips0  = ips
+c             print*,'ig=',j,'  ips0=',ips
+              factp = (paprs(j,1)         -psurfve(ips0,lat))
+     .               /(psurfve(ips0+1,lat)-psurfve(ips0,lat))
               exit
          endif
        enddo
-       if (mat0.eq.0) then
+       isza0=0
+       if (nbszave(lat).gt.1) then
+        do isza=1,nbszave(lat)-1
+         if (  (szave(isza,lat).ge.zrmu0)
+     .    .and.(szave(isza+1,lat).lt.zrmu0) ) then
+              isza0  = isza
+c             print*,'ig=',j,'  isza0=',isza
+              factz = (zrmu0             -szave(isza0,lat))
+     .               /(szave(isza0+1,lat)-szave(isza0,lat))
+              exit
+         endif
+        enddo
+       else            ! Only one sza, no interpolation
+        isza0=-99
+       endif
+       
+       if ((ips0.eq.0).or.(isza0.eq.0)) then
          write(*,*) 'Finding the right matrix in radlwsw'
-         print*,'Probleme pour interpolation au point ig=',j
-         print*,'psurf = ',paprs(j,1),' ztop = ',ztop(j)
+         print*,'Interpolation problem, grid point ig=',j
+         print*,'psurf = ',paprs(j,1),' mu0 = ',zrmu0
          stop
        endif
-
+       
+       if (isza0.eq.-99) then
+         mat0 = indexve(lat)+ips0
+       else
+         mat0 = indexve(lat)+(isza0-1)*nbpsve(lat)+ips0
+       endif
+       
 c interpolation of ksi and computation of psi,deltapsi
 c ----------------------------------------------------
@@ -217,11 +239,17 @@
         do k=0,klev+1
          do i=0,klev+1
-          ksi = ksive(i,k,band,mat0)*(1-factz)*(1-factp)
-     .         +ksive(i,k,band,mat0+1)*factz  *(1-factp)
-     .         +ksive(i,k,band,mat0+nbztopve)*(1-factz)*factp
-     .         +ksive(i,k,band,mat0+nbztopve+1)*factz  *factp
+          if (isza0.eq.-99) then
+           ksi = ksive(i,k,band,mat0)*(1-factp)
+     .          +ksive(i,k,band,mat0+1)*factp
+          else
+           ksi = ksive(i,k,band,mat0)*(1-factp)*(1-factz)
+     .          +ksive(i,k,band,mat0+1)*factp  *(1-factz)
+     .          +ksive(i,k,band,mat0+nbpsve(lat))*(1-factp)*factz
+     .          +ksive(i,k,band,mat0+nbpsve(lat)+1)*factp  *factz
+          endif
+     
           psi(i,k) = psi(i,k) + 
-     .               ksi*(bplck(i,band)-bplck(k,band))
-          deltapsi(i,k) = deltapsi(i,k) + ksi*zdblay(i,band)
+     .               RPI*ksi*(bplck(i,band)-bplck(k,band))
+          deltapsi(i,k) = deltapsi(i,k) + RPI*ksi*zdblay(i,band)
          enddo
         enddo
@@ -241,4 +269,10 @@
 c SW call
 c---------
+      znivs=zzlev(j,:)
+c      CALL SW_venus_ve(zrmu0, zfract,
+c     S        PPB,temp,znivs,
+c     S        zheat, 
+c     S        ztopsw,zsolsw,ZFSNET)
+
       CALL SW_venus_dc(zrmu0, zfract,
      S        PPB,temp, 
@@ -283,11 +317,11 @@
 c     j = 1
 c     print*,'mu0=',rmu0(j)
-c     print*,'   net flux vis   HEAT(K/day)'
+c     print*,'   net flux vis   HEAT(K/Eday)'
 c     do k=1,klev
-c     print*,k,ZFSNET(k),heat(j,k)*8.56548e-3
+c     print*,k,ZFSNET(k),heat(j,k)*86400.
 c     enddo
-c     print*,'   net flux IR    COOL(K/day)'
+c     print*,'   net flux IR    COOL(K/Eday)'
 c     do k=1,klev
-c     print*,k,ZFLNET(k),cool(j,k)*8.56548e-3
+c     print*,k,ZFLNET(k),cool(j,k)*86400.
 c     enddo
 
Index: /trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.NewtonCool
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.NewtonCool	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.NewtonCool	(revision 1301)
@@ -2,11 +2,7 @@
 ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/radlwsw.F,v 1.2 2004/10/27 10:14:46 lmdzadmin Exp $
 !
-      SUBROUTINE radlwsw(dist, rmu0, fract, 
-     .                  paprs, pplay,tsol, t,
-     .                  heat,cool,radsol,
-     .                  topsw,toplw,solsw,sollw,
-     .                  sollwdown,
-     .                  lwnet, swnet)
-c      
+      SUBROUTINE radlwsw(dist, rmu0, fract, zzlev,
+     .                  paprs, pplay,tsol, pt, nq, nmicro, pq,qaer)
+      
 c======================================================================
 c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
@@ -16,43 +12,36 @@
 c rmu0-----input-R- cosinus de l'angle zenithal
 c fract----input-R- duree d'ensoleillement normalisee
-c solaire--input-R- constante solaire (W/m**2) (dans clesphys.h)
 c paprs----input-R- pression a inter-couche (Pa)
 c pplay----input-R- pression au milieu de couche (Pa)
 c tsol-----input-R- temperature du sol (en K)
-c t--------input-R- temperature (K)
-c heat-----output-R- echauffement atmospherique (visible) (K/jour)
-c cool-----output-R- refroidissement dans l'IR (K/jour)
-c radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas)
-c topsw----output-R- flux solaire net au sommet de l'atm. (+ vers le bas)
-c toplw----output-R- ray. IR net au sommet de l'atmosphere (+ vers le haut)
-c solsw----output-R- flux solaire net a la surface (+ vers le bas)
-c sollw----output-R- ray. IR net a la surface (+ vers le bas)
-c sollwdown-output-R- ray. IR descendant a la surface (+ vers le bas)
-c lwnet____output-R- flux IR net (+ vers le haut)
-c swnet____output-R- flux solaire net (+ vers le bas)
+c pt-------input-R- temperature (K)
 c
       
 c   S. Lebonnois    12/04/2007
 c  VERSION NEWTONIAN COOLING pour Venus (no diurnal cycle)
+c  update 01/2014
 
 c======================================================================
       use dimphy
       USE comgeomphy
+      USE phys_state_var_mod, only: falbe,heat,cool,radsol,
+     .      topsw,toplw,solsw,sollw,sollwdown,lwnet,swnet
+      USE write_field_phy
       IMPLICIT none
 #include "dimensions.h"
 #include "YOMCST.h"
 #include "clesphys.h" 
-c
+
+c ARGUMENTS
+      INTEGER nq,nmicro
       real rmu0(klon), fract(klon), dist
-c
-      real paprs(klon,klev+1), pplay(klon,klev)
+ 
+      real zzlev(klon,klev+1),paprs(klon,klev+1), pplay(klon,klev)
       real tsol(klon)
-      real t(klon,klev)
-      real heat(klon,klev), cool(klon,klev)
-      real radsol(klon), topsw(klon), toplw(klon)
-      real solsw(klon), sollw(klon)
-      real sollwdown(klon)
-      REAL swnet(klon,klev+1),lwnet(klon,klev+1)
-c
+      real pt(klon,klev)
+      real pq(klon,klev,nq)
+      REAL qaer(klon,klev,nq)
+ 
+c LOCAL VARIABLES
       INTEGER i,j,k
       integer    nlevCLee,level
@@ -64,6 +53,4 @@
       real   ztemp,zdt,fact
       real   dTsdt(klev),zt_eq(klon,klev)
-      real   dureejour
-      parameter (dureejour=10.087e6)
       save   zt_eq
       
@@ -144,5 +131,5 @@
 
         DO k = 1, klev
-         heat (j,k) = dTsdt(k)*dureejour ! K/Venusday
+         heat (j,k) = dTsdt(k)    ! K/s
          cool (j,k) = 0.
         ENDDO
Index: /trunk/LMDZ.VENUS/libf/phyvenus/rcm1d.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/rcm1d.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/rcm1d.F	(revision 1301)
@@ -6,4 +6,6 @@
       USE phys_state_var_mod
       use cpdet_mod, only: ini_cpdet
+      use moyzon_mod, only: tmoy
+
       IMPLICIT NONE
 
@@ -288,4 +290,7 @@
         print*,ilayer,play(ilayer),zlay(ilayer),temp(ilayer)
       ENDDO
+
+      allocate(tmoy(llm))
+      tmoy(:)=temp(:)
       
 c     temperature du sous-sol
@@ -465,4 +470,3 @@
 #include "../dyn3d_common/disvert_noterre.F"
 #include "../dyn3d/abort_gcm.F"
-!#include "../dyn3d/dump2d.F"
-
+
Index: /trunk/LMDZ.VENUS/libf/phyvenus/readstartphy.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/readstartphy.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/readstartphy.F	(revision 1301)
@@ -7,5 +7,5 @@
      .            rlat,rlon, tsol,tsoil,
      .           albe, solsw, sollw,
-     .           fder,radsol,
+     .           fder,dlw, sollwdown, radsol,
      .    zmea, zstd, zsig, zgam, zthe, zpic, zval,
      .           tabcntr0)
@@ -29,11 +29,10 @@
       REAL tsoil(ngridmx,nsoilmx)
       REAL albe(ngridmx)
-cIM BEG alblw
-      REAL alblw(ngridmx)
-cIM END alblw
       REAL radsol(ngridmx)
       REAL sollw(ngridmx)
       real solsw(ngridmx)
       real fder(ngridmx)
+      real dlw(ngridmx)
+      real sollwdown(ngridmx)
       REAL zmea(ngridmx), zstd(ngridmx)
       REAL zsig(ngridmx), zgam(ngridmx), zthe(ngridmx)
@@ -275,4 +274,58 @@
 
 c
+c Lecture derive des flux IR:
+c
+      ierr = NF_INQ_VARID (nid, "dlw", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Le champ <dlw> est absent'
+         PRINT*, 'mis a zero'
+         dlw = 0.
+      ELSE
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, dlw)
+#else
+        ierr = NF_GET_VAR_REAL(nid, nvarid, dlw)
+#endif
+        IF (ierr.NE.NF_NOERR) THEN
+          PRINT*, 'phyetat0: Lecture echouee pour <dlw>'
+          CALL abort
+        ENDIF
+      ENDIF
+      xmin = 1.0E+20
+      xmax = -1.0E+20
+      DO i = 1, ngridmx
+         xmin = MIN(dlw(i),xmin)
+         xmax = MAX(dlw(i),xmax)
+      ENDDO
+      PRINT*,'Derive des flux IR dlw:', xmin, xmax
+
+c
+c Lecture rayonnement IR vers le bas au sol:
+c
+      ierr = NF_INQ_VARID (nid, "sollwdown", nvarid)
+      IF (ierr.NE.NF_NOERR) THEN
+         PRINT*, 'phyetat0: Le champ <sollwdown> est absent'
+         PRINT*, 'mis a zero'
+         sollwdown = 0.
+      ELSE
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, sollwdown)
+#else
+        ierr = NF_GET_VAR_REAL(nid, nvarid, sollwdown)
+#endif
+        IF (ierr.NE.NF_NOERR) THEN
+          PRINT*, 'phyetat0: Lecture echouee pour <sollwdown>'
+          CALL abort
+        ENDIF
+      ENDIF
+      xmin = 1.0E+20
+      xmax = -1.0E+20
+      DO i = 1, ngridmx
+         xmin = MIN(sollwdown(i),xmin)
+         xmax = MAX(sollwdown(i),xmax)
+      ENDDO
+      PRINT*,'Flux IR vers le bas au sol sollwdown:', xmin, xmax
+
+c
 c Lecture du rayonnement net au sol:
 c
Index: /trunk/LMDZ.VENUS/libf/phyvenus/start2archive.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/start2archive.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/start2archive.F	(revision 1301)
@@ -18,5 +18,5 @@
       USE infotrac
       USE control_mod
-      use cpdet_mod, only: tpot2t
+      use cpdet_mod, only: tpot2t,ini_cpdet
 
       implicit none
@@ -42,5 +42,5 @@
 c -----------------------------
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
-      REAL teta(ip1jmp1,llm)                    ! temperature potentielle 
+      REAL teta(ip1jmp1,llm)                 ! temperature potentielle 
       REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
       REAL pks(ip1jmp1)                      ! exner (f pour filtre)
@@ -51,5 +51,5 @@
       REAL masse(ip1jmp1,llm)                ! masse de l'atmosphere
       REAL ps(ip1jmp1)                       ! pression au sol
-      REAL p3d(iip1, jjp1, llm+1)            ! pression aux interfaces
+      REAL p3d(iip1,jjp1,llm+1)              ! pression aux interfaces
       
 c Variable Physiques (grille physique)
@@ -62,5 +62,6 @@
       REAL tsurf(ngridmx),tsoil(ngridmx,nsoilmx)
       REAL albe(ngridmx),radsol(ngridmx),sollw(ngridmx)
-      real solsw(ngridmx),dlw(ngridmx)
+      real solsw(ngridmx),fder(ngridmx)
+      real sollwdown(ngridmx),dlw(ngridmx)
       REAL zmea(ngridmx), zstd(ngridmx)
       REAL zsig(ngridmx), zgam(ngridmx), zthe(ngridmx)
@@ -69,5 +70,5 @@
       INTEGER start,length
       PARAMETER (length = 100)
-      REAL tab_cntrl_fi(length) ! tableau des parametres de startfi
+      REAL tab_cntrl_fi(length)  ! tableau des parametres de startfi
       REAL tab_cntrl_dyn(length) ! tableau des parametres de start
       INTEGER*4 day_ini_fi
@@ -79,5 +80,6 @@
       real rlatS(ip1jmp1),rlonS(ip1jmp1)
       real albeS(ip1jmp1),radsolS(ip1jmp1),sollwS(ip1jmp1)
-      real solswS(ip1jmp1),dlwS(ip1jmp1)
+      real solswS(ip1jmp1),fderS(ip1jmp1)
+      real dlwS(ip1jmp1),sollwdownS(ip1jmp1)
       real zmeaS(ip1jmp1),zstdS(ip1jmp1),zsigS(ip1jmp1)
       real zgamS(ip1jmp1),ztheS(ip1jmp1),zpicS(ip1jmp1)
@@ -132,5 +134,5 @@
      .       rlat,rlon,tsurf,tsoil,
      .       albe, solsw, sollw,
-     .       dlw,radsol,
+     .       fder,dlw,sollwdown,radsol,
      .       zmea,zstd,zsig,zgam,zthe,zpic,zval,
      .       tab_cntrl_fi)
@@ -140,7 +142,10 @@
 c-----------------------------------------------------------------------
 
+      CALL conf_gcm( 99, .TRUE. )
       call iniconst
       call inigeom
       call inifilr
+      call ini_cpdet
+
       CALL pression(ip1jmp1, ap, bp, ps, p3d)
          if (disvert_type==1) then
@@ -198,5 +203,7 @@
       call gr_fi_dyn(1,ngridmx,iip1,jjp1,sollw,sollwS)
       call gr_fi_dyn(1,ngridmx,iip1,jjp1,solsw,solswS)
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,fder,fderS)
       call gr_fi_dyn(1,ngridmx,iip1,jjp1,dlw,dlwS)
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,sollwdown,sollwdownS)
       call gr_fi_dyn(1,ngridmx,iip1,jjp1,zmea,zmeaS)
       call gr_fi_dyn(1,ngridmx,iip1,jjp1,zstd,zstdS)
@@ -295,5 +302,8 @@
       call write_archive(nid,'solsw',
      .             'SW flux at surface','W m-2',2,solswS)
+      call write_archive(nid,'fder','derive','?',2,fderS)
       call write_archive(nid,'dlw','LW derive','?',2,dlwS)
+      call write_archive(nid,'sollwdown',
+     .             'LW dwn flux at surface','?',2,sollwdownS)
       call write_archive(nid,'zmea','param oro sous-maille','m',2,zmeaS)
       call write_archive(nid,'zstd','param oro sous-maille','m',2,zstdS)
Index: /trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_dc.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_dc.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_dc.F	(revision 1301)
@@ -40,5 +40,5 @@
 c output
 
-      REAL   PHEAT(klev) ! SHORTWAVE HEATING (K/VENUSDAY) within each layer
+      REAL   PHEAT(klev)  ! SHORTWAVE HEATING (K/s) within each layer
       REAL   PTOPSW       ! SHORTWAVE FLUX AT T.O.A. (net)
       REAL   PSOLSW       ! SHORTWAVE FLUX AT SURFACE (net)
@@ -49,8 +49,6 @@
 C
       integer nldc,nszadc
-      real    dureejour
       parameter (nldc=49)  ! fichiers Crisp
       parameter (nszadc=8) ! fichiers Crisp
-      parameter (dureejour=10.087e6)
       
       integer i,j,nsza,nsza0,nl0
@@ -234,5 +232,4 @@
         PHEAT(j) = (ZFSNET(j+1)-ZFSNET(j))
      .            *RG/cpdet(pt(j)) / ((PPB(j)-PPB(j+1))*1.e5)
-        PHEAT(j) = PHEAT(j)*dureejour  ! K/venus_day
       enddo
 
Index: /trunk/LMDZ.VENUS/libf/phyvenus/write_histday.h
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/write_histday.h	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/write_histday.h	(revision 1301)
@@ -86,8 +86,8 @@
 c en K/s      
       call histwrite_phy(nid_day,.false.,"dtajs",itau_w,d_t_ajs)
-c K/day ==> K/s
-      call histwrite_phy(nid_day,.false.,"dtswr",itau_w,heat/RDAY)
-c K/day ==> K/s      
-      call histwrite_phy(nid_day,.false.,"dtlwr",itau_w,-1.*cool/RDAY)
+c en K/s      
+      call histwrite_phy(nid_day,.false.,"dtswr",itau_w,heat)
+c en K/s      
+      call histwrite_phy(nid_day,.false.,"dtlwr",itau_w,-1.*cool)
 c en K/s      
 c     call histwrite_phy(nid_day,.false.,"dtec",itau_w,d_t_ec)
Index: /trunk/LMDZ.VENUS/libf/phyvenus/write_histins.h
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/write_histins.h	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/write_histins.h	(revision 1301)
@@ -86,8 +86,8 @@
 c en K/s      
       call histwrite_phy(nid_ins,.false.,"dtajs",itau_w,d_t_ajs)
-c K/day ==> K/s
-      call histwrite_phy(nid_ins,.false.,"dtswr",itau_w,heat/RDAY)
-c K/day ==> K/s      
-      call histwrite_phy(nid_ins,.false.,"dtlwr",itau_w,-1.*cool/RDAY)
+c en K/s      
+      call histwrite_phy(nid_ins,.false.,"dtswr",itau_w,heat)
+c en K/s      
+      call histwrite_phy(nid_ins,.false.,"dtlwr",itau_w,-1.*cool)
 c en K/s      
 c     call histwrite_phy(nid_ins,.false.,"dtec",itau_w,d_t_ec)
Index: /trunk/LMDZ.VENUS/libf/phyvenus/write_histmth.h
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/write_histmth.h	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/write_histmth.h	(revision 1301)
@@ -86,8 +86,8 @@
 c en K/s      
       call histwrite_phy(nid_mth,.false.,"dtajs",itau_w,d_t_ajs)
-c K/day ==> K/s
-      call histwrite_phy(nid_mth,.false.,"dtswr",itau_w,heat/RDAY)
-c K/day ==> K/s      
-      call histwrite_phy(nid_mth,.false.,"dtlwr",itau_w,-1.*cool/RDAY)
+c en K/s      
+      call histwrite_phy(nid_mth,.false.,"dtswr",itau_w,heat)
+c en K/s      
+      call histwrite_phy(nid_mth,.false.,"dtlwr",itau_w,-1.*cool)
 c en K/s      
 c     call histwrite_phy(nid_mth,.false.,"dtec",itau_w,d_t_ec)
Index: /trunk/LMDZ.VENUS/libf/phyvenus/writerestartphy.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/writerestartphy.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/writerestartphy.F	(revision 1301)
@@ -2,6 +2,6 @@
      .           rlat,rlon, tsol,tsoil,
      .           albedo, 
-     .           solsw, sollw,fder,
-     .           radsol,
+     .           solsw, sollw,fder,dlw,
+     .           sollwdown,radsol,
      .    zmea, zstd, zsig, zgam, zthe, zpic, zval,
      .           t_ancien)
@@ -27,4 +27,6 @@
       real sollw(klon)
       real fder(klon)
+      real dlw(klon)
+      real sollwdown(klon)
       REAL radsol(klon)
       REAL zmea(klon), zstd(klon)
@@ -182,5 +184,5 @@
 #endif
       ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 27,
-     .                        "Rayonnement IF a la surface")
+     .                        "Rayonnement IR a la surface")
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
@@ -207,4 +209,34 @@
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "dlw", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "dlw", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 14,
+     .                        "Derivee flux IR")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,dlw)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,dlw)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "sollwdown", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "sollwdown", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 14,
+     .                        "Flux IR vers le bas a la surface")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,sollwdown)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,sollwdown)
+#endif
+c
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
       ierr = NF_DEF_VAR (nid, "RADS", NF_DOUBLE, 1, idim2,nvarid)
 #else
Index: /trunk/LMDZ.VENUS/libf/phyvenus/zenang.F
===================================================================
--- /trunk/LMDZ.VENUS/libf/phyvenus/zenang.F	(revision 1300)
+++ /trunk/LMDZ.VENUS/libf/phyvenus/zenang.F	(revision 1301)
@@ -21,6 +21,6 @@
 c lat------INPUT : latitude en degres
 c long-----INPUT : longitude en degres
-c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
-c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
+c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad/RDAY
+c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad/RDAY
 c================================================================
       use dimphy
@@ -98,9 +98,12 @@
       IF (omega1.LE.omega2) THEN  !--on est dans la meme journee locale
 c
-      IF (omega2.LE.-omega .OR. omega1.GE.omega
-     .                     .OR. omega.LT.1e-5) THEN   !--nuit
-         frac(i)=0.0
-         pmu0(i)=0.0
-      ELSE                                              !--jour+nuit/jour
+      IF (omega2.LE.-omega .OR. omega1.GE.omega      !--nuit
+     .           .OR. omega.LT.1e-5) THEN            !--nuit polaire
+        frac(i)=0.0
+        pmu0(i)=SIN(latr)*SIN(lat_sun) + 
+     .          COS(latr)*COS(lat_sun)*
+     .          (SIN(omega2)-SIN(omega1))/
+     .          (omega2-omega1)        
+      ELSE                                           !--jour+nuit/jour
         omegadeb=MAX(-omega,omega1)
         omegafin=MIN(omega,omega2)
@@ -117,5 +120,8 @@
       IF (omega1.GE.omega) THEN  !--nuit
          zfrac1=0.0
-         z1_mu =0.0
+         z1_mu =SIN(latr)*SIN(lat_sun) + 
+     .          COS(latr)*COS(lat_sun)*
+     .          (-SIN(omega1))/       
+     .          (RPI-omega1)
       ELSE                       !--jour+nuit
         omegadeb=MAX(-omega,omega1)
@@ -130,5 +136,8 @@
       IF (omega2.LE.-omega) THEN   !--nuit
          zfrac2=0.0
-         z2_mu =0.0
+         z2_mu =SIN(latr)*SIN(lat_sun) + 
+     .          COS(latr)*COS(lat_sun)*
+     .          (SIN(omega2))/       
+     .          (omega2+RPI) 
       ELSE                         !--jour+nuit
          omegadeb=-omega
@@ -143,5 +152,10 @@
 c-----------------------moyenne 
       frac(i)=(zfrac1+zfrac2)/(omega2+2*RPI-omega1)
-      pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2,1.E-10)
+      if (frac(i).ne.0.) then
+       pmu0(i)=(zfrac1*z1_mu+zfrac2*z2_mu)/MAX(zfrac1+zfrac2,1.E-10)
+      else
+       pmu0(i)=((RPI-omega1)*z1_mu+(omega2+RPI)*z2_mu)/
+     .                   (omega2+2*RPI-omega1)
+      endif
 c
       ENDIF   !---comparaison omega1 et omega2
