Index: LMDZ5/trunk/libf/phylmd/limit_netcdf.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/limit_netcdf.F90	(revision 2239)
+++ LMDZ5/trunk/libf/phylmd/limit_netcdf.F90	(revision 2293)
@@ -1,8 +1,3 @@
-!
-! $Id: limit_netcdf.F90 1508 2011-04-15 13:05:34Z jghattas $
-!-------------------------------------------------------------------------------
-!
 SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque)
-#ifndef CPP_1D
 !
 !-------------------------------------------------------------------------------
@@ -21,4 +16,5 @@
 !  *    12/2009: D. Cugnet   (f77->f90, calendars, files from coupled runs)
 !-------------------------------------------------------------------------------
+#ifndef CPP_1D
   USE control_mod
   USE indice_sol_mod
@@ -32,12 +28,13 @@
 		   NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED, NF90_FLOAT
   USE inter_barxy_m, only: inter_barxy
-  use netcdf95, only: nf95_def_var, nf95_put_att, nf95_put_var
+  USE netcdf95,    ONLY: nf95_def_var, nf95_put_att, nf95_put_var
+  USE grid_atob_m, ONLY: grille_m, rugosite, sea_ice
 #endif
   IMPLICIT NONE
 !-------------------------------------------------------------------------------
 ! Arguments:
-#include "dimensions.h"
-#include "paramet.h"
-#include "iniprint.h"
+  include "dimensions.h"
+  include "paramet.h"
+  include "iniprint.h"
   LOGICAL,                    INTENT(IN) :: interbar ! barycentric interpolation
   LOGICAL,                    INTENT(IN) :: extrap   ! SST extrapolation flag
@@ -49,33 +46,33 @@
 !-------------------------------------------------------------------------------
 ! Local variables:
-#include "logic.h"
-#include "comvert.h"
-#include "comgeom2.h"
-#include "comconst.h"
+  include "logic.h"
+  include "comvert.h"
+  include "comgeom2.h"
+  include "comconst.h"
 
 !--- INPUT NETCDF FILES NAMES --------------------------------------------------
-  CHARACTER(LEN=25) :: icefile, sstfile, dumstr
-  CHARACTER(LEN=25), PARAMETER :: famipsst ='amipbc_sst_1x1.nc        ',        &
-                                  famipsic ='amipbc_sic_1x1.nc        ',        &
-                                  fcpldsst ='cpl_atm_sst.nc           ',        &
-                                  fcpldsic ='cpl_atm_sic.nc           ',        &
-                                  fhistsst ='histmth_sst.nc           ',        &
-                                  fhistsic ='histmth_sic.nc           ',        &
-                                  frugo    ='Rugos.nc                 ',        &
-                                  falbe    ='Albedo.nc                ',        &
-                                  feraisstk='sstk.nc                  ',        &
-                                  feraici  ='ci.nc                    '
+  CHARACTER(LEN=20) :: icefile, sstfile, dumstr, fnam
+  CHARACTER(LEN=20), PARAMETER :: &
+  fsst(4)=['amipbc_sst_1x1.nc   ','cpl_atm_sst.nc      ','histmth_sst.nc      '&
+          ,'sstk.nc             ']
+  CHARACTER(LEN=20), PARAMETER :: &
+  fsic(4)=['amipbc_sic_1x1.nc   ','cpl_atm_sic.nc      ','histmth_sic.nc      '&
+          ,'ci.nc               ']
+  CHARACTER(LEN=10), PARAMETER :: &
+  vsst(4)=['tosbcs    ','SISUTESW  ','tsol_oce  ','sstk      '], &
+  vsic(4)=['sicbcs    ','SIICECOV  ','pourc_sic ','ci        ']
+  CHARACTER(LEN=20), PARAMETER :: frugo='Rugos.nc            ',                &
+                                  falbe='Albedo.nc           '
+  CHARACTER(LEN=10), PARAMETER :: vrug='RUGOS     ', valb='ALBEDO    '
   CHARACTER(LEN=10) :: varname
 !--- 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
+  REAL               :: fi_ice(klon), verif(klon)
+  REAL, POINTER      :: phy_rug(:,:)=>NULL(), phy_ice(:,:)=>NULL()
+  REAL, POINTER      :: phy_sst(:,:)=>NULL(), phy_alb(:,:)=>NULL()
+  REAL, ALLOCATABLE  :: phy_bil(:,:), pctsrf_t(:,:,:)
+  INTEGER            :: nbad
 
 !--- VARIABLES FOR OUTPUT FILE WRITING -----------------------------------------
-  INTEGER :: ierr, nid, ndim, ntim, k
-  INTEGER, DIMENSION(2) :: dims
+  INTEGER :: ierr, nid, ndim, ntim, k, dims(2), ix_sic, ix_sst
   INTEGER :: id_tim,  id_SST,  id_BILS, id_RUG, id_ALB
   INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC, varid_longitude, varid_latitude
@@ -89,13 +86,4 @@
   NF90_FORMAT=NF90_FLOAT
 #endif
-
-  pi    = 4.*ATAN(1.)
-  rad   = 6371229.
-  daysec= 86400.
-  omeg  = 2.*pi/daysec
-  g     = 9.8
-  kappa = 0.2857143
-  cpp   = 1004.70885
-  dtvr  = daysec/REAL(day_step)
   CALL inigeom
 
@@ -104,33 +92,25 @@
 
 !--- RUGOSITY TREATMENT --------------------------------------------------------
-  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la rugosite'
-  varname='RUGOS'
-  CALL get_2Dfield(frugo,varname,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
+  CALL msg(1,'Traitement de la rugosite')
+  CALL get_2Dfield(frugo,vrug,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
 
 !--- OCEAN TREATMENT -----------------------------------------------------------
-  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la glace oceanique'
+  CALL msg(1,'Traitement de la glace oceanique')
 
 ! Input SIC file selection
 ! Open file only to test if available
-  IF ( NF90_OPEN(TRIM(famipsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     icefile=TRIM(famipsic)
-     varname='sicbcs'
-  ELSE IF( NF90_OPEN(TRIM(fcpldsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     icefile=TRIM(fcpldsic)
-     varname='SIICECOV'
-  ELSE IF ( NF90_OPEN(TRIM(fhistsic),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     icefile=TRIM(fhistsic)
-     varname='pourc_sic'
-  ELSE IF ( NF90_OPEN(TRIM(feraici),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     icefile=TRIM(feraici)
-     varname='ci'
-  ELSE
+  DO ix_sic=1,SIZE(fsic)
+     IF ( NF90_OPEN(TRIM(fsic(ix_sic)),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+        icefile=fsic(ix_sic); varname=vsic(ix_sic); EXIT
+     END IF
+  END DO
+  IF(ix_sic==SIZE(fsic)+1) THEN
      WRITE(lunout,*) 'ERROR! No sea-ice input file was found.'
-     WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsic),', ',trim(fcpldsic),', ', &
-                      trim(fhistsic), trim(feraici)
+     WRITE(lunout,*) 'One of following files must be available : '
+     DO k=1,SIZE(fsic); WRITE(lunout,*) TRIM(fsic(k)); END DO
      CALL abort_gcm('limit_netcdf','No sea-ice file was found',1)
   END IF
-  ierr=NF90_CLOSE(nid)
-  IF (prt_level>=0) WRITE(lunout,*)'Pour la glace de mer a ete choisi le fichier '//TRIM(icefile)
+  CALL ncerr(NF90_CLOSE(nid),icefile)
+  CALL msg(-1,'Fichier choisi pour la glace de mer:'//TRIM(icefile))
 
   CALL get_2Dfield(icefile,varname, 'SIC',interbar,ndays,phy_ice,flag=oldice)
@@ -143,11 +123,12 @@
      pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
      pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
-     IF (icefile==trim(fcpldsic)) THEN           ! SIC=pICE*(1-LIC-TER)
+     SELECT CASE(ix_sic)
+        CASE(2)                                   ! SIC=pICE*(1-LIC-TER) (CPL)
         pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
-     ELSE IF (icefile==trim(fhistsic)) THEN      ! SIC=pICE
+        CASE(3)                                   ! SIC=pICE            (HIST)
         pctsrf_t(:,is_sic,k)=fi_ice(:)
-     ELSE ! icefile==famipsic                    ! SIC=pICE-LIC
+        CASE DEFAULT                              ! SIC=pICE-LIC   (AMIP,ERAI)
         pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
-     END IF
+     END SELECT
      WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
      WHERE(1.0-zmasq<EPSFRA)
@@ -167,6 +148,6 @@
      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 maille pour nb points = ',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
   END DO
@@ -174,34 +155,27 @@
 
 !--- SST TREATMENT -------------------------------------------------------------
-  IF (prt_level>1) WRITE(lunout,*) 'Traitement de la sst'
+  CALL msg(1,'Traitement de la sst')
 
 ! Input SST file selection
 ! Open file only to test if available
-  IF ( NF90_OPEN(TRIM(famipsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     sstfile=TRIM(famipsst)
-     varname='tosbcs'
-  ELSE IF ( NF90_OPEN(TRIM(fcpldsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     sstfile=TRIM(fcpldsst)
-     varname='SISUTESW'
-  ELSE IF ( NF90_OPEN(TRIM(fhistsst),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     sstfile=TRIM(fhistsst)
-     varname='tsol_oce'
-  ELSE IF ( NF90_OPEN(TRIM(feraisstk),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
-     sstfile=TRIM(feraisstk)
-     varname='sstk'
-  ELSE
+  DO ix_sst=1,SIZE(fsst)
+     IF ( NF90_OPEN(TRIM(fsst(ix_sst)),NF90_NOWRITE,nid)==NF90_NOERR ) THEN
+       sstfile=fsst(ix_sst); varname=vsst(ix_sst); EXIT
+     END IF
+  END DO
+  IF(ix_sst==SIZE(fsst)+1) THEN
      WRITE(lunout,*) 'ERROR! No sst input file was found.'
-     WRITE(lunout,*) 'One of following files must be availible : ',trim(famipsst),trim(fcpldsst),trim(fhistsst),trim(feraisstk)
+     WRITE(lunout,*) 'One of following files must be available : '
+     DO k=1,SIZE(fsst); WRITE(lunout,*) TRIM(fsst(k)); END DO
      CALL abort_gcm('limit_netcdf','No sst file was found',1)
   END IF
-  ierr=NF90_CLOSE(nid)
-  IF (prt_level>=0) WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(sstfile)
+  CALL ncerr(NF90_CLOSE(nid),sstfile)
+  CALL msg(-1,'Fichier choisi pour la temperature de mer: '//TRIM(sstfile))
 
   CALL get_2Dfield(sstfile,varname,'SST',interbar,ndays,phy_sst,flag=extrap)
 
 !--- ALBEDO TREATMENT ----------------------------------------------------------
-  IF (prt_level>1) WRITE(lunout,*) 'Traitement de l albedo'
-  varname='ALBEDO'
-  CALL get_2Dfield(falbe,varname,'ALB',interbar,ndays,phy_alb)
+  CALL msg(1,'Traitement de l albedo')
+  CALL get_2Dfield(falbe,valb,'ALB',interbar,ndays,phy_alb)
 
 !--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
@@ -209,39 +183,40 @@
 
 !--- OUTPUT FILE WRITING -------------------------------------------------------
-  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : debut'
+  CALL msg(5,'Ecriture du fichier limit : debut')
+  fnam="limit.nc"
 
   !--- File creation
-  ierr=NF90_CREATE("limit.nc",NF90_CLOBBER,nid)
-  ierr=NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier conditions aux limites")
+  CALL ncerr(NF90_CREATE(fnam,NF90_CLOBBER,nid),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier conditions aux limites"),fnam)
 
   !--- Dimensions creation
-  ierr=NF90_DEF_DIM(nid,"points_physiques",klon,ndim)
-  ierr=NF90_DEF_DIM(nid,"time",NF90_UNLIMITED,ntim)
-
-  dims=(/ndim,ntim/)
+  CALL ncerr(NF90_DEF_DIM(nid,"points_physiques",klon,ndim),fnam)
+  CALL ncerr(NF90_DEF_DIM(nid,"time",NF90_UNLIMITED,ntim),fnam)
+
+  dims=[ndim,ntim]
 
   !--- Variables creation
-  ierr=NF90_DEF_VAR(nid,"TEMPS",NF90_FORMAT,(/ntim/),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)
+  CALL ncerr(NF90_DEF_VAR(nid,"TEMPS",NF90_FORMAT,[ntim],id_tim),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"FOCE", NF90_FORMAT,dims,id_FOCE),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"FSIC", NF90_FORMAT,dims,id_FSIC),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"FTER", NF90_FORMAT,dims,id_FTER),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"FLIC", NF90_FORMAT,dims,id_FLIC),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"SST",  NF90_FORMAT,dims,id_SST),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"BILS", NF90_FORMAT,dims,id_BILS),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"ALB",  NF90_FORMAT,dims,id_ALB),fnam)
+  CALL ncerr(NF90_DEF_VAR(nid,"RUG",  NF90_FORMAT,dims,id_RUG),fnam)
   call nf95_def_var(nid, "longitude", NF90_FLOAT, ndim, varid_longitude)
-  call nf95_def_var(nid, "latitude", NF90_FLOAT, ndim, varid_latitude)
+  call nf95_def_var(nid, "latitude",  NF90_FLOAT, ndim, varid_latitude)
 
   !--- 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")
+  CALL ncerr(NF90_PUT_ATT(nid,id_tim, "title","Jour dans l annee"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_FOCE,"title","Fraction ocean"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_FSIC,"title","Fraction glace de mer"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_FTER,"title","Fraction terre"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_FLIC,"title","Fraction land ice"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_SST ,"title","Temperature superficielle de la mer"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_BILS,"title","Reference flux de chaleur au sol"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_ALB, "title","Albedo a la surface"),fnam)
+  CALL ncerr(NF90_PUT_ATT(nid,id_RUG, "title","Rugosite"),fnam)
 
   call nf95_put_att(nid, varid_longitude, "standard_name", "longitude")
@@ -251,22 +226,22 @@
   call nf95_put_att(nid, varid_latitude, "units", "degrees_north")
 
-  ierr=NF90_ENDDEF(nid)
+  CALL ncerr(NF90_ENDDEF(nid),fnam)
 
   !--- Variables saving
-  ierr=NF90_PUT_VAR(nid,id_tim,(/(REAL(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/))
+  CALL ncerr(NF90_PUT_VAR(nid,id_tim,[(REAL(k),k=1,ndays)]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_FOCE,pctsrf_t(:,is_oce,:),[1,1],[klon,ndays]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_FSIC,pctsrf_t(:,is_sic,:),[1,1],[klon,ndays]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_FTER,pctsrf_t(:,is_ter,:),[1,1],[klon,ndays]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_FLIC,pctsrf_t(:,is_lic,:),[1,1],[klon,ndays]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_SST ,phy_sst(:,:),[1,1],[klon,ndays]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_BILS,phy_bil(:,:),[1,1],[klon,ndays]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_ALB ,phy_alb(:,:),[1,1],[klon,ndays]),fnam)
+  CALL ncerr(NF90_PUT_VAR(nid,id_RUG ,phy_rug(:,:),[1,1],[klon,ndays]),fnam)
   call nf95_put_var(nid, varid_longitude, rlon)
   call nf95_put_var(nid, varid_latitude, rlat)
 
-  ierr=NF90_CLOSE(nid)
-
-  IF (prt_level>5) WRITE(lunout,*) 'Ecriture du fichier limit : fin'
+  CALL ncerr(NF90_CLOSE(nid),fnam)
+
+  CALL msg(6,'Ecriture du fichier limit : fin')
 
   DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
@@ -296,4 +271,5 @@
   USE dimphy, ONLY : klon
   USE phys_state_var_mod, ONLY : pctsrf
+  USE conf_dat_m, ONLY: conf_dat2d
   USE control_mod
   USE pchsp_95_m, only: pchsp_95
@@ -303,8 +279,8 @@
 
   IMPLICIT NONE
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom2.h"
-#include "iniprint.h"
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom2.h"
+  include "iniprint.h"
 !-----------------------------------------------------------------------------
 ! Arguments:
@@ -320,29 +296,29 @@
 ! Local variables:
 !--- NetCDF
-  INTEGER :: ncid, varid                  ! NetCDF identifiers
-  CHARACTER(LEN=30)               :: dnam       ! dimension name
+  INTEGER           :: ncid, varid        ! NetCDF identifiers
+  CHARACTER(LEN=30) :: dnam               ! dimension 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
+  INTEGER           :: dids(4)            ! NetCDF dimensions identifiers
+  REAL, ALLOCATABLE :: dlon_ini(:)        ! initial longitudes vector
+  REAL, ALLOCATABLE :: dlat_ini(:)        ! initial latitudes  vector
+  REAL, POINTER     :: 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
+  INTEGER :: imdep, jmdep, lmdep          ! dimensions of 'champ'
+  REAL, ALLOCATABLE :: champ(:,:)         ! wanted field on initial grid
+  REAL, ALLOCATABLE :: yder(:), timeyear(:)
+  REAL              :: champint(iim,jjp1) ! interpolated field
+  REAL, ALLOCATABLE :: champtime(:,:,:)
+  REAL, ALLOCATABLE :: champan(:,:,:)
 !--- input files
-  CHARACTER(LEN=20)                 :: cal_in   ! calendar
-  CHARACTER(LEN=20)                 :: unit_sic ! attribute unit in sea-ice file
-  INTEGER                           :: ndays_in ! number of days
+  CHARACTER(LEN=20) :: cal_in             ! calendar
+  CHARACTER(LEN=20) :: unit_sic           ! attribute unit in sea-ice file
+  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
-  LOGICAL                           :: oldice   ! flag for old way ice computation 
-  REAL                              :: chmin, chmax
+  INTEGER           :: i, j, k, l         ! loop counters
+  REAL, ALLOCATABLE :: work(:,:)          ! used for extrapolation
+  CHARACTER(LEN=25) :: title              ! for messages
+  LOGICAL           :: extrp              ! flag for extrapolation
+  LOGICAL           :: oldice             ! flag for old way ice computation 
+  REAL              :: chmin, chmax
   INTEGER ierr
   integer n_extrap ! number of extrapolated points
@@ -369,51 +345,47 @@
 
 !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
-  IF (prt_level>5) WRITE(lunout,*) ' Now reading file : ',fnam
-  ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid);             CALL ncerr(ierr, fnam)
-  ierr=NF90_INQ_VARID(ncid, trim(varname), varid);            CALL ncerr(ierr, fnam)
-  ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam)
+  CALL msg(5,' Now reading file : '//TRIM(fnam))
+  CALL ncerr(NF90_OPEN(fnam, NF90_NOWRITE, ncid),fnam)
+  CALL ncerr(NF90_INQ_VARID(ncid, trim(varname), varid),fnam)
+  CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam)
 
 !--- Read unit for sea-ice variable only
   IF (mode=='SIC') THEN
-     ierr=NF90_GET_ATT(ncid, varid, 'units', unit_sic)
-     IF(ierr/=NF90_NOERR) THEN
-        IF (prt_level>5) WRITE(lunout,*) 'No unit was given in sea-ice file. Take percentage as default value'
-        unit_sic='X'
-     ELSE
-        IF (prt_level>5) WRITE(lunout,*) ' Sea-ice cover has unit=',unit_sic
-     END IF
+    IF(NF90_GET_ATT(ncid, varid, 'units', unit_sic)/=NF90_NOERR) THEN
+      CALL msg(5,'No unit in sea-ice file. Take percentage as default value')
+      unit_sic='X'
+    ELSE
+      CALL msg(5,'Sea-ice cover has unit='//TRIM(unit_sic))
+    END IF
   END IF
 
 !--- 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)
-  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
+  CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep),fnam)
+  ALLOCATE(dlon_ini(imdep), dlon(imdep))
+  CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
+  CALL ncerr(NF90_GET_VAR(ncid, varid, dlon_ini), fnam)
+  CALL msg(5,'variable '//TRIM(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)
-  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
+  CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep),fnam)
+  ALLOCATE(dlat_ini(jmdep), dlat(jmdep))
+  CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
+  CALL ncerr(NF90_GET_VAR(ncid, varid, dlat_ini), fnam)
+  CALL msg(5,'variable '//TRIM(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)
+  CALL ncerr(NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep), fnam)
+  ALLOCATE(timeyear(lmdep))
+  CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
   cal_in=' '
-  ierr=NF90_GET_ATT(ncid, varid, 'calendar', cal_in)
-  IF(ierr/=NF90_NOERR) THEN
+  IF(NF90_GET_ATT(ncid, varid, 'calendar', cal_in)/=NF90_NOERR) THEN
     SELECT CASE(mode)
       CASE('RUG', 'ALB'); cal_in='360d'
       CASE('SIC', 'SST'); cal_in='gregorian'
     END SELECT
-    IF (prt_level>5) WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
-         // 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
-  END IF
-  IF (prt_level>5) WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
-       cal_in
-
+  CALL msg(5,'WARNING: missing "calendar" attribute for "time" in '&
+     &//TRIM(fnam)//'. Choosing default value.')
+  END IF
+  CALL msg(5,'var, calendar, dim: '//TRIM(dnam)//' '//TRIM(cal_in), lmdep)
   
 !--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
@@ -430,26 +402,20 @@
   ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep))
   IF(extrp) ALLOCATE(work(imdep, jmdep))
-
-  IF (prt_level>5) WRITE(lunout, *)
-  IF (prt_level>5) WRITE(lunout,*)'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.'
-  ierr=NF90_INQ_VARID(ncid, varname, varid);             CALL ncerr(ierr, fnam)
+  CALL msg(5,'')
+  CALL msg(5,'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, ' CHAMPS.')
+  CALL ncerr(NF90_INQ_VARID(ncid, varname, varid), 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)
+    CALL ncerr(NF90_GET_VAR(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam)
+    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, ibar)
+    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
 
     IF(ibar .AND. .NOT.oldice) THEN
-      IF(l==1 .AND. prt_level>5) THEN
-        WRITE(lunout, *) '-------------------------------------------------------------------------'
-        WRITE(lunout, *) 'Utilisation de l''interpolation barycentrique pour ',TRIM(title),' $$$'
-        WRITE(lunout, *) '-------------------------------------------------------------------------'
+      IF(l==1) THEN
+      CALL msg(5,"----------------------------------------------------------")
+      CALL msg(5,"$$$ Interpolation barycentrique pour "//TRIM(title)//" $$$")
+      CALL msg(5,"----------------------------------------------------------")
       END IF
       IF(mode=='RUG') champ=LOG(champ)
-      CALL inter_barxy(dlon, dlat(:jmdep-1), champ, rlonu(:iim),     &
-                         rlatv, champint)
+      CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
       IF(mode=='RUG') THEN
         champint=EXP(champint)
@@ -458,15 +424,12 @@
     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)
+        CASE('RUG');       CALL rugosite(dlon,dlat,champ,rlonv,rlatu,champint,mask)
+        CASE('SIC');       CALL sea_ice (dlon,dlat,champ,rlonv,rlatu,champint)
+        CASE('SST','ALB'); CALL grille_m(dlon,dlat,champ,rlonv,rlatu,champint)
       END SELECT
     END IF
     champtime(:, :, l)=champint
   END DO
-  ierr=NF90_CLOSE(ncid);                                 CALL ncerr(ierr, fnam)
+  CALL ncerr(NF90_CLOSE(ncid), fnam)
 
   DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
@@ -474,5 +437,5 @@
 
 !--- TIME INTERPOLATION ------------------------------------------------------
-  IF (prt_level>5) THEN
+  IF(prt_level>5) THEN
      WRITE(lunout, *)
      WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
@@ -508,5 +471,5 @@
 !--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
   IF(mode=='SST') THEN
-    IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
+    CALL msg(5,'Filtrage de la SST: SST >= 271.38')
     WHERE(champan<271.38) champan=271.38
   END IF
@@ -514,13 +477,13 @@
 !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
   IF(mode=='SIC') THEN
-    IF (prt_level>5) WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
+    CALL msg(5,'Filtrage de la SIC: 0.0 < Sea-ice < 1.0')
 
     IF (unit_sic=='1') THEN
        ! Nothing to be done for sea-ice field is already in fraction of 1
        ! This is the case for sea-ice in file cpl_atm_sic.nc
-       IF (prt_level>5) WRITE(lunout,*) 'Sea-ice field already in fraction of 1'
+       CALL msg(5,'Sea-ice field already in fraction of 1')
     ELSE
        ! Convert sea ice from percentage to fraction of 1
-       IF (prt_level>5) WRITE(lunout,*) 'Transformt sea-ice field from percentage to fraction of 1.' 
+       CALL msg(5,'Transformt sea-ice field from percentage to fraction of 1.')
        champan(:, :, :)=champan(:, :, :)/100.
     END IF
@@ -541,5 +504,4 @@
 !
 !-------------------------------------------------------------------------------
-
 
 
@@ -620,5 +582,5 @@
 
   ELSE
-    mnth=(/(m,m=1,nm,nd/nm)/)
+    mnth=[(m,m=1,nm,nd/nm)]
   END IF
 
@@ -634,4 +596,31 @@
 
 
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE msg(lev,str1,i,str2)
+!
+!-------------------------------------------------------------------------------
+! Arguments:
+  INTEGER,                    INTENT(IN) :: lev
+  CHARACTER(LEN=*),           INTENT(IN) :: str1
+  INTEGER,          OPTIONAL, INTENT(IN) :: i
+  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: str2
+!-------------------------------------------------------------------------------
+  IF(prt_level>lev) THEN
+    IF(PRESENT(str2)) THEN
+      WRITE(lunout,*) TRIM(str1), i, TRIM(str2)
+    ELSE IF(PRESENT(i)) THEN
+      WRITE(lunout,*) TRIM(str1), i
+    ELSE
+      WRITE(lunout,*) TRIM(str1)
+    END IF
+  END IF
+
+END SUBROUTINE msg
+!
+!-------------------------------------------------------------------------------
+
+
 !-------------------------------------------------------------------------------
 !
