ch_emission_fluxn.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 ch_emission_flux_n (YSC, &
                                     hprogram,psimtime,psfsv, prhoa, ptstep, knbts_max)
!     ######################################################################
!!
!!***  *CH_EMISSION_FLUX_n* - 
!!
!!    PURPOSE
!!    -------
!!      Return a time-dependent emission flux based on tabulated values
!!
!!**  METHOD
!!    ------
!!
!!    AUTHOR
!!    ------
!!    D. Gazen
!!
!!    MODIFICATIONS
!!    -------------
!!    Original 08/02/00
!!    C. Mari  30/10/00 call to MODD_TYPE_EFUTIL and MODD_CST
!!    D.Gazen  01/12/03  change emissions handling for surf. externalization
!!    P.Tulet  01/01/04  change emission conversion factor
!!    P.Tulet  01/01/05  add dust, orilam
!!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!
!
USE modd_surfex_n, ONLY : surfex_t
!
USE modd_type_efutil,      ONLY: emissvar_t, pronosvar_t
USE modd_csts,             ONLY: ndaysec
!
USE modi_read_surf
USE modi_init_io_surf_n
USE modi_end_io_surf_n
USE modi_get_luout
!UPG*AERO1
USE modd_chs_aerosol, ONLY: lch_aero_flux
USE modi_ch_aer_emission
!UPG*AERO1
!!
!------------------------------------------------------------------------------
!
!*       0.   DECLARATIONS
!        -----------------
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
!*       0.1  declaration of arguments
!
TYPE(surfex_t), INTENT(INOUT) :: ysc
!
 CHARACTER(LEN=6),   INTENT(IN)  :: hprogram    ! program calling surf. schemes
REAL,               INTENT(IN)  :: psimtime    ! time of simulation in sec UTC
                                               ! (counting from midnight of
                                               ! the current day)
REAL,DIMENSION(:,:),  INTENT(INOUT) :: psfsv   ! emission flux in ppp*m/s
REAL, DIMENSION(:),   INTENT(IN)    :: prhoa     ! air density (kg/m3)
REAL,                 INTENT(IN)    :: ptstep    ! atmospheric time-step                 (s)
INTEGER,              INTENT(IN)    :: knbts_max !max size of TEMISS%NETIMES

!
!*       0.2  declaration of local variables
!
INTEGER       :: iverb   ! verbosity level
INTEGER       :: ksize1d ! 1D size = X*Y physical domain 
INTEGER       :: ji      ! loop control
REAL          :: zalpha  ! interpolation weight
!
INTEGER :: inbts       ! Number of emission times for a species
INTEGER :: itim1,itim2 ! first/last time for interpolation
INTEGER :: indx1,indx2 ! first/next index for data interpolation
INTEGER :: isimtime, itperiod
 CHARACTER (LEN=16)  :: yrecfm          ! LFI article name
TYPE(pronosvar_t),POINTER :: curpronos !Current pronostic variable
!
!*       0.3  declaration of saved local variables
!
 CHARACTER(LEN=6), DIMENSION(:), POINTER :: cnames
REAL,DIMENSION(SIZE(PSFSV,1),KNBTS_MAX)     :: zwork ! temporary array for reading data
REAL,DIMENSION(SIZE(PSFSV,1),SIZE(PSFSV,2)) :: zemis ! interpolated in time emission flux
REAL,DIMENSION(SIZE(PSFSV,1),SIZE(PSFSV,2)) :: zdepot! interpolated in time deposition flux
REAL,DIMENSION(SIZE(PSFSV,1))               :: zfco  ! CO flux
INTEGER                          :: ineq  ! number of chemical var
                                          !(=NEQ (chimie gaz) + NSV_AER (chimie aerosol)
INTEGER                          :: iws   ! window size
INTEGER                          :: iresp ! return code for I/O
INTEGER                          :: iluout ! Outputlisting unit
LOGICAL                          :: lioinit ! True if I/O init done
INTEGER                          :: jw
INTEGER                          :: itime
LOGICAL                          :: gco = .false. ! switch if CO emission are available
REAL(KIND=JPRB) :: zhook_handle
!
!------------------------------------------------------------------------------
!
!*    EXECUTABLE STATEMENTS
!     ---------------------
!
IF (lhook) CALL dr_hook('CH_EMISSION_FLUX_N',0,zhook_handle)
 CALL get_luout(hprogram,iluout)
lioinit = .false.
iverb   = 5
ksize1d = SIZE(psfsv,1)
ineq    = SIZE(psfsv,2)
!
!------------------------------------------------------------------------------
!
!*    3.  INTERPOLATE SURFACE FLUXES IN TIME IF NEEDED
!     ------------------------------------------------
!
IF (ysc%CHE%XTIME_SIMUL == 0.) THEN
   ysc%CHE%XTIME_SIMUL = psimtime
ELSE
   ysc%CHE%XTIME_SIMUL = ysc%CHE%XTIME_SIMUL + ptstep
END IF

IF (iverb >= 5) WRITE(iluout,*) '******** CH_EMISSION_FLUX  ********'
DO ji=1,SIZE(ysc%CHE%TSEMISS)
! Simulation time (counting from midnight) is saved
  isimtime = ysc%CHE%XTIME_SIMUL
!
  inbts = SIZE(ysc%CHE%TSEMISS(ji)%NETIMES) ! 
  iws   = ysc%CHE%TSEMISS(ji)%NWS           ! Window Size for I/O
  indx1 = ysc%CHE%TSEMISS(ji)%NDX           ! Current data index
!
  IF (inbts == 1) THEN
!   Time Constant Flux
!   XFWORK already points on data (see ch_buildemiss.f90)
    IF (iverb >= 6) THEN
      WRITE(iluout,*) 'NO interpolation for ',trim(ysc%CHE%TSEMISS(ji)%CNAME)
      IF (iverb >= 10 ) WRITE(iluout,*) ysc%CHE%TSEMISS(ji)%XFWORK
    END IF
  ELSE
    IF (iverb >= 6) THEN
      WRITE(iluout,*) 'Interpolation (T =',isimtime,') : ',ysc%CHE%TSEMISS(ji)%CNAME
    END IF
    IF (isimtime < ysc%CHE%TSEMISS(ji)%NETIMES(1)) THEN
!     Tsim < T(1)=Tmin should not happen but who knows ?
      ysc%CHE%TSEMISS(ji)%NTX = 1
    ELSE
!     Check for periodicity when ISIMTIME is beyond last emission time
!     and probably correct ISIMTIME
      IF (isimtime > ysc%CHE%TSEMISS(ji)%NETIMES(inbts)) THEN 
!       Tsim > T(INBTS)=Tmax
        itperiod = (1+(ysc%CHE%TSEMISS(ji)%NETIMES(inbts)-&
                ysc%CHE%TSEMISS(ji)%NETIMES(ysc%CHE%TSEMISS(ji)%NPX))/ndaysec)*ndaysec  
        isimtime = modulo(isimtime-ysc%CHE%TSEMISS(ji)%NETIMES(ysc%CHE%TSEMISS(ji)%NPX),itperiod)+&
                ysc%CHE%TSEMISS(ji)%NETIMES(ysc%CHE%TSEMISS(ji)%NPX)  
        IF (iverb >= 6) THEN
          WRITE(iluout,*) '  ITPERIOD = ', itperiod
          WRITE(iluout,*) '  ISIMTIME modifie = ', isimtime
        END IF
        IF (ysc%CHE%TSEMISS(ji)%NTX == inbts .AND. isimtime<ysc%CHE%TSEMISS(ji)%NETIMES(inbts)) THEN
!         Update time index NTX 
          ysc%CHE%TSEMISS(ji)%NTX = ysc%CHE%TSEMISS(ji)%NPX
!         Increment data index NDX : NDX correction will occur later
!                                    to assure 1 <= NDX <= IWS
          indx1 = indx1 + 1
        END IF
      END IF
!
!     search NTX such that : ETIMES(NTX) < ISIMTIME <= ETIMES(NTX+1)
!     and make NDX follow NTX : NDX correction will occur later
!                               to assure 1 <= NDX <= IWS
      DO WHILE (ysc%CHE%TSEMISS(ji)%NTX < inbts)
        IF (isimtime >= ysc%CHE%TSEMISS(ji)%NETIMES(ysc%CHE%TSEMISS(ji)%NTX+1)) THEN
          ysc%CHE%TSEMISS(ji)%NTX = ysc%CHE%TSEMISS(ji)%NTX + 1
          indx1 = indx1 + 1
          indx2 = indx1 + 1
        ELSE
          EXIT
        END IF
      END DO
    END IF
!
!   Check availability of data within memory Window (XEMISDATA(:,1:IWS))
    IF (indx1 >= iws) THEN
!
!     Data index reached the memory window limits
!
      IF (ysc%CHE%TSEMISS(ji)%LREAD) THEN 
!
!       File must be read to update XEMISDATA array for this species 
!
        IF (.NOT. lioinit) THEN
!         Must be done once before reading
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        hprogram,'FULL  ','SURF  ','READ ')
          IF (iverb >= 6) WRITE(iluout,*) 'INIT des I/O DONE.'
          lioinit=.true.
        END IF
        yrecfm='E_'//trim(ysc%CHE%TSEMISS(ji)%CNAME)
        IF (iverb >= 6)&
               WRITE (iluout,*) 'READ emission :',trim(yrecfm),&
               ', SIZE(ZWORK)=',SIZE(zwork,1),inbts 
        CALL read_surf(hprogram,yrecfm,zwork(:,1:inbts),iresp)
!
! Correction : Replace 999. with 0. value in the Emission FLUX
        WHERE(zwork(:,1:inbts) == 999.)
          zwork(:,1:inbts) = 0. 
        END WHERE
        WHERE(zwork(:,1:inbts) == 1.e20)
          zwork(:,1:inbts) = 0. 
        END WHERE
        DO itime=1,inbts
        zwork(:,itime) = zwork(:,itime)*ysc%CHU%XCONVERSION(:)
        END DO
!
!       
        IF ((ysc%CHE%TSEMISS(ji)%NTX+iws-1) > inbts) THEN
!
!         ===== Periodic CASE =====
!
          IF (iverb >= 6)&
                 WRITE (iluout,*) 'Periodic CASE : NPX =',ysc%CHE%TSEMISS(ji)%NPX  
          IF (iws <  (inbts-ysc%CHE%TSEMISS(ji)%NPX+1)) THEN
!           Window size is smaller then number of periodical times
!
!           example : IWS=5, NPX=2, INBTS=11, NTX=9
!                               NTX       NPX
!                                |         |
!           time index :      ...9 10 11 # 2 3 4...11 #
!       old data index :[1 2 3 4 5] 
!       new data index :        [1  2  3   4 5]
!                                |  
!                               NDX    
!
            ysc%CHE%TSEMISS(ji)%XEMISDATA(:,1:inbts-ysc%CHE%TSEMISS(ji)%NTX+1) = &
                   zwork(:,ysc%CHE%TSEMISS(ji)%NTX:inbts)  
!
            IF (iverb >= 6) THEN
              WRITE(iluout,*) 'Window SIZE smaller than INBTS !'
              WRITE(iluout,*) 'Window index, Time index'
              DO jw=1,inbts-ysc%CHE%TSEMISS(ji)%NTX+1
                WRITE(iluout,*) jw,ysc%CHE%TSEMISS(ji)%NTX+jw-1
              END DO
            END IF
!
            ysc%CHE%TSEMISS(ji)%XEMISDATA(:,inbts-ysc%CHE%TSEMISS(ji)%NTX+2:iws) = &
                   zwork(:,ysc%CHE%TSEMISS(ji)%NPX:ysc%CHE%TSEMISS(ji)%NPX+iws-inbts+ysc%CHE%TSEMISS(ji)%NTX-2)  
!
            IF (iverb >= 6) THEN
              DO jw=inbts-ysc%CHE%TSEMISS(ji)%NTX+2,iws
                WRITE(iluout,*) jw,ysc%CHE%TSEMISS(ji)%NPX+jw-(inbts-ysc%CHE%TSEMISS(ji)%NTX+2)
              END DO
            END IF
            indx1 = 1
            indx2 = 2
          ELSE
!           Window size may get smaller AND it will be the last reading
!
!           example : IWS=6, NPX=7, INBTS=11, NTX=9
!
!                         NTX       NPX NTX
!                          |         |   |
!           time index: ...9 10 11 # 7 8 9 10 11 #
!       old data index: ...6]
!       new data index:             [1 2 3  4  5]
!                                        |
!                                       NDX=NTX-NPX+1
!
            iws = inbts-ysc%CHE%TSEMISS(ji)%NPX+1
            ysc%CHE%TSEMISS(ji)%NWS = iws
            ysc%CHE%TSEMISS(ji)%XEMISDATA(:,1:iws) = zwork(:,ysc%CHE%TSEMISS(ji)%NPX:inbts)
            IF (iverb >= 6) THEN
              WRITE(iluout,*) 'Window SIZE equal or greater than INBTS !'
              WRITE(iluout,*) 'Window index, Time index'
              DO jw=1,iws
                WRITE(iluout,*) jw,ysc%CHE%TSEMISS(ji)%NPX+jw-1
              END DO
            END IF
            indx1 = ysc%CHE%TSEMISS(ji)%NTX-ysc%CHE%TSEMISS(ji)%NPX+1
            indx2 = mod((indx1+1),iws)
            ysc%CHE%TSEMISS(ji)%LREAD = .false. ! no more reading
          END IF
        ELSE
!
!         ===== NON periodic (normal) CASE =====
!
! example : with IWS=5, the window moves forward
!                             NTX
!                              | 
!         time index : 1 2 3 4 5 6 7 8 9 10 11 ... INBTS # 
!     old data index :[1 2 3 4 5] 
!     new data index :        [1 2 3 4 5] 
!                              |
!                             NDX
!
          ysc%CHE%TSEMISS(ji)%XEMISDATA(:,1:iws) = zwork(:,ysc%CHE%TSEMISS(ji)%NTX:ysc%CHE%TSEMISS(ji)%NTX+iws-1)
          IF (iverb >= 6) THEN
            WRITE(iluout,*) 'Window index, Time index'
            DO jw=1,iws
              WRITE(iluout,*) jw,ysc%CHE%TSEMISS(ji)%NTX+jw-1
            END DO
          END IF
          indx1 = 1
          indx2 = 2
        END IF
      ELSE
!       Data is already in memory because window size is sufficient 
!       to hold INBTS emission times => simply update NDX according to NTX
!       
        IF (iws==inbts) THEN 
!
!         'Window size' = 'Nb emis times' at INIT (ch_init_emission)
!         so NDX must be set equal to NTX (the window does not move)
! example :
!                         NPX    NTX
!                          |      | 
!         time index :  1  2  3  ... INBTS
!         data index : [1  2  3  ... INBTS]
!                                 |
!                                NDX

          indx1 = ysc%CHE%TSEMISS(ji)%NTX
          indx2 = indx1+1
          IF (indx2 > iws) indx2=ysc%CHE%TSEMISS(ji)%NPX
        ELSE
!          
!         Windows size changed during periodic case
!         NDX must be equal to NTX - NPX + 1
!         (the window does not move)
! example :
!                                NTX
!                                 | 
!         time index : NPX NPX+1 NPX+2 ... INBTS
!         data index : [1    2    3    ...   IWS]
!                                 |
!                                NDX
          indx1 = ysc%CHE%TSEMISS(ji)%NTX-ysc%CHE%TSEMISS(ji)%NPX+1
          indx2 = mod((indx1+1),iws)
        END IF
      END IF
    ELSE ! (INDX1 < IWS)
      indx2 = indx1+1
    END IF
!
!   Don't forget to update NDX with new value INDX1
    ysc%CHE%TSEMISS(ji)%NDX = indx1
!
!   Compute both times for interpolation
    IF (ysc%CHE%TSEMISS(ji)%NTX < inbts) THEN 
      itim1 = ysc%CHE%TSEMISS(ji)%NETIMES(ysc%CHE%TSEMISS(ji)%NTX)
      itim2 = ysc%CHE%TSEMISS(ji)%NETIMES(ysc%CHE%TSEMISS(ji)%NTX+1)
    ELSE
      itim1 = ysc%CHE%TSEMISS(ji)%NETIMES(inbts)
      itim2 = ysc%CHE%TSEMISS(ji)%NETIMES(ysc%CHE%TSEMISS(ji)%NPX)+itperiod
    END IF
!
! Interpolate variables in time -> update XFWORK
!
!
!  time       :  ITIM1...Tsim...ITIM2
!                  |              |  
!  data index :  INDX1          INDX2
!
!
    zalpha = (REAL(ISIMTIME) - itim1) / (itim2-itim1)
    ysc%CHE%TSEMISS(ji)%XFWORK(:) = zalpha*ysc%CHE%TSEMISS(ji)%XEMISDATA(:,indx2) +&
            (1.-zalpha)*ysc%CHE%TSEMISS(ji)%XEMISDATA(:,indx1)  
    IF (iverb >= 6) THEN
      WRITE(iluout,*) '  Current time INDEX : ',ysc%CHE%TSEMISS(ji)%NTX
      WRITE(iluout,*) '  TIME : ',isimtime, ' (',itim1,',',itim2,')'
      WRITE(iluout,*) '  Window size : ',ysc%CHE%TSEMISS(ji)%NWS
      WRITE(iluout,*) '  Current data INDEX : ',indx1,indx2
      IF (iverb >= 10) WRITE(iluout,*) '  FLUX : ',ysc%CHE%TSEMISS(ji)%XFWORK
    END IF
  END IF
END DO
! 
! Agregation : flux computation
!
zemis(:,:) = 0.
!
! Point on head of Pronostic variable list
! to cover the entire list.
IF (ysc%SV%NSV_AEREND > 0) THEN
 cnames=>ysc%SV%CSV(ysc%SV%NSV_CHSBEG:ysc%SV%NSV_AEREND)
ELSE
 cnames=>ysc%SV%CSV(ysc%SV%NSV_CHSBEG:ysc%SV%NSV_CHSEND)
END IF
 curpronos=>ysc%CHE%TSPRONOSLIST
DO WHILE(ASSOCIATED(curpronos))
  IF (curpronos%NAMINDEX > ineq) THEN
    WRITE(iluout,*) 'FATAL ERROR in CH_EMISSION_FLUXN : SIZE(ZEMIS,2) =',&
           ineq,', INDEX bugge =',curpronos%NAMINDEX  
    CALL abor1_sfx('CH_EMISSION_FLUXN: FATAL ERROR')
  END IF
  
  zemis(:,curpronos%NAMINDEX) = 0.
!
! Loop on the number of agreg. coeff.
  DO ji=1,curpronos%NBCOEFF
!   Compute agregated flux    
    zemis(:,curpronos%NAMINDEX) = zemis(:,curpronos%NAMINDEX)+&
            curpronos%XCOEFF(ji)*ysc%CHE%TSEMISS(curpronos%NEFINDEX(ji))%XFWORK(:)  
  END DO

  IF (iverb >= 6) THEN
    WRITE(iluout,*) 'Agregation for ',cnames(curpronos%NAMINDEX)
    IF (iverb >= 10) WRITE(iluout,*) 'ZEMIS = ',zemis(:,curpronos%NAMINDEX)
  END IF
  IF ((cnames(curpronos%NAMINDEX) == "CO") .AND. any(zemis(:,curpronos%NAMINDEX).GT.0.)) THEN
  zfco(:) = zemis(:,curpronos%NAMINDEX)
  gco = .true.
  END IF

  curpronos=>curpronos%NEXT
!
END DO
!
zdepot(:,:) = 0.
WHERE (psfsv(:,:) >= 0.) 
  zemis(:,:) = zemis(:,:) + psfsv(:,:)
ELSE WHERE
  zdepot(:,:) = psfsv(:,:)
END WHERE
!
IF ((lch_aero_flux).AND.(ysc%SV%NSV_AERBEG > 0)) THEN
  IF (gco) THEN
    CALL ch_aer_emission(zemis, prhoa, ysc%SV%CSV, ysc%SV%NSV_CHSBEG, pfco=zfco)
  ELSE
    CALL ch_aer_emission(zemis, prhoa, ysc%SV%CSV, ysc%SV%NSV_CHSBEG)
  ENDIF
END IF
!
psfsv(:,:) = psfsv(:,:) + zemis(:,:)
!
IF (lioinit) CALL end_io_surf_n(hprogram)
!
IF (iverb >= 6) WRITE(iluout,*) '******** END CH_EMISSION_FLUX  ********'
IF (lhook) CALL dr_hook('CH_EMISSION_FLUX_N',1,zhook_handle)
!
END SUBROUTINE ch_emission_flux_n