oi_control.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.
SUBROUTINE oi_control (YSC, &
 & ldinline,                 &
 & p__surftemperature,       &
 & p__surfprec_eau_con,      &
 & p__surfprec_eau_gec,      &
 & p__surfprec_nei_con,      &
 & p__surfprec_nei_gec,      &
 & p__atmonebul_basse,       &
 & p__surfxevapotransp,      &
 & p__surfflu_lat_meva,      &
 & p__surfaccpluie,          &
 & p__surfaccneige,          &
 & p__surfaccgraupel,        &
 & p__clstemperature,        &
 & p__clshumi_relative,      &
 & p__clsvent_zonal,         &
 & p__clsvent_meridien,      &
 & p__surfind_terremer,      &
 & p__surfreserv_neige,      &
 & p__lon,                   &
 & p__lat,                   &
 & ld_maskext)

! ------------------------------------------------------------------------------------------
!  *****************************************************************************************
!
!  Program to perform within SURFEX 
!  a soil analysis for water content and temperature 
!  using the Meteo-France optimum interpolation technique of Giard and Bazile (2000)
!
!  Derived from CANARI subroutines externalized by Lora Taseva (Dec. 2007)
!
!  Author : Jean-Francois Mahfouf (01/2008)
!
!  Modifications : 
!   (05/2008)  : The I/O of this version follow the newly available LFI format in SURFEX  
!   (01/2009)  : Read directly atmospheric FA files using XRD library instead of using "edf"
!   (06/2009)  : Modifications to allow the assimilation of ASCAT superficial soil moisture
!   (09/2010)  : More parameters to goto_surfex
!   (03/2011)  : Initialization of ZEVAPTR (F.Bouyssel)
!   (03/2013)  : Use 10m wind from upperair instead surfex one (F.Taillefer)
!
! ******************************************************************************************
! ------------------------------------------------------------------------------------------
!
!
!
!
USE modd_surfex_n, ONLY : surfex_t
USE mode_modeln_surfex_handler
!
USE modd_type_date_surf
USE modd_csts,       ONLY : xday, xpi, xrholw, xlvtt, ndaysec
USE modd_surf_par,   ONLY : xundef
USE modd_assim
USE modd_ol_fileid

USE modd_surfex_mpi, ONLY : nrank, npio
USE modd_surfex_omp, ONLY : nindx2sfx, nwork, nwork2, xwork, xwork2, xwork3, &
                            nwork_full, nwork2_full, xwork_full, xwork2_full

USE modn_io_offline, ONLY : nam_io_offline, cnamelist, csurf_filetype


#ifdef SFX_LFI
USE modd_io_surf_lfi,ONLY : cfilein_lfi, cfileout_lfi, cfilepgd_lfi, cfilein_lfi_save
#endif
#ifdef SFX_FA
USE modd_io_surf_fa, ONLY : cfilein_fa, cfilein_fa_save, cdnomc, &
                            ndgux,  ndlux,  perpk,  pelon0, pelat0, &
                            pedelx, pedely, pelon1, pelat1, pebeta
#endif
#ifdef SFX_ARO
USE modd_io_surf_aro,ONLY : ngptot, ngptot_cap, nproma, nindx1, nindx2, nblock, nkproma, &
                                                             ysurfex_cache_out,          &
                            surfex_field_buf_prealloc, surfex_field_buf_set_record
USE modd_surfex_aro, ONLY : ysurfex_aro_all, ysurfex_aro_cur
#endif

USE mode_pos_surf,  ONLY : posnam

USE modi_open_namelist
USE modi_close_namelist
USE modi_read_all_namelists
USE modi_ini_data_cover
USE modi_init_io_surf_n
USE modi_read_surf
USE modi_set_surfex_filein
USE modi_get_size_full_n
USE modi_read_cover_n
USE modi_convert_cover_frac
USE modi_get_1d_mask
USE modi_end_io_surf_n
USE modi_io_buff_clean
USE modi_oi_bc_soil_moisture
USE modi_oi_latlon_conf_proj
USE modi_oi_cacsts
USE modi_oi_hor_extrapol_surf
USE modi_flag_update
USE modi_write_surf

USE yomhook ,ONLY : lhook, dr_hook
USE parkind1  ,ONLY : jprb

IMPLICIT NONE
!
TYPE(surfex_t), INTENT(INOUT) :: ysc
!
LOGICAL, INTENT (IN) :: ldinline
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__surftemperature
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__surfprec_eau_con
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__surfprec_eau_gec
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__surfprec_nei_con
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__surfprec_nei_gec
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__atmonebul_basse
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__surfxevapotransp
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__surfflu_lat_meva
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__surfaccpluie
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__surfaccneige
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__surfaccgraupel
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__clstemperature
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__clshumi_relative
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__clsvent_zonal
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__clsvent_meridien
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__surfind_terremer
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__surfreserv_neige
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__lon
REAL(KIND=JPRB), OPTIONAL, DIMENSION (:) ::  p__lat
LOGICAL,         OPTIONAL, DIMENSION (:) ::  ld_maskext

INTEGER :: igpcomp
INTEGER :: idat

 CHARACTER(LEN=28) :: ynamelist = 'OPTIONS.nam                 '

!    Declarations of local variables

 CHARACTER(LEN=6)            :: yprogram
 CHARACTER(LEN=6), PARAMETER :: yprogram2 = 'FA    '
 CHARACTER(LEN=2)            :: cmonth
INTEGER                     :: iyear                ! current year (UTC)
INTEGER                     :: imonth               ! current month (UTC)
INTEGER                     :: iday                 ! current day (UTC)
INTEGER                     :: nsssss               ! current time since start of the run (s)
INTEGER                     :: iresp                ! return code
TYPE (date_time)            :: ttime                ! Current date and time
INTEGER                     :: isize
INTEGER                     :: isize1
LOGICAL                     :: llkeepextzone

! Arrays for soil OI analysis
REAL, DIMENSION (:,:), ALLOCATABLE :: pws, pwp, pts, ptp, ptl, psns, prsmin, pd2, plai, pveg
REAL, DIMENSION (:),   ALLOCATABLE :: psst, psab, parg, plat, plon, ptcls, phcls, pucls, pvcls,  &
                                     & PEVAP, PEVAPTR, PT2M_O, PHU2M_O, PTS_O, ZT2INC, ZH2INC,   &
                                     & ZWS, ZWP, ZTL, ZTS, ZTP, ZTCLS, ZHCLS, ZUCLS, ZVCLS,      &
                                     & PSSTC, PWPINC1, PWPINC2, PWPINC3, PT2MBIAS, PH2MBIAS,     &
                                     & PRRCN, PRRCL, PRRSN, PRRSL, PATMNEB, PITM, PALBF, PEMISF, &
                                     & PZ0F, PIVEG, PZ0H, PTSC, PTPC, PWSC, PWPC, PSNC, ZEVAP,   &
                                     & ZEVAPTR, PGELAT, PGELAM, PGEMU, ZWSINC, ZWPINC, ZTSINC,   &
                                     & ZTPINC, ZTLINC, ZSNINC, ZSNS, ZPX, ZPY, PSM_O, PSIG_SMO,  &
                                     & PLSM_O, PWS_O, ZWGINC, PLST, PTRD3, ZSST, ZLST, ZALT
REAL, DIMENSION (:),   ALLOCATABLE :: zsst1, zlst1, psst1, plst1, plat1, plon1, zalt1

INTEGER                            :: iversion, ibugfix
INTEGER                            :: jj,j1
 CHARACTER(LEN=10)                  :: yvar    ! Name of the prognostic variable (in LFI file)
 CHARACTER(LEN=100)                 :: yprefix ! Prefix of the prognostic variable  (in LFI file)
INTEGER                            :: iluout  ! ascii output unit number
INTEGER                            :: inobs   ! number of observations
INTEGER :: ilunam
LOGICAL :: gfound

REAL                               :: plat0,plon0,prpk,plator,plonor,delx,dely,pbeta,zthres
REAL(KIND=JPRB)                    :: z1s2pi, zpis180

LOGICAL, DIMENSION(:), ALLOCATABLE :: ointerp_lst, ointerp_sst
LOGICAL, DIMENSION(:), ALLOCATABLE :: ointerp_lst1, ointerp_sst1

REAL(KIND=JPRB) :: zhook_handle
! ----------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('OI_CONTROL', 0, zhook_handle)

print *,'--------------------------------------------------------------------------'
print *,'|                                                                        |'
print *,'|                             ENTER OI_ASSIM                             |'
print *,'|                                                                        |'
print *,'--------------------------------------------------------------------------'

 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)

!
 CALL read_all_namelists(ysc,csurf_filetype,'ALL',.false.)
!
IF (ldinline) THEN

  yprogram = 'AROME'
#ifdef SFX_ARO
  igpcomp = min(ngptot, ngptot_cap)
  
  nblock   = 1
  nindx1   = 1 + (nblock - 1) * nproma
  nindx2   = min(nblock * nproma, igpcomp)
  nkproma  = nindx2 - nindx1 + 1
  ysurfex_aro_cur => ysurfex_aro_all(nblock)
#endif

ELSE

  yprogram = 'LFI'
  
ENDIF

icurrent_model = 1

iluout = 111
llkeepextzone = .false.

z1s2pi=1.0_jprb/(2.0_jprb*xpi)
zpis180=xpi/180.0_jprb

!   Update some constants dependant from NACVEG

!  scaling of soil moisture increments when assimilation window is different from 6 hours
xrscaldw = REAL(nechgu)/6.0_jprb
!  half assimilation window in sec
nitrad   = nechgu*1800

 CALL ini_data_cover(ysc%DTCO, ysc%U)

!   File handling definition

IF (.NOT. ldinline) THEN
#ifdef SFX_LFI
  cfilepgd_lfi = 'PGD'
  cfilein_lfi = 'PREP'        ! input PREP file (surface fields) 
  cfilein_lfi_save = cfilein_lfi
#endif
ENDIF

!   Read grid dimension for allocation

 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        yprogram,'FULL  ','SURF  ','READ ')

!   Find current time

 CALL read_surf(&
               yprogram,'DTCUR',ttime,iresp)

!   Time initializations

iyear  = ttime%TDATE%YEAR
imonth = ttime%TDATE%MONTH
iday   = ttime%TDATE%DAY
nsssss = ttime%TIME
IF (nsssss > ndaysec) nsssss = nsssss - ndaysec

!   Reading grid characteristics to perform nature mask

 CALL end_io_surf_n(yprogram)
 CALL set_surfex_filein(yprogram,'PGD ') ! change input file name to pgd name
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        yprogram,'FULL  ','SURF  ','READ ')

 CALL read_surf(&
               yprogram,'SEA   ',ysc%U%CSEA   ,iresp)
 CALL read_surf(&
               yprogram,'WATER ',ysc%U%CWATER ,iresp)
 CALL read_surf(&
               yprogram,'NATURE',ysc%U%CNATURE,iresp)
 CALL read_surf(&
               yprogram,'TOWN  ',ysc%U%CTOWN  ,iresp)

 CALL read_surf(&
               yprogram,'DIM_FULL  ',ysc%U%NDIM_FULL,  iresp)
nindx2sfx = ysc%U%NDIM_FULL
 CALL end_io_surf_n(yprogram)
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        yprogram,'FULL  ','SURF  ','READ ')

 CALL read_surf(&
               yprogram,'DIM_SEA   ',ysc%U%NDIM_SEA,   iresp)
 CALL read_surf(&
               yprogram,'DIM_NATURE',ysc%U%NDIM_NATURE,iresp)
 CALL read_surf(&
               yprogram,'DIM_WATER ',ysc%U%NDIM_WATER, iresp)
 CALL read_surf(&
               yprogram,'DIM_TOWN  ',ysc%U%NDIM_TOWN,  iresp)

ALLOCATE(nwork(ysc%U%NDIM_FULL))
ALLOCATE(xwork(ysc%U%NDIM_FULL))
ALLOCATE(nwork2(ysc%U%NDIM_FULL,10))
ALLOCATE(xwork2(ysc%U%NDIM_FULL,10))
ALLOCATE(xwork3(ysc%U%NDIM_FULL,10,10))
IF (nrank==npio) THEN
  ALLOCATE(nwork_full(ysc%U%NDIM_FULL))
  ALLOCATE(xwork_full(ysc%U%NDIM_FULL))
  ALLOCATE(nwork2_full(ysc%U%NDIM_FULL,10))
  ALLOCATE(xwork2_full(ysc%U%NDIM_FULL,10))
ELSE
  ALLOCATE(nwork_full(0))
  ALLOCATE(xwork_full(0))
  ALLOCATE(nwork2_full(0,0))
  ALLOCATE(xwork2_full(0,0))
ENDIF
!
!   Get total dimension of domain (excluding extension zone)

 CALL get_size_full_n(ysc%U, &
                     yprogram,ysc%U%NDIM_FULL,ysc%U%NSIZE_FULL)

IF (ldinline) THEN
  isize = ysc%U%NSIZE_FULL
ELSE
  isize = ysc%U%NDIM_FULL
ENDIF

ALLOCATE (psab(isize)) 
ALLOCATE (parg(isize))
ALLOCATE (zalt(isize))

 CALL read_surf(&
               yprogram,'SAND',      psab,  iresp)
 CALL read_surf(&
               yprogram,'CLAY',      parg,  iresp)
 CALL read_surf(&
               yprogram,'ZS',        zalt,  iresp)

 CALL read_cover_n(ysc%DTCO, ysc%U, &
                  yprogram)

!   Perform masks (only nature used)

ALLOCATE(ysc%U%XSEA   (isize))
ALLOCATE(ysc%U%XNATURE(isize))
ALLOCATE(ysc%U%XWATER (isize))
ALLOCATE(ysc%U%XTOWN  (isize))

 CALL convert_cover_frac(ysc%DTCO, &
                        ysc%U%XCOVER,ysc%U%LCOVER,ysc%U%XSEA,ysc%U%XNATURE,ysc%U%XTOWN,ysc%U%XWATER)

ysc%U%NSIZE_NATURE = count(ysc%U%XNATURE(:) > 0.0)
ysc%U%NSIZE_TOWN   = count(ysc%U%XTOWN(:)   > 0.0)
ysc%U%NSIZE_WATER  = count(ysc%U%XWATER(:)  > 0.0)
ysc%U%NSIZE_SEA    = count(ysc%U%XSEA(:)    > 0.0)

ALLOCATE(ysc%U%NR_NATURE (ysc%U%NSIZE_NATURE))
ALLOCATE(ysc%U%NR_TOWN   (ysc%U%NSIZE_TOWN  ))
ALLOCATE(ysc%U%NR_WATER  (ysc%U%NSIZE_WATER ))
ALLOCATE(ysc%U%NR_SEA    (ysc%U%NSIZE_SEA   ))

 CALL get_1d_mask( ysc%U%NSIZE_SEA,    isize, ysc%U%XSEA   , ysc%U%NR_SEA   )
 CALL get_1d_mask( ysc%U%NSIZE_WATER,  isize, ysc%U%XWATER , ysc%U%NR_WATER )
 CALL get_1d_mask( ysc%U%NSIZE_TOWN,   isize, ysc%U%XTOWN  , ysc%U%NR_TOWN  )
 CALL get_1d_mask( ysc%U%NSIZE_NATURE, isize, ysc%U%XNATURE, ysc%U%NR_NATURE)

! Allocate arrays

ysc%IM%I%NPATCH = 1

ALLOCATE (pws(isize,1))
ALLOCATE (pwp(isize,1))
ALLOCATE (pts(isize,1))
ALLOCATE (ptp(isize,1))
ALLOCATE (ptl(isize,1))
ALLOCATE (psst(isize))
ALLOCATE (psns(isize,1))
ALLOCATE (plai(isize,1))
ALLOCATE (pveg(isize,1))
ALLOCATE (prsmin(isize,1))
ALLOCATE (pd2(isize,1))
ALLOCATE (ptcls(isize))
ALLOCATE (phcls(isize))
ALLOCATE (pucls(isize))
ALLOCATE (pvcls(isize))
ALLOCATE (pevap(isize))
ALLOCATE (plst(isize))
ALLOCATE (ptrd3(isize))

ALLOCATE (ointerp_lst(isize))
ALLOCATE (ointerp_sst(isize))
ALLOCATE (zlst(isize))
ALLOCATE (zsst(isize))

!  Read prognostic variables

 CALL end_io_surf_n(yprogram)
 CALL set_surfex_filein(yprogram,'PREP') ! change input file name to pgd name
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        yprogram,'FULL  ','SURF  ','READ ')

IF (ysc%U%NSIZE_NATURE>0 .AND. ysc%U%CNATURE/='NONE') THEN
  CALL read_surf(&
               yprogram,'WG1',       pws,   iresp)
  CALL read_surf(&
               yprogram,'WG2',       pwp,   iresp)
  CALL read_surf(&
               yprogram,'TG1',       pts,   iresp)
  CALL read_surf(&
               yprogram,'TG2',       ptp,   iresp)
  CALL read_surf(&
               yprogram,'WGI2',      ptl,   iresp)

  CALL read_surf(&
               yprogram,'VERSION',iversion,iresp)
  CALL read_surf(&
               yprogram,'BUG',ibugfix,iresp)
  IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=3) THEN
    CALL read_surf(&
               yprogram,'WSN_VEG1',psns,  iresp)
  ELSE
    CALL read_surf(&
               yprogram,'WSNOW_VEG1',psns,  iresp)
  ENDIF
ENDIF

IF (ysc%U%NSIZE_SEA>0 .AND. ysc%U%CSEA/='NONE') THEN
  CALL read_surf(&
               yprogram,'SST',       psst,  iresp)
ENDIF

IF (ysc%U%NSIZE_WATER>0 .AND. ysc%U%CWATER/='NONE') THEN
  CALL read_surf(&
               yprogram,'TS_WATER',  plst,  iresp)
ENDIF

IF (ysc%U%NSIZE_TOWN>0 .AND. ysc%U%CTOWN/='NONE' .AND. larome) THEN
  CALL read_surf(&
               yprogram,'VERSION',iversion,iresp)
  CALL read_surf(&
               yprogram,'BUG',ibugfix,iresp)
  IF (iversion>7 .OR. iversion==7 .AND. ibugfix>=3) THEN
    CALL read_surf(&
               yprogram,'TROAD3',   ptrd3,  iresp)
  ELSE
    CALL read_surf(&
               yprogram,'T_ROAD3',   ptrd3,  iresp)
  ENDIF
ELSE
  ptrd3(:) = xundef
ENDIF

 CALL read_surf(&
               yprogram,'T2M',       ptcls, iresp)
 CALL read_surf(&
               yprogram,'HU2M',      phcls, iresp)

! Read constant surface fields

 CALL read_surf(&
               yprogram,'RSMIN',     prsmin,iresp)
 CALL read_surf(&
               yprogram,'DG2',       pd2,   iresp)
 CALL read_surf(&
               yprogram,'LAI',       plai,  iresp)
 CALL read_surf(&
               yprogram,'VEG',       pveg,  iresp)

IF (nprintlev>0) THEN
  jj = ysc%U%NR_NATURE(1)
  print *,'value in PREP file => WG1       ',pws(jj,1)
  print *,'value in PREP file => WG2       ',pwp(jj,1)
  print *,'value in PREP file => TG1       ',pts(jj,1)
  print *,'value in PREP file => TG2       ',ptp(jj,1)
  print *,'value in PREP file => WGI2      ',ptl(jj,1)
  print *,'value in PREP file => WSNOW_VEG1',psns(jj,1)
  print *,'value in PREP file => SST       ',psst(jj)
  print *,'value in PREP file => LAI       ',plai(jj,1)
  print *,'value in PREP file => VEG       ',pveg(jj,1)
  print *,'value in PREP file => RSMIN     ',prsmin(jj,1)
  print *,'value in PREP file => DATA_DG2  ',pd2(jj,1)
  print *,'value in PREP file => SAND      ',psab(jj)
  print *,'value in PREP file => CLAY      ',parg(jj)
  print *,'value in PREP file => ZS        ',zalt(jj)
ENDIF

 CALL end_io_surf_n(yprogram)
 CALL io_buff_clean

!  Interface (allocate arrays)

ALLOCATE (plat(isize))
ALLOCATE (plon(isize))
ALLOCATE (zpx(isize))
ALLOCATE (zpy(isize))
ALLOCATE (pevaptr(isize))
ALLOCATE (zwp(isize))
ALLOCATE (zws(isize))
ALLOCATE (ztl(isize))
ALLOCATE (zts(isize))
ALLOCATE (ztp(isize))
ALLOCATE (zsns(isize))
ALLOCATE (ztcls(isize))
ALLOCATE (zhcls(isize))
ALLOCATE (zucls(isize))
ALLOCATE (zvcls(isize))
ALLOCATE (psstc(isize))
ALLOCATE (pwpinc1(isize))
ALLOCATE (pwpinc2(isize))
ALLOCATE (pwpinc3(isize))
ALLOCATE (pt2mbias(isize))
ALLOCATE (ph2mbias(isize))
ALLOCATE (prrcn(isize))
ALLOCATE (prrcl(isize))
ALLOCATE (prrsn(isize))
ALLOCATE (prrsl(isize))
ALLOCATE (patmneb(isize))
ALLOCATE (pitm(isize))
ALLOCATE (palbf(isize))
ALLOCATE (pemisf(isize))
ALLOCATE (pz0f(isize))
ALLOCATE (piveg(isize))
ALLOCATE (pz0h(isize))
ALLOCATE (ptsc(isize))
ALLOCATE (ptpc(isize))
ALLOCATE (pwsc(isize))
ALLOCATE (pwpc(isize))
ALLOCATE (psnc(isize))
ALLOCATE (zevap(isize))
ALLOCATE (zevaptr(isize))
ALLOCATE (pgelat(isize))
ALLOCATE (pgelam(isize))
ALLOCATE (pgemu(isize))
ALLOCATE (pt2m_o(isize))
ALLOCATE (phu2m_o(isize))
ALLOCATE (pts_o(isize))
ALLOCATE (psm_o(isize))
ALLOCATE (psig_smo(isize))
ALLOCATE (plsm_o(isize))
ALLOCATE (pws_o(isize))
ALLOCATE (zwginc(isize))

IF (.NOT. ldinline) THEN

!  Read atmospheric forecast fields from FA files 
#ifdef SFX_FA
  cfilein_fa = 'FG_OI_MAIN'        ! input forecast
  cfilein_fa_save  = cfilein_fa
#endif
!  Open FA file (LAM version with extension zone)
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                     yprogram2,'EXTZON','SURF  ','READ ')
ENDIF

!  Read model forecast quantities

IF (larome) THEN
  IF (ldinline) THEN
    prrsl(:) = p__surfaccpluie(1:isize)
    prrsn(:) = p__surfaccneige(1:isize)
    prrcn(:) = p__surfaccgraupel(1:isize)
  ELSE
    CALL read_surf(&
               yprogram2,'SURFACCPLUIE',    prrsl  ,iresp)
    CALL read_surf(&
               yprogram2,'SURFACCNEIGE',    prrsn  ,iresp)
    CALL read_surf(&
               yprogram2,'SURFACCGRAUPEL',  prrcn  ,iresp)
  ENDIF
  prrcl(:) = 0.0
!   CALL READ_SURF(YPROGRAM2,'SURFIND.VEG.DOMI',PIVEG  ,IRESP) 
  piveg(:) = 0.0
ELSE
  IF (ldinline) THEN
    prrcl(:) = p__surfprec_eau_con(1:isize)
    prrsl(:) = p__surfprec_eau_gec(1:isize)
    prrcn(:) = p__surfprec_nei_con(1:isize)
    prrsn(:) = p__surfprec_nei_gec(1:isize)
  ELSE
    CALL read_surf(&
               yprogram2,'SURFPREC.EAU.CON',prrcl  ,iresp)
    CALL read_surf(&
               yprogram2,'SURFPREC.EAU.GEC',prrsl  ,iresp)
    CALL read_surf(&
               yprogram2,'SURFPREC.NEI.CON',prrcn  ,iresp)
    CALL read_surf(&
               yprogram2,'SURFPREC.NEI.GEC',prrsn  ,iresp)
  ENDIF
  piveg(:) = 0.0
ENDIF
IF (ldinline) THEN
  patmneb(:) = p__atmonebul_basse(1:isize)
  pitm(:)    = p__surfind_terremer(1:isize)
  pevap(:)   = p__surfflu_lat_meva(1:isize)
ELSE
  CALL read_surf(&
               yprogram2,'ATMONEBUL.BASSE ',patmneb,iresp)
  CALL read_surf(&
               yprogram2,'SURFIND.TERREMER',pitm   ,iresp) 
  CALL read_surf(&
               yprogram2,'SURFFLU.LAT.MEVA',pevap  ,iresp) ! accumulated fluxes (not available in LFI)
ENDIF
IF (.NOT.laladsurf) THEN
  IF (ldinline) THEN
    pevaptr(:) = p__surfxevapotransp(1:isize)
  ELSE
    CALL read_surf(&
               yprogram2,'SURFXEVAPOTRANSP',pevaptr,iresp) ! not in ALADIN SURFEX
  ENDIF
ELSE
  pevaptr(:) = 0.0
ENDIF

IF (.NOT. ldinline) THEN
!  Close FA file
  CALL end_io_surf_n(yprogram2)
  CALL io_buff_clean
  print *,'READ FG_OI_MAIN OK'
ENDIF

IF (.NOT. ldinline) THEN
!  Define FA file name for CANARI analysis
#ifdef SFX_FA
  cfilein_fa = 'CANARI'        ! input CANARI analysis
  cfilein_fa_save  = cfilein_fa
#endif
!  Open FA file 
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                     yprogram2,'EXTZON','SURF  ','READ ')
ENDIF

IF (ldinline) THEN
  pt2m_o(:)  = p__clstemperature(1:isize)
  phu2m_o(:) = p__clshumi_relative(1:isize)
  pts_o(:)   = p__surftemperature(1:isize)
  pucls(:)   = p__clsvent_zonal(1:isize)
  pvcls(:)   = p__clsvent_meridien(1:isize)
ELSE
!  Read CANARI analysis
  CALL read_surf(&
               yprogram2,'CLSTEMPERATURE  ',pt2m_o ,iresp)
  CALL read_surf(&
               yprogram2,'CLSHUMI.RELATIVE',phu2m_o,iresp)
  CALL read_surf(&
               yprogram2,'SURFTEMPERATURE ',pts_o  ,iresp)
  CALL read_surf(&
               yprogram2,'CLSVENT.ZONAL   ',pucls  ,iresp)
  CALL read_surf(&
               yprogram2,'CLSVENT.MERIDIEN',pvcls  ,iresp)
ENDIF

IF (.NOT. ldinline) THEN
!  Close CANARI file
  CALL end_io_surf_n(yprogram2)
  CALL io_buff_clean
  print *,'READ CANARI OK'
ENDIF

!  Read ASCAT SM observations (in percent)

inobs = 0
IF (lobswg) THEN
  OPEN(unit=111,file='ASCAT_SM.DAT')
  DO jj=1,ysc%U%NDIM_FULL
    READ(111,*,end=990) psm_o(jj),psig_smo(jj),plsm_o(jj)
    IF (plsm_o(jj) < 1.0)          psm_o(jj) = 999.0 ! data rejection if not on land
    IF (psig_smo(jj) > xsigwgo_max) psm_o(jj) = 999.0 ! data rejection of error too large
    IF (psm_o(jj) /= 999.0) inobs = inobs + 1
  ENDDO
990 CONTINUE
  CLOSE(unit=111)
  print *,'READ ASCAT SM OK'
ELSE
  psm_o(:)    = 999.0
  psig_smo(:) = 999.0
  plsm_o(:)   = 0.0
ENDIF
print *,' NUMBER OF ASCAT OBSERVATIONS AFTER INITIAL CHECKS  :: ',inobs
inobs = 0

! Perform bias correction of SM observations

 CALL oi_bc_soil_moisture(isize,psm_o,psab,pws_o)

IF (.NOT. ldinline) THEN
!  Define FA file name for surface climatology
#ifdef SFX_FA
  cfilein_fa = 'clim_isba'               ! input climatology
  cfilein_fa_save  = cfilein_fa
  cdnomc     = 'climat'                  ! new frame name
#endif
!  Open FA file 
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                     yprogram2,'EXTZON','SURF  ','READ ')
ENDIF

IF (ldinline) THEN
  psnc(:) = p__surfreserv_neige(1:isize)
ELSE
!  Read climatology file (snow water equivalent)
  CALL read_surf(&
               yprogram2,'SURFRESERV.NEIGE',psnc  ,iresp)
ENDIF

IF (.NOT. ldinline) THEN
!  Close climatology file
  CALL end_io_surf_n(yprogram2)
  CALL io_buff_clean
  print *,'READ CLIMATOLOGY OK'
ENDIF

IF (.NOT. ldinline) THEN
#ifdef SFX_FA
  plat0  = pelat0 
  plon0  = pelon0 
  plator = pelat1 
  plonor = pelon1 
  prpk   = perpk  
  pbeta  = pebeta 
  delx   = pedelx 
  dely   = pedely 
  IF (plonor > 180.0) plonor = plonor - 360.0
  IF (plon0  > 180.0) plon0  = plon0  - 360.0
  DO jj=1,ndgux
    DO j1=1,ndlux
      zpx((jj-1)*ndlux + j1) = delx*REAL(j1-1)
      zpy((jj-1)*ndlux + j1) = dely*REAL(jj-1)
    ENDDO
  ENDDO
#endif
  CALL oi_latlon_conf_proj(isize,plat0,plon0,prpk,pbeta,plator,plonor,zpx,zpy,plat,plon)
ELSE
  plat(:) = p__lat(1:isize)
  plon(:) = p__lon(1:isize)
ENDIF

!  Allocate arrays to produce analysis increments  

ALLOCATE (zt2inc(isize))
ALLOCATE (zh2inc(isize))
ALLOCATE (zwsinc(isize))
ALLOCATE (zwpinc(isize))
ALLOCATE (ztlinc(isize))
ALLOCATE (ztsinc(isize))
ALLOCATE (ztpinc(isize))
ALLOCATE (zsninc(isize))

! Screen-level innovations

zt2inc(:) = pt2m_o(:) - ptcls(:)
zh2inc(:) = phu2m_o(:) - phcls(:)

! Threshold for background check

zthres=xrthr_qc*sqrt(xsigwgo**2 + xsigwgb**2)

! Superficial soil moisture innovations in (m3/m3)

DO jj = 1, isize
  IF (pws_o(jj) /= 999.0) THEN
    zwginc(jj) = pws_o(jj) - pws(jj,1)
    IF (abs(zwginc(jj)) > zthres) THEN 
      zwginc(jj) = 0.0 ! background check
    ELSE
      inobs = inobs + 1
    ENDIF
  ELSE
    zwginc(jj) = 0.0
  ENDIF
ENDDO
print *,' NUMBER OF ASCAT OBSERVATIONS AFTER BACKGROUND CHECK  :: ',inobs

! Interface (define arrays and perform unit conversions)

parg(:) = parg(:)*100.0
psab(:) = psab(:)*100.0

zws(:) = xundef
zwp(:) = xundef
ztl(:) = xundef

WHERE (pws(:,1)/=xundef) 
  zws(:)      = pws(:,1)*xrd1*xrholw     ! conversion of m3/m3 -> mm
  zwp(:)      = pwp(:,1)*pd2(:,1)*xrholw  ! conversion of m3/m3 -> mm
  ztl(:)      = ptl(:,1)*pd2(:,1)*xrholw  ! conversion of m3/m3 -> mm
END WHERE

ztcls(:)    = ptcls(:)
zhcls(:)    = phcls(:)
zucls(:)    = pucls(:)
zvcls(:)    = pvcls(:)
psstc(:)    = pts_o(:)
pwpinc1(:)  = xundef
pwpinc2(:)  = xundef
pwpinc3(:)  = xundef
pt2mbias(:) = xundef
ph2mbias(:) = xundef

! Sea-ice surface properties

palbf(:)    = xundef
pemisf(:)   = xundef
pz0f(:)     = xundef
pz0h(:)     = xundef

! Climatological arrays set to missing values

psnc(:)     =  psns(:,1) ! need to read the snow climatology
pwsc(:)     =  xundef
pwpc(:)     =  xundef
ptsc(:)     =  xundef
ptpc(:)     =  xundef

DO jj = 1, isize
  pgelat(jj)   = plat(jj)
  pgelam(jj)   = plon(jj)
  pgemu(jj)    = sin(plat(jj)*zpis180)
ENDDO

zevap(:)   =  (pevap(:)/xlvtt*xday)/(nechgu*3600.) ! conversion W/m2 -> mm/day
zevaptr(:) =  pevaptr(:)*xday
zsns(:)    =  psns(:,1)

DO jj = 1, isize
  zts(jj) = pts(jj,1)
  ztp(jj) = ptp(jj,1)
ENDDO

idat = iyear*10000. + imonth*100. + iday

IF (ldinline) THEN
! Avoid division by zero in next WHERE statement; 
! this may occur in the extension zone
  WHERE (ld_maskext(1:isize))
    pd2(:,1) = 1.0
    zt2inc(:) = 0.0_jprb
    zh2inc(:) = 0.0_jprb
  END WHERE
ENDIF

print *,'           '
print *,'Mean T2m increments  ',sum(zt2inc)/ysc%U%NDIM_FULL
print *,'Mean HU2m increments ',sum(zh2inc)/ysc%U%NDIM_FULL
print *,'           '

!  Soil analysis based on optimal interpolation

 CALL oi_cacsts(isize,zt2inc,zh2inc,zwginc,pws_o,                      &
                idat,nsssss,                                           &
                ztp,zwp,ztl,zsns,zts,zws,                              &
                ztcls,zhcls,zucls,zvcls,psstc,pwpinc1,pwpinc2,pwpinc3, &
                pt2mbias,ph2mbias,                                     &
                prrcl,prrsl,prrcn,prrsn,patmneb,zevap,zevaptr,         &
                pitm,pveg(:,1),palbf,pemisf,pz0f,                      &
                piveg,parg,pd2(:,1),psab,plai(:,1),prsmin(:,1),pz0h,   &
                ptsc,ptpc,pwsc,pwpc,psnc,                              &
                pgelat,pgelam,pgemu) 

!  Store increments

zwsinc(:) = 0.0_jprb
zwpinc(:) = 0.0_jprb
ztlinc(:) = 0.0_jprb
zsninc(:) = 0.0_jprb

WHERE (pws(:,1)/=xundef)
  zwsinc(:) = zws(:) - pws(:,1)*(xrd1*xrholw)    
  zwpinc(:) = zwp(:) - pwp(:,1)*(pd2(:,1)*xrholw) 
  ztlinc(:) = ztl(:) - ptl(:,1)*(pd2(:,1)*xrholw) 
  zsninc(:) = zsns(:) - psns(:,1)
END WHERE

IF (ldinline) THEN
! Avoid division by zero in next WHERE statement; 
! this may occur in the extension zone
  WHERE (ld_maskext(1:isize))
    pd2(:,1) = 1.0
  END WHERE
ENDIF

!  Define soil moiture analyses over NATURE points

WHERE (pws(:,1)/=xundef)
  pws(:,1)  = zws(:)/(xrd1*xrholw)
  pwp(:,1)  = zwp(:)/(pd2(:,1)*xrholw)
  ptl(:,1)  = ztl(:)/(pd2(:,1)*xrholw)
  psns(:,1) = zsns(:)
END WHERE

!  Perform temperature analysis according to surface types

ointerp_lst(:) = .false.
ointerp_sst(:) = .false.

ztsinc(:) = 0.0_jprb
ztpinc(:) = 0.0_jprb

! a) Temperature analysis of NATURE points

WHERE (pts(:,1)/=xundef)
  ztsinc(:) = zts(:) - pts(:,1)
  ztpinc(:) = ztp(:) - ptp(:,1)
  pts(:,1)  = zts(:)
  ptp(:,1)  = ztp(:)
END WHERE

! b) Temperature analysis of SEA and LAKE points

DO jj = 1, isize
  IF (pitm(jj) < 0.5_jprb) THEN
    IF (psst(jj)/=xundef) THEN
      ztsinc(jj) = pts_o(jj) - psst(jj)
      psst(jj) = pts_o(jj)   ! canari
    ENDIF
    IF (plst(jj)/=xundef) THEN
      plst(jj) = pts_o(jj)   ! canari
    ENDIF
  ELSE
    IF (psst(jj)/=xundef) THEN
      psst(jj) = xundef
      ointerp_sst(jj) = .true.
    ENDIF
    IF (plst(jj)/=xundef) THEN
      plst(jj) = xundef
      ointerp_lst(jj) = .true.
    ENDIF
  ENDIF
ENDDO

! c) Temperature analysis of TOWN points

WHERE (ptrd3(:)/=xundef)
  ptrd3(:) = ptrd3(:) + zt2inc(:)*z1s2pi
END WHERE

! Search for the nearest grid point values for lake and sea points
! at locations where the water fraction is less than 50 % 
! and therefore no useful information is given from the SST analysis
! A standard temperature gradient is applied to account for the atitude differences

IF (ldinline) THEN

  IF (llkeepextzone) THEN

    zlst(:) = plst(:)
     
    IF (ldinline) THEN
      WHERE (ld_maskext(1:isize))
        zlst = xundef
      END WHERE
    ENDIF
     
    CALL oi_hor_extrapol_surf(isize,plat,plon,zlst,plat,plon,plst,ointerp_lst,zalt)
     
    zsst(:) = psst(:)
     
    IF (ldinline) THEN
      WHERE (ld_maskext(1:isize))
        zsst = xundef
      END WHERE
    ENDIF
     
    CALL oi_hor_extrapol_surf(isize,plat,plon,zsst,plat,plon,psst,ointerp_sst,zalt)

  ELSE

    isize1 = count(.NOT. ld_maskext)
   
    ALLOCATE (psst1(isize1), plst1(isize1), zsst1(isize1), zlst1(isize1), plat1(isize1), &
            & plon1(isize1), zalt1(isize1), ointerp_lst1(isize1), ointerp_sst1(isize1))
   
    ! remove extension zone
    jj = 1
    DO j1 = 1, isize
      IF (.NOT. ld_maskext(j1)) THEN
        psst1(jj) = psst(j1)
        plst1(jj) = plst(j1)
        plat1(jj) = plat(j1)
        plon1(jj) = plon(j1)
        zalt1(jj) = zalt(j1)
        ointerp_lst1(jj) = ointerp_lst(j1)
        ointerp_sst1(jj) = ointerp_sst(j1)
        jj = jj + 1
      ENDIF
    ENDDO
   
    zlst1(:) = plst1(:)
    CALL oi_hor_extrapol_surf(isize1,plat1,plon1,zlst1,plat1,plon1,plst1,ointerp_lst1,zalt1)
     
    zsst1(:) = psst1(:)
    CALL oi_hor_extrapol_surf(isize1,plat1,plon1,zsst1,plat1,plon1,psst1,ointerp_sst1,zalt1)
   
    ! copy back
    jj = 1
    DO j1 = 1, isize
      IF (.NOT. ld_maskext(j1)) THEN
        psst(j1) = psst1(jj)
        plst(j1) = plst1(jj) 
        jj = jj + 1
      ENDIF
    ENDDO
   
    DEALLOCATE (psst1, plst1, zsst1, zlst1, plat1, plon1, &
              & zalt1, ointerp_lst1, ointerp_sst1)

  ENDIF

ELSE
 
  zlst(:) = plst(:)
  CALL oi_hor_extrapol_surf(isize,plat,plon,zlst,plat,plon,plst,ointerp_lst,zalt)

  zsst(:) = psst(:)
  CALL oi_hor_extrapol_surf(isize,plat,plon,zsst,plat,plon,psst,ointerp_sst,zalt)

ENDIF

! PRINT values produced by OI_HO_EXTRAPOL_SURF

IF (nprintlev>0) THEN
  DO jj = 1, isize
    IF (ointerp_lst(jj)) THEN
      print *,'Lake surface temperature set to ',plst(jj),'from nearest neighbour at J=',jj
    ENDIF
    IF (ointerp_sst(jj)) THEN
      print *,'Sea surface temperature set to ',psst(jj),'from nearest neighbour at J=',jj
    ENDIF
  ENDDO
ENDIF

! PRINT statistics of the soil analysis

print *,'---------------------------------------------------------------'
print *,'Mean WS increments over NATURE ',sum(zwsinc,ysc%U%XNATURE > 0.)/ysc%U%NDIM_NATURE
print *,'Mean WP increments over NATURE ',sum(zwpinc,ysc%U%XNATURE > 0.)/ysc%U%NDIM_NATURE
print *,'Mean TS increments over NATURE ',sum(ztsinc,ysc%U%XNATURE > 0.)/ysc%U%NDIM_NATURE
print *,'Mean TP increments over NATURE ',sum(ztpinc,ysc%U%XNATURE > 0.)/ysc%U%NDIM_NATURE
print *,'Mean TL increments over NATURE ',sum(ztlinc,ysc%U%XNATURE > 0.)/ysc%U%NDIM_NATURE
print *,'Mean SN increments over NATURE ',sum(zsninc,ysc%U%XNATURE > 0.)/ysc%U%NDIM_NATURE
print *,'---------------------------------------------------------------'
print *,'Mean WS increments over SEA    ',sum(zwsinc,ysc%U%XSEA > 0.)/ysc%U%NDIM_SEA
print *,'Mean WP increments over SEA    ',sum(zwpinc,ysc%U%XSEA > 0.)/ysc%U%NDIM_SEA
print *,'Mean TS increments over SEA    ',sum(ztsinc,ysc%U%XSEA > 0.)/ysc%U%NDIM_SEA
print *,'Mean TP increments over SEA    ',sum(ztpinc,ysc%U%XSEA > 0.)/ysc%U%NDIM_SEA
print *,'Mean TL increments over SEA    ',sum(ztlinc,ysc%U%XSEA > 0.)/ysc%U%NDIM_SEA
print *,'Mean SN increments over SEA    ',sum(zsninc,ysc%U%XSEA > 0.)/ysc%U%NDIM_SEA
print *,'---------------------------------------------------------------'

IF (.NOT. ldinline) THEN 
!   Write analysis in LFI file PREP
#ifdef SFX_LFI
  cfileout_lfi='PREP'
#endif
ENDIF

 CALL flag_update(ysc%IM%DGI, ysc%DGU, &
                  .false.,.false.,.true.,.false.)
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        yprogram,'FULL  ','SURF  ','WRITE')

IF (ldinline) THEN
#ifdef SFX_ARO
! Count 2D fields in MSE
! NINDX1, NINDX2, NKPROMA already set 
  CALL surfex_field_buf_set_record(ysurfex_cache_out, .false.)

  CALL WRITE
   
  CALL surfex_field_buf_prealloc(ysurfex_cache_out)
  CALL surfex_field_buf_set_record(ysurfex_cache_out, .true.)

#endif

ENDIF

 CALL WRITE

 CALL end_io_surf_n(yprogram)
 CALL io_buff_clean

DEALLOCATE(nwork)
DEALLOCATE(xwork)
DEALLOCATE(nwork2)
DEALLOCATE(xwork2)
DEALLOCATE(xwork3)
DEALLOCATE(nwork_full)
DEALLOCATE(xwork_full)
DEALLOCATE(nwork2_full)
DEALLOCATE(xwork2_full)

IF (.NOT. ldinline) THEN
  print *,'after write in PREP file'
ENDIF

! -------------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('OI_CONTROL', 1, zhook_handle)

 CONTAINS

SUBROUTINE write 

REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('OI_CONTROL:WRITE', 0, zhook_handle)

 CALL dd('WG1', pws(:,1))

yvar='WG1'
yprefix='X_Y_WG1 (m3/m3)                                   '
 CALL write_surf(ysc%DGU, ysc%U, &
                yprogram,yvar,pws,iresp,hcomment=yprefix)

 CALL dd('WG2', pwp(:,1))

yvar='WG2'
yprefix='X_Y_WG2 (m3/m3)                                   '
 CALL write_surf(ysc%DGU, ysc%U, &
                yprogram,yvar,pwp,iresp,hcomment=yprefix)

 CALL dd('WGI2', ptl(:,1))

yvar='WGI2'
yprefix='X_Y_WGI2 (m3/m3)                                  '
 CALL write_surf(ysc%DGU, ysc%U, &
                yprogram,yvar,ptl,iresp,hcomment=yprefix)

 CALL dd('TG1', pts(:,1))

yvar='TG1'
yprefix='X_Y_TG1 (K)                                       '
 CALL write_surf(ysc%DGU, ysc%U, &
                yprogram,yvar,pts,iresp,hcomment=yprefix)

 CALL dd('TG2', ptp(:,1))

yvar='TG2'
yprefix='X_Y_TG2 (K)                                       '
 CALL write_surf(ysc%DGU, ysc%U, &
                yprogram,yvar,ptp,iresp,hcomment=yprefix)

 CALL dd('SST', psst)

yvar='SST'
yprefix='X_Y_SST (K)                                       '
 CALL write_surf(ysc%DGU, ysc%U, &
                yprogram,yvar,psst,iresp,hcomment=yprefix)

 CALL dd('TS_WATER', plst)

yvar='TS_WATER'
yprefix='X_Y_TS_WATER (K)                                  '
 CALL write_surf(ysc%DGU, ysc%U, &
                yprogram,yvar,plst,iresp,hcomment=yprefix)

IF (ysc%U%NSIZE_TOWN > 0 .AND. larome) THEN
  CALL dd('T_ROAD3', ptrd3)

  yvar='TROAD3'
  yprefix='X_Y_T_ROAD3 (K)                                   '
  CALL write_surf(ysc%DGU, ysc%U, &
                yprogram,yvar,ptrd3,iresp,hcomment=yprefix)
ENDIF

 CALL dd('WSNOW_VEG1', psns(:,1))

yvar='WSN_VEG1'
yprefix='X_Y_WSNOW_VEG1 (kg/m2)                            '
 CALL write_surf(ysc%DGU, ysc%U, &
                yprogram,yvar,psns,iresp,hcomment=yprefix)

IF (lhook) CALL dr_hook('OI_CONTROL:WRITE', 1, zhook_handle)

END SUBROUTINE write

SUBROUTINE dd (CDN, PX)
 CHARACTER(LEN=*), INTENT (IN) :: cdn
REAL, INTENT (IN) :: px (:)

REAL :: zx (size (px))
INTEGER :: ji, jn
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('OI_CONTROL:DD', 0, zhook_handle)

IF (ldinline) THEN
#ifdef SFX_ARO
  IF (lhook) CALL dr_hook('OI_CONTROL:DD', 1, zhook_handle)
  RETURN
#endif
  jn = count(.NOT. ld_maskext)
  zx(1:jn) = pack(px, .NOT. ld_maskext)
ELSE
  zx = px
  jn = SIZE (px)
ENDIF

WRITE (0, *) trim(cdn)//" = " 
WRITE (0, *) jn, minval(zx(1:jn)), maxval(zx(1:jn))
!WRITE (0, '(10(E14.6,", "))') ZX (1:N)

IF (lhook) CALL dr_hook('OI_CONTROL:DD', 1, zhook_handle)

END SUBROUTINE dd

END SUBROUTINE oi_control