assim_isban.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_isba_n (DGMI, IG, I, U, &
                         hprogram,ki,                                   &
                        pcon_rain, pstrat_rain, pcon_snow, pstrat_snow,&
                        pclouds,   plsm,        pevaptr,   pevap,      &
                        pswec,     ptsc,        pucls,     pvcls,      &
                        pts,       pt2m,        phu2m,     pswe,       &
                        htest,     od_maskext,  plon_in,   plat_in )

!     ###############################################################################
!
!!****  *ASSIM_ISBA_n * - Chooses the surface assimilation schemes for ISBA
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     T. Aspelien
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    04/2012
!!      Trygve Aspelien, Separating IO  06/2013
!!--------------------------------------------------------------------
!
!
!
USE modd_diag_misc_isba_n, ONLY : diag_misc_isba_t
USE modd_isba_grid_n, ONLY : isba_grid_t
USE modd_isba_n, ONLY : isba_t
USE modd_surf_atm_n, ONLY : surf_atm_t
!
USE modd_surfex_mpi, ONLY : nrank, npio
USE modd_surf_par,       ONLY : xundef
USE modd_assim,          ONLY : cassim_isba,laesnm,lextrap_nature,nprintlev
!
!
USE yomhook,             ONLY : lhook,   dr_hook
USE parkind1,            ONLY : jprb
!
USE modi_abor1_sfx
USE modi_oi_hor_extrapol_surf
USE modi_assim_isba_update_snow
USE modi_assim_nature_isba_ekf
USE modi_assim_nature_isba_enkf
USE modi_assim_nature_isba_oi
USE modi_average_diag_misc_isba_n
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
TYPE(diag_misc_isba_t), INTENT(INOUT) :: dgmi
TYPE(isba_grid_t), INTENT(INOUT) :: ig
TYPE(isba_t), INTENT(INOUT) :: i
TYPE(surf_atm_t), INTENT(INOUT) :: u
!
 CHARACTER(LEN=6),    INTENT(IN) :: hprogram  ! program calling surf. schemes
INTEGER,             INTENT(IN) :: ki
REAL, DIMENSION(KI), INTENT(IN) :: pcon_rain
REAL, DIMENSION(KI), INTENT(IN) :: pstrat_rain
REAL, DIMENSION(KI), INTENT(IN) :: pcon_snow
REAL, DIMENSION(KI), INTENT(IN) :: pstrat_snow
REAL, DIMENSION(KI), INTENT(IN) :: pclouds
REAL, DIMENSION(KI), INTENT(IN) :: plsm
REAL, DIMENSION(KI), INTENT(IN) :: pevaptr
REAL, DIMENSION(KI), INTENT(IN) :: pevap
REAL, DIMENSION(KI), INTENT(IN) :: pswec
REAL, DIMENSION(KI), INTENT(IN) :: ptsc
REAL, DIMENSION(KI), INTENT(IN) :: pucls
REAL, DIMENSION(KI), INTENT(IN) :: pvcls
REAL, DIMENSION(KI), INTENT(IN) :: pts
REAL, DIMENSION(KI), INTENT(IN) :: pt2m
REAL, DIMENSION(KI), INTENT(IN) :: phu2m
REAL, DIMENSION(KI), INTENT(IN) :: pswe
 CHARACTER(LEN=2),    INTENT(IN) :: htest ! must be equal to 'OK'
LOGICAL,  DIMENSION (KI), INTENT(IN) ::  od_maskext
REAL(KIND=JPRB), DIMENSION (:), INTENT(IN) ::  plon_in
REAL(KIND=JPRB), DIMENSION (:), INTENT(IN) ::  plat_in
!
!*      0.2    declarations of local variables
!
!-------------------------------------------------------------------------------------
!
LOGICAL, DIMENSION(:), ALLOCATABLE :: ginterp_nature
LOGICAL, DIMENSION(:), ALLOCATABLE :: ginterp_sn
REAL,    DIMENSION(:), ALLOCATABLE :: zts_ep,zts_ep0
REAL,    DIMENSION(:), ALLOCATABLE :: ztp_ep,ztp_ep0
REAL,    DIMENSION(:), ALLOCATABLE :: zws_ep,zws_ep0
REAL,    DIMENSION(:), ALLOCATABLE :: zwp_ep,zwp_ep0
REAL,    DIMENSION(:), ALLOCATABLE :: ztl_ep,ztl_ep0
REAL,    DIMENSION(:), ALLOCATABLE :: zswe_ep,zswe_ep0
REAL,    DIMENSION(:), ALLOCATABLE :: zsnr_ep,zsnr_ep0
REAL,    DIMENSION(:), ALLOCATABLE :: zsna_ep,zsna_ep0
REAL,    DIMENSION(KI) :: zswe
REAL,    DIMENSION(KI) :: zswe_orig
INTEGER :: ji,jl,jp
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('ASSIM_ISBA_N',0,zhook_handle)
!
IF (htest/='OK') THEN
  CALL abor1_sfx('ASSIM_ISBA_n: FATAL ERROR DURING ARGUMENT TRANSFER')
END IF
!
! Set snow layers and patches
jp = 1
jl = 1
!
zswe = pswe
!
! General snow update
IF ( cassim_isba /= 'OI   ' ) THEN
  !
  ! Snow analysis/update
  IF (laesnm) THEN
    IF (nrank==npio) WRITE(*,*) 'UPDATE SNOW FROM ANALYSED VALUES'
    CALL assim_isba_update_snow(i, &
                                hprogram,ki,zswe,zswe_orig,.true.,.true.,htest)
  ELSE
    IF (nrank==npio) WRITE(*,*) 'SNOW IS NOT UPDATED FROM ANALYSED VALUES'
  ENDIF
ENDIF
!
! Soil assimilation
IF ( cassim_isba == 'EKF  ' ) THEN
  !
  ! Run EKF for soil
  CALL assim_nature_isba_ekf(i, &
                             hprogram, ki, pt2m, phu2m, htest)
  !
ELSEIF ( cassim_isba == 'ENKF ') THEN
  !
  CALL assim_nature_isba_enkf(i, hprogram, ki, pt2m, phu2m, htest)
  !  
ELSEIF ( cassim_isba == 'OI   ' ) THEN
  !
  ! Snow analysis/update. Store the original field in the surfex file
  IF (laesnm) THEN
    IF (nrank==npio) WRITE(*,*) 'UPDATE SNOW FROM ANALYSED VALUES'
    CALL assim_isba_update_snow(i, &
                                hprogram,ki,zswe,zswe_orig,.true.,.false.,htest)
  ELSE
    IF (nrank==npio) WRITE(*,*) 'SNOW IS NOT UPDATED FROM ANALYSED VALUES'
  ENDIF
  !
  ! Run OI for soil
  CALL assim_nature_isba_oi(i, &
                            hprogram, ki,                                  &
                            pcon_rain, pstrat_rain, pcon_snow, pstrat_snow,&
                            pclouds,   plsm,        pevaptr,   pevap,      &
                            pswec,     ptsc,        pucls,     pvcls,      &
                            pts,       pt2m,        phu2m,     zswe,       &
                            htest,     od_maskext,  plon_in,   plat_in )
  !
  ! Snow analysis/update (changed in oi_cacsts). Get the full increment
  IF (laesnm) THEN
    IF (nrank==npio) WRITE(*,*) 'UPDATE SNOW FROM ANALYSED OI_CACSTS VALUES'
    CALL assim_isba_update_snow(i, &
                                hprogram,ki,zswe,zswe_orig,.false.,.true.,htest)
  ELSE
    IF (nrank==npio) WRITE(*,*) 'SNOW IS NOT UPDATED FROM ANALYSED OI_CACSTS VALUES'
  ENDIF
  !
ELSE
  CALL abor1_sfx(cassim_isba//' is not a defined scheme for ASSIM_ISBA_N')
ENDIF

! Extrapolation if requested
IF ( lextrap_nature ) THEN
  !
  ALLOCATE(zws_ep(ki),zwp_ep(ki),zts_ep(ki),ztp_ep(ki),&
           ztl_ep(ki),zswe_ep(ki),zsnr_ep(ki),zsna_ep(ki))
  !  
  zws_ep  = i%XWG(:,1,jp)
  zwp_ep  = i%XWG(:,2,jp)
  zts_ep  = i%XTG(:,1,jp)
  ztp_ep  = i%XTG(:,2,jp)
  ztl_ep  = i%XWGI(:,2,jp)
  zswe_ep = i%TSNOW%WSNOW(:,jl,jp)
  zsnr_ep = i%TSNOW%RHO  (:,jl,jp)
  zsna_ep = i%TSNOW%ALB  (:,   jp)
  !
  ALLOCATE(ginterp_nature(ki),ginterp_sn(ki))
  !
  ! Search for the nearest grid point values for land surface fields
  ! at locations where the CANARI land fraction is less than 50%
  ! and therefore useless values MIGTH be given
  ginterp_nature = .false.
  ginterp_sn     = .false.
  !
  ! Snow albedo and density are also extrapolated in points 
  ! which get initial snow in the snow analysis
  WHERE ( zswe_ep(:) < 1.0e-10 .AND. pswe(:)>= 1.0e-10 )
    ginterp_sn(:) = .true.
    zsna_ep(:)    = xundef
    zsnr_ep(:)    = xundef
  END WHERE
  zswe_ep(:) = pswe(:)
  !
  WHERE ( plsm(:) < 0.5 )
    ginterp_nature(:) = .true.
    ginterp_sn(:) = .true.
    zts_ep(:)     = xundef
    ztp_ep(:)     = xundef
    zws_ep(:)     = xundef
    zwp_ep(:)     = xundef
    ztl_ep(:)     = xundef
    zswe_ep(:)    = xundef
    zsna_ep(:)    = xundef
    zsnr_ep(:)    = xundef
  END WHERE
  !
  ALLOCATE(zws_ep0(ki),zwp_ep0(ki),zts_ep0(ki),ztp_ep0(ki),&
           ztl_ep0(ki),zswe_ep0(ki),zsnr_ep0(ki),zsna_ep0(ki))
  !
  zws_ep0(:) = zws_ep(:)
  CALL oi_hor_extrapol_surf(ki,ig%XLAT,ig%XLON,zws_ep0,ig%XLAT,ig%XLON,zws_ep,ginterp_nature)
  zwp_ep0(:) = zwp_ep(:)
  CALL oi_hor_extrapol_surf(ki,ig%XLAT,ig%XLON,zwp_ep0,ig%XLAT,ig%XLON,zwp_ep,ginterp_nature)
  zts_ep0(:) = zts_ep(:)
  CALL oi_hor_extrapol_surf(ki,ig%XLAT,ig%XLON,zts_ep0,ig%XLAT,ig%XLON,zts_ep,ginterp_nature,i%XZS)
  ztp_ep0(:) = ztp_ep(:)
  CALL oi_hor_extrapol_surf(ki,ig%XLAT,ig%XLON,ztp_ep0,ig%XLAT,ig%XLON,ztp_ep,ginterp_nature,i%XZS)
  ztl_ep0(:) = ztl_ep(:)
  CALL oi_hor_extrapol_surf(ki,ig%XLAT,ig%XLON,ztl_ep0,ig%XLAT,ig%XLON,ztl_ep,ginterp_nature)
  zswe_ep0(:) = zswe_ep(:)
  CALL oi_hor_extrapol_surf(ki,ig%XLAT,ig%XLON,zswe_ep0,ig%XLAT,ig%XLON,zswe_ep,ginterp_sn)
  zsnr_ep0(:) = zsnr_ep(:)
  CALL oi_hor_extrapol_surf(ki,ig%XLAT,ig%XLON,zsnr_ep0,ig%XLAT,ig%XLON,zsnr_ep,ginterp_sn)
  zsna_ep0(:) = zsna_ep(:)
  CALL oi_hor_extrapol_surf(ki,ig%XLAT,ig%XLON,zsna_ep0,ig%XLAT,ig%XLON,zsna_ep,ginterp_sn)  
  !
  DEALLOCATE(zws_ep0,zwp_ep0,zts_ep0,ztp_ep0,ztl_ep0,zswe_ep0,zsnr_ep0,zsna_ep0)
  !
  ! PRINT values produced by OI_HO_EXTRAPOL_SURF for TS
  IF ( nprintlev > 2 ) THEN
    DO ji=1,ki
     IF (ginterp_nature(ji)) THEN
       print *,'Surface temperature set to ',zts_ep(ji),'from nearest neighbour at I=',u%NR_NATURE(ji)
     ENDIF
    ENDDO
  ENDIF
  !
  DEALLOCATE(ginterp_nature,ginterp_sn)
  !
  ! Set extrpolated fields to global
  i%XWG (:,1,jp) = zws_ep(:)
  i%XWG (:,2,jp) = zwp_ep(:)
  i%XTG (:,1,jp) = zts_ep(:)
  i%XTG (:,2,jp) = ztp_ep(:)
  i%XWGI(:,2,jp) = ztl_ep(:)
  i%TSNOW%WSNOW(:,jl,jp) = zswe_ep(:)
  i%TSNOW%RHO  (:,jl,jp) = zsnr_ep(:)
  i%TSNOW%ALB  (:,   jp) = zsna_ep(:)
  !
  DEALLOCATE(zws_ep,zwp_ep,zts_ep,ztp_ep,ztl_ep,zswe_ep,zsnr_ep,zsna_ep)
  !
ENDIF

! Snow analysis/update security
IF (laesnm) THEN

  ! removes very small values due to computation precision
  WHERE( i%TSNOW%WSNOW(:,jl,jp) < 1.0e-10 ) i%TSNOW%WSNOW(:,jl,jp) = 0.0

  ! No SNOW
  WHERE ( i%TSNOW%WSNOW(:,jl,jp) == 0.0 )
    i%TSNOW%RHO(:,jl,jp) = xundef
    i%TSNOW%ALB(:,jp)    = xundef
  END WHERE
  !
ENDIF
!
!to be improved later - needed for surfex course
 CALL average_diag_misc_isba_n(dgmi, i)
 !
IF (lhook) CALL dr_hook('ASSIM_ISBA_N',1,zhook_handle)
!
!-------------------------------------------------------------------------------------
!
END SUBROUTINE assim_isba_n