ecearth3/sources/tm5mp/proj/cbm4/sources_sinks.F90

#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') 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:      SOURCES_SINKS
!
! !DESCRIPTION: Perform all calculations needed for CBM4 chemistry simulation
!               in TM5: this is mainly emissions, process updates after
!               changes in meteo, boundary, sedimentation, photolysis,...
!
!FD: all emission are converted to kg X (fmw) /month exceptions are mentioned in the code
!
!\\
!\\
! !INTERFACE: 
!
MODULE SOURCES_SINKS
  !
  ! !USES:
  !
  use GO, only : gol, goErr, goPr, goBug, goLabel

  implicit none
  private  
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  PUBLIC  ::  SOURCES_SINKS_INIT, SOURCES_SINKS_DONE  ! Init and Done methods
  PUBLIC  ::  SS_MONTHLY_UPDATE   ! monthly initialization (photolysis,..)
  PUBLIC  ::  SS_AFTER_READ_METEO_UPDATE ! Update SS after met fields are updated. Called from modelIntegration/Proces_update
  PUBLIC  ::  SOURCES_SINKS_APPLY ! apply SS
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'sources_sinks'
  !
  ! !REVISION HISTORY: 
  !   19 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------------

contains
  
  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    SOURCES_SINKS_INIT
  !
  ! !DESCRIPTION: switch ON required meteo; init emissions
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE SOURCES_SINKS_INIT( status )
    !
    ! !USES:
    !
    use meteo,     only : Set
    use meteodata, only : temper_dat, humid_dat, oro_dat, gph_dat
    use Meteodata, only : cp_dat, lsp_dat
    use Meteodata, only : cvl_dat, cvh_dat, tv_dat
    use Meteodata, only : ci_dat, sd_dat, swvl1_dat
    use Meteodata, only : t2m_dat, d2m_dat
    use Meteodata, only : u10m_dat, v10m_dat, lsmask_dat,ci_dat
#if defined(with_online_bvoc) || defined(with_online_nox)
    use Meteodata, only : skt_dat
#endif
#ifdef with_online_bvoc
    use emission_bvoc_data, only : online_bvoc_skt
#endif
#ifdef with_online_nox
    use online_nox_data, only : online_nox_skt
    use Meteodata,       only : src_dat, lsmask_dat
#endif
    use Meteodata,      only : ssr_dat,  sshf_dat, slhf_dat, ewss_dat,   nsss_dat
    use Meteodata,      only : u10m_dat, v10m_dat, src_dat,  albedo_dat, nveg
    use chem_rates,     only : rates    
#ifndef without_emission
    use emission,       only : Emission_Init
#endif
#ifndef without_sedimentation
    use sedimentation , only : Sedimentation_Init
#endif
#ifndef without_photolysis
    use photolysis    , only : Photolysis_Init
#endif
#ifndef without_boundary
    use boundary      , only : Boundary_Init, o3du, use_o3du
#endif
#ifdef with_m7
    use emission_dust, only : emission_dust_init
    use emission_ss,   only : emission_ss_init
#endif
    
    use GO,             only : TrcFile, Init, Done, ReadRc
    use global_data,    only : rcfile 
    use dims,           only : iglbsfc, nregions 
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)    ::  status
    !
    ! !REVISION HISTORY: 
    !    3 Oct 2012 - P. Le Sager - get Henry coeff (call rates)
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    integer, parameter          :: kr=31 ! standard unit to read auxiliary files
    character(len=*), parameter :: rname = mname//'/Sources_Sinks_Init'
    type(TrcFile)               :: rcF
    integer                     :: region, iveg 

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

    !--------------------------------------------------
    ! ** select meteo (cases not accounted for in the **_init procedures)    
    !--------------------------------------------------    
    do region = 1, nregions

#ifndef without_emission
       call Set( temper_dat(region), status, used=.true. )
       call Set(  humid_dat(region), status, used=.true. )
       call Set(    oro_dat(region), status, used=.true. )
       call Set(    gph_dat(region), status, used=.true. )
#endif

       ! other
       call Set(    cvl_dat(region), status, used=.true. )
       call Set(    cvh_dat(region), status, used=.true. )
       do iveg=1,nveg
          call Set(     tv_dat(region,nveg), status, used=.true. )
       enddo
       call Set(     ci_dat(region), status, used=.true. )
       call Set(     sd_dat(region), status, used=.true. )
       call Set(  swvl1_dat(region), status, used=.true. )
       call Set(    t2m_dat(region), status, used=.true. )
       call Set(    d2m_dat(region), status, used=.true. )
       call Set(    ssr_dat(region), status, used=.true. )
       call Set(   sshf_dat(region), status, used=.true. )
       call Set(   slhf_dat(region), status, used=.true. )
       call Set(   ewss_dat(region), status, used=.true. )
       call Set(   nsss_dat(region), status, used=.true. )
       call Set(   u10m_dat(region), status, used=.true. )
       call Set(   v10m_dat(region), status, used=.true. )
       call Set(    src_dat(region), status, used=.true. )
       call Set( albedo_dat(region), status, used=.true. )
    enddo

    ! special set for DMS
    call Set(  t2m_dat(iglbsfc), status, used=.true. ) 
    call Set( u10m_dat(iglbsfc), status, used=.true. )
    call Set( v10m_dat(iglbsfc), status, used=.true. )
    
    !--------------------------------------------------
    ! ** Henry coefficients (must be before sedimentation_init)
    !--------------------------------------------------    
    call rates(status)
    IF_NOTOK_RETURN(status=1)

    !--------------------------------------------------
    ! ** Sedimentation
    !--------------------------------------------------    
#ifndef without_sedimentation
    call Sedimentation_Init( status )
    IF_NOTOK_RETURN(status=1)    
#endif    

    !--------------------------------------------------
    ! ** Stratospheric boundary (must be before photolysis)
    !--------------------------------------------------    
#ifndef without_boundary
    call Boundary_Init( .true., status )
    IF_NOTOK_RETURN(status=1)
#endif
    
    !--------------------------------------------------
    ! ** Photolysis
    !--------------------------------------------------    
#ifndef without_photolysis
#ifdef without_boundary
    call photolysis_init(.true., kr )
#else
    if (use_o3du) then
       call Photolysis_Init(.true., kr, o3du )
    else
       call Photolysis_Init(.true., kr )
    end if
#endif
#endif

    !--------------------------------------------------
    ! ** Emissions    
    !--------------------------------------------------    
#ifndef without_emission
    call Emission_Init( status )
    IF_NOTOK_RETURN(status=1)

#ifdef with_online_bvoc
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'online.bvoc.skt', online_bvoc_skt, status )
    IF_NOTOK_RETURN(status=1)
    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)
    if (online_bvoc_skt) then
       call Set(    skt_dat(iglbsfc), status, used=.true. )
    else
       call Set(    t2m_dat(iglbsfc), status, used=.true. )
    endif
    call Set(    ssr_dat(iglbsfc), status, used=.true. )
#endif

#ifdef with_online_nox
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'online.nox.skt', online_nox_skt, status )
    IF_NOTOK_RETURN(status=1)
    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)
    if (online_nox_skt) then
       call Set(    skt_dat(iglbsfc), status, used=.true. )
    else
       call Set(    t2m_dat(iglbsfc), status, used=.true. )
    endif
    call Set(    lsp_dat(iglbsfc), status, used=.true. )
    call Set(     cp_dat(iglbsfc), status, used=.true. )
    call Set(    src_dat(iglbsfc), status, used=.true. )
    call Set( lsmask_dat(iglbsfc), status, used=.true. )
#endif

#endif /* EMISSIONS */

#ifdef with_m7
    call emission_dust_init( status )
    call emission_ss_init( status )
#endif
    
    !--------------------------------------------------
    ! ** Done
    !--------------------------------------------------    
    status = 0

  END SUBROUTINE SOURCES_SINKS_INIT
  !EOC

  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    SOURCES_SINKS_DONE
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE SOURCES_SINKS_DONE( status )
    !
    ! !USES:
    !
#ifndef without_photolysis
    use photolysis,    only: photolysis_done
#endif
#ifndef without_sedimentation
    use sedimentation, only: Sedimentation_Done
#endif
#ifndef without_boundary
    use Boundary,      only: Boundary_Done
#endif
#ifndef without_emission
    use emission,      only: Emission_Done
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)  ::  status
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

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

    ! --- begin --------------------------------
#ifndef without_photolysis
    call photolysis_done ( status )
    IF_NOTOK_RETURN(status=1)
#endif

#ifndef without_boundary
    call Boundary_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif

#ifndef without_sedimentation
    call Sedimentation_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif

#ifndef without_emission
    call Emission_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif

    status = 0

  END SUBROUTINE SOURCES_SINKS_DONE
  !EOC

  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    SS_MONTHLY_UPDATE
  !
  ! !DESCRIPTION: monthly (re)initialisation of sources/sinks 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE SS_MONTHLY_UPDATE( status )
    !
    ! !USES:
    !
    use dims,           only : mlen, sec_day, sec_month, sec_year
    use datetime,       only : calc_sm
#ifndef without_photolysis
    use photolysis    , only : photolysis_init
#endif
#ifndef without_emission
    use emission,       only : declare_emission
#endif
#ifndef without_boundary
    use boundary,       only : Boundary_Init, o3du, use_o3du
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)     ::  status
    !
    ! !REVISION HISTORY: 
    !   19 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !   (1) routine is called at the beginning of every month
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    character(len=*), parameter :: rname = mname//'/ss_monthly_update'
    integer,          parameter :: kr=31 ! standard unit to read auxiliary files

    ! --- begin ------------------------------------
    
    ! calculate some conversion factors related to time...
    call calc_sm( mlen, sec_day, sec_month, sec_year )

    ! Read monthly emissions
#ifndef without_emission
    call declare_emission( status )
    IF_NOTOK_RETURN(status=1)
#endif
    
    ! Monthly update for stratospheric boundary
#ifndef without_boundary
    call Boundary_Init( .false., status )
    IF_NOTOK_RETURN(status=1)
#endif
 
    ! Monthly update for photolysis    
#ifndef without_photolysis   
#ifndef without_boundary
    if (use_o3du) then
       call Photolysis_Init( .false., kr, o3du )
    end if
#endif
#endif
    
    status = 0

  END SUBROUTINE SS_MONTHLY_UPDATE
  !EOC
  
  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    SS_AFTER_READ_METEO_UPDATE
  !
  ! !DESCRIPTION: subroutine that is called after reading new met fields (clouds,
  !               surface winds, etc.).
  !               In this routine, 'chemistry' fields that depend on these
  !               data are calculated. Called from modelIntegration/Proces_update.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE SS_AFTER_READ_METEO_UPDATE( status )
    !
    ! !USES:
    !
    use dims,          only : nregions, sec_month
    use tm5_distgrid,  only : dgrid, Get_DistGrid
#ifndef without_photolysis
    use photolysis,    only : ozone_info_online, slingo, aerosol_info, update_csqy
#endif

#ifndef without_emission
    use emission_nox,  only : eminox_lightning
#ifndef without_convection
    use emission_nox,  only : lightningNOX
#endif
    use emission_dms,  only : getDMS
#if defined(with_online_bvoc) || defined(with_online_nox)
    use dims,          only : itau, ndyn_max
#endif
#ifdef with_online_bvoc
    use emission_bvoc, only : getBVOC
#endif
#ifdef with_online_nox
    use online_nox,    only : getNOx
    use emission_nox,  only : nat_nox
#endif
#ifdef with_m7
    use emission_dust, only : calc_emission_dust, read_emission_dust
    use emission_ss,   only : calc_emission_ss
#endif
#endif /* EMISSIONS */

#ifndef without_sedimentation    
    use sedimentation, only : calculate_rh
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)      ::  status
    !
    ! !REVISION HISTORY: 
    !   19 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = 'ss_after_read_meteo_update'
    integer :: region, i1, j1

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

    call goLabel()
    
#ifndef without_emission

#ifdef with_online_nox
    if (mod(itau, ndyn_max) == 0) then
       call getNOx(status) ! returns nat_nox in kg(N)/month
       IF_NOTOK_RETURN(status=1)
    endif
#endif

    !get dms_emissions dms_sea, flagged by new_surface meteo fields  get in kgS/month
    call getDMS( status )
    IF_NOTOK_RETURN(status=1)

#ifdef with_online_bvoc
    if (mod(itau, ndyn_max) == 0) then
       ! TvN: added fixed time step ndyn_max (1 hour)
       ! because the averages over the past 24 hours are calculated from hourly averages
       call getBVOC (status)
       IF_NOTOK_RETURN(status=1)
    endif
#endif

    ! Lightning NOx (defined only if convection is turned on)
#ifdef without_convection
    do region = 1, nregions
       eminox_lightning(region)%d3=0.
    end do
#else
    do region = 1, nregions

       call Get_DistGrid( dgrid(region), I_STRT=i1, J_STRT=j1 )

       call lightningNOX(region, i1, j1, eminox_lightning(region)%d3, status)
       IF_NOTOK_RETURN(status=1)
       eminox_lightning(region)%d3(:,:,:) = eminox_lightning(region)%d3(:,:,:)*sec_month   !from kg N/s ----> kg N/month

    end do
#endif
        
#ifdef with_m7
    call calc_emission_ss( status )
    IF_NOTOK_RETURN(status=1)
    call read_emission_dust( status ) ! this is active if (input.emis.dust /= ONLINE)
    IF_NOTOK_RETURN(status=1)
    call calc_emission_dust( status ) ! this is active if (input.emis.dust == ONLINE)
    IF_NOTOK_RETURN(status=1)
#endif
#endif /* EMISSIONS */
    
#ifndef without_photolysis
    do region = 1, nregions

       ! cloud optical depth
       call slingo(region)

       ! calculate optical depth ozone from current ozone field
       ! note: this routine does not depend on clouds but on rm ...
       call ozone_info_online(region)
       
       call update_csqy( region ) ! t/p-dependent cross-sections and quantum yields

    end do
#endif

#ifndef without_sedimentation
    call calculate_rh
#endif

    ! ok
    call goLabel()
    status = 0

  END SUBROUTINE SS_AFTER_READ_METEO_UPDATE
  !EOC
  
  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    SOURCES_SINKS_APPLY
  !
  ! !DESCRIPTION: this subroutine changes the tracer mass and its
  !               slopes by chemical sources. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE SOURCES_SINKS_APPLY( region, tr, status )
    !
    ! !USES:
    !
    use GO, only : TDate
    
#ifndef without_emission
    use emission, only: emission_apply
#endif
#ifndef without_sedimentation
    use sedimentation, only: Sedimentation_Apply
#endif    
#ifndef without_boundary
    use boundary, only : Boundary_Apply
#endif
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)        ::  region
    type(TDate)                ::  tr(2)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)       ::  status
    !
    ! !REVISION HISTORY: 
    !
    ! !REMARKS:
    ! - called each time step, during "source" step, by modelIntegration/do_steps
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC
    
    character(len=*), parameter ::  rname = mname//'/Sources_sinks_apply'
    
    ! --- begin ----------------------------------

#ifndef without_sedimentation
    call Sedimentation_Apply ( region, status )
    IF_NOTOK_RETURN(status=1) 
#endif

#ifndef without_emission
    ! note dust/ss emissions are ported to sedimentation routine
    call emission_apply( region, status )
    IF_NOTOK_RETURN(status=1)
#endif

    ! Apply boundary conditions to selected tracers
#ifndef without_chemistry
#ifndef without_boundary
    call Boundary_Apply( region, status )
    IF_NOTOK_RETURN(status=1)
#endif
#endif
    
    status = 0
    
  END SUBROUTINE SOURCES_SINKS_APPLY
  !EOC

END MODULE SOURCES_SINKS