Index: trunk/LMDZ.MARS/libf/phymars/writediagsoil.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/writediagsoil.F90	(revision 1525)
+++ trunk/LMDZ.MARS/libf/phymars/writediagsoil.F90	(revision 1528)
@@ -12,16 +12,14 @@
 ! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4
 
-use comsoil_h, only: nsoilmx
+use comsoil_h, only: nsoilmx, inertiedat
+use comgeomphy, only: airephy
 use time_phylmdz_mod, only: ecritphy, day_step, iphysiq
 use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
-use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
+use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo, &
+                              nbp_lon, nbp_lat
 
 implicit none
 
-#include"dimensions.h"
-#include"paramet.h"
-!#include"control.h"
-!#include"comsoil.h"
-#include"netcdf.inc"
+include"netcdf.inc"
 
 ! Arguments:
@@ -33,10 +31,8 @@
 integer,intent(in) :: dimpx ! dimension of the variable (3,2 or 0)
 real,dimension(ngrid,nsoilmx),intent(in) :: px ! variable
-! Note: nsoilmx is a parameter set in 'comsoil_h'
 
 ! Local variables:
-real*4,dimension(iip1,jjp1,nsoilmx) :: data3 ! to store 3D data
-! Note iip1,jjp1 known from paramet.h; nsoilmx known from comsoil_h
-real*4,dimension(iip1,jjp1) :: data2 ! to store 2D data
+real*4,dimension(nbp_lon+1,nbp_lat,nsoilmx) :: data3 ! to store 3D data
+real*4,dimension(nbp_lon+1,nbp_lat) :: data2 ! to store 2D data
 real*4 :: data0 ! to store 0D data
 integer :: i,j,l ! for loops
@@ -44,4 +40,11 @@
 
 real*4,save :: date ! time counter (in elapsed days)
+
+real :: inertia((nbp_lon+1),nbp_lat,nsoilmx)
+real :: area((nbp_lon+1),nbp_lat)
+
+real :: inertiafi_glo(klon_glo,nsoilmx)
+real :: areafi_glo(klon_glo)
+
 integer,save :: isample ! sample rate at which data is to be written to output
 integer,save :: ntime=0 ! counter to internally store time steps
@@ -61,7 +64,7 @@
 ! Added to work in parallel mode
 real dx3_glop(klon_glo,nsoilmx)
-real dx3_glo(iim,jjp1,nsoilmx) ! to store a global 3D data set
+real dx3_glo(nbp_lon,nbp_lat,nsoilmx) ! to store a global 3D data set
 real dx2_glop(klon_glo)
-real dx2_glo(iim,jjp1)     ! to store a global 2D (surface) data set
+real dx2_glo(nbp_lon,nbp_lat)     ! to store a global 2D (surface) data set
 real px2(ngrid)
 #endif
@@ -92,8 +95,38 @@
     stop
    endif
+
+#ifdef CPP_PARA
+   ! Gather inertiedat() soil thermal inertia on physics grid
+   call Gather(inertiedat,inertiafi_glo)
+   ! Gather airephy() mesh area on physics grid
+   call Gather(airephy,areafi_glo)
+#else
+         inertiafi_glo(:,:)=inertiedat(:,:)
+         areafi_glo(:)=airephy(:)
+#endif
+
+   ! build inertia() and area()
+   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)
+     ! for area, divide at the poles by nbp_lon
+     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
+     ig0= 1+(j-2)*nbp_lon
+     do i=1,nbp_lon
+        inertia(i,j,1:nsoilmx)=inertiafi_glo(ig0+i,1:nsoilmx)
+        area(i,j)=areafi_glo(ig0+i)
+     enddo
+     ! handle redundant point in longitude
+     inertia(nbp_lon+1,j,1:nsoilmx)=inertia(1,j,1:nsoilmx)
+     area(nbp_lon+1,j)=area(1,j)
+   enddo
+   
+   ! write "header" of file (longitudes, latitudes, geopotential, ...)
+   call iniwritesoil(nid,ngrid,inertia,area)
+
   endif ! of if (is_master)
-
-  ! Define dimensions and axis attributes
-  call iniwritesoil(nid,ngrid)
   
   ! set zitau to -1 to be compatible with zitau incrementation step below
@@ -149,6 +182,6 @@
     call Grid1Dto2D_glo(dx3_glop,dx3_glo)
     ! copy dx3_glo() to dx3(:) and add redundant longitude
-    data3(1:iim,:,:)=dx3_glo(1:iim,:,:)
-    data3(iip1,:,:)=data3(1,:,:)
+    data3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
+    data3(nbp_lon+1,:,:)=data3(1,:,:)
   endif
 !$OMP END MASTER
@@ -157,15 +190,15 @@
   do l=1,nsoilmx
     ! handle the poles
-    do i=1,iip1
+    do i=1,nbp_lon+1
       data3(i,1,l)=px(1,l)
-      data3(i,jjp1,l)=px(ngrid,l)
+      data3(i,nbp_lat,l)=px(ngrid,l)
     enddo
     ! rest of the grid
-    do j=2,jjm
-      ig0=1+(j-2)*iim
-      do i=1,iim
+    do j=2,nbp_lat-1
+      ig0=1+(j-2)*nbp_lon
+      do i=1,nbp_lon
         data3(i,j,l)=px(ig0+i,l)
       enddo
-      data3(iip1,j,l)=data3(1,j,l) ! extra (modulo) longitude
+      data3(nbp_lon+1,j,l)=data3(1,j,l) ! extra (modulo) longitude
     enddo
   enddo
@@ -196,6 +229,6 @@
   corners(4)=ntime
   
-  edges(1)=iip1
-  edges(2)=jjp1
+  edges(1)=nbp_lon+1
+  edges(2)=nbp_lat
   edges(3)=nsoilmx
   edges(4)=1
@@ -224,6 +257,6 @@
     call Grid1Dto2D_glo(dx2_glop,dx2_glo)
     ! copy dx3_glo() to dx3(:) and add redundant longitude
-    data2(1:iim,:)=dx2_glo(1:iim,:)
-    data2(iip1,:)=data2(1,:)
+    data2(1:nbp_lon,:)=dx2_glo(1:nbp_lon,:)
+    data2(nbp_lon+1,:)=data2(1,:)
   endif
 !$OMP END MASTER
@@ -231,15 +264,15 @@
 #else
   ! handle the poles
-  do i=1,iip1
+  do i=1,nbp_lon+1
     data2(i,1)=px(1,1)
-    data2(i,jjp1)=px(ngrid,1)
+    data2(i,nbp_lat)=px(ngrid,1)
   enddo
   ! rest of the grid
-  do j=2,jjm
-    ig0=1+(j-2)*iim
-    do i=1,iim
+  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(iip1,j)=data2(1,j) ! extra (modulo) longitude
+    data2(nbp_lon+1,j)=data2(1,j) ! extra (modulo) longitude
   enddo
 #endif
@@ -267,6 +300,6 @@
   corners(3)=ntime
   
-  edges(1)=iip1
-  edges(2)=jjp1
+  edges(1)=nbp_lon+1
+  edges(2)=nbp_lat
   edges(3)=1
   
