Index: LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90	(revision 2542)
+++ LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90	(revision 2594)
@@ -12,12 +12,7 @@
                      prad,pg,pr,pcpp,iflag_phys)
   USE dimphy, ONLY: init_dimphy
-  USE mod_grid_phy_lmdz, ONLY: klon_glo,  & ! number of atmospheric columns (on full grid)
-                               regular_lonlat, &  ! regular longitude-latitude grid type
-                               nbp_lon, nbp_lat, nbp_lev
-  USE mod_phys_lmdz_para, ONLY: klon_omp, & ! number of columns (on local omp grid)
-                                klon_omp_begin, & ! start index of local omp subgrid
-                                klon_omp_end, & ! end index of local omp subgrid
-                                klon_mpi_begin ! start indes of columns (on local mpi grid)
-  USE geometry_mod, ONLY : init_geometry
+  USE inigeomphy_mod, ONLY: inigeomphy
+  USE mod_grid_phy_lmdz, ONLY: klon_glo ! number of atmospheric columns (on full grid)
+  USE mod_phys_lmdz_para, ONLY: klon_omp ! number of columns (on local omp grid)
   USE vertical_layers_mod, ONLY : init_vertical_layers
   USE infotrac, ONLY: nqtot,nqo,nbtr,tname,ttext,type_trac,&
@@ -39,10 +34,4 @@
   USE phystokenc_mod, ONLY: init_phystokenc
   USE phyaqua_mod, ONLY: iniaqua
-  USE physics_distribution_mod, ONLY : init_physics_distribution
-  USE regular_lonlat_mod, ONLY : init_regular_lonlat, &
-                                 east, west, north, south, &
-                                 north_east, north_west, &
-                                 south_west, south_east
-  USE mod_interface_dyn_phys, ONLY :  init_interface_dyn_phys
 #ifdef INCA
   USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic
@@ -92,28 +81,5 @@
   CHARACTER (LEN=20) :: modname = 'iniphysiq'
   CHARACTER (LEN=80) :: abort_message
-  REAL :: total_area_phy, total_area_dyn
 
-  ! boundaries, on global grid
-  REAL,ALLOCATABLE :: boundslon_reg(:,:)
-  REAL,ALLOCATABLE :: boundslat_reg(:,:)
-
-  ! global array, on full physics grid:
-  REAL,ALLOCATABLE :: latfi_glo(:)
-  REAL,ALLOCATABLE :: lonfi_glo(:)
-  REAL,ALLOCATABLE :: cufi_glo(:)
-  REAL,ALLOCATABLE :: cvfi_glo(:)
-  REAL,ALLOCATABLE :: airefi_glo(:)
-  REAL,ALLOCATABLE :: boundslonfi_glo(:,:)
-  REAL,ALLOCATABLE :: boundslatfi_glo(:,:)
-
-  ! local arrays, on given MPI/OpenMP domain:
-  REAL,ALLOCATABLE,SAVE :: latfi(:)
-  REAL,ALLOCATABLE,SAVE :: lonfi(:)
-  REAL,ALLOCATABLE,SAVE :: cufi(:)
-  REAL,ALLOCATABLE,SAVE :: cvfi(:)
-  REAL,ALLOCATABLE,SAVE :: airefi(:)
-  REAL,ALLOCATABLE,SAVE :: boundslonfi(:,:)
-  REAL,ALLOCATABLE,SAVE :: boundslatfi(:,:)
-!$OMP THREADPRIVATE (latfi,lonfi,cufi,cvfi,airefi,boundslonfi,boundslatfi)
 
 #ifndef CPP_PARA
@@ -123,145 +89,15 @@
 #endif
 
-  ! Initialize Physics distibution and parameters and interface with dynamics
-  IF (ii*jj>1) THEN ! general 3D case
-    CALL init_physics_distribution(regular_lonlat,4, &
-                                 nbp,ii,jj+1,nlayer,communicator)
-  ELSE ! For 1D model
-    CALL init_physics_distribution(regular_lonlat,4, &
-                                 1,1,1,nlayer,communicator)  
-  ENDIF
-  CALL init_interface_dyn_phys
+  ! --> initialize physics distribution, global fields and geometry
+  ! (i.e. things in phy_common or dynphy_lonlat)
+  CALL inigeomphy(ii,jj,nlayer, &
+               nbp, communicator, &
+               rlatu,rlatv, &
+               rlonu,rlonv, &
+               aire,cu,cv)
+
+  ! --> now initialize things specific to the phylmd physics package
   
-  ! init regular global longitude-latitude grid points and boundaries
-  ALLOCATE(boundslon_reg(ii,2))
-  ALLOCATE(boundslat_reg(jj+1,2))
-  
-  DO i=1,ii
-   boundslon_reg(i,east)=rlonu(i) 
-   boundslon_reg(i,west)=rlonu(i+1) 
-  ENDDO
-
-  boundslat_reg(1,north)= PI/2 
-  boundslat_reg(1,south)= rlatv(1)
-  DO j=2,jj
-   boundslat_reg(j,north)=rlatv(j-1) 
-   boundslat_reg(j,south)=rlatv(j) 
-  ENDDO
-  boundslat_reg(jj+1,north)= rlatv(jj) 
-  boundslat_reg(jj+1,south)= -PI/2
-
-  ! Write values in module regular_lonlat_mod
-  CALL init_regular_lonlat(ii,jj+1, rlonv(1:ii), rlatu, &
-                           boundslon_reg, boundslat_reg)
-
-  ! Generate global arrays on full physics grid
-  ALLOCATE(latfi_glo(klon_glo),lonfi_glo(klon_glo))
-  ALLOCATE(cufi_glo(klon_glo),cvfi_glo(klon_glo))
-  ALLOCATE(airefi_glo(klon_glo))
-  ALLOCATE(boundslonfi_glo(klon_glo,4))
-  ALLOCATE(boundslatfi_glo(klon_glo,4))
-  
-
-  IF (klon_glo>1) THEN ! general case
-    ! North pole
-    latfi_glo(1)=rlatu(1)
-    lonfi_glo(1)=0.
-    cufi_glo(1) = cu(1)
-    cvfi_glo(1) = cv(1)
-    boundslonfi_glo(1,north_east)=0
-    boundslatfi_glo(1,north_east)=PI/2
-    boundslonfi_glo(1,north_west)=2*PI
-    boundslatfi_glo(1,north_west)=PI/2
-    boundslonfi_glo(1,south_west)=2*PI
-    boundslatfi_glo(1,south_west)=rlatv(1)
-    boundslonfi_glo(1,south_east)=0
-    boundslatfi_glo(1,south_east)=rlatv(1)
-    DO j=2,jj
-      DO i=1,ii
-        k=(j-2)*ii+1+i
-        latfi_glo(k)= rlatu(j)
-        lonfi_glo(k)= rlonv(i)
-        cufi_glo(k) = cu((j-1)*(ii+1)+i)
-        cvfi_glo(k) = cv((j-1)*(ii+1)+i)
-        boundslonfi_glo(k,north_east)=rlonu(i)
-        boundslatfi_glo(k,north_east)=rlatv(j-1)
-        boundslonfi_glo(k,north_west)=rlonu(i+1)
-        boundslatfi_glo(k,north_west)=rlatv(j-1)
-        boundslonfi_glo(k,south_west)=rlonu(i+1)
-        boundslatfi_glo(k,south_west)=rlatv(j)
-        boundslonfi_glo(k,south_east)=rlonu(i)
-        boundslatfi_glo(k,south_east)=rlatv(j)
-      ENDDO
-    ENDDO
-    ! South pole
-    latfi_glo(klon_glo)= rlatu(jj+1)
-    lonfi_glo(klon_glo)= 0.
-    cufi_glo(klon_glo) = cu((ii+1)*jj+1)
-    cvfi_glo(klon_glo) = cv((ii+1)*jj-ii)
-    boundslonfi_glo(klon_glo,north_east)= 0
-    boundslatfi_glo(klon_glo,north_east)= rlatv(jj)
-    boundslonfi_glo(klon_glo,north_west)= 2*PI
-    boundslatfi_glo(klon_glo,north_west)= rlatv(jj)
-    boundslonfi_glo(klon_glo,south_west)= 2*PI
-    boundslatfi_glo(klon_glo,south_west)= -PI/2
-    boundslonfi_glo(klon_glo,south_east)= 0
-    boundslatfi_glo(klon_glo,south_east)= -Pi/2
-
-    ! build airefi_glo(), mesh area on physics grid
-    CALL gr_dyn_fi(1,ii+1,jj+1,klon_glo,aire,airefi_glo)
-    ! Poles are single points on physics grid
-    airefi_glo(1)=sum(aire(1:ii,1))
-    airefi_glo(klon_glo)=sum(aire(1:ii,jj+1))
-
-    ! Sanity check: do total planet area match between physics and dynamics?
-    total_area_dyn=sum(aire(1:ii,1:jj+1))
-    total_area_phy=sum(airefi_glo(1:klon_glo))
-    IF (total_area_dyn/=total_area_phy) THEN
-      WRITE (lunout, *) 'iniphysiq: planet total surface discrepancy !!!'
-      WRITE (lunout, *) '     in the dynamics total_area_dyn=', total_area_dyn
-      WRITE (lunout, *) '  but in the physics total_area_phy=', total_area_phy
-      IF (abs(total_area_dyn-total_area_phy)>0.00001*total_area_dyn) THEN
-        ! stop here if the relative difference is more than 0.001%
-        abort_message = 'planet total surface discrepancy'
-        CALL abort_gcm(modname, abort_message, 1)
-      ENDIF
-    ENDIF
-  ELSE ! klon_glo==1, running the 1D model
-    ! just copy over input values
-    latfi_glo(1)=rlatu(1)
-    lonfi_glo(1)=rlonv(1)
-    cufi_glo(1)=cu(1)
-    cvfi_glo(1)=cv(1)
-    airefi_glo(1)=aire(1,1)
-    boundslonfi_glo(1,north_east)=rlonu(1)
-    boundslatfi_glo(1,north_east)=PI/2
-    boundslonfi_glo(1,north_west)=rlonu(2)
-    boundslatfi_glo(1,north_west)=PI/2
-    boundslonfi_glo(1,south_west)=rlonu(2)
-    boundslatfi_glo(1,south_west)=rlatv(1)
-    boundslonfi_glo(1,south_east)=rlonu(1)
-    boundslatfi_glo(1,south_east)=rlatv(1)
-  ENDIF ! of IF (klon_glo>1)
-
 !$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/)
-  ! Now generate local lon/lat/cu/cv/area/bounds arrays
-  ALLOCATE(latfi(klon_omp),lonfi(klon_omp),cufi(klon_omp),cvfi(klon_omp))
-  ALLOCATE(airefi(klon_omp))
-  ALLOCATE(boundslonfi(klon_omp,4))
-  ALLOCATE(boundslatfi(klon_omp,4))
-
-
-  offset = klon_mpi_begin - 1
-  airefi(1:klon_omp) = airefi_glo(offset+klon_omp_begin:offset+klon_omp_end)
-  cufi(1:klon_omp) = cufi_glo(offset+klon_omp_begin:offset+klon_omp_end)
-  cvfi(1:klon_omp) = cvfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
-  lonfi(1:klon_omp) = lonfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
-  latfi(1:klon_omp) = latfi_glo(offset+klon_omp_begin:offset+klon_omp_end)
-  boundslonfi(1:klon_omp,:) = boundslonfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
-  boundslatfi(1:klon_omp,:) = boundslatfi_glo(offset+klon_omp_begin:offset+klon_omp_end,:)
-
-  ! copy over local grid longitudes and latitudes
-  CALL init_geometry(klon_omp,lonfi,latfi,boundslonfi,boundslatfi, &
-                     airefi,cufi,cvfi)
 
   ! copy over preff , ap(), bp(), etc 
Index: LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/limit_netcdf.F90
===================================================================
--- LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/limit_netcdf.F90	(revision 2542)
+++ LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/limit_netcdf.F90	(revision 2594)
@@ -117,4 +117,7 @@
     ALLOCATE(pctsrf(klon,nbsrf))
     CALL start_init_subsurf(.FALSE.)
+  !--- TO MATCH EXACTLY WHAT WOULD BE DONE IN etat0phys_netcdf
+    WHERE(   masque(:,:)<EPSFRA) masque(:,:)=0.
+    WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
   END IF
 
@@ -336,4 +339,5 @@
   REAL, ALLOCATABLE :: champan(:,:,:)
 !--- input files
+  CHARACTER(LEN=20) :: fnam_m, fnam_p     ! previous/next files names
   CHARACTER(LEN=20) :: cal_in             ! calendar
   CHARACTER(LEN=20) :: unit_sic           ! attribute unit in sea-ice file
@@ -345,5 +349,5 @@
   LOGICAL           :: extrp              ! flag for extrapolation
   REAL              :: chmin, chmax
-  INTEGER ierr
+  INTEGER ierr, idx
   integer n_extrap ! number of extrapolated points
   logical skip
@@ -360,4 +364,5 @@
   END SELECT
   extrp=.FALSE.; IF(PRESENT(flag).AND.mode=='SST') extrp=flag
+  idx=INDEX(fnam,'.nc')-1
 
 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
@@ -393,5 +398,5 @@
 !--- 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(MAX(lmdep,14)))
+  ALLOCATE(timeyear(lmdep+2))
   CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
   cal_in=' '
@@ -412,19 +417,16 @@
   IF(lmdep==12) THEN
     timeyear=mid_month(anneeref, cal_in, ndays_in)
-    CALL msg(0,'Monthly periodic input file (perpetual run).')
-  ELSE IF(lmdep==14) THEN
-    timeyear=mid_month(anneeref, cal_in, ndays_in)
-    CALL msg(0,'Monthly 14-records input file (interannual run).')
+    CALL msg(0,'Monthly input file(s) for '//TRIM(title)//'.')
   ELSE IF(lmdep==ndays_in) THEN
-    timeyear=[(REAL(k)-0.5,k=1,ndays_in)]
+    timeyear=[(REAL(k)-0.5,k=0,ndays_in+1)]
     CALL msg(0,'Daily input file (no time interpolation).')
   ELSE
     WRITE(mess,'(a,i3,a,i3,a)')'Mismatching input file: found',lmdep,        &
-      ' records, 12/14/',ndays_in,' (periodic/interannual/daily) needed'
+      ' records, 12/',ndays_in,' (monthly/daily needed).'
     CALL abort_physic('mid_months',TRIM(mess),1)
   END IF
 
 !--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
-  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, MAX(lmdep,14)))
+  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep+2))
   IF(extrp) ALLOCATE(work(imdep, jmdep))
   CALL msg(5,'')
@@ -446,14 +448,51 @@
       WHERE(NINT(mask)/=1) champint=0.001
     END IF
-    !--- Special case for periodic input file: index shifted
-    ll = l; IF(lmdep==12) ll = l + 1
-    champtime(:, :, ll)=champint
+    champtime(:, :, l+1)=champint
   END DO
-  IF(lmdep==12) THEN
-    champtime(:,:, 1)=champtime(:,:,13)
-    champtime(:,:,14)=champtime(:,:, 2)
-  END IF
   CALL ncerr(NF90_CLOSE(ncid), fnam)
 
+!--- FIRST RECORD: LAST ONE OF PREVIOUS YEAR (CURRENT YEAR IF UNAVAILABLE)
+  fnam_m=fnam(1:idx)//'_m.nc'
+  IF(NF90_OPEN(fnam_m,NF90_NOWRITE,ncid)==NF90_NOERR) THEN
+    CALL msg(0,'Reading previous year file ("'//TRIM(fnam_m)//'") last record for '//TRIM(title))
+    CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_m)
+    CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam_m)
+    CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), len=l), fnam_m)
+    CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam_m)
+    CALL ncerr(NF90_CLOSE(ncid), fnam_m)
+    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(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(:, :, 1)=champint
+  ELSE
+    CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") last record for '//TRIM(title))
+    champtime(:, :, 1)=champtime(:, :, lmdep+1)
+  END IF
+
+!--- LAST RECORD: FIRST ONE OF NEXT YEAR (CURRENT YEAR IF UNAVAILABLE)
+  fnam_p=fnam(1:idx)//'_p.nc'
+  IF(NF90_OPEN(fnam_p,NF90_NOWRITE,ncid)==NF90_NOERR) THEN
+    CALL msg(0,'Reading previous year file ("'//TRIM(fnam_p)//'") first record for '//TRIM(title))
+    CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_p)
+    CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,1],[imdep,jmdep,1]),fnam_p)
+    CALL ncerr(NF90_CLOSE(ncid), fnam_p)
+    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(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(:, :, lmdep+2)=champint
+  ELSE
+    CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") first record for '//TRIM(title))
+    champtime(:, :, lmdep+2)=champtime(:, :, 2)
+  END IF
   DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
   IF(extrp) DEALLOCATE(work)
@@ -477,5 +516,5 @@
      END IF
   END IF
-  ALLOCATE(yder(MAX(lmdep,14)), champan(iip1, jjp1, ndays))
+  ALLOCATE(yder(lmdep+2), champan(iip1, jjp1, ndays))
   IF(lmdep==ndays_in) THEN
      champan(1:iim,:,:)=champtime
@@ -546,4 +585,5 @@
 !
 !-------------------------------------------------------------------------------
+  USE grid_noro_m, ONLY: grid_noro0
   IMPLICIT NONE
 !===============================================================================
@@ -599,147 +639,4 @@
 
 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
 !
 !-------------------------------------------------------------------------------
