C
C $Header$
C
      SUBROUTINE limit_netcdf ( interbar, extrap, oldice )
c
      IMPLICIT none
c
c-------------------------------------------------------------
C Author : L. Fairhead
C Date   : 27/01/94
C Objet  : Construction des fichiers de conditions aux limites
C          pour le nouveau
C          modele a partir de fichiers de climatologie. Les deux
C          grilles doivent etre regulieres
c
c Modifie par z.x.li (le23mars1994)
c Modifie par L. Fairhead (fairhead@lmd.jussieu.fr) septembre 1999
c                         pour lecture netcdf dans LMDZ.3.3
c Modifie par P;Le Van  ,  juillet 2001
c-------------------------------------------------------------
c
#include "dimensions.h"
#include "paramet.h"
#include "control.h"
#include "logic.h"
#include "netcdf.inc"
#include "comvert.h"
#include "comgeom2.h"
#include "comconst.h"
c
c-----------------------------------------------------------------------
      LOGICAL interbar, extrap, oldice

      INTEGER KIDIA, KFDIA, KLON, KLEV
      PARAMETER (KIDIA=1,KFDIA=iim*(jjm-1)+2,
     .           KLON=KFDIA-KIDIA+1,KLEV=llm)
c-----------------------------------------------------------------------
      REAL phy_nat(klon,360), phy_nat0(klon)
      REAL phy_alb(klon,360)
      REAL phy_sst(klon,360)
      REAL phy_bil(klon,360)
      REAL phy_rug(klon,360)
      REAL phy_ice(klon,360)
c
      REAL masque(iip1,jjp1)
      REAL mask(iim,jjp1)

C Declarations pour le champ de depart
      INTEGER imdep, jmdep,lmdep
      INTEGER  tbid
      PARAMETER ( tbid = 60 )        ! >52 semaines
      REAL  timecoord(tbid)
c
      REAL , ALLOCATABLE :: dlon_msk(:), dlat_msk(:)
      REAL , ALLOCATABLE :: lonmsk_ini(:), latmsk_ini(:)
      REAL , ALLOCATABLE :: dlon(:), dlat(:)
      REAL , ALLOCATABLE :: dlon_ini(:), dlat_ini(:)
      REAL , ALLOCATABLE :: champ_msk(:), champ(:)
      REAL , ALLOCATABLE :: work(:,:)

      CHARACTER*25 title

C Declarations pour le champ interpole 2D
      REAL champint(iim,jjp1)
      real chmin,chmax

C Declarations pour le champ interpole 3D
      REAL champtime(iim,jjp1,tbid)
      REAL timeyear(tbid)
      REAL champan(iip1,jjp1,366)

C Declarations pour l'inteprolation verticale
      REAL ax(tbid), ay(tbid)
      REAL by
      REAL yder(tbid)


      INTEGER ierr
      INTEGER dimfirst(3)
      INTEGER dimlast(3)
c
      INTEGER nid, ndim, ntim
      INTEGER dims(2), debut(2), epais(2)
      INTEGER id_tim
      INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB

      INTEGER i, j, k, l
c Diverses variables locales
      REAL time

      INTEGER          longcles
      PARAMETER      ( longcles = 20 )
      REAL  clesphy0 ( longcles      )
#include "serre.h"
      INTEGER ncid,varid,ndimid(4),dimid
      character*30 namedim

      pi     = 4. * ATAN(1.)
      rad    = 6 371 229.
      omeg   = 4.* ASIN(1.)/(24.*3600.)
      g      = 9.8
      daysec = 86400.
      kappa  = 0.2857143
      cpp    = 1004.70885
      dtvr    = daysec/FLOAT(day_step)
c
C Traitement du relief au sol
c
      write(*,*) 'Traitement du relief au sol pour fabriquer masque'
      ierr = NF_OPEN('Relief.nc', NF_NOWRITE, ncid)

      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF

      ierr = NF_INQ_VARID(ncid,'RELIEF',varid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
      print*,'variable ', namedim, 'dimension ', imdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF

      ALLOCATE( lonmsk_ini(imdep) )
      ALLOCATE(   dlon_msk(imdep) )

#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,lonmsk_ini)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,lonmsk_ini)
#endif

c
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
      print*,'variable ', namedim, 'dimension ', jmdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF

      ALLOCATE( latmsk_ini(jmdep) )
      ALLOCATE(   dlat_msk(jmdep) )
      ALLOCATE(  champ_msk(imdep*jmdep) )

#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,latmsk_ini)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,latmsk_ini)
#endif
c
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,varid,champ_msk)
#else
      ierr = NF_GET_VAR_REAL(ncid,varid,champ_msk)
#endif
c
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
c
      title='RELIEF'

      CALL conf_dat2d(title,imdep, jmdep, lonmsk_ini, latmsk_ini,
     . dlon_msk, dlat_msk, champ_msk, interbar  )

      CALL mask_c_o(imdep, jmdep, dlon_msk, dlat_msk,champ_msk,
     .             iim, jjp1, rlonv, rlatu, champint)
      CALL gr_int_dyn(champint, masque, iim, jjp1)
      DO i = 1, iim
         masque(i,1) = FLOAT(NINT(masque(i,1)))
         masque(i,jjp1) = FLOAT(NINT(masque(i,jjp1)))
      ENDDO
      DO i = 1, iim
      DO j = 1, jjp1
         mask(i,j) = champint(i,j)
      ENDDO
      ENDDO
      CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, phy_nat0)
      ierr = NF_CLOSE(ncid)
c
c
C Traitement de la rugosite
c
      PRINT*, 'Traitement de la rugosite'
      ierr = NF_OPEN('Rugos.nc', NF_NOWRITE, ncid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF

      ierr = NF_INQ_VARID(ncid,'RUGOS',varid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      print*,'variable ', namedim, 'dimension ', imdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 

      ALLOCATE( dlon_ini(imdep) )
      ALLOCATE(     dlon(imdep) )

#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
#endif
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      print*,'variable ', namedim, 'dimension ', jmdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 

      ALLOCATE( dlat_ini(jmdep) )
      ALLOCATE(     dlat(jmdep) )

#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
#endif
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
      print*,'variable ', namedim, 'dimension ', lmdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
#endif
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
c
      ALLOCATE( champ(imdep*jmdep) )

      DO  200 l = 1, lmdep
         dimfirst(1) = 1
         dimfirst(2) = 1
         dimfirst(3) = l
c
         dimlast(1) = imdep
         dimlast(2) = jmdep
         dimlast(3) = 1
c
         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
         print*,dimfirst,dimlast
#ifdef NC_DOUBLE
         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
#else
         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
#endif
         if (ierr.ne.0) then
           print *, NF_STRERROR(ierr)
           STOP
         ENDIF 
   
        title = 'Rugosite Amip '
c
        CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
     .                      dlon, dlat, champ, interbar          )

       IF ( interbar )   THEN
         DO j = 1, imdep * jmdep
           champ(j) = LOG(champ(j))
         ENDDO

        IF( l.EQ.1 )  THEN
         WRITE(6,*) '-------------------------------------------------',
     ,'------------------------'
         WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
     , ' pour la rugosite $$$ '
         WRITE(6,*) '-------------------------------------------------',
     ,'------------------------'
        ENDIF
        CALL inter_barxy ( imdep,jmdep -1,dlon,dlat,champ ,
     ,                  iim,jjm,rlonu,rlatv, jjp1,champint )
         DO j=1,jjp1
          DO i=1,iim
           champint(i,j)=EXP(champint(i,j))
          ENDDO
         ENDDO

         DO j = 1, jjp1
           DO i = 1, iim
             IF(NINT(mask(i,j)).NE.1)  THEN
               champint( i,j ) = 0.001
             ENDIF
           ENDDO
         ENDDO
      ELSE
         CALL rugosite(imdep, jmdep, dlon, dlat, champ,
     .             iim, jjp1, rlonv, rlatu, champint, mask)
      ENDIF
         DO j = 1,jjp1
         DO i = 1, iim
            champtime (i,j,l) = champint(i,j)
         ENDDO
         ENDDO
200      CONTINUE
c
      DO l = 1, lmdep
         timeyear(l) = timecoord(l)
      ENDDO

      PRINT 222, timeyear
222   FORMAT(2x,' Time year ',10f6.1)
c
        
      PRINT*, 'Interpolation temporelle dans l annee'

      DO j = 1, jjp1
      DO i = 1, iim
          DO l = 1, lmdep
            ax(l) = timeyear(l)
            ay(l) = champtime (i,j,l)
          ENDDO
          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
          DO k = 1, 360
            time = FLOAT(k-1)
            CALL SPLINT(ax,ay,yder,lmdep,time,by)
            champan(i,j,k) = by
          ENDDO
      ENDDO
      ENDDO
      DO k = 1, 360
      DO j = 1, jjp1
         champan(iip1,j,k) = champan(1,j,k)
      ENDDO
        IF ( k.EQ.10 )  THEN
          DO j = 1, jjp1
            CALL minmax( iip1,champan(1,j,10),chmin,chmax )
            PRINT *,' Rugosite au temps 10 ', chmin,chmax,j
          ENDDO
        ENDIF
      ENDDO
c
      DO k = 1, 360
         CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k), phy_rug(1,k))
      ENDDO
c
      ierr = NF_CLOSE(ncid)

       DEALLOCATE( dlon      )
       DEALLOCATE( dlon_ini  )
       DEALLOCATE( dlat      )
       DEALLOCATE( dlat_ini  )
       DEALLOCATE( champ     )
c
c
C Traitement de la glace oceanique
c
      PRINT*, 'Traitement de la glace oceanique'
      ierr = NF_OPEN('AMIP.nc', NF_NOWRITE, ncid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF

      ierr = NF_INQ_VARID(ncid,'SEA_ICE',varid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      print*,'variable ', namedim, 'dimension ', imdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 

      ALLOCATE ( dlon_ini(imdep) )
      ALLOCATE (     dlon(imdep) )

#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
#endif
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      print*,'variable ', namedim, jmdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 

      ALLOCATE ( dlat_ini(jmdep) )
      ALLOCATE (     dlat(jmdep) )

#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
#endif
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      print*,'variable ', namedim, lmdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
#endif
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
c
      ALLOCATE ( champ(imdep*jmdep) )

      DO l = 1, lmdep
         dimfirst(1) = 1
         dimfirst(2) = 1
         dimfirst(3) = l
c
         dimlast(1) = imdep
         dimlast(2) = jmdep
         dimlast(3) = 1
c
         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
#ifdef NC_DOUBLE
         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
#else
         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
#endif
         if (ierr.ne.0) then
           print *, NF_STRERROR(ierr)
           STOP
         ENDIF
 
         title = 'Sea-ice Amip '
c
         IF( oldice )  THEN
          CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
     .                        dlon, dlat, champ, .FALSE.           )
         ELSE
          CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
     .                        dlon, dlat, champ, interbar          )
         ENDIF
c
      IF( oldice )  THEN
                 CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
     .             iim, jjp1, rlonv, rlatu, champint )
      ELSEIF ( interbar )  THEN
       IF( l.EQ.1 )  THEN
        WRITE(6,*) '-------------------------------------------------',
     ,'------------------------'
        WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
     , ' pour Sea-ice Amip  $$$ '
        WRITE(6,*) '-------------------------------------------------',
     ,'------------------------'
       ENDIF

         CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
     ,     champ, iim, jjm, rlonu, rlatv, jjp1, champint )
      ELSE
         CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
     .             iim, jjp1, rlonv, rlatu, champint )
      ENDIF
         DO j = 1,jjp1
         DO i = 1, iim
            champtime (i,j,l) = champint(i,j)
         ENDDO
         ENDDO
      ENDDO
c
      DO l = 1, lmdep
         timeyear(l) = timecoord(l)
      ENDDO
      PRINT 222,  timeyear
c
      PRINT*, 'Interpolation temporelle'
      DO j = 1, jjp1
      DO i = 1, iim
          DO l = 1, lmdep
            ax(l) = timeyear(l)
            ay(l) = champtime (i,j,l)
          ENDDO
          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
          DO k = 1, 360
            time = FLOAT(k-1)
            CALL SPLINT(ax,ay,yder,lmdep,time,by)
            champan(i,j,k) = by
          ENDDO
      ENDDO
      ENDDO
      DO k = 1, 360
      DO j = 1, jjp1
         champan(iip1, j, k) = champan(1, j, k)
      ENDDO
        IF ( k.EQ.10 )  THEN
          DO j = 1, jjp1
            CALL minmax( iip1,champan(1,j,10),chmin,chmax )
            PRINT *,' Sea ice au temps 10 ', chmin,chmax,j
          ENDDO
        ENDIF
      ENDDO
c

      DO k = 1, 360
         CALL gr_dyn_fi(1, iip1, jjp1, klon,
     .                  champan(1,1,k), phy_ice(1,k))
         DO i = 1, klon
            phy_nat(i,k) = phy_nat0(i)
            IF ( (phy_ice(i,k) - 0.5).GE.1.e-5 ) THEN
               IF (NINT(phy_nat0(i)).EQ.0) THEN
                  phy_nat(i,k) = 3.0
               ELSE
                  phy_nat(i,k) = 2.0
               ENDIF
            ENDIF
            IF( NINT(phy_nat(i,k)).EQ.0 ) THEN
              IF ( phy_rug(i,k).NE.0.001 ) phy_rug(i,k) = 0.001
            ENDIF
         ENDDO

      ENDDO
c

      ierr = NF_CLOSE(ncid)
c
       DEALLOCATE( dlon      )
       DEALLOCATE( dlon_ini  )
       DEALLOCATE( dlat      )
       DEALLOCATE( dlat_ini  )
       DEALLOCATE( champ     )

477    continue
c
C Traitement de la sst
c
      PRINT*, 'Traitement de la sst'
      ierr = NF_OPEN('AMIP.nc', NF_NOWRITE, ncid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF

      ierr = NF_INQ_VARID(ncid,'SST',varid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      print*,'variable SST ', namedim,'dimension ', imdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
 
      ALLOCATE( dlon_ini(imdep) )
      ALLOCATE(     dlon(imdep) )

#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
#endif

      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      print*,'variable SST ', namedim, 'dimension ', jmdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 

      ALLOCATE( dlat_ini(jmdep) )
      ALLOCATE(     dlat(jmdep) )

#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
#endif
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      print*,'variable ', namedim, 'dimension ', lmdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
#endif
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 

       ALLOCATE( champ(imdep*jmdep) )
       IF( extrap )   THEN
         ALLOCATE ( work(imdep,jmdep) )
       ENDIF
c
      DO l = 1, lmdep
         dimfirst(1) = 1
         dimfirst(2) = 1
         dimfirst(3) = l
c
         dimlast(1) = imdep
         dimlast(2) = jmdep
         dimlast(3) = 1
c
         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
#ifdef NC_DOUBLE
         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
#else
         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
#endif
         if (ierr.ne.0) then
           print *, NF_STRERROR(ierr)
           STOP
         ENDIF

         title='Sst Amip'
c
         CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
     .                            dlon, dlat, champ, interbar     )
       IF ( extrap )  THEN
        CALL extrapol(champ, imdep, jmdep, 999999.,.TRUE.,.TRUE.,2,work)
       ENDIF
c

      IF ( interbar )  THEN
        IF( l.EQ.1 )  THEN
         WRITE(6,*) '-------------------------------------------------',
     ,'------------------------'
         WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
     , ' pour la Sst Amip $$$ '
         WRITE(6,*) '-------------------------------------------------',
     ,'------------------------'
        ENDIF
       CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
     , champ, iim, jjm, rlonu, rlatv, jjp1, champint )
      ELSE
       CALL grille_m(imdep, jmdep, dlon, dlat, champ,
     .          iim, jjp1, rlonv, rlatu, champint   )
      ENDIF

         DO j = 1,jjp1
         DO i = 1, iim
            champtime (i,j,l) = champint(i,j)
         ENDDO
         ENDDO
      ENDDO
c
      DO l = 1, lmdep
         timeyear(l) = timecoord(l)
      ENDDO
      print 222,  timeyear
c
C interpolation temporelle
      DO j = 1, jjp1
      DO i = 1, iim
          DO l = 1, lmdep
            ax(l) = timeyear(l)
            ay(l) = champtime (i,j,l)
          ENDDO
          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
          DO k = 1, 360
            time = FLOAT(k-1)
            CALL SPLINT(ax,ay,yder,lmdep,time,by)
            champan(i,j,k) = by
          ENDDO
      ENDDO
      ENDDO
      DO k = 1, 360
      DO j = 1, jjp1
         champan(iip1,j,k) = champan(1,j,k)
      ENDDO
        IF ( k.EQ.10 )  THEN
          DO j = 1, jjp1
            CALL minmax( iip1,champan(1,j,10),chmin,chmax )
            PRINT *,' SST au temps 10 ', chmin,chmax,j
          ENDDO
        ENDIF
      ENDDO
c
      DO k = 1, 360
         CALL gr_dyn_fi(1, iip1, jjp1, klon, 
     .                  champan(1,1,k), phy_sst(1,k))
      ENDDO
c
      ierr = NF_CLOSE(ncid)
c
       DEALLOCATE( dlon      )
       DEALLOCATE( dlon_ini  )
       DEALLOCATE( dlat      )
       DEALLOCATE( dlat_ini  )
       DEALLOCATE( champ     )
c
C Traitement de l'albedo
c
      PRINT*, 'Traitement de l albedo'
      ierr = NF_OPEN('Albedo.nc', NF_NOWRITE, ncid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
      ierr = NF_INQ_VARID(ncid,'ALBEDO',varid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
      ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      print*,'variable ', namedim, 'dimension ', imdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 

      ALLOCATE ( dlon_ini(imdep) )
      ALLOCATE (     dlon(imdep) )

#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_ini)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_ini)
#endif
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      print*,'variable ', namedim, 'dimension ', jmdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 

      ALLOCATE ( dlat_ini(jmdep) )
      ALLOCATE (     dlat(jmdep) )

#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_ini)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_ini)
#endif
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      print*,'variable ', namedim, 'dimension ', lmdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
#endif
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
c
      ALLOCATE ( champ(imdep*jmdep) )

      DO l = 1, lmdep
         dimfirst(1) = 1
         dimfirst(2) = 1
         dimfirst(3) = l
c
         dimlast(1) = imdep
         dimlast(2) = jmdep
         dimlast(3) = 1
c
         PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
#ifdef NC_DOUBLE
         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
#else
         ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
#endif
         if (ierr.ne.0) then
           print *, NF_STRERROR(ierr)
           STOP
         ENDIF

         title='Albedo Amip'
c
         CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
     .                            dlon, dlat, champ, interbar      )
c
c
      IF ( interbar )  THEN
        IF( l.EQ.1 )  THEN
         WRITE(6,*) '-------------------------------------------------',
     ,'------------------------'
         WRITE(6,*) '$$$ Utilisation de l interpolation barycentrique ',
     , ' pour l Albedo Amip $$$ '
         WRITE(6,*) '-------------------------------------------------',
     ,'------------------------'
        ENDIF

       CALL inter_barxy ( imdep,jmdep -1,dlon, dlat ,
     , champ, iim, jjm, rlonu, rlatv, jjp1, champint )
      ELSE
       CALL grille_m(imdep, jmdep, dlon, dlat, champ,
     .          iim, jjp1, rlonv, rlatu, champint   )
      ENDIF
c
         DO j = 1,jjp1
         DO i = 1, iim
            champtime (i, j, l) = champint(i, j)
         ENDDO
         ENDDO
      ENDDO
c
      DO l = 1, lmdep
         timeyear(l) = timecoord(l)
      ENDDO
      print 222,  timeyear
c
C interpolation temporelle
      DO j = 1, jjp1
      DO i = 1, iim
          DO l = 1, lmdep
            ax(l) = timeyear(l)
            ay(l) = champtime (i, j, l)
          ENDDO
          CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
          DO k = 1, 360
            time = FLOAT(k-1)
            CALL SPLINT(ax,ay,yder,lmdep,time,by)
            champan(i,j,k) = by
          ENDDO
      ENDDO
      ENDDO
      DO k = 1, 360
      DO j = 1, jjp1
         champan(iip1, j, k) = champan(1, j, k)
      ENDDO
        IF ( k.EQ.10 )  THEN
          DO j = 1, jjp1
            CALL minmax( iip1,champan(1,j,10),chmin,chmax )
            PRINT *,' Albedo au temps 10 ', chmin,chmax,j
          ENDDO
        ENDIF
      ENDDO
c
      DO k = 1, 360
         CALL gr_dyn_fi(1, iip1, jjp1, klon,
     .                  champan(1,1,k), phy_alb(1,k))
      ENDDO
c
      ierr = NF_CLOSE(ncid)
c
c
      DO k = 1, 360
      DO i = 1, klon
         phy_bil(i,k) = 0.0
      ENDDO
      ENDDO
c
      PRINT*, 'Ecriture du fichier limit'
c
      ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
c
      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
     .                       "Fichier conditions aux limites")
      ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
      ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
c
      dims(1) = ndim
      dims(2) = ntim
c
ccc      ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
      ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
      ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
     .                        "Jour dans l annee")
ccc      ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
      ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
      ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
     .                        "Nature du sol (0,1,2,3)")
ccc      ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
      ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
      ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
     .                        "Temperature superficielle de la mer")
ccc      ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
      ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
      ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
     .                        "Reference flux de chaleur au sol")
ccc      ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
      ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
      ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
     .                        "Albedo a la surface")
ccc      ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
      ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
      ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
     .                        "Rugosite")
c
      ierr = NF_ENDDEF(nid)
c
      DO k = 1, 360
c
      debut(1) = 1
      debut(2) = k
      epais(1) = klon
      epais(2) = 1
c
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
      ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais,phy_nat(1,k))
      ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k))
      ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k))
      ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k))
      ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k))
#else
      ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
      ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais,phy_nat(1,k))
      ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k))
      ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k))
      ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k))
      ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k))
#endif
c
      ENDDO
c
      ierr = NF_CLOSE(nid)
c
      STOP
      END