Index: trunk/libf/phylmd/physiq.F
===================================================================
--- trunk/libf/phylmd/physiq.F	(revision 1)
+++ trunk/libf/phylmd/physiq.F	(revision 66)
@@ -1,3 +1,3 @@
-! $Id: physiq.F 1428 2010-09-13 08:43:37Z fairhead $
+! $Id: physiq.F 1486 2011-02-11 12:07:39Z fairhead $
 c#define IO_DEBUG
 
@@ -158,4 +158,8 @@
       save ok_LES                            
 c$OMP THREADPRIVATE(ok_LES)                  
+c
+      LOGICAL callstats ! sortir le fichier stats 
+      save callstats                            
+c$OMP THREADPRIVATE(callstats)                  
 c
       LOGICAL ok_region ! sortir le fichier regional
@@ -606,5 +610,5 @@
       real, save :: ale_bl_prescr=0.
 
-      real, save :: ale_max=100.
+      real, save :: ale_max=1000.
       real, save :: alp_max=2.
 
@@ -1150,4 +1154,5 @@
 !     and 360
 
+      INTEGER ierr
 #include "YOMCST.h"
 #include "YOETHF.h"
@@ -1222,4 +1227,5 @@
      .     ok_instan, ok_hf,
      .     ok_LES,
+     .     callstats,
      .     solarlong0,seuil_inversion,
      .     fact_cldcon, facttemps,ok_newmicro,iflag_radia,
@@ -1250,5 +1256,4 @@
 cym Attention pbase pas initialise dans concvl !!!!
           pbase=0
-          paire_ter(:)=0.    
 cIM 180608
 c         pmflxr=0.
@@ -1847,21 +1852,22 @@
 ! doit donc etre placé avant radlwsw et pbl_surface
 
+!!!   jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
+      day_since_equinox = (jD_cur + jH_cur) - jD_eq 
+!
+!   choix entre calcul de la longitude solaire vraie ou valeur fixee a 
+!   solarlong0
+      if (solarlong0<-999.) then
+       if (new_orbit) then
 ! calcul selon la routine utilisee pour les planetes
-      if (new_orbit) then
-        call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
-        day_since_equinox = (jD_cur + jH_cur) - jD_eq 
-!        day_since_equinox = (jD_cur) - jD_eq 
         call solarlong(day_since_equinox, zlongi, dist)
-      else      
+       else    
 ! calcul selon la routine utilisee pour l'AR4
-!   choix entre calcul de la longitude solaire vraie ou valeur fixee a
-!   solarlong0
-        if (solarlong0<-999.) then
-           CALL orbite(REAL(days_elapsed+1),zlongi,dist)
-        else
+        CALL orbite(REAL(days_elapsed+1),zlongi,dist)
+       endif   
+      else    
            zlongi=solarlong0  ! longitude solaire vraie
-           dist=1.            ! distance au soleil / moyenne
-        endif
-      endif
+           dist=1.            ! distance au soleil / moyenne 
+      endif   
       if(prt_level.ge.1)                                                &
      &    write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
@@ -3376,5 +3382,6 @@
      I                   cdragh,coefh,u1,v1,ftsol,pctsrf,
      I                   frac_impa, frac_nucl,
-     I                   pphis,airephy,dtime,itap)
+     I                   pphis,airephy,dtime,itap,
+     I                   rlon,rlat,qx(:,:,ivap),da,phi,mp,upwd,dnwd)
 
 
@@ -3694,5 +3701,58 @@
 c====================================================================
 c
-      
+
+c        -----------------------------------------------------------------
+c        WSTATS: Saving statistics
+c        -----------------------------------------------------------------
+c        ("stats" stores and accumulates 8 key variables in file "stats.nc"
+c        which can later be used to make the statistic files of the run:
+c        "stats")          only possible in 3D runs !
+
+         
+         IF (callstats) THEN
+
+           call wstats(klon,o_psol%name,"Surface pressure","Pa"
+     &                 ,2,paprs(:,1))
+           call wstats(klon,o_tsol%name,"Surface temperature","K",
+     &                 2,zxtsol)
+           zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
+           call wstats(klon,o_precip%name,"Precip Totale liq+sol",
+     &                 "kg/(s*m2)",2,zx_tmp_fi2d)
+           zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
+           call wstats(klon,o_plul%name,"Large-scale Precip",
+     &                 "kg/(s*m2)",2,zx_tmp_fi2d)
+           zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
+           call wstats(klon,o_pluc%name,"Convective Precip",
+     &                 "kg/(s*m2)",2,zx_tmp_fi2d)
+           call wstats(klon,o_sols%name,"Solar rad. at surf.",
+     &                 "W/m2",2,solsw)
+           call wstats(klon,o_soll%name,"IR rad. at surf.",
+     &                 "W/m2",2,sollw)
+          zx_tmp_fi2d(:) = topsw(:)-toplw(:)
+          call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA",
+     &                 "W/m2",2,zx_tmp_fi2d)
+
+
+
+           call wstats(klon,o_temp%name,"Air temperature","K",
+     &                 3,t_seri)
+           call wstats(klon,o_vitu%name,"Zonal wind","m.s-1",
+     &                 3,u_seri)
+           call wstats(klon,o_vitv%name,"Meridional wind",
+     &                "m.s-1",3,v_seri)
+           call wstats(klon,o_vitw%name,"Vertical wind",
+     &                "m.s-1",3,omega)
+           call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg",
+     &                 3,q_seri)
+ 
+
+
+           IF(lafin) THEN
+             write (*,*) "Writing stats..."
+             call mkstats(ierr)
+           ENDIF
+
+         ENDIF !if callstats
+     
 
       IF (lafin) THEN
