offline.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
! -------------------------------------------------
PROGRAM offline
!
! -------------------------------------------------
! Driver structure
! ----------------
! 1. Initializations
! 2. Temporal loops
!   2.a Read forcing
!   2.b Interpolate forcing in time
!   2.c Run surface
!   2.d Write prognostics and diagnostics variables
!
! modifications 
! 09/2012 G. Pigeon: coherence between radiation and zenith angle because of
!                    trouble with radiation received by wall in TEB
! 03/2014 E. Martin change indices names in OMP module according to GMAP changes
! 05/2014 B. Decharme delete trip
!                     I5 format to print DAY
! 04/2013 P. Lemoigne Add XDELTA_OROG to fix the maximum difference allowed between
!                     forcing and surface file orographies if LSET_FORC_ZS=.F
! 12/2013 S.Senesi    Add call to Gelato diag files init and close
!
! 02/2016: replace DOUBLE PRECISION by REAL to handle problem for promotion of real with GMKPACK or IBM SP
!
! -------------------------------------------------
!
USE modd_off_surfex_n
!
USE modd_forc_atm,  ONLY: csv         ,&! name of all scalar variables
                            xdir_alb    ,&! direct albedo for each band
                            xsca_alb    ,&! diffuse albedo for each band
                            xemis       ,&! emissivity
                            xtsrad      ,&! radiative temperature
                            xtsun       ,&! solar time                    (s from midnight)
                            xzs         ,&! orography                             (m)
                            xzref       ,&! height of T,q forcing                 (m)
                            xuref       ,&! height of wind forcing                (m)
                            xta         ,&! air temperature forcing               (K)
                            xqa         ,&! air humidity forcing                  (kg/m3)
                            xsv         ,&! scalar variables
                            xu          ,&! zonal wind                            (m/s)
                            xv          ,&! meridian wind                         (m/s)
                            xdir_sw     ,&! direct  solar radiation (on horizontal surf.)
                            xsca_sw     ,&! diffuse solar radiation (on horizontal surf.)
                            xsw_bands   ,&! mean wavelength of each shortwave band (m)
                            xzenith     ,&! zenithal angle       (radian from the vertical)
                            xzenith2    ,&! zenithal angle       (radian from the vertical)
                            xazim       ,&! azimuthal angle      (radian from North, clockwise)
                            xlw         ,&! longwave radiation (on horizontal surf.)
                            xps         ,&! pressure at atmospheric model surface (Pa)
                            xpa         ,&! pressure at forcing level             (Pa)
                            xrhoa       ,&! density at forcing level              (kg/m3)
                            xco2        ,&! CO2 concentration in the air          (kg/m3)
                            xsnow       ,&! snow precipitation                    (kg/m2/s)
                            xrain       ,&! liquid precipitation                  (kg/m2/s)
                            xsfth       ,&! flux of heat                          (W/m2)
                            xsftq       ,&! flux of water vapor                   (kg/m2/s)
                            xsfu        ,&! zonal momentum flux                   (m/s)
                            xsfv        ,&! meridian momentum flux                (m/s)
                            xsfco2      ,&! flux of CO2                           (kg/m2/s)
                            xsfts       ,&! flux of scalar var.                   (kg/m2/s)
                            xpew_a_coef ,&! implicit coefficients
                            xpew_b_coef ,&! needed if HCOUPLING='I'
                            xpet_a_coef ,&
                            xpeq_a_coef ,&
                            xpet_b_coef ,&
                            xpeq_b_coef ,&
                            xtsurf      ,&! effective temperature                  (K)
                            xz0         ,&! surface roughness length for momentum  (m)
                            xz0h        ,&! surface roughness length for heat      (m)
                            xqsurf        ! specific humidity at surface           (kg/kg)
!
USE modd_surf_conf,  ONLY : cprogname, csoftware
USE modd_csts,       ONLY : xpi, xday, xrv, xrd, xg
USE modd_io_surf_asc,ONLY : cfilein,cfilein_save,cfileout,cfilepgd
USE modd_surf_par
USE modd_io_surf_fa, ONLY : cfilein_fa, cfilein_fa_save,       &
                            cfileout_fa, nunit_fa, cdnomc,     &
                            iverbfa, lfanocompact, cfilepgd_fa  
USE modd_io_surf_lfi,ONLY : cfilein_lfi, cfilein_lfi_save, cluout_lfi, cfileout_lfi, &
                            lmnh_compatible, cfilepgd_lfi  
USE modd_io_surf_nc, ONLY : cfilein_nc, cfilein_nc_save, cfileout_nc, cluout_nc, &
                            cfilepgd_nc, ldef
USE modd_io_surf_ol, ONLY : xstart, xcount, xstride,           &
                              ldefined_nature, ldefined_sea,    &
                              ldefined_water,  ldefined_town,   &
                              ldefined_surf_atm, lpartw,        &
                              xstartw, xcountw, ltime_written,  &
                              nstep_output  
USE modd_write_bin,  ONLY : nwrite
!
USE modd_surfex_mpi, ONLY : ncomm, nproc, nrank, npio, wlog_mpi, prep_log_mpi,   &
                            nindex, nsize_task, xtime_npio_read, xtime_npio_write, &
                            xtime_comm_read, xtime_comm_write, xtime_sea, &
                            xtime_nature, xtime_water, xtime_town, &
                            xtime_init_sea, xtime_init_water, &
                            xtime_init_nature, xtime_init_town, &
                            xtime_write, xtime_calc, xtime_omp_barr, &
                            xtime_calc_write, idx_w, end_log_mpi
!
USE modd_surfex_omp, ONLY :  nindx1sfx, nindx2sfx, nblock, nblocktot, &
                             init_dim, reset_dim, nwork, nwork2, xwork, xwork2, xwork3, &
                             nwork_full, nwork2_full, xwork_full, xwork2_full
!
USE modd_coupling_topd, ONLY : nnb_topd, nnb_stp_restart, lbudget_topd, ltopd_step, &
                               lcoupl_topd, ntopd_step, nyear, nmonth, nday, nh, nm
USE modd_topodyn, ONLY : xtopd_step, nnb_topd_step, xqtot, xqb_run, xqb_dr
!
USE modd_slope_effect, ONLY: xzs_thread,xzs_xy_thread,xslopang_thread,&
                             xslopazi_thread,xsurf_triangle_thread
!
USE modd_sfx_oasis, ONLY : loasis, xruntime
!
USE mode_pos_surf
!
USE mode_crodebug
!
USE modn_io_offline
!
USE modi_get_luout
USE modi_open_namelist
USE modi_test_nam_var_surf
USE modi_close_namelist
USE modi_read_all_namelists
USE modi_open_close_bin_asc_forc
USE modi_open_filein_ol
USE modi_ol_read_atm_conf
USE modi_abor1_sfx
USE modi_ol_alloc_atm
USE modi_compare_orography
USE modi_sunpos
USE modi_init_index_mpi
USE modi_ol_read_atm
USE modi_get_sizes_parallel
USE modi_io_buff_clean
USE modi_init_surf_atm_n
USE modi_init_surf_landuse_n
USE modi_ol_time_interp_atm
USE modi_coupling_surf_atm_n
USE modi_add_forecast_to_date_surf
USE modi_write_surf_atm_n
USE modi_write_header_mnh
USE modi_flag_update
USE modi_flag_diag_update
USE modi_diag_surf_atm_n
USE modi_write_diag_surf_atm_n
USE modi_get_surf_var_n
USE modi_gather_and_write_mpi
USE modi_close_filein_ol
USE modi_close_fileout_ol
USE modi_init_output_ol_n
USE modi_init_output_nc_n
!
USE modi_write_header_fa
USE modi_abor1_sfx
!
USE modi_write_discharge_file
USE modi_write_budget_coupl_rout
USE modi_prep_restart_coupl_topd
!
USE modi_init_slope_param
USE modi_slope_radiative_effect
!
USE modi_sfx_oasis_init
USE modi_sfx_oasis_read_nam
USE modi_sfx_oasis_def_ol
USE modi_sfx_oasis_recv_ol
USE modi_sfx_oasis_send_ol
USE modi_sfx_oasis_end
!RJ: missing modi
USE modi_local_slope_param
!
USE mode_glt_dia_lu
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
#ifdef AIX64
!$ USE OMP_LIB
#endif
!
IMPLICIT NONE
!
#ifdef SFX_MPI
include 'mpif.h'
#endif
!
#ifndef AIX64
!$ INCLUDE 'omp_lib.h'
#endif
!
!*      0.    declarations of local variables
!
 CHARACTER(LEN=3), PARAMETER       :: yinit     = 'ALL'
!
 CHARACTER(LEN=28)                 :: yluout    = 'LISTING_OFFLINE             '
!
INTEGER                           :: iyear               ! current year (UTC)
INTEGER                           :: imonth              ! current month (UTC)
INTEGER                           :: iday                ! current day (UTC)
REAL                              :: ztime               ! current time since start of the run (s)
REAL                              :: ztimec              ! current duration since start of the run (s)
!
INTEGER                           :: iyear_out           ! output year name
INTEGER                           :: imonth_out          ! output month name
INTEGER                           :: iday_out            ! output day name
REAL                              :: ztime_out           ! output time since start of the run (s)
!
INTEGER, DIMENSION(11)  :: idatef
!
 CHARACTER(LEN=28), PARAMETER      :: yatmfile     = '                            '
 CHARACTER(LEN=6),  PARAMETER      :: yatmfiletype = '      '
 CHARACTER(LEN=2),  PARAMETER      :: ytest        = 'OK'          ! must be equal to 'OK'
!
REAL, DIMENSION(:), POINTER       :: zlat                ! latitude                         (rad)
REAL, DIMENSION(:), POINTER       :: zlon                ! longitude                        (rad)
REAL, DIMENSION(:), POINTER       :: zzs_forc            ! orography                        (m)  
REAL, DIMENSION(:), POINTER       :: zzref               ! Forcing level for T
REAL, DIMENSION(:), POINTER       :: zuref               ! Forcing level for U
!
REAL                              :: ztstep              ! atmospheric time-step            (s)
!
INTEGER                           :: ini                 ! grid dimension
INTEGER                           :: jloop               ! loop counter
INTEGER                           :: ibands              ! Number of radiative bands 
INTEGER                           :: inb_step_atm        ! Number of atmospheric time-steps
INTEGER                           :: inb_atm             ! Number of Isba time-steps 
                                                         ! within a forcing time-step
INTEGER                           :: id_forc             ! indice of forcing in the file
INTEGER                           :: inb_lines           ! nb of lines to read in the forcing file
INTEGER                           :: idmax               ! nb of lines to read in the forcing file at last 
INTEGER                           :: jforc_step          ! atmospheric loop index
INTEGER                           :: jsurf_step          ! isba loop index
INTEGER                           :: icount              ! day counter 
REAL                              :: zduration           ! duration of run                     (s)
!
REAL, DIMENSION(:,:), ALLOCATABLE :: zta                 ! air temperature forcing               (K)
REAL, DIMENSION(:,:), ALLOCATABLE :: zqa                 ! air humidity forcing                  (kg/m3)
REAL, DIMENSION(:,:), ALLOCATABLE :: zwind               ! wind speed                            (m/s)
REAL, DIMENSION(:,:), ALLOCATABLE :: zsca_sw             ! diffuse solar radiation (on horizontal surf.)
REAL, DIMENSION(:,:), ALLOCATABLE :: zdir_sw             ! direct  solar radiation (on horizontal surf.)
REAL, DIMENSION(:,:), ALLOCATABLE :: zlw                 ! longwave radiation (on horizontal surf.)
REAL, DIMENSION(:,:), ALLOCATABLE :: zsnow               ! snow precipitation                    (kg/m2/s)
REAL, DIMENSION(:,:), ALLOCATABLE :: zrain               ! liquid precipitation                  (kg/m2/s)
REAL, DIMENSION(:,:), ALLOCATABLE :: zps                 ! pressure at forcing level             (Pa)
REAL, DIMENSION(:,:), ALLOCATABLE :: zco2                ! CO2 concentration in the air          (kg/m3)
REAL, DIMENSION(:,:), ALLOCATABLE :: zdir                ! wind direction
INTEGER                           :: iluout              ! ascii output unit number
INTEGER                           :: ilunam              ! namelist unit number
INTEGER                           :: iret                ! error return code
INTEGER                           :: inb 
INTEGER                           :: inw, jnw
 CHARACTER(LEN=14)                :: ytag      
LOGICAL                           :: gfound              ! return logical when reading namelist
LOGICAL                           :: gshadows    
REAL, DIMENSION(:),   ALLOCATABLE :: zsw                 ! total solar radiation (on horizontal surf.)
REAL, DIMENSION(:),   ALLOCATABLE :: zcoef               ! coefficient for solar radiation interpolation near sunset/sunrise
!
! Flag diag :
!
INTEGER                           :: i2m, ibeq, idsteq
LOGICAL                           :: gfrac, gdiag_grid, gsurf_budget, grad_budget, gcoef,    &
                                     gsurf_vars, gdiag_ocean, gdiag_seaice, gwater_profile,  &
                                     ginterpol_ts, gsurf_evap_budget, gflood, gpgd_isba,     &
                                     gch_no_flux_isba, gsurf_misc_budget_isba, gpgd_teb,     &
                                     gsurf_misc_budget_teb
!
! Inquiry mode arrays:
!
REAL, DIMENSION(:), ALLOCATABLE   :: zsea, zwater, znature, ztown
REAL, DIMENSION(:), ALLOCATABLE   :: zsea_full, zwater_full, znature_full, ztown_full
REAL, DIMENSION(:), ALLOCATABLE   :: zt2m, zq2m
REAL, DIMENSION(:), ALLOCATABLE   :: zz0, zz0h, zqs
REAL, DIMENSION(:), ALLOCATABLE   :: zqs_sea, zqs_water, zqs_nature, zqs_town
REAL, DIMENSION(:), ALLOCATABLE   :: zpsng, zpsnv
REAL, DIMENSION(:), ALLOCATABLE   :: zz0eff
REAL, DIMENSION(:), ALLOCATABLE   :: zzs
REAL, DIMENSION(:), ALLOCATABLE   :: zz0_full, zz0eff_full, zzs_full
INTEGER :: iseries, isize
!
! MPI variables
!
 CHARACTER(LEN=100) :: yname
 CHARACTER(LEN=10)  :: yrank
INTEGER :: ilevel, infompi, j, inkproma, jblock
INTEGER, DIMENSION(:), ALLOCATABLE :: isize_omp
REAL :: xtime0, xtime1, xtime
!
! SFX - OASIS coupling variables
!
INTEGER  :: ilocal_comm  ! Local communicator
!
REAL(KIND=JPRB) :: zhook_handle
!
! --------------------------------------------------------------------------------------
!
!*     0.1.   MPI, OASIS, and dr_hook initializations
!
infompi=1
!
#ifdef CPLOASIS
!Must be call before DRHOOK !
 CALL sfx_oasis_init(cnamelist,ilocal_comm)
#else
loasis   = .false.
xruntime = 0.0
#endif
!
#ifdef SFX_MPI
IF(.NOT.loasis)THEN
 CALL mpi_init_thread(mpi_thread_multiple,ilevel,infompi)
 IF (infompi /= mpi_success) THEN 
    CALL abor1_sfx('OFFLINE: ERROR WHEN INITIALIZING MPI')
 ENDIF
ENDIF
#endif
!
IF (lhook) CALL dr_hook('OFFLINE',0,zhook_handle)
!
 csoftware='OFFLINE'
!
#ifdef SFX_MPI
IF(loasis)THEN
  ncomm=ilocal_comm
ELSE
  ncomm=mpi_comm_world
ENDIF
 CALL mpi_comm_size(ncomm,nproc,infompi)
 CALL mpi_comm_rank(ncomm,nrank,infompi)
#endif
!
!RJ: init modd_surefx_omp
!$OMP PARALLEL
!$ NBLOCKTOT = OMP_GET_NUM_THREADS()
!$ NBLOCK = OMP_GET_THREAD_NUM()
!$OMP END PARALLEL
!
 CALL prep_log_mpi
!
 CALL wlog_mpi(' ')
!
 CALL wlog_mpi('NBLOCKTOT ',klog=nblocktot)
!
#ifdef SFX_MPI
xtime0 = mpi_wtime()
#endif
!
!
!*      0.3.   Open ascii file for writing
!
WRITE(yrank,fmt='(I10)') nrank
yname=trim(yluout)//adjustl(yrank)
!
 cluout_lfi =  adjustl(adjustr(yname)//'.txt')
 cluout_nc  =  adjustl(adjustr(yname)//'.txt')
!
 CALL get_luout('ASCII ',iluout)
OPEN(unit=iluout,file=adjustl(adjustr(yname)//'.txt'),form='FORMATTED',action='WRITE')
!
!
IF ( nrank==npio ) THEN
  !
!$OMP SINGLE
  !
!RJ: be verbose just for openmp
  IF(nblocktot==1) THEN
!$  WRITE(*,*) "CAUTION: DID YOU THINK TO SET OMP_NUM_THREADS=1?"
!$  WRITE(*,*) "PLEASE VERIFY OMP_NUM_THREADS IS INITIALIZED : TYPE ECHO $OMP_NUM_THREADS IN A TERMINAL" 
  !
!$  WRITE(ILUOUT,*) "CAUTION: DID YOU THINK TO SET OMP_NUM_THREADS=1?"
!$  WRITE(ILUOUT,*) "PLEASE VERIFY OMP_NUM_THREADS IS INITIALIZED : TYPE ECHO $OMP_NUM_THREADS IN A TERMINAL" 
  ENDIF
  !
!$OMP END SINGLE
  !
ENDIF
!
!*      0.4.   Reads namelists
!
 CALL open_namelist('ASCII ',ilunam,cnamelist)
!
 CALL posnam(ilunam,'NAM_IO_OFFLINE',gfound,iluout)
IF (gfound) READ (unit=ilunam,nml=nam_io_offline)
 CALL close_namelist('ASCII ',ilunam)
!
IF (nproc==1) THEN 
  xio_frac=1.
ELSE
  xio_frac = max(min(xio_frac,1.),0.)
ENDIF
!
 CALL test_nam_var_surf(iluout,'CSURF_FILETYPE',csurf_filetype,'ASCII ','LFI   ','FA    ','NC    ')
 CALL test_nam_var_surf(iluout,'CTIMESERIES_FILETYPE',ctimeseries_filetype,'NETCDF','TEXTE ','BINARY',&
                                                                            'ASCII ','LFI   ','FA    ',&
                                                                            'NONE  ','OFFLIN','NC    ')  
 CALL test_nam_var_surf(iluout,'CFORCING_FILETYPE',cforcing_filetype,'NETCDF','ASCII ','BINARY')
!
IF (nscal>59) CALL abor1_sfx("OFFLINE: NSCAL MUST BE LOWER THAN OR EQUAL TO 59")
!
!
IF (ctimeseries_filetype=='NETCDF') ctimeseries_filetype='OFFLIN'
!
IF (nrank==npio) THEN
  !
  cfilepgd = adjustl(adjustr(cpgdfile)//'.txt')
  cfilein  = adjustl(adjustr(cprepfile)//'.txt')
  cfilein_save = cfilein
  !
  cfilepgd_lfi = cpgdfile
  cfilein_lfi  = cprepfile
  cfilein_lfi_save = cfilein_lfi
  !
  cfilepgd_fa = adjustl(adjustr(cpgdfile)//'.fa')
  cfilein_fa  = adjustl(adjustr(cprepfile)//'.fa')
  cfilein_fa_save  = cfilein_fa
  !
  cfilepgd_nc = adjustl(adjustr(cpgdfile)//'.nc')
  cfilein_nc  = adjustl(adjustr(cprepfile)//'.nc')
  cfilein_nc_save  = cfilein_nc
  !
ENDIF
!
!     Allocations of Surfex Types
 CALL surfex_alloc_list(1)
 ysurf_cur => ysurf_list(1)
!
!     Reading all namelist (also assimilation)
 CALL read_all_namelists(ysurf_cur, &
                         csurf_filetype,'ALL',.false.)
!
!
!*      0.5.   Reads SFX - OASIS coupling namelists
!
 CALL sfx_oasis_read_nam(csurf_filetype,xtstep_surf)
!
!*      0.6   Assume FA filetype consistency 
!
 cprogname = csurf_filetype
!
! --------------------------------------------------------------------------------------
!
!*      1.    Initializations
!
!       netcdf file handling
!
IF (nrank==npio) THEN
  !
  xstart            = nundef
  xstride           = nundef
  xcount            = nundef
  xstartw           = 0
  xcountw           = 1
  lpartw            = .true.
  ldefined_surf_atm = .false.
  ldefined_nature   = .false.
  ldefined_town     = .false.
  ldefined_water    = .false.
  ldefined_sea      = .false.
  !
ENDIF
!
#ifdef SFX_MPI
xtime = (mpi_wtime() - xtime0)
#endif
 CALL wlog_mpi('READ NAMELISTS ',plog=xtime)
#ifdef SFX_MPI
xtime0 = mpi_wtime()
#endif
!
!       forcing file handling
!
IF (cforcing_filetype=='ASCII ' .OR. cforcing_filetype=='BINARY') CALL open_close_bin_asc_forc('CONF ',cforcing_filetype,'R')
IF (cforcing_filetype=='NETCDF') CALL open_filein_ol
!
!
!       splitting of the grid
!
gshadows = lshadows_slope .OR. lshadows_other
 CALL init_index_mpi(ysurf_cur, &
                    csurf_filetype, yalg_mpi, xio_frac, gshadows)
!
 CALL wlog_mpi(' ')
 CALL wlog_mpi('TIME_NPIO_READ init_index ',plog=xtime_npio_read)
 CALL wlog_mpi('TIME_COMM_READ init_index ',plog=xtime_comm_read)
xtime_npio_read = 0.
xtime_comm_read = 0.
!
#ifdef SFX_MPI
xtime = (mpi_wtime() - xtime0)
#endif
 CALL wlog_mpi(' ')
 CALL wlog_mpi('INIT_INDEX_MPI ',plog=xtime)
 CALL wlog_mpi(' ')
#ifdef SFX_MPI
xtime0 = mpi_wtime()
#endif
!
!       configuration of run
!
 CALL ol_read_atm_conf(ysurf_cur, &
                       csurf_filetype, cforcing_filetype,            &
                      zduration, ztstep, ini, iyear, imonth, iday,  &
                      ztime, zlat, zlon, zzs_forc, zzref, zuref     )
!
 CALL wlog_mpi(' ')
 CALL wlog_mpi('TIME_NPIO_READ forc conf ',plog=xtime_npio_read)
 CALL wlog_mpi('TIME_COMM_READ forc conf ',plog=xtime_comm_read)
xtime_npio_read = 0.
xtime_comm_read = 0.
!
#ifdef SFX_MPI
xtime = (mpi_wtime() - xtime0)
#endif
 CALL wlog_mpi('OL_READ_ATM_CONF ',plog=xtime)
 CALL wlog_mpi(' ')
#ifdef SFX_MPI
xtime0 = mpi_wtime()
#endif
!
!*     time steps coherence check 
!
IF ( (mod(xtstep_output,ztstep)*mod(ztstep,xtstep_output) /= 0) .OR. (mod(ztstep,xtstep_surf) /= 0) ) THEN
   WRITE(iluout,*)' FORCING  AND OUTPUT/SURFACE TIME STEP SHOULD BE MULTIPLE', &
     nint(ztstep),nint(xtstep_output),nint(xtstep_surf)    
   CALL abor1_sfx('OFFLINE: FORCING  AND OUTPUT/SURFACE TIME STEP SHOULD BE MULTIPLE')
ENDIF
!
IF ( ztime /= 0. .AND. mod(ztime,xtstep_surf) /= 0  ) THEN
   WRITE(iluout,*)' INITIAL AND SURFACE TIME STEP SHOULD BE MULTIPLE', &
   nint(ztime),nint(xtstep_surf)  
   CALL abor1_sfx('OFFLINE: INITIAL AND SURFACE TIME STEP SHOULD BE MULTIPLE')
ENDIF
!
IF(loasis.AND.zduration/=xruntime)THEN
   WRITE(iluout,*)'Total simulated time given by Forcing field and OASIS namcouple are different'
   WRITE(iluout,*)'From Forcing (s) : ',zduration, 'From OASIS   (s) : ',xruntime
   CALL abor1_sfx('OFFLINE: TOTAL SIMULATED TIME DIFFERENT BETWEEN FORCING AND OASIS')
ENDIF
!
inb_step_atm  = int(zduration / ztstep)
inb_atm       = int(ztstep / xtstep_surf)
nstep_output  = int(zduration / xtstep_output)
!
xtopd_step = 0
nnb_topd_step = 0
ntopd_step = 0
IF ( lcoupl_topd ) THEN
  !
  xtopd_step = float(nnb_topd)* xtstep_surf
  nnb_topd_step = int( zduration / xtopd_step )
  !
  IF ( nnb_stp_restart==0 .AND. .NOT.lrestart ) nnb_stp_restart = -1
  !
  ntopd_step = 1
  !
ENDIF
!
!       allocation of variables
!
ibands = 1
!
 CALL ol_alloc_atm(ini,ibands,nscal)
!
xzs   = zzs_forc
xzref = zzref
xuref = zuref
!
!       compare orography
!
 CALL compare_orography(ysurf_cur, &
                         csurf_filetype, lset_forc_zs, xdelta_orog)
!
!       miscellaneous initialization
!
icount = 0
ztimec = 0.
!
ALLOCATE(isize_omp(0:nblocktot-1))
 CALL get_sizes_parallel(ysurf_cur%DTCO, ysurf_cur%DGU, ysurf_cur%UG, ysurf_cur%U, &
                         nblocktot,ini,0,isize_omp, gshadows)
 CALL sunpos(isize_omp, iyear, imonth, iday, ztime, zlon, zlat, xtsun, xzenith, xazim)
DEALLOCATE(isize_omp)
!
!number of lines read in forcing files
inb_lines=1
IF (nb_read_forc.EQ.1) THEN
  inb_lines=inb_step_atm
ELSEIF (nb_read_forc.NE.0) THEN
  !to be sure the number of readings will be NB_READ_FORC as a maximum
  inb_lines=ceiling(1.*(inb_step_atm+1)/nb_read_forc)
ENDIF
!number of lines to be read effectively
idmax=inb_lines+1
!effective number of readings of the forcing files
nb_read_forc=ceiling(1.*(inb_step_atm+1)/inb_lines)
!
!     open Gelato specific diagnostic files (if requested by 
!     Gelato wizzard user)
!
#if ! defined in_arpege
 CALL opndia()
#endif
!
!       allocate local atmospheric variables
!
IF (.NOT.ALLOCATED(zta)) ALLOCATE(zta(ini,inb_lines+1)) 
IF (.NOT.ALLOCATED(zqa))ALLOCATE(zqa(ini,inb_lines+1))
IF (.NOT.ALLOCATED(zwind))ALLOCATE(zwind(ini,inb_lines+1))
IF (.NOT.ALLOCATED(zdir_sw))ALLOCATE(zdir_sw(ini,inb_lines+1))
IF (.NOT.ALLOCATED(zsca_sw))ALLOCATE(zsca_sw(ini,inb_lines+1))
IF (.NOT.ALLOCATED(zlw))ALLOCATE(zlw(ini,inb_lines+1))
IF (.NOT.ALLOCATED(zsnow))ALLOCATE(zsnow(ini,inb_lines+1))
IF (.NOT.ALLOCATED(zrain))ALLOCATE(zrain(ini,inb_lines+1))
IF (.NOT.ALLOCATED(zps))ALLOCATE(zps(ini,inb_lines+1))
IF (.NOT.ALLOCATED(zco2))ALLOCATE(zco2(ini,inb_lines+1))
IF (.NOT.ALLOCATED(zdir))ALLOCATE(zdir(ini,inb_lines+1))
IF (.NOT.ALLOCATED(zcoef))ALLOCATE(zcoef(ini))
!
IF (.NOT.ALLOCATED(zsw))ALLOCATE(zsw(ini))
!
!      computes initial air co2 concentration and  density
!
#ifdef SFX_MPI
xtime = (mpi_wtime() - xtime0)
#endif
 CALL wlog_mpi('COMPARE_OROGRAPHY SUNPOS ',plog=xtime)
#ifdef SFX_MPI
xtime0 = mpi_wtime()
#endif
!
!* opens forcing files (if ASCII or BINARY)
!
IF (cforcing_filetype=='ASCII ' .OR. cforcing_filetype=='BINARY') CALL open_close_bin_asc_forc('OPEN ',cforcing_filetype,'R')
!
 CALL ol_read_atm(&
                  csurf_filetype, cforcing_filetype, 1,             &
                   zta,zqa,zwind,zdir_sw,zsca_sw,zlw,zsnow,zrain,zps,&
                   zco2,zdir,llimit_qair                           ) 
!
 CALL wlog_mpi(' ')
 CALL wlog_mpi('TIME_NPIO_READ forc ',plog=xtime_npio_read)
 CALL wlog_mpi('TIME_COMM_READ forc ',plog=xtime_comm_read)
xtime_npio_read = 0.
xtime_comm_read = 0.
!
#ifdef SFX_MPI
xtime = (mpi_wtime() - xtime0)
#endif
 CALL wlog_mpi(' ')
 CALL wlog_mpi('OL_READ_ATM0 ',plog=xtime)
 CALL wlog_mpi(' ')
#ifdef SFX_MPI
xtime0 = mpi_wtime()
#endif
!
xco2(:)  = zco2(:,1)
xrhoa(:) = zps(:,1) / (xrd * zta(:,1) * ( 1.+((xrv/xrd)-1.)*zqa(:,1) ) + xg * xzref )
!                 
!       surface Initialisation     
!
#ifdef SFX_MPI
xtime = (mpi_wtime() - xtime0)
#endif
 CALL wlog_mpi('CO2 RHOA ',plog=xtime)
!
 CALL io_buff_clean
!
 CALL surfex_deallo_list
 CALL surfex_alloc_list(nblocktot)
!
 CALL goto_model(1)
ALLOCATE(isize_omp(0:nblocktot-1))
 CALL get_sizes_parallel(ysurf_cur%DTCO, ysurf_cur%DGU, ysurf_cur%UG, ysurf_cur%U, &
                         nblocktot,ini,0,isize_omp, gshadows)
DO j=0,nblocktot-1
 CALL wlog_mpi("SIZES_OMP ",klog=j,klog2=isize_omp(j))
ENDDO
!
!$OMP PARALLEL PRIVATE(INKPROMA,XTIME,XTIME0)
!
#ifdef SFX_MPI
xtime0 = mpi_wtime()
#endif
!
!$ NBLOCK = OMP_GET_THREAD_NUM()
!
IF (nblock==nblocktot) THEN
  CALL init_dim(isize_omp,0,inkproma,nindx1sfx,nindx2sfx)
ELSE
  CALL init_dim(isize_omp,nblock,inkproma,nindx1sfx,nindx2sfx)
ENDIF
!
IF (nblock==0) THEN
  CALL goto_model(nblocktot)
ELSE
  CALL goto_model(nblock)
ENDIF
!
 CALL init_surf_atm_n(ysurf_cur, &
                            csurf_filetype, yinit, lland_use,                      &
                     inkproma, nscal, ibands,                               &
                     csv,xco2(nindx1sfx:nindx2sfx),xrhoa(nindx1sfx:nindx2sfx),          &
                     xzenith(nindx1sfx:nindx2sfx),xazim(nindx1sfx:nindx2sfx),xsw_bands, &
                     xdir_alb(nindx1sfx:nindx2sfx,:), xsca_alb(nindx1sfx:nindx2sfx,:),  &
                     xemis(nindx1sfx:nindx2sfx), xtsrad(nindx1sfx:nindx2sfx),           &
                     xtsurf(nindx1sfx:nindx2sfx),                           &
                     iyear, imonth, iday, ztime,                            &
                     yatmfile, yatmfiletype, ytest                          )
!
! initialization routines to compute shadows
IF (gshadows) THEN
  IF (nblock==0) THEN
    CALL init_slope_param(ysurf_cur%UG, &
                          zzs_forc,ini,zlat)
  END IF 
  !$OMP BARRIER
  CALL local_slope_param(nindx1sfx,nindx2sfx)
END IF
!
 CALL reset_dim(ini,inkproma,nindx1sfx,nindx2sfx)
!
#ifdef SFX_MPI
xtime = (mpi_wtime() - xtime0)
#endif
 CALL wlog_mpi(' ')
 CALL wlog_mpi('INIT_SURF_ATM ',plog=xtime)
 CALL wlog_mpi(' ')
!
 CALL wlog_mpi('TIME_NPIO_READ init ',plog=xtime_npio_read)
 CALL wlog_mpi('TIME_COMM_READ init ',plog=xtime_comm_read)
 CALL wlog_mpi(' ')
!
!$OMP END PARALLEL
!
xtime_npio_read = 0.
xtime_comm_read = 0.
!
!   Land use or/and vegetation dynamic
!                  
 CALL init_surf_landuse_n(ysurf_cur, &
                          csurf_filetype,yinit,lland_use,             &
                       ini, nscal, ibands,                           &
                       csv,xco2(nindx1sfx:nindx2sfx),xrhoa(nindx1sfx:nindx2sfx), &
                       xzenith(nindx1sfx:nindx2sfx),xazim(nindx1sfx:nindx2sfx),  &
                       xsw_bands,xdir_alb(nindx1sfx:nindx2sfx,:),          &
                       xsca_alb(nindx1sfx:nindx2sfx,:),                    &
                       xemis(nindx1sfx:nindx2sfx),xtsrad(nindx1sfx:nindx2sfx),   &
                       xtsurf(nindx1sfx:nindx2sfx),                           &
                       iyear, imonth, iday, ztime,                   &
                       yatmfile, yatmfiletype, ytest                 )
!
#ifdef SFX_MPI
xtime0 = mpi_wtime()
#endif
!
 CALL init_crodebug(ysurf_cur%IM%I)
!
! * SURFEX - OASIS  grid, partitions and local field definitions
!
IF(loasis)THEN
  CALL sfx_oasis_def_ol(ysurf_cur%IM%I, &
                        ysurf_cur%U, &
                        csurf_filetype,yalg_mpi)
ENDIF
!
IF (ctimeseries_filetype=="OFFLIN") THEN
  CALL goto_model(1)
  CALL init_output_ol_n(ysurf_cur)
ENDIF
! --------------------------------------------------------------------------------------
!
inw = 1
IF (ctimeseries_filetype=="NC    ") inw = 2
!
nwrite = 0
!
#ifdef SFX_MPI
xtime = (mpi_wtime() - xtime0)
#endif
 CALL wlog_mpi('INIT FINISHED ',plog=xtime)
#ifdef SFX_MPI
xtime0 = mpi_wtime()
#endif
!*      2.    Temporal loops
!
xtime_calc(:) = 0.
xtime_write(:) = 0.
!
DO jforc_step=1,inb_step_atm
  !
#ifdef SFX_MPI
  xtime1 = mpi_wtime()
#endif
  ! read Forcing
  !
  !indice of forcing line in forcing arrays
  id_forc=jforc_step-int(jforc_step/inb_lines)*inb_lines
  IF (id_forc==0) id_forc=inb_lines
  !new forcings to read
  IF (id_forc==1 .AND. jforc_step.NE.1) THEN
    !if last part of forcing, the last point has to be adjusted on the end of
    !files
    IF (jforc_step/inb_lines==nb_read_forc-1) THEN 
      idmax=inb_step_atm-jforc_step+1+1
      !for ascii and binary forcing files
      zta(:,idmax)=zta(:,SIZE(zta,2))
      zqa(:,idmax)=zqa(:,SIZE(zta,2))
      zwind(:,idmax)=zwind(:,SIZE(zta,2))
      zdir_sw(:,idmax)=zdir_sw(:,SIZE(zta,2))
      zsca_sw(:,idmax)=zsca_sw(:,SIZE(zta,2))
      zlw(:,idmax)=zlw(:,SIZE(zta,2))
      zsnow(:,idmax)=zsnow(:,SIZE(zta,2))
      zrain(:,idmax)=zrain(:,SIZE(zta,2))
      zps(:,idmax)=zps(:,SIZE(zta,2))
      zco2(:,idmax)=zco2(:,SIZE(zta,2))
      zdir(:,idmax)=zdir(:,SIZE(zta,2))
    ENDIF
    CALL ol_read_atm(&
                  csurf_filetype, cforcing_filetype, jforc_step,    &
                     zta(:,1:idmax),zqa(:,1:idmax),zwind(:,1:idmax), &
                     zdir_sw(:,1:idmax),zsca_sw(:,1:idmax),zlw(:,1:idmax), &
                     zsnow(:,1:idmax),zrain(:,1:idmax),zps(:,1:idmax),&
                     zco2(:,1:idmax),zdir(:,1:idmax),llimit_qair         )
  ENDIF

#ifdef SFX_MPI
  xtime_calc(1) = xtime_calc(1) + (mpi_wtime() - xtime1)
  xtime1 = mpi_wtime()
#endif
  !
  DO jsurf_step=1,inb_atm
    !
    ! time interpolation of the forcing
    !
#ifdef SFX_MPI
    xtime1 = mpi_wtime()
#endif
    CALL sunpos(isize_omp, iyear, imonth, iday, ztime, &
                zlon, zlat, xtsun, xzenith, xazim)
    CALL sunpos(isize_omp, iyear, imonth, iday, ztime+xtstep_surf, &
                zlon, zlat, xtsun, xzenith2, xazim)
#ifdef SFX_MPI
    xtime_calc(2) = xtime_calc(2) + (mpi_wtime() - xtime1)
    xtime1 = mpi_wtime()
#endif
    !interpolation between beginning and end of current forcing time step
    CALL ol_time_interp_atm(jsurf_step,inb_atm,isize_omp,            &
                            zta(:,id_forc),zta(:,id_forc+1),         &
                            zqa(:,id_forc),zqa(:,id_forc+1),         &
                            zwind(:,id_forc),zwind(:,id_forc+1),     &
                            zdir_sw(:,id_forc),zdir_sw(:,id_forc+1), &
                            zsca_sw(:,id_forc),zsca_sw(:,id_forc+1), &
                            zlw(:,id_forc),zlw(:,id_forc+1),         &
                            zsnow(:,id_forc+1),zrain(:,id_forc+1),   &
                            zps(:,id_forc),zps(:,id_forc+1),         &
                            zco2(:,id_forc), zco2(:,id_forc+1),      &
                            zdir(:,id_forc) ,zdir(:,id_forc+1)       )
#ifdef SFX_MPI
    xtime_calc(3) = xtime_calc(3) + (mpi_wtime() - xtime1)
    xtime1 = mpi_wtime()
#endif
    !
    IF(ladapt_sw)THEN
      !
      ! coherence between solar zenithal angle and radiation
      ! when solar beam close to horizontal -> reduction of direct radiation to
      ! the benefit of scattered radiation
      ! when pi/2 - 0.1 < ZENITH < pi/2 - 0.05 => weight of direct to scattered radiation decreases linearly with zenith 
      ! when pi/2 - 0.05 < ZENITH => all the direct radiation is converted to scattered radiation
      ! coherence between solar zenithal angle and radiation
      !
      zcoef(:) = (xpi/2. - xzenith(:) - 0.05) / 0.05
      zcoef(:) = max(min(zcoef,1.),0.)
      DO jloop=1,SIZE(xdir_sw,2)
        xsca_sw(:,jloop) = xsca_sw(:,jloop) + xdir_sw(:,jloop) * (1 - zcoef)
        xdir_sw(:,jloop) = xdir_sw(:,jloop) * zcoef(:)
      ENDDO
      !
    ELSE
      !
      zsw(:) = 0.
      DO jloop=1,SIZE(xdir_sw,2)
        zsw(:) = zsw(:) + xdir_sw(:,jloop) + xsca_sw(:,jloop)
      END DO
      WHERE (zsw(:)>0.)
        xzenith  = min(xzenith ,xpi/2.-0.01)
        xzenith2 = min(xzenith2,xpi/2.-0.01)
      ELSEWHERE
        xzenith  = max(xzenith ,xpi/2.)
        xzenith2 = max(xzenith2,xpi/2.)
      END WHERE
      !
    ENDIF       
    !
    ! updates time
    ztimec= ztimec+xtstep_surf
    IF (lcoupl_topd) ltopd_step = ( mod((((jforc_step-1)*inb_atm)+jsurf_step),nnb_topd) == 0 )
    !
    ! run Surface
    !
#ifdef SFX_MPI
    xtime_calc(4) = xtime_calc(4) + (mpi_wtime() - xtime1)
#endif
    !
    CALL io_buff_clean
    !
    IF(loasis)THEN
     ! Receive fields to other models proc by proc
     CALL sfx_oasis_recv_ol(ysurf_cur%FM%F, ysurf_cur%IM%I, ysurf_cur%SM%S, ysurf_cur%U, ysurf_cur%WM%W, &
                            csurf_filetype,ini,ibands,ztimec,       &
                            xtstep_surf,isize_omp,xzenith,xsw_bands,&
                            xtsrad,xdir_alb,xsca_alb,xemis,xtsurf   )
    ENDIF
    !
!$OMP PARALLEL PRIVATE(INKPROMA,XTIME1)
    !
#ifdef SFX_MPI
    xtime1 = mpi_wtime()
#endif
    !
!$ NBLOCK = OMP_GET_THREAD_NUM()
    !
    IF (nblock==nblocktot) THEN
      CALL init_dim(isize_omp,0,inkproma,nindx1sfx,nindx2sfx)
    ELSE
      CALL init_dim(isize_omp,nblock,inkproma,nindx1sfx,nindx2sfx)
    ENDIF
    !
    IF (nblock==0) THEN
      CALL goto_model(nblocktot)
    ELSE
      CALL goto_model(nblock)
    ENDIF
    !
    IF(gshadows) THEN 
      CALL slope_radiative_effect(xtstep_surf,xzenith(nindx1sfx:nindx2sfx),&
      xazim(nindx1sfx:nindx2sfx),xps(nindx1sfx:nindx2sfx),xta(nindx1sfx:nindx2sfx),&
      xrain(nindx1sfx:nindx2sfx),xdir_sw(nindx1sfx:nindx2sfx,:),xlw(nindx1sfx:nindx2sfx),&
      xzs_thread,xzs_xy_thread, xslopang_thread,xslopazi_thread,xsurf_triangle_thread)
    END IF
    !
    CALL coupling_surf_atm_n(ysurf_cur, &
                             csurf_filetype, 'E', ztimec,                    &
           xtstep_surf, iyear, imonth, iday, ztime, inkproma, nscal, ibands, &
           xtsun(nindx1sfx:nindx2sfx), xzenith(nindx1sfx:nindx2sfx),         &
           xzenith2(nindx1sfx:nindx2sfx), xazim(nindx1sfx:nindx2sfx),        &
           xzref(nindx1sfx:nindx2sfx), xuref(nindx1sfx:nindx2sfx),           &
           xzs(nindx1sfx:nindx2sfx), xu(nindx1sfx:nindx2sfx),                &
           xv(nindx1sfx:nindx2sfx), xqa(nindx1sfx:nindx2sfx),                &
           xta(nindx1sfx:nindx2sfx), xrhoa(nindx1sfx:nindx2sfx),             &
           xsv(nindx1sfx:nindx2sfx,:), xco2(nindx1sfx:nindx2sfx), csv,       &
           xrain(nindx1sfx:nindx2sfx),  xsnow(nindx1sfx:nindx2sfx),          &
           xlw(nindx1sfx:nindx2sfx), xdir_sw(nindx1sfx:nindx2sfx,:),          &
           xsca_sw(nindx1sfx:nindx2sfx,:), xsw_bands, xps(nindx1sfx:nindx2sfx),&
           xpa(nindx1sfx:nindx2sfx), xsftq(nindx1sfx:nindx2sfx),             &
           xsfth(nindx1sfx:nindx2sfx), xsfts(nindx1sfx:nindx2sfx,:),         &
           xsfco2(nindx1sfx:nindx2sfx), xsfu(nindx1sfx:nindx2sfx),           &
           xsfv(nindx1sfx:nindx2sfx), xtsrad(nindx1sfx:nindx2sfx),           &
           xdir_alb(nindx1sfx:nindx2sfx,:), xsca_alb(nindx1sfx:nindx2sfx,:), &
           xemis(nindx1sfx:nindx2sfx),                                       &
           xtsurf(nindx1sfx:nindx2sfx), xz0(nindx1sfx:nindx2sfx),            &
           xz0h(nindx1sfx:nindx2sfx), xqsurf(nindx1sfx:nindx2sfx),           &          
           xpew_a_coef(nindx1sfx:nindx2sfx),xpew_b_coef(nindx1sfx:nindx2sfx),&
           xpet_a_coef(nindx1sfx:nindx2sfx),xpeq_a_coef(nindx1sfx:nindx2sfx),&
           xpet_b_coef(nindx1sfx:nindx2sfx),xpeq_b_coef(nindx1sfx:nindx2sfx),&
           ytest                                                             )
    !
    CALL reset_dim(ini,inkproma,nindx1sfx,nindx2sfx)
    !
#ifdef SFX_MPI
    xtime_calc(5) = xtime_calc(5) + (mpi_wtime() - xtime1)
#endif
    !
!$OMP END PARALLEL
    !
#ifdef SFX_MPI
    xtime1 = mpi_wtime()
#endif
    !
    IF(loasis)THEN
     ! Send fields to other models proc by proc
     CALL sfx_oasis_send_ol(ysurf_cur%FM%F, ysurf_cur%IM%I, ysurf_cur%SM%S, ysurf_cur%U, ysurf_cur%WM%W, &
                            csurf_filetype,ini,ztimec,xtstep_surf,isize_omp)
    ENDIF
    !
    ztime = ztime + xtstep_surf
    CALL add_forecast_to_date_surf(iyear, imonth, iday, ztime)
#ifdef SFX_MPI
    xtime_calc(6) = xtime_calc(6) + (mpi_wtime() - xtime1)
    !
    xtime1 =  mpi_wtime()
#endif
    ! ecrit Surface
    !
    IF ( lcoupl_topd .AND. ltopd_step ) THEN
      !
      IF (.NOT.ALLOCATED(nyear))  ALLOCATE(nyear(nnb_topd_step))
      IF (.NOT.ALLOCATED(nmonth)) ALLOCATE(nmonth(nnb_topd_step))
      IF (.NOT.ALLOCATED(nday))   ALLOCATE(nday(nnb_topd_step))
      IF (.NOT.ALLOCATED(nh))     ALLOCATE(nh(nnb_topd_step))
      IF (.NOT.ALLOCATED(nm))     ALLOCATE(nm(nnb_topd_step))
      !
      nyear(ntopd_step) = iyear
      nmonth(ntopd_step) = imonth
      nday(ntopd_step) = iday
      nh(ntopd_step) = int(ztime/3600.)
      nm(ntopd_step) = int((ztime-nh(ntopd_step)*3600.)/60.)
      !
      IF ( nm(ntopd_step)==60 ) THEN
        !
        nm(ntopd_step) = 0
        nh(ntopd_step) = nh(ntopd_step)+1
        !
      ENDIF
      !
      IF ( nh(ntopd_step)==24 ) THEN
        !
        nh(ntopd_step) = 0
        nday(ntopd_step) = nday(ntopd_step)+1
        !
        !!AJOUT BEC 
        SELECT CASE (nmonth(ntopd_step))
          CASE(4,6,9,11)
            IF ( nday(ntopd_step)==31 ) THEN
              nmonth(ntopd_step) = nmonth(ntopd_step)+1
              nday(ntopd_step) = 1
            ENDIF
          CASE(1,3,5,7:8,10)
            IF ( nday(ntopd_step)==32 ) THEN
              nmonth(ntopd_step) = nmonth(ntopd_step)+1
              nday(ntopd_step) = 1
            ENDIF
          CASE(12)
            IF ( nday(ntopd_step)==32 ) THEN
              nyear(ntopd_step) = nyear(ntopd_step)+1
              nmonth(ntopd_step) = 1
              nday(ntopd_step) = 1
            ENDIF
          CASE(2)
            IF( mod(nyear(ntopd_step),4)==0 .AND. mod(nyear(ntopd_step),100)/=0 .OR. mod(nyear(ntopd_step),400)==0 ) THEN 
              IF (nday(ntopd_step)==30) THEN
                nmonth(ntopd_step) = nmonth(ntopd_step)+1
                nday(ntopd_step) = 1
              ENDIF
            ELSE
              IF (nday(ntopd_step)==29) THEN
                nmonth(ntopd_step) = nmonth(ntopd_step)+1
                nday(ntopd_step) = 1
              ENDIF
            ENDIF
        END SELECT
        !
      ENDIF
      !
      ! * 2. Stocking date of each time step
      !
      ntopd_step = ntopd_step + 1 
      !     
    ENDIF
    !
    IF (mod(ztimec,xtstep_output) == 0. .AND. ctimeseries_filetype/='NONE  ') THEN
      !
      IF (nrank==npio) THEN
        !
        !* name of the file
        IF (ctimeseries_filetype=="ASCII " .OR. &
            ctimeseries_filetype=="LFI   " .OR. &
            ctimeseries_filetype=="FA    " .OR. &
            ctimeseries_filetype=="NC    "    ) THEN  
          !
          !
          ztime_out  = ztime
          iday_out   = iday
          imonth_out = imonth
          iyear_out  = iyear
          !
          IF(lout_timename)THEN
            ! if true, change the name of output file at the end of a day
            ! (ex: 19860502_00h00 -> 19860501_24h00)                     
            IF(ztime==0.0)THEN
              ztime_out = 86400.
              iday_out   = iday-1
              IF(iday_out==0)THEN
                imonth_out = imonth - 1
                IF(imonth_out==0)THEN
                  imonth_out=12
                  iyear_out = iyear - 1
                ENDIF
                SELECT CASE (imonth_out)
                  CASE(4,6,9,11)
                    iday_out=30
                  CASE(1,3,5,7:8,10,12)
                    iday_out=31
                  CASE(2)
                    IF( ((mod(iyear_out,4)==0).AND.(mod(iyear_out,100)/=0)) .OR. (mod(iyear_out,400)==0))THEN 
                      iday_out=29
                    ELSE
                      iday_out=28
                    ENDIF
                END SELECT
              ENDIF
            ENDIF
            !
          ENDIF
          !
          WRITE(ytag,fmt='(I4.4,I2.2,I2.2,A1,I2.2,A1,I2.2)') iyear_out,imonth_out,iday_out,&
                '_',int(ztime_out/3600.),'h',nint(ztime_out)/60-60*int(ztime_out/3600.)  
          cfileout    = adjustl(adjustr(csurffile)//'.'//ytag//'.txt')
          cfileout_lfi= adjustl(adjustr(csurffile)//'.'//ytag)
          cfileout_fa = adjustl(adjustr(csurffile)//'.'//ytag//'.fa')
          cfileout_nc = adjustl(adjustr(csurffile)//'.'//ytag//'.nc')
          !
          IF (ctimeseries_filetype=='FA    ') THEN
#ifdef SFX_FA
            lfanocompact = ldiag_fa_nocompact
            idatef(1)= iyear!_OUT
            idatef(2)= imonth!_OUT
            idatef(3)= iday!_OUT
            !ZTIME instead of ZTIME_OUT (FA XRD39 do not like 24h)
            idatef(4)= floor(ztime/3600.)
            idatef(5)= floor(ztime/60.) - idatef(4) * 60 
            idatef(6)= nint(ztime) - idatef(4) * 3600 - idatef(5) * 60
            idatef(7:11) = 0
            IF (csurf_filetype/='FA    ') THEN
              CALL write_header_fa(ysurf_cur%UG, &
                                   csurf_filetype,'ALL')
            ELSE
              CALL faitou(iret,nunit_fa,.true.,cfileout_fa,'UNKNOWN',.true.,.false.,iverbfa,0,inb,cdnomc)
            ENDIF
            CALL fandar(iret,nunit_fa,idatef)
#endif
          END IF
          !
        END IF
        !
        xstartw = xstartw + 1
        nwrite  = nwrite  + 1
        ltime_written(:)=.false.
        !
      ENDIF
      !
#ifdef SFX_MPI
      xtime_write(1) = xtime_write(1) + (mpi_wtime() - xtime1)
#endif
      !
      ldef = .true.
      !
      IF (ctimeseries_filetype=="NC    ") THEN
        CALL goto_model(1)
        CALL init_output_nc_n(ysurf_cur%TM%BDD, ysurf_cur%CHE, ysurf_cur%CHN, ysurf_cur%CHU, &
                               ysurf_cur%SM%DTS, ysurf_cur%TM%DTT, ysurf_cur%DTZ, ysurf_cur%IM%I, &
                               ysurf_cur%UG, ysurf_cur%U, ysurf_cur%DGU)                   
      ENDIF
      !
      idx_w = 0
      !
      DO jnw = 1,inw
        ! 
        CALL io_buff_clean
        !
!$OMP PARALLEL PRIVATE(INKPROMA,XTIME1)
        !
!$ NBLOCK = OMP_GET_THREAD_NUM()
        !
        IF (nblock==nblocktot) THEN
          CALL init_dim(isize_omp,0,inkproma,nindx1sfx,nindx2sfx)
        ELSE
          CALL init_dim(isize_omp,nblock,inkproma,nindx1sfx,nindx2sfx)
        ENDIF
        !
        IF (nblock==0) THEN
          CALL goto_model(nblocktot)
        ELSE
          CALL goto_model(nblock)
        ENDIF
        !
#ifdef SFX_MPI
        xtime1 =  mpi_wtime()
#endif
        CALL write_surf_atm_n(ysurf_cur, &
                              ctimeseries_filetype,'ALL',lland_use)
#ifdef SFX_MPI
        xtime_write(2) = xtime_write(2) + (mpi_wtime() - xtime1)
        xtime1 =  mpi_wtime()
#endif
        CALL diag_surf_atm_n(ysurf_cur%IM%DGEI, ysurf_cur%FM%DGF, ysurf_cur%DGL, ysurf_cur%IM%DGI, &
                             ysurf_cur%SM%DGS, ysurf_cur%DGU, ysurf_cur%TM%DGT, ysurf_cur%WM%DGW, &
                             ysurf_cur%U, ysurf_cur%USS, &
                             ctimeseries_filetype)
#ifdef SFX_MPI
        xtime_write(3) = xtime_write(3) + (mpi_wtime() - xtime1)
        xtime1 =  mpi_wtime()
#endif
        CALL write_diag_surf_atm_n(ysurf_cur, &
                                   ctimeseries_filetype,'ALL')
#ifdef SFX_MPI
        xtime_write(4) = xtime_write(4) + (mpi_wtime() - xtime1)
#endif
        !
        CALL reset_dim(ini,inkproma,nindx1sfx,nindx2sfx)
        !
!$OMP BARRIER
        !
!$OMP END PARALLEL
        !
        ldef = .false.
        !
      ENDDO
      !
      IF (lcoupl_topd .AND. ntopd_step > nnb_topd_step) THEN
        !
        ! Writing of file resulting of coupling with TOPMODEL or routing ****
        CALL write_discharge_file(csurf_filetype,'q_total.txt','FORMATTED',&
                                  nyear,nmonth,nday,nh,nm,xqtot)
        CALL write_discharge_file(csurf_filetype,'q_runoff.txt','FORMATTED',&
                                  nyear,nmonth,nday,nh,nm,xqb_run)
        CALL write_discharge_file(csurf_filetype,'q_drainage.txt','FORMATTED',&
                                  nyear,nmonth,nday,nh,nm,xqb_dr)
        ! Writing of budget files 
        IF (lbudget_topd) CALL write_budget_coupl_rout
        !
      ENDIF
      !
#ifdef SFX_MPI
      xtime1 =  mpi_wtime()
#endif
      !
      IF (nrank==npio) THEN
        IF (ctimeseries_filetype=='FA    ') THEN
#ifdef SFX_FA
          CALL fairme(iret,nunit_fa,'UNKNOWN')
#endif
        END IF
        !* add informations in the file
        IF (ctimeseries_filetype=='LFI   ' .AND. lmnh_compatible) CALL write_header_mnh
      ENDIF
#ifdef SFX_MPI
      xtime_write(5) = xtime_write(5) + (mpi_wtime() - xtime1)
#endif
      !
    ENDIF
    !
  END DO
  !
  
  IF (nrank==npio) THEN
    IF (lprint) THEN
      IF (mod(ztimec,xday) == 0.) THEN
!$OMP SINGLE              
        icount = icount + 1
        CALL wlog_mpi('SFX DAY :',klog=icount,klog2=int(zduration/xday))
        WRITE(*,'(A10,I5,A2,I5)')'SFX  DAY :',icount,' /',int(zduration/xday)
!$OMP END SINGLE        
      ENDIF
    ENDIF
  ENDIF
  !
END DO
!
!$OMP PARALLEL PRIVATE(XTIME)
!
#ifdef SFX_MPI
xtime = (mpi_wtime() - xtime0)
#endif
!
 CALL wlog_mpi(' ')
 CALL wlog_mpi('OL_READ_ATM ',plog=xtime_calc(1))
 CALL wlog_mpi('SUNPOS ',plog=xtime_calc(2))
 CALL wlog_mpi('OL_TIME_INTERP_ATM ',plog=xtime_calc(3))
 CALL wlog_mpi('')
 CALL wlog_mpi('ZENITH ',plog=xtime_calc(4))
 CALL wlog_mpi('')
 CALL wlog_mpi('COUPLING_SURF_ATM ',plog=xtime_calc(5))
 CALL wlog_mpi('')
 CALL wlog_mpi('ADD_FORECAST_TO_DATE_SURF ',plog=xtime_calc(6))
 CALL wlog_mpi('DEF_DATE ',plog=xtime_write(1))
 CALL wlog_mpi('')
 CALL wlog_mpi('WRITE_SURF_ATM ',plog=xtime_write(2))
 CALL wlog_mpi('DIAG_SURF_ATM ',plog=xtime_write(3))
 CALL wlog_mpi('WRITE_DIAG_SURF_ATM ',plog=xtime_write(4))
 CALL wlog_mpi('')
 CALL wlog_mpi('CLOSE FILES ',plog=xtime_write(5))
 CALL wlog_mpi('')
 CALL wlog_mpi('END LOOP ',plog=xtime)
 CALL wlog_mpi('')
 CALL wlog_mpi('TIME_NPIO_WRITE ',plog=xtime_npio_write)
 CALL wlog_mpi('TIME_COMM_WRITE ',plog=xtime_comm_write)
 CALL wlog_mpi('TIME_OMP_BARR ',plog=xtime_omp_barr)
 CALL wlog_mpi('TIME_CALC_WRITE ',plog=xtime_calc_write)
 CALL wlog_mpi('')
 CALL wlog_mpi('TIME_INIT_SEA ',plog=xtime_init_sea)
 CALL wlog_mpi('TIME_INIT_WATER ',plog=xtime_init_water)
 CALL wlog_mpi('TIME_INIT_NATURE ',plog=xtime_init_nature)
 CALL wlog_mpi('TIME_INIT_TOWN ',plog=xtime_init_town)
 CALL wlog_mpi('')
 CALL wlog_mpi('TIME_SEA ',plog=xtime_sea)
 CALL wlog_mpi('TIME_WATER ',plog=xtime_water)
 CALL wlog_mpi('TIME_NATURE ',plog=xtime_nature)
 CALL wlog_mpi('TIME_TOWN ',plog=xtime_town)
!
!$OMP END PARALLEL
!
!
IF (cforcing_filetype=='ASCII ' .OR. cforcing_filetype=='BINARY') CALL open_close_bin_asc_forc('CLOSE',cforcing_filetype,'R')
!
IF (cforcing_filetype=='NETCDF') CALL close_filein_ol
IF (ctimeseries_filetype=='OFFLIN') CALL close_fileout_ol
!
! --------------------------------------------------------------------------------------
!
!*    3.     write restart file
!            ------------------
!
IF ( lrestart ) THEN
  !
  IF (nrank==npio) THEN
    !* name of the file
    cfileout    = adjustl(adjustr(csurffile)//'.txt')
    cfileout_lfi= csurffile
    cfileout_fa = adjustl(adjustr(csurffile)//'.fa')
    cfileout_nc = adjustl(adjustr(csurffile)//'.nc')

    !* opens the file
    IF (csurf_filetype=='FA    ') THEN
#ifdef SFX_FA
      lfanocompact = .true.
      idatef(1)= iyear
      idatef(2)= imonth
      idatef(3)= iday
      idatef(4)= floor(ztime/3600.)
      idatef(5)= floor(ztime/60.) - idatef(4) * 60 
      idatef(6)= nint(ztime) - idatef(4) * 3600 - idatef(5) * 60
      idatef(7:11) = 0    
      CALL faitou(iret,nunit_fa,.true.,cfileout_fa,'UNKNOWN',.true.,.false.,iverbfa,0,inb,cdnomc)
      CALL fandar(iret,nunit_fa,idatef)
#endif
    END IF
    !
  ENDIF
  !
  inw = 1
  IF (csurf_filetype=="NC    ") inw = 2
  !
  ldef = .true.
  !
  IF (csurf_filetype=="NC    ") THEN
    CALL goto_model(1)
    CALL init_output_nc_n(ysurf_cur%TM%BDD, ysurf_cur%CHE, ysurf_cur%CHN, ysurf_cur%CHU, &
                         ysurf_cur%SM%DTS, ysurf_cur%TM%DTT, ysurf_cur%DTZ, ysurf_cur%IM%I, &
                         ysurf_cur%UG, ysurf_cur%U, ysurf_cur%DGU)
  ENDIF
  !
  DO jnw = 1,inw
    !
    CALL io_buff_clean
    !
!$OMP PARALLEL PRIVATE(INKPROMA,JNW)
  !
!$ NBLOCK = OMP_GET_THREAD_NUM()
    !
    IF (nblock==nblocktot) THEN
      CALL init_dim(isize_omp,0,inkproma,nindx1sfx,nindx2sfx)
    ELSE
      CALL init_dim(isize_omp,nblock,inkproma,nindx1sfx,nindx2sfx)
    ENDIF
    !
    IF (nblock==0) THEN
      CALL goto_model(nblocktot)
    ELSE
      CALL goto_model(nblock)
    ENDIF
    !  
    CALL flag_update(ysurf_cur%IM%DGI, ysurf_cur%DGU, &
                     .false.,.true.,.false.,.false.)
    !
    IF (lrestart_2m) THEN
      i2m       = 1
      gpgd_isba = .true.
    ELSE
      i2m       = 0
      gpgd_isba = .false.
    ENDIF
    gfrac                  = .true.
    gdiag_grid             = .true.
    gsurf_budget           = .false.
    grad_budget            = .false.
    gcoef                  = .false.
    gsurf_vars             = .false.
    ibeq                   = 0
    idsteq                 = 0
    gdiag_ocean            = .false.
    gdiag_seaice           = .false.
    gwater_profile         = .false.
    gsurf_evap_budget      = .false.
    gflood                 = .false.
    gch_no_flux_isba       = .false.
    gsurf_misc_budget_isba = .false.
    gpgd_teb               = .false.
    gsurf_misc_budget_teb  = .false.  
    !
    CALL flag_diag_update(ysurf_cur%FM%CHF, ysurf_cur%IM%CHI, ysurf_cur%SM%CHS, ysurf_cur%TM%CHT, &
                          ysurf_cur%WM%CHW, ysurf_cur%IM%DGEI, ysurf_cur%FM%DGF, ysurf_cur%IM%DGI, &
                          ysurf_cur%FM%DGMF, ysurf_cur%IM%DGMI, ysurf_cur%TM%DGMTO, ysurf_cur%SM%DGO, &
                          ysurf_cur%SM%DGS, ysurf_cur%SM%DGSI, ysurf_cur%DGU, ysurf_cur%TM%DGT, &
                          ysurf_cur%WM%DGW, ysurf_cur%IM%I, ysurf_cur%U, &
                          gfrac, gdiag_grid, i2m, gsurf_budget, grad_budget, gcoef,  &
                          gsurf_vars, ibeq, idsteq, gdiag_ocean, gdiag_seaice,       &
                          gwater_profile,                                            &                          
                          gsurf_evap_budget, gflood,  gpgd_isba, gch_no_flux_isba,   &
                          gsurf_misc_budget_isba, gpgd_teb, gsurf_misc_budget_teb    )
    !
    !* writes into the file
    CALL write_surf_atm_n(ysurf_cur, &
                              csurf_filetype,'ALL',lland_use)
    IF(csurf_filetype/='FA    ' .OR. lrestart_2m) THEN
       CALL write_diag_surf_atm_n(ysurf_cur, &
                                   csurf_filetype,'ALL')
    ENDIF
    !
    CALL reset_dim(ini,inkproma,nindx1sfx,nindx2sfx)
    !
!$OMP END PARALLEL     
    !
    ldef = .false.
    !
  ENDDO
  !
  !* closes the file
  IF (nrank==0 ) THEN
    IF (csurf_filetype=='FA    ') THEN
#ifdef SFX_FA
      CALL fairme(iret,nunit_fa,'UNKNOWN')
#endif
    END IF
    !* add informations in the file
    IF (csurf_filetype=='LFI   ' .AND. lmnh_compatible) CALL write_header_mnh
    !
  ENDIF
  !
  IF (lcoupl_topd .AND. ntopd_step > nnb_topd_step) CALL prep_restart_coupl_topd(ysurf_cur%UG, ysurf_cur%U, &
                                                                                 csurf_filetype,ini)
  !
END IF
!
! --------------------------------------------------------------------------------------
!
!*    4.     inquiry mode
!            ------------
!
IF ( linquire ) THEN
  !
!$OMP PARALLEL PRIVATE(ZSEA,ZWATER,ZNATURE,ZTOWN,ZT2M,ZQ2M,ZZ0,ZZ0H,ZQS_SEA, &
!$OMP ZQS_WATER,ZQS_NATURE,ZQS_TOWN,ZQS,ZPSNG,ZPSNV,ZZ0EFF,ZZS ,&
!$OMP INKPROMA,INI)
!
!$ NBLOCK = OMP_GET_THREAD_NUM()
!
  IF (nblock==nblocktot) THEN
    CALL init_dim(isize_omp,0,inkproma,nindx1sfx,nindx2sfx)
  ELSE
    CALL init_dim(isize_omp,nblock,inkproma,nindx1sfx,nindx2sfx)
  ENDIF
  ini = nindx2sfx-nindx1sfx+1
  !
  ALLOCATE( zsea( ini ) )
  ALLOCATE( zwater( ini ) )
  ALLOCATE( znature( ini ) )
  ALLOCATE( ztown( ini ) )
  ALLOCATE( zt2m( ini ) )
  ALLOCATE( zq2m( ini ) )
  ALLOCATE( zz0( ini ) )
  ALLOCATE( zz0h( ini ) )
  ALLOCATE( zqs_sea( ini ) )
  ALLOCATE( zqs_water( ini ) )
  ALLOCATE( zqs_nature( ini ) )
  ALLOCATE( zqs_town( ini ) )
  ALLOCATE( zqs( ini ) )
  ALLOCATE( zpsng( ini ) )
  ALLOCATE( zpsnv( ini ) )
  ALLOCATE( zz0eff( ini ) )
  ALLOCATE( zzs( ini ) )
  !
  iseries = 0
  CALL get_surf_var_n(ysurf_cur%FM%DGF, ysurf_cur%IM%I, ysurf_cur%IM%DGI, ysurf_cur%IM%DGMI, &
                      ysurf_cur%SM%DGS, ysurf_cur%DGU, ysurf_cur%TM%DGT, ysurf_cur%WM%DGW, &
                      ysurf_cur%FM%F, ysurf_cur%UG, ysurf_cur%U, ysurf_cur%USS, &
                      csurf_filetype,ini,iseries,psea=zsea,pwater=zwater,pnature=znature,ptown=ztown, &
                        pt2m=zt2m,pq2m=zq2m,pqs=zqs,pz0=zz0,pz0h=zz0h,pz0eff=zz0eff,pqs_sea=zqs_sea,  &
                        pqs_water=zqs_water,pqs_nature=zqs_nature,pqs_town=zqs_town,                  &
                        ppsng=zpsng,ppsnv=zpsnv,pzs=zzs                                         )  
  !
  isize = SIZE(nindex)
  IF (nrank==npio) THEN
!$OMP SINGLE
    ALLOCATE(zsea_full(isize))
    ALLOCATE(zwater_full(isize))
    ALLOCATE(znature_full(isize))
    ALLOCATE(ztown_full(isize))
    ALLOCATE(zz0_full(isize))
    ALLOCATE(zz0eff_full(isize))
    ALLOCATE(zzs_full(isize))
!$OMP END SINGLE
  ELSE
!$OMP SINGLE          
    ALLOCATE(zsea_full(0))
    ALLOCATE(zwater_full(0))
    ALLOCATE(znature_full(0))
    ALLOCATE(ztown_full(0))
    ALLOCATE(zz0_full(0))
    ALLOCATE(zz0eff_full(0))
    ALLOCATE(zzs_full(0))
!$OMP END SINGLE    
  ENDIF
  CALL gather_and_write_mpi(zsea,zsea_full) 
  CALL gather_and_write_mpi(zwater,zwater_full)
  CALL gather_and_write_mpi(znature,znature_full)
  CALL gather_and_write_mpi(ztown,ztown_full)
  CALL gather_and_write_mpi(zz0,zz0_full)
  CALL gather_and_write_mpi(zz0eff,zz0eff_full)
  CALL gather_and_write_mpi(zzs,zzs_full)

  IF (nrank==npio) THEN
!$OMP SINGLE
    WRITE(iluout,'(A32,I4,A3,I4)') ' GRID BOXES CONTAINING SEA    : ',count( zsea_full(:) > 0. ),' / ',isize
    WRITE(iluout,'(A32,I4,A3,I4)') ' GRID BOXES CONTAINING WATER  : ',count( zwater_full(:) > 0. ),' / ',isize
    WRITE(iluout,'(A32,I4,A3,I4)') ' GRID BOXES CONTAINING NATURE : ',count( znature_full(:) > 0. ),' / ',isize
    WRITE(iluout,'(A32,I4,A3,I4)') ' GRID BOXES CONTAINING TOWN   : ',count( ztown_full(:) > 0. ),' / ',isize
    WRITE(iluout,*)'ZZ0    = ',zz0_full
    WRITE(iluout,*)'ZZ0EFF = ',zz0eff_full
    WRITE(iluout,*)'ZZS = ',zzs_full
    WRITE(iluout,*)'MINVAL(ZZS) = ',minval(zzs_full),' MAXVAL(ZZS) = ',maxval(zzs_full)
!$OMP END SINGLE
  ENDIF
  !
  DEALLOCATE( zsea       )
  DEALLOCATE( zwater     )
  DEALLOCATE( znature    )
  DEALLOCATE( ztown      )
  DEALLOCATE( zt2m       )
  DEALLOCATE( zq2m       )
  DEALLOCATE( zz0        )
  DEALLOCATE( zz0h       )
  DEALLOCATE( zqs_sea    )
  DEALLOCATE( zqs_water  )
  DEALLOCATE( zqs_nature )
  DEALLOCATE( zqs_town   )
  DEALLOCATE( zqs        )
  DEALLOCATE( zpsng      )
  DEALLOCATE( zpsnv      )
  DEALLOCATE( zz0eff     )
  DEALLOCATE( zzs        )
  !
  IF (nrank==npio) THEN
!$OMP SINGLE
    DEALLOCATE(zsea_full   )
    DEALLOCATE(zwater_full )
    DEALLOCATE(znature_full)
    DEALLOCATE(ztown_full  )
    DEALLOCATE(zz0_full    )
    DEALLOCATE(zz0eff_full )
    DEALLOCATE(zzs_full    )
!$OMP END SINGLE
  ENDIF
  !
!$OMP END PARALLEL
  !
ENDIF   
!
! --------------------------------------------------------------------------------------
!
!    4'    Close Gelato specific diagnostic 
#if ! defined in_arpege
 CALL clsdia()
#endif
!
!
!*    5.     Close parallelized I/O
!            ----------------------
!
IF (nrank==npio) THEN
  WRITE(iluout,*) ' '
  WRITE(iluout,*) '    --------------------------'
  WRITE(iluout,*) '    | OFFLINE ENDS CORRECTLY |'
  WRITE(iluout,*) '    --------------------------'
  WRITE(iluout,*) ' '
!$OMP SINGLE
  CLOSE(iluout)
!$OMP END SINGLE  
  WRITE(*,*) ' '
  WRITE(*,*) '    --------------------------'
  WRITE(*,*) '    | OFFLINE ENDS CORRECTLY |'
  WRITE(*,*) '    --------------------------'
  WRITE(*,*) ' '
ENDIF
!
 CALL surfex_deallo_list
!
IF (ALLOCATED(nindex)) DEALLOCATE(nindex)
IF (ALLOCATED(nsize_task)) DEALLOCATE(nsize_task)
!
IF (ASSOCIATED(nwork)) DEALLOCATE(nwork)
IF (ASSOCIATED(xwork)) DEALLOCATE(xwork)
IF (ASSOCIATED(nwork2)) DEALLOCATE(nwork2)
IF (ASSOCIATED(xwork2)) DEALLOCATE(xwork2)
IF (ASSOCIATED(xwork3)) DEALLOCATE(xwork3)
IF (ASSOCIATED(nwork_full)) DEALLOCATE(nwork_full)
IF (ASSOCIATED(xwork_full)) DEALLOCATE(xwork_full)
IF (ASSOCIATED(nwork2_full)) DEALLOCATE(nwork2_full)
IF (ASSOCIATED(xwork2_full)) DEALLOCATE(xwork2_full)
!
 CALL end_log_mpi
!
IF (lhook) CALL dr_hook('OFFLINE',1,zhook_handle)
!
! * MPI and OASIS must be finalized after the last DR_HOOK call
!
IF(loasis)THEN
  CALL sfx_oasis_end
ELSE
#ifdef SFX_MPI
 CALL mpi_finalize(infompi)
#endif
ENDIF
!
! --------------------------------------------------------------------------------------
!
END PROGRAM offline