ch_emission_snapn.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_snap_n (CHN, &
                                     hprogram,ksize,psimtime,psuntime, &
                                    kyear,kmonth,kday,prhoa,plon      )
!     ######################################################################
!!
!!***  *CH_EMISSION_SNAP_n* - 
!!
!!    PURPOSE
!!    -------
!!      Return a time-dependent emission flux based on tabulated values
!!
!!**  METHOD
!!    ------
!!
!!    AUTHOR
!!    ------
!!    S. Queguiner
!!
!!    MODIFICATIONS
!!    -------------
!!    Original 10/2011
!!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!
USE modd_ch_snap_n, ONLY : ch_emis_snap_t
!
USE modd_csts,        ONLY: xday
!
USE modi_add_forecast_to_date_surf
USE modi_substract_to_date_surf
USE modi_ch_conversion_factor
!------------------------------------------------------------------------------
!
!*       0.   DECLARATIONS
!        -----------------
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
!
IMPLICIT NONE
!
!*       0.1  declaration of arguments
!
!
TYPE(ch_emis_snap_t), INTENT(INOUT) :: chn
!
 CHARACTER(LEN=6),       INTENT(IN)  :: hprogram! program calling surf. schemes
INTEGER,                INTENT(IN)  :: ksize   ! number of points
REAL,                   INTENT(IN)  :: psimtime! time of simulation in sec UTC
                                               ! (counting from midnight of
                                               ! the current day)
REAL, DIMENSION(KSIZE), INTENT(IN)  :: psuntime! Solar time (s since midnight)
INTEGER,                INTENT(IN)  :: kyear,kmonth,kday ! UTC year, month, day
REAL, DIMENSION(KSIZE), INTENT(IN)  :: prhoa   ! Air density
REAL, DIMENSION(KSIZE), INTENT(IN)  :: plon    ! Longitude (deg, from Greenwich)
!                                              ! (must be between -180deg and 180deg)
!
!*       0.2  declaration of local variables
!
REAL,   DIMENSION(KSIZE) :: zlon   ! Longitude centered in Greenwich meridian
REAL,   DIMENSION(KSIZE) :: ztime0
INTEGER,DIMENSION(KSIZE,2) :: iyear ! Year        at the begining of current hour
INTEGER,DIMENSION(KSIZE,2) :: imonth! Month       at the begining of current hour
INTEGER,DIMENSION(KSIZE,2) :: iday  ! Day         at the begining of current hour
INTEGER,DIMENSION(KSIZE,2) :: idow  ! Day of Week at the begining of current hour
INTEGER,DIMENSION(KSIZE,2) :: ihour ! Entire hour at the begining of current hour
REAL,   DIMENSION(KSIZE,2) :: ztime ! time (s)    at the begining of current hour
INTEGER                  :: jspec ! Loop counter on chemical species
INTEGER                  :: jsnap ! Loop counter on snap categories
INTEGER                  :: ji, jt    ! Loop counter on points
REAL,DIMENSION(KSIZE,2)    :: ze  ! Emissions at beginning and end of the hour
!
REAL(KIND=JPRB) :: zhook_handle
!
!------------------------------------------------------------------------------
!
!*    EXECUTABLE STATEMENTS
!     ---------------------
!
IF (lhook) CALL dr_hook('CH_EMISSION_SNAP_N',0,zhook_handle)
!
!------------------------------------------------------------------------------
!
!*  1.  Updates Conversion Factor (may depends on air density)
!       ------------------------------------------------------
!
 CALL ch_conversion_factor(chn, &
                           chn%CCONVERSION,prhoa(:))
!
!------------------------------------------------------------------------------
!
!*  2.  time and date for each point depending on Time reference
!       --------------------------------------------------------
!
!*  2.1 time at the beginning of current hour
!
iday(:,1)=kday
imonth(:,1)=kmonth
iyear(:,1)=kyear
!
SELECT CASE (chn%CSNAP_TIME_REF)
  CASE ('UTC  ')
    ztime0(:)=psimtime
  CASE ('SOLAR')
    zlon(:)=plon(:)
    WHERE(plon(:)>  180.) zlon(:)=plon(:)-360.
    WHERE(plon(:)<=-180.) zlon(:)=plon(:)+360.
    !*  retrieves solar date and time
    ztime0(:)=psimtime + zlon(:)*240. ! first guess is approximated solar time.
                                     ! The suntime should be close to this.
    DO ji=1,ksize
      IF (ztime0(ji)>psuntime(ji)+xday/2.) THEN
        ztime0(ji) = psuntime(ji) + xday
      ELSEIF (ztime0(ji)<psuntime(ji)-xday/2.) THEN
        ztime0(ji) = psuntime(ji) - xday
      ELSE
        ztime0(ji) = psuntime(ji)
      END IF
      CALL add_forecast_to_date_surf(iyear(ji,1),imonth(ji,1),iday(ji,1),ztime0(ji))
      CALL substract_to_date_surf(iyear(ji,1),imonth(ji,1),iday(ji,1),ztime0(ji))
    ENDDO
    
  CASE ('LEGAL')
    ztime0(:)=psimtime + chn%XDELTA_LEGAL_TIME(:) * 3600.
    DO ji=1,ksize
      CALL add_forecast_to_date_surf(iyear(ji,1),imonth(ji,1),iday(ji,1),ztime0(ji))
      CALL substract_to_date_surf(iyear(ji,1),imonth(ji,1),iday(ji,1),ztime0(ji))
    ENDDO

END SELECT
!
 CALL day_of_week(iday(:,1), imonth(:,1), iyear(:,1), idow(:,1))
!
ihour(:,1) = int((ztime0(:)+1.e-10)/3600.)! 1.E-10 and the where condition after are
WHERE (ihour(:,1)==24) ihour(:,1)=23      ! set to avoid computer precision problems
ztime(:,1) =  ihour(:,1)    * 3600.
!
!*   2.2 time at the end       of current hour
!
iday(:,2)=iday(:,1)
imonth(:,2)=imonth(:,1)
iyear(:,2)=iyear(:,1)
!
ztime(:,2) = (ihour(:,1)+1) * 3600.
!
DO ji=1,ksize
  CALL add_forecast_to_date_surf(iyear(ji,2),imonth(ji,2),iday(ji,2),ztime(ji,2))
ENDDO
!
 CALL day_of_week(iday(:,2), imonth(:,2), iyear(:,2), idow(:,2))
!
ihour(:,2)=nint(ztime(:,2))/3600
!
!------------------------------------------------------------------------------
!
!*  3.  Emission at the begining of the current hour
!       --------------------------------------------
!
 chn%XEMIS_FIELDS(:,:)=0.
!
DO jspec=1,chn%NEMIS_NBR
  !
  ze(:,:) = 0.
  !
  DO jsnap=1,chn%NEMIS_SNAP
    !
    DO jt=1,2
      !
      DO ji=1,ksize
        !
        ze(ji,jt) = ze(ji,jt) +  chn%XEMIS_FIELDS_SNAP(ji,jsnap,jspec) &
                      *chn%XSNAP_MONTHLY(imonth(ji,jt)  ,jsnap,jspec) &
                      *chn%XSNAP_DAILY  (idow(ji,jt)  ,jsnap,jspec) &
                      *chn%XSNAP_HOURLY (ihour(ji,jt)+1,jsnap,jspec) &
                      *chn%XCONVERSION(ji)
      ENDDO
      !
    ENDDO
    !
  ENDDO
!
!*  5.  Temporal interpolation within the current hour
!       ----------------------------------------------
!
  chn%XEMIS_FIELDS(:,jspec) = ze(:,1) + (ze(:,2)-ze(:,1))/3600.*(ztime0(:)-ihour(:,1)*3600.)

END DO
!
IF (lhook) CALL dr_hook('CH_EMISSION_SNAP_N',1,zhook_handle)
!
!-------------------------------------------------------------------------------
 CONTAINS
!
SUBROUTINE day_of_week(DATE, MONTH, YEAR, DOW)
!!    AUTHOR
!!    ------
!!    J.Arteta 
!!    Original   August 2010
!!
!!
!!    MODifICATIONS
!!    -------------
!!    S. Queguiner 10/2011  DAY:Monday->Sunday => DOW:1->7
!
IMPLICIT NONE
INTEGER, DIMENSION(:), INTENT(IN) :: date, month, year
INTEGER, DIMENSION(:), INTENT(OUT):: dow
INTEGER, DIMENSION(SIZE(DOW))     :: day, yr, mn, n1, n2
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('CH_EMISSION_SNAP_N:DAY_OF_WEEK',0,zhook_handle)
!
yr = year
mn = month
!
WHERE (mn.LE.2)
  mn = mn + 12
  yr = yr -1
END WHERE
!
n1 = (26 * (mn + 1)) /10
n2 = (125 * yr) / 100
day = (date + n1 + n2 - (yr / 100) + (yr / 400) - 1)
!
dow = mod(day,7) + 7
WHERE (dow.GT.7) dow = dow - 7
!
IF (lhook) CALL dr_hook('CH_EMISSION_SNAP_N:DAY_OF_WEEK',1,zhook_handle)
END SUBROUTINE day_of_week
!
END SUBROUTINE ch_emission_snap_n