Index: /LMDZ6/branches/Amaury_dev/libf/dyn3d/dynetat0.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3d/dynetat0.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3d/dynetat0.F90	(revision 5100)
@@ -8,6 +8,6 @@
   USE infotrac,    ONLY: nqtot, tracers, niso, iqIsoPha, iH2O, isoName
   USE strings_mod, ONLY: maxlen, msg, strStack, real2str, int2str
-  USE netcdf,      ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQ_VARID, &
-                         NF90_CLOSE, nf90_get_var, nf90_noerr
+  USE netcdf,      ONLY: nf90_open,  nf90_nowrite, NF90_INQ_VARID, &
+                         nf90_close, nf90_get_var, nf90_noerr
   USE readTracFiles_mod, ONLY: new2oldH2O, newHNO3, oldHNO3, getKey
   USE control_mod, ONLY: planet_type
@@ -50,5 +50,5 @@
 !--- Initial state file opening
   var=fichnom
-  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
+  CALL err(nf90_open(var,nf90_nowrite,fID),"open",var)
   CALL get_var1("controle",tab_cntrl)
 
@@ -188,5 +188,5 @@
   END DO
 
-  CALL err(NF90_CLOSE(fID),"close",fichnom)
+  CALL err(nf90_close(fID),"close",fichnom)
   day_ini=day_ini+INT(time)
   time=time-INT(time)
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3d/dynredem.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3d/dynredem.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3d/dynredem.F90	(revision 5100)
@@ -9,7 +9,7 @@
   USE strings_mod, ONLY: maxlen
   USE infotrac, ONLY: nqtot, tracers
-  USE netcdf, ONLY: NF90_CREATE, NF90_DEF_DIM, NF90_INQ_VARID, NF90_GLOBAL,    &
-                    NF90_CLOSE,  NF90_PUT_ATT, NF90_UNLIMITED, NF90_CLOBBER,   &
-                    NF90_64BIT_OFFSET
+  USE netcdf, ONLY: nf90_create, nf90_def_dim, NF90_INQ_VARID, nf90_global,    &
+                    nf90_close,  nf90_put_att, nf90_unlimited, nf90_clobber,   &
+                    nf90_64bit_offset
   USE dynredem_mod, ONLY: cre_var, put_var1, put_var2, err, modname, fil
   USE comvert_mod,  ONLY: ap, bp, presnivs, pa, preff, nivsig, nivsigs
@@ -105,18 +105,18 @@
 
 !--- File creation
-  CALL err(NF90_CREATE(fichnom,IOR(NF90_CLOBBER,NF90_64BIT_OFFSET),nid))
+  CALL err(nf90_create(fichnom,IOR(nf90_clobber,nf90_64bit_offset),nid))
 
 !--- Some global attributes
-  CALL err(NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier demarrage dynamique"))
+  CALL err(nf90_put_att(nid,nf90_global,"title","Fichier demarrage dynamique"))
 
 !--- Dimensions
-  CALL err(NF90_DEF_DIM(nid,"index", length, indexID))
-  CALL err(NF90_DEF_DIM(nid,"rlonu", iip1,   rlonuID))
-  CALL err(NF90_DEF_DIM(nid,"rlatu", jjp1,   rlatuID))
-  CALL err(NF90_DEF_DIM(nid,"rlonv", iip1,   rlonvID))
-  CALL err(NF90_DEF_DIM(nid,"rlatv", jjm,    rlatvID))
-  CALL err(NF90_DEF_DIM(nid,"sigs",  llm,        sID))
-  CALL err(NF90_DEF_DIM(nid,"sig",   llmp1,    sigID))
-  CALL err(NF90_DEF_DIM(nid,"temps", NF90_UNLIMITED, timID))
+  CALL err(nf90_def_dim(nid,"index", length, indexID))
+  CALL err(nf90_def_dim(nid,"rlonu", iip1,   rlonuID))
+  CALL err(nf90_def_dim(nid,"rlatu", jjp1,   rlatuID))
+  CALL err(nf90_def_dim(nid,"rlonv", iip1,   rlonvID))
+  CALL err(nf90_def_dim(nid,"rlatv", jjm,    rlatvID))
+  CALL err(nf90_def_dim(nid,"sigs",  llm,        sID))
+  CALL err(nf90_def_dim(nid,"sig",   llmp1,    sigID))
+  CALL err(nf90_def_dim(nid,"temps", nf90_unlimited, timID))
 
 !--- Define and save invariant fields
@@ -149,5 +149,5 @@
   CALL cre_var(nid,"masse","Masse d air"    ,[rlonvID,rlatuID,sID,timID])
   CALL cre_var(nid,"ps"   ,"Pression au sol",[rlonvID,rlatuID    ,timID])
-  CALL err(NF90_CLOSE (nid))
+  CALL err(nf90_close (nid))
 
   WRITE(lunout,*)TRIM(modname)//': iim,jjm,llm,iday_end',iim,jjm,llm,iday_end
@@ -169,6 +169,6 @@
   USE infotrac, ONLY: nqtot, tracers, type_trac
   USE control_mod
-  USE netcdf,   ONLY: NF90_OPEN,  NF90_NOWRITE, nf90_get_var, NF90_INQ_VARID,  &
-                      NF90_CLOSE, NF90_WRITE,   NF90_PUT_VAR, nf90_noerr
+  USE netcdf,   ONLY: nf90_open,  nf90_nowrite, nf90_get_var, NF90_INQ_VARID,  &
+                      nf90_close, NF90_WRITE,   nf90_put_var, nf90_noerr
   USE dynredem_mod, ONLY: dynredem_write_u, dynredem_write_v, dynredem_read_u, &
                           err, modname, fil, msg
@@ -202,5 +202,5 @@
 
   modname='dynredem1'; fil=fichnom
-  CALL err(NF90_OPEN(fil,NF90_WRITE,nid),"open",fil)
+  CALL err(nf90_open(fil,NF90_WRITE,nid),"open",fil)
 
 !--- Write/extend time coordinate
@@ -208,5 +208,5 @@
   var="temps"
   CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
-  CALL err(NF90_PUT_VAR(nid,vID,[time]),"put",var)
+  CALL err(nf90_put_var(nid,vID,[time]),"put",var)
   WRITE(lunout,*)TRIM(modname)//": Saving for ", nb, time
 
@@ -217,5 +217,5 @@
   tab_cntrl(31)=DBLE(itau_dyn + itaufin)
   CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
-  CALL err(NF90_PUT_VAR(nid,vID,tab_cntrl),"put",var)
+  CALL err(nf90_put_var(nid,vID,tab_cntrl),"put",var)
 
 !--- Save fields
@@ -229,5 +229,5 @@
   lread_inca=.FALSE.; fil="start_trac.nc"
   IF(ANY(type_trac == ['inca','inco'])) INQUIRE(FILE=fil,EXIST=lread_inca)
-  IF(lread_inca) CALL err(NF90_OPEN(fil,NF90_NOWRITE,nid_trac),"open")
+  IF(lread_inca) CALL err(nf90_open(fil,nf90_nowrite,nid_trac),"open")
 
 !--- Save tracers
@@ -245,7 +245,7 @@
     CALL dynredem_write_u(nid,var,q(:,:,:,iq),llm)
   END DO
-  CALL err(NF90_CLOSE(nid),"close")
+  CALL err(nf90_close(nid),"close")
   fil="start_trac.nc"
-  IF(lread_inca) CALL err(NF90_CLOSE(nid_trac),"close")
+  IF(lread_inca) CALL err(nf90_close(nid_trac),"close")
 
 END SUBROUTINE dynredem1
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3d/dynredem_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3d/dynredem_mod.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3d/dynredem_mod.F90	(revision 5100)
@@ -32,5 +32,5 @@
   start(:)=[1,1,1,1]; count(:)=[iip1,jjp1,ll,1]
   CALL err(NF90_INQ_VARID(ncid,id,nvarid),"inq",id)
-  CALL err(NF90_PUT_VAR(ncid,nvarid,var,start,count),"put",id)
+  CALL err(nf90_put_var(ncid,nvarid,var,start,count),"put",id)
   
 END SUBROUTINE dynredem_write_u
@@ -55,5 +55,5 @@
   start(:)=[1,1,1,1]; count(:)=[iip1,jjm,ll,1]
   CALL err(NF90_INQ_VARID(ncid,id,nvarid),"inq",id)
-  CALL err(NF90_PUT_VAR(ncid,nvarid,var,start,count),"put",id)
+  CALL err(nf90_put_var(ncid,nvarid,var,start,count),"put",id)
   
 END SUBROUTINE dynredem_write_v
@@ -96,7 +96,7 @@
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: units
 !===============================================================================
-  CALL err(NF90_DEF_VAR(ncid,var,nf90_format,did,nvarid),"inq",var)
-  IF(title/="")      CALL err(NF90_PUT_ATT(ncid,nvarid,"title",title),var)
-  IF(PRESENT(units)) CALL err(NF90_PUT_ATT(ncid,nvarid,"units",units),var)
+  CALL err(nf90_def_var(ncid,var,nf90_format,did,nvarid),"inq",var)
+  IF(title/="")      CALL err(nf90_put_att(ncid,nvarid,"title",title),var)
+  IF(PRESENT(units)) CALL err(nf90_put_att(ncid,nvarid,"units",units),var)
 
 END SUBROUTINE cre_var
@@ -119,6 +119,6 @@
   IF(     PRESENT(units)) CALL cre_var(ncid,var,title,did,units)
   IF(.NOT.PRESENT(units)) CALL cre_var(ncid,var,title,did)
-  CALL err(NF90_ENDDEF(ncid))
-  CALL err(NF90_PUT_VAR(ncid,nvarid,v),"put",var)
+  CALL err(nf90_enddef(ncid))
+  CALL err(nf90_put_var(ncid,nvarid,v),"put",var)
   CALL err(NF90_REDEF(ncid))
 
@@ -142,6 +142,6 @@
   IF(     PRESENT(units)) CALL cre_var(ncid,var,title,did,units)
   IF(.NOT.PRESENT(units)) CALL cre_var(ncid,var,title,did)
-  CALL err(NF90_ENDDEF(ncid))
-  CALL err(NF90_PUT_VAR(ncid,nvarid,v),"put",var)
+  CALL err(nf90_enddef(ncid))
+  CALL err(nf90_put_var(ncid,nvarid,v),"put",var)
   CALL err(NF90_REDEF(ncid))
 
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3d/guide_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3d/guide_mod.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3d/guide_mod.F90	(revision 5100)
@@ -175,5 +175,5 @@
     if (guide_plevs==1) then
        if (ncidpl==-99) then
-          rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
+          rcod=nf90_open('apbp.nc',nf90_nowrite, ncidpl)
           if (rcod/=nf90_noerr) THEN
              abort_message=' Nudging error -> no file apbp.nc'
@@ -183,5 +183,5 @@
     elseif (guide_plevs==2) then
        if (ncidpl==-99) then
-          rcod=nf90_open('P.nc',Nf90_NOWRITe,ncidpl)
+          rcod=nf90_open('P.nc',nf90_nowrite,ncidpl)
           if (rcod/=nf90_noerr) THEN
              abort_message=' Nudging error -> no file P.nc'
@@ -192,5 +192,5 @@
     elseif (guide_u) then
            if (ncidpl==-99) then
-               rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
+               rcod=nf90_open('u.nc',nf90_nowrite,ncidpl)
                if (rcod/=nf90_noerr) THEN
                   CALL abort_gcm(modname, &
@@ -1690,13 +1690,13 @@
 
 ! Creation des variables dimensions
-        ierr=NF90_DEF_VAR(nid,"LONU",NF90_FLOAT,id_lonu,vid_lonu)
-        ierr=NF90_DEF_VAR(nid,"LONV",NF90_FLOAT,id_lonv,vid_lonv)
-        ierr=NF90_DEF_VAR(nid,"LATU",NF90_FLOAT,id_latu,vid_latu)
-        ierr=NF90_DEF_VAR(nid,"LATV",NF90_FLOAT,id_latv,vid_latv)
-        ierr=NF90_DEF_VAR(nid,"LEVEL",NF90_FLOAT,id_lev,vid_lev)
-        ierr=NF90_DEF_VAR(nid,"cu",NF90_FLOAT,(/id_lonu,id_latu/),vid_cu)
-        ierr=NF90_DEF_VAR(nid,"cv",NF90_FLOAT,(/id_lonv,id_latv/),vid_cv)
-        ierr=NF90_DEF_VAR(nid,"au",NF90_FLOAT,(/id_lonu,id_latu/),vid_au)
-        ierr=NF90_DEF_VAR(nid,"av",NF90_FLOAT,(/id_lonv,id_latv/),vid_av)
+        ierr=nf90_def_var(nid,"LONU",nf90_float,id_lonu,vid_lonu)
+        ierr=nf90_def_var(nid,"LONV",nf90_float,id_lonv,vid_lonv)
+        ierr=nf90_def_var(nid,"LATU",nf90_float,id_latu,vid_latu)
+        ierr=nf90_def_var(nid,"LATV",nf90_float,id_latv,vid_latv)
+        ierr=nf90_def_var(nid,"LEVEL",nf90_float,id_lev,vid_lev)
+        ierr=nf90_def_var(nid,"cu",nf90_float,(/id_lonu,id_latu/),vid_cu)
+        ierr=nf90_def_var(nid,"cv",nf90_float,(/id_lonv,id_latv/),vid_cv)
+        ierr=nf90_def_var(nid,"au",nf90_float,(/id_lonu,id_latu/),vid_au)
+        ierr=nf90_def_var(nid,"av",nf90_float,(/id_lonv,id_latv/),vid_av)
         call nf95_def_var(nid, "alpha_T", nf90_float, (/id_lonv, id_latu/), &
              varid_alpha_t)
@@ -1724,33 +1724,33 @@
 ! Pressure (GCM)
         dim4=(/id_lonv,id_latu,id_lev,id_tim/)
-        ierr = NF90_DEF_VAR(nid,"SP",NF90_FLOAT,dim4,varid)
+        ierr = nf90_def_var(nid,"SP",nf90_float,dim4,varid)
 ! Surface pressure (guidage)
         IF (guide_P) THEN
             dim3=(/id_lonv,id_latu,id_tim/)
-            ierr = NF90_DEF_VAR(nid,"ps",NF90_FLOAT,dim3,varid)
+            ierr = nf90_def_var(nid,"ps",nf90_float,dim3,varid)
         ENDIF
 ! Zonal wind
         IF (guide_u) THEN
             dim4=(/id_lonu,id_latu,id_lev,id_tim/)
-            ierr = NF90_DEF_VAR(nid,"u",NF90_FLOAT,dim4,varid)
-            ierr = NF90_DEF_VAR(nid,"ua",NF90_FLOAT,dim4,varid)
-            ierr = NF90_DEF_VAR(nid,"ucov",NF90_FLOAT,dim4,varid)
+            ierr = nf90_def_var(nid,"u",nf90_float,dim4,varid)
+            ierr = nf90_def_var(nid,"ua",nf90_float,dim4,varid)
+            ierr = nf90_def_var(nid,"ucov",nf90_float,dim4,varid)
         ENDIF
 ! Merid. wind
         IF (guide_v) THEN
             dim4=(/id_lonv,id_latv,id_lev,id_tim/)
-            ierr = NF90_DEF_VAR(nid,"v",NF90_FLOAT,dim4,varid)
-            ierr = NF90_DEF_VAR(nid,"va",NF90_FLOAT,dim4,varid)
-            ierr = NF90_DEF_VAR(nid,"vcov",NF90_FLOAT,dim4,varid)
+            ierr = nf90_def_var(nid,"v",nf90_float,dim4,varid)
+            ierr = nf90_def_var(nid,"va",nf90_float,dim4,varid)
+            ierr = nf90_def_var(nid,"vcov",nf90_float,dim4,varid)
         ENDIF
 ! Pot. Temperature
         IF (guide_T) THEN
             dim4=(/id_lonv,id_latu,id_lev,id_tim/)
-            ierr = NF90_DEF_VAR(nid,"teta",NF90_FLOAT,dim4,varid)
+            ierr = nf90_def_var(nid,"teta",nf90_float,dim4,varid)
         ENDIF
 ! Specific Humidity
         IF (guide_Q) THEN
             dim4=(/id_lonv,id_latu,id_lev,id_tim/)
-            ierr = NF90_DEF_VAR(nid,"q",NF90_FLOAT,dim4,varid)
+            ierr = nf90_def_var(nid,"q",nf90_float,dim4,varid)
         ENDIF
         
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90	(revision 5100)
@@ -22,5 +22,5 @@
   USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
   USE readTracFiles_mod, ONLY: addPhase
-  use netcdf, only : NF90_NOWRITE,NF90_OPEN,nf90_noerr,NF90_INQ_VARID,NF90_CLOSE,nf90_get_var
+  use netcdf, only : nf90_nowrite,nf90_open,nf90_noerr,NF90_INQ_VARID,nf90_close,nf90_get_var
 
   !   Author:    Frederic Hourdin      original: 15/01/93
@@ -140,5 +140,5 @@
 
      relief=0.
-     ierr = NF90_OPEN ('relief_in.nc', NF90_NOWRITE,nid_relief)
+     ierr = nf90_open ('relief_in.nc', nf90_nowrite,nid_relief)
      if (ierr==nf90_noerr) THEN
          ierr=NF90_INQ_VARID(nid_relief,'RELIEF',varid)
@@ -150,5 +150,5 @@
          endif
      endif
-     ierr = NF90_CLOSE (nid_relief)
+     ierr = nf90_close (nid_relief)
 
      !------------------------------------------------------------------
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90	(revision 5100)
@@ -122,20 +122,20 @@
   !-------------------------
   !   Longitudes on "u" dynamical grid
-  status=NF90_DEF_VAR(ncid_out,'lonu',NF90_FLOAT,out_lonudim, out_lonuid)
+  status=nf90_def_var(ncid_out,'lonu',nf90_float,out_lonudim, out_lonuid)
   CALL handle_err(status)
   status=nf90_put_att(ncid_out,out_lonuid,'units','degrees_east')
   status=nf90_put_att(ncid_out,out_lonuid,'long_name','Longitude on u grid')
   !   Longitudes on "v" dynamical grid
-  status=NF90_DEF_VAR(ncid_out,'lonv',NF90_FLOAT,out_lonvdim, out_lonvid)
+  status=nf90_def_var(ncid_out,'lonv',nf90_float,out_lonvdim, out_lonvid)
   CALL handle_err(status)
   status=nf90_put_att(ncid_out,out_lonvid,'units','degrees_east')
   status=nf90_put_att(ncid_out,out_lonvid,'long_name','Longitude on v grid')
   !   Latitudes on "u" dynamical grid
-  status=NF90_DEF_VAR(ncid_out,'latu',NF90_FLOAT,out_latudim, out_latuid)
+  status=nf90_def_var(ncid_out,'latu',nf90_float,out_latudim, out_latuid)
   CALL handle_err(status)
   status=nf90_put_att(ncid_out,out_latuid,'units','degrees_north')
   status=nf90_put_att(ncid_out,out_latuid,'long_name','Latitude on u grid')
   !  Latitudes on "v" dynamical grid
-  status=NF90_DEF_VAR(ncid_out,'latv',NF90_FLOAT,out_latvdim, out_latvid)
+  status=nf90_def_var(ncid_out,'latv',nf90_float,out_latvdim, out_latvid)
   CALL handle_err(status)
   status=nf90_put_att(ncid_out,out_latvid,'units','degrees_north')
@@ -144,5 +144,5 @@
   out_dim(1)=out_lonudim
   out_dim(2)=out_latudim
-  status=NF90_DEF_VAR(ncid_out,'grille_u',NF90_FLOAT,out_dim, out_uid)
+  status=nf90_def_var(ncid_out,'grille_u',nf90_float,out_dim, out_uid)
   CALL handle_err(status)
   status=nf90_put_att(ncid_out,out_uid,'units','m/s')
@@ -151,5 +151,5 @@
   out_dim(1)=out_lonvdim
   out_dim(2)=out_latvdim
-  status=NF90_DEF_VAR(ncid_out,'grille_v',NF90_FLOAT,out_dim, out_vid)
+  status=nf90_def_var(ncid_out,'grille_v',nf90_float,out_dim, out_vid)
   CALL handle_err(status)
   status=nf90_put_att(ncid_out,out_vid,'units','m/s')
@@ -158,5 +158,5 @@
   out_dim(1)=out_lonvdim
   out_dim(2)=out_latudim
-  status=NF90_DEF_VAR(ncid_out,'grille_s',NF90_FLOAT,out_dim, out_tempid)
+  status=nf90_def_var(ncid_out,'grille_s',nf90_float,out_dim, out_tempid)
   CALL handle_err(status)
   status=nf90_put_att(ncid_out,out_tempid,'units','Kelvin')
@@ -165,5 +165,5 @@
   ! for INCA :
   ! vertical levels "presnivs"
-  status=NF90_DEF_VAR(ncid_out,'presnivs',NF90_FLOAT,out_levdim, presnivs_id)
+  status=nf90_def_var(ncid_out,'presnivs',nf90_float,out_levdim, presnivs_id)
   CALL handle_err(status)
   status=nf90_put_att(ncid_out,presnivs_id,'units','Pa')
@@ -172,10 +172,10 @@
   out_dim(1)=out_lonvdim
   out_dim(2)=out_latudim
-  status = nf90_def_var(ncid_out,'phis',NF90_FLOAT,out_dim,phis_id)
+  status = nf90_def_var(ncid_out,'phis',nf90_float,out_dim,phis_id)
   CALL handle_err(status)
   status=nf90_put_att(ncid_out,phis_id,'units','m')
   status=nf90_put_att(ncid_out,phis_id,'long_name','surface geopotential height')
   ! gridcell area
-  status = nf90_def_var(ncid_out,'aire',NF90_FLOAT,out_dim,area_id)
+  status = nf90_def_var(ncid_out,'aire',nf90_float,out_dim,area_id)
   CALL handle_err(status)
   status=nf90_put_att(ncid_out,area_id,'units','m2')
@@ -193,13 +193,13 @@
   !-------------------------
   ! 1D : lonu, lonv,latu,latv ; INCA : presnivs
-  status=NF90_PUT_VAR(ncid_out,out_lonuid,rlonudeg,[1],[iip1])
-  CALL handle_err(status)
-  status=NF90_PUT_VAR(ncid_out,out_lonvid,rlonvdeg,[1],[iip1])
-  CALL handle_err(status)
-  status=NF90_PUT_VAR(ncid_out,out_latuid,rlatudeg,[1],[jjp1])
-  CALL handle_err(status)
-  status=NF90_PUT_VAR(ncid_out,out_latvid,rlatvdeg,[1],[jjm])
-  CALL handle_err(status)
-  status=NF90_PUT_VAR(ncid_out,presnivs_id,rlev,[1],[llm])
+  status=nf90_put_var(ncid_out,out_lonuid,rlonudeg,[1],[iip1])
+  CALL handle_err(status)
+  status=nf90_put_var(ncid_out,out_lonvid,rlonvdeg,[1],[iip1])
+  CALL handle_err(status)
+  status=nf90_put_var(ncid_out,out_latuid,rlatudeg,[1],[jjp1])
+  CALL handle_err(status)
+  status=nf90_put_var(ncid_out,out_latvid,rlatvdeg,[1],[jjm])
+  CALL handle_err(status)
+  status=nf90_put_var(ncid_out,presnivs_id,rlev,[1],[llm])
   CALL handle_err(status)
 
@@ -209,11 +209,11 @@
   
   COUNT(2)=jjp1  ! for "u" and "s" grids
-  status=NF90_PUT_VAR(ncid_out,out_uid,uwnd,start, count)
+  status=nf90_put_var(ncid_out,out_uid,uwnd,start, count)
   CALL handle_err(status)
   COUNT(2)=jjm  ! for "v" grid
-  status=NF90_PUT_VAR(ncid_out,out_vid,vwnd,start, count)
+  status=nf90_put_var(ncid_out,out_vid,vwnd,start, count)
   CALL handle_err(status) 
   COUNT(2)=jjp1  ! as "s" grid, for all the following vars
-  status=NF90_PUT_VAR(ncid_out,out_tempid,temp,start, count)
+  status=nf90_put_var(ncid_out,out_tempid,temp,start, count)
   CALL handle_err(status)
   status = nf90_put_var(ncid_out, phis_id, phis_loc,start,count)
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/dynetat0_loc.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/dynetat0_loc.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/dynetat0_loc.F90	(revision 5100)
@@ -9,6 +9,6 @@
   USE infotrac,    ONLY: nqtot, tracers, niso, iqIsoPha, iH2O, isoName
   USE strings_mod, ONLY: maxlen, msg, strStack, real2str, int2str, strIdx
-  USE netcdf,      ONLY: NF90_OPEN,  NF90_NOWRITE, nf90_inquire_dimension, NF90_INQ_VARID, &
-                         NF90_CLOSE, nf90_get_var, NF90_INQUIRE_VARIABLE,  nf90_noerr
+  USE netcdf,      ONLY: nf90_open,  nf90_nowrite, nf90_inquire_dimension, NF90_INQ_VARID, &
+                         nf90_close, nf90_get_var, NF90_INQUIRE_VARIABLE,  nf90_noerr
   USE readTracFiles_mod, ONLY: new2oldH2O, newHNO3, oldHNO3, getKey
   USE control_mod, ONLY: planet_type
@@ -55,5 +55,5 @@
 !--- Initial state file opening
   var=fichnom
-  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
+  CALL err(nf90_open(var,nf90_nowrite,fID),"open",var)
   CALL get_var1("controle",tab_cntrl)
 
@@ -215,5 +215,5 @@
   END DO
   DEALLOCATE(q_glo)
-  CALL err(NF90_CLOSE(fID),"close",fichnom)
+  CALL err(nf90_close(fID),"close",fichnom)
   day_ini=day_ini+INT(time)
   time=time-INT(time)
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/dynredem_loc.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/dynredem_loc.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/dynredem_loc.F90	(revision 5100)
@@ -11,7 +11,7 @@
   USE strings_mod, ONLY: maxlen
   USE infotrac, ONLY: nqtot, tracers
-  USE netcdf, ONLY: NF90_CREATE, NF90_DEF_DIM, NF90_INQ_VARID, NF90_GLOBAL,    &
-                    NF90_CLOSE,  NF90_PUT_ATT, NF90_UNLIMITED, NF90_CLOBBER,   &
-                    NF90_64BIT_OFFSET
+  USE netcdf, ONLY: nf90_create, nf90_def_dim, NF90_INQ_VARID, nf90_global,    &
+                    nf90_close,  nf90_put_att, nf90_unlimited, nf90_clobber,   &
+                    nf90_64bit_offset
   USE dynredem_mod, ONLY: cre_var, put_var, err, modname, fil
   USE comvert_mod,  ONLY: ap, bp, presnivs, pa, preff, nivsig, nivsigs,&
@@ -112,18 +112,18 @@
 
 !--- File creation
-  CALL err(NF90_CREATE(fichnom,IOR(NF90_CLOBBER,NF90_64BIT_OFFSET),nid))
+  CALL err(nf90_create(fichnom,IOR(nf90_clobber,nf90_64bit_offset),nid))
 
 !--- Some global attributes
-  CALL err(NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier demarrage dynamique"))
+  CALL err(nf90_put_att(nid,nf90_global,"title","Fichier demarrage dynamique"))
 
 !--- Dimensions
-  CALL err(NF90_DEF_DIM(nid,"index", length, indexID))
-  CALL err(NF90_DEF_DIM(nid,"rlonu", iip1,   rlonuID))
-  CALL err(NF90_DEF_DIM(nid,"rlatu", jjp1,   rlatuID))
-  CALL err(NF90_DEF_DIM(nid,"rlonv", iip1,   rlonvID))
-  CALL err(NF90_DEF_DIM(nid,"rlatv", jjm,    rlatvID))
-  CALL err(NF90_DEF_DIM(nid,"sigs",  llm,        sID))
-  CALL err(NF90_DEF_DIM(nid,"sig",   llmp1,    sigID))
-  CALL err(NF90_DEF_DIM(nid,"temps", NF90_UNLIMITED, timID))
+  CALL err(nf90_def_dim(nid,"index", length, indexID))
+  CALL err(nf90_def_dim(nid,"rlonu", iip1,   rlonuID))
+  CALL err(nf90_def_dim(nid,"rlatu", jjp1,   rlatuID))
+  CALL err(nf90_def_dim(nid,"rlonv", iip1,   rlonvID))
+  CALL err(nf90_def_dim(nid,"rlatv", jjm,    rlatvID))
+  CALL err(nf90_def_dim(nid,"sigs",  llm,        sID))
+  CALL err(nf90_def_dim(nid,"sig",   llmp1,    sigID))
+  CALL err(nf90_def_dim(nid,"temps", nf90_unlimited, timID))
 
 !--- Define and save invariant fields
@@ -156,5 +156,5 @@
   CALL cre_var(nid,"masse","Masse d air"    ,[rlonvID,rlatuID,sID,timID])
   CALL cre_var(nid,"ps"   ,"Pression au sol",[rlonvID,rlatuID    ,timID])
-  CALL err(NF90_CLOSE (nid))
+  CALL err(nf90_close (nid))
 
   WRITE(lunout,*)TRIM(modname)//': iim,jjm,llm,iday_end',iim,jjm,llm,iday_end
@@ -178,6 +178,6 @@
   USE infotrac, ONLY: nqtot, tracers, type_trac
   USE control_mod
-  USE netcdf,   ONLY: NF90_OPEN,  NF90_NOWRITE, nf90_get_var, NF90_INQ_VARID,  &
-                      NF90_CLOSE, NF90_WRITE,   NF90_PUT_VAR, nf90_noerr
+  USE netcdf,   ONLY: nf90_open,  nf90_nowrite, nf90_get_var, NF90_INQ_VARID,  &
+                      nf90_close, NF90_WRITE,   nf90_put_var, nf90_noerr
   USE dynredem_mod, ONLY: dynredem_write_u, dynredem_write_v, dynredem_read_u, &
                           err, modname, fil, msg
@@ -213,5 +213,5 @@
   IF(mpi_rank==0) THEN !++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   modname='dynredem1_loc'; fil=fichnom
-  CALL err(NF90_OPEN(fil,NF90_WRITE,nid),"open",fil)
+  CALL err(nf90_open(fil,NF90_WRITE,nid),"open",fil)
 
 !--- Write/extend time coordinate
@@ -219,5 +219,5 @@
   var="temps"
   CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
-  CALL err(NF90_PUT_VAR(nid,vID,[time]),"put",var)
+  CALL err(nf90_put_var(nid,vID,[time]),"put",var)
   WRITE(lunout,*)TRIM(modname)//": Saving for ", nb, time
 
@@ -228,5 +228,5 @@
   tab_cntrl(31)=DBLE(itau_dyn + itaufin)
   CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
-  CALL err(NF90_PUT_VAR(nid,vID,tab_cntrl),"put",var)
+  CALL err(nf90_put_var(nid,vID,tab_cntrl),"put",var)
   END IF               !++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 !$OMP END MASTER
@@ -244,5 +244,5 @@
   fil="start_trac.nc"
   IF(ANY(type_trac == ['inca','inco'])) INQUIRE(FILE=fil,EXIST=lread_inca)
-  IF(lread_inca) CALL err(NF90_OPEN(fil,NF90_NOWRITE,nid_trac),"open")
+  IF(lread_inca) CALL err(nf90_open(fil,nf90_nowrite,nid_trac),"open")
 !$OMP END MASTER
 !$OMP BARRIER
@@ -266,7 +266,7 @@
 !$OMP MASTER
   IF(mpi_rank==0) THEN !++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  CALL err(NF90_CLOSE(nid),"close")
+  CALL err(nf90_close(nid),"close")
   fil="start_trac.nc"
-  IF(lread_inca) CALL err(NF90_CLOSE(nid_trac),"close")
+  IF(lread_inca) CALL err(nf90_close(nid_trac),"close")
   END IF               !++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 !$OMP END MASTER
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/dynredem_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/dynredem_mod.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/dynredem_mod.F90	(revision 5100)
@@ -52,5 +52,5 @@
     !$OMP MASTER
       start(3)=l
-      CALL err(NF90_PUT_VAR(ncid,nvarid,var_glo,start,count),"put",id)
+      CALL err(nf90_put_var(ncid,nvarid,var_glo,start,count),"put",id)
     !$OMP END MASTER
     END IF
@@ -102,5 +102,5 @@
     !$OMP MASTER
       start(3)=l
-      CALL err(NF90_PUT_VAR(ncid,nvarid,var_glo,start,count),"put",id)
+      CALL err(nf90_put_var(ncid,nvarid,var_glo,start,count),"put",id)
     !$OMP END MASTER
     END IF
@@ -182,7 +182,7 @@
   CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: units
 !===============================================================================
-  CALL err(NF90_DEF_VAR(ncid,var,nf90_format,did,nvarid),"inq",var)
-  IF(title/="")      CALL err(NF90_PUT_ATT(ncid,nvarid,"title",title),var)
-  IF(PRESENT(units)) CALL err(NF90_PUT_ATT(ncid,nvarid,"units",units),var)
+  CALL err(nf90_def_var(ncid,var,nf90_format,did,nvarid),"inq",var)
+  IF(title/="")      CALL err(nf90_put_att(ncid,nvarid,"title",title),var)
+  IF(PRESENT(units)) CALL err(nf90_put_att(ncid,nvarid,"units",units),var)
 
 END SUBROUTINE cre_var
@@ -208,9 +208,9 @@
   IF(     PRESENT(units)) CALL cre_var(ncid,var,title,did,units)
   IF(.NOT.PRESENT(units)) CALL cre_var(ncid,var,title,did)
-  CALL err(NF90_ENDDEF(ncid))
+  CALL err(nf90_enddef(ncid))
   nd=SIZE(did)
   DO k=1,nd; CALL err(nf90_inquire_dimension(ncid,did(k),len=nn(k))); END DO
-  IF(nd==1) CALL err(NF90_PUT_VAR(ncid,nvarid,RESHAPE(v,nn(1:1))),var)
-  IF(nd==2) CALL err(NF90_PUT_VAR(ncid,nvarid,RESHAPE(v,nn(1:2))),var)
+  IF(nd==1) CALL err(nf90_put_var(ncid,nvarid,RESHAPE(v,nn(1:1))),var)
+  IF(nd==2) CALL err(nf90_put_var(ncid,nvarid,RESHAPE(v,nn(1:2))),var)
   CALL err(NF90_REDEF(ncid))
 END SUBROUTINE put_var
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/guide_loc_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/guide_loc_mod.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/guide_loc_mod.F90	(revision 5100)
@@ -175,5 +175,5 @@
     if (guide_plevs==1) then
        if (ncidpl==-99) then
-          rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
+          rcod=nf90_open('apbp.nc',nf90_nowrite, ncidpl)
           if (rcod/=nf90_noerr) THEN
              abort_message=' Nudging error -> no file apbp.nc'
@@ -183,5 +183,5 @@
     elseif (guide_plevs==2) then
        if (ncidpl==-99) then
-          rcod=nf90_open('P.nc',Nf90_NOWRITe,ncidpl)
+          rcod=nf90_open('P.nc',nf90_nowrite,ncidpl)
           if (rcod/=nf90_noerr) THEN
              abort_message=' Nudging error -> no file P.nc'
@@ -192,5 +192,5 @@
     elseif (guide_u) then
        if (ncidpl==-99) then
-          rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
+          rcod=nf90_open('u.nc',nf90_nowrite,ncidpl)
           if (rcod/=nf90_noerr) THEN
              abort_message=' Nudging error -> no file u.nc'
Index: /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/iniacademic_loc.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/iniacademic_loc.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/dyn3dmem/iniacademic_loc.F90	(revision 5100)
@@ -23,5 +23,5 @@
   USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
   USE readTracFiles_mod, ONLY: addPhase
-  use netcdf, only : NF90_NOWRITE,NF90_OPEN,nf90_noerr,NF90_INQ_VARID,NF90_CLOSE, nf90_get_var
+  use netcdf, only : nf90_nowrite,nf90_open,nf90_noerr,NF90_INQ_VARID,nf90_close, nf90_get_var
 
   !   Author:    Frederic Hourdin      original: 15/01/93
@@ -152,5 +152,5 @@
 
      relief=0.
-     ierr = NF90_OPEN ('relief_in.nc', NF90_NOWRITE,nid_relief)
+     ierr = nf90_open ('relief_in.nc', nf90_nowrite,nid_relief)
      if (ierr==nf90_noerr) THEN
          ierr=NF90_INQ_VARID(nid_relief,'RELIEF',varid)
@@ -162,5 +162,5 @@
          endif
      endif
-     ierr = NF90_CLOSE (nid_relief)
+     ierr = nf90_close (nid_relief)
 
 
Index: /LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/ce0l.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/ce0l.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/ce0l.F90	(revision 5100)
@@ -21,5 +21,5 @@
   USE etat0phys,      ONLY: etat0phys_netcdf
   USE limit,          ONLY: limit_netcdf
-  USE netcdf,         ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, nf90_noerr,    &
+  USE netcdf,         ONLY: nf90_open, nf90_nowrite, nf90_close, nf90_noerr,    &
          nf90_inquire_dimension, nf90_inq_dimid, NF90_INQ_VARID, nf90_get_var
   USE infotrac,       ONLY: init_infotrac
@@ -154,6 +154,6 @@
 ! weights to ensure ocean fractions are the same for atmosphere and ocean.
 !*******************************************************************************
-  IF(NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)==nf90_noerr) THEN
-    iret=NF90_CLOSE(nid_o2a)
+  IF(nf90_open("o2a.nc", nf90_nowrite, nid_o2a)==nf90_noerr) THEN
+    iret=nf90_close(nid_o2a)
     WRITE(lunout,*)'BEWARE !! Ocean mask "o2a.nc" file found'
     WRITE(lunout,*)'Coupled run.'
@@ -188,5 +188,5 @@
     masque(iip1 ,:)=masque(1,:)
     DEALLOCATE(ocemask)
-  ELSE IF(NF90_OPEN("startphy0.nc", NF90_NOWRITE, nid_sta)==nf90_noerr) THEN
+  ELSE IF(nf90_open("startphy0.nc", nf90_nowrite, nid_sta)==nf90_noerr) THEN
     WRITE(lunout,*)'BEWARE !! File "startphy0.nc" found.'
     WRITE(lunout,*)'Getting the land mask from a previous run.'
@@ -196,5 +196,5 @@
       WRITE(lunout,*)'Mismatching dimensions for land mask'
       WRITE(lunout,*)'nphys  = ',nphys ,' klon = ',klon
-      iret=NF90_CLOSE(nid_sta)
+      iret=nf90_close(nid_sta)
       CALL abort_gcm(modname,'',1)
     END IF
@@ -202,5 +202,5 @@
     iret=NF90_INQ_VARID(nid_sta,'masque',nid_msk)
     iret=nf90_get_var(nid_sta,nid_msk,masktmp)
-    iret=NF90_CLOSE(nid_sta)
+    iret=nf90_close(nid_sta)
     CALL gr_fi_dyn(1,klon,iip1,jjp1,masktmp,masque)
     IF(prt_level>=1) THEN
Index: DZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/limit_netcdf.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/limit_netcdf.F90	(revision 5099)
+++ 	(revision )
@@ -1,834 +1,0 @@
-MODULE limit
-
-!*******************************************************************************
-! 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)
-!-------------------------------------------------------------------------------
-
-  USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo
-  USE assert_eq_m,        ONLY: assert_eq
-  USE cal_tools_m,        ONLY: year_len, mid_month
-  USE conf_dat_m,         ONLY: conf_dat2d, conf_dat3d
-  USE dimphy,             ONLY: klon, zmasq
-  USE geometry_mod,       ONLY: longitude_deg, latitude_deg
-  USE phys_state_var_mod, ONLY: pctsrf
-  USE control_mod,        ONLY: anneeref
-  USE init_ssrf_m,        ONLY: start_init_subsurf
-
-  INTEGER,           PARAMETER :: ns=256
-  CHARACTER(LEN=ns), PARAMETER :: &
-  fsst(5)=['amipbc_sst_1x1.nc   ','amip_sst_1x1.nc     ','cpl_atm_sst.nc      '&
-          ,'histmth_sst.nc      ','sstk.nc             '],                     &
-  fsic(5)=['amipbc_sic_1x1.nc   ','amip_sic_1x1.nc     ','cpl_atm_sic.nc      '&
-          ,'histmth_sic.nc      ','ci.nc               '],                     &
-  vsst(5)=['tosbcs    ','tos       ','SISUTESW  ','tsol_oce  ','sstk      '],  &
-  vsic(5)=['sicbcs    ','sic       ','SIICECOV  ','pourc_sic ','ci        '],  &
-  frugo='Rugos.nc  ', falbe='Albedo.nc ', frelf='Relief.nc ',                  &
-   vrug='RUGOS     ',  valb='ALBEDO    ',  vrel='RELIEF    ',                  &
-  DegK(11)=['degK          ','degree_K      ','degreeK       ','deg_K         '&
-           ,'degsK         ','degrees_K     ','degreesK      ','degs_K        '&
-           ,'degree_kelvin ','degrees_kelvin','K             '],               &
-  DegC(10)=['degC          ','degree_C      ','degreeC       ','deg_C         '&
-           ,'degsC         ','degrees_C     ','degreesC      ','degs_C        '&
-           ,'degree_Celsius','celsius       '], &
-  Perc(2) =['%             ','percent       '], &
-  Frac(2) =['1.0           ','1             ']
-
-CONTAINS
-
-!-------------------------------------------------------------------------------
-
-SUBROUTINE limit_netcdf(masque, phis, extrap)
-
-!-------------------------------------------------------------------------------
-! 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)
-!  *    04/2016: D. Cugnet   (12/14 recs SST/SIC files: cyclic/interannual runs)
-!  *    05/2017: D. Cugnet   (linear time interpolation for BCS files)
-!-------------------------------------------------------------------------------
-#ifndef CPP_1D
-  USE indice_sol_mod
-  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_GLOBAL,       &
-                  NF90_CLOBBER, NF90_ENDDEF,  NF90_UNLIMITED, NF90_FLOAT,      &
-                  NF90_64BIT_OFFSET
-  USE lmdz_cppkeys_wrapper, ONLY: nf90_format
-  USE inter_barxy_m,      ONLY: inter_barxy
-  USE netcdf95,           ONLY: nf95_def_var, nf95_put_att, nf95_put_var
-  USE comconst_mod, ONLY: pi
-  USE phys_cal_mod, ONLY: calend
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  include "iniprint.h"
-  include "dimensions.h"
-  include "paramet.h"
-  REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
-  REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: phis   ! ground geopotential
-  LOGICAL,                    INTENT(IN)    :: extrap ! SST extrapolation flag
-!-------------------------------------------------------------------------------
-! Local variables:
-  include "comgeom2.h"
-
-!--- INPUT NETCDF FILES AND VARIABLES NAMES ------------------------------------
-  CHARACTER(LEN=ns) :: icefile, sstfile, fnam, varname
-
-!--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
-  REAL               :: fi_ice(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 :: 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
-  INTEGER :: ndays                   !--- Depending on the output calendar
-  CHARACTER(LEN=ns) :: str
-
-!--- INITIALIZATIONS -----------------------------------------------------------
-  CALL inigeom
-
-!--- MASK, GROUND GEOPOT. & SUBSURFACES COMPUTATION (IN CASE ok_etat0==.FALSE.)
-   IF(ALL(masque==-99999.)) THEN
-    CALL start_init_orog0(rlonv,rlatu,phis,masque)
-    CALL gr_dyn_fi(1,iip1,jjp1,klon,masque,zmasq)          !--- To physical grid
-    ALLOCATE(pctsrf(klon,nbsrf))
-    CALL start_init_subsurf(.FALSE.)
-  !--- TO MATCH EXACTLY WHAT WOULD BE DONE IN etat0phys_netcdf
-    WHERE(   masque(:,:)<EPSFRA) masque(:,:)=0.
-    WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
-  END IF
-
-!--- Beware: anneeref (from gcm.def) is used to determine output time sampling
-  ndays=year_len(anneeref)
-
-!--- RUGOSITY TREATMENT --------------------------------------------------------
-  CALL msg(0,""); CALL msg(0," *** TRAITEMENT DE LA RUGOSITE ***")
-  CALL get_2Dfield(frugo,vrug,'RUG',ndays,phy_rug,mask=masque(1:iim,:))
-
-!--- OCEAN TREATMENT -----------------------------------------------------------
-  CALL msg(0,""); CALL msg(0," *** TRAITEMENT DE LA GLACE OCEANIQUE ***")
-
-! Input SIC file selection
-! Open file only to test if available
-  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 available : '
-     DO k=1,SIZE(fsic); WRITE(lunout,*) TRIM(fsic(k)); END DO
-     CALL abort_physic('limit_netcdf','No sea-ice file was found',1)
-  END IF
-  CALL ncerr(NF90_CLOSE(nid),icefile)
-  CALL msg(0,'Fichier choisi pour la glace de mer:'//TRIM(icefile))
-
-  CALL get_2Dfield(icefile,varname, 'SIC',ndays,phy_ice)
-
-  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
-     pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)       ! land soil
-     pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic)       ! land ice
-     SELECT CASE(ix_sic)
-        CASE(3)                                   ! SIC=pICE*(1-LIC-TER) (CPL)
-        pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
-        CASE(4)                                   ! SIC=pICE            (HIST)
-        pctsrf_t(:,is_sic,k)=fi_ice(:)
-        CASE DEFAULT                              ! SIC=pICE-LIC   (AMIP,ERAI)
-        pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
-     END SELECT
-     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 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
-  DEALLOCATE(phy_ice)
-
-!--- SST TREATMENT -------------------------------------------------------------
-  CALL msg(0,""); CALL msg(0," *** TRAITEMENT DE LA SST ***")
-
-! Input SST file selection
-! Open file only to test if available
-  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 available : '
-     DO k=1,SIZE(fsst); WRITE(lunout,*) TRIM(fsst(k)); END DO
-     CALL abort_physic('limit_netcdf','No sst file was found',1)
-  END IF
-  CALL ncerr(NF90_CLOSE(nid),sstfile)
-  CALL msg(0,'Fichier choisi pour la temperature de mer: '//TRIM(sstfile))
-
-  CALL get_2Dfield(sstfile,varname,'SST',ndays,phy_sst,flag=extrap)
-
-!--- ALBEDO TREATMENT ----------------------------------------------------------
-  CALL msg(0,""); CALL msg(0," *** TRAITEMENT DE L'ALBEDO ***")
-  CALL get_2Dfield(falbe,valb,'ALB',ndays,phy_alb)
-
-!--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
-  ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0
-
-!--- OUTPUT FILE WRITING -------------------------------------------------------
-  CALL msg(0,""); CALL msg(0,' *** Ecriture du fichier limit : debut ***')
-  fnam="limit.nc"
-
-  !--- File creation
-  CALL ncerr(NF90_CREATE(fnam,IOR(NF90_CLOBBER,NF90_64BIT_OFFSET),nid),fnam)
-  CALL ncerr(NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier conditions aux limites"),fnam)
-  str='File produced using ce0l executable.'
-  str=TRIM(str)//NEW_LINE(' ')//'Sea Ice Concentration built from'
-  SELECT CASE(ix_sic)
-    CASE(1); str=TRIM(str)//' Amip mid-month boundary condition (BCS).'
-    CASE(2); str=TRIM(str)//' Amip monthly mean observations.'
-    CASE(3); str=TRIM(str)//' IPSL coupled model outputs.'
-    CASE(4); str=TRIM(str)//' LMDZ model outputs.'
-    CASE(5); str=TRIM(str)//' ci.nc file.'
-  END SELECT
-  str=TRIM(str)//NEW_LINE(' ')//'Sea Surface Temperature built from'
-  SELECT CASE(ix_sst)
-    CASE(1); str=TRIM(str)//' Amip mid-month boundary condition (BCS).'
-    CASE(2); str=TRIM(str)//' Amip monthly mean observations.'
-    CASE(3); str=TRIM(str)//' IPSL coupled model outputs.'
-    CASE(4); str=TRIM(str)//' LMDZ model outputs.'
-    CASE(5); str=TRIM(str)//' sstk.nc file.'
-  END SELECT
-  CALL ncerr(NF90_PUT_ATT(nid,NF90_GLOBAL,"history",TRIM(str)),fnam)
-
-  !--- Dimensions creation
-  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
-  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)
-
-  !--- Attributes creation
-  CALL ncerr(NF90_PUT_ATT(nid,id_tim, "title","Jour dans l annee"),fnam)
-  CALL ncerr(NF90_PUT_ATT(nid,id_tim, "calendar",calend),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")
-  call nf95_put_att(nid, varid_longitude, "units", "degrees_east")
-
-  call nf95_put_att(nid, varid_latitude, "standard_name", "latitude")
-  call nf95_put_att(nid, varid_latitude, "units", "degrees_north")
-
-  CALL ncerr(NF90_ENDDEF(nid),fnam)
-
-  !--- Variables saving
-  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, longitude_deg)
-  call nf95_put_var(nid, varid_latitude, latitude_deg)
-
-  CALL ncerr(NF90_CLOSE(nid),fnam)
-
-  CALL msg(0,""); CALL msg(0,' *** Ecriture du fichier limit : fin ***')
-
-  DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
-
-
-!===============================================================================
-
-  CONTAINS
-
-!===============================================================================
-
-
-!-------------------------------------------------------------------------------
-
-SUBROUTINE get_2Dfield(fnam, varname, mode, 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 pchsp_95_m, only: pchsp_95
-  USE pchfe_95_m, only: pchfe_95
-  USE arth_m, only: arth
-  USE indice_sol_mod
-
-  IMPLICIT NONE
-  include "dimensions.h"
-  include "paramet.h"
-  include "comgeom2.h"
-!-----------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*),  INTENT(IN)     :: fnam     ! NetCDF file name
-  CHARACTER(LEN=*),  INTENT(IN)     :: varname  ! NetCDF variable name
-  CHARACTER(LEN=*),  INTENT(IN)     :: mode     ! RUG, SIC, SST or ALB
-  INTEGER,           INTENT(IN)     :: ndays    ! current year number of days
-  REAL,    POINTER,  DIMENSION(:, :) :: champo  ! output field = f(t)
-  LOGICAL, OPTIONAL, INTENT(IN)     :: flag     ! extrapol. (SST) old ice (SIC)
-  REAL,    OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
-!------------------------------------------------------------------------------
-! Local variables:
-!--- NetCDF
-  INTEGER           :: ncid, varid        ! NetCDF identifiers
-  CHARACTER(LEN=ns) :: dnam               ! dimension name
-!--- dimensions
-  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 :: 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=ns) :: fnam_m, fnam_p     ! previous/next files names
-  CHARACTER(LEN=ns) :: cal_in             ! calendar
-  CHARACTER(LEN=ns) :: units              ! attribute "units" in sic/sst file
-  INTEGER           :: ndays_in           ! number of days
-  REAL              :: value              ! mean/max value near equator
-!--- misc
-  INTEGER           :: i, j, k, l         ! loop counters
-  REAL, ALLOCATABLE :: work(:,:)          ! used for extrapolation
-  CHARACTER(LEN=ns) :: title, mess        ! for messages
-  LOGICAL           :: is_bcs             ! flag for BCS data
-  LOGICAL           :: extrp              ! flag for extrapolation
-  LOGICAL           :: ll
-  REAL              :: chmin, chmax, timeday, al
-  INTEGER ierr, idx
-  integer n_extrap ! number of extrapolated points
-  logical skip
-
-!------------------------------------------------------------------------------
-!---Variables depending on keyword 'mode' -------------------------------------
-  NULLIFY(champo)
-
-  SELECT CASE(mode)
-  CASE('RUG'); title='Rugosite'
-  CASE('SIC'); title='Sea-ice'
-  CASE('SST'); title='SST'
-  CASE('ALB'); title='Albedo'
-  END SELECT
-  extrp=.FALSE.; IF(PRESENT(flag).AND.mode=='SST') extrp=flag
-  is_bcs=(mode=='SIC'.AND.ix_sic==1).OR.(mode=='SST'.AND.ix_sst==1)
-  idx=INDEX(fnam,'.nc')-1
-
-!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
-  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)
-
-!--- Longitude
-  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
-  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)
-  CALL ncerr(nf90_inquire_dimension(ncid, dids(3), name=dnam, len=lmdep), fnam)
-  ALLOCATE(timeyear(lmdep+2))
-  CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
-  cal_in=' '
-  IF(NF90_GET_ATT(ncid, varid, 'calendar', cal_in)/=nf90_noerr) THEN
-    SELECT CASE(mode)
-      CASE('RUG', 'ALB'); cal_in='360_day'
-      CASE('SIC', 'SST'); cal_in='gregorian'
-    END SELECT
-    CALL msg(0,'WARNING: missing "calendar" attribute for "time" in '&
-      //TRIM(fnam)//'. Choosing default value.')
-  END IF
-  CALL strclean(cal_in)                     !--- REMOVE (WEIRD) NULL CHARACTERS
-  CALL msg(0,'var, calendar, dim: '//TRIM(dnam)//' '//TRIM(cal_in), lmdep)
-  
-!--- Determining input file number of days, depending on calendar
-  ndays_in=year_len(anneeref, cal_in)
-
-!--- Rebuilding input time vector (field from input file might be unreliable)
-  IF(lmdep==12) THEN
-    timeyear=mid_month(anneeref, cal_in)
-    CALL msg(0,'Monthly input file(s) for '//TRIM(title)//'.')
-  ELSE IF(lmdep==ndays_in) THEN
-    timeyear=[(REAL(k)-0.5,k=0,ndays_in+1)]
-    CALL msg(0,'Daily input file (no time interpolation).')
-  ELSE
-    WRITE(mess,'(a,i3,a,i3,a)')'Mismatching input file: found',lmdep,        &
-      ' records, 12/',ndays_in,' (monthly/daily needed).'
-    CALL abort_physic('mid_month',TRIM(mess),1)
-  END IF
-
-!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
-  ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep+2))
-  IF(extrp) ALLOCATE(work(imdep, jmdep))
-  CALL msg(5,'')
-  CALL msg(5,'READ AND INTERPOLATE HORIZONTALLY ', lmdep, ' FIELDS.')
-  CALL ncerr(NF90_INQ_VARID(ncid, varname, varid), fnam)
-  DO l=1, lmdep
-    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, .TRUE.)
-
-    !--- FOR SIC/SST FIELDS ONLY
-    IF(l==1.AND.is_in(mode,['SIC','SST'])) THEN
-
-      !--- DETERMINE THE UNIT: READ FROM FILE OR ASSUMED USING FIELD VALUES
-      ierr=NF90_GET_ATT(ncid, varid, 'units', units)
-      IF(ierr==nf90_noerr) THEN !--- ATTRIBUTE "units" FOUND IN THE FILE
-        CALL strclean(units)
-        IF(mode=='SIC'.AND.is_in(units,Perc)) units="%"
-        IF(mode=='SIC'.AND.is_in(units,Frac)) units="1"
-        IF(mode=='SST'.AND.is_in(units,DegC)) units="C"
-        IF(mode=='SST'.AND.is_in(units,DegK)) units="K"
-      ELSE                      !--- CHECK THE FIELD VALUES
-        IF(mode=='SIC') value=MAXVAL(champ(:,:))
-        IF(mode=='SST') value=   SUM(champ(:,jmdep/2),DIM=1)/REAL(imdep)
-        IF(mode=='SIC') THEN; units="1"; IF(value>= 10.) units="%"; END IF
-        IF(mode=='SST') THEN; units="C"; IF(value>=100.) units="K"; END IF
-      END IF
-      CALL msg(0,'INPUT FILE '//TRIM(title)//' UNIT IS: "'//TRIM(units)//'".')
-      IF(ierr/=nf90_noerr) CALL msg(0,'WARNING ! UNIT TO BE CHECKED ! '      &
-        //'No "units" attribute, so only based on the fields values.')
-
-      !--- CHECK VALUES ARE IN THE EXPECTED RANGE
-      SELECT CASE(units)
-        CASE('%'); ll=ANY(champ>100.0+EPSFRA); str='percentages > 100.'
-        CASE('1'); ll=ANY(champ>  1.0+EPSFRA); str='fractions > 1.'
-        CASE('C'); ll=ANY(champ<-100.).OR.ANY(champ> 60.); str='<-100 or >60 DegC'
-        CASE('K'); ll=ANY(champ< 180.).OR.ANY(champ>330.); str='<180 or >330 DegK'
-        CASE DEFAULT; CALL abort_physic(mode, 'Unrecognized '//TRIM(title)   &
-                                                  //' unit: '//TRIM(units),1)
-      END SELECT
-
-      !--- DROPPED FOR BCS DATA (FRACTIONS CAN BE HIGHER THAN 1)
-      IF(ll.AND.ix_sic/=1.AND.mode=='SIC') &
-        CALL abort_physic(mode,'unrealistic '//TRIM(mode)//' found: '//TRIM(str), 1)
-
-    END IF
-
-    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
-    IF(l==1) THEN
-      CALL msg(5,"--------------------------------------------------------")
-      CALL msg(5,"$$$ Barycentric interpolation for "//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)
-    IF(mode=='RUG') THEN
-      champint=EXP(champint)
-      WHERE(NINT(mask)/=1) champint=0.001
-    END IF
-    champtime(:, :, l+1)=champint
-  END DO
-  CALL ncerr(NF90_CLOSE(ncid), fnam)
-
-!--- FIRST RECORD: LAST ONE OF PREVIOUS YEAR (CURRENT YEAR IF UNAVAILABLE)
-  fnam_m=fnam(1:idx)//'_m.nc'
-  IF(NF90_OPEN(fnam_m,NF90_NOWRITE,ncid)==nf90_noerr) THEN
-    CALL msg(0,'Reading previous year file ("'//TRIM(fnam_m)//'") last record for '//TRIM(title))
-    CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_m)
-    CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids),fnam_m)
-    CALL ncerr(nf90_inquire_dimension(ncid, dids(3), len=l), fnam_m)
-    CALL ncerr(nf90_get_var(ncid,varid,champ,[1,1,l],[imdep,jmdep,1]),fnam_m)
-    CALL ncerr(NF90_CLOSE(ncid), fnam_m)
-    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
-    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
-    IF(mode=='RUG') champ=LOG(champ)
-    CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
-    IF(mode=='RUG') THEN
-      champint=EXP(champint)
-      WHERE(NINT(mask)/=1) champint=0.001
-    END IF
-    champtime(:, :, 1)=champint
-  ELSE
-    CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") last record for '//TRIM(title))
-    champtime(:, :, 1)=champtime(:, :, lmdep+1)
-  END IF
-
-!--- LAST RECORD: FIRST ONE OF NEXT YEAR (CURRENT YEAR IF UNAVAILABLE)
-  fnam_p=fnam(1:idx)//'_p.nc'
-  IF(NF90_OPEN(fnam_p,NF90_NOWRITE,ncid)==nf90_noerr) THEN
-    CALL msg(0,'Reading next year file ("'//TRIM(fnam_p)//'") first record for '//TRIM(title))
-    CALL ncerr(NF90_INQ_VARID(ncid, varname, varid),fnam_p)
-    CALL ncerr(nf90_get_var(ncid,varid,champ,[1,1,1],[imdep,jmdep,1]),fnam_p)
-    CALL ncerr(NF90_CLOSE(ncid), fnam_p)
-    CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
-    IF(extrp) CALL extrapol(champ,imdep,jmdep,999999.,.TRUE.,.TRUE.,2,work)
-    IF(mode=='RUG') champ=LOG(champ)
-    CALL inter_barxy(dlon,dlat(:jmdep-1),champ,rlonu(:iim),rlatv,champint)
-    IF(mode=='RUG') THEN
-      champint=EXP(champint)
-      WHERE(NINT(mask)/=1) champint=0.001
-    END IF
-    champtime(:, :, lmdep+2)=champint
-  ELSE
-    CALL msg(0,'Using current year file ("'//TRIM(fnam)//'") first record for '//TRIM(title))
-    champtime(:, :, lmdep+2)=champtime(:, :, 2)
-  END IF
-  DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
-  IF(extrp) DEALLOCATE(work)
-
-!--- TIME INTERPOLATION ------------------------------------------------------
-  IF(prt_level>0) THEN
-     IF(ndays/=ndays_in) THEN
-        WRITE(lunout,*)'DIFFERENT YEAR LENGTHS:'
-        WRITE(lunout,*)' In the  input file: ',ndays_in
-        WRITE(lunout,*)' In the output file: ',ndays
-     END IF
-     IF(lmdep==ndays_in) THEN
-        WRITE(lunout, *)'NO TIME INTERPOLATION.'
-        WRITE(lunout, *)' Daily input file.'
-     ELSE
-        IF(     is_bcs) WRITE(lunout, *)'LINEAR TIME INTERPOLATION.'
-        IF(.NOT.is_bcs) WRITE(lunout, *)'SPLINES TIME INTERPOLATION.'
-        WRITE(lunout, *)' Input time vector: ', timeyear
-        WRITE(lunout, *)' Output time vector: from 0.5 to ', ndays-0.5
-     END IF
-  END IF
-  ALLOCATE(champan(iip1, jjp1, ndays))
-
-  IF(lmdep==ndays_in) THEN  !--- DAILY DATA: NO     TIME INTERPOLATION
-    DO l=1,lmdep
-      champan(1:iim,:,l)=champtime(:,:,l+1)
-    END DO
-  ELSE IF(is_bcs) THEN      !--- BCS   DATA: LINEAR TIME INTERPOLATION
-    l=1
-    DO k=1, ndays
-      timeday = (REAL(k)-0.5)*REAL(ndays_in)/ndays
-      IF(timeyear(l+1)<timeday) l=l+1
-      al=(timeday-timeyear(l))/(timeyear(l+1)-timeyear(l))
-      champan(1:iim,:,k) = champtime(1:iim,:,l)+al*(champtime(1:iim,:,l+1)-champtime(1:iim,:,l))
-    END DO
-  ELSE                      !--- AVE   DATA: SPLINE TIME INTERPOLATION
-     skip = .false.
-     n_extrap = 0
-     ALLOCATE(yder(lmdep+2))
-     DO j=1, jjp1
-       DO i=1, iim
-         yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
-              vc_beg=0., vc_end=0.)
-         CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
-              arth(0.5, real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
-         if (ierr < 0) call abort_physic("get_2Dfield", "", 1)
-         n_extrap = n_extrap + ierr
-       END DO
-     END DO
-     IF(n_extrap /= 0) WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
-     DEALLOCATE(yder)
-  END IF
-  champan(iip1, :, :)=champan(1, :, :)
-  DEALLOCATE(champtime, timeyear)
-
-!--- Checking the result
-  DO j=1, jjp1
-    CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
-    IF (prt_level>5) WRITE(lunout, *)' ',TRIM(title),' at time 10 ', chmin, chmax, j
-  END DO
-
-!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
-  IF(mode=='SST') THEN
-    SELECT CASE(units)
-      CASE("K"); CALL msg(0,'SST field is already in kelvins.')
-      CASE("C"); CALL msg(0,'SST field converted from celcius degrees to kelvins.')
-      champan(:, :, :)=champan(:, :, :)+273.15
-    END SELECT
-    CALL msg(0,'Filtering SST: Sea Surface Temperature >= 271.38')
-    WHERE(champan<271.38) champan=271.38
-  END IF
-
-!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
-  IF(mode=='SIC') THEN
-    SELECT CASE(units)
-      CASE("1"); CALL msg(0,'SIC field already in fraction of 1')
-      CASE("%"); CALL msg(0,'SIC field converted from percentage to fraction of 1.')
-       champan(:, :, :)=champan(:, :, :)/100.
-    END SELECT
-    CALL msg(0,'Filtering SIC: 0.0 <= Sea-ice <=1.0')
-    WHERE(champan>1.0) champan=1.0
-    WHERE(champan<0.0) champan=0.0
- END IF
-
-!--- 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(:, k))
-  END DO
-  DEALLOCATE(champan)
-
-END SUBROUTINE get_2Dfield
-
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-
-SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque)
-
-!-------------------------------------------------------------------------------
-  USE grid_noro_m, ONLY: grid_noro0
-  IMPLICIT NONE
-!===============================================================================
-! Purpose:  Compute "phis" just like it would be in start_init_orog.
-!===============================================================================
-! Arguments:
-  REAL,             INTENT(IN)    :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
-  REAL,             INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
-!-------------------------------------------------------------------------------
-! Local variables:
-  CHARACTER(LEN=ns)  :: modname="start_init_orog0"
-  INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
-  REAL               :: lev(1), date, dt, deg2rad
-  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
-  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:), lat_rel(:,:)
-!-------------------------------------------------------------------------------
-  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
-  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
-  IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1)
-  IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1)
-  pi=2.0*ASIN(1.0); deg2rad=pi/180.0
-  IF(ANY(phis/=-99999.)) RETURN                  !--- phis ALREADY KNOWN
-
-!--- HIGH RESOLUTION OROGRAPHY
-  CALL flininfo(frelf, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
-
-  ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
-  CALL flinopen(frelf, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,   &
-                lev, ttm_tmp, itau, date, dt, fid)
-  ALLOCATE(relief_hi(iml_rel,jml_rel))
-  CALL flinget(fid, vrel, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
-  CALL flinclo(fid)
-
-!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
-  ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
-  lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
-  lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
-
-!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
-  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
-  CALL conf_dat2d(vrel, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
-  DEALLOCATE(lon_ini,lat_ini)
-
-!--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0
-  WRITE(lunout,*)
-  WRITE(lunout,*)'*** Compute surface geopotential ***'
-
-!--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS
-  CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque)
-  phis = phis * 9.81
-  phis(iml,:) = phis(1,:)
-  DEALLOCATE(relief_hi,lon_rad,lat_rad)
-
-END SUBROUTINE start_init_orog0
-
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-
-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
-
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-
-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
-!-------------------------------------------------------------------------------
-  IF(ncres/=nf90_noerr) THEN
-    WRITE(lunout,*)'Problem with file '//TRIM(fnam)//' in routine limit_netcdf.'
-    CALL abort_physic('limit_netcdf',NF90_STRERROR(ncres),1)
-  END IF
-
-END SUBROUTINE ncerr
-
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-
-SUBROUTINE strclean(s)
-
-!-------------------------------------------------------------------------------
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Purpose: Remove tail null characters from the input string.
-!-------------------------------------------------------------------------------
-! Parameters:
-  CHARACTER(LEN=*), INTENT(INOUT) :: s
-!-------------------------------------------------------------------------------
-! Local variable:
-  INTEGER :: k
-!-------------------------------------------------------------------------------
-  k=LEN_TRIM(s); DO WHILE(ICHAR(s(k:k))==0); s(k:k)=' '; k=LEN_TRIM(s); END DO
-
-END SUBROUTINE strclean
-
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-
-FUNCTION is_in(s1,s2) RESULT(res)
-
-!-------------------------------------------------------------------------------
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Purpose: Check wether s1 is present in the s2(:) list (case insensitive).
-!-------------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*), INTENT(IN) :: s1, s2(:)
-  LOGICAL                      :: res
-!-------------------------------------------------------------------------------
-  res=.FALSE.; DO k=1,SIZE(s2); res=res.OR.strLow(s1)==strLow(s2(k)); END DO
-
-END FUNCTION is_in
-
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-
-ELEMENTAL FUNCTION strLow(s) RESULT(res)
-
-!-------------------------------------------------------------------------------
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Purpose: Lower case conversion.
-!-------------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*), INTENT(IN) :: s
-  CHARACTER(LEN=ns)            :: res
-!-------------------------------------------------------------------------------
-! Local variable:
-  INTEGER :: k, ix
-!-------------------------------------------------------------------------------
-  res=s
-  DO k=1,LEN(s); ix=IACHAR(s(k:k))
-    IF(64<ix.AND.ix<91) res(k:k)=ACHAR(ix+32)
-  END DO
-
-END FUNCTION strLow
-
-!-------------------------------------------------------------------------------
-
-#endif
-! of #ifndef CPP_1D
-END SUBROUTINE limit_netcdf
-
-END MODULE limit
-
-!*******************************************************************************
-
Index: /LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/limit_netcdf.f90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/limit_netcdf.f90	(revision 5100)
+++ /LMDZ6/branches/Amaury_dev/libf/dynphy_lonlat/phylmd/limit_netcdf.f90	(revision 5100)
@@ -0,0 +1,837 @@
+MODULE limit
+
+  !*******************************************************************************
+  ! 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)
+  !-------------------------------------------------------------------------------
+
+  USE ioipsl, ONLY : flininfo, flinopen, flinget, flinclo
+  USE assert_eq_m, ONLY : assert_eq
+  USE cal_tools_m, ONLY : year_len, mid_month
+  USE conf_dat_m, ONLY : conf_dat2d, conf_dat3d
+  USE dimphy, ONLY : klon, zmasq
+  USE geometry_mod, ONLY : longitude_deg, latitude_deg
+  USE phys_state_var_mod, ONLY : pctsrf
+  USE control_mod, ONLY : anneeref
+  USE init_ssrf_m, ONLY : start_init_subsurf
+
+  INTEGER, PARAMETER :: ns = 256
+  CHARACTER(LEN = ns), PARAMETER :: &
+          fsst(5) = ['amipbc_sst_1x1.nc   ', 'amip_sst_1x1.nc     ', 'cpl_atm_sst.nc      '&
+                  , 'histmth_sst.nc      ', 'sstk.nc             '], &
+          fsic(5) = ['amipbc_sic_1x1.nc   ', 'amip_sic_1x1.nc     ', 'cpl_atm_sic.nc      '&
+                  , 'histmth_sic.nc      ', 'ci.nc               '], &
+          vsst(5) = ['tosbcs    ', 'tos       ', 'SISUTESW  ', 'tsol_oce  ', 'sstk      '], &
+          vsic(5) = ['sicbcs    ', 'sic       ', 'SIICECOV  ', 'pourc_sic ', 'ci        '], &
+          frugo = 'Rugos.nc  ', falbe = 'Albedo.nc ', frelf = 'Relief.nc ', &
+          vrug = 'RUGOS     ', valb = 'ALBEDO    ', vrel = 'RELIEF    ', &
+          DegK(11) = ['degK          ', 'degree_K      ', 'degreeK       ', 'deg_K         '&
+                  , 'degsK         ', 'degrees_K     ', 'degreesK      ', 'degs_K        '&
+                  , 'degree_kelvin ', 'degrees_kelvin', 'K             '], &
+          DegC(10) = ['degC          ', 'degree_C      ', 'degreeC       ', 'deg_C         '&
+                  , 'degsC         ', 'degrees_C     ', 'degreesC      ', 'degs_C        '&
+                  , 'degree_Celsius', 'celsius       '], &
+          Perc(2) = ['%             ', 'percent       '], &
+          Frac(2) = ['1.0           ', '1             ']
+
+CONTAINS
+
+  !-------------------------------------------------------------------------------
+
+  SUBROUTINE limit_netcdf(masque, phis, extrap)
+
+    !-------------------------------------------------------------------------------
+    ! 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)
+    !  *    04/2016: D. Cugnet   (12/14 recs SST/SIC files: cyclic/interannual runs)
+    !  *    05/2017: D. Cugnet   (linear time interpolation for BCS files)
+    !-------------------------------------------------------------------------------
+    USE indice_sol_mod
+    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_global, &
+            nf90_clobber, nf90_enddef, nf90_unlimited, nf90_float, &
+            nf90_64bit_offset
+    USE lmdz_cppkeys_wrapper, ONLY : nf90_format
+    USE inter_barxy_m, ONLY : inter_barxy
+    USE netcdf95, ONLY : nf95_def_var, nf95_put_att, nf95_put_var
+    USE comconst_mod, ONLY : pi
+    USE phys_cal_mod, ONLY : calend
+    IMPLICIT NONE
+    !-------------------------------------------------------------------------------
+    ! Arguments:
+    include "iniprint.h"
+    include "dimensions.h"
+    include "paramet.h"
+    REAL, DIMENSION(iip1, jjp1), INTENT(INOUT) :: masque ! land mask
+    REAL, DIMENSION(iip1, jjp1), INTENT(INOUT) :: phis   ! ground geopotential
+    LOGICAL, INTENT(IN) :: extrap ! SST extrapolation flag
+    !-------------------------------------------------------------------------------
+    ! Local variables:
+    include "comgeom2.h"
+
+    !--- INPUT NETCDF FILES AND VARIABLES NAMES ------------------------------------
+    CHARACTER(LEN = ns) :: icefile, sstfile, fnam, varname
+
+    !--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
+    REAL :: fi_ice(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 :: 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
+    INTEGER :: ndays                   !--- Depending on the output calendar
+    CHARACTER(LEN = ns) :: str
+
+    !--- INITIALIZATIONS -----------------------------------------------------------
+    CALL inigeom
+
+    !--- MASK, GROUND GEOPOT. & SUBSURFACES COMPUTATION (IN CASE ok_etat0==.FALSE.)
+    IF(ALL(masque==-99999.)) THEN
+      CALL start_init_orog0(rlonv, rlatu, phis, masque)
+      CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq)          !--- To physical grid
+      ALLOCATE(pctsrf(klon, nbsrf))
+      CALL start_init_subsurf(.FALSE.)
+      !--- TO MATCH EXACTLY WHAT WOULD BE DONE IN etat0phys_netcdf
+      WHERE(masque(:, :)<EPSFRA) masque(:, :) = 0.
+      WHERE(1. - masque(:, :)<EPSFRA) masque(:, :) = 1.
+    END IF
+
+    !--- Beware: anneeref (from gcm.def) is used to determine output time sampling
+    ndays = year_len(anneeref)
+
+    !--- RUGOSITY TREATMENT --------------------------------------------------------
+    CALL msg(0, ""); CALL msg(0, " *** TRAITEMENT DE LA RUGOSITE ***")
+    CALL get_2Dfield(frugo, vrug, 'RUG', ndays, phy_rug, mask = masque(1:iim, :))
+
+    !--- OCEAN TREATMENT -----------------------------------------------------------
+    CALL msg(0, ""); CALL msg(0, " *** TRAITEMENT DE LA GLACE OCEANIQUE ***")
+
+    ! Input SIC file selection
+    ! Open file only to test if available
+    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 available : '
+      DO k = 1, SIZE(fsic); WRITE(lunout, *) TRIM(fsic(k));
+      END DO
+      CALL abort_physic('limit_netcdf', 'No sea-ice file was found', 1)
+    END IF
+    CALL ncerr(nf90_close(nid), icefile)
+    CALL msg(0, 'Fichier choisi pour la glace de mer:' // TRIM(icefile))
+
+    CALL get_2Dfield(icefile, varname, 'SIC', ndays, phy_ice)
+
+    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
+      pctsrf_t(:, is_ter, k) = pctsrf(:, is_ter)       ! land soil
+      pctsrf_t(:, is_lic, k) = pctsrf(:, is_lic)       ! land ice
+      SELECT CASE(ix_sic)
+      CASE(3)                                   ! SIC=pICE*(1-LIC-TER) (CPL)
+        pctsrf_t(:, is_sic, k) = fi_ice(:) * (1. - pctsrf(:, is_lic) - pctsrf(:, is_ter))
+      CASE(4)                                   ! SIC=pICE            (HIST)
+        pctsrf_t(:, is_sic, k) = fi_ice(:)
+      CASE DEFAULT                              ! SIC=pICE-LIC   (AMIP,ERAI)
+        pctsrf_t(:, is_sic, k) = fi_ice - pctsrf_t(:, is_lic, k)
+      END SELECT
+      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 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
+    DEALLOCATE(phy_ice)
+
+    !--- SST TREATMENT -------------------------------------------------------------
+    CALL msg(0, ""); CALL msg(0, " *** TRAITEMENT DE LA SST ***")
+
+    ! Input SST file selection
+    ! Open file only to test if available
+    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 available : '
+      DO k = 1, SIZE(fsst); WRITE(lunout, *) TRIM(fsst(k));
+      END DO
+      CALL abort_physic('limit_netcdf', 'No sst file was found', 1)
+    END IF
+    CALL ncerr(nf90_close(nid), sstfile)
+    CALL msg(0, 'Fichier choisi pour la temperature de mer: ' // TRIM(sstfile))
+
+    CALL get_2Dfield(sstfile, varname, 'SST', ndays, phy_sst, flag = extrap)
+
+    !--- ALBEDO TREATMENT ----------------------------------------------------------
+    CALL msg(0, ""); CALL msg(0, " *** TRAITEMENT DE L'ALBEDO ***")
+    CALL get_2Dfield(falbe, valb, 'ALB', ndays, phy_alb)
+
+    !--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
+    ALLOCATE(phy_bil(klon, ndays)); phy_bil = 0.0
+
+    !--- OUTPUT FILE WRITING -------------------------------------------------------
+    CALL msg(0, ""); CALL msg(0, ' *** Ecriture du fichier limit : debut ***')
+    fnam = "limit.nc"
+
+    !--- File creation
+    CALL ncerr(nf90_create(fnam, IOR(nf90_clobber, nf90_64bit_offset), nid), fnam)
+    CALL ncerr(nf90_put_att(nid, nf90_global, "title", "Fichier conditions aux limites"), fnam)
+    str = 'File produced using ce0l executable.'
+    str = TRIM(str) // NEW_LINE(' ') // 'Sea Ice Concentration built from'
+    SELECT CASE(ix_sic)
+    CASE(1); str = TRIM(str) // ' Amip mid-month boundary condition (BCS).'
+    CASE(2); str = TRIM(str) // ' Amip monthly mean observations.'
+    CASE(3); str = TRIM(str) // ' IPSL coupled model outputs.'
+    CASE(4); str = TRIM(str) // ' LMDZ model outputs.'
+    CASE(5); str = TRIM(str) // ' ci.nc file.'
+    END SELECT
+    str = TRIM(str) // NEW_LINE(' ') // 'Sea Surface Temperature built from'
+    SELECT CASE(ix_sst)
+    CASE(1); str = TRIM(str) // ' Amip mid-month boundary condition (BCS).'
+    CASE(2); str = TRIM(str) // ' Amip monthly mean observations.'
+    CASE(3); str = TRIM(str) // ' IPSL coupled model outputs.'
+    CASE(4); str = TRIM(str) // ' LMDZ model outputs.'
+    CASE(5); str = TRIM(str) // ' sstk.nc file.'
+    END SELECT
+    CALL ncerr(nf90_put_att(nid, nf90_global, "history", TRIM(str)), fnam)
+
+    !--- Dimensions creation
+    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
+    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)
+
+    !--- Attributes creation
+    CALL ncerr(nf90_put_att(nid, id_tim, "title", "Jour dans l annee"), fnam)
+    CALL ncerr(nf90_put_att(nid, id_tim, "calendar", calend), 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")
+    call nf95_put_att(nid, varid_longitude, "units", "degrees_east")
+
+    call nf95_put_att(nid, varid_latitude, "standard_name", "latitude")
+    call nf95_put_att(nid, varid_latitude, "units", "degrees_north")
+
+    CALL ncerr(nf90_enddef(nid), fnam)
+
+    !--- Variables saving
+    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, longitude_deg)
+    call nf95_put_var(nid, varid_latitude, latitude_deg)
+
+    CALL ncerr(nf90_close(nid), fnam)
+
+    CALL msg(0, ""); CALL msg(0, ' *** Ecriture du fichier limit : fin ***')
+
+    DEALLOCATE(pctsrf_t, phy_sst, phy_bil, phy_alb, phy_rug)
+
+
+    !===============================================================================
+
+  CONTAINS
+
+    !===============================================================================
+
+
+    !-------------------------------------------------------------------------------
+
+    SUBROUTINE get_2Dfield(fnam, varname, mode, 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 pchsp_95_m, only : pchsp_95
+      USE pchfe_95_m, only : pchfe_95
+      USE arth_m, only : arth
+      USE indice_sol_mod
+
+      IMPLICIT NONE
+      include "dimensions.h"
+      include "paramet.h"
+      include "comgeom2.h"
+      !-----------------------------------------------------------------------------
+      ! Arguments:
+      CHARACTER(LEN = *), INTENT(IN) :: fnam     ! NetCDF file name
+      CHARACTER(LEN = *), INTENT(IN) :: varname  ! NetCDF variable name
+      CHARACTER(LEN = *), INTENT(IN) :: mode     ! RUG, SIC, SST or ALB
+      INTEGER, INTENT(IN) :: ndays    ! current year number of days
+      REAL, POINTER, DIMENSION(:, :) :: champo  ! output field = f(t)
+      LOGICAL, OPTIONAL, INTENT(IN) :: flag     ! extrapol. (SST) old ice (SIC)
+      REAL, OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
+      !------------------------------------------------------------------------------
+      ! Local variables:
+      !--- NetCDF
+      INTEGER :: ncid, varid        ! NetCDF identifiers
+      CHARACTER(LEN = ns) :: dnam               ! dimension name
+      !--- dimensions
+      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 :: 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 = ns) :: fnam_m, fnam_p     ! previous/next files names
+      CHARACTER(LEN = ns) :: cal_in             ! calendar
+      CHARACTER(LEN = ns) :: units              ! attribute "units" in sic/sst file
+      INTEGER :: ndays_in           ! number of days
+      REAL :: value              ! mean/max value near equator
+      !--- misc
+      INTEGER :: i, j, k, l         ! loop counters
+      REAL, ALLOCATABLE :: work(:, :)          ! used for extrapolation
+      CHARACTER(LEN = ns) :: title, mess        ! for messages
+      LOGICAL :: is_bcs             ! flag for BCS data
+      LOGICAL :: extrp              ! flag for extrapolation
+      LOGICAL :: ll
+      REAL :: chmin, chmax, timeday, al
+      INTEGER ierr, idx
+      integer n_extrap ! number of extrapolated points
+      logical skip
+
+      !------------------------------------------------------------------------------
+      !---Variables depending on keyword 'mode' -------------------------------------
+      NULLIFY(champo)
+
+      SELECT CASE(mode)
+      CASE('RUG'); title = 'Rugosite'
+      CASE('SIC'); title = 'Sea-ice'
+      CASE('SST'); title = 'SST'
+      CASE('ALB'); title = 'Albedo'
+      END SELECT
+      extrp = .FALSE.; IF(PRESENT(flag).AND.mode=='SST') extrp = flag
+      is_bcs = (mode=='SIC'.AND.ix_sic==1).OR.(mode=='SST'.AND.ix_sst==1)
+      idx = INDEX(fnam, '.nc') - 1
+
+      !--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
+      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)
+
+      !--- Longitude
+      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
+      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)
+      CALL ncerr(nf90_inquire_dimension(ncid, dids(3), name = dnam, len = lmdep), fnam)
+      ALLOCATE(timeyear(lmdep + 2))
+      CALL ncerr(NF90_INQ_VARID(ncid, dnam, varid), fnam)
+      cal_in = ' '
+      IF(nf90_get_att(ncid, varid, 'calendar', cal_in)/=nf90_noerr) THEN
+        SELECT CASE(mode)
+        CASE('RUG', 'ALB'); cal_in = '360_day'
+        CASE('SIC', 'SST'); cal_in = 'gregorian'
+        END SELECT
+        CALL msg(0, 'WARNING: missing "calendar" attribute for "time" in '&
+                // TRIM(fnam) // '. Choosing default value.')
+      END IF
+      CALL strclean(cal_in)                     !--- REMOVE (WEIRD) NULL CHARACTERS
+      CALL msg(0, 'var, calendar, dim: ' // TRIM(dnam) // ' ' // TRIM(cal_in), lmdep)
+
+      !--- Determining input file number of days, depending on calendar
+      ndays_in = year_len(anneeref, cal_in)
+
+      !--- Rebuilding input time vector (field from input file might be unreliable)
+      IF(lmdep==12) THEN
+        timeyear = mid_month(anneeref, cal_in)
+        CALL msg(0, 'Monthly input file(s) for ' // TRIM(title) // '.')
+      ELSE IF(lmdep==ndays_in) THEN
+        timeyear = [(REAL(k) - 0.5, k = 0, ndays_in + 1)]
+        CALL msg(0, 'Daily input file (no time interpolation).')
+      ELSE
+        WRITE(mess, '(a,i3,a,i3,a)')'Mismatching input file: found', lmdep, &
+                ' records, 12/', ndays_in, ' (monthly/daily needed).'
+        CALL abort_physic('mid_month', TRIM(mess), 1)
+      END IF
+
+      !--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
+      ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep + 2))
+      IF(extrp) ALLOCATE(work(imdep, jmdep))
+      CALL msg(5, '')
+      CALL msg(5, 'READ AND INTERPOLATE HORIZONTALLY ', lmdep, ' FIELDS.')
+      CALL ncerr(NF90_INQ_VARID(ncid, varname, varid), fnam)
+      DO l = 1, lmdep
+        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, .TRUE.)
+
+        !--- FOR SIC/SST FIELDS ONLY
+        IF(l==1.AND.is_in(mode, ['SIC', 'SST'])) THEN
+
+          !--- DETERMINE THE UNIT: READ FROM FILE OR ASSUMED USING FIELD VALUES
+          ierr = nf90_get_att(ncid, varid, 'units', units)
+          IF(ierr==nf90_noerr) THEN !--- ATTRIBUTE "units" FOUND IN THE FILE
+            CALL strclean(units)
+            IF(mode=='SIC'.AND.is_in(units, Perc)) units = "%"
+            IF(mode=='SIC'.AND.is_in(units, Frac)) units = "1"
+            IF(mode=='SST'.AND.is_in(units, DegC)) units = "C"
+            IF(mode=='SST'.AND.is_in(units, DegK)) units = "K"
+          ELSE                      !--- CHECK THE FIELD VALUES
+            IF(mode=='SIC') value = MAXVAL(champ(:, :))
+            IF(mode=='SST') value = SUM(champ(:, jmdep / 2), DIM = 1) / REAL(imdep)
+            IF(mode=='SIC') THEN; units = "1"; IF(value>= 10.) units = "%";
+            END IF
+            IF(mode=='SST') THEN; units = "C"; IF(value>=100.) units = "K";
+            END IF
+          END IF
+          CALL msg(0, 'INPUT FILE ' // TRIM(title) // ' UNIT IS: "' // TRIM(units) // '".')
+          IF(ierr/=nf90_noerr) CALL msg(0, 'WARNING ! UNIT TO BE CHECKED ! '      &
+                  // 'No "units" attribute, so only based on the fields values.')
+
+          !--- CHECK VALUES ARE IN THE EXPECTED RANGE
+          SELECT CASE(units)
+          CASE('%'); ll = ANY(champ>100.0 + EPSFRA); str = 'percentages > 100.'
+          CASE('1'); ll = ANY(champ>  1.0 + EPSFRA); str = 'fractions > 1.'
+          CASE('C'); ll = ANY(champ<-100.).OR.ANY(champ> 60.); str = '<-100 or >60 DegC'
+          CASE('K'); ll = ANY(champ< 180.).OR.ANY(champ>330.); str = '<180 or >330 DegK'
+          CASE DEFAULT; CALL abort_physic(mode, 'Unrecognized ' // TRIM(title)   &
+                  // ' unit: ' // TRIM(units), 1)
+          END SELECT
+
+          !--- DROPPED FOR BCS DATA (FRACTIONS CAN BE HIGHER THAN 1)
+          IF(ll.AND.ix_sic/=1.AND.mode=='SIC') &
+                  CALL abort_physic(mode, 'unrealistic ' // TRIM(mode) // ' found: ' // TRIM(str), 1)
+
+        END IF
+
+        IF(extrp) CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, work)
+        IF(l==1) THEN
+          CALL msg(5, "--------------------------------------------------------")
+          CALL msg(5, "$$$ Barycentric interpolation for " // 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)
+        IF(mode=='RUG') THEN
+          champint = EXP(champint)
+          WHERE(NINT(mask)/=1) champint = 0.001
+        END IF
+        champtime(:, :, l + 1) = champint
+      END DO
+      CALL ncerr(nf90_close(ncid), fnam)
+
+      !--- FIRST RECORD: LAST ONE OF PREVIOUS YEAR (CURRENT YEAR IF UNAVAILABLE)
+      fnam_m = fnam(1:idx) // '_m.nc'
+      IF(nf90_open(fnam_m, nf90_nowrite, ncid)==nf90_noerr) THEN
+        CALL msg(0, 'Reading previous year file ("' // TRIM(fnam_m) // '") last record for ' // TRIM(title))
+        CALL ncerr(NF90_INQ_VARID(ncid, varname, varid), fnam_m)
+        CALL ncerr(NF90_INQUIRE_VARIABLE(ncid, varid, dimids = dids), fnam_m)
+        CALL ncerr(nf90_inquire_dimension(ncid, dids(3), len = l), fnam_m)
+        CALL ncerr(nf90_get_var(ncid, varid, champ, [1, 1, l], [imdep, jmdep, 1]), fnam_m)
+        CALL ncerr(nf90_close(ncid), fnam_m)
+        CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
+        IF(extrp) CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, work)
+        IF(mode=='RUG') champ = LOG(champ)
+        CALL inter_barxy(dlon, dlat(:jmdep - 1), champ, rlonu(:iim), rlatv, champint)
+        IF(mode=='RUG') THEN
+          champint = EXP(champint)
+          WHERE(NINT(mask)/=1) champint = 0.001
+        END IF
+        champtime(:, :, 1) = champint
+      ELSE
+        CALL msg(0, 'Using current year file ("' // TRIM(fnam) // '") last record for ' // TRIM(title))
+        champtime(:, :, 1) = champtime(:, :, lmdep + 1)
+      END IF
+
+      !--- LAST RECORD: FIRST ONE OF NEXT YEAR (CURRENT YEAR IF UNAVAILABLE)
+      fnam_p = fnam(1:idx) // '_p.nc'
+      IF(nf90_open(fnam_p, nf90_nowrite, ncid)==nf90_noerr) THEN
+        CALL msg(0, 'Reading next year file ("' // TRIM(fnam_p) // '") first record for ' // TRIM(title))
+        CALL ncerr(NF90_INQ_VARID(ncid, varname, varid), fnam_p)
+        CALL ncerr(nf90_get_var(ncid, varid, champ, [1, 1, 1], [imdep, jmdep, 1]), fnam_p)
+        CALL ncerr(nf90_close(ncid), fnam_p)
+        CALL conf_dat2d(title, dlon_ini, dlat_ini, dlon, dlat, champ, .TRUE.)
+        IF(extrp) CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, work)
+        IF(mode=='RUG') champ = LOG(champ)
+        CALL inter_barxy(dlon, dlat(:jmdep - 1), champ, rlonu(:iim), rlatv, champint)
+        IF(mode=='RUG') THEN
+          champint = EXP(champint)
+          WHERE(NINT(mask)/=1) champint = 0.001
+        END IF
+        champtime(:, :, lmdep + 2) = champint
+      ELSE
+        CALL msg(0, 'Using current year file ("' // TRIM(fnam) // '") first record for ' // TRIM(title))
+        champtime(:, :, lmdep + 2) = champtime(:, :, 2)
+      END IF
+      DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
+      IF(extrp) DEALLOCATE(work)
+
+      !--- TIME INTERPOLATION ------------------------------------------------------
+      IF(prt_level>0) THEN
+        IF(ndays/=ndays_in) THEN
+          WRITE(lunout, *)'DIFFERENT YEAR LENGTHS:'
+          WRITE(lunout, *)' In the  input file: ', ndays_in
+          WRITE(lunout, *)' In the output file: ', ndays
+        END IF
+        IF(lmdep==ndays_in) THEN
+          WRITE(lunout, *)'NO TIME INTERPOLATION.'
+          WRITE(lunout, *)' Daily input file.'
+        ELSE
+          IF(is_bcs) WRITE(lunout, *)'LINEAR TIME INTERPOLATION.'
+          IF(.NOT.is_bcs) WRITE(lunout, *)'SPLINES TIME INTERPOLATION.'
+          WRITE(lunout, *)' Input time vector: ', timeyear
+          WRITE(lunout, *)' Output time vector: from 0.5 to ', ndays - 0.5
+        END IF
+      END IF
+      ALLOCATE(champan(iip1, jjp1, ndays))
+
+      IF(lmdep==ndays_in) THEN  !--- DAILY DATA: NO     TIME INTERPOLATION
+        DO l = 1, lmdep
+          champan(1:iim, :, l) = champtime(:, :, l + 1)
+        END DO
+      ELSE IF(is_bcs) THEN      !--- BCS   DATA: LINEAR TIME INTERPOLATION
+        l = 1
+        DO k = 1, ndays
+          timeday = (REAL(k) - 0.5) * REAL(ndays_in) / ndays
+          IF(timeyear(l + 1)<timeday) l = l + 1
+          al = (timeday - timeyear(l)) / (timeyear(l + 1) - timeyear(l))
+          champan(1:iim, :, k) = champtime(1:iim, :, l) + al * (champtime(1:iim, :, l + 1) - champtime(1:iim, :, l))
+        END DO
+      ELSE                      !--- AVE   DATA: SPLINE TIME INTERPOLATION
+        skip = .false.
+        n_extrap = 0
+        ALLOCATE(yder(lmdep + 2))
+        DO j = 1, jjp1
+          DO i = 1, iim
+            yder = pchsp_95(timeyear, champtime(i, j, :), ibeg = 2, iend = 2, &
+                    vc_beg = 0., vc_end = 0.)
+            CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
+                    arth(0.5, real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
+            if (ierr < 0) call abort_physic("get_2Dfield", "", 1)
+            n_extrap = n_extrap + ierr
+          END DO
+        END DO
+        IF(n_extrap /= 0) WRITE(lunout, *) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
+        DEALLOCATE(yder)
+      END IF
+      champan(iip1, :, :) = champan(1, :, :)
+      DEALLOCATE(champtime, timeyear)
+
+      !--- Checking the result
+      DO j = 1, jjp1
+        CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
+        IF (prt_level>5) WRITE(lunout, *)' ', TRIM(title), ' at time 10 ', chmin, chmax, j
+      END DO
+
+      !--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
+      IF(mode=='SST') THEN
+        SELECT CASE(units)
+        CASE("K"); CALL msg(0, 'SST field is already in kelvins.')
+        CASE("C"); CALL msg(0, 'SST field converted from celcius degrees to kelvins.')
+        champan(:, :, :) = champan(:, :, :) + 273.15
+        END SELECT
+        CALL msg(0, 'Filtering SST: Sea Surface Temperature >= 271.38')
+        WHERE(champan<271.38) champan = 271.38
+      END IF
+
+      !--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
+      IF(mode=='SIC') THEN
+        SELECT CASE(units)
+        CASE("1"); CALL msg(0, 'SIC field already in fraction of 1')
+        CASE("%"); CALL msg(0, 'SIC field converted from percentage to fraction of 1.')
+        champan(:, :, :) = champan(:, :, :) / 100.
+        END SELECT
+        CALL msg(0, 'Filtering SIC: 0.0 <= Sea-ice <=1.0')
+        WHERE(champan>1.0) champan = 1.0
+        WHERE(champan<0.0) champan = 0.0
+      END IF
+
+      !--- 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(:, k))
+      END DO
+      DEALLOCATE(champan)
+
+    END SUBROUTINE get_2Dfield
+
+    !-------------------------------------------------------------------------------
+
+
+    !-------------------------------------------------------------------------------
+
+    SUBROUTINE start_init_orog0(lon_in, lat_in, phis, masque)
+
+      !-------------------------------------------------------------------------------
+      USE grid_noro_m, ONLY : grid_noro0
+      IMPLICIT NONE
+      !===============================================================================
+      ! Purpose:  Compute "phis" just like it would be in start_init_orog.
+      !===============================================================================
+      ! Arguments:
+      REAL, INTENT(IN) :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
+      REAL, INTENT(INOUT) :: phis(:, :), masque(:, :) ! dim (iml,jml)
+      !-------------------------------------------------------------------------------
+      ! Local variables:
+      CHARACTER(LEN = ns) :: modname = "start_init_orog0"
+      INTEGER :: fid, llm_tmp, ttm_tmp, iml, jml, iml_rel, jml_rel, itau(1)
+      REAL :: lev(1), date, dt, deg2rad
+      REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), lon_rel(:, :), relief_hi(:, :)
+      REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:), lat_rel(:, :)
+      !-------------------------------------------------------------------------------
+      iml = assert_eq(SIZE(lon_in), SIZE(phis, 1), SIZE(masque, 1), TRIM(modname) // " iml")
+      jml = assert_eq(SIZE(lat_in), SIZE(phis, 2), SIZE(masque, 2), TRIM(modname) // " jml")
+      IF(iml/=iip1) CALL abort_gcm(TRIM(modname), 'iml/=iip1', 1)
+      IF(jml/=jjp1) CALL abort_gcm(TRIM(modname), 'jml/=jjp1', 1)
+      pi = 2.0 * ASIN(1.0); deg2rad = pi / 180.0
+      IF(ANY(phis/=-99999.)) RETURN                  !--- phis ALREADY KNOWN
+
+      !--- HIGH RESOLUTION OROGRAPHY
+      CALL flininfo(frelf, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
+
+      ALLOCATE(lat_rel(iml_rel, jml_rel), lon_rel(iml_rel, jml_rel))
+      CALL flinopen(frelf, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel, &
+              lev, ttm_tmp, itau, date, dt, fid)
+      ALLOCATE(relief_hi(iml_rel, jml_rel))
+      CALL flinget(fid, vrel, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
+      CALL flinclo(fid)
+
+      !--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
+      ALLOCATE(lon_ini(iml_rel), lat_ini(jml_rel))
+      lon_ini(:) = lon_rel(:, 1); IF(MAXVAL(lon_rel)>pi) lon_ini = lon_ini * deg2rad
+      lat_ini(:) = lat_rel(1, :); IF(MAXVAL(lat_rel)>pi) lat_ini = lat_ini * deg2rad
+
+      !--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
+      ALLOCATE(lon_rad(iml_rel), lat_rad(jml_rel))
+      CALL conf_dat2d(vrel, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
+      DEALLOCATE(lon_ini, lat_ini)
+
+      !--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0
+      WRITE(lunout, *)
+      WRITE(lunout, *)'*** Compute surface geopotential ***'
+
+      !--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS
+      CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque)
+      phis = phis * 9.81
+      phis(iml, :) = phis(1, :)
+      DEALLOCATE(relief_hi, lon_rad, lat_rad)
+
+    END SUBROUTINE start_init_orog0
+
+    !-------------------------------------------------------------------------------
+
+
+    !-------------------------------------------------------------------------------
+
+    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
+
+    !-------------------------------------------------------------------------------
+
+
+    !-------------------------------------------------------------------------------
+
+    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
+      !-------------------------------------------------------------------------------
+      IF(ncres/=nf90_noerr) THEN
+        WRITE(lunout, *)'Problem with file ' // TRIM(fnam) // ' in routine limit_netcdf.'
+        CALL abort_physic('limit_netcdf', NF90_STRERROR(ncres), 1)
+      END IF
+
+    END SUBROUTINE ncerr
+
+    !-------------------------------------------------------------------------------
+
+
+    !-------------------------------------------------------------------------------
+
+    SUBROUTINE strclean(s)
+
+      !-------------------------------------------------------------------------------
+      IMPLICIT NONE
+      !-------------------------------------------------------------------------------
+      ! Purpose: Remove tail null characters from the input string.
+      !-------------------------------------------------------------------------------
+      ! Parameters:
+      CHARACTER(LEN = *), INTENT(INOUT) :: s
+      !-------------------------------------------------------------------------------
+      ! Local variable:
+      INTEGER :: k
+      !-------------------------------------------------------------------------------
+      k = LEN_TRIM(s); DO WHILE(ICHAR(s(k:k))==0); s(k:k) = ' '; k = LEN_TRIM(s);
+      END DO
+
+    END SUBROUTINE strclean
+
+    !-------------------------------------------------------------------------------
+
+
+    !-------------------------------------------------------------------------------
+
+    FUNCTION is_in(s1, s2) RESULT(res)
+
+      !-------------------------------------------------------------------------------
+      IMPLICIT NONE
+      !-------------------------------------------------------------------------------
+      ! Purpose: Check wether s1 is present in the s2(:) list (case insensitive).
+      !-------------------------------------------------------------------------------
+      ! Arguments:
+      CHARACTER(LEN = *), INTENT(IN) :: s1, s2(:)
+      LOGICAL :: res
+      !-------------------------------------------------------------------------------
+      res = .FALSE.; DO k = 1, SIZE(s2); res = res.OR.strLow(s1)==strLow(s2(k));
+      END DO
+
+    END FUNCTION is_in
+
+    !-------------------------------------------------------------------------------
+
+
+    !-------------------------------------------------------------------------------
+
+    ELEMENTAL FUNCTION strLow(s) RESULT(res)
+
+      !-------------------------------------------------------------------------------
+      IMPLICIT NONE
+      !-------------------------------------------------------------------------------
+      ! Purpose: Lower case conversion.
+      !-------------------------------------------------------------------------------
+      ! Arguments:
+      CHARACTER(LEN = *), INTENT(IN) :: s
+      CHARACTER(LEN = ns) :: res
+      !-------------------------------------------------------------------------------
+      ! Local variable:
+      INTEGER :: k, ix
+      !-------------------------------------------------------------------------------
+      res = s
+      DO k = 1, LEN(s); ix = IACHAR(s(k:k))
+      IF(64<ix.AND.ix<91) res(k:k) = ACHAR(ix + 32)
+      END DO
+
+    END FUNCTION strLow
+
+    !-------------------------------------------------------------------------------
+
+  END SUBROUTINE limit_netcdf
+
+END MODULE limit
+
+!*******************************************************************************
+
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/create_etat0_unstruct_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/create_etat0_unstruct_mod.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/create_etat0_unstruct_mod.F90	(revision 5100)
@@ -23,5 +23,5 @@
   SUBROUTINE init_create_etat0_unstruct
   USE lmdz_xios
-  USE netcdf, ONLY: NF90_NOWRITE,nf90_close,nf90_noerr,nf90_open
+  USE netcdf, ONLY: nf90_nowrite,nf90_close,nf90_noerr,nf90_open
   USE mod_phys_lmdz_para
   IMPLICIT NONE
@@ -31,20 +31,20 @@
     IF (is_omp_master) THEN
 
-      IF (NF90_OPEN("ocean_fraction.nc", NF90_NOWRITE, file_id)==nf90_noerr) THEN
+      IF (nf90_open("ocean_fraction.nc", nf90_nowrite, file_id)==nf90_noerr) THEN
         CALL xios_set_file_attr("frac_ocean",enabled=.TRUE.)
         CALL xios_set_field_attr("mask",field_ref="frac_ocean_read")
-        iret=NF90_CLOSE(file_id)
-      ELSE IF (NF90_OPEN("land_water_0.05.nc", NF90_NOWRITE, file_id)==nf90_noerr) THEN
+        iret=nf90_close(file_id)
+      ELSE IF (nf90_open("land_water_0.05.nc", nf90_nowrite, file_id)==nf90_noerr) THEN
         CALL xios_set_file_attr("land_water",name="land_water_0.05",enabled=.TRUE.)
         CALL xios_set_field_attr("mask",field_ref="land_water")
-        iret=NF90_CLOSE(file_id)
-      ELSE IF (NF90_OPEN("land_water_0.25.nc", NF90_NOWRITE, file_id)==nf90_noerr) THEN
+        iret=nf90_close(file_id)
+      ELSE IF (nf90_open("land_water_0.25.nc", nf90_nowrite, file_id)==nf90_noerr) THEN
         CALL xios_set_file_attr("land_water",name="land_water_0.25",enabled=.TRUE.)
         CALL xios_set_field_attr("mask",field_ref="land_water")
-        iret=NF90_CLOSE(file_id)
-      ELSE IF (NF90_OPEN("land_water_0.50.nc", NF90_NOWRITE, file_id)==nf90_noerr) THEN
+        iret=nf90_close(file_id)
+      ELSE IF (nf90_open("land_water_0.50.nc", nf90_nowrite, file_id)==nf90_noerr) THEN
         CALL xios_set_file_attr("land_water",name="land_water_0.50",enabled=.TRUE.)
         CALL xios_set_field_attr("mask",field_ref="land_water")
-        iret=NF90_CLOSE(file_id)
+        iret=nf90_close(file_id)
       ENDIF
 
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz1d.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz1d.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz1d.F90	(revision 5100)
@@ -2,5 +2,4 @@
 ! $Id$
 
-!#ifdef CPP_1D
 !#include "../dyn3d/mod_const_mpi.F90"
 !#include "../dyn3d_common/control_mod.F90"
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/grid_noro_m.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/grid_noro_m.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/grid_noro_m.F90	(revision 5100)
@@ -435,7 +435,7 @@
 ! Purpose: Read parameters usually determined with grid_noro from a file.
 !===============================================================================
-  USE netcdf, ONLY: NF90_OPEN,  nf90_inq_dimid, nf90_inquire_dimension,        &
-        nf90_noerr, NF90_CLOSE, NF90_INQ_VARID, nf90_get_var, NF90_STRERROR,   &
-        NF90_NOWRITE
+  USE netcdf, ONLY: nf90_open,  nf90_inq_dimid, nf90_inquire_dimension,        &
+        nf90_noerr, nf90_close, NF90_INQ_VARID, nf90_get_var, NF90_STRERROR,   &
+        nf90_nowrite
   IMPLICIT NONE 
 !-------------------------------------------------------------------------------
@@ -471,5 +471,5 @@
   masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
   WRITE(lunout,*)'Masque lu: ',masque_lu
-  CALL ncerr(NF90_OPEN(fname,NF90_NOWRITE,fid))
+  CALL ncerr(nf90_open(fname,nf90_nowrite,fid))
   CALL check_dim('x','longitude',x(1:imar))
   CALL check_dim('y','latitude' ,y(1:jmar))
@@ -483,5 +483,5 @@
   zpic=zmea+2*zstd
   zval=MAX(0.,zmea-2.*zstd)
-  CALL ncerr(NF90_CLOSE(fid))
+  CALL ncerr(nf90_close(fid))
   WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
   WRITE(lunout,*)'  ST. DEV.:' ,MAXVAL(zstd)
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/iostart.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/iostart.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/iostart.F90	(revision 5100)
@@ -37,5 +37,5 @@
 
     IF (is_mpi_root .AND. is_omp_root) THEN
-      ierr = NF90_OPEN (filename, NF90_NOWRITE,nid_start)
+      ierr = nf90_open (filename, nf90_nowrite,nid_start)
       IF (ierr/=nf90_noerr) THEN
         write(6,*)' Pb d''ouverture du fichier '//filename
@@ -54,5 +54,5 @@
 
     IF (is_mpi_root .AND. is_omp_root) THEN
-        ierr = NF90_CLOSE (nid_start)
+        ierr = nf90_close (nid_start)
     ENDIF
 
@@ -312,5 +312,5 @@
     
     IF (is_master) THEN
-      ierr = NF90_CREATE(filename, IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), &
+      ierr = nf90_create(filename, IOR(nf90_clobber,nf90_64bit_offset), &
                           nid_restart)
       IF (ierr/=nf90_noerr) THEN
@@ -320,12 +320,12 @@
       ENDIF
 
-      ierr = NF90_PUT_ATT (nid_restart, NF90_GLOBAL, "title","Fichier redemmarage physique")
-
-      ierr = NF90_DEF_DIM (nid_restart, "index", length, idim1)
-      ierr = NF90_DEF_DIM (nid_restart, "points_physiques", klon_glo, idim2)
-      ierr = NF90_DEF_DIM (nid_restart, "horizon_vertical", klon_glo*klev, idim3)
-      ierr = NF90_DEF_DIM (nid_restart, "horizon_klevp1", klon_glo*klevp1, idim4)
-
-!      ierr = NF90_ENDDEF(nid_restart)
+      ierr = nf90_put_att (nid_restart, nf90_global, "title","Fichier redemmarage physique")
+
+      ierr = nf90_def_dim (nid_restart, "index", length, idim1)
+      ierr = nf90_def_dim (nid_restart, "points_physiques", klon_glo, idim2)
+      ierr = nf90_def_dim (nid_restart, "horizon_vertical", klon_glo*klev, idim3)
+      ierr = nf90_def_dim (nid_restart, "horizon_klevp1", klon_glo*klevp1, idim4)
+
+!      ierr = nf90_enddef(nid_restart)
     ENDIF
 
@@ -338,5 +338,5 @@
     INTEGER          :: ierr
 
-    IF (is_master) ierr = NF90_ENDDEF(nid_restart)
+    IF (is_master) ierr = nf90_enddef(nid_restart)
  
   END SUBROUTINE enddef_restartphy
@@ -348,5 +348,5 @@
     INTEGER          :: ierr
 
-    IF (is_master) ierr = NF90_CLOSE (nid_restart)
+    IF (is_master) ierr = nf90_close (nid_restart)
  
   END SUBROUTINE close_restartphy
@@ -426,7 +426,7 @@
          
 !      ierr = NF90_REDEF (nid_restart)
-      ierr = NF90_DEF_VAR (nid_restart, field_name, nf90_format,(/ idim /),nvarid)
-      IF (LEN_TRIM(title) > 0) ierr = NF90_PUT_ATT (nid_restart,nvarid,"title", title)
-!      ierr = NF90_ENDDEF(nid_restart)
+      ierr = nf90_def_var (nid_restart, field_name, nf90_format,(/ idim /),nvarid)
+      IF (LEN_TRIM(title) > 0) ierr = nf90_put_att (nid_restart,nvarid,"title", title)
+!      ierr = nf90_enddef(nid_restart)
      ENDIF
 
@@ -454,5 +454,5 @@
 
        ierr = NF90_INQ_VARID(nid_restart, field_name, nvarid)
-       ierr = NF90_PUT_VAR(nid_restart,nvarid,RESHAPE(field_glo,(/klon_glo*field_size/)))
+       ierr = nf90_put_var(nid_restart,nvarid,RESHAPE(field_glo,(/klon_glo*field_size/)))
       ENDIF
    ENDIF
@@ -536,12 +536,12 @@
 !      ierr = NF90_REDEF (nid_restart)
 
-        ierr = NF90_DEF_VAR (nid_restart, var_name, nf90_format,(/ idim1 /),nvarid)
-        IF (LEN_TRIM(title)>0) ierr = NF90_PUT_ATT (nid_restart,nvarid,"title", title)
-!      ierr = NF90_ENDDEF(nid_restart)
+        ierr = nf90_def_var (nid_restart, var_name, nf90_format,(/ idim1 /),nvarid)
+        IF (LEN_TRIM(title)>0) ierr = nf90_put_att (nid_restart,nvarid,"title", title)
+!      ierr = nf90_enddef(nid_restart)
 
     ! second pass : write     
       ELSE IF (pass==2) THEN 
         ierr = NF90_INQ_VARID(nid_restart, var_name, nvarid)
-        ierr = NF90_PUT_VAR(nid_restart,nvarid,var)
+        ierr = nf90_put_var(nid_restart,nvarid,var)
       ENDIF
     ENDIF
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/iotd_ecrit.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/iotd_ecrit.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/iotd_ecrit.F90	(revision 5100)
@@ -109,5 +109,5 @@
 
 
-           ierr= NF90_PUT_VAR(nid,varid,date,[ntime])
+           ierr= nf90_put_var(nid,varid,date,[ntime])
 
 !          print*,'date ',date,ierr,nid
@@ -173,5 +173,5 @@
 
 
-      ierr= NF90_PUT_VAR(nid,varid,zx,corner,edges)
+      ierr= nf90_put_var(nid,varid,zx,corner,edges)
 
       if (ierr/=nf90_noerr) then
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/limit_read_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/limit_read_mod.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/limit_read_mod.F90	(revision 5100)
@@ -234,5 +234,5 @@
 !$OMP MASTER  ! Only master thread
        IF (is_mpi_root) THEN ! Only master processus
-          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
+          ierr = nf90_open ('limit.nc', nf90_nowrite, nid)
           IF (ierr /= nf90_noerr) CALL abort_physic(modname,&
                'Pb d''ouverture du fichier de conditions aux limites',1)
@@ -240,5 +240,5 @@
           !--- WARNING IF CALENDAR IS KNOWN AND DOES NOT MATCH THE ONE OF LMDZ
           ierr=NF90_INQ_VARID(nid, 'TEMPS', nvarid)
-          ierr=NF90_GET_ATT(nid, nvarid, 'calendar', calendar)
+          ierr=nf90_get_att(nid, nvarid, 'calendar', calendar)
           IF(ierr==nf90_noerr.AND.calendar/=calend.AND.prt_level>=1) THEN
              WRITE(lunout,*)'BEWARE: gcm and limit.nc calendars differ: '
@@ -269,5 +269,5 @@
           IF(nn/=klon_glo) CALL abort_physic(modname,abort_message,1)
 
-          ierr = NF90_CLOSE(nid)
+          ierr = nf90_close(nid)
           IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Pb when closing file', 1)
        END IF ! is_mpi_root
@@ -332,5 +332,5 @@
        IF (is_mpi_root) THEN ! Only master processus!
 
-          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
+          ierr = nf90_open ('limit.nc', nf90_nowrite, nid)
           IF (ierr /= nf90_noerr) CALL abort_physic(modname,&
                'Pb d''ouverture du fichier de conditions aux limites',1)
@@ -426,5 +426,5 @@
 
 !****************************************************************************************
-          ierr = NF90_CLOSE(nid)
+          ierr = nf90_close(nid)
           IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Pb when closing file', 1)
        ENDIF ! is_mpi_root
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/limit_slab.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/limit_slab.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/limit_slab.F90	(revision 5100)
@@ -61,5 +61,5 @@
         read_siv=.TRUE.
         
-        ierr = NF90_OPEN ('limit_slab.nc', NF90_NOWRITE, nid)
+        ierr = nf90_open ('limit_slab.nc', nf90_nowrite, nid)
         IF (ierr /= nf90_noerr) THEN
             PRINT *,'LIMIT_SLAB file not found'
@@ -145,5 +145,5 @@
 
 !****************************************************************************************
-        ierr = NF90_CLOSE(nid)
+        ierr = nf90_close(nid)
         IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Pb when closing file', 1)
         END IF ! Read File 
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/mo_simple_plumes.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/mo_simple_plumes.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/mo_simple_plumes.F90	(revision 5100)
@@ -88,5 +88,5 @@
     IF (is_mpi_root.AND.is_omp_root) THEN
 
-       iret = nf90_open("MACv2.0-SP_v1.nc", NF90_NOWRITE, ncid)
+       iret = nf90_open("MACv2.0-SP_v1.nc", nf90_nowrite, ncid)
        IF (iret /= nf90_noerr) THEN
           abort_message='NetCDF File not opened'
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/open_climoz_m.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/open_climoz_m.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/open_climoz_m.F90	(revision 5100)
@@ -53,5 +53,5 @@
     press_in_cen = press_in_cen * 100.
     nlev = SIZE(press_in_cen)
-    CALL NF95_INQ_VARID(ncID, "time", varID)
+    CALL nf95_inq_varid(ncID, "time", varID)
     CALL NF95_GW_VAR(ncid, varid, time_in)
     ntim = SIZE(time_in)
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/read_map2D.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/read_map2D.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/read_map2D.F90	(revision 5100)
@@ -30,5 +30,5 @@
 ! Read variable from file. Done by master process MPI and master thread OpenMP
   IF (is_mpi_root .AND. is_omp_root) THEN
-     ierr = NF90_OPEN(trim(filename), NF90_NOWRITE, nid)
+     ierr = nf90_open(trim(filename), nf90_nowrite, nid)
      IF (ierr /= nf90_noerr) CALL write_err_mess('Problem in opening file')
 
@@ -41,5 +41,5 @@
      IF (ierr /= nf90_noerr) CALL write_err_mess('Problem in reading varaiable')
 
-     ierr = NF90_CLOSE(nid)
+     ierr = nf90_close(nid)
      IF (ierr /= nf90_noerr) CALL write_err_mess('Problem in closing file')
 
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/readaerosol_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/readaerosol_mod.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/readaerosol_mod.F90	(revision 5100)
@@ -200,6 +200,6 @@
         IF (is_mpi_root) THEN
     
-          IF (nf90_open(TRIM(file_aerosol), NF90_NOWRITE, ncid) /= nf90_noerr) THEN
-            CALL check_err( nf90_open(TRIM(file_so4), NF90_NOWRITE, ncid), "pb open "//trim(file_so4) )
+          IF (nf90_open(TRIM(file_aerosol), nf90_nowrite, ncid) /= nf90_noerr) THEN
+            CALL check_err( nf90_open(TRIM(file_so4), nf90_nowrite, ncid), "pb open "//trim(file_so4) )
           ENDIF
 
@@ -351,5 +351,5 @@
   
        WRITE(lunout,*) 'reading variable ',TRIM(varname),' in file ', TRIM(fname)
-       CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid), "pb open "//trim(fname) )
+       CALL check_err( nf90_open(TRIM(fname), nf90_nowrite, ncid), "pb open "//trim(fname) )
 
 
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/regr_horiz_time_climoz_m.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/regr_horiz_time_climoz_m.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/regr_horiz_time_climoz_m.F90	(revision 5100)
@@ -1,723 +1,722 @@
 MODULE regr_horiz_time_climoz_m
 
-  USE interpolation,     ONLY: locate
-  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, grid_type, unstructured
-  USE nrtype,            ONLY: pi
-  USE netcdf,   ONLY: NF90_CLOBBER, NF90_FLOAT,     NF90_OPEN,   &
-                      NF90_NOWRITE, nf90_noerr,     NF90_GET_ATT, NF90_GLOBAL
-  USE netcdf95, ONLY: NF95_DEF_DIM, NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION,    &
-       NF95_DEF_VAR, NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, &
-       NF95_OPEN,  NF95_CREATE,  NF95_GET_ATT,   NF95_GW_VAR,  nf95_get_var,  &
-       NF95_CLOSE, NF95_ENDDEF,  NF95_PUT_ATT,   NF95_PUT_VAR, NF95_COPY_ATT
-  USE print_control_mod, ONLY: lunout
+  USE interpolation, ONLY : locate
+  USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat, grid_type, unstructured
+  USE nrtype, ONLY : pi
+  USE netcdf, ONLY : nf90_clobber, nf90_float, nf90_open, &
+          nf90_nowrite, nf90_noerr, nf90_get_att, nf90_global
+  USE netcdf95, ONLY : nf95_def_dim, nf95_inq_dimid, nf95_inquire_dimension, &
+          nf95_def_var, nf95_inq_varid, NF95_INQUIRE_VARIABLE, &
+          NF95_OPEN, NF95_CREATE, NF95_GET_ATT, NF95_GW_VAR, nf95_get_var, &
+          NF95_CLOSE, NF95_ENDDEF, NF95_PUT_ATT, NF95_PUT_VAR, NF95_COPY_ATT
+  USE print_control_mod, ONLY : lunout
   USE dimphy
   IMPLICIT NONE
   PRIVATE
   PUBLIC :: regr_horiz_time_climoz
-  REAL, PARAMETER :: deg2rad=pi/180.
-  CHARACTER(LEN=13), PARAMETER :: vars_in(2)=['tro3         ','tro3_daylight']
+  REAL, PARAMETER :: deg2rad = pi / 180.
+  CHARACTER(LEN = 13), PARAMETER :: vars_in(2) = ['tro3         ', 'tro3_daylight']
 
   INTEGER :: nlat_ou, nlon_ou
   REAL, ALLOCATABLE :: latitude_glo(:)
-!$OMP THREADPRIVATE(latitude_glo)
+  !$OMP THREADPRIVATE(latitude_glo)
   INTEGER, ALLOCATABLE :: ind_cell_glo_glo(:)
-!$OMP THREADPRIVATE(ind_cell_glo_glo)
+  !$OMP THREADPRIVATE(ind_cell_glo_glo)
 
 CONTAINS
 
-!-------------------------------------------------------------------------------
-
-SUBROUTINE regr_horiz_time_climoz(read_climoz,interpt)
-
-!-------------------------------------------------------------------------------
-! Purpose: Regrid horizontally and in time zonal or 3D ozone climatologies.
-!   * Read ozone climatology from netcdf file
-!   * Regrid it horizontaly to LMDZ grid (quasi-conservative method)
-!   * If interpt=T, interpolate linearly in time (one record each day)
-!     If interpt=F, keep original time sampling  (14 months).
-!   * Save it to a new netcdf file.
-!-------------------------------------------------------------------------------
-! Remarks:
-!   * Up to 2 variables treated: "tro3" and "tro3_daylight" (if read_climoz=2)
-!   * Input fields coordinates: (longitudes, latitudes, pressure_levels, time)
-!   * Output grid cells centers coordinates given by [rlonv,] rlatu.
-!   * Output grid cells edges   coordinates given by [rlonu,] rlatv.
-!   * Input file [longitudes and] latitudes given in degrees.
-!   * Input file pressure levels are given in Pa or hPa.
-!   * All coordinates variables are stricly monotonic.
-!   * Monthly fields are interpolated linearly in time to get daily values.
-!   * Fields are known at the middle of the months, so interpolation requires an
-!     additional record both for 1st half of january and 2nd half of december:
-!     - For a 14-records "climoz.nc": records 1 and 14.
-!     - For 12-records files:
-!       record 12 of "climoz_m.nc" if available, or record 1  of "climoz.nc".
-!       record 1  of "climoz_p.nc" if available, or record 12 of "climoz.nc".
-!   * Calendar is taken into account to get one record each day (not 360 always).
-!   * Missing values are filled in from sky to ground by copying lowest valid one.
-!     Attribute "missing_value" or "_FillValue" must be present in input file.
-!-------------------------------------------------------------------------------
-  USE assert_m,           ONLY: assert
-  USE cal_tools_m,        ONLY: year_len, mid_month
-!!  USE control_mod,        ONLY: anneeref
-  USE time_phylmdz_mod,   ONLY: annee_ref
-  USE ioipsl,             ONLY: ioget_year_len, ioget_calendar
-  USE regr_conserv_m,     ONLY: regr_conserv
-  USE regr_lint_m,        ONLY: regr_lint
-  USE regular_lonlat_mod, ONLY: boundslon_reg, boundslat_reg, south, west, east
-  USE slopes_m,           ONLY: slopes
-  USE lmdz_xios
-  USE mod_phys_lmdz_para, ONLY: is_mpi_root, is_master, is_omp_master, gather, gather_mpi, bcast_mpi, klon_mpi
-  USE geometry_mod, ONLY : latitude_deg, ind_cell_glo
-  USE mod_grid_phy_lmdz, ONLY: klon_glo
-
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER, INTENT(IN) :: read_climoz ! read ozone climatology, 1 or 2
-!                         1: read a single ozone climatology used day and night
-!                         2: same + read also a daylight climatology
-  LOGICAL, INTENT(IN) :: interpt     ! TRUE  => daily interpolation
-                                     ! FALSE => no interpolation (14 months)
-!-------------------------------------------------------------------------------
-! Local variables:
-
-!--- Input files variables
-  INTEGER :: nlon_in                       ! Number of longitudes
-  INTEGER :: nlat_in                       ! Number of latitudes
-  INTEGER :: nlev_in                       ! Number of pressure levels
-  INTEGER :: nmth_in                       ! Number of months
-  REAL, ALLOCATABLE:: lon_in(:)           ! Longitudes   (ascending order, rad)
-  REAL, ALLOCATABLE:: lat_in(:)           ! Latitudes    (ascending order, rad)
-  REAL, ALLOCATABLE:: lev_in(:)           ! Pressure levels (ascen. order, hPa)
-  REAL, ALLOCATABLE :: lon_in_edge(:)      ! Longitude intervals edges
-                                           !              (ascending order,  / )
-  REAL, ALLOCATABLE :: sinlat_in_edge(:)   ! Sinus of latitude intervals edges
-                                           !              (ascending order,  / )
-  LOGICAL :: ldec_lon, ldec_lat, ldec_lev  ! Decreasing order in input file
-  CHARACTER(LEN=20) :: cal_in              ! Calendar
-  REAL, ALLOCATABLE :: o3_in3(:,:,:,:,:)   ! Ozone climatologies
-  REAL, ALLOCATABLE :: o3_in3bis(:,:,:,:,:)   ! Ozone climatologies
-  REAL, ALLOCATABLE :: o3_in2  (:,:,:,:)   ! Ozone climatologies
-  REAL, ALLOCATABLE :: o3_in2bis(:,:,:,:,:)   ! Ozone climatologies
-  ! last index: 1 for the day-night average, 2 for the daylight field.
-  REAL :: NaN
-
-!--- Partially or totally regridded variables      (:,:,nlev_in,:,read_climoz)
-  REAL, ALLOCATABLE :: o3_regr_lon   (:,:,:,:,:) ! (nlon_ou,nlat_in,:,0:13   ,:)
-  REAL, ALLOCATABLE :: o3_regr_lonlat(:,:,:,:,:) ! (nlon_ou,nlat_ou,:,0:13   ,:)
-  REAL, ALLOCATABLE :: o3_out3       (:,:,:,:,:) ! (nlon_ou,nlat_ou,:,ntim_ou,:)
-  REAL, ALLOCATABLE :: o3_out3_glo   (:,:,:,:) !   (nbp_lat,:,ntim_ou,:)
-  REAL, ALLOCATABLE :: o3_regr_lat     (:,:,:,:) !         (nlat_in,:,0:13   ,:)
-  REAL, ALLOCATABLE :: o3_out2         (:,:,:,:) !         (nlat_ou,:,ntim_ou,:)
-  REAL, ALLOCATABLE :: o3_out2_glo     (:,:,:,:) !         (nbp_lat,:,ntim_ou,:)
-  REAL, ALLOCATABLE :: o3_out          (:,:,:,:) !         (nbp_lat,:,ntim_ou,:)
-! Dimension number  | Interval                | Contains  | For variables:
-!   1 (longitude)   | [rlonu(i-1), rlonu(i)]  | rlonv(i)  | all
-!   2 (latitude)    | [rlatv(j), rlatv(j-1)]  | rlatu(j)  | all but o3_regr_lon
-!   3 (press level) |                         |   lev(k)  | all
-! Note that rlatv(0)=pi/2 and rlatv(nlat_ou)=-pi/2.
-! Dimension 4 is: month number                             (all vars but o3_out)
-!                 days elapsed since Jan. 1st 0h at mid-day (o3_out only)
-  REAL, ALLOCATABLE :: v1(:)
-
-!--- For NetCDF:
-  INTEGER :: fID_in_m, fID_in, levID_ou, dimid, vID_in(read_climoz), ntim_ou
-  INTEGER :: fID_in_p, fID_ou, timID_ou, varid, vID_ou(read_climoz), ndims, ncerr
-  INTEGER, ALLOCATABLE :: dIDs(:)
-  CHARACTER(LEN=20) :: cal_ou     !--- Calendar; no time inter => same as input
-  CHARACTER(LEN=80) :: press_unit !--- Pressure unit
-  REAL    :: tmidmonth(0:13)      !--- Elapsed days since Jan-1 0h at mid-months
-                                  ! Additional records 0, 13 for interpolation
-  REAL, ALLOCATABLE :: tmidday(:) !--- Output times (mid-days since Jan 1st 0h)
-  LOGICAL :: lprev, lnext         !--- Flags: previous/next files are present
-  LOGICAL :: l3D, l2D             !--- Flag:  input fields are 3D or zonal
-  INTEGER :: ii, i, j, k, l, m, dln, ib, ie, iv, dx1, dx2
-  INTEGER, ALLOCATABLE :: sta(:), cnt(:)
-  CHARACTER(LEN=80) :: sub, dim_nam, msg
-  REAL :: null_array(0)
-  LOGICAL,SAVE :: first=.TRUE.
-!$OMP THREADPRIVATE(first)  
-  REAL, ALLOCATABLE :: test_o3_in(:,:)
-  REAL, ALLOCATABLE :: test_o3_out(:)
-
-
-  IF (grid_type==unstructured) THEN 
-    IF (first) THEN
-      IF (is_master) THEN
-        ALLOCATE(latitude_glo(klon_glo))
-        ALLOCATE(ind_cell_glo_glo(klon_glo))
-      ELSE
-        ALLOCATE(latitude_glo(0))
-        ALLOCATE(ind_cell_glo_glo(0))
+  !-------------------------------------------------------------------------------
+
+  SUBROUTINE regr_horiz_time_climoz(read_climoz, interpt)
+
+    !-------------------------------------------------------------------------------
+    ! Purpose: Regrid horizontally and in time zonal or 3D ozone climatologies.
+    !   * Read ozone climatology from netcdf file
+    !   * Regrid it horizontaly to LMDZ grid (quasi-conservative method)
+    !   * If interpt=T, interpolate linearly in time (one record each day)
+    !     If interpt=F, keep original time sampling  (14 months).
+    !   * Save it to a new netcdf file.
+    !-------------------------------------------------------------------------------
+    ! Remarks:
+    !   * Up to 2 variables treated: "tro3" and "tro3_daylight" (if read_climoz=2)
+    !   * Input fields coordinates: (longitudes, latitudes, pressure_levels, time)
+    !   * Output grid cells centers coordinates given by [rlonv,] rlatu.
+    !   * Output grid cells edges   coordinates given by [rlonu,] rlatv.
+    !   * Input file [longitudes and] latitudes given in degrees.
+    !   * Input file pressure levels are given in Pa or hPa.
+    !   * All coordinates variables are stricly monotonic.
+    !   * Monthly fields are interpolated linearly in time to get daily values.
+    !   * Fields are known at the middle of the months, so interpolation requires an
+    !     additional record both for 1st half of january and 2nd half of december:
+    !     - For a 14-records "climoz.nc": records 1 and 14.
+    !     - For 12-records files:
+    !       record 12 of "climoz_m.nc" if available, or record 1  of "climoz.nc".
+    !       record 1  of "climoz_p.nc" if available, or record 12 of "climoz.nc".
+    !   * Calendar is taken into account to get one record each day (not 360 always).
+    !   * Missing values are filled in from sky to ground by copying lowest valid one.
+    !     Attribute "missing_value" or "_FillValue" must be present in input file.
+    !-------------------------------------------------------------------------------
+    USE assert_m, ONLY : assert
+    USE cal_tools_m, ONLY : year_len, mid_month
+    !!  USE control_mod,        ONLY: anneeref
+    USE time_phylmdz_mod, ONLY : annee_ref
+    USE ioipsl, ONLY : ioget_year_len, ioget_calendar
+    USE regr_conserv_m, ONLY : regr_conserv
+    USE regr_lint_m, ONLY : regr_lint
+    USE regular_lonlat_mod, ONLY : boundslon_reg, boundslat_reg, south, west, east
+    USE slopes_m, ONLY : slopes
+    USE lmdz_xios
+    USE mod_phys_lmdz_para, ONLY : is_mpi_root, is_master, is_omp_master, gather, gather_mpi, bcast_mpi, klon_mpi
+    USE geometry_mod, ONLY : latitude_deg, ind_cell_glo
+    USE mod_grid_phy_lmdz, ONLY : klon_glo
+
+    !-------------------------------------------------------------------------------
+    ! Arguments:
+    INTEGER, INTENT(IN) :: read_climoz ! read ozone climatology, 1 or 2
+    !                         1: read a single ozone climatology used day and night
+    !                         2: same + read also a daylight climatology
+    LOGICAL, INTENT(IN) :: interpt     ! TRUE  => daily interpolation
+    ! FALSE => no interpolation (14 months)
+    !-------------------------------------------------------------------------------
+    ! Local variables:
+
+    !--- Input files variables
+    INTEGER :: nlon_in                       ! Number of longitudes
+    INTEGER :: nlat_in                       ! Number of latitudes
+    INTEGER :: nlev_in                       ! Number of pressure levels
+    INTEGER :: nmth_in                       ! Number of months
+    REAL, ALLOCATABLE :: lon_in(:)           ! Longitudes   (ascending order, rad)
+    REAL, ALLOCATABLE :: lat_in(:)           ! Latitudes    (ascending order, rad)
+    REAL, ALLOCATABLE :: lev_in(:)           ! Pressure levels (ascen. order, hPa)
+    REAL, ALLOCATABLE :: lon_in_edge(:)      ! Longitude intervals edges
+    !              (ascending order,  / )
+    REAL, ALLOCATABLE :: sinlat_in_edge(:)   ! Sinus of latitude intervals edges
+    !              (ascending order,  / )
+    LOGICAL :: ldec_lon, ldec_lat, ldec_lev  ! Decreasing order in input file
+    CHARACTER(LEN = 20) :: cal_in              ! Calendar
+    REAL, ALLOCATABLE :: o3_in3(:, :, :, :, :)   ! Ozone climatologies
+    REAL, ALLOCATABLE :: o3_in3bis(:, :, :, :, :)   ! Ozone climatologies
+    REAL, ALLOCATABLE :: o3_in2  (:, :, :, :)   ! Ozone climatologies
+    REAL, ALLOCATABLE :: o3_in2bis(:, :, :, :, :)   ! Ozone climatologies
+    ! last index: 1 for the day-night average, 2 for the daylight field.
+    REAL :: NaN
+
+    !--- Partially or totally regridded variables      (:,:,nlev_in,:,read_climoz)
+    REAL, ALLOCATABLE :: o3_regr_lon   (:, :, :, :, :) ! (nlon_ou,nlat_in,:,0:13   ,:)
+    REAL, ALLOCATABLE :: o3_regr_lonlat(:, :, :, :, :) ! (nlon_ou,nlat_ou,:,0:13   ,:)
+    REAL, ALLOCATABLE :: o3_out3       (:, :, :, :, :) ! (nlon_ou,nlat_ou,:,ntim_ou,:)
+    REAL, ALLOCATABLE :: o3_out3_glo   (:, :, :, :) !   (nbp_lat,:,ntim_ou,:)
+    REAL, ALLOCATABLE :: o3_regr_lat     (:, :, :, :) !         (nlat_in,:,0:13   ,:)
+    REAL, ALLOCATABLE :: o3_out2         (:, :, :, :) !         (nlat_ou,:,ntim_ou,:)
+    REAL, ALLOCATABLE :: o3_out2_glo     (:, :, :, :) !         (nbp_lat,:,ntim_ou,:)
+    REAL, ALLOCATABLE :: o3_out          (:, :, :, :) !         (nbp_lat,:,ntim_ou,:)
+    ! Dimension number  | Interval                | Contains  | For variables:
+    !   1 (longitude)   | [rlonu(i-1), rlonu(i)]  | rlonv(i)  | all
+    !   2 (latitude)    | [rlatv(j), rlatv(j-1)]  | rlatu(j)  | all but o3_regr_lon
+    !   3 (press level) |                         |   lev(k)  | all
+    ! Note that rlatv(0)=pi/2 and rlatv(nlat_ou)=-pi/2.
+    ! Dimension 4 is: month number                             (all vars but o3_out)
+    !                 days elapsed since Jan. 1st 0h at mid-day (o3_out only)
+    REAL, ALLOCATABLE :: v1(:)
+
+    !--- For NetCDF:
+    INTEGER :: fID_in_m, fID_in, levID_ou, dimid, vID_in(read_climoz), ntim_ou
+    INTEGER :: fID_in_p, fID_ou, timID_ou, varid, vID_ou(read_climoz), ndims, ncerr
+    INTEGER, ALLOCATABLE :: dIDs(:)
+    CHARACTER(LEN = 20) :: cal_ou     !--- Calendar; no time inter => same as input
+    CHARACTER(LEN = 80) :: press_unit !--- Pressure unit
+    REAL :: tmidmonth(0:13)      !--- Elapsed days since Jan-1 0h at mid-months
+    ! Additional records 0, 13 for interpolation
+    REAL, ALLOCATABLE :: tmidday(:) !--- Output times (mid-days since Jan 1st 0h)
+    LOGICAL :: lprev, lnext         !--- Flags: previous/next files are present
+    LOGICAL :: l3D, l2D             !--- Flag:  input fields are 3D or zonal
+    INTEGER :: ii, i, j, k, l, m, dln, ib, ie, iv, dx1, dx2
+    INTEGER, ALLOCATABLE :: sta(:), cnt(:)
+    CHARACTER(LEN = 80) :: sub, dim_nam, msg
+    REAL :: null_array(0)
+    LOGICAL, SAVE :: first = .TRUE.
+    !$OMP THREADPRIVATE(first)
+    REAL, ALLOCATABLE :: test_o3_in(:, :)
+    REAL, ALLOCATABLE :: test_o3_out(:)
+
+    IF (grid_type==unstructured) THEN
+      IF (first) THEN
+        IF (is_master) THEN
+          ALLOCATE(latitude_glo(klon_glo))
+          ALLOCATE(ind_cell_glo_glo(klon_glo))
+        ELSE
+          ALLOCATE(latitude_glo(0))
+          ALLOCATE(ind_cell_glo_glo(0))
+        ENDIF
+        CALL gather(latitude_deg, latitude_glo)
+        CALL gather(ind_cell_glo, ind_cell_glo_glo)
       ENDIF
-      CALL gather(latitude_deg,  latitude_glo)
-      CALL gather(ind_cell_glo,  ind_cell_glo_glo)
-    ENDIF 
-  ENDIF
-    
-  IF (is_omp_master) THEN
-    nlat_ou=nbp_lat
-    nlon_ou=nbp_lon
-    
-   !-------------------------------------------------------------------------------
-    IF (is_mpi_root) THEN
-      sub="regr_horiz_time_climoz"
-      WRITE(lunout,*)"Call sequence information: "//TRIM(sub)
-      CALL assert(read_climoz == 1 .OR. read_climoz == 2, "regr_lat_time_climoz")
-
-      CALL  NF95_OPEN("climoz.nc"  , NF90_NOWRITE, fID_in)
-      lprev=NF90_OPEN("climoz_m.nc", NF90_NOWRITE, fID_in_m)==nf90_noerr
-      lnext=NF90_OPEN("climoz_p.nc", NF90_NOWRITE, fID_in_p)==nf90_noerr
-
-      !--- Get coordinates from the input file. Converts lon/lat in radians.
-      !    Few inversions because "regr_conserv" and gcm need ascending vectors.
-      CALL NF95_INQ_VARID(fID_in, vars_in(1), varid)
-      CALL NF95_INQUIRE_VARIABLE(fID_in, varid, dimids=dIDs, ndims=ndims)
-      l3D=ndims==4; l2D=ndims==3
-      IF(l3D) WRITE(lunout,*)"Input files contain full 3D ozone fields."
-      IF(l2D) WRITE(lunout,*)"Input files contain zonal 2D ozone fields."
-      DO i=1,ndims
-        CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), name=dim_nam, nclen=dln)
-        CALL NF95_INQ_VARID(fID_in, dim_nam, varid)
-        ii=i; IF(l2D) ii=i+1                              !--- ndims==3:NO LONGITUDE 
-        SELECT CASE(ii)
+    ENDIF
+
+    IF (is_omp_master) THEN
+      nlat_ou = nbp_lat
+      nlon_ou = nbp_lon
+
+      !-------------------------------------------------------------------------------
+      IF (is_mpi_root) THEN
+        sub = "regr_horiz_time_climoz"
+        WRITE(lunout, *)"Call sequence information: " // TRIM(sub)
+        CALL assert(read_climoz == 1 .OR. read_climoz == 2, "regr_lat_time_climoz")
+
+        CALL  NF95_OPEN("climoz.nc", nf90_nowrite, fID_in)
+        lprev = nf90_open("climoz_m.nc", nf90_nowrite, fID_in_m)==nf90_noerr
+        lnext = nf90_open("climoz_p.nc", nf90_nowrite, fID_in_p)==nf90_noerr
+
+        !--- Get coordinates from the input file. Converts lon/lat in radians.
+        !    Few inversions because "regr_conserv" and gcm need ascending vectors.
+        CALL nf95_inq_varid(fID_in, vars_in(1), varid)
+        CALL NF95_INQUIRE_VARIABLE(fID_in, varid, dimids = dIDs, ndims = ndims)
+        l3D = ndims==4; l2D = ndims==3
+        IF(l3D) WRITE(lunout, *)"Input files contain full 3D ozone fields."
+        IF(l2D) WRITE(lunout, *)"Input files contain zonal 2D ozone fields."
+        DO i = 1, ndims
+          CALL nf95_inquire_dimension(fID_in, dIDs(i), name = dim_nam, nclen = dln)
+          CALL nf95_inq_varid(fID_in, dim_nam, varid)
+          ii = i; IF(l2D) ii = i + 1                              !--- ndims==3:NO LONGITUDE
+          SELECT CASE(ii)
           CASE(1)                                         !--- LONGITUDE
             CALL NF95_GW_VAR(fID_in, varid, lon_in)
-            ldec_lon=lon_in(1)>lon_in(dln); IF(ldec_lon) lon_in=lon_in(dln:1:-1)
-            nlon_in=dln; lon_in=lon_in*deg2rad
+            ldec_lon = lon_in(1)>lon_in(dln); IF(ldec_lon) lon_in = lon_in(dln:1:-1)
+            nlon_in = dln; lon_in = lon_in * deg2rad
           CASE(2)                                         !--- LATITUDE
             CALL NF95_GW_VAR(fID_in, varid, lat_in)
-            ldec_lat=lat_in(1)>lat_in(dln); IF(ldec_lat) lat_in=lat_in(dln:1:-1)
-            nlat_in=dln; lat_in=lat_in*deg2rad
+            ldec_lat = lat_in(1)>lat_in(dln); IF(ldec_lat) lat_in = lat_in(dln:1:-1)
+            nlat_in = dln; lat_in = lat_in * deg2rad
           CASE(3)                                         !--- PRESSURE LEVELS
             CALL NF95_GW_VAR(fID_in, varid, lev_in)
-            ldec_lev=lev_in(1)>lev_in(dln); IF(ldec_lev) lev_in=lev_in(dln:1:-1)
-            nlev_in=dln
+            ldec_lev = lev_in(1)>lev_in(dln); IF(ldec_lev) lev_in = lev_in(dln:1:-1)
+            nlev_in = dln
             CALL NF95_GET_ATT(fID_in, varid, "units", press_unit)
-            k=LEN_TRIM(press_unit)
+            k = LEN_TRIM(press_unit)
             DO WHILE(ICHAR(press_unit(k:k))==0)
-              press_unit(k:k)=' '; k=LEN_TRIM(press_unit) !--- REMOVE NULL END CHAR
+              press_unit(k:k) = ' '; k = LEN_TRIM(press_unit) !--- REMOVE NULL END CHAR
             END DO
             IF(press_unit ==  "Pa") THEN
-              lev_in = lev_in/100.                        !--- CONVERT TO hPa
+              lev_in = lev_in / 100.                        !--- CONVERT TO hPa
             ELSE IF(press_unit /= "hPa") THEN
-              CALL abort_physic(sub, "the only recognized units are Pa and hPa.",1)
+              CALL abort_physic(sub, "the only recognized units are Pa and hPa.", 1)
             END IF
           CASE(4)                                         !--- TIME
-            CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), nclen=nmth_in)
-            cal_in='gregorian'
-            IF(NF90_GET_ATT(fID_in, varid, 'calendar', cal_in)/=nf90_noerr)        &
-            WRITE(lunout,*)'WARNING: missing "calendar" attribute for "'//       &
-            TRIM(dim_nam)//'" in "climoz.nc". Choosing default: "gregorian".'
-            k=LEN_TRIM(cal_in)
+            CALL nf95_inquire_dimension(fID_in, dIDs(i), nclen = nmth_in)
+            cal_in = 'gregorian'
+            IF(nf90_get_att(fID_in, varid, 'calendar', cal_in)/=nf90_noerr)        &
+                    WRITE(lunout, *)'WARNING: missing "calendar" attribute for "' // &
+                            TRIM(dim_nam) // '" in "climoz.nc". Choosing default: "gregorian".'
+            k = LEN_TRIM(cal_in)
             DO WHILE(ICHAR(cal_in(k:k))==0)
-              cal_in(k:k)=' '; k=LEN_TRIM(cal_in)         !--- REMOVE NULL END CHAR
-            END DO
-        END SELECT
-      END DO
-
-      !--- Prepare quantities for time interpolation
-      tmidmonth=mid_month(annee_ref, cal_in)
-      IF(interpt) THEN
-        ntim_ou=ioget_year_len(annee_ref)
-        ALLOCATE(tmidday(ntim_ou))
-        tmidday=[(REAL(k)-0.5,k=1,ntim_ou)]
-        CALL ioget_calendar(cal_ou)
+              cal_in(k:k) = ' '; k = LEN_TRIM(cal_in)         !--- REMOVE NULL END CHAR
+            END DO
+          END SELECT
+        END DO
+
+        !--- Prepare quantities for time interpolation
+        tmidmonth = mid_month(annee_ref, cal_in)
+        IF(interpt) THEN
+          ntim_ou = ioget_year_len(annee_ref)
+          ALLOCATE(tmidday(ntim_ou))
+          tmidday = [(REAL(k) - 0.5, k = 1, ntim_ou)]
+          CALL ioget_calendar(cal_ou)
+        ELSE
+          ntim_ou = 14
+          cal_ou = cal_in
+        END IF
+      ENDIF
+
+      IF (grid_type==unstructured) THEN
+        CALL bcast_mpi(nlon_in)
+        CALL bcast_mpi(nlat_in)
+        CALL bcast_mpi(nlev_in)
+        CALL bcast_mpi(l3d)
+        CALL bcast_mpi(tmidmonth)
+        IF(interpt) CALL bcast_mpi(tmidday)
+        CALL bcast_mpi(ntim_ou)
+
+        IF (is_mpi_root) THEN
+          CALL xios_set_domain_attr("domain_climoz", nj_glo = nlat_in, nj = nlat_in, jbegin = 0, latvalue_1d = lat_in / deg2rad)
+          IF (l3D) THEN
+            CALL xios_set_domain_attr("domain_climoz", ni_glo = nlon_in, ni = nlon_in, ibegin = 0, lonvalue_1d = lon_in / deg2rad)
+          ELSE
+            CALL xios_set_domain_attr("domain_climoz", ni_glo = 8, ni = 8, ibegin = 0, lonvalue_1d = (/ 0., 45., 90., 135., 180., 225., 270., 315. /))
+          ENDIF
+        ELSE
+          CALL xios_set_domain_attr("domain_climoz", nj_glo = nlat_in, nj = 0, jbegin = 0, latvalue_1d = null_array)
+          IF (l3D) THEN
+            CALL xios_set_domain_attr("domain_climoz", ni_glo = nlon_in, ni = 0, ibegin = 0, lonvalue_1d = null_array)
+          ELSE
+            CALL xios_set_domain_attr("domain_climoz", ni_glo = 8, ni = 0, ibegin = 0, lonvalue_1d = null_array)
+          ENDIF
+        ENDIF
+        CALL  xios_set_axis_attr("axis_climoz", n_glo = nlev_in)
+        CALL  xios_set_axis_attr("time_axis_climoz", n_glo = ntim_ou)
+        CALL  xios_set_axis_attr("time_axis_climoz", n_glo = ntim_ou)
+        CALL  xios_set_axis_attr("tr_climoz", n_glo = read_climoz)
+        CALL  xios_set_field_attr("tro3_out", enabled = .TRUE.)
+        CALL  xios_set_field_attr("tro3_out", enabled = .TRUE.)
+
+        IF (first) THEN
+          first = .FALSE.
+          RETURN
+        ENDIF
+      ENDIF
+
+      IF (is_mpi_root) THEN
+        !--- Longitudes management:
+        !    * Need to shift data if the origin of input file longitudes /= -pi
+        !    * Need to add some margin in longitude to ensure input interval contains
+        !      all the output intervals => at least one longitudes slice has to be
+        !      duplicated, possibly more for undersampling.
+        IF(l3D) THEN
+          IF (grid_type==unstructured) THEN
+            dx2 = 0
+          ELSE
+            !--- Compute input edges longitudes vector (no end point yet)
+            ALLOCATE(v1(nlon_in + 1))
+            v1(1) = (lon_in(nlon_in) + lon_in(1)) / 2. - pi
+            FORALL(i = 2:nlon_in) v1(i) = (lon_in(i - 1) + lon_in(i)) / 2.
+            v1(nlon_in + 1) = v1(1) + 2. * pi
+            DEALLOCATE(lon_in)
+
+            !--- Shift input longitudes vector until it contains first output point boundslon_reg(1,west)
+            v1 = v1 + 2 * pi * REAL(FLOOR((boundslon_reg(1, west) - v1(1)) / (2. * pi)))
+
+            !--- Ensure first input longitudes interval contains first output point boundslon_reg(1,west)
+            dx1 = locate(v1, boundslon_reg(1, west)) - 1
+            v1 = CSHIFT(v1, SHIFT = dx1, DIM = 1)
+            v1(nlon_in - dx1 + 2:) = v1(nlon_in - dx1 + 2:) + 2. * pi
+
+            !--- Extend input longitudes vector until last interval contains boundslon_reg(nlat_ou,east)
+            dx2 = 0; DO WHILE(v1(1 + dx2) + 2. * pi<=boundslon_reg(nlon_ou, east)); dx2 = dx2 + 1;
+            END DO
+
+            !--- Final edges longitudes vector (with margin and end point)
+            ALLOCATE(lon_in_edge(nlon_in + dx2 + 1)); lon_in_edge = [v1, v1(2:1 + dx2) + 2. * pi]
+            DEALLOCATE(v1)
+          ENDIF
+        END IF
+
+        !--- Compute sinus of intervals edges latitudes:
+        ALLOCATE(sinlat_in_edge(nlat_in + 1))
+        sinlat_in_edge(1) = -1. ; sinlat_in_edge(nlat_in + 1) = 1.
+        FORALL(j = 2:nlat_in) sinlat_in_edge(j) = SIN((lat_in(j - 1) + lat_in(j)) / 2.)
+        DEALLOCATE(lat_in)
+
+
+
+        !--- Check for contiguous years:
+        ib = 0; ie = 13
+        IF(nmth_in == 14) THEN; lprev = .FALSE.; lnext = .FALSE.
+        WRITE(lunout, *)'Using 14 months ozone climatology "climoz.nc"...'
+        ELSE
+          IF(lprev) WRITE(lunout, *)'Using "climoz_m.nc" last record (previous year).'
+          IF(.NOT.lprev) WRITE(lunout, *)"No previous year file ; assuming periodicity."
+          IF(lnext) WRITE(lunout, *)'Using "climoz_p.nc" first record (next year).'
+          IF(.NOT.lnext) WRITE(lunout, *)"No next year file ; assuming periodicity."
+          IF(.NOT.lprev) ib = 1
+          IF(.NOT.lnext) ie = 12
+        END IF
+        ALLOCATE(sta(ndims), cnt(ndims)); sta(:) = 1
+        IF(l3D) cnt = [nlon_in, nlat_in, nlev_in, 1]
+        IF(l2D) cnt = [        nlat_in, nlev_in, 1]
+        IF(l3D) ALLOCATE(o3_in3(nlon_in + dx2, nlat_in, nlev_in, ib:ie, read_climoz))
+        IF(l2D) ALLOCATE(o3_in2(nlat_in, nlev_in, ib:ie, read_climoz))
+
+        !--- Read full current file and one record each available contiguous file
+        DO iv = 1, read_climoz
+          CALL nf95_inq_varid(fID_in, vars_in(1), vID_in(iv))
+          IF(l3D) call NF95_GET_VAR(fID_in, vID_in(iv), o3_in3(1:nlon_in, :, :, 1:12, iv))
+          IF(l2D) call NF95_GET_VAR(fID_in, vID_in(iv), o3_in2(:, :, 1:12, iv))
+          IF(lprev) THEN; sta(ndims) = 12
+          CALL nf95_inq_varid(fID_in_m, vars_in(1), vID_in(iv))
+          IF(l3D) call NF95_GET_VAR(fID_in_m, vID_in(iv), o3_in3(1:nlon_in, :, :, 0, iv), sta, cnt)
+          IF(l2d) call NF95_GET_VAR(fID_in_m, vID_in(iv), o3_in2(:, :, 0, iv), sta, cnt)
+          END IF
+          IF(lnext) THEN; sta(ndims) = 1
+          CALL nf95_inq_varid(fID_in_p, vars_in(1), vID_in(iv))
+          IF(l3D) call NF95_GET_VAR(fID_in_p, vID_in(iv), o3_in3(1:nlon_in, :, :, 13, iv), sta, cnt)
+          IF(l2D) call NF95_GET_VAR(fID_in_p, vID_in(iv), o3_in2(:, :, 13, iv), sta, cnt)
+          END IF
+        END DO
+        IF(lprev.OR.lnext) DEALLOCATE(sta, cnt)
+        IF(lprev) CALL NF95_CLOSE(fID_in_m)
+        IF(lnext) CALL NF95_CLOSE(fID_in_p)
+
+        !--- Revert decreasing coordinates vector
+        IF(l3D) THEN
+          IF(ldec_lon) o3_in3(1:nlon_in, :, :, :, :) = o3_in3(nlon_in:1:-1, :, :, :, :)
+          IF(ldec_lat) o3_in3 = o3_in3(:, nlat_in:1:-1, :, :, :)
+          IF(ldec_lev) o3_in3 = o3_in3(:, :, nlev_in:1:-1, :, :)
+
+          IF (grid_type /= unstructured) THEN
+            !--- Shift values for longitude and duplicate some longitudes slices
+            o3_in3(1:nlon_in, :, :, :, :) = CSHIFT(o3_in3(1:nlon_in, :, :, :, :), SHIFT = dx1, DIM = 1)
+            o3_in3(nlon_in + 1:nlon_in + dx2, :, :, :, :) = o3_in3(1:dx2, :, :, :, :)
+          ENDIF
+        ELSE
+          IF(ldec_lat) o3_in2 = o3_in2(nlat_in:1:-1, :, :, :)
+          IF(ldec_lev) o3_in2 = o3_in2(:, nlev_in:1:-1, :, :)
+        END IF
+
+        !--- Deal with missing values
+        DO m = 1, read_climoz
+          WRITE(msg, '(a,i0)')"regr_lat_time_climoz: field Nr.", m
+          IF(nf90_get_att(fID_in, vID_in(m), "missing_value", NaN)/= nf90_noerr) THEN
+            IF(nf90_get_att(fID_in, vID_in(m), "_FillValue", NaN)/= nf90_noerr) THEN
+              WRITE(lunout, *)TRIM(msg) // ": no missing value attribute found."; CYCLE
+            END IF
+          END IF
+          WRITE(lunout, *)TRIM(msg) // ": missing value attribute found."
+          WRITE(lunout, *)"Trying to fill in NaNs ; a full field would be better."
+
+          !--- Check top layer contains no NaNs & search NaNs from top to ground
+          msg = TRIM(sub) // ": NaNs in top layer !"
+          IF(l3D) THEN
+            IF(ANY(o3_in3(:, :, 1, :, m)==NaN)) CALL abort_physic(sub, msg, 1)
+            DO k = 2, nlev_in
+              WHERE(o3_in3(:, :, k, :, m)==NaN) o3_in3(:, :, k, :, m) = o3_in3(:, :, k - 1, :, m)
+            END DO
+          ELSE
+            IF(ANY(o3_in2(:, 1, :, m)==NaN)) THEN
+              WRITE(lunout, *)msg
+              !--- Fill in latitudes where all values are missing
+              DO l = 1, nmth_in
+                !--- Next to south pole
+                j = 1;       DO WHILE(o3_in2(j, 1, l, m)==NaN); j = j + 1;
+                END DO
+                IF(j>1) &
+                        o3_in2(:j - 1, :, l, m) = SPREAD(o3_in2(j, :, l, m), DIM = 1, ncopies = j - 1)
+                !--- Next to north pole
+                j = nlat_in; DO WHILE(o3_in2(j, 1, l, m)==NaN); j = j + 1;
+                END DO
+                IF(j<nlat_in) &
+                        o3_in2(j + 1:, :, l, m) = SPREAD(o3_in2(j, :, l, m), DIM = 1, ncopies = nlat_in - j)
+              END DO
+            END IF
+
+            !--- Fill in high latitudes missing values
+            !--- Highest level been filled-in, so has always valid values.
+            DO k = 2, nlev_in
+              WHERE(o3_in2(:, k, :, m)==NaN) o3_in2(:, k, :, m) = o3_in2(:, k - 1, :, m)
+            END DO
+          END IF
+        END DO
+
+      ENDIF
+
+      !=============================================================================
+      IF(l3D) THEN                                                   !=== 3D FIELDS
+        !=============================================================================
+        IF (grid_type==unstructured) THEN
+          nlat_ou = klon_mpi
+
+          IF (is_mpi_root) THEN
+            ALLOCATE(o3_in3bis(nlon_in, nlat_in, nlev_in, 0:13, read_climoz))
+            o3_in3bis(:, :, :, ib:ie, :) = o3_in3(1:nlon_in, :, :, ib:ie, :)
+          ELSE
+            ALLOCATE(o3_in3bis(0, 0, 0, 0, read_climoz))
+          ENDIF
+          ALLOCATE(o3_regr_lonlat(1, nlat_ou, nlev_in, 0:13, read_climoz))
+
+          CALL xios_send_field("tro3_in", o3_in3bis(:, :, :, :, :))
+          CALL xios_recv_field("tro3_out", o3_regr_lonlat(1, :, :, :, :))
+        ELSE
+
+          !--- Regrid in longitude
+          ALLOCATE(o3_regr_lon(nlon_ou, nlat_in, nlev_in, ie - ib + 1, read_climoz))
+          CALL regr_conserv(1, o3_in3, xs = lon_in_edge, &
+                  xt = [boundslon_reg(1, west), boundslon_reg(:, east)], &
+                  vt = o3_regr_lon, slope = slopes(1, o3_in3, lon_in_edge))
+          DEALLOCATE(o3_in3)
+
+          !--- Regrid in latitude: averaging with respect to SIN(lat) is
+          !                        equivalent to weighting by COS(lat)
+          !--- (inverted indices in "o3_regr_lonlat" because "rlatu" is decreasing)
+          ALLOCATE(o3_regr_lonlat(nlon_ou, nlat_ou, nlev_in, 0:13, read_climoz))
+          CALL regr_conserv(2, o3_regr_lon, xs = sinlat_in_edge, &
+                  xt = [- 1., SIN(boundslat_reg(nlat_ou - 1:1:-1, south)), 1.], &
+                  vt = o3_regr_lonlat(:, nlat_ou:1:- 1, :, ib:ie, :), &
+                  slope = slopes(2, o3_regr_lon, sinlat_in_edge))
+          DEALLOCATE(o3_regr_lon)
+
+        ENDIF
+
+        !--- Duplicate previous/next record(s) if they are not available
+        IF(.NOT.lprev) o3_regr_lonlat(:, :, :, 0, :) = o3_regr_lonlat(:, :, :, 12, :)
+        IF(.NOT.lnext) o3_regr_lonlat(:, :, :, 13, :) = o3_regr_lonlat(:, :, :, 1, :)
+
+        !--- Regrid in time by linear interpolation:
+        ALLOCATE(o3_out3(nlon_ou, nlat_ou, nlev_in, ntim_ou, read_climoz))
+        IF(interpt) CALL regr_lint(4, o3_regr_lonlat, tmidmonth, tmidday, o3_out3)
+        IF(.NOT.interpt) o3_out3 = o3_regr_lonlat
+        DEALLOCATE(o3_regr_lonlat)
+
+        nlat_ou = nbp_lat
+        IF (grid_type==unstructured) THEN
+          CALL xios_send_field('o3_out', o3_out3)
+          ndims = 3
+          ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
+          CALL gather_mpi(o3_out3(1, :, :, :, :), o3_out3_glo)
+        ENDIF
+
+        !--- Create the output file and get the variable IDs:
+        CALL prepare_out(fID_in, nlev_in, ntim_ou, fID_ou, levID_ou, timID_ou, vID_ou, &
+                ndims, cal_ou)
+
+        IF (is_mpi_root) THEN
+          !--- Write remaining coordinate variables:
+          CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
+          IF(interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
+          IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
+
+          !--- Write to file (the order of "rlatu" is inverted in the output file):
+          IF (grid_type==unstructured) THEN
+
+            ALLOCATE(o3_out(nlat_ou, nlev_in, ntim_ou, read_climoz))
+            DO i = 1, klon_glo
+              o3_out(ind_cell_glo_glo(i), :, :, :) = o3_out3_glo(i, :, :, :)
+            ENDDO
+
+            DO m = 1, read_climoz
+              CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1, :, :, m))
+            END DO
+
+          ELSE
+            DO m = 1, read_climoz
+              CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out3(:, nlat_ou:1:-1, :, :, m))
+            END DO
+          ENDIF
+          CALL NF95_CLOSE(fID_ou)
+
+        ENDIF
+
+
+        !=============================================================================
+      ELSE                                                         !=== ZONAL FIELDS
+        !=============================================================================
+
+        IF (grid_type==unstructured) THEN
+          nlat_ou = klon_mpi
+
+          IF (is_mpi_root) THEN
+            ALLOCATE(o3_in2bis(8, nlat_in, nlev_in, 0:13, read_climoz))
+            o3_in2bis(:, :, :, ib:ie, :) = SPREAD(o3_in2, 1, 8)
+          ELSE
+            ALLOCATE(o3_in2bis(0, 0, 0, 0, read_climoz))
+          ENDIF
+          ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
+          CALL xios_send_field("tro3_in", o3_in2bis(:, :, :, :, :))
+          CALL xios_recv_field("tro3_out", o3_regr_lat(:, :, :, :))
+          IF(.NOT.lprev) o3_regr_lat(:, :, 0, :) = o3_regr_lat(:, :, 12, :)
+          IF(.NOT.lnext) o3_regr_lat(:, :, 13, :) = o3_regr_lat(:, :, 1, :)
+
+        ELSE
+          !--- Regrid in latitude: averaging with respect to SIN(lat) is
+          !                        equivalent to weighting by COS(lat)
+          !--- (inverted indices in "o3_regr_lat" because "rlatu" is decreasing)
+          ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
+          CALL regr_conserv(1, o3_in2, xs = sinlat_in_edge, &
+                  xt = [- 1., SIN(boundslat_reg(nlat_ou - 1:1:-1, south)), 1.], &
+                  vt = o3_regr_lat(nlat_ou:1:- 1, :, ib:ie, :), &
+                  slope = slopes(1, o3_in2, sinlat_in_edge))
+          DEALLOCATE(o3_in2)
+
+          !--- Duplicate previous/next record(s) if they are not available
+          IF(.NOT.lprev) o3_regr_lat(:, :, 0, :) = o3_regr_lat(:, :, 12, :)
+          IF(.NOT.lnext) o3_regr_lat(:, :, 13, :) = o3_regr_lat(:, :, 1, :)
+
+        ENDIF
+
+        !--- Regrid in time by linear interpolation:
+        ALLOCATE(o3_out2(nlat_ou, nlev_in, ntim_ou, read_climoz))
+        IF(interpt) CALL regr_lint(3, o3_regr_lat, tmidmonth, tmidday, o3_out2)
+        IF(.NOT.interpt) o3_out2 = o3_regr_lat
+        DEALLOCATE(o3_regr_lat)
+
+        nlat_ou = nbp_lat
+
+        IF (grid_type==unstructured) THEN
+          ndims = 3
+          ALLOCATE(o3_out2_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
+          CALL gather_mpi(o3_out2, o3_out2_glo)
+        ENDIF
+
+        !--- Create the output file and get the variable IDs:
+        CALL prepare_out(fID_in, nlev_in, ntim_ou, fID_ou, levID_ou, timID_ou, vID_ou, &
+                ndims, cal_ou)
+
+        IF (is_mpi_root) THEN
+
+          !--- Write remaining coordinate variables:
+          CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
+          IF(interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
+          IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
+
+          IF (grid_type==unstructured) THEN
+
+            ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
+            DO i = 1, klon_glo
+              o3_out(ind_cell_glo_glo(i), :, :, :) = o3_out2_glo(i, :, :, :)
+            ENDDO
+
+            DO m = 1, read_climoz
+              CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1, :, :, m))
+            END DO
+          ELSE
+            !--- Write to file (the order of "rlatu" is inverted in the output file):
+            DO m = 1, read_climoz
+              CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out2(nlat_ou:1:-1, :, :, m))
+            END DO
+          ENDIF
+
+          CALL NF95_CLOSE(fID_ou)
+
+        ENDIF
+
+        !=============================================================================
+      END IF
+      !=============================================================================
+
+      IF (is_mpi_root) CALL NF95_CLOSE(fID_in)
+
+    ENDIF ! is_omp_master
+
+    first = .FALSE.
+  END SUBROUTINE regr_horiz_time_climoz
+
+  !-------------------------------------------------------------------------------
+
+
+  !-------------------------------------------------------------------------------
+
+  SUBROUTINE prepare_out(fID_in, nlev_in, ntim_ou, fID_ou, vlevID, vtimID, &
+          vID_ou, ndims, cal_ou)
+    !-------------------------------------------------------------------------------
+    ! Purpose:  This subroutine creates the NetCDF output file, defines
+    !     dimensions and variables, and writes some of the coordinate variables.
+    !-------------------------------------------------------------------------------
+    USE regular_lonlat_mod, ONLY : lon_reg, lat_reg
+    USE regular_lonlat_mod, ONLY : lon_reg, lat_reg
+    USE mod_phys_lmdz_para, ONLY : is_mpi_root
+    USE mod_grid_phy_lmdz, ONLY : klon_glo
+
+    !-------------------------------------------------------------------------------
+    ! Arguments:
+    INTEGER, INTENT(IN) :: fID_in, nlev_in, ntim_ou
+    INTEGER, INTENT(OUT) :: fID_ou, vlevID, vtimID
+    INTEGER, INTENT(OUT) :: vID_ou(:)      ! dim(1/2) 1: O3day&night 2: O3daylight
+    INTEGER, INTENT(IN) :: ndims          ! fields rank (3 or 4)
+    CHARACTER(LEN = *), INTENT(IN) :: cal_ou ! calendar
+    !-------------------------------------------------------------------------------
+    ! Local variables:
+    INTEGER :: dlonID, dlatID, dlevID, dtimID, dIDs(4)
+    INTEGER :: vlonID, vlatID, ncerr, is
+    REAL, ALLOCATABLE :: latitude_glo_(:)
+    CHARACTER(LEN = 80) :: sub
+    INTEGER :: i
+
+
+    !-------------------------------------------------------------------------------
+
+    IF (is_mpi_root) THEN
+      sub = "prepare_out"
+      WRITE(lunout, *)"CALL sequence information: " // TRIM(sub)
+      CALL NF95_CREATE("climoz_LMDZ.nc", nf90_clobber, fID_ou)
+
+      !--- Dimensions:
+      IF(ndims==4) &
+              CALL nf95_def_dim(fID_ou, "rlonv", nlon_ou, dlonID)
+      CALL nf95_def_dim(fID_ou, "rlatu", nlat_ou, dlatID)
+      CALL nf95_def_dim(fID_ou, "plev", nlev_in, dlevID)
+      CALL nf95_def_dim(fID_ou, "time", ntim_ou, dtimID)
+
+      !--- Define coordinate variables:
+      IF(ndims==4) &
+              CALL nf95_def_var(fID_ou, "rlonv", nf90_float, dlonID, vlonID)
+      CALL nf95_def_var(fID_ou, "rlatu", nf90_float, dlatID, vlatID)
+      CALL nf95_def_var(fID_ou, "plev", nf90_float, dlevID, vlevID)
+      CALL nf95_def_var(fID_ou, "time", nf90_float, dtimID, vtimID)
+      IF(ndims==4) &
+              CALL NF95_PUT_ATT(fID_ou, vlonID, "units", "degrees_east")
+      CALL NF95_PUT_ATT(fID_ou, vlatID, "units", "degrees_north")
+      CALL NF95_PUT_ATT(fID_ou, vlevID, "units", "millibar")
+      CALL NF95_PUT_ATT(fID_ou, vtimID, "units", "days since 2000-1-1")
+      IF(ndims==4) &
+              CALL NF95_PUT_ATT(fID_ou, vlonID, "standard_name", "longitude")
+      CALL NF95_PUT_ATT(fID_ou, vlatID, "standard_name", "latitude")
+      CALL NF95_PUT_ATT(fID_ou, vlevID, "standard_name", "air_pressure")
+      CALL NF95_PUT_ATT(fID_ou, vtimID, "standard_name", "time")
+      CALL NF95_PUT_ATT(fID_ou, vlevID, "long_name", "air pressure")
+      CALL NF95_PUT_ATT(fID_ou, vtimID, "calendar", cal_ou)
+
+      !--- Define the main variables:
+      IF(ndims==3) dIDs(1:3) = [ dlatID, dlevID, dtimID]
+      IF(ndims==4) dIDs = [dlonID, dlatID, dlevID, dtimID]
+      CALL nf95_def_var(fID_ou, vars_in(1), nf90_float, dIDs(1:ndims), vID_ou(1))
+      CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "long_name", "ozone mole fraction")
+      CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "standard_name", "mole_fraction_of_ozone&
+              _in_air")
+      IF(SIZE(vID_ou) == 2) THEN
+        CALL nf95_def_var(fID_ou, vars_in(2), nf90_float, dIDs(1:ndims), vID_ou(2))
+        CALL NF95_PUT_ATT(fID_ou, vID_ou(2), "long_name", "ozone mole fraction in da&
+                ylight")
+      END IF
+
+      !--- Global attributes:
+      ! The following commands, copying attributes, may fail. That is OK.
+      ! It should just mean that the attribute is not defined in the input file.
+      CALL NF95_COPY_ATT(fID_in, nf90_global, "Conventions", fID_ou, nf90_global, ncerr)
+      CALL handle_err_copy_att("Conventions")
+      CALL NF95_COPY_ATT(fID_in, nf90_global, "title", fID_ou, nf90_global, ncerr)
+      CALL handle_err_copy_att("title")
+      CALL NF95_COPY_ATT(fID_in, nf90_global, "institution", fID_ou, nf90_global, ncerr)
+      CALL handle_err_copy_att("institution")
+      CALL NF95_COPY_ATT(fID_in, nf90_global, "source", fID_ou, nf90_global, ncerr)
+      CALL handle_err_copy_att("source")
+      CALL NF95_PUT_ATT (fID_ou, nf90_global, "comment", "Regridded for LMDZ")
+      CALL NF95_ENDDEF(fID_ou)
+
+      IF (grid_type==unstructured) THEN
+        ALLOCATE(latitude_glo_(klon_glo))
+        DO i = 1, klon_glo
+          latitude_glo_(ind_cell_glo_glo(i)) = latitude_glo(i)
+        ENDDO
+        CALL NF95_PUT_VAR(fID_ou, vlatID, latitude_glo_)
       ELSE
-        ntim_ou=14
-        cal_ou=cal_in
-      END IF
-    ENDIF
-
-    IF (grid_type==unstructured) THEN
-      CALL bcast_mpi(nlon_in)
-      CALL bcast_mpi(nlat_in)
-      CALL bcast_mpi(nlev_in)
-      CALL bcast_mpi(l3d)
-      CALL bcast_mpi(tmidmonth)
-      IF(interpt) CALL bcast_mpi(tmidday)
-      CALL bcast_mpi(ntim_ou)
-
-      IF (is_mpi_root) THEN
-        CALL xios_set_domain_attr("domain_climoz",nj_glo=nlat_in, nj=nlat_in, jbegin=0, latvalue_1d=lat_in/deg2rad)
-        IF (l3D) THEN
-          CALL xios_set_domain_attr("domain_climoz",ni_glo=nlon_in, ni=nlon_in, ibegin=0, lonvalue_1d=lon_in/deg2rad)
-        ELSE
-          CALL xios_set_domain_attr("domain_climoz",ni_glo=8, ni=8, ibegin=0, lonvalue_1d = (/ 0.,45.,90.,135.,180.,225.,270., 315. /))
-        ENDIF
-      ELSE
-        CALL xios_set_domain_attr("domain_climoz",nj_glo=nlat_in, nj=0, jbegin=0, latvalue_1d=null_array )
-        IF (l3D) THEN
-          CALL xios_set_domain_attr("domain_climoz",ni_glo=nlon_in, ni=0, ibegin=0, lonvalue_1d=null_array)
-        ELSE
-          CALL xios_set_domain_attr("domain_climoz",ni_glo=8, ni=0, ibegin=0, lonvalue_1d=null_array)
-        ENDIF
-      ENDIF
-      CALL  xios_set_axis_attr("axis_climoz", n_glo=nlev_in)
-      CALL  xios_set_axis_attr("time_axis_climoz", n_glo=ntim_ou)
-      CALL  xios_set_axis_attr("time_axis_climoz", n_glo=ntim_ou)
-      CALL  xios_set_axis_attr("tr_climoz", n_glo=read_climoz)
-      CALL  xios_set_field_attr("tro3_out", enabled=.TRUE.)
-      CALL  xios_set_field_attr("tro3_out", enabled=.TRUE.)
-      
-      IF (first) THEN
-        first=.FALSE.
-        RETURN
+        !--- Write one of the coordinate variables:
+        IF(ndims==4) CALL NF95_PUT_VAR(fID_ou, vlonID, lon_reg / deg2rad)
+        CALL NF95_PUT_VAR(fID_ou, vlatID, lat_reg(nlat_ou:1:-1) / deg2rad)
+        !    (convert from rad to degrees and sort in ascending order)
       ENDIF
     ENDIF
-    
-    
-    IF (is_mpi_root) THEN     
-      !--- Longitudes management:
-      !    * Need to shift data if the origin of input file longitudes /= -pi
-      !    * Need to add some margin in longitude to ensure input interval contains
-      !      all the output intervals => at least one longitudes slice has to be
-      !      duplicated, possibly more for undersampling.
-      IF(l3D) THEN
-        IF (grid_type==unstructured) THEN
-          dx2=0
-        ELSE
-          !--- Compute input edges longitudes vector (no end point yet)
-          ALLOCATE(v1(nlon_in+1))
-          v1(1)=(lon_in(nlon_in)+lon_in(1))/2.-pi
-          FORALL(i=2:nlon_in) v1(i)=(lon_in(i-1)+lon_in(i))/2.
-          v1(nlon_in+1)=v1(1)+2.*pi
-          DEALLOCATE(lon_in)
-
-          !--- Shift input longitudes vector until it contains first output point boundslon_reg(1,west)
-          v1=v1+2*pi*REAL(FLOOR((boundslon_reg(1,west)-v1(1))/(2.*pi)))
-
-          !--- Ensure first input longitudes interval contains first output point boundslon_reg(1,west)
-          dx1=locate(v1,boundslon_reg(1,west))-1
-          v1=CSHIFT(v1,SHIFT=dx1,DIM=1) 
-          v1(nlon_in-dx1+2:)=v1(nlon_in-dx1+2:)+2.*pi
-    
-          !--- Extend input longitudes vector until last interval contains boundslon_reg(nlat_ou,east)
-          dx2=0; DO WHILE(v1(1+dx2)+2.*pi<=boundslon_reg(nlon_ou,east)); dx2=dx2+1; END DO
-
-          !--- Final edges longitudes vector (with margin and end point)
-          ALLOCATE(lon_in_edge(nlon_in+dx2+1)); lon_in_edge=[v1,v1(2:1+dx2)+2.*pi]
-          DEALLOCATE(v1)
-        ENDIF
-      END IF
-
-      !--- Compute sinus of intervals edges latitudes:
-      ALLOCATE(sinlat_in_edge(nlat_in+1))
-      sinlat_in_edge(1) = -1. ; sinlat_in_edge(nlat_in+1) = 1.
-      FORALL(j=2:nlat_in) sinlat_in_edge(j)=SIN((lat_in(j-1)+lat_in(j))/2.)
-      DEALLOCATE(lat_in)
-
-
-
-      !--- Check for contiguous years:
-      ib=0; ie=13
-      IF(nmth_in == 14) THEN; lprev=.FALSE.; lnext=.FALSE.
-        WRITE(lunout,*)'Using 14 months ozone climatology "climoz.nc"...'
-      ELSE  
-        IF(     lprev) WRITE(lunout,*)'Using "climoz_m.nc" last record (previous year).'
-        IF(.NOT.lprev) WRITE(lunout,*)"No previous year file ; assuming periodicity."
-        IF(     lnext) WRITE(lunout,*)'Using "climoz_p.nc" first record (next year).'
-        IF(.NOT.lnext) WRITE(lunout,*)"No next year file ; assuming periodicity."
-        IF(.NOT.lprev) ib=1
-        IF(.NOT.lnext) ie=12
-      END IF
-      ALLOCATE(sta(ndims),cnt(ndims)); sta(:)=1  
-      IF(l3D) cnt=[nlon_in,nlat_in,nlev_in,1]
-      IF(l2D) cnt=[        nlat_in,nlev_in,1]  
-      IF(l3D) ALLOCATE(o3_in3(nlon_in+dx2,nlat_in,nlev_in,ib:ie,read_climoz))
-      IF(l2D) ALLOCATE(o3_in2(            nlat_in,nlev_in,ib:ie,read_climoz))
-
-      !--- Read full current file and one record each available contiguous file
-      DO iv=1,read_climoz
-        CALL NF95_INQ_VARID(fID_in, vars_in(1), vID_in(iv))
-        IF(l3D) call NF95_GET_VAR(fID_in, vID_in(iv), o3_in3(1:nlon_in,:,:,1:12,iv))
-        IF(l2D) call NF95_GET_VAR(fID_in, vID_in(iv), o3_in2(          :,:,1:12,iv))
-        IF(lprev) THEN; sta(ndims)=12  
-          CALL NF95_INQ_VARID(fID_in_m, vars_in(1), vID_in(iv))
-          IF(l3D) call NF95_GET_VAR(fID_in_m,vID_in(iv),o3_in3(1:nlon_in,:,:, 0,iv),sta,cnt)
-          IF(l2d) call NF95_GET_VAR(fID_in_m,vID_in(iv),o3_in2(          :,:, 0,iv),sta,cnt)
-        END IF
-        IF(lnext) THEN; sta(ndims)=1  
-          CALL NF95_INQ_VARID(fID_in_p, vars_in(1), vID_in(iv))
-          IF(l3D) call NF95_GET_VAR(fID_in_p,vID_in(iv),o3_in3(1:nlon_in,:,:,13,iv),sta,cnt)
-          IF(l2D) call NF95_GET_VAR(fID_in_p,vID_in(iv),o3_in2(          :,:,13,iv),sta,cnt)
-        END IF
-      END DO
-      IF(lprev.OR.lnext) DEALLOCATE(sta,cnt)
-      IF(lprev) CALL NF95_CLOSE(fID_in_m)
-      IF(lnext) CALL NF95_CLOSE(fID_in_p)
-
-      !--- Revert decreasing coordinates vector
-      IF(l3D) THEN
-        IF(ldec_lon) o3_in3(1:nlon_in,:,:,:,:) = o3_in3(nlon_in:1:-1,:,:,:,:)
-        IF(ldec_lat) o3_in3 = o3_in3(:,nlat_in:1:-1,:,:,:)
-        IF(ldec_lev) o3_in3 = o3_in3(:,:,nlev_in:1:-1,:,:)
-        
-        IF (grid_type /= unstructured) THEN
-          !--- Shift values for longitude and duplicate some longitudes slices
-          o3_in3(1:nlon_in,:,:,:,:)=CSHIFT(o3_in3(1:nlon_in,:,:,:,:),SHIFT=dx1,DIM=1)
-          o3_in3(nlon_in+1:nlon_in+dx2,:,:,:,:)=o3_in3(1:dx2,:,:,:,:)
-        ENDIF
-      ELSE
-        IF(ldec_lat) o3_in2 = o3_in2(  nlat_in:1:-1,:,:,:)
-        IF(ldec_lev) o3_in2 = o3_in2(  :,nlev_in:1:-1,:,:)
-      END IF
-
-     !--- Deal with missing values 
-      DO m=1, read_climoz
-        WRITE(msg,'(a,i0)')"regr_lat_time_climoz: field Nr.",m
-        IF(NF90_GET_ATT(fID_in,vID_in(m),"missing_value",NaN)/= nf90_noerr) THEN
-          IF(NF90_GET_ATT(fID_in, vID_in(m),"_FillValue",NaN)/= nf90_noerr) THEN
-            WRITE(lunout,*)TRIM(msg)//": no missing value attribute found."; CYCLE
-          END IF
-        END IF
-        WRITE(lunout,*)TRIM(msg)//": missing value attribute found."
-        WRITE(lunout,*)"Trying to fill in NaNs ; a full field would be better."
-
-        !--- Check top layer contains no NaNs & search NaNs from top to ground
-        msg=TRIM(sub)//": NaNs in top layer !"
-        IF(l3D) THEN
-          IF(ANY(o3_in3(:,:,1,:,m)==NaN)) CALL abort_physic(sub,msg,1)
-          DO k = 2,nlev_in
-            WHERE(o3_in3(:,:,k,:,m)==NaN) o3_in3(:,:,k,:,m)=o3_in3(:,:,k-1,:,m)
-          END DO
-        ELSE
-          IF(ANY(o3_in2(  :,1,:,m)==NaN)) THEN
-            WRITE(lunout,*)msg
-            !--- Fill in latitudes where all values are missing
-            DO l=1,nmth_in
-              !--- Next to south pole
-              j=1;       DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
-              IF(j>1) &
-                o3_in2(:j-1,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=j-1)
-              !--- Next to north pole
-              j=nlat_in; DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
-              IF(j<nlat_in) &
-                o3_in2(j+1:,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=nlat_in-j)
-            END DO
-          END IF
-
-          !--- Fill in high latitudes missing values
-          !--- Highest level been filled-in, so has always valid values.
-          DO k = 2,nlev_in
-            WHERE(o3_in2(:,k,:,m)==NaN) o3_in2(:,k,:,m)=o3_in2(:,k-1,:,m)
-          END DO
-        END IF
-      END DO
-
-    ENDIF
-    
-    !=============================================================================
-    IF(l3D) THEN                                                   !=== 3D FIELDS
-    !=============================================================================
-     IF (grid_type==unstructured) THEN
-       nlat_ou=klon_mpi
-       
-       IF (is_mpi_root) THEN 
-         ALLOCATE(o3_in3bis(nlon_in,nlat_in,nlev_in,0:13,read_climoz)) 
-         o3_in3bis(:,:,:,ib:ie,:)=o3_in3(1:nlon_in,:,:,ib:ie,:)
-       ELSE
-         ALLOCATE(o3_in3bis(0,0,0,0,read_climoz))
-       ENDIF
-       ALLOCATE(o3_regr_lonlat(1, nlat_ou, nlev_in, 0:13, read_climoz))
-       
-       CALL xios_send_field("tro3_in",o3_in3bis(:,:,:,:,:))
-       CALL xios_recv_field("tro3_out",o3_regr_lonlat(1,:,:,:,:))
-     ELSE
-          
-       !--- Regrid in longitude
-        ALLOCATE(o3_regr_lon(nlon_ou, nlat_in, nlev_in, ie-ib+1, read_climoz))
-        CALL regr_conserv(1, o3_in3, xs = lon_in_edge,                             &
-                            xt = [boundslon_reg(1,west),boundslon_reg(:,east)],    &
-                            vt = o3_regr_lon, slope = slopes(1,o3_in3, lon_in_edge))
-        DEALLOCATE(o3_in3)
-
-        !--- Regrid in latitude: averaging with respect to SIN(lat) is
-        !                        equivalent to weighting by COS(lat)
-        !--- (inverted indices in "o3_regr_lonlat" because "rlatu" is decreasing)
-        ALLOCATE(o3_regr_lonlat(nlon_ou, nlat_ou, nlev_in, 0:13, read_climoz))
-        CALL regr_conserv(2, o3_regr_lon, xs = sinlat_in_edge,                     &
-                        xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
-                        vt = o3_regr_lonlat(:,nlat_ou:1:- 1,:,ib:ie,:),            &
-                   slope = slopes(2,o3_regr_lon, sinlat_in_edge))
-        DEALLOCATE(o3_regr_lon)
-
-     ENDIF
-
-     !--- Duplicate previous/next record(s) if they are not available
-     IF(.NOT.lprev) o3_regr_lonlat(:,:,:, 0,:) = o3_regr_lonlat(:,:,:,12,:)
-     IF(.NOT.lnext) o3_regr_lonlat(:,:,:,13,:) = o3_regr_lonlat(:,:,:, 1,:)
-      
-     !--- Regrid in time by linear interpolation:
-     ALLOCATE(o3_out3(nlon_ou, nlat_ou, nlev_in, ntim_ou, read_climoz))
-     IF(     interpt) CALL regr_lint(4,o3_regr_lonlat,tmidmonth,tmidday,o3_out3)
-     IF(.NOT.interpt) o3_out3=o3_regr_lonlat
-     DEALLOCATE(o3_regr_lonlat)
-
-     nlat_ou=nbp_lat
-     IF (grid_type==unstructured) THEN
-       CALL xios_send_field('o3_out',o3_out3)
-       ndims=3 
-       ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
-       CALL gather_mpi(o3_out3(1,:,:,:,:), o3_out3_glo)
-     ENDIF
-
-    !--- Create the output file and get the variable IDs:
-    CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
-                     ndims, cal_ou)
-
-    IF (is_mpi_root) THEN
-      !--- Write remaining coordinate variables:
-      CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
-      IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
-      IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
-
-      !--- Write to file (the order of "rlatu" is inverted in the output file):
-        IF (grid_type==unstructured) THEN
-
-          ALLOCATE(o3_out(nlat_ou, nlev_in, ntim_ou, read_climoz))
-          DO i=1,klon_glo
-            o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out3_glo(i,:,:,:)
-          ENDDO
-
-          DO m = 1, read_climoz
-            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
-          END DO
-          
-        ELSE
-          DO m = 1, read_climoz
-            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out3(:,nlat_ou:1:-1,:,:,m))
-          END DO
-      ENDIF
-      CALL NF95_CLOSE(fID_ou)
-
-
-    ENDIF
-
-
-    !=============================================================================
-    ELSE                                                         !=== ZONAL FIELDS
-    !=============================================================================
-   
-     IF (grid_type==unstructured) THEN
-       nlat_ou=klon_mpi
-
-       IF (is_mpi_root) THEN 
-         ALLOCATE(o3_in2bis(8,nlat_in,nlev_in,0:13,read_climoz)) 
-         o3_in2bis(:,:,:,ib:ie,:)=SPREAD(o3_in2,1,8)
-       ELSE
-         ALLOCATE(o3_in2bis(0,0,0,0,read_climoz))
-       ENDIF
-       ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
-       CALL xios_send_field("tro3_in",o3_in2bis(:,:,:,:,:))
-       CALL xios_recv_field("tro3_out",o3_regr_lat(:,:,:,:))
-       IF(.NOT.lprev) o3_regr_lat(:,:, 0, :) = o3_regr_lat(:,:,12,:)
-       IF(.NOT.lnext) o3_regr_lat(:,:,13, :) = o3_regr_lat(:,:, 1,:)
-     
-     ELSE 
-        !--- Regrid in latitude: averaging with respect to SIN(lat) is
-        !                        equivalent to weighting by COS(lat)
-        !--- (inverted indices in "o3_regr_lat" because "rlatu" is decreasing)
-        ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
-        CALL regr_conserv(1, o3_in2, xs = sinlat_in_edge,                          &
-                        xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
-                        vt = o3_regr_lat(nlat_ou:1:- 1,:,ib:ie,:),                 &
-                     slope = slopes(1,o3_in2, sinlat_in_edge))
-        DEALLOCATE(o3_in2)
-
-        !--- Duplicate previous/next record(s) if they are not available
-        IF(.NOT.lprev) o3_regr_lat(:,:, 0,:) = o3_regr_lat(:,:,12,:)
-        IF(.NOT.lnext) o3_regr_lat(:,:,13,:) = o3_regr_lat(:,:, 1,:)
-
-     ENDIF
-     
-      !--- Regrid in time by linear interpolation:
-      ALLOCATE(o3_out2(nlat_ou, nlev_in, ntim_ou, read_climoz))
-      IF(     interpt) CALL regr_lint(3,o3_regr_lat, tmidmonth, tmidday, o3_out2)
-      IF(.NOT.interpt) o3_out2=o3_regr_lat
-      DEALLOCATE(o3_regr_lat)
-
-      nlat_ou=nbp_lat
-    
-      IF (grid_type==unstructured) THEN
-        ndims=3 
-        ALLOCATE(o3_out2_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
-        CALL gather_mpi(o3_out2, o3_out2_glo)
-      ENDIF
-      
-      !--- Create the output file and get the variable IDs:
-      CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
-                         ndims, cal_ou)
-
-      IF (is_mpi_root) THEN
-      
-        !--- Write remaining coordinate variables:
-        CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
-        IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
-        IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
-
-        IF (grid_type==unstructured) THEN
-
-          ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
-          DO i=1,klon_glo
-            o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out2_glo(i,:,:,:)
-          ENDDO
-
-
-          DO m = 1, read_climoz
-            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
-          END DO
-        ELSE
-          !--- Write to file (the order of "rlatu" is inverted in the output file):
-          DO m = 1, read_climoz
-            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out2(nlat_ou:1:-1,:,:,m))
-          END DO
-        ENDIF
-        
-        CALL NF95_CLOSE(fID_ou)
-      
-      ENDIF
-
-    !=============================================================================
-    END IF
-    !=============================================================================
-
-    IF (is_mpi_root) CALL NF95_CLOSE(fID_in)
-
-  ENDIF ! is_omp_master
-
-  first=.FALSE.
-END SUBROUTINE regr_horiz_time_climoz
+
+  CONTAINS
+
+    !-------------------------------------------------------------------------------
+
+    SUBROUTINE handle_err_copy_att(att_name)
+
+      !-------------------------------------------------------------------------------
+      USE netcdf, ONLY : nf90_noerr, NF90_strerror
+      !-------------------------------------------------------------------------------
+      ! Arguments:
+      CHARACTER(LEN = *), INTENT(IN) :: att_name
+      !-------------------------------------------------------------------------------
+      IF(ncerr /= nf90_noerr) &
+              WRITE(lunout, *)TRIM(sub) // " prepare_out NF95_COPY_ATT " // TRIM(att_name) // &
+                      " -- " // TRIM(NF90_strerror(ncerr))
+
+    END SUBROUTINE handle_err_copy_att
+
+    !-------------------------------------------------------------------------------
+
+  END SUBROUTINE prepare_out
+
+  !-------------------------------------------------------------------------------
+
+END MODULE regr_horiz_time_climoz_m
 
 !-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-
-SUBROUTINE prepare_out(fID_in, nlev_in, ntim_ou, fID_ou, vlevID, vtimID, &
-                       vID_ou, ndims, cal_ou)
-!-------------------------------------------------------------------------------
-! Purpose:  This subroutine creates the NetCDF output file, defines
-!     dimensions and variables, and writes some of the coordinate variables.
-!-------------------------------------------------------------------------------
-  USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
-  USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
-  USE mod_phys_lmdz_para, ONLY: is_mpi_root
-  USE mod_grid_phy_lmdz, ONLY: klon_glo
-
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER, INTENT(IN)  :: fID_in, nlev_in, ntim_ou
-  INTEGER, INTENT(OUT) :: fID_ou, vlevID,  vtimID
-  INTEGER, INTENT(OUT) :: vID_ou(:)      ! dim(1/2) 1: O3day&night 2: O3daylight
-  INTEGER, INTENT(IN)  :: ndims          ! fields rank (3 or 4)
-  CHARACTER(LEN=*), INTENT(IN) :: cal_ou ! calendar
-!-------------------------------------------------------------------------------
-! Local variables:
-  INTEGER :: dlonID, dlatID, dlevID, dtimID, dIDs(4)
-  INTEGER :: vlonID, vlatID, ncerr,  is
-  REAL,ALLOCATABLE    :: latitude_glo_(:)
-  CHARACTER(LEN=80) :: sub
-  INTEGER :: i
-
-
-!-------------------------------------------------------------------------------
-  
-  IF (is_mpi_root) THEN  
-    sub="prepare_out"
-    WRITE(lunout,*)"CALL sequence information: "//TRIM(sub)
-    CALL NF95_CREATE("climoz_LMDZ.nc", NF90_clobber, fID_ou)
-
-  !--- Dimensions:
-    IF(ndims==4) &
-    CALL NF95_DEF_DIM(fID_ou, "rlonv", nlon_ou, dlonID)
-    CALL NF95_DEF_DIM(fID_ou, "rlatu", nlat_ou, dlatID)
-    CALL NF95_DEF_DIM(fID_ou, "plev",  nlev_in, dlevID)
-    CALL NF95_DEF_DIM(fID_ou, "time",  ntim_ou, dtimID)
-
-    !--- Define coordinate variables:
-    IF(ndims==4) &
-    CALL NF95_DEF_VAR(fID_ou, "rlonv", NF90_FLOAT, dlonID, vlonID)
-    CALL NF95_DEF_VAR(fID_ou, "rlatu", NF90_FLOAT, dlatID, vlatID)
-    CALL NF95_DEF_VAR(fID_ou, "plev",  NF90_FLOAT, dlevID, vlevID)
-    CALL NF95_DEF_VAR(fID_ou, "time",  NF90_FLOAT, dtimID, vtimID)
-    IF(ndims==4) &
-    CALL NF95_PUT_ATT(fID_ou, vlonID, "units", "degrees_east")
-    CALL NF95_PUT_ATT(fID_ou, vlatID, "units", "degrees_north")
-    CALL NF95_PUT_ATT(fID_ou, vlevID, "units", "millibar")
-    CALL NF95_PUT_ATT(fID_ou, vtimID, "units", "days since 2000-1-1")
-    IF(ndims==4) &
-    CALL NF95_PUT_ATT(fID_ou, vlonID, "standard_name", "longitude")
-    CALL NF95_PUT_ATT(fID_ou, vlatID, "standard_name", "latitude")
-    CALL NF95_PUT_ATT(fID_ou, vlevID, "standard_name", "air_pressure")
-    CALL NF95_PUT_ATT(fID_ou, vtimID, "standard_name", "time")
-    CALL NF95_PUT_ATT(fID_ou, vlevID, "long_name",     "air pressure")
-    CALL NF95_PUT_ATT(fID_ou, vtimID, "calendar",      cal_ou)
-
-  !--- Define the main variables:
-    IF(ndims==3) dIDs(1:3) = [ dlatID, dlevID, dtimID]
-    IF(ndims==4) dIDs=[dlonID, dlatID, dlevID, dtimID]
-    CALL NF95_DEF_VAR(fID_ou, vars_in(1), NF90_FLOAT, dIDs(1:ndims), vID_ou(1))
-    CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "long_name", "ozone mole fraction")
-    CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "standard_name", "mole_fraction_of_ozone&
-  _in_air")
-    IF(SIZE(vID_ou) == 2) THEN
-      CALL NF95_DEF_VAR(fID_ou, vars_in(2), NF90_FLOAT, dIDs(1:ndims), vID_ou(2))
-      CALL NF95_PUT_ATT(fID_ou, vID_ou(2), "long_name","ozone mole fraction in da&
-  ylight")
-    END IF
-
-  !--- Global attributes:
-  ! The following commands, copying attributes, may fail. That is OK.
-  ! It should just mean that the attribute is not defined in the input file.
-    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"Conventions",fID_ou,NF90_GLOBAL, ncerr)
-    CALL handle_err_copy_att("Conventions")
-    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"title",      fID_ou,NF90_GLOBAL, ncerr)
-    CALL handle_err_copy_att("title")
-    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"institution",fID_ou,NF90_GLOBAL, ncerr)
-    CALL handle_err_copy_att("institution")
-    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"source",     fID_ou,NF90_GLOBAL, ncerr)
-    CALL handle_err_copy_att("source")
-    CALL NF95_PUT_ATT (fID_ou,NF90_GLOBAL,"comment", "Regridded for LMDZ")
-    CALL NF95_ENDDEF(fID_ou)
-
-    IF (grid_type==unstructured) THEN
-      ALLOCATE(latitude_glo_(klon_glo))
-      DO i=1,klon_glo
-        latitude_glo_(ind_cell_glo_glo(i))=latitude_glo(i)
-      ENDDO
-      CALL NF95_PUT_VAR(fID_ou, vlatID, latitude_glo_)
-    ELSE
-      !--- Write one of the coordinate variables:
-      IF(ndims==4) CALL NF95_PUT_VAR(fID_ou, vlonID, lon_reg/deg2rad)
-      CALL NF95_PUT_VAR(fID_ou, vlatID, lat_reg(nlat_ou:1:-1)/deg2rad)
-    !    (convert from rad to degrees and sort in ascending order)
-    ENDIF
-  ENDIF
-  
-CONTAINS
-
-!-------------------------------------------------------------------------------
-
-SUBROUTINE handle_err_copy_att(att_name)
-
-!-------------------------------------------------------------------------------
-  USE netcdf, ONLY: nf90_noerr, NF90_strerror
-!-------------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*), INTENT(IN) :: att_name
-!-------------------------------------------------------------------------------
-  IF(ncerr /= nf90_noerr) &
-    WRITE(lunout,*)TRIM(sub)//" prepare_out NF95_COPY_ATT "//TRIM(att_name)//  &
-                      " -- "//TRIM(NF90_strerror(ncerr))
-
-END SUBROUTINE handle_err_copy_att
-
-!-------------------------------------------------------------------------------
-
-END SUBROUTINE prepare_out
-
-!-------------------------------------------------------------------------------
-
-END MODULE regr_horiz_time_climoz_m
-
-!-------------------------------------------------------------------------------
Index: /LMDZ6/branches/Amaury_dev/libf/phylmd/regr_pr_time_av_m.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmd/regr_pr_time_av_m.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmd/regr_pr_time_av_m.F90	(revision 5100)
@@ -113,6 +113,6 @@
 !-------------------------------------------------------------------------------
   USE dimphy,         ONLY: klon
-  USE netcdf95,       ONLY: NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, &
-                            NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION, nf95_get_var
+  USE netcdf95,       ONLY: nf95_inq_varid, NF95_INQUIRE_VARIABLE, &
+                            nf95_inq_dimid, nf95_inquire_dimension, nf95_get_var
   USE netcdf,         ONLY: NF90_INQ_VARID, nf90_noerr
   USE assert_m,       ONLY: assert
@@ -213,6 +213,6 @@
       lPrTfile=lAdjTro.AND.NF90_INQ_VARID(fID,"tropopause_air_pressure",vID)==nf90_noerr
       lO3Tfile=lAdjTro.AND.NF90_INQ_VARID(fID,"tro3_at_tropopause"     ,vID)==nf90_noerr
-      CALL NF95_INQ_DIMID(fID,"time",vID)
-      CALL NF95_INQUIRE_DIMENSION(fID,vID,nclen=ntim_in)
+      CALL nf95_inq_dimid(fID,"time",vID)
+      CALL nf95_inquire_dimension(fID,vID,nclen=ntim_in)
       linterp=PRESENT(time_in).AND.ntim_in==14
       ALLOCATE(v1(nlon,nlat,nlev_in,n_var))
@@ -480,5 +480,5 @@
   CHARACTER(LEN=*), INTENT(IN)    :: var
 !-------------------------------------------------------------------------------
-  CALL NF95_INQ_VARID(fID, TRIM(var), vID)
+  CALL nf95_inq_varid(fID, TRIM(var), vID)
   CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
   IF(n_dim==2) call NF95_GET_VAR(fID,vID,v(1,:), start=[  1,irec])
@@ -511,5 +511,5 @@
 !-------------------------------------------------------------------------------
   DO i=1,SIZE(nam)
-    CALL NF95_INQ_VARID(fID, TRIM(nam(i)), vID)
+    CALL nf95_inq_varid(fID, TRIM(nam(i)), vID)
     CALL NF95_INQUIRE_VARIABLE(fID, vID, ndims=n_dim)
     IF(n_dim==3) call NF95_GET_VAR(fID,vID,v(1,:,:,i), start=[  1,1,irec])
Index: /LMDZ6/branches/Amaury_dev/libf/phylmdiso/limit_read_mod.F90
===================================================================
--- /LMDZ6/branches/Amaury_dev/libf/phylmdiso/limit_read_mod.F90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/libf/phylmdiso/limit_read_mod.F90	(revision 5100)
@@ -352,5 +352,5 @@
 !$OMP MASTER  ! Only master thread
        IF (is_mpi_root) THEN ! Only master processus
-          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
+          ierr = nf90_open ('limit.nc', nf90_nowrite, nid)
           IF (ierr /= nf90_noerr) CALL abort_physic(modname,&
                'Pb d''ouverture du fichier de conditions aux limites',1)
@@ -358,5 +358,5 @@
           !--- WARNING IF CALENDAR IS KNOWN AND DOES NOT MATCH THE ONE OF LMDZ
           ierr=NF90_INQ_VARID(nid, 'TEMPS', nvarid)
-          ierr=NF90_GET_ATT(nid, nvarid, 'calendar', calendar)
+          ierr=nf90_get_att(nid, nvarid, 'calendar', calendar)
           IF(ierr==nf90_noerr.AND.calendar/=calend.AND.prt_level>=1) THEN
              WRITE(lunout,*)'BEWARE: gcm and limit.nc calendars differ: '
@@ -387,5 +387,5 @@
           IF(nn/=klon_glo) CALL abort_physic(modname,abort_message,1)
 
-          ierr = NF90_CLOSE(nid)
+          ierr = nf90_close(nid)
           IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Pb when closing file', 1)
        END IF ! is_mpi_root
@@ -450,5 +450,5 @@
        IF (is_mpi_root) THEN ! Only master processus!
 
-          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
+          ierr = nf90_open ('limit.nc', nf90_nowrite, nid)
           IF (ierr /= nf90_noerr) CALL abort_physic(modname,&
                'Pb d''ouverture du fichier de conditions aux limites',1)
@@ -564,5 +564,5 @@
 
 !****************************************************************************************
-          ierr = NF90_CLOSE(nid)
+          ierr = nf90_close(nid)
           IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Pb when closing file', 1)
        ENDIF ! is_mpi_root
Index: /LMDZ6/branches/Amaury_dev/tools/make_sso/make_sso_SpherePack.f90
===================================================================
--- /LMDZ6/branches/Amaury_dev/tools/make_sso/make_sso_SpherePack.f90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/tools/make_sso/make_sso_SpherePack.f90	(revision 5100)
@@ -117,6 +117,6 @@
     SELECT CASE(arg)
       CASE('-i');   f_in=args(k); msg='Missing file "'//TRIM(f_in)//'".'
-                    ll=NF90_OPEN(f_in,NF90_NOWRITE,fID)/=nf90_noerr
-                    IF(.NOT.ll) n=NF90_CLOSE(fID)
+                    ll=nf90_open(f_in,nf90_nowrite,fID)/=nf90_noerr
+                    IF(.NOT.ll) n=nf90_close(fID)
       CASE('-res'); nlon_ou=str2int(args(k)); k=k+1
                     nlat_ou=str2int(args(k))
@@ -135,10 +135,10 @@
   IF(ALL(['noro','spec']/=fmsk)) THEN
     msg='Missing or wrong "-m" option ; can be "noro", "spec" or a mask file'
-    CALL err(NF90_OPEN(fmsk,NF90_NOWRITE,fID)/=nf90_noerr,msg)
+    CALL err(nf90_open(fmsk,nf90_nowrite,fID)/=nf90_noerr,msg)
     CALL nc(NF90_INQ_VARID(fID,"MaskOcean",vID),"MaskOcean")    !--- MASK ID
     CALL nc(NF90_INQUIRE_VARIABLE(fID,vID,dimids=dIDs))         !--- DIMS IDS
     CALL nc(nf90_inquire_dimension(fID,dIDs(1),len=nlon_ou),'x')!--- NB LONG
     CALL nc(nf90_inquire_dimension(fID,dIDs(2),len=nlat_ou),'y')!--- NB LAT
-    CALL nc(NF90_CLOSE(fID))
+    CALL nc(nf90_close(fID))
   END IF
   IF(nlon_ou<=0.OR.nlat_ou<=0) THEN
@@ -147,5 +147,5 @@
 
 !=== READ THE INPUT FIELD =======================================================
-  CALL nc(NF90_OPEN(f_in,NF90_NOWRITE,fID))
+  CALL nc(nf90_open(f_in,nf90_nowrite,fID))
   WRITE(*,*)'>> Reading variable "'//TRIM(vnam)//'" from "'//TRIM(f_in)//'"...'
 
@@ -163,8 +163,8 @@
   CALL nc(nf90_get_var  (fID,loID,lon_in)      ,lonn)
   CALL nc(nf90_get_var  (fID,laID,lat_in)      ,latn)
-  CALL nc(NF90_GET_ATT  (fID,loID,'units',lonu),lonn)
-  CALL nc(NF90_GET_ATT  (fID,laID,'units',latu),latn)
+  CALL nc(nf90_get_att  (fID,loID,'units',lonu),lonn)
+  CALL nc(nf90_get_att  (fID,laID,'units',latu),latn)
   CALL nc(nf90_get_var  (fID,vID,h(:,:),[1,1],[nlon_in+1,nlat_in]),vnam)
-  CALL nc(NF90_CLOSE(fID))
+  CALL nc(nf90_close(fID))
 
 !--- CHECK WETHER GRID IS CORRECT (GLOBAL DOMAIN, IDENTIFIED UNITS...)
@@ -368,8 +368,8 @@
       ALLOCATE(msko(nlon_ou,nlat_ou))
       msg='Missing or wrong "-m" option ; can be "noro", "spec" or a mask file'
-      CALL err(NF90_OPEN(fmsk,NF90_NOWRITE,fID)/=nf90_noerr,msg)
+      CALL err(nf90_open(fmsk,nf90_nowrite,fID)/=nf90_noerr,msg)
       CALL nc(NF90_INQ_VARID(fID,"MaskOcean",vID),"MaskOcean")    !--- MASK ID
       CALL nc(nf90_get_var(fID,vID,msko(:,:)))                    !--- MASK
-      CALL nc(NF90_CLOSE(fID))
+      CALL nc(nf90_close(fID))
       msko(:,:)=1.0-msko(:,:)
   END SELECT
@@ -457,67 +457,67 @@
   f_ou='make_sso_'//TRIM(res_ou)//'_'//TRIM(f_in)
   fnam=f_ou
-  CALL nc(NF90_CREATE(f_ou,NF90_CLOBBER,fID))
-
-  CALL nc(NF90_DEF_DIM(fID,'x',nlon_ou,xID))
-  CALL nc(NF90_DEF_VAR(fID,'x',NF90_REAL,xID,loID)      ,'x')
-  CALL nc(NF90_PUT_ATT(fID,loID,'long_name','Longitude'),'x')
-  CALL nc(NF90_PUT_ATT(fID,loID,'units','degrees_east') ,'x')
-
-  CALL nc(NF90_DEF_DIM(fID,'y',nlat_ou,yID))
-  CALL nc(NF90_DEF_VAR(fID,'y',NF90_REAL,yID,laID)      ,'y')
-  CALL nc(NF90_PUT_ATT(fID,laID,'long_name','Latitude') ,'y')
-  CALL nc(NF90_PUT_ATT(fID,laID,'units','degrees_north'),'y')
-
-  CALL nc(NF90_DEF_VAR(fID,'mask',NF90_REAL,[xID,yID],mskID),'mask')
+  CALL nc(nf90_create(f_ou,nf90_clobber,fID))
+
+  CALL nc(nf90_def_dim(fID,'x',nlon_ou,xID))
+  CALL nc(nf90_def_var(fID,'x',NF90_REAL,xID,loID)      ,'x')
+  CALL nc(nf90_put_att(fID,loID,'long_name','Longitude'),'x')
+  CALL nc(nf90_put_att(fID,loID,'units','degrees_east') ,'x')
+
+  CALL nc(nf90_def_dim(fID,'y',nlat_ou,yID))
+  CALL nc(nf90_def_var(fID,'y',NF90_REAL,yID,laID)      ,'y')
+  CALL nc(nf90_put_att(fID,laID,'long_name','Latitude') ,'y')
+  CALL nc(nf90_put_att(fID,laID,'units','degrees_north'),'y')
+
+  CALL nc(nf90_def_var(fID,'mask',NF90_REAL,[xID,yID],mskID),'mask')
   IF(fmsk=='noro') &
-  CALL nc(NF90_DEF_VAR(fID,'Zphi',NF90_REAL,[xID,yID],phiID),'Zphi')
-  CALL nc(NF90_DEF_VAR(fID,'Zmea',NF90_REAL,[xID,yID],meaID),'Zmea')
-  CALL nc(NF90_DEF_VAR(fID,'mu'  ,NF90_REAL,[xID,yID], muID),'mu'  )
-  CALL nc(NF90_DEF_VAR(fID,'Zsig',NF90_REAL,[xID,yID],sigID),'Zsig')
-  CALL nc(NF90_DEF_VAR(fID,'Zgam',NF90_REAL,[xID,yID],gamID),'Zgam')
-  CALL nc(NF90_DEF_VAR(fID,'Zthe',NF90_REAL,[xID,yID],theID),'Zthe')
-  CALL nc(NF90_DEF_VAR(fID,'Zpic',NF90_REAL,[xID,yID],picID),'Zpic')
-  CALL nc(NF90_DEF_VAR(fID,'Zval',NF90_REAL,[xID,yID],valID),'Zval')
-
-  CALL nc(NF90_PUT_ATT(fID,mskID,'long_name','Fractional land mask'                      ),'mask')
+  CALL nc(nf90_def_var(fID,'Zphi',NF90_REAL,[xID,yID],phiID),'Zphi')
+  CALL nc(nf90_def_var(fID,'Zmea',NF90_REAL,[xID,yID],meaID),'Zmea')
+  CALL nc(nf90_def_var(fID,'mu'  ,NF90_REAL,[xID,yID], muID),'mu'  )
+  CALL nc(nf90_def_var(fID,'Zsig',NF90_REAL,[xID,yID],sigID),'Zsig')
+  CALL nc(nf90_def_var(fID,'Zgam',NF90_REAL,[xID,yID],gamID),'Zgam')
+  CALL nc(nf90_def_var(fID,'Zthe',NF90_REAL,[xID,yID],theID),'Zthe')
+  CALL nc(nf90_def_var(fID,'Zpic',NF90_REAL,[xID,yID],picID),'Zpic')
+  CALL nc(nf90_def_var(fID,'Zval',NF90_REAL,[xID,yID],valID),'Zval')
+
+  CALL nc(nf90_put_att(fID,mskID,'long_name','Fractional land mask'                      ),'mask')
   IF(fmsk=='noro') &
-  CALL nc(NF90_PUT_ATT(fID,phiID,'long_name','Geopotential'                              ),'Zphi')
-  CALL nc(NF90_PUT_ATT(fID,meaID,'long_name','Mean orography'                            ),'Zmea')
-  CALL nc(NF90_PUT_ATT(fID, muID,'long_name','Std deviation of sub-cell scales orography'),'mu'  )
-  CALL nc(NF90_PUT_ATT(fID,sigID,'long_name','Slope along principal axis'                ),'Zsig')
-  CALL nc(NF90_PUT_ATT(fID,gamID,'long_name','Anisotropy (aspect ratio)'                 ),'Zgam')
-  CALL nc(NF90_PUT_ATT(fID,theID,'long_name','Orientation (principal axis)'              ),'Zthe')
-  CALL nc(NF90_PUT_ATT(fID,picID,'long_name','Maximum height'                            ),'Zpic')
-  CALL nc(NF90_PUT_ATT(fID,valID,'long_name','Minimum height'                            ),'Zval')
-
-  CALL nc(NF90_PUT_ATT(fID,mskID,'units','none'   ),'mask')
+  CALL nc(nf90_put_att(fID,phiID,'long_name','Geopotential'                              ),'Zphi')
+  CALL nc(nf90_put_att(fID,meaID,'long_name','Mean orography'                            ),'Zmea')
+  CALL nc(nf90_put_att(fID, muID,'long_name','Std deviation of sub-cell scales orography'),'mu'  )
+  CALL nc(nf90_put_att(fID,sigID,'long_name','Slope along principal axis'                ),'Zsig')
+  CALL nc(nf90_put_att(fID,gamID,'long_name','Anisotropy (aspect ratio)'                 ),'Zgam')
+  CALL nc(nf90_put_att(fID,theID,'long_name','Orientation (principal axis)'              ),'Zthe')
+  CALL nc(nf90_put_att(fID,picID,'long_name','Maximum height'                            ),'Zpic')
+  CALL nc(nf90_put_att(fID,valID,'long_name','Minimum height'                            ),'Zval')
+
+  CALL nc(nf90_put_att(fID,mskID,'units','none'   ),'mask')
   IF(fmsk=='noro') &
-  CALL nc(NF90_PUT_ATT(fID,phiID,'units','m'      ),'Zphi')
-  CALL nc(NF90_PUT_ATT(fID,meaID,'units','m'      ),'Zmea')
-  CALL nc(NF90_PUT_ATT(fID, muID,'units','m'      ),'mu'  )
-  CALL nc(NF90_PUT_ATT(fID,sigID,'units','m/m'    ),'Zsig')
-  CALL nc(NF90_PUT_ATT(fID,gamID,'units','none'   ),'Zgam')
-  CALL nc(NF90_PUT_ATT(fID,theID,'units','degrees'),'Zthe')
-  CALL nc(NF90_PUT_ATT(fID,picID,'units','m'      ),'Zpic')
-  CALL nc(NF90_PUT_ATT(fID,valID,'units','m'      ),'Zval')
-
-  CALL nc(NF90_PUT_ATT(fID,NF90_GLOBAL,'Conventions','COARDS/CF-1.0'))
-  CALL nc(NF90_PUT_ATT(fID,NF90_GLOBAL,'Initial_Grid',TRIM(res_in)))
-  CALL nc(NF90_PUT_ATT(fID,NF90_GLOBAL,'history',TRIM(call_seq)))
-  CALL nc(NF90_ENDDEF(fID))
-
-  CALL nc(NF90_PUT_VAR(fID, loID,lon_ou),'x' )
-  CALL nc(NF90_PUT_VAR(fID, laID,lat_ou),'y' )
-  CALL nc(NF90_PUT_VAR(fID,mskID,msko),'mask')
+  CALL nc(nf90_put_att(fID,phiID,'units','m'      ),'Zphi')
+  CALL nc(nf90_put_att(fID,meaID,'units','m'      ),'Zmea')
+  CALL nc(nf90_put_att(fID, muID,'units','m'      ),'mu'  )
+  CALL nc(nf90_put_att(fID,sigID,'units','m/m'    ),'Zsig')
+  CALL nc(nf90_put_att(fID,gamID,'units','none'   ),'Zgam')
+  CALL nc(nf90_put_att(fID,theID,'units','degrees'),'Zthe')
+  CALL nc(nf90_put_att(fID,picID,'units','m'      ),'Zpic')
+  CALL nc(nf90_put_att(fID,valID,'units','m'      ),'Zval')
+
+  CALL nc(nf90_put_att(fID,nf90_global,'Conventions','COARDS/CF-1.0'))
+  CALL nc(nf90_put_att(fID,nf90_global,'Initial_Grid',TRIM(res_in)))
+  CALL nc(nf90_put_att(fID,nf90_global,'history',TRIM(call_seq)))
+  CALL nc(nf90_enddef(fID))
+
+  CALL nc(nf90_put_var(fID, loID,lon_ou),'x' )
+  CALL nc(nf90_put_var(fID, laID,lat_ou),'y' )
+  CALL nc(nf90_put_var(fID,mskID,msko),'mask')
   IF(fmsk=='noro') &
-  CALL nc(NF90_PUT_VAR(fID,phiID,Zphi),'Zphi')
-  CALL nc(NF90_PUT_VAR(fID,meaID,h0  ),'Zmea')
-  CALL nc(NF90_PUT_VAR(fID, muID,mu  ),'mu'  )
-  CALL nc(NF90_PUT_VAR(fID,sigID,Zsig),'Zsig')
-  CALL nc(NF90_PUT_VAR(fID,gamID,Zgam),'Zgam')
-  CALL nc(NF90_PUT_VAR(fID,theID,Zthe),'Zthe')
-  CALL nc(NF90_PUT_VAR(fID,picID,       h0+2*mu ),'Zpic')
-  CALL nc(NF90_PUT_VAR(fID,valID,MAX(0.,h0-2*mu)),'Zval')
-  CALL nc(NF90_CLOSE(fID))
+  CALL nc(nf90_put_var(fID,phiID,Zphi),'Zphi')
+  CALL nc(nf90_put_var(fID,meaID,h0  ),'Zmea')
+  CALL nc(nf90_put_var(fID, muID,mu  ),'mu'  )
+  CALL nc(nf90_put_var(fID,sigID,Zsig),'Zsig')
+  CALL nc(nf90_put_var(fID,gamID,Zgam),'Zgam')
+  CALL nc(nf90_put_var(fID,theID,Zthe),'Zthe')
+  CALL nc(nf90_put_var(fID,picID,       h0+2*mu ),'Zpic')
+  CALL nc(nf90_put_var(fID,valID,MAX(0.,h0-2*mu)),'Zval')
+  CALL nc(nf90_close(fID))
   WRITE(*,*)'Finished.'
 
Index: /LMDZ6/branches/Amaury_dev/tools/netcdf95/Datasets/nf95_create_single.f90
===================================================================
--- /LMDZ6/branches/Amaury_dev/tools/netcdf95/Datasets/nf95_create_single.f90	(revision 5099)
+++ /LMDZ6/branches/Amaury_dev/tools/netcdf95/Datasets/nf95_create_single.f90	(revision 5100)
@@ -19,5 +19,5 @@
     ! Shortcut to create a file containing a single primary variable.
 
-    use netcdf, only: NF90_CLOBBER, NF90_FLOAT
+    use netcdf, only: nf90_clobber, nf90_float
 
     use nf95_create_m, only: nf95_create
@@ -36,10 +36,10 @@
     !----------------------------------------------------------------------
 
-    call nf95_create(name // ".nc", NF90_CLOBBER, ncid)
+    call nf95_create(name // ".nc", nf90_clobber, ncid)
 
     do i = 1, size(coordinates)
        call nf95_def_dim(ncid, coordinates(i)%name, coordinates(i)%nclen, &
             dimids(i))
-       call nf95_def_var(ncid, coordinates(i)%name, NF90_FLOAT, dimids(i), &
+       call nf95_def_var(ncid, coordinates(i)%name, nf90_float, dimids(i), &
             varid_coord(i))
 
@@ -50,5 +50,5 @@
     END DO
 
-    call nf95_def_var(ncid, name, NF90_FLOAT, dimids, varid)
+    call nf95_def_var(ncid, name, nf90_float, dimids, varid)
 
   end subroutine nf95_create_single
