!-------------------------------------------------------------------
! manages reading and interpolation of prescribed volcanic aerosol
! Created by: Francis Vitt
!-------------------------------------------------------------------
module prescribed_volcaero 7,5
use shr_kind_mod
, only : r8 => shr_kind_r8
use abortutils, only : endrun
use spmd_utils, only : masterproc
use tracer_data, only : trfld, trfile
use cam_logfile, only : iulog
implicit none
private
save
type(trfld), pointer :: fields(:)
type(trfile) :: file
public :: prescribed_volcaero_readnl
public :: prescribed_volcaero_register
public :: prescribed_volcaero_init
public :: prescribed_volcaero_adv
public :: write_prescribed_volcaero_restart
public :: read_prescribed_volcaero_restart
public :: has_prescribed_volcaero
public :: init_prescribed_volcaero_restart
logical :: has_prescribed_volcaero = .false.
character(len=8), parameter :: volcaero_name = 'VOLC_MMR'
character(len=13), parameter :: volcrad_name = 'VOLC_RAD_GEOM'
character(len=9), parameter :: volcmass_name = 'VOLC_MASS'
character(len=11), parameter :: volcmass_column_name = 'VOLC_MASS_C'
! These variables are settable via the namelist (with longer names)
character(len=16) :: fld_name = 'MMRVOLC'
character(len=256) :: filename = ''
character(len=256) :: filelist = ''
character(len=256) :: datapath = ''
character(len=32) :: data_type = 'SERIAL'
logical :: rmv_file = .false.
integer :: start_ymd = 0
integer :: start_tod = 0
integer :: rad_chnk_offset
real(r8), dimension(:,:,:), pointer :: radius
contains
!-------------------------------------------------------------------
!-------------------------------------------------------------------
subroutine prescribed_volcaero_readnl(nlfile) 1,15
use namelist_utils, only: find_group_name
use units, only: getunit, freeunit
use mpishorthand
character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
! Local variables
integer :: unitn, ierr
character(len=*), parameter :: subname = 'prescribed_volcaero_readnl'
character(len=16) :: prescribed_volcaero_name
character(len=256) :: prescribed_volcaero_file
character(len=256) :: prescribed_volcaero_filelist
character(len=256) :: prescribed_volcaero_datapath
character(len=32) :: prescribed_volcaero_type
logical :: prescribed_volcaero_rmfile
integer :: prescribed_volcaero_ymd
integer :: prescribed_volcaero_tod
namelist /prescribed_volcaero_nl/ &
prescribed_volcaero_name, &
prescribed_volcaero_file, &
prescribed_volcaero_filelist, &
prescribed_volcaero_datapath, &
prescribed_volcaero_type, &
prescribed_volcaero_rmfile, &
prescribed_volcaero_ymd, &
prescribed_volcaero_tod
!-----------------------------------------------------------------------------
! Initialize namelist variables from local module variables.
prescribed_volcaero_name = fld_name
prescribed_volcaero_file = filename
prescribed_volcaero_filelist = filelist
prescribed_volcaero_datapath = datapath
prescribed_volcaero_type = data_type
prescribed_volcaero_rmfile = rmv_file
prescribed_volcaero_ymd = start_ymd
prescribed_volcaero_tod = start_tod
! Read namelist
if (masterproc) then
unitn = getunit()
open( unitn, file=trim(nlfile), status='old' )
call find_group_name(unitn, 'prescribed_volcaero_nl', status=ierr)
if (ierr == 0) then
read(unitn, prescribed_volcaero_nl, iostat=ierr)
if (ierr /= 0) then
call endrun(subname // ':: ERROR reading namelist')
end if
end if
close(unitn)
call freeunit(unitn)
end if
#ifdef SPMD
! Broadcast namelist variables
call mpibcast(prescribed_volcaero_name, len(prescribed_volcaero_name), mpichar, 0, mpicom)
call mpibcast(prescribed_volcaero_file, len(prescribed_volcaero_file), mpichar, 0, mpicom)
call mpibcast(prescribed_volcaero_filelist, len(prescribed_volcaero_filelist), mpichar, 0, mpicom)
call mpibcast(prescribed_volcaero_datapath, len(prescribed_volcaero_datapath), mpichar, 0, mpicom)
call mpibcast(prescribed_volcaero_type, len(prescribed_volcaero_type), mpichar, 0, mpicom)
call mpibcast(prescribed_volcaero_rmfile, 1, mpilog, 0, mpicom)
call mpibcast(prescribed_volcaero_ymd, 1, mpiint, 0, mpicom)
call mpibcast(prescribed_volcaero_tod, 1, mpiint, 0, mpicom)
#endif
! Update module variables with user settings.
fld_name = prescribed_volcaero_name
filename = prescribed_volcaero_file
filelist = prescribed_volcaero_filelist
datapath = prescribed_volcaero_datapath
data_type = prescribed_volcaero_type
rmv_file = prescribed_volcaero_rmfile
start_ymd = prescribed_volcaero_ymd
start_tod = prescribed_volcaero_tod
! Turn on prescribed volcanics if user has specified an input dataset.
if (len_trim(filename) > 0 ) has_prescribed_volcaero = .true.
end subroutine prescribed_volcaero_readnl
!-------------------------------------------------------------------
!-------------------------------------------------------------------
subroutine prescribed_volcaero_register() 1,4
use ppgrid, only: pver
use phys_buffer, only: pbuf_add
integer :: idx
if (has_prescribed_volcaero) then
call pbuf_add(volcaero_name,'physpkg',1,pver,1,idx)
call pbuf_add(volcrad_name, 'physpkg',1,pver,1,idx)
endif
endsubroutine prescribed_volcaero_register
!-------------------------------------------------------------------
!-------------------------------------------------------------------
subroutine prescribed_volcaero_init() 1,12
use tracer_data, only : trcdata_init
use cam_history, only : addfld, phys_decomp
use ppgrid, only : pver
use error_messages, only: handle_err
use ppgrid, only: pcols, pver, begchunk, endchunk
use phys_buffer, only : pbuf, pbuf_get_fld_idx
implicit none
integer :: ndx, istat
character(len=32) :: specifier(1)
if ( has_prescribed_volcaero ) then
if ( masterproc ) then
write(iulog,*) 'volcanic aerosol is prescribed in :'//trim(filename)
endif
else
return
endif
specifier(1) = trim(volcaero_name)//':'//trim(fld_name)
file%in_pbuf = .true.
call trcdata_init( specifier, filename, filelist, datapath, fields, file, &
rmv_file, start_ymd, start_tod, data_type )
call addfld(volcaero_name,'kg/kg', pver, 'I', 'prescribed volcanic aerosol dry mass mixing ratio', phys_decomp )
call addfld(volcrad_name,'m', pver, 'I', 'volcanic aerosol geometric-mean radius', phys_decomp )
call addfld(volcmass_name,'kg/m^2', pver, 'I', 'volcanic aerosol vertical mass path in layer', phys_decomp )
call addfld(volcmass_column_name,'kg/m^2', 1, 'I', 'volcanic aerosol column mass', phys_decomp )
ndx = pbuf_get_fld_idx(volcrad_name, failcode=-1 )
radius => pbuf(ndx)%fld_ptr(1,:,:,:,1)
rad_chnk_offset = -begchunk+1
end subroutine prescribed_volcaero_init
!-------------------------------------------------------------------
!-------------------------------------------------------------------
subroutine prescribed_volcaero_adv( state ) 1,17
use tracer_data, only : advance_trcdata
use physics_types,only : physics_state
use ppgrid, only : begchunk, endchunk
use ppgrid, only : pcols, pver
use string_utils, only : to_lower, GLC
use cam_history, only : outfld
use physconst, only : amass => mwdry ! molecular weight dry air ~ kg/kmole
use physconst, only : boltz, gravit ! J/K/molecule
use tropopause, only : tropopause_find, TROP_ALG_TWMO, TROP_ALG_CLIMATE
implicit none
type(physics_state), intent(in):: state(begchunk:endchunk)
integer :: c,ncol,i,k
real(r8) :: to_mmr(pcols,pver)
real(r8), parameter :: molmass = 47.9981995_r8
real(r8) :: ptrop
real(r8) :: concvolc ! micrograms of wetted aerosol per cubic centimeter
real(r8) :: volcmass(pcols,pver)
real(r8) :: columnmass(pcols)
real(r8) :: mmrvolc
integer :: tropLev(pcols)
!WACCM-derived relation between mass concentration and wet aerosol radius in meters
real(r8),parameter :: radius_conversion = 1.9e-4_r8
if( .not. has_prescribed_volcaero ) return
call advance_trcdata( fields, file, state )
do c = begchunk,endchunk
radius(:,:,c+rad_chnk_offset) = 0._r8
enddo
! copy prescribed tracer fields into state svariable with the correct units
do c = begchunk,endchunk
ncol = state(c)%ncol
select case ( to_lower(trim(fields(1)%units(:GLC(fields(1)%units)))) )
case ("molec/cm3","/cm3","molecules/cm3","cm^-3","cm**-3")
to_mmr(:ncol,:) = (molmass*1.e6_r8*boltz*state(c)%t(:ncol,:))/(amass*state(c)%pmiddry(:ncol,:))
case ('kg/kg','mmr','kg kg-1')
to_mmr(:ncol,:) = 1._r8
case ('mol/mol','mole/mole','vmr','fraction')
to_mmr(:ncol,:) = molmass/amass
case default
write(iulog,*) 'prescribed_volcaero_adv: units = ',trim(fields(1)%units) ,' are not recognized'
call endrun('prescribed_volcaero_adv: units are not recognized')
end select
fields(1)%data(:ncol,:,c+fields(1)%chnk_offset) = to_mmr(:ncol,:) * fields(1)%data(:ncol,:,c+fields(1)%chnk_offset) ! mmr
call tropopause_find(state(c), tropLev, primary=TROP_ALG_TWMO, backup=TROP_ALG_CLIMATE)
do i = 1,ncol
do k = 1,pver
! set to zero below tropopause
if ( k >= tropLev(i) ) then
fields(1)%data(i,k,c+fields(1)%chnk_offset) = 0._r8
endif
mmrvolc = fields(1)%data(i,k,c+fields(1)%chnk_offset)
if (mmrvolc > 0._r8) then
concvolc = (mmrvolc * state(c)%pdel(i,k))/(gravit * state(c)%zm(i,k))
radius(i,k,c+rad_chnk_offset) = radius_conversion*(concvolc**(1._r8/3._r8))
endif
enddo
enddo
volcmass(:ncol,:) = fields(1)%data(:ncol,:,c+fields(1)%chnk_offset)*state(c)%pdel(:ncol,:)/gravit
columnmass(:ncol) = sum(volcmass(:ncol,:), 2)
call outfld( volcaero_name, fields(1)%data(:,:,c+fields(1)%chnk_offset), pcols, state(c)%lchnk)
call outfld( volcrad_name, radius(:,:,c+rad_chnk_offset), pcols, state(c)%lchnk)
call outfld( volcmass_name, volcmass, pcols, state(c)%lchnk)
call outfld( volcmass_column_name, columnmass, pcols, state(c)%lchnk)
enddo
end subroutine prescribed_volcaero_adv
!-------------------------------------------------------------------
subroutine init_prescribed_volcaero_restart( piofile ) 1,2
use pio, only : file_desc_t
use tracer_data, only : init_trc_restart
implicit none
type(file_desc_t),intent(inout) :: pioFile ! pio File pointer
call init_trc_restart( 'prescribed_volcaero', piofile, file )
end subroutine init_prescribed_volcaero_restart
!-------------------------------------------------------------------
subroutine write_prescribed_volcaero_restart( piofile ) 1,2
use tracer_data, only : write_trc_restart
use pio, only : file_desc_t
implicit none
type(file_desc_t) :: piofile
call write_trc_restart( piofile, file )
end subroutine write_prescribed_volcaero_restart
!-------------------------------------------------------------------
!-------------------------------------------------------------------
subroutine read_prescribed_volcaero_restart( pioFile ) 1,2
use tracer_data, only : read_trc_restart
use pio, only : file_desc_t
implicit none
type(file_desc_t) :: piofile
call read_trc_restart( 'prescribed_volcaero', piofile, file )
end subroutine read_prescribed_volcaero_restart
end module prescribed_volcaero