Index: dynamico_lmdz/aquaplanet/LMDZ5/libf/dyn3d_common/conf_dat_m.F90
===================================================================
--- dynamico_lmdz/aquaplanet/LMDZ5/libf/dyn3d_common/conf_dat_m.F90	(revision 3883)
+++ dynamico_lmdz/aquaplanet/LMDZ5/libf/dyn3d_common/conf_dat_m.F90	(revision 3883)
@@ -0,0 +1,235 @@
+MODULE conf_dat_m
+!
+!*******************************************************************************
+
+  PRIVATE
+  PUBLIC :: conf_dat2d, conf_dat3d
+
+
+CONTAINS
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE conf_dat2d(title, xd, yd, xf, yf, champd, interbar)
+!
+!-------------------------------------------------------------------------------
+! Author: P. Le Van
+!-------------------------------------------------------------------------------
+! Purpose: Configure the 2D data field "champd" so that:
+!     - Longitudes are in [ -pi     pi   ]
+!     - Latitudes  are in [  pi/2. -pi/2 ]
+!  * xd / yd are initial lon / lats.
+!  * xf / yf are output  lon / lats, possibly modified to satisfy configuration.
+!  * interbar is TRUE for barycentric interpolation.
+!-------------------------------------------------------------------------------
+  USE assert_eq_m, ONLY: assert_eq
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN)    :: title
+  REAL,             INTENT(IN)    :: xd(:), yd(:)  ! dim (lons) (lats)
+  REAL,             INTENT(INOUT) :: xf(:), yf(:)  ! dim (lons) (lats)
+  REAL,             INTENT(INOUT) :: champd(:,:)   ! dim (lons,lats)
+  LOGICAL,          INTENT(IN)    :: interbar
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="conf_dat2d"
+  INTEGER :: i, j, ip180, ind, lons, lats
+  LOGICAL :: radlon, invlon ,radlat, invlat
+  REAL    :: pi, pis2, depi, rlatmin, rlatmax, oldxd1
+  REAL, ALLOCATABLE :: xtemp(:) , ytemp(:), champf(:,:)
+!-------------------------------------------------------------------------------
+  lons=assert_eq(SIZE(xd),SIZE(xf),SIZE(champd,1),TRIM(modname)//" lons")
+  lats=assert_eq(SIZE(yd),SIZE(yf),SIZE(champd,2),TRIM(modname)//" lats")
+  ALLOCATE(xtemp(lons),ytemp(lats)); xtemp(:)=xd(:); ytemp(:)=yd(:)
+  ALLOCATE(champf(lons,lats))
+  pi = 2. * ASIN(1.) 
+  pis2 = pi/2.
+  depi = 2. * pi
+  radlon=.FALSE.; invlon=.FALSE.
+  IF     (xtemp(1)>=-pi-0.5.AND.xtemp(lons)<=  pi+0.5) THEN
+    radlon = .TRUE.;  invlon = .FALSE.
+  ELSE IF(xtemp(1)>=   -0.5.AND.xtemp(lons)<=depi+0.5) THEN
+    radlon = .TRUE.;  invlon = .TRUE.
+  ELSE IF(xtemp(1)>= -180.5.AND.xtemp(lons)<=   180.5) THEN
+    radlon = .FALSE.; invlon = .FALSE.
+  ELSE IF(xtemp(1)>=   -0.5.AND.xtemp(lons)<=   360.5) THEN
+    radlon = .FALSE.; invlon = .TRUE.
+  ELSE; WRITE(6,*) 'Problems with data longitudes for file '//TRIM(title)
+  END IF
+  invlat = ytemp(1)<ytemp(lats)
+  rlatmin = MIN( ytemp(1), ytemp(lats) )
+  rlatmax = MAX( ytemp(1), ytemp(lats) )
+      
+  IF     (rlatmin>=-pis2-0.5.AND.rlatmax<=pis2+0.5) THEN; radlat = .TRUE.
+  ELSE IF(rlatmin>= -90.-0.5.AND.rlatmax<= 90.+0.5) THEN; radlat = .FALSE.
+  ELSE; WRITE(6,*) 'Problems with data latitudes for file '//TRIM(title)
+  END IF
+
+  IF(.NOT.radlon) xtemp(:)=xtemp(:)*pi/180.
+  IF(.NOT.radlat) ytemp(:)=ytemp(:)*pi/180.
+
+!--- FLIPPED LONGITUDES
+  IF(invlon) THEN
+    champf(:,:)=champd(:,:); xf(:)=xtemp(:)
+
+  !--- Longitudes rotated to get them in [ -pi pi] interval.
+    DO i=1,lons;      IF(xf(i)>pi) EXIT;                   END DO; ip180 = i
+    DO i=1,lons;      IF(xf(i)>pi) xf(i)=xf(i)-depi;       END DO
+    DO i= ip180,lons; ind=i-ip180+1; xtemp(ind)=xf(i);     END DO
+    DO i= ind+1,lons;                xtemp(i  )=xf(i-ind); END DO
+
+  !--- Longitudes rotated in champf
+    DO j=1,lats
+      DO i=ip180,lons; ind=i-ip180+1; champd(ind,j)=champf(i    ,j); END DO
+      DO i=ind+1,lons;                champd(i  ,j)=champf(i-ind,j); END DO
+    END DO
+  END IF
+
+!--- FLIPPED LATITUDES
+  IF(invlat) THEN
+    yf(:)=ytemp(:)
+    champf(:,:)=champd(:,:)
+    ytemp(lats:1:-1)=yf(:)
+    DO j=1,lats; champd(:,lats-j+1)=champf(:,j); END DO
+  END IF
+
+!--- FOR BARYCENTRIC INTERPOLATION
+  IF(interbar) THEN
+    oldxd1 = xtemp(1)
+    xtemp(1:lons-1)=0.5*(xtemp(1:lons-1)+xtemp(2:lons))
+    xtemp(  lons  )=0.5*(xtemp(  lons  )+oldxd1+depi)
+    ytemp(1:lats-1)=0.5*(ytemp(1:lats-1)+ytemp(2:lats))
+  END IF
+  DEALLOCATE(champf)
+  xf(:)=xtemp(:); DEALLOCATE(xtemp)
+  yf(:)=ytemp(:); DEALLOCATE(ytemp)
+
+END SUBROUTINE conf_dat2d
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE conf_dat3d(title, xd, yd, zd, xf, yf, zf, champd, interbar)
+!
+!-------------------------------------------------------------------------------
+! Author: P. Le Van
+!-------------------------------------------------------------------------------
+! Purpose: Configure the 3D data field "champd" so that:
+!     - Longitudes are in [ -pi     pi   ]
+!     - Latitudes  are in [  pi/2. -pi/2 ]
+!     - Vertical levels from ground to model top (in Pascals)
+!  * xd / yd are initial lon / lats.
+!  * xf / yf are output  lon / lats, possibly modified to satisfy configuration.
+!  * zf      are output pressures
+!  * interbar is TRUE for barycentric interpolation.
+!-------------------------------------------------------------------------------
+  USE assert_eq_m, ONLY: assert_eq
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN)    :: title
+  REAL,             INTENT(IN)    :: xd(:), yd(:), zd(:) ! (lons) (lats) (levs)
+  REAL,             INTENT(INOUT) :: xf(:), yf(:), zf(:) ! (lons) (lats) (levs)
+  REAL,             INTENT(INOUT) :: champd(:,:,:)       ! (lons,lats,levs)
+  LOGICAL,          INTENT(IN)    :: interbar
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="conf_dat3d"
+  INTEGER :: i, j, l, ip180, ind, lons, lats, levs
+  LOGICAL :: radlon, invlon ,radlat, invlat, invlev
+  REAL    :: pi, pis2, depi, presmax, rlatmin, rlatmax, oldxd1
+  REAL, ALLOCATABLE :: xtemp(:) , ytemp(:), ztemp(:), champf(:,:,:)
+!-------------------------------------------------------------------------------
+  lons=assert_eq(SIZE(xd),SIZE(xf),SIZE(champd,1),TRIM(modname)//" lons")
+  lats=assert_eq(SIZE(yd),SIZE(yf),SIZE(champd,2),TRIM(modname)//" lats")
+  levs=assert_eq(SIZE(zd),SIZE(zf),SIZE(champd,3),TRIM(modname)//" levs")
+  ALLOCATE(xtemp(lons),ytemp(lats),ztemp(levs),champf(lons,lats,levs))
+  xtemp(:)=xd(:); ytemp(:)=yd(:); ztemp(:)=zd(:)
+  pi = 2. * ASIN(1.) 
+  pis2 = pi/2.
+  depi = 2. * pi
+  radlon=.FALSE.; invlon=.FALSE.
+  IF     (xtemp(1)>=-pi-0.5.AND.xtemp(lons)<=  pi+0.5) THEN
+    radlon = .TRUE.;  invlon = .FALSE.
+  ELSE IF(xtemp(1)>=   -0.5.AND.xtemp(lons)<=depi+0.5) THEN
+    radlon = .TRUE.;  invlon = .TRUE.
+  ELSE IF(xtemp(1)>= -180.5.AND.xtemp(lons)<=   180.5) THEN
+    radlon = .FALSE.; invlon = .FALSE.
+  ELSE IF(xtemp(1)>=   -0.5.AND.xtemp(lons)<=   360.5) THEN
+    radlon = .FALSE.; invlon = .TRUE.
+  ELSE; WRITE(6,*) 'Problems with data longitudes for file '//TRIM(title)
+  END IF
+  invlat = ytemp(1)<ytemp(lats)
+  rlatmin = MIN( ytemp(1), ytemp(lats) )
+  rlatmax = MAX( ytemp(1), ytemp(lats) )
+      
+  IF     (rlatmin>=-pis2-0.5.AND.rlatmax<=pis2+0.5) THEN; radlat = .TRUE.
+  ELSE IF(rlatmin>= -90.-0.5.AND.rlatmax<= 90.+0.5) THEN; radlat = .FALSE.
+  ELSE; WRITE(6,*) 'Problems with data latitudes for file '//TRIM(title)
+  END IF
+
+  IF(.NOT.radlon) xtemp(:)=xtemp(:)*pi/180.
+  IF(.NOT.radlat) ytemp(:)=ytemp(:)*pi/180.
+
+!--- FLIPPED LONGITUDES
+  IF(invlon) THEN
+    champf(:,:,:)=champd(:,:,:); xf(:)=xtemp(:)
+
+  !--- Longitudes rotated to get them in [ -pi pi] interval.
+    DO i=1,lons;      IF(xf(i)>pi) EXIT;                   END DO; ip180 = i
+    DO i=1,lons;      IF(xf(i)>pi) xf(i)=xf(i)-depi;       END DO
+    DO i= ip180,lons; ind=i-ip180+1; xtemp(ind)=xf(i);     END DO
+    DO i= ind+1,lons;                xtemp(i  )=xf(i-ind); END DO
+
+  !--- Longitudes rotated in champf
+    DO l=1,levs
+    DO j=1,lats
+      DO i=ip180,lons; ind=i-ip180+1; champd(ind,j,l)=champf(i    ,j,l); END DO
+      DO i=ind+1,lons;                champd(i  ,j,l)=champf(i-ind,j,l); END DO
+    END DO
+    END DO
+  END IF
+
+!--- FLIPPED LATITUDES
+  IF(invlat) THEN
+    yf(:)=ytemp(:)
+    champf(:,:,:)=champd(:,:,:)
+    ytemp(lats:1:-1)=yf(:)
+    DO l=1,levs
+      DO j=1,lats; champd(:,lats-j+1,l)=champf(:,j,l); END DO
+    END DO
+  END IF
+
+!--- FOR BARYCENTRIC INTERPOLATION
+  IF(interbar) THEN
+    oldxd1 = xtemp(1)
+    xtemp(1:lons-1)=0.5*(xtemp(1:lons-1)+xtemp(2:lons))
+    xtemp(  lons  )=0.5*(xtemp(  lons  )+oldxd1+depi)
+    ytemp(1:lats-1)=0.5*(ytemp(1:lats-1)+ytemp(2:lats))
+  END IF
+
+!--- FLIPPED LEVELS
+  invlev=ztemp(1)<ztemp(levs)
+  IF(MAX(ztemp(1),ztemp(levs))<1200.) ztemp(:)=ztemp(:)*100.
+  IF(invlev) THEN
+    zf(:)=ztemp(:)
+    champf(:,:,:)=champd(:,:,:)
+    ztemp(levs:1:-1)=zf(:)
+    DO l=1,levs; champd(:,:,levs+1-l)=champf(:,:,l); END DO
+  END IF
+
+  DEALLOCATE(champf)
+  xf(:)=xtemp(:); DEALLOCATE(xtemp)
+  yf(:)=ytemp(:); DEALLOCATE(ytemp)
+  zf(:)=ztemp(:); DEALLOCATE(ztemp)
+
+END SUBROUTINE conf_dat3d
+!
+!-------------------------------------------------------------------------------
+
+END MODULE conf_dat_m
+
Index: dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/grid_atob_m.f90
===================================================================
--- dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/grid_atob_m.f90	(revision 3883)
+++ dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/grid_atob_m.f90	(revision 3883)
@@ -0,0 +1,292 @@
+!*******************************************************************************
+!
+MODULE grid_atob_m
+!
+!*******************************************************************************
+
+  USE assert_eq_m, ONLY: assert_eq
+
+  PRIVATE
+  PUBLIC :: grille_m, rugosite, sea_ice, rugsoro
+
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE fine2coarse(x_i, y_i, x_o, y_o, d_o1, d_i, msk, d_o2)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL,              INTENT(IN)  :: x_i(:), y_i(:) !-- INPUT  X&Y COOR. (mi)(ni)
+  REAL,              INTENT(IN)  :: x_o(:), y_o(:) !-- OUTPUT X&Y COOR. (mo)(no)
+  REAL,              INTENT(OUT) :: d_o1(:,:)      !-- OUTPUT FIELD     (mo,no)
+  REAL,    OPTIONAL, INTENT(IN)  :: d_i (:,:)      !-- INPUT FIELD      (mi,ni)
+  LOGICAL, OPTIONAL, INTENT(IN)  :: msk (:,:)      !-- MASK             (mo,no)
+  REAL,    OPTIONAL, INTENT(OUT) :: d_o2(:,:)      !-- OUTPUT FOR d_i^2 (mo,no)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="fine2coarse"
+  INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1, n1, m2, n2, mx, my, nn
+  LOGICAL :: li, lo
+  REAL    :: inc
+  INTEGER, DIMENSION(SIZE(x_o),SIZE(y_o)) :: num_tot
+  LOGICAL, DIMENSION(SIZE(x_o),SIZE(y_o)) :: found, mask
+  REAL,    DIMENSION(SIZE(x_o),SIZE(y_o)) :: dist
+  REAL,    DIMENSION(SIZE(x_o))           :: a, b
+  REAL,    DIMENSION(SIZE(y_o))           :: c, d
+  REAL,    PARAMETER :: thresh=1.E-5
+!-------------------------------------------------------------------------------
+  mi=SIZE(x_i); ni=SIZE(y_i); mo=SIZE(x_o); no=SIZE(y_o)
+  m1=m1; m2=mo; mx=mo; IF(PRESENT(msk)) mx=SIZE(msk,1)
+  n1=ni; n2=no; my=no; IF(PRESENT(msk)) my=SIZE(msk,2)
+  li=PRESENT(d_i ); IF(li) THEN; m1=SIZE(d_i ,1); n1=SIZE(d_i ,2); END IF
+  lo=PRESENT(d_o2); IF(lo) THEN; m2=SIZE(d_o2,1); n2=SIZE(d_o2,2); END IF
+  mi=assert_eq(mi,m1,TRIM(modname)//" mi")
+  ni=assert_eq(ni,n1,TRIM(modname)//" ni")
+  mo=assert_eq(mo,m2,mx,SIZE(d_o1,1),TRIM(modname)//" mo")
+  no=assert_eq(no,n2,my,SIZE(d_o1,2),TRIM(modname)//" no")
+  mask(:,:)=.TRUE.; IF(PRESENT(msk)) mask(:,:)=msk(:,:)
+
+!--- COMPUTE CELLS INTERFACES COORDINATES OF OUTPUT GRID
+  b(mo)=x_o(mo)+(x_o(mo)-x_o(mo-1))/2.0; b(1:mo-1)=(x_o(1:mo-1)+x_o(2:mo))/2.0
+  d(no)=y_o(no)+(y_o(no)-y_o(no-1))/2.0; d(1:no-1)=(y_o(1:no-1)+y_o(2:no))/2.0
+  a(1 )=x_o(1 )-(x_o(2 )-x_o(1   ))/2.0; a(2:mo  )=   b(1:mo-1)
+  c(1 )=y_o(1 )-(y_o(2 )-y_o(1   ))/2.0; c(2:no  )=   d(1:no-1)
+
+!--- ACCUMULATE INPUT POINTS ON OUTPUT GRID
+  d_o1(:,:)=0.; num_tot(:,:)=0; inc=1.0
+  IF(lo) d_o2(:,:)=0.
+  DO ji = 1, ni
+    DO ii = 1, mi
+      IF(li) inc=d_i(ii,ji)
+      DO jo = 1, no
+        IF((y_i(ji)-c(jo)<thresh.OR.y_i(ji)-d(jo)>thresh).AND.   &
+           (y_i(ji)-c(jo)>thresh.OR.y_i(ji)-d(jo)<thresh)) CYCLE
+        DO io = 1, mo
+          IF((x_i(ii)-a(io)<thresh.OR.x_i(ii)-b(io)>thresh).AND. &
+             (x_i(ii)-a(io)>thresh.OR.x_i(ii)-b(io)<thresh)) CYCLE
+          num_tot(io,jo)=num_tot(io,jo)+1
+          IF(mask(io,jo)) d_o1(io,jo)=d_o1(io,jo)+inc
+          IF(.NOT.lo) CYCLE
+          IF(mask(io,jo)) d_o2(io,jo)=d_o2(io,jo)+inc*inc
+        END DO
+      END DO
+    END DO
+  END DO
+
+!--- CHECK INPUT POINTS HAVE BEEN FOUND IN EACH OUTPUT CELL
+  found(:,:)=num_tot(:,:)/=0
+  WHERE(found.AND.mask) d_o1(:,:)=d_o1(:,:)/REAL(num_tot(:,:))
+  IF(PRESENT(d_o2)) THEN
+    WHERE(found.AND.mask) d_o2(:,:)=d_o2(:,:)/REAL(num_tot(:,:))
+    RETURN
+  END IF
+  nn=COUNT(found); IF(nn==0) RETURN
+
+!--- MISSING POINTS ; USE DISTANCE ON THE SPHERE TO FIND NEAREST POINT nr(2)
+  DO io = 1, mo
+    DO jo = 1, no
+      IF(found(io,jo)) CYCLE
+!      IF(prt_level>=1) PRINT*, "Problem: point out of domain (i,j)=", io,jo
+      CALL dist_sphe(x_o(io),y_o(jo),x_i,y_i,dist(:,:))
+      nr=MINLOC(dist(:,:))!; IF(prt_level>=1) PRINT*, "Solution: ", nr
+      inc=1.0; IF(li) inc=d_i(nr(1),nr(2))
+      IF(mask(io,jo)) d_o1(io,jo)=inc
+    END DO
+  END DO
+
+END SUBROUTINE fine2coarse
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE grille_m(xdata, ydata, entree, x, y, sortie)
+!
+!-------------------------------------------------------------------------------
+! Author:  Z.X. Li (april 1st 1994) (see also A. Harzallah and L. Fairhead)
+!-------------------------------------------------------------------------------
+! Purpose: Naive method to regrid on a coarser grid.
+!   Value at a new point is an average of the old points lcoated in a zone
+!   surrounding the considered new point.
+!   No ponderation is used (see grille_p)
+!
+!           (c)
+!        ----d-----
+!        | . . . .|
+!        |        |
+!     (b)a . * . .b(a)
+!        |        |
+!        | . . . .|
+!        ----c-----
+!           (d)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)  :: xdata(:), ydata(:)       !--- INPUT  FIELD X AND Y COORD.
+  REAL, INTENT(IN)  :: entree(SIZE(xdata),SIZE(ydata)) !--- INPUT  FIELD
+  REAL, INTENT(IN)  :: x(:), y(:)               !--- OUTPUT FIELD X AND Y COORD.
+  REAL, INTENT(OUT) :: sortie(SIZE(x),SIZE(y))  !--- OUTPUT FIELD 
+!-------------------------------------------------------------------------------
+  CALL fine2coarse(xdata,ydata,x,y,sortie,entree)
+
+END SUBROUTINE grille_m
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE rugosite(xdata, ydata, entree, x, y, sortie, mask)
+!
+!-------------------------------------------------------------------------------
+! Author:  Z.X. Li (april 1st 1994)
+!-------------------------------------------------------------------------------
+! Purpose: Remap rugosity length ; constant value (0.001) on oceans.
+! Naive method  (see grille_m)
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)  :: xdata(:), ydata(:)      !--- INPUT  FIELD X AND Y COORD.
+  REAL, INTENT(IN)  :: entree(SIZE(xdata),SIZE(ydata)) !--- INPUT  FIELD
+  REAL, INTENT(IN)  :: x(:), y(:)              !--- OUTPUT FIELD X AND Y COORD.
+  REAL, INTENT(OUT) :: sortie(SIZE(x),SIZE(y)) !--- OUTPUT FIELD 
+  REAL, INTENT(IN)  :: mask  (SIZE(x),SIZE(y)) !--- MASK
+!-------------------------------------------------------------------------------
+  CALL fine2coarse(xdata,ydata,x,y,sortie,LOG(entree))
+  WHERE(NINT(mask)==1)
+    sortie(:,:)=EXP(sortie(:,:))
+  ELSE WHERE
+    sortie(:,:)=0.001
+  END WHERE
+
+END SUBROUTINE rugosite
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE sea_ice(xdata, ydata, glace01, x, y, frac_ice)
+!
+!-------------------------------------------------------------------------------
+! Author:  Z.X. Li (april 1st 1994)
+! Purpose: Regrid ice indicator (0 or 1) on a coarser grid to get an ice fract.
+! field (between 0 and 1).
+! Naive method  (see grille_m)
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)  :: xdata(:), ydata(:)      !--- INPUT  FIELD X AND Y COORD.
+  REAL, INTENT(IN)  :: glace01(:,:)            !--- INPUT  FIELD
+  REAL, INTENT(IN)  :: x(:), y(:)              !--- OUTPUT FIELD X AND Y COORD.
+  REAL, INTENT(OUT) :: frac_ice(SIZE(x),SIZE(y)) !--- OUTPUT FIELD 
+!-------------------------------------------------------------------------------
+  CALL fine2coarse(xdata,ydata,x,y,frac_ice,msk=NINT(glace01)==1)
+
+END SUBROUTINE sea_ice
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE rugsoro(xrel, yrel, relief, xmod, ymod, rugs)
+!
+!-------------------------------------------------------------------------------
+! Author:  Z.X. Li (april 1st 1994) ; Modifications: august 23rd 1995.
+! Purpose: Compute rugosity due to orography by using standard dev in a 1x1 cell
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)  :: xrel(:), yrel(:)        !--- INPUT  FIELD X AND Y COORD.
+  REAL, INTENT(IN)  :: relief(:,:)             !--- INPUT  FIELD
+  REAL, INTENT(IN)  :: xmod(:), ymod(:)        !--- OUTPUT FIELD X AND Y COORD.
+  REAL, INTENT(OUT) :: rugs(SIZE(xmod),SIZE(ymod)) !--- OUTPUT FIELD 
+!-------------------------------------------------------------------------------
+! Local variable:
+  INTEGER           :: k, nn
+  INTEGER, PARAMETER:: itmp=360, jtmp=180
+  REAL  :: out(SIZE(xmod),SIZE(xmod)), amin, amax
+  REAL  :: cham1tmp(itmp,jtmp), cham2tmp(itmp,jtmp), xtmp(itmp), ytmp(jtmp)
+!-------------------------------------------------------------------------------
+
+!--- TEST INPUT FILE FITS FOR ONE DEGREE RESOLUTION
+  xtmp(:)=4.0*ATAN(1.0)*[(-1.0+REAL(2*k-1)/REAL(itmp),k=1,itmp)]
+  ytmp(:)=2.0*ATAN(1.0)*[(-1.0+REAL(2*k-1)/REAL(jtmp),k=1,jtmp)]
+  CALL fine2coarse(xrel,yrel,xtmp,ytmp,cham1tmp,relief,d_o2=cham2tmp)
+  cham2tmp(:,:)=cham2tmp(:,:)-cham1tmp(:,:)**2
+  nn=COUNT(cham2tmp<=-7.5)
+  IF(nn/=0) THEN
+    PRINT*,'Problem for rugsoro ; std**2 < -7.5 for several points: ',nn
+    CALL ABORT_GCM("", "", 1)
+  END IF
+  nn=COUNT(cham2tmp<0)
+  IF(nn/=0) PRINT*,'Problem for rugsoro ; std**2 < 0. for several points: ',nn
+  WHERE(cham2tmp<0.0) cham2tmp=0.0
+  cham2tmp(:,:)=SQRT(cham2tmp(:,:))
+  amin=MINVAL(cham2tmp); amax=MAXVAL(cham2tmp)
+  PRINT*, 'Ecart-type 1x1:', amin, amax
+
+!--- COMPUTE RUGOSITY AT REQUIRED SCALE
+  WHERE(cham2tmp<0.001) cham2tmp=0.001
+  CALL fine2coarse(xtmp,ytmp,xmod,ymod,out,REAL(LOG(cham2tmp)))
+  out=EXP(out)
+  amin=MINVAL(out); amax=MAXVAL(out)
+  PRINT*, 'Ecart-type du modele:', amin, amax
+  out=out/amax*20.0
+  amin=MINVAL(out); amax=MAXVAL(out)
+  PRINT*, 'Longueur de rugosite du modele:', amin, amax
+  rugs=REAL(out)
+
+END SUBROUTINE rugsoro
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,distance)
+!
+!-------------------------------------------------------------------------------
+! Author:  Laurent Li (december 30th 1996).
+! Purpose: Compute min. distance (along big circle) between 2 points in radians.
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN) :: rf_lon, rf_lat  !--- Reference point coordinates (radians)
+  REAL, INTENT(IN) :: rlon(:), rlat(:)!--- Points longitudes/latitudes (radians)
+  REAL, INTENT(OUT):: distance(SIZE(rlon),SIZE(rlat)) !--- Distance    (radians)
+!-------------------------------------------------------------------------------
+! Local variables:
+  LOGICAL, SAVE :: first=.TRUE.
+  REAL,    SAVE :: pi, hpi
+  REAL    :: pa, pb, cpa, spa, crlo(SIZE(rlon))
+  INTEGER :: i, j
+!-------------------------------------------------------------------------------
+  IF(first) THEN; pi=4.0*ATAN(1.0); hpi=pi/2.0; first=.FALSE.; END IF
+  crlo(:)=COS(rf_lon-rlon(:))         !--- COS of points 1 and 2 angle
+  pa=hpi-rf_lat                       !--- North Pole - Point 1 distance
+  cpa=COS(pa); spa=SIN(pa)
+  DO j=1,SIZE(rlat)
+    pb=hpi-rlat(j)                    !--- North Pole - Point 2 distance
+    distance(:,j)=ACOS(cpa*COS(pb)+spa*SIN(pb)*crlo(:))
+  END DO
+
+END SUBROUTINE dist_sphe
+!
+!-------------------------------------------------------------------------------
+
+END MODULE grid_atob_m
+!
+!*******************************************************************************
+
Index: dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/ce0l.F90
===================================================================
--- dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/ce0l.F90	(revision 3882)
+++ dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/ce0l.F90	(revision 3883)
@@ -1,109 +1,214 @@
+PROGRAM ce0l
 !
-! $Id: ce0l.F90 2248 2015-03-25 18:04:54Z lguez $
-!
-!-------------------------------------------------------------------------------
-!
-PROGRAM ce0l
-!-------------------------------------------------------------------------------
-! Purpose: Calls etat0, creates initial states and limit_netcdf
-!
-!     interbar=.T. for barycentric interpolation inter_barxy
-!     extrap  =.T. for data extrapolation, like for the SSTs when file does not
-!                  contain ocean points only.
-!     oldice  =.T. for old-style ice, obtained using grille_m (grid_atob).
-!     masque is created in etat0, passed to limit to ensure consistancy.
-!-------------------------------------------------------------------------------
-  USE control_mod
-#ifdef CPP_EARTH
-! This prog. is designed to work for Earth
-  USE dimphy
-  USE comgeomphy
-  USE infotrac
-  USE indice_sol_mod
-
+!-------------------------------------------------------------------------------
+! Purpose: Initial states and boundary conditions files creation:
+!     * start.nc    for dynamics    (using etat0dyn     routine)
+!     * startphy.nc for physics     (using etat0phys    routine)
+!     * limit.nc    for forced runs (using limit_netcdf routine)
+!-------------------------------------------------------------------------------
+! Notes:
+!     * extrap=.T. (default) for data extrapolation, like for the SSTs when file
+!                   does contain ocean points only.
+!     * "masque" can be:
+!       - read from file "o2a.nc"          (for coupled runs).
+!       - created in etat0phys or etat0dyn (for forced  runs).
+!     It is then passed to limit_netcdf to ensure consistancy.
+!-------------------------------------------------------------------------------
+  USE ioipsl, ONLY: ioconf_calendar, getin, flininfo, flinopen, flinget, flinclo
+  USE control_mod,    ONLY: day_step, dayref, nsplit_phys
+  USE etat0dyn,       ONLY: etat0dyn_netcdf
+  USE etat0phys,      ONLY: etat0phys_netcdf
+  USE limit,          ONLY: limit_netcdf
+  USE netcdf,         ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
+  USE infotrac,       ONLY: type_trac, infotrac_init
+  USE dimphy,         ONLY: klon
+  USE test_disvert_m, ONLY: test_disvert
+  USE filtreg_mod,    ONLY: inifilr
+!  USE iniphysiq_mod,  ONLY: iniphysiq
+  USE mod_const_mpi,  ONLY: comm_lmdz
+#ifdef inca
+  USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic
+#endif
+#ifdef CPP_PARA
+  USE mod_const_mpi,  ONLY: init_const_mpi
+  USE parallel_lmdz,  ONLY: init_parallel, mpi_rank, omp_rank, mpi_size
+  USE bands,          ONLY: read_distrib, distrib_phys
+  USE mod_hallo,      ONLY: init_mod_hallo
+  USE mod_interface_dyn_phys, ONLY: init_interface_dyn_phys
+#endif
+
+  IMPLICIT NONE
+
+!-------------------------------------------------------------------------------
+! Local variables:
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom2.h"
+  include "comconst.h"
+  include "comvert.h"
+  include "iniprint.h"
+  include "temps.h"
+  include "logic.h"
+  REAL               :: masque(iip1,jjp1)             !--- CONTINENTAL MASK
+  REAL               :: phis  (iip1,jjp1)             !--- GROUND GEOPOTENTIAL
+  CHARACTER(LEN=256) :: modname, fmt, calnd           !--- CALENDAR TYPE
+  LOGICAL            :: use_filtre_fft
+  LOGICAL, PARAMETER :: extrap=.FALSE.
+
+!--- Local variables for ocean mask reading:
+  INTEGER            :: nid_o2a, iml_omask, jml_omask, j
+  INTEGER            :: fid, iret, llm_tmp, ttm_tmp, itaul(1)
+  REAL, ALLOCATABLE  :: lon_omask(:,:), dlon_omask(:), ocemask(:,:)
+  REAL, ALLOCATABLE  :: lat_omask(:,:), dlat_omask(:), ocetmp (:,:)
+  REAL               :: date, lev(1)
+#ifndef CPP_PARA
+! for iniphysiq in serial mode
+  INTEGER,PARAMETER :: mpi_rank=0
+  INTEGER :: distrib_phys(mpi_rank:mpi_rank)=(jjm-1)*iim+2
+#endif
+!-------------------------------------------------------------------------------
+  modname="ce0l"
+
+!--- Constants
+  pi     = 4. * ATAN(1.)
+  rad    = 6371229.
+  daysec = 86400.
+  omeg   = 2.*pi/daysec
+  g      = 9.8
+  kappa  = 0.2857143
+  cpp    = 1004.70885
+  jmp1   = jjm + 1
+  preff   = 101325.
+  pa      = 50000.
+
+  CALL conf_gcm( 99, .TRUE. )
+  dtvr = daysec/REAL(day_step)
+  WRITE(lunout,*)'dtvr',dtvr
+  CALL iniconst()
+  CALL inigeom()
+
+!--- Calendar choice
 #ifdef CPP_IOIPSL
-  USE ioipsl, ONLY: ioconf_calendar, getin
-#else
-  ! if not using IOIPSL, we still need to use (a local version of) getin
-  use ioipsl_getincom, only: getin
-#endif
-
-#endif
-  IMPLICIT NONE
-#ifndef CPP_EARTH
-#include "iniprint.h"
-  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
-#else
-!-------------------------------------------------------------------------------
-! Local variables:
-  LOGICAL, PARAMETER :: interbar=.TRUE., extrap=.FALSE., oldice=.FALSE.
-#include "dimensions.h"
-#include "paramet.h"
-!#include "indicesol.h"
-#include "iniprint.h"
-#include "temps.h"
-#include "logic.h"
-  REAL,    DIMENSION(iip1,jjp1) :: masque
-  CHARACTER(LEN=15)             :: calnd
-  REAL,    DIMENSION(iip1,jjp1) :: phis ! geopotentiel au sol
-  logical use_filtre_fft
-!-------------------------------------------------------------------------------
-  CALL conf_gcm( 99, .TRUE. )
-
-  use_filtre_fft=.FALSE.
-  CALL getin('use_filtre_fft',use_filtre_fft)
-  IF (use_filtre_fft) THEN
-     write(lunout, fmt = *) 'FFT filter is not available in the ' &
-          // 'sequential version of the dynamics.'
-     write(lunout, fmt = *) &
-          "Your setting of variable use_filtre_fft is not used."
-  ENDIF
-
-  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
-  WRITE(lunout,*)'---> klon=',klon
-  CALL InitComgeomphy
-
-#ifdef CPP_IOIPSL
+  calnd='gregorian'
   SELECT CASE(calend)
-    CASE('earth_360d');CALL ioconf_calendar('360d');      calnd='a 360 jours/an'
-    CASE('earth_365d');CALL ioconf_calendar('noleap');    calnd='a 365 jours/an'
-    CASE('earth_366d');CALL ioconf_calendar('366d');      calnd='bissextile'
-    CASE('gregorian'); CALL ioconf_calendar('gregorian'); calnd='gregorien'
-    CASE('standard');  CALL ioconf_calendar('gregorian'); calnd='gregorien'
-    CASE('julian');    CALL ioconf_calendar('julian');    calnd='julien'
+    CASE('earth_360d');CALL ioconf_calendar('360d');   calnd='with 360 days/year'
+    CASE('earth_365d');CALL ioconf_calendar('noleap'); calnd='with no leap year'
+    CASE('earth_366d');CALL ioconf_calendar('366d');   calnd='with leap years only'
+    CASE('gregorian'); CALL ioconf_calendar('gregorian')
+    CASE('standard');  CALL ioconf_calendar('gregorian')
+    CASE('julian');    CALL ioconf_calendar('julian'); calnd='julian'
     CASE('proleptic_gregorian'); CALL ioconf_calendar('gregorian')
   !--- DC Bof...  => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian
     CASE DEFAULT
-      CALL abort_gcm('ce0l','Mauvais choix de calendrier',1)
+      CALL abort_gcm('ce0l','Bad choice for calendar',1)
   END SELECT
-  WRITE(lunout,*)'CALENDRIER CHOISI: Terrestre '//TRIM(calnd)
-#endif
-
+  WRITE(lunout,*)'CHOSEN CALENDAR: Earth '//TRIM(calnd)
+#endif
+
+#ifdef CPP_PARA
+!--- Physical grid + parallel initializations
+  CALL init_const_mpi()
+  CALL init_parallel()
+  CALL read_distrib()
+  CALL init_mod_hallo()
+#endif
+  WRITE(lunout,*)'---> klon=',klon
+
+!--- Tracers initializations
   IF (type_trac == 'inca') THEN
 #ifdef INCA
-    CALL init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday)
-    CALL init_inca_para(iim,jjm+1,klon,1,klon_mpi_para_nb,0)
+    CALL init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,&
+                         nbsrf,is_oce,is_sic,is_ter,is_lic,calend)
+    CALL init_inca_para(iim,jjp1,llm,klon_glo,mpi_size,distrib_phys,&
+                        COMM_LMDZ)
     WRITE(lunout,*)'nbtr =' , nbtr 
 #endif
   END IF
-
-  WRITE(lunout,'(//)')
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,*) '  ***  etat0_netcdf ***  '
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,'(//)')
-  WRITE(lunout,*) ' interbar = ',interbar
-  CALL etat0_netcdf(interbar,masque,phis,ok_etat0)
+  CALL infotrac_init()
+
+  CALL inifilr()
+  CALL iniphysiq(iim,jjm, &
+                 distrib_phys(mpi_rank),comm_lmdz, llm, &
+                 daysec,dtphys/nsplit_phys, &
+                 rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
+  IF(pressure_exner) CALL test_disvert
+
+#ifdef CPP_PARA
+  IF (mpi_rank==0.AND.omp_rank==0) THEN
+#endif
+  use_filtre_fft=.FALSE.
+  CALL getin('use_filtre_fft',use_filtre_fft)
+  IF(use_filtre_fft) THEN
+     WRITE(lunout,*)"FFT filter not available for sequential dynamics."
+     WRITE(lunout,*)"Your setting of variable use_filtre_fft is not used."
+  ENDIF
+
+!--- LAND MASK. TWO CASES:
+!   1) read from ocean model    file "o2a.nc"    (coupled runs)
+!   2) computed from topography file "Relief.nc" (masque(:,:)=-99999.)
+! Coupled simulations (case 1) use the ocean model mask to compute the
+! weights to ensure ocean fractions are the same for atmosphere and ocean.
+!*******************************************************************************
+  IF(NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)/=NF90_NOERR) THEN
+    WRITE(lunout,*)'BEWARE !! No ocean mask "o2a.nc" file found'
+    WRITE(lunout,*)'Forced run.'
+    masque(:,:)=-99999.
+  ELSE
+    iret=NF90_CLOSE(nid_o2a)
+    WRITE(lunout,*)'BEWARE !! Ocean mask "o2a.nc" file found'
+    WRITE(lunout,*)'Coupled run.'
+    CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
+    IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
+      WRITE(lunout,*)'Mismatching dimensions for ocean mask'
+      WRITE(lunout,*)'iim  = ',iim ,' iml_omask = ',iml_omask
+      WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
+      CALL abort_gcm(modname,'',1)
+    END IF
+    ALLOCATE(ocemask(iim,jjp1),lon_omask(iim,jjp1),dlon_omask(iim ))
+    ALLOCATE(ocetmp (iim,jjp1),lat_omask(iim,jjp1),dlat_omask(jjp1))
+    CALL flinopen("o2a.nc", .FALSE.,iml_omask,jml_omask,llm_tmp,               &
+                  lon_omask,lat_omask,lev,ttm_tmp,itaul,date,dt,fid)
+    CALL flinget(fid, "OceMask",    iim,jjp1,llm_tmp,ttm_tmp,1,1,ocetmp)
+    CALL flinclo(fid)
+    dlon_omask(1:iim ) = lon_omask(1:iim,1)
+    dlat_omask(1:jjp1) = lat_omask(1,1:jjp1)
+    ocemask = ocetmp
+    IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
+      DO j=1,jjp1; ocemask(:,j) = ocetmp(:,jjp1-j+1); END DO
+    END IF
+    DEALLOCATE(ocetmp,lon_omask,lat_omask,dlon_omask,dlat_omask)
+    IF(prt_level>=1) THEN
+      WRITE(fmt,"(i4,'i1)')")iim ; fmt='('//ADJUSTL(fmt)
+      WRITE(lunout,*)'OCEAN MASK :'
+      WRITE(lunout,fmt) NINT(ocemask)
+    END IF
+    masque(1:iim,:)=1.-ocemask(:,:)
+    masque(iip1 ,:)=masque(1,:)
+    DEALLOCATE(ocemask)
+  END IF
+  phis(:,:)=-99999.
+
+  IF(ok_etat0) THEN
+    WRITE(lunout,'(//)')
+    WRITE(lunout,*) '  ************************  '
+    WRITE(lunout,*) '  ***  etat0phy_netcdf ***  '
+    WRITE(lunout,*) '  ************************  '
+    CALL etat0phys_netcdf(masque,phis)
+    WRITE(lunout,'(//)')
+    WRITE(lunout,*) '  ************************  '
+    WRITE(lunout,*) '  ***  etat0dyn_netcdf ***  '
+    WRITE(lunout,*) '  ************************  '
+    CALL etat0dyn_netcdf(masque,phis)
+  END IF
 
   IF(ok_limit) THEN
-  WRITE(lunout,'(//)')
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,*) '  ***  Limit_netcdf ***  '
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,'(//)')
-  CALL limit_netcdf(interbar,extrap,oldice,masque)
-  END IF
-
-  
+    WRITE(lunout,'(//)')
+    WRITE(lunout,*) '  *********************  '
+    WRITE(lunout,*) '  ***  Limit_netcdf ***  '
+    WRITE(lunout,*) '  *********************  '
+    WRITE(lunout,'(//)')
+    CALL limit_netcdf(masque,phis,extrap)
+  END IF
+
   WRITE(lunout,'(//)')
   WRITE(lunout,*) '  ***************************  '
@@ -113,6 +218,7 @@
   CALL grilles_gcm_netcdf_sub(masque,phis)
 
-#endif
-! of #ifndef CPP_EARTH #else
+#ifdef CPP_PARA
+  END IF
+#endif
 
 END PROGRAM ce0l
Index: dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90
===================================================================
--- dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90	(revision 3883)
+++ dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90	(revision 3883)
@@ -0,0 +1,503 @@
+MODULE etat0dyn
+!
+!*******************************************************************************
+! Purpose: Create dynamical initial state using atmospheric fields from a
+!          database of atmospheric to initialize the model.
+!-------------------------------------------------------------------------------
+! Comments:
+!
+!    *  This module is designed to work for Earth (and with ioipsl)
+!
+!    *  etat0dyn_netcdf routine can access to NetCDF data through the following
+!  routine (to be called after restget):
+!    CALL startget_dyn3d(varname, lon_in,  lat_in, pls, workvar,&
+!                          champ, lon_in2, lat_in2)
+!
+!    *  Variables should have the following names in the NetCDF files:
+!            'U'      : East ward wind              (in "ECDYN.nc")
+!            'V'      : Northward wind              (in "ECDYN.nc")
+!            'TEMP'   : Temperature                 (in "ECDYN.nc")
+!            'R'      : Relative humidity           (in "ECDYN.nc")
+!            'RELIEF' : High resolution orography   (in "Relief.nc") 
+!
+!    * The land mask and corresponding weights can be:
+!      1) already known (in particular if etat0dyn has been called before) ;
+!         in this case, ANY(masque(:,:)/=-99999.) = .TRUE.
+!      2) computed using the ocean mask from the ocean model (to ensure ocean
+!         fractions are the same for atmosphere and ocean) for coupled runs.
+!         File name: "o2a.nc"  ;  variable name: "OceMask"
+!      3) computed from topography file "Relief.nc" for forced runs.
+!
+!   *   There is a big mess with the longitude size. Should it be iml or iml+1 ?
+!  I have chosen to use the iml+1 as an argument to this routine and we declare
+!  internaly smaller fields when needed. This needs to be cleared once and for
+!  all in LMDZ. A convention is required.
+!-------------------------------------------------------------------------------
+  USE ioipsl,         ONLY: flininfo, flinopen, flinget, flinclo, histclo
+  USE assert_eq_m,    ONLY: assert_eq
+  IMPLICIT NONE
+
+  PRIVATE
+  PUBLIC :: etat0dyn_netcdf
+
+  include "iniprint.h"
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom2.h"
+  include "comvert.h"
+  include "comconst.h"
+  include "temps.h"
+  include "comdissnew.h"
+  include "serre.h"
+  REAL, SAVE :: deg2rad
+  INTEGER,            SAVE      :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn
+  REAL, ALLOCATABLE,  SAVE      :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:)
+  CHARACTER(LEN=120), PARAMETER :: dynfname='ECDYN.nc'
+
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE etat0dyn_netcdf(masque, phis)
+!
+!-------------------------------------------------------------------------------
+! Purpose: Create dynamical initial states.
+!-------------------------------------------------------------------------------
+! Notes:  1) This routine is designed to work for Earth
+!         2) If masque(:,:)/=-99999., masque and phis are already known.
+!         Otherwise: compute it.
+!-------------------------------------------------------------------------------
+  USE control_mod
+  USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz
+  USE regr_pr_o3_m,   ONLY: regr_pr_o3
+  USE press_coefoz_m, ONLY: press_coefoz
+  USE exner_hyb_m,    ONLY: exner_hyb
+  USE exner_milieu_m, ONLY: exner_milieu
+  USE infotrac,       ONLY: nqtot, tname
+  USE filtreg_mod
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL,    INTENT(INOUT) :: masque(iip1,jjp1)   !--- Land-ocean mask
+  REAL,    INTENT(INOUT) :: phis  (iip1,jjp1)   !--- Ground geopotential
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname, fmt
+  INTEGER            :: i, j, l, ji, itau, iday
+  REAL               :: xpn, xps, time, phystep
+  REAL, DIMENSION(iip1,jjp1)       :: psol
+  REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d
+  REAL, DIMENSION(iip1,jjp1,llm)   :: uvent, t3d, tpot, qsat, qd
+  REAL, DIMENSION(iip1,jjp1,llm)   :: pk, pls, y, masse
+  REAL, DIMENSION(iip1,jjm ,llm)   :: vvent
+  REAL, DIMENSION(ip1jm    ,llm)   :: pbarv
+  REAL, DIMENSION(ip1jmp1  ,llm)   :: pbaru, phi, w
+  REAL, DIMENSION(ip1jmp1)         :: pks
+  REAL, DIMENSION(iim)             :: xppn, xpps
+  REAL, ALLOCATABLE                :: q3d(:,:,:,:)
+!-------------------------------------------------------------------------------
+  modname='etat0dyn_netcdf'
+
+  deg2rad = pi/180.0
+
+! Compute psol AND tsol, knowing phis.
+!*******************************************************************************
+  CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, phis, psol)
+
+! Mid-levels pressure computation
+!*******************************************************************************
+  CALL pression(ip1jmp1, ap, bp, psol, p3d)             !--- Update p3d
+  IF(pressure_exner) THEN                               !--- Update pk, pks
+    CALL exner_hyb   (ip1jmp1,psol,p3d,pks,pk)
+  ELSE
+    CALL exner_milieu(ip1jmp1,psol,p3d,pks,pk)
+  END IF
+  pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa)          !--- Update pls
+
+! Update uvent, vvent, t3d and tpot
+!*******************************************************************************
+  uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0
+  CALL startget_dyn3d('u'   ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv)
+  CALL startget_dyn3d('v'   ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent,      &
+ &                           rlonu,rlatu(:jjm))
+  CALL startget_dyn3d('t'   ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv)
+  tpot(:,:,:)=t3d(:,:,:)
+  CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv)
+
+  WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:))
+  WRITE(lunout,*) 'PLS min,max:',MINVAL(pls(:,:,:)),MAXVAL(pls(:,:,:))
+
+! Humidity at saturation computation
+!*******************************************************************************
+  WRITE(lunout,*) 'avant q_sat'
+  CALL q_sat(llm*jjp1*iip1, t3d, pls, qsat)
+  WRITE(lunout,*) 'apres q_sat'
+  WRITE(lunout,*) 'QSAT min,max:',MINVAL(qsat(:,:,:)),MAXVAL(qsat(:,:,:))
+!  WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
+  qd (:,:,:) = 0.0
+  CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv)
+  ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:)
+  CALL flinclo(fid_dyn)
+
+! Parameterization of ozone chemistry:
+!*******************************************************************************
+! Look for ozone tracer:
+  DO i=1,nqtot; IF(ANY(["O3","o3"]==tname(i))) EXIT; END DO
+  IF(i/=nqtot+1) THEN
+    CALL regr_lat_time_coefoz
+    CALL press_coefoz
+    CALL regr_pr_o3(p3d, q3d(:,:,:,i))
+    q3d(:,:,:,i)=q3d(:,:,:,i)*48./ 29.                  !--- Mole->mass fraction         
+  END IF
+  q3d(iip1,:,:,:)=q3d(1,:,:,:)
+
+! Writing
+!*******************************************************************************
+  CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, tetagrot,   &
+       tetatemp, vert_prof_dissip)
+  WRITE(lunout,*)'sortie inidissip'
+  itau=0
+  itau_dyn=0
+  itau_phy=0
+  iday=dayref+itau/day_step
+  time=FLOAT(itau-(iday-dayref)*day_step)/day_step
+  IF(time>1.) THEN
+   time=time-1
+   iday=iday+1
+  END IF
+  day_ref=dayref
+  annee_ref=anneeref
+  CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
+  WRITE(lunout,*)'sortie geopot'
+  CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis,               &
+                phi,  w, pbaru, pbarv, time+iday-dayref)
+  WRITE(lunout,*)'sortie caldyn0'     
+#ifdef CPP_PARA
+  CALL dynredem0_loc( "start.nc", dayref, phis)
+#else
+  CALL dynredem0( "start.nc", dayref, phis)
+#endif
+  WRITE(lunout,*)'sortie dynredem0'
+#ifdef CPP_PARA
+  CALL dynredem1_loc( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
+#else
+  CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
+#endif
+  WRITE(lunout,*)'sortie dynredem1' 
+  CALL histclo()
+
+END SUBROUTINE etat0dyn_netcdf
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE startget_dyn3d(var, lon_in,  lat_in,  pls,  workvar,&
+                        champ, lon_in2, lat_in2)
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!===============================================================================
+! Purpose: Compute some quantities (u,v,t,q,tpot) using variables U,V,TEMP and R
+!     (3D fields) of file dynfname.
+!-------------------------------------------------------------------------------
+! Note: An input auxilliary field "workvar" has to be specified in two cases:
+!     * for "q":    the saturated humidity.
+!     * for "tpot": the Exner function.
+!===============================================================================
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN)    :: var
+  REAL,             INTENT(IN)    :: lon_in(:)        ! dim (iml)
+  REAL,             INTENT(IN)    :: lat_in(:)        ! dim (jml)
+  REAL,             INTENT(IN)    :: pls    (:, :, :) ! dim (iml, jml, lml)
+  REAL,             INTENT(IN)    :: workvar(:, :, :) ! dim (iml, jml, lml)
+  REAL,             INTENT(INOUT) :: champ  (:, :, :) ! dim (iml, jml, lml)
+  REAL,             INTENT(IN)    :: lon_in2(:)       ! dim (iml)
+  REAL,             INTENT(IN)    :: lat_in2(:)       ! dim (jml2)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=10)  :: vname
+  CHARACTER(LEN=256) :: msg, modname="startget_dyn3d"
+  INTEGER            :: iml, jml, jml2, lml, il
+  REAL               :: xppn, xpps
+!-------------------------------------------------------------------------------
+  iml=assert_eq([SIZE(lon_in),SIZE(pls,1),SIZE(workvar,1),SIZE(champ,1),       &
+     &                                    SIZE(lon_in2)], TRIM(modname)//" iml")
+  jml=assert_eq( SIZE(lat_in),SIZE(pls,2),SIZE(workvar,2),SIZE(champ,2),       &
+     &                                                    TRIM(modname)//" jml")
+  lml=assert_eq(              SIZE(pls,3),SIZE(workvar,3),SIZE(champ,3),       &
+     &                                                    TRIM(modname)//" lml")
+  jml2=SIZE(lat_in2)
+
+!--- CHECK IF THE FIELD IS KNOWN
+   SELECT CASE(var)
+    CASE('u');    vname='U'
+    CASE('v');    vname='V'
+    CASE('t');    vname='TEMP'
+    CASE('q');    vname='R';    msg='humidity as the saturated humidity'
+    CASE('tpot'); msg='potential temperature as the Exner function'
+    CASE DEFAULT; msg='No rule to extract variable '//TRIM(var)
+      CALL abort_gcm(modname,TRIM(msg)//' from any data set',1)
+  END SELECT
+
+!--- CHECK IF SOMETHING IS MISSING
+  IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN
+    msg='Could not compute '//TRIM(msg)//' is missing or constant.'
+    CALL abort_gcm(modname,TRIM(msg),1)
+  END IF
+
+!--- INTERPOLATE 3D FIELD IF NEEDED
+  IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2,      &
+                                                  lat_in2,pls,champ)
+
+!--- COMPUTE THE REQUIRED FILED
+  SELECT CASE(var)
+    CASE('u'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
+      champ(iml,:,:)=champ(1,:,:)                   !--- Eastward wind
+
+    CASE('v'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
+      champ(iml,:,:)=champ(1,:,:)                   !--- Northward wind
+
+    CASE('tpot','q')
+      IF(var=='tpot') THEN; champ=champ*cpp/workvar !--- Potential temperature
+      ELSE;                 champ=champ*.01*workvar !--- Relative humidity
+        WHERE(champ<0.) champ=1.0E-10
+      END IF
+      DO il=1,lml
+        xppn = SUM(aire(:,1  )*champ(:,1  ,il))/apoln
+        xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
+        champ(:,1  ,il) = xppn
+        champ(:,jml,il) = xpps
+      END DO
+  END SELECT
+
+END SUBROUTINE startget_dyn3d
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,zs,psol)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!===============================================================================
+! Purpose:   Compute psol, knowing phis.
+!===============================================================================
+! Arguments:
+  REAL,    INTENT(IN)  :: lon_in (:),  lat_in (:)    ! dim (iml) (jml)
+  REAL,    INTENT(IN)  :: lon_in2(:),  lat_in2(:)    ! dim (iml) (jml2)
+  REAL,    INTENT(IN)  :: zs  (:,:)                  ! dim (iml,jml)
+  REAL,    INTENT(OUT) :: psol(:,:)                  ! dim (iml,jml)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname='start_init_dyn'
+  REAL               :: date, dt
+  INTEGER            :: iml, jml, jml2, itau(1)
+  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), var_ana(:,:)
+  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:)
+  REAL, ALLOCATABLE  :: z(:,:), ps(:,:), ts(:,:)
+!-------------------------------------------------------------------------------
+  iml=assert_eq(SIZE(lon_in),SIZE(zs,1),SIZE(psol,1),SIZE(lon_in2),            &
+      &                                              TRIM(modname)//" iml")
+  jml=assert_eq(SIZE(lat_in),SIZE(zs,2),SIZE(psol,2),TRIM(modname)//" jml")
+  jml2=SIZE(lat_in2)
+
+  WRITE(lunout,*) 'Opening the surface analysis'
+  CALL flininfo(dynfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
+  WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
+
+  ALLOCATE(lon_dyn(iml_dyn,jml_dyn), lat_dyn(iml_dyn,jml_dyn))
+  ALLOCATE(levdyn_ini(llm_dyn))
+  CALL flinopen(dynfname, .FALSE., iml_dyn, jml_dyn, llm_dyn,                  &
+                lon_dyn,lat_dyn,levdyn_ini,ttm_dyn,itau,date,dt,fid_dyn)
+
+!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
+  ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
+  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
+  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
+
+  ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn))
+  CALL get_var_dyn('Z',z)                        !--- SURFACE GEOPOTENTIAL
+  CALL get_var_dyn('SP',ps)                      !--- SURFACE PRESSURE
+  CALL get_var_dyn('ST',ts)                      !--- SURFACE TEMPERATURE
+!  CALL flinclo(fid_dyn)
+  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
+
+!--- PSOL IS COMPUTED IN PASCALS
+  psol(:iml-1,:) = ps(:iml-1,:)*(1.0+(z(:iml-1,:)-zs(:iml-1,:))/287.0          &
+     &            /ts(:iml-1,:))
+  psol(iml,:)=psol(1,:)
+  DEALLOCATE(z,ps,ts)
+  psol(:,1  )=SUM(aire(1:iml-1,1  )*psol(1:iml-1,1  ))/apoln  !--- NORTH POLE
+  psol(:,jml)=SUM(aire(1:iml-1,jml)*psol(1:iml-1,jml))/apols  !--- SOUTH POLE
+
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE get_var_dyn(title,field)
+!
+!-------------------------------------------------------------------------------
+  USE conf_dat_m, ONLY: conf_dat2d
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*),  INTENT(IN)    :: title
+  REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: msg
+  INTEGER :: tllm
+!-------------------------------------------------------------------------------
+  SELECT CASE(title)
+    CASE('Z');     tllm=0;       msg='geopotential'
+    CASE('SP');    tllm=0;       msg='surface pressure'
+    CASE('ST');    tllm=llm_dyn; msg='temperature'
+  END SELECT
+  IF(.NOT.ALLOCATED(field)) THEN
+    ALLOCATE(field(iml,jml))
+    CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana)
+    CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
+    CALL interp_startvar(title, .TRUE., lon_rad,lat_rad, var_ana,              &
+                                        lon_in, lat_in, lon_in2, lat_in2, field)
+  ELSE IF(SIZE(field)/=SIZE(z)) THEN
+    msg='The '//TRIM(msg)//' field we have does not have the right size'
+    CALL abort_gcm(TRIM(modname),msg,1)
+  END IF
+
+END SUBROUTINE get_var_dyn
+!
+!-------------------------------------------------------------------------------
+
+END SUBROUTINE start_init_dyn
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d)
+!
+!-------------------------------------------------------------------------------
+  USE conf_dat_m, ONLY: conf_dat3d
+  USE pchsp_95_m, ONLY: pchsp_95
+  USE pchfe_95_m, ONLY: pchfe_95
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN) :: var
+  REAL,    INTENT(IN)  :: lon_in(:),  lat_in(:)   ! dim (iml) (jml)
+  REAL,    INTENT(IN)  :: lon_in2(:), lat_in2(:)  ! dim (iml) (jml2)
+  REAL,    INTENT(IN)  :: pls_in(:,:,:)           ! dim (iml,jml,lml)
+  REAL,    INTENT(OUT) :: var3d (:,:,:)           ! dim (iml,jml,lml)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname='start_inter_3d'
+  LOGICAL :: skip
+  REAL    :: chmin, chmax
+  INTEGER :: iml, jml, lml, jml2, ii, ij, il, ierr
+  INTEGER :: n_extrap                             ! Extrapolated points number
+  REAL, ALLOCATABLE :: ax(:), lon_rad(:), lon_ini(:), lev_dyn(:), yder(:)
+  REAL, ALLOCATABLE :: ay(:), lat_rad(:), lat_ini(:), var_tmp3d(:,:,:)
+  REAL, ALLOCATABLE, SAVE :: var_ana3d(:,:,:)
+!-------------------------------------------------------------------------------
+  iml=assert_eq(SIZE(lon_in),SIZE(lon_in2),SIZE(pls_in,1),SIZE(var3d,1),TRIM(modname)//" iml")
+  jml=assert_eq(SIZE(lat_in),              SIZE(pls_in,2),SIZE(var3d,2),TRIM(modname)//" jml")
+  lml=assert_eq(SIZE(pls_in,3),SIZE(var3d,3),TRIM(modname)//" lml"); jml2=SIZE(lat_in2)
+
+  WRITE(lunout, *)'Going into flinget to extract the 3D field.'
+  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
+  CALL flinget(fid_dyn,var,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
+
+!--- ANGLES IN DEGREES ARE CONVERTED INTO RADIANS
+  ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
+  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
+  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
+
+!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
+  ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn))
+  CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini,                           &
+                       lon_rad, lat_rad, lev_dyn, var_ana3d, .TRUE.)
+  DEALLOCATE(lon_ini, lat_ini)
+
+!--- COMPUTE THE REQUIRED FIELDS USING ROUTINE grid_noro
+  ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
+  DO il = 1,llm_dyn
+    CALL interp_startvar(var,il==1,lon_rad,lat_rad,var_ana3d(:,:,il),          &
+                          lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il))
+  END DO
+  DEALLOCATE(lon_rad, lat_rad)
+
+!--- VERTICAL INTERPOLATION FROM TOP OF ATMOSPHERE TO GROUND
+  ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
+  ax = lev_dyn(llm_dyn:1:-1) 
+  skip = .FALSE.
+  n_extrap = 0
+  DO ij=1, jml
+    DO ii=1, iml-1
+      ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
+      yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
+      CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1),              &
+           var3d(ii, ij, lml:1:-1), ierr)
+      IF(ierr<0) CALL abort_gcm(TRIM(modname),'error in pchfe_95',1)
+      n_extrap = n_extrap + ierr
+    END DO
+  END DO
+  IF(n_extrap/=0) WRITE(lunout,*)TRIM(modname)//" pchfe_95: n_extrap=", n_extrap
+  var3d(iml, :, :) = var3d(1, :, :) 
+
+  DO il=1, lml
+    CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
+    WRITE(lunout, *)' '//TRIM(var)//'  min max l ', il, chmin, chmax
+  END DO
+
+END SUBROUTINE start_inter_3d
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
+!
+!-------------------------------------------------------------------------------
+  USE inter_barxy_m, ONLY: inter_barxy
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN)  :: nam
+  LOGICAL,          INTENT(IN)  :: ibeg
+  REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
+  REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
+  REAL,             INTENT(IN)  :: lon1(:), lat1(:) ! dim (i1) (j1)
+  REAL,             INTENT(IN)  :: lon2(:), lat2(:) ! dim (i1) (j2)
+  REAL,             INTENT(OUT) :: varo(:,:)        ! dim (i1) (j1)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="interp_startvar"
+  INTEGER            :: ii, jj, i1, j1, j2
+  REAL, ALLOCATABLE  :: vtmp(:,:)
+!-------------------------------------------------------------------------------
+  ii=assert_eq(SIZE(lon),            SIZE(vari,1),TRIM(modname)//" ii")
+  jj=assert_eq(SIZE(lat),            SIZE(vari,2),TRIM(modname)//" jj")
+  i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
+  j1=assert_eq(SIZE(lat1),           SIZE(varo,2),TRIM(modname)//" j1")
+  j2=SIZE(lat2)
+  ALLOCATE(vtmp(i1-1,j1))
+  IF(ibeg.AND.prt_level>1) THEN
+    WRITE(lunout,*)"---------------------------------------------------------"
+    WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
+    WRITE(lunout,*)"---------------------------------------------------------"
+  END IF
+  CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
+  CALL gr_int_dyn(vtmp, varo, i1-1, j1)
+
+END SUBROUTINE interp_startvar
+!
+!-------------------------------------------------------------------------------
+
+END MODULE etat0dyn
+!
+!*******************************************************************************
Index: dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90
===================================================================
--- dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90	(revision 3883)
+++ dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90	(revision 3883)
@@ -0,0 +1,454 @@
+MODULE etat0phys
+!
+!*******************************************************************************
+! Purpose: Create physical initial state using atmospheric fields from a
+!          database of atmospheric to initialize the model.
+!-------------------------------------------------------------------------------
+! Comments:
+!
+!    *  This module is designed to work for Earth (and with ioipsl)
+!
+!    *  etat0phys_netcdf routine can access to NetCDF data through subroutines:
+!         "start_init_phys" for variables contained in file "ECPHY.nc":
+!            'ST'     : Surface temperature
+!            'CDSW'   : Soil moisture
+!         "start_init_orog" for variables contained in file "Relief.nc":
+!            'RELIEF' : High resolution orography 
+!
+!    * The land mask and corresponding weights can be:
+!      1) computed using the ocean mask from the ocean model (to ensure ocean
+!         fractions are the same for atmosphere and ocean) for coupled runs.
+!         File name: "o2a.nc"  ;  variable name: "OceMask"
+!      2) computed from topography file "Relief.nc" for forced runs.
+!
+!    * Allowed values for read_climoz flag are 0, 1 and 2:
+!      0: do not read an ozone climatology
+!      1: read a single ozone climatology that will be used day and night
+!      2: read two ozone climatologies, the average day and night climatology
+!         and the daylight climatology
+!-------------------------------------------------------------------------------
+!    * There is a big mess with the longitude size. Should it be iml or iml+1 ?
+!  I have chosen to use the iml+1 as an argument to this routine and we declare
+!  internaly smaller fields when needed. This needs to be cleared once and for
+!  all in LMDZ. A convention is required.
+!-------------------------------------------------------------------------------
+
+  USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo
+  USE assert_eq_m,        ONLY: assert_eq
+  USE dimphy,             ONLY: klon
+  USE conf_dat_m,         ONLY: conf_dat2d
+  USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, &
+          solsw, radsol, t_ancien, wake_deltat, wake_s,  rain_fall, qsol, z0h, &
+          sollw, rugoro, q_ancien, wake_deltaq, wake_pe, snow_fall, ratqs,w01, &
+    sig1, ftsol, clwcon, fm_therm, wake_Cstar,  pctsrf,  entr_therm,radpas, f0,&
+    zmax0,fevap, rnebcon,falb_dir, wake_fip,    agesno,  detr_therm, pbl_tke,  &
+    phys_state_var_init
+
+  PRIVATE
+  PUBLIC :: etat0phys_netcdf
+
+  include "iniprint.h"
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom2.h"
+  include "comconst.h"
+  include "dimsoil.h"
+  include "temps.h"
+  include "clesphys.h"
+  REAL, SAVE :: deg2rad
+  REAL, SAVE, ALLOCATABLE :: tsol(:)
+  INTEGER,            SAVE      :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys
+  REAL, ALLOCATABLE,  SAVE      :: lon_phys(:,:), lat_phys(:,:), levphys_ini(:)
+  CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc", orogvar="RELIEF"
+  CHARACTER(LEN=256), PARAMETER :: phyfname="ECPHY.nc",  psrfvar="SP"
+  CHARACTER(LEN=256), PARAMETER :: qsolvar="CDSW",       tsrfvar="ST"
+
+
+CONTAINS
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE etat0phys_netcdf(masque, phis)
+!
+!-------------------------------------------------------------------------------
+! Purpose: Creates initial states
+!-------------------------------------------------------------------------------
+! Notes:  1) This routine is designed to work for Earth
+!         2) If masque(:,:)/=-99999., masque and phis are already known.
+!         Otherwise: compute it.
+!-------------------------------------------------------------------------------
+  USE control_mod
+  USE fonte_neige_mod
+  USE pbl_surface_mod
+  USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz
+  USE indice_sol_mod
+  USE conf_phys_m, ONLY: conf_phys
+  USE init_ssrf_m, ONLY: start_init_subsurf
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL,    INTENT(INOUT) :: masque(:,:) !--- Land mask           dim(iip1,jjp1)
+  REAL,    INTENT(INOUT) :: phis  (:,:) !--- Ground geopotential dim(iip1,jjp1)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="etat0phys_netcdf", fmt
+  INTEGER            :: i, j, l, ji, iml, jml
+  LOGICAL            :: read_mask
+  REAL               :: phystep, dummy
+  REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp
+  REAL, DIMENSION(klon)               :: sn, rugmer, run_off_lic_0, fder
+  REAL, DIMENSION(klon,nbsrf)         :: qsolsrf, snsrf
+  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil
+
+!--- Arguments for conf_phys
+  LOGICAL :: ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, callstats
+  REAL    :: solarlong0, seuil_inversion, fact_cldcon, facttemps
+  LOGICAL :: ok_newmicro
+  INTEGER :: iflag_radia, iflag_cldcon, iflag_ratqs
+  REAL    :: ratqsbas, ratqshaut, tau_ratqs
+  LOGICAL :: ok_ade, ok_aie, ok_cdnc, aerosol_couple
+  INTEGER :: flag_aerosol
+  LOGICAL :: flag_aerosol_strat
+  LOGICAL :: new_aod
+  REAL    :: bl95_b0, bl95_b1
+  INTEGER :: read_climoz                        !--- Read ozone climatology
+  REAL    :: alp_offset
+
+  deg2rad= pi/180.0
+  iml=assert_eq(SIZE(masque,1),SIZE(phis,1),TRIM(modname)//" iml")
+  jml=assert_eq(SIZE(masque,2),SIZE(phis,2),TRIM(modname)//" jml")
+
+! Physics configuration
+!*******************************************************************************
+  CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,     &
+                   callstats,                                           &
+                   solarlong0,seuil_inversion,                          &
+                   fact_cldcon, facttemps,ok_newmicro,iflag_radia,      &
+                   iflag_cldcon,                                        &
+                   iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,            &
+                   ok_ade, ok_aie, ok_cdnc, aerosol_couple,             &
+                   flag_aerosol, flag_aerosol_strat, new_aod,           &
+                   bl95_b0, bl95_b1,                                    &
+                   read_climoz,                                         &
+                   alp_offset)
+  CALL phys_state_var_init(read_climoz)
+
+!--- Initial atmospheric CO2 conc. from .def file
+  co2_ppm0 = co2_ppm
+
+! Compute ground geopotential, sub-cells quantities and possibly the mask.
+!*******************************************************************************
+  read_mask=ANY(masque/=-99999.); masque_tmp=masque
+  CALL start_init_orog(rlonv, rlatu, phis, masque_tmp)
+  WRITE(fmt,"(i4,'i1)')")iml ; fmt='('//ADJUSTL(fmt)
+  IF(.NOT.read_mask) THEN                       !--- Keep mask form orography
+    masque=masque_tmp
+    IF(prt_level>=1) THEN
+      WRITE(lunout,*)'BUILT MASK :'
+      WRITE(lunout,fmt) NINT(masque)
+    END IF
+    WHERE(   masque(:,:)<EPSFRA) masque(:,:)=0.
+    WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
+  END IF
+  CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid
+
+! Compute tsol and qsol on physical grid, knowing phis on 2D grid.
+!*******************************************************************************
+  CALL start_init_phys(rlonu, rlatv, phis)
+
+! Some initializations.
+!*******************************************************************************
+  sn    (:) = 0.0                               !--- Snow
+  radsol(:) = 0.0                               !--- Net radiation at ground
+  rugmer(:) = 0.001                             !--- Ocean rugosity
+  IF(read_climoz>=1) &                          !--- Ozone climatology
+    CALL regr_lat_time_climoz(read_climoz)
+
+! Sub-surfaces initialization.
+!*******************************************************************************
+  CALL start_init_subsurf(read_mask)
+
+! Write physical initial state
+!*******************************************************************************
+  WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad
+  phystep = dtvr * FLOAT(iphysiq)
+  radpas  = NINT (86400./phystep/ FLOAT(nbapp_rad) )
+  WRITE(lunout,*)'phystep =', phystep, radpas
+
+! Init: ftsol, snsrf, qsolsrf, tsoil, rain_fall, snow_fall, solsw, sollw, z0
+!*******************************************************************************
+  DO i=1,nbsrf; ftsol(:,i) = tsol; END DO
+  DO i=1,nbsrf; snsrf(:,i) = sn;   END DO
+  falb_dir(:,is_ter,:) = 0.08
+  falb_dir(:,is_lic,:) = 0.6
+  falb_dir(:,is_oce,:) = 0.5
+  falb_dir(:,is_sic,:) = 0.6
+  fevap(:,:) = 0.
+  DO i=1,nbsrf; qsolsrf(:,i)=150.; END DO
+  DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO
+  rain_fall  = 0.
+  snow_fall  = 0.
+  solsw      = 165.
+  sollw      = -53.
+  t_ancien   = 273.15
+  q_ancien   = 0.
+  agesno     = 0.
+
+  z0m(:,is_oce) = rugmer(:)
+  z0m(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
+  z0m(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
+  z0m(:,is_sic) = 0.001
+  z0h(:,:)=z0m(:,:)
+
+  fder    = 0.0
+  clwcon  = 0.0
+  rnebcon = 0.0
+  ratqs   = 0.0
+  run_off_lic_0 = 0.0 
+  rugoro  = 0.0
+
+! Before phyredem calling, surface modules and values to be saved in startphy.nc
+! are initialized
+!*******************************************************************************
+  dummy            = 1.0
+  pbl_tke(:,:,:)   = 1.e-8 
+  zmax0(:)         = 40.
+  f0(:)            = 1.e-5
+  sig1(:,:)        = 0.
+  w01(:,:)         = 0.
+  wake_deltat(:,:) = 0.
+  wake_deltaq(:,:) = 0.
+  wake_s(:)        = 0.
+  wake_cstar(:)    = 0.
+  wake_fip(:)      = 0.
+  wake_pe          = 0.
+  fm_therm         = 0.
+  entr_therm       = 0.
+  detr_therm       = 0.
+
+  CALL fonte_neige_init(run_off_lic_0)
+  CALL pbl_surface_init( fder, snsrf, qsolsrf, tsoil )
+  CALL phyredem( "startphy.nc" )
+
+!  WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
+!  WRITE(lunout,*)'entree histclo'
+  CALL histclo()
+
+END SUBROUTINE etat0phys_netcdf
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE start_init_orog(lon_in,lat_in,phis,masque)
+!
+!===============================================================================
+! Comment:
+!   This routine launch grid_noro, which computes parameters for SSO scheme as
+!   described in LOTT & MILLER (1997) and LOTT(1999).
+!===============================================================================
+  USE grid_noro_m, ONLY: grid_noro
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL,    INTENT(IN)    :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
+  REAL,    INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname
+  INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
+  REAL               :: lev(1), date, dt
+  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
+  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:), lat_rel(:,:), tmp_var  (:,:)
+  REAL, ALLOCATABLE  :: zmea0(:,:), zstd0(:,:), zsig0(:,:)
+  REAL, ALLOCATABLE  :: zgam0(:,:), zthe0(:,:), zpic0(:,:), zval0(:,:)
+!-------------------------------------------------------------------------------
+  modname="start_init_orog"
+  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
+  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
+
+!--- HIGH RESOLUTION OROGRAPHY
+  CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
+
+  ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
+  CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
+                lev, ttm_tmp, itau, date, dt, fid)
+  ALLOCATE(relief_hi(iml_rel,jml_rel))
+  CALL flinget(fid, orogvar, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1,1, relief_hi)
+  CALL flinclo(fid)
+
+!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
+  ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
+  lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
+  lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
+
+!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
+  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
+  CALL conf_dat2d(orogvar, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi,.FALSE.)
+  DEALLOCATE(lon_ini,lat_ini)
+
+!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
+  WRITE(lunout,*)
+  WRITE(lunout,*)'*** Compute parameters needed for gravity wave drag code ***'
+
+!--- ALLOCATIONS OF SUB-CELL SCALES QUANTITIES
+  ALLOCATE(zmea0(iml,jml),zstd0(iml,jml)) !--- Mean orography and std deviation
+  ALLOCATE(zsig0(iml,jml),zgam0(iml,jml)) !--- Slope and nisotropy
+  ALLOCATE(zthe0(iml,jml))                !--- Highest slope orientation
+  ALLOCATE(zpic0(iml,jml),zval0(iml,jml)) !--- Peaks and valley heights
+
+!--- CALL OROGRAPHY MODULE TO COMPUTE FIELDS
+  CALL grid_noro(lon_rad,lat_rad,relief_hi,lon_in,lat_in,phis,zmea0,zstd0,     &
+                                      zsig0,zgam0,zthe0,zpic0,zval0,masque)
+  phis = phis * 9.81
+  phis(iml,:) = phis(1,:)
+  DEALLOCATE(relief_hi,lon_rad,lat_rad)
+
+!--- PUT QUANTITIES TO PHYSICAL GRID
+  CALL gr_dyn_fi(1,iml,jml,klon,zmea0,zmea); DEALLOCATE(zmea0)
+  CALL gr_dyn_fi(1,iml,jml,klon,zstd0,zstd); DEALLOCATE(zstd0)
+  CALL gr_dyn_fi(1,iml,jml,klon,zsig0,zsig); DEALLOCATE(zsig0)
+  CALL gr_dyn_fi(1,iml,jml,klon,zgam0,zgam); DEALLOCATE(zgam0)
+  CALL gr_dyn_fi(1,iml,jml,klon,zthe0,zthe); DEALLOCATE(zthe0)
+  CALL gr_dyn_fi(1,iml,jml,klon,zpic0,zpic); DEALLOCATE(zpic0)
+  CALL gr_dyn_fi(1,iml,jml,klon,zval0,zval); DEALLOCATE(zval0)
+
+
+END SUBROUTINE start_init_orog
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE start_init_phys(lon_in,lat_in,phis)
+!
+!===============================================================================
+! Purpose:   Compute tsol and qsol, knowing phis.
+!===============================================================================
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL,    INTENT(IN) :: lon_in(:),  lat_in(:)       ! dim (iml) (jml2)
+  REAL,    INTENT(IN) :: phis(:,:)                   ! dim (iml,jml)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname
+  REAL               :: date, dt
+  INTEGER            :: iml, jml, jml2, itau(1)
+  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), var_ana(:,:)
+  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:)
+  REAL, ALLOCATABLE  :: ts(:,:), qs(:,:)
+!-------------------------------------------------------------------------------
+  modname="start_init_phys"
+  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),TRIM(modname)//" iml")
+  jml=SIZE(phis,2); jml2=SIZE(lat_in)
+
+  WRITE(lunout,*)'Opening the surface analysis'
+  CALL flininfo(phyfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)
+  WRITE(lunout,*) 'Values read: ',  iml_phys, jml_phys, llm_phys, ttm_phys
+
+  ALLOCATE(lat_phys(iml_phys,jml_phys),lon_phys(iml_phys,jml_phys))
+  ALLOCATE(levphys_ini(llm_phys))
+  CALL flinopen(phyfname, .FALSE., iml_phys, jml_phys, llm_phys,              &
+                lon_phys,lat_phys,levphys_ini,ttm_phys,itau,date,dt,fid_phys)
+
+!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
+  ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys))
+  lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad
+  lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad
+
+  ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys))
+  CALL get_var_phys(tsrfvar,ts)                   !--- SURFACE TEMPERATURE
+  CALL get_var_phys(qsolvar,qs)                   !--- SOIL MOISTURE
+  CALL flinclo(fid_phys)
+  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
+
+!--- TSOL AND QSOL ON PHYSICAL GRID
+  ALLOCATE(tsol(klon))
+  CALL gr_dyn_fi(1,iml,jml,klon,ts,tsol)
+  CALL gr_dyn_fi(1,iml,jml,klon,qs,qsol)
+  DEALLOCATE(ts,qs)
+
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE get_var_phys(title,field)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*),  INTENT(IN)    :: title
+  REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: tllm
+!-------------------------------------------------------------------------------
+  SELECT CASE(title)
+    CASE(psrfvar);         tllm=0
+    CASE(tsrfvar,qsolvar); tllm=llm_phys
+  END SELECT
+  IF(ALLOCATED(field)) RETURN
+  ALLOCATE(field(iml,jml)); field(:,:)=0.
+  CALL flinget(fid_phys,title,iml_phys,jml_phys,tllm,ttm_phys,1,1,var_ana)
+  CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
+  CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana,               &
+                                      lon_in,  lat_in,  field)
+
+END SUBROUTINE get_var_phys
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE start_init_phys
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon2,lat2,varo)
+!
+!-------------------------------------------------------------------------------
+  USE inter_barxy_m, ONLY: inter_barxy
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN)  :: nam
+  LOGICAL,          INTENT(IN)  :: ibeg
+  REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
+  REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
+  REAL,             INTENT(IN)  :: lon2(:), lat2(:) ! dim (i1) (j2)
+  REAL,             INTENT(OUT) :: varo(:,:)        ! dim (i1) (j1)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname
+  INTEGER            :: ii, jj, i1, j1, j2
+  REAL, ALLOCATABLE  :: vtmp(:,:)
+!-------------------------------------------------------------------------------
+  modname="interp_startvar"
+  ii=assert_eq(SIZE(lon), SIZE(vari,1),TRIM(modname)//" ii")
+  jj=assert_eq(SIZE(lat), SIZE(vari,2),TRIM(modname)//" jj")
+  i1=assert_eq(SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
+  j1=SIZE(varo,2); j2=SIZE(lat2)
+  ALLOCATE(vtmp(i1-1,j1))
+  IF(ibeg.AND.prt_level>1) THEN
+    WRITE(lunout,*)"--------------------------------------------------------"
+    WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
+    WRITE(lunout,*)"--------------------------------------------------------"
+  END IF
+  CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
+  CALL gr_int_dyn(vtmp, varo, i1-1, j1)
+
+END SUBROUTINE interp_startvar
+!
+!-------------------------------------------------------------------------------
+
+
+END MODULE etat0phys
+!
+!*******************************************************************************
+
Index: dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/init_ssrf_m.F90
===================================================================
--- dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/init_ssrf_m.F90	(revision 3883)
+++ dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/init_ssrf_m.F90	(revision 3883)
@@ -0,0 +1,114 @@
+MODULE init_ssrf_m
+!
+!*******************************************************************************
+
+  USE indice_sol_mod, ONLY: is_ter, is_oce, is_oce, is_lic, epsfra
+  USE dimphy,             ONLY: klon, zmasq
+  USE phys_state_var_mod, ONLY: pctsrf
+  USE geometry_mod, ONLY : longitude_deg=>lon_degrees, latitude_deg=>lat_degrees
+  USE grid_atob_m,        ONLY: grille_m
+  USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo
+
+  CHARACTER(LEN=256), PARAMETER :: icefname="landiceref.nc", icevar="landice"
+  PRIVATE
+  PUBLIC :: start_init_subsurf
+  include "iniprint.h"
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom2.h"
+  include "comconst.h"
+
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE start_init_subsurf(known_mask)
+!
+!-------------------------------------------------------------------------------
+! Purpose: Subsurfaces initialization.
+!-------------------------------------------------------------------------------
+! Comment: Called by etat0phys_netcdf ; also called by limit_netcdf in case
+!          no starting states are required (ok_etat0==.FALSE.).
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  LOGICAL, INTENT(IN) :: known_mask
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER           :: iml_lic, jml_lic
+  INTEGER           :: fid, llm_tmp, ttm_tmp, itaul(1), ji, j
+  REAL, ALLOCATABLE :: dlon_lic(:), lon_lic(:,:), fraclic (:,:)
+  REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:,:), flic_tmp(:,:), vtmp(:,:)
+  REAL              :: date, lev(1), dt, deg2rad
+!-------------------------------------------------------------------------------
+  deg2rad= pi/180.0
+
+!--- Physical grid points coordinates
+  DO j=2,jjm; latitude_deg((j-2)*iim+2:(j-1)*iim+1)=rlatu(j);    END DO
+  DO j=2,jjm; longitude_deg((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:im); END DO
+  latitude_deg(1) = pi/2.; latitude_deg(klon) = - pi/2.
+  latitude_deg(:)=latitude_deg(:)/deg2rad
+  longitude_deg(1) = 0.0;   longitude_deg(klon) = 0.0;
+  longitude_deg(:)=longitude_deg(:)/deg2rad
+
+! Compute ground geopotential, sub-cells quantities and possibly the mask.
+! Sub-surfaces initialization
+!*******************************************************************************
+!--- Read and interpolate on model T-grid soil fraction and soil ice fraction.
+  CALL flininfo(icefname, iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
+  ALLOCATE(lat_lic(iml_lic,jml_lic),lon_lic(iml_lic,jml_lic))
+  ALLOCATE(fraclic(iml_lic,jml_lic))
+  CALL flinopen(icefname, .FALSE., iml_lic, jml_lic, llm_tmp,  &
+ &               lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
+  CALL flinget(fid, icevar, iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic)
+  CALL flinclo(fid)
+  WRITE(lunout,*)'landice dimensions: iml_lic, jml_lic : ',iml_lic,jml_lic
+
+  ALLOCATE(dlon_lic(iml_lic),dlat_lic(jml_lic))
+  dlon_lic(:)=lon_lic(:,1); IF(MAXVAL(dlon_lic)>pi) dlon_lic=dlon_lic*pi/180.
+  dlat_lic(:)=lat_lic(1,:); IF(MAXVAL(dlat_lic)>pi) dlat_lic=dlat_lic*pi/180.
+  DEALLOCATE(lon_lic,lat_lic); ALLOCATE(flic_tmp(iip1,jjp1))
+  CALL grille_m(dlon_lic,dlat_lic,fraclic,rlonv(1:iim),rlatu,flic_tmp(1:iim,:))
+  flic_tmp(iip1,:)=flic_tmp(1,:)
+
+!--- To the physical grid
+  pctsrf(:,:) = 0.
+  CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, pctsrf(:,is_lic))
+  DEALLOCATE(flic_tmp)
+
+!--- Adequation with soil/sea mask
+  WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0. 
+  WHERE(zmasq(:)<EPSFRA)         pctsrf(:,is_lic)=0.
+  pctsrf(:,is_ter)=zmasq(:)
+  DO ji=1,klon
+    IF(zmasq(ji)>EPSFRA) THEN 
+      IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN
+        pctsrf(ji,is_lic)=zmasq(ji)
+        pctsrf(ji,is_ter)=0.
+      ELSE
+        pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic)
+        IF(pctsrf(ji,is_ter)<EPSFRA) THEN
+          pctsrf(ji,is_ter)=0.
+          pctsrf(ji,is_lic)=zmasq(ji)
+        END IF 
+      END IF 
+    END IF 
+  END DO 
+
+!--- Sub-surface ocean and sea ice (sea ice set to zero for start).
+  pctsrf(:,is_oce)=(1.-zmasq(:))
+  WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0.
+  IF(known_mask) pctsrf(:,is_oce)=1-zmasq(:)
+
+!--- It is checked that the sub-surfaces sum is equal to 1.
+  ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA)
+  IF(ji/=0) WRITE(lunout,*) 'Sub-cell distribution problem for ',ji,' points'
+
+END SUBROUTINE start_init_subsurf
+!
+!-------------------------------------------------------------------------------
+
+END MODULE init_ssrf_m
+!
+!*******************************************************************************
Index: dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/limit_netcdf.F90
===================================================================
--- dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/limit_netcdf.F90	(revision 3883)
+++ dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/limit_netcdf.F90	(revision 3883)
@@ -0,0 +1,864 @@
+MODULE limit
+!
+!*******************************************************************************
+! Author : L. Fairhead, 27/01/94
+!-------------------------------------------------------------------------------
+! Purpose: Boundary conditions files building for new model using climatologies.
+!          Both grids have to be regular.
+!-------------------------------------------------------------------------------
+! Note: This routine is designed to work for Earth
+!-------------------------------------------------------------------------------
+! Modification history:
+!  * 23/03/1994: Z. X. Li
+!  *    09/1999: L. Fairhead (netcdf reading in LMDZ.3.3)
+!  *    07/2001: P. Le Van
+!  *    11/2009: L. Guez     (ozone day & night climatos, see etat0_netcdf.F90)
+!  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
+!-------------------------------------------------------------------------------
+
+  USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo,          &
+   ioconf_calendar, ioget_calendar, lock_calendar, ioget_mon_len, ioget_year_len
+  USE assert_eq_m,        ONLY: assert_eq
+  USE conf_dat_m,         ONLY: conf_dat2d, conf_dat3d
+  USE dimphy,             ONLY: klon, zmasq
+  USE geometry_mod, ONLY: longitude_deg=>lon_degrees, latitude_deg=>lat_degrees
+  USE phys_state_var_mod, ONLY: pctsrf
+  USE control_mod, ONLY: anneeref
+  USE init_ssrf_m,        ONLY: start_init_subsurf
+
+  CHARACTER(LEN=20), PARAMETER :: &
+  fsst(4)=['amipbc_sst_1x1.nc   ','cpl_atm_sst.nc      ','histmth_sst.nc      '&
+          ,'sstk.nc             ']
+  CHARACTER(LEN=20), PARAMETER :: &
+  fsic(4)=['amipbc_sic_1x1.nc   ','cpl_atm_sic.nc      ','histmth_sic.nc      '&
+          ,'ci.nc               ']
+  CHARACTER(LEN=10), PARAMETER :: &
+  vsst(4)=['tosbcs    ','SISUTESW  ','tsol_oce  ','sstk      '], &
+  vsic(4)=['sicbcs    ','SIICECOV  ','pourc_sic ','ci        ']
+  CHARACTER(LEN=10), PARAMETER :: &
+  frugo='Rugos.nc  ', falbe='Albedo.nc ', frelf='Relief.nc ',    &
+   vrug='RUGOS     ',  valb='ALBEDO    ',  vrel='RELIEF    '
+
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE limit_netcdf(masque, phis, extrap)
+!
+!-------------------------------------------------------------------------------
+! Author : L. Fairhead, 27/01/94
+!-------------------------------------------------------------------------------
+! Purpose: Boundary conditions files building for new model using climatologies.
+!          Both grids have to be regular.
+!-------------------------------------------------------------------------------
+! Note: This routine is designed to work for Earth
+!-------------------------------------------------------------------------------
+! Modification history:
+!  * 23/03/1994: Z. X. Li
+!  *    09/1999: L. Fairhead (netcdf reading in LMDZ.3.3)
+!  *    07/2001: P. Le Van
+!  *    11/2009: L. Guez     (ozone day & night climatos, see etat0_netcdf.F90)
+!  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
+!-------------------------------------------------------------------------------
+#ifndef CPP_1D
+  USE indice_sol_mod
+  USE netcdf,             ONLY: NF90_OPEN,    NF90_CREATE,  NF90_CLOSE,        &
+                  NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT,      &
+                  NF90_NOERR,   NF90_NOWRITE, NF90_DOUBLE,  NF90_GLOBAL,       &
+		  NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED, NF90_FLOAT
+  USE inter_barxy_m,      ONLY: inter_barxy
+  USE netcdf95,           ONLY: nf95_def_var, nf95_put_att, nf95_put_var
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  include "iniprint.h"
+  include "dimensions.h"
+  include "paramet.h"
+  REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
+  REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: phis   ! ground geopotential
+  LOGICAL,                    INTENT(IN)    :: extrap ! SST extrapolation flag
+!-------------------------------------------------------------------------------
+! Local variables:
+  include "logic.h"
+  include "comgeom2.h"
+  include "comconst.h"
+
+!--- INPUT NETCDF FILES NAMES --------------------------------------------------
+  CHARACTER(LEN=20) :: icefile, sstfile, dumstr, fnam
+  CHARACTER(LEN=10) :: varname
+
+!--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
+  REAL               :: fi_ice(klon), verif(klon)
+  REAL, POINTER      :: phy_rug(:,:)=>NULL(), phy_ice(:,:)=>NULL()
+  REAL, POINTER      :: phy_sst(:,:)=>NULL(), phy_alb(:,:)=>NULL()
+  REAL, ALLOCATABLE  :: phy_bil(:,:), pctsrf_t(:,:,:)
+  INTEGER            :: nbad
+
+!--- VARIABLES FOR OUTPUT FILE WRITING -----------------------------------------
+  INTEGER :: ierr, nid, ndim, ntim, k, dims(2), ix_sic, ix_sst
+  INTEGER :: id_tim,  id_SST,  id_BILS, id_RUG, id_ALB
+  INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC, varid_longitude, varid_latitude
+  INTEGER :: NF90_FORMAT
+  INTEGER :: ndays                   !--- Depending on the output calendar
+
+!--- INITIALIZATIONS -----------------------------------------------------------
+#ifdef NC_DOUBLE
+  NF90_FORMAT=NF90_DOUBLE
+#else
+  NF90_FORMAT=NF90_FLOAT
+#endif
+  CALL inigeom
+
+!--- MASK, GROUND GEOPOT. & SUBSURFACES COMPUTATION (IN CASE ok_etat0==.FALSE.)
+   IF(ALL(masque==-99999.)) THEN
+    CALL start_init_orog0(rlonv,rlatu,phis,masque)
+    CALL gr_dyn_fi(1,iip1,jjp1,klon,masque,zmasq)          !--- To physical grid
+    ALLOCATE(pctsrf(klon,nbsrf))
+    CALL start_init_subsurf(.FALSE.)
+  END IF
+
+!--- Beware: anneeref (from gcm.def) is used to determine output time sampling
+  ndays=ioget_year_len(anneeref)
+
+!--- RUGOSITY TREATMENT --------------------------------------------------------
+  CALL msg(1,'Traitement de la rugosite')
+  CALL get_2Dfield(frugo,vrug,'RUG',ndays,phy_rug,mask=masque(1:iim,:))
+
+!--- OCEAN TREATMENT -----------------------------------------------------------
+  CALL msg(1,'Traitement de la glace oceanique')
+
+! Input SIC file selection
+! Open file only to test if available
+  DO ix_sic=1,SIZE(fsic)
+     IF ( NF90_OPEN(TRIM(fsic(ix_sic)),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+        icefile=fsic(ix_sic); varname=vsic(ix_sic); EXIT
+     END IF
+  END DO
+  IF(ix_sic==SIZE(fsic)+1) THEN
+     WRITE(lunout,*) 'ERROR! No sea-ice input file was found.'
+     WRITE(lunout,*) 'One of following files must be available : '
+     DO k=1,SIZE(fsic); WRITE(lunout,*) TRIM(fsic(k)); END DO
+     CALL abort_physic('limit_netcdf','No sea-ice file was found',1)
+  END IF
+  CALL ncerr(NF90_CLOSE(nid),icefile)
+  CALL msg(-1,'Fichier choisi pour la glace de mer:'//TRIM(icefile))
+
+  CALL get_2Dfield(icefile,varname, 'SIC',ndays,phy_ice)
+
+  ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
+  DO k=1,ndays
+     fi_ice=phy_ice(:,k)
+     WHERE(fi_ice>=1.0  ) fi_ice=1.0
+     WHERE(fi_ice<EPSFRA) fi_ice=0.0
+     pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
+     pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
+     SELECT CASE(ix_sic)
+        CASE(2)                                   ! SIC=pICE*(1-LIC-TER) (CPL)
+        pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
+        CASE(3)                                   ! SIC=pICE            (HIST)
+        pctsrf_t(:,is_sic,k)=fi_ice(:)
+        CASE DEFAULT                              ! SIC=pICE-LIC   (AMIP,ERAI)
+        pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
+     END SELECT
+     WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
+     WHERE(1.0-zmasq<EPSFRA)
+        pctsrf_t(:,is_sic,k)=0.0
+        pctsrf_t(:,is_oce,k)=0.0
+     ELSEWHERE
+        WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
+           pctsrf_t(:,is_sic,k)=1.0-zmasq
+           pctsrf_t(:,is_oce,k)=0.0
+        ELSEWHERE
+           pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
+           WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
+              pctsrf_t(:,is_oce,k)=0.0
+              pctsrf_t(:,is_sic,k)=1.0-zmasq
+           END WHERE
+        END WHERE
+     END WHERE
+     nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
+     IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb points = ',nbad
+     nbad=COUNT(ABS(SUM(pctsrf_t(:,:,k),DIM=2)-1.0)>EPSFRA)
+     IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
+  END DO
+  DEALLOCATE(phy_ice)
+
+!--- SST TREATMENT -------------------------------------------------------------
+  CALL msg(1,'Traitement de la sst')
+
+! Input SST file selection
+! Open file only to test if available
+  DO ix_sst=1,SIZE(fsst)
+     IF ( NF90_OPEN(TRIM(fsst(ix_sst)),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+       sstfile=fsst(ix_sst); varname=vsst(ix_sst); EXIT
+     END IF
+  END DO
+  IF(ix_sst==SIZE(fsst)+1) THEN
+     WRITE(lunout,*) 'ERROR! No sst input file was found.'
+     WRITE(lunout,*) 'One of following files must be available : '
+     DO k=1,SIZE(fsst); WRITE(lunout,*) TRIM(fsst(k)); END DO
+     CALL abort_physic('limit_netcdf','No sst file was found',1)
+  END IF
+  CALL ncerr(NF90_CLOSE(nid),sstfile)
+  CALL msg(-1,'Fichier choisi pour la temperature de mer: '//TRIM(sstfile))
+
+  CALL get_2Dfield(sstfile,varname,'SST',ndays,phy_sst,flag=extrap)
+
+!--- ALBEDO TREATMENT ----------------------------------------------------------
+  CALL msg(1,'Traitement de l albedo')
+  CALL get_2Dfield(falbe,valb,'ALB',ndays,phy_alb)
+
+!--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
+  ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0
+
+!--- OUTPUT FILE WRITING -------------------------------------------------------
+  CALL msg(5,'Ecriture du fichier limit : debut')
+  fnam="limit.nc"
+
+  !--- File creation
+  CALL ncerr(NF90_CREATE(fnam,NF90_CLOBBER,nid),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier conditions aux limites"),fnam)
+
+  !--- Dimensions creation
+  CALL ncerr(NF90_DEF_DIM(nid,"points_physiques",klon,ndim),fnam)
+  CALL ncerr(NF90_DEF_DIM(nid,"time",NF90_UNLIMITED,ntim),fnam)
+
+  dims=[ndim,ntim]
+
+  !--- Variables creation
+  CALL ncerr(NF90_DEF_VAR(nid,"TEMPS",NF90_FORMAT,[ntim],id_tim),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"FOCE", NF90_FORMAT,dims,id_FOCE),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"FSIC", NF90_FORMAT,dims,id_FSIC),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"FTER", NF90_FORMAT,dims,id_FTER),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"FLIC", NF90_FORMAT,dims,id_FLIC),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"SST",  NF90_FORMAT,dims,id_SST),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"BILS", NF90_FORMAT,dims,id_BILS),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"ALB",  NF90_FORMAT,dims,id_ALB),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"RUG",  NF90_FORMAT,dims,id_RUG),fnam)
+  call nf95_def_var(nid, "longitude", NF90_FLOAT, ndim, varid_longitude)
+  call nf95_def_var(nid, "latitude",  NF90_FLOAT, ndim, varid_latitude)
+
+  !--- Attributes creation
+  CALL ncerr(NF90_PUT_ATT(nid,id_tim, "title","Jour dans l annee"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_FOCE,"title","Fraction ocean"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_FSIC,"title","Fraction glace de mer"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_FTER,"title","Fraction terre"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_FLIC,"title","Fraction land ice"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_SST ,"title","Temperature superficielle de la mer"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_BILS,"title","Reference flux de chaleur au sol"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_ALB, "title","Albedo a la surface"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_RUG, "title","Rugosite"),fnam)
+
+  call nf95_put_att(nid, varid_longitude, "standard_name", "longitude")
+  call nf95_put_att(nid, varid_longitude, "units", "degrees_east")
+
+  call nf95_put_att(nid, varid_latitude, "standard_name", "latitude")
+  call nf95_put_att(nid, varid_latitude, "units", "degrees_north")
+
+  CALL ncerr(NF90_ENDDEF(nid),fnam)
+
+  !--- Variables saving
+  CALL ncerr(NF90_PUT_VAR(nid,id_tim,[(REAL(k),k=1,ndays)]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_FOCE,pctsrf_t(:,is_oce,:),[1,1],[klon,ndays]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_FSIC,pctsrf_t(:,is_sic,:),[1,1],[klon,ndays]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_FTER,pctsrf_t(:,is_ter,:),[1,1],[klon,ndays]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_FLIC,pctsrf_t(:,is_lic,:),[1,1],[klon,ndays]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_SST ,phy_sst(:,:),[1,1],[klon,ndays]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_BILS,phy_bil(:,:),[1,1],[klon,ndays]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_ALB ,phy_alb(:,:),[1,1],[klon,ndays]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_RUG ,phy_rug(:,:),[1,1],[klon,ndays]),fnam)
+  call nf95_put_var(nid, varid_longitude, longitude_deg)
+  call nf95_put_var(nid, varid_latitude, latitude_deg)
+
+  CALL ncerr(NF90_CLOSE(nid),fnam)
+
+  CALL msg(6,'Ecriture du fichier limit : fin')
+
+  DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
+
+
+!===============================================================================
+!
+  CONTAINS
+!
+!===============================================================================
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE get_2Dfield(fnam, varname, mode, ndays, champo, flag, mask)
+!
+!-----------------------------------------------------------------------------
+! Comments:
+!   There are two assumptions concerning the NetCDF files, that are satisfied
+!   with files that are conforming NC convention:
+!     1) The last dimension of the variables used is the time record.
+!     2) Dimensional variables have the same names as corresponding dimensions.
+!-----------------------------------------------------------------------------
+  USE netcdf, ONLY: NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, &
+       NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_GET_VAR, &
+       NF90_GET_ATT
+  USE pchsp_95_m, only: pchsp_95
+  USE pchfe_95_m, only: pchfe_95
+  USE arth_m, only: arth
+  USE indice_sol_mod
+
+  IMPLICIT NONE
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom2.h"
+!-----------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
+  CHARACTER(LEN=10), INTENT(IN)     :: varname  ! NetCDF variable name
+  CHARACTER(LEN=3),  INTENT(IN)     :: mode     ! RUG, SIC, SST or ALB
+  INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
+  REAL,    POINTER,  DIMENSION(:, :) :: champo  ! output field = f(t)
+  LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST) old ice (SIC)
+  REAL,    OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
+!------------------------------------------------------------------------------
+! Local variables:
+!--- NetCDF
+  INTEGER           :: ncid, varid        ! NetCDF identifiers
+  CHARACTER(LEN=30) :: dnam               ! dimension name
+!--- dimensions
+  INTEGER           :: dids(4)            ! NetCDF dimensions identifiers
+  REAL, ALLOCATABLE :: dlon_ini(:)        ! initial longitudes vector
+  REAL, ALLOCATABLE :: dlat_ini(:)        ! initial latitudes  vector
+  REAL, POINTER     :: dlon(:), dlat(:)   ! reordered lon/lat  vectors
+!--- fields
+  INTEGER :: imdep, jmdep, lmdep          ! dimensions of 'champ'
+  REAL, ALLOCATABLE :: champ(:,:)         ! wanted field on initial grid
+  REAL, ALLOCATABLE :: yder(:), timeyear(:)
+  REAL              :: champint(iim,jjp1) ! interpolated field
+  REAL, ALLOCATABLE :: champtime(:,:,:)
+  REAL, ALLOCATABLE :: champan(:,:,:)
+!--- input files
+  CHARACTER(LEN=20) :: cal_in             ! calendar
+  CHARACTER(LEN=20) :: unit_sic           ! attribute unit in sea-ice file
+  INTEGER           :: ndays_in           ! number of days
+!--- misc
+  INTEGER           :: i, j, k, l         ! loop counters
+  REAL, ALLOCATABLE :: work(:,:)          ! used for extrapolation
+  CHARACTER(LEN=25) :: title              ! for messages
+  LOGICAL           :: extrp              ! flag for extrapolation
+  REAL              :: chmin, chmax
+  INTEGER ierr
+  integer n_extrap ! number of extrapolated points
+  logical skip
+
+!------------------------------------------------------------------------------
+!---Variables depending on keyword 'mode' -------------------------------------
+  NULLIFY(champo)
+
+  SELECT CASE(mode)
+  CASE('RUG'); title='Rugosite'
+  CASE('SIC'); title='Sea-ice'
+  CASE('SST'); title='SST'
+  CASE('ALB'); title='Albedo'
+  END SELECT
+  extrp=.FALSE.; IF(PRESENT(flag).AND.mode=='SST') extrp=flag
+
+!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
+  CALL msg(5,' Now reading file : '//TRIM(fnam))
+  CALL ncerr(NF90_OPEN(fnam, NF90_NOWRITE, ncid),fnam)
+  CALL ncerr(NF90_INQ_VARID(ncid, trim(varname), varid),fnam)
+  CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam)
+
+!--- Read unit for sea-ice variable only
+  IF (mode=='SIC') THEN
+    IF(NF90_GET_ATT(ncid, varid, 'units', unit_sic)/=NF90_NOERR) THEN
+      CALL msg(5,'No unit in sea-ice file. Take percentage as default value')
+      unit_sic='X'
+    ELSE
+      CALL msg(5,'Sea-ice cover has unit='//TRIM(unit_sic))
+    END IF
+  END IF
+
+!--- Longitude
+  CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep),fnam)
+  ALLOCATE(dlon_ini(imdep), dlon(imdep))
+  CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
+  CALL ncerr(NF90_GET_VAR(ncid, varid, dlon_ini), fnam)
+  CALL msg(5,'variable '//TRIM(dnam)//' dimension ', imdep)
+
+!--- Latitude
+  CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep),fnam)
+  ALLOCATE(dlat_ini(jmdep), dlat(jmdep))
+  CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
+  CALL ncerr(NF90_GET_VAR(ncid, varid, dlat_ini), fnam)
+  CALL msg(5,'variable '//TRIM(dnam)//' dimension ', jmdep)
+
+!--- Time (variable is not needed - it is rebuilt - but calendar is)
+  CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep), fnam)
+  ALLOCATE(timeyear(lmdep))
+  CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
+  cal_in=' '
+  IF(NF90_GET_ATT(ncid, varid, 'calendar', cal_in)/=NF90_NOERR) THEN
+    SELECT CASE(mode)
+      CASE('RUG', 'ALB'); cal_in='360d'
+      CASE('SIC', 'SST'); cal_in='gregorian'
+    END SELECT
+  CALL msg(5,'WARNING: missing "calendar" attribute for "time" in '&
+     &//TRIM(fnam)//'. Choosing default value.')
+  END IF
+  CALL msg(5,'var, calendar, dim: '//TRIM(dnam)//' '//TRIM(cal_in), lmdep)
+  
+!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
+  !--- Determining input file number of days, depending on calendar
+  ndays_in=year_len(anneeref, cal_in)
+
+!--- Time vector reconstruction (time vector from file is not trusted)
+!--- If input records are not monthly, time sampling has to be constant !
+  timeyear=mid_months(anneeref, cal_in, lmdep)
+  IF (lmdep /= 12) WRITE(lunout,*) 'Note : les fichiers de ', TRIM(mode), &
+       ' ne comportent pas 12, mais ', lmdep, ' enregistrements.'
+
+!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
+  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep))
+  IF(extrp) ALLOCATE(work(imdep, jmdep))
+  CALL msg(5,'')
+  CALL msg(5,'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.')
+  CALL ncerr(NF90_INQ_VARID(ncid, varname, varid), fnam)
+  DO l=1, lmdep
+    CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam)
+    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
+    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
+    IF(l==1) THEN
+      CALL msg(5,"----------------------------------------------------------")
+      CALL msg(5,"$$$ Interpolation barycentrique pour "//TRIM(title)//" $$$")
+      CALL msg(5,"----------------------------------------------------------")
+    END IF
+    IF(mode=='RUG') champ=LOG(champ)
+    CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
+    IF(mode=='RUG') THEN
+      champint=EXP(champint)
+      WHERE(NINT(mask)/=1) champint=0.001
+    END IF
+    champtime(:, :, l)=champint
+  END DO
+  CALL ncerr(NF90_CLOSE(ncid), fnam)
+
+  DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
+  IF(extrp) DEALLOCATE(work)
+
+!--- TIME INTERPOLATION ------------------------------------------------------
+  IF(prt_level>5) THEN
+     WRITE(lunout, *)
+     WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
+     WRITE(lunout, *)' Vecteur temps en entree: ', timeyear
+     WRITE(lunout, *)' Vecteur temps en sortie de 0 a ', ndays
+  END IF
+
+  ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
+  skip = .false.
+  n_extrap = 0
+  DO j=1, jjp1
+    DO i=1, iim
+      yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
+           vc_beg=0., vc_end=0.)
+      CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
+           arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
+      if (ierr < 0) stop 1
+      n_extrap = n_extrap + ierr
+    END DO
+  END DO
+  if (n_extrap /= 0) then
+     WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
+  end if
+  champan(iip1, :, :)=champan(1, :, :)
+  DEALLOCATE(yder, champtime, timeyear)
+
+!--- Checking the result
+  DO j=1, jjp1
+    CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
+    IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' au temps 10 ', chmin, chmax, j
+  END DO
+
+!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
+  IF(mode=='SST') THEN
+    CALL msg(5,'Filtrage de la SST: SST >= 271.38')
+    WHERE(champan<271.38) champan=271.38
+  END IF
+
+!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
+  IF(mode=='SIC') THEN
+    CALL msg(5,'Filtrage de la SIC: 0.0 < Sea-ice < 1.0')
+
+    IF (unit_sic=='1') THEN
+       ! Nothing to be done for sea-ice field is already in fraction of 1
+       ! This is the case for sea-ice in file cpl_atm_sic.nc
+       CALL msg(5,'Sea-ice field already in fraction of 1')
+    ELSE
+       ! Convert sea ice from percentage to fraction of 1
+       CALL msg(5,'Transformt sea-ice field from percentage to fraction of 1.')
+       champan(:, :, :)=champan(:, :, :)/100.
+    END IF
+
+    champan(iip1, :, :)=champan(1, :, :)
+    WHERE(champan>1.0) champan=1.0
+    WHERE(champan<0.0) champan=0.0
+ END IF
+
+!--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
+  ALLOCATE(champo(klon, ndays))
+  DO k=1, ndays
+    CALL gr_dyn_fi(1, iip1, jjp1, klon, champan(1, 1, k), champo(1, k))
+  END DO
+  DEALLOCATE(champan)
+
+END SUBROUTINE get_2Dfield
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!===============================================================================
+! Purpose:  Compute "phis" just like it would be in start_init_orog.
+!===============================================================================
+! Arguments:
+  REAL,             INTENT(IN)    :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
+  REAL,             INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="start_init_orog0"
+  INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
+  REAL               :: lev(1), date, dt, deg2rad
+  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
+  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:), lat_rel(:,:)
+!-------------------------------------------------------------------------------
+  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
+  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
+  IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1)
+  IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1)
+  pi=2.0*ASIN(1.0); deg2rad=pi/180.0
+  IF(ANY(phis/=-99999.)) RETURN                  !--- phis ALREADY KNOWN
+
+!--- HIGH RESOLUTION OROGRAPHY
+  CALL flininfo(frelf, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
+
+  ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
+  CALL flinopen(frelf, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,   &
+                lev, ttm_tmp, itau, date, dt, fid)
+  ALLOCATE(relief_hi(iml_rel,jml_rel))
+  CALL flinget(fid, vrel, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
+  CALL flinclo(fid)
+
+!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
+  ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
+  lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
+  lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
+
+!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
+  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
+  CALL conf_dat2d(vrel, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
+  DEALLOCATE(lon_ini,lat_ini)
+
+!--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0
+  WRITE(lunout,*)
+  WRITE(lunout,*)'*** Compute surface geopotential ***'
+
+!--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS
+  CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque)
+  phis = phis * 9.81
+  phis(iml,:) = phis(1,:)
+  DEALLOCATE(relief_hi,lon_rad,lat_rad)
+
+END SUBROUTINE start_init_orog0
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
+!
+!===============================================================================
+! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
+!          without any call to physics subroutines.
+!===============================================================================
+  IMPLICIT NONE 
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)   :: xd(:), yd(:) !--- INPUT  COORDINATES     (imdp) (jmdp)
+  REAL, INTENT(IN)   :: zd(:,:)      !--- INPUT  FIELD           (imdp,jmdp)
+  REAL, INTENT(IN)   :: x(:), y(:)   !--- OUTPUT COORDINATES     (imar+1) (jmar)
+  REAL, INTENT(OUT)  :: zphi(:,:)    !--- GEOPOTENTIAL           (imar+1,jmar)
+  REAL, INTENT(INOUT):: mask(:,:)    !--- MASK                   (imar+1,jmar)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="grid_noro0"
+  REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
+  REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
+  REAL, ALLOCATABLE :: weight(:,:)                ! dim (imar+1,jmar)
+  REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:)   ! dim (imar+1,jmar)
+  REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
+  REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
+  REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
+  LOGICAL :: masque_lu
+  INTEGER :: i, ii, imdp, imar, iext
+  INTEGER :: j, jj, jmdp, jmar, nn
+  REAL    :: xpi, zlenx, weighx, xincr,  zbordnor, zmeanor, zweinor, zbordest
+  REAL    :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
+!-------------------------------------------------------------------------------
+  imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
+  jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
+  imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
+  jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
+  IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
+  IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
+  iext=imdp/10
+  xpi = ACOS(-1.)
+  rad = 6371229.
+
+!--- ARE WE USING A READ MASK ?
+  masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
+  WRITE(lunout,*)'Masque lu: ',masque_lu
+
+!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
+  ALLOCATE(xusn(imdp+2*iext))
+  xusn(1     +iext:imdp  +iext)=xd(:)
+  xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
+  xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
+
+  ALLOCATE(yusn(jmdp+2))
+  yusn(1       )=yd(1)   +(yd(1)   -yd(2))
+  yusn(2:jmdp+1)=yd(:)
+  yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
+
+  ALLOCATE(zusn(imdp+2*iext,jmdp+2))
+  zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
+  zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
+  zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
+  zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
+  zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
+  zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
+  zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
+
+!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
+  ALLOCATE(a(imar+1),b(imar+1))
+  b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
+  b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
+  a(1)=x(1)-(x(2)-x(1))/2.0
+  a(2:imar+1)= b(1:imar)
+
+  ALLOCATE(c(jmar),d(jmar))
+  d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
+  d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
+  c(1)=y(1)-(y(2)-y(1))/2.0
+  c(2:jmar)=d(1:jmar-1)
+
+!--- INITIALIZATIONS:
+  ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
+  ALLOCATE(zmea  (imar+1,jmar)); zmea  (:,:)= 0.0
+
+!--- SUMMATION OVER GRIDPOINT AREA
+  zleny=xpi/REAL(jmdp)*rad
+  xincr=xpi/REAL(jmdp)/2.
+  ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
+  ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
+  DO ii = 1, imar+1
+    DO jj = 1, jmar
+      DO j = 2,jmdp+1 
+        zlenx  =zleny  *COS(yusn(j))
+        zbordnor=(xincr+c(jj)-yusn(j))*rad
+        zbordsud=(xincr-d(jj)+yusn(j))*rad
+        weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
+        IF(weighy/=0) THEN
+          DO i = 2, imdp+2*iext-1
+            zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
+            zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
+            weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
+            IF(weighx/=0)THEN
+              num_tot(ii,jj)=num_tot(ii,jj)+1.0
+              IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
+              weight(ii,jj)=weight(ii,jj)+weighx*weighy
+              zmea  (ii,jj)=zmea  (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
+            END IF
+          END DO
+        END IF
+      END DO
+    END DO
+  END DO
+
+!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
+  IF(.NOT.masque_lu) THEN
+    WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
+  END IF
+  nn=COUNT(weight(:,1:jmar-1)==0.0)
+  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
+  WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
+
+!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
+  ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
+  WHERE(mask>=0.1) mask_tmp = 1.
+  WHERE(weight(:,:)/=0.0)
+    zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
+    zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
+  END WHERE
+  WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
+
+!--- Values at poles
+  zphi(imar+1,:)=zphi(1,:)
+
+  zweinor=SUM(weight(1:imar,   1),DIM=1)
+  zweisud=SUM(weight(1:imar,jmar),DIM=1)
+  zmeanor=SUM(weight(1:imar,   1)*zmea(1:imar,   1),DIM=1)
+  zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
+  zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
+
+END SUBROUTINE grid_noro0
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+FUNCTION year_len(y,cal_in)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  INTEGER                       :: year_len
+  INTEGER,           INTENT(IN) :: y
+  CHARACTER(LEN=*),  INTENT(IN) :: cal_in
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=20)             :: cal_out              ! calendar (for outputs)
+!-------------------------------------------------------------------------------
+!--- Getting the input calendar to reset at the end of the function
+  CALL ioget_calendar(cal_out)
+
+!--- Unlocking calendar and setting it to wanted one
+  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
+
+!--- Getting the number of days in this year
+  year_len=ioget_year_len(y)
+
+!--- Back to original calendar
+  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
+
+END FUNCTION year_len
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+FUNCTION mid_months(y,cal_in,nm)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  INTEGER,                INTENT(IN) :: y               ! year
+  CHARACTER(LEN=*),       INTENT(IN) :: cal_in          ! calendar
+  INTEGER,                INTENT(IN) :: nm              ! months/year number
+  REAL,    DIMENSION(nm)             :: mid_months      ! mid-month times
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=99)                  :: mess            ! error message
+  CHARACTER(LEN=20)                  :: cal_out         ! calendar (for outputs)
+  INTEGER, DIMENSION(nm)             :: mnth            ! months lengths (days)
+  INTEGER                            :: m               ! months counter
+  INTEGER                            :: nd              ! number of days
+!-------------------------------------------------------------------------------
+  nd=year_len(y,cal_in)
+
+  IF(nm==12) THEN
+
+  !--- Getting the input calendar to reset at the end of the function
+    CALL ioget_calendar(cal_out)
+
+  !--- Unlocking calendar and setting it to wanted one
+    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
+
+  !--- Getting the length of each month
+    DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
+
+  !--- Back to original calendar
+    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
+
+  ELSE IF(MODULO(nd,nm)/=0) THEN
+    WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
+      nm,' months/year. Months number should divide days number.'
+    CALL abort_physic('mid_months',TRIM(mess),1)
+
+  ELSE
+    mnth=[(m,m=1,nm,nd/nm)]
+  END IF
+
+!--- Mid-months times
+  mid_months(1)=0.5*REAL(mnth(1))
+  DO k=2,nm
+    mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
+  END DO
+
+END FUNCTION mid_months
+!
+!-------------------------------------------------------------------------------
+
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE msg(lev,str1,i,str2)
+!
+!-------------------------------------------------------------------------------
+! Arguments:
+  INTEGER,                    INTENT(IN) :: lev
+  CHARACTER(LEN=*),           INTENT(IN) :: str1
+  INTEGER,          OPTIONAL, INTENT(IN) :: i
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: str2
+!-------------------------------------------------------------------------------
+  IF(prt_level>lev) THEN
+    IF(PRESENT(str2)) THEN
+      WRITE(lunout,*) TRIM(str1), i, TRIM(str2)
+    ELSE IF(PRESENT(i)) THEN
+      WRITE(lunout,*) TRIM(str1), i
+    ELSE
+      WRITE(lunout,*) TRIM(str1)
+    END IF
+  END IF
+
+END SUBROUTINE msg
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE ncerr(ncres,fnam)
+!
+!-------------------------------------------------------------------------------
+! Purpose: NetCDF errors handling.
+!-------------------------------------------------------------------------------
+  USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  INTEGER,          INTENT(IN) :: ncres
+  CHARACTER(LEN=*), INTENT(IN) :: fnam
+!-------------------------------------------------------------------------------
+  IF(ncres/=NF90_NOERR) THEN
+    WRITE(lunout,*)'Problem with file '//TRIM(fnam)//' in routine limit_netcdf.'
+    CALL abort_physic('limit_netcdf',NF90_STRERROR(ncres),1)
+  END IF
+
+END SUBROUTINE ncerr
+!
+!-------------------------------------------------------------------------------
+
+#endif
+! of #ifndef CPP_1D
+END SUBROUTINE limit_netcdf
+
+END MODULE limit
+!
+!*******************************************************************************
Index: dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/test_disvert_m.F90
===================================================================
--- dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/test_disvert_m.F90	(revision 3883)
+++ dynamico_lmdz/aquaplanet/LMDZ5/libf/dynlonlat_phylonlat/phylmd/test_disvert_m.F90	(revision 3883)
@@ -0,0 +1,68 @@
+module test_disvert_m
+
+  implicit none
+
+contains
+
+  subroutine test_disvert
+
+    ! Author: Lionel GUEZ
+
+    ! This procedure tests the order of pressure values at half-levels
+    ! and full levels. We arbitrarily choose to test ngrid values of
+    ! the surface pressure, which sample possible values on Earth.
+
+    use exner_hyb_m, only: exner_hyb
+    use vertical_layers_mod, only: ap,bp,preff
+
+    ! For llm:
+    include "dimensions.h"
+
+    ! For kappa, cpp:
+    include "comconst.h"
+
+    ! Local:
+    integer l, i
+    integer, parameter:: ngrid = 7
+    real p(ngrid, llm + 1) ! pressure at half-level, in Pa
+    real pks(ngrid) ! exner function at the surface, in J K-1 kg-1 
+    real pk(ngrid, llm) ! exner function at full level, in J K-1 kg-1 
+    real ps(ngrid) ! surface pressure, in Pa
+    real p_lay(ngrid, llm) ! pressure at full level, in Pa
+    real delta_ps ! in Pa
+
+    !---------------------
+
+    print *, "Call sequence information: test_disvert"
+
+    delta_ps = 6e4 / (ngrid - 1)
+    ps = (/(5e4 + delta_ps * i, i = 0, ngrid - 1)/)
+    forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
+    call exner_hyb(ngrid, ps, p, pks, pk)
+    p_lay = preff * (pk / cpp)**(1. / kappa)
+
+    ! Are pressure values in the right order?
+    if (any(p(:, :llm) <= p_lay .or. p_lay <= p(:, 2:))) then
+       ! List details and stop:
+       do l = 1, llm
+          do i = 1, ngrid
+             if (p(i, l) <= p_lay(i, l)) then
+                print 1000, "ps = ", ps(i) / 100., "hPa, p(level ",  l, &
+                     ") = ", p(i, l) / 100., " hPa <= p(layer ", l, ") = ", &
+                     p_lay(i, l) / 100., " hPa"
+             end if
+             if (p_lay(i, l) <= p(i, l + 1)) then
+                print 1000, "ps = ", ps(i) / 100., "hPa, p(layer ", l, ") = ", &
+                     p_lay(i, l) / 100., " hPa <= p(level ", l + 1, ") = ", &
+                     p(i, l + 1) / 100., " hPa"
+             end if
+          end do
+       end do
+       call abort_physic("test_disvert", "bad order of pressure values", 1)
+    end if
+
+1000 format (3(a, g10.4, a, i0))
+
+  end subroutine test_disvert
+
+end module test_disvert_m
Index: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/grid_noro_m.F90
===================================================================
--- dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/grid_noro_m.F90	(revision 3883)
+++ dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/grid_noro_m.F90	(revision 3883)
@@ -0,0 +1,353 @@
+MODULE grid_noro_m
+!
+!*******************************************************************************
+
+  PRIVATE
+  PUBLIC :: grid_noro
+
+
+CONTAINS
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE grid_noro(xd,yd,zd,x,y,zphi,zmea,zstd,zsig,zgam,zthe,zpic,zval,mask)
+!
+!-------------------------------------------------------------------------------
+! Author: F. Lott (see also Z.X. Li, A. Harzallah et L. Fairhead)
+!-------------------------------------------------------------------------------
+! Purpose: Compute the Parameters of the SSO scheme as described in LOTT &MILLER
+!         (1997) and LOTT(1999).
+!-------------------------------------------------------------------------------
+! Comments:
+!  * Target points are on a rectangular grid:
+!      iim+1 latitudes including North and South Poles;
+!      jjm+1 longitudes, with periodicity jjm+1=1.
+!  * At the poles, the fields value is repeated jjm+1 time.
+!  * The parameters a,b,c,d represent the limits of the target gridpoint region.
+!    The means over this region are calculated from USN data, ponderated by a
+!    weight proportional to the surface occupated by the data inside the model
+!    gridpoint area. In most circumstances, this weight is the ratio between the
+!    surfaces of the USN gridpoint area and the model gridpoint area. 
+!
+!           (c)
+!        ----d-----
+!        | . . . .|
+!        |        |
+!     (b)a . * . .b(a)
+!        |        |
+!        | . . . .|
+!        ----c-----
+!           (d)
+!  * Hard-coded US Navy dataset dimensions (imdp=2160 ; jmdp=1080) have been
+!    removed (ALLOCATABLE used).
+!  * iext (currently 10% of imdp) represents the margin to ensure output cells
+!    on the edge are contained in input cells.
+!===============================================================================
+  USE assert_eq_m, ONLY: assert_eq
+  USE print_control_mod, ONLY: lunout
+  IMPLICIT NONE 
+  REAL, PARAMETER :: epsfra = 1.e-5
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)  :: xd(:), yd(:)  !--- INPUT  COORDINATES     (imdp) (jmdp)
+  REAL, INTENT(IN)  :: zd(:,:)       !--- INPUT  FIELD           (imdp,jmdp)
+  REAL, INTENT(IN)  :: x(:), y(:)    !--- OUTPUT COORDINATES     (imar+1) (jmar)
+  REAL, INTENT(OUT) :: zphi(:,:)     !--- GEOPOTENTIAL           (imar+1,jmar)
+  REAL, INTENT(OUT) :: zmea(:,:)     !--- MEAN OROGRAPHY         (imar+1,jmar)
+  REAL, INTENT(OUT) :: zstd(:,:)     !--- STANDARD DEVIATION     (imar+1,jmar)
+  REAL, INTENT(OUT) :: zsig(:,:)     !--- SLOPE                  (imar+1,jmar)
+  REAL, INTENT(OUT) :: zgam(:,:)     !--- ANISOTROPY             (imar+1,jmar)
+  REAL, INTENT(OUT) :: zthe(:,:)     !--- SMALL AXIS ORIENTATION (imar+1,jmar)
+  REAL, INTENT(OUT) :: zpic(:,:)     !--- MAXIMUM ALTITITUDE     (imar+1,jmar)
+  REAL, INTENT(OUT) :: zval(:,:)     !--- MINIMUM ALTITITUDE     (imar+1,jmar)
+  REAL, INTENT(OUT) :: mask(:,:)     !--- MASK                   (imar+1,jmar)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="grid_noro"
+  REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
+  REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
+! CORRELATIONS OF OROGRAPHY GRADIENT              ! dim (imar+1,jmar)
+  REAL, ALLOCATABLE :: ztz(:,:), zxtzx(:,:), zytzy(:,:), zxtzy(:,:), weight(:,:)
+! CORRELATIONS OF USN OROGRAPHY GRADIENTS         ! dim (imar+2*iext,jmdp+2)
+  REAL, ALLOCATABLE :: zxtzxusn(:,:), zytzyusn(:,:), zxtzyusn(:,:)
+  REAL, ALLOCATABLE :: mask_tmp(:,:), zmea0(:,:)  ! dim (imar+1,jmar)
+  REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
+  REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
+  REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
+  LOGICAL :: masque_lu
+  INTEGER :: i, ii, imdp, imar, iext
+  INTEGER :: j, jj, jmdp, jmar, nn
+  REAL    :: xpi, zdeltax, zlenx, weighx, xincr,  zmeanor0
+  REAL    :: rad, zdeltay, zleny, weighy, masque, zmeasud0
+  REAL    :: zbordnor, zmeanor, zstdnor, zsignor, zweinor, zpicnor, zvalnor
+  REAL    :: zbordsud, zmeasud, zstdsud, zsigsud, zweisud, zpicsud, zvalsud
+  REAL    :: zbordest, zbordoue, xk, xl, xm, xp, xq, xw
+!-------------------------------------------------------------------------------
+  imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
+  jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
+  imar=assert_eq([SIZE(x),SIZE(zphi,1),SIZE(zmea,1),SIZE(zstd,1),SIZE(zsig,1), &
+                          SIZE(zgam,1),SIZE(zthe,1),SIZE(zpic,1),SIZE(zval,1), &
+                          SIZE(mask,1)],TRIM(modname)//" imar")-1
+  jmar=assert_eq([SIZE(y),SIZE(zphi,2),SIZE(zmea,2),SIZE(zstd,2),SIZE(zsig,2), &
+                          SIZE(zgam,2),SIZE(zthe,2),SIZE(zpic,2),SIZE(zval,2), &
+                          SIZE(mask,2)],TRIM(modname)//" jmar")
+!  IF(imar/=iim)   CALL abort_physic(TRIM(modname),'imar/=iim'  ,1)
+!  IF(jmar/=jjm+1) CALL abort_physic(TRIM(modname),'jmar/=jjm+1',1)
+  iext=imdp/10                                !--- OK up to 36 degrees cell
+  xpi = ACOS(-1.)
+  rad = 6371229.
+  zdeltay=2.*xpi/REAL(jmdp)*rad
+  WRITE(lunout,*)"*** Orography parameters at sub-cell scale ***"
+
+!--- ARE WE USING A READ MASK ?
+  masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
+  WRITE(lunout,*)'Masque lu: ',masque_lu
+
+!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
+  ALLOCATE(xusn(imdp+2*iext))
+  xusn(1     +iext:imdp  +iext)=xd(:)
+  xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
+  xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
+
+  ALLOCATE(yusn(jmdp+2))
+  yusn(1       )=yd(1)   +(yd(1)   -yd(2))
+  yusn(2:jmdp+1)=yd(:)
+  yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
+
+  ALLOCATE(zusn(imdp+2*iext,jmdp+2))
+  zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
+  zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
+  zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
+  zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
+  zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
+  zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
+  zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
+
+!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
+  ALLOCATE(a(imar+1),b(imar+1))
+  b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
+  b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
+  a(1)=x(1)-(x(2)-x(1))/2.0
+  a(2:imar+1)= b(1:imar)
+
+  ALLOCATE(c(jmar),d(jmar))
+  d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
+  d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
+  c(1)=y(1)-(y(2)-y(1))/2.0
+  c(2:jmar)=d(1:jmar-1)
+
+!--- INITIALIZATIONS:
+  ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
+  ALLOCATE(zxtzx (imar+1,jmar)); zxtzx (:,:)= 0.0
+  ALLOCATE(zytzy (imar+1,jmar)); zytzy (:,:)= 0.0
+  ALLOCATE(zxtzy (imar+1,jmar)); zxtzy (:,:)= 0.0
+  ALLOCATE(ztz   (imar+1,jmar)); ztz   (:,:)= 0.0
+  zmea(:,:)= 0.0
+  zpic(:,:)=-1.E+10
+  zval(:,:)= 1.E+10
+
+!--- COMPUTE SLOPES CORRELATIONS ON USN GRID
+! CORRELATIONS OF USN OROGRAPHY GRADIENTS  ! dim (imdp+2*iext,jmdp+2)
+  ALLOCATE(zytzyusn(imdp+2*iext,jmdp+2)); zytzyusn(:,:)=0.0
+  ALLOCATE(zxtzxusn(imdp+2*iext,jmdp+2)); zxtzxusn(:,:)=0.0
+  ALLOCATE(zxtzyusn(imdp+2*iext,jmdp+2)); zxtzyusn(:,:)=0.0
+  DO j = 2, jmdp+1
+    zdeltax=zdeltay*cos(yusn(j))
+    DO i = 2, imdp+2*iext-1
+      zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))**2/zdeltay**2
+      zxtzxusn(i,j)=(zusn(i+1,j)-zusn(i-1,j))**2/zdeltax**2
+      zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))   /zdeltay  &
+     &             *(zusn(i+1,j)-zusn(i-1,j))   /zdeltax
+    END DO
+  END DO
+
+!--- SUMMATION OVER GRIDPOINT AREA
+  zleny=xpi/REAL(jmdp)*rad
+  xincr=xpi/REAL(jmdp)/2.
+  ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
+  ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
+  DO ii = 1, imar+1
+    DO jj = 1, jmar
+      DO j = 2,jmdp+1 
+        zlenx  =zleny  *COS(yusn(j))
+        zdeltax=zdeltay*COS(yusn(j))
+        zbordnor=(xincr+c(jj)-yusn(j))*rad
+        zbordsud=(xincr-d(jj)+yusn(j))*rad
+        weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
+        IF(weighy==0.) CYCLE
+        DO i = 2, imdp+2*iext-1
+          zbordest=(xusn(i)-a(ii)+xincr)*rad*COS(yusn(j))
+          zbordoue=(b(ii)+xincr-xusn(i))*rad*COS(yusn(j))
+          weighx=AMAX1(0.,AMIN1(zbordest,zbordoue,zlenx))
+          IF(weighx==0.) CYCLE
+          num_tot(ii,jj)=num_tot(ii,jj)+1.0
+          IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
+          weight(ii,jj)=weight(ii,jj)+weighx*weighy
+          zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)*weighx*weighy
+          zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy
+          zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy
+          ztz  (ii,jj)=  ztz(ii,jj)+zusn(i,j)*zusn(i,j)*weighx*weighy
+          zmea (ii,jj)= zmea(ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
+          zpic (ii,jj)=AMAX1(zpic(ii,jj),zusn(i,j))         !--- PEAKS
+          zval (ii,jj)=AMIN1(zval(ii,jj),zusn(i,j))         !--- VALLEYS
+        END DO
+      END DO
+    END DO
+  END DO
+
+!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
+  IF(.NOT.masque_lu) THEN
+    WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
+  END IF
+  nn=COUNT(weight(:,1:jmar-1)==0.0)
+  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
+  WHERE(weight(:,:)/=0.0)
+    zmea (:,:)=zmea (:,:)/weight(:,:)
+    zxtzx(:,:)=zxtzx(:,:)/weight(:,:)
+    zytzy(:,:)=zytzy(:,:)/weight(:,:)
+    zxtzy(:,:)=zxtzy(:,:)/weight(:,:)
+    ztz  (:,:)=ztz  (:,:)/weight(:,:)
+    zstd (:,:)=ztz  (:,:)-zmea(:,:)**2
+  END WHERE
+  WHERE(zstd(:,:)<0) zstd(:,:)=0.
+  zstd (:,:)=SQRT(zstd(:,:))
+
+!--- CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:
+  zxtzx(:,   1)=zxtzx(:,     2)
+  zxtzx(:,jmar)=zxtzx(:,jmar-1)
+  zxtzy(:,   1)=zxtzy(:,     2)
+  zxtzy(:,jmar)=zxtzy(:,jmar-1)
+  zytzy(:,   1)=zytzy(:,     2)
+  zytzy(:,jmar)=zytzy(:,jmar-1)
+
+!=== FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.
+!--- FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
+!-------------------------------------------------------------------------------
+  ALLOCATE(zmea0(imar+1,jmar))
+  zmea0(:,:)=zmea(:,:)                           ! GK211005 (CG) UNSMOOTHED TOPO
+  CALL MVA9(zmea);  CALL MVA9(zstd);  CALL MVA9(zpic);  CALL MVA9(zval)
+  CALL MVA9(zxtzx); CALL MVA9(zxtzy); CALL MVA9(zytzy)
+
+!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD. (SURFACE PARAMS MEANINGLESS)
+  ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
+  WHERE(mask>=0.1) mask_tmp = 1.
+  WHERE(weight(:,:)/=0.0)
+!   zphi (:,:)= mask_tmp(:,:)*zmea (:,:) ! GK211005 (CG) not necessarly smoothed
+    zphi (:,:)= mask_tmp(:,:)*zmea0(:,:)
+    zmea0(:,:)= mask_tmp(:,:)*zmea0(:,:)
+    zmea (:,:)= mask_tmp(:,:)*zmea (:,:)
+    zpic (:,:)= mask_tmp(:,:)*zpic (:,:)
+    zval (:,:)= mask_tmp(:,:)*zval (:,:)
+    zstd (:,:)= mask_tmp(:,:)*zstd (:,:)
+  END WHERE
+  DO ii = 1, imar
+    DO jj = 1, jmar
+      IF (weight(ii,jj)/=0.0) THEN
+      !--- Coefficients K, L et M:
+        xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2.
+        xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2.
+        xm=zxtzy(ii,jj)
+        xp=xk-SQRT(xl**2+xm**2)
+        xq=xk+SQRT(xl**2+xm**2)
+        xw=1.e-8
+        IF(xp<=xw) xp=0.
+        IF(xq<=xw) xq=xw
+        IF(ABS(xm)<=xw) xm=xw*SIGN(1.,xm)
+      !--- SLOPE
+        zsig(ii,jj)=SQRT(xq)*mask_tmp(ii,jj)
+      !---ISOTROPY
+        zgam(ii,jj)=xp/xq*mask_tmp(ii,jj)
+      !--- THETA ANGLE
+        zthe(ii,jj)=57.29577951*ATAN2(xm,xl)/2.*mask_tmp(ii,jj)
+      END IF
+    END DO
+  END DO
+  WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
+  WRITE(lunout,*)'  ST. DEV.:' ,MAXVAL(zstd)
+  WRITE(lunout,*)'  PENTE:'    ,MAXVAL(zsig)
+  WRITE(lunout,*)'  ANISOTROP:',MAXVAL(zgam)
+  WRITE(lunout,*)'  ANGLE:'    ,MINVAL(zthe),MAXVAL(zthe)
+  WRITE(lunout,*)'  pic:'      ,MAXVAL(zpic)
+  WRITE(lunout,*)'  val:'      ,MAXVAL(zval)
+      
+!--- Values at poles
+  zmea0(imar+1,:)=zmea0(1,:)
+  zmea (imar+1,:)=zmea (1,:)
+  zphi (imar+1,:)=zphi (1,:)
+  zpic (imar+1,:)=zpic (1,:)
+  zval (imar+1,:)=zval (1,:)
+  zstd (imar+1,:)=zstd (1,:)
+  zsig (imar+1,:)=zsig (1,:)
+  zgam (imar+1,:)=zgam (1,:)
+  zthe (imar+1,:)=zthe (1,:)
+
+  zweinor =SUM(weight(1:imar,   1),DIM=1)
+  zweisud =SUM(weight(1:imar,jmar),DIM=1)
+  zmeanor0=SUM(weight(1:imar,   1)*zmea0(1:imar,   1),DIM=1)
+  zmeasud0=SUM(weight(1:imar,jmar)*zmea0(1:imar,jmar),DIM=1)
+  zmeanor =SUM(weight(1:imar,   1)*zmea (1:imar,   1),DIM=1)
+  zmeasud =SUM(weight(1:imar,jmar)*zmea (1:imar,jmar),DIM=1)
+  zstdnor =SUM(weight(1:imar,   1)*zstd (1:imar,   1),DIM=1)
+  zstdsud =SUM(weight(1:imar,jmar)*zstd (1:imar,jmar),DIM=1)
+  zsignor =SUM(weight(1:imar,   1)*zsig (1:imar,   1),DIM=1)
+  zsigsud =SUM(weight(1:imar,jmar)*zsig (1:imar,jmar),DIM=1)
+  zpicnor =SUM(weight(1:imar,   1)*zpic (1:imar,   1),DIM=1)
+  zpicsud =SUM(weight(1:imar,jmar)*zpic (1:imar,jmar),DIM=1)
+  zvalnor =SUM(weight(1:imar,   1)*zval (1:imar,   1),DIM=1)
+  zvalsud =SUM(weight(1:imar,jmar)*zval (1:imar,jmar),DIM=1)
+
+  zmea(:,1)=zmeanor /zweinor; zmea(:,jmar)=zmeasud /zweisud
+!  zphi(:,1)=zmeanor0/zweinor; zphi(:,jmar)=zmeasud0/zweisud   TO COMMIT
+  zphi(:,1)=zmeanor /zweinor; zphi(:,jmar)=zmeasud /zweisud
+  zpic(:,1)=zpicnor /zweinor; zpic(:,jmar)=zpicsud /zweisud
+  zval(:,1)=zvalnor /zweinor; zval(:,jmar)=zvalsud /zweisud
+  zstd(:,1)=zstdnor /zweinor; zstd(:,jmar)=zstdsud /zweisud
+  zsig(:,1)=zsignor /zweinor; zsig(:,jmar)=zsigsud /zweisud
+  zgam(:,1)=1.;               zgam(:,jmar)=1.
+  zthe(:,1)=0.;               zthe(:,jmar)=0.
+
+END SUBROUTINE grid_noro
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE MVA9(x)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+! MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(INOUT) :: x(:,:)
+!-------------------------------------------------------------------------------
+! Local variables:
+  REAL    :: xf(SIZE(x,DIM=1),SIZE(x,DIM=2)), WEIGHTpb(-1:1,-1:1)
+  INTEGER :: i, j, imar, jmar
+!-------------------------------------------------------------------------------
+  WEIGHTpb=RESHAPE([((1./REAL((1+i**2)*(1+j**2)),i=-1,1),j=-1,1)],SHAPE=[3,3])
+  WEIGHTpb=WEIGHTpb/SUM(WEIGHTpb)
+  imar=SIZE(X,DIM=1); jmar=SIZE(X,DIM=2)
+  DO j=2,jmar-1
+    DO i=2,imar-1
+      xf(i,j)=SUM(x(i-1:i+1,j-1:j+1)*WEIGHTpb(:,:))
+    END DO
+  END DO
+  DO j=2,jmar-1
+    xf(1,j)=SUM(x(imar-1,j-1:j+1)*WEIGHTpb(-1,:))
+    xf(1,j)=xf(1,j)+SUM(x(1:2,j-1:j+1)*WEIGHTpb(0:1,-1:1))
+    xf(imar,j)=xf(1,j)
+  END DO
+  xf(:,   1)=xf(:,     2)
+  xf(:,jmar)=xf(:,jmar-1)
+  x(:,:)=xf(:,:)
+
+END SUBROUTINE MVA9
+!
+!-------------------------------------------------------------------------------
+
+
+END MODULE grid_noro_m
+
