build_emisstabn.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 build_emisstab_n (PCONVERSION, HPROGRAM,KCH,HEMIS_GR_NAME, KNBTIMES,&
              KEMIS_GR_TIME,KOFFNDX,TPEMISS,KSIZE,KLUOUT, KVERB,PRHODREF)  
!!    #####################################################################
!!
!!*** *BUILD_EMISSTAB*
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    AUTHOR
!!    ------
!!    D. Gazen
!!
!!    MODIFICATIONS
!!    -------------
!!    Original 01/02/00
!!    C. Mari  30/10/00  call of MODD_TYPE_EFUTIL and MODD_CST
!!    D.Gazen  01/12/03  change emissions handling for surf. externalization!!
!!    P.Tulet  01/01/04  change conversion for externalization (flux unit is
!!                        molec./m2/s)
!!    M.Leriche  04/14   apply conversion factor if lead = f
!!                       change emissions name EMIS_ -> E_ name for coherence with PGD
!!    M.Moge    01/2016  using READ_SURF_FIELD2D for 2D surfex fields reads
!!
!!    EXTERNAL
!!    --------
!
USE modi_ch_open_inputb
USE modi_read_surf_field2d
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
USE modd_type_efutil, ONLY : emissvar_t
USE modd_csts,        ONLY : ndaysec, xmd, xavogadro
!------------------------------------------------------------------------------
!
!*       0.   DECLARATIONS
!        -----------------
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
!*       0.1  declaration of arguments
!
!
!
REAL, DIMENSION(:), POINTER :: PCONVERSION
!
 CHARACTER(LEN=6),                INTENT(IN) :: HPROGRAM   ! Program name
INTEGER,                         INTENT(IN) :: KCH
 CHARACTER(LEN=*),DIMENSION(:),   INTENT(IN) :: HEMIS_GR_NAME ! Offline species name
INTEGER, DIMENSION(:),           INTENT(IN) :: KNBTIMES ! nb of emis times array
INTEGER, DIMENSION(:),           INTENT(IN) :: KEMIS_GR_TIME
INTEGER, DIMENSION(:),           INTENT(IN) :: KOFFNDX ! index of offline species
TYPE(emissvar_t),DIMENSION(:),   INTENT(OUT):: TPEMISS ! emission struct array to fill
INTEGER,                         INTENT(IN) :: KSIZE   ! size X*Y (1D) of physical domain
INTEGER,                         INTENT(IN) :: KLUOUT  ! output listing channel
INTEGER,                         INTENT(IN) :: KVERB   ! verbose level
REAL, DIMENSION(:),              INTENT(IN) :: PRHODREF ! dry density for ref. state
!
!
!*       0.2  declaration of local variables
!
 CHARACTER(LEN=3):: YUNIT       ! unit of the flux
INTEGER         :: INBTS       ! Number of emis times for a species
INTEGER         :: IRESP       ! I/O return value
INTEGER         :: IIND1, IIND2
INTEGER         :: JSPEC       ! loop index
INTEGER         :: ITIME       ! loop index
INTEGER         :: IWS_DEFAULT ! Default Memory window size for emission reading
 CHARACTER (LEN=16):: YRECFM    ! LFI article name
REAL(KIND=JPRB) :: ZHOOK_HANDLE

!
!------------------------------------------------------------------------------
!
!*    EXECUTABLE STATEMENTS
!     ---------------------
!

IF (lhook) CALL dr_hook('BUILD_EMISSTAB_N',0,zhook_handle)
IF (kverb >= 5) THEN
  WRITE(kluout,*) '********     SUBROUTINE (CHIMIE): BUILD_EMISSTAB_n     ********'
END IF
!
!*       1.   READ DATA 
!        --------------
!
 CALL ch_open_inputb("EMISUNIT", kch, kluout)
!
! read unit identifier
READ(kch,'(A3)') yunit
!
!*       2.   MAP DATA ONTO PROGNOSTIC VARIABLES
!        ---------------------------------------
!
ALLOCATE (pconversion(SIZE(prhodref,1)))
! determine the conversion factor
  pconversion(:) = 1.
SELECT CASE (yunit)
CASE ('MIX') ! flux given ppp*m/s,  conversion to molec/m2/s
! where 1 molecule/cm2/s = (224.14/6.022136E23) ppp*m/s
  pconversion(:) = xavogadro * prhodref(:) / xmd
CASE ('CON') ! flux given in molecules/cm2/s, conversion to molec/m2/s 
  pconversion(:) =  1e4
CASE ('MOL') ! flux given in microMol/m2/day, conversion to molec/m2/s  
! where 1 microMol/m2/day = (22.414/86.400)*1E-12 ppp*m/s
  !XCONVERSION(:) = (22.414/86.400)*1E-12 * XAVOGADRO * PRHODREF(:) / XMD
  pconversion(:) = 1e-6 * xavogadro / 86400.

CASE DEFAULT
  CALL abor1_sfx('CH_BUILDEMISSN: UNKNOWN CONVERSION FACTOR')
END SELECT
!
! Read Window size default value >= 2
iws_default = 5 ! Should be set by namelist
IF (iws_default < 2) iws_default = 2
!
iind1 = 0
iind2 = 0
DO jspec=1,SIZE(tpemiss) ! loop on offline emission species
!
  inbts = knbtimes(jspec)
!
! Fill %CNAME
  tpemiss(jspec)%CNAME = hemis_gr_name(koffndx(jspec))
! Allocate and Fill %NETIMES 
  ALLOCATE(tpemiss(jspec)%NETIMES(inbts))
  iind1 = iind2+1
  iind2 = iind2+inbts
  tpemiss(jspec)%NETIMES(:) = kemis_gr_time(iind1:iind2)
! 
! Update %NWS, %NDX, %NTX, %LREAD, %XEMISDATA
  IF (inbts <= iws_default) THEN
! Number of times smaller than read window size allowed
! Read emis data once and for all
    tpemiss(jspec)%NWS = inbts
    tpemiss(jspec)%NDX = 1
    tpemiss(jspec)%NTX = 1
    tpemiss(jspec)%LREAD = .false. ! to prevent future reading
    ALLOCATE(tpemiss(jspec)%XEMISDATA(ksize,inbts))
! Read file for emission data
    yrecfm='E_'//trim(tpemiss(jspec)%CNAME)
    CALL read_surf_field2d(hprogram,tpemiss(jspec)%XEMISDATA(:,:),yrecfm)
!
! Correction : Replace 999. with 0. value in the Emission FLUX
! and apply conversion
    WHERE(tpemiss(jspec)%XEMISDATA(:,:) == 999.)
      tpemiss(jspec)%XEMISDATA(:,:) = 0. 
    END WHERE
    WHERE(tpemiss(jspec)%XEMISDATA(:,:) == 1.e20)
      tpemiss(jspec)%XEMISDATA(:,:) = 0. 
    END WHERE
      DO itime=1,inbts
      ! XCONVERSION HAS BEEN ALREADY APPLY IN CH_EMISSION_FLUXN ONLY FOR LREAD = T
      tpemiss(jspec)%XEMISDATA(:,itime) = tpemiss(jspec)%XEMISDATA(:,itime) * pconversion(:)
      !TPEMISS(JSPEC)%XEMISDATA(:,ITIME) = TPEMISS(JSPEC)%XEMISDATA(:,ITIME)
      END DO
    ELSE
! Read window size is smaller than number of emission times
    tpemiss(jspec)%NWS = iws_default
    tpemiss(jspec)%NDX = iws_default
    tpemiss(jspec)%NTX = 0
    tpemiss(jspec)%LREAD = .true.
    ALLOCATE(tpemiss(jspec)%XEMISDATA(ksize,iws_default))
  END IF
 
  IF (inbts == 1) THEN
    tpemiss(jspec)%XFWORK=>tpemiss(jspec)%XEMISDATA(:,1)
  ELSE
    ALLOCATE(tpemiss(jspec)%XFWORK(ksize))
  END IF
! Compute index for periodic case
  tpemiss(jspec)%NPX = maxval(minloc(tpemiss(jspec)%NETIMES(:)+&
         (1+(tpemiss(jspec)%NETIMES(inbts)-&
         tpemiss(jspec)%NETIMES(:))/ndaysec)*ndaysec))  
!
! Some di###ay
  IF (kverb >= 6) THEN
    WRITE(kluout,*) '====== Species ',trim(tpemiss(jspec)%CNAME), ' ======'
    WRITE(kluout,*) '  Emission Times :' ,tpemiss(jspec)%NETIMES
    WRITE(kluout,*) '  Current time index :' ,tpemiss(jspec)%NTX
    WRITE(kluout,*) '  Current data index :' ,tpemiss(jspec)%NDX
    WRITE(kluout,*) '  Periodic index = ',tpemiss(jspec)%NPX,&
            ' at time :',tpemiss(jspec)%NETIMES(tpemiss(jspec)%NPX)  
    WRITE(kluout,*) '  Read window size :', tpemiss(jspec)%NWS
    IF (tpemiss(jspec)%LREAD) THEN
      WRITE(kluout,*) ' -> Data must be read during simulation.'
    ELSE
      WRITE(kluout,*) ' -> Data already in memory.'
    END IF
  END IF
END DO

IF (kverb >= 5) THEN
  WRITE(kluout,*) '******** END SUBROUTINE (CHIMIE) : BUILD_EMISSTAB_n     ********'
END IF
IF (lhook) CALL dr_hook('BUILD_EMISSTAB_N',1,zhook_handle)

END SUBROUTINE build_emisstab_n