PROGRAM create_limit
      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-------------------------------------------------------------
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-----------------------------------------------------------------------
      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 ibid, jbid, tbid
      PARAMETER (ibid = 400,       ! >360 pts
     .           jbid = 200,       ! >181 pts
     .           tbid = 60)        ! >52 semaines
      REAL champ(ibid*jbid)
      REAL dlon(ibid), dlat(jbid), timecoord(tbid)
c
      INTEGER ibid_msk, jbid_msk
      PARAMETER(ibid_msk=2200,jbid_msk=1100)
      REAL champ_msk(ibid_msk*jbid_msk)
      REAL dlon_msk(ibid_msk), dlat_msk(jbid_msk)

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

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

c initialisations:
      OPEN (8,file='run.def',form='formatted')
      CALL defrun_new(8,.TRUE.,clesphy0)
      CLOSE(8)

      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
ccc      CALL iniconst ( non  indispensable )

      CALL inigeom
c
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
c     ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_msk)
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon_msk)
c$$$      zbidon=0.0
c$$$      ierr = NF_GET_VAR_REAL(ncid,dimid,zbidon(1:imdep))
c$$$      dlon_msk(1 : imdep) = dble(zbidon(1:imdep))
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
c      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat_msk)
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat_msk)
c$$$      zbidon=0.
c$$$      ierr = NF_GET_VAR_REAL(ncid,dimid,zbidon(1:jmdep))
c$$$      dlat_msk=dble(zbidon(1:jmdep))
c
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
c      ierr = NF_GET_VAR_REAL(ncid,varid,champ_msk)
      ierr = NF_GET_VAR_DOUBLE(ncid,varid,champ_msk)
c$$$      zbidon=0.
c$$$      ierr = NF_GET_VAR_REAL(ncid,varid,zbidon(1:imdep*jmdep))
c$$$      champ_msk(1: imdep*jmdep) = zbidon(1:imdep*jmdep)
c
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF

c
      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 
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
      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 
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
      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
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      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)
         print*,dimfirst,dimlast
         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
         if (ierr.ne.0) then
           print *, NF_STRERROR(ierr)
           STOP
         ENDIF 
   
         CALL rugosite(imdep, jmdep, dlon, dlat, champ,
     .             iim, jjp1, rlonv, rlatu, champint, mask)
         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
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
      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)
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 
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
      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 
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
      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 
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      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)
         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
         if (ierr.ne.0) then
           print *, NF_STRERROR(ierr)
           STOP
         ENDIF

         CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
     .             iim, jjp1, rlonv, rlatu, champint)
         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
      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
         ENDDO
      ENDDO
c
      ierr = NF_CLOSE(ncid)
c
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 ', namedim,'dimension ', imdep
      ierr = NF_INQ_VARID(ncid,namedim,dimid)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF 
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
      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 
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
      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 
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      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)
         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
         if (ierr.ne.0) then
           print *, NF_STRERROR(ierr)
           STOP
         ENDIF
         CALL grille_m(imdep, jmdep, dlon, dlat, champ,
     .             iim, jjp1, rlonv, rlatu, champint)

         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
      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
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 
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
      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 
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
      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 
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      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)
         ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
         if (ierr.ne.0) then
           print *, NF_STRERROR(ierr)
           STOP
         ENDIF
         CALL grille_m(imdep, jmdep, dlon, dlat, champ,
     .             iim, jjp1, rlonv, rlatu, champint)
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
      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