Index: LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90
===================================================================
--- LMDZ4/trunk/libf/dyn3d/limit_netcdf.F	(revision 1279)
+++ LMDZ4/trunk/libf/dyn3d/limit_netcdf.F90	(revision 1319)
@@ -1,1334 +1,617 @@
 !
 ! $Id$
-!
-C
-C
-      SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque)
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque)
+!
+!-------------------------------------------------------------------------------
+! Author : L. Fairhead, 27/01/94
+!-------------------------------------------------------------------------------
+! Purpose: Boundary conditions files building for new model using climatologies.
+!          Both grids have to be regular.
+!-------------------------------------------------------------------------------
+! Note: This routine is designed to work for Earth
+!-------------------------------------------------------------------------------
+! Modification history:
+!  * 23/03/1994: Z. X. Li
+!  *    09/1999: L. Fairhead (netcdf reading in LMDZ.3.3)
+!  *    07/2001: P. Le Van
+!  *    11/2009: L. Guez     (ozone day & night climatos, see etat0_netcdf.F90)
+!  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
+!-------------------------------------------------------------------------------
 #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
+  USE dimphy
+  USE ioipsl,             ONLY : ioget_year_len
+  USE phys_state_var_mod, ONLY : pctsrf
+  USE netcdf,             ONLY : NF90_OPEN,    NF90_CREATE,  NF90_CLOSE,       &
+                   NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT,     &
+                   NF90_NOERR,   NF90_NOWRITE, NF90_DOUBLE,  NF90_GLOBAL,      &
+		   NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED
+#endif
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
 #include "dimensions.h"
 #include "paramet.h"
+#include "iniprint.h"
+  LOGICAL,                    INTENT(IN) :: interbar ! barycentric interpolation
+  LOGICAL,                    INTENT(IN) :: extrap   ! SST extrapolation flag
+  LOGICAL,                    INTENT(IN) :: oldice   ! old way ice computation
+  REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque   ! land mask
+#ifndef CPP_EARTH
+  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
+#else
+!-------------------------------------------------------------------------------
+! Local variables:
 #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"
+
+!--- For fractionary sub-cell use (old coding used soil index: 0,1,2,3) --------
+  LOGICAL, PARAMETER :: fracterre=.TRUE.
+
+!--- INPUT NETCDF FILES NAMES --------------------------------------------------
+  CHARACTER(LEN=25) :: icefile, sstfile, dumstr
+  CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc        ',        &
+                                  famipsic='amipbc_sic_1x1.nc        ',        &
+                                  fclimsst='amipbc_sst_1x1_clim.nc   ',        &
+                                  fclimsic='amipbc_sic_1x1_clim.nc   ',        &
+                                  fcpldsst='cpl_atm_sst.nc           ',        &
+                                  fcpldsic='cpl_atm_sic.nc           ',        &
+                                  frugo   ='Rugos.nc                 ',        &
+                                  falbe   ='Albedo.nc                '
+
+!--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
+  REAL,   DIMENSION(klon)                :: fi_ice, verif
+  REAL,   DIMENSION(:,:),   POINTER      :: phy_rug=>NULL(), phy_ice=>NULL()
+  REAL,   DIMENSION(:,:),   POINTER      :: phy_sst=>NULL(), phy_alb=>NULL()
+  REAL,   DIMENSION(:,:),   ALLOCATABLE  :: phy_bil
+  REAL,   DIMENSION(:,:,:), ALLOCATABLE  :: pctsrf_t
+  INTEGER                                :: nbad
+
+!--- VARIABLES FOR OUTPUT FILE WRITING -----------------------------------------
+  INTEGER :: ierr, nid, ndim, ntim, k
+  INTEGER, DIMENSION(2) :: dims
+  INTEGER :: id_tim,  id_SST,  id_BILS, id_RUG, id_ALB
+  INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC
+  INTEGER :: NF90_FORMAT
+  LOGICAL :: lCPL                    !--- T: IPCC-IPSL cpl model output files
+
+!--- MICS (PHYSICAL KEYS, TIME) ------------------------------------------------
+!  INTEGER, PARAMETER        :: longcles=20
+!  REAL, DIMENSION(longcles) :: clesphy0
+  INTEGER                   :: ndays
+  
+!--- INITIALIZATIONS -----------------------------------------------------------
+#ifdef NC_DOUBLE
+  NF90_FORMAT=NF90_DOUBLE
+#else
+  NF90_FORMAT=NF90_FLOAT
+#endif
+
+!  CALL conf_gcm(99, .TRUE., clesphy0)
+  pi    = 4.*ATAN(1.)
+  rad   = 6371229.
+  daysec= 86400.
+  omeg  = 2.*pi/daysec
+  g     = 9.8
+  kappa = 0.2857143
+  cpp   = 1004.70885
+  dtvr  = daysec/FLOAT(day_step)
+  CALL inigeom
+
+!--- Beware: anneeref (from gcm.def) is used to determine output time sampling
+  ndays=ioget_year_len(anneeref)
+
+!--- RUGOSITY TREATMENT --------------------------------------------------------
+  WRITE(lunout,*) 'Traitement de la rugosite'
+  CALL get_2Dfield(frugo,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
+
+!--- OCEAN TREATMENT -----------------------------------------------------------
+  PRINT*, 'Traitement de la glace oceanique' ; icefile=''; lCPL=.FALSE.
+
+! Input SIC file selection
+  icefile='fake'
+  IF(NF90_OPEN(famipsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(famipsic)
+  IF(NF90_OPEN(fclimsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fclimsic)
+  IF(NF90_OPEN(fcpldsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fcpldsic)
+  SELECT CASE(icefile)
+    CASE(famipsic); dumstr='Amip.'
+    CASE(fclimsic); dumstr='Amip climatologique.'
+    CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE.
+    CASE('fake');   CALL abort_gcm('limit_netcdf','Fichier SIC non reconnu.',1)
+  END SELECT
+  ierr=NF90_CLOSE(nid)
+  WRITE(lunout,*)'Pour la glace de mer a ete choisi un fichier '//TRIM(dumstr)
+
+  CALL get_2Dfield(icefile,'SIC',interbar,ndays,phy_ice,flag=oldice)
+
+  ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
+  DO k=1,ndays
+    fi_ice=phy_ice(:,k)
+    WHERE(fi_ice>=1.0  ) fi_ice=1.0
+    WHERE(fi_ice<EPSFRA) fi_ice=0.0
+    IF(fracterre) THEN
+      pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
+      pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
+      IF(lCPL) THEN                              ! SIC=pICE*(1-LIC-TER)
+        pctsrf_t(:,is_sic,k)=fi_ice*(1-pctsrf(:,is_lic)-pctsrf(:,is_ter))
+      ELSE                                        ! SIC=pICE-LIC
+        pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
+      END IF
+      WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
+      WHERE(1.0-zmasq<EPSFRA)
+        pctsrf_t(:,is_sic,k)=0.0
+        pctsrf_t(:,is_oce,k)=0.0
+      ELSEWHERE
+        WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
+          pctsrf_t(:,is_sic,k)=1.0-zmasq
+          pctsrf_t(:,is_oce,k)=0.0
+        ELSEWHERE
+          pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
+          WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
+            pctsrf_t(:,is_oce,k)=0.0
+            pctsrf_t(:,is_sic,k)=1.0-zmasq
+          END WHERE
+        END WHERE
+      END WHERE
+      nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
+      IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
+      nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
+      IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
+    ELSE 
+      pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)
+      WHERE(NINT(pctsrf(:,is_ter))==1)
+        pctsrf_t(:,is_sic,k)=0.
+        pctsrf_t(:,is_oce,k)=0.                  
+        WHERE(fi_ice>=1.e-5)
+          pctsrf_t(:,is_lic,k)=fi_ice
+          pctsrf_t(:,is_ter,k)=pctsrf_t(:,is_ter,k)-pctsrf_t(:,is_lic,k)
+        ELSEWHERE
+          pctsrf_t(:,is_lic,k)=0.0
+        END WHERE
+      ELSEWHERE
+        pctsrf_t(:,is_lic,k) = 0.0 
+        WHERE(fi_ice>=1.e-5)
+          pctsrf_t(:,is_ter,k)=0.0
+          pctsrf_t(:,is_sic,k)=fi_ice
+          pctsrf_t(:,is_oce,k)=1.0-pctsrf_t(:,is_sic,k)
+        ELSEWHERE
+          pctsrf_t(:,is_sic,k)=0.0
+          pctsrf_t(:,is_oce,k)=1.0
+        END WHERE
+      END WHERE
+      verif=sum(pctsrf_t(:,:,k),dim=2)
+      nbad=COUNT(verif<1.0-1.e-5.OR.verif>1.0+1.e-5)
+      IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
+    END IF 
+  END DO
+  DEALLOCATE(phy_ice)
+
+!--- SST TREATMENT -------------------------------------------------------------
+  WRITE(lunout,*) 'Traitement de la sst' ; sstfile=''; lCPL=.FALSE.
+
+! Input SST file selection
+  sstfile='fake'
+  IF(NF90_OPEN(famipsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(famipsst)
+  IF(NF90_OPEN(fclimsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fclimsst)
+  IF(NF90_OPEN(fcpldsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fcpldsst)
+  SELECT CASE(icefile)
+    CASE(famipsic); dumstr='Amip.'
+    CASE(fclimsic); dumstr='Amip climatologique.'
+    CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE.
+    CASE('fake');   CALL abort_gcm('limit_netcdf','Fichier SST non reconnu',1)
+  END SELECT
+  ierr=NF90_CLOSE(nid)
+  WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(dumstr)
+
+  CALL get_2Dfield(trim(sstfile),'SST',interbar,ndays,phy_sst,flag=extrap)
+
+!--- ALBEDO TREATMENT ----------------------------------------------------------
+  WRITE(lunout,*) 'Traitement de l albedo'
+  CALL get_2Dfield(falbe,'ALB',interbar,ndays,phy_alb)
+
+!--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
+  ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0
+
+!--- OUTPUT FILE WRITING -------------------------------------------------------
+  WRITE(lunout,*) 'Ecriture du fichier limit'
+
+  !--- File creation
+  ierr=NF90_CREATE("limit.nc",NF90_CLOBBER,nid)
+  ierr=NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier conditions aux limites")
+
+  !--- Dimensions creation
+  ierr=NF90_DEF_DIM(nid,"points_physiques",klon,ndim)
+  ierr=NF90_DEF_DIM(nid,"time",NF90_UNLIMITED,ntim)
+
+  dims=(/ndim,ntim/)
+
+  !--- Variables creation
+  ierr=NF90_DEF_VAR(nid,"TEMPS",NF90_FORMAT,dims,id_tim)
+  ierr=NF90_DEF_VAR(nid,"FOCE", NF90_FORMAT,dims,id_FOCE)
+  ierr=NF90_DEF_VAR(nid,"FSIC", NF90_FORMAT,dims,id_FSIC)
+  ierr=NF90_DEF_VAR(nid,"FTER", NF90_FORMAT,dims,id_FTER)
+  ierr=NF90_DEF_VAR(nid,"FLIC", NF90_FORMAT,dims,id_FLIC)
+  ierr=NF90_DEF_VAR(nid,"SST",  NF90_FORMAT,dims,id_SST)
+  ierr=NF90_DEF_VAR(nid,"BILS", NF90_FORMAT,dims,id_BILS)
+  ierr=NF90_DEF_VAR(nid,"ALB",  NF90_FORMAT,dims,id_ALB)
+  ierr=NF90_DEF_VAR(nid,"RUG",  NF90_FORMAT,dims,id_RUG)
+
+  !--- Attributes creation
+  ierr=NF90_PUT_ATT(nid,id_tim, "title","Jour dans l annee")
+  ierr=NF90_PUT_ATT(nid,id_FOCE,"title","Fraction ocean")
+  ierr=NF90_PUT_ATT(nid,id_FSIC,"title","Fraction glace de mer")
+  ierr=NF90_PUT_ATT(nid,id_FTER,"title","Fraction terre")
+  ierr=NF90_PUT_ATT(nid,id_FLIC,"title","Fraction land ice")
+  ierr=NF90_PUT_ATT(nid,id_SST ,"title","Temperature superficielle de la mer")
+  ierr=NF90_PUT_ATT(nid,id_BILS,"title","Reference flux de chaleur au sol")
+  ierr=NF90_PUT_ATT(nid,id_ALB, "title","Albedo a la surface")
+  ierr=NF90_PUT_ATT(nid,id_RUG, "title","Rugosite")
+
+  ierr=NF90_ENDDEF(nid)
+
+  !--- Variables saving
+  ierr=NF90_PUT_VAR(nid,id_tim,(/(DBLE(k),k=1,ndays)/))
+  ierr=NF90_PUT_VAR(nid,id_FOCE,pctsrf_t(:,is_oce,:),(/1,1/),(/klon,ndays/))
+  ierr=NF90_PUT_VAR(nid,id_FSIC,pctsrf_t(:,is_sic,:),(/1,1/),(/klon,ndays/))
+  ierr=NF90_PUT_VAR(nid,id_FTER,pctsrf_t(:,is_ter,:),(/1,1/),(/klon,ndays/))
+  ierr=NF90_PUT_VAR(nid,id_FLIC,pctsrf_t(:,is_lic,:),(/1,1/),(/klon,ndays/))
+  ierr=NF90_PUT_VAR(nid,id_SST ,phy_sst(:,:),(/1,1/),(/klon,ndays/))
+  ierr=NF90_PUT_VAR(nid,id_BILS,phy_bil(:,:),(/1,1/),(/klon,ndays/))
+  ierr=NF90_PUT_VAR(nid,id_ALB ,phy_alb(:,:),(/1,1/),(/klon,ndays/))
+  ierr=NF90_PUT_VAR(nid,id_RUG ,phy_rug(:,:),(/1,1/),(/klon,ndays/))
+
+  ierr=NF90_CLOSE(nid)
+
+  DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
+#endif
+! of #ifdef CPP_EARTH
+  STOP
+
+
+!===============================================================================
+!
+  CONTAINS
+!
+!===============================================================================
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask)
+!
+!-------------------------------------------------------------------------------
+! Comments:
+!   There are two assumptions concerning the NetCDF files, that are satisfied
+!   with files that are conforming NC convention:
+!     1) The last dimension of the variables used is the time record.
+!     2) Dimensional variables have the same names as corresponding dimensions.
+!-------------------------------------------------------------------------------
+  USE netcdf, ONLY : NF90_OPEN,    NF90_INQ_VARID, NF90_INQUIRE_VARIABLE,      &
+                     NF90_CLOSE,   NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION,     &
+                     NF90_GET_VAR, NF90_GET_ATT
+  USE dimphy, ONLY : klon
+  USE phys_state_var_mod, ONLY : pctsrf
+  IMPLICIT NONE
+#include "dimensions.h"
+#include "paramet.h"
+#include "comgeom2.h"
+#include "control.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
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
+  CHARACTER(LEN=3),  INTENT(IN)     :: mode     ! RUG, SIC, SST or ALB
+  LOGICAL,           INTENT(IN)     :: ibar     ! interp on pressure levels
+  INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
+  REAL,    POINTER,  DIMENSION(:,:) :: champo   ! output field = f(t)
+  LOGICAL, OPTIONAL,                      INTENT(IN) :: flag
+! flag=T means:   extrapolation (SST case)  or  old ice (SIC case)
+  REAL,    OPTIONAL, DIMENSION(iim,jjp1), INTENT(IN) :: mask
+!-------------------------------------------------------------------------------
+! Local variables:
+!--- NetCDF
+  INTEGER :: ierr, ncid, varid                  ! NetCDF identifiers
+  CHARACTER(LEN=30)               :: dnam       ! dimension name
+  CHARACTER(LEN=80)               :: varname    ! NetCDF variable name
+!--- dimensions
+  INTEGER,           DIMENSION(4) :: dids       ! NetCDF dimensions identifiers
+  REAL, ALLOCATABLE, DIMENSION(:) :: dlon_ini   ! initial longitudes vector
+  REAL, ALLOCATABLE, DIMENSION(:) :: dlat_ini   ! initial latitudes  vector
+  REAL, POINTER,     DIMENSION(:) :: dlon, dlat ! reordered lon/lat  vectors
+!--- fields
+  INTEGER :: imdep, jmdep, lmdep                ! dimensions of 'champ'
+  REAL, ALLOCATABLE, DIMENSION(:) :: champ      ! wanted field on initial grid
+  REAL, ALLOCATABLE, DIMENSION(:) :: yder, timeyear
+  REAL,              DIMENSION(iim,jjp1) :: champint   ! interpolated field
+  REAL, ALLOCATABLE, DIMENSION(:,:,:)    :: champtime
+  REAL, ALLOCATABLE, DIMENSION(:,:,:)    :: champan
+  REAL                              :: time, by
+!--- input files
+  CHARACTER(LEN=20)                 :: cal_in   ! calendar
+  INTEGER                           :: ndays_in ! number of days
+!--- misc
+  INTEGER :: i, j, k, l                         ! loop counters
+  REAL, ALLOCATABLE, DIMENSION(:,:) :: work     ! used for extrapolation
+  CHARACTER(LEN=25)                 :: title    ! for messages
+  LOGICAL                           :: extrp    ! flag for extrapolation
+  REAL                              :: chmin, chmax
+!-------------------------------------------------------------------------------
+!---Variables depending on keyword 'mode' --------------------------------------
+  NULLIFY(champo)
+  SELECT CASE(mode)
+    CASE('RUG'); varname='RUGOS';  title='Rugosite'
+    CASE('SIC'); varname='sicbcs'; title='Sea-ice'
+    CASE('SST'); varname='tosbcs'; title='SST'
+    CASE('ALB'); varname='ALBEDO'; title='Albedo'
+  END SELECT
+  extrp=(flag.AND.mode=='SST')
+
+!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE ------------------------------
+  ierr=NF90_OPEN(fnam,NF90_NOWRITE,ncid);                  CALL ncerr(ierr,fnam)
+  ierr=NF90_INQ_VARID(ncid,varname,varid);                 CALL ncerr(ierr,fnam)
+  ierr=NF90_INQUIRE_VARIABLE(ncid,varid,dimids=dids);      CALL ncerr(ierr,fnam)
+
+!--- Longitude
+  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(1),name=dnam,len=imdep)
+  CALL ncerr(ierr,fnam); ALLOCATE(dlon_ini(imdep),dlon(imdep))
+  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
+  ierr=NF90_GET_VAR(ncid,varid,dlon_ini);                  CALL ncerr(ierr,fnam)
+  WRITE(lunout,*) 'variable ', dnam, 'dimension ', imdep
+
+!--- Latitude
+  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(2),name=dnam,len=jmdep)
+  CALL ncerr(ierr,fnam); ALLOCATE(dlat_ini(jmdep),dlat(jmdep))
+  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
+  ierr=NF90_GET_VAR(ncid,varid,dlat_ini);                  CALL ncerr(ierr,fnam)
+  WRITE(lunout,*) 'variable ', dnam, 'dimension ', jmdep
+
+!--- Time (variable is not needed - it is rebuilt - but calendar is)
+  ierr=NF90_INQUIRE_DIMENSION(ncid,dids(3),name=dnam,len=lmdep)
+  CALL ncerr(ierr,fnam); ALLOCATE(timeyear(lmdep))
+  ierr=NF90_INQ_VARID(ncid,dnam,varid);                    CALL ncerr(ierr,fnam)
+  cal_in=' '
+  ierr=NF90_GET_ATT(ncid,varid,'calendar',cal_in)
+  IF(ierr/=NF90_NOERR) THEN
+    SELECT CASE(mode)
+      CASE('RUG','ALB'); cal_in='360d'
+      CASE('SIC','SST'); cal_in='gregorian'
+    END SELECT
+    WRITE(lunout,*)'ATTENTION: variable ''time'' sans attribut ''calendrier'' d&
+     &ans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
+  END IF
+  WRITE(lunout,*) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', cal_in
+
+!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION ----------------------
+  !--- Determining input file number of days, depending on calendar
+  ndays_in=year_len(anneeref,cal_in)
+
+!--- Time vector reconstruction (time vector from file is not trusted)
+!--- If input records are not monthly, time sampling has to be constant !
+  timeyear=mid_months(anneeref,cal_in,lmdep)
+  IF(lmdep/=12) WRITE(lunout,'(a,i3,a)')'Note: les fichiers de '//TRIM(mode)   &
+    //' ne comportent pas 12, mais ',lmdep,' renregistrements.'
+
+!--- GETTING THE FIELD AND INTERPOLATING IT ------------------------------------
+  ALLOCATE(champ(imdep*jmdep),champtime(iim,jjp1,lmdep))
+  IF(extrp) ALLOCATE(work(imdep,jmdep))
+
+  WRITE(lunout,*)
+  WRITE(lunout,'(a,i3,a)')'LECTURE ET INTERPOLATION HORIZ. DE ',lmdep,' CHAMPS.'
+  ierr=NF90_INQ_VARID(ncid,varname,varid);                 CALL ncerr(ierr,fnam)
+  DO l=1,lmdep
+    ierr=NF90_GET_VAR(ncid,varid,champ,(/1,1,l/),(/imdep,jmdep,1/))
+    CALL ncerr(ierr,fnam)
+    CALL conf_dat2d(title,imdep,jmdep,dlon_ini,dlat_ini,dlon,dlat,champ,ibar)
+
+    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
+
+    IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN
+      IF(l==1) THEN
+        WRITE(lunout,*)                                                        &
+  '-------------------------------------------------------------------------'
+        WRITE(lunout,*)                                                        &
+  '$$$ Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
+        WRITE(lunout,*)                                                        &
+  '-------------------------------------------------------------------------'
+      END IF
+      IF(mode=='RUG') champ=LOG(champ)
+      CALL inter_barxy(imdep, jmdep-1, dlon, dlat, champ, iim, jjm, rlonu,     &
+                         rlatv, jjp1, champint)
+      IF(mode=='RUG') THEN
+        champint=EXP(champint)
+        WHERE(NINT(mask)/=1) champint=0.001
+      END IF
+    ELSE
+      SELECT CASE(mode)
+        CASE('RUG');       CALL rugosite(imdep, jmdep, dlon, dlat, champ,      &
+                                    iim, jjp1, rlonv, rlatu, champint, mask)
+        CASE('SIC');       CALL sea_ice (imdep, jmdep, dlon, dlat, champ,      &
+                                    iim, jjp1, rlonv, rlatu, champint)
+        CASE('SST','ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ,      &
+                                    iim, jjp1, rlonv, rlatu, champint)
+      END SELECT
+    END IF
+    champtime(:,:,l)=champint
+  END DO
+  ierr=NF90_CLOSE(ncid);                                   CALL ncerr(ierr,fnam)
+
+  DEALLOCATE(dlon_ini,dlat_ini,dlon,dlat,champ)
+  IF(extrp) DEALLOCATE(work)
+
+!--- TIME INTERPOLATION --------------------------------------------------------
+  WRITE(lunout,*)
+  WRITE(lunout,*)'INTERPOLATION TEMPORELLE.'
+  WRITE(lunout,"(2x,' Vecteur temps en entree: ',10f6.1)") timeyear
+  WRITE(lunout,"(2x,' Vecteur temps en sortie de 0 a ',i3)") ndays
+  ALLOCATE(yder(lmdep),champan(iip1,jjp1,ndays))
+  DO j=1,jjp1
+    DO i=1,iim
+      CALL spline(timeyear,champtime(i,j,:),lmdep,1.e30,1.e30,yder)
+      DO k=1,ndays
+        time=FLOAT((k-1)*ndays_in)/FLOAT(ndays)
+        CALL splint(timeyear,champtime(i,j,:),yder,lmdep,time,by)
+        champan(i,j,k) = by
+      END DO
+    END DO
+  END DO
+  champan(iip1,:,:)=champan(1,:,:)
+  DEALLOCATE(yder,champtime,timeyear)
+
+!--- Checking the result
+  DO j=1,jjp1
+    CALL minmax(iip1,champan(1,j,10),chmin,chmax)
+    WRITE(lunout,*)' '//TRIM(title)//' au temps 10 ',chmin,chmax,j
+  END DO
+
+!--- SPECIAL FILTER FOR SST: SST>271.38 ----------------------------------------
+  IF(mode=='SST') THEN
+    WRITE(lunout,*) 'Filtrage de la SST: SST >= 271.38'
+    WHERE(champan<271.38) champan=271.38
+  END IF
+
+!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 ---------------------------------------
+  IF(mode=='SIC') THEN
+    WRITE(lunout,*) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
+    champan(:,:,:)=champan(:,:,:)/100.
+    champan(iip1,:,:)=champan(1,:,:)
+    WHERE(champan>1.0) champan=1.0
+    WHERE(champan<0.0) champan=0.0
+  END IF
+
+write(*,*)'coin1'
+!--- DYNAMICAL TO PHYSICAL GRID ------------------------------------------------
+  ALLOCATE(champo(klon,ndays))
+  DO k=1,ndays
+    CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k),champo(1,k))
+  END DO
+write(*,*)'coin2'
+  DEALLOCATE(champan)
+write(*,*)'coin3'
+END SUBROUTINE get_2Dfield
+!
+!-------------------------------------------------------------------------------
+
+
+
+!-------------------------------------------------------------------------------
+!
+FUNCTION year_len(y,cal_in)
+!
+!-------------------------------------------------------------------------------
+  USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  INTEGER                       :: year_len
+  INTEGER,           INTENT(IN) :: y
+  CHARACTER(LEN=*),  INTENT(IN) :: cal_in
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=20)             :: cal_out              ! calendar (for outputs)
+!-------------------------------------------------------------------------------
+!--- Getting the input calendar to reset at the end of the function
+  CALL ioget_calendar(cal_out)
+
+!--- Unlocking calendar and setting it to wanted one
+  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
+
+!--- Getting the number of days in this year
+  year_len=ioget_year_len(y)
+
+!--- Back to original calendar
+  CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
+
+END FUNCTION year_len
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+FUNCTION mid_months(y,cal_in,nm)
+!
+!-------------------------------------------------------------------------------
+  USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  INTEGER,                INTENT(IN) :: y               ! year
+  CHARACTER(LEN=*),       INTENT(IN) :: cal_in          ! calendar
+  INTEGER,                INTENT(IN) :: nm              ! months/year number
+  REAL,    DIMENSION(nm)             :: mid_months      ! mid-month times
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=99)                  :: mess            ! error message
+  CHARACTER(LEN=20)                  :: cal_out         ! calendar (for outputs)
+  INTEGER, DIMENSION(nm)             :: mnth            ! months lengths (days)
+  INTEGER                            :: m               ! months counter
+  INTEGER                            :: nd              ! number of days
+!-------------------------------------------------------------------------------
+  nd=year_len(y,cal_in)
+
+  IF(nm==12) THEN
+
+  !--- Getting the input calendar to reset at the end of the function
+    CALL ioget_calendar(cal_out)
+
+  !--- Unlocking calendar and setting it to wanted one
+    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
+
+  !--- Getting the length of each month
+    DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
+
+  !--- Back to original calendar
+    CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
+
+  ELSE IF(MODULO(nd,nm)/=0) THEN
+    WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
+      nm,' months/year. Months number should divide days number.'
+    CALL abort_gcm('mid_months',TRIM(mess),1)
+
+  ELSE
+    mnth=(/(m,m=1,nm,nd/nm)/)
+  END IF
+
+!--- Mid-months times
+  mid_months(1)=0.5*FLOAT(mnth(1))
+  DO k=2,nm
+    mid_months(k)=mid_months(k-1)+0.5*FLOAT(mnth(k-1)+mnth(k))
+  END DO
+
+END FUNCTION mid_months
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE ncerr(ncres,fnam)
+!
+!-------------------------------------------------------------------------------
+! Purpose: NetCDF errors handling.
+!-------------------------------------------------------------------------------
+  USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  INTEGER,          INTENT(IN) :: ncres
+  CHARACTER(LEN=*), INTENT(IN) :: fnam
+!-------------------------------------------------------------------------------
+#include "iniprint.h"
+  IF(ncres/=NF90_NOERR) THEN
+    WRITE(lunout,*)'Problem with file '//TRIM(fnam)//' in routine limit_netcdf.'
+    CALL abort_gcm('limit_netcdf',NF90_STRERROR(ncres),1)
+  END IF
+
+END SUBROUTINE ncerr
+!
+!-------------------------------------------------------------------------------
+
+
+END SUBROUTINE limit_netcdf
