restart_coupl_topd.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 restart_coupl_topd (UG, U, &
                               hprogram,ki)
!###################################################################
!
!!****  *RESTART_COUPL_TOPDn*  
!!
!!    PURPOSE
!!    -------
!!  Read all files needed in case of restart 
!!      
!!    REFERENCE
!!    ---------
!!      
!!    AUTHOR
!!    ------
!!      B. Vincendon
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    07/06/11 
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!
!
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
!
USE modd_surf_par,        ONLY : xundef,nundef
USE modd_isba_par,        ONLY : xwgmin
!
USE modd_topd_par, ONLY : nunit
USE modd_topodyn,       ONLY : nncat, ccat, nnpt, nline, nnmc, npmax,&
                               nnb_topd_step

USE modd_coupling_topd, ONLY : xas_nature, &
                                 nnb_stp_stock,nnb_stp_restart,xwtopt,&
                                 xrun_torout,xdr_torout
!
USE modi_read_topd_file
USE modi_read_file_isbamap
!
USE modi_get_luout
USE modi_abor1_sfx
USE modi_unpack_same_rank
USE modi_pack_same_rank
!
USE modi_open_file
USE modi_close_file
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
TYPE(surf_atm_grid_t), INTENT(INOUT) :: ug
TYPE(surf_atm_t), INTENT(INOUT) :: u
!
 CHARACTER(LEN=6), INTENT(IN)         :: hprogram ! program calling surf. schemes
INTEGER,          INTENT(IN)         :: ki       ! Surfex grid dimension
!
!
!*      0.2    declarations of local variables
!
INTEGER :: jj
INTEGER                                     :: iluout   ! unit of output listing file
INTEGER                                     :: jstp,jcat,jpix! loop control indexes
REAL, DIMENSION(:),ALLOCATABLE              :: zas      ! Saturated area fraction for each Isba meshes
REAL, DIMENSION(:),ALLOCATABLE              :: zwtopt   ! Initial water content in case of restart
 CHARACTER(LEN=50), DIMENSION(:),ALLOCATABLE :: yfiletop ! File names
LOGICAL                                     :: lstock, lwg, lasat
REAL                                        :: zcorr_stock ! used to avoid to lose stock
REAL                                        :: zcnt_undef,zsum1,zsum2, zdenom ! used to correct budget
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('RESTART_COUPL_TOPD',0,zhook_handle)
!
 CALL get_luout(hprogram,iluout)
!
! * 1. Read stock files
!       
WRITE(*,*) 'Read STOCK file ',nnb_stp_stock
nnb_stp_stock = min(nnb_stp_stock, nnb_topd_step + nnb_stp_restart)
!
INQUIRE(file='stocks_init.txt', exist=lstock)
INQUIRE(file='surfcont_init.map', exist=lasat)
!
IF (.NOT.lstock) THEN
  WRITE(iluout,*) 'You asked to run in restart mode but stock file is missing'
  CALL abor1_sfx("RESTART_COUPL_TOPD_n: stock file is missing")
ELSEIF (.NOT.lasat) THEN
  WRITE(iluout,*) 'You asked to run in restart mode but contributive area file is missing'
  CALL abor1_sfx("RESTART_COUPL_TOPD_n: contributive area file is missing")
ELSE
  !
  CALL open_file('ASCII ',nunit,'stocks_init.txt','FORMATTED',haction='READ ')
  DO jstp=1,nnb_stp_stock
    READ(nunit,*)  xrun_torout(1:nncat,jstp),xdr_torout(1:nncat,jstp)
  ENDDO
  CALL close_file('ASCII ',nunit)
  !  
  ! * 2. Read pixels water content
  !
  DO jcat=1,nncat
    !
    IF (xrun_torout(jcat,nnb_stp_stock)/=0.) THEN
      !
      zcorr_stock = xrun_torout(jcat,nnb_stp_stock) - xrun_torout(jcat,nnb_stp_stock-1)
      DO jstp = nnb_stp_stock+1,nnb_topd_step
        xrun_torout(jcat,jstp) = max(0.,xrun_torout(jcat,jstp-1)+zcorr_stock)
      ENDDO
      !
    ENDIF
    !
    IF (xdr_torout(jcat,nnb_stp_stock)/=0.) THEN
      !
      zcorr_stock = xdr_torout(jcat,nnb_stp_stock) - xdr_torout(jcat,nnb_stp_stock-1)
      DO jstp = nnb_stp_stock+1,nnb_topd_step
        xdr_torout(jcat,jstp) = max(0.,xdr_torout(jcat,jstp-1)+zcorr_stock)
      ENDDO
      !
    ENDIF
    !
  ENDDO
  !
  WRITE(*,*) 'Write pixels water content files'
  !
  ALLOCATE(zwtopt(npmax))
  ALLOCATE(yfiletop(nncat))
  !
  DO jcat=1,nncat
    !
    yfiletop(jcat)=trim(ccat(jcat))//'_xwtop_init.map'
    INQUIRE(file=yfiletop(jcat), exist=lwg)
    IF (.NOT.lwg) THEN
      !
      WRITE(iluout,*) 'You asked to run in restart mode but pixels water content file is missing'
      WRITE(iluout,*) 'for catchment : ',ccat(jcat)
      CALL abor1_sfx("RESTART_COUPL_TOPD_n: pixels water content file is missing")
      !
    ELSE
      !
      zsum1=sum(xwtopt(jcat,:),mask=xwtopt(jcat,:)/=xundef)
      !
      CALL read_topd_file('ASCII ',yfiletop(jcat),'FORMATTED',nnpt(jcat),zwtopt)
      !
      DO jpix=1,SIZE(nline(jcat,:))
        IF ( nline(jcat,jpix)/=0 .AND. nline(jcat,jpix)/=xundef ) THEN
          IF (zwtopt(jpix) /= xundef) xwtopt(jcat,nline(jcat,jpix)) = zwtopt(jpix)
        ENDIF
      ENDDO
      !
      zsum2=sum(xwtopt(jcat,:),mask=xwtopt(jcat,:)<xundef)
      !
      IF ( abs(zsum2-zsum1)>100. ) THEN
        !
        zcnt_undef = 0.
        DO jj=1,SIZE(nline,2)
          IF ( nline(jcat,jj)/=0 ) THEN
            IF ( xwtopt(jcat,nline(jcat,jj))/=xundef ) zcnt_undef = zcnt_undef + 1.
          ENDIF
        ENDDO
        !
        IF (zcnt_undef/=0.) THEN
          zdenom = (zsum2-zsum1)/zcnt_undef
          DO jj=1,SIZE(nline,2)
            IF ( nline(jcat,jj)/=0 ) THEN
              IF ( xwtopt(jcat,nline(jcat,jj))/=xundef ) THEN
                xwtopt(jcat,nline(jcat,jj)) = xwtopt(jcat,nline(jcat,jj)) - zdenom
              ENDIF
            ENDIF
          ENDDO
        ENDIF
        zsum2=sum(xwtopt(jcat,:),mask=xwtopt(jcat,:)<xundef)
        !
      ENDIF
      !
    ENDIF
    !
  ENDDO
  ! 
  ! * 3. Read Asat files
  ! 
  WRITE(*,*) 'Read Asat files'
  ALLOCATE(zas(ki))
  CALL open_file('ASCII ',nunit,'surfcont_init.map','FORMATTED',haction='READ ')
  CALL read_file_isbamap(ug, &
                         nunit,zas,ki)
  CALL close_file('ASCII ',nunit)
  CALL pack_same_rank(u%NR_NATURE,zas,xas_nature)
  !
ENDIF
!
IF (lhook) CALL dr_hook('RESTART_COUPL_TOPD',1,zhook_handle)
!-------------------------------------------------------------------------------
END SUBROUTINE restart_coupl_topd