ecearth3/sources/tm5mp/proj/cb05/emission_rn222.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_RN222
!
! !DESCRIPTION: Perform RN222 emissions needed for TM5 CBM4 version. 
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION_RN222
  !
  ! !USES:
  !
  use GO,            only : gol, goErr, goPr
  use tm5_distgrid,  only : dgrid, get_distgrid, scatter
  use partools,      only : isRoot
  use global_types,  only : emis_data, d3_data
  use Dims,          only : nregions

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public   :: emission_rn222_init      ! allocate
  public   :: emission_rn222_declare   ! read monthly input
  public   :: emission_rn222_apply     ! distribute & add emissions to tracer array
  public   :: emission_rn222_done      ! deallocate
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter                   :: mname = 'emission_rn222'
  type(emis_data),  dimension(nregions), target :: emis_rn222
  type(emis_data),  dimension(nregions), target :: emis_rn222_sec
  !
  ! !REVISION HISTORY: 
  !    1 Oct 2010 - Achim Strunk - standardized routines name, new apply method
  !   28 Mar 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_RN222_INIT
  !
  ! !DESCRIPTION: Allocate space needed to handle the emissions
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_RN222_INIT( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)          ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - extracted from old 'declare' routine
    !   22 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/Emission_RN222_Init'
    integer                     ::  region
    integer                     ::  i1, i2, j1, j2

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

    do region = 1, nregions

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

       allocate(emis_rn222    (region)%surf(i1:i2,j1:j2) )
       allocate(emis_rn222_sec(region)%surf(i1:i2,j1:j2) )

    enddo

    status = 0

  END SUBROUTINE EMISSION_RN222_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_RN222_DONE
  !
  ! !DESCRIPTION: Free space after handling of the emissions
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_RN222_DONE( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)           ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - rename old 'free_emission_rn222'
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    do region = 1, nregions
       deallocate(emis_rn222    (region)%surf)
       deallocate(emis_rn222_sec(region)%surf)
    end do

    status = 0

  END SUBROUTINE EMISSION_RN222_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_RN222_DECLARE
  !
  ! !DESCRIPTION: Opens, reads and evaluates input files (per month). 
  !               Provides emissions on 2d/3d-arrays which are then added
  !               to tracers in routine *apply. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_RN222_DECLARE( status )
    !
    ! !USES:
    !
    use dims,            only : im, jm, sec_month, sec_year, newsrun, iglbsfc, nlat180, nlon360
    USE MDF,             ONLY : MDF_Open, MDF_NETCDF, MDF_READ, MDF_Inq_VarID, MDF_Get_Var, MDF_Close
    use toolbox,         only : coarsen_emission
    use chem_param,      only : irn222, xmrn222
    use emission_data,   only : emis_input_dir_rn222, msg_emis
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - update for new emission standard; rid of GAIA default
    !   28 Mar 2012 - P. Le Sager  - adapted for lon-lat MPI domain decomposition
    !                              - read only once
    !                              - use MDF
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    real, dimension(:,:), allocatable    :: RN222_temp     ! global emiss for reading
    type(emis_data), dimension(nregions) :: emis_rn222_glb ! global emiss for coarsening
    integer                              :: region, nlatsrc, nlonsrc
    integer, parameter                   :: add_field = 0
    integer, parameter                   :: amonth=1

    character(len=256)                   :: fname, vname
    integer                              :: varid, hid

    ! --- begin -----------------------------------------
    
    ! constant emission : read & regrid once only
    if(newsrun) then
       
       write(gol,'(" EMISS-INFO ------------- read Rn222 emissions -------------")'); call goPr
           
       nlatsrc = nlat180 ! or dgrid(iglbsfc)%jm_region, accessible thru get_distgrid(...)
       nlonsrc = nlon360 ! or dgrid(iglbsfc)%im_region, accessible thru get_distgrid(...)

       ! global array for reading
       if(isRoot)then
          allocate(RN222_temp(nlonsrc,nlatsrc))
       else
          allocate(RN222_temp(1,1))
       end if

       ! global arrays for coarsening    
       do region = 1, nregions
          
          if (isRoot)then
             allocate(emis_rn222_glb(region)%surf(im(region),jm(region)))
          else
             allocate(emis_rn222_glb(region)%surf(1,1))
          end if
       enddo
        
       if (isRoot) then  ! read data. Units= kg/s/ 1x1grid 

          fname=trim(emis_input_dir_rn222)//'/RN222_WMO2004.nc4'
          vname='Rn222_emis'
          
          CALL MDF_Open( TRIM(fname), MDF_NETCDF, MDF_READ, hid, status )
          IF_NOTOK_RETURN(status=1)
          CALL MDF_Inq_VarID( hid, TRIM(vname), varid, status )
          IF_NOTOK_RETURN(status=1)
          CALL MDF_Get_Var( hid, varid, RN222_temp, status )
          IF_NOTOK_RETURN(status=1)
          CALL MDF_Close( hid, status )
          IF_NOTOK_RETURN(status=1)

          ! info for full year
          call msg_emis( amonth, 'WMO','(2004)', 'Rn222', xmrn222, sum(RN222_temp*sec_year) )

          call coarsen_emission('Rn222_emis', nlonsrc, nlatsrc, RN222_temp, emis_rn222_glb, add_field, status)
          IF_NOTOK_RETURN(status=1)
          deallocate(RN222_temp)
          
       endif

       do region = 1, nregions
          call scatter(dgrid(region), emis_rn222_sec(region)%surf, emis_rn222_glb(region)%surf, 0, status)
          IF_NOTOK_RETURN(status=1)

          deallocate(emis_rn222_glb(region)%surf)
       end do

    end if

    do region = 1, nregions
       emis_rn222(region)%surf = emis_rn222_sec(region)%surf *sec_month   ! convert to [kg/month/box]
    end do
    
    ! ok
    status = 0

  END SUBROUTINE EMISSION_RN222_DECLARE
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_RN222_APPLY
  !
  ! !DESCRIPTION: Take monthly emissions, and
  !               - split them vertically
  !               - apply time splitting factors 
  !               - add them up (add_3d) 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_RN222_APPLY( region, status )
    !
    ! !USES:
    !
    use dims,          only : okdebug, itaur, nsrce, tref, lm
    use datetime,      only : tau2date
    use emission_data, only : emission_vdist_by_sector
    use emission_data, only : do_add_3d
    use chem_param,    only : irn222, xmrn222
    use emission_data, only : vd_class_name_len
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - updated to new emission methods
    !   28 Mar 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/emission_rn222_apply'
    integer, dimension(6)       :: idater
    real                        :: dtime, fraction
    integer                     :: imr, jmr, lmr, i1, i2, j1, j2
    type(d3_data)               :: emis3d

    character(len=vd_class_name_len) :: splittype

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

    if( okdebug ) then
       write(gol,*) 'start of emission_rn222_apply'; call goPr
    endif

    call tau2date(itaur(region),idater)
    dtime=float(nsrce)/(2*tref(region))    ! emissions are added in two steps...XYZECCEZYX.

    if( okdebug ) then
       write(gol,*) 'emission_rn222_apply in region ',region,' at date: ',idater, ' with time step:', dtime ; call goPr
    endif

    ! get a working structure for 3d emissions
    call get_distgrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    allocate( emis3d%d3(i1:i2,j1:j2,lm(region)) ) ; emis3d%d3 = 0.0

    ! default: no additional splitting
    fraction = 1.0

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

    splittype = 'surface'

    ! vertically distribute according to sector
    call emission_vdist_by_sector( splittype, 'RN', region, emis_rn222(region), emis3d, status )
    IF_NOTOK_RETURN(status=1;deallocate(emis3d%d3))

    ! add dataset
    call do_add_3d( region, irn222, i1, j1, emis3d%d3, xmrn222, xmrn222, status, fraction )
    IF_NOTOK_RETURN(status=1)

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

    deallocate( emis3d%d3 )

    status=0
    
  END SUBROUTINE EMISSION_RN222_APPLY
  !EOC
  

END MODULE EMISSION_RN222