!
! $Id: limit_netcdf.F 1239 2009-09-08 09:16:12Z evignon $
!
C
C
      SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque)
#ifdef CPP_EARTH
! This routine is designed to work for Earth
      USE dimphy
      use phys_state_var_mod , ONLY : pctsrf
      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"
cy#include "dimphy.h"
#include "indicesol.h"
#include "iniprint.h"
c
c-----------------------------------------------------------------------
      LOGICAL interbar, extrap, oldice

      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)
c
      real pctsrf_t(klon,nbsrf,360)

      REAL verif

      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.) 

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
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

cIM28/02/2002 <== PM
      REAL tmidmonth(12)
      SAVE tmidmonth
      DATA tmidmonth/15,45,75,105,135,165,195,225,255,285,315,345/

c initialisations:
      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)
      CALL inigeom
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  )

      DO i = 1, iim
      DO j = 1, jjp1
         mask(i,j) = masque(i,j)
      ENDDO
      ENDDO
      WRITE(*,*) 'MASK:'
      WRITE(*,'(96i1)')INT(mask)     
      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(:lmdep)
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('amipbc_sic_1x1.nc', NF_NOWRITE, ncid)
      if (ierr.ne.0) THEN
        ierr = NF_OPEN('amipbc_sic_1x1_clim.nc', NF_NOWRITE, ncid)
        if (ierr.ne.0) THEN
          print *, NF_STRERROR(ierr)
          STOP
        endif
      ENDIF

cIM22/02/2002
cIM07/03/2002 AMIP.nc & amip79to95.nc
cIM   ierr = NF_INQ_VARID(ncid,'SEA_ICE',varid)
cIM07/03/2002 amipbc_sic_1x1_clim.nc & amipbc_sic_1x1.nc
      ierr = NF_INQ_VARID(ncid,'sicbcs',varid)
cIM22/02/2002
      if (ierr.ne.0) then
        print *, NF_STRERROR(ierr),'sicbcs'
        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
cIM28/02/2002
cPM28/02/2002 : nouvelle coord temporelle fichiers AMIP pas en jours
c               Ici on suppose qu'on a 12 mois (de 30 jours).
      IF (lmdep.NE.12) THEN
          print *, 'Unknown AMIP file: not 12 months ?'
          STOP
       ENDIF
cIM28/02/2002
      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
         CALL conf_dat2d(title,imdep, jmdep, dlon_ini, dlat_ini,
     .                        dlon, dlat, champ, interbar          )
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
cIM07/03/2002 
cIM22/02/2002 : Sea-ice Amip entre 0. et 1.
cIM    PRINT*,'SUB. limit_netcdf.F IM : Sea-ice et SST Amip_new clim' 
cIM   DO j = 1, imdep * jmdep
cIM28/02/2002 <==PM         champ(j) = champ(j)/100.
cIM14/03/2002      champ(j) = max(0.0,(min(1.0, (champ(j)/100.) )))
cIM      champ(j) = amax1(0.0,(amin1(1.0, (champ(j)/100.) )))
cIM   ENDDO
cIM22/02/2002
         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
cIM28/02/2002 <== PM  timeyear(l) = timecoord(l)
cIM      timeyear(l) = timecoord(l)
cIM07/03/2002      
         timeyear(l) = tmidmonth(l)
      ENDDO
      PRINT 222,  timeyear(:lmdep)
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
cIM14/03/2002 : Sea-ice Amip entre 0. et 1.
      PRINT*,'SUB. limit_netcdf.F IM : Sea-ice Amipbc '
      DO k = 1, 360
      DO j = 1, jjp1
      DO i = 1, iim
        champan(i, j, k) = 
     $ amax1(0.0,(amin1(1.0,(champan(i, j, k)/100.))))
      ENDDO
        champan(iip1, j, k) = champan(1, j, k)
      ENDDO
      ENDDO
cIM14/03/2002

      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
              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
            IF( NINT(phy_nat(i,k)).EQ.0 ) THEN
              IF ( phy_rug(i,k).NE.0.001 ) phy_rug(i,k) = 0.001
            ENDIF
          END DO
        ENDIF
      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'
c     ierr = NF_OPEN('AMIP_SST.nc', NF_NOWRITE, ncid)
      ierr = NF_OPEN('amipbc_sst_1x1.nc', NF_NOWRITE, ncid)
      if (ierr.ne.0) THEN
        ierr = NF_OPEN('amipbc_sst_1x1_clim.nc', NF_NOWRITE, ncid)
        if (ierr.ne.0) THEN
          print *, NF_STRERROR(ierr)
          STOP
        endif
      ENDIF

cIM22/02/2002
cIM   ierr = NF_INQ_VARID(ncid,'SST',varid)
      ierr = NF_INQ_VARID(ncid,'tosbcs',varid)
cIM22/02/2002
      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
cIM28/02/2002
cPM28/02/2002 : nouvelle coord temporelle fichiers AMIP pas en jours
c               Ici on suppose qu'on a 12 mois (de 30 jours).
      IF (lmdep.NE.12) THEN
          print *, 'Unknown AMIP file: not 12 months ?'
          STOP
       ENDIF
cIM28/02/2002
      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
cIM28/02/2002 <==PM  timeyear(l) = timecoord(l)
         timeyear(l) = tmidmonth(l)
      ENDDO
      print 222,  timeyear(:lmdep)
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
cIM14/03/2002 : SST amipbc greater then 271.38
      PRINT*,'SUB. limit_netcdf.F IM : SST Amipbc >= 271.38 '
      DO k = 1, 360
      DO j = 1, jjp1
      DO i = 1, iim
         champan(i, j, k) = amax1(champan(i, j, k), 271.38)
      ENDDO
         champan(iip1, j, k) = champan(1, j, k)
      ENDDO
      ENDDO
cIM14/03/2002
      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(:lmdep)
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
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
#else
      ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
     .                        "Jour dans l annee")
      IF (newlmt) THEN
c
#ifdef NC_DOUBLE
        ierr = NF_DEF_VAR (nid, "FOCE", NF_DOUBLE, 2,dims, id_FOCE)
#else
        ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
#endif
        ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 14,
     .                      "Fraction ocean")
c
#ifdef NC_DOUBLE
        ierr = NF_DEF_VAR (nid, "FSIC", NF_DOUBLE, 2,dims, id_FSIC)
#else
        ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
#endif
        ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 21,
     .                      "Fraction glace de mer")
c
#ifdef NC_DOUBLE
        ierr = NF_DEF_VAR (nid, "FTER", NF_DOUBLE, 2,dims, id_FTER)
#else
        ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
#endif
        ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 14,
     .                      "Fraction terre")
c
#ifdef NC_DOUBLE
        ierr = NF_DEF_VAR (nid, "FLIC", NF_DOUBLE, 2,dims, id_FLIC)
#else
        ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
#endif
        ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 17,
     .                      "Fraction land ice")
c
      ELSE 
#ifdef NC_DOUBLE
        ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
#else
        ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
#endif
        ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
     .                      "Nature du sol (0,1,2,3)")
      ENDIF 
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
#else
      ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
     .                      "Temperature superficielle de la mer")
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
#else
      ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
     .                        "Reference flux de chaleur au sol")
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
#else
      ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
     .                        "Albedo a la surface")
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
#else
      ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
#endif
      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
#else
      WRITE(lunout,*)
     & 'limit_netcdf: Earth-specific routine, needs Earth physics'
#endif
! of #ifdef CPP_EARTH
      STOP
      END