soda.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 soda

#ifdef USE_SODA
!
! ------------------------------------------------------------------------------------------
!!
!!    SODA: SURFEX Offline Data Assimilation
!!
!!    PURPOSE
!!    -------
!!    Program to perform surface data assimilation within SURFEX 
!!
!!
!!    METHOD
!!    ------
!!    Different methods for different tiles
!!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    T. Aspelien                  met.no
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original         04/2012
!!
!! 03/2014 E. Martin change indices names in OMP module according to GMAP changes
!  05/2013 B. Decharme New coupling variables XTSURF (for AGCM)
!  02/2016 B. Decharme MODD_IO_SURF_ARO not used
!----------------------------------------------------------------------------
!
USE modd_off_surfex_n
USE mode_modeln_surfex_handler
!
USE modd_surfex_mpi, ONLY : nrank, npio, wlog_mpi, prep_log_mpi, nproc, ncomm,   &
                            nindex, nsize_task, end_log_mpi, nsize
USE modd_surfex_omp, ONLY : nindx2sfx, nwork, xwork, xwork2, xwork3, &
                            nwork_full, xwork_full, xwork2_full, nwork2, &
                            nwork2_full
!
USE modd_mask, ONLY: nmask_full
!
USE modd_type_date_surf, ONLY : date_time
!
USE modd_surf_conf, ONLY : cprogname, csoftware
USE modd_surf_par,  ONLY : xundef,nundef
!
USE modd_assim, ONLY : lassim, larome, laladsurf, cassim_isba, nvar, xf, xf_patch, &
                       nobstype, xat2m_isba, xahu2m_isba, cvar, cobs, nechgu, xi,  &
                       xlai_pass, xbio_pass, cbio, nivar, xyo, nific, nprintlev, &
                       nobs, nprintlev, lread_all, nens,lbias_correction, &
                       lextrap_sea,lextrap_water,lextrap_nature,lobsheader,nobsmax, &
                       cfile_format_fg,cfile_format_lsm,cfile_format_obs,cfile_format_clim,&
                       nnco,lwatertg2,lobsnat
!
USE modd_forc_atm,       ONLY : csv, xdir_alb, xsca_alb, xemis, xtsrad, xtsun, xzs, &
                                xzref, xuref, xta, xqa, xsv, xu, xv, xsw_bands,     &
                                xzenith, xazim, xco2, xrhoa, xtsurf
!
USE modd_write_bin,  ONLY : nwrite
!
#ifdef SFX_OL
USE modd_io_surf_ol, ONLY : xstart, xcount, xstride, xstartw, xcountw, ltime_written,  &
                            ldefined_nature, ldefined_sea, ldefined_water, &
                            ldefined_town, ldefined_surf_atm, lpartw
                            
#endif
!
#ifdef SFX_NC
USE modd_io_surf_nc,   ONLY : cfilein_nc, cfilein_nc_save, cfilepgd_nc, cfileout_nc, ldef, &
                              cluout_nc
#endif
#ifdef SFX_ASC
USE modd_io_surf_asc,  ONLY : cfilein, cfilein_save, cfilepgd, cfileout, lcreated
#endif
#ifdef SFX_FA
USE modd_io_surf_fa,   ONLY : cfilein_fa, cfilein_fa_save, cfilepgd_fa, cdnomc, cfileout_fa, &
                              nunit_fa, iverbfa, lfanocompact
#endif
#ifdef SFX_LFI
USE modd_io_surf_lfi,    ONLY : cfilein_lfi, cfilein_lfi_save, &
                                cfilepgd_lfi, cfile_lfi, cluout_lfi, cfileout_lfi 
#endif
!
USE modn_io_offline,     ONLY : nam_io_offline, cnamelist, cpgdfile, cprepfile, csurffile, &
                                csurf_filetype, ctimeseries_filetype, lland_use, yalg_mpi, &
                                ldiag_fa_nocompact, lout_timename, xio_frac, lrestart_2m
!
USE mode_pos_surf,  ONLY : posnam
!
USE modi_dealloc_surf_atm_n
USE modi_init_output_ol_n
USE modi_set_surfex_filein
USE modi_init_index_mpi
USE modi_gather_and_write_mpi
USE modi_read_and_send_mpi
USE modi_abor1_sfx
USE modi_get_luout
USE modi_open_namelist
USE modi_close_namelist
USE modi_read_all_namelists
USE modi_init_io_surf_n
USE modi_end_io_surf_n
USE modi_read_surf
USE modi_io_buff_clean
USE modi_get_size_full_n
USE modi_init_surf_atm_n
USE modi_assim_surf_atm_n
USE modi_write_surf_atm_n
USE modi_write_diag_surf_atm_n
USE modi_assim_set_sst
USE modi_add_forecast_to_date_surf
USE modi_flag_update
USE modi_flag_diag_update
USE modi_test_nam_var_surf
USE modi_init_output_nc_n
USE modi_close_fileout_ol
!
USE mode_ekf, ONLY : get_file_name, set_filein
!
USE yomhook,             ONLY : lhook,dr_hook
USE parkind1,            ONLY : jprb
!
IMPLICIT NONE
!
#ifdef SFX_MPI
include 'mpif.h'
#endif
!
!*    0.     Declaration of local variables
!            ------------------------------
!
 CHARACTER(LEN=200) :: ymfile     ! Name of the observation, perturbed or reference file!
 CHARACTER(LEN=3)  :: yinit
 CHARACTER(LEN=2), PARAMETER  :: ytest        = 'OK'          ! must be equal to 'OK'
 CHARACTER(LEN=28)            :: yatmfile  ='   '  ! name of the Atmospheric file
 CHARACTER(LEN=6)             :: yatmfiletype ='      '                     ! type of the Atmospheric file
 CHARACTER(LEN=28)            :: yluout    ='LISTING_SODA                '  ! name of listing
 CHARACTER(LEN=6)             :: yprogram2 = 'FA    '
 CHARACTER(LEN=28)            :: yfilein
 CHARACTER(LEN=3)             :: yvar
!
 CHARACTER(LEN=100) :: yname
 CHARACTER(LEN=10)  :: yrank
 CHARACTER(LEN=3) :: yens
 CHARACTER(LEN=10),DIMENSION(:), ALLOCATABLE  :: cobsinfile     ! Identifier for simulated observations in file
!
REAL, ALLOCATABLE, DIMENSION(:,:) :: zyo_nat
REAL, ALLOCATABLE, DIMENSION(:) :: znature
!
REAL,ALLOCATABLE, DIMENSION(:,:) :: zwork
REAL,ALLOCATABLE, DIMENSION(:)   :: zlsm                ! Land-Sea mask
REAL,ALLOCATABLE, DIMENSION(:)   :: zcon_rain           ! Amount of convective liquid precipitation
REAL,ALLOCATABLE, DIMENSION(:)   :: zstrat_rain         ! Amount of stratiform liquid precipitation
REAL,ALLOCATABLE, DIMENSION(:)   :: zcon_snow           ! Amount of convective solid precipitation
REAL,ALLOCATABLE, DIMENSION(:)   :: zstrat_snow         ! Amount of stratiform solid precipitation
REAL,ALLOCATABLE, DIMENSION(:)   :: zclouds             ! Cloudcover
REAL,ALLOCATABLE, DIMENSION(:)   :: zevaptr             ! Evaporation
REAL,ALLOCATABLE, DIMENSION(:)   :: zevap               ! Evaporation
REAL,ALLOCATABLE, DIMENSION(:)   :: ztsc                ! Climatological surface temperature
REAL,ALLOCATABLE, DIMENSION(:)   :: zts                 ! Surface temperature
REAL,ALLOCATABLE, DIMENSION(:)   :: zt2m                ! Screen level temperature
REAL,ALLOCATABLE, DIMENSION(:)   :: zhu2m               ! Screen level relative humidity
REAL,ALLOCATABLE, DIMENSION(:)   :: zsnc
REAL,ALLOCATABLE, DIMENSION(:)   :: zswe                ! Snow water equvivalent (amount of snow on the ground)
REAL,ALLOCATABLE, DIMENSION(:)   :: zswec               ! Climatological snow water equvivalent (amount of snow on the ground)
REAL,ALLOCATABLE, DIMENSION(:)   :: zucls
REAL,ALLOCATABLE, DIMENSION(:)   :: zvcls
REAL,ALLOCATABLE, DIMENSION(:)   :: zsst                ! SST from external file
REAL,ALLOCATABLE, DIMENSION(:)   :: zsic                ! SIC from external file
REAL,ALLOCATABLE, DIMENSION(:)   :: zlat
REAL,ALLOCATABLE, DIMENSION(:)   :: zlon
!
REAL    :: ztime
REAL    :: ztime_out           ! output time since start of the run (s)
!
LOGICAL, ALLOCATABLE, DIMENSION(:) :: gd_maskext
LOGICAL :: glkeepextzone
LOGICAL :: gfound
!
TYPE (date_time)                 :: ttime               ! Current date and time  
!
 CHARACTER(LEN=14)                :: ytag
!
INTEGER, DIMENSION(11)  :: idatef
INTEGER :: idim_full
INTEGER :: isv                 ! Number of scalar species
INTEGER :: isw                 ! Number of radiative bands 
INTEGER :: iyear, imonth, iday, ihour
INTEGER :: iyear_out, imonth_out, iday_out
INTEGER :: jl,ji,jj,inb,icpt
INTEGER :: inw, jnw
INTEGER :: istep
INTEGER :: iobs
INTEGER :: igpcomp
INTEGER :: iluout
INTEGER :: ilunam
INTEGER :: iret, inbfa
INTEGER :: iresp, istat               ! Response value
INTEGER :: infompi, ilevel
INTEGER :: isize, iens, isize_full
!
! Flag diag :
!
INTEGER :: i2m, ibeq, idsteq
LOGICAL :: gfrac, gdiag_grid, gsurf_budget, grad_budget, gcoef,    &
           gsurf_vars, 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
!
REAL(KIND=JPRB)                  :: zhook_handle
! ******************************************************************************************
!
infompi=1
!
#ifdef SFX_MPI
 CALL mpi_init_thread(mpi_thread_multiple,ilevel,infompi)
#endif
!
IF (lhook) CALL dr_hook('SODA',0,zhook_handle)
!
 csoftware = 'SODA'
!
#ifdef SFX_MPI
ncomm = mpi_comm_world
 CALL mpi_comm_size(ncomm,nproc,infompi)
 CALL mpi_comm_rank(ncomm,nrank,infompi)
#endif
!
 CALL prep_log_mpi
!
!--------------------------------------
!
IF (nrank==npio) THEN
  WRITE(*,*)
  WRITE(*,*) '   ------------------------------------'
  WRITE(*,*) '   |               SODA               |'
  WRITE(*,*) '   | SURFEX OFFLINE DATA ASSIMILATION |'
  WRITE(*,*) '   ------------------------------------'
  WRITE(*,*)
ENDIF
!
WRITE(yrank,fmt='(I10)') nrank
yname=trim(yluout)//adjustl(yrank)
!
! Open ascii outputfile for writing
#ifdef SFX_LFI
 cluout_lfi =  adjustl(adjustr(yluout)//'.txt')
#endif
#ifdef SFX_NC
 cluout_nc = adjustl(adjustr(yluout)//'.txt')
#endif
 CALL get_luout('ASCII ',iluout)
OPEN(unit=iluout,file=adjustl(adjustr(yname)//'.txt'),form='FORMATTED',action='WRITE')
!
! Read offline specific things
 CALL open_namelist('ASCII ',ilunam,cnamelist)
 CALL posnam(ilunam,'NAM_IO_OFFLINE',gfound)
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
!
! Check validity of NAM_IO_OFFLINE settings
 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    ')  
!
IF (ctimeseries_filetype=='NETCDF') ctimeseries_filetype='OFFLIN'
!
! Setting input files read from namelist
IF ( csurf_filetype == "LFI   " ) THEN
#ifdef SFX_LFI
  cfilein_lfi      = cprepfile
  cfile_lfi        = cprepfile
  cfilein_lfi_save = cprepfile
  cfilepgd_lfi     = cpgdfile
#endif
ELSEIF ( csurf_filetype == "FA    " ) THEN
#ifdef SFX_FA
  cfilein_fa      = adjustl(adjustr(cprepfile)//'.fa')
  cfilein_fa_save = adjustl(adjustr(cprepfile)//'.fa')
  cfilepgd_fa     = adjustl(adjustr(cpgdfile)//'.fa')
#endif
ELSEIF ( csurf_filetype == "ASCII " ) THEN
#ifdef SFX_ASC
  cfilein      = adjustl(adjustr(cprepfile)//'.txt')
  cfilein_save = adjustl(adjustr(cprepfile)//'.txt')
  cfilepgd     = adjustl(adjustr(cpgdfile)//'.txt')
#endif
ELSEIF ( csurf_filetype == "NC    " ) THEN
#ifdef SFX_NC
  cfilein_nc      = adjustl(adjustr(cprepfile)//'.nc')
  cfilein_nc_save = adjustl(adjustr(cprepfile)//'.nc')
  cfilepgd_nc     = adjustl(adjustr(cpgdfile)//'.nc')
#endif
ELSE
  CALL abor1_sfx(trim(csurf_filetype)//" is not implemented!")
ENDIF
!
! Allocation 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.2.    Goto model of Surfex Types 
!
icurrent_model = 1
!
 cprogname = csurf_filetype
lread_all = .true.
!
! Initialization 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
!
 CALL init_index_mpi(ysurf_cur, &
                    csurf_filetype, yalg_mpi, xio_frac, .false.)
!
!
! Initialize time information
iyear    = nundef
imonth   = nundef
iday     = nundef
ztime    = xundef
 CALL set_surfex_filein(csurf_filetype,'PREP')
 CALL init_io_surf_n(ysurf_cur%DTCO, ysurf_cur%DGU, ysurf_cur%U, &
                        csurf_filetype,'FULL  ','SURF  ','READ ') 
 CALL read_surf(&
                csurf_filetype,'DIM_FULL  ',idim_full,  iresp)
 CALL read_surf(&
                csurf_filetype,'DTCUR     ',ttime,  iresp)
 CALL end_io_surf_n(csurf_filetype)
!
 CALL get_size_full_n(ysurf_cur%U,csurf_filetype,idim_full,isize_full)
!
IF (ALLOCATED(nmask_full)) DEALLOCATE(nmask_full)
!
isw = 0
isv = 0
ALLOCATE(csv(isv))
ALLOCATE(xco2(isize_full))
ALLOCATE(xrhoa(isize_full))
ALLOCATE(xzenith(isize_full))
ALLOCATE(xazim(isize_full))
ALLOCATE(xsw_bands(isw))
ALLOCATE(xdir_alb(isize_full,isw))
ALLOCATE(xsca_alb(isize_full,isw))
ALLOCATE(xemis(isize_full))
ALLOCATE(xtsrad(isize_full))
ALLOCATE(xtsurf(isize_full))
!
! Indicate that zenith and azimuth angles are not initialized
xzenith = xundef
xazim   = xundef
xco2    = 0.
xrhoa   = 1.
!
! Sanity check
IF ( .NOT. lassim ) CALL abor1_sfx("YOU CAN'T RUN SODA WITHOUT SETTING LASSIM=.TRUE. IN THE ASSIM NAMELIST")
!
! Set the number of initializations to be done
! Default is one
inb = 1
IF ( cassim_isba == 'EKF  ' ) THEN
  ! Has to do initialization for all the perturbations + 
  ! control + the real run at last
  inb = nvar + 2
  isize = nvar
ELSEIF ( cassim_isba == 'ENKF ' ) THEN
  inb = nens
  IF (lbias_correction) inb = inb + 1
  isize = nens
ENDIF
!
IF (nrank==npio) WRITE(*,*) "INITIALIZING SURFEX..."
!
yinit = 'ALL'
!
iyear  = ttime%TDATE%YEAR
imonth = ttime%TDATE%MONTH
iday   = ttime%TDATE%DAY
ztime  = ttime%TIME
ihour  = nint(ztime/3600.)
!
nobs = 0
!
lread_all = .false.
!
DO nific = inb,1,-1
  !
  ! If we have more than one initialization to do
  ! For last initialization, we must re-do the first.
  ! Could be avoided by introducing knowlegde of LASSIM on this level
  IF ( cassim_isba == 'EKF  ' .OR. cassim_isba == 'ENKF ' ) THEN
    !
    IF (cassim_isba == 'EKF  ') THEN
      IF ( nific<inb ) THEN
        ymfile = "PREP_"
        CALL get_file_name(iyear,imonth,iday,ihour,ymfile)              
        WRITE(yvar,'(I1.1)') nific-1
        yfilein = trim(ymfile)//"_EKF_PERT"//adjustl(yvar)
      ELSE
        yfilein = "PREP_INIT"
      ENDIF
    ELSEIF (cassim_isba == 'ENKF ') THEN
      ymfile = "PREP_"
      CALL get_file_name(iyear,imonth,iday,ihour,ymfile)              
      WRITE(yvar,'(I1.1)') nific
      yfilein = trim(ymfile)//"_EKF_ENS"//adjustl(yvar)
    ENDIF
    !
    CALL set_filein(yfilein)
    !
  ENDIF
  !
  CALL dealloc_surf_atm_n(ysurf_cur)
  !
  IF ( (cassim_isba=='EKF  '.OR.cassim_isba=='ENKF ') .AND. nific==1 ) lread_all = .true.
  !    
  ! Initialize the SURFEX interface
  CALL io_buff_clean
  CALL init_surf_atm_n(ysurf_cur, &
                       csurf_filetype,yinit, lland_use, isize_full, isv, isw,  &
                       csv, xco2, xrhoa, xzenith, xazim, xsw_bands,     &
                       xdir_alb, xsca_alb, xemis, xtsrad, xtsurf,       &
                       iyear, imonth, iday, ztime,                      &
                       yatmfile, yatmfiletype, ytest              )
  !
  IF ( cassim_isba=='EKF  ' .OR. cassim_isba=='ENKF ' ) THEN
    !
    IF ( nific==inb ) THEN
      ALLOCATE(xlai_pass(ysurf_cur%U%NSIZE_NATURE,ysurf_cur%IM%I%NPATCH)) 
      ALLOCATE(xbio_pass(ysurf_cur%U%NSIZE_NATURE,ysurf_cur%IM%I%NPATCH))     
      IF (cassim_isba=='EKF  ') ALLOCATE(xi(ysurf_cur%U%NSIZE_NATURE,ysurf_cur%IM%I%NPATCH,isize))
      ALLOCATE(xf(ysurf_cur%U%NSIZE_NATURE,ysurf_cur%IM%I%NPATCH,isize+1,nvar))
      ALLOCATE(xf_patch(ysurf_cur%U%NSIZE_NATURE,ysurf_cur%IM%I%NPATCH,isize+1,nobstype))
    ENDIF
    !
    IF ( cassim_isba=='EKF  ' .AND. nific<inb .OR. cassim_isba=='ENKF ') THEN
      !
      ! Set the global state values for this control value
      DO iobs = 1,nobstype
        DO ji=1,ysurf_cur%IM%I%NPATCH
          SELECT CASE (trim(cobs(iobs)))
            CASE("T2M")
              xf_patch(:,ji,nific,iobs) = xat2m_isba(:,1)
            CASE("HU2M")
              xf_patch(:,ji,nific,iobs) = xahu2m_isba(:,1)
            CASE("WG1")
              xf_patch(:,ji,nific,iobs) = ysurf_cur%IM%I%XWG(:,1,ji)
            CASE("LAI")
              xf_patch(:,ji,nific,iobs) = ysurf_cur%IM%I%XLAI(:,ji)
            CASE default
              CALL abor1_sfx("Mapping of "//cobs(iobs)//" is not defined in SODA!")
          END SELECT
        ENDDO
      ENDDO
      !
      ! Prognostic fields for assimilation (Control vector)
      DO jl = 1,nvar
        SELECT CASE (trim(cvar(jl)))
          CASE("TG1")
            xf(:,:,nific,jl) = ysurf_cur%IM%I%XTG(:,1,:)
          CASE("TG2")
            xf(:,:,nific,jl) = ysurf_cur%IM%I%XTG(:,2,:)
          CASE("WG1")
            xf(:,:,nific,jl) = ysurf_cur%IM%I%XWG(:,1,:)
          CASE("WG2")
            xf(:,:,nific,jl) = ysurf_cur%IM%I%XWG(:,2,:)
          CASE("WG3")
            xf(:,:,nific,jl) = ysurf_cur%IM%I%XWG(:,3,:)              
          CASE("LAI")
            xf(:,:,nific,jl) = ysurf_cur%IM%I%XLAI(:,:)
          CASE default
            CALL abor1_sfx("Mapping of "//trim(cvar(jl))//" is not defined in SODA!")
        END SELECT
      ENDDO
      !
      IF ( nific==1 ) THEN
        !
        DO jl = 1,nvar
          IF (trim(cvar(jl))=="LAI") THEN
            IF ( ysurf_cur%IM%I%NPATCH==1 .AND. trim(cbio)/="LAI" ) THEN
              CALL abor1_sfx("Mapping of "//cbio//" is not defined in EKF with NPATCH=1!")
            ENDIF
            SELECT CASE (trim(cbio))
              CASE("BIOMA1","BIOMASS1")
                xbio_pass(:,:) = ysurf_cur%IM%I%XBIOMASS(:,1,:)
              CASE("BIOMA2","BIOMASS2")
                xbio_pass(:,:) = ysurf_cur%IM%I%XBIOMASS(:,2,:)
              CASE("RESPI1","RESP_BIOM1")
                xbio_pass(:,:) = ysurf_cur%IM%I%XRESP_BIOMASS(:,1,:)
              CASE("RESPI2","RESP_BIOM2")
                xbio_pass(:,:) = ysurf_cur%IM%I%XRESP_BIOMASS(:,2,:)
              CASE("LAI")
                xbio_pass(:,:) = ysurf_cur%IM%I%XLAI(:,:)
              CASE default
                CALL abor1_sfx("Mapping of "//cbio//" is not defined in EKF!")
            END SELECT
            !
            xlai_pass(:,:) = ysurf_cur%IM%I%XLAI(:,:)          
            !
          ENDIF
          !
        ENDDO
      ENDIF
      !
    ELSE
      !
      DO jl = 1,nvar
        SELECT CASE (trim(cvar(jl)))
          CASE("TG1")
            xi(:,:,jl) = ysurf_cur%IM%I%XTG(:,1,:)
          CASE("TG2")
            xi(:,:,jl) = ysurf_cur%IM%I%XTG(:,2,:)
          CASE("WG1")
            xi(:,:,jl) = ysurf_cur%IM%I%XWG(:,1,:)
          CASE("WG2")
            xi(:,:,jl) = ysurf_cur%IM%I%XWG(:,2,:)
          CASE("WG3")
            xi(:,:,jl) = ysurf_cur%IM%I%XWG(:,3,:)               
          CASE("LAI")
            xi(:,:,jl) = ysurf_cur%IM%I%XLAI(:,:)
          CASE default
            CALL abor1_sfx("Mapping of "//trim(cvar(jl))//" is not defined in SODA!")
        END SELECT
      ENDDO        
      !
    ENDIF
    !
  ENDIF
  !
ENDDO
!
! Allocate input fields to the assimilation interface
ALLOCATE(zlsm(isize_full))
ALLOCATE(zcon_rain(isize_full))
ALLOCATE(zstrat_rain(isize_full))
ALLOCATE(zcon_snow(isize_full))
ALLOCATE(zstrat_snow(isize_full))
ALLOCATE(zclouds(isize_full))
ALLOCATE(zevaptr(isize_full))
ALLOCATE(zevap(isize_full))
ALLOCATE(ztsc(isize_full))
ALLOCATE(zswec(isize_full))
ALLOCATE(zts(isize_full))
ALLOCATE(zucls(isize_full))
ALLOCATE(zvcls(isize_full))
ALLOCATE(zsst(isize_full))
ALLOCATE(zsic(isize_full))
zts(:) = xundef

! Allocate observations
ALLOCATE(zt2m(isize_full))
ALLOCATE(zhu2m(isize_full))
ALLOCATE(zswe(isize_full))
!
! OI needs first guess values used in oi_cacsts
IF (cassim_isba=="OI   ") THEN
  !
  IF ( trim(cfile_format_fg) == "ASCII" ) THEN

    ALLOCATE(zwork(ysurf_cur%U%NDIM_FULL,8))

    IF (nrank==npio) THEN

      ymfile = 'FIRST_GUESS_'
      CALL get_file_name(iyear,imonth,iday,ihour,ymfile)
      WRITE(*,*) "READING first guess from file "//trim(ymfile)//".DAT"
      istat = 0
      OPEN(unit=55,file=trim(ymfile)//".DAT",form='FORMATTED',status='OLD',iostat=istat)
      IF ( istat /= 0 ) CALL abor1_sfx("Can not open "//trim(ymfile))

      zwork(:,:) = xundef

      ! Read first guess
      DO ji = 1,ysurf_cur%U%NDIM_FULL
        READ (55,*,iostat=istat)  (zwork(ji,jj),jj=1,8)
        IF ( istat /= 0 ) CALL abor1_sfx("Error reading file "//trim(ymfile))
      ENDDO
      CLOSE(55)
    ENDIF

    ! Distribute values on processors
    IF (nproc>1) THEN
      CALL read_and_send_mpi(zwork(:,1),zcon_rain(:))
      CALL read_and_send_mpi(zwork(:,2),zstrat_snow(:))
      CALL read_and_send_mpi(zwork(:,3),zcon_snow(:))
      CALL read_and_send_mpi(zwork(:,4),zstrat_snow(:))
      CALL read_and_send_mpi(zwork(:,5),zclouds(:))
      CALL read_and_send_mpi(zwork(:,6),zlsm(:))
      CALL read_and_send_mpi(zwork(:,7),zevap(:))
      CALL read_and_send_mpi(zwork(:,8),zevaptr(:))
    ELSE
      ! Set First-Guess variables
      zcon_rain(:) = zwork(:,1)
      zstrat_snow(:) = zwork(:,2)
      zcon_snow(:) = zwork(:,3)
      zstrat_snow(:) = zwork(:,4)
      zclouds(:) = zwork(:,5)
      zlsm(:) = zwork(:,6)
      zevap(:) = zwork(:,7)
      zevaptr(:) = zwork(:,8)
    ENDIF

    DEALLOCATE(zwork)

  ELSEIF ( trim(cfile_format_fg) == "FA" ) THEN
    !
    !  Read atmospheric forecast fields from FA files 
#ifdef SFX_FA
    cfilein_fa = 'FG_OI_MAIN'
    cdnomc     = 'oimain'                  ! new frame name

    !  Open FA file (LAM version with extension zone)
    CALL init_io_surf_n(ysurf_cur%DTCO, ysurf_cur%DGU, ysurf_cur%U, &
                        yprogram2,'EXTZON','SURF  ','READ ') 
    !
    !  Read model forecast quantities
    IF (larome) THEN  
      CALL read_surf(yprogram2,'SURFACCPLUIE',  zstrat_rain,iresp)
      CALL read_surf(yprogram2,'SURFACCNEIGE',  zstrat_snow,iresp)
      CALL read_surf(yprogram2,'SURFACCGRAUPEL',zcon_snow,iresp)
      ! So far graupel has not been used
      !ZCON_SNOW=ZCON_SNOW+ZCON_GRAUPEL
      zcon_rain(:) = 0.0
    ELSE    
      CALL read_surf(yprogram2,'SURFPREC.EAU.CON',zcon_rain    ,iresp)
      CALL read_surf(yprogram2,'SURFPREC.EAU.GEC',zstrat_rain  ,iresp)
      CALL read_surf(yprogram2,'SURFPREC.NEI.CON',zcon_snow    ,iresp)
      CALL read_surf(yprogram2,'SURFPREC.NEI.GEC',zstrat_snow  ,iresp)
    ENDIF
    !
    CALL read_surf(yprogram2,'ATMONEBUL.BASSE ',zclouds,iresp)
    CALL read_surf(yprogram2,'SURFIND.TERREMER',zlsm   ,iresp)
    CALL read_surf(yprogram2,'SURFFLU.LAT.MEVA',zevap  ,iresp) ! accumulated fluxes (not available in LFI)
    !
    IF (.NOT.laladsurf) THEN    
      CALL read_surf(yprogram2,'SURFXEVAPOTRANSP',zevaptr,iresp) ! not in ALADIN SURFEX
    ELSE
      zevaptr(:) = 0.0
    ENDIF
    !
    !  Close FA file
    CALL end_io_surf_n(yprogram2)
    CALL io_buff_clean
#else
    CALL abor1_sfx("The first guess is supposed to be an FA file. You must compile with FA support enabled: -DSFX_FA")
#endif
  ELSE
    CALL abor1_sfx("CFILE_FORMAT_FG="//trim(cfile_format_fg)//" not implemented!")
  ENDIF
  IF (nrank==npio .AND. nprintlev>0) WRITE(*,*)'READ FIRST GUESS OK'
ENDIF
!
! If we want to extrapolate values, we need to have a land-sea-mask available
IF (lextrap_sea .OR. lextrap_water .OR. lextrap_nature .OR. .NOT.lwatertg2) THEN

  IF (nrank==npio .AND. nprintlev>0) WRITE(*,*) "READING Land-Sea mask"

  IF ( trim(cfile_format_lsm) == "ASCII" ) THEN

    ALLOCATE(zwork(ysurf_cur%U%NDIM_FULL,1))

    IF (nrank==npio) THEN
      ymfile = 'LSM.DAT'
      IF (nprintlev>0) WRITE(*,*) "READING LSM from file "//trim(ymfile)
      OPEN(unit=55,file=trim(ymfile),form='FORMATTED',status='OLD',iostat=istat)
      IF ( istat /= 0 ) CALL abor1_sfx("Can not open "//trim(ymfile))

      zwork(:,:) = xundef
      ! Read LSM
      DO ji = 1,ysurf_cur%U%NDIM_FULL
        READ (55,*,iostat=istat)  zwork(ji,1)
        IF ( istat /= 0 ) CALL abor1_sfx("Error reading file "//trim(ymfile))
      ENDDO
      CLOSE(55)
    ENDIF

    ! Distribute values on processors
    IF (nproc>1) THEN
      CALL read_and_send_mpi(zwork(:,1),zlsm(:))
    ELSE
      ! Set First-Guess variables
      zlsm(:)=zwork(:,1)
    ENDIF
    DEALLOCATE(zwork)

  ELSEIF ( trim(cfile_format_lsm) == "FA" ) THEN
    !  Read atmospheric forecast fields from FA files 
#ifdef SFX_FA
    cfilein_fa = 'FG_OI_MAIN'
    cdnomc     = 'oimain'                  ! new frame name

    !  Open FA file (LAM version with extension zone)
    CALL init_io_surf_n(ysurf_cur%DTCO, ysurf_cur%DGU, ysurf_cur%U, &
                                yprogram2,'EXTZON','SURF  ','READ ')
    CALL read_surf(yprogram2,'SURFIND.TERREMER',zlsm   ,iresp)
    IF (iresp /=0) CALL abor1_sfx("Can not read Land-Sea mask from FA file "//trim(cfilein_fa))
    !  Close FA file
    CALL end_io_surf_n(yprogram2)
    CALL io_buff_clean
#else
    CALL abor1_sfx("The Land-sea mask file is assumed to be an FA file. You must compile with FA support enabled: -DSFX_FA")
#endif
  ELSE
    CALL abor1_sfx("CFILE_FORMAT_LSM="//trim(cfile_format_lsm)//" not implemented!")
  ENDIF
  IF (nrank==npio.AND.nprintlev>0) WRITE(*,*)'READ LSM OK'
ENDIF
!
! The observations used in the analysis is read.
! The options are either from a CANARI analysis or from a file
IF ( cassim_isba=="EKF  " .OR. cassim_isba=="ENKF " ) THEN
  ALLOCATE(xyo(ysurf_cur%U%NSIZE_NATURE,nobstype))
ENDIF

IF ( trim(cfile_format_obs) == "ASCII") THEN

  ALLOCATE(zwork(ysurf_cur%U%NDIM_FULL,nobstype))

  IF (nrank==npio) THEN

    ymfile = 'OBSERVATIONS_'
    CALL get_file_name(iyear,imonth,iday,ihour,ymfile)
    ymfile=trim(ymfile)//".DAT"
    istat = 0
    OPEN(unit=55,file=trim(ymfile),form='FORMATTED',status='OLD',iostat=istat)
    IF ( istat /=0 ) CALL abor1_sfx("Error opening file "//trim(ymfile))

    ! Get the nature points from processors
    ! If the file has an header, we check for consistency
    IF ( lobsheader ) THEN
      ! Read in first line and check if variables are consistent
      ALLOCATE(cobsinfile(nobstype))
      READ (55,*,iostat=istat)  (cobsinfile(jj),jj=1,nobstype)
      IF ( istat /= 0 ) CALL abor1_sfx("Error reading header in "//trim(ymfile))

      DO iobs = 1,nobstype
        IF ( trim(cobs(iobs)) /= trim(cobsinfile(iobs))) THEN
          CALL abor1_sfx("Mapping of observations in "//trim(ymfile)//&
               " is not consistent with setup! "//trim(cobs(iobs))//" /= "//trim(cobsinfile(iobs)))
        ENDIF
      ENDDO
      DEALLOCATE(cobsinfile)
    ENDIF

    zwork(:,:) = xundef

  ENDIF

  !   Read all observations (NDIM_FULL)
  IF (lobsnat) THEN

    IF (nrank==npio) THEN
      ALLOCATE(zyo_nat(ysurf_cur%U%NDIM_NATURE,nobstype))
      DO ji = 1,ysurf_cur%U%NDIM_NATURE
        READ (55,*,iostat=istat)  (zyo_nat(ji,jj),jj=1,nobstype)
        IF ( istat /= 0 ) CALL abor1_sfx("Error reading file "//trim(ymfile))
      ENDDO
      ALLOCATE(znature(ysurf_cur%U%NDIM_FULL))
    ENDIF

    IF (nproc>1) THEN    
      CALL gather_and_write_mpi(ysurf_cur%U%XNATURE,znature)
    ELSEIF (nrank==npio) THEN
      znature(:) = ysurf_cur%U%XNATURE
    ENDIF

    IF (nrank==npio) THEN
      icpt = 0
      DO ji = 1,ysurf_cur%U%NDIM_FULL
        IF (znature(ji)>0.) THEN
          icpt = icpt + 1
          zwork(ji,:) = zyo_nat(icpt,:)
        ENDIF
      ENDDO
      DEALLOCATE(znature,zyo_nat) 
    ENDIF

  ELSEIF (nrank==npio) THEN

    DO ji = 1,ysurf_cur%U%NDIM_FULL
      READ (55,*,iostat=istat)  (zwork(ji,jj),jj=1,nobstype)
      IF ( istat /= 0 ) CALL abor1_sfx("Error reading file "//trim(ymfile))
    ENDDO
  ENDIF

  IF (nrank==npio) THEN

    CLOSE(55)
    IF (nprintlev>0) WRITE(*,*) 'Read observation file OK'

  ENDIF

  ! Initialize possible observations
  zt2m  = 999.
  zhu2m = 999.
  zswe  = 999.

  IF (nproc>1) THEN

    ! Running on more CPU's
    ! For EKF/EnKF we must distribute variables for nature tile
    IF ( cassim_isba=="EKF  " .OR. cassim_isba=="ENKF " ) THEN
      DO jj=1,nobstype
        CALL read_and_send_mpi(zwork(:,jj),xyo(:,jj),ysurf_cur%U%NR_NATURE)
      ENDDO
    ENDIF

    ! Set observations used for possibly other tiles than nature
    ! Distribute read variables
    jj = 1
    DO ji = 1,nobsmax
      IF (nnco(ji) == 1 .AND. jj <= nobstype ) THEN
        SELECT CASE (ji)
          CASE (1)
            CALL read_and_send_mpi(zwork(:,jj),zt2m(:))
          CASE (2)
            CALL read_and_send_mpi(zwork(:,jj),zhu2m(:))
          CASE (5)
            CALL read_and_send_mpi(zwork(:,jj),zswe(:))
        END SELECT
        jj = jj + 1
      ENDIF
    ENDDO

  ELSE

    ! Running on one CPU
    IF ( cassim_isba=="EKF  " .OR. cassim_isba=="ENKF " ) THEN
      DO ji = 1,ysurf_cur%U%NSIZE_NATURE
        xyo(ji,:) = zwork(ysurf_cur%U%NR_NATURE(ji),:)
      ENDDO
    ENDIF
    ! Set observations used for possibly other tiles than nature
    jj = 1
    DO ji = 1,nobsmax
      IF (nnco(ji) == 1 .AND. jj <= nobstype ) THEN 
        SELECT CASE (ji)
          CASE (1)
            zt2m(:)=zwork(:,jj)
          CASE (2)
            zhu2m(:)=zwork(:,jj)
          CASE (5)
            zswe(:)=zwork(:,jj)
         END SELECT
         jj = jj + 1
       ENDIF
     ENDDO
  ENDIF
  DEALLOCATE(zwork)

  nobs = nobs + nobstype
  IF (( cassim_isba=="EKF  " .OR. cassim_isba=="ENKF " ) .AND. ( nprintlev > 2 )) WRITE(iluout,*) 'read in obs: ', xyo(1,:), nobs

ELSEIF ( trim(cfile_format_obs) == "FA") THEN
  !      
  nobs = nobstype
  !
  !
  !  Define FA file name for CANARI analysis
#ifdef SFX_FA
  cfilein_fa = 'CANARI'        ! input CANARI analysis
  cdnomc     = 'canari'                  ! new frame name

!  Open FA file 
  CALL init_io_surf_n(ysurf_cur%DTCO, ysurf_cur%DGU, ysurf_cur%U, &
                        yprogram2,'EXTZON','SURF  ','READ ')
!
!  Read CANARI analysis
  CALL read_surf(yprogram2,'CLSTEMPERATURE  ',zt2m ,iresp)
  CALL read_surf(yprogram2,'CLSHUMI.RELATIVE',zhu2m,iresp)
  CALL read_surf(yprogram2,'SURFTEMPERATURE ',zts  ,iresp)
  CALL read_surf(yprogram2,'SURFRESERV.NEIGE',zswe ,iresp)
  CALL read_surf(yprogram2,'CLSVENT.ZONAL   ',zucls,iresp)
  CALL read_surf(yprogram2,'CLSVENT.MERIDIEN',zvcls,iresp)  

!  Close CANARI file
  CALL end_io_surf_n(yprogram2)
  CALL io_buff_clean
  IF (nrank==npio.AND.nprintlev>0) WRITE(*,*) 'READ CANARI OK'
#else
  CALL abor1_sfx("CANARI analyis is supposed to be an FA file. You must compile with FA support enabled: -DSFX_FA")       
#endif
ELSE
  CALL abor1_sfx("CFILE_FORMAT_OBS="//trim(cfile_format_obs)//" not implemented!")  
ENDIF
!
! Climatological fields are only used in OI
IF (cassim_isba=="OI   ") THEN
  !
  IF (trim(cfile_format_clim) == "ASCII" ) THEN

    ALLOCATE(zwork(ysurf_cur%U%NDIM_FULL,2))

    IF (nrank==npio) THEN
      ymfile = 'CLIMATE.DAT'
      WRITE(*,*) "READING CLIM from file "//trim(ymfile)
      OPEN(unit=55,file=trim(ymfile),form='FORMATTED',status='OLD',iostat=istat)
      IF ( istat /= 0 ) CALL abor1_sfx("Can not open "//trim(ymfile))

      zwork(:,:) = xundef
      ! Read CLIMATE file
      DO ji = 1,ysurf_cur%U%NDIM_FULL
        READ (55,*,iostat=istat) (zwork(ji,jj),jj=1,2)
        IF ( istat /= 0 ) CALL abor1_sfx("Error reading file "//trim(ymfile))
      ENDDO
      CLOSE(55)
    ENDIF

    ! Distribute values on processors
    IF (nproc>1) THEN
      CALL read_and_send_mpi(zwork(:,1),zswec(:))
      CALL read_and_send_mpi(zwork(:,2),ztsc(:))
    ELSE
      ! Set First-Guess variables
      zswec(:)=zwork(:,1)
      ztsc=zwork(:,2)
    ENDIF
    DEALLOCATE(zwork)

  ELSEIF (trim(cfile_format_clim) == "FA" ) THEN  

  !  Define FA file name for surface climatology
#ifdef SFX_FA
  cfilein_fa = 'clim_isba'               ! input climatology
  cdnomc     = 'climat'                  ! new frame name

!  Open FA file 
  CALL init_io_surf_n(ysurf_cur%DTCO, ysurf_cur%DGU, ysurf_cur%U, &
                        yprogram2,'EXTZON','SURF  ','READ ')
!
!  Read climatology file
  CALL read_surf(yprogram2,'SURFRESERV.NEIGE',zswec,iresp)
  CALL read_surf(yprogram2,'SURFTEMPERATURE' ,ztsc ,iresp)

!  Close climatology file
  CALL end_io_surf_n(yprogram2)
  CALL io_buff_clean
#else
    CALL abor1_sfx("The climate file is supposed to be an FA file. You must compile with FA support enabled: -DSFX_FA")
#endif
  ELSE
    CALL abor1_sfx("CFILE_FORMAT_CLIM="//trim(cfile_format_clim)//" not implemented!")
  ENDIF  
  IF (nrank==npio.AND.nprintlev>0) WRITE(*,*) 'READ CLIMATOLOGY OK'
  !
ENDIF
!
 CALL assim_set_sst(ysurf_cur%DTCO, ysurf_cur%DGU, ysurf_cur%SM%S, ysurf_cur%U, &
                    isize_full,zlsm,zsst,zsic,ytest)

IF ( .NOT. lassim ) CALL abor1_sfx("YOU CAN'T RUN SODA WITHOUT SETTING LASSIM=.TRUE. IN THE ASSIM NAMELIST")
!
ALLOCATE(gd_maskext(isize_full))
gd_maskext(:) = .false.
!
ALLOCATE(zlon(isize_full))
ALLOCATE(zlat(isize_full))
zlon(:) = ysurf_cur%UG%XLON(:)
zlat(:) = ysurf_cur%UG%XLAT(:)        
!
glkeepextzone = .true.
!
IF (nrank==npio) WRITE(*,*) 'PERFORMIMG OFFLINE SURFEX DATA ASSIMILATION...'
 CALL assim_surf_atm_n(ysurf_cur%IM%DGMI, ysurf_cur%IM%IG, ysurf_cur%IM%I, ysurf_cur%SM%S, &
                       ysurf_cur%U, ysurf_cur%TM%T, ysurf_cur%TM%TOP, ysurf_cur%WM%W, &
                       csurf_filetype,isize_full,   &
                      zcon_rain,  zstrat_rain, zcon_snow, zstrat_snow, &
                      zclouds,    zlsm,        zevaptr,   zevap,       &
                      zswec,      ztsc,       &
                      zts,        zt2m,        zhu2m,     zswe,        &
                      zsst,       zsic,  zucls,     zvcls,       &
                      ytest, gd_maskext, zlon, zlat, glkeepextzone )
!
DEALLOCATE(zcon_rain)
DEALLOCATE(zstrat_rain)
DEALLOCATE(zcon_snow)
DEALLOCATE(zstrat_snow)
DEALLOCATE(zclouds)
DEALLOCATE(zlsm)
DEALLOCATE(zevaptr)
DEALLOCATE(zevap)
DEALLOCATE(ztsc)
DEALLOCATE(zswec)
DEALLOCATE(zts)
DEALLOCATE(zt2m)
DEALLOCATE(zhu2m)
DEALLOCATE(zswe)
DEALLOCATE(zucls)
DEALLOCATE(zvcls)
DEALLOCATE(zsst)
DEALLOCATE(zsic)
!
!
IF (ctimeseries_filetype=="OFFLIN") THEN
  CALL init_output_ol_n(ysurf_cur)
ENDIF
!
ztime_out  = ztime
iday_out   = iday
imonth_out = imonth
iyear_out  = iyear
!
IF (nrank==npio) THEN
  !
  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')
#ifdef SFX_LFI
  cfileout_lfi= adjustl(adjustr(csurffile)//'.'//ytag)
#endif
#ifdef SFX_FA  
  cfileout_fa = adjustl(adjustr(csurffile)//'.'//ytag//'.fa')
#endif
#ifdef SFX_NC
  cfileout_nc = adjustl(adjustr(csurffile)//'.'//ytag//'.nc')
#endif
  !
  IF (csurf_filetype=='FA    ') THEN
#ifdef SFX_FA
    cdnomc = 'header'
    lfanocompact = ldiag_fa_nocompact
    idatef(1)= iyear_out
    idatef(2)= imonth_out
    idatef(3)= iday_out
    idatef(4)= floor(ztime_out/3600.)
    idatef(5)= floor(ztime_out/60.) - idatef(4) * 60 
    idatef(6)= nint(ztime_out) - 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
!
isize = 1
IF (cassim_isba=="ENKF ") THEN
  isize = nens
  IF (lbias_correction) isize = isize + 1
ENDIF
!
nwrite = 1
xstartw = 1
ltime_written(:)=.false.
!
DO iens = 1,isize
  !
  IF (cassim_isba=="ENKF ") THEN
    !
    ymfile = "PREP_"
    CALL get_file_name(iyear,imonth,iday,ihour,ymfile)              
    WRITE(yvar,'(I3)') iens
    yfilein = trim(ymfile)//"_EKF_ENS"//adjustl(yvar)
    CALL set_filein(yfilein)
    !
    lread_all = .true.
    !
    CALL dealloc_surf_atm_n(ysurf_cur)    
    !
    ! Initialize the SURFEX interface
    CALL io_buff_clean
    CALL init_surf_atm_n(ysurf_cur, &
                            csurf_filetype,yinit, lland_use, isize_full, isv, isw,  &
                         csv, xco2, xrhoa, xzenith, xazim, xsw_bands,     &
                         xdir_alb, xsca_alb, xemis, xtsrad, xtsurf,       &
                         iyear, imonth, iday, ztime,                      &
                          yatmfile, yatmfiletype, ytest              )
    !
    !
    DO jl=1,nvar
      !
      ! Update the modified values
      SELECT CASE (trim(cvar(jl)))
        CASE("TG1")
          ysurf_cur%IM%I%XTG(:,1,:) = xf(:,:,iens,jl)
        CASE("TG2")
          ysurf_cur%IM%I%XTG(:,2,:) = xf(:,:,iens,jl)
        CASE("WG1")
          ysurf_cur%IM%I%XWG(:,1,:) = xf(:,:,iens,jl)
        CASE("WG2")
          ysurf_cur%IM%I%XWG(:,2,:) = xf(:,:,iens,jl)
        CASE("WG3")
          ysurf_cur%IM%I%XWG(:,3,:) = xf(:,:,iens,jl)          
        CASE("LAI") 
          ysurf_cur%IM%I%XLAI(:,:) = xf(:,:,iens,jl)
        CASE default
          CALL abor1_sfx("Mapping of "//trim(cvar(jl))//" is not defined in EKF!")
      END SELECT
    ENDDO
    !
  ENDIF
  !
#ifdef SFX_NC  
  ldef = .true.
#endif
  !
  IF (ctimeseries_filetype=="NC    ") THEN
    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
  !
  inw = 1
  IF (ctimeseries_filetype=="NC    ") inw = 2
  !
  DO jnw = 1,inw
    CALL io_buff_clean
    CALL write_surf_atm_n(ysurf_cur, &
                        ctimeseries_filetype,'ALL',lland_use)
    CALL write_diag_surf_atm_n(ysurf_cur, &
                             ctimeseries_filetype,'ALL')
#ifdef SFX_NC
      ldef = .false.
#endif
  ENDDO
  !
  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.
  gpgd_isba              = .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    )
  ! 
  yens = '   '
  IF (isize>1) WRITE(yens,'(I3)') iens
  !
  IF ( csurf_filetype == "LFI   " ) THEN
#ifdef SFX_LFI
    cfileout_lfi     = trim(trim(csurffile)//adjustl(yens))
#endif
  ELSEIF ( csurf_filetype == "FA    " ) THEN
#ifdef SFX_FA
    cfileout_fa  = adjustl(trim(adjustr(csurffile)//adjustl(yens))//'.fa')
#endif
  ELSEIF ( csurf_filetype == "ASCII " ) THEN
#ifdef SFX_ASC
    cfileout = adjustl(trim(adjustr(csurffile)//adjustl(yens))//'.txt')
    lcreated = .false.
#endif
  ELSEIF ( csurf_filetype == "NC    " ) THEN
#ifdef SFX_NC
    cfileout_nc = adjustl(trim(adjustr(csurffile)//adjustl(yens))//'.nc')
#endif
  ELSE
    CALL abor1_sfx(trim(csurf_filetype)//" is not implemented!")
  ENDIF
  !
#ifdef SFX_NC
  ldef = .true.
#endif
  !  
  IF (csurf_filetype=="NC    ") THEN
    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
  !  
  inw = 1
  IF (csurf_filetype=="NC    ") inw = 2
  !  
  DO jnw = 1,inw
    !
    CALL io_buff_clean
    !  
    ! Store results from assimilation
    CALL write_surf_atm_n(ysurf_cur, &
                        csurf_filetype,'ALL',lland_use)
    IF (ysurf_cur%DGU%LREAD_BUDGETC.AND..NOT.ysurf_cur%IM%DGEI%LRESET_BUDGETC) THEN
      CALL write_diag_surf_atm_n(ysurf_cur, &
                             csurf_filetype,'ALL')
    ENDIF
#ifdef SFX_NC
    ldef = .false.
#endif
    !
  ENDDO  
  !
ENDDO
!
IF (nrank==npio .AND. csurf_filetype=='FA    ') THEN
#ifdef SFX_FA
  CALL fairme(iret,nunit_fa,'UNKNOWN')
#endif
END IF
!
!*    3.     Close parallelized I/O
!            ----------------------
!
IF (ctimeseries_filetype=='OFFLIN') CALL close_fileout_ol
!
IF (nrank==npio) THEN
  WRITE(iluout,*) ' '
  WRITE(iluout,*) '    -----------------------'
  WRITE(iluout,*) '    | SODA ENDS CORRECTLY |'
  WRITE(iluout,*) '    -----------------------'
  !
  WRITE(*,*) ' '
  WRITE(*,*) '    -----------------------'
  WRITE(*,*) '    | SODA ENDS CORRECTLY |'
  WRITE(*,*) '    -----------------------'
  !
ENDIF
!
 CLOSE(iluout)
!
 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('SODA',1,zhook_handle)
!
#ifdef SFX_MPI
 CALL mpi_finalize(infompi)
#endif
!
!-------------------------------------------------------------------------------

#endif
END PROGRAM soda