assim_nature_isba_oi.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 assim_nature_isba_oi (I, &
                                 hprogram, ki,                     &
                                prrcl,    prrsl,  prrcn,   prrsn, &
                                patmneb,  pitm,   pevaptr, pevap, &
                                psnc,     ptsc,   pucls, pvcls,   &
                                pts_o,    pt2m_o, phu2m_o, pswe,  &
                                htest, od_maskext,      &
                                plon_in, plat_in )

! ------------------------------------------------------------------------------------------
!  *****************************************************************************************
!
!  Routine to perform OI 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)
!   (07/2011)  : Read pgd+prep (B. Decharme)
!   (04/2012)  : Made as a subroutine (T. Aspelien)
!   (06/2013)  : Separating IO (T. Aspelien)
! ******************************************************************************************
! ------------------------------------------------------------------------------------------
!
USE modd_isba_n, ONLY : isba_t
!
USE modd_csts,            ONLY : xday, xpi, xrholw, xlvtt, ndaysec
USE modd_surf_par,        ONLY : xundef 
!
USE modd_assim,           ONLY : lobswg, nitrad, nprintlev, nechgu, xrd1, xrscaldw,   &
                                 xrthr_qc, xsigwgb, xsigwgo, xsigwgo_max, xat2m_isba, &
                                 xahu2m_isba, xazon10m_isba, xamer10m_isba
!

USE yomhook,              ONLY : lhook,   dr_hook
USE parkind1,             ONLY : jprb
!
USE modi_abor1_sfx
USE modi_oi_bc_soil_moisture
USE modi_oi_cacsts
!
IMPLICIT NONE
!
TYPE(isba_t), INTENT(INOUT) :: i
!
 CHARACTER(LEN=6),    INTENT(IN) :: hprogram  ! program calling surf. schemes
INTEGER,             INTENT(IN) :: ki
REAL, DIMENSION(KI), INTENT(IN) :: prrcl
REAL, DIMENSION(KI), INTENT(IN) :: prrsl
REAL, DIMENSION(KI), INTENT(IN) :: prrcn
REAL, DIMENSION(KI), INTENT(IN) :: prrsn
REAL, DIMENSION(KI), INTENT(IN) :: patmneb
REAL, DIMENSION(KI), INTENT(IN) :: pitm
REAL, DIMENSION(KI), INTENT(IN) :: pevaptr
REAL, DIMENSION(KI), INTENT(IN) :: pevap
REAL, DIMENSION(KI), INTENT(IN) :: psnc
REAL, DIMENSION(KI), INTENT(IN) :: ptsc
REAL, DIMENSION(KI), INTENT(IN) :: pucls
REAL, DIMENSION(KI), INTENT(IN) :: pvcls
REAL, DIMENSION(KI), INTENT(IN) :: pts_o
REAL, DIMENSION(KI), INTENT(IN) :: pt2m_o
REAL, DIMENSION(KI), INTENT(IN) :: phu2m_o
REAL, DIMENSION(KI), INTENT(OUT):: pswe
 CHARACTER(LEN=2),    INTENT(IN) :: htest ! must be equal to 'OK'
LOGICAL,  DIMENSION (KI) ::  od_maskext
REAL(KIND=JPRB), DIMENSION (:), INTENT(IN) ::  plon_in
REAL(KIND=JPRB), DIMENSION (:), INTENT(IN) ::  plat_in

!    Declarations of local variables
!
! Arrays for soil OI analysis
!
REAL, DIMENSION (KI) :: zsab, zarg, zws, zwp, ztl, zws0, zwp0, ztl0,         &
                        zts, ztp, zsns, ztcls, zhcls, zd2, zucls, zvcls,     &
                        zrsmin, zlai, zveg, zsns0, zts0, ztp0,               &
                        ziveg,  zsm_o, zsig_smo, zlsm_o, zws_o, zlon, zlat,  &           
                        zt2inc, zh2inc, zwginc, zwpinc1, zwpinc2, zwpinc3,   &
                        zt2mbias, zh2mbias, zalbf, zemisf, zz0f, zz0h,       &
                        zwsc, zwpc, ztpc, zevap, zevaptr, zsstc,             &
                        zgelat, zgelam, zgemu,  zwsinc, zwpinc, ztlinc,      &
                        zsninc, ztsinc, ztpinc
!                 
REAL :: zthres
INTEGER  :: idat
INTEGER  :: iyear                      ! current year (UTC)
INTEGER  :: imonth                     ! current month (UTC)
INTEGER  :: iday                       ! current day (UTC)
INTEGER  :: isssss                     ! current time since start of the run (s)
INTEGER  :: ji, jj, jp, jl
INTEGER  :: inobs   ! number of observations
REAL(KIND=JPRB)      :: zhook_handle
!
! ----------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('ASSIM_NATURE_ISBA_OI',0,zhook_handle)

IF (htest/='OK') THEN
  CALL abor1_sfx('ASSIM_NATURE_ISBA_OI: FATAL ERROR DURING ARGUMENT TRANSFER')
END IF

IF ( nprintlev > 0 ) THEN 
  WRITE(*,*) '--------------------------------------------------------------------------'
  WRITE(*,*) '|                                                                        |'
  WRITE(*,*) '|                             ENTER OI_ASSIM                             |'
  WRITE(*,*) '|                                                                        |'
  WRITE(*,*) '--------------------------------------------------------------------------'
ENDIF
!
!   Update some constants dependant from NACVEG
!
!  scaling of soil moisture increments when assimilation window is different
!  from 6 hours
xrscaldw = REAL(nechgu)/6.0
!  half assimilation window in sec
nitrad = nechgu*1800
!
!   Time initializations 
!
iyear  = i%TTIME%TDATE%YEAR
imonth = i%TTIME%TDATE%MONTH
iday   = i%TTIME%TDATE%DAY
isssss = i%TTIME%TIME
IF ( isssss>ndaysec ) isssss = isssss - ndaysec
idat = iyear*10000. + imonth*100. + iday
!
!
DO ji=1,ki     
  zgelam(ji) = plon_in(ji)  
  zgelat(ji) = plat_in(ji)
ENDDO
!
DO ji=1,ki
  zgemu(ji) = sin(zgelat(ji)*xpi/180.)
ENDDO
!
!
jp = 1
jl = 1
!
!
zsab(:) = i%XSAND(:,jp)*100.
zarg(:) = i%XCLAY(:,jp)*100.
!
zts0(:) = i%XTG  (:,1,jp)
ztp0(:) = i%XTG  (:,2,jp)
!
zws0(:) = i%XWG  (:,1,jp)
zwp0(:) = i%XWG  (:,2,jp)
ztl0(:) = i%XWGI (:,2,jp)
!
zsns0(:) = i%TSNOW%WSNOW(:,jl,jp)
zsns(:) = zsns0(:)
!
ztcls(:) = xat2m_isba(:,jp)
zhcls(:) = xahu2m_isba(:,jp)
!
zd2(:) = i%XDG   (:,2,jp)
zrsmin(:) = i%XRSMIN(:,jp)
zlai(:) = i%XLAI  (:,jp)
zveg(:) = i%XVEG  (:,jp)
!
zucls(:) = pucls(:)
zvcls(:) = pvcls(:)
!
zwpinc1(:) = xundef
zwpinc2(:) = xundef
zwpinc3(:) = xundef
zt2mbias(:) = xundef
zh2mbias(:) = xundef
!
! Sea-ice surface properties
zalbf(:) = xundef
zemisf(:) = xundef
zz0f(:) = xundef
zz0h(:) = xundef
!
! Climatological arrays set to missing values
zwsc(:) = xundef
zwpc(:) = xundef
ztpc(:) = xundef
!
! PRINT 
!
IF ( nprintlev > 1 ) THEN
  WRITE(*,*) 'value in PREP file => WG1       ',sum(zws0)/ki
  WRITE(*,*) 'value in PREP file => WG2       ',sum(zwp0)/ki
  WRITE(*,*) 'value in PREP file => TG1       ',sum(zts0)/ki
  WRITE(*,*) 'value in PREP file => TG2       ',sum(ztp0)/ki
  WRITE(*,*) 'value in PREP file => WGI2      ',sum(ztl0)/ki
  WRITE(*,*) 'value in PREP file => WSNOW_VEG1',sum(zsns)/ki
  WRITE(*,*) 'value in PREP file => LAI       ',sum(zlai)/ki
  WRITE(*,*) 'value in PREP file => VEG       ',sum(zveg)/ki
  WRITE(*,*) 'value in PREP file => RSMIN     ',sum(zrsmin)/ki
  WRITE(*,*) 'value in PREP file => DATA_DG2  ',sum(zd2)/ki
  WRITE(*,*) 'value in PREP file => SAND      ',sum(zsab)/ki
  WRITE(*,*) 'value in PREP file => CLAY      ',sum(zarg)/ki
ENDIF
!
! SST not used in cacsts
zsstc(:)    = 0.
!
WHERE ( zws0(:)/=xundef )
  zws(:) = zws0(:) * xrd1   * xrholw  ! conversion of m3/m3 -> mm
  zwp(:) = zwp0(:) * zd2(:) * xrholw  ! conversion of m3/m3 -> mm
  ztl(:) = ztl0(:) * zd2(:) * xrholw  ! conversion of m3/m3 -> mm
ENDWHERE
!
zevap(:) =  (pevap(:)/xlvtt*xday)/(nechgu*3600.) ! conversion W/m2 -> mm/day
zevaptr(:) =  pevaptr(:)*xday 
!
! Set PIVEG (SURFIND.VEG.DOMI) since it is not available
ziveg(:) = 0.0
!
!
!  Read ASCAT SM observations (in percent)
!
inobs = 0
IF ( lobswg ) THEN
  OPEN(unit=111,file='ASCAT_SM.DAT')
  DO ji=1,ki
    READ(111,*) zsm_o(ji),zsig_smo(ji),zlsm_o(ji)
    ! data rejection if not on land
    ! data rejection of error too large
    IF ( zlsm_o(ji)<1.0 .OR. zsig_smo(ji)>xsigwgo_max ) zsm_o(ji) = 999.0 
    IF ( zsm_o(ji)/=999.0 ) inobs = inobs + 1
  ENDDO
  CLOSE(unit=111)
  IF ( nprintlev > 0 ) WRITE(*,*) 'READ ASCAT SM OK'
ELSE
  zsm_o(:) = 999.0
  zsig_smo(:) = 999.0
  zlsm_o(:) = 0.0
ENDIF
IF ( nprintlev > 0 ) WRITE(*,*) ' NUMBER OF ASCAT OBSERVATIONS AFTER INITIAL CHECKS  :: ',inobs
!
! Perform bias correction of SM observations
 CALL oi_bc_soil_moisture(ki,zsm_o,zsab,zws_o)
!
!
! Screen-level innovations
zt2inc=0.
WHERE ( pt2m_o(:) /= 999.0 ) 
  zt2inc(:) = pt2m_o(:) - ztcls(:)
END WHERE
zh2inc=0.
WHERE ( phu2m_o(:) /= 999.0 ) 
  zh2inc(:) = phu2m_o(:) - zhcls(:)
END WHERE
!
! Avoid division by zero in next WHERE statement; 
! this may occur in the extension zone
WHERE (od_maskext(1:ki))
  zd2(:) = 1.0
  zt2inc(:) = 0.0
  zh2inc(:) = 0.0
END WHERE
!
! Threshold for background check
zthres = xrthr_qc*sqrt(xsigwgo**2 + xsigwgb**2)
! Superficial soil moisture innovations in (m3/m3)
inobs = 0
DO ji=1,ki
  IF ( zws_o(ji)/=999.0 ) THEN
    zwginc(ji) = zws_o(ji) - zws(ji)
    IF ( abs(zwginc(ji))>zthres ) THEN 
      zwginc(ji) = 0.0 ! background check
    ELSE
      inobs = inobs + 1
    ENDIF
  ELSE
    zwginc(ji) = 0.0
  ENDIF
ENDDO
IF ( nprintlev > 0 ) THEN
  WRITE(*,*) ' NUMBER OF ASCAT OBSERVATIONS AFTER BACKGROUND CHECK  :: ',inobs
  WRITE(*,*) 'Mean T2m increments  ',sum(zt2inc)/ki
  WRITE(*,*) 'Mean HU2m increments ',sum(zh2inc)/ki
ENDIF
!
zts(:) = zts0(:)
ztp(:) = ztp0(:)
write(*,*) 'PERFORMING OI SOIL ANALYSIS'
 CALL oi_cacsts(ki, zt2inc, zh2inc, zwginc, zws_o,                   &
               idat, isssss,                                        &
               ztp, zwp, ztl, zsns, zts, zws,                       &
               ztcls, zhcls, zucls, zvcls, zsstc,                   &
               zwpinc1, zwpinc2, zwpinc3, zt2mbias, zh2mbias,       &
               prrcl, prrsl, prrcn, prrsn, patmneb, zevap, zevaptr, &
               pitm, zveg, zalbf, zemisf, zz0f,                     &
               ziveg, zarg, zd2, zsab, zlai, zrsmin, zz0h,          &
               ptsc, ztpc, zwsc, zwpc, psnc, zgelat, zgelam, zgemu )  
!
!  Perform soil moiture analyses
pswe(:) = zsns(:)
!
zwsinc(:) = 0.0
zwpinc(:) = 0.0
ztlinc(:) = 0.0
zsninc(:) = 0.0
WHERE ( zws(:)/=xundef )
  zwsinc(:) = zws(:) - zws0(:) * xrd1 * xrholw
  zwpinc(:) = zwp(:) - zwp0(:) * zd2(:) * xrholw
  ztlinc(:) = ztl(:) - ztl0(:) * zd2(:) * xrholw
  zsninc(:) = zsns(:) - zsns0(:)
END WHERE

! Avoid division by zero in next WHERE statement; 
! this may occur in the extension zone
WHERE (od_maskext(1:ki)) zd2(:) = 1.0
!
WHERE (zws0(:)/=xundef)
  zws0(:)  = zws(:)/ (xrd1*xrholw)
  zwp0(:)  = zwp(:)/ (zd2(:)*xrholw)
  ztl0(:)  = ztl(:)/ (zd2(:)*xrholw)
  zsns0(:) = zsns(:)
ENDWHERE
!  Perform temperature analyses
!
ztsinc(:) = 0.0
ztpinc(:) = 0.0
WHERE (zts(:)/=xundef)
  ztsinc(:) = zts(:) - zts0(:)
  ztpinc(:) = ztp(:) - ztp0(:)
  zts0(:)   = zts(:)
  ztp0(:)   = ztp(:)
END WHERE


! PRINT statistics of the soil analysis

WRITE(*,*) '---------------------------------------------------------------'
WRITE(*,*) 'Mean WS increments over NATURE ',sum(zwsinc)/ki
WRITE(*,*) 'Mean WP increments over NATURE ',sum(zwpinc)/ki
WRITE(*,*) 'Mean TS increments over NATURE ',sum(ztsinc)/ki
WRITE(*,*) 'Mean TP increments over NATURE ',sum(ztpinc)/ki
WRITE(*,*) 'Mean TL increments over NATURE ',sum(ztlinc)/ki
WRITE(*,*) '---------------------------------------------------------------'

! Update modified variables
i%XWG (:,1,jp) = zws0(:)
i%XWG (:,2,jp) = zwp0(:)
i%XTG (:,1,jp) = zts0(:)
i%XTG (:,2,jp) = ztp0(:)
i%XWGI(:,2,jp) = ztl0(:)

!
! -------------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('ASSIM_NATURE_ISBA_OI',1,zhook_handle)
END SUBROUTINE assim_nature_isba_oi