c $Header$
      PROGRAM create_limit
      USE startvar
      USE ioipsl
      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"
#include "dimphy.h"
#include "indicesol.h"
c-----------------------------------------------------------------------
      REAL phy_nat(klon,360)
      real 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)
      logical interbar
      parameter (interbar=.false.)
CPB
c      REAL phy_icet(klon,360)
c      REAL phy_oce(klon,360)
      real pctsrf_t(klon,nbsrf,360)
      real pctsrf(klon,nbsrf)
      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, ji
c declarations pour lecture glace de mer
      INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
      INTEGER :: itaul(1), fid
      REAL :: lev(1), date, dt
      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_lic, dlat_lic
      REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
      REAL :: flic_tmp(iip1, jjp1)
c Diverses variables locales
      REAL time
! pour la lecture du fichier masque ocean
      integer :: nid_o2a
      logical :: couple = .false.
      INTEGER :: iml_omask, jml_omask
      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_omask, dlat_omask
      REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
      real, dimension(klon) :: ocemask_fi


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

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

      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(*,*) 'Fabrication masque'

      varname = 'masque'
      masque(:,:) = 0.0
      CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0, jjm
     ,     ,rlonu,rlatv , interbar )
      pctsrf=0.
      varname = 'zmasq'
      zmasq(:) = 0.
      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmasq,0.0,
     .jjm ,rlonu,rlatv , interbar)
      WHERE (zmasq(1 : klon) .LT. EPSFRA)
          zmasq(1 : klon) = 0.
      END WHERE 
      WHERE (1 - zmasq(1 : klon) .LT. EPSFRA)
          zmasq(1 : klon) = 1.
      END WHERE 
!      WRITE(*,*)zmasq

      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 
c$$$      DO i = 1, iim
c$$$      DO j = 1, jjp1
c$$$         mask(i,j) = masque(i,j)
c$$$      ENDDO
c$$$      ENDDO
c$$$      CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, phy_nat0)
      phy_nat0(1:klon) = zmasq(1:klon)
      mask = 0.
      DO j = 1, jjp1
        DO i = 1, iim
          IF ( masque(i,j) .GE. EPSFRA) mask (i,j) = 1
        END DO 
      END DO  
C
C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
C utilise pour calculer les poids et pour assurer l'adequation entre les
C fractions d'ocean vu par l'atmosphere et l'ocean
C

      write(*,*)'Essai de lecture masque ocean'
      iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
      if (iret .ne. 0) then
        write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve'
        write(*,*)'Run force'
      else
        couple = .true.
        iret = nf_close(nid_o2a)
        call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
     $    , nid_o2a)
        if (iml_omask /= iim .or. jml_omask /= jjp1) then
          write(*,*)'Dimensions non compatibles pour masque ocean'
          write(*,*)'iim = ',iim,' iml_omask = ',iml_omask
          write(*,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
          stop
        endif
        ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
        ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
        ALLOCATE(dlon_omask(iml_omask), stat=iret)
        ALLOCATE(dlat_omask(jml_omask), stat=iret)
        ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
        ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
        CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
     $    , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
        CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp, 
     $      ttm_tmp, 1, 1, ocetmp)
        CALL flinclo(fid)
        dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
        dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
        ocemask = ocetmp
        if (dlat_omask(1) < dlat_omask(jml_omask)) then
          do j = 1, jml_omask
            ocemask(:,j) = ocetmp(:,jml_omask-j+1)
          enddo
        endif 
C
C passage masque ocean a la grille physique
C
        ocemask_fi(1) = ocemask(1,1)
        do j = 2, jjm
          do i = 1, iim
            ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j)
          enddo
        enddo
        ocemask_fi(klon) = ocemask(1,jjp1)
        zmasq = 1. - ocemask_fi
      endif


C
C lecture du fichier glace de terre pour fixer la fraction de terre 
C et de glace de terre
C
      CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
     $    , fid)
      ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
      ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
      ALLOCATE(dlon_lic(iml_lic), stat=iret)
      ALLOCATE(dlat_lic(jml_lic), stat=iret)
      ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
      CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
     $    , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
      CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
     $    , 1, 1, fraclic)
      CALL flinclo(fid)
C
C interpolation sur la grille T du modele
C
      WRITE(*,*) 'dimensions de landice iml_lic, jml_lic : ', 
     $    iml_lic, jml_lic
c
C sil les coordonnees sont en degres, on les transforme
C
      IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) )  THEN
          lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
      ENDIF 
      IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN 
          lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
      ENDIF 

      dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
      dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic) 
C
      CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
     $    ,iim, jjp1,
     $    rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
c$$$      flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
      flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
C
C passage sur la grille physique
C
      CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
     $    pctsrf(1:klon, is_lic))
C adequation avec le maque terre/mer
      WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA ) 
          pctsrf(1 : klon, is_lic) = 0. 
      END WHERE
      WHERE (zmasq( 1 : klon) .LT. EPSFRA) 
          pctsrf(1 : klon, is_lic) = 0.
      END WHERE 
      pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
      DO ji = 1, klon
        IF (zmasq(ji) .GT. EPSFRA) THEN 
            IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
                pctsrf(ji, is_lic) = zmasq(ji)
                pctsrf(ji, is_ter) = 0.
            ELSE 
                pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
                IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
                    pctsrf(ji,is_ter) = 0.
                    pctsrf(ji, is_lic) = zmasq(ji)
                ENDIF 
            ENDIF 
        ENDIF 
      END DO 
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 
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
#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 
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
#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
      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
#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 
   
         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      write(70,*)champtime
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 
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
#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 
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
#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
      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

         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
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)
        IF ( newlmt) THEN

CPB  en attendant de mettre fraction de terre
c
          WHERE(phy_ice(1:klon) .GE. 1.) phy_ice(1 : klon) = 1.
          WHERE(phy_ice(1:klon) .LT. EPSFRA) phy_ice(1 : klon) = 0.
c 
          IF (fracterre ) THEN
c            WRITE(*,*) 'passe dans cas fracterre' 
            pctsrf_t(:,is_ter,k) = pctsrf(:,is_ter)
            pctsrf_t(:,is_lic,k) = pctsrf(:,is_lic)
            pctsrf_t(1:klon,is_sic,k) =   phy_ice(1:klon) 
     $            - pctsrf_t(1:klon,is_lic,k)
c§§ Il y a des cas ou il y a de la glace dans landiceref et pas dans AMIP
            WHERE (pctsrf_t(1:klon,is_sic,k) .LE. 0)
                pctsrf_t(1:klon,is_sic,k) = 0.
            END WHERE 
            WHERE( 1. - zmasq(1:klon) .LT. EPSFRA)
                pctsrf_t(1:klon,is_sic,k) = 0.
                pctsrf_t(1:klon,is_oce,k) = 0.
            END WHERE
            DO i = 1, klon
c$$              pctsrf_t(i,is_sic,k) = (1. - pctsrf_t(i,is_lic,k) -
c$$     .                               pctsrf_t(i,is_ter,k)) * phy_ice(i)
c$$              pctsrf_t(i,is_oce,k) = 1. - pctsrf_t(i,is_lic,k) -
c$$     .                      pctsrf_t(i,is_ter,k) - pctsrf_t(i,is_sic,k)
              IF ( 1. - zmasq(i) .GT. EPSFRA) THEN 
                  IF ( pctsrf_t(i,is_sic,k) .GE. 1 - zmasq(i)) THEN
                      pctsrf_t(i,is_sic,k) = 1 - zmasq(i)
                      pctsrf_t(i,is_oce,k) = 0.
                  ELSE 
                      pctsrf_t(i,is_oce,k) = 1 - zmasq(i) 
     $                    - pctsrf_t(i,is_sic,k)
                      IF (pctsrf_t(i,is_oce,k) .LT. EPSFRA) THEN
                          pctsrf_t(i,is_oce,k) = 0.
                          pctsrf_t(i,is_sic,k) = 1 - zmasq(i)
                      ENDIF 
                  ENDIF
              ENDIF  
              if (pctsrf_t(i,is_oce,k) .lt. 0.) then
                  WRITE(*,*) 'pb sous maille au point : i,k '
     $                , i,k,pctsrf_t(:,is_oce,k)
              ENDIF
              IF ( abs( pctsrf_t(i, is_ter,k) + pctsrf_t(i, is_lic,k) + 
     $            pctsrf_t(i, is_oce,k) + pctsrf_t(i, is_sic,k)  - 1.) 
     $            .GT. EPSFRA) THEN 
                  WRITE(*,*) 'physiq : pb sous surface au point ', i, 
     $                pctsrf_t(i, 1 : nbsrf,k), phy_ice(i)
              ENDIF 
            END DO
        ELSE 
            DO i = 1, klon
              pctsrf_t(i,is_ter,k) = pctsrf(i,is_ter)
              IF (NINT(pctsrf(i,is_ter)).EQ.1 ) THEN
                pctsrf_t(i,is_sic,k) = 0.
                pctsrf_t(i,is_oce,k) = 0.                  
                IF(phy_ice(i) .GE. 1.e-5) THEN
                  pctsrf_t(i,is_lic,k) = phy_ice(i)
                  pctsrf_t(i,is_ter,k) = pctsrf_t(i,is_ter,k) 
     .                                   - pctsrf_t(i,is_lic,k)
                ELSE
                  pctsrf_t(i,is_lic,k) = 0.
                ENDIF 
              ELSE
                pctsrf_t(i,is_lic,k) = 0. 
                IF(phy_ice(i) .GE. 1.e-5) THEN 
                  pctsrf_t(i,is_ter,k) = 0.
                  pctsrf_t(i,is_sic,k) = phy_ice(i)
                  pctsrf_t(i,is_oce,k) = 1. - pctsrf_t(i,is_sic,k)
                ELSE
                  pctsrf_t(i,is_sic,k) = 0.
                  pctsrf_t(i,is_oce,k) = 1.                      
                ENDIF 
              ENDIF
              verif = pctsrf_t(i,is_ter,k) +
     .                pctsrf_t(i,is_oce,k) + 
     .                pctsrf_t(i,is_sic,k) +
     .                pctsrf_t(i,is_lic,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) - 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 
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
#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 
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
#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
      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
         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
      WHERE(phy_sst .LT. 271.35) phy_sst = 271.35
      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 
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
#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 
#ifdef NC_DOUBLE
      ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
#else
      ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
#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
      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
         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
     $        ,pctsrf_t(1,is_oce,k))
          ierr = NF_PUT_VARA_DOUBLE (nid,id_FSIC,debut,epais
     $        ,pctsrf_t(1,is_sic,k))
          ierr = NF_PUT_VARA_DOUBLE (nid,id_FTER,debut,epais
     $        ,pctsrf_t(1,is_ter,k))
          ierr = NF_PUT_VARA_DOUBLE (nid,id_FLIC,debut,epais
     $        ,pctsrf_t(1,is_lic,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
     $        ,pctsrf_t(1,is_oce,k))
          ierr = NF_PUT_VARA_REAL (nid,id_FSIC,debut,epais
     $        ,pctsrf_t(1,is_sic,k))
          ierr = NF_PUT_VARA_REAL (nid,id_FTER,debut,epais
     $        ,pctsrf_t(1,is_ter,k))
          ierr = NF_PUT_VARA_REAL (nid,id_FLIC,debut,epais
     $        ,pctsrf_t(1,is_lic,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