Index: LMDZ4/trunk/libf/dyn3dpar/startvar.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3dpar/startvar.F90	(revision 1323)
+++ LMDZ4/trunk/libf/dyn3dpar/startvar.F90	(revision 1425)
@@ -643,8 +643,10 @@
 !-------------------------------------------------------------------------------
 !
-SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2,lon_in2,&
-                          lat_in2, pls_in, var3d, ibar)
-!
-!-------------------------------------------------------------------------------
+SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2, &
+     lon_in2, lat_in2, pls_in, var3d, ibar)
+
+  use pchsp_95_m, only: pchsp_95
+  use pchfe_95_m, only: pchfe_95
+
 ! Arguments:
   CHARACTER(LEN=*),             INTENT(IN)    :: varname
@@ -655,73 +657,71 @@
   REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in2
   REAL, DIMENSION(jml2),        INTENT(IN)    :: lat_in2
-  REAL, DIMENSION(iml,jml,lml), INTENT(IN)    :: pls_in
-  REAL, DIMENSION(iml,jml,lml), INTENT(OUT)   :: var3d
+  REAL, DIMENSION(iml, jml, lml), INTENT(IN)    :: pls_in
+  REAL, DIMENSION(iml, jml, lml), INTENT(OUT)   :: var3d
   LOGICAL,                      INTENT(IN)    :: ibar
-!-------------------------------------------------------------------------------
+!----------------------------------------------------------------------------
 ! Local variables:
 #include "iniprint.h"
-  LOGICAL               :: check=.TRUE.
-  REAL                  :: bx, by, chmin, chmax
-  INTEGER               :: ii, ij, il
-  REAL,    DIMENSION(:,:,:), ALLOCATABLE :: var_tmp3d
+  LOGICAL:: check=.TRUE., skip
+  REAL                  chmin, chmax
+  INTEGER ii, ij, il, ierr
+  integer n_extrap ! number of extrapolated points
+  REAL, DIMENSION(iml, jml, llm_dyn):: var_tmp3d
   REAL,    DIMENSION(:),     ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
-  REAL,    DIMENSION(:),     ALLOCATABLE :: lev_dyn, ax, ay, yder
-  INTEGER, DIMENSION(:),     ALLOCATABLE :: lind
-!-------------------------------------------------------------------------------
-  IF(check) WRITE(lunout,*)'Going into flinget to extract the 3D  field.'
-  IF(check) WRITE(lunout,*)fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn
-  IF(check) WRITE(lunout,*)'Allocating space for interpolation',iml,jml,llm_dyn
-
-  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn,jml_dyn,llm_dyn))
-  CALL flinget(fid_dyn,varname,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
+  REAL, DIMENSION(llm_dyn):: lev_dyn, ax, ay, yder
+
+!---------------------------------------------------------------------------
+  IF(check) WRITE(lunout, *)'Going into flinget to extract the 3D  field.'
+  IF(check) WRITE(lunout, *) fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, &
+       ttm_dyn
+  IF(check) WRITE(lunout, *) 'Allocating space for interpolation', iml, jml, &
+       llm_dyn
+
+  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
+  CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, &
+       var_ana3d)
 
 !--- 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(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 (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,       &
+  ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn))
+  CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,      &
                    levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
-  DEALLOCATE(lon_ini,lat_ini)
+  DEALLOCATE(lon_ini, lat_ini)
 
 !--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
-  ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
-  DO il=1,llm_dyn
-    CALL interp_startvar(varname, ibar, il==1,                                 &
-      iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana3d(:,:,il), iml, jml, jml2,   &
-      lon_in,  lat_in,  lon_in2, lat_in2, var_tmp3d(:,:,il))
+  DO il=1, llm_dyn
+    CALL interp_startvar(varname, ibar, il==1, iml_dyn, jml_dyn, lon_rad, &
+         lat_rad, var_ana3d(:, :, il), iml, jml, jml2, lon_in, lat_in, &
+         lon_in2, lat_in2, var_tmp3d(:, :, il))
   END DO
-  DEALLOCATE(lon_rad,lat_rad)
-
-  ALLOCATE(lind(llm_dyn))
-  DO il=1,llm_dyn
-    lind(il) = llm_dyn-il+1
+  DEALLOCATE(lon_rad, lat_rad)
+
+!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
+  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) stop 1
+      n_extrap = n_extrap + ierr
+    END DO
   END DO
-
-!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
-  ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
-  DO ij=1,jml
-    DO ii=1,iml-1
-      ax(:)=lev_dyn(lind(:)) 
-      ay(:)=var_tmp3d(ii,ij,lind(:))
-      CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
-      DO il=1,lml
-        bx=pls_in(ii,ij,il)
-        CALL SPLINT(ax, ay, yder, llm_dyn, bx, by)
-        var3d(ii,ij,il) = by
-      END DO
-    END DO
-    var3d(iml,ij,:) = var3d(1,ij,:) 
+  if (n_extrap /= 0) then
+     print *, "start_inter_3d pchfe_95: n_extrap = ", n_extrap
+  end if
+  var3d(iml, :, :) = var3d(1, :, :) 
+
+  DO il=1, lml
+    CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
+    WRITE(lunout, *)' '//TRIM(varname)//'  min max l ', il, chmin, chmax
   END DO
-  DEALLOCATE(var_tmp3d,lev_dyn,ax,ay,yder,lind)
-
-  DO il=1,lml
-    CALL minmax(iml*jml,var3d(1,1,il),chmin,chmax)
-    WRITE(lunout,*)' '//TRIM(varname)//'  min max l ',il,chmin,chmax
-  END DO
-
-  RETURN
 
 END SUBROUTINE start_inter_3d
