Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 1528)
+++ /trunk/LMDZ.GENERIC/README	(revision 1529)
@@ -1152,2 +1152,13 @@
 - Got rid of references to "control_mod" from the physics. Added a couple
   of relevent variables for outputs in time_phylmdz_mod.
+
+== 05/04/2016 == EM
+- Got rid of references to "dimensions.h" from physics packages:
+  use nbp_lon (=iim), nbp_lat (=jjp1) and nbp_lev 
+  from module mod_grid_phy_lmdz (in phy_common) instead.
+- Removed module "comhdiff_mod.F90", as it is only used by module
+  surf_heat_transp_mod.F90, moved module variables there.
+- Addedin "surf_heat_transp_mod" local versions of some arrays
+  and routines (from dyn3d) required to compute gradient, divergence, etc. 
+  on the global dynamics grid.
+  As before, the slab ocean only works in serial.
Index: /trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/iniphysiq_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/iniphysiq_mod.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/iniphysiq_mod.F90	(revision 1529)
@@ -4,5 +4,5 @@
 
 subroutine iniphysiq(ii,jj,nlayer,punjours, pdayref,ptimestep,           &
-                     rlatu,rlatv,rlonu,rlonv,aire,cu,cv,                 &
+               rlatudyn,rlatvdyn,rlonudyn,rlonvdyn,airedyn,cudyn,cvdyn,  &
                      prad,pg,pr,pcpp,iflag_phys)
 
@@ -21,4 +21,5 @@
                        rlond, & ! longitudes
                        rlatd ! latitudes
+use surf_heat_transp_mod, only: ini_surf_heat_transp
 use infotrac, only : nqtot ! number of advected tracers
 use planete_mod, only: ini_planete_mod
@@ -29,7 +30,10 @@
                               north_east, north_west, &
                               south_west, south_east
+use ioipsl_getin_p_mod, only: getin_p
 
 implicit none
 include "dimensions.h"
+include "paramet.h"
+include "comgeom.h"
 include "iniprint.h"
 
@@ -43,11 +47,11 @@
 integer,intent(in) :: ii ! number of atmospheric coulumns along longitudes
 integer,intent(in) :: jj  ! number of atompsheric columns along latitudes
-real,intent(in) :: rlatu(jj+1) ! latitudes of the physics grid
-real,intent(in) :: rlatv(jj) ! latitude boundaries of the physics grid
-real,intent(in) :: rlonv(ii+1) ! longitudes of the physics grid
-real,intent(in) :: rlonu(ii+1) ! longitude boundaries of the physics grid
-real,intent(in) :: aire(ii+1,jj+1) ! area of the dynamics grid (m2)
-real,intent(in) :: cu((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
-real,intent(in) :: cv((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
+real,intent(in) :: rlatudyn(jj+1) ! latitudes of the physics grid
+real,intent(in) :: rlatvdyn(jj) ! latitude boundaries of the physics grid
+real,intent(in) :: rlonvdyn(ii+1) ! longitudes of the physics grid
+real,intent(in) :: rlonudyn(ii+1) ! longitude boundaries of the physics grid
+real,intent(in) :: airedyn(ii+1,jj+1) ! area of the dynamics grid (m2)
+real,intent(in) :: cudyn((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
+real,intent(in) :: cvdyn((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
 integer,intent(in) :: pdayref ! reference day of for the simulation
 real,intent(in) :: ptimestep !physics time step (s)
@@ -60,4 +64,5 @@
 real :: total_area_phy, total_area_dyn
 real :: pi
+logical :: ok_slab_ocean
 
 ! boundaries, on global grid
@@ -99,49 +104,49 @@
   
 DO i=1,ii
-   boundslon_reg(i,east)=rlonu(i) 
-   boundslon_reg(i,west)=rlonu(i+1) 
+   boundslon_reg(i,east)=rlonudyn(i) 
+   boundslon_reg(i,west)=rlonudyn(i+1) 
 ENDDO
 
 boundslat_reg(1,north)= PI/2 
-boundslat_reg(1,south)= rlatv(1)
+boundslat_reg(1,south)= rlatvdyn(1)
 DO j=2,jj
-   boundslat_reg(j,north)=rlatv(j-1) 
-   boundslat_reg(j,south)=rlatv(j) 
+   boundslat_reg(j,north)=rlatvdyn(j-1) 
+   boundslat_reg(j,south)=rlatvdyn(j) 
 ENDDO
-boundslat_reg(jj+1,north)= rlatv(jj) 
+boundslat_reg(jj+1,north)= rlatvdyn(jj) 
 boundslat_reg(jj+1,south)= -PI/2
 
 ! Write values in module regular_lonlat_mod
-CALL init_regular_lonlat(ii,jj+1, rlonv(1:ii), rlatu, &
+CALL init_regular_lonlat(ii,jj+1, rlonvdyn(1:ii), rlatudyn, &
                          boundslon_reg, boundslat_reg)
 
 ! Generate global arrays on full physics grid
 allocate(latfi(klon_glo),lonfi(klon_glo),cufi(klon_glo),cvfi(klon_glo))
-latfi(1)=rlatu(1)
+latfi(1)=rlatudyn(1)
 lonfi(1)=0.
-cufi(1) = cu(1)
-cvfi(1) = cv(1)
+cufi(1) = cudyn(1)
+cvfi(1) = cvdyn(1)
 DO j=2,jj
   DO i=1,ii
-    latfi((j-2)*ii+1+i)= rlatu(j)
-    lonfi((j-2)*ii+1+i)= rlonv(i)
-    cufi((j-2)*ii+1+i) = cu((j-1)*(ii+1)+i)
-    cvfi((j-2)*ii+1+i) = cv((j-1)*(ii+1)+i)
+    latfi((j-2)*ii+1+i)= rlatudyn(j)
+    lonfi((j-2)*ii+1+i)= rlonvdyn(i)
+    cufi((j-2)*ii+1+i) = cudyn((j-1)*(ii+1)+i)
+    cvfi((j-2)*ii+1+i) = cvdyn((j-1)*(ii+1)+i)
   ENDDO
 ENDDO
-latfi(klon_glo)= rlatu(jj+1)
+latfi(klon_glo)= rlatudyn(jj+1)
 lonfi(klon_glo)= 0.
-cufi(klon_glo) = cu((ii+1)*jj+1)
-cvfi(klon_glo) = cv((ii+1)*jj-ii)
+cufi(klon_glo) = cudyn((ii+1)*jj+1)
+cvfi(klon_glo) = cvdyn((ii+1)*jj-ii)
 
 ! build airefi(), mesh area on physics grid
 allocate(airefi(klon_glo))
-CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,aire,airefi)
+CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,airedyn,airefi)
 ! Poles are single points on physics grid
-airefi(1)=sum(aire(1:ii,1))
-airefi(klon_glo)=sum(aire(1:ii,jj+1))
+airefi(1)=sum(airedyn(1:ii,1))
+airefi(klon_glo)=sum(airedyn(1:ii,jj+1))
 
 ! Sanity check: do total planet area match between physics and dynamics?
-total_area_dyn=sum(aire(1:ii,1:jj+1))
+total_area_dyn=sum(airedyn(1:ii,1:jj+1))
 total_area_phy=sum(airefi(1:klon_glo))
 IF (total_area_dyn/=total_area_phy) THEN
@@ -174,4 +179,13 @@
 call ini_planete_mod(nlayer,preff,ap,bp)
 
+! for slab ocean, copy over some arrays
+ok_slab_ocean=.false. ! default value
+call getin_p("ok_slab_ocean",ok_slab_ocean)
+if (ok_slab_ocean) then
+  call ini_surf_heat_transp(ip1jm,ip1jmp1,unsairez,fext,unsaire, &
+                            cu,cuvsurcv,cv,cvusurcu,aire,apoln,apols, &
+                            aireu,airev)
+endif
+
 ! copy some fundamental parameters to physics 
 ! and do some initializations 
Index: /trunk/LMDZ.GENERIC/libf/phystd/aeroptproperties.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/aeroptproperties.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/aeroptproperties.F90	(revision 1529)
@@ -31,7 +31,4 @@
 !     Varying nueff section removed by R. Wordsworth for simplicity
 !     ==============================================================
-
-!#include "dimensions.h"
-!#include "dimphys.h"
 
 !     Local variables 
Index: /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 1529)
@@ -19,6 +19,6 @@
       USE tracer_h
       use comcstfi_mod, only: pi, mugaz, cpp
-      use callkeys_mod, only: varactive,diurnal,tracer,water,nosurf,varfixed,satval,	&
-		              kastprof,strictboundcorrk,specOLR,CLFvarying
+      use callkeys_mod, only: varactive,diurnal,tracer,water,nosurf,varfixed,satval,        &
+                              kastprof,strictboundcorrk,specOLR,CLFvarying
 
       implicit none
@@ -144,14 +144,11 @@
 
       ! Local aerosol optical properties for each column on RADIATIVE grid.
-      real*8,save ::  QXVAER(L_LEVELS+1,L_NSPECTV,naerkind)
-      real*8,save ::  QSVAER(L_LEVELS+1,L_NSPECTV,naerkind)
-      real*8,save ::  GVAER(L_LEVELS+1,L_NSPECTV,naerkind)
-      real*8,save ::  QXIAER(L_LEVELS+1,L_NSPECTI,naerkind)
-      real*8,save ::  QSIAER(L_LEVELS+1,L_NSPECTI,naerkind)
-      real*8,save ::  GIAER(L_LEVELS+1,L_NSPECTI,naerkind)
-
-      !REAL :: QREFvis3d(ngrid,nlayer,naerkind)
-      !REAL :: QREFir3d(ngrid,nlayer,naerkind)
-      !save QREFvis3d, QREFir3d 
+      real*8,save,allocatable ::  QXVAER(:,:,:)
+      real*8,save,allocatable ::  QSVAER(:,:,:)
+      real*8,save,allocatable ::  GVAER(:,:,:)
+      real*8,save,allocatable ::  QXIAER(:,:,:)
+      real*8,save,allocatable ::  QSIAER(:,:,:)
+      real*8,save,allocatable ::  GIAER(:,:,:)
+
       real, dimension(:,:,:), save, allocatable :: QREFvis3d
       real, dimension(:,:,:), save, allocatable :: QREFir3d
@@ -188,13 +185,13 @@
 
 
-      qxvaer(:,:,:)=0.0
-      qsvaer(:,:,:)=0.0
-      gvaer(:,:,:) =0.0
-
-      qxiaer(:,:,:)=0.0
-      qsiaer(:,:,:)=0.0
-      giaer(:,:,:) =0.0
-
       if(firstcall) then
+
+        ! test on allocated necessary because of CLFvarying (two calls to callcorrk in physiq)
+        if(.not.allocated(QXVAER)) allocate(QXVAER(L_LEVELS+1,L_NSPECTV,naerkind))
+        if(.not.allocated(QSVAER)) allocate(QSVAER(L_LEVELS+1,L_NSPECTV,naerkind))
+        if(.not.allocated(GVAER)) allocate(GVAER(L_LEVELS+1,L_NSPECTV,naerkind))
+        if(.not.allocated(QXIAER)) allocate(QXIAER(L_LEVELS+1,L_NSPECTI,naerkind))
+        if(.not.allocated(QSIAER)) allocate(QSIAER(L_LEVELS+1,L_NSPECTI,naerkind))
+        if(.not.allocated(GIAER)) allocate(GIAER(L_LEVELS+1,L_NSPECTI,naerkind))
 
          !!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...)
@@ -212,5 +209,5 @@
          endif
          call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
-	 
+         
          
 !--------------------------------------------------
@@ -264,10 +261,16 @@
       end if ! of if (firstcall)
 
-
 !=======================================================================
 !          I.b  Initialization on every call   
 !=======================================================================
  
- 
+      qxvaer(:,:,:)=0.0
+      qsvaer(:,:,:)=0.0
+      gvaer(:,:,:) =0.0
+
+      qxiaer(:,:,:)=0.0
+      qsiaer(:,:,:)=0.0
+      giaer(:,:,:) =0.0
+
 !--------------------------------------------------
 !     Effective radius and variance of the aerosols
@@ -277,11 +280,11 @@
 
          if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! Treat condensed co2 particles.
-	    call co2_reffrad(ngrid,nlayer,nq,pq,reffrad(1,1,iaero_co2))
+            call co2_reffrad(ngrid,nlayer,nq,pq,reffrad(1,1,iaero_co2))
             print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
             print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
-	 end if
+         end if
          
          if ((iaer.eq.iaero_h2o).and.water) then ! Treat condensed water particles. To be generalized for other aerosols ...
-	    call h2o_reffrad(ngrid,nlayer,pq(1,1,igcm_h2o_ice),pt, &
+            call h2o_reffrad(ngrid,nlayer,pq(1,1,igcm_h2o_ice),pt, &
                              reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
             print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
@@ -290,15 +293,15 @@
          
          if(iaer.eq.iaero_dust)then
-	    call dust_reffrad(ngrid,nlayer,reffrad(1,1,iaero_dust))
+            call dust_reffrad(ngrid,nlayer,reffrad(1,1,iaero_dust))
             print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
          endif
          
          if(iaer.eq.iaero_h2so4)then
-	    call h2so4_reffrad(ngrid,nlayer,reffrad(1,1,iaero_h2so4))
+            call h2so4_reffrad(ngrid,nlayer,reffrad(1,1,iaero_h2so4))
             print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
          endif
          
           if(iaer.eq.iaero_back2lay)then
-	    call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
+            call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
          endif
          
@@ -321,5 +324,5 @@
            reffrad,QREFvis3d,QREFir3d,                             & 
            tau_col,cloudfrac,totcloudfrac,clearsky)                
-	  
+          
 
 
@@ -475,14 +478,14 @@
          DO nw=1,L_NSPECTV ! Short Wave loop.
             albv(nw)=albedo(ig,nw)
-	 ENDDO
+         ENDDO
 
          if (nosurf) then ! Case with no surface.
             DO nw=1,L_NSPECTV
                if(albv(nw).gt.0.0) then
-	          print*,'For open lower boundary in callcorrk must'
+                  print*,'For open lower boundary in callcorrk must'
                   print*,'have spectral surface band albedos all set to zero!'
                   call abort
-	       endif
-	    ENDDO         
+               endif
+            ENDDO         
          endif
 
@@ -770,5 +773,5 @@
 
          Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
-	 glat_ig=glat(ig)
+         glat_ig=glat(ig)
 
 !-----------------------------------------------------------------------
@@ -798,5 +801,5 @@
          else ! During the night, fluxes = 0.
             nfluxtopv       = 0.0d0
-	    fluxtopvdn      = 0.0d0
+            fluxtopvdn      = 0.0d0
             nfluxoutv_nu(:) = 0.0d0
             nfluxgndv_nu(:) = 0.0d0
@@ -813,15 +816,15 @@
             surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))      
             if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive.
-	       DO nw=1,L_NSPECTV	          
-	          albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
+               DO nw=1,L_NSPECTV                  
+                  albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
                ENDDO
-	       albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
+               albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
                albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
             else
                albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
             endif
-	 else
-	    albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
-	 endif
+         else
+            albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
+         endif
 
 
@@ -956,5 +959,5 @@
         IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
 !$OMP END MASTER
-!$OMP BARRIER	
+!$OMP BARRIER
         IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
         IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
Index: unk/LMDZ.GENERIC/libf/phystd/comhdiff_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/comhdiff_mod.F90	(revision 1528)
+++ 	(revision )
@@ -1,10 +1,0 @@
-MODULE comhdiff_mod
-IMPLICIT NONE  
-#include "dimensions.h"
-#include "paramet.h"
-
-  LOGICAL alpha_var,slab_upstream
-  REAL zmasqu(ip1jmp1),zmasqv(ip1jm),	&
-	unsfv(ip1jm),unsev(ip1jm),unsfu(ip1jmp1),unseu(ip1jmp1)
-
-END MODULE comhdiff_mod
Index: unk/LMDZ.GENERIC/libf/phystd/dyn1d/diverg.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/diverg.F	(revision 1528)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dyn3d/diverg.F
Index: unk/LMDZ.GENERIC/libf/phystd/dyn1d/gr_dyn_fi.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/gr_dyn_fi.F	(revision 1528)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dynlonlat_phylonlat/gr_dyn_fi.F
Index: unk/LMDZ.GENERIC/libf/phystd/dyn1d/gr_fi_dyn.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/gr_fi_dyn.F	(revision 1528)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dynlonlat_phylonlat/gr_fi_dyn.F
Index: unk/LMDZ.GENERIC/libf/phystd/dyn1d/gr_u_scal.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/gr_u_scal.F	(revision 1528)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dyn3d/gr_u_scal.F
Index: unk/LMDZ.GENERIC/libf/phystd/dyn1d/gr_v_scal.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/gr_v_scal.F	(revision 1528)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dyn3d/gr_v_scal.F
Index: unk/LMDZ.GENERIC/libf/phystd/dyn1d/grad.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/grad.F	(revision 1528)
+++ 	(revision )
@@ -1,1 +1,0 @@
-link ../../dyn3d/grad.F
Index: /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/dyn1d/rcm1d.F	(revision 1529)
@@ -907,22 +907,3 @@
       end                       !rcm1d
  
-c***********************************************************************
-c***********************************************************************
-c     Subroutines Bidons utilise seulement en 3D, mais
-c     necessaire a la compilation de rcm1d en 1D
-
-!      subroutine gr_fi_dyn
-!      RETURN
-!      END
- 
-c***********************************************************************
-c***********************************************************************
-
-!#include "../dyn3d/disvert.F"
-!#include "../dyn3d/abort_gcm.F"
-!#include "../dyn3d/diverg.F"
-!#include "../dyn3d/grad.F"
-!#include "../dyn3d/gr_u_scal.F"
-!#include "../dyn3d/gr_v_scal.F"
-!#include "../dyn3d/gr_dyn_fi.F"
-
+
Index: /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 1529)
@@ -9,5 +9,6 @@
              prad,pg,pr,pcpp)
 
-  use radinc_h, only : naerkind
+  use radinc_h, only: ini_radinc_h, naerkind
+  use radcommon_h, only: ini_radcommon_h
   use datafile_mod, only: datadir
   use comdiurn_h, only: sinlat, coslat, sinlon, coslon
@@ -21,4 +22,5 @@
   use planetwide_mod, only: planetwide_sumval
   use callkeys_mod
+  use mod_phys_lmdz_para, only : is_parallel
 
 !=======================================================================
@@ -74,11 +76,4 @@
   REAL SSUM
  
-  CHARACTER ch1*12
-  CHARACTER ch80*80
-
-  logical chem, h2o
-
-  real psurf,pN2 ! added by RW for Gliese 581d N2+CO2
-
   ! initialize constants in comcstfi_mod
   rad=prad
@@ -333,4 +328,9 @@
      call getin_p("ok_slab_ocean",ok_slab_ocean)
      write(*,*) "ok_slab_ocean = ",ok_slab_ocean
+     ! Sanity check: for now slab oncean only works in serial mode
+     if (ok_slab_ocean.and.is_parallel) then
+       write(*,*) " Error: slab ocean should only be used in serial mode!"
+       call abort
+     endif
 
      write(*,*) "Use slab-sea-ice ?"
@@ -718,4 +718,10 @@
   ENDDO
 
+  ! initialize variables in radinc_h
+  call ini_radinc_h(nlayer)
+ 
+  ! allocate "radcommon_h" arrays
+  call ini_radcommon_h()
+
   ! allocate "comsoil_h" arrays
   call ini_comsoil_h(ngrid)
Index: /trunk/LMDZ.GENERIC/libf/phystd/inistats.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/inistats.F	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/inistats.F	(revision 1529)
@@ -1,19 +1,20 @@
       subroutine inistats(ierr)
 
+      use statto_mod, only: istats,istime
       use mod_phys_lmdz_para, only : is_master
-      use statto_mod, only: istats,istime
       USE comvert_mod, ONLY: ap,bp,aps,bps,preff,pseudoalt,presnivs
-      USE comconst_mod, ONLY: daysec,dtphys,pi
+      USE comconst_mod, ONLY: pi
+      USE time_phylmdz_mod, ONLY: daysec,dtphys
+      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
       implicit none
 
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "netcdf.inc"
+      include "netcdf.inc"
 
       integer,intent(out) :: ierr
       integer :: nid
       integer :: l,nsteppd
-      real, dimension(llm) ::  sig_s
+      real, dimension(nbp_lev) ::  sig_s
+      real :: lon_reg_ext(nbp_lon+1) ! extended longitudes
       integer :: idim_lat,idim_lon,idim_llm,idim_llmp1,idim_time
       real, dimension(istime) :: lt
@@ -39,8 +40,12 @@
       write (*,*) 
 
-      do l= 1, llm
+      do l= 1, nbp_lev
          sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2.
          pseudoalt(l)=-10.*log(presnivs(l)/preff)   
       enddo
+      
+      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 (is_master) then
@@ -53,8 +58,8 @@
       endif
 
-      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_lat)
-      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_lon)
-      ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
-      ierr = NF_DEF_DIM (nid, "llmp1", llm+1, idim_llmp1)
+      ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_lat)
+      ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_lon)
+      ierr = NF_DEF_DIM (nid, "altitude", nbp_lev, idim_llm)
+      ierr = NF_DEF_DIM (nid, "llmp1", nbp_lev+1, idim_llmp1)
       ierr = NF_DEF_DIM (nid, "Time", NF_UNLIMITED, idim_time)
 
@@ -68,14 +73,14 @@
      &            "degrees_north",1,idim_lat,nvarid,ierr)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
 #else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
 #endif
       call def_var_stats(nid,"longitude","East longitude",
      &            "degrees_east",1,idim_lon,nvarid,ierr)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
 #else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
 #endif
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/initracer.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/initracer.F	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/initracer.F	(revision 1529)
@@ -22,7 +22,4 @@
 c            Ehouarn Millour (oct. 2008) identify tracers by their names
 c=======================================================================
-
-
-#include "dimensions.h"
 
       integer :: ngrid,nq
Index: /trunk/LMDZ.GENERIC/libf/phystd/iniwrite.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/iniwrite.F	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/iniwrite.F	(revision 1529)
@@ -1,11 +1,13 @@
-      SUBROUTINE iniwrite(nid,idayref,phis)
+      SUBROUTINE iniwrite(nid,idayref,phis,area)
 
       use comsoil_h, only: mlayer, nsoilmx
-      use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, pi
-      use time_phylmdz_mod, only: daysec, dtphys
+      USE comcstfi_mod, only: g, mugaz, omeg, rad, rcp, pi 
       USE comvert_mod, ONLY: ap,bp,aps,bps,pseudoalt
       USE logic_mod, ONLY: fxyhypb,ysinus
       USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy
+      USE time_phylmdz_mod, ONLY: daysec, dtphys
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
       IMPLICIT NONE
 
@@ -26,8 +28,5 @@
 c   -------------
 
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "netcdf.inc"
+      include "netcdf.inc"
 
 c   Arguments:
@@ -36,5 +35,6 @@
       integer,intent(in) :: nid        ! NetCDF file ID
       INTEGER*4,intent(in) :: idayref  ! date (initial date for this run)
-      real,intent(in) :: phis(ip1jmp1) ! surface geopotential
+      real,intent(in) :: phis(nbp_lon+1,nbp_lat) ! surface geopotential
+      real,intent(in) :: area(nbp_lon+1,nbp_lat) ! mesh area (m2)
 
 c   Local:
@@ -44,7 +44,8 @@
       REAL tab_cntrl(length) ! run parameters are stored in this array
       INTEGER ierr
-
-      integer :: nvarid,idim_index,idim_rlonu,idim_rlonv
-      integer :: idim_rlatu,idim_rlatv,idim_llmp1,idim_llm
+      REAl :: lon_reg_ext(nbp_lon+1) ! extended longitudes
+
+      integer :: nvarid,idim_index,idim_rlonv
+      integer :: idim_rlatu,idim_llmp1,idim_llm
       integer :: idim_nsoilmx ! "subsurface_layers" dimension ID #
       integer, dimension(2) :: id  
@@ -54,7 +55,7 @@
          tab_cntrl(l)=0.
       ENDDO
-      tab_cntrl(1)  = real(iim)
-      tab_cntrl(2)  = real(jjm)
-      tab_cntrl(3)  = real(llm)
+      tab_cntrl(1)  = real(nbp_lon)
+      tab_cntrl(2)  = real(nbp_lat-1)
+      tab_cntrl(3)  = real(nbp_lev)
       tab_cntrl(4)  = real(idayref)
       tab_cntrl(5)  = rad
@@ -99,10 +100,10 @@
 
       ierr = NF_DEF_DIM (nid, "index", length, idim_index)
-      ierr = NF_DEF_DIM (nid, "rlonu", iip1, idim_rlonu)
-      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_rlatu)
-      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_rlonv)
-      ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv)
-      ierr = NF_DEF_DIM (nid, "interlayer", (llm+1), idim_llmp1)
-      ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
+!      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)
+!      ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv)
+      ierr = NF_DEF_DIM (nid, "interlayer", (nbp_lev+1), idim_llmp1)
+      ierr = NF_DEF_DIM (nid, "altitude", nbp_lev, idim_llm)
       ierr = NF_DEF_DIM (nid,"subsurface_layers",nsoilmx,idim_nsoilmx)
 c
@@ -129,18 +130,19 @@
 c --------------------------
 c  longitudes and latitudes
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "rlonu", NF_DOUBLE, 1, idim_rlonu,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "rlonu", NF_FLOAT, 1, idim_rlonu,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 21,
-     .                       "Longitudes at u nodes")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonu/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonu/pi*180)
-#endif
+!
+!      ierr = NF_REDEF (nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_DEF_VAR (nid, "rlonu", NF_DOUBLE, 1, idim_rlonu,nvarid)
+!#else
+!      ierr = NF_DEF_VAR (nid, "rlonu", NF_FLOAT, 1, idim_rlonu,nvarid)
+!#endif
+!      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 21,
+!     .                       "Longitudes at u nodes")
+!      ierr = NF_ENDDEF(nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonu/pi*180)
+!#else
+!      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonu/pi*180)
+!#endif
 c
 c --------------------------
@@ -156,10 +158,14 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
-#endif
-c
-c --------------------------
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
+#endif
+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)
+
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
@@ -173,7 +179,7 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
 #endif
 c
@@ -199,18 +205,18 @@
 c
 c --------------------------
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "rlatv", NF_DOUBLE, 1, idim_rlatv,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "rlatv", NF_FLOAT, 1, idim_rlatv,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
-     .                       "Latitudes at v nodes")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatv/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatv/pi*180)
-#endif
+!      ierr = NF_REDEF (nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_DEF_VAR (nid, "rlatv", NF_DOUBLE, 1, idim_rlatv,nvarid)
+!#else
+!      ierr = NF_DEF_VAR (nid, "rlatv", NF_FLOAT, 1, idim_rlatv,nvarid)
+!#endif
+!      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
+!     .                       "Latitudes at v nodes")
+!      ierr = NF_ENDDEF(nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatv/pi*180)
+!#else
+!      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatv/pi*180)
+!#endif
 c
 c --------------------------
@@ -274,40 +280,40 @@
 c  Mesh area and conversion coefficients cov. <-> contra. <--> natural
 
-      id(1)=idim_rlonu
-      id(2)=idim_rlatu
-c
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "cu", NF_DOUBLE, 2, id,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "cu", NF_FLOAT, 2, id,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
-     .             "Conversion coefficients cov <--> natural")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cu)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,cu)
-#endif
-c
-      id(1)=idim_rlonv
-      id(2)=idim_rlatv
-c
-c --------------------------
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "cv", NF_DOUBLE, 2, id,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "cv", NF_FLOAT, 2, id,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
-     .             "Conversion coefficients cov <--> natural")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cv)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,cv)
-#endif
+!      id(1)=idim_rlonu
+!      id(2)=idim_rlatu
+c
+!      ierr = NF_REDEF (nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_DEF_VAR (nid, "cu", NF_DOUBLE, 2, id,nvarid)
+!#else
+!      ierr = NF_DEF_VAR (nid, "cu", NF_FLOAT, 2, id,nvarid)
+!#endif
+!      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
+!     .             "Conversion coefficients cov <--> natural")
+!      ierr = NF_ENDDEF(nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cu)
+!#else
+!      ierr = NF_PUT_VAR_REAL (nid,nvarid,cu)
+!#endif
+c
+!      id(1)=idim_rlonv
+!      id(2)=idim_rlatv
+c
+c --------------------------
+!      ierr = NF_REDEF (nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_DEF_VAR (nid, "cv", NF_DOUBLE, 2, id,nvarid)
+!#else
+!      ierr = NF_DEF_VAR (nid, "cv", NF_FLOAT, 2, id,nvarid)
+!#endif
+!      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
+!     .             "Conversion coefficients cov <--> natural")
+!      ierr = NF_ENDDEF(nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cv)
+!#else
+!      ierr = NF_PUT_VAR_REAL (nid,nvarid,cv)
+!#endif
 c
       id(1)=idim_rlonv
@@ -325,7 +331,7 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,aire)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,area)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,area)
 #endif
 c
@@ -350,7 +356,8 @@
 c
 
-      write(*,*)'iniwrite: iim,jjm,llm,idayref',iim,jjm,llm,idayref
+      write(*,*)'iniwrite: nbp_lon,nbp_lat,nbp_lev,idayref',
+     & nbp_lon,nbp_lat,nbp_lev,idayref
       write(*,*)'iniwrite: rad,omeg,g,mugaz,rcp',
-     s rad,omeg,g,mugaz,rcp
+     & rad,omeg,g,mugaz,rcp
       write(*,*)'iniwrite: daysec,dtphys',daysec,dtphys
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specIR.F	(revision 1529)
@@ -1,3 +1,3 @@
-      SUBROUTINE iniwrite_specIR(nid,idayref)
+      SUBROUTINE iniwrite_specIR(nid,idayref,area)
 
       use radinc_h, only: L_NSPECTI
@@ -8,4 +8,6 @@
       USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
 
       implicit none
@@ -27,9 +29,5 @@
 c   -------------
 
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "netcdf.inc"
-!#include"dimphys.h"
+      include "netcdf.inc"
 
 c   Arguments:
@@ -38,4 +36,5 @@
       integer,intent(in) :: nid        ! NetCDF file ID
       INTEGER*4,intent(in) :: idayref  ! date (initial date for this run)
+      real,intent(in) :: area(nbp_lon+1,nbp_lat) ! mesh area (m2)
 
 c   Local:
@@ -45,4 +44,5 @@
       REAL tab_cntrl(length) ! run parameters are stored in this array
       INTEGER ierr
+      REAl :: lon_reg_ext(nbp_lon+1) ! extended longitudes
 
       integer :: nvarid,idim_index,idim_rlonu,idim_rlonv
@@ -57,7 +57,7 @@
          tab_cntrl(l)=0.
       ENDDO
-      tab_cntrl(1)  = FLOAT(iim)
-      tab_cntrl(2)  = FLOAT(jjm)
-      tab_cntrl(3)  = FLOAT(llm)
+      tab_cntrl(1)  = FLOAT(nbp_lon)
+      tab_cntrl(2)  = FLOAT(nbp_lat-1)
+      tab_cntrl(3)  = FLOAT(nbp_lev)
       tab_cntrl(4)  = FLOAT(idayref)
       tab_cntrl(5)  = rad
@@ -102,6 +102,6 @@
 
       ierr = NF_DEF_DIM (nid, "index", length, idim_index)
-      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_rlatu)
-      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_rlonv)
+      ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_rlatu)
+      ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_rlonv)
       ierr = NF_DEF_DIM (nid, "IR Wavenumber",L_NSPECTI,idim_bandsIR)
 
@@ -141,10 +141,14 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
-#endif
-c
-c --------------------------
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
+#endif
+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)
+      
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
@@ -158,7 +162,7 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
 #endif
 c
@@ -172,8 +176,8 @@
 #ifdef NC_DOUBLE
       ierr=NF_DEF_VAR(nid,"IR Wavenumber",NF_DOUBLE,1,
-     .				idim_bandsIR,nvarid)
+     .                                idim_bandsIR,nvarid)
 #else
       ierr=NF_DEF_VAR(nid,"IR Wavenumber",NF_FLOAT,1,
-     .				idim_bandsIR,nvarid)
+     .                                idim_bandsIR,nvarid)
 #endif
       ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 34,
@@ -196,8 +200,8 @@
 #ifdef NC_DOUBLE
       ierr=NF_DEF_VAR(nid,"IR Bandwidth",NF_DOUBLE,1,
-     .				idim_bandsIR,nvarid)
+     .                                idim_bandsIR,nvarid)
 #else
       ierr=NF_DEF_VAR(nid,"IR Bandwidth",NF_FLOAT,1,
-     .				idim_bandsIR,nvarid)
+     .                                idim_bandsIR,nvarid)
 #endif
       ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 25,
@@ -230,7 +234,7 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,aire)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,area)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,area)
 #endif
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/iniwrite_specVI.F	(revision 1529)
@@ -1,7 +1,6 @@
-      SUBROUTINE iniwrite_specVI(nid,idayref)
+      SUBROUTINE iniwrite_specVI(nid,idayref,area)
 
       use radinc_h, only: L_NSPECTV
       use radcommon_h, only: WNOV,DWNV
-      use comsoil_h
       use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, pi
       use time_phylmdz_mod, only: daysec, dtphys
@@ -9,4 +8,6 @@
       USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
 
       implicit none
@@ -28,15 +29,12 @@
 c   -------------
 
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "netcdf.inc"
-!#include"dimphys.h"
+      include "netcdf.inc"
 
 c   Arguments:
 c   ----------
 
-      integer nid        ! NetCDF file ID
-      INTEGER*4 idayref  ! date (initial date for this run)
+      integer,intent(in) :: nid        ! NetCDF file ID
+      INTEGER*4,INTENT(IN) :: idayref  ! date (initial date for this run)
+      real,intent(in) :: area(nbp_lon+1,nbp_lat) ! mesh area (m2)
 
 c   Local:
@@ -46,4 +44,5 @@
       REAL tab_cntrl(length) ! run parameters are stored in this array
       INTEGER ierr
+      REAl :: lon_reg_ext(nbp_lon+1) ! extended longitudes
 
       integer :: nvarid,idim_index,idim_rlonu,idim_rlonv
@@ -58,7 +57,7 @@
          tab_cntrl(l)=0.
       ENDDO
-      tab_cntrl(1)  = FLOAT(iim)
-      tab_cntrl(2)  = FLOAT(jjm)
-      tab_cntrl(3)  = FLOAT(llm)
+      tab_cntrl(1)  = FLOAT(nbp_lon)
+      tab_cntrl(2)  = FLOAT(nbp_lat-1)
+      tab_cntrl(3)  = FLOAT(nbp_lev)
       tab_cntrl(4)  = FLOAT(idayref)
       tab_cntrl(5)  = rad
@@ -103,6 +102,6 @@
 
       ierr = NF_DEF_DIM (nid, "index", length, idim_index)
-      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_rlatu)
-      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_rlonv)
+      ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_rlatu)
+      ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_rlonv)
       ierr = NF_DEF_DIM (nid, "VI Wavenumber",L_NSPECTV,idim_bandsVI)
 
@@ -141,10 +140,14 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
-#endif
-c
-c --------------------------
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
+#endif
+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)
+
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
@@ -158,7 +161,7 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
 #endif
 c
@@ -172,8 +175,8 @@
 #ifdef NC_DOUBLE
       ierr=NF_DEF_VAR(nid,"VI Wavenumber",NF_DOUBLE,1,
-     .				idim_bandsVI,nvarid)
+     .                                idim_bandsVI,nvarid)
 #else
       ierr=NF_DEF_VAR(nid,"VI Wavenumber",NF_FLOAT,1,
-     .				idim_bandsVI,nvarid)
+     .                                idim_bandsVI,nvarid)
 #endif
       ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 33,
@@ -197,8 +200,8 @@
 #ifdef NC_DOUBLE
       ierr=NF_DEF_VAR(nid,"VI Bandwidth",NF_DOUBLE,1,
-     .				idim_bandsVI,nvarid)
+     .                                idim_bandsVI,nvarid)
 #else
       ierr=NF_DEF_VAR(nid,"VI Bandwidth",NF_FLOAT,1,
-     .				idim_bandsVI,nvarid)
+     .                                idim_bandsVI,nvarid)
 #endif
       ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 24,
@@ -231,7 +234,7 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,aire)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,area)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,area)
 #endif
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/iniwritesoil.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/iniwritesoil.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/iniwritesoil.F90	(revision 1529)
@@ -1,3 +1,3 @@
-subroutine iniwritesoil(nid,ngrid)
+subroutine iniwritesoil(nid,ngrid,inertia,area)
 
 ! initialization routine for 'writediagoil'. Here we create/define
@@ -5,17 +5,18 @@
 ! (time-independent) parameters.
 
-use comsoil_h, only: mlayer, inertiedat, nsoilmx
-use comcstfi_mod, only: pi
+use comsoil_h, only: mlayer, nsoilmx
+USE comcstfi_mod, only: pi
+USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+use mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
 
 implicit none
 
-#include"dimensions.h"
-#include"paramet.h"
-#include"comgeom.h"
-#include"netcdf.inc"
+include"netcdf.inc"
 
 ! Arguments:
 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)
 
 ! Local variables:
@@ -30,6 +31,7 @@
 integer,dimension(3) :: dimids ! to store IDs of dimensions of a variable
 character(len=60) :: text ! to store some text
-real,dimension(iip1,jjp1,nsoilmx) :: data3 ! to store 3D data
+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
 
 ! 1. Define the dimensions
@@ -38,16 +40,14 @@
 
 ! Define the dimensions
-ierr=NF_DEF_DIM(nid,"longitude",iip1,idim_rlonv)
-! iip1 known from paramet.h
+ierr=NF_DEF_DIM(nid,"longitude",nbp_lon+1,idim_rlonv)
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not define longitude dimension"
 endif
-ierr=NF_DEF_DIM(nid,"latitude",jjp1,idim_rlatu)
-! jjp1 known from paramet.h
+ierr=NF_DEF_DIM(nid,"latitude",nbp_lat,idim_rlatu)
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not define latitude dimension"
 endif
 ierr=NF_DEF_DIM(nid,"depth",nsoilmx,idim_depth)
-! nsoilmx known from dimphys.h
+! nsoilmx known from comsoil_h
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not define depth dimension"
@@ -81,13 +81,17 @@
 ierr=NF_PUT_ATT_TEXT(nid,varid,"units",len_trim(text),text)
 
+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)
+
 ! Write longitude to file
 ierr=NF_ENDDEF(nid) ! switch out of NetCDF define mode
 ! Write
 #ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(nid,varid,rlonv*(180./pi))
-#else
-ierr=NF_PUT_VAR_REAL(nid,varid,rlonv*(180./pi))
-#endif
-! Note: rlonv is known from comgeom.h
+ierr=NF_PUT_VAR_DOUBLE(nid,varid,lon_reg_ext*(180./pi))
+#else
+ierr=NF_PUT_VAR_REAL(nid,varid,lon_reg_ext*(180./pi))
+#endif
+! Note: rlonv is known from comgeom.h and pi from comcstfi.h
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not write longitude variable"
@@ -117,9 +121,8 @@
 ! Write
 #ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(nid,varid,rlatu*(180./pi))
-#else
-ierr=NF_PUT_VAR_REAL(nid,varid,rlatu*(180./pi))
-#endif
-! Note: rlatu is known from comgeom.h
+ierr=NF_PUT_VAR_DOUBLE(nid,varid,lat_reg*(180./pi))
+#else
+ierr=NF_PUT_VAR_REAL(nid,varid,lat_reg*(180./pi))
+#endif
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not write longitude variable"
@@ -209,9 +212,8 @@
 ! Write
 #ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(nid,varid,aire)
-#else
-ierr=NF_PUT_VAR_REAL(nid,varid,aire)
-#endif
-! Note: aire is known from comgeom.h
+ierr=NF_PUT_VAR_DOUBLE(nid,varid,area)
+#else
+ierr=NF_PUT_VAR_REAL(nid,varid,area)
+#endif
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not write area variable"
@@ -240,30 +242,11 @@
 ierr=NF_PUT_ATT_TEXT(nid,varid,"units",len_trim(text),text)
 
-! Recast data along 'dynamics' grid
-! Note: inertiedat is known from comsoil_h
-
-do l=1,nsoilmx
-  ! handle the poles
-  do i=1,iip1
-    data3(i,1,l)=inertiedat(1,l)
-    data3(i,jjp1,l)=inertiedat(ngrid,l)
-  enddo
-  ! rest of the grid
-  do j=2,jjm
-    ig0=1+(j-2)*iim
-    do i=1,iim
-      data3(i,j,l)=inertiedat(ig0+i,l)
-    enddo
-    data3(iip1,j,l)=data3(1,j,l) ! extra (modulo) longitude
-  enddo
-enddo ! of do l=1,nsoilmx
-
 ! Write data2 to file
 ierr=NF_ENDDEF(nid) ! switch out of NetCDF define mode
 ! Write
 #ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(nid,varid,data3)
-#else
-ierr=NF_PUT_VAR_REAL(nid,varid,data3)
+ierr=NF_PUT_VAR_DOUBLE(nid,varid,inertia)
+#else
+ierr=NF_PUT_VAR_REAL(nid,varid,inertia)
 #endif
 if (ierr.ne.NF_NOERR) then
Index: /trunk/LMDZ.GENERIC/libf/phystd/mass_redistribution.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/mass_redistribution.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/mass_redistribution.F90	(revision 1529)
@@ -1,3 +1,3 @@
-      SUBROUTINE mass_redistribution(ngrid,nlayer,nq,ptimestep,   &
+SUBROUTINE mass_redistribution(ngrid,nlayer,nq,ptimestep,   &
                        rnat,pcapcal,pplay,pplev,pt,ptsrf,pq,pqs,     &
 		       pu,pv,pdt,pdtsrf,pdq,pdu,pdv,pdmassmr,  &
@@ -224,7 +224,7 @@
 !        Van Leer scheme:
          w(1:nlayer+1)=-zmflux(1:nlayer+1)*ptimestep
-         call vl1d(zzt,2.,zzmass,w,ztm) 
-         call vl1d(zzu ,2.,zzmass,w,zum) 
-         call vl1d(zzv ,2.,zzmass,w,zvm) 
+         call vl1d(nlayer,zzt,2.,zzmass,w,ztm) 
+         call vl1d(nlayer,zzu,2.,zzmass,w,zum) 
+         call vl1d(nlayer,zzv,2.,zzmass,w,zvm) 
          do iq=1,nq
            zq1(1:nlayer)=zzq(1:nlayer,iq)
@@ -232,5 +232,5 @@
 !		print*,iq
 !		print*,zq1
-           call vl1d(zq1,2.,zzmass,w,zqm1)
+           call vl1d(nlayer,zq1,2.,zzmass,w,zqm1)
            do l=2,nlayer
               zzq(l,iq)=zq1(l)
@@ -281,11 +281,8 @@
       END DO  ! loop on ig 
 
-      return
-      end
-
-
+CONTAINS
 
 ! *****************************************************************
-      SUBROUTINE vl1d(q,pente_max,zzmass,w,qm)
+      SUBROUTINE vl1d(llm,q,pente_max,zzmass,w,qm)
 !
 !    
@@ -299,10 +296,10 @@
 !
 !   --------------------------------------------------------------------
+      
       IMPLICIT NONE
-
-#include "dimensions.h"
 
 !   Arguments:
 !   ----------
+      integer,intent(in) :: llm
       real zzmass(llm),pente_max
       REAL q(llm),qm(llm+1)
@@ -406,4 +403,5 @@
 !         end if
 
-      return
-      end
+       END SUBROUTINE vl1d
+
+END SUBROUTINE mass_redistribution
Index: /trunk/LMDZ.GENERIC/libf/phystd/mkstat.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/mkstat.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/mkstat.F90	(revision 1529)
@@ -10,24 +10,22 @@
 !  Yann W. july 2003
 
+use statto_mod, only: istime,count
 use mod_phys_lmdz_para, only : is_master
-use statto_mod, only: istime,count
+use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev
 
 implicit none
 
-#include "dimensions.h"
-#include "netcdf.inc"
+include "netcdf.inc"
 
-integer,parameter :: iip1=iim+1
-integer,parameter :: jjp1=jjm+1
 integer :: ierr,nid,nbvar,i,ndims,lt,nvarid
 integer, dimension(4) :: id,varid,start,size
 integer, dimension(5) :: dimids
 character (len=50) :: name,nameout,units,title
-real, dimension(iip1,jjp1,llm) :: sum3d,square3d,mean3d,sd3d
-real, dimension(iip1,jjp1) :: sum2d,square2d,mean2d,sd2d
+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, dimension(istime) :: time
-real, dimension(jjp1) :: lat
-real, dimension(iip1) :: lon
-real, dimension(llm) :: alt
+real, dimension(nbp_lat) :: lat
+real, dimension(nbp_lon+1) :: lon
+real, dimension(nbp_lev) :: alt
 logical :: lcopy=.true.
 !integer :: latid,lonid,altid,timeid
@@ -35,9 +33,9 @@
 !integer, dimension(4) :: dimout
 
+! Incrementation of count for the last step, which is not done in wstats
+count(istime)=count(istime)+1
+
 if (is_master) then
 ! only the master needs do this
-
-! Incrementation of count for the last step, which is not done in wstats
-count(istime)=count(istime)+1
 
 ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
@@ -107,5 +105,5 @@
 !      dimout(4)=timeid
 
-      size=(/iip1,jjp1,llm,1/)
+      size=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
       do lt=1,istime
          start=(/1,1,1,lt/)
@@ -137,5 +135,5 @@
 !      dimout(3)=timeid
 
-      size=(/iip1,jjp1,1,0/)
+      size=(/nbp_lon+1,nbp_lat,1,0/)
       do lt=1,istime
          start=(/1,1,lt,0/)
Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 1529)
@@ -26,6 +26,8 @@
       use time_phylmdz_mod, only: ecritphy, iphysiq, nday
       use phyredem, only: physdem0, physdem1
-      use slab_ice_h
+      use slab_ice_h, only: capcalocean, capcalseaice,capcalsno, &
+                            noceanmx
       use ocean_slab_mod, only :ocean_slab_init, ocean_slab_ice, &
+                                ini_surf_heat_transp_mod, &
                                 ocean_slab_get_vars,ocean_slab_final
       use surf_heat_transp_mod,only :init_masquv
@@ -558,4 +560,6 @@
                                 sea_ice, pctsrf_sic)
 
+           call ini_surf_heat_transp_mod()
+           
            knindex(:) = 0
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90	(revision 1529)
@@ -1,5 +1,5 @@
-      module radcommon_h
-
-      use radinc_h
+module radcommon_h
+      use radinc_h, only: L_NSPECTI, L_NSPECTV, L_NGAUSS, NTstar, NTstop, &
+                          naerkind, nsizemax
       implicit none
 
@@ -121,13 +121,14 @@
       REAL :: omegaREFir(naerkind,nsizemax)
 
-      REAL tstellar ! Stellar brightness temperature (SW)
+      REAL,SAVE :: tstellar ! Stellar brightness temperature (SW)
 
-      real*8 planckir(L_NSPECTI,NTstop-NTstar+1)
+      real*8,save :: planckir(L_NSPECTI,NTstop-NTstar+1)
 
-      real*8 PTOP, TAUREF(L_LEVELS+1)
+      real*8,save :: PTOP
+      real*8,save,allocatable :: TAUREF(:)
 
-      real*8, parameter :: UBARI = 0.5D0
+      real*8,parameter :: UBARI = 0.5D0
 
-      real*8 gweight(L_NGAUSS)
+      real*8,save :: gweight(L_NGAUSS)
 !$OMP THREADPRIVATE(QREFvis,QREFir,omegaREFvis,omegaREFir,& 	! gweight read by master in sugas_corrk
 		!$OMP tstellar,planckir,PTOP,TAUREF)
@@ -144,16 +145,24 @@
       real*8, parameter :: grav = 6.672E-11
 
-      real*8 Cmk
-      save Cmk
-      real*8 glat_ig
-      save glat_ig
+      real*8,save :: Cmk
+      real*8,save :: glat_ig
 !$OMP THREADPRIVATE(Cmk,glat_ig)
 
       ! extinction of incoming sunlight (Saturn's rings, eclipses, etc...)
-      REAL, DIMENSION(:), ALLOCATABLE :: eclipse
+      REAL, DIMENSION(:), ALLOCATABLE ,SAVE :: eclipse
 
       !Latitude-dependent gravity
-      REAL, DIMENSION(:), ALLOCATABLE :: glat
+      REAL, DIMENSION(:), ALLOCATABLE , SAVE :: glat
 !$OMP THREADPRIVATE(glat,eclipse)
 
-      end module radcommon_h
+contains
+
+      subroutine ini_radcommon_h
+      use radinc_h, only: L_LEVELS
+      implicit none
+      
+        allocate(TAUREF(L_LEVELS+1))
+      
+      end subroutine ini_radcommon_h
+
+end module radcommon_h
Index: /trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90	(revision 1529)
@@ -8,7 +8,7 @@
 !     water cloud optical properties
 
-      use callkeys_mod, only: radfixed,Nmix_co2,			&
-		pres_bottom_tropo,pres_top_tropo,size_tropo,	&
-		pres_bottom_strato,size_strato
+      use callkeys_mod, only: radfixed,Nmix_co2,                    &
+                pres_bottom_tropo,pres_top_tropo,size_tropo,        &
+                pres_bottom_strato,size_strato
       
       real, save ::  rad_h2o
@@ -21,5 +21,5 @@
 
 
-      contains
+contains
 
 
@@ -38,10 +38,7 @@
       use ioipsl_getin_p_mod, only: getin_p
       use radinc_h, only: naerkind
-      use aerosol_mod
-!      USE tracer_h
-      Implicit none
-
-!      include "dimensions.h"
-!      include "dimphys.h"
+      use aerosol_mod, only: iaero_back2lay, iaero_co2, iaero_dust, &
+                             iaero_h2o, iaero_h2so4
+      Implicit none
 
       integer,intent(in) :: ngrid
@@ -81,6 +78,6 @@
                nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
             endif
-	    
-	    if(iaer.eq.iaero_back2lay)then ! Two-layer aerosols
+            
+            if(iaer.eq.iaero_back2lay)then ! Two-layer aerosols
                reffrad(1:ngrid,1:nlayer,iaer) = 2.e-6
                nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
@@ -101,15 +98,15 @@
          if (radfixed) then
 
-	    write(*,*)"radius of H2O water particles:"
+            write(*,*)"radius of H2O water particles:"
             rad_h2o=13. ! default value
             call getin_p("rad_h2o",rad_h2o)
             write(*,*)" rad_h2o = ",rad_h2o
 
-	    write(*,*)"radius of H2O ice particles:"
+            write(*,*)"radius of H2O ice particles:"
             rad_h2o_ice=35. ! default value
             call getin_p("rad_h2o_ice",rad_h2o_ice)
             write(*,*)" rad_h2o_ice = ",rad_h2o_ice
 
-	 else
+         else
 
             write(*,*)"Number mixing ratio of H2O water particles:"
@@ -122,5 +119,5 @@
             call getin_p("Nmix_h2o_ice",Nmix_h2o_ice)
             write(*,*)" Nmix_h2o_ice = ",Nmix_h2o_ice
-	 endif
+         endif
 
       print*,'exit su_aer_radii'
@@ -173,10 +170,10 @@
                zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
                zfice = MIN(MAX(zfice,0.0),1.0)
-	       zrad_liq  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o*pi*rhowater) )
-	       zrad_ice  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o_ice*pi*rhowaterice) )
+               zrad_liq  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o*pi*rhowater) )
+               zrad_ice  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o_ice*pi*rhowaterice) )
                nueffrad(ig,l) = coef_chaud * (1.-zfice) + coef_froid * zfice
                zrad = zrad_liq * (1.-zfice) + zrad_ice * zfice
 
-	       reffrad(ig,l) = min(max(zrad,1.e-6),1000.e-6)
+               reffrad(ig,l) = min(max(zrad,1.e-6),1000.e-6)
                enddo
             enddo      
@@ -213,16 +210,16 @@
 
       if (radfixed) then
-	 reffliq(1:ngrid,1:nlayer)= rad_h2o
-	 reffice(1:ngrid,1:nlayer)= rad_h2o_ice
+         reffliq(1:ngrid,1:nlayer)= rad_h2o
+         reffice(1:ngrid,1:nlayer)= rad_h2o_ice
       else
-	 do k=1,nlayer
-	   do i=1,ngrid
-	     reffliq(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o*pi*rhowater))
-	     reffliq(i,k) = min(max(reffliq(i,k),1.e-6),1000.e-6)
-	   
-	     reffice(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o_ice*pi*rhowaterice))
-	     reffice(i,k) = min(max(reffice(i,k),1.e-6),1000.e-6)
-	   enddo
-	 enddo
+         do k=1,nlayer
+           do i=1,ngrid
+             reffliq(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o*pi*rhowater))
+             reffliq(i,k) = min(max(reffliq(i,k),1.e-6),1000.e-6)
+           
+             reffice(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o_ice*pi*rhowaterice))
+             reffice(i,k) = min(max(reffice(i,k),1.e-6),1000.e-6)
+           enddo
+         enddo
       endif
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90	(revision 1529)
@@ -1,9 +1,8 @@
-      module radinc_h
+module radinc_h
 
-      implicit none
+  implicit none
 
-#include "dimensions.h"
-#include "bands.h"
-#include "scatterers.h"
+  include "bands.h"
+  include "scatterers.h"
 
 !======================================================================
@@ -57,7 +56,8 @@
 !----------------------------------------------------------------------
 
-      integer, parameter :: L_NLAYRAD  = llm
-      integer, parameter :: L_LEVELS   = 2*(llm-1)+3
-      integer, parameter :: L_NLEVRAD  = llm+1
+      integer,save :: L_NLAYRAD  ! = nbp_lev ! set by ini_radinc_h
+      integer,save :: L_LEVELS   ! = 2*(nbp_lev-1)+3 ! set by ini_radinc_h
+      integer,save :: L_NLEVRAD  ! = nbp_lev+1 ! set by ini_radinc_h
+!$OMP THREADPRIVATE(L_NLAYRAD,L_LEVELS,L_NLEVRAD)
 
       ! These are set in sugas_corrk
@@ -96,3 +96,16 @@
 !$OMP THREADPRIVATE(banddir)
 
-      end module radinc_h
+contains
+
+  subroutine ini_radinc_h(nbp_lev)
+  ! Initialize module variables
+  implicit none
+  integer,intent(in) :: nbp_lev
+  
+  L_NLAYRAD = nbp_lev
+  L_LEVELS = 2*(nbp_lev-1)+3
+  L_NLEVRAD = nbp_lev+1
+  
+  end subroutine
+
+end module radinc_h
Index: /trunk/LMDZ.GENERIC/libf/phystd/soil.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/soil.F	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/soil.F	(revision 1529)
@@ -20,8 +20,4 @@
 !        heat capacity are commons in comsoil.h
 !-----------------------------------------------------------------------
-
-!      include "dimensions.h"
-!      include "dimphys.h"
-
 
 c-----------------------------------------------------------------------
Index: /trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90	(revision 1529)
@@ -2,5 +2,5 @@
 
       ! inputs
-      use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,iim,jjm,naerkind
+      use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind
       use radcommon_h, only: blamv,blami,lamrefir,lamrefvis
       use datafile_mod, only: datadir, aerdir
Index: /trunk/LMDZ.GENERIC/libf/phystd/surf_heat_transp_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/surf_heat_transp_mod.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/surf_heat_transp_mod.F90	(revision 1529)
@@ -4,25 +4,123 @@
 MODULE surf_heat_transp_mod
 
-
+IMPLICIT NONE  
+
+  ! Variables copied over from dyn3d dynamics:
+  REAL,SAVE,ALLOCATABLE :: fext(:)
+  REAL,SAVE,ALLOCATABLE :: unsairez(:)
+  REAL,SAVE,ALLOCATABLE :: unsaire(:)
+  REAL,SAVE,ALLOCATABLE :: cu(:)
+  REAL,SAVE,ALLOCATABLE :: cv(:)
+  REAL,SAVE,ALLOCATABLE :: cuvsurcv(:)
+  REAL,SAVE,ALLOCATABLE :: cvusurcu(:)
+  REAL,SAVE,ALLOCATABLE :: aire(:)
+  REAL,SAVE :: apoln
+  REAL,SAVE :: apols
+  REAL,SAVE,ALLOCATABLE :: aireu(:)
+  REAL,SAVE,ALLOCATABLE :: airev(:) 
+
+  LOGICAL,SAVE :: alpha_var
+  LOGICAL,SAVE :: slab_upstream
+  REAL,SAVE,ALLOCATABLE :: zmasqu(:)
+  REAL,SAVE,ALLOCATABLE :: zmasqv(:)
+  REAL,SAVE,ALLOCATABLE :: unsfv(:)
+  REAL,SAVE,ALLOCATABLE :: unsev(:)
+  REAL,SAVE,ALLOCATABLE :: unsfu(:)
+  REAL,SAVE,ALLOCATABLE :: unseu(:)
+
+  ! Routines usable only by routines within this module:
+  PRIVATE :: gr_fi_dyn, gr_dyn_fi
+  ! Routines from dyn3d, valid on global dynamics grid only:
+  PRIVATE :: grad,diverg,gr_u_scal,gr_v_scal
+  
 CONTAINS
-
-      SUBROUTINE divgrad_phy(ngrid,nlevs,temp,delta)
-
-      USE comhdiff_mod, ONLY: zmasqu,zmasqv
-
+  
+  SUBROUTINE ini_surf_heat_transp(ip1jm,ip1jmp1,unsairez_,fext_,unsaire_,&
+                                  cu_,cuvsurcv_,cv_,cvusurcu_, &
+                                  aire_,apoln_,apols_, &
+                                  aireu_,airev_)
+    USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
+    IMPLICIT NONE
+    ! Transfer some variables from dyn3d dynamics
+    INTEGER,INTENT(IN) :: ip1jm
+    INTEGER,INTENT(IN) :: ip1jmp1
+    REAL,INTENT(IN) :: unsairez_(ip1jm)
+    REAL,INTENT(IN) :: fext_(ip1jm)
+    REAL,INTENT(IN) :: unsaire_(ip1jmp1)
+    REAL,INTENT(IN) :: cu_(ip1jmp1)
+    REAL,INTENT(IN) :: cuvsurcv_(ip1jm)
+    REAL,INTENT(IN) :: cv_(ip1jm)
+    REAL,INTENT(IN) :: cvusurcu_(ip1jmp1)
+    REAL,INTENT(IN) :: aire_(ip1jmp1)
+    REAL,INTENT(IN) :: apoln_
+    REAL,INTENT(IN) :: apols_
+    REAL,INTENT(IN) :: aireu_(ip1jmp1)
+    REAL,INTENT(IN) :: airev_(ip1jm)
+    
+    ! Sanity check
+    if ((ip1jm.ne.((nbp_lon+1)*(nbp_lat-1))).or. &
+        (ip1jmp1.ne.((nbp_lon+1)*nbp_lat))) then
+      write(*,*) "ini_surf_heat_transp Error: wrong array sizes"
+      stop
+    endif
+    
+    allocate(unsairez(ip1jm))
+    unsairez(:)=unsairez_(:)
+    allocate(fext(ip1jm))
+    fext(:)=fext_(:)
+    allocate(unsaire(ip1jmp1))
+    unsaire(:)=unsaire_(:)
+    allocate(cu(ip1jmp1))
+    cu(:)=cu_(:)
+    allocate(cuvsurcv(ip1jm))
+    cuvsurcv(:)=cuvsurcv_(:)
+    allocate(cv(ip1jm))
+    cv(:)=cv_(:)
+    allocate(cvusurcu(ip1jmp1))
+    cvusurcu(:)=cvusurcu_(:)
+    allocate(aire(ip1jmp1))
+    aire(:)=aire_(:)
+    apoln=apoln_
+    apols=apols_
+    allocate(aireu(ip1jmp1))
+    aireu(:)=aireu_(:)
+    allocate(airev(ip1jm))
+    airev(:)=airev_(:)
+    
+  END SUBROUTINE ini_surf_heat_transp
+  
+  SUBROUTINE ini_surf_heat_transp_mod
+    USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
+    IMPLICIT NONE
+    INTEGER :: ip1jm, ip1jmp1
+    
+    ip1jm=(nbp_lon+1)*(nbp_lat-1)
+    ip1jmp1=(nbp_lon+1)*nbp_lat
+  
+    allocate(zmasqu(ip1jmp1))
+    allocate(zmasqv(ip1jm))
+    allocate(unsfv(ip1jm))
+    allocate(unsev(ip1jm))
+    allocate(unsfu(ip1jmp1))
+    allocate(unseu(ip1jmp1))
+
+  END SUBROUTINE ini_surf_heat_transp_mod
+
+  SUBROUTINE divgrad_phy(ngrid,nlevs,temp,delta)
+
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
       IMPLICIT NONE
 
-#include "dimensions.h"
-!#include "dimphys.h"
-#include "paramet.h"
-#include "comgeom.h"
-      
       INTEGER,INTENT(IN) :: ngrid, nlevs
       REAL,INTENT(IN) :: temp(ngrid,nlevs)
       REAL,INTENT(OUT) :: delta(ngrid,nlevs)
-      REAL delta_2d(ip1jmp1,nlevs)
+      REAL delta_2d((nbp_lon+1)*nbp_lat,nlevs)
       INTEGER :: ll
-      REAL ghx(ip1jmp1,nlevs), ghy(ip1jm,nlevs)
-
+      REAL ghx((nbp_lon+1)*nbp_lat,nlevs)
+      REAL ghy((nbp_lon+1)*(nbp_lat-1),nlevs)
+      INTEGER :: iip1,jjp1
+      
+      iip1=nbp_lon+1
+      jjp1=nbp_lat
 
       CALL gr_fi_dyn(nlevs,ngrid,iip1,jjp1,temp,delta_2d)
@@ -42,23 +140,23 @@
 
 
-      SUBROUTINE init_masquv(ngrid,zmasq)
-
-      USE comhdiff_mod, ONLY: zmasqu,zmasqv,unsfu,unsfv,unseu,unsev
-     
+  SUBROUTINE init_masquv(ngrid,zmasq)
+
+      USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
       IMPLICIT NONE
-
-#include "dimensions.h"
-!#include "dimphys.h"
-#include "paramet.h"
-#include "comgeom.h"
 
 
       INTEGER,INTENT(IN) :: ngrid
       REAL zmasq(ngrid)
-      REAL zmasq_2d(ip1jmp1)
-      REAL ff(ip1jm)
+      REAL zmasq_2d((nbp_lon+1)*nbp_lat)
+      REAL ff((nbp_lon+1)*(nbp_lat-1))
       REAL eps
       INTEGER i
-
+      INTEGER :: iim,iip1,jjp1,ip1jm,ip1jmp1
+      
+      iim=nbp_lon
+      iip1=nbp_lon+1
+      jjp1=nbp_lat
+      ip1jm=(nbp_lon+1)*(nbp_lat-1)
+      ip1jmp1=(nbp_lon+1)*nbp_lat
 
 ! Masques u,v
@@ -104,27 +202,32 @@
 
 
-      SUBROUTINE slab_ekman2(ngrid,tx_phy,ty_phy,ts_phy,dt_phy)
-  
-          use slab_ice_h 
-          USE comhdiff_mod, ONLY: zmasqu,zmasqv,unsfu,unsfv,unseu,unsev
+  SUBROUTINE slab_ekman2(ngrid,tx_phy,ty_phy,ts_phy,dt_phy)
+  
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
+      USE slab_ice_h, ONLY: noceanmx
 
       IMPLICIT NONE
       
-#include "dimensions.h"
-!#include "dimphys.h"
-#include "paramet.h"
-#include "comgeom.h"
-
       INTEGER,INTENT(IN) :: ngrid
       INTEGER ij
-      REAL txv(ip1jm),fluxm(ip1jm),tyv(ip1jm)
-      REAL fluxtm(ip1jm,noceanmx),fluxtz(ip1jmp1,noceanmx)
-      REAL tyu(ip1jmp1),txu(ip1jmp1),fluxz(ip1jmp1),fluxv(ip1jmp1)
-      REAL dt(ip1jmp1,noceanmx),ts(ip1jmp1,noceanmx)
+      REAL txv((nbp_lon+1)*(nbp_lat-1)),fluxm((nbp_lon+1)*(nbp_lat-1))
+      REAL tyv((nbp_lon+1)*(nbp_lat-1))
+      REAL fluxtm((nbp_lon+1)*(nbp_lat-1),noceanmx)
+      REAL fluxtz((nbp_lon+1)*nbp_lat,noceanmx)
+      REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat)
+      REAL fluxz((nbp_lon+1)*nbp_lat),fluxv((nbp_lon+1)*nbp_lat)
+      REAL dt((nbp_lon+1)*nbp_lat,noceanmx),ts((nbp_lon+1)*nbp_lat,noceanmx)
       REAL tx_phy(ngrid),ty_phy(ngrid)
       REAL dt_phy(ngrid,noceanmx),ts_phy(ngrid,noceanmx)
 
-
-
+      INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
+
+      iim=nbp_lon
+      iip1=nbp_lon+1
+      iip2=nbp_lon+2
+      jjp1=nbp_lon
+      ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+      ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
+      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
 
 ! Passage taux,y sur grilles 2d
@@ -226,5 +329,223 @@
   END SUBROUTINE slab_ekman2
 
- 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
+  ! Transfer a variable on global "physics" grid to global "dynamics" grid
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: im,jm,ngrid,nfield
+  REAL,INTENT(IN) :: pfi(ngrid,nfield)
+  REAL,INTENT(OUT) :: pdyn(im,jm,nfield)
+  
+  INTEGER :: i,j,ifield,ig
+  
+  DO ifield=1,nfield
+    ! Handle poles
+    DO i=1,im
+      pdyn(i,1,ifield)=pfi(1,ifield)
+      pdyn(i,jm,ifield)=pfi(ngrid,ifield)
+    ENDDO
+    ! Other points
+    DO j=2,jm-1
+      ig=2+(j-2)*(im-1)
+      CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
+      pdyn(im,j,ifield)=pdyn(1,j,ifield)
+    ENDDO
+  ENDDO ! of DO ifield=1,nfield
+  
+  END SUBROUTINE gr_fi_dyn
+  
+
+
+  SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
+  ! Transfer a variable on global "dynamics" grid to global "physics" grid
+  IMPLICIT NONE
+  
+  INTEGER,INTENT(IN) :: im,jm,ngrid,nfield
+  REAL,INTENT(IN) :: pdyn(im,jm,nfield)
+  REAL,INTENT(OUT) :: pfi(ngrid,nfield)
+  
+  INTEGER j,ifield,ig
+
+  ! Sanity check:
+  IF(ngrid.NE.2+(jm-2)*(im-1)) THEN
+    WRITE(*,*) "gr_dyn_fi error, wrong sizes"
+    STOP
+  ENDIF
+  
+  ! Handle poles
+  CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
+  CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
+  ! Other points
+  DO ifield=1,nfield
+    DO j=2,jm-1
+      ig=2+(j-2)*(im-1)
+      CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
+    ENDDO
+  ENDDO
+
+  END SUBROUTINE gr_dyn_fi
+
+
+
+  SUBROUTINE  grad(klevel,pg,pgx,pgy)
+  ! compute the covariant components x,y of the gradient of pg
+  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
+  IMPLICIT NONE
+  
+  INTEGER,INTENT(IN) :: klevel
+  REAL,INTENT(IN) :: pg((nbp_lon+1)*nbp_lat,klevel)
+  REAL,INTENT(OUT) :: pgx((nbp_lon+1)*nbp_lat,klevel)
+  REAL,INTENT(OUT) :: pgy((nbp_lon+1)*(nbp_lat-1),klevel)
+
+  INTEGER :: l,ij
+  INTEGER :: iim,iip1,ip1jm,ip1jmp1
+  
+  iim=nbp_lon
+  iip1=nbp_lon+1
+  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+  
+  DO l=1,klevel
+    DO ij=1,ip1jmp1-1
+      pgx(ij,l)=pg(ij+1,l)-pg(ij,l)
+    ENDDO
+    ! correction for pgx(ip1,j,l) ...
+    ! ... pgx(iip1,j,l)=pgx(1,j,l) ...
+    DO ij=iip1,ip1jmp1,iip1
+      pgx(ij,l)=pgx(ij-iim,l)
+    ENDDO
+    DO ij=1,ip1jm
+      pgy(ij,l)=pg(ij,l)-pg(ij+iip1,l)
+    ENDDO
+  ENDDO
+  
+  END SUBROUTINE grad
+
+
+
+  SUBROUTINE diverg(klevel,x,y,div)
+  ! compute the divergence of a vector field of components
+  ! x,y. y and y being covriant components
+  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
+  IMPLICIT NONE
+  
+  INTEGER,INTENT(IN) :: klevel
+  REAL,INTENT(IN) :: x((nbp_lon+1)*nbp_lat,klevel)
+  REAL,INTENT(IN) :: y((nbp_lon+1)*(nbp_lat-1),klevel)
+  REAL,INTENT(OUT) :: div((nbp_lon+1)*nbp_lat,klevel)
+  
+  INTEGER :: l,ij
+  INTEGER :: iim,iip1,iip2,ip1jm,ip1jmp1,ip1jmi1
+  
+  REAL :: aiy1(nbp_lon+1),aiy2(nbp_lon+1)
+  REAL :: sumypn,sumyps
+  REAL,EXTERNAL :: SSUM
+  
+  iim=nbp_lon
+  iip1=nbp_lon+1
+  iip2=nbp_lon+2
+  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+  ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
+  
+  DO l=1,klevel
+    DO ij=iip2,ip1jm-1
+      div(ij+1,l)= &
+        cvusurcu(ij+1)*x(ij+1,l)-cvusurcu(ij)*x(ij,l)+ &
+        cuvsurcv(ij-iim)*y(ij-iim,l)-cuvsurcv(ij+1)*y(ij+1,l) 
+    ENDDO
+    ! correction for div(1,j,l) ...
+    ! ... div(1,j,l)= div(iip1,j,l) ...
+    DO ij=iip2,ip1jm,iip1
+      div(ij,l)=div(ij+iim,l)
+    ENDDO
+    ! at the poles
+    DO ij=1,iim
+      aiy1(ij)=cuvsurcv(ij)*y(ij,l)
+      aiy2(ij)=cuvsurcv(ij+ip1jmi1)*y(ij+ip1jmi1,l)
+    ENDDO
+    sumypn=SSUM(iim,aiy1,1)/apoln
+    sumyps=SSUM(iim,aiy2,1)/apols
+    DO ij=1,iip1
+      div(ij,l)=-sumypn
+      div(ij+ip1jm,l)=sumyps
+    ENDDO
+  ENDDO ! of DO l=1,klevel
+  
+  !!! CALL filtreg( div, jjp1, klevel, 2, 2, .TRUE., 1 )
+  
+  DO l=1,klevel
+    DO ij=iip2,ip1jm
+      div(ij,l)=div(ij,l)*unsaire(ij) 
+    ENDDO
+  ENDDO
+  
+  END SUBROUTINE diverg
+
+
+
+  SUBROUTINE gr_u_scal(nx,x_u,x_scal)
+  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
+  IMPLICIT NONE
+  
+  INTEGER,INTENT(IN) :: nx
+  REAL,INTENT(IN) :: x_u((nbp_lon+1)*nbp_lat,nx)
+  REAL,INTENT(OUT) :: x_scal((nbp_lon+1)*nbp_lat,nx)
+  
+  INTEGER :: l,ij
+  INTEGER :: iip1,jjp1,ip1jmp1
+  
+  iip1=nbp_lon+1
+  jjp1=nbp_lat
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+  
+  DO l=1,nx
+    DO ij=ip1jmp1,2,-1
+      x_scal(ij,l)= &
+                   (aireu(ij)*x_u(ij,l)+aireu(ij-1)*x_u(ij-1,l)) &
+                  /(aireu(ij)+aireu(ij-1))
+    ENDDO
+  ENDDO
+  
+  CALL SCOPY(nx*jjp1,x_scal(iip1,1),iip1,x_scal(1,1),iip1)
+
+  END SUBROUTINE gr_u_scal
+
+
+
+  SUBROUTINE gr_v_scal(nx,x_v,x_scal)
+  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
+  IMPLICIT NONE
+  
+  INTEGER,INTENT(IN) :: nx
+  REAL,INTENT(IN) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
+  REAL,INTENT(OUT) :: x_scal((nbp_lon+1)*nbp_lat,nx)
+  
+  INTEGER :: l,ij
+  INTEGER :: iip1,iip2,ip1jm,ip1jmp1
+  
+  iip1=nbp_lon+1
+  iip2=nbp_lon+2
+  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+  
+  DO l=1,nx
+    DO ij=iip2,ip1jm
+      x_scal(ij,l)= &
+                   (airev(ij-iip1)*x_v(ij-iip1,l)+airev(ij)*x_v(ij,l)) &
+                  /(airev(ij-iip1)+airev(ij))
+    ENDDO
+    DO ij=1,iip1
+      x_scal(ij,l)=0.
+    ENDDO
+    DO ij=ip1jm+1,ip1jmp1
+      x_scal(ij,l)=0.
+    ENDDO
+  ENDDO
+
+  END SUBROUTINE gr_v_scal
+  
 END MODULE surf_heat_transp_mod
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/vdif_kc.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/vdif_kc.F	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/vdif_kc.F	(revision 1529)
@@ -1,7 +1,4 @@
       SUBROUTINE vdif_kc(ngrid,nlay,dt,g,zlev,zlay,u,v,teta,cd,q2,km,kn)
       IMPLICIT NONE
-c.......................................................................
-!#include "dimensions.h"
-!#include "dimphys.h"
 c.......................................................................
 c
Index: /trunk/LMDZ.GENERIC/libf/phystd/vlz_fi.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/vlz_fi.F	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/vlz_fi.F	(revision 1529)
@@ -15,8 +15,4 @@
 c   --------------------------------------------------------------------
       IMPLICIT NONE
-c
-!#include "dimensions.h"
-!#include "dimphys.h"
-
 c
 c
Index: /trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/writediagfi.F	(revision 1529)
@@ -40,14 +40,13 @@
 !=================================================================
       use surfdat_h, only: phisfi
+      use comgeomphy, only: airephy
       use time_phylmdz_mod, only: ecritphy, day_step, iphysiq, day_ini
       USE mod_phys_lmdz_para, only : is_parallel, 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, nbp_lev
       implicit none
 
 ! Commons
-      include "dimensions.h"
-      include "paramet.h"
-      include "comgeom.h"
       include "netcdf.inc"
 
@@ -56,11 +55,11 @@
       character (len=*),intent(in) :: nom,titre,unite
       integer,intent(in) :: dim
-      real,intent(in) :: px(ngrid,llm)
+      real,intent(in) :: px(ngrid,nbp_lev)
 
 ! Local variables:
 
-      real*4 dx3(iip1,jjp1,llm) ! to store a 3D data set
-      real*4 dx2(iip1,jjp1)     ! to store a 2D (surface) data set
-      real*4 dx1(llm)           ! to store a 1D (column) data set
+      real*4 dx3(nbp_lon+1,nbp_lat,nbp_lev) ! to store a 3D data set
+      real*4 dx2(nbp_lon+1,nbp_lat)     ! to store a 2D (surface) data set
+      real*4 dx1(nbp_lev)           ! to store a 1D (column) data set
       real*4 dx0
 
@@ -68,10 +67,10 @@
 !$OMP THREADPRIVATE(date)
 
-      REAL phis(ip1jmp1)
+      REAL phis((nbp_lon+1),nbp_lat)
+      REAL area((nbp_lon+1),nbp_lat)
 
       integer irythme
       integer ierr,ierr2
-      integer iq
-      integer i,j,l,zmax , ig0
+      integer i,j,l, ig0
 
       integer,save :: zitau=0
@@ -102,14 +101,16 @@
 #ifdef CPP_PARA
 ! Added to work in parallel mode
-      real dx3_glop(klon_glo,llm)
-      real dx3_glo(iim,jjp1,llm) ! to store a global 3D data set
+      real dx3_glop(klon_glo,nbp_lev)
+      real dx3_glo(nbp_lon,nbp_lat,nbp_lev) ! 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)
-!      real dx1_glo(llm)          ! to store a 1D (column) data set
+!      real dx1_glo(nbp_lev)          ! to store a 1D (column) data set
 !      real dx0_glo
       real phisfi_glo(klon_glo) ! surface geopotential on global physics grid
+      real areafi_glo(klon_glo) ! mesh area on global physics grid
 #else
       real phisfi_glo(ngrid) ! surface geopotential on global physics grid
+      real areafi_glo(ngrid) ! mesh area on global physics grid
 #endif
 
@@ -117,5 +118,5 @@
 !Sortie des variables au rythme voulu
 
-      irythme = ecritphy ! sortie au rythme de ecritphy
+      irythme = int(ecritphy) ! sortie au rythme de ecritphy
 !     irythme = iconser  ! sortie au rythme des variables de controle
 !     irythme = iphysiq  ! sortie a tous les pas physique
@@ -188,6 +189,9 @@
           ! Gather phisfi() geopotential on physics grid
           call Gather(phisfi,phisfi_glo)
+          ! Gather airephy() mesh area on physics grid
+          call Gather(airephy,areafi_glo)
 #else
          phisfi_glo(:)=phisfi(:)
+         areafi_glo(:)=airephy(:)
 #endif
 
@@ -216,7 +220,25 @@
          ierr = NF_ENDDEF(nid)
 
+         ! Build phis() and area()
+         do i=1,nbp_lon+1 ! poles
+           phis(i,1)=phisfi_glo(1)
+           phis(i,nbp_lat)=phisfi_glo(klon_glo)
+           ! 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
+              phis(i,j)=phisfi_glo(ig0+i)
+              area(i,j)=areafi_glo(ig0+i)
+           enddo
+           ! handle redundant point in longitude
+           phis(nbp_lon+1,j)=phis(1,j)
+           area(nbp_lon+1,j)=area(1,j)
+         enddo
+         
          ! write "header" of file (longitudes, latitudes, geopotential, ...)
-         call gr_fi_dyn(1,size(phisfi_glo),iip1,jjp1,phisfi_glo,phis)
-         call iniwrite(nid,day_ini,phis)
+         call iniwrite(nid,day_ini,phis,area)
 
          endif ! of if (is_master)
@@ -234,7 +256,8 @@
 
       if (ngrid.eq.1) then
-        ! in testphys1d, for the 1d version of the GCM, iphysiq 
-        ! should be most likely 1 (because no dyn!)
+        ! 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
@@ -296,6 +319,6 @@
             call Grid1Dto2D_glo(dx3_glop,dx3_glo)
             ! copy dx3_glo() to dx3(:) and add redundant longitude
-            dx3(1:iim,:,:)=dx3_glo(1:iim,:,:)
-            dx3(iip1,:,:)=dx3(1,:,:)
+            dx3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
+            dx3(nbp_lon+1,:,:)=dx3(1,:,:)
           endif
 !$OMP END MASTER
@@ -304,15 +327,15 @@
 !         Passage variable physique -->  variable dynamique
 !         recast (copy) variable from physics grid to dynamics grid
-           DO l=1,llm
-             DO i=1,iip1
+           DO l=1,nbp_lev
+             DO i=1,nbp_lon+1
                 dx3(i,1,l)=px(1,l)
-                dx3(i,jjp1,l)=px(ngrid,l)
+                dx3(i,nbp_lat,l)=px(ngrid,l)
              ENDDO
-             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
                    dx3(i,j,l)=px(ig0+i,l)
                 ENDDO
-                dx3(iip1,j,l)=dx3(1,j,l)
+                dx3(nbp_lon+1,j,l)=dx3(1,j,l)
              ENDDO
            ENDDO
@@ -344,7 +367,7 @@
            corner(4)=ntime
 
-           edges(1)=iip1
-           edges(2)=jjp1
-           edges(3)=llm
+           edges(1)=nbp_lon+1
+           edges(2)=nbp_lat
+           edges(3)=nbp_lev
            edges(4)=1
 !#ifdef NC_DOUBLE
@@ -380,6 +403,6 @@
             call Grid1Dto2D_glo(dx2_glop,dx2_glo)
             ! copy dx2_glo() to dx2(:) and add redundant longitude
-            dx2(1:iim,:)=dx2_glo(1:iim,:)
-            dx2(iip1,:)=dx2(1,:)
+            dx2(1:nbp_lon,:)=dx2_glo(1:nbp_lon,:)
+            dx2(nbp_lon+1,:)=dx2(1,:)
           endif
 !$OMP END MASTER
@@ -390,14 +413,14 @@
 !         recast (copy) variable from physics grid to dynamics grid
 
-             DO i=1,iip1
+             DO i=1,nbp_lon+1
                 dx2(i,1)=px(1,1)
-                dx2(i,jjp1)=px(ngrid,1)
+                dx2(i,nbp_lat)=px(ngrid,1)
              ENDDO
-             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
                    dx2(i,j)=px(ig0+i,1)
                 ENDDO
-                dx2(iip1,j)=dx2(1,j)
+                dx2(nbp_lon+1,j)=dx2(1,j)
              ENDDO
 #endif
@@ -426,6 +449,6 @@
            corner(2)=1
            corner(3)=ntime
-           edges(1)=iip1
-           edges(2)=jjp1
+           edges(1)=nbp_lon+1
+           edges(2)=nbp_lat
            edges(3)=1
 
@@ -457,5 +480,5 @@
 !         Passage variable physique -->  physique dynamique
 !         recast (copy) variable from physics grid to dynamics grid
-          do l=1,llm
+          do l=1,nbp_lev
             dx1(l)=px(1,l)
           enddo
@@ -479,5 +502,5 @@
            corner(2)=ntime
            
-           edges(1)=llm
+           edges(1)=nbp_lev
            edges(2)=1
 !#ifdef NC_DOUBLE
@@ -543,3 +566,4 @@
 
 #endif
+! of #ifndef MESOSCALE
       end
Index: /trunk/LMDZ.GENERIC/libf/phystd/writediagsoil.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/writediagsoil.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/writediagsoil.F90	(revision 1529)
@@ -12,14 +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"netcdf.inc"
+include"netcdf.inc"
 
 ! Arguments:
@@ -31,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
@@ -42,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
@@ -60,7 +65,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
@@ -81,5 +86,5 @@
   
   ! Set output sample rate
-  isample=ecritphy ! same as for diagfi outputs
+  isample=int(ecritphy) ! same as for diagfi outputs
   ! Note ecritphy is known from control.h
   
@@ -91,8 +96,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
@@ -148,6 +183,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
@@ -156,15 +191,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
@@ -195,6 +230,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
@@ -223,6 +258,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
@@ -230,15 +265,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
@@ -266,6 +301,6 @@
   corners(3)=ntime
   
-  edges(1)=iip1
-  edges(2)=jjp1
+  edges(1)=nbp_lon+1
+  edges(2)=nbp_lat
   edges(3)=1
   
Index: /trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/writediagspecIR.F	(revision 1529)
@@ -43,9 +43,8 @@
 ! Addition by RW (2010) to allow OLR to be saved in .nc format
       use radinc_h, only : L_NSPECTI
-!      USE surfdat_h, only : phisfi
-#ifdef CPP_PARA
+      use comgeomphy, only: airephy
       use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
-      use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
-#endif
+      use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo,
+     &                              nbp_lon, nbp_lat
       use time_phylmdz_mod, only: ecritphy, iphysiq, day_step, day_ini
       use callkeys_mod, only: iradia
@@ -53,11 +52,5 @@
       implicit none
 
-! Commons
-#include "dimensions.h"
-!#include "dimphys.h"
-#include "paramet.h"
-!#include "control.h"
-#include "comgeom.h"
-#include "netcdf.inc"
+      include "netcdf.inc"
 
 ! Arguments on input:
@@ -73,8 +66,4 @@
 !      real dx0
 
-      real date
-
-!      REAL phis(ip1jmp1)
-
       integer irythme
       integer ierr
@@ -82,11 +71,7 @@
       integer i,j,l,zmax , ig0
 
-      integer zitau
-      character firstnom*20
-      SAVE firstnom
-      SAVE zitau
-      SAVE date
-      data firstnom /'1234567890'/
-      data zitau /0/
+      integer,save :: zitau=0
+      character(len=20),save :: firstnom='1234567890'
+      real,save :: date
 !$OMP THREADPRIVATE(firstnom,zitau,date)
 
@@ -100,14 +85,15 @@
       integer, dimension(4) :: edges,corner
 
+      real area((nbp_lon+1),nbp_lat)
 ! added by RDW for OLR output
-       real dx3(iip1,jjp1,L_NSPECTI) ! to store the data set
+       real dx3(nbp_lon+1,nbp_lat,L_NSPECTI) ! to store the data set
 
 #ifdef CPP_PARA
 ! Added to work in parallel mode
       real dx3_glop(klon_glo,L_NSPECTI)
-      real dx3_glo(iim,jjp1,L_NSPECTI) ! to store a global 3D data set
-#else
-      logical,parameter :: is_master=.true.
-      logical,parameter :: is_mpi_root=.true.
+      real dx3_glo(nbp_lon,nbp_lat,L_NSPECTI) ! to store a global 3D data set
+      real areafi_glo(klon_glo) ! mesh area on global physics grid
+#else
+      real areafi_glo(ngrid) ! mesh area on global physics grid
 #endif
 
@@ -139,4 +125,10 @@
          endif
 
+#ifdef CPP_PARA
+          ! Gather airephy() mesh area on physics grid
+          call Gather(airephy,areafi_glo)
+#else
+         areafi_glo(:)=airephy(:)
+#endif
          ! Create the NetCDF file
          if (is_master) then
@@ -159,7 +151,21 @@
          ierr = NF_ENDDEF(nid)
 
-!         call gr_fi_dyn(1,size(phisfi_glo),iip1,jjp1,phisfi_glo,phis)
+         ! Build area()
+         do i=1,nbp_lon+1 ! poles
+           ! 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
+              area(i,j)=areafi_glo(ig0+i)
+           enddo
+           ! handle redundant point in longitude
+           area(nbp_lon+1,j)=area(1,j)
+         enddo
+
          ! write "header" of file (longitudes, latitudes, area, ...)
-         call iniwrite_specIR(nid,day_ini)
+         call iniwrite_specIR(nid,day_ini,area)
          endif ! of if (is_master)
 
@@ -215,5 +221,5 @@
              endif
 
-             write(6,*)'WRITEDIAGSPEC: date= ', date
+             write(6,*)'WRITEDIAGSPECIR: date= ', date
            endif ! of if (is_master)
         end if ! of if (nom.eq.firstnom)
@@ -233,6 +239,6 @@
             call Grid1Dto2D_glo(dx3_glop,dx3_glo)
             ! copy dx3_glo() to dx3(:) and add redundant longitude
-            dx3(1:iim,:,:)=dx3_glo(1:iim,:,:)
-            dx3(iip1,:,:)=dx3(1,:,:)
+            dx3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
+            dx3(nbp_lon+1,:,:)=dx3(1,:,:)
           endif
 !$OMP END MASTER
@@ -240,14 +246,14 @@
 #else
            DO l=1,L_NSPECTI
-             DO i=1,iip1
+             DO i=1,nbp_lon+1
                 dx3(i,1,l)=px(1,l)
-                dx3(i,jjp1,l)=px(ngrid,l)
+                dx3(i,nbp_lat,l)=px(ngrid,l)
              ENDDO
-             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
                    dx3(i,j,l)=px(ig0+i,l)
                 ENDDO
-                dx3(iip1,j,l)=dx3(1,j,l)
+                dx3(nbp_lon+1,j,l)=dx3(1,j,l)
              ENDDO
            ENDDO
@@ -279,6 +285,6 @@
            corner(4)=ntime
 
-           edges(1)=iip1
-           edges(2)=jjp1
+           edges(1)=nbp_lon+1
+           edges(2)=nbp_lat
            edges(3)=L_NSPECTI
            edges(4)=1
Index: /trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/writediagspecVI.F	(revision 1529)
@@ -43,9 +43,8 @@
 ! Addition by RW (2010) to allow OSR to be saved in .nc format
       use radinc_h, only : L_NSPECTV
-!      USE surfdat_h
-#ifdef CPP_PARA
+      use comgeomphy, only: airephy
       use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
-      use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
-#endif
+      use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo,
+     &                              nbp_lon, nbp_lat
       use time_phylmdz_mod, only: ecritphy, iphysiq, day_step, day_ini
       use callkeys_mod, only: iradia
@@ -53,11 +52,5 @@
       implicit none
 
-! Commons
-#include "dimensions.h"
-!#include "dimphys.h"
-#include "paramet.h"
-!#include "control.h"
-#include "comgeom.h"
-#include "netcdf.inc"
+      include "netcdf.inc"
 
 ! Arguments on input:
@@ -73,8 +66,4 @@
 !      real dx0
 
-      real date
-
-!      REAL phis(ip1jmp1)
-
       integer irythme
       integer ierr
@@ -82,11 +71,7 @@
       integer i,j,l,zmax , ig0
 
-      integer zitau
-      character firstnom*20
-      SAVE firstnom
-      SAVE zitau
-      SAVE date
-      data firstnom /'1234567890'/
-      data zitau /0/
+      integer,save :: zitau=0
+      character(len=20),save :: firstnom='1234567890'
+      real,save :: date
 !$OMP THREADPRIVATE(firstnom,zitau,date)
 
@@ -100,14 +85,15 @@
       integer, dimension(4) :: edges,corner
 
+      real area((nbp_lon+1),nbp_lat)
 ! added by RDW for OSR output
-       real dx3(iip1,jjp1,L_NSPECTV) ! to store the data set
+       real dx3(nbp_lon+1,nbp_lat,L_NSPECTV) ! to store the data set
 
 #ifdef CPP_PARA
 ! Added to work in parallel mode
       real dx3_glop(klon_glo,L_NSPECTV)
-      real dx3_glo(iim,jjp1,L_NSPECTV) ! to store a global 3D data set
-#else
-      logical,parameter :: is_master=.true.
-      logical,parameter :: is_mpi_root=.true.
+      real dx3_glo(nbp_lon,nbp_lat,L_NSPECTV) ! to store a global 3D data set
+      real areafi_glo(klon_glo) ! mesh area on global physics grid
+#else
+      real areafi_glo(ngrid) ! mesh area on global physics grid
 #endif
 
@@ -138,4 +124,10 @@
          endif
 
+#ifdef CPP_PARA
+          ! Gather airephy() mesh area on physics grid
+          call Gather(airephy,areafi_glo)
+#else
+         areafi_glo(:)=airephy(:)
+#endif
          ! Create the NetCDF file
          if (is_master) then
@@ -158,8 +150,21 @@
          ierr = NF_ENDDEF(nid)
 
+         ! Build area()
+         do i=1,nbp_lon+1 ! poles
+           ! 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
+              area(i,j)=areafi_glo(ig0+i)
+           enddo
+           ! handle redundant point in longitude
+           area(nbp_lon+1,j)=area(1,j)
+         enddo
+
          ! write "header" of file (longitudes, latitudes, geopotential, ...)
-!         call gr_fi_dyn(1,ngrid,iip1,jjp1,phisfi,phis)
-!         call iniwrite(nid,day_ini,phis)
-         call iniwrite_specVI(nid,day_ini)
+         call iniwrite_specVI(nid,day_ini,area)
          endif ! of if (is_master)
 
@@ -215,5 +220,5 @@
              endif
 
-             write(6,*)'WRITEDIAGSPEC: date= ', date
+             write(6,*)'WRITEDIAGSPECVI: date= ', date
            endif ! of if (is_master)
         end if ! of if (nom.eq.firstnom)
@@ -233,6 +238,6 @@
             call Grid1Dto2D_glo(dx3_glop,dx3_glo)
             ! copy dx3_glo() to dx3(:) and add redundant longitude
-            dx3(1:iim,:,:)=dx3_glo(1:iim,:,:)
-            dx3(iip1,:,:)=dx3(1,:,:)
+            dx3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
+            dx3(nbp_lon+1,:,:)=dx3(1,:,:)
           endif
 !$OMP END MASTER
@@ -240,14 +245,14 @@
 #else
            DO l=1,L_NSPECTV
-             DO i=1,iip1
+             DO i=1,nbp_lon+1
                 dx3(i,1,l)=px(1,l)
-                dx3(i,jjp1,l)=px(ngrid,l)
+                dx3(i,nbp_lat,l)=px(ngrid,l)
              ENDDO
-             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
                    dx3(i,j,l)=px(ig0+i,l)
                 ENDDO
-                dx3(iip1,j,l)=dx3(1,j,l)
+                dx3(nbp_lon+1,j,l)=dx3(1,j,l)
              ENDDO
            ENDDO
@@ -279,6 +284,6 @@
            corner(4)=ntime
 
-           edges(1)=iip1
-           edges(2)=jjp1
+           edges(1)=nbp_lon+1
+           edges(2)=nbp_lat
            edges(3)=L_NSPECTV
            edges(4)=1
Index: /trunk/LMDZ.GENERIC/libf/phystd/wstats.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/wstats.F90	(revision 1528)
+++ /trunk/LMDZ.GENERIC/libf/phystd/wstats.F90	(revision 1529)
@@ -1,22 +1,18 @@
 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
-use statto_mod, only: istats,istime,count
-
+use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo, &
+                              nbp_lon, nbp_lat, nbp_lev
 implicit none
 
-#include "dimensions.h"
-!#include "dimphys.h"
-#include "netcdf.inc"
+include "netcdf.inc"
 
 integer,intent(in) :: ngrid
 character (len=*),intent(in) :: nom,titre,unite
 integer,intent(in) :: dim
-integer,parameter :: iip1=iim+1
-integer,parameter :: jjp1=jjm+1
-real,intent(in) :: px(ngrid,llm)
-real, dimension(iip1,jjp1,llm) :: mean3d,sd3d,dx3
-real, dimension(iip1,jjp1) :: mean2d,sd2d,dx2
+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
 character (len=50) :: namebis
 character (len=50), save :: firstvar
@@ -34,16 +30,16 @@
 ! Added to work in parallel mode
 #ifdef CPP_PARA
-real px3_glop(klon_glo,llm) ! to store a 3D data set on global physics grid
-real px3_glo(iim,jjp1,llm) ! to store a global 3D data set on global lonxlat grid
+real px3_glop(klon_glo,nbp_lev) ! to store a 3D data set on global physics grid
+real px3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set on global lonxlat grid
 real px2_glop(klon_glo) ! to store a 2D data set on global physics grid
-real px2_glo(iim,jjp1) ! to store a 2D data set on global lonxlat grid
+real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
 real px2(ngrid)
-real px3(ngrid,llm)
+real px3(ngrid,nbp_lev)
 #else
 ! When not running in parallel mode:
-real px3_glop(ngrid,llm) ! to store a 3D data set on global physics grid
-real px3_glo(iim,jjp1,llm) ! to store a global 3D data set on global lonxlat grid
+real px3_glop(ngrid,nbp_lev) ! to store a 3D data set on global physics grid
+real px3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set on global lonxlat grid
 real px2_glop(ngrid) ! to store a 2D data set on global physics grid
-real px2_glo(iim,jjp1) ! to store a 2D data set on global lonxlat grid
+real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
 #endif
 
@@ -67,5 +63,5 @@
 #ifdef CPP_PARA
  if (dim.eq.3) then
-  px3(1:ngrid,1:llm)=px(1:ngrid,1:llm)
+  px3(1:ngrid,1:nbp_lev)=px(1:ngrid,1:nbp_lev)
   ! Gather fieds on a "global" (without redundant longitude) array
   call Gather(px3,px3_glop)
@@ -74,6 +70,6 @@
     call Grid1Dto2D_glo(px3_glop,px3_glo)
     ! copy dx3_glo() to dx3(:) and add redundant longitude
-    dx3(1:iim,:,:)=px3_glo(1:iim,:,:)
-    dx3(iip1,:,:)=dx3(1,:,:)
+    dx3(1:nbp_lon,:,:)=px3_glo(1:nbp_lon,:,:)
+    dx3(nbp_lon+1,:,:)=dx3(1,:,:)
   endif
 !$OMP END MASTER
@@ -87,6 +83,6 @@
             call Grid1Dto2D_glo(px2_glop,px2_glo)
             ! copy px2_glo() to dx2(:) and add redundant longitude
-            dx2(1:iim,:)=px2_glo(1:iim,:)
-            dx2(iip1,:)=dx2(1,:)
+            dx2(1:nbp_lon,:)=px2_glo(1:nbp_lon,:)
+            dx2(nbp_lon+1,:)=dx2(1,:)
           endif
 !$OMP END MASTER
@@ -95,14 +91,14 @@
 #else
   if (dim.eq.3) then
-    px3_glop(:,1:llm)=px(:,1:llm)
+    px3_glop(:,1:nbp_lev)=px(:,1:nbp_lev)
 !  Passage variable physique -->  variable dynamique
-    DO l=1,llm
-      DO i=1,iim
+    DO l=1,nbp_lev
+      DO i=1,nbp_lon
          px3_glo(i,1,l)=px(1,l)
-         px3_glo(i,jjp1,l)=px(ngrid,l)
+         px3_glo(i,nbp_lat,l)=px(ngrid,l)
       ENDDO
-      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
             px3_glo(i,j,l)=px(ig0+i,l)
          ENDDO
@@ -112,11 +108,11 @@
     px2_glop(:)=px(:,1)
 !    Passage variable physique -->  physique dynamique
-   DO i=1,iim
+   DO i=1,nbp_lon
      px2_glo(i,1)=px(1,1)
-     px2_glo(i,jjp1)=px(ngrid,1)
+     px2_glo(i,nbp_lat)=px(ngrid,1)
    ENDDO
-   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
         px2_glo(i,j)=px(ig0+i,1)
      ENDDO
@@ -198,10 +194,10 @@
    if (dim.eq.3) then
       start=(/1,1,1,indx/)
-      sizes=(/iip1,jjp1,llm,1/)
+      sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
       mean3d(:,:,:)=0
       sd3d(:,:,:)=0
    else if (dim.eq.2) then
       start=(/1,1,indx,0/)
-      sizes=(/iip1,jjp1,1,0/)
+      sizes=(/nbp_lon+1,nbp_lev,1,0/)
       mean2d(:,:)=0
       sd2d(:,:)=0
@@ -211,5 +207,5 @@
    if (dim.eq.3) then
       start=(/1,1,1,indx/)
-      sizes=(/iip1,jjp1,llm,1/)
+      sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
 #ifdef NC_DOUBLE
       ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
@@ -226,5 +222,5 @@
    else if (dim.eq.2) then
       start=(/1,1,indx,0/)
-      sizes=(/iip1,jjp1,1,0/)
+      sizes=(/nbp_lon+1,nbp_lat,1,0/)
 #ifdef NC_DOUBLE
       ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
@@ -245,9 +241,9 @@
 
 if (dim.eq.3) then
-  dx3(1:iim,:,:)=px3_glo(:,:,:)
-  dx3(iip1,:,:)=dx3(1,:,:)
+  dx3(1:nbp_lon,:,:)=px3_glo(:,:,:)
+  dx3(nbp_lon+1,:,:)=dx3(1,:,:)
 else ! dim.eq.2
-  dx2(1:iim,:)=px2_glo(:,:)
-  dx2(iip1,:)=dx2(1,:)
+  dx2(1:nbp_lon,:)=px2_glo(:,:)
+  dx2(nbp_lon+1,:)=dx2(1,:)
 endif
 
@@ -290,41 +286,37 @@
 !======================================================
 subroutine inivar(nid,varid,ngrid,dim,indx,px,ierr)
+use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev
 
 implicit none
 
-include "dimensions.h"
-!include "dimphys.h"
 include "netcdf.inc"
 
 integer, intent(in) :: nid,varid,dim,indx,ngrid
-real, dimension(ngrid,llm), intent(in) :: px
+real, dimension(ngrid,nbp_lev), intent(in) :: px
 integer, intent(out) :: ierr
-
-integer,parameter :: iip1=iim+1
-integer,parameter :: jjp1=jjm+1
 
 integer :: l,i,j,ig0
 integer, dimension(4) :: start,sizes
-real, dimension(iip1,jjp1,llm) :: dx3
-real, dimension(iip1,jjp1) :: dx2
+real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: dx3
+real, dimension(nbp_lon+1,nbp_lat) :: dx2
 
 if (dim.eq.3) then
 
    start=(/1,1,1,indx/)
-   sizes=(/iip1,jjp1,llm,1/)
+   sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
 
 !  Passage variable physique -->  variable dynamique
 
-   DO l=1,llm
-      DO i=1,iip1
+   DO l=1,nbp_lev
+      DO i=1,nbp_lon+1
          dx3(i,1,l)=px(1,l)
-         dx3(i,jjp1,l)=px(ngrid,l)
+         dx3(i,nbp_lat,l)=px(ngrid,l)
       ENDDO
-      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
             dx3(i,j,l)=px(ig0+i,l)
          ENDDO
-         dx3(iip1,j,l)=dx3(1,j,l)
+         dx3(nbp_lon+1,j,l)=dx3(1,j,l)
       ENDDO
    ENDDO
@@ -339,18 +331,18 @@
 
       start=(/1,1,indx,0/)
-      sizes=(/iip1,jjp1,1,0/)
+      sizes=(/nbp_lon+1,nbp_lat,1,0/)
 
 !    Passage variable physique -->  physique dynamique
 
-  DO i=1,iip1
+  DO i=1,nbp_lon+1
      dx2(i,1)=px(1,1)
-     dx2(i,jjp1)=px(ngrid,1)
+     dx2(i,nbp_lat)=px(ngrid,1)
   ENDDO
-  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
         dx2(i,j)=px(ig0+i,1)
      ENDDO
-     dx2(iip1,j)=dx2(1,j)
+     dx2(nbp_lon+1,j)=dx2(1,j)
   ENDDO
 
@@ -380,5 +372,5 @@
 implicit none
 
-#include "netcdf.inc"
+include "netcdf.inc"
 
 integer,intent(in) :: nid ! NetCDF file ID
