 MODULE sflxread

   !!======================================================================
   !!                   ***  MODULE  sflxread  ***
   !! Atmospherical fluxes :  read of data and fluxes
   !!
   !!   History :
   !!        !  2011-01 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
   !!
   !!======================================================================

#include "toolcpp.h"

   !! * Modules used
   USE parameters, ONLY : rsh,rlg,imin,imax,jmin,jmax,limin,limax, &
                          ljmin,ljmax,limaxp1,ljmaxp1,pi,valmanq
   USE ionc4, ONLY : ionc4_openr,ionc4_corres,ionc4_flimax,ionc4_fljmax, &
                     ionc4_read_lat,ionc4_read_lon,ionc4_read_dimt,ionc4_read_time,&
                     ionc4_read_subxyt,ionc4_read_xy,ionc4_read_var_dim
   USE comvars2d,  ONLY : h0,hx,hy
   USE sflxatm,    ONLY : rad,senslos,dsenslos,irlos,dirlos,latlos,dlatlos
#ifdef key_sflx_analyse
   USE sflxatm,    ONLY : windx,windy
#endif
#ifdef key_grid_curv
  USE parameters, ONLY : valmanq_curv
#endif

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC sflx_read                                 ! function called by step

   !! * Shared module variables
#ifdef key_sflx_analyse
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:), PUBLIC ::                  &
                                    windx1,     &   ! wind at u-location and time 1
                                    windx2,     &   ! wind at u-location and time 2
                                    windy1,     &   ! wind at v-location and time 1
                                    windy2          ! wind at v-location and time 2
#endif

   !! * Private variables
   REAL(kind=rlg),DIMENSION(:,:),ALLOCATABLE :: lonmet2d,latmet2d
   REAL(kind=rlg),DIMENSION(:,:),ALLOCATABLE :: lonmet2dtmp,latmet2dtmp
#ifdef key_sflx_interp_scrip
   LOGICAL,DIMENSION(:,:),ALLOCATABLE        :: maskdst,masksrc
   REAL(kind=rlg),DIMENSION(:,:),ALLOCATABLE :: matrix,matrix_u,matrix_v
   REAL(kind=rlg),DIMENSION(:,:),ALLOCATABLE :: tabvar2d,tabinterp2d
   INTEGER,DIMENSION(:),ALLOCATABLE          :: src_add,dst_add,src_add_u,dst_add_u,src_add_v,dst_add_v
   INTEGER :: nxfin,nyfin
#endif
   INTEGER           :: ixmet,jxmet

 CONTAINS

   !!======================================================================

  SUBROUTINE sflx_read

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE sflx_read  ***
   !&E
   !&E ** Purpose : Read met. data generated from met. models analyses.
   !&E              These data are provided in a NetCDF file with an AO structure :
   !&E              the grid is oriented from South to North and from West to East
   !&E              and spatial steps are constant.
   !&E
   !&E ** Description : Read the following variables :
   !&E              air temperature at 2m, 
   !&E              relative humidity at 2m, 
   !&E              2 wind components at 10m (in m/s),
   !&E              air pressure at 2m.
   !&E              These data are subsequently interpolated onto the computation grid.
   !&E
   !&E ** Called by : step
   !&E
   !&E ** External calls : tool_interp2d_sflx, ionc4_openr, ionc4_corres, ionc4_flimax
   !&E                     ionc4_fljmax, ionc4_read_lat, ionc4_read_lon, ionc4_read_dimt
   !&E                     ionc4_read_time, ionc4_read_subxyt
   !&E
   !&E ** Used ij-arrays : 
   !&E
   !&E ** Modified variables : pa, sat, rh, tosx, tosy, ventvent
   !&E                         pa1, sat1, rh1, tosx1, tosy1, vv1
   !&E                in/out : pa2, sat2, rh2, tosx2, tosy2, vv2
   !&E                         xdpxe, ydpxe, i1xe, j1xe, xdpu, ydpu
   !&E                         i1u, j1u, xdpv, ydpv, i1v, j1v, ilec, ilecmem
   !&E                         tnc1, tnc2
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  1987-02  (J.-C. Salomon, P. Lazure)  Original code
   !&E       !  1990-10  (F. Dumas ) new formulation
   !&E       !  2004-08-19
   !&E       !  2008-07  (T. Odaka) MPI
   !&E       !  2009-04  (V. Garnier) MPI2D
   !&E       !  2009-11  (S. Michel, HOCER) Check meteorological fields units,
   !&E                                      correct the precipitation unit,
   !&E                                      deal with increasing or decreasing lon/lat in met.model
   !&E       !  2011-01  (M. Caillaud, V. Garnier) Integration of SCRIP interpolation
   !&E       !  01-2011 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
   !&E       !  2012-02  (V. Garnier) Reading of the correct record done by the MASTER only
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters,ONLY : rsh,imin,imax,jmin,jmax,ljmin,ljmax,ljmaxp2,limaxp2,  &
                         riosh,rlg,kmax,tref,epsilon,ljminp1,liminp1,liminm2,ljminm2
   USE comvars2d, ONLY : t,nbouc,iscreenlog,ierrorlog,iwarnlog, &
                         lon2d,lat2d,lon2du,lon2dv,lat2du,lat2dv,  &
                         lonmin,latmin,lonmax,latmax
#if defined key_grid_rotated || defined key_sflx_interp_scrip
   USE comvars2d, ONLY : grid_angle
#endif
#ifdef key_grid_curv
   USE comvars2d, ONLY : grid_angle2D
#endif
   USE cometdyn,  ONLY : xdpxe,ydpxe,i1xe,j1xe,xdpu,ydpu,i1u,j1u,xdpv,ydpv
   USE cometdyn,  ONLY : i1v,j1v,idimt,ilec,tnc1,tnc2,ilecmem,cdf_id
   USE cometdyn,  ONLY : name_sflx_rh,name_sflx_pa,name_sflx_rain,name_sflx_ir,name_sflx_windy,name_sflx_fcc
   USE cometdyn,  ONLY : name_sflx_windx,name_sflx_solar,name_sflx_sat,file_meteo
   USE cometdyn,  ONLY : longitude_name_meteo,latitude_name_meteo
   USE cometdyn,  ONLY : imetmin,imetmax,jmetmin,jmetmax,pa1,pa2,rad2,rad1
   USE cometdyn,  ONLY : rain2,rain1,tosx1,tosx2,irlos1,irlos2
   USE cometdyn,  ONLY : l_sflx_rpa,l_sflx_rrh,l_sflx_rrain,l_sflx_rfcc,l_sflx_rwind,l_sflx_rir,l_sflx_rsolar
   USE cometdyn,  ONLY : l_sflx_rsat,l_sflx_rflux,tosy1,tosy2,pa,tosx,tosy
   USE cometemp,  ONLY : rh1,rh2,sat1,sat2,fcc1,fcc2,vv1,vv2,sat,rh,fcc
   USE cometemp,  ONLY : ventvent
   USE comfluxdo, ONLY : rain
   USE_MPI toolmpi, ONLY : ex_i_rsh,ex_j_rsh,MASTER,BCAST_MPI
#ifdef key_sflx_interp_scrip
   USE bicubic_interp,ONLY : get_remap_bicub,make_bicubic_remap
#endif

   IMPLICIT NONE

   !! * Arguments

   !! * Local declarations
   LOGICAL           :: l_pb
   CHARACTER(LEN=19) :: tool_sectodat,tdate
   INTEGER           :: i,j,ncerror,IERR_MPI,dimlon,imetmin_ref,imetmax_ref,jmetmin_ref,jmetmax_ref
   REAL(kind=rsh)    :: dgmet,dfimet,xm,ym,tdif
   REAL(kind=rsh)    :: tabmin,tabmax
   REAL(kind=rlg)    :: tnc0
   REAL(kind=rlg),  ALLOCATABLE,DIMENSION(:)   :: lonmet,latmet
   REAL(kind=riosh),ALLOCATABLE,DIMENSION(:,:) :: uvent1,vvent1,tab1,tosi,tosi2,tab2,uvent2,vvent2
   REAL(kind=riosh),ALLOCATABLE,DIMENSION(:,:) :: tosi_2,tosi2_2
   REAL(kind=rlg)    :: torigin,tmet,dt1,dt2
#ifdef key_grid_rotated
   REAL(kind=riosh),ALLOCATABLE,DIMENSION(:,:) :: uvent1new,vvent1new,uvent2new,vvent2new
#endif
#if defined key_grid_curv
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:) :: tosx1tmp,tosy1tmp,tosx2tmp,tosy2tmp
#endif
   REAL(kind=rsh) :: radeg
   REAL(kind=rlg)                              :: tool_datosec

   !!----------------------------------------------------------------------
   !! * Executable part

   PRINT_DBG*, 'enter SFLX_READ routine, nbouc=',nbouc

!$OMP PARALLEL DEFAULT(none) PRIVATE (i,j,dt1,dt2,tdif,dgmet,dfimet,   &
!$OMP& IERR_MPI,l_pb,tdate,dimlon,                                     &
!$OMP& imetmin_ref,imetmax_ref,jmetmin_ref,jmetmax_ref,                &
!$OMP& xm,ym,torigin,radeg,ncerror,tabmin,tabmax,tnc0)                 &
!$OMP& SHARED(t,nbouc,tmet,h0,iscreenlog,ierrorlog,iwarnlog,           &
!$OMP& lonmin,latmin,lonmax,latmax,                                    &
#ifdef key_MPI_2D
!$OMP& MASTER,liminm2,ljminm2,ljmaxp2,limaxp2,                         &
#endif
!$OMP& lonmet,latmet,uvent1,vvent1,tab1,tosi,tosi2,tab2,uvent2,vvent2, &
#ifdef key_sflx_analyse
!$OMP& windx1,windy1,windx2,windy2,windx,windy,                        &
#endif
!$OMP& lon2d,lat2d,lon2du,lat2du,lon2dv,lat2dv,                        &
#if defined key_grid_rotated || defined key_sflx_interp_scrip
!$OMP& grid_angle,                                                     &
#endif
#ifdef key_grid_rotated
!$OMP& uvent1new,vvent1new,uvent2new,vvent2new,                        &
#endif
#ifdef key_grid_curv
!$OMP& tosx1tmp,tosy1tmp,tosx2tmp,tosy2tmp,grid_angle2D,               &
#endif
!$OMP& lonmet2d,latmet2d,lonmet2dtmp,latmet2dtmp,                      &
#ifdef key_sflx_interp_scrip
!$OMP& src_add,dst_add,src_add_u,dst_add_u,src_add_v,dst_add_v,        &
!$OMP& nxfin,nyfin,                                                    &
#endif
#ifdef key_sflx_interp_scrip
!$OMP& maskdst,masksrc,matrix,matrix_u,matrix_v,tabvar2d,tabinterp2d,  &
#endif
!$OMP& tosi_2,tosi2_2,                                                 &
!$OMP& ixmet,jxmet,cdf_id,i1v,j1v,idimt,ilec,tnc1,tnc2,ilecmem,        &
!$OMP& xdpxe,ydpxe,i1xe,j1xe,xdpu,ydpu,i1u,j1u,xdpv,ydpv,              &
!$OMP& name_sflx_rh,name_sflx_pa,name_sflx_rain,name_sflx_ir,          &
!$OMP& name_sflx_windy,name_sflx_fcc,                                  &
!$OMP& name_sflx_windx,name_sflx_solar,name_sflx_sat,file_meteo,       &
!$OMP& longitude_name_meteo,latitude_name_meteo,                       &
!$OMP& imetmin,imetmax,jmetmin,jmetmax,pa1,pa2,rad2,rad1,              &
!$OMP& rain2,rain1,tosx1,tosx2,irlos1,irlos2,                          &
!$OMP& l_sflx_rpa,l_sflx_rrh,l_sflx_rrain,l_sflx_rfcc,l_sflx_rwind,    &
!$OMP& l_sflx_rir,l_sflx_rsolar,l_sflx_rflux,                          &
!$OMP& l_sflx_rsat,tosy1,tosy2,pa,tosx,tosy,                           &
!$OMP& limin,limax,ljminp1,ljmax,ljmin,liminp1,                        &
!$OMP& rh1,rh2,sat1,sat2,fcc1,fcc2,vv1,vv2,sat,rh,fcc,                 &
!$OMP& ventvent,rain,rad,irlos)

   l_pb = .false. ! No problem

   radeg = pi/180.0_rsh

   tmet=t
   torigin=tref
   torigin=0.0_rlg
   ! include time lag if time origine is not 01/01/1900 in meteorological file
   !torigin=tool_datosec('01/01/1900 00:00:00')-tool_datosec('01/01/0000 00:00:00')
#ifdef key_menor_chain
   torigin=tool_datosec('01/01/1900 00:00:00')-tool_datosec('01/01/0000 00:00:00')
#endif

   IF(nbouc==1) THEN

!$OMP single
     CALL ionc4_openr(file_meteo,l_out_nc4par= .false.)
     CALL ionc4_corres(file_meteo,cdf_id)
     ! global dimension of the atmospherical file
     ixmet=ionc4_flimax(cdf_id)
     jxmet=ionc4_fljmax(cdf_id)

     ! Get the dimensions of the longitude or latitude variables
     CALL ionc4_read_var_dim(file_meteo,TRIM(longitude_name_meteo),dimlon)

     ALLOCATE(latmet2dtmp(1:ixmet,1:jxmet),lonmet2dtmp(1:ixmet,1:jxmet))
     IF (dimlon == 2) THEN
       CALL ionc4_read_xy(file_meteo,TRIM(latitude_name_meteo),latmet2dtmp(1:ixmet,1:jxmet),1,ixmet,1,jxmet)
       CALL ionc4_read_xy(file_meteo,TRIM(longitude_name_meteo),lonmet2dtmp(1:ixmet,1:jxmet),1,ixmet,1,jxmet)
     ELSE
       ALLOCATE(latmet(1:jxmet),lonmet(1:ixmet))
       CALL ionc4_read_lat(file_meteo,latmet(1:jxmet),1,jxmet)
       CALL ionc4_read_lon(file_meteo,lonmet(1:ixmet),1,ixmet)
       DO j=1,jxmet
         lonmet2dtmp(1:ixmet,j)=lonmet(1:ixmet)
       ENDDO
       DO i=1,ixmet
         latmet2dtmp(i,1:jxmet)=latmet(1:jxmet)
       ENDDO
     END IF

     dgmet =ABS(lonmet2dtmp(2,1)-lonmet2dtmp(1,1))
     dfimet=ABS(latmet2dtmp(1,2)-latmet2dtmp(1,1))


#ifdef key_grid_rotated
     ! to estimate min and max indexes relative to the rotated grid
     ! to be used in atmospherical file
     ! Work indifferently if met. model longitudes decrease or increase.
     imetmin=1               ! global dimensions of the meteorological file by default
     imetmax=ixmet
     jmetmin=1
     jmetmax=jxmet
     DO i=1,ixmet-1
       IF((lonmin-lonmet(i))*(lonmin-lonmet(i+1)) <= 0.0_rlg) imetmin=i
       IF((lonmax-lonmet(i)+epsilon)*(lonmax-lonmet(i+1)-epsilon) <= 0.0_rlg) imetmax=i+1
     ENDDO
     ! Work indifferently if met. model latitudes decrease or increase.
     DO j=1,jxmet-1
       IF((latmin-latmet(j))*(latmin-latmet(j+1)) <= 0.0_rlg) jmetmin=j
       IF((latmax-latmet(j)+epsilon)*(latmax-latmet(j+1)-epsilon) <= 0.0_rlg) jmetmax=j+1
     ENDDO

#else

!toseelater (if pb) cval warning before it was xmin=REAL(imin+1,rsh) and ymin=REAL(jmin+1,rsh)
     ! lon2d_g(liminp1,ljminp1) ! p1 because the first rowx and column (imin,jmin)
     ! llatd_g(liminp1,ljminp1) ! are outside the domain of calculation
                                ! p1 because we interpolate on u and v-locations (i+0.5 and j+0.5)
     ! lon2d_g(limaxp2,ljmaxp2) ! + p1 because of the definition of ii2 in toolinterp2dsflx.F90
     ! lat2d_g(limaxp2,ljmaxp2) ! total = p2 but maybe too safe

     ! dimensions of the meteorological file usefull for each cpu
     imetmin=ixmet           ! global dimensions of the meteorological file by default
     imetmax=1
     jmetmin=jxmet
     jmetmax=1
     imetmin_ref=imetmin
     imetmax_ref=imetmax
     jmetmin_ref=jmetmin
     jmetmax_ref=jmetmax
     ! check the useful indexes of the meteorological file
     ! Work indifferently if met. model longitudes decrease or increase.
     DO j=1,jxmet
     DO i=1,ixmet-1
       IF((lonmin-lonmet2dtmp(i,j)+epsilon) * (lonmin-lonmet2dtmp(i+1,j)-epsilon) <= 0.0_rlg) THEN
         IF (i<imetmin .AND. MAXVAL(latmet2dtmp(i:ixmet,j))>latmin-dfimet .AND. MINVAL(latmet2dtmp(i:ixmet,j))<latmax+dfimet) THEN
           imetmin=i
           IF (j<jmetmin) jmetmin=j
         END IF
       END IF
     ENDDO
     ENDDO
     DO i=imetmin,ixmet-1
     DO j=jmetmin,jxmet-1
       IF((latmax-latmet2dtmp(i,j)  +epsilon) * (latmax-latmet2dtmp(i,j+1)-epsilon) <= 0.0_rlg) THEN
         IF (j+1>=jmetmax  .AND. MAXVAL(lonmet2dtmp(i+1,jmetmin:j+1))>lonmin-dgmet .AND. &
                                 MINVAL(lonmet2dtmp(i+1,jmetmin:j+1))<lonmax+dgmet ) THEN
           jmetmax=j+1
           IF (i+1>imetmax) imetmax=i+1
         END IF
       END IF
     ENDDO
     ENDDO

#endif

     IF (imetmin==imetmin_ref) THEN
       WRITE(iwarnlog,*) 'None index has been found for imetmin'
       WRITE(iwarnlog,*) 'imetmin =',imetmin,' imetmax =',imetmax
       WRITE(iwarnlog,*) 'jmetmin =',jmetmin,' jmetmax =',jmetmax
       WRITE(iwarnlog,*) 'For this cpu :'
       WRITE(iwarnlog,*) 'lonmin =',lonmin,' lonmax =',lonmax
       WRITE(iwarnlog,*) 'latmin =',latmin,' latmax =',latmax
       WRITE(iwarnlog,*) 'You use too many cpus for a so coarse meteorological grid'
       WRITE(iwarnlog,*) 'or the file is not large enough'
       WRITE(iwarnlog,*) 'imetmin is set to 1'
!       WRITE(ierrorlog,*) 'Simulation stopped.'
!       CALL_MPI MPI_FINALIZE(IERR_MPI)
!       STOP
       imetmin=1
     END IF

     IF (imetmax==imetmax_ref) THEN
       WRITE(iwarnlog,*) 'None index has been found for imetmax'
       WRITE(iwarnlog,*) 'imetmin =',imetmin,' imetmax =',imetmax
       WRITE(iwarnlog,*) 'jmetmin =',jmetmin,' jmetmax =',jmetmax
       WRITE(iwarnlog,*) 'For this cpu :'
       WRITE(iwarnlog,*) 'lonmin =',lonmin,' lonmax =',lonmax
       WRITE(iwarnlog,*) 'latmin =',latmin,' latmax =',latmax
       WRITE(iwarnlog,*) 'You use too many cpus for a so coarse meteorological grid'
       WRITE(iwarnlog,*) 'or the file is not large enough'
       WRITE(iwarnlog,*) 'imetmax is set to ixmet'
!       WRITE(ierrorlog,*) 'Simulation stopped.'
!       CALL_MPI MPI_FINALIZE(IERR_MPI)
!       STOP
       imetmax=ixmet
     END IF

     IF (jmetmin==jmetmin_ref) THEN
       WRITE(iwarnlog,*) 'None index has been found for jmetmin'
       WRITE(iwarnlog,*) 'imetmin =',imetmin,' imetmax =',imetmax
       WRITE(iwarnlog,*) 'jmetmin =',jmetmin,' jmetmax =',jmetmax
       WRITE(iwarnlog,*) 'For this cpu :'
       WRITE(iwarnlog,*) 'lonmin =',lonmin,' lonmax =',lonmax
       WRITE(iwarnlog,*) 'latmin =',latmin,' latmax =',latmax
       WRITE(iwarnlog,*) 'You use too many cpus for a so coarse meteorological grid'
       WRITE(iwarnlog,*) 'or the file is not large enough'
       WRITE(iwarnlog,*) 'jmetmin is set to 1'
!       WRITE(ierrorlog,*) 'Simulation stopped.'
!       CALL_MPI MPI_FINALIZE(IERR_MPI)
!       STOP
       jmetmin=1
     END IF

     IF (jmetmax==jmetmax_ref) THEN
       WRITE(iwarnlog,*) 'None index has been found for jmetmax'
       WRITE(iwarnlog,*) 'imetmin =',imetmin,' imetmax =',imetmax
       WRITE(iwarnlog,*) 'jmetmin =',jmetmin,' jmetmax =',jmetmax
       WRITE(iwarnlog,*) 'For this cpu :'
       WRITE(iwarnlog,*) 'lonmin =',lonmin,' lonmax =',lonmax
       WRITE(iwarnlog,*) 'latmin =',latmin,' latmax =',latmax
       WRITE(iwarnlog,*) 'You use too many cpus for a so coarse meteorological grid'
       WRITE(iwarnlog,*) 'or the file is not large enough'
       WRITE(iwarnlog,*) 'jmetmax is set to jxmet'
!       WRITE(ierrorlog,*) 'Simulation stopped.'
!       CALL_MPI MPI_FINALIZE(IERR_MPI)
!       STOP
       jmetmax=jxmet
     END IF

     IF (imetmin>=imetmax .OR. jmetmin>=jmetmax) THEN
       WRITE(ierrorlog,*) 'The size of the domain read from the meteorological field is not large enough'
       WRITE(ierrorlog,*) 'imetmin =',imetmin,' imetmax =',imetmax
       WRITE(ierrorlog,*) 'jmetmin =',jmetmin,' jmetmax =',jmetmax
       WRITE(ierrorlog,*) 'You use too many cpus for a so coarse meteorological grid'
       WRITE(ierrorlog,*) 'Simulation stopped.'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     END IF

     IF_MPI (MASTER) THEN
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) 'Meteorological file : spatial limits'
     WRITE(iscreenlog,*) 'imetmin = ',imetmin
     WRITE(iscreenlog,*) 'imetmax = ',imetmax
     WRITE(iscreenlog,*) 'jmetmin = ',jmetmin
     WRITE(iscreenlog,*) 'jmetmax = ',jmetmax
     WRITE(iscreenlog,*) ' '
     ENDIF_MPI

#ifdef key_sflx_interp_scrip
     IF (imetmin <= 2) THEN
       WRITE(ierrorlog,*) 'Meteorological file too small'
       WRITE(ierrorlog,*) 'At the Western boundary'
       WRITE(ierrorlog,*) 'Simulation stopped.'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     ENDIF
     IF (jmetmin <= 2) THEN
       WRITE(ierrorlog,*) 'Meteorological file too small'
       WRITE(ierrorlog,*) 'At the Southern boundary'
       WRITE(ierrorlog,*) 'Simulation stopped.'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     ENDIF
     IF (imetmax >= ixmet - 1) THEN
       WRITE(ierrorlog,*) 'Meteorological file too small'
       WRITE(ierrorlog,*) 'At the Eastern boundary'
       WRITE(ierrorlog,*) 'Simulation stopped.'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     ENDIF
     IF (jmetmax >= jxmet - 1) THEN
       WRITE(ierrorlog,*) 'Meteorological file too small'
       WRITE(ierrorlog,*) 'At the Northern boundary'
       WRITE(ierrorlog,*) 'Simulation stopped.'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     ENDIF
#endif

     imetmin = MAX0(1 , imetmin - 2)
     jmetmin = MAX0(1 , jmetmin - 2)
     imetmax = MIN0(ixmet, imetmax + 2)
     jmetmax = MIN0(jxmet, jmetmax + 2)

     IF_MPI (MASTER) THEN
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) 'Meteorological file : spatial limits'
     WRITE(iscreenlog,*) 'imetmin = ',imetmin
     WRITE(iscreenlog,*) 'imetmax = ',imetmax
     WRITE(iscreenlog,*) 'jmetmin = ',jmetmin
     WRITE(iscreenlog,*) 'jmetmax = ',jmetmax
     WRITE(iscreenlog,*) ' '
     ENDIF_MPI

     ! Save the longitude latitude in 2D vraiables
     ALLOCATE(lonmet2d(imetmin:imetmax,jmetmin:jmetmax))
     ALLOCATE(latmet2d(imetmin:imetmax,jmetmin:jmetmax))

     DO j=jmetmin,jmetmax
       lonmet2d(imetmin:imetmax,j)=lonmet2dtmp(imetmin:imetmax,j)
     ENDDO
     DO i=imetmin,imetmax
       latmet2d(i,jmetmin:jmetmax)=latmet2dtmp(i,jmetmin:jmetmax)
     ENDDO

     DO j=jmetmin,jmetmax
     DO i=imetmin,imetmax
       IF(lonmet2d(i,j) >= 180.0_rlg) lonmet2d(i,j)=lonmet2d(i,j)-360.0_rlg
     END DO
     END DO

#ifdef key_grid_rotated
     ALLOCATE(uvent1new(imetmin:imetmax,jmetmin:jmetmax))
     ALLOCATE(vvent1new(imetmin:imetmax,jmetmin:jmetmax))
     ALLOCATE(uvent2new(imetmin:imetmax,jmetmin:jmetmax))
     ALLOCATE(vvent2new(imetmin:imetmax,jmetmin:jmetmax))
#endif
#ifdef key_grid_curv
     ALLOCATE(tosx1tmp(limin:limax,ljmin:ljmax))
     ALLOCATE(tosy1tmp(limin:limax,ljmin:ljmax))
     ALLOCATE(tosx2tmp(limin:limax,ljmin:ljmax))
     ALLOCATE(tosy2tmp(limin:limax,ljmin:ljmax))
#endif
     ALLOCATE( uvent1(imetmin:imetmax,jmetmin:jmetmax))
     ALLOCATE( vvent1(imetmin:imetmax,jmetmin:jmetmax))
     ALLOCATE( uvent2(imetmin:imetmax,jmetmin:jmetmax))
     ALLOCATE( vvent2(imetmin:imetmax,jmetmin:jmetmax))
     ALLOCATE(   tosi(imetmin:imetmax,jmetmin:jmetmax))
     ALLOCATE(  tosi2(imetmin:imetmax,jmetmin:jmetmax))
     ALLOCATE( tosi_2(imetmin:imetmax,jmetmin:jmetmax))
     ALLOCATE(tosi2_2(imetmin:imetmax,jmetmin:jmetmax))
     ALLOCATE(   tab1(imetmin:imetmax,jmetmin:jmetmax))
     ALLOCATE(   tab2(imetmin:imetmax,jmetmin:jmetmax))

#ifdef key_sflx_interp_scrip
     ALLOCATE(masksrc(imetmin:imetmax,jmetmin:jmetmax))
     ALLOCATE(tabvar2d(imetmin:imetmax,jmetmin:jmetmax))

     ALLOCATE(tabinterp2d(limin:limax,ljmin:ljmax))
     ALLOCATE(maskdst(limin:limax,ljmin:ljmax))

     masksrc(:,:)=.true.
     maskdst(:,:)=.true.
     tabvar2d(:,:)=0.0_rlg
#endif

!$OMP end single

   ! -2- Ponderation weights for the spatial interpolation
   ! -----------------------------------------------------

#ifdef key_sflx_interp_scrip
     nxfin=limax-limin+1
     nyfin=ljmax-ljmin+1

!$OMP single
     maskdst(:,:)=.true.
#ifdef key_grid_curv
     where(lon2d(limin:limax,ljmin:ljmax)==valmanq_curv) maskdst=.False.
#endif
     CALL get_remap_bicub(latmet2d,REAL(lat2d(limin:limax,ljmin:ljmax),rlg),lonmet2d,REAL(lon2d(limin:limax,ljmin:ljmax),rlg), &
                          masksrc,maskdst,matrix,src_add,dst_add)
     maskdst(:,:)=.true.
#ifdef key_grid_curv
     where(lon2du==valmanq_curv) maskdst=.False.
#endif
     CALL get_remap_bicub(latmet2d,REAL(lat2du(limin:limax,ljmin:ljmax),rlg),lonmet2d,REAL(lon2du(limin:limax,ljmin:ljmax),rlg), &
                          masksrc,maskdst,matrix_u,src_add_u,dst_add_u)
     maskdst(:,:)=.true.
#ifdef key_grid_curv
     where(lon2dv==valmanq_curv) maskdst=.False.
#endif
     CALL get_remap_bicub(latmet2d,REAL(lat2dv(limin:limax,ljmin:ljmax),rlg),lonmet2d,REAL(lon2dv(limin:limax,ljmin:ljmax),rlg), &
                          masksrc,maskdst,matrix_v,src_add_v,dst_add_v)
!$OMP end single

#else

     dgmet =ABS(lonmet2dtmp(2,1)-lonmet2dtmp(1,1))
     dfimet=ABS(latmet2dtmp(1,2)-latmet2dtmp(1,1))

!$OMP do schedule(static)
     DO j=ljminp1,ljmax
       DO i=liminp1,limax
         xm = (lon2d(i,j)-lonmet2dtmp(1,1))/dgmet+1        !pour les points xe 
         ym = (lat2d(i,j)-latmet2dtmp(1,1))/dfimet+1
         xdpxe(i,j) = xm - REAL(INT(xm))
         ydpxe(i,j) = ym - REAL(INT(ym))
         i1xe(i,j) = INT(xm)
         j1xe(i,j) = INT(ym)

         xm = (lon2du(i,j)-lonmet2dtmp(1,1))/dgmet+1        !pour les points u 
         ym = (lat2du(i,j)-latmet2dtmp(1,1))/dfimet+1
         xdpu(i,j) = xm - REAL(INT(xm))
         ydpu(i,j) = ym - REAL(INT(ym))
         i1u(i,j) = INT(xm)
         j1u(i,j) = INT(ym)

         xm = (lon2dv(i,j)-lonmet2dtmp(1,1))/dgmet+1        !pour les points v 
         ym = (lat2dv(i,j)-latmet2dtmp(1,1))/dfimet+1
         xdpv(i,j) = xm - REAL(INT(xm))
         ydpv(i,j) = ym - REAL(INT(ym))
         i1v(i,j) = INT(xm)
         j1v(i,j) = INT(ym)
       ENDDO
     ENDDO
!$OMP end do

#endif

   ! Time index in the meteorological file
   ! -------------------------------------

     IF_MPI (MASTER) THEN

!$OMP single
     idimt=ionc4_read_dimt(file_meteo)
     ilec=0 
     tnc1=0.0_rlg
     tnc2=0.0_rlg
     DO WHILE((tmet-(tnc1-torigin))*(tmet-(tnc2-torigin)) > 0.0_rlg .AND. ilec < idimt)
       ilec=ilec+1 
       IF (ilec == idimt) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine SFLXREAD.F90'
         WRITE(ierrorlog,*) 'All records are read, dates are not compatible with simulation time.'
         WRITE(ierrorlog,*) 'Make sure your meteorological file is correct: ',TRIM(file_meteo)
         WRITE(ierrorlog,*) 'Make sure the origin time is 01/01/1900 00:00:00 in the met. file.'
         WRITE(ierrorlog,*) 'If not, modify variable "torigin" in the routine SFLXREAD.F90.'
         WRITE(ierrorlog,*) ''
         tdate = tool_sectodat(tmet)
         WRITE(ierrorlog,*) 't in MARS model =',tdate
         tdate = tool_sectodat(tnc1-torigin)
         WRITE(ierrorlog,*) 't1 in atmospherical file =',tdate
         tdate = tool_sectodat(tnc2-torigin)
         WRITE(ierrorlog,*) 't2 in atmospherical file =',tdate
         WRITE(ierrorlog,*) 'Simulation stopped.'
         WRITE(ierrorlog,*) ''
         l_pb = .true.
       END IF
       IF (.NOT. l_pb) THEN
         CALL ionc4_read_time(file_meteo,ilec,tnc1)
         CALL ionc4_read_time(file_meteo,ilec+1,tnc2)
       END IF
     ENDDO
     ilecmem=ilec
!$OMP end single

     ENDIF_MPI

     ! the MASTER sends informations to every cpus
OMPMPI barrier
OMPMPI master
     CALL_MPI BCAST_MPI(l_pb)
     CALL_MPI BCAST_MPI(idimt)
     CALL_MPI BCAST_MPI(ilec)
     CALL_MPI BCAST_MPI(ilecmem)
     CALL_MPI BCAST_MPI(tnc1)
     CALL_MPI BCAST_MPI(tnc2)
OMPMPI end master
OMPMPI barrier
OMPMPI flush(l_pb,idimt,ilec,ilecmem,tnc1,tnc2)

     IF (l_pb) THEN
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     END IF

     ! -3- lecture des tableaux deux indices qui encadrent le pas de temps courant

     IF(l_sflx_rpa) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_pa),tab1,imetmin,imetmax,jmetmin,jmetmax,ilec,1,1)
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_pa),tab2,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       pa1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       tabvar2d(:,:)=tab2(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       pa2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,pa1)
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab2,pa2)
#endif
     END IF

     IF(l_sflx_rrh) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_rh),tab1,imetmin,imetmax,jmetmin,jmetmax,ilec,1,1)
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_rh),tab2,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       rh1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       tabvar2d(:,:)=tab2(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       rh2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,rh1)
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab2,rh2)
#endif
     END IF

     IF(l_sflx_rsat) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_sat),tab1,imetmin,imetmax,jmetmin,jmetmax,ilec,1,1)
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_sat),tab2,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       sat1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       tabvar2d(:,:)=tab2(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       sat2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,sat1)
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab2,sat2)
#endif
     END IF

     IF(l_sflx_rsolar) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_solar),tab1,imetmin,imetmax,jmetmin,jmetmax,ilec,1,1)
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_solar),tab2,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       rad1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       tabvar2d(:,:)=tab2(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       rad2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,rad1)
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab2,rad2)
#endif
     END IF

     IF(l_sflx_rir) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_ir),tab1,imetmin,imetmax,jmetmin,jmetmax,ilec,1,1)
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_ir),tab2,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       irlos1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       tabvar2d(:,:)=tab2(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       irlos2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,irlos1)
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab2,irlos2)
#endif
     END IF

     IF(l_sflx_rrain) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_rain),tab1,imetmin,imetmax,jmetmin,jmetmax,ilec,1,1)
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_rain),tab2,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       rain1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       IF (.NOT. l_sflx_rflux) WHERE( rain1(:,:) < 0.0_rsh ) rain1(:,:) = 0.0_rsh
       tabvar2d(:,:)=tab2(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       rain2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       IF (.NOT. l_sflx_rflux) WHERE( rain2(:,:) < 0.0_rsh ) rain2(:,:) = 0.0_rsh
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,rain1)
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab2,rain2)
#endif
       ! Convert the precipitation amount (cumulated between 2 meteorological steps, in kg.m-2=mm)
       ! into a precipitation flux (per unit of time, i.e. mm.s-1).
!$OMP single
       IF (ilec>1) THEN
         CALL ionc4_read_time(file_meteo,ilec-1,tnc0)
         rain1(:,:)=rain1(:,:)/(tnc1-tnc0)
       ELSE
         rain1(:,:)=rain1(:,:)/(tnc2-tnc1)
       ENDIF
       rain2(:,:)=rain2(:,:)/(tnc2-tnc1)
!$OMP end single
     END IF

     IF(l_sflx_rfcc) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_fcc),tab1,imetmin,imetmax,jmetmin,jmetmax,ilec,1,1)
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_fcc),tab2,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       fcc1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       tabvar2d(:,:)=tab2(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       fcc2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,fcc1)
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab2,fcc2)
#endif
     END IF

     IF(l_sflx_rwind) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_windy),vvent1,imetmin,imetmax,jmetmin,jmetmax,ilec,1,1)
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_windx),uvent1,imetmin,imetmax,jmetmin,jmetmax,ilec,1,1)
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_windy),vvent2,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_windx),uvent2,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
!$OMP do schedule(static)
       DO j= jmetmin,jmetmax
#ifdef key_grid_rotated 
         uvent1new(:,j)=uvent1(:,j)*COS(grid_angle*radeg)-vvent1(:,j)*SIN(grid_angle*radeg)
         vvent1new(:,j)=vvent1(:,j)*COS(grid_angle*radeg)+uvent1(:,j)*SIN(grid_angle*radeg)
         uvent2new(:,j)=uvent2(:,j)*COS(grid_angle*radeg)-vvent2(:,j)*SIN(grid_angle*radeg)
         vvent2new(:,j)=vvent2(:,j)*COS(grid_angle*radeg)+uvent2(:,j)*SIN(grid_angle*radeg)
         uvent1(:,j)=uvent1new(:,j)
         vvent1(:,j)=vvent1new(:,j) 
         uvent2(:,j)=uvent2new(:,j)
         vvent2(:,j)=vvent2new(:,j)
#endif
         tab1(:,j)=SQRT(uvent1(:,j)*uvent1(:,j)+vvent1(:,j)*vvent1(:,j))
         tab2(:,j)=SQRT(uvent2(:,j)*uvent2(:,j)+vvent2(:,j)*vvent2(:,j))
       ENDDO
!$OMP end do
       IF (.NOT. l_sflx_rflux) THEN
!$OMP do schedule(static)
         DO j= jmetmin,jmetmax
           tosi(:,j)=tab1(:,j)*uvent1(:,j)
           tosi2(:,j)=tab1(:,j)*vvent1(:,j)
           tosi_2(:,j)=tab2(:,j)*uvent2(:,j)
           tosi2_2(:,j)=tab2(:,j)*vvent2(:,j)
         ENDDO
!$OMP end do
       ELSE
!$OMP do schedule(static)
         DO j= jmetmin,jmetmax
           tosi(:,j)=uvent1(:,j)
           tosi2(:,j)=vvent1(:,j)
           tosi_2(:,j)=uvent2(:,j)
           tosi2_2(:,j)=vvent2(:,j)
         ENDDO
!$OMP end do
       END IF
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       vv1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       WHERE(vv1(:,:)<0.0_rsh) vv1(:,:)=0.0_rsh
       tabvar2d(:,:)=tab2(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       vv2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       WHERE(vv2(:,:)<0.0_rsh) vv2(:,:)=0.0_rsh
       tabvar2d(:,:)=tosi(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       tosx1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       tabvar2d(:,:)=tosi_2(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       tosx2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       tabvar2d(:,:)=tosi2(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       tosy1(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       tabvar2d(:,:)=tosi2_2(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       tosy2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,vv1)
       CALL tool_interp2d_sflx(i1u,j1u,xdpu,ydpu,tosi,tosx1)
       CALL tool_interp2d_sflx(i1v,j1v,xdpv,ydpv,tosi2,tosy1)
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab2,vv2)
       CALL tool_interp2d_sflx(i1u,j1u,xdpu,ydpu,tosi_2,tosx2)
       CALL tool_interp2d_sflx(i1v,j1v,xdpv,ydpv,tosi2_2,tosy2)
#endif

#ifdef key_grid_curv
!$OMP do schedule(static)
     DO j=ljmin,ljmax
       tosx1tmp(limin:limax,j)=tosx1(limin:limax,j)*COS(grid_angle2D(limin:limax,j)*radeg)+tosy1(limin:limax,j)*SIN(grid_angle2D(limin:limax,j)*radeg)
       tosy1tmp(limin:limax,j)=tosy1(limin:limax,j)*COS(grid_angle2D(limin:limax,j)*radeg)-tosx1(limin:limax,j)*SIN(grid_angle2D(limin:limax,j)*radeg)
       tosx2tmp(limin:limax,j)=tosx2(limin:limax,j)*COS(grid_angle2D(limin:limax,j)*radeg)+tosy2(limin:limax,j)*SIN(grid_angle2D(limin:limax,j)*radeg)
       tosy2tmp(limin:limax,j)=tosy2(limin:limax,j)*COS(grid_angle2D(limin:limax,j)*radeg)-tosx2(limin:limax,j)*SIN(grid_angle2D(limin:limax,j)*radeg)
       tosx1(limin:limax,j)=tosx1tmp(limin:limax,j) 
       tosy1(limin:limax,j)=tosy1tmp(limin:limax,j) 
       tosx2(limin:limax,j)=tosx2tmp(limin:limax,j) 
       tosy2(limin:limax,j)=tosy2tmp(limin:limax,j) 
     ENDDO
!$OMP end do
#endif

#ifdef key_sflx_analyse
!$OMP single
       ALLOCATE( windx1(limin:limax,ljmin:ljmax))
       ALLOCATE( windy1(limin:limax,ljmin:ljmax))
       ALLOCATE( windx2(limin:limax,ljmin:ljmax))
       ALLOCATE( windy2(limin:limax,ljmin:ljmax))
       windx1(:,:)=0.0_rsh ; windx2(:,:)=0.0_rsh
       windy1(:,:)=0.0_rsh ; windy2(:,:)=0.0_rsh
       CALL tool_interp2d_sflx(i1u,j1u,xdpu,ydpu,uvent1,windx1)
       CALL tool_interp2d_sflx(i1u,j1u,xdpu,ydpu,uvent2,windx2)
       CALL tool_interp2d_sflx(i1v,j1v,xdpv,ydpv,vvent1,windy1)
       CALL tool_interp2d_sflx(i1v,j1v,xdpv,ydpv,vvent2,windy2)
!$OMP end single
#endif

!$OMP single
#ifdef key_grid_rotated 
       DEALLOCATE(uvent1new,vvent1new)
#endif
#ifdef key_grid_curv
       DEALLOCATE(tosx1tmp,tosy1tmp,tosx2tmp,tosy2tmp)
#endif
       DEALLOCATE(uvent1,vvent1,tosi,tosi2,uvent2,vvent2,tosi_2,tosi2_2,tab1,tab2)
       IF (dimlon == 1) DEALLOCATE(latmet,lonmet)
       DEALLOCATE(lonmet2d,latmet2d,lonmet2dtmp,latmet2dtmp)
!$OMP end single
     ENDIF

     ! Check meteorological fields units.
!$OMP single
     ncerror=0

     IF(l_sflx_rpa) THEN
       tabmin=1000000.0_rsh
       tabmax=0.0_rsh
       DO j=ljminp1,ljmax
         DO i=liminp1,limax
           IF (h0(i,j)/=-valmanq) THEN
             tabmin=MIN(tabmin,MIN(pa1(i,j),pa2(i,j)))
             tabmax=MAX(tabmax,MAX(pa1(i,j),pa2(i,j)))
           ENDIF
         ENDDO
       ENDDO
       IF ((tabmin<90000.0_rsh).OR.(tabmax>109000.0_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine SFLXREAD.F90'
         WRITE(ierrorlog,*) 'In the meteorological file: ',TRIM(file_meteo)
         WRITE(ierrorlog,*) 'the atmospheric pressure is not expressed in Pa (=1e-5 bar).'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_sflx_pa)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       ENDIF
     ENDIF

     IF(l_sflx_rrh) THEN
       tabmin=MIN(MINVAL(rh1),MINVAL(rh2))
       tabmax=MAX(MAXVAL(rh1),MAXVAL(rh2))
       IF ((tabmin<0.0_rsh).OR.(tabmax>1.1_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine SFLXREAD.F90'
         WRITE(ierrorlog,*) 'In the meteorological file: ',TRIM(file_meteo)
         WRITE(ierrorlog,*) 'the relative humidity is not expressed in fraction (0.0 to 1.1).'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_sflx_rh)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       ENDIF
     ENDIF

     IF(l_sflx_rsat) THEN
       tabmin=MIN(MINVAL(sat1),MINVAL(sat2))
       tabmax=MAX(MAXVAL(sat1),MAXVAL(sat2))
       IF ((tabmin<-50.0_rsh).OR.(tabmax>50.0_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine SFLXREAD.F90'
         WRITE(ierrorlog,*) 'In the meteorological file: ',TRIM(file_meteo)
         WRITE(ierrorlog,*) 'the atmospheric temperature is not expressed in degrees Celcius.'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_sflx_sat)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       ENDIF
     ENDIF

     IF(l_sflx_rsolar) THEN
       tabmin=MIN(MINVAL(rad1),MINVAL(rad2))
       tabmax=MAX(MAXVAL(rad1),MAXVAL(rad2))
       IF ((tabmin<0.0_rsh).OR.(tabmax>1400.0_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine SFLXREAD.F90'
         WRITE(ierrorlog,*) 'In the meteorological file: ',TRIM(file_meteo)
         WRITE(ierrorlog,*) 'the solar radiation flux is not expressed in W/m2.'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_sflx_solar)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       ENDIF
     ENDIF

     IF(l_sflx_rir) THEN
       tabmin=MIN(MINVAL(irlos1),MINVAL(irlos2))
       tabmax=MAX(MAXVAL(irlos1),MAXVAL(irlos2))
       IF ((tabmin<-200.0_rsh).OR.(tabmax>200.0_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine SFLXREAD.F90'
         WRITE(ierrorlog,*) 'In the meteorological file: ',TRIM(file_meteo)
         WRITE(ierrorlog,*) 'the infra-red radiation flux (atm - sea) is not expressed in W/m2.'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_sflx_ir)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       ENDIF
     ENDIF

     IF(l_sflx_rrain) THEN
       tabmin=MIN(MINVAL(rain1),MINVAL(rain2))*(tnc2-tnc1)
       tabmax=MAX(MAXVAL(rain1),MAXVAL(rain2))*(tnc2-tnc1)
       IF ((tabmin<0.0_rsh).OR.(tabmax>1000.0_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine SFLXREAD.F90'
         WRITE(ierrorlog,*) 'In the meteorological file: ',TRIM(file_meteo)
         WRITE(ierrorlog,*) 'the precipitation amount is not expressed in kg.m-2 (=1 mm).'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_sflx_rain)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       ENDIF
     ENDIF

     IF(l_sflx_rfcc) THEN
       tabmin=MIN(MINVAL(fcc1),MINVAL(fcc2))
       tabmax=MAX(MAXVAL(fcc1),MAXVAL(fcc2))
#ifdef key_sflx_interp_scrip
       IF ((tabmin<-0.2_rsh).OR.(tabmax>1.2_rsh)) THEN
#else
       IF ((tabmin<0.0_rsh).OR.(tabmax>1.1_rsh)) THEN
#endif
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine SFLXREAD.F90'
         WRITE(ierrorlog,*) 'In the meteorological file: ',TRIM(file_meteo)
         WRITE(ierrorlog,*) 'the cloud cover fraction is not expressed in fraction (0.0 to 1.1).'
         WRITE(ierrorlog,*) 'variable name: ',TRIM(name_sflx_fcc)
         WRITE(ierrorlog,*) 'min. value = ',tabmin,', max. value = ',tabmax
         ncerror=ncerror+1
       ENDIF
     ENDIF

     IF(l_sflx_rwind) THEN
       tabmin=SQRT(MIN(MINVAL(vv1),MINVAL(vv2)))
       tabmax=SQRT(MAX(MAXVAL(vv1),MAXVAL(vv2)))
       IF ((tabmin<0.0_rsh).OR.(tabmax>50.0_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: routine SFLXREAD.F90'
         WRITE(ierrorlog,*) 'In the meteorological file: ',TRIM(file_meteo)
         WRITE(ierrorlog,*) 'the wind velocities are not expressed in m.s-1.'
         WRITE(ierrorlog,*) 'variable names: ',TRIM(name_sflx_windx),TRIM(name_sflx_windy)
         WRITE(ierrorlog,*) 'min. module = ',tabmin,', max. module = ',tabmax
         ncerror=ncerror+1
       ENDIF
     ENDIF

     IF (ncerror>0) THEN
       WRITE(ierrorlog,*) ' '
       WRITE(ierrorlog,*) 'Simulation stopped.'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     ENDIF
!$OMP end single

   END IF !endif for nbouc=1?

   ! deroulement normal

!$OMP single
   DO WHILE((tmet-(tnc1-torigin))*(tmet-(tnc2-torigin)) > 0.0_rlg .AND. ilec < idimt)
     ilec=ilec+1 
     IF (ilec == idimt) THEN
       WRITE(ierrorlog,*) ' '
       WRITE(ierrorlog,*) 'ERROR: routine SFLXREAD.F90'
       WRITE(ierrorlog,*) 'All records are already read.'
       WRITE(ierrorlog,*) 'The meteorological file is not long enough: ',TRIM(file_meteo)
       WRITE(ierrorlog,*) 'Simulation stopped.'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     END IF
     CALL ionc4_read_time(file_meteo,ilec,tnc1)
     CALL ionc4_read_time(file_meteo,ilec+1,tnc2)
   ENDDO
!$OMP end single

   IF(ilecmem/=ilec) THEN
     IF(nbouc/=1) THEN
!$OMP single
#ifdef key_grid_rotated 
       ALLOCATE(uvent1new(imetmin:imetmax,jmetmin:jmetmax))
       ALLOCATE(vvent1new(imetmin:imetmax,jmetmin:jmetmax))
#endif
#ifdef key_grid_curv
       ALLOCATE(tosx1tmp(limin:limax,ljmin:ljmax))
       ALLOCATE(tosy1tmp(limin:limax,ljmin:ljmax))
       ALLOCATE(tosx2tmp(limin:limax,ljmin:ljmax))
       ALLOCATE(tosy2tmp(limin:limax,ljmin:ljmax))
#endif
       ALLOCATE(uvent1(imetmin:imetmax,jmetmin:jmetmax))
       ALLOCATE(vvent1(imetmin:imetmax,jmetmin:jmetmax))
       ALLOCATE(  tosi(imetmin:imetmax,jmetmin:jmetmax))
       ALLOCATE( tosi2(imetmin:imetmax,jmetmin:jmetmax))
       ALLOCATE(  tab1(imetmin:imetmax,jmetmin:jmetmax))
!$OMP end single
     END IF

     ilecmem=ilec

     ! -1- on transfert les tableaux 2 en 1

!$OMP do schedule(static)
     DO j=ljmin,ljmax
       IF(l_sflx_rpa)    pa1(:,j)   =pa2(:,j)
       IF(l_sflx_rsat)   sat1(:,j)  =sat2(:,j)
       IF(l_sflx_rrh)    rh1(:,j)   =rh2(:,j)
       IF(l_sflx_rsolar) rad1(:,j)  =rad2(:,j)
       IF(l_sflx_rir)    irlos1(:,j)=irlos2(:,j)
       IF(l_sflx_rrain)  rain1(:,j) =rain2(:,j)
       IF(l_sflx_rfcc)   fcc1(:,j)  =fcc2(:,j)
       IF(l_sflx_rwind) THEN
         vv1(:,j)   =vv2(:,j)
         tosx1(:,j) =tosx2(:,j)
         tosy1(:,j) =tosy2(:,j)
#ifdef key_sflx_analyse
         windx1(:,j)=windx2(:,j)
         windy1(:,j)=windy2(:,j)
#endif
       ENDIF
     ENDDO
!$OMP end do
     IF(l_sflx_rpa) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_pa),tab1,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
#ifdef key_sflx_interp_scrip
!$OMP single
       IF (imetmin==-1) THEN
         tab1(imetmin,:)=tab1(imetmin+2,:) ; tab1(imetmin+1,:)=tab1(imetmin+2,:)
       END IF
       IF (jmetmin==-1) THEN
         tab1(:,jmetmin)=tab1(:,jmetmin+2) ; tab1(:,jmetmin+1)=tab1(:,jmetmin+2)
       END IF
       IF (imetmax==ixmet+2) THEN
         tab1(imetmax,:)=tab1(imetmax-2,:) ; tab1(imetmax-1,:)=tab1(imetmax-2,:)
       END IF
       IF (jmetmax==jxmet+2) THEN
         tab1(:,jmetmax)=tab1(:,jmetmax-2) ; tab1(:,jmetmax-1)=tab1(:,jmetmax-2)
       END IF
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       pa2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,pa2)
#endif
     END IF

     IF(l_sflx_rrh) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_rh),tab1,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       rh2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,rh2)
#endif
     END IF

     IF(l_sflx_rsat) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_sat),tab1,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       sat2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,sat2)
#endif
     END IF

     IF(l_sflx_rsolar) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_solar),tab1,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       rad2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,rad2)
#endif
     END IF

     IF(l_sflx_rir) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_ir),tab1,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       irlos2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,irlos2)
#endif
     END IF

     IF(l_sflx_rrain) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_rain),tab1,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       rain2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       IF (.NOT. l_sflx_rflux) WHERE( rain2(:,:) < 0.0_rsh ) rain2(:,:) = 0.0_rsh
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,rain2)
#endif
!$OMP single
       rain2(:,:)=rain2(:,:)/(tnc2-tnc1)
!$OMP end single
     END IF

     IF(l_sflx_rfcc) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_fcc),tab1,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       fcc2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,fcc2)
#endif
     END IF

     IF(l_sflx_rwind) THEN
!$OMP single
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_windy),vvent1,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
       CALL ionc4_read_subxyt(file_meteo,TRIM(name_sflx_windx),uvent1,imetmin,imetmax,jmetmin,jmetmax,ilec+1,1,1)
!$OMP end single
!$OMP do schedule(static)
       DO j=jmetmin,jmetmax

#ifdef key_grid_rotated  
         uvent1new(:,j)=uvent1(:,j)*COS(grid_angle*radeg)-vvent1(:,j)*SIN(grid_angle*radeg)
         vvent1new(:,j)=vvent1(:,j)*COS(grid_angle*radeg)+uvent1(:,j)*SIN(grid_angle*radeg) 
         uvent1(:,j)=uvent1new(:,j)
         vvent1(:,j)=vvent1new(:,j)
#endif

         tab1(:,j)=SQRT(uvent1(:,j)*uvent1(:,j)+vvent1(:,j)*vvent1(:,j))
       ENDDO
!$OMP enddo
       IF (.NOT. l_sflx_rflux) THEN
!$OMP do schedule(static)
         DO j= jmetmin,jmetmax
           tosi(:,j)=tab1(:,j)*uvent1(:,j)
           tosi2(:,j)=tab1(:,j)*vvent1(:,j)
         ENDDO
!$OMP end do
       ELSE
!$OMP do schedule(static)
         DO j= jmetmin,jmetmax
           tosi(:,j)=uvent1(:,j)
           tosi2(:,j)=vvent1(:,j)
         ENDDO
!$OMP end do
       END IF
#ifdef key_sflx_interp_scrip
!$OMP single
       tabvar2d(:,:)=tab1(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       vv2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       tabvar2d(:,:)=tosi(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       tosx2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
       tabvar2d(:,:)=tosi2(:,:)
       CALL make_bicubic_remap(tabvar2d,masksrc,tabinterp2d(limin:limax,ljmin:ljmax),nxfin,nyfin,matrix,src_add,dst_add)
       tosy2(:,:)=tabinterp2d(limin:limax,ljmin:ljmax)
!$OMP end single
#else
       CALL tool_interp2d_sflx(i1xe,j1xe,xdpxe,ydpxe,tab1,vv2)
       CALL tool_interp2d_sflx(i1u,j1u,xdpu,ydpu,tosi,tosx2)
       CALL tool_interp2d_sflx(i1v,j1v,xdpv,ydpv,tosi2,tosy2)
#ifdef key_sflx_analyse
       CALL tool_interp2d_sflx(i1u,j1u,xdpu,ydpu,uvent1,windx2)
       CALL tool_interp2d_sflx(i1v,j1v,xdpv,ydpv,vvent1,windy2)
#endif
#endif

#ifdef key_grid_curv
!$OMP do schedule(static)
     DO j=ljmin,ljmax
         tosx2tmp(limin:limax,j)=tosx2(limin:limax,j)*COS(grid_angle2D(limin:limax,j)*radeg)+tosy2(limin:limax,j)*SIN(grid_angle2D(limin:limax,j)*radeg)
         tosy2tmp(limin:limax,j)=tosy2(limin:limax,j)*COS(grid_angle2D(limin:limax,j)*radeg)-tosx2(limin:limax,j)*SIN(grid_angle2D(limin:limax,j)*radeg)
         tosx2(limin:limax,j)=tosx2tmp(limin:limax,j)
         tosy2(limin:limax,j)=tosy2tmp(limin:limax,j)
     ENDDO
!$OMP end do
#endif

     END IF

!$OMP single
#ifdef key_grid_rotated 
     DEALLOCATE(uvent1new,vvent1new)
#endif
#ifdef key_grid_curv
     DEALLOCATE(tosx1tmp,tosy1tmp,tosx2tmp,tosy2tmp)
#endif
     DEALLOCATE(uvent1,vvent1,tosi,tosi2,tab1)
!$OMP end single

   END IF

   dt1=(tnc2-torigin-tmet)/(tnc2-tnc1)
   dt2=1.0_rlg-dt1

!$OMP do schedule(static)
   DO j=ljmin,ljmax
     IF(l_sflx_rpa)   pa(limin:limax,j)       = pa1(limin:limax,j)*dt1   + pa2(limin:limax,j)*dt2
     IF(l_sflx_rsat)  sat(limin:limax,j)      = sat1(limin:limax,j)*dt1  + sat2(limin:limax,j)*dt2
     IF(l_sflx_rrh)   rh(limin:limax,j)       = rh1(limin:limax,j)*dt1   + rh2(limin:limax,j)*dt2
     IF(l_sflx_rsolar)rad(limin:limax,j)      = rad1(limin:limax,j)*dt1  + rad2(limin:limax,j)*dt2
     IF(l_sflx_rir)   irlos(limin:limax,j)    = irlos1(limin:limax,j)*dt1+ irlos2(limin:limax,j)*dt2
     IF(l_sflx_rrain) rain(limin:limax,j)     = rain1(limin:limax,j)*dt1 + rain2(limin:limax,j)*dt2
     IF(l_sflx_rfcc)  fcc(limin:limax,j)      = fcc1(limin:limax,j)*dt1  + fcc2(limin:limax,j)*dt2
     IF(l_sflx_rwind) THEN
#ifdef key_sflx_analyse
       windx(limin:limax,j) = windx1(limin:limax,j)*dt1   + windx2(limin:limax,j)*dt2
       windy(limin:limax,j) = windy1(limin:limax,j)*dt1   + windy2(limin:limax,j)*dt2
#endif
       ventvent(limin:limax,j) = vv1(limin:limax,j)*dt1   + vv2(limin:limax,j)*dt2
       ! the wind stress drag coefficient estimated at the same time that turbulent fluxes
       ! so, the wind stress estimate is completed in routines of turbulent fluxes
#if ! defined key_palm_ww3
       tosx(limin:limax,j)   = (tosx1(limin:limax,j)*dt1 + tosx2(limin:limax,j)*dt2)
       tosy(limin:limax,j)   = (tosy1(limin:limax,j)*dt1 + tosy2(limin:limax,j)*dt2)
#endif
     ENDIF
   ENDDO
!$OMP end do
!$OMP end parallel

   IF_MPI (l_sflx_rpa) THEN
     CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,pa(liminm2:limaxp2,ljminm2:ljmaxp2))
     CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,pa(liminm2:limaxp2,ljminm2:ljmaxp2))
   ENDIF_MPI

   PRINT_DBG*, 'exit SFLX_READ routine'

  END SUBROUTINE sflx_read

 END MODULE
