Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 1531)
+++ /trunk/LMDZ.MARS/README	(revision 1532)
@@ -2243,2 +2243,7 @@
 - Added "ioipsl_getin_p_mod.F90" (getin_p routine) in phy_common to
   correctly read in parameters from *.def files in a parallel environment.
+
+== 07/04/2016 == EM
+- Some fixes for buggy outputs in 1D introduced by previous code modifications.
+- made statto.h a module.
+- ecritphy in dyn3d/control_mod.F90 should be an integer, not a real.
Index: /trunk/LMDZ.MARS/libf/dyn3d/control_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/dyn3d/control_mod.F90	(revision 1531)
+++ /trunk/LMDZ.MARS/libf/dyn3d/control_mod.F90	(revision 1532)
@@ -12,5 +12,5 @@
   integer,save :: anneeref ! reference year # ! not used
   real,save :: periodav
-  real,save :: ecritphy ! output data in "diagfi.nc" every ecritphy dynamical steps 
+  integer,save :: ecritphy ! output data in "diagfi.nc" every ecritphy dynamical steps 
   integer,save :: ecritstart ! output data in "start.nc" every ecritstart dynamical steps 
   real,save :: timestart ! time start for run in "start.nc"
Index: /trunk/LMDZ.MARS/libf/phymars/inistats.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/inistats.F	(revision 1531)
+++ /trunk/LMDZ.MARS/libf/phymars/inistats.F	(revision 1532)
@@ -1,4 +1,5 @@
       subroutine inistats(ierr)
 
+      use statto_mod, only: istats,istime
       use mod_phys_lmdz_para, only : is_master
       USE comvert_mod, ONLY: ap,bp,aps,bps,preff,pseudoalt,presnivs
@@ -9,5 +10,4 @@
       implicit none
 
-      include "statto.h"
       include "netcdf.inc"
 
@@ -16,9 +16,18 @@
       integer :: l,nsteppd
       real, dimension(nbp_lev) ::  sig_s
-      real :: lon_reg_ext(nbp_lon+1) ! extended longitudes
+      real,allocatable :: lon_reg_ext(:) ! extended longitudes
       integer :: idim_lat,idim_lon,idim_llm,idim_llmp1,idim_time
       real, dimension(istime) :: lt
       integer :: nvarid
 
+
+      IF (nbp_lon*nbp_lat==1) THEN
+        ! 1D model
+        ALLOCATE(lon_reg_ext(1))
+      ELSE
+        ! 3D model
+        ALLOCATE(lon_reg_ext(nbp_lon+1))
+      ENDIF
+      
       write (*,*) 
       write (*,*) '                        || STATS ||'
@@ -46,6 +55,9 @@
       
       lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
-      !add extra redundant point (180 degrees, since lon_reg starts at -180
-      lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
+      IF (nbp_lon*nbp_lat/=1) THEN
+        ! In 3D, add extra redundant point (180 degrees,
+        ! since lon_reg starts at -180)
+        lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
+      ENDIF
 
       if (is_master) then
@@ -59,5 +71,9 @@
 
       ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_lat)
-      ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_lon)
+      IF (nbp_lon*nbp_lat==1) THEN
+        ierr = NF_DEF_DIM (nid, "longitude", 1, idim_lon)
+      ELSE
+        ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_lon)
+      ENDIF
       ierr = NF_DEF_DIM (nid, "altitude", nbp_lev, idim_llm)
       ierr = NF_DEF_DIM (nid, "llmp1", nbp_lev+1, idim_llmp1)
Index: /trunk/LMDZ.MARS/libf/phymars/iniwrite.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/iniwrite.F	(revision 1531)
+++ /trunk/LMDZ.MARS/libf/phymars/iniwrite.F	(revision 1532)
@@ -1,3 +1,3 @@
-      SUBROUTINE iniwrite(nid,idayref,phis,area)
+      SUBROUTINE iniwrite(nid,idayref,phis,area,nbplon,nbplat)
 
       use comsoil_h, only: mlayer, nsoilmx
@@ -35,6 +35,7 @@
       integer,intent(in) :: nid        ! NetCDF file ID
       INTEGER*4,intent(in) :: idayref  ! date (initial date for this run)
-      real,intent(in) :: phis(nbp_lon+1,nbp_lat) ! surface geopotential
-      real,intent(in) :: area(nbp_lon+1,nbp_lat) ! mesh area (m2)
+      real,intent(in) :: phis(nbplon,nbp_lat) ! surface geopotential
+      real,intent(in) :: area(nbplon,nbp_lat) ! mesh area (m2)
+      integer,intent(in) :: nbplon,nbplat ! sizes of area and phis arrays
 
 c   Local:
@@ -44,5 +45,5 @@
       REAL tab_cntrl(length) ! run parameters are stored in this array
       INTEGER ierr
-      REAl :: lon_reg_ext(nbp_lon+1) ! extended longitudes
+      REAl,ALLOCATABLE :: lon_reg_ext(:) ! extended longitudes
 
       integer :: nvarid,idim_index,idim_rlonv
@@ -51,4 +52,12 @@
       integer, dimension(2) :: id  
 c-----------------------------------------------------------------------
+
+      IF (nbp_lon*nbp_lat==1) THEN
+        ! 1D model
+        ALLOCATE(lon_reg_ext(1))
+      ELSE
+        ! 3D model
+        ALLOCATE(lon_reg_ext(nbp_lon+1))
+      ENDIF
 
       DO l=1,length
@@ -104,5 +113,9 @@
 !      ierr = NF_DEF_DIM (nid, "rlonu", iip1, idim_rlonu)
       ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_rlatu)
-      ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_rlonv)
+      IF (nbp_lon*nbp_lat==1) THEN
+        ierr = NF_DEF_DIM (nid, "longitude", 1, idim_rlonv)
+      ELSE
+        ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_rlonv)
+      ENDIF
 !      ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv)
       ierr = NF_DEF_DIM (nid, "interlayer", (nbp_lev+1), idim_llmp1)
@@ -166,7 +179,11 @@
 c
 c --------------------------
+      
       lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
-      !add extra redundant point (180 degrees, since lon_reg starts at -180
-      lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
+      IF (nbp_lon*nbp_lat/=1) THEN
+        ! In 3D, add extra redundant point (180 degrees,
+        ! since lon_reg starts at -180)
+        lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
+      ENDIF
 
       ierr = NF_REDEF (nid)
Index: /trunk/LMDZ.MARS/libf/phymars/iniwritesoil.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/iniwritesoil.F90	(revision 1531)
+++ /trunk/LMDZ.MARS/libf/phymars/iniwritesoil.F90	(revision 1532)
@@ -1,3 +1,3 @@
-subroutine iniwritesoil(nid,ngrid,inertia,area)
+subroutine iniwritesoil(nid,ngrid,inertia,area,nbplon,nbplat)
 
 ! initialization routine for 'writediagoil'. Here we create/define
@@ -17,6 +17,7 @@
 integer,intent(in) :: ngrid
 integer,intent(in) :: nid ! NetCDF output file ID
-real,intent(in) :: inertia(nbp_lon+1,nbp_lat,nsoilmx)
-real,intent(in) :: area(nbp_lon+1,nbp_lat) ! mesh area (m2)
+real,intent(in) :: inertia(nbplon,nbplat,nsoilmx)
+real,intent(in) :: area(nbplon,nbp_lat) ! mesh area (m2)
+integer,intent(in) :: nbplon,nbplat ! sizes of area
 
 ! Local variables:
@@ -33,5 +34,14 @@
 real,dimension(nbp_lon+1,nbp_lat,nsoilmx) :: data3 ! to store 3D data
 integer :: i,j,l,ig0
-real :: lon_reg_ext(nbp_lon+1) ! extended longitudes
+real,allocatable :: lon_reg_ext(:) ! extended longitudes
+
+
+if (nbp_lon*nbp_lat==1) then
+  ! 1D model
+  allocate(lon_reg_ext(1))
+else
+  ! 3D model
+  allocate(lon_reg_ext(nbp_lon+1))
+endif
 
 ! 1. Define the dimensions
@@ -40,5 +50,9 @@
 
 ! Define the dimensions
-ierr=NF_DEF_DIM(nid,"longitude",nbp_lon+1,idim_rlonv)
+if (nbp_lon*nbp_lat==1) then
+  ierr=NF_DEF_DIM(nid,"longitude",1,idim_rlonv)
+else
+  ierr=NF_DEF_DIM(nid,"longitude",nbp_lon+1,idim_rlonv)
+endif
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not define longitude dimension"
@@ -82,6 +96,8 @@
 
 lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
-!add extra redundant point (180 degrees, since lon_reg starts at -180
-lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
+if (nbp_lon*nbp_lat/=1) then ! in 3D only:
+  ! add extra redundant point (180 degrees, since lon_reg starts at -180
+  lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
+endif
 
 ! Write longitude to file
Index: /trunk/LMDZ.MARS/libf/phymars/mkstat.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/mkstat.F90	(revision 1531)
+++ /trunk/LMDZ.MARS/libf/phymars/mkstat.F90	(revision 1532)
@@ -10,10 +10,10 @@
 !  Yann W. july 2003
 
+use statto_mod, only: istime,count
 use mod_phys_lmdz_para, only : is_master
-use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev
+use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev, klon_glo
 
 implicit none
 
-include "statto.h"
 include "netcdf.inc"
 
@@ -22,9 +22,9 @@
 integer, dimension(5) :: dimids
 character (len=50) :: name,nameout,units,title
-real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: sum3d,square3d,mean3d,sd3d
-real, dimension(nbp_lon+1,nbp_lat) :: sum2d,square2d,mean2d,sd2d
+real,allocatable :: sum3d(:,:,:),square3d(:,:,:),mean3d(:,:,:),sd3d(:,:,:)
+real,allocatable :: sum2d(:,:),square2d(:,:),mean2d(:,:),sd2d(:,:)
 real, dimension(istime) :: time
 real, dimension(nbp_lat) :: lat
-real, dimension(nbp_lon+1) :: lon
+real,allocatable :: lon(:)
 real, dimension(nbp_lev) :: alt
 logical :: lcopy=.true.
@@ -38,4 +38,25 @@
 if (is_master) then
 ! only the master needs do this
+if (klon_glo>1) then
+  allocate(lon(nbp_lon+1))
+  allocate(sum3d(nbp_lon+1,nbp_lat,nbp_lev))
+  allocate(square3d(nbp_lon+1,nbp_lat,nbp_lev))
+  allocate(mean3d(nbp_lon+1,nbp_lat,nbp_lev))
+  allocate(sd3d(nbp_lon+1,nbp_lat,nbp_lev))
+  allocate(sum2d(nbp_lon+1,nbp_lat))
+  allocate(square2d(nbp_lon+1,nbp_lat))
+  allocate(mean2d(nbp_lon+1,nbp_lat))
+  allocate(sd2d(nbp_lon+1,nbp_lat))
+else ! 1D model case
+  allocate(lon(1))
+  allocate(sum3d(1,1,nbp_lev))
+  allocate(square3d(1,1,nbp_lev))
+  allocate(mean3d(1,1,nbp_lev))
+  allocate(sd3d(1,1,nbp_lev))
+  allocate(sum2d(1,1))
+  allocate(square2d(1,1))
+  allocate(mean2d(1,1))
+  allocate(sd2d(1,1))
+endif
 
 ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
@@ -105,5 +126,9 @@
 !      dimout(4)=timeid
 
-      size=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
+      if (klon_glo>1) then ! general case
+        size=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
+      else ! 1D model
+        size=(/1,1,nbp_lev,1/)
+      endif
       do lt=1,istime
          start=(/1,1,1,lt/)
@@ -135,5 +160,9 @@
 !      dimout(3)=timeid
 
-      size=(/nbp_lon+1,nbp_lat,1,0/)
+      if (klon_glo>1) then ! general case
+        size=(/nbp_lon+1,nbp_lat,1,0/)
+      else
+        size=(/1,1,1,0/)
+      endif
       do lt=1,istime
          start=(/1,1,lt,0/)
Index: unk/LMDZ.MARS/libf/phymars/statto.h
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/statto.h	(revision 1531)
+++ 	(revision )
@@ -1,34 +1,0 @@
-
-!  statto:
-!     This include file controls the production of statistics.
-!     Some variables could be set in a namelist, but it is easier to
-!     do it here since arrays can then be dimensioned using parameters
-!     and values shouldn't have to change too often.   SRL
-
-!     Calculate stats every istats physics timesteps, starting at first
-!     call.  If istats=0 then don't do statistics at all.  Check value
-!     if number of physics timesteps changes.
-	integer istats
-
-!     Calculate itime independent sums and sums of squares,
-!     example, istat=1,istime=1 gives a single time mean
-	integer, parameter :: istime=12
-
-!     Number of 2D and 3D variables on which to do statistics.
-	integer n2dvar, n3dvar
-	parameter (n2dvar = 8, n3dvar = 5)
-
-!     Units for writing stats header and data
-	integer usdata
-
-!     count tab to know the variable record
-        integer count(istime)
-
-!     Record of the number of stores made for each time.
-	integer nstore(istime)
-
-! Size of the "controle" array
-        integer, parameter :: cntrlsize=15
-
-!       common /sttcom/ dummy,nstore,istats,usdata
-        common /sttcom/ nstore,istats,usdata,count
Index: /trunk/LMDZ.MARS/libf/phymars/statto_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/statto_mod.F90	(revision 1532)
+++ /trunk/LMDZ.MARS/libf/phymars/statto_mod.F90	(revision 1532)
@@ -0,0 +1,35 @@
+MODULE statto_mod
+IMPLICIT NONE  
+
+!  statto_mod:
+!     This include file controls the production of statistics.
+!     Some variables could be set in a namelist, but it is easier to
+!     do it here since arrays can then be dimensioned using parameters
+!     and values shouldn't have to change too often.   SRL
+
+!     Calculate stats every istats physics timesteps, starting at first
+!     call.  If istats=0 then don't do statistics at all.  Check value
+!     if number of physics timesteps changes.
+	integer istats
+
+!     Calculate itime independent sums and sums of squares,
+!     example, istat=1,istime=1 gives a single time mean
+	integer, parameter :: istime=12
+
+!     Number of 2D and 3D variables on which to do statistics.
+	integer n2dvar, n3dvar
+	parameter (n2dvar = 8, n3dvar = 5)
+
+!     Units for writing stats header and data
+	integer usdata
+
+!     count tab to know the variable record
+        integer count(istime)
+
+!     Record of the number of stores made for each time.
+	integer nstore(istime)
+
+! Size of the "controle" array
+        integer, parameter :: cntrlsize=15
+
+END MODULE statto_mod
Index: /trunk/LMDZ.MARS/libf/phymars/writediagfi.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/writediagfi.F	(revision 1531)
+++ /trunk/LMDZ.MARS/libf/phymars/writediagfi.F	(revision 1532)
@@ -63,6 +63,9 @@
       real*4 dx1(nbp_lev)           ! to store a 1D (column) data set
       real*4 dx0
+      real*4 dx3_1d(1,nbp_lev) ! to store a profile with 1D model
+      real*4 dx2_1d ! to store a surface value with 1D model
 
       real*4,save :: date
+!$OMP THREADPRIVATE(date)
 
       REAL phis((nbp_lon+1),nbp_lat)
@@ -75,7 +78,9 @@
       integer,save :: zitau=0
       character(len=20),save :: firstnom='1234567890'
+!$OMP THREADPRIVATE(zitau,firstnom)
 
 ! Ajouts
       integer, save :: ntime=0
+!$OMP THREADPRIVATE(ntime)
       integer :: idim,varid
       integer :: nid
@@ -92,7 +97,6 @@
       character(len=120),save :: nom_def(n_nom_def_max)
       logical,save :: firstcall=.true.
+!$OMP THREADPRIVATE(firstcall) 	!diagfi_def,n_nom_def,nom_def read in diagfi.def
       
-#ifndef MESOSCALE
-
 #ifdef CPP_PARA
 ! Added to work in parallel mode
@@ -127,4 +131,5 @@
          firstcall=.false.
 
+!$OMP MASTER
   !      Open diagfi.def definition file if there is one:
          open(99,file="diagfi.def",status='old',form='formatted',
@@ -150,4 +155,6 @@
             diagfi_def=.false.
          endif
+!$OMP END MASTER
+!$OMP BARRIER
       END IF ! of IF (firstcall)
 
@@ -214,5 +221,6 @@
 
          ! Build phis() and area()
-         do i=1,nbp_lon+1 ! poles
+         IF (klon_glo>1) THEN
+          do i=1,nbp_lon+1 ! poles
            phis(i,1)=phisfi_glo(1)
            phis(i,nbp_lat)=phisfi_glo(klon_glo)
@@ -220,6 +228,6 @@
            area(i,1)=areafi_glo(1)/nbp_lon
            area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
-         enddo
-         do j=2,nbp_lat-1
+          enddo
+          do j=2,nbp_lat-1
            ig0= 1+(j-2)*nbp_lon
            do i=1,nbp_lon
@@ -230,8 +238,13 @@
            phis(nbp_lon+1,j)=phis(1,j)
            area(nbp_lon+1,j)=area(1,j)
-         enddo
+          enddo
+         ENDIF
          
          ! write "header" of file (longitudes, latitudes, geopotential, ...)
-         call iniwrite(nid,day_ini,phis,area)
+         IF (klon_glo>1) THEN ! general 3D case
+           call iniwrite(nid,day_ini,phis,area,nbp_lon+1,nbp_lat)
+         ELSE ! 1D model
+           call iniwrite(nid,day_ini,phisfi_glo(1),areafi_glo(1),1,1)
+         ENDIF
 
          endif ! of if (is_master)
@@ -248,10 +261,9 @@
       endif ! if (firstnom.eq.'1234567890')
 
-      if (ngrid.eq.1) then
+      if (klon_glo.eq.1) then
         ! in testphys1d, for the 1d version of the GCM, iphysiq and irythme
         ! are undefined; so set them to 1
         iphysiq=1 
         irythme=1
-        ! NB: 
       endif
 
@@ -320,4 +332,5 @@
 !         Passage variable physique -->  variable dynamique
 !         recast (copy) variable from physics grid to dynamics grid
+          IF (klon_glo>1) THEN ! General case
            DO l=1,nbp_lev
              DO i=1,nbp_lon+1
@@ -333,4 +346,7 @@
              ENDDO
            ENDDO
+          ELSE ! 1D model case
+           dx3_1d(1,1:nbp_lev)=px(1,1:nbp_lev)
+          ENDIF
 #endif
 !         Ecriture du champs
@@ -360,5 +376,9 @@
            corner(4)=ntime
 
-           edges(1)=nbp_lon+1
+           IF (klon_glo==1) THEN
+             edges(1)=1
+           ELSE
+             edges(1)=nbp_lon+1
+           ENDIF
            edges(2)=nbp_lat
            edges(3)=nbp_lev
@@ -371,12 +391,16 @@
 !           write(*,*)"       edges()=",edges
 !           write(*,*)"       dx3()=",dx3
-           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
+           IF (klon_glo>1) THEN ! General case
+             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
+           ELSE
+             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3_1d)
+           ENDIF
 !#endif
 
            if (ierr.ne.NF_NOERR) then
               write(*,*) "***** PUT_VAR problem in writediagfi"
-              write(*,*) "***** with ",nom
+              write(*,*) "***** with dx3: ",nom
               write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
-c             call abort
+              stop
            endif 
 
@@ -405,5 +429,5 @@
 !         Passage variable physique -->  physique dynamique
 !         recast (copy) variable from physics grid to dynamics grid
-
+          IF (klon_glo>1) THEN ! General case
              DO i=1,nbp_lon+1
                 dx2(i,1)=px(1,1)
@@ -417,4 +441,7 @@
                 dx2(nbp_lon+1,j)=dx2(1,j)
              ENDDO
+          ELSE ! 1D model case
+            dx2_1d=px(1,1)
+          ENDIF
 #endif
 
@@ -442,5 +469,9 @@
            corner(2)=1
            corner(3)=ntime
-           edges(1)=nbp_lon+1
+           IF (klon_glo==1) THEN
+             edges(1)=1
+           ELSE
+             edges(1)=nbp_lon+1
+           ENDIF
            edges(2)=nbp_lat
            edges(3)=1
@@ -450,12 +481,16 @@
 !           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx2) 
 !#else         
-           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2)
+           IF (klon_glo>1) THEN ! General case
+             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2)
+           ELSE
+             ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2_1d)
+           ENDIF
 !#endif     
 
            if (ierr.ne.NF_NOERR) then
               write(*,*) "***** PUT_VAR matter in writediagfi"
-              write(*,*) "***** with ",nom
+              write(*,*) "***** with dx2: ",nom
               write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
-c             call abort
+              stop
            endif 
 
@@ -505,7 +540,7 @@
            if (ierr.ne.NF_NOERR) then
               write(*,*) "***** PUT_VAR problem in writediagfi"
-              write(*,*) "***** with ",nom
+              write(*,*) "***** with dx1: ",nom
               write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
-c             call abort
+              stop
            endif 
 
@@ -543,7 +578,7 @@
            if (ierr.ne.NF_NOERR) then
               write(*,*) "***** PUT_VAR matter in writediagfi"
-              write(*,*) "***** with ",nom
+              write(*,*) "***** with dx0: ",nom
               write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
-c             call abort
+              stop
            endif 
 
@@ -558,4 +593,3 @@
       endif
 
-#endif
       end
Index: /trunk/LMDZ.MARS/libf/phymars/writediagsoil.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/writediagsoil.F90	(revision 1531)
+++ /trunk/LMDZ.MARS/libf/phymars/writediagsoil.F90	(revision 1532)
@@ -36,4 +36,6 @@
 real*4,dimension(nbp_lon+1,nbp_lat) :: data2 ! to store 2D data
 real*4 :: data0 ! to store 0D data
+real*4 :: data3_1d(1,nsoilmx) ! to store a profile in 1D model
+real*4 :: data2_1d ! to store surface value with 1D model
 integer :: i,j,l ! for loops
 integer :: ig0
@@ -51,4 +53,5 @@
 character(len=20),save :: firstname="1234567890"
 integer,save :: zitau=0
+!$OMP THREADPRIVATE(date,isample,ntime,firstname,zitau)
 
 character(len=30) :: filename="diagsoil.nc"
@@ -107,5 +110,6 @@
 
    ! build inertia() and area()
-   do i=1,nbp_lon+1 ! poles
+   if (klon_glo>1) then
+    do i=1,nbp_lon+1 ! poles
      inertia(i,1,1:nsoilmx)=inertiafi_glo(1,1:nsoilmx)
      inertia(i,nbp_lat,1:nsoilmx)=inertiafi_glo(klon_glo,1:nsoilmx)
@@ -113,6 +117,6 @@
      area(i,1)=areafi_glo(1)/nbp_lon
      area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
-   enddo
-   do j=2,nbp_lat-1
+    enddo
+    do j=2,nbp_lat-1
      ig0= 1+(j-2)*nbp_lon
      do i=1,nbp_lon
@@ -123,8 +127,13 @@
      inertia(nbp_lon+1,j,1:nsoilmx)=inertia(1,j,1:nsoilmx)
      area(nbp_lon+1,j)=area(1,j)
-   enddo
+    enddo
+   endif
    
    ! write "header" of file (longitudes, latitudes, geopotential, ...)
-   call iniwritesoil(nid,ngrid,inertia,area)
+   if (klon_glo>1) then ! general 3D case
+     call iniwritesoil(nid,ngrid,inertia,area,nbp_lon+1,nbp_lat)
+   else ! 1D model
+     call iniwritesoil(nid,ngrid,inertiafi_glo(1,:),areafi_glo(1),1,1)
+   endif
 
   endif ! of if (is_master)
@@ -188,5 +197,6 @@
 !$OMP BARRIER
 #else
-  do l=1,nsoilmx
+  if (klon_glo>1) then ! General case
+   do l=1,nsoilmx
     ! handle the poles
     do i=1,nbp_lon+1
@@ -202,5 +212,8 @@
       data3(nbp_lon+1,j,l)=data3(1,j,l) ! extra (modulo) longitude
     enddo
-  enddo
+   enddo
+  else ! 1D model case
+   data3_1d(1,1:nsoilmx)=px(1,1:nsoilmx)
+  endif
 #endif
   
@@ -229,5 +242,9 @@
   corners(4)=ntime
   
-  edges(1)=nbp_lon+1
+  if (klon_glo==1) then
+    edges(1)=1
+  else
+    edges(1)=nbp_lon+1
+  endif
   edges(2)=nbp_lat
   edges(3)=nsoilmx
@@ -238,5 +255,9 @@
 !  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data3)
 !#else
-  ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data3)
+  if (klon_glo>1) then
+    ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data3)
+  else
+    ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data3_1d)
+  endif
 !#endif
   if (ierr.ne.NF_NOERR) then
@@ -263,17 +284,21 @@
 !$OMP BARRIER
 #else
-  ! handle the poles
-  do i=1,nbp_lon+1
-    data2(i,1)=px(1,1)
-    data2(i,nbp_lat)=px(ngrid,1)
-  enddo
-  ! rest of the grid
-  do j=2,nbp_lat-1
-    ig0=1+(j-2)*nbp_lon
-    do i=1,nbp_lon
-      data2(i,j)=px(ig0+i,1)
-    enddo
-    data2(nbp_lon+1,j)=data2(1,j) ! extra (modulo) longitude
-  enddo
+  if (klon_glo>1) then ! general case
+    ! handle the poles
+    do i=1,nbp_lon+1
+      data2(i,1)=px(1,1)
+      data2(i,nbp_lat)=px(ngrid,1)
+    enddo
+    ! rest of the grid
+    do j=2,nbp_lat-1
+      ig0=1+(j-2)*nbp_lon
+      do i=1,nbp_lon
+        data2(i,j)=px(ig0+i,1)
+      enddo
+      data2(nbp_lon+1,j)=data2(1,j) ! extra (modulo) longitude
+    enddo
+  else ! 1D model case
+    data2_1d=px(1,1)
+  endif
 #endif
 
@@ -300,5 +325,9 @@
   corners(3)=ntime
   
-  edges(1)=nbp_lon+1
+  if (klon_glo==1) then
+    edges(1)=1
+  else
+    edges(1)=nbp_lon+1
+  endif
   edges(2)=nbp_lat
   edges(3)=1
@@ -308,5 +337,9 @@
 !  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data2)
 !#else
-  ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data2)
+  if (klon_glo>1) then ! General case
+    ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data2)
+  else
+    ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data2_1d)
+  endif
 !#endif
   if (ierr.ne.NF_NOERR) then
Index: /trunk/LMDZ.MARS/libf/phymars/wstats.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/wstats.F90	(revision 1531)
+++ /trunk/LMDZ.MARS/libf/phymars/wstats.F90	(revision 1532)
@@ -1,4 +1,5 @@
 subroutine wstats(ngrid,nom,titre,unite,dim,px)
 
+use statto_mod, only: istats,istime,count
 use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather, klon_mpi_begin
 use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo, &
@@ -6,5 +7,4 @@
 implicit none
 
-include "statto.h"
 include "netcdf.inc"
 
@@ -13,8 +13,9 @@
 integer,intent(in) :: dim
 real,intent(in) :: px(ngrid,nbp_lev)
-real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: mean3d,sd3d,dx3
-real, dimension(nbp_lon+1,nbp_lat) :: mean2d,sd2d,dx2
+real,allocatable,save :: mean3d(:,:,:),sd3d(:,:,:),dx3(:,:,:)
+real,allocatable,save :: mean2d(:,:),sd2d(:,:),dx2(:,:)
 character (len=50) :: namebis
 character (len=50), save :: firstvar
+!$OMP THREADPRIVATE(firstvar)
 integer :: ierr,varid,nbdim,nid
 integer :: meanid,sdid
@@ -25,4 +26,5 @@
 
 integer, save :: step=0
+!$OMP THREADPRIVATE(firstcall,indx,step)
 
 ! Added to work in parallel mode
@@ -47,4 +49,19 @@
    firstvar=trim((nom))
    call inistats(ierr)
+   if (klon_glo>1) then ! general case, 3D GCM
+     allocate(mean3d(nbp_lon+1,nbp_lat,nbp_lev))
+     allocate(sd3d(nbp_lon+1,nbp_lat,nbp_lev))
+     allocate(dx3(nbp_lon+1,nbp_lat,nbp_lev))
+     allocate(mean2d(nbp_lon+1,nbp_lat))
+     allocate(sd2d(nbp_lon+1,nbp_lat))
+     allocate(dx2(nbp_lon+1,nbp_lat))
+   else ! 1D model
+     allocate(mean3d(1,1,nbp_lev))
+     allocate(sd3d(1,1,nbp_lev))
+     allocate(dx3(1,1,nbp_lev))
+     allocate(mean2d(1,1))
+     allocate(sd2d(1,1))
+     allocate(dx2(1,1))
+   endif
 endif
 
@@ -192,10 +209,18 @@
    if (dim.eq.3) then
       start=(/1,1,1,indx/)
-      sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
+      if (klon_glo>1) then !general case
+        sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
+      else
+        sizes=(/1,1,nbp_lev,1/)
+      endif
       mean3d(:,:,:)=0
       sd3d(:,:,:)=0
    else if (dim.eq.2) then
       start=(/1,1,indx,0/)
-      sizes=(/nbp_lon+1,nbp_lev,1,0/)
+      if (klon_glo>1) then !general case
+        sizes=(/nbp_lon+1,nbp_lev,1,0/)
+      else
+        sizes=(/1,1,1,0/)
+      endif
       mean2d(:,:)=0
       sd2d(:,:)=0
@@ -205,5 +230,9 @@
    if (dim.eq.3) then
       start=(/1,1,1,indx/)
-      sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
+      if (klon_glo>1) then !general case
+        sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
+      else ! 1D model case
+        sizes=(/1,1,nbp_lev,1/)
+      endif
 #ifdef NC_DOUBLE
       ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
@@ -214,4 +243,5 @@
 #endif
       if (ierr.ne.NF_NOERR) then
+         write (*,*) "wstats error reading :",trim(nom)
          write (*,*) NF_STRERROR(ierr)
          stop ""
@@ -220,5 +250,9 @@
    else if (dim.eq.2) then
       start=(/1,1,indx,0/)
-      sizes=(/nbp_lon+1,nbp_lat,1,0/)
+      if (klon_glo>1) then ! general case
+        sizes=(/nbp_lon+1,nbp_lat,1,0/)
+      else
+        sizes=(/1,1,1,0/)
+      endif
 #ifdef NC_DOUBLE
       ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
@@ -229,4 +263,5 @@
 #endif
       if (ierr.ne.NF_NOERR) then
+         write (*,*) "wstats error reading :",trim(nom)
          write (*,*) NF_STRERROR(ierr)
          stop ""
@@ -240,8 +275,12 @@
 if (dim.eq.3) then
   dx3(1:nbp_lon,:,:)=px3_glo(:,:,:)
-  dx3(nbp_lon+1,:,:)=dx3(1,:,:)
+  IF (klon_glo>1) THEN ! in 3D, add redundant longitude point
+    dx3(nbp_lon+1,:,:)=dx3(1,:,:)
+  ENDIF
 else ! dim.eq.2
   dx2(1:nbp_lon,:)=px2_glo(:,:)
-  dx2(nbp_lon+1,:)=dx2(1,:)
+  IF (klon_glo>1) THEN ! in 3D, add redundant longitude point
+    dx2(nbp_lon+1,:)=dx2(1,:)
+  ENDIF
 endif
 
@@ -261,4 +300,9 @@
    ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd3d)
 #endif
+  if (ierr.ne.NF_NOERR) then
+     write (*,*) "wstats error writing :",trim(nom)
+     write (*,*) NF_STRERROR(ierr)
+     stop ""
+  endif
 
 else if (dim.eq.2) then
@@ -274,4 +318,13 @@
    ierr = NF_PUT_VARA_REAL(nid,sdid,start,sizes,sd2d)
 #endif
+  if (ierr.ne.NF_NOERR) then
+     write (*,*) "wstats error writing :",trim(nom)
+     write(*,*) "start:",start
+     write(*,*) "sizes:",sizes
+     write(*,*) "mean2d:",mean2d
+     write(*,*) "sd2d:",sd2d
+     write (*,*) NF_STRERROR(ierr)
+     stop ""
+  endif
 
 endif ! of if (dim.eq.3) elseif (dim.eq.2)
@@ -284,5 +337,5 @@
 !======================================================
 subroutine inivar(nid,varid,ngrid,dim,indx,px,ierr)
-use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev
+use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev, klon_glo
 
 implicit none
@@ -298,13 +351,20 @@
 real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: dx3
 real, dimension(nbp_lon+1,nbp_lat) :: dx2
+real :: dx3_1d(nbp_lev) ! for 1D outputs
+real :: dx2_1d ! for 1D outputs
 
 if (dim.eq.3) then
 
    start=(/1,1,1,indx/)
-   sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
+   if (klon_glo>1) then ! general 3D case
+     sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
+   else
+     sizes=(/1,1,nbp_lev,1/)
+   endif
 
 !  Passage variable physique -->  variable dynamique
 
-   DO l=1,nbp_lev
+   if (klon_glo>1) then ! general case
+    DO l=1,nbp_lev
       DO i=1,nbp_lon+1
          dx3(i,1,l)=px(1,l)
@@ -318,19 +378,40 @@
          dx3(nbp_lon+1,j,l)=dx3(1,j,l)
       ENDDO
-   ENDDO
-
-#ifdef NC_DOUBLE
-   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3)
-#else
-   ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3)
-#endif
+    ENDDO
+   else ! 1D model case
+     dx3_1d(1:nbp_lev)=px(1,1:nbp_lev)
+   endif
+
+#ifdef NC_DOUBLE
+   if (klon_glo>1) then
+     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3)
+   else
+     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx3_1d)
+   endif
+#else
+   if (klon_glo>1) then
+     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3)
+   else
+     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx3_1d)
+   endif
+#endif
+  if (ierr.ne.NF_NOERR) then
+     write (*,*) "inivar error writing variable"
+     write (*,*) NF_STRERROR(ierr)
+     stop ""
+  endif
 
 else if (dim.eq.2) then
 
       start=(/1,1,indx,0/)
-      sizes=(/nbp_lon+1,nbp_lat,1,0/)
+      if (klon_glo>1) then ! general 3D case
+        sizes=(/nbp_lon+1,nbp_lat,1,0/)
+      else
+        sizes=(/1,1,1,0/)
+      endif
 
 !    Passage variable physique -->  physique dynamique
 
+  if (klon_glo>1) then ! general case
   DO i=1,nbp_lon+1
      dx2(i,1)=px(1,1)
@@ -344,10 +425,26 @@
      dx2(nbp_lon+1,j)=dx2(1,j)
   ENDDO
-
-#ifdef NC_DOUBLE
-   ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2)
-#else
-   ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2)
-#endif
+  else ! 1D model case
+    dx2_1d=px(1,1)
+  endif
+  
+#ifdef NC_DOUBLE
+   if (klon_glo>1) then
+     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2)
+   else
+     ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,sizes,dx2_1d)
+   endif
+#else
+   if (klon_glo>1) then
+     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2)
+   else
+     ierr = NF_PUT_VARA_REAL(nid,varid,start,sizes,dx2_1d)
+   endif
+#endif
+  if (ierr.ne.NF_NOERR) then
+     write (*,*) "inivar error writing variable"
+     write (*,*) NF_STRERROR(ierr)
+     stop ""
+  endif
 
 endif
