ecearth3/sources/tm5mp/proj/cb05/emission_isop.F90

!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      EMISSION_ISOP
!
! !DESCRIPTION: data and methods for Isoprene emissions.
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION_ISOP
  !
  ! !USES:
  !
  use GO,            only : gol, goPr, goErr, goBug
  use tm5_distgrid,  only : dgrid, get_distgrid, scatter, gather
  use partools,      only : isRoot,  par_broadcast
  use dims,          only : nregions, idate, okdebug
  use global_types,  only : emis_data, d3_data
  use global_types,  only : isop_data
  use emission_read, only : used_providers_isop, has_isop_emis

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public   :: Emission_isop_Init
  public   :: Emission_isop_Done
  public   :: emission_isop_Declare
  public   :: Emission_isop_Apply
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'emission_isop'

  !isoprene emission distribution function
  type(isop_data),dimension(nregions),target       :: isop_dat          

  !flags the daily calculation of isoprene distribution func.
  integer,dimension(nregions)                      :: day_isop

  type( emis_data ), dimension(:,:), allocatable :: isop_emis_2d
  type( d3_data   ), dimension(:,:), allocatable :: isop_emis_3d
  
  integer :: j_start, j_end     ! indices for tropics
  
  integer              :: isop_3dsec,     isop_2dsec
  logical, allocatable :: has_data_3d(:), has_data_2d(:)
  !
  ! !REVISION HISTORY: 
  !    1 Oct 2010 - Achim Strunk - overhaul for AR5
  !   28 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

contains


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_ISOP_INIT
  !
  ! !DESCRIPTION: Allocate memory
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_ISOP_INIT( status )
    !
    ! !USES:
    !
    use dims,           only : ndyn_max, tref, im, jm, lm, idate, nlon360, nlat180
    use dims,           only : sec_month, dy, yref, ybeg
    use global_data,    only : rcfile
    use emission_read,  only : providers_def, numb_providers, sectors_def
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)     ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - adapted for AR5
    !   28 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/emission_isop_init'

    integer                   :: region, i1, i2, j1, j2
    integer                   :: lmr, lsec, ntim, lprov
    real                      :: dtime, dy1

    ! --- begin -----------------------------------------

    status = 0
    if(.not. has_isop_emis) return
    
    isop_2dsec = 0
    isop_3dsec = 0
    
    do lprov = 1, numb_providers
       if (count(used_providers_isop == providers_def(lprov)%name)/=0) then

          ! limit AR5 to the 2D biomass burning sector
          if ( trim(providers_def(lprov)%name) == 'AR5' ) then
             isop_2dsec = isop_2dsec + count( sectors_def%prov == 'AR5     ' &
                                        .and. sectors_def%catname == 'biomassburning')
          else
             ! use all sectors
             isop_2dsec = isop_2dsec + providers_def(lprov)%nsect2d
             isop_3dsec = isop_3dsec + providers_def(lprov)%nsect3d
          end if
       endif
    enddo

    allocate( isop_emis_2d( nregions, isop_2dsec ) )
    allocate( isop_emis_3d( nregions, isop_3dsec ) )
    allocate( has_data_2d(isop_2dsec)) ; has_data_2d=.false.
    allocate( has_data_3d(isop_3dsec)) ; has_data_3d=.false.

    ! allocate information arrays (2d and 3d)
    do region=1,nregions

       
       CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )          
       lmr = lm(region)

       do lsec=1,isop_2dsec
          allocate( isop_emis_2d(region,lsec)%surf(i1:i2, j1:j2) )
       end do

       do lsec=1,isop_3dsec
          allocate( isop_emis_3d(region,lsec)%d3(i1:i2, j1:j2,lmr) )
       end do

       dtime=float(ndyn_max)/(2*tref(region))    ! timestep emissions
       ntim=86400/nint(dtime)                    ! number of timesteps in 24 hours for this region

       allocate(isop_dat(region)%scalef_isop(ntim,j1:j2))   ! allocate enough memory!

    enddo
    
    ! find indices of tropics for splitting emissions
    dy1=dy/yref(1)
    j_start = floor((-20. - ybeg(1))/dy1) + 1
    j_end   = floor(( 20. - ybeg(1))/dy1) + 1
    
    ! ok
    status = 0

  end subroutine Emission_isop_Init
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_ISOP_DONE
  !
  ! !DESCRIPTION: Free memory
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Emission_isop_Done( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)          ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - adapted to new structures
    !   28 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/Emission_isop_Done'
    integer                     ::  region, lsec

    ! --- begin --------------------------------------

    status = 0
    if(.not. has_isop_emis) return

    do region = 1, nregions
       do lsec=1,isop_2dsec
          deallocate( isop_emis_2d(region,lsec)%surf )
       end do
       do lsec=1,isop_3dsec
          deallocate( isop_emis_3d(region,lsec)%d3 )
       end do
       deallocate(isop_dat(region)%scalef_isop )
    end do

    deallocate( isop_emis_2d, isop_emis_3d )
    deallocate( has_data_2d,  has_data_3d  )

    status = 0

  end subroutine Emission_isop_Done
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_ISOP_DECLARE
  !
  ! !DESCRIPTION: Opens, reads and evaluates input files (per month). 
  !               Provides emissions on 2d/3d-arrays which are then added
  !               to mixing ratios in routine *apply. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_ISOP_DECLARE( status )
    !
    ! !USES:
    !
    use toolbox,       only : coarsen_emission
    use dims,          only : im, jm, lm, idate, sec_month, nlon360, nlat180, iglbsfc
    use chem_param,    only : xmisop
    use emission_data, only : msg_emis, LCMIP6BMB, LAR5BMB, LMEGAN
    ! ---------------- CMIP6 - AR5 - MACC --------------------
    use emission_data, only : emis_input_year_nmvoc, emis_input_year_nat
    use emission_data, only : emis_input_dir_mac
    use emission_data, only : emis_input_dir_megan
    use emission_data, only : emis_input_dir_retro
    use emission_data, only : emis_input_dir_gfed
    use emission_read, only : emission_ar5_regrid_aircraft

    use emission_read, only : emission_cmip6_ReadSector
    use emission_read, only : emission_cmip6bmb_ReadSector
    use emission_read, only : emission_ar5_ReadSector
    use emission_read, only : emission_macc_ReadSector
    use emission_read, only : emission_megan_ReadSector
    use emission_read, only : emission_gfed_ReadSector
    use emission_read, only : emission_retro_ReadSector
    
    use emission_read, only : sectors_def, numb_sectors
    use emission_read, only : ar5_dim_3ddata
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - adapted for AR5
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter :: rname = mname//'/emission_isop_declare'
    integer                     :: region, hasData
    integer, parameter          :: add_field=0
    integer, parameter          :: amonth=2
    integer                     :: imr, jmr, lmr, lsec

    ! AR5
    real,dimension(:,:,:),allocatable  :: field3d, field3d2
    type(d3_data)                      :: emis3d,  work(nregions)
    type(emis_data)                    :: wrk2D(nregions)
    integer                            :: seccount2d, seccount3d

    ! --- begin -----------------------------------------

    status = 0
    if(.not. has_isop_emis) return
    
    write(gol,'(" EMISS-INFO ------------- read ISOPRENE emissions -------------")'); call goPr

    ! flag calculation of daily isoprene emission distribution    
    day_isop(1:nregions) = -1 

    do region = 1, nregions
       do lsec=1,isop_2dsec
          isop_emis_2d(region,lsec)%surf = 0.0
       end do
       do lsec=1,isop_3dsec
          isop_emis_3d(region,lsec)%d3   = 0.0
       end do
    end do

    ! global arrays for coarsening    
    do region = 1, nregions
       if (isRoot)then
          allocate(work(region)%d3(im(region),jm(region),ar5_dim_3ddata))
       else
          allocate(work(region)%d3(1,1,1))
       end if
    enddo
    
    do region = 1, nregions
       wrk2D(region)%surf => work(region)%d3(:,:,1)
    end do

    ! --------------------------------
    ! do a loop over available sectors
    ! --------------------------------

    ! count 2d and 3d sectors
    seccount2d = 0
    seccount3d = 0

    ! always allocate here 3d data set (for 2d sectors it will be filled in first layer only)
    if (isRoot) then
       allocate( field3d( nlon360, nlat180, ar5_dim_3ddata ) ) ; field3d = 0.0
    else
       allocate( field3d( 1, 1, 1 ) )
    end if

    
    sec : do lsec = 1, numb_sectors

       ! screen out unused providers
       if (count(used_providers_isop.eq.sectors_def(lsec)%prov).eq.0) cycle

       ! screen out CMIP6 and AR5 non-fires
       if ( sectors_def(lsec)%prov == 'AR5     ' .and. sectors_def(lsec)%catname /= 'biomassburning') cycle

       field3d = 0.0

       if( sectors_def(lsec)%f3d ) then 
          seccount3d = seccount3d + 1
       else 
          seccount2d = seccount2d + 1
       end if

       if (isRoot) then ! READ

          select case( trim(sectors_def(lsec)%prov) )

          case( 'CMIP6BMB' )

            call emission_cmip6bmb_ReadSector( 'NMVOC-'//'C5H8', emis_input_year_nmvoc, idate(2), lsec, field3d, status )
            IF_NOTOK_RETURN(status=1;deallocate(field3d))

          case( 'AR5' )

             ! here only if we select AR5, and if sector is BMB
             if (LAR5BMB) then
                call emission_ar5_ReadSector( 'isoprene', emis_input_year_nmvoc, idate(2), lsec, field3d, status )
                IF_NOTOK_RETURN(status=1)
             endif

          case( 'MACC' )

             ! skip if MEGAN is used, screen out sectors w/o isop else
             if( (trim(sectors_def(lsec)%name) .eq. 'emiss_bio') .and. (.not. LMEGAN) ) then

                if (trim(sectors_def(lsec)%catname) .eq. 'natural') then
                   call emission_macc_ReadSector( emis_input_dir_mac, 'ISOPRENE', emis_input_year_nat, idate(2), &
                        '0.5x0.5_kg.nc', sectors_def(lsec)%name, 'kg / s', field3d, status )
                   IF_NOTOK_RETURN(status=1)
                else
                   call emission_macc_ReadSector( emis_input_dir_mac, 'ISOPRENE', emis_input_year_nmvoc, idate(2), &
                        '0.5x0.5_kg.nc', sectors_def(lsec)%name, 'kg / s', field3d, status )
                   IF_NOTOK_RETURN(status=1)
                endif

             endif

          case('GFEDv3')

             call emission_gfed_ReadSector( emis_input_dir_gfed, 'isoprene', emis_input_year_nmvoc, idate(2), &
                  sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status )
             IF_NOTOK_RETURN(status=1)


          case('RETRO')

             call emission_retro_ReadSector( emis_input_dir_retro, 'ISOPRENE', emis_input_year_nmvoc, idate(2), &
                  sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status )
             IF_NOTOK_RETURN(status=1)
	     
          case('MEGAN')
	  	     
             call emission_megan_ReadSector( emis_input_dir_megan, 'isoprene', emis_input_year_nat, idate(2), &
                  sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status )
             IF_NOTOK_RETURN(status=1)
	     
          case('DUMMY')

          case default
             write(gol,*) "Error in buidling list of providers USED_PROVIDERS_ISOP" ; call goErr
             status=1; TRACEBACK; return

          end select
          
          ! nothing found??? 
          if( sum(field3d) < 100.*TINY(1.0) ) then
             if (okdebug) then
                write(gol,'("EMISS-INFO - no ISOP emissions found for ",a," ",a," for month ",i2 )') &
                     trim(sectors_def(lsec)%prov), trim(sectors_def(lsec)%name), idate(2) ; call goPr
             endif
             hasData=0
          else
             if (okdebug) then
                write(gol,'("EMISS-INFO - found ISOP emissions for ",a," ",a," for month ",i2 )') &
                     trim(sectors_def(lsec)%prov), trim(sectors_def(lsec)%name), idate(2) ; call goPr
             endif
             ! scale from kg/s to kg/month
             field3d = field3d * sec_month ! kg / month
             hasData=1
          end if

       end if

       call Par_broadcast(hasData, status)
       IF_NOTOK_RETURN(status=1)

       if (hasData == 0) then
          cycle sec
       else
          if ( sectors_def(lsec)%f3d ) then
             has_data_3d(seccount3d)=.true.
          else
             has_data_2d(seccount2d)=.true.
          end if
       end if

       ! Distinguish 2d/3d sectors
       if( sectors_def(lsec)%f3d ) then 

          write(gol,'("EMISS-ERROR - Unexpected 3D data: Uncomment code below ")'); call goErr
          status=1; TRACEBACK; return

          ! ---------------------------
          ! 3d data (AIRCRAFT)
          ! ---------------------------          
!           if (isRoot) then ! regrid

!              ! helper array for regridding
!              allocate( field3d2( nlon360,nlat180,lm(1) ) ) ; field3d2 = 0.0

!              ! aircraft data: regrid vertically to model layers
!              call emission_ar5_regrid_aircraft( iglbsfc, field3d, nlon360, nlat180, ar5_dim_3ddata, lm(1), field3d2, status )
!              IF_NOTOK_RETURN(status=1;deallocate(field3d,field3d2))

!              ! write some numbers
!              call msg_emis( amonth, 'aircraft. mass month', 'ISOP', xmisop, sum(field3d2) )

!              ! distribute to isop_emis in regions
!              call Coarsen_Emission( 'ISOP '//trim(sectors_def(lsec)%name), nlon360,  nlat180, lm(1),            & 
!                                                                            field3d2, work,    add_field, status )
!              IF_NOTOK_RETURN(status=1;deallocate(field3d,field3d2))

!              deallocate( field3d2 )

!           end if
          
!           do region = 1, nregions
!              call scatter(dgrid(region), isop_emis_3d(region,seccount3d)%d3, work(region)%d3, 0, status)
!              IF_NOTOK_RETURN(status=1)
!           end do

       else

          ! ---------------------------
          ! 2d data (Anthropogenic, Ships, Biomassburning)
          ! ---------------------------

          if (isRoot) then  ! print total & regrid
             
             call msg_emis( amonth, trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, 'ISOP', xmisop, sum(field3d(:,:,1)) )
             
             call coarsen_emission( 'ISOP '//sectors_def(lsec)%name, nlon360, nlat180, field3d(:,:,1), wrk2D, add_field, status )
             IF_NOTOK_RETURN(status=1)

          end if

          do region = 1, nregions
             call scatter(dgrid(region), isop_emis_2d(region,seccount2d)%surf, work(region)%d3(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)
          end do

       end if

       
    end do sec ! sectors


    deallocate( field3d )
    do region = 1, nregions
       if (associated(wrk2D(region)%surf)) nullify(wrk2D(region)%surf)
       deallocate( work(region)%d3 )
    end do

    status = 0

  END SUBROUTINE EMISSION_ISOP_DECLARE
  !EOC  
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_ISOP_APPLY
  !
  ! !DESCRIPTION: Take monthly emissions, and
  !               - split them vertically
  !               - apply time splitting factors 
  !               - add them up (add_isop) 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_ISOP_APPLY( region, status )
    !
    ! !USES:
    !
    use dims,          only : itaur, nsrce, tref
    use dims,          only : im, jm, lm, mlen, ndyn_max
    use dims,          only : yref, dy, ybeg
    use datetime,      only : tau2date, get_day
    use emission_data, only : emission_vdist_by_sector, bb_cycle
    use emission_data, only : emis_bb_trop_cycle
    use chem_param,    only : iisop, xmisop
    use emission_data, only : do_add_3d,   do_add_3d_cycle
    use emission_read, only : sectors_def, numb_sectors
    use tracer_data,   only : tracer_print
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)         :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        :: status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - adapted for AR5
    !   27 Apr 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter :: rname = mname//'/emission_isop_apply'
    integer, dimension(6)       :: idater
    real                        :: dtime, fraction
    integer                     :: imr, jmr, lmr, lsec
    integer                     :: seccount2d, seccount3d
    type(d3_data)               :: emis3d
    real                        :: dlatr
    integer                     :: day, ntim, i1, i2, j1, j2, j01, j02
    
    ! --- begin -----------------------------------------
    status = 0
    if(.not. has_isop_emis) return

    if( okdebug ) then
       write(gol,*) 'start of emission_apply_isop'; call goPr
    end if

    ! --- WHERE, WHEN
    
    CALL TAU2DATE( itaur(region), idater)
    CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    
    day   = GET_DAY(idater(2),idater(3),mlen)
    dtime = float(nsrce)/(2*tref(region))     !emissions are added in two steps...XYZECCEZYX.

    if(okdebug) then
       write(gol,*)'emission_isop_apply in region ',region,' at date: ',idater, ' with time step:', dtime
       call goPr
    end if
    
    ! --- DIURNAL CYCLE

    if ( day /= day_isop(region) ) then

       dlatr = dy/yref(region)  ! should be the same as lli%dlat_deg

       dtime=float(ndyn_max)/(2*tref(region))   ! timestep emissions
       ntim=86400/nint(dtime)                   ! number of timesteps in 24 hours.

       if ( okdebug ) then
          write (gol,*) 'day',day,'ntim',ntim,'jm',jm(region),'ybeg',&
               ybeg(region),'dlatr',dlatr,'sum(mlen)',sum(mlen),nsrce,dtime
          call goPr
       end if
       
       ! NOTE : latitude_start must be an integer. This is ok as long as
       ! the regions ybeg are integer.       
       call scale_isop(day, ntim, j1, isop_dat(region)%scalef_isop, ybeg(region), dlatr, sum(mlen))
       day_isop(region) = day   !avoid recalculation during this day!
    endif

    
    ! --- ADD EMISSIONS
    ! default: no additional splitting
    fraction = 1.0
    
    ! get a working structure for 3d emissions
    allocate( emis3d%d3(i1:i2,j1:j2,lm(region)) ) ; emis3d%d3 = 0.0

    ! count 2d and 3d sectors
    seccount2d = 0
    seccount3d = 0

    ! cycle over sectors
    do lsec = 1, numb_sectors

       ! screen out unused providers
       if (count(used_providers_isop.eq.sectors_def(lsec)%prov).eq.0) cycle

       ! screen out AR5 non-fires
       if ( sectors_def(lsec)%prov == 'AR5     ' .and. sectors_def(lsec)%catname /= 'biomassburning') cycle


       ! default: no additional splitting
       fraction = 1.0

       ! ----------------------------------------------------------------------------------------
       ! distinguish here between sectors and whether they should have additional splitting 
       ! if( sectors_def(lsec)%catname == 'biomassburning' ) fraction = fraction * bb_frac etc...
       ! ----------------------------------------------------------------------------------------

       ! distinguish between 2d/3d sectors
       if( sectors_def(lsec)%f3d ) then 

          seccount3d = seccount3d + 1
          if (.not.has_data_3d(seccount3d)) cycle
          emis3d%d3 = isop_emis_3d(region,seccount3d)%d3

       else

          seccount2d = seccount2d + 1
          if (.not.has_data_2d(seccount2d)) cycle
          emis3d%d3 = 0.0

          ! vertically distribute according to sector
          call emission_vdist_by_sector( sectors_def(lsec)%vdisttype, 'ISOP', region, isop_emis_2d(region,seccount2d), &
                                         emis3d, status)
          IF_NOTOK_RETURN(status=1)
                    
          ! is tile in tropics
          if ((j_end .lt. j1 ).or.(j_start .gt. j2) ) then
             continue
          else

             ! limit range
             j02=j_end
             j01=j_start
             if (j_end .gt. j2 ) j02 = j2
             if (j_start .lt. j1 ) j01 = j1
             
             emis3d%d3(:,j01:j02,1)=0.5*emis3d%d3(:,j01:j02,1)
             emis3d%d3(:,j01:j02,2)=emis3d%d3(:,j01:j02,1)
                          
          end if
          
       end if
              
       ! add dataset according to sector and category
       if( trim(sectors_def(lsec)%catname) == "natural" ) then 
          call do_add_isop( region, i1, j1, emis3d%d3, status, fraction )   ! special routine with additional emission scaling...
          IF_NOTOK_RETURN(status=1)

       elseif (trim(sectors_def(lsec)%catname) == "biomassburning") then

          if( emis_bb_trop_cycle ) then 
             call do_add_3d_cycle( region, iisop, i1, j1, emis3d%d3, bb_cycle(region)%scalef, xmisop, xmisop, status, fraction )
             IF_NOTOK_RETURN(status=1)
          else
             call do_add_3d( region, iisop, i1, j1, emis3d%d3, xmisop, xmisop, status, fraction )
             IF_NOTOK_RETURN(status=1)
          end if
       endif
       
    end do

    deallocate( emis3d%d3 )

    if(okdebug) then
       write(gol,*) 'end of emission_isop_apply' ; call goPr
    end if

    ! OK 
    status = 0

  END SUBROUTINE EMISSION_ISOP_APPLY
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:     SCALE_ISOP
  !
  ! !DESCRIPTION:  calculates factors to be multiplied with an emission field
  !                (e.g. isop) in order to simulate a diurnal cycle in
  !                emissions (caused by solar dependency), e.g. :
  !                
  !                rm(i,j,1,isop) = rm(i,j,1,isop) + e_isop(i,j)*scalef(ipos,j)*dt
  !                
  !                with ipos depending on time and longitude.
  !                
  !                The scalefactor is calculated for -180 longitude and the
  !                mean value for ntim timesteps during a day is scaled to 1.
  !                The routine should be called every day, since the position
  !                relative to the sun changes. The scaling is based on the
  !                cos(zenith) which is 1 for overhead sun and 0 at sunset.
  !                For the polar night, the values are set to 1.0 (even distribution)
  !                One thing is sure afterwards: we know that e_isop kg/s is added.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE SCALE_ISOP(day, ntim, j1, scalef, lat_start, dlat, idayy)
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: day       ! day number based on 365 days a year
    integer, intent(in) :: ntim      ! number of timesteps during one day
    integer, intent(in) :: j1        ! start index for latitudes
    integer, intent(in) :: lat_start ! southern edge of the domain (integer!)
    real,    intent(in) :: dlat      ! increment   (real!)
    integer, intent(in) :: idayy     ! # days this year
    !
    ! !OUTPUT PARAMETERS:
    !
    real, dimension(:,j1:), intent(out) :: scalef
    !
    ! !REVISION HISTORY: 
    !   10 May 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real              :: pi, piby, obliq, deday, delta, dt, lon, hr, lat, houra, smm
    integer           :: i, j

    pi    = acos(-1.0)
    piby  = pi/180.0
    obliq = 23.45 * piby
    deday = 4.88 + 2*pi/idayy  * day
    delta = asin(sin(obliq)*sin(deday))
    dt    = 24./ntim          ! timestep in hours
    lon   = -180.0*piby       ! calculate for dateline

    do j=j1,ubound(scalef,2)   ! over latitudes

       hr  = - 0.5*dt   ! shift half a interval back
       lat = (lat_start + (j-0.5)*dlat)*piby
       smm = 0.0

       do i=1,ntim
          hr = hr + dt
          houra = lon - pi + hr * (2.*pi/24.)
          scalef(i,j) = &
               max((sin(delta)*sin(lat)  + cos(delta)*cos(lat)*cos(houra)),0.0)
          smm = smm+scalef(i,j)/ntim
       end do

       if ( smm > 1e-5 ) then 
          scalef(1:ntim,j) = scalef(1:ntim,j)/smm 
       else 
          scalef(1:ntim,j) = 1.0
       end if

    end do

  END SUBROUTINE SCALE_ISOP
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DO_ADD_ISOP
  !
  ! !DESCRIPTION: Add local time splitted emissions to chemical array
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DO_ADD_ISOP(region, i1, j1, em_isop, status, fraction)
    !
    ! !USES:
    !
    use dims,          only : tref, nsrce, itaur, ndyn_max
    use dims,          only : sec_month, lm, dx, xref, xbeg
    use global_data,   only : mass_dat
    use chem_param,    only : iisop, xmisop
#ifdef with_budgets
    use budget_global, only : budg_dat, nzon_vg
    use emission_data, only : sum_emission, budemi_dat
#endif
    use datetime,      only : tau2date
    !
    ! !INPUT PARAMETERS:
    !
    integer,        intent(in) :: region, i1, j1
    real,           intent(in) :: em_isop(i1:,j1:,:)
    real, optional, intent(in) :: fraction
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)  :: status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - added fraction input
    !                              - switch from 2D to 3D input => remove tropics case
    !   10 May 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !     (1) scalef_isop has been calculated w/ ndyn_max. 
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/do_add_isop'

    real                              :: x, xlon, dtime, month2dt, dtime2
    integer                           :: i, j, l, ipos, nzone, nzone_v 
    integer                           :: sec_in_day, ntim, itim, i2, j2
    integer, dimension(6)             :: idater
    real, dimension(:,:,:,:), pointer :: rm
    
    ! --- begin -----------------------------------------
    
    rm  => mass_dat(region)%rm

    CALL GET_DISTGRID( dgrid(region), I_STOP=i2, J_STOP=j2 )

    dtime    = float(nsrce)    / (2*tref(region)) ! timestep emissions
    dtime2   = float(ndyn_max) / (2*tref(region)) ! timestep emissions based on non-CFL interference (CMK 05/2006)
    month2dt = dtime/sec_month                    ! conversion to emission per timestep (CFl induced timestep CMK 05/2006)
    ntim     = 86400/nint(dtime2)                 ! number of timesteps in 24 hours based on non-CFL interference (CMK 05/2006)

    call tau2date(itaur(region),idater)
    sec_in_day = idater(4)*3600 + idater(5)*60 + idater(6) ! elapsd seconds this day
    itim       = sec_in_day/nint(dtime2)+1                 ! # time interval  CMK 05/2006


    do l= 1,lm(region)
    do j=j1,j2
    do i=i1,i2

       xlon = xbeg(region) + (i-0.5)*dx/xref(region)
       ! itim = 1 and lon = -180 --->position 1
       ! itim = ntim ant lon = -180 --->position ntim, etc.
       ! itim = 1 and lon = 0 ---->position ntim/2
       ipos = 1 + mod(int((xlon+180.)*ntim/360.) + (itim-1),ntim)  !position in array depending on time and lon.
       x=em_isop(i,j,l)*isop_dat(region)%scalef_isop(ipos,j)*month2dt

       rm(i,j,l,iisop)=rm(i,j,l,iisop)+x

#ifdef with_budgets
       nzone=budg_dat(region)%nzong(i,j)    !global budget
       nzone_v=nzon_vg(l)
       budemi_dat(region)%emi(i,j,nzone_v,iisop) = &
            budemi_dat(region)%emi(i,j,nzone_v,iisop) + x/xmisop*1e3
       if(iisop ==1) sum_emission(region) = sum_emission(region) + x  !in kg
#endif
       !--
       !AJS: strange statement: injection of emission decreases vertical slope ?
       !AJS: this gives negative slopes in case of large emissions ...
       !if(adv_scheme.eq.'slope') rzm(i,j,L,iisop)=rzm(i,j,L,iisop)-x
       !--
       
    enddo !i
    enddo !j
    enddo !l   

    nullify(rm)
    status=0

  END SUBROUTINE DO_ADD_ISOP
  !EOC

END MODULE EMISSION_ISOP