Index: /utils/tools/ABL_TOOLS/Makefile
===================================================================
--- /utils/tools/ABL_TOOLS/Makefile	(revision 11589)
+++ /utils/tools/ABL_TOOLS/Makefile	(revision 11589)
@@ -0,0 +1,62 @@
+SHELL = /bin/bash
+#-------------------------------------------------------------------------------------------------------------
+EXEC_BIN   = vinterp_abl_frc.exe
+EXEC_BIN2  = uvg_hpg_abl_frc.exe 
+EXEC_BIN3  = drown_abl_frc.exe
+
+NETCDF     := $(NETCDF_DIR)
+NETCDF_LIB := -L$(NETCDF)/lib -lnetcdf -lnetcdff
+NETCDF_INC := -I$(NETCDF)/include
+
+#FC         := gfortran
+#FFLAGS     := -fdefault-real-8 -fdefault-double-8 -O3 $(NETCDF_INC)
+#FFLAGS     := -g -O0 -fcheck=bounds -Wall -Wextra -fbacktrace -finit-real=snan -finit-integer=8888 -finit-character=90 $(NETCDF_INC)
+
+FC         := ifort
+FFLAGS     := -check noarg_temp_created -autodouble -O3 $(NETCDF_INC)
+#FFLAGS     := -check noarg_temp_created -autodouble -O0 -g -traceback -check all $(NETCDF_INC)
+
+#---------------------------------------------------------------------  
+SRC_DIR        =  $(PWD)
+OBJ_DIR        =  $(PWD)
+FINAL_OBJS     =  $(OBJ_DIR)/main_vinterp.o $(OBJ_DIR)/module_io.o $(OBJ_DIR)/module_interp.o $(OBJ_DIR)/module_grid.o
+FINAL_OBJS2    =  $(OBJ_DIR)/main_uvg_hpg.o $(OBJ_DIR)/module_io.o $(OBJ_DIR)/module_grid.o
+FINAL_OBJS3    =  $(OBJ_DIR)/main_drown.o $(OBJ_DIR)/module_io.o  $(OBJ_DIR)/module_grid.o
+#---------------------------------------------------------------------
+
+all: | $(EXEC_BIN) $(EXEC_BIN2) $(EXEC_BIN3)
+	@echo 
+	@echo ===================================================
+	@echo     ABL Preprocessing tools for ECMWF data: OK
+	@echo ===================================================
+	@echo
+
+clean:
+	$(RM)    $(OBJ_DIR)/*.o   *.mod 
+	$(RM)    $(EXEC_BIN) $(EXEC_BIN2) $(EXEC_BIN3)
+	
+$(EXEC_BIN): $(FINAL_OBJS) 
+	$(FC) -o $(EXEC_BIN) $(FINAL_OBJS) $(NETCDF_LIB)
+
+$(EXEC_BIN2): $(FINAL_OBJS2) 
+	$(FC) -o $(EXEC_BIN2) $(FINAL_OBJS2) $(NETCDF_LIB)
+
+$(EXEC_BIN3): $(FINAL_OBJS3) 
+	$(FC) -o $(EXEC_BIN3) $(FINAL_OBJS3) $(NETCDF_LIB)
+
+# Rules
+.DEFAULT:
+
+.SUFFIXES:
+
+.PRECIOUS: $(SRC_DIR)/%.F90
+	
+$(OBJ_DIR)/%.o : $(SRC_DIR)/%.F90
+	$(FC) $(FFLAGS) -o $@ -c $<
+
+$(OBJ_DIR)/main_vinterp.o:	$(OBJ_DIR)/module_grid.o $(OBJ_DIR)/module_io.o $(OBJ_DIR)/module_interp.o
+$(OBJ_DIR)/main_uvg_hpg.o:	$(OBJ_DIR)/module_grid.o $(OBJ_DIR)/module_io.o
+$(OBJ_DIR)/main_drown.o:	$(OBJ_DIR)/module_grid.o $(OBJ_DIR)/module_io.o
+$(OBJ_DIR)/module_io.o:
+$(OBJ_DIR)/module_interp.o:
+$(OBJ_DIR)/module_grid.o: $(OBJ_DIR)/module_io.o
Index: /utils/tools/ABL_TOOLS/main_drown.F90
===================================================================
--- /utils/tools/ABL_TOOLS/main_drown.F90	(revision 11589)
+++ /utils/tools/ABL_TOOLS/main_drown.F90	(revision 11589)
@@ -0,0 +1,347 @@
+PROGRAM main
+   !!======================================================================
+   !!                     ***  PROGRAM main  ***
+   !!
+   !! ** Purpose : horizontal extrapolation ("drowning") of ECMWF dataset
+   !!
+   !! 
+   !!======================================================================
+   !! History : 2016-10  (F. Lemarié)  Original code largely inspired by SOSIE (L. Brodeau & J.M. Molines)
+   !!   
+   !!----------------------------------------------------------------------
+   USE module_io       ! I/O routines
+   USE module_grid     ! compute input and output grids 
+   !!
+   IMPLICIT NONE
+   !!----------------------------------------------------------------------
+   !! 
+   !!  
+   !! 
+   !!----------------------------------------------------------------------
+   !
+   INTEGER                                  :: ji,jj,jk,kt,kv,jjp1,jjm1,jip1,jim1
+   INTEGER                                  :: jpka,jpvar            ! number of vertical levels for input and target grids 
+   INTEGER                                  :: jpi , jpj             ! number of grid points in x and y directions                 
+   INTEGER                                  :: jptime,ctrl,niter,status,ncid
+   INTEGER                                  :: ioerr
+   INTEGER, PARAMETER                       :: stdout   = 6
+   INTEGER, PARAMETER                       :: max_iter = 50
+   !!
+   REAL(8)                                  :: cff,cff1,cff2
+   !!  
+   REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: varin
+   REAL(8), ALLOCATABLE, DIMENSION(:,:    ) :: tmask, tmask2, ff_t
+   REAL(8), ALLOCATABLE, DIMENSION(:,:    ) :: mask_coast, maskv              ! land-sea mask   
+   REAL(8), ALLOCATABLE, DIMENSION(:,:    ) :: varold,varnew   
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: tmp1d, lat
+   !!
+   CHARACTER(len=500)                       :: file_u,file_v,file_hpg,file_geos ! ERAi files containing wind components 
+   CHARACTER(len=500)                       :: file_t,file_q, file_m   ! ERAi files containing tracers and mask 
+   CHARACTER(len=500)                       :: file_z,file_p,cn_dir    ! ERAi files containing surface geopot and pressure
+   CHARACTER(len=500)                       :: abl_file, grd_file, drwn_file, out_file
+   CHARACTER(len=500)                       :: namelistf              
+   CHARACTER(len=500)                       :: argument   
+   CHARACTER(len= 20), DIMENSION(4)         :: dimnames
+   CHARACTER(:), ALLOCATABLE, DIMENSION(:)  :: varnames
+   CHARACTER(6)                             :: mask_var            ! name of mask variable in file_m file  
+   CHARACTER(6)                             :: var_name
+   !!
+   LOGICAL                                  :: ln_read_zsurf      ! read surface geopotential or not
+   LOGICAL                                  :: ln_read_mask       ! read land-sea mask or not
+   LOGICAL                                  :: ln_perio_latbc     ! use periodic BC along the domain latitudinal edges (for global data) or use zero-gradient BC (for regional data)
+   LOGICAL                                  :: ln_c1d             ! output only a single column in output file
+   LOGICAL                                  :: ln_hpg_frc         ! compute horizontal pressure gradient
+   LOGICAL                                  :: ln_geo_wnd         ! compute goestrophic wind components 
+   LOGICAL                                  :: ln_slp_smth        ! apply gibbs oscillation filetring on mean sea level pressure
+   LOGICAL                                  :: ln_drw_smth        ! apply gibbs oscillation filetring on mean sea level pressure
+   LOGICAL                                  :: ln_slp_log         ! log(sea-level pressure) or sea-level pressure
+   LOGICAL                                  :: ln_lsm_land        ! if T mask is 1 over land and 0 over ocean if F it is the other way around
+   INTEGER                                  :: ptemp_method       ! way to compute potential temperature 
+   !!
+   REAL(8), PARAMETER     :: omega =   7.292116e-05
+   REAL(8), PARAMETER     :: rad   =   3.141592653589793 / 180.
+
+
+   !!---------------------------------------------------------------------
+   !! List of variables read in the namelist file 
+   NAMELIST/nml_out/    grd_file, abl_file, drwn_file, var_name
+   NAMELIST/nml_opt/    ptemp_method, ln_slp_log, ln_slp_smth, ln_read_mask, ln_perio_latbc, &
+      &                 ln_hpg_frc, ln_geo_wnd, ln_c1d, ln_read_zsurf, ln_lsm_land, ln_drw_smth
+   NAMELIST/nml_fld/    cn_dir, file_u, file_v, file_t,              &
+      &                 file_q, file_z, file_p, file_hpg, file_geos, &
+      &                 file_m, mask_var
+   !!
+   !! get the namelist file name
+   CALL get_command_argument(1, argument, ctrl, status)
+   !
+   SELECT CASE(status)
+   CASE(0)
+      namelistf = trim(argument)
+   CASE(-1)
+      WRITE(stdout,*) "### Error: file name too long"
+      STOP
+   CASE DEFAULT
+      namelistf = 'namelist_abl_tools'
+   END SELECT
+   !!---------------------------------------------------------------------
+
+
+   !!---------------------------------------------------------------------   
+   !! read namelist variables
+   ctrl = 0
+   OPEN(50, file=namelistf, status='old', form='formatted', access='sequential', iostat=ioerr)
+   IF (ioerr /= 0) ctrl = ctrl + 1   
+   READ(50,nml_opt, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
+   READ(50,nml_fld, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
+   READ(50,nml_out, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
+   
+   IF (ctrl > 0) then
+      WRITE(stdout,*) "### E R R O R while reading namelist file '",trim(namelistf),"'"
+      WRITE(stdout,*) " ctrl = ",ctrl
+      STOP
+   ELSE 
+      WRITE(stdout,*) " Namelist file ",trim(namelistf)," OK "    
+   END IF
+   !!-------------------------------------------------------------------------------------   
+   !! list of variables to treat
+   !!
+   !! get the variable name
+   CALL get_command_argument( 2, argument, ctrl, status)
+   SELECT CASE(status)
+   CASE(0)
+      var_name = trim(argument)
+   CASE(-1)
+      WRITE(stdout,*) "### Error: file name too long"
+      STOP
+   END SELECT
+   WRITE(stdout,*) "var_name: ", trim(var_name), Len_Trim(var_name)
+   !
+   IF (Len_Trim(var_name) == 0) THEN
+     jpvar    = 7
+     allocate(character(len=20) :: varnames(jpvar) )
+     varnames(1:5) = [character(len=4) :: 'slp', 'uwnd', 'vwnd', 'tair', 'humi' ]
+     IF (ln_hpg_frc) varnames(6:7) = [character(len=4) :: 'uhpg', 'vhpg' ]
+     IF (ln_geo_wnd) varnames(6:7) = [character(len=4) :: 'ugeo', 'vgeo' ]
+     IF( .NOT. Var_Existence( varnames(5), trim(cn_dir)//'/'//trim(abl_file) ) ) jpvar = jpvar - 1
+     IF( .NOT. Var_Existence( varnames(6), trim(cn_dir)//'/'//trim(abl_file) ) ) jpvar = jpvar - 1 
+   ELSE
+     jpvar     = 1 
+     allocate(character(len=20) :: varnames(jpvar) )
+     varnames  = var_name
+     abl_file  = trim(var_name)//'_'//trim( abl_file)
+     drwn_file = trim(var_name)//'_'//trim(drwn_file)
+   END IF 
+   WRITE(stdout,*) "Number of variables to treat : ",jpvar    
+   !!-------------------------------------------------------------------------------------    
+   !! read the dimensions for the input files
+   CALL Read_Ncdf_dim ( 'jpka' , trim(cn_dir)//'/'//trim(abl_file), jpka    ) 
+   CALL Read_Ncdf_dim ( 'time' , trim(cn_dir)//'/'//trim(abl_file), jptime  ) 
+   CALL Read_Ncdf_dim ( 'lon'  , trim(cn_dir)//'/'//trim(abl_file), jpi     )     
+   CALL Read_Ncdf_dim ( 'lat'  , trim(cn_dir)//'/'//trim(abl_file), jpj     )        
+   WRITE(stdout,*) "jpka, jptime, jpi, jpj: ", jpka, jptime, jpi, jpj
+   !
+   !!---------------------------------------------------------------------
+
+
+   !!---------------------------------------------------------------------
+   !! allocate arrays  
+   ALLOCATE(varin (1:jpi,1:jpj,1:jpka,1))
+   ALLOCATE(tmask (1:jpi,1:jpj    ))
+   ALLOCATE(tmask2(1:jpi,1:jpj    ))
+   ALLOCATE(maskv      (1:jpi,1:jpj )) 
+   ALLOCATE(mask_coast (1:jpi,1:jpj ))      
+   ALLOCATE(varold     (1:jpi,1:jpj ))  
+   ALLOCATE(varnew     (1:jpi,1:jpj ))     
+
+   !!---------------------------------------------------------------------
+   !! Read the mask and remove some closed seas 
+   IF (ln_read_mask) THEN
+     CALL init_atm_mask(jpi,jpj,trim(cn_dir)//'/'//trim(file_m),trim(mask_var),ln_lsm_land, tmask )
+     CALL Read_Ncdf_var ( 'tmask' , trim(cn_dir)//'/'//trim(abl_file), tmask2(:,:) )
+   ELSE
+     tmask(:,:)  = 1.
+     tmask2(:,:) = 1.
+   END IF
+   !!   
+
+
+
+   !!---------------------------------------------------------------------   
+   !! create output file
+   !!
+   out_file = trim(cn_dir)//'/'//trim(drwn_file)
+   CALL Init_output_File ( jpi, jpj, jpka-1, trim(cn_dir)//'/'//trim(abl_file), out_file, tmask(:,:) )
+
+   !!---------------------------------------------------------------------   
+   !! Initialize the name of the dimensions for the result of the drowning
+   !!
+   dimnames(1) = 'lon'
+   dimnames(2) = 'lat'       
+   dimnames(3) = 'jpka'
+   dimnames(4) = 'time'  
+
+   CALL Write_Ncdf_var( 'tmask', dimnames(1:2), trim(out_file), tmask2, 'float' )
+
+   !!---------------------------------------------------------------------
+   ! Read time variable
+   ALLOCATE(tmp1d (1:jptime))    
+   ALLOCATE(lat   (1:jpj))    
+   CALL Read_Ncdf_var ( 'time', trim(cn_dir)//'/'//trim(file_t), tmp1d )    
+   CALL Read_Ncdf_var ( 'lat' , trim(cn_dir)//'/'//trim(file_t), lat )      !<-- latitude  
+   !!---------------------------------------------------------------------     
+
+   ! force drowning in the equatorial band with geo wind
+   IF (ln_geo_wnd) THEN
+     ALLOCATE( ff_t(1:jpi,1:jpj) )
+     DO jj = 1, jpj
+        DO ji = 1, jpi
+           ff_t(ji,jj) = 2. * omega * SIN( rad * lat(jj) )
+           IF (abs(ff_t(ji,jj)) < 2.5e-5) tmask2(ji,jj) = 0. 
+        END DO
+     END DO   
+   END IF
+
+   !===========
+   DO kv = 1,jpvar
+   !===========
+      DO kt = 1,jptime
+      !===========
+
+         IF( kv == 1 ) THEN
+           CALL Write_Ncdf_var( 'time', dimnames(4:4), trim(out_file), tmp1d(kt:kt), kt, 'double' )
+         ENDIF
+
+         ! Read variable to treat
+         print*,'Treat variable ',trim(varnames(kv)),' at time level ',kt
+         IF (trim(varnames(kv)).EQ."slp") THEN
+           CALL Read_Ncdf_var ( trim(varnames(kv)), trim(cn_dir)//'/'//trim(file_hpg), varin(:,:,1,1),  kt )
+         ELSE
+           CALL Read_Ncdf_var ( trim(varnames(kv)), trim(cn_dir)//'/'//trim(abl_file), varin,  kt )
+         END IF
+
+
+         !===========
+         DO jk = 1,jpka
+         !===========
+
+            IF ( (trim(varnames(kv))=="uwnd").OR.(trim(varnames(kv))=="vwnd").OR. &
+               & (trim(varnames(kv))=="ugeo").OR.(trim(varnames(kv))=="vgeo").OR. &
+               & (trim(varnames(kv))=="uhpg").OR.(trim(varnames(kv))=="vhpg") ) THEN
+                 maskv (1:jpi,1:jpj     ) = tmask2(1:jpi,1:jpj )
+            ELSE
+                 maskv (1:jpi,1:jpj     ) = tmask(1:jpi,1:jpj )
+            END IF
+
+
+            varold (1:jpi,1:jpj) = varin(1:jpi,1:jpj,jk,1)
+            varnew (1:jpi,1:jpj) = varin(1:jpi,1:jpj,jk,1)            
+            mask_coast(1:jpi,1:jpj) = 0.
+
+            !$$$$$$$$$$$$$$$$
+            DO niter = 1, max_iter
+            !$$$$$$$$$$$$$$$$                 
+
+               varold(1:jpi,1:jpj) = varnew (1:jpi,1:jpj)
+
+               IF ( .NOT. (ANY(maskv == 0))  ) THEN
+                  EXIT
+               END IF
+
+               !+++++++++++++++++++++++++++++++++++
+               ! Build mask_coast
+               DO jj = 1,jpj
+                  DO ji = 1,jpi
+
+                     IF (ln_perio_latbc) THEN
+                       jip1 = ji + 1; if(ji==jpi) jip1 = 1
+                       jim1 = ji - 1; if(ji==  1) jim1 = jpi
+                     ELSE
+                       jip1 = ji + 1; if(ji==jpi) jip1 = jpi-1
+                       jim1 = ji - 1; if(ji==  1) jim1 = 2
+                     END IF
+                     jjp1 = jj + 1; if(jj==jpj) jjp1 = jpj-1
+                     jjm1 = jj - 1; if(jj==  1) jjm1 = 2
+
+                     cff  = (1.-maskv(ji,jj))
+                     cff1 = maskv(ji,jjp1)+maskv(jip1,jj)+maskv(jim1,jj)+maskv(ji,jjm1) 
+                     mask_coast(ji,jj) = cff * cff1 
+                     IF( mask_coast(ji,jj) > 0. ) mask_coast(ji,jj) = 1.
+
+                  END DO
+               END DO   
+               !+++++++++++++++++++++++++++++++++++          
+
+               DO jj=1,jpj
+                  DO ji=1,jpi
+
+                     IF ( mask_coast(ji,jj) == 1. ) THEN
+
+                        IF (ln_perio_latbc) THEN
+                          jip1 = ji + 1; if(ji==jpi) jip1 = 1
+                          jim1 = ji - 1; if(ji==  1) jim1 = jpi
+                        ELSE
+                          jip1 = ji + 1; if(ji==jpi) jip1 = jpi-1
+                          jim1 = ji - 1; if(ji==  1) jim1 = 2
+                        END IF
+                        jjp1 = jj + 1; if(jj==jpj) jjp1 = jpj-1
+                        jjm1 = jj - 1; if(jj==  1) jjm1 = 2
+
+                        cff  = maskv(jim1,jjm1)+maskv(jim1,jj)+maskv(jim1,jjp1) &
+                           & + maskv(ji  ,jjm1)+               maskv(ji  ,jjp1) &
+                           & + maskv(jip1,jjm1)+maskv(jip1,jj)+maskv(jip1,jjp1)
+
+                        varnew(ji,jj) = (1./cff)*(                                   &
+                           &                    varold(jim1,jjm1)*maskv(jim1,jjm1)   &
+                           &              +     varold(jim1,jj  )*maskv(jim1,jj  )   &
+                           &              +     varold(jim1,jjp1)*maskv(jim1,jjp1)   &
+                           &              +     varold(ji  ,jjm1)*maskv(ji  ,jjm1)   &
+                           &              +     varold(ji  ,jjp1)*maskv(ji  ,jjp1)   &                
+                           &              +     varold(jip1,jjm1)*maskv(jip1,jjm1)   &
+                           &              +     varold(jip1,jj  )*maskv(jip1,jj  )   &
+                           &              +     varold(jip1,jjp1)*maskv(jip1,jjp1)   )
+
+                     END IF
+
+                  END DO
+               END DO   
+               !+++++++++++++++++++++++++++++++++++
+               maskv(1:jpi,1:jpj) = maskv(1:jpi,1:jpj) + mask_coast(1:jpi,1:jpj)
+               !+++++++++++++++++++++++++++++++++++                   
+           !$$$$$$$$$$$$$$$$
+           END DO
+           !$$$$$$$$$$$$$$$$
+
+           ! SMOOTHING AFTER DROWNING
+           IF (ln_drw_smth) THEN
+             maskv (1:jpi,1:jpj     ) = 1.
+           ELSE
+             !!only over continent
+             maskv (1:jpi,1:jpj     ) = 1. - tmask(1:jpi,1:jpj  )
+           END IF        
+
+           CALL smooth_field( jpi, jpj, varnew, maskv, 3 )
+           varin(1:jpi,1:jpj,jk,1) = varnew(1:jpi,1:jpj)
+
+           IF (trim(varnames(kv)).EQ."slp") EXIT
+
+        !===========
+        END DO ! jk
+        !===========
+
+        IF (trim(varnames(kv)).EQ."slp") THEN
+          CALL Write_Ncdf_var ( trim(varnames(kv)), dimnames(1:4), trim(cn_dir)//'/'//drwn_file, varin(:,:,1,1), kt, 'float' )
+        ELSE
+          CALL Write_Ncdf_var ( trim(varnames(kv)), dimnames(1:4), trim(cn_dir)//'/'//drwn_file, varin(:,:,:,:), kt, 'float' )
+        END IF
+
+      !===========
+      END DO ! time
+      !===========
+
+   !===========
+   END DO ! variable
+   !===========
+
+   STOP
+   !
+END PROGRAM main
Index: /utils/tools/ABL_TOOLS/main_uvg_hpg.F90
===================================================================
--- /utils/tools/ABL_TOOLS/main_uvg_hpg.F90	(revision 11589)
+++ /utils/tools/ABL_TOOLS/main_uvg_hpg.F90	(revision 11589)
@@ -0,0 +1,647 @@
+PROGRAM main
+   !!======================================================================
+   !!                     ***  PROGRAM main  ***
+   !!
+   !! ** Purpose : compute geostrophic wind or horizontal pressure gradient
+   !!              from ECMWF atmospheric model vertical levels altitude, 
+   !!              temperature and humidity 3D fields
+   !! 
+   !!======================================================================
+   !! History : 2016-10  (F. Lemarié)  Original code
+   !!   
+   !!----------------------------------------------------------------------
+   USE module_io       ! I/O routines
+   USE module_grid     ! compute input and output grids 
+   !!
+   IMPLICIT NONE
+   !!----------------------------------------------------------------------
+   !! 
+   !!  
+   !! 
+   !!----------------------------------------------------------------------
+   !
+   INTEGER                                  :: ji,jj,jk,kt, nhym, nhyi
+   INTEGER                                  :: jpka                  ! number of vertical levels for input and target grids 
+   INTEGER                                  :: jpi , jpj             ! number of grid points in x and y directions     
+   INTEGER                                  :: status             
+   INTEGER                                  :: jptime,ctrl
+   INTEGER                                  :: ioerr,ncid
+   INTEGER, ALLOCATABLE, DIMENSION(:,:  )   :: ind  
+   INTEGER, PARAMETER                       :: stdout  = 6
+   !!
+   REAL(8)                                  :: cff,tv
+   !!  
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: lev                ! A coefficients to reconstruct ECMWF grid
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: A_w                ! A coefficients to reconstruct ECMWF grid
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: B_w                ! B coefficients to reconstruct ECMWF grid
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: B_r                ! B coefficients to reconstruct ECMWF grid   
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: A_r                ! A coefficients to reconstruct ECMWF grid
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: ph,lon,lat   
+   REAL(8), ALLOCATABLE, DIMENSION(:,:,:  ) :: e3t,ghw                ! thickness of vertical layers in target grid
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: tmp1d, tmp_fullw, tmp_fullm ! temporary/working 1D arrays
+   REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: humi
+   REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: temp
+   REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: uhpg
+   REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: vhpg
+   REAL(8), ALLOCATABLE, DIMENSION(:,:    ) :: slp, zsurf, zsurf_smth, slp_smth, ghw_smth !, slp_mask, ghw_mask
+   REAL(8), ALLOCATABLE, DIMENSION(:,:    ) :: dx,dy,ff_t,tmask, tmask2
+   REAL(8), ALLOCATABLE, DIMENSION(:,:    ) :: FX,FE,wrkx,wrke
+   REAL(8), ALLOCATABLE, DIMENSION(:,:    ) :: dZx,dZe
+   !!
+   CHARACTER(len=500)                       :: file_u,file_v,file_hpg, file_geos ! ECMWF files containing wind components 
+   CHARACTER(len=500)                       :: file_t,file_q, file_m             ! ECMWF files containing tracers and mask 
+   CHARACTER(len=500)                       :: file_z,file_p,cn_dir              ! ECMWF files containing surface geopot and pressure
+   CHARACTER(len=500)                       :: out_file      
+   CHARACTER(len=500)                       :: namelistf
+   CHARACTER(len=500)                       :: argument   
+   CHARACTER(len= 20),DIMENSION(4)          :: dimnames
+   CHARACTER(len= 20),DIMENSION(4)          :: varnames
+   CHARACTER(6)                             :: mask_var           ! name of mask variable in file_m file  
+   !!
+   LOGICAL                                  :: ln_read_zsurf      ! read surface geopotential or not
+   LOGICAL                                  :: ln_read_mask       ! read land-sea mask or not
+   LOGICAL                                  :: ln_perio_latbc     ! use periodic BC along the domain latitudinal edges (for global data) or use zero-gradient BC (for regional data)
+   LOGICAL                                  :: ln_c1d             ! output only a single column in output file
+   LOGICAL                                  :: ln_hpg_frc         ! compute horizontal pressure gradient
+   LOGICAL                                  :: ln_geo_wnd         ! compute goestrophic wind components 
+   LOGICAL                                  :: ln_slp_smth        ! apply gibbs oscillation filetring on mean sea level pressure
+   LOGICAL                                  :: ln_drw_smth        ! apply gibbs oscillation filetring on mean sea level pressure
+   LOGICAL                                  :: ln_slp_log         ! log(sea-level pressure) or sea-level pressure
+   LOGICAL                                  :: ln_lsm_land        ! if T mask is 1 over land and 0 over ocean if F it is the other way around
+   LOGICAL                                  :: ln_impose_z1       ! impose the altitude of the first level in target grid
+   INTEGER                                  :: ptemp_method       ! way to compute potential temperature 
+                                                                  ! = 0  (absolute temperature)
+                                                                  ! = 1  (potential temperature with local ref pressure)
+                                                                  ! = 2  (potential temperature with global ref pressure on temperature perturbation)
+                                                                  ! = 3  (potential temperature with global ref pressure)
+   !!   
+   REAL(8), PARAMETER     :: grav  =   9.80665
+   REAL(8), PARAMETER     :: Rd    = 287.058   
+   REAL(8), PARAMETER     :: zvir  =   0.609133 
+   REAL(8), PARAMETER     :: omega =   7.292116e-05
+   REAL(8), PARAMETER     :: rad   =   3.141592653589793 / 180.
+   REAL(8), PARAMETER     :: rt    = 6371229.
+   REAL(8), PARAMETER     :: lat_smth = 0. !65.   !!GS: possibility to smooth only above a selected latitude
+
+
+   !!---------------------------------------------------------------------
+   !! List of variables read in the namelist file 
+   NAMELIST/nml_opt/    ptemp_method, ln_slp_log, ln_slp_smth, ln_read_mask, ln_perio_latbc, & 
+      &                 ln_hpg_frc, ln_geo_wnd, ln_c1d, ln_read_zsurf, ln_lsm_land, ln_drw_smth
+   NAMELIST/nml_fld/    cn_dir, file_u, file_v, file_t,               &
+      &                 file_q, file_z, file_p, file_hpg, file_geos,  &
+      &                 file_m, mask_var
+   !!
+   !! get the namelist file name
+   CALL get_command_argument( 1, argument, ctrl, status)
+   !
+   SELECT CASE(status)
+   CASE(0)
+      namelistf = trim(argument)
+   CASE(-1)
+      WRITE(stdout,*) "### Error: file name too long"
+      STOP
+   CASE DEFAULT
+      namelistf = 'namelist_abl_tools'
+   END SELECT
+   !!---------------------------------------------------------------------
+
+
+   !!---------------------------------------------------------------------
+   !! read namelist variables
+   ctrl = 0
+   OPEN(50, file=namelistf, status='old', form='formatted', access='sequential', iostat=ioerr)
+   IF (ioerr /= 0) ctrl = ctrl + 1   
+   READ(50,nml_opt, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1   
+   READ(50,nml_fld, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
+   
+   IF (ctrl > 0) then
+      WRITE(stdout,*) "### E R R O R while reading namelist file '",trim(namelistf),"'"
+      WRITE(stdout,*) " ctrl = ",ctrl
+      STOP
+   ELSE 
+      WRITE(stdout,*) " Namelist file ",trim(namelistf), " OK "    
+   END IF
+
+   IF(ln_hpg_frc) THEN
+      out_file = trim(cn_dir)//'/'//trim(file_hpg)
+   ELSE IF(ln_geo_wnd) THEN
+      out_file = trim(cn_dir)//'/'//trim(file_geos)   
+   ELSE
+      WRITE(stdout,*) "### E R R O R in namelist variable "
+      WRITE(stdout,*) "either ln_hpg_frc or ln_geo_wnd should be set to True"
+      STOP       
+   END IF
+   !!---------------------------------------------------------------------   
+       
+       
+   !!---------------------------------------------------------------------
+   ! check files content
+   CALL Read_Ncdf_dim('lev',trim(cn_dir)//'/'//trim(file_t),jpka)   
+   !                          geop_surf  temp  humi  pressure
+   varnames = [ character(len=3) :: 'Z', 'T', 'Q', 'MSL' ]
+   IF (ln_read_zsurf) varnames(4) = "SP"
+   IF (ln_slp_log)    varnames(4) = "LNSP"
+
+   ctrl = 0
+   IF (ln_read_zsurf) THEN
+     IF( .not. VAR_EXISTENCE( trim(varnames(1)) , trim(cn_dir)//'/'//trim(file_z) ) ) ctrl = ctrl + 1 
+   END IF
+   IF ( .not. VAR_EXISTENCE( trim(varnames(2)) , trim(cn_dir)//'/'//trim(file_t) ) ) ctrl = ctrl + 1 
+   IF ( .not. VAR_EXISTENCE( trim(varnames(3)) , trim(cn_dir)//'/'//trim(file_q) ) ) ctrl = ctrl + 1   
+   IF ( .not. VAR_EXISTENCE( trim(varnames(4)) , trim(cn_dir)//'/'//trim(file_p) ) ) ctrl = ctrl + 1
+   WRITE(*,*) " pressure variable name: ", varnames(4)
+ 
+   IF ( ctrl > 0 ) THEN  
+      WRITE(stdout,*) "### E R R O R while reading ECMWF atmospheric files "
+      STOP
+   ELSE
+      WRITE(stdout,*) " ECMWF atmospheric files OK "        
+   END IF
+   !!---------------------------------------------------------------------
+   
+
+   !!---------------------------------------------------------------------
+   !! read the dimensions for the input files
+   CALL Read_Ncdf_dim ( 'time', trim(cn_dir)//'/'//trim(file_t), jptime  ) 
+   CALL Read_Ncdf_dim ( 'lon' , trim(cn_dir)//'/'//trim(file_t), jpi     )     
+   CALL Read_Ncdf_dim ( 'lat' , trim(cn_dir)//'/'//trim(file_t), jpj     )     
+   CALL Read_Ncdf_dim ( 'nhym' , trim(cn_dir)//'/'//trim(file_t), nhym   )
+   CALL Read_Ncdf_dim ( 'nhyi' , trim(cn_dir)//'/'//trim(file_t), nhyi   )
+   !
+   !!---------------------------------------------------------------------
+
+
+   !!---------------------------------------------------------------------
+   !! allocate arrays  
+   ALLOCATE( A_w    (               0:jpka   ) )         
+   ALLOCATE( B_w    (               0:jpka   ) )     
+   ALLOCATE( B_r    (               1:jpka   ) )
+   ALLOCATE( A_r    (               1:jpka   ) )   
+
+   ALLOCATE( e3t    ( 1:jpi, 1:jpj, 1:jpka   ) )
+   ALLOCATE( ghw    ( 1:jpi, 1:jpj, 0:jpka   ) ) 
+   ALLOCATE( slp_smth(1:jpi, 1:jpj           ) )
+   ALLOCATE( ghw_smth(1:jpi, 1:jpj           ) )
+   ALLOCATE( slp    ( 1:jpi, 1:jpj           ) )
+   ALLOCATE( zsurf  ( 1:jpi, 1:jpj           ) )
+   ALLOCATE( zsurf_smth  ( 1:jpi, 1:jpj           ) )
+   ALLOCATE( temp   ( 1:jpi, 1:jpj, 1:jpka, 1) )
+   ALLOCATE( humi   ( 1:jpi, 1:jpj, 1:jpka, 1) )   
+   ALLOCATE( uhpg   ( 1:jpi, 1:jpj, 1:jpka, 1) )
+   ALLOCATE( vhpg   ( 1:jpi, 1:jpj, 1:jpka, 1) )   
+   ALLOCATE( ph     (               0:jpka   ) ) 
+   ALLOCATE( dx     ( 1:jpi, 1:jpj           ) ) 
+   ALLOCATE( dy     ( 1:jpi, 1:jpj           ) )         
+   ALLOCATE( FX     ( 1:jpi, 1:jpj           ) ) 
+   ALLOCATE( FE     ( 1:jpi, 1:jpj           ) ) 
+   ALLOCATE( dzx    ( 0:jpi, 1:jpj           ) ) 
+   ALLOCATE( dze    ( 1:jpi, 0:jpj           ) ) 
+   ALLOCATE( wrkx   ( 1:jpi, 1:jpj           ) )      
+   ALLOCATE( wrke   ( 1:jpi, 1:jpj           ) )  
+   ALLOCATE( tmask  ( 1:jpi, 1:jpj           ) ) 
+   ALLOCATE( tmask2 ( 1:jpi, 1:jpj           ) ) 
+   IF (jpka.NE.nhym) THEN
+     ALLOCATE( tmp_fullw(1:nhyi) )
+     ALLOCATE( tmp_fullm(1:nhym) )
+   END IF
+
+   !!---------------------------------------------------------------------
+   !! Read the mask and remove some closed seas 
+   IF (ln_read_mask) THEN
+     CALL init_atm_mask( jpi, jpj, trim(cn_dir)//'/'//trim(file_m), trim(mask_var), ln_lsm_land, tmask)
+   ELSE
+     tmask(:,:) = 1.
+   END IF
+   tmask2 = tmask
+   !!   
+   
+   !! Read the static A and B coefficients for the ECMWF vertical grid
+   IF (jpka.EQ.nhym) THEN
+     CALL Read_Ncdf_var    ( 'hyai', trim(cn_dir)//'/'//trim(file_t), A_w )
+     CALL Read_Ncdf_var    ( 'hybi', trim(cn_dir)//'/'//trim(file_t), B_w )   
+     CALL Read_Ncdf_var    ( 'hyam', trim(cn_dir)//'/'//trim(file_t), A_r ) 
+     CALL Read_Ncdf_var    ( 'hybm', trim(cn_dir)//'/'//trim(file_t), B_r ) 
+   ELSE
+     CALL Read_Ncdf_var    ( 'hyai', trim(cn_dir)//'/'//trim(file_t), tmp_fullw )
+     A_w(0:jpka) = tmp_fullw(nhyi-(jpka+1)+1:nhyi)
+     CALL Read_Ncdf_var    ( 'hybi', trim(cn_dir)//'/'//trim(file_t), tmp_fullw )
+     B_w(0:jpka) = tmp_fullw(nhyi-(jpka+1)+1:nhyi)
+     CALL Read_Ncdf_var    ( 'hyam', trim(cn_dir)//'/'//trim(file_t), tmp_fullm )
+     A_r(1:jpka) = tmp_fullm(nhym-jpka+1:nhym)
+     CALL Read_Ncdf_var    ( 'hybm', trim(cn_dir)//'/'//trim(file_t), tmp_fullm )
+     B_r(1:jpka) = tmp_fullm(nhym-jpka+1:nhym)
+   END IF
+
+   ALLOCATE(lat(1:jpj),lon(1:jpi),ff_t(1:jpi,1:jpj),lev(1:jpka))            
+   CALL Read_Ncdf_var ( 'lon' , trim(cn_dir)//'/'//trim(file_t), lon )      !<-- longitude
+   CALL Read_Ncdf_var ( 'lat' , trim(cn_dir)//'/'//trim(file_t), lat )      !<-- latitude  
+   CALL Read_Ncdf_var ( 'lev' , trim(cn_dir)//'/'//trim(file_t), lev )    
+
+
+   !!---------------------------------------------------------------------
+   !++ Compute Coriolis frequency at cell centers
+   !    
+   DO jj = 1, jpj
+      DO ji = 1, jpi
+         ff_t(ji,jj) = 2. * omega * SIN( rad * lat(jj) )
+      END DO
+   END DO   
+
+
+   !!---------------------------------------------------------------------
+   !++ Compute dx and dy at cell centers
+   !
+   dx(:,:) = 0.
+   dy(:,:) = 0.
+   DO jj = 2, jpj-1
+      DO ji = 1, jpi-1
+         dx(ji,jj) = rt * rad * abs( lon(ji+1) - lon(ji) ) * COS( rad * lat(jj) )
+      END DO
+   END DO 
+   dx(  jpi,1:jpj) = dx(  jpi-1,1:jpj  )
+   dx(1:jpi,  jpj) = dx(1:jpi  ,  jpj-1)
+   dx(1:jpi,1    ) = dx(1:jpi  ,2      )
+   !++ 
+   DO jj = 1, jpj-1
+      DO ji = 1, jpi   
+         dy(ji,jj) = rt * rad * abs( lat(jj+1) - lat(jj) )
+      END DO
+   END DO   
+   dy(1:jpi,jpj) = dy(1:jpi,jpj-1) 
+
+
+   !!---------------------------------------------------------------------   
+   !! create output file
+   !!
+   status   = nf90_create( trim(out_file), NF90_WRITE, ncid )
+   status   = nf90_close ( ncid )      
+
+   CALL Write_Ncdf_dim ( 'lon'    , trim(out_file), jpi    )
+   CALL Write_Ncdf_dim ( 'lat'    , trim(out_file), jpj    )
+   CALL Write_Ncdf_dim ( 'lev'    , trim(out_file), jpka   )
+   CALL Write_Ncdf_dim ( 'nhym'   , trim(out_file), jpka   )  
+   CALL Write_Ncdf_dim ( 'nhyi'   , trim(out_file), jpka+1 )  
+   CALL Write_Ncdf_dim ( 'time'   , trim(out_file), 0      )    
+
+
+   !!---------------------------------------------------------------------   
+   !! Initialize the name of the dimensions for geostrophic winds in the output file
+   !!
+   dimnames(1) = 'lon'
+   dimnames(2) = 'lat'       
+   dimnames(3) = 'lev'
+   dimnames(4) = 'time'
+
+   CALL Write_Ncdf_var(  'lon',  'lon', trim(out_file), lon, 'double' ) 
+   CALL Write_Ncdf_var(  'lat',  'lat', trim(out_file), lat, 'double' ) 
+   CALL Write_Ncdf_var(  'lev',  'lev', trim(out_file), lev, 'double' )
+   CALL Write_Ncdf_var( 'hyai', 'nhyi', trim(out_file), A_w, 'double' ) 
+   CALL Write_Ncdf_var( 'hybi', 'nhyi', trim(out_file), B_w, 'double' ) 
+   CALL Write_Ncdf_var( 'hyam', 'nhym', trim(out_file), A_r, 'double' )
+   CALL Write_Ncdf_var( 'hybm', 'nhym', trim(out_file), B_r, 'double' )
+
+   CALL Write_Ncdf_var( 'tmask', dimnames(1:2), trim(out_file), tmask, 'float' )
+
+
+   !!---------------------------------------------------------------------
+   ! Read time variable
+   ALLOCATE(tmp1d (1:jptime))    
+   CALL Read_Ncdf_var ( 'time', trim(cn_dir)//'/'//trim(file_t), tmp1d )    
+   !!---------------------------------------------------------------------     
+
+
+   DO kt=1,jptime  
+   !
+      WRITE(stdout,*) '======================'
+      WRITE(stdout,*) 'time = ',kt,'/',jptime   
+      !
+      CALL Write_Ncdf_var( 'time', dimnames(4:4), trim(out_file), tmp1d(kt:kt), kt, 'double' )
+      !
+      IF( kt == 1 ) THEN
+         CALL Duplicate_lon_lat_time( trim(cn_dir)//'/'//trim(file_t), out_file )
+         CALL Duplicate_lev_hyb     ( trim(cn_dir)//'/'//trim(file_t), out_file ) 
+      ENDIF      
+      !
+      IF ( varnames(4) == "LNSP" ) THEN
+         CALL Read_Ncdf_var( varnames(4) , trim(cn_dir)//'/'//trim(file_p), slp(:,:),  kt, 1 )      !<-- log of surface pressure
+      ELSE
+         CALL Read_Ncdf_var( varnames(4) , trim(cn_dir)//'/'//trim(file_p), slp(:,:),  kt )
+      END IF
+      !
+      IF (ln_slp_log) THEN
+        DO jj = 1, jpj
+           DO ji = 1, jpi      
+             slp(ji,jj) = exp( slp(ji,jj) )
+           END DO
+        END DO
+      ENDIF
+      !
+      IF (ln_read_zsurf) THEN
+        CALL Read_Ncdf_var( varnames(1) , trim(cn_dir)//'/'//trim(file_z), zsurf(:,:),  kt, 1 )      !<-- surface geopotential 
+      ELSE
+        zsurf(:,:) = 0.
+      END IF
+      !    
+      CALL Read_Ncdf_var ( varnames(2), trim(cn_dir)//'/'//trim(file_t), temp(:,:,:,:), kt )      !<-- temperature   
+      CALL Read_Ncdf_var ( varnames(3), trim(cn_dir)//'/'//trim(file_q), humi(:,:,:,:), kt )      !<-- humidity        
+      WHERE(humi.LT.1.E-08) humi = 1.E-08 !<-- negative values in ECMWF
+      !
+      ! Smoothing of surface fields to remove gibbs oscillations (must be done on both fields or none of them)
+      !IF( ln_slp_smth ) CALL smooth_field( jpi, jpj, slp  (:,:), tmask(:,:), 3 )
+      !IF (ln_read_zsurf.AND.ln_slp_smth) CALL smooth_field( jpi, jpj, zsurf(:,:), tmask(:,:), 3 )
+      !IF( ln_slp_smth ) CALL DTV_Filter( jpi, jpj, slp(:,:), tmask(:,:,1), 25, kt )  !<-- not yet robust enough
+
+      !!GS: DO NOT USE SMOOTH + LAND MASK FOR NOW (BUG) 
+      IF( ln_slp_smth ) THEN
+        slp_smth(:,:) = slp(:,:)
+        wrke(:,:) = 1.
+        CALL smooth_field( jpi, jpj, slp_smth(:,1:jpj),  wrke(:,1:jpj), 3 )
+        IF (ln_read_zsurf) THEN
+          zsurf_smth(:,:) = zsurf(:,:)
+          CALL smooth_field( jpi, jpj, zsurf_smth(:,:), wrke(:,:), 3 )
+        END IF
+        IF( ABS(lat_smth).GT.0.1) THEN
+          DO jj = 1, jpj
+             DO ji = 1, jpi
+               IF ((lat(jj).GE.lat_smth).OR.(lat(jj).LE.-1.*lat_smth)) THEN
+                 slp(ji,jj) = slp_smth(ji,jj)
+                 IF (ln_read_zsurf) zsurf(ji,jj) = zsurf_smth(ji,jj)
+               END IF
+             END DO
+          END DO
+        ELSE
+          slp(:,:) = slp_smth(:,:)
+          IF (ln_read_zsurf) zsurf(:,:) = zsurf_smth(:,:)
+        END IF
+      END IF
+      CALL Write_Ncdf_var( 'slp', dimnames(1:2), trim(out_file), slp(:,:), kt, 'float' )
+      IF (ln_read_zsurf) CALL Write_Ncdf_var( 'zsurf', dimnames(1:2), trim(out_file), zsurf(:,:), kt, 'float' )
+
+      !
+      ! Compute the altitude at layer interfaces
+      ghw(:,:,1:jpka)  = 0.
+      ghw(:,:,  jpka)  = zsurf(:,:) * (1. / grav) ! * tmask(:,:)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            DO jk = 0, jpka
+               ph(jk) = A_w( jk ) + B_w( jk ) * slp( ji, jj )    !<-- Pa
+            END DO
+            !ph(0) = 0.1
+            IF ( nhym .EQ. jpka) ph(0) = 1.
+            DO jk = jpka,1,-1
+               tv      = temp( ji, jj, jk, 1 ) * ( 1. + zvir*humi( ji, jj, jk, 1 ) )  !<-- Virtual temperature
+               e3t ( ji, jj, jk   ) =  (1./grav)*( Rd * tv * log( ph( jk ) / ph( jk-1 ) ) ) !* tmask(ji, jj)
+               ghw ( ji, jj, jk-1 ) =  e3t( ji, jj, jk ) + ghw( ji, jj, jk ) 
+            END DO  
+         END DO
+      END DO
+
+      IF( ln_slp_smth ) THEN
+        wrke(:,:) = 1.
+        DO jk = 0, jpka
+          ghw_smth(:,:) = ghw(:,:,jk)
+          CALL smooth_field( jpi, jpj, ghw_smth(:,1:jpj),  wrke(:,1:jpj), 3 )
+          IF( ABS(lat_smth).GT.0.1) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                 IF ((lat(jj).GE.lat_smth).OR.(lat(jj).LE.-1.*lat_smth)) ghw(ji,jj,jk) = ghw_smth(ji,jj)
+               END DO
+            END DO
+          ELSE
+            ghw(:,:,jk) = ghw_smth(:,:)
+          END IF
+        END DO
+      END IF
+
+      !
+      ! Compute horizontal gradient of slp in x-direction (FX = dslp / dx)
+      FX(:,:) = 0.
+      DO jj = 1, jpj
+         DO ji = 1, jpi-1     
+            IF ((tmask(ji,jj) .gt. 0.5).AND.(tmask(ji+1,jj) .gt. 0.5)) THEN
+              cff            =  2. / (  dx( ji+1,jj ) +  dx( ji,jj ) )
+              FX( ji ,jj )   = cff * ( slp( ji+1,jj ) - slp( ji,jj ) )
+            ELSE
+              tmask2(ji:ji+1,jj) = 0.
+            END IF
+         END DO
+      END DO
+
+      IF (ln_perio_latbc) THEN
+         ! apply periodicity       
+         DO jj = 1, jpj
+            IF ((tmask(1,jj) .gt. 0.5).AND.(tmask(jpi,jj) .gt. 0.5)) THEN
+               cff              =  2. / (  dx( 1, jj) +  dx( jpi, jj) )
+               FX( jpi , jj)    = cff * ( slp( 1, jj) - slp( jpi, jj) )
+            ELSE
+               tmask2(   1, jj) = 0.
+               tmask2( jpi, jj) = 0.
+            END IF
+         END DO
+      ELSE
+         ! apply no-gradient
+         DO jj = 1, jpj
+            FX( jpi ,jj )   = FX( jpi-1 ,jj )
+         END DO
+      ENDIF    
+
+      !
+      ! Compute horizontal gradient of slp in y-direction (FE = dslp / dy)
+      FE(:,:) = 0.
+      DO jj = 1, jpj-1
+         DO ji = 1, jpi     
+            IF ((tmask(ji,jj) .gt. 0.5).AND.(tmask(ji,jj+1) .gt. 0.5)) THEN
+               cff            =  2. / (  dy( ji,jj+1 ) +  dy( ji,jj ) )
+               FE( ji ,jj )   = cff * ( slp( ji,jj+1 ) - slp( ji,jj ) )
+            ELSE
+               tmask2(ji,jj:jj+1) = 0.
+            END IF
+         END DO
+      END DO
+
+      ! apply no-gradient
+      DO ji = 1, jpi
+         FE( ji ,jpj ) = FE( ji ,jpj-1 )          
+      END DO              
+
+      !
+      !++ Compute the geostrophic winds
+      !
+      dZx(:,:)  = 0.
+      dZe(:,:)  = 0.
+      wrkX(:,:) = 0.
+      wrkE(:,:) = 0.
+      !////////////   
+      DO jk=1,jpka
+      !////////////           
+
+         !
+         ! Compute horizontal gradient of altitude in x-direction along the coordinate dZx = (dz / dx)s
+         DO jj = 1, jpj
+            DO ji = 1, jpi-1
+               IF ((tmask(ji,jj) .gt. 0.5).AND.(tmask(ji+1,jj) .gt. 0.5)) THEN
+                  cff          = 1. / (dx(ji+1,jj)+dx(ji,jj))
+                  dZx(ji,jj)   = cff * (  (ghw( ji+1, jj, jk-1) - ghw( ji, jj, jk-1))   &
+                     &                   +(ghw( ji+1, jj, jk  ) - ghw( ji, jj, jk  )) )
+               ELSE
+                 tmask2(ji:ji+1,jj) = 0.
+               END IF
+            END DO
+         END DO   
+
+         IF (ln_perio_latbc) THEN
+           ! apply periodicity
+           DO jj = 1, jpj
+              IF ((tmask(1,jj) .gt. 0.5).AND.(tmask(jpi,jj) .gt. 0.5)) THEN
+                 cff          = 1. / (dx(1,jj)+dx(jpi,jj)) 
+                 dZx(jpi,jj)  = cff * ( (ghw( 1, jj, jk-1) - ghw( jpi, jj, jk-1))   &
+                    &                  +(ghw( 1, jj, jk  ) - ghw( jpi, jj, jk  )) )
+                 dZx(  0,jj)  = dZx(jpi,jj)
+              ELSE
+                 tmask2(  1, jj) = 0.
+                 tmask2(jpi, jj) = 0.
+              END IF
+           END DO 
+         ELSE
+           ! apply no-gradient
+           DO jj = 1, jpj
+              dZx(jpi,jj)     = dZx(jpi-1,jj)
+              dZx(  0,jj)     = dZx(    1,jj)            
+           END DO         
+         END IF
+
+         !
+         ! Compute horizontal gradient of altitude in y-direction along the coordinate dZy = (dz / dy)s
+         DO jj = 1, jpj-1
+            DO ji = 1, jpi
+               IF ((tmask(ji,jj) .gt. 0.5).AND.(tmask(ji,jj+1) .gt. 0.5)) THEN
+                  cff          = 1. / (dy(ji,jj)+dy(ji,jj+1))
+                  dZe(ji,jj)   = cff * (  (ghw( ji, jj+1, jk-1) - ghw( ji, jj, jk-1))  &
+                     &                   +(ghw( ji, jj+1, jk  ) - ghw( ji, jj, jk  )) )   
+               ELSE
+                  tmask2(ji,jj:jj+1) = 0.
+               END IF
+            END DO
+         END DO 
+
+         ! apply no-gradient 
+         DO ji = 1, jpi
+            dZe(ji,jpj)     = dZe(ji,jpj-1)
+            dZe(ji,  0)     = dZe(ji,    1)            
+         END DO              
+
+         
+         ! Compute horizontal pressure gradient in x-direction along the coordinate wrkX = (dp/dx)s
+         DO jj = 1, jpj
+            DO ji = 2, jpi
+               IF ((tmask(ji,jj) .gt. 0.5).AND.(tmask(ji-1,jj) .gt. 0.5)) THEN
+                  cff = slp(ji,jj) * (B_w(jk)-B_w(jk-1)) + (A_w(jk)-A_w(jk-1))
+                  wrkX(ji,jj) = B_r(jk) * 0.5 * (FX(ji,jj)+FX(ji-1,jj))  &
+                     &               * (ghw(ji,jj,jk)-ghw(ji,jj,jk-1)) / cff
+               ELSE
+                  tmask2(ji-1:ji,jj) = 0.
+               END IF
+            END DO
+         END DO
+
+         IF (ln_perio_latbc) THEN
+           ! apply periodicity
+           ji = 1
+           DO jj = 1, jpj
+              IF ((tmask(1,jj) .gt. 0.5).AND.(tmask(jpi,jj) .gt. 0.5)) THEN
+                 cff = slp(ji,jj) * (B_w(jk)-B_w(jk-1)) + (A_w(jk)-A_w(jk-1))
+                 wrkX(ji,jj) = B_r(jk) * 0.5 * (FX(ji,jj)+FX(jpi,jj))  &
+                    &               * (ghw(ji,jj,jk)-ghw(ji,jj,jk-1)) / cff            
+              ELSE
+                 tmask2(  1,jj) = 0.
+                 tmask2(jpi,jj) = 0.
+              END IF
+           END DO
+         ELSE
+           ! apply no gradient  
+           DO jj = 1, jpj
+              wrkX(1,jj) = wrkX(2,jj)
+           END DO     
+         END IF 
+
+
+         ! Compute horizontal pressure gradient in y-direction along the coordinate wrkE = (dp/dy)s
+         DO jj = 2, jpj
+            DO ji = 1, jpi
+               IF ((tmask(ji,jj) .gt. 0.5).AND.(tmask(ji,jj-1) .gt. 0.5)) THEN
+                  cff = slp(ji,jj) * (B_w(jk)-B_w(jk-1)) + (A_w(jk)-A_w(jk-1))
+                  wrkE(ji,jj) = B_r(jk) * 0.5 * (FE(ji,jj)+FE(ji,jj-1))  &
+                     &               * (ghw(ji,jj,jk)-ghw(ji,jj,jk-1)) / cff
+               ELSE
+                  tmask2(ji,jj-1:jj) = 0.
+               END IF
+            END DO
+         END DO    
+
+         ! apply no gradient  
+         jj = 1
+         DO ji = 1, jpi       
+            wrkE(ji,1) = wrkE(ji,2) 
+         END DO   
+
+
+         !+++ Finalize pressure gradient/geostrophic wind computation
+         IF(ln_hpg_frc) THEN
+            DO jj=1,jpj
+               DO ji=1,jpi                     
+                 IF (tmask2(ji,jj).GT.0.5) THEN
+                   uhpg(ji,jj,jk,1) = - grav*( wrkX(ji,jj) - 0.5*( dZx(ji,jj)+dZx(ji-1,jj  ) ) ) 
+                   if (lat(1).GT.0.) vhpg(ji,jj,jk,1) =   grav*( wrkE(ji,jj) - 0.5*( dZe(ji,jj)+dZe(ji  ,jj-1) ) )
+                   if (lat(1).LT.0.) vhpg(ji,jj,jk,1) = - grav*( wrkE(ji,jj) - 0.5*( dZe(ji,jj)+dZe(ji  ,jj-1) ) )
+                 ELSE
+                   uhpg(ji,jj,jk,1) = 0.
+                   vhpg(ji,jj,jk,1) = 0.
+                 END IF
+               END DO   
+            END DO   
+         ELSE
+            DO jj=1,jpj               
+               DO ji=1,jpi    
+                 IF (tmask2(ji,jj).GT.0.5) THEN
+                   cff = grav / ff_t(ji,jj)
+                   ! geostrophic wind computed only where Coriolis .ge. 3.e-5 (~12deg)
+                   if(abs(ff_t(ji,jj)) < 2.5e-5) cff = 0.    
+                   ! minus sign for uhpg because y-derivatives are inverted
+                   vhpg(ji,jj,jk,1) = - cff*( wrkX(ji,jj) - 0.5*( dZx(ji,jj)+dZx(ji-1,jj  ) ) ) 
+                   if (lat(1).GT.0.) uhpg(ji,jj,jk,1) = - cff*( wrkE(ji,jj) - 0.5*( dZe(ji,jj)+dZe(ji  ,jj-1) ) )
+                   if (lat(1).LT.0.) uhpg(ji,jj,jk,1) =   cff*( wrkE(ji,jj) - 0.5*( dZe(ji,jj)+dZe(ji  ,jj-1) ) )
+                 ELSE
+                   uhpg(ji,jj,jk,1) = 0.
+                   vhpg(ji,jj,jk,1) = 0.
+                 END IF
+               END DO   
+            END DO          
+         ENDIF
+
+      !////////////        
+      END DO ! jk
+      !////////////        
+
+      IF (ln_geo_wnd) THEN
+        CALL Write_Ncdf_var ( 'ugeo', dimnames(1:4), trim(out_file), uhpg, kt, 'float' )
+        CALL Write_Ncdf_var ( 'vgeo', dimnames(1:4), trim(out_file), vhpg, kt, 'float' )
+      ELSE
+        CALL Write_Ncdf_var ( 'uhpg', dimnames(1:4), trim(out_file), uhpg, kt, 'float' )
+        CALL Write_Ncdf_var ( 'vhpg', dimnames(1:4), trim(out_file), vhpg, kt, 'float' )
+      END IF
+
+      IF (kt .EQ. 1) THEN
+        tmask2(:,1)   = 0. ! force northern v line drowning
+        tmask2(:,jpj) = 0. ! force northern v line drowning
+        CALL Write_Ncdf_var( 'tmask', dimnames(1:2), trim(out_file), tmask2, 'float' )
+      END IF
+
+   END DO ! kt
+   !
+   DEALLOCATE( zsurf, zsurf_smth, slp, slp_smth, temp,humi,uhpg,vhpg)     
+   IF (jpka.NE.nhym) DEALLOCATE(tmp_fullw,tmp_fullm)
+   !
+   STOP
+   !
+END PROGRAM main
Index: /utils/tools/ABL_TOOLS/main_vinterp.F90
===================================================================
--- /utils/tools/ABL_TOOLS/main_vinterp.F90	(revision 11589)
+++ /utils/tools/ABL_TOOLS/main_vinterp.F90	(revision 11589)
@@ -0,0 +1,590 @@
+PROGRAM main
+   !!======================================================================
+   !!                     ***  PROGRAM main  ***
+   !!
+   !! ** Purpose : Vertical interpolation of ECMWF dataset on a given fixed 
+   !!              vertical grid
+   !! 
+   !!======================================================================
+   !! History : 2016-10  (F. Lemarié)  Original code
+   !!   
+   !!----------------------------------------------------------------------
+   USE module_io       ! I/O routines
+   USE module_interp   ! vertical interpolation routines
+   USE module_grid     ! compute input and output grids 
+   !!
+   IMPLICIT NONE
+   !!----------------------------------------------------------------------
+   !! 
+   !!  
+   !! 
+   !!----------------------------------------------------------------------
+   !
+   INTEGER                                  :: ji,jj,jk,kt, jk_in, nhym, nhyi
+   INTEGER                                  :: jpka_in, jpka      ! number of vertical levels for input and target grids 
+   INTEGER                                  :: jpi , jpj          ! number of grid points in x and y directions     
+   INTEGER                                  :: iloc, jloc         ! grid indexes for c1d case
+   INTEGER                                  :: status             
+   INTEGER                                  :: jptime,ctrl
+   INTEGER                                  :: ioerr
+   INTEGER, ALLOCATABLE, DIMENSION(:,:  )   :: ind  
+   INTEGER, PARAMETER                       :: stdout  = 6
+   INTEGER, PARAMETER                       :: jp_weno = 1
+   INTEGER, PARAMETER                       :: jp_spln = 2   
+   !!
+   REAL(8)                                  :: hc,hmax,theta_s,z1 ! parameters related to the target vertical grid
+   REAL(8)                                  :: cff
+   !!  
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: A_w                ! A coefficients to reconstruct ECMWF grid
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: A_wa               ! A coefficients to reconstruct ECMWF grid
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: B_w                ! B coefficients to reconstruct ECMWF grid
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: B_wa               ! B coefficients to reconstruct ECMWF grid
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: tmp1d, tmp_fullw, tmp_fullm ! temporary/working 1D arrays
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: e3t,e3w            ! thickness of vertical layers in target grid
+   REAL(8), ALLOCATABLE, DIMENSION(:      ) :: ght,ghw            ! altitude of vertical grid points 
+   REAL(8), ALLOCATABLE, DIMENSION(:,:    ) :: e3_bak
+   REAL(8), ALLOCATABLE, DIMENSION(:,:,:  ) :: ghw_in             ! altitude of cell interfaces  of ECMWF grid
+   REAL(8), ALLOCATABLE, DIMENSION(:,:,:  ) :: e3t_in             ! thickness of vertical layers in ECMWF grid
+   REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: humi
+   REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: tair, tpot
+   REAL(8), ALLOCATABLE, DIMENSION(:,:,:,:) :: varout, varc1d 
+   REAL(8), ALLOCATABLE, DIMENSION(:,:    ) :: slp, zsurf
+   REAL(8), ALLOCATABLE, DIMENSION(:,:    ) :: tmask, tmask2      ! land-sea mask  
+   !!
+   CHARACTER(len=500)                       :: file_u,file_v,file_hpg,file_geos  ! ECMWF files containing wind components 
+   CHARACTER(len=500)                       :: file_t,file_q, file_m             ! ECMWF files containing tracers and mask 
+   CHARACTER(len=500)                       :: file_z,file_p,cn_dir,file_in      ! ECMWF files containing surface geopot and pressure
+   CHARACTER(len=500)                       :: grd_file, abl_file, drwn_file, out_file
+   CHARACTER(len=500)                       :: namelistf,stmp
+   CHARACTER(len=500)                       :: argument, var_file
+   CHARACTER(len= 20),DIMENSION(4)          :: dimnames
+   CHARACTER(len= 20),DIMENSION(11)         :: varnames, outnames
+   CHARACTER(len=500),DIMENSION(11)         :: filnames
+   CHARACTER(6)                             :: mask_var           ! name of mask variable in file_m file            
+   CHARACTER(6)                             :: var_name
+   !!
+   LOGICAL                                  :: ln_read_zsurf      ! read surface geopotential or not
+   LOGICAL                                  :: ln_read_mask       ! read land-sea mask or not
+   LOGICAL                                  :: ln_perio_latbc     ! use periodic BC along the domain latitudinal edges (for global data) or use zero-gradient BC (for regional data)
+   LOGICAL                                  :: ln_c1d             ! output only a single column in output file
+   LOGICAL                                  :: ln_hpg_frc         ! compute horizontal pressure gradient
+   LOGICAL                                  :: ln_geo_wnd         ! compute goestrophic wind components 
+   LOGICAL                                  :: ln_slp_smth        ! apply gibbs oscillation filetring on mean sea level pressure
+   LOGICAL                                  :: ln_drw_smth        ! apply gibbs oscillation filetring on mean sea level pressure
+   LOGICAL                                  :: ln_slp_log         ! log(sea-level pressure) or sea-level pressure
+   LOGICAL                                  :: ln_lsm_land        ! if T mask is 1 over land and 0 over ocean if F it is the other way around
+   LOGICAL                                  :: ln_impose_z1       ! impose the altitude of the first level in target grid
+   INTEGER                                  :: ptemp_method       ! way to compute potential temperature 
+                                                                  ! = 0  (absolute temperature)
+                                                                  ! = 1  (potential temperature with local ref pressure)
+                                                                  ! = 2  (potential temperature with global ref pressure on temperature perturbation)
+                                                                  ! = 3  (potential temperature with global ref pressure)
+   !!   
+   REAL(8), PARAMETER     :: grav  =   9.80665
+
+   !!---------------------------------------------------------------------
+   !! List of variables read in the namelist file 
+   NAMELIST/nml_dom/    jpka, hmax, theta_s, hc, ln_impose_z1, z1   
+   NAMELIST/nml_opt/    ptemp_method, ln_slp_log, ln_slp_smth, ln_read_mask, ln_perio_latbc, &
+      &                 ln_hpg_frc, ln_geo_wnd, ln_c1d, ln_read_zsurf, ln_lsm_land, ln_drw_smth
+   NAMELIST/nml_fld/    cn_dir, file_u, file_v, file_t,               &
+      &                 file_q, file_z, file_p, file_hpg, file_geos,  &
+      &                 file_m, mask_var
+   NAMELIST/nml_out/    grd_file, abl_file, drwn_file, var_name
+   NAMELIST/nml_c1d/    iloc, jloc
+   !!
+   !! get the namelist file name
+   CALL get_command_argument( 1, argument, ctrl, status)
+   !
+   SELECT CASE(status)
+   CASE(0)
+      namelistf = trim(argument)
+   CASE(-1)
+      WRITE(stdout,*) "### Error: file name too long"
+      STOP
+   CASE DEFAULT
+      namelistf = 'namelist_abl_tools'
+   END SELECT
+   !!---------------------------------------------------------------------
+
+
+   !!---------------------------------------------------------------------
+   !! read namelist variables
+   ctrl = 0
+   OPEN(50, file=namelistf, status='old', form='formatted', access='sequential', iostat=ioerr)
+   IF (ioerr /= 0) ctrl = ctrl + 1   
+   READ(50,nml_dom, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1   
+   READ(50,nml_opt, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
+   READ(50,nml_fld, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
+   READ(50,nml_out, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1  
+   IF( ln_c1d ) THEN
+      print*,'c1d is activated'
+      READ(50,nml_c1d, iostat=ioerr); IF (ioerr /= 0) ctrl = ctrl + 1
+   ENDIF
+
+   IF (ctrl > 0) THEN
+      WRITE(stdout,*) "### E R R O R while reading namelist file '",trim(namelistf),"'"
+      WRITE(stdout,*) " ctrl = ",ctrl
+      STOP
+   ELSE 
+      WRITE(stdout,*) " Namelist file ",trim(namelistf)," OK "    
+   ENDIF
+   IF( ln_hpg_frc .AND. ln_geo_wnd ) THEN
+      WRITE(stdout,*) "### E R R O R conflicting options "
+      WRITE(stdout,*) "ln_hpg_frc and ln_geo_wnd can not both be set to True"
+      STOP   
+   ENDIF
+   
+   SELECT CASE (ptemp_method) 
+   CASE(0)
+      WRITE(stdout,*) "Absolute temperature option is activated"
+   CASE(1)   
+      WRITE(stdout,*) "Potential temperature option with  local reference pressure is activated"   
+   CASE(2)    
+      WRITE(stdout,*) "Potential temperature option with global reference pressure on temperature perturbation is activated"   
+   CASE(3)    
+      WRITE(stdout,*) "Potential temperature option with global reference pressure is activated"  
+   END SELECT
+   
+   IF(ln_slp_smth) WRITE(stdout,*) "MSLP smoothing option is activated" 
+   IF(ln_hpg_frc ) WRITE(stdout,*) "Large-scale pressure gradient will be interpolated"      
+   IF(ln_geo_wnd ) WRITE(stdout,*) "Geostrophic winds will be interpolated" 
+   !!---------------------------------------------------------------------   
+       
+   !!-------------------------------------------------------------------------------------   
+   !! list of variables to treat
+   !!
+   !! get the variable name
+   CALL get_command_argument( 2, argument, ctrl, status)
+   SELECT CASE(status)
+   CASE(0)
+      var_name = trim(argument)
+   CASE(-1)
+      WRITE(stdout,*) "### Error: file name too long"
+      STOP
+   END SELECT
+   WRITE(stdout,*) "var_name: ", trim(var_name), "lenght: ", Len_Trim(var_name)
+   !!
+   IF( ln_hpg_frc ) THEN 
+      file_in = trim(file_hpg )
+   ELSE
+      file_in = trim(file_geos)
+   ENDIF
+   varnames = [character(len=4  ) :: 'Z'    ,    'T',    'Q',    'U',    'V', 'MSL', 'LSM', 'uhpg', 'vhpg', 'ugeo', 'vgeo' ]
+   outnames = [character(len=5  ) :: 'zsurf', 'tair', 'humi', 'uwnd', 'vwnd', 'slp',    '', 'uhpg', 'vhpg', 'ugeo', 'vgeo' ]
+   filnames = [character(len=500) ::  file_z, file_t, file_q, file_u, file_v, file_p, file_m, file_in, file_in, file_in, file_in ]
+   IF( ln_slp_log ) varnames(6) = 'LNSP'
+ 
+   !!---------------------------------------------------------------------
+   ! check files content
+   ctrl = 0
+   CALL Read_Ncdf_dim('lev',trim(cn_dir)//'/'//trim(file_t),jpka_in)   
+   !
+   IF (ln_read_zsurf) THEN
+     IF( .not. VAR_EXISTENCE( trim(varnames(1)) , trim(cn_dir)//'/'//trim(file_z) )      &
+        & .or. jpka_in == 1 ) ctrl = ctrl + 1 
+   ENDIF
+   WRITE(stdout,*) trim(varnames(1)), ctrl
+   IF ( .not. VAR_EXISTENCE( trim(varnames(2)) , trim(cn_dir)//'/'//trim(file_t) ) ) ctrl = ctrl + 1 
+   WRITE(stdout,*) trim(varnames(2)), ctrl
+   IF ( .not. VAR_EXISTENCE( trim(varnames(3)) , trim(cn_dir)//'/'//trim(file_q) ) ) ctrl = ctrl + 1   
+   WRITE(stdout,*) trim(varnames(3)), ctrl
+   IF ( .not. VAR_EXISTENCE( trim(varnames(4)) , trim(cn_dir)//'/'//trim(file_u) ) ) ctrl = ctrl + 1
+   WRITE(stdout,*) trim(varnames(4)), ctrl
+   IF ( .not. VAR_EXISTENCE( trim(varnames(5)) , trim(cn_dir)//'/'//trim(file_v) ) ) ctrl = ctrl + 1 
+   WRITE(stdout,*) trim(varnames(5)), ctrl
+   IF ( .not. VAR_EXISTENCE( trim(varnames(6)) , trim(cn_dir)//'/'//trim(file_p) ) ) ctrl = ctrl + 1 
+   WRITE(stdout,*) trim(varnames(6)), ctrl
+   IF (ln_read_mask) THEN
+     IF ( .not. VAR_EXISTENCE( trim(varnames(7)) , trim(cn_dir)//'/'//trim(file_m) ) ) ctrl = ctrl + 1 
+     WRITE(stdout,*) trim(varnames(7)), ctrl
+     varnames(7) = TRIM(mask_var)
+   ENDIF
+   IF (ln_hpg_frc) THEN
+      IF ( .not. VAR_EXISTENCE( trim(varnames(8)) , trim(cn_dir)//'/'//trim(file_hpg) ) ) ctrl = ctrl + 1 
+      WRITE(stdout,*) trim(varnames(8)), ctrl
+      IF ( .not. VAR_EXISTENCE( trim(varnames(9)) , trim(cn_dir)//'/'//trim(file_hpg) ) ) ctrl = ctrl + 1    
+      WRITE(stdout,*) trim(varnames(9)), ctrl
+   ENDIF
+   IF (ln_geo_wnd) THEN
+      IF ( .not. VAR_EXISTENCE( trim(varnames(10)) , trim(cn_dir)//'/'//trim(file_geos) ) ) ctrl = ctrl + 1 
+      WRITE(stdout,*) trim(varnames(10)), ctrl
+      IF ( .not. VAR_EXISTENCE( trim(varnames(11)) , trim(cn_dir)//'/'//trim(file_geos) ) ) ctrl = ctrl + 1    
+      WRITE(stdout,*) trim(varnames(11)), ctrl
+   ENDIF
+          
+   IF ( ctrl > 0 ) THEN  
+      WRITE(stdout,*) "### E R R O R while reading ECMWF atmospheric files "
+      STOP
+   ELSE
+      WRITE(stdout,*) " ECMWF atmospheric files OK "        
+   ENDIF
+   !!---------------------------------------------------------------------
+
+
+   !!---------------------------------------------------------------------
+   !! read the dimensions for the input files
+   CALL Read_Ncdf_dim ( 'time', trim(cn_dir)//'/'//trim(file_t), jptime  ) 
+   CALL Read_Ncdf_dim ( 'lon' , trim(cn_dir)//'/'//trim(file_t), jpi     )     
+   CALL Read_Ncdf_dim ( 'lat' , trim(cn_dir)//'/'//trim(file_t), jpj     )     
+   CALL Read_Ncdf_dim ( 'nhym' , trim(cn_dir)//'/'//trim(file_t), nhym   )
+   CALL Read_Ncdf_dim ( 'nhyi' , trim(cn_dir)//'/'//trim(file_t), nhyi   )
+   !WRITE(stdout,*) "jpka_in, jptime, jpi, jpj, nhym, nhyi: ", jpka_in, jptime, jpi, jpj, nhym, nhyi
+   !
+   !!---------------------------------------------------------------------
+
+
+   !!---------------------------------------------------------------------
+   !! allocate arrays  
+   ALLOCATE( A_w    (               0:jpka_in) )         
+   ALLOCATE( B_w    (               0:jpka_in) )
+   ALLOCATE( e3t_in ( 1:jpi, 1:jpj, 1:jpka_in) )
+   ALLOCATE( ghw_in ( 1:jpi, 1:jpj, 0:jpka_in) ) 
+   ALLOCATE( slp    ( 1:jpi, 1:jpj              ) )
+   ALLOCATE( zsurf  ( 1:jpi, 1:jpj              ) )
+   ALLOCATE( tair   ( 1:jpi, 1:jpj, 1:jpka_in, 1) )
+   ALLOCATE( tpot   ( 1:jpi, 1:jpj, 1:jpka_in, 1) )
+   ALLOCATE( humi   ( 1:jpi, 1:jpj, 1:jpka_in, 1) )       
+   ALLOCATE( varout ( 1:jpi, 1:jpj, 1:jpka+1 , 1) )
+   ALLOCATE( ind    ( 1:jpi, 1:jpj              ) )
+   ALLOCATE( e3_bak ( 1:jpi, 1:jpj              ) )  
+   ALLOCATE( ght    (               1:jpka+1    ) )
+   ALLOCATE( ghw    (               1:jpka+1    ) )   
+   ALLOCATE( e3t    (               1:jpka+1    ) )
+   ALLOCATE( e3w    (               1:jpka+1    ) )      
+   ALLOCATE( tmask  ( 1:jpi, 1:jpj              ) )
+   ALLOCATE( tmask2 ( 1:jpi, 1:jpj              ) )
+   IF( ln_c1d ) ALLOCATE( varc1d( 1:3, 1:3, 1:jpka+1 , 1 ) )
+   IF( ln_c1d ) varc1d( 1:3, 1:3, 1:jpka+1 , 1 ) = 0.
+   IF (jpka_in.NE.nhym) THEN
+     ALLOCATE( tmp_fullw(1:nhyi) )
+     ALLOCATE( tmp_fullm(1:nhym) )
+   ENDIF
+   !
+   varout(:,:,:,1)    = 0.
+      
+   !!---------------------------------------------------------------------
+   !! Read the mask and remove some closed seas 
+   IF (ln_read_mask) THEN
+     CALL init_atm_mask(jpi,jpj,trim(cn_dir)//'/'//trim(file_m),trim(mask_var),ln_lsm_land,tmask)
+     IF (ln_geo_wnd) CALL Read_Ncdf_var ( 'tmask' , trim(cn_dir)//'/'//trim(file_geos), tmask2(:,:) )
+     IF (ln_hpg_frc) CALL Read_Ncdf_var ( 'tmask' , trim(cn_dir)//'/'//trim(file_hpg) , tmask2(:,:) )
+   ELSE
+     tmask(:,:)  = 1.
+     tmask2(:,:) = 1.
+   ENDIF
+   !!   
+   
+   !!---------------------------------------------------------------------
+   !! Compute the altitude and layer thickness of the target grid
+   CALL init_target_grid ( jpka, ght, ghw, e3t, e3w, hmax, hc, theta_s,   &
+      &                                              ln_impose_z1, z1  )
+
+   !! Write the grid file for the target grid
+   CALL Write_Grid_File  ( jpka, ght, ghw, e3t, e3w, trim(cn_dir)//'/'//trim(grd_file) )
+
+   
+   !! Read the static A and B coefficients for the ECMWF vertical grid
+   IF (jpka_in.EQ.nhym) THEN
+     CALL Read_Ncdf_var    ( 'hyai', trim(cn_dir)//'/'//trim(file_t), A_w )
+     CALL Read_Ncdf_var    ( 'hybi', trim(cn_dir)//'/'//trim(file_t), B_w )   
+   ELSE
+     CALL Read_Ncdf_var    ( 'hyai', trim(cn_dir)//'/'//trim(file_t), tmp_fullw )
+     A_w(0:jpka_in) = tmp_fullw(nhyi-(jpka_in+1)+1:nhyi)
+     CALL Read_Ncdf_var    ( 'hybi', trim(cn_dir)//'/'//trim(file_t), tmp_fullw )
+     B_w(0:jpka_in) = tmp_fullw(nhyi-(jpka_in+1)+1:nhyi)
+   ENDIF
+
+
+   !!---------------------------------------------------------------------   
+   !! create output file
+   !!
+   IF (Len_Trim(var_name) == 0) THEN
+      out_file = trim(cn_dir)//'/'//trim(abl_file)
+   ELSE
+      out_file = trim(cn_dir)//'/'//trim(var_name)//'_'//trim(abl_file)
+   ENDIF
+
+   IF(ln_c1d) THEN
+      CALL Init_output_File_c1d ( jpi, jpj, jpka, trim(cn_dir)//'/'//trim(file_t), out_file, tmask(:,:), iloc, jloc )
+   ELSE 
+      CALL Init_output_File ( jpi, jpj, jpka, trim(cn_dir)//'/'//trim(file_t), out_file, tmask(:,:) )
+   ENDIF
+
+   !!---------------------------------------------------------------------   
+   !! Initialize the name of the dimensions for the result of the interpolation
+   !!
+   dimnames(1) = 'lon'
+   dimnames(2) = 'lat'       
+   dimnames(3) = 'jpka'
+   dimnames(4) = 'time'
+   CALL Write_Ncdf_var( 'tmask', dimnames(1:2), trim(out_file), tmask2, 'float' )
+
+   !!---------------------------------------------------------------------
+   ! Read time variable
+   ALLOCATE(tmp1d (1:jptime))    
+   CALL Read_Ncdf_var ( 'time', trim(cn_dir)//'/'//trim(file_t), tmp1d )    
+   !!---------------------------------------------------------------------     
+ 
+   DO kt=1,jptime
+   !
+      WRITE(stdout,*) '======================'
+      WRITE(stdout,*) 'time = ',kt,'/',jptime   
+      !
+      CALL Write_Ncdf_var( 'time', dimnames(4:4), trim(out_file), tmp1d(kt:kt), kt, 'double' )
+      !
+      IF( kt == 1 ) THEN
+         CALL Duplicate_lon_lat_time( trim(cn_dir)//'/'//trim(file_t), out_file ) 
+         CALL add_globatt_real( out_file, "jpka"   , REAL(jpka) )
+         CALL add_globatt_real( out_file, "hmax"   , hmax       )
+         CALL add_globatt_real( out_file, "theta_s", theta_s    )
+         CALL add_globatt_real( out_file, "hc"     , hc         )
+         IF (ln_impose_z1) CALL add_globatt_real( out_file, "z1", z1 )
+      ENDIF
+      !
+      ! read SLP
+      !CALL Read_Ncdf_var ( varnames(6) , trim(cn_dir)//'/'//trim(file_p), slp,  kt )      !<-- (log of) surface pressure
+      !IF (ln_slp_log) THEN
+      !  DO jj = 1, jpj
+      !     DO ji = 1, jpi      
+      !       slp( ji ,jj )   = exp( slp(ji ,jj ) )
+      !     END DO
+      !  END DO
+      !ENDIF
+      CALL Read_Ncdf_var ( outnames(6) , trim(cn_dir)//'/'//trim(file_in), slp,  kt )      !<-- (log of) surface pressure
+      !
+      ! read ZSURF
+      IF (ln_read_zsurf) THEN
+        !CALL Read_Ncdf_var ( varnames(1) , trim(cn_dir)//'/'//trim(file_z), zsurf,  kt )      !<-- surface geopotential 
+        CALL Read_Ncdf_var ( outnames(1) , trim(cn_dir)//'/'//trim(file_in), zsurf,  kt )      !<-- surface geopotential 
+      ELSE
+        zsurf(:,:) = 0.
+      ENDIF
+      !    
+      ! Smoothing of SLP and ZSURF to remove gibbs oscillations (must be done on both fields or none of them)
+      !IF( ln_slp_smth ) CALL smooth_field( jpi, jpj, slp(:,:), tmask(:,:), 3 )
+      !IF (ln_read_zsurf.AND.ln_slp_smth) CALL smooth_field( jpi, jpj, zsurf(:,:), tmask(:,:), 3 )
+      !IF( ln_slp_smth ) CALL DTV_Filter( jpi, jpj, slp(:,:), tmask(:,:), 25, kt )  !<-- not yet robust enough
+      !
+      ! read tair and HUMI
+      CALL Read_Ncdf_var ( varnames(2) , trim(cn_dir)//'/'//trim(file_t), tair, kt    )      !<-- temperature   
+      CALL Read_Ncdf_var ( varnames(3) , trim(cn_dir)//'/'//trim(file_q), humi, kt    )      !<-- humidity        
+      WHERE(humi.LT.1.E-08) varout = 1.E-08 !<-- negative values in ECMWF
+      !
+      ! Reconstruct the ERA-Interim vertical grid in terms of altitude 
+      ghw_in(:,:,1:jpka_in)  = 0.
+      ghw_in(:,:,  jpka_in)  = zsurf(:,:) * (1. / grav)
+      CALL get_atm_grid( jpi, jpj, jpka_in, slp, tair,  &
+         &                          humi, A_w, B_w, e3t_in, ghw_in )
+      !
+      ! Compute potential temperature
+      tpot = tair ! save tpot
+      CALL get_pot_temp    ( jpi, jpj, jpka_in, slp, tpot, A_w, B_w,   &
+         &                   tmask(:,:), ptemp_method, humi(:,:,jpka_in,1), 0.5*ghw_in(:,:,jpka_in-1) )
+
+      ! Flip the vertical axis to go from k=0 at the bottom to k=N_in at the top of the atmosphere
+      CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, e3t_in          )       
+      CALL flip_vert_dim ( 0, jpka_in, jpi, jpj, ghw_in          ) 
+      CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, tpot(:,:,:,1)   )       
+      CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, humi(:,:,:,1)   ) 
+      
+      ! Correct the layer thickness to match hmax  
+      DO jj = 1, jpj
+         DO ji = 1, jpi          
+            cff = 0.
+            DO jk=1,jpka_in
+               cff = cff + e3t_in( ji, jj, jk )
+               IF ( cff > hmax ) THEN
+                  jk_in = jk
+                  EXIT
+               ENDIF 
+            END DO  
+            ind    ( ji, jj       ) = jk_in    
+            e3_bak ( ji, jj       ) = e3t_in ( ji, jj, jk_in )  ! store the value of the original layer thickness 
+            e3t_in ( ji, jj, jk_in) = e3t_in ( ji, jj, jk_in ) - ( cff - hmax )
+         END DO  
+      END DO 
+      !
+      IF (Len_Trim(var_name) == 0) THEN
+        !/
+        ! Interpolation of potential temperature TPOT
+        CALL zinterp ( jpi, jpj, jpka, jpka_in, ind,  &
+            &                         tpot, e3t_in, e3_bak, e3t, varout, jp_weno )  
+
+        varout(:,:,1,1) = varout(:,:,2,1)
+
+        IF (ln_c1d) THEN
+           DO jj=1,3
+              DO ji=1,3
+                 varc1d( ji, jj , 1:jpka+1 , 1 ) = varout(iloc,jloc, 1:jpka+1 , 1 )
+              END DO
+           END DO
+           CALL Write_Ncdf_var ('tair', dimnames(1:4), out_file, varc1d, kt, 'float' )
+        ELSE
+           CALL Write_Ncdf_var ('tair', dimnames(1:4), out_file, varout, kt, 'float' ) 
+        ENDIF
+        !
+        ! Interpolation of HUMI
+        CALL zinterp ( jpi, jpj, jpka, jpka_in, ind,  &
+            &                        humi, e3t_in, e3_bak, e3t, varout, jp_weno  )
+
+        !FL: dirty loop (possible issue in boundary conditions for WENO scheme)
+        DO jj = 1, jpj
+           DO ji = 1, jpi
+              DO jk = 2, jpka+1
+                 varout(ji,jj,jk,1) = MAX(varout(ji,jj,jk,1),1.E-08)                       !<-- negative values in ECMWF
+              END DO
+           END DO
+        END DO
+
+        varout(:,:,1,1) = varout(:,:,2,1) 
+
+        IF (ln_c1d) THEN
+           DO jj=1,3
+              DO ji=1,3
+                 varc1d( ji, jj , 1:jpka+1 , 1 ) = varout(iloc,jloc, 1:jpka+1 , 1 )
+              END DO
+           END DO
+           CALL Write_Ncdf_var  ('humi', dimnames(1:4), out_file, varc1d, kt, 'float' )
+        ELSE
+           CALL Write_Ncdf_var ('humi', dimnames(1:4), out_file, varout, kt, 'float' )  
+        ENDIF
+
+        !
+        ! Interpolate large-scale HPG or geostrophic wind
+        !
+        ! Read HPG
+        IF (ln_hpg_frc) THEN
+           CALL Read_Ncdf_var ( varnames(8) , trim(cn_dir)//'/'//trim(file_in), tair, kt )        
+           CALL Read_Ncdf_var ( varnames(9) , trim(cn_dir)//'/'//trim(file_in), humi, kt ) 
+        ENDIF
+        ! Read geostrophic wind
+        IF (ln_geo_wnd) THEN
+           CALL Read_Ncdf_var ( varnames(10) , trim(cn_dir)//'/'//trim(file_in), tair, kt )        
+           CALL Read_Ncdf_var ( varnames(11) , trim(cn_dir)//'/'//trim(file_in), humi, kt ) 
+        ENDIF
+        !
+        CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, tair( :,:,:,1 ) )       
+        CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, humi( :,:,:,1 ) ) 
+        !      
+        ! Interpolation of geostrophic U
+        CALL zinterp ( jpi, jpj, jpka, jpka_in, ind,  &
+            &                         tair, e3t_in, e3_bak, e3t, varout, jp_spln  )      
+        varout(:,:,1,1) = varout(:,:,2,1)
+           
+           IF (ln_c1d) THEN
+              DO jj=1,3
+                 DO ji=1,3
+                    varc1d( ji, jj , 1:jpka+1 , 1 ) = varout(iloc,jloc, 1:jpka+1 , 1 )
+                 END DO
+              END DO
+              CALL Write_Ncdf_var ( 'uhpg', dimnames(1:4), out_file, varc1d, kt, 'float' )
+           ELSE
+              IF (ln_hpg_frc) CALL Write_Ncdf_var ( varnames( 8), dimnames(1:4), out_file, varout, kt, 'float' )
+              IF (ln_geo_wnd) CALL Write_Ncdf_var ( varnames(10), dimnames(1:4), out_file, varout, kt, 'float' )
+           ENDIF
+           !
+           ! Interpolation of geostrophic V
+           CALL zinterp ( jpi, jpj, jpka, jpka_in, ind,  &
+               &                         humi, e3t_in, e3_bak, e3t, varout, jp_spln  )       
+           varout(:,:,1,1) = varout(:,:,2,1)
+
+           IF (ln_c1d) THEN
+              DO jj=1,3
+                 DO ji=1,3
+                    varc1d( ji, jj , 1:jpka+1 , 1 ) = varout(iloc,jloc, 1:jpka+1 , 1 )
+                 END DO
+              END DO
+              CALL Write_Ncdf_var ( 'vhpg', dimnames(1:4), out_file, varc1d, kt, 'float' )
+           ELSE
+              IF (ln_hpg_frc) CALL Write_Ncdf_var ( varnames( 9), dimnames(1:4), out_file, varout, kt, 'float' )
+              IF (ln_geo_wnd) CALL Write_Ncdf_var ( varnames(11), dimnames(1:4), out_file, varout, kt, 'float' )
+           ENDIF
+        !
+        ! Interpolation of total winds
+        !      
+        ! Read wind
+        CALL Read_Ncdf_var ( varnames(4) , trim(cn_dir)//'/'//trim(file_u), tair, kt )        
+        CALL Read_Ncdf_var ( varnames(5) , trim(cn_dir)//'/'//trim(file_v), humi, kt ) 
+        CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, tair( :,:,:,1 ) )       
+        CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, humi( :,:,:,1 ) ) 
+        !      
+        ! Interpolation of total U
+        CALL zinterp ( jpi, jpj, jpka, jpka_in, ind,  &
+            &                         tair, e3t_in, e3_bak, e3t, varout, jp_spln  )      
+        varout(:,:,1,1) = varout(:,:,2,1)
+
+        IF(ln_c1d) THEN
+           DO jj=1,3
+              DO ji=1,3
+                 varc1d( ji, jj , 1:jpka+1 , 1 ) = varout(iloc,jloc, 1:jpka+1 , 1 )
+              END DO
+           END DO
+           CALL Write_Ncdf_var ( 'uwnd', dimnames(1:4), out_file, varc1d, kt, 'float' )
+        ELSE
+           CALL Write_Ncdf_var ( 'uwnd', dimnames(1:4), out_file, varout, kt, 'float' )
+        ENDIF
+        !
+        ! Interpolation of total V
+        CALL zinterp ( jpi, jpj, jpka, jpka_in, ind,  &
+            &                         humi, e3t_in, e3_bak, e3t, varout, jp_spln  )       
+        varout(:,:,1,1) = varout(:,:,2,1)
+
+        IF(ln_c1d) THEN
+           DO jj=1,3
+              DO ji=1,3
+                 varc1d( ji, jj , 1:jpka+1 , 1 ) = varout(iloc,jloc, 1:jpka+1 , 1 )
+              END DO
+           END DO
+           CALL Write_Ncdf_var ( 'vwnd', dimnames(1:4), out_file, varc1d, kt, 'float' )
+        ELSE
+           CALL Write_Ncdf_var ( 'vwnd', dimnames(1:4), out_file, varout, kt, 'float' )
+        ENDIF
+
+      ELSE ! var_name
+        !
+        ctrl     =  minval(pack([(ji,ji=1,size(outnames))],outnames==var_name))
+        var_file = filnames(ctrl)
+        !
+        !
+        ! Interpolation of var_name
+        !      
+        IF ( (var_name.NE."tair").AND.(var_name.NE."humi") ) THEN
+          ! Read var_name
+          CALL Read_Ncdf_var ( varnames(ctrl) , trim(cn_dir)//'/'//trim(var_file), tair, kt )
+          CALL flip_vert_dim ( 1, jpka_in, jpi, jpj, tair( :,:,:,1 ) )       
+          ! Interpolation of var_name
+          CALL zinterp ( jpi, jpj, jpka, jpka_in, ind,  &
+              &                         tair, e3t_in, e3_bak, e3t, varout, jp_spln  )
+        ELSE
+          ! humi and tpot already read     
+          IF (var_name.EQ."humi") tair = humi
+          IF (var_name.EQ."tair") tair = tpot
+          ! Interpolation of var_name
+          CALL zinterp ( jpi, jpj, jpka, jpka_in, ind,  &
+              &                         tair, e3t_in, e3_bak, e3t, varout, jp_weno )
+        END IF
+
+        varout(:,:,1,1) = varout(:,:,2,1)
+
+        IF(ln_c1d) THEN
+           DO jj=1,3
+              DO ji=1,3
+                 varc1d( ji, jj , 1:jpka+1 , 1 ) = varout(iloc,jloc, 1:jpka+1 , 1 )
+              END DO
+           END DO
+           CALL Write_Ncdf_var ( outnames(ctrl), dimnames(1:4), out_file, varc1d, kt, 'float' )
+        ELSE
+           CALL Write_Ncdf_var ( outnames(ctrl), dimnames(1:4), out_file, varout, kt, 'float' )
+        ENDIF
+        !
+      ENDIF ! var_name
+   !
+   END DO ! kt
+   !
+   DEALLOCATE(zsurf,slp,tair,humi,varout)     
+   IF (jpka_in.NE.nhym) DEALLOCATE(tmp_fullw,tmp_fullm)
+   !
+   STOP
+   !
+END PROGRAM main
Index: /utils/tools/ABL_TOOLS/module_grid.F90
===================================================================
--- /utils/tools/ABL_TOOLS/module_grid.F90	(revision 11589)
+++ /utils/tools/ABL_TOOLS/module_grid.F90	(revision 11589)
@@ -0,0 +1,904 @@
+MODULE module_grid
+   !!======================================================================
+   !!                   ***  MODULE  module_grid  ***
+   !! ABL utilities to define and store vertical grids 
+   !!=====================================================================
+   !! History : 2016-10  (F. Lemarié)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   FUNCTIONS   : None    
+   !!   SUBROUTINES : get_atm_grid, init_atm_mask, get_pot_temp
+   !!                 Write_Grid_File, Init_output_File, init_target_grid
+   !!                 flip_vert_dim, smooth_field
+   !!----------------------------------------------------------------------
+   IMPLICIT NONE
+
+CONTAINS
+
+   SUBROUTINE get_atm_grid( jpi, jpj, jpka, slp, temp, humi, Aw, Bw, &   ! in
+      &                                                      e3t, ghw )   ! out
+      !!--------------------------------------------------------------------------
+      !!                    ***  ROUTINE get_atm_grid  ***
+      !!                   
+      !! ** Purpose : compute layer thickness and altitude of interfaces   
+      !!              of the ECMWF atmospheric grid
+      !!
+      !! ** Method  : 
+      !!        (1) recompute the pressure levels thanks to the sea-level pressure
+      !!        (2) use the hydrostatic relation to convert pressure into altitudes
+      !!
+      !!---------------------------------------------------------------------------      
+      INTEGER, INTENT(in   ) ::   jpi, jpj
+      INTEGER, INTENT(in   ) ::   jpka
+      REAL(8), INTENT(in   ) ::   slp   ( 1:jpi, 1:jpj         )
+      REAL(8), INTENT(in   ) ::   temp  ( 1:jpi, 1:jpj, 1:jpka )      
+      REAL(8), INTENT(in   ) ::   humi  ( 1:jpi, 1:jpj, 1:jpka ) 
+      REAL(8), INTENT(in   ) ::   Aw    (               0:jpka )
+      REAL(8), INTENT(in   ) ::   Bw    (               0:jpka )                 
+      REAL(8), INTENT(  out) ::   e3t   ( 1:jpi, 1:jpj, 1:jpka )  
+      REAL(8)                ::   ghw   ( 1:jpi, 1:jpj, 0:jpka )                           
+      REAL(8), PARAMETER     ::   g    =   9.80665
+      REAL(8), PARAMETER     ::   Rd   = 287.058   
+      REAL(8), PARAMETER     ::   zvir =   0.609133 
+      REAL(8), PARAMETER     ::   ig   =   1./g
+      !!
+      INTEGER                :: ji,jj,jk
+      REAL(8)                :: tv,ph(0:jpka)
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+
+            DO jk=0,jpka
+               ph(jk) = Aw( jk ) + Bw( jk ) * slp( ji, jj )    !<-- Pa
+            END DO
+            
+            DO jk=jpka,1,-1
+              tv      = temp( ji, jj, jk ) * ( 1. + zvir*humi( ji, jj, jk ) )  !<-- Virtual temperature
+              e3t ( ji, jj, jk   ) =  ig*( Rd * tv * log( ph( jk ) / ph( jk-1 ) ) )
+              ghw ( ji, jj, jk-1 ) =  e3t( ji, jj, jk ) + ghw( ji, jj, jk )
+            END DO  
+
+         END DO
+      END DO
+      !
+   END SUBROUTINE get_atm_grid
+
+
+
+
+   SUBROUTINE init_atm_mask( jpi, jpj, mask_file, mask_name, ln_lsm_land, tmask )
+      USE module_io
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE INIT_atm_MASK  ***
+      !!                   
+      !! ** Purpose : extract the land/sea mask and remove isolated sea points 
+      !!
+      !! ** Method  : mask is 1 over the ocean and 0 over land 
+      !!
+      !!----------------------------------------------------------------------         
+      INTEGER, INTENT(in   ) ::   jpi, jpj 
+      LOGICAL, INTENT(in   ) ::   ln_lsm_land
+      REAL(8)                ::   tmask  ( 1:jpi, 1:jpj, 1 ) 
+      INTEGER                ::   jj, ji        
+      INTEGER                ::   status,ncid,varid  
+      CHARACTER(len = *  )   ::   mask_file, mask_name
+      REAL(8)                ::   cff
+        
+      ! Read land-sea mask variable
+      status    = nf90_open(trim(mask_file),NF90_NOWRITE,ncid)      
+      status    = nf90_inq_varid(ncid,mask_name,varid)        
+      status    = nf90_get_var(ncid,varid,tmask,start=(/1,1,1/))                                             
+      status    = nf90_close(ncid)  
+      ! invert the mask (1 over the ocean, 0 over land if ln_lsm_land)
+      IF(ln_lsm_land) THEN
+         DO jj=1,jpj
+            DO ji=1,jpi
+               IF( tmask(ji,jj,1) <= 0. ) THEN
+                  tmask(ji,jj,1) = 1.     !! Ocean points
+               ELSE
+                  tmask(ji,jj,1) = 0.     !! Land points     
+               END IF
+            END DO
+         END DO  
+      ENDIF
+      ! remove some closed seas
+      DO jj=2,jpj-1
+         DO ji=2,jpi-1
+            cff = MAX( tmask(ji+1,jj  ,1),tmask(ji-1,jj  ,1), &
+               &       tmask(ji  ,jj+1,1),tmask(ji  ,jj-1,1), &
+               &       tmask(ji+1,jj+1,1),tmask(ji-1,jj-1,1), &
+               &       tmask(ji+1,jj-1,1),tmask(ji-1,jj+1,1) )     
+            IF( tmask( ji, jj, 1 ) .gt. 0.5 .and. cff .lt. 0.5 ) THEN
+               tmask( ji, jj, 1 ) = 0.
+            END IF
+         END DO
+      END DO
+      WRITE(*,*)' init_atm_mask: ',mask_name,' in ',mask_file, ' OK'
+      !!---------------------------------------------------------------------- 
+      !
+   END SUBROUTINE init_atm_mask
+
+
+       
+   SUBROUTINE get_pot_temp( jpi, jpj, jpka, slp, temp, Aw, Bw, tmask, method, hum1, z1 )         
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE get_pot_temp  ***
+      !!                   
+      !! ** Purpose :  compute the potential temperature based on the 
+      !!               absolute temperature and the sea level pressure
+      !!
+      !! ** Method  : five different ways are implemented depending on the 
+      !!              value of 'method'
+      !!       (0) potential temperature = absolute temperature (not recommended)
+      !!       (1) potential temperature is computed using a local reference 
+      !!                                pressure equal to the sea-level-pressure
+      !!       (2) potential temperature is computed only on a perturbation of 
+      !!                                  the absolute temperature around t0
+      !!       (3) a local reference pressure is used consistently with AEROBULK gamma_moist
+      !!       (4) a constant global reference pressure is used
+      !!----------------------------------------------------------------------  
+      INTEGER, INTENT(in   ) ::   jpi, jpj
+      INTEGER, INTENT(in   ) ::   jpka, method
+      REAL(8), INTENT(in   ) ::   slp   ( 1:jpi, 1:jpj         )
+      REAL(8), INTENT(in   ) ::   tmask ( 1:jpi, 1:jpj         )   
+      REAL(8), INTENT(in   ) ::   z1    ( 1:jpi, 1:jpj         )            
+      REAL(8), INTENT(inout) ::   temp  ( 1:jpi, 1:jpj, 1:jpka ) 
+      REAL(8), INTENT(in   ) ::   hum1  ( 1:jpi, 1:jpj         )       
+      REAL(8), INTENT(in   ) ::   Aw    (               0:jpka )
+      REAL(8), INTENT(in   ) ::   Bw    (               0:jpka )                         
+      REAL(8), PARAMETER     ::   grav   =   9.80665
+      REAL(8), PARAMETER     ::   R_dry  =   287.058 
+      REAL(8), PARAMETER     ::   R_vap  =   461.495
+      REAL(8), PARAMETER     ::   Cp_dry =   1005.             
+      REAL(8), PARAMETER     ::   cevap  =   2.5E+06       
+      REAL(8), PARAMETER     ::   reps0  =   R_dry / R_vap        
+      REAL(8), PARAMETER     ::   gamma  =   2./7.
+      REAL(8)                ::   pres0, gamma_moist, zrv, zirt     
+      REAL(8), PARAMETER     ::   t0     =      288.     !<-- K
+      !!
+      INTEGER                :: ji,jj,jk
+      REAL(8)                :: pres,ph(0:jpka), cff         
+      !!
+      SELECT CASE ( method ) 
+      CASE(0) 
+         RETURN
+      CASE(1)         
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+               IF( tmask(ji,jj) .gt. 0.5 ) THEN
+                  DO jk=0,jpka
+                     ph(jk) = Aw( jk ) + Bw( jk ) * slp( ji, jj )
+                  END DO            
+
+                  pres0 = 0.5*(ph(jpka)+ph(jpka-1))
+            
+                  DO jk=1,jpka
+                     pres = 0.5*(ph(jk)+ph(jk-1))
+                     cff  = ( pres0/pres )**gamma
+                     temp( ji, jj, jk ) = cff * temp( ji, jj, jk )
+                  END DO
+               ELSE
+                  temp( ji, jj, 1:jpka ) = 273.15
+               END IF
+            END DO
+         END DO        
+      CASE(2)
+         pres0 =   100900.      
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( tmask(ji,jj) .gt. 0.5 ) THEN
+                  DO jk=0,jpka
+                     ph(jk) = Aw( jk ) + Bw( jk ) * slp( ji, jj )
+                  END DO            
+            
+                  DO jk=1,jpka
+                     pres = 0.5*(ph(jk)+ph(jk-1))
+                     cff  = ( pres0/pres )**gamma
+                     temp( ji, jj, jk ) = cff * ( temp( ji, jj, jk ) - t0 ) + t0
+                  END DO
+               ELSE
+                  temp( ji, jj, 1:jpka ) = 273.15
+               END IF
+            END DO
+         END DO  
+      CASE(3)
+         DO jj = 1, jpj
+            DO ji = 1, jpi      
+               IF( tmask(ji,jj) .gt. 0.5 ) THEN
+                  DO jk=0,jpka
+                     ph(jk) = Aw( jk ) + Bw( jk ) * slp( ji, jj )
+                  END DO 
+                  !! compute gamma_moist consistently with AEROBULK
+                  zrv         = hum1( ji, jj ) / ( 1. - hum1( ji, jj ) )
+                  zirt        = 1. / ( R_dry * temp( ji, jj, jpka ) )
+                  gamma_moist = grav * ( 1. + cevap*zrv*ziRT )   & 
+                     &     / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/temp( ji, jj, jpka ) )       
+                  !! pressure at z = z1
+                  pres        = 0.5*(ph(jpka)+ph(jpka-1))                  
+                  pres0       = pres * ( 1. + gamma_moist * z1(ji,jj) / temp( ji, jj, jpka ) )**(1./gamma)
+                  DO jk=1,jpka
+                     pres = 0.5*(ph(jk)+ph(jk-1))
+                     cff  = ( pres0/pres )**gamma
+                     temp( ji, jj, jk ) = cff * temp( ji, jj, jk )
+                  END DO
+               ELSE
+                  temp( ji, jj, 1:jpka ) = 273.15
+               END IF
+            END DO
+         END DO            
+      CASE(4)
+         pres0 =   100900.      
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( tmask(ji,jj) .gt. 0.5 ) THEN
+                  DO jk=0,jpka
+                     ph(jk) = Aw( jk ) + Bw( jk ) * slp( ji, jj )
+                  END DO            
+            
+                  DO jk=1,jpka
+                     pres = 0.5*(ph(jk)+ph(jk-1))
+                     cff  = ( pres0/pres )**gamma
+                     temp( ji, jj, jk ) = cff * temp( ji, jj, jk )
+                  END DO
+               ELSE
+                  temp( ji, jj, 1:jpka ) = 273.15
+               END IF
+            END DO
+         END DO       
+      END SELECT
+      !
+   END SUBROUTINE get_pot_temp
+
+
+   SUBROUTINE Write_Grid_File( jpka, ght, ghw, e3t, e3w, grd_file  )  
+      !!
+      USE module_io
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE write_Grid_File  ***
+      !!                   
+      !! ** Purpose : write the ABL grid file  
+      !!
+      !! ** Method  : store the layer thickness and altitude of grid points 
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER,      INTENT(in   )  ::   jpka
+      REAL(8),      INTENT(in   )  ::   ght ( 1:jpka+1 )
+      REAL(8),      INTENT(in   )  ::   ghw ( 1:jpka+1 )      
+      REAL(8),      INTENT(in   )  ::   e3t ( 1:jpka+1 )
+      REAL(8),      INTENT(in   )  ::   e3w ( 1:jpka+1 )            
+      CHARACTER(*), INTENT(in   )  ::   grd_file
+      !!
+      INTEGER                      :: status, ncid
+      !!
+      status   = nf90_create( grd_file, NF90_WRITE, ncid )
+      status   = nf90_close ( ncid )    
+      !!  
+      Call Write_Ncdf_dim ( 'jpka' , grd_file, jpka+1 )       
+      Call Write_Ncdf_var ( 'ghw', 'jpka', grd_file, ghw, 'double' ) 
+      Call Write_Ncdf_var ( 'ght', 'jpka', grd_file, ght, 'double' )
+      Call Write_Ncdf_var ( 'e3t', 'jpka', grd_file, e3t, 'double' )
+      Call Write_Ncdf_var ( 'e3w', 'jpka', grd_file, e3w, 'double' )
+      !
+   END SUBROUTINE Write_Grid_File
+
+       
+
+
+
+   SUBROUTINE Init_output_File ( jpi, jpj, jpka, atm_file, abl_file, tmask  )  
+      !!
+      USE module_io
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE Init_output_File  ***
+      !!                   
+      !! ** Purpose :  write longitude, latitude and mask in the output file
+      !!
+      !! ** Method  : define dimensions in the netcdf file
+      !!----------------------------------------------------------------------
+      INTEGER,                    INTENT(in   )  ::   jpi,jpj,jpka
+      CHARACTER(*),               INTENT(in   )  ::   atm_file      
+      CHARACTER(*),               INTENT(in   )  ::   abl_file
+      REAL(8)                                    ::   tmask(1:jpi,1:jpj)
+      !!
+      REAL(8), ALLOCATABLE, DIMENSION(:      )   ::   tmp1d
+      INTEGER                                    ::   status, ncid
+      CHARACTER(len= 20),DIMENSION(4)            :: dimnames
+      !!
+      status   = nf90_create( abl_file, NF90_WRITE, ncid )
+      status   = nf90_close ( ncid )    
+      !!     
+      CALL Write_Ncdf_dim ( 'lon'    , abl_file, jpi     )
+      CALL Write_Ncdf_dim ( 'lat'    , abl_file, jpj     )
+      CALL Write_Ncdf_dim ( 'jpka'   , abl_file, jpka+1  )         
+      CALL Write_Ncdf_dim ( 'time'   , abl_file, 0       )     
+      !
+      ALLOCATE( tmp1d( 1:jpi ) )
+      !
+      CALL Read_Ncdf_var  ( 'lon', atm_file, tmp1d )
+      CALL Write_Ncdf_var ( 'lon', 'lon', abl_file, tmp1d, 'double' ) 
+      !
+      DEALLOCATE( tmp1d )
+      !
+      ALLOCATE( tmp1d( 1:jpj ) )           
+      !
+      CALL Read_Ncdf_var ( 'lat', atm_file, tmp1d )   
+      CALL Write_Ncdf_var( 'lat', 'lat', abl_file, tmp1d, 'double' )     
+      !
+      DEALLOCATE( tmp1d )   
+      !
+      dimnames(1) = 'lon'
+      dimnames(2) = 'lat'
+      CALL Write_Ncdf_var( 'lsm', dimnames , abl_file, tmask, 'float' )
+      !
+   END SUBROUTINE Init_output_File
+       
+       
+   SUBROUTINE Init_output_File_c1d ( jpi, jpj, jpka, atm_file, abl_file, tmask, iloc, jloc )
+      !!
+      USE module_io
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE Init_output_File  ***
+      !!
+      !! ** Purpose :  write longitude, latitude and mask in the output file
+      !!
+      !! ** Method  : define dimensions in the netcdf file
+      !!----------------------------------------------------------------------
+      INTEGER,                    INTENT(in   )  ::   iloc,jloc,jpi,jpj,jpka
+      CHARACTER(*),               INTENT(in   )  ::   atm_file
+      CHARACTER(*),               INTENT(in   )  ::   abl_file
+      REAL(8)                                    ::   tmask(1:jpi,1:jpj)
+      !!
+      REAL(8), ALLOCATABLE, DIMENSION(:      )   ::   tmp1d
+      REAL(8)                                    ::   tmp1d_loc(3),tmp2d_loc(3,3)
+      INTEGER                                    ::   status, ncid
+      CHARACTER(len= 20),DIMENSION(4)            :: dimnames
+      !!
+      status   = nf90_create( abl_file, NF90_WRITE, ncid )
+      status   = nf90_close ( ncid )
+      !!
+      CALL Write_Ncdf_dim ( 'lon'    , abl_file, 3       )
+      CALL Write_Ncdf_dim ( 'lat'    , abl_file, 3       )
+      CALL Write_Ncdf_dim ( 'jpka'   , abl_file, jpka+1  )
+      CALL Write_Ncdf_dim ( 'time'   , abl_file, 0       )
+      !
+      ALLOCATE( tmp1d    ( 1:jpi ) )
+      !
+      CALL Read_Ncdf_var  ( 'lon', atm_file, tmp1d )
+      tmp1d_loc( 1 ) = tmp1d(iloc)
+      tmp1d_loc( 2 ) = tmp1d(iloc)
+      tmp1d_loc( 3 ) = tmp1d(iloc)
+      CALL Write_Ncdf_var ( 'lon', 'lon', abl_file, tmp1d_loc, 'double' )
+      print*,'1D column located at : '
+      print*,'longitude = ',tmp1d(iloc),' degree_east'
+      !
+      DEALLOCATE( tmp1d )
+      !
+      ALLOCATE( tmp1d( 1:jpj ) )
+      !
+      CALL Read_Ncdf_var ( 'lat', atm_file, tmp1d )
+      tmp1d_loc( 1 ) = tmp1d(jloc)
+      tmp1d_loc( 2 ) = tmp1d(jloc)
+      tmp1d_loc( 3 ) = tmp1d(jloc)
+      CALL Write_Ncdf_var( 'lat', 'lat', abl_file, tmp1d_loc, 'double' )
+      print*,'latitude = ',tmp1d(jloc),' degree_north'
+      !
+      DEALLOCATE( tmp1d )
+      !
+      tmp2d_loc(:,:) = tmask(iloc,jloc)       
+      dimnames(1) = 'lon'
+      dimnames(2) = 'lat'
+      CALL Write_Ncdf_var( 'lsm', dimnames , abl_file, tmp2d_loc, 'float' )
+      !
+   END SUBROUTINE Init_output_File_c1d       
+       
+
+   SUBROUTINE init_target_grid( jpka, ght, ghw, e3t, e3w, hmax, hc, theta_s,  &
+      &                                                        force_z1, z1 )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE init_target_grid  ***
+      !!                   
+      !! ** Purpose : compute the layer thickness and altitude of grid points 
+      !!              for the ABL model based on the namelist parameter values  
+      !!
+      !! ** Method  : depending on the logical 'force_z1' two methods are used
+      !!              (1) true -> the user chooses the value of the first vertical
+      !!                  grid point. A few Newton iterations are used to correct  
+      !!                  the value of the parameter theta_s to satisfy this constraint
+      !!              (2) false -> use the parameter values in the namelist to 
+      !!                  compute the vertical grid
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, intent(in   )    :: jpka
+      REAL(8), intent(in   )    :: hmax,hc
+      REAL(8), intent(inout)    :: theta_s      
+      LOGICAL, intent(in   )    :: force_z1
+      REAL(8), intent(in   )    :: z1
+      !!
+      REAL(8), intent(  out)    :: ghw( 1:jpka+1 )
+      REAL(8), intent(  out)    :: ght( 1:jpka+1 )      
+      REAL(8), intent(  out)    :: e3w( 1:jpka+1 )
+      REAL(8), intent(  out)    :: e3t( 1:jpka+1 ) 
+      !! 
+      REAL(8)                   ::  ds,cff,sc_w,sc_r,alpha,x
+      REAL(8)                   ::  fx,fxp
+      INTEGER                   ::  jk,maxiter,jiter
+      REAL(8), PARAMETER        :: tol = 1.E-12
+      !!
+      IF(force_z1) THEN
+        IF(z1.LT.10.) THEN
+          WRITE(*,*) " ERROR: z1 < 1st ECMWF level height (~10m)"
+          STOP
+        ELSE
+          !! Newton iterations to find the appropriate value of theta_s
+          maxiter = 1000
+          x       = theta_s
+          sc_r    = (float(1)-0.5)/float(jpka)
+          alpha   = (z1 - hc*sc_r) / (hmax - hc)
+          !
+          DO jiter=1,maxiter
+             fx   = (sinh(sc_r*x)/sinh(x))-alpha
+             fxp  = (sc_r*cosh(sc_r*x)-sinh(sc_r*x)*cosh(x)/sinh(x))/sinh(x)
+             IF( abs(fx) .lt. tol ) THEN
+                   exit
+             ENDIF
+             cff  = fx / fxp
+             x    = x - cff      
+          ENDDO
+          !
+          theta_s = x
+        END IF
+        !   
+      ENDIF
+      !
+      ds =1./float(jpka)
+      cff=(hmax-hc)/sinh(theta_s)
+      !
+      DO jk = jpka,1,-1
+        sc_w      = ds*float(jk)
+        ghw(jk+1) = hc*sc_w + cff*sinh(theta_s*sc_w)
+        sc_r      = ds*(float(jk)-0.5)
+        ght(jk+1) = hc*sc_r + cff*sinh(theta_s*sc_r)
+      END DO
+      !
+      ghw(1) = 0.
+      e3t(1) = 0.
+      ght(1) = 0.
+      !
+      DO jk = 2,jpka+1
+        e3t(jk) = ghw(jk)-ghw(jk-1)
+      END DO
+      !
+      DO jk=1,jpka
+        e3w(jk)  = ght(jk+1)-ght(jk)
+      END DO
+      !
+      e3w(jpka+1) = ghw(jpka+1) - ght(jpka+1)          
+      !
+      IF(force_z1) THEN   !++ print the new parameter values
+         print*,'*** Updated grid parameters'
+         print*,'theta_s = ',theta_s
+         print*,'hc      = ',hc      
+         print*,'hmax    = ',hmax         
+         print*,'ght(2)  = ',ght(2)
+      ENDIF      
+      !
+   END SUBROUTINE init_target_grid
+   
+   
+   
+   
+   
+
+   SUBROUTINE flip_vert_dim( kstr, kend, jpi, jpj, tabin )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE flip_vert_dim  ***
+      !!                   
+      !! ** Purpose : flip the vertical axis of the array tabin so that    
+      !!              the vertical grid goes from k=kstr at the bottom 
+      !!              of the ABL to k=kend at the top
+      !!----------------------------------------------------------------------
+      INTEGER, intent(in   )   :: kstr,kend
+      INTEGER, intent(in   )   :: jpi, jpj 
+      REAL(8), intent(inout)   :: tabin( 1:jpi, 1:jpj, kstr:kend )
+      !!
+      INTEGER                  :: ji,jj,jk,ks
+      REAL(8)                  :: tabfl(kstr:kend)
+      !
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            DO jk=kstr,kend
+               ks=(kend-jk)+kstr
+               tabfl(ks) = tabin( ji, jj, jk ) 
+            END DO
+            tabin( ji, jj, kstr:kend ) = tabfl( kstr:kend )
+         END DO 
+      END DO 
+      !
+   END SUBROUTINE flip_vert_dim
+   
+  
+
+
+
+
+
+   SUBROUTINE smooth_field( jpi, jpj, varin, tmask, niter )         
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE smooth_field  ***
+      !!                   
+      !! ** Purpose : smooth the sea level pressure over the ocean
+      !!              to attenuate Gibbs oscillation 
+      !! ** Method  : a 9-point isotropic laplacian filter is applied
+      !!              iteratively on ocean grid points only
+      !!
+      !! Proper treatment of the periodicity is still not yet implemented
+      !!
+      !!----------------------------------------------------------------------     
+      INTEGER, INTENT(in   ) ::   jpi, jpj, niter  
+      REAL(8), INTENT(in   ) ::   tmask  ( 1:jpi, 1:jpj )        
+      REAL(8), INTENT(inout) ::   varin  ( 1:jpi, 1:jpj )
+      INTEGER                :: ji,jj,nit
+      REAL(8)                :: smth_a,smth_b,umask,vmask 
+      REAL(8)                :: FX ( 0:jpi  , 0:jpj+1 )
+      REAL(8)                :: FE1( 0:jpi+1, 0:jpj   )
+      REAL(8)                :: FE ( 1:jpi  , 0:jpj   )       
+      !!  
+      !!=========================================================
+      !!    
+      ! Hanning filter
+      !smth_a = 1./8.
+      !smth_b = 1./4.     
+      ! 9-point isotropic laplacian filter
+      smth_a = 1./12.
+      smth_b = 3./16.      
+
+      FX ( 0:jpi  , 0:jpj+1 ) = 0.
+      FE1( 0:jpi+1, 0:jpj   ) = 0.
+      FE ( 1:jpi  , 0:jpj   ) = 0.
+
+      !!+++++++++
+      DO nit = 1,niter
+      !!+++++++++
+         DO jj=1,jpj              
+            DO ji=1,jpi-1          
+            umask = tmask(ji,jj)*tmask(ji+1,jj)  
+            FX ( ji, jj ) = ( varin( ji+1,jj ) - varin( ji  ,jj ) ) * umask  
+            END DO                    
+         END DO
+         FX(   0  , 1:jpj ) = FX( jpi-1, 1:jpj  )
+         FX  ( jpi, 1:jpj ) = FX(     1, 1:jpj  )      
+         FX( 0:jpi, 0     ) = FX( 0:jpi,     1  )
+         FX( 0:jpi, jpj+1 ) = FX( 0:jpi,     jpj)           
+         !!
+         DO jj=1,jpj-1                 
+            DO ji=1,jpi
+               vmask = tmask(ji,jj)*tmask(ji,jj+1)  
+               FE1( ji, jj ) = ( varin( ji, jj+1 ) - varin( ji  ,jj ) ) * vmask 
+           END DO
+         END DO
+         !! 
+         FE1( 0      , 1:jpj-1 ) = FE1( jpi , 1:jpj-1 )
+         FE1( jpi+1  , 1:jpj-1 ) = FE1(   1 , 1:jpj-1 )      
+         FE1( 0:jpi+1,    0    ) = 0.
+         FE1( 0:jpi+1,    jpj  ) = 0. 
+         !!
+         DO jj=0,jpj                 
+            DO ji=1,jpi
+               FE ( ji, jj ) =   FE1( ji, jj   )                   &
+               & +  smth_a*(  FX ( ji, jj+1 ) + FX( ji-1, jj   )   &
+               &           -  FX ( ji, jj   ) - FX( ji-1, jj+1 ) )
+            END DO
+         END DO
+         !!
+         DO jj = 1, jpj
+            DO ji = 0,jpi
+               FX( ji, jj ) = FX( ji, jj )         &
+              &    + smth_a*( FE1( ji+1, jj ) + FE1( ji  , jj-1 )    &
+              &              -FE1( ji  , jj ) - FE1( ji+1, jj-1 ) )
+            END DO
+            DO ji = 1,jpi
+               varin( ji  ,jj ) = varin( ji  ,jj ) &
+     &         + tmask(ji,jj) * smth_b * (         &
+     &                  FX( ji, jj ) - FX( ji-1, jj )          &
+     &                 +FE( ji, jj ) - FE( ji, jj-1 )  )
+            END DO
+         END DO      
+         !!
+      !!+++++++++      
+      END DO
+      !!+++++++++   
+            
+   END SUBROUTINE smooth_field
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+   
+
+  
+
+
+
+
+!   SUBROUTINE DTV_Filter( jpi, jpj, varin, tmask, niter, time )         
+!      USE module_io
+!      !!---------------------------------------------------------------------
+!      !!                    ***  ROUTINE DTV_Filter  ***
+!      !!                   
+!      !! ** Purpose :   
+!      !!
+!      !! ** Method  : 
+!      !!
+!      !! ** Action  : 
+!      !!----------------------------------------------------------------------   
+!      INTEGER, INTENT(in   ) ::   jpi, jpj, niter, time 
+!      REAL(8), INTENT(in   ) ::   tmask  ( 1:jpi, 1:jpj )        
+!      REAL(8), INTENT(inout) ::   varin  ( 1:jpi, 1:jpj )
+!      INTEGER                ::   jn, jj, ji, nit, nt_n, nt_a, nl2      
+!      REAL(8)                ::   var0,mean,sigma2,cff,lambda
+!      REAL(8)                ::   FX(0:jpi+1,0:jpj+1)
+!      REAL(8)                ::   FE(0:jpi+1,0:jpj+1)
+!      REAL(8)                ::   FL(0:jpi+1,0:jpj+1)
+!      REAL(8)                ::   FR(0:jpi+1,0:jpj+1)  
+!      REAL(8)                ::   wrk ( 1:jpi, 1:jpj, 2 ) 
+!      REAL(8)                ::   div ( 0:jpi+1, 0:jpj+1 ), diag               
+!      REAL(8), PARAMETER     ::   rsmall = 1.E-08
+!      REAL(8), PARAMETER     ::   rbig   = 1.E+14
+!      REAL(8)                ::   umask, vmask, fmask, wght(8), dt
+!      REAL(8)                ::   L2norm_n, L2norm_a
+!      REAL(8), ALLOCATABLE, DIMENSION(:      )   ::   tmp1d      
+!      INTEGER                :: status,ncid,varid1,varid2,dimid1,dimid2  
+!      !!=========================================================
+!      CHARACTER(len=  3  )   :: nn
+!      CHARACTER(len = 500)   :: erai_file, smth_file
+!      LOGICAL                :: ln_diag_smoothing
+!      ln_diag_smoothing = .true.
+!      IF( ln_diag_smoothing ) THEN
+!         erai_file = '/Users/florianlemarie/Documents/INRIA/SIMBAD/ERA_INTERIM_DECEMBRE_2015/phi_ml1_6h_erai_201512.nc'
+!         WRITE(nn,'(I3)') time
+!         DO ji=1,3
+!            IF(nn(ji:ji)==' ') nn(ji:ji)='0'
+!         END DO
+!      
+!         smth_file = 'smoothing_results_DTV_'//nn//'.nc'  
+!         status    = nf90_create( trim(smth_file) , NF90_WRITE, ncid )
+!         status    = nf90_close ( ncid )    
+!         !   
+!         CALL Write_Ncdf_dim ( 'lon'    , trim(smth_file) , jpi     )
+!         CALL Write_Ncdf_dim ( 'lat'    , trim(smth_file) , jpj     )
+!         !
+!         ALLOCATE( tmp1d( 1:jpi ) )
+!         CALL Read_Ncdf_var  ( 'lon', erai_file, tmp1d )
+!         CALL Write_Ncdf_var ( 'lon', 'lon', trim(smth_file), tmp1d, 'double' ) 
+!         DEALLOCATE( tmp1d )
+!         !
+!         ALLOCATE( tmp1d( 1:jpj ) )           
+!         CALL Read_Ncdf_var ( 'lat', erai_file, tmp1d )   
+!         CALL Write_Ncdf_var( 'lat', 'lat', trim(smth_file), tmp1d, 'double' )     
+!         DEALLOCATE( tmp1d ) 
+!         !         
+!         status = nf90_open(trim(smth_file),NF90_WRITE,ncid)     
+!         status = nf90_inq_dimid(ncid, 'lon', dimid1)
+!         status = nf90_inq_dimid(ncid, 'lat', dimid2)
+!         status = nf90_redef(ncid)      
+!         status = nf90_def_var(ncid,'varinp',nf90_double,(/dimid1,dimid2/),varid1)        
+!         status = nf90_def_var(ncid,'varout',nf90_double,(/dimid1,dimid2/),varid2)   
+!         status = nf90_enddef(ncid)
+!         ! 
+!         !
+!         status = nf90_put_var(ncid,varid1,varin)
+!      END IF
+!      !!=========================================================
+!      nit    = 0
+!      nl2    = 0
+!      L2norm_a = 0.
+!      L2norm_n = rbig      
+!      nt_n   = 1 + MOD( nit  , 2 )
+!      nt_a   = 1 + MOD( nit+1, 2 )      
+!      var0   = varin(  nint(0.5*jpi) , nint(0.5*jpj) )  
+!      varin( 1:jpi, 1:jpj       ) = varin( 1:jpi, 1:jpj ) - var0
+!      wrk  ( 1:jpi, 1:jpj, nt_n ) = varin( 1:jpi, 1:jpj )
+!      !! Compute the mean first
+!      mean   = Get_Mean ( jpi , jpj,  varin, tmask       )
+!      sigma2 = Get_Vari ( jpi , jpj,  varin, tmask, mean )
+!      !!
+!      print*,'mean value over the ocean = ',mean + var0
+!      print*,'variance of the input field = ',sigma2      
+!      !!
+!      lambda = (1. / sigma2      )
+!      dt     =  1. / (16.+lambda)
+!      diag   =  0.
+!      
+!      !!>>>>>>>>>>>>>>>>
+!      DO nit = 1,niter
+!      !!>>>>>>>>>>>>>>>>
+!      FX(0:jpi+1,0:jpj+1) = 0.
+!      FE(0:jpi+1,0:jpj+1) = 0.
+!      FL(0:jpi+1,0:jpj+1) = 0.
+!      FR(0:jpi+1,0:jpj+1) = 0.
+!      !!
+!      DO jj = 1,jpj
+!         DO ji = 1,jpi-1
+!            umask        = tmask(ji,jj)*tmask(ji+1,jj)  
+!            FX( ji, jj ) = umask * ( wrk( ji+1, jj, nt_n ) - wrk( ji, jj, nt_n ) )
+!         END DO
+!      END DO  
+!      !!
+!      DO jj = 1,jpj-1
+!         DO ji = 1,jpi
+!            vmask        = tmask(ji,jj)*tmask(ji,jj+1)
+!            FE( ji, jj ) = vmask * ( wrk( ji, jj+1, nt_n ) - wrk( ji, jj, nt_n ) )
+!         END DO
+!      END DO       
+!      !!
+!      DO jj = 1,jpj-1
+!         DO ji = 1,jpi-1
+!            fmask        = tmask(ji,jj)*tmask(ji+1,jj+1)*diag
+!            FL( ji, jj ) = fmask * ( wrk( ji+1, jj+1, nt_n ) - wrk( ji  , jj, nt_n   ) )
+!            fmask        = tmask(ji+1,jj)*tmask(ji,jj+1)*diag       
+!            FR( ji, jj ) = fmask * ( wrk( ji  , jj+1, nt_n ) - wrk( ji+1, jj, nt_n ) )
+!         END DO
+!      END DO 
+!      !!
+!      div( 0:jpi+1, 0:jpj+1 ) = rsmall
+!      !!
+!      DO jj = 1,jpj+1
+!         DO ji = 1,jpi+1           
+!            div(ji,jj) = MAX( sqrt(  FX(ji,jj)**2 + FX(ji-1,jj  )**2 + FE(ji,jj  )**2 + FX(ji,jj-1)**2          &
+!            &                      + FL(ji,jj)**2 + FL(ji-1,jj-1)**2 + FR(ji,jj-1)**2 + FR(ji-1,jj)**2 ), rsmall ) 
+!         END DO
+!      END DO
+!      
+!      
+!       
+!      
+!      DO jj=1,jpj
+!         DO ji=1,jpi
+!            IF( div(ji,jj) .eq. rsmall .or. tmask(ji,jj).lt.0.5 ) THEN
+!               wrk(ji,jj,nt_a) = wrk(ji,jj,nt_n)
+!            ELSE
+!               wght(1)     = 1. +  div(ji,jj) / div(ji+1,jj  ) 
+!               wght(2)     = 1. +  div(ji,jj) / div(ji-1,jj  ) 
+!               wght(3)     = 1. +  div(ji,jj) / div(ji  ,jj+1)                        
+!               wght(4)     = 1. +  div(ji,jj) / div(ji  ,jj-1)           
+!               wght(5)     = 1. +  div(ji,jj) / div(ji+1,jj+1)   
+!               wght(6)     = 1. +  div(ji,jj) / div(ji-1,jj+1)
+!               wght(7)     = 1. +  div(ji,jj) / div(ji-1,jj-1)                       
+!               wght(8)     = 1. +  div(ji,jj) / div(ji+1,jj-1)             
+!              
+!            wrk(ji,jj,nt_a) = wrk(ji,jj,nt_n) +  dt*(                   &
+!               &            + FX(ji  ,jj  ) * wght(1)                   &
+!               &            - FX(ji-1,jj  ) * wght(2)                   &
+!               &            + FE(ji  ,jj  ) * wght(3)                   &
+!               &            - FE(ji  ,jj-1) * wght(4)                   &
+!               &            + FL(ji  ,jj  ) * wght(5)                   &
+!               &            + FR(ji-1,jj  ) * wght(6)                   &
+!               &            - FL(ji-1,jj-1) * wght(7)                   &
+!               &            - FR(ji  ,jj-1) * wght(8)                   &
+!               &            - lambda * div(ji,jj) * ( wrk(ji,jj,nt_n) - varin(ji,jj) ) )  
+!               
+!              IF( isnan(wrk(ji,jj,nt_a)) ) THEN 
+!                  print*,'Nan in smoothing at iteration ',nit
+!                  print*,'at grid point ',ji,jj
+!              END IF    
+!              IF( abs(wrk(ji,jj,nt_a)) .gt. rbig ) THEN
+!                  print*,'Inf in smoothing at iteration ',nit
+!                  print*,'at grid point ',ji,jj                                    
+!              END IF
+!                 L2norm_a = L2norm_a + ( wrk(ji,jj,nt_n)-wrk(ji,jj,nt_a) )**2
+!              END IF                 
+!         END DO
+!      END DO
+!
+!
+!      mean   = Get_Mean ( jpi , jpj,  wrk(:,:,nt_a), tmask       )
+!      sigma2 = Get_Vari ( jpi , jpj,  wrk(:,:,nt_a), tmask, mean )
+!      !!
+!      print*,'mean value over the ocean = ',mean + var0
+!      print*,'variance of the input field = ',sigma2   
+!
+!      
+!      IF( L2norm_a .gt. L2norm_n ) THEN
+!         print*,'convergence after ',nit,' iterations'
+!!         EXIT
+!      END IF     
+!      L2norm_n = L2norm_a 
+!      L2norm_a = 0.     
+!      nt_n = 1 + MOD( nit  , 2 )
+!      nt_a = 1 + MOD( nit+1, 2 ) 
+!           
+!      !!>>>>>>>>>>>>>>>>
+!      END DO
+!      !!>>>>>>>>>>>>>>>>
+!
+!       
+!      DO jj=1,jpj
+!         DO ji=1,jpi
+!            varin(ji,jj) = ( wrk(ji,jj,nt_n) + var0 )*tmask(ji,jj)
+!         END DO
+!      END DO
+!
+!      IF( ln_diag_smoothing ) THEN
+!         status = nf90_put_var(ncid,varid2,varin)
+!         status = nf90_close(ncid)
+!      END IF   
+!      !
+!   END SUBROUTINE dtv_filter
+!
+!
+!
+!   REAL(8) FUNCTION Get_Mean ( jpi , jpj, tabvar, tmask )
+!      !!---------------------------------------------------------------------
+!      !!                    ***  FUNCTION Get_Mean  ***
+!      !!                   
+!      !! ** Purpose : get the mean of the input field 
+!      !!
+!      !!----------------------------------------------------------------------
+!      INTEGER,   INTENT(in)      :: jpi,jpj
+!      REAL(8),   INTENT(in)      :: tabvar(jpi,jpj)
+!      REAL(8),   INTENT(in)      :: tmask (jpi,jpj)
+!      INTEGER                    :: ji,jj,jn    
+!      
+!      !! Compute the mean first
+!      jn       = 0
+!      Get_Mean = 0.
+!      DO jj = 1, jpj
+!         DO ji = 1, jpi
+!            IF(tmask(ji,jj) .gt. 0.5) THEN
+!               Get_Mean = Get_Mean + tabvar( ji,jj )  
+!               jn=jn+1
+!            END IF  
+!         END DO
+!      END DO      
+!      Get_Mean   = (1./jn)   * Get_Mean  
+!      !
+!   END FUNCTION Get_Mean
+!      
+!
+!
+!
+!
+!   REAL(8) FUNCTION Get_Vari ( jpi , jpj, tabvar, tmask, mean )
+!      !!---------------------------------------------------------------------
+!      !!                    ***  FUNCTION Get_Mean  ***
+!      !!                   
+!      !! ** Purpose : get the mean of the input field 
+!      !!
+!      !!----------------------------------------------------------------------
+!      INTEGER,   INTENT(in)      :: jpi,jpj
+!      REAL(8),   INTENT(in)      :: tabvar(jpi,jpj)
+!      REAL(8),   INTENT(in)      :: tmask (jpi,jpj), mean
+!      INTEGER                    :: ji,jj,jn    
+!      
+!      jn       = 0
+!      Get_Vari = 0
+!      !!
+!      DO jj = 1, jpj
+!         DO ji = 1, jpi
+!            IF(tmask(ji,jj) .gt. 0.5) THEN
+!              Get_Vari = Get_Vari + ( tabvar(ji,jj) - mean )**2 
+!              jn = jn+1 
+!            ENDIF  
+!         END DO
+!      END DO 
+!      Get_Vari   = (1./jn)   * Get_Vari 
+!      !
+!   END FUNCTION Get_Vari
+
+
+   
+END MODULE module_grid
Index: /utils/tools/ABL_TOOLS/module_interp.F90
===================================================================
--- /utils/tools/ABL_TOOLS/module_interp.F90	(revision 11589)
+++ /utils/tools/ABL_TOOLS/module_interp.F90	(revision 11589)
@@ -0,0 +1,608 @@
+MODULE module_interp
+   !!======================================================================
+   !!                   ***  MODULE  module_interp  ***
+   !! Ocean forcing:  bulk thermohaline forcing of the ocean (or ice)
+   !!=====================================================================
+   !! History : 2016-10  (F. Lemarié)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   zinterp                : 
+   !!   reconstructandremap    : 
+   !!   reconstructandremap_ps :
+   !! 
+   !!----------------------------------------------------------------------
+   IMPLICIT NONE
+
+CONTAINS
+
+   SUBROUTINE zinterp( jpi, jpj, jpka, jpka_in, ind, tab_in, e3t_in, e3_bak, &   
+      &                              e3t_out, tab_out, interp_type           )   
+      
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE zinterp  ***
+      !!                   
+      !! ** Purpose :   
+      !!
+      !! ** Method  :
+      !!
+      !! ** Action  :  
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in   ) ::   jpi, jpj
+      INTEGER, INTENT(in   ) ::   jpka, jpka_in, interp_type
+      INTEGER, INTENT(in   ) ::   ind     ( 1:jpi, 1:jpj            )
+      REAL(8), INTENT(inout) ::   tab_in  ( 1:jpi, 1:jpj, 1:jpka_in )
+      REAL(8), INTENT(in   ) ::   e3t_in  ( 1:jpi, 1:jpj, 1:jpka_in )
+      REAL(8), INTENT(in   ) ::   e3_bak  ( 1:jpi, 1:jpj            )            
+      REAL(8), INTENT(in   ) ::   e3t_out (               1:jpka+1  )                           
+      REAL(8), INTENT(  out) ::   tab_out ( 1:jpi, 1:jpj, 1:jpka+1  )  
+      !!
+      INTEGER                :: ji,jj,k_in
+      REAL(8)                :: val1,val2,cff
+      
+      SELECT CASE(interp_type)
+      CASE(1)      ! WENO
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               k_in = ind( ji, jj )
+               val1 = tab_in ( ji, jj, k_in-1 )
+               val2 = tab_in ( ji, jj, k_in   )
+               cff  =        val1 * e3_bak( ji, jj         )     &
+                  & +        val2 * e3t_in( ji, jj, k_in-1 )     &
+                  & + (val2-val1) * e3t_in( ji, jj, k_in   )             
+               tab_in( ji, jj, k_in ) = cff / ( e3_bak( ji, jj ) + e3t_in ( ji, jj, k_in-1 ) )
+               !
+               CALL reconstructandremap( tab_in ( ji, jj, 1:k_in   ), e3t_in( ji, jj, 1:k_in ),   &
+                  &                      tab_out( ji, jj, 2:jpka+1 ), e3t_out (    2:jpka+1  ),   &
+                  &                      k_in, jpka )  
+               !
+            END DO
+         END DO   
+      CASE(2)      ! SPLINES
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               k_in = ind( ji, jj )
+               val1 = tab_in ( ji, jj, k_in-1 )
+               val2 = tab_in ( ji, jj, k_in   )
+               cff  =        val1 * e3_bak( ji, jj         )     &
+                  & +        val2 * e3t_in( ji, jj, k_in-1 )     &
+                  & + (val2-val1) * e3t_in( ji, jj, k_in   )             
+               tab_in( ji, jj, k_in ) = cff / ( e3_bak( ji, jj ) + e3t_in ( ji, jj, k_in-1 ) )
+               !
+               CALL reconstructandremap_ps( tab_in ( ji, jj, 1:k_in   ), e3t_in( ji, jj, 1:k_in ),   &
+                  &                         tab_out( ji, jj, 2:jpka+1 ), e3t_out (    2:jpka+1  ),   &
+                  &                         k_in, jpka )  
+               !
+            END DO
+         END DO  
+      CASE DEFAULT
+         WRITE(*,*) "### Error: problem in zinterp, interp_type not set properly"
+         STOP
+      END SELECT      
+      !
+   END SUBROUTINE zinterp     
+
+
+
+
+
+!
+!===================================================================================================
+subroutine reconstructandremap(tabin,hin,tabout,hout,N,Nout)
+!---------------------------------------------------------------------------------------------------
+      implicit none
+      integer            :: N, Nout
+      real(8)            :: tabin(N), tabout(Nout)
+      real(8)            :: hin(N), hout(Nout)
+      real(8)            :: coeffremap(N,3),zwork(N,3)
+      real(8)            :: zwork2(N+1,3)
+      integer            :: k
+      real(8), parameter :: dsmll=1.0d-8  
+      real(8)            :: q,q01,q02,q001,q002,q0
+      real(8)            :: z_win(1:N+1), z_wout(1:Nout+1)
+      real(8),parameter  :: dpthin = 1.D-3
+      integer            :: k1, kbox, ktop, ka, kbot
+      real(8)            :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
+!-----
+
+!----------------------
+      z_win(1)=0.; z_wout(1)= 0.
+      do k=1,N
+       z_win(k+1)=z_win(k)+hin(k)
+      enddo 
+      
+      do k=1,Nout
+       z_wout(k+1)=z_wout(k)+hout(k)       
+      enddo       
+
+        do k=2,N
+          zwork(k,1)=1./(hin(k-1)+hin(k))
+        enddo
+        
+        do k=2,N-1
+          q0 = 1./(hin(k-1)+hin(k)+hin(k+1))
+          zwork(k,2)=hin(k-1)+2.*hin(k)+hin(k+1)
+          zwork(k,3)=q0
+        enddo       
+      
+        do k= 2,N
+        zwork2(k,1)=zwork(k,1)*(tabin(k)-tabin(k-1))
+        enddo
+        
+        coeffremap(:,1) = tabin(:)
+ 
+         do k=2,N-1
+        q001 = hin(k)*zwork2(k+1,1)
+        q002 = hin(k)*zwork2(k,1)        
+        if (q001*q002 < 0) then
+          q001 = 0.
+          q002 = 0.
+        endif
+        q=zwork(k,2)
+        q01=q*zwork2(k+1,1)
+        q02=q*zwork2(k,1)
+        if (abs(q001) > abs(q02)) q001 = q02
+        if (abs(q002) > abs(q01)) q002 = q01
+        
+        q=(q001-q002)*zwork(k,3)
+        q001=q001-q*hin(k+1)
+        q002=q002+q*hin(k-1)
+        
+        coeffremap(k,3)=coeffremap(k,1)+q001
+        coeffremap(k,2)=coeffremap(k,1)-q002
+        
+        zwork2(k,1)=(2.*q001-q002)**2
+        zwork2(k,2)=(2.*q002-q001)**2
+        enddo
+        
+        do k=1,N
+        if     (k.eq.1 .or. k.eq.N .or.   hin(k).le.dpthin) then
+        coeffremap(k,3) = coeffremap(k,1)
+        coeffremap(k,2) = coeffremap(k,1)
+        zwork2(k,1) = 0.
+        zwork2(k,2) = 0.
+        endif
+        enddo
+        
+        do k=2,N
+        q002=max(zwork2(k-1,2),dsmll)
+        q001=max(zwork2(k,1),dsmll)
+        zwork2(k,3)=(q001*coeffremap(k-1,3)+q002*coeffremap(k,2))/(q001+q002)
+        enddo
+        
+        zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)
+        zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)
+ 
+        do k=1,N
+        q01=zwork2(k+1,3)-coeffremap(k,1)
+        q02=coeffremap(k,1)-zwork2(k,3)
+        q001=2.*q01
+        q002=2.*q02
+        if (q01*q02<0) then
+          q01=0.
+          q02=0.
+        elseif (abs(q01)>abs(q002)) then
+          q01=q002
+        elseif (abs(q02)>abs(q001)) then
+          q02=q001
+        endif
+        coeffremap(k,2)=coeffremap(k,1)-q02
+        coeffremap(k,3)=coeffremap(k,1)+q01
+        enddo
+
+      zbot=0.0
+      kbot=1
+      do k=1,Nout
+        ztop=zbot  !top is bottom of previous layer
+        ktop=kbot
+        if     (ztop.ge.z_win(ktop+1)) then
+          ktop=ktop+1
+        endif
+        
+        zbot=z_wout(k+1)
+        zthk=zbot-ztop
+
+        if     (zthk.gt.dpthin .and. ztop.lt.z_wout(Nout+1)) then
+
+          kbot=ktop
+          do while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)
+            kbot=kbot+1
+          enddo
+          zbox=zbot
+          do k1= k+1,Nout
+            if     (z_wout(k1+1)-z_wout(k1).gt.dpthin) then
+              exit !thick layer
+            else
+              zbox=z_wout(k1+1)  !include thin adjacent layers
+              if     (zbox.eq.z_wout(Nout+1)) then
+                exit !at bottom
+              endif
+            endif
+          enddo
+          zthk=zbox-ztop
+
+          kbox=ktop
+          do while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)
+            kbox=kbox+1
+          enddo
+          
+          if     (ktop.eq.kbox) then
+
+
+            if     (z_wout(k)  .ne.z_win(kbox)   .or.z_wout(k+1).ne.z_win(kbox+1)     ) then
+
+              if     (hin(kbox).gt.dpthin) then
+                q001 = (zbox-z_win(kbox))/hin(kbox)
+                q002 = (ztop-z_win(kbox))/hin(kbox)
+                q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)
+                q02=q01-1.+(q001+q002)
+                q0=1.-q01-q02
+              else
+                q0 = 1.0
+                q01 = 0.
+                q02 = 0.
+              endif
+          tabout(k)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)
+              
+            else
+            tabout(k) = tabin(kbox)
+              
+            endif 
+
+          else
+
+            if     (ktop.le.k .and. kbox.ge.k) then
+              ka = k
+            elseif (kbox-ktop.ge.3) then
+              ka = (kbox+ktop)/2
+            elseif (hin(ktop).ge.hin(kbox)) then
+              ka = ktop
+            else
+              ka = kbox
+            endif !choose ka
+
+            offset=coeffremap(ka,1)
+
+            qtop = z_win(ktop+1)-ztop !partial layer thickness
+            if     (hin(ktop).gt.dpthin) then
+              q=(ztop-z_win(ktop))/hin(ktop)
+              q01=q*(q-1.)
+              q02=q01+q
+              q0=1-q01-q02            
+            else
+              q0 = 1.
+              q01 = 0.
+              q02 = 0.
+            endif
+            
+            tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop
+
+            do k1= ktop+1,kbox-1
+              tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)
+            enddo !k1
+
+            
+            qbot = zbox-z_win(kbox) !partial layer thickness
+            if     (hin(kbox).gt.dpthin) then
+              q=qbot/hin(kbox)
+              q01=(q-1.)**2
+              q02=q01-1.+q
+              q0=1-q01-q02                            
+            else
+              q0 = 1.0
+              q01 = 0.
+              q02 = 0.
+            endif
+           
+            tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot
+            
+            rpsum=1.0d0/zthk
+              tabout(k)=offset+tsum*rpsum
+              
+          endif !single or multiple layers
+        else
+        if (k==1) then
+        print *,'problem = ',zthk,z_wout(k+1),hout(1)
+        endif
+         tabout(k) = tabout(k-1)
+          
+        endif !normal:thin layer
+      enddo !k
+            
+      return
+      
+!---------------------------------------------------------------------------------------------------
+end subroutine reconstructandremap
+!===================================================================================================
+!
+
+
+
+
+
+
+
+
+
+!
+!===================================================================================================
+      subroutine reconstructandremap_ps(tabin,hin,tabout,hout,N,Nout)      ! parabloc spline
+!---------------------------------------------------------------------------------------------------      
+      implicit none
+      integer N, Nout
+      real(8) tabin(N), tabout(Nout)
+      real(8) hin(N), hout(Nout)
+      real(8) coeffremap(N,3),zwork(N,3)
+      real(8) zwork2(N+1,3)
+      
+      real(8) my_zwork(0:N,3)
+      real(8) my_zwork2(0:N,3)
+      
+      integer k
+      double precision, parameter :: dsmll=1.0d-8  
+      real(8) q,q01,q02,q001,q002,q0
+      real(8) z_win(1:N+1), z_wout(1:Nout+1)
+      real(8),parameter :: dpthin = 1.D-3
+      integer :: k1, kbox, ktop, ka, kbot
+      real(8) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
+      real(8) :: p
+      real(8) :: qtri(0:N)
+
+      z_win(1)=0.; z_wout(1)= 0.
+      do k=1,N
+       z_win(k+1)=z_win(k)+hin(k)
+      enddo 
+      
+      do k=1,Nout
+       z_wout(k+1)=z_wout(k)+hout(k)       
+      enddo       
+
+        do k=2,N
+          zwork(k,1)=1./(hin(k-1)+hin(k))
+        enddo
+        
+        do k=2,N-1
+          q0 = 1./(hin(k-1)+hin(k)+hin(k+1))
+          zwork(k,2)=hin(k-1)+2.*hin(k)+hin(k+1)
+          zwork(k,3)=q0
+        enddo       
+      
+        do k= 2,N
+        zwork2(k,1)=zwork(k,1)*(tabin(k)-tabin(k-1))
+        enddo
+        
+        coeffremap(:,1) = tabin(:)
+ 
+         do k=2,N-1
+        q001 = hin(k)*zwork2(k+1,1)
+        q002 = hin(k)*zwork2(k,1)        
+ !       if (q001*q002 < 0) then
+ !         q001 = 0.
+ !         q002 = 0.
+ !       endif
+        q=zwork(k,2)
+        q01=q*zwork2(k+1,1)
+        q02=q*zwork2(k,1)
+ !       if (abs(q001) > abs(q02)) q001 = q02
+ !       if (abs(q002) > abs(q01)) q002 = q01
+        
+        q=(q001-q002)*zwork(k,3)
+        q001=q001-q*hin(k+1)
+        q002=q002+q*hin(k-1)
+        
+        coeffremap(k,3)=coeffremap(k,1)+q001
+        coeffremap(k,2)=coeffremap(k,1)-q002
+        
+        zwork2(k,1)=(2.*q001-q002)**2
+        zwork2(k,2)=(2.*q002-q001)**2
+        enddo
+        
+        do k=1,N
+        if     (k.eq.1 .or. k.eq.N .or.   hin(k).le.dpthin) then
+        coeffremap(k,3) = coeffremap(k,1)
+        coeffremap(k,2) = coeffremap(k,1)
+        zwork2(k,1) = 0.
+        zwork2(k,2) = 0.
+        endif
+        enddo
+        
+        do k=2,N
+        q002=max(zwork2(k-1,2),dsmll)
+        q001=max(zwork2(k,1),dsmll)
+        zwork2(k,3)=(q001*coeffremap(k-1,3)+q002*coeffremap(k,2))/(q001+q002)
+        enddo
+        
+        zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)
+        zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)
+ 
+        do k=1,N
+        q01=zwork2(k+1,3)-coeffremap(k,1)
+        q02=coeffremap(k,1)-zwork2(k,3)
+!        q001=2.*q01
+!        q002=2.*q02
+!        if (q01*q02<0) then
+!          q01=0.
+!          q02=0.
+!        elseif (abs(q01)>abs(q002)) then
+!          q01=q002
+!        elseif (abs(q02)>abs(q001)) then
+!          q02=q001
+!        endif
+        coeffremap(k,2)=coeffremap(k,1)-q02
+        coeffremap(k,3)=coeffremap(k,1)+q01
+        enddo
+        
+        
+        do k=0,N
+        if (k==0) then
+        my_zwork(k,1)=0.
+        my_zwork(k,2)=1.
+        my_zwork(k,3)=0.5
+        my_zwork2(k,1)=1.5*tabin(1)
+        elseif (k==N) then
+        my_zwork(k,1)=0.5
+        my_zwork(k,2)=1.
+        my_zwork(k,3)=0.
+        my_zwork2(k,1)=1.5*tabin(k)        
+        else
+        my_zwork(k,1)=hin(k+1)
+        my_zwork(k,2)=2.*(hin(k)+hin(k+1))
+        my_zwork(k,3)=hin(k)
+        my_zwork2(k,1)=3.*(hin(k+1)*tabin(k)+hin(k)*tabin(k+1))
+        my_zwork2(k,2)=my_zwork2(k,1)
+        endif
+        enddo
+        
+        qtri(0)=-my_zwork(0,3)/my_zwork(0,2)
+        my_zwork2(0,1)=my_zwork2(0,1)/my_zwork(0,2)
+        
+        do k=1,N
+        p=1.0/(my_zwork(k,2)+my_zwork(k,1)*qtri(k-1))
+        qtri(k)=-my_zwork(k,3)*p
+        my_zwork2(k,1)=(my_zwork2(k,1)-my_zwork(k,1)*my_zwork2(k-1,1))*p
+        enddo
+        
+        do k=N-1,0,-1
+        my_zwork2(k,1)=my_zwork2(k,1)+qtri(k)*my_zwork2(k+1,1)
+        enddo
+        
+        do k=1,N
+        coeffremap(k,2)=my_zwork2(k-1,1)
+        coeffremap(k,3)=my_zwork2(k,1)
+        enddo
+        
+        do k=2,N-1
+!        print *,'VAL22 = ',my_zwork(k,1)*my_zwork2(k-1,1)
+!     &+my_zwork(k,2)*my_zwork2(k,1)
+!     &  +my_zwork(k,3)*my_zwork2(k+1,1),my_zwork2(k,2)
+        enddo
+
+      zbot=0.0
+      kbot=1
+      do k=1,Nout
+        ztop=zbot  !top is bottom of previous layer
+        ktop=kbot
+        if     (ztop.ge.z_win(ktop+1)) then
+          ktop=ktop+1
+        endif
+        
+        zbot=z_wout(k+1)
+        zthk=zbot-ztop
+
+        if     (zthk.gt.dpthin .and. ztop.lt.z_wout(Nout+1)) then
+
+          kbot=ktop
+          do while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)
+            kbot=kbot+1
+          enddo
+          zbox=zbot
+          do k1= k+1,Nout
+            if     (z_wout(k1+1)-z_wout(k1).gt.dpthin) then
+              exit !thick layer
+            else
+              zbox=z_wout(k1+1)  !include thin adjacent layers
+              if     (zbox.eq.z_wout(Nout+1)) then
+                exit !at bottom
+              endif
+            endif
+          enddo
+          zthk=zbox-ztop
+
+          kbox=ktop
+          do while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)
+            kbox=kbox+1
+          enddo
+          
+          if     (ktop.eq.kbox) then
+
+
+            if     (z_wout(k)  .ne.z_win(kbox).or.z_wout(k+1).ne.z_win(kbox+1)     ) then
+
+              if     (hin(kbox).gt.dpthin) then
+                q001 = (zbox-z_win(kbox))/hin(kbox)
+                q002 = (ztop-z_win(kbox))/hin(kbox)
+                q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)
+                q02=q01-1.+(q001+q002)
+                q0=1.-q01-q02
+              else
+                q0 = 1.0
+                q01 = 0.
+                q02 = 0.
+              endif
+          tabout(k)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2) &
+       +q02*coeffremap(kbox,3)
+            else
+            tabout(k) = tabin(kbox)
+              
+            endif 
+
+          else
+
+            if     (ktop.le.k .and. kbox.ge.k) then
+              ka = k
+            elseif (kbox-ktop.ge.3) then
+              ka = (kbox+ktop)/2
+            elseif (hin(ktop).ge.hin(kbox)) then
+              ka = ktop
+            else
+              ka = kbox
+            endif !choose ka
+
+            offset=coeffremap(ka,1)
+
+            qtop = z_win(ktop+1)-ztop !partial layer thickness
+            if     (hin(ktop).gt.dpthin) then
+              q=(ztop-z_win(ktop))/hin(ktop)
+              q01=q*(q-1.)
+              q02=q01+q
+              q0=1-q01-q02            
+            else
+              q0 = 1.
+              q01 = 0.
+              q02 = 0.
+            endif
+            
+            tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+ &
+       q02*coeffremap(ktop,3))-offset)*qtop
+
+            do k1= ktop+1,kbox-1
+              tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)
+            enddo !k1
+
+            
+            qbot = zbox-z_win(kbox) !partial layer thickness
+            if     (hin(kbox).gt.dpthin) then
+              q=qbot/hin(kbox)
+              q01=(q-1.)**2
+              q02=q01-1.+q
+              q0=1-q01-q02                            
+            else
+              q0 = 1.0
+              q01 = 0.
+              q02 = 0.
+            endif
+           
+        tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+ &
+       q02*coeffremap(kbox,3))-offset)*qbot
+            
+            rpsum=1.0d0/zthk
+              tabout(k)=offset+tsum*rpsum
+              
+          endif !single or multiple layers
+        else
+        if (k==1) then
+        print *,'problem = ',zthk,z_wout(k+1),hout(1),hin(1),hin(2)
+        stop
+        endif
+         tabout(k) = tabout(k-1)
+          
+        endif !normal:thin layer
+      enddo !k
+            
+      return
+!---------------------------------------------------------------------------------------------------
+      end subroutine reconstructandremap_ps
+!===================================================================================================
+!
+
+end module module_interp
Index: /utils/tools/ABL_TOOLS/module_io.F90
===================================================================
--- /utils/tools/ABL_TOOLS/module_io.F90	(revision 11589)
+++ /utils/tools/ABL_TOOLS/module_io.F90	(revision 11589)
@@ -0,0 +1,823 @@
+MODULE module_io
+   !!======================================================================
+   !!                   ***  MODULE  module_io  ***
+   !! ABL preprocessing tool netcdf I/O subroutines
+   !!=====================================================================
+   !! History : 2016-10  (F. Lemarié)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!
+   !!   FUNCTIONS   : Var_Existence    
+   !!   SUBROUTINES : Read_Ncdf_dim, Write_Ncdf_dim, Read_Ncdf_var1d_Real,
+   !!                 Read_Ncdf_var4d_Real_nt,       Read_Ncdf_var4d_Real_t,    
+   !!                 Write_Ncdf_var1d_Real       , Write_Ncdf_var2d_Real  ,
+   !!                 Write_Ncdf_var4d_Real_t     , Write_Ncdf_var1d_Real_t,
+   !!                 Duplicate_lon_lat_time
+   !!----------------------------------------------------------------------
+   USE netcdf
+   IMPLICIT NONE
+    
+   INTERFACE Read_Ncdf_var
+          MODULE PROCEDURE Read_Ncdf_var1d_Real   ,       &
+             &             Read_Ncdf_var2d_Real   ,       &
+             &             Read_Ncdf_var2d_Real_t ,       &
+             &             Read_Ncdf_var3d_Real_t ,       &
+             &             Read_Ncdf_var4d_Real_t ,       &
+             &             Read_Ncdf_var2d_Real_nt,       &
+             &             Read_Ncdf_var4d_Real_nt
+   END INTERFACE       
+!
+   INTERFACE Write_Ncdf_var
+          MODULE PROCEDURE Write_Ncdf_var1d_Real  ,       &
+             &             Write_Ncdf_var2d_Real  ,       &          
+             &             Write_Ncdf_var1d_Real_t,       &
+             &             Write_Ncdf_var2d_Real_t,       &
+             &             Write_Ncdf_var4d_Real_t
+   END INTERFACE      
+!      
+CONTAINS
+
+
+
+
+
+   LOGICAL FUNCTION Var_Existence( varname , filename )
+      !!---------------------------------------------------------------------
+      !!                    ***  FUNCTION Var_Existence  ***
+      !!                   
+      !! ** Purpose : check if variables varname exists in filename file  
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(*),   intent(in) :: varname,filename
+      INTEGER                    :: status,ncid,varid        
+      
+      status = nf90_open(TRIM(filename),NF90_NOWRITE,ncid)
+      
+      IF ( status/=nf90_noerr ) THEN    
+            WRITE(*,*) "*** Var_Existence: unable to open netcdf file : ",TRIM(filename)
+            stop
+      END IF 
+      
+      status = nf90_inq_varid( ncid, varname, varid )
+      
+      IF ( status/=nf90_noerr ) THEN
+         Var_Existence = .false.
+      ELSE
+         Var_Existence = .true.
+      END IF   
+
+   END FUNCTION Var_Existence
+   
+   
+   
+   SUBROUTINE Read_Ncdf_dim ( dimname, file, dimval )         
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Read_Ncdf_dim  ***
+      !!                   
+      !! ** Purpose : read the integer dimension dimname in input file and   
+      !!              store it in dimval
+      !!
+      !!----------------------------------------------------------------------
+
+      CHARACTER(*),INTENT(in) :: dimname,file    
+      INTEGER                 :: dimval
+      INTEGER                 :: ncid,status,dimid
+      !     
+      status = nf90_open(file,NF90_NOWRITE,ncid)
+      
+      IF ( status/=nf90_noerr ) THEN    
+         WRITE(*,*)"*** Read_Ncdf_dim: unable to open netcdf file : ",trim(file)
+         STOP
+      END IF     
+      !     
+      status = nf90_inq_dimid        ( ncid, dimname,     dimid  ) 
+      status = nf90_inquire_dimension( ncid, dimid,   len=dimval )
+      status = nf90_close( ncid )
+
+   END SUBROUTINE Read_Ncdf_dim  
+
+
+
+
+
+
+
+   SUBROUTINE Write_Ncdf_dim( dimname, file, dimval )   
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Write_Ncdf_dim  ***
+      !!                   
+      !! ** Purpose : write the dimension dimname in the output file   
+      !!
+      !!----------------------------------------------------------------------               
+      CHARACTER(*), INTENT(in) :: dimname,file    
+      INTEGER                  :: dimval
+      INTEGER                  :: ncid,status,dimid
+      !     
+      status = nf90_open(file,NF90_WRITE,ncid)
+      IF ( status/=nf90_noerr ) THEN    
+         WRITE(*,*)"*** Write_Ncdf_dim: to open netcdf file : ",file
+         STOP
+      END IF          
+      
+      status = nf90_redef(ncid)
+      
+      If( dimval.eq.0 ) THEN
+         status = nf90_def_dim(ncid,dimname,nf90_unlimited,dimid)       
+      ELSE
+         status = nf90_def_dim(ncid,dimname,dimval,dimid)  
+      END If    
+      !!
+      status = nf90_enddef(ncid)
+      status = nf90_close(ncid)
+      !
+   END SUBROUTINE Write_Ncdf_dim   
+
+
+
+
+
+
+   SUBROUTINE Read_Ncdf_var1d_Real( varname, file, tabvar )
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Read_Ncdf_var1d_Real  ***
+      !!                   
+      !! ** Purpose : read the 1D variable varname in the input file   
+      !!              and store it in tabvar
+      !!
+      !!----------------------------------------------------------------------  
+      CHARACTER(*),              INTENT(in)    :: varname,file
+      REAL(8), DIMENSION(:),ALLOCATABLE        :: tabvar
+      INTEGER, DIMENSION(1)                    :: dimID
+      INTEGER                                  :: dim1
+      INTEGER                                  :: status,ncid
+      INTEGER                                  :: varid             
+      !
+      status = nf90_open(file,NF90_NOWRITE,ncid)    
+      IF ( status/=nf90_noerr ) THEN    
+         WRITE(*,*) "*** Read_Ncdf_var1d_Real: unable to open netcdf file : ",file
+         STOP
+      END IF          
+      status = nf90_inq_varid(ncid,varname,varid)
+      status = nf90_inquire_variable(ncid,varid,dimids=dimID)
+      status = nf90_inquire_dimension(ncid,dimID(1),len=dim1)                
+      IF( .not. allocated(tabvar) ) THEN
+         ALLOCATE( tabvar ( dim1 ) )  
+      ELSE
+         IF ( any(shape(tabvar)/=(/dim1/)) ) THEN
+            DEALLOCATE ( tabvar          )   
+            ALLOCATE   ( tabvar ( dim1 ) )   
+            WRITE(*,*) 'Warning change shape of array for ',trim(varname)   
+         END IF      
+      END IF        
+      status = nf90_get_var(ncid,varid,tabvar)     
+      status = nf90_close(ncid)
+      !
+   END SUBROUTINE Read_Ncdf_var1d_Real
+
+
+   SUBROUTINE Read_Ncdf_var2d_Real( varname, file, tabvar )    
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Read_Ncdf_var2d_Real  ***
+      !!                   
+      !! ** Purpose : read the 2D variable varname in the input file   
+      !!              and store it in 2D tabvar
+      !!
+      !!----------------------------------------------------------------------          
+      CHARACTER(*),               INTENT(in)    :: varname, file
+      REAL(8), DIMENSION(:,:),    INTENT(inout) :: tabvar
+      INTEGER, DIMENSION(4)                     :: dimIDS
+      INTEGER                                   :: dim1,dim2
+      INTEGER                                   :: status,ncid
+      INTEGER                                   :: varid             
+      !
+      status = nf90_open(file,NF90_NOWRITE,ncid)      
+      IF (status/=nf90_noerr) THEN    
+         WRITE(*,*)"*** Read_Ncdf_var2d_Real: unable to open netcdf file : ",file
+         STOP
+      END IF         
+      status = nf90_inq_varid        (ncid , varname, varid)       
+      status = nf90_inquire_variable (ncid , varid, dimids=dimIDS)
+      status = nf90_inquire_dimension(ncid , dimIDS(1), len=dim1)
+      status = nf90_inquire_dimension(ncid , dimIDS(2), len=dim2)
+      status = nf90_get_var( ncid, varid, tabvar(:,:),      &
+        &                start = (/1,1/), count=(/dim1,dim2/))                                            
+      IF (status/=nf90_noerr) THEN    
+         WRITE(*,*)"unable to retrieve netcdf variable : ",trim(varname)
+         WRITE(*,*)"in file : ",trim(file)
+         WRITE(*,*) "error code: ", status
+         STOP
+      END IF         
+      status = nf90_close(ncid)
+      !    
+   END SUBROUTINE Read_Ncdf_var2d_Real
+
+   
+   SUBROUTINE Read_Ncdf_var2d_Real_nt( varname, file, tabvar, time, level )    
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Read_Ncdf_var2d_Real_nt  ***
+      !!                   
+      !! ** Purpose : read the 4D variable varname in the input file   
+      !!              for a given time and vertical level, and store it in 2D tabvar
+      !!
+      !!----------------------------------------------------------------------          
+      CHARACTER(*),               INTENT(in)    :: varname, file
+      INTEGER     ,               INTENT(in)    :: time, level
+      !REAL(8), DIMENSION(:,:),    ALLOCATABLE   :: tabvar
+      REAL(8), DIMENSION(:,:),    INTENT(inout) :: tabvar
+      INTEGER, DIMENSION(4)                     :: dimIDS
+      INTEGER                                   :: dim1,dim2
+      INTEGER                                   :: status,ncid
+      INTEGER                                   :: varid             
+      !
+      status = nf90_open(file,NF90_NOWRITE,ncid)      
+      IF (status/=nf90_noerr) THEN    
+         WRITE(*,*)"*** Read_Ncdf_var2d_Real_nt: unable to open netcdf file : ",file
+         STOP
+      END IF         
+      status = nf90_inq_varid        (ncid , varname, varid)       
+      status = nf90_inquire_variable (ncid , varid, dimids=dimIDS)
+      status = nf90_inquire_dimension(ncid , dimIDS(1), len=dim1)
+      status = nf90_inquire_dimension(ncid , dimIDS(2), len=dim2)
+      !IF( .not. allocated( tabvar ) ) then
+      !   ALLOCATE ( tabvar( dim1, dim2 ) )  
+      !ELSE
+      !   IF ( (size(tabvar,1) /= dim1) .OR. (size(tabvar,2) /= dim2) ) THEN
+      !      DEALLOCATE( tabvar )
+      !      ALLOCATE  ( tabvar (dim1, dim2 ) )
+      !   END IF     
+      !END IF       
+      status = nf90_get_var( ncid, varid, tabvar(:,:),      &
+        &                start = (/1,1,level,time/), count=(/dim1,dim2,1,1/))                                            
+      IF (status/=nf90_noerr) THEN    
+         WRITE(*,*)"unable to retrieve netcdf variable : ",trim(varname)
+         WRITE(*,*)"in file : ",trim(file)
+         WRITE(*,*) "error code: ", status
+         STOP
+      END IF         
+      status = nf90_close(ncid)
+      !    
+   END SUBROUTINE Read_Ncdf_var2d_Real_nt
+
+
+ 
+ 
+   SUBROUTINE Read_Ncdf_var4d_Real_nt( varname, file, tabvar, time, level )    
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Read_Ncdf_var4d_Real_nt  ***
+      !!                   
+      !! ** Purpose : read the 4D variable varname in the input file   
+      !!              for a given time and vertical level, and store it in 4D tabvar
+      !!
+      !!----------------------------------------------------------------------          
+      CHARACTER(*),               INTENT(in)    :: varname,file
+      INTEGER     ,               INTENT(in)    :: time,level
+      !REAL(8), DIMENSION(:,:,:,:),ALLOCATABLE   :: tabvar
+      REAL(8), DIMENSION(:,:,:,:),INTENT(inout) :: tabvar
+      INTEGER, DIMENSION(4)                     :: dimIDS
+      INTEGER                                   :: dim1,dim2
+      INTEGER                                   :: status,ncid
+      INTEGER                                   :: varid             
+      !
+      status = nf90_open(file,NF90_NOWRITE,ncid)      
+      IF (status/=nf90_noerr) THEN    
+         WRITE(*,*)"*** Read_Ncdf_var4d_Real_nt: unable to open netcdf file : ",file
+         STOP
+      END IF         
+      status = nf90_inq_varid        (ncid , varname, varid)       
+      status = nf90_inquire_variable (ncid , varid, dimids=dimIDS)
+      status = nf90_inquire_dimension(ncid , dimIDS(1), len=dim1)
+      status = nf90_inquire_dimension(ncid , dimIDS(2), len=dim2)
+      !IF( .not. allocated( tabvar ) ) then
+      !   ALLOCATE ( tabvar( dim1, dim2, 1, 1 ) )  
+      !ELSE
+      !   IF ( (size(tabvar,1) /= dim1) .OR. (size(tabvar,2) /= dim2) ) THEN
+      !      DEALLOCATE( tabvar )
+      !      ALLOCATE  ( tabvar (dim1, dim2, 1, 1 ) )
+      !   END IF     
+      !END IF       
+      status = nf90_get_var( ncid, varid, tabvar(:,:,:,:),      &
+        &                start = (/1,1,level,time/), count=(/dim1,dim2,1,1/))                                            
+      IF (status/=nf90_noerr) THEN    
+         WRITE(*,*)"unable to retrieve netcdf variable : ",trim(varname)
+         WRITE(*,*)"in file : ",trim(file)
+         WRITE(*,*) "error code: ", status
+         STOP
+      END IF         
+      status = nf90_close(ncid)
+      !    
+   END SUBROUTINE Read_Ncdf_var4d_Real_nt
+      
+ 
+ 
+ 
+ 
+ 
+   SUBROUTINE Read_Ncdf_var3d_Real_t( varname, file, tabvar, time)   
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Read_Ncdf_var3d_Real_t  ***
+      !!                   
+      !! ** Purpose : read the 3D variable varname in the input file   
+      !!              for a given time level, and store it in tabvar
+      !!
+      !!----------------------------------------------------------------------           
+      CHARACTER(*),               INTENT(in)    :: varname,file
+      INTEGER     ,               INTENT(in)    :: time
+      !REAL(8), DIMENSION(:,:,:),ALLOCATABLE     :: tabvar
+      REAL(8), DIMENSION(:,:,:), INTENT(inout)  :: tabvar
+      INTEGER, DIMENSION(3)                     :: dimIDS
+      INTEGER                                   :: dim1,dim2
+      INTEGER                                   :: status,ncid
+      INTEGER                                   :: varid             
+      !
+      status = nf90_open(file,NF90_NOWRITE,ncid)      
+      IF (status/=nf90_noerr) then    
+         WRITE(*,*)"*** Read_Ncdf_var3d_Real_t: unable to open netcdf file : ",file
+         STOP
+      END IF
+
+      status = nf90_inq_varid        (ncid, varname,varid)        
+      status = nf90_inquire_variable (ncid, varid,dimids=dimIDS)
+      status = nf90_inquire_dimension(ncid, dimIDS(1),len=dim1)
+      status = nf90_inquire_dimension(ncid, dimIDS(2),len=dim2)
+
+      !IF( .not. allocated(tabvar) ) ALLOCATE(tabvar(dim1,dim2,1))  
+
+      status=nf90_get_var(ncid,varid,tabvar,start=(/1,1,time/))                                             
+      IF ( status/=nf90_noerr ) THEN    
+         WRITE(*,*)"unable to retrieve netcdf variable : ",trim(varname)
+         WRITE(*,*)"in file : ",trim(file)
+         WRITE(*,*) "error code: ", status
+         STOP
+      END IF           
+      status = nf90_close(ncid)
+      !
+   END SUBROUTINE Read_Ncdf_var3d_Real_t        
+ 
+      
+   SUBROUTINE Read_Ncdf_var2d_Real_t(varname,file,tabvar,time)   
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Read_Ncdf_var2d_Real_t  ***
+      !!                   
+      !! ** Purpose : read the 3D variable varname in the input file   
+      !!              for a given time level, and store it in tabvar
+      !!
+      !!----------------------------------------------------------------------           
+      CHARACTER(*),               INTENT(in)    :: varname,file
+      INTEGER     ,               INTENT(in)    :: time
+      !REAL(8), DIMENSION(:,:),ALLOCATABLE       :: tabvar
+      REAL(8), DIMENSION(:,:), INTENT(inout)    :: tabvar
+      INTEGER, DIMENSION(3)                     :: dimIDS
+      INTEGER                                   :: dim1,dim2
+      INTEGER                                   :: status,ncid
+      INTEGER                                   :: varid             
+      !
+      status = nf90_open(file,NF90_NOWRITE,ncid)      
+      IF (status/=nf90_noerr) then    
+         WRITE(*,*)"*** Read_Ncdf_var2d_Real_t: unable to open netcdf file : ",file
+         STOP
+      END IF
+
+      status = nf90_inq_varid        (ncid, varname,varid)        
+      status = nf90_inquire_variable (ncid, varid,dimids=dimIDS)
+      status = nf90_inquire_dimension(ncid, dimIDS(1),len=dim1)
+      status = nf90_inquire_dimension(ncid, dimIDS(2),len=dim2)
+
+      !IF( .not. allocated(tabvar) ) ALLOCATE(tabvar(dim1,dim2))  
+
+      status=nf90_get_var(ncid,varid,tabvar,start=(/1,1,time/))                                             
+      IF ( status/=nf90_noerr ) THEN    
+         WRITE(*,*)"unable to retrieve netcdf variable : ",trim(varname)
+         WRITE(*,*)"in file : ",trim(file)
+         STOP
+      END IF           
+      status = nf90_close(ncid)
+      !
+   END SUBROUTINE Read_Ncdf_var2d_Real_t
+
+
+   SUBROUTINE Read_Ncdf_var4d_Real_t(varname,file,tabvar,time)   
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Read_Ncdf_var4d_Real_t  ***
+      !!                   
+      !! ** Purpose : read the 4D variable varname in the input file   
+      !!              for a given time level, and store it in tabvar
+      !!
+      !!----------------------------------------------------------------------           
+      CHARACTER(*),               INTENT(in)    :: varname,file
+      INTEGER     ,               INTENT(in)    :: time
+      !REAL(8), DIMENSION(:,:,:,:),ALLOCATABLE   :: tabvar
+      REAL(8), DIMENSION(:,:,:,:), INTENT(inout):: tabvar
+      INTEGER, DIMENSION(4)                     :: dimIDS
+      INTEGER                                   :: dim1,dim2,dim3
+      INTEGER                                   :: status,ncid
+      INTEGER                                   :: varid             
+      !
+      status = nf90_open(file,NF90_NOWRITE,ncid)      
+      IF (status/=nf90_noerr) then    
+         WRITE(*,*)"*** Read_Ncdf_var4d_Real_t: unable to open netcdf file : ",file
+         STOP
+      END IF         
+
+      status = nf90_inq_varid        (ncid, varname,varid)        
+      status = nf90_inquire_variable (ncid, varid,dimids=dimIDS)
+      status = nf90_inquire_dimension(ncid, dimIDS(1),len=dim1)
+      status = nf90_inquire_dimension(ncid, dimIDS(2),len=dim2)
+      status = nf90_inquire_dimension(ncid, dimIDS(3),len=dim3)
+
+      !IF( .not. allocated(tabvar) ) ALLOCATE(tabvar(dim1,dim2,dim3,1))
+
+      status=nf90_get_var(ncid,varid,tabvar,start=(/1,1,1,time/))                                             
+      IF ( status/=nf90_noerr ) THEN    
+         WRITE(*,*) "unable to retrieve netcdf variable : ",trim(varname)
+         WRITE(*,*) "in file : ",trim(file)
+         WRITE(*,*) "error code: ", status
+         STOP
+      END IF           
+      status = nf90_close(ncid)
+      !
+   END SUBROUTINE Read_Ncdf_var4d_Real_t        
+   
+
+
+
+
+
+   SUBROUTINE Write_Ncdf_var1d_Real( varname, dimname, file, tabvar, typevar )
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Write_Ncdf_var1d_Real  ***
+      !!                   
+      !! ** Purpose : write the 1D variable varname stored in tabvar 
+      !!              in the output file using typevar type (float or double)             
+      !!
+      !!----------------------------------------------------------------------             
+      CHARACTER(*),               INTENT(in)    :: varname,file,dimname,typevar
+      REAL(8), DIMENSION(:), INTENT(in)         :: tabvar
+      INTEGER                                   :: dimid
+      INTEGER                                   :: status,ncid
+      INTEGER                                   :: varid             
+!
+      status = nf90_open(file,NF90_WRITE,ncid)       
+      IF (status/=nf90_noerr) then    
+         WRITE(*,*)"*** Write_Ncdf_var1d_Real: unable to open netcdf file : ",file
+         STOP
+      END IF         
+      status = nf90_inq_dimid(ncid,dimname, dimid)
+      status = nf90_redef(ncid)
+      SELECT CASE( TRIM(typevar) )
+      CASE('double')
+         status = nf90_def_var(ncid,varname,nf90_double,(/dimid/),varid)
+      CASE('float')
+         status = nf90_def_var(ncid,varname,nf90_float,(/dimid/),varid)      
+      END SELECT
+      status = nf90_enddef ( ncid )
+      status = nf90_put_var( ncid ,varid,tabvar)        
+      status = nf90_close  ( ncid )
+      !
+   END SUBROUTINE Write_Ncdf_var1d_Real 
+      
+      
+      
+
+
+   SUBROUTINE Write_Ncdf_var2d_Real( varname, dimname, file, tabvar, typevar )
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Write_Ncdf_var2d_Real  ***
+      !!                   
+      !! ** Purpose : write the 2D variable varname stored in tabvar 
+      !!              in the output file using typevar type (float or double)             
+      !!
+      !!----------------------------------------------------------------------          
+      CHARACTER(*),               INTENT(in)    :: varname, file, typevar
+      CHARACTER(*), DIMENSION(2), INTENT(in)    :: dimname      
+      REAL(8), DIMENSION(:,:), INTENT(in)       :: tabvar
+      INTEGER                                   :: dimid1, dimid2
+      INTEGER                                   :: status, ncid
+      INTEGER                                   :: varid             
+!
+      status = nf90_open(file,NF90_WRITE,ncid)       
+      IF (status/=nf90_noerr) then    
+         WRITE(*,*)"*** Write_Ncdf_var2d_Real: unable to open netcdf file : ",file
+         WRITE(*,*)"*** Write_Ncdf_var2d_Real: variable : ", varname
+         STOP
+      END IF
+         status = nf90_inq_dimid(ncid,dimname(1), dimid1)
+         status = nf90_inq_dimid(ncid,dimname(2), dimid2)
+         status = nf90_inq_varid(ncid,varname,varid)
+         IF (status /= nf90_noerr) THEN
+            status = nf90_redef(ncid)    
+            SELECT CASE( TRIM(typevar) )
+            CASE('double')
+               status = nf90_def_var(ncid,varname,nf90_double,     &
+               &                    (/dimid1,dimid2/),varid)
+            CASE('float')
+               status = nf90_def_var(ncid,varname,nf90_float,      &
+               &                    (/dimid1,dimid2/),varid)     
+            END SELECT
+            status = nf90_enddef(ncid)
+         END IF
+      status = nf90_put_var(ncid,varid,tabvar)     
+      status = nf90_close(ncid)
+      ! 
+   END SUBROUTINE Write_Ncdf_var2d_Real 
+
+
+
+
+
+
+   SUBROUTINE Write_Ncdf_var4d_Real_t( varname, dimname, file, tabvar, time, typevar )
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Write_Ncdf_var4d_Real_t  ***
+      !!                   
+      !! ** Purpose : write the 4D variable varname stored in tabvar 
+      !!              in the output file at time level time using typevar type (float or double)             
+      !!
+      !!----------------------------------------------------------------------  
+      CHARACTER(*),                INTENT(in)   :: varname,file,typevar
+      CHARACTER(*), DIMENSION(4) , INTENT(in)   :: dimname
+      INTEGER, INTENT(in)                       :: time
+      REAL(8), DIMENSION(:,:,:,:), INTENT(in)   :: tabvar
+      INTEGER                                   :: dimid1,dimid2,dimid3,dimid4
+      INTEGER                                   :: status,ncid
+      INTEGER                                   :: varid             
+!
+      status = nf90_open(file,NF90_WRITE,ncid)       
+      IF ( status/=nf90_noerr ) THEN    
+         WRITE(*,*)"*** Write_Ncdf_var4d_Real_t: unable to open netcdf file : ",file
+         STOP
+      END IF            
+      IF ( time==1 .and. (TRIM(typevar)=='double'  .or. TRIM(typevar)=='float') ) THEN      
+         status = nf90_inq_dimid(ncid,dimname(1), dimid1)
+         status = nf90_inq_dimid(ncid,dimname(2), dimid2)
+         status = nf90_inq_dimid(ncid,dimname(3), dimid3)
+         status = nf90_inq_dimid(ncid,dimname(4), dimid4)
+         status = nf90_redef(ncid)      
+         SELECT CASE(TRIM(typevar))
+         CASE('double')
+            status = nf90_def_var(ncid,TRIM(varname),nf90_double,     &
+            &                     (/dimid1,dimid2,dimid3,dimid4/),varid)  
+         CASE('float')
+            status = nf90_def_var(ncid,TRIM(varname),nf90_float,     &
+            &                     (/dimid1,dimid2,dimid3,dimid4/),varid)    
+         END SELECT
+         status = nf90_enddef(ncid)
+      ELSE
+         status = nf90_inq_varid(ncid, varname, varid)
+      END IF    
+      status = nf90_put_var(ncid,varid,tabvar,start=(/1,1,1,time/))
+      IF (status/=nf90_noerr) THEN    
+         WRITE(*,*)"unable to store variable ",varname, &
+	                                       " in file ",file
+         WRITE(*,*)"erorr code: ", status
+         STOP
+      END IF        
+      status = nf90_close(ncid)
+      !
+   END SUBROUTINE Write_Ncdf_var4d_Real_t           
+
+
+
+
+
+
+
+    
+   SUBROUTINE Write_Ncdf_var1d_Real_t(varname,dimname,file,tabvar,time,typevar)
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Write_Ncdf_var1d_Real_t  ***
+      !!                   
+      !! ** Purpose : write the 1D variable varname stored in tabvar 
+      !!              in the output file at time level time using typevar type (float or double)             
+      !!
+      !!----------------------------------------------------------------------  
+      CHARACTER(*),               INTENT(in)    :: varname,file,typevar
+      CHARACTER(*), DIMENSION(1) ,INTENT(in)    :: dimname
+      INTEGER, INTENT(in)                       :: time
+      REAL(8), DIMENSION(:), INTENT(in)         :: tabvar
+      INTEGER                                   :: dimid1
+      INTEGER                                   :: status,ncid
+      INTEGER                                   :: varid             
+      !
+      status = nf90_open(file,NF90_WRITE,ncid)       
+      IF (status/=nf90_noerr) THEN    
+         WRITE(*,*) "*** Write_Ncdf_var1d_Real_t: unable to open netcdf file : ",file
+         STOP
+      END IF          
+            
+      IF( time==1 ) THEN
+         status = nf90_inq_dimid(ncid,dimname(1), dimid1) 
+         status = nf90_redef(ncid)
+         SELECT CASE ( TRIM(typevar) )
+         CASE ('double')
+            status = nf90_def_var(ncid,varname,nf90_double,     &
+             &                                (/dimid1/),varid) 
+         CASE ('float')
+            status = nf90_def_var(ncid,varname,nf90_float,     &
+             &                                (/dimid1/),varid)    
+         END SELECT      
+         status = nf90_enddef(ncid)
+      ELSE
+         status = nf90_inq_varid(ncid, varname, varid)
+      END IF  
+         
+      status = nf90_put_var(ncid,varid,tabvar,start=(/time/))
+      
+      IF (status/=nf90_noerr) THEN    
+         WRITE(*,*) "unable to store variable ",varname, &
+	                                     " in file ",file
+         STOP
+      END IF        
+      status = nf90_close(ncid)
+      !  
+   END SUBROUTINE Write_Ncdf_var1d_Real_t 
+
+
+
+   SUBROUTINE Write_Ncdf_var2d_Real_t( varname, dimname, file, tabvar, time, typevar)
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Write_Ncdf_var2d_Real_t  ***
+      !!                   
+      !! ** Purpose : write the 1D variable varname stored in tabvar 
+      !!              in the output file at time level time using typevar type (float or double)             
+      !!
+      !!----------------------------------------------------------------------  
+      CHARACTER(*),               INTENT(in)    :: varname, file, typevar
+      CHARACTER(*), DIMENSION(2), INTENT(in)    :: dimname
+      INTEGER, INTENT(in)                       :: time
+      REAL(8), DIMENSION(:,:), INTENT(in)       :: tabvar
+      INTEGER                                   :: dimid1, dimid2, dimidt
+      INTEGER                                   :: status,ncid
+      INTEGER                                   :: varid
+      INTEGER, DIMENSION(2)                     :: dim_size
+      !
+      status = nf90_open( file, NF90_WRITE, ncid)       
+      IF (status/=nf90_noerr) THEN    
+         WRITE(*,*) "*** Write_Ncdf_var2d_Real_t: unable to open netcdf file : ",file
+         STOP
+      END IF          
+
+      IF( time==1 ) THEN
+         status = nf90_inq_dimid( ncid, dimname(1), dimid1) 
+         status = nf90_inq_dimid( ncid, dimname(2), dimid2)
+         status = nf90_inq_dimid( ncid, "time"    , dimidt)
+         !status = nf90_def_dim( ncid, "time", NF90_UNLIMITED, dimidt)
+         status = nf90_redef(ncid)
+         SELECT CASE ( TRIM(typevar) )
+         CASE ('double')
+            status = nf90_def_var( ncid, varname, nf90_double,     &
+             &                                (/dimid1,dimid2,dimidt/), varid) 
+         CASE ('float')
+            status = nf90_def_var( ncid, varname, nf90_float,     &
+             &                                (/dimid1,dimid2,dimidt/), varid)    
+         END SELECT
+         status = nf90_enddef(ncid)
+      ELSE
+         status = nf90_inq_varid( ncid, varname, varid )
+      END IF  
+
+      dim_size = SHAPE(tabvar)
+      status = nf90_put_var( ncid, varid, tabvar(:,:), start=(/1,1,time/), count=(/dim_size(1),dim_size(2),1/) )
+      
+      IF (status/=nf90_noerr) THEN    
+         WRITE(*,*) "unable to store variable ",varname, &
+	                                     " in file ",file
+         WRITE(*,*) "error code: ", status
+         STOP
+      END IF        
+      status = nf90_close(ncid)
+      !  
+   END SUBROUTINE Write_Ncdf_var2d_Real_t 
+
+
+
+   SUBROUTINE add_globatt_real( file, att_name, att_value )
+
+      INTEGER      :: ncid, status
+      CHARACTER(*) :: file, att_name
+      REAL(8)      :: att_value
+
+      status = nf90_open( file, nf90_write, ncid)
+
+      ! Enter define mode so we can add the attribute
+      status = nf90_redef( ncid )
+
+      ! ...  put the range attribute, setting it to eight byte reals...
+      status = nf90_put_att( ncid, NF90_GLOBAL, att_name, att_value )
+
+      ! Leave define mode
+      status = nf90_enddef(ncid)
+
+   END SUBROUTINE add_globatt_real
+
+
+
+   SUBROUTINE Duplicate_lon_lat_time( file_in, file_out )
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Duplicate_lon_lat_time  ***
+      !!                   
+      !! ** Purpose : duplicate the attribute of lon, lat and time variables            
+      !!              from file_in in file_out
+      !!----------------------------------------------------------------------  
+      CHARACTER(*),               INTENT(in)    ::  file_in, file_out
+      INTEGER                                   ::  status
+      INTEGER                                   ::  ncid_in, ncid_out    
+      INTEGER                                   :: varid_in,varid_out 
+      
+      status = nf90_open(file_in ,NF90_NOWRITE,ncid_in )
+      status = nf90_open(file_out,NF90_WRITE  ,ncid_out)
+        
+      status = nf90_inq_varid(ncid_in,'lon',varid_in)
+      status = nf90_inq_varid(ncid_out,'lon',varid_out)
+      status = nf90_redef(ncid_out)     
+      status = nf90_copy_att(ncid_in,varid_in,'standard_name',ncid_out,varid_out)
+      status = nf90_copy_att(ncid_in,varid_in,'long_name',ncid_out,varid_out)
+      status = nf90_copy_att(ncid_in,varid_in,'units',ncid_out,varid_out) 
+      status = nf90_copy_att(ncid_in,varid_in,'axis',ncid_out,varid_out)     
+      status = nf90_enddef(ncid_out)  
+      
+      status = nf90_inq_varid(ncid_in,'lat',varid_in)
+      status = nf90_inq_varid(ncid_out,'lat',varid_out)
+      status = nf90_redef(ncid_out)     
+      status = nf90_copy_att(ncid_in,varid_in,'standard_name',ncid_out,varid_out)
+      status = nf90_copy_att(ncid_in,varid_in,'long_name',ncid_out,varid_out)
+      status = nf90_copy_att(ncid_in,varid_in,'units',ncid_out,varid_out) 
+      status = nf90_copy_att(ncid_in,varid_in,'axis',ncid_out,varid_out)     
+      status = nf90_enddef(ncid_out)        
+      
+      status = nf90_inq_varid(ncid_in,'time',varid_in)
+      status = nf90_inq_varid(ncid_out,'time',varid_out)
+      status = nf90_redef(ncid_out)     
+      status = nf90_copy_att(ncid_in,varid_in,'standard_name',ncid_out,varid_out)
+      status = nf90_copy_att(ncid_in,varid_in,'units',ncid_out,varid_out)
+      status = nf90_copy_att(ncid_in,varid_in,'calendar',ncid_out,varid_out)    
+      status = nf90_copy_att(ncid_in,varid_in,'axis',ncid_out,varid_out)  
+      status = nf90_enddef(ncid_out)        
+      
+      status = nf90_close(ncid_in)
+      status = nf90_close(ncid_out)
+      !   
+   END SUBROUTINE Duplicate_lon_lat_time 
+
+
+
+   SUBROUTINE Duplicate_lev_hyb( file_in, file_out )
+      !!---------------------------------------------------------------------
+      !!                    ***  SUBROUTINE Duplicate_lon_lat_time  ***
+      !!                   
+      !! ** Purpose : duplicate the attribute of lon, lat and time variables            
+      !!              from file_in in file_out
+      !!----------------------------------------------------------------------  
+      CHARACTER(*),               INTENT(in)    ::  file_in, file_out
+      INTEGER                                   ::  status
+      INTEGER                                   ::  ncid_in, ncid_out    
+      INTEGER                                   :: varid_in,varid_out 
+      
+      status = nf90_open(file_in ,NF90_NOWRITE,ncid_in )
+      status = nf90_open(file_out,NF90_WRITE  ,ncid_out)
+        
+      status = nf90_inq_varid(ncid_in,'hyai',varid_in)
+      status = nf90_inq_varid(ncid_out,'hyai',varid_out)
+      status = nf90_redef(ncid_out)     
+      status = nf90_copy_att(ncid_in,varid_in,'long_name',ncid_out,varid_out)
+      status = nf90_copy_att(ncid_in,varid_in,'units',ncid_out,varid_out)     
+      status = nf90_enddef(ncid_out)  
+      
+      status = nf90_inq_varid(ncid_in,'hybi',varid_in)
+      status = nf90_inq_varid(ncid_out,'hybi',varid_out)
+      status = nf90_redef(ncid_out)     
+      status = nf90_copy_att(ncid_in,varid_in,'long_name',ncid_out,varid_out)
+      status = nf90_copy_att(ncid_in,varid_in,'units',ncid_out,varid_out)   
+      status = nf90_enddef(ncid_out)        
+
+      status = nf90_inq_varid(ncid_in,'hyam',varid_in)
+      status = nf90_inq_varid(ncid_out,'hyam',varid_out)
+      status = nf90_redef(ncid_out)     
+      status = nf90_copy_att(ncid_in,varid_in,'long_name',ncid_out,varid_out)
+      status = nf90_copy_att(ncid_in,varid_in,'units',ncid_out,varid_out)     
+      status = nf90_enddef(ncid_out)  
+      
+      status = nf90_inq_varid(ncid_in,'hybm',varid_in)
+      status = nf90_inq_varid(ncid_out,'hybm',varid_out)
+      status = nf90_redef(ncid_out)     
+      status = nf90_copy_att(ncid_in,varid_in,'long_name',ncid_out,varid_out)
+      status = nf90_copy_att(ncid_in,varid_in,'units',ncid_out,varid_out)   
+      status = nf90_enddef(ncid_out)  
+      
+      status = nf90_inq_varid(ncid_in,'lev',varid_in)
+      status = nf90_inq_varid(ncid_out,'lev',varid_out)
+      status = nf90_redef(ncid_out)     
+      status = nf90_copy_att(ncid_in,varid_in,'standard_name',ncid_out,varid_out)
+      status = nf90_copy_att(ncid_in,varid_in,'long_name',ncid_out,varid_out)      
+      status = nf90_copy_att(ncid_in,varid_in,'formula',ncid_out,varid_out)    
+      status = nf90_copy_att(ncid_in,varid_in,'formula_terms',ncid_out,varid_out)                
+      status = nf90_copy_att(ncid_in,varid_in,'units',ncid_out,varid_out)
+      status = nf90_copy_att(ncid_in,varid_in,'positive',ncid_out,varid_out)    
+      status = nf90_enddef(ncid_out)        
+      
+      status = nf90_close(ncid_in)
+      status = nf90_close(ncid_out)
+      !   
+   END SUBROUTINE Duplicate_lev_hyb
+
+END MODULE module_io
Index: /utils/tools/ABL_TOOLS/namelist_abl_tools
===================================================================
--- /utils/tools/ABL_TOOLS/namelist_abl_tools	(revision 11589)
+++ /utils/tools/ABL_TOOLS/namelist_abl_tools	(revision 11589)
@@ -0,0 +1,53 @@
+:
+:-----------------------------------------------------------------------------
+: Atmospheric Boundary Layer preprocessing tool 
+:-----------------------------------------------------------------------------
+:
+&nml_dom
+    jpka         =      50,
+    hmax         =   2000.,
+    theta_s      =      2.,
+    hc           =    100.,
+    ln_impose_z1 =  .true.,
+    z1           =     10.,
+/
+
+&nml_opt
+    ptemp_method    =        3  ,  ! potential temperature computation method
+    ln_slp_smth     =   .true.  ,  ! smooth slp and ghw at high latitudes only
+    ln_drw_smth     =   .false. ,  ! smooth after drowning
+    ln_slp_log      =   .false. ,  ! read log(slp)
+    ln_read_zsurf   =   .false. ,  ! read surface geopotential
+    ln_hpg_frc      =   .true.  ,  ! compute horizontal pressure gradient
+    ln_geo_wnd      =   .false. ,  ! compute geostrophic wind
+    ln_c1d          =   .false. ,  ! 1D case
+    ln_read_mask    =   .true.  ,  ! read mask file
+    ln_lsm_land     =   .false. ,  ! inverse land & sea masks
+    ln_perio_latbc  =   .true.  ,  ! periodic lateral boundary conditions
+/
+
+&nml_fld
+    cn_dir      = '',
+    mask_var    = 'LSM',
+    file_m      = 'MASK.nc',
+    file_u      = 'U3D.nc',
+    file_v      = 'V3D.nc',
+    file_t      = 'T3D.nc',
+    file_q      = 'Q3D.nc',
+    file_p      = 'P2D.nc',
+    file_z      = 'Z2D.nc',
+    file_geos   = 'UVG_OUT.nc',
+    file_hpg    = 'HPG_OUT.nc',
+/
+
+&nml_out
+    grd_file  = 'dom_cfg_abl_L50Z10.nc',
+    abl_file  = 'ABL_L50Z10_OUT.nc',
+    drwn_file = 'ABL_DRWN_L50Z10_OUT.nc',
+    var_name  = '',
+/
+
+&nml_c1d
+    iloc = 283,
+    jloc =  52,
+/
