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 modifie par P. Braconnot pour utiliser la version sous-surfaces
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)
CPB
      REAL phy_icet(klon,360)
      REAL phy_oce(klon,360)
      REAL verif
c
      REAL masque(iip1,jjp1)
      REAL mask(iim,jjp1)
CPB
C newlmt indique l'utilisation de la sous-maille fractionnelle
C tandis que l'ancien codage utilise l'indicateur du sol (0,1,2,3)
      LOGICAL newlmt, fracterre
      PARAMETER(newlmt=.TRUE.)
      PARAMETER(fracterre = .TRUE.) 
CPB
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)
      REAL*4 zbidon(ibid_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
CPB
      INTEGER id_FOCE, id_FSIC, id_FTER, id_FLIC

      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
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon_msk)
      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_REAL(ncid,dimid,dlat_msk)
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr)
        STOP
      ENDIF
      ierr = NF_GET_VAR_REAL(ncid,varid,champ_msk)
      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)
      IF ( fracterre ) THEN
          DO i = 1, iim
            masque(i,1) = masque(i,1)
            masque(i,jjp1) = masque(i,jjp1)
          END DO
      ELSE 
          DO i = 1, iim
            masque(i,1) = FLOAT(NINT(masque(i,1)))
            masque(i,jjp1) = FLOAT(NINT(masque(i,jjp1)))
          END DO
      ENDIF 
      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_REAL(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_REAL(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_REAL(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_REAL(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_REAL(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_REAL(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_REAL(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_REAL(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
      WRITE(*,*) 'phy_nat'
c     WRITE(*,'(72f4.1)') phy_nat0(1:klon)
c
      DO k = 1, 360
        CALL gr_dyn_fi(1, iip1, jjp1, klon,
     .      champan(1,1,k), phy_ice(1,k))
        IF ( newlmt) THEN

CPB  en attendant de mettre fraction de terre
c
            WHERE(phy_ice(1:klon, k) .GT. 1.) phy_ice(1 : klon, k) = 1.
            WHERE(phy_ice(1:klon, k) .LT. 0.) phy_ice(1 : klon, k) = 0.
            WRITE(*,*) 'phy_ice : ', k
c           WRITE(*,'(72f4.1)') phy_ice(1 : klon, k)
c 
            IF (fracterre ) THEN
                WRITE(*,*) 'passe dans cas fracterre' 
                DO i = 1, klon
                  phy_nat(i,k) = phy_nat0(i)
                  IF (phy_nat0(i) .GE. 0.5 ) THEN
                      IF(phy_ice(i,k) .GE. 1.e-5) THEN 
                          IF ( phy_ice(i,k) .LE. phy_nat(i,k)) THEN
                              phy_oce(i,k) = 1. - phy_nat(i,k)
                              phy_icet(i,k) = phy_ice(i,k)
                              phy_ice(i,k) = 0.
                              phy_nat(i,k)= phy_nat(i,k)- phy_icet(i,k)
                          ELSE 
                              phy_oce(i,k) = 1. - phy_ice(i,k)
                              phy_icet(i,k) = phy_nat(i,k)
                              phy_nat(i,k) = 0.
                              phy_ice(i,k) = phy_ice(i,k) 
     $                            - phy_icet(i,k)
                          ENDIF 
                      ELSE
                          phy_icet(i,k) = 0.
                          phy_ice(i,k) = 0.
                          phy_oce(i,k) = 1. - phy_nat(i,k)
                      ENDIF 
                  ELSE
                      phy_oce(i,k) = 1. - phy_nat(i,k)
                      IF(phy_ice(i,k) .GE. 1.e-5) THEN 
                          IF( phy_ice(i,k) .LE. phy_oce(i,k) ) THEN 
                              phy_icet(i,k) = 0.
                              phy_oce(i,k) = phy_oce(i,k) - phy_ice(i,k)
                          ELSE
                              phy_icet(i,k)=phy_ice(i,k) - phy_oce(i,k)
                              phy_ice(i,k) = phy_oce(i,k)
                              phy_oce(i,k) = 0.
                              phy_nat(i,k) = phy_nat(i,k)-phy_icet(i,k)
                          ENDIF 
                      ELSE
                          phy_icet(i,k) = 0.
                          phy_ice(i,k) = 0.
                      ENDIF 
                  ENDIF
                  verif = phy_nat(i,k) + phy_icet(i,k)+ phy_ice(i,k) 
     $                + phy_oce(i,k)
                  IF ( verif .LT. 1. - 1.e-5 .OR. 
     $                verif .GT. 1 + 1.e-5) THEN  
                      WRITE(*,*) 'pb sous maille au point : i,k,verif '
     $                    , i,k,verif
                  ENDIF 
                END DO
            ELSE 
                DO i = 1, klon
                  phy_nat(i,k) = phy_nat0(i)
                  IF (NINT(phy_nat0(i)).EQ.1 ) THEN
                      IF(phy_ice(i,k) .GE. 1.e-5) THEN 
                          phy_icet(i,k) = phy_ice(i,k)
                          phy_ice(i,k) = 0.
                          phy_nat(i,k) = phy_nat(i,k) - phy_icet(i,k)
                          phy_oce(i,k) = 0.
                      ELSE
                          phy_icet(i,k) = 0.
                          phy_ice(i,k) = 0.
                          phy_oce(i,k) = 0.                      
                      ENDIF 
                  ELSE 
                      IF(phy_ice(i,k) .GE. 1.e-5) THEN 
                          phy_icet(i,k) = 0.
                          phy_nat(i,k) = 0.
                          phy_oce(i,k) = 1. - phy_ice(i,k)
                      ELSE
                          phy_icet(i,k) = 0.
                          phy_ice(i,k) = 0.
                          phy_oce(i,k) = 1.                      
                      ENDIF 
                  ENDIF
                  verif = phy_nat(i,k) + phy_icet(i,k)+ phy_ice(i,k) 
     $                + phy_oce(i,k)
                  IF ( verif .LT. 1. - 1.e-5 .OR. 
     $                verif .GT. 1 + 1.e-5) THEN  
                      WRITE(*,*) 'pb sous maille au point : i,k,verif '
     $                    , i,k,verif
                  ENDIF 
                END DO
            ENDIF 
        ELSE  
            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
            END DO
        ENDIF 
      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_REAL(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_REAL(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_REAL(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_REAL(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_REAL(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_REAL(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_REAL(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_REAL(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
      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")
      IF (newlmt) THEN
c
          ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
          ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 14,
     .                        "Fraction ocean")
c
          ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
          ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 21,
     .                        "Fraction glace de mer")
c
          ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
          ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 14,
     .                        "Fraction terre")
c
          ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
          ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 17,
     .                        "Fraction land ice")
c
      ELSE 
          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)")
      ENDIF 
C 
      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")
      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")
      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")
      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))
c
      IF (newlmt ) THEN
          ierr = NF_PUT_VARA_DOUBLE (nid,id_FOCE,debut,epais
     $        ,phy_oce(1,k))
          ierr = NF_PUT_VARA_DOUBLE (nid,id_FSIC,debut,epais
     $        ,phy_ice(1,k))
          ierr = NF_PUT_VARA_DOUBLE (nid,id_FTER,debut,epais
     $        ,phy_nat(1,k))
          ierr = NF_PUT_VARA_DOUBLE (nid,id_FLIC,debut,epais
     $        ,phy_icet(1,k))
      ELSE 
          ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais
     $        ,phy_nat(1,k))
      ENDIF 
c
      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))
      IF (newlmt ) THEN
          ierr = NF_PUT_VARA_REAL (nid,id_FOCE,debut,epais
     $        ,phy_oce(1,k))
          ierr = NF_PUT_VARA_REAL (nid,id_FSIC,debut,epais
     $        ,phy_ice(1,k))
          ierr = NF_PUT_VARA_REAL (nid,id_FTER,debut,epais
     $        ,phy_nat(1,k))
          ierr = NF_PUT_VARA_REAL (nid,id_FLIC,debut,epais
     $        ,phy_icet(1,k))
      ELSE 
          ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais
     $        ,phy_nat(1,k))
      ENDIF 
      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