Index: LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90	(revision 1404)
+++ LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90	(revision 1425)
@@ -267,5 +267,4 @@
 #endif
 ! of #ifdef CPP_EARTH
-  STOP
 
 
@@ -281,5 +280,5 @@
 SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL)
 !
-!-------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
 ! Comments:
 !   There are two assumptions concerning the NetCDF files, that are satisfied
@@ -287,11 +286,15 @@
 !     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 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 dimphy, ONLY : klon
   USE phys_state_var_mod, ONLY : pctsrf
   USE control_mod
+  use pchsp_95_m, only: pchsp_95
+  use pchfe_95_m, only: pchfe_95
+  use arth_m, only: arth
+
   IMPLICIT NONE
 #include "dimensions.h"
@@ -300,5 +303,5 @@
 #include "indicesol.h"
 #include "iniprint.h"
-!-------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
 ! Arguments:
   CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
@@ -306,12 +309,12 @@
   LOGICAL,           INTENT(IN)     :: ibar     ! interp on pressure levels
   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
+  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
   LOGICAL, OPTIONAL, INTENT(IN)     :: lCPL     ! Coupled model flag (for ICE)
-!-------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
 ! Local variables:
 !--- NetCDF
-  INTEGER :: ierr, ncid, varid                  ! NetCDF identifiers
+  INTEGER :: ncid, varid                  ! NetCDF identifiers
   CHARACTER(LEN=30)               :: dnam       ! dimension name
   CHARACTER(LEN=80)               :: varname    ! NetCDF variable name
@@ -323,10 +326,9 @@
 !--- fields
   INTEGER :: imdep, jmdep, lmdep                ! dimensions of 'champ'
-  REAL, ALLOCATABLE, DIMENSION(:,:) :: champ      ! wanted field on initial grid
+  REAL, ALLOCATABLE, DIMENSION(:, :) :: champ   ! wanted field on initial grid
   REAL, ALLOCATABLE, DIMENSION(:) :: yder, timeyear
-  REAL,              DIMENSION(iim,jjp1) :: champint   ! interpolated field
-  REAL, ALLOCATABLE, DIMENSION(:,:,:)    :: champtime
-  REAL, ALLOCATABLE, DIMENSION(:,:,:)    :: champan
-  REAL                              :: time, by
+  REAL,              DIMENSION(iim, jjp1) :: champint   ! interpolated field
+  REAL, ALLOCATABLE, DIMENSION(:, :, :)    :: champtime
+  REAL, ALLOCATABLE, DIMENSION(:, :, :)    :: champan
 !--- input files
   CHARACTER(LEN=20)                 :: cal_in   ! calendar
@@ -334,10 +336,13 @@
 !--- misc
   INTEGER :: i, j, k, l                         ! loop counters
-  REAL, ALLOCATABLE, DIMENSION(:,:) :: work     ! used for extrapolation
+  REAL, ALLOCATABLE, DIMENSION(:, :) :: work     ! used for extrapolation
   CHARACTER(LEN=25)                 :: title    ! for messages
   LOGICAL                           :: extrp    ! flag for extrapolation
   REAL                              :: chmin, chmax
-!-------------------------------------------------------------------------------
-!---Variables depending on keyword 'mode' --------------------------------------
+  INTEGER ierr
+  integer n_extrap ! number of extrapolated points
+  logical skip
+!------------------------------------------------------------------------------
+!---Variables depending on keyword 'mode' -------------------------------------
   NULLIFY(champo)
   SELECT CASE(mode)
@@ -352,70 +357,75 @@
   END IF
 
-!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ------------------------------
-  ierr=NF90_OPEN(fnam,NF90_NOWRITE,ncid);                  CALL ncerr(ierr,fnam)
-  ierr=NF90_INQ_VARID(ncid,varname,varid);                 CALL ncerr(ierr,fnam)
-  ierr=NF90_INQUIRE_VARIABLE(ncid,varid,dimids=dids);      CALL ncerr(ierr,fnam)
+!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
+  ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid);             CALL ncerr(ierr, fnam)
+  ierr=NF90_INQ_VARID(ncid, varname, varid);            CALL ncerr(ierr, fnam)
+  ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam)
 
 !--- Longitude
-  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(1),name=dnam,len=imdep)
-  CALL ncerr(ierr,fnam); ALLOCATE(dlon_ini(imdep),dlon(imdep))
-  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
-  ierr=NF90_GET_VAR(ncid,varid,dlon_ini);                  CALL ncerr(ierr,fnam)
-  WRITE(lunout,*) 'variable ', dnam, 'dimension ', imdep
+  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep)
+  CALL ncerr(ierr, fnam); ALLOCATE(dlon_ini(imdep), dlon(imdep))
+  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
+  ierr=NF90_GET_VAR(ncid, varid, dlon_ini);              CALL ncerr(ierr, fnam)
+  WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
 
 !--- Latitude
-  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(2),name=dnam,len=jmdep)
-  CALL ncerr(ierr,fnam); ALLOCATE(dlat_ini(jmdep),dlat(jmdep))
-  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
-  ierr=NF90_GET_VAR(ncid,varid,dlat_ini);                  CALL ncerr(ierr,fnam)
-  WRITE(lunout,*) 'variable ', dnam, 'dimension ', jmdep
+  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep)
+  CALL ncerr(ierr, fnam); ALLOCATE(dlat_ini(jmdep), dlat(jmdep))
+  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
+  ierr=NF90_GET_VAR(ncid, varid, dlat_ini);              CALL ncerr(ierr, fnam)
+  WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
 
 !--- Time (variable is not needed - it is rebuilt - but calendar is)
-  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(3),name=dnam,len=lmdep)
-  CALL ncerr(ierr,fnam); ALLOCATE(timeyear(lmdep))
-  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
+  ierr=NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep)
+  CALL ncerr(ierr, fnam); ALLOCATE(timeyear(lmdep))
+  ierr=NF90_INQ_VARID(ncid, dnam, varid);                CALL ncerr(ierr, fnam)
   cal_in=' '
-  ierr=NF90_GET_ATT(ncid,varid,'calendar',cal_in)
+  ierr=NF90_GET_ATT(ncid, varid, 'calendar', cal_in)
   IF(ierr/=NF90_NOERR) THEN
     SELECT CASE(mode)
-      CASE('RUG','ALB'); cal_in='360d'
-      CASE('SIC','SST'); cal_in='gregorian'
+      CASE('RUG', 'ALB'); cal_in='360d'
+      CASE('SIC', 'SST'); cal_in='gregorian'
     END SELECT
-    WRITE(lunout,*)'ATTENTION: variable ''time'' sans attribut ''calendrier'' d&
-     &ans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
+    WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
+         // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
   END IF
-  WRITE(lunout,*) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', cal_in
-
-!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION ----------------------
+  WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
+       cal_in
+
+!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
   !--- Determining input file number of days, depending on calendar
-  ndays_in=year_len(anneeref,cal_in)
+  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,'(a,i3,a)')'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))
-
-  WRITE(lunout,*)
-  WRITE(lunout,'(a,i3,a)')'LECTURE ET INTERPOLATION HORIZ. DE ',lmdep,' CHAMPS.'
-  ierr=NF90_INQ_VARID(ncid,varname,varid);                 CALL ncerr(ierr,fnam)
-  DO l=1,lmdep
-    ierr=NF90_GET_VAR(ncid,varid,champ,(/1,1,l/),(/imdep,jmdep,1/))
-    CALL ncerr(ierr,fnam)
-    CALL conf_dat2d(title,imdep,jmdep,dlon_ini,dlat_ini,dlon,dlat,champ,ibar)
-
-    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
+  timeyear=mid_months(anneeref, cal_in, lmdep)
+  IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') '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))
+
+  WRITE(lunout, *)
+  WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, &
+       ' CHAMPS.'
+  ierr=NF90_INQ_VARID(ncid, varname, varid);             CALL ncerr(ierr, fnam)
+  DO l=1, lmdep
+    ierr=NF90_GET_VAR(ncid, varid, champ, (/1, 1, l/), (/imdep, jmdep, 1/))
+    CALL ncerr(ierr, fnam)
+    CALL conf_dat2d(title, imdep, jmdep, dlon_ini, dlat_ini, dlon, dlat, &
+         champ, ibar)
+
+    IF (extrp) CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, &
+         work)
 
     IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN
       IF(l==1) THEN
-        WRITE(lunout,*)                                                        &
+        WRITE(lunout, *)                                                      &
   '-------------------------------------------------------------------------'
-        WRITE(lunout,*)                                                        &
-  '$$$ Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
-        WRITE(lunout,*)                                                        &
+        WRITE(lunout, *)                                                     &
+  'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
+        WRITE(lunout, *)                                                      &
   '-------------------------------------------------------------------------'
       END IF
@@ -433,59 +443,64 @@
         CASE('SIC');       CALL sea_ice (imdep, jmdep, dlon, dlat, champ,    &
                                     iim, jjp1, rlonv, rlatu, champint)
-        CASE('SST','ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ,    &
+        CASE('SST', 'ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ,    &
                                     iim, jjp1, rlonv, rlatu, champint)
       END SELECT
     END IF
-    champtime(:,:,l)=champint
+    champtime(:, :, l)=champint
   END DO
-  ierr=NF90_CLOSE(ncid);                                   CALL ncerr(ierr,fnam)
-
-  DEALLOCATE(dlon_ini,dlat_ini,dlon,dlat,champ)
+  ierr=NF90_CLOSE(ncid);                                 CALL ncerr(ierr, fnam)
+
+  DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
   IF(extrp) DEALLOCATE(work)
 
-!--- TIME INTERPOLATION --------------------------------------------------------
-  WRITE(lunout,*)
-  WRITE(lunout,*)'INTERPOLATION TEMPORELLE.'
-  WRITE(lunout,"(2x,' Vecteur temps en entree: ',10f6.1)") timeyear
-  WRITE(lunout,"(2x,' Vecteur temps en sortie de 0 a ',i3)") ndays
-  ALLOCATE(yder(lmdep),champan(iip1,jjp1,ndays))
-  DO j=1,jjp1
-    DO i=1,iim
-      CALL spline(timeyear,champtime(i,j,:),lmdep,1.e30,1.e30,yder)
-      DO k=1,ndays
-        time=FLOAT((k-1)*ndays_in)/FLOAT(ndays)
-        CALL splint(timeyear,champtime(i,j,:),yder,lmdep,time,by)
-        champan(i,j,k) = by
-      END DO
+!--- TIME INTERPOLATION ------------------------------------------------------
+  WRITE(lunout, *)
+  WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
+  WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear
+  WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays
+  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
-  champan(iip1,:,:)=champan(1,:,:)
-  DEALLOCATE(yder,champtime,timeyear)
+  if (n_extrap /= 0) then
+     print *, "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)
-    WRITE(lunout,*)' '//TRIM(title)//' au temps 10 ',chmin,chmax,j
+  DO j=1, jjp1
+    CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
+    WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j
   END DO
 
-!--- SPECIAL FILTER FOR SST: SST>271.38 ----------------------------------------
+!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
   IF(mode=='SST') THEN
-    WRITE(lunout,*) 'Filtrage de la SST: SST >= 271.38'
+    WRITE(lunout, *) '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 ---------------------------------------
+!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
   IF(mode=='SIC') THEN
-    WRITE(lunout,*) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
-    IF(.NOT.lCPL) champan(:,:,:)=champan(:,:,:)/100.
-    champan(iip1,:,:)=champan(1,:,:)
+    WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
+    IF(.NOT.lCPL) champan(:, :, :)=champan(:, :, :)/100.
+    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))
+!--- 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)
