module cam3_aero_data 4,13
!-----------------------------------------------------------------------
!
! Purposes:
! read, store, interpolate, and return fields
! of aerosols to CAM. The initialization
! file (mass.nc) is assumed to be a monthly climatology
! of aerosols from MATCH (on a sigma pressure
! coordinate system).
! also provide a "background" aerosol field to correct
! for any deficiencies in the physical parameterizations
! This fields is a "tuning" parameter.
! Public methods:
! (1) - initialization
! read aerosol masses from external file
! also pressure coordinates
! convert from monthly average values to mid-month values
! (2) - interpolation (time and vertical)
! interpolate onto pressure levels of CAM
! interpolate to time step of CAM
! return mass of aerosols
!
!-----------------------------------------------------------------------
use shr_kind_mod
, only: r8 => shr_kind_r8
use shr_scam_mod, only: shr_scam_GetCloseLatLon
use spmd_utils, only: masterproc
use ppgrid, only: pcols, pver, pverp, begchunk, endchunk
use phys_grid, only: get_ncols_p, scatter_field_to_chunk
use time_manager, only: get_curr_calday
use infnan, only: inf, bigint
use abortutils, only: endrun
use scamMod, only: scmlon,scmlat,single_column
use error_messages, only: handle_ncerr
use physics_types, only: physics_state
use boundarydata, only: boundarydata_init, boundarydata_type
use perf_mod, only: t_startf, t_stopf
use cam_logfile, only: iulog
use netcdf
implicit none
private
save
public :: &
cam3_aero_data_readnl, & ! read namelist
cam3_aero_data_register, & ! register these aerosols with pbuf
cam3_aero_data_init, & ! read from file, interpolate onto horiz grid
cam3_aero_data_timestep_init ! update data-aerosols to this timestep
! namelist variables
logical, public :: cam3_aero_data_on = .false.
character(len=256) :: bndtvaer = 'bndtvaer' ! full pathname for time-variant aerosol mass climatology dataset
! naer is number of species in climatology
integer, parameter :: naer = 11
real, parameter :: wgt_sscm = 6.0 / 7.0 ! Fraction of total seasalt mass in coarse mode
! indices to aerosol array (species portion)
integer, parameter :: &
idxSUL = 1, &
idxSSLTA = 2, & ! accumulation mode
idxSSLTC = 3, & ! coarse mode
idxOCPHO = 8, &
idxBCPHO = 9, &
idxOCPHI = 10, &
idxBCPHI = 11
! indices to sections of array that represent
! groups of aerosols
integer, parameter :: &
idxSSLTfirst = 2, numSSLT = 2, &
idxDUSTfirst = 4, &
numDUST = 4, &
idxCARBONfirst = 8, &
numCARBON = 4
! names of aerosols are they are represented in
! the climatology file.
! Appended '_V' indicates field has been vertically summed.
character(len=8), parameter :: aerosol_name(naer) = &
(/"MSUL_V "&
,"MSSLTA_V"&
,"MSSLTC_V"&
,"MDUST1_V"&
,"MDUST2_V"&
,"MDUST3_V"&
,"MDUST4_V"&
,"MOCPHO_V"&
,"MBCPHO_V"&
,"MOCPHI_V"&
,"MBCPHI_V"/)
! number of different "groups" of aerosols
integer, parameter :: num_aer_groups=4
! which group does each bin belong to?
integer, dimension(naer), parameter :: &
group =(/1,2,2,3,3,3,3,4,4,4,4/)
! name of each group
character(len=10), dimension(num_aer_groups), parameter :: &
aerosol_names = (/'sul ','sslt ','dust ','car '/)
! this boundarydata_type is used for datasets in the ncols format only.
type(boundarydata_type) :: aerosol_datan
integer :: aernid = -1 ! netcdf id for aerosol file (init to invalid)
integer :: species_id(naer) = -1 ! netcdf_id of each aerosol species (init to invalid)
integer :: Mpsid ! netcdf id for MATCH PS
integer :: nm = 1 ! index to prv month in array. init to 1 and toggle between 1 and 2
integer :: np = 2 ! index to nxt month in array. init to 2 and toggle between 1 and 2
integer :: mo_nxt = bigint ! index to nxt month in file
real(r8) :: cdaym = inf ! calendar day of prv month
real(r8) :: cdayp = inf ! calendar day of next month
! aerosol mass
real(r8), allocatable :: aer_mass(:, :, :, :)
! Days into year for mid month date
! This variable is dumb, the dates are in the dataset to be read in but they are
! slightly different than this so getting rid of it causes a change which
! exceeds roundoff.
real(r8) :: Mid(12) = (/16.5_r8, 46.0_r8, 75.5_r8, 106.0_r8, 136.5_r8, 167.0_r8, &
197.5_r8, 228.5_r8, 259.0_r8, 289.5_r8, 320.0_r8, 350.5_r8 /)
! values read from file and temporary values used for interpolation
!
! aerosolc is:
! Cumulative Mass at midpoint of each month
! on CAM's horizontal grid (col)
! on MATCH's levels (lev)
! aerosolc
integer, parameter :: paerlev = 28 ! number of levels for aerosol fields (MUST = naerlev)
integer :: naerlev ! size of level dimension in MATCH data
integer :: naerlon
integer :: naerlat
real(r8), pointer :: M_hybi(:) ! MATCH hybi
real(r8), pointer :: M_ps(:,:) ! surface pressure from MATCH file
real(r8), pointer :: aerosolc(:,:,:,:,:) ! Aerosol cumulative mass from MATCH
real(r8), pointer :: M_ps_cam_col(:,:,:) ! PS from MATCH on Cam Columns
! indices for fields in the physics buffer
integer :: cam3_sul_idx, cam3_ssam_idx, cam3_sscm_idx, &
cam3_dust1_idx, cam3_dust2_idx, cam3_dust3_idx, cam3_dust4_idx,&
cam3_ocpho_idx, cam3_bcpho_idx, cam3_ocphi_idx, cam3_bcphi_idx
!================================================================================================
contains
!================================================================================================
subroutine cam3_aero_data_readnl(nlfile) 1,9
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 = 'cam3_aero_data_readnl'
namelist /cam3_aero_data_nl/ cam3_aero_data_on, bndtvaer
!-----------------------------------------------------------------------------
if (masterproc) then
unitn = getunit()
open( unitn, file=trim(nlfile), status='old' )
call find_group_name(unitn, 'cam3_aero_data_nl', status=ierr)
if (ierr == 0) then
read(unitn, cam3_aero_data_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(cam3_aero_data_on, 1, mpilog, 0, mpicom)
call mpibcast(bndtvaer, len(bndtvaer), mpichar, 0, mpicom)
#endif
end subroutine cam3_aero_data_readnl
!================================================================================================
subroutine cam3_aero_data_register 1,12
! register old prescribed aerosols with physics buffer
use phys_buffer, only: pbuf_add
call pbuf_add('cam3_sul', 'physpkg', 1, pver, 1, cam3_sul_idx)
call pbuf_add('cam3_ssam', 'physpkg', 1, pver, 1, cam3_ssam_idx)
call pbuf_add('cam3_sscm', 'physpkg', 1, pver, 1, cam3_sscm_idx)
call pbuf_add('cam3_dust1', 'physpkg', 1, pver, 1, cam3_dust1_idx)
call pbuf_add('cam3_dust2', 'physpkg', 1, pver, 1, cam3_dust2_idx)
call pbuf_add('cam3_dust3', 'physpkg', 1, pver, 1, cam3_dust3_idx)
call pbuf_add('cam3_dust4', 'physpkg', 1, pver, 1, cam3_dust4_idx)
call pbuf_add('cam3_ocpho', 'physpkg', 1, pver, 1, cam3_ocpho_idx)
call pbuf_add('cam3_bcpho', 'physpkg', 1, pver, 1, cam3_bcpho_idx)
call pbuf_add('cam3_ocphi', 'physpkg', 1, pver, 1, cam3_ocphi_idx)
call pbuf_add('cam3_bcphi', 'physpkg', 1, pver, 1, cam3_bcphi_idx)
end subroutine cam3_aero_data_register
!================================================================================================
subroutine cam3_aero_data_init(phys_state) 1,38
!------------------------------------------------------------------
! Reads in:
! file from which to read aerosol Masses on CAM grid. Currently
! assumed to be MATCH ncep runs, averaged by month.
! NOTE (Data have been externally interpolated onto CAM grid
! and backsolved to provide Mid-month values)
!
! Populates:
! module variables:
! aerosolc(pcols,paerlev+1,begchunk:endchunk,naer,2))
! aerosolc( column_index
! , level_index (match levels)
! , chunk_index
! , species_index
! , month = 1:2 )
! M_hybi(level_index = Lev_MATCH) = pressure at mid-level.
! M_ps_cam_col(column,chunk,month) ! PS from MATCH on Cam Columns
!
! Method:
! read data from file
! allocate memory for storage of aerosol data on CAM horizontal grid
! distribute data to remote nodes
! populates the module variables
!
!------------------------------------------------------------------
use ioFileMod, only: getfil
#if ( defined SPMD )
use mpishorthand
#endif
type(physics_state), intent(in) :: phys_state(begchunk:endchunk)
! local variables
integer :: naerlev
integer dateid ! netcdf id for date variable
integer secid ! netcdf id for seconds variable
integer londimid ! netcdf id for longitude dimension
integer latdimid ! netcdf id for latitude dimension
integer levdimid ! netcdf id for level dimension
integer timesiz ! number of time samples (=12) in netcdf file
integer latid ! netcdf id for latitude variable
integer Mhybiid ! netcdf id for MATCH hybi
integer timeid ! netcdf id for time variable
integer dimids(nf90_max_var_dims) ! variable shape
integer :: start(4) ! start vector for netcdf calls
integer :: kount(4) ! count vector for netcdf calls
integer mo ! month index
integer m ! constituent index
integer :: n ! loop index
integer :: i,j,k ! spatial indices
integer :: date_aer(12) ! Date on aerosol dataset (YYYYMMDD)
integer :: attnum ! attribute number
integer :: ierr ! netcdf return code
real(r8) :: coldata(paerlev) ! aerosol field read in from dataset
integer :: ret
integer mo_prv ! index to previous month
integer latidx,lonidx
character(len=8) :: aname ! temporary aerosol name
character(len=8) :: tmp_aero_name(naer) ! name for input to boundary data
character(len=256) :: locfn ! netcdf local filename to open
!
! aerosol_data will be read in from the aerosol boundary dataset, then scattered to chunks
! after filling in the bottom level with zeros
!
real(r8), allocatable :: aerosol_data(:,:,:) ! aerosol field read in from dataset
real(r8), allocatable :: aerosol_field(:,:,:) ! (plon,paerlev+1,plat) aerosol field to be scattered
real(r8) :: caldayloc ! calendar day of current timestep
real(r8) :: closelat,closelon
character(len=*), parameter :: subname = 'cam3_aero_data_init'
!------------------------------------------------------------------
call t_startf(subname)
allocate (aer_mass(pcols, pver, naer, begchunk:endchunk) )
! set new aerosol names because input file has 1 seasalt bin
do m = 1, naer
tmp_aero_name(m)=aerosol_name(m)
if (aerosol_name(m)=='MSSLTA_V') tmp_aero_name(m) = 'MSSLT_V'
if (aerosol_name(m)=='MSSLTC_V') tmp_aero_name(m) = 'MSSLT_V'
end do
allocate (aerosolc(pcols,paerlev+1,begchunk:endchunk,naer,2))
aerosolc(:,:,:,:,:) = 0._r8
caldayloc = get_curr_calday ()
if (caldayloc < Mid(1)) then
mo_prv = 12
mo_nxt = 1
else if (caldayloc >= Mid(12)) then
mo_prv = 12
mo_nxt = 1
else
do i = 2 , 12
if (caldayloc < Mid(i)) then
mo_prv = i-1
mo_nxt = i
exit
end if
end do
end if
! Set initial calendar day values
cdaym = Mid(mo_prv)
cdayp = Mid(mo_nxt)
if (masterproc) &
write(iulog,*) subname//': CAM3 prescribed aerosol dataset is: ', trim(bndtvaer)
call getfil (bndtvaer, locfn, 0)
call handle_ncerr( nf90_open (locfn, 0, aernid),&
subname, __LINE__)
if (single_column) &
call shr_scam_GetCloseLatLon(aernid,scmlat,scmlon,closelat,closelon,latidx,lonidx)
! Check to see if this dataset is in ncol format.
aerosol_datan%isncol=.false.
ierr = nf90_inq_dimid( aernid, 'ncol', londimid )
if ( ierr==NF90_NOERR ) then
aerosol_datan%isncol=.true.
call handle_ncerr(nf90_close(aernid),subname, __LINE__)
call boundarydata_init(bndtvaer, phys_state, tmp_aero_name, naer, &
aerosol_datan, 3)
aerosolc(:,1:paerlev,:,:,:)=aerosol_datan%fields
M_ps_cam_col=>aerosol_datan%ps
M_hybi=>aerosol_datan%hybi
else
! Allocate memory for dynamic arrays local to this module
allocate (M_ps_cam_col(pcols,begchunk:endchunk,2))
allocate (M_hybi(paerlev+1))
! TBH: HACK to avoid use of uninitialized values when ncols < pcols
M_ps_cam_col(:,:,:) = 0._r8
if (masterproc) then
! First ensure dataset is CAM-ready
call handle_ncerr(nf90_inquire_attribute (aernid, nf90_global, 'cam-ready', attnum=attnum),&
subname//': interpaerosols needs to be run to create a cam-ready aerosol dataset')
! Get and check dimension info
call handle_ncerr( nf90_inq_dimid( aernid, 'lon', londimid ),&
subname, __LINE__)
call handle_ncerr( nf90_inq_dimid( aernid, 'lev', levdimid ),&
subname, __LINE__)
call handle_ncerr( nf90_inq_dimid( aernid, 'time', timeid ),&
subname, __LINE__)
call handle_ncerr( nf90_inq_dimid( aernid, 'lat', latdimid ),&
subname, __LINE__)
call handle_ncerr( nf90_inquire_dimension( aernid, londimid, len=naerlon ),&
subname, __LINE__)
call handle_ncerr( nf90_inquire_dimension( aernid, levdimid, len=naerlev ),&
subname, __LINE__)
call handle_ncerr( nf90_inquire_dimension( aernid, latdimid, len=naerlat ),&
subname, __LINE__)
call handle_ncerr( nf90_inquire_dimension( aernid, timeid, len=timesiz ),&
subname, __LINE__)
call handle_ncerr( nf90_inq_varid( aernid, 'date', dateid ),&
subname, __LINE__)
call handle_ncerr( nf90_inq_varid( aernid, 'datesec', secid ),&
subname, __LINE__)
do m = 1, naer
aname=aerosol_name(m)
! rename because file has only one seasalt field
if (aname=='MSSLTA_V') aname = 'MSSLT_V'
if (aname=='MSSLTC_V') aname = 'MSSLT_V'
call handle_ncerr( nf90_inq_varid( aernid, TRIM(aname), species_id(m)), &
subname, __LINE__)
end do
call handle_ncerr( nf90_inq_varid( aernid, 'lat', latid ),&
subname, __LINE__)
! quick sanity check on one field
call handle_ncerr( nf90_inquire_variable (aernid, species_id(1), dimids=dimids),&
subname, __LINE__)
if ( (dimids(4) /= timeid) .or. &
(dimids(3) /= levdimid) .or. &
(dimids(2) /= latdimid) .or. &
(dimids(1) /= londimid) ) then
write(iulog,*) subname//': Data must be ordered time, lev, lat, lon'
write(iulog,*) 'data are ordered as', dimids(4), dimids(3), dimids(2), dimids(1)
write(iulog,*) 'data should be ordered as', timeid, levdimid, latdimid, londimid
call endrun ()
end if
! use hybi,PS from MATCH
call handle_ncerr( nf90_inq_varid( aernid, 'hybi', Mhybiid ),&
subname, __LINE__)
call handle_ncerr( nf90_inq_varid( aernid, 'PS', Mpsid ),&
subname, __LINE__)
! check dimension order for MATCH's surface pressure
call handle_ncerr( nf90_inquire_variable (aernid, Mpsid, dimids=dimids),&
subname, __LINE__)
if ( (dimids(3) /= timeid) .or. &
(dimids(2) /= latdimid) .or. &
(dimids(1) /= londimid) ) then
write(iulog,*) subname//': Pressure must be ordered time, lat, lon'
write(iulog,*) 'data are ordered as', dimids(3), dimids(2), dimids(1)
write(iulog,*) 'data should be ordered as', timeid, levdimid, latdimid, londimid
call endrun ()
end if
! read in hybi from MATCH
call handle_ncerr( nf90_get_var (aernid, Mhybiid, M_hybi),&
subname, __LINE__)
! Retrieve date and sec variables.
call handle_ncerr( nf90_get_var (aernid, dateid, date_aer),&
subname, __LINE__)
if (timesiz < 12) then
write(iulog,*) subname//': When cycling aerosols, dataset must have 12 consecutive ', &
'months of data starting with Jan'
write(iulog,*) 'Current dataset has only ',timesiz,' months'
call endrun ()
end if
do mo = 1,12
if (mod(date_aer(mo),10000)/100 /= mo) then
write(iulog,*) subname//': When cycling aerosols, dataset must have 12 consecutive ', &
'months of data starting with Jan'
write(iulog,*)'Month ',mo,' of dataset says date=',date_aer(mo)
call endrun ()
end if
end do
if (single_column) then
naerlat=1
naerlon=1
endif
kount(:) = (/naerlon,naerlat,paerlev,1/)
end if ! masterproc
! broadcast hybi to nodes
#if ( defined SPMD )
call mpibcast (M_hybi, paerlev+1, mpir8, 0, mpicom)
call mpibcast (kount, 3, mpiint, 0, mpicom)
naerlon = kount(1)
naerlat = kount(2)
#endif
allocate(aerosol_field(kount(1),kount(3)+1,kount(2)))
allocate(M_ps(kount(1),kount(2)))
if (masterproc) allocate(aerosol_data(kount(1),kount(2),kount(3)))
! Retrieve Aerosol Masses (kg/m^2 in each layer), transpose to model order (lon,lev,lat),
! then scatter to slaves.
if (nm /= 1 .or. np /= 2) call endrun (subname//': bad nm or np value')
do n=nm,np
if (n == 1) then
mo = mo_prv
else
mo = mo_nxt
end if
do m=1,naer
if (masterproc) then
if (single_column) then
start(:) = (/lonidx,latidx,1,mo/)
else
start(:) = (/1,1,1,mo/)
endif
kount(:) = (/naerlon,naerlat,paerlev,1/)
call handle_ncerr( nf90_get_var (aernid, species_id(m),aerosol_data, start, kount),&
subname, __LINE__)
do j=1,naerlat
do k=1,paerlev
aerosol_field(:,k,j) = aerosol_data(:,j,k)
end do
aerosol_field(:,paerlev+1,j) = 0._r8 ! value at bottom
end do
end if
call scatter_field_to_chunk (1, paerlev+1, 1, naerlon, aerosol_field, &
aerosolc(:,:,:,m,n))
end do
! Retrieve PS from Match
if (masterproc) then
if (single_column) then
start(:) = (/lonidx,latidx,mo,-1/)
else
start(:) = (/1,1,mo,-1/)
endif
kount(:) = (/naerlon,naerlat,1,-1/)
call handle_ncerr( nf90_get_var(aernid, Mpsid, M_ps,start,kount),&
subname, __LINE__)
end if
call scatter_field_to_chunk (1, 1, 1, naerlon, M_ps(:,:), M_ps_cam_col(:,:,n))
end do ! n=nm,np (=1,2)
if(masterproc) deallocate(aerosol_data)
deallocate(aerosol_field)
end if ! Check to see if this dataset is in ncol format.
call t_stopf(subname)
end subroutine cam3_aero_data_init
!================================================================================================
subroutine cam3_aero_data_timestep_init(pbuf, phys_state) 1,12
!------------------------------------------------------------------
!
! Input:
! time at which aerosol masses are needed (get_curr_calday())
! chunk index
! CAM's vertical grid (pint)
!
! Output:
! values for Aerosol Mass at time specified by get_curr_calday
! on vertical grid specified by pint (aer_mass) :: aerosol at time t
!
! Method:
! first determine which indexs of aerosols are the bounding data sets
! interpolate both onto vertical grid aerm(),aerp().
! from those two, interpolate in time.
!
!------------------------------------------------------------------
use interpolate_data, only: get_timeinterp_factors
use phys_buffer, only: pbuf_size_max, pbuf_fld
use cam_logfile, only: iulog
use ppgrid, only: begchunk,endchunk
use physconst, only: gravit
!
! aerosol fields interpolated to current time step
! on pressure levels of this time step.
! these should be made read-only for other modules
! Is allocation done correctly here?
!
type(pbuf_fld), intent(inout), dimension(pbuf_size_max) :: pbuf
type(physics_state), intent(in), dimension(begchunk:endchunk) :: phys_state
!
! Local workspace
!
real(r8) :: pint(pcols,pverp) ! interface pres.
integer :: c ! chunk index
real(r8) caldayloc ! calendar day of current timestep
real(r8) fact1, fact2 ! time interpolation factors
integer i, k, j ! spatial indices
integer m ! constituent index
integer lats(pcols),lons(pcols) ! latitude and longitudes of column
integer ncol ! number of columns
real(r8) speciesmin(naer) ! minimal value for each species
!
! values before current time step "the minus month"
! aerosolm(pcols,pver) is value of preceeding month's aerosol masses
! aerosolp(pcols,pver) is value of next month's aerosol masses
! (think minus and plus or values to left and right of point to be interpolated)
!
real(r8) aerosolm(pcols,pver,naer,begchunk:endchunk) ! aerosol mass from MATCH in column,level at previous (minus) month
!
! values beyond (or at) current time step "the plus month"
!
real(r8) aerosolp(pcols,pver,naer,begchunk:endchunk) ! aerosol mass from MATCH in column,level at next (plus) month
real(r8) :: mass_to_mmr(pcols,pver)
character(len=*), parameter :: subname = 'cam3_aero_data_timestep_init'
logical error_found
!------------------------------------------------------------------
call aerint(phys_state)
caldayloc = get_curr_calday ()
! Determine time interpolation factors. 1st arg says we are cycling 1 year of data
call get_timeinterp_factors (.true., mo_nxt, cdaym, cdayp, caldayloc, &
fact1, fact2, 'GET_AEROSOL:')
! interpolate (prv and nxt month) bounding datasets onto cam vertical grid.
! compute mass mixing ratios on CAMS's pressure coordinate
! for both the "minus" and "plus" months
!
! This loop over chunk could probably be removed by working with the whole
! begchunk:endchunk group at once. It would require a slight generalization
! in vert_interpolate.
do c = begchunk,endchunk
pint = phys_state(c)%pint
ncol = get_ncols_p(c)
call vert_interpolate (M_ps_cam_col(:,c,nm), pint, nm, aerosolm(:,:,:,c), ncol, c)
call vert_interpolate (M_ps_cam_col(:,c,np), pint, np, aerosolp(:,:,:,c), ncol, c)
! Time interpolate.
do m=1,naer
do k=1,pver
do i=1,ncol
aer_mass(i,k,m,c) = aerosolm(i,k,m,c)*fact1 + aerosolp(i,k,m,c)*fact2
end do
end do
! Partition seasalt aerosol mass
if (m .eq. idxSSLTA) then
aer_mass(:ncol,:,m,c) = (1.-wgt_sscm)*aer_mass(:ncol,:,m,c) ! fraction of seasalt mass in accumulation mode
elseif (m .eq. idxSSLTC) then
aer_mass(:ncol,:,m,c) = wgt_sscm*aer_mass(:ncol,:,m,c) ! fraction of seasalt mass in coarse mode
endif
end do
! exit if mass is negative (we have previously set
! cumulative mass to be a decreasing function.)
speciesmin(:) = 0._r8 ! speciesmin(m) = 0 is minimum mass for each species
error_found = .false.
do m=1,naer
do k=1,pver
do i=1,ncol
if (aer_mass(i, k, m,c) < speciesmin(m)) error_found = .true.
end do
end do
end do
if (error_found) then
do m=1,naer
do k=1,pver
do i=1,ncol
if (aer_mass(i, k, m,c) < speciesmin(m)) then
write(iulog,*) subname//': negative mass mixing ratio, exiting'
write(iulog,*) 'm, column, pver',m, i, k ,aer_mass(i, k, m,c)
call endrun ()
end if
end do
end do
end do
end if
do k = 1, pver
mass_to_mmr(1:ncol,k) = gravit/(pint(1:ncol,k+1)-pint(1:ncol,k))
enddo
pbuf(cam3_sul_idx )%fld_ptr(1,1:ncol,:,c,1)=aer_mass(1:ncol,:, idxSUL,c)*mass_to_mmr(:ncol,:)
pbuf(cam3_ssam_idx )%fld_ptr(1,1:ncol,:,c,1)=aer_mass(1:ncol,:, idxSSLTA,c)*mass_to_mmr(:ncol,:)
pbuf(cam3_sscm_idx )%fld_ptr(1,1:ncol,:,c,1)=aer_mass(1:ncol,:, idxSSLTC,c)*mass_to_mmr(:ncol,:)
pbuf(cam3_dust1_idx)%fld_ptr(1,1:ncol,:,c,1)=aer_mass(1:ncol,:, idxDUSTfirst,c)*mass_to_mmr(:ncol,:)
pbuf(cam3_dust2_idx)%fld_ptr(1,1:ncol,:,c,1)=aer_mass(1:ncol,:,idxDUSTfirst+1,c)*mass_to_mmr(:ncol,:)
pbuf(cam3_dust3_idx)%fld_ptr(1,1:ncol,:,c,1)=aer_mass(1:ncol,:,idxDUSTfirst+2,c)*mass_to_mmr(:ncol,:)
pbuf(cam3_dust4_idx)%fld_ptr(1,1:ncol,:,c,1)=aer_mass(1:ncol,:,idxDUSTfirst+3,c)*mass_to_mmr(:ncol,:)
pbuf(cam3_ocpho_idx)%fld_ptr(1,1:ncol,:,c,1)=aer_mass(1:ncol,:, idxOCPHO,c)*mass_to_mmr(:ncol,:)
pbuf(cam3_bcpho_idx)%fld_ptr(1,1:ncol,:,c,1)=aer_mass(1:ncol,:, idxBCPHO,c)*mass_to_mmr(:ncol,:)
pbuf(cam3_ocphi_idx)%fld_ptr(1,1:ncol,:,c,1)=aer_mass(1:ncol,:, idxOCPHI,c)*mass_to_mmr(:ncol,:)
pbuf(cam3_bcphi_idx)%fld_ptr(1,1:ncol,:,c,1)=aer_mass(1:ncol,:, idxBCPHI,c)*mass_to_mmr(:ncol,:)
enddo ! c = begchunk:endchunk
end subroutine cam3_aero_data_timestep_init
!================================================================================================
subroutine vert_interpolate (Match_ps, pint, n, aerosol_mass, ncol, c) 2,2
!--------------------------------------------------------------------
! Input: match surface pressure, cam interface pressure,
! month index, number of columns, chunk index
!
! Output: Aerosol mass mixing ratio (aerosol_mass)
!
! Method:
! interpolate column mass (cumulative) from match onto
! cam's vertical grid (pressure coordinate)
! convert back to mass mixing ratio
!
!--------------------------------------------------------------------
real(r8), intent(out) :: aerosol_mass(pcols,pver,naer) ! aerosol mass from MATCH
real(r8), intent(in) :: Match_ps(pcols) ! surface pressure at a particular month
real(r8), intent(in) :: pint(pcols,pverp) ! interface pressure from CAM
integer, intent(in) :: ncol,c ! chunk index and number of columns
integer, intent(in) :: n ! prv or nxt month index
!
! Local workspace
!
integer m ! index to aerosol species
integer kupper(pcols) ! last upper bound for interpolation
integer i, k, kk, kkstart, kount ! loop vars for interpolation
integer isv, ksv, msv ! loop indices to save
logical bad ! indicates a bad point found
logical lev_interp_comp ! interpolation completed for a level
logical error_found
real(r8) aerosol(pcols,pverp,naer) ! cumulative mass of aerosol in column beneath upper
! interface of level in column at particular month
real(r8) dpl, dpu ! lower and upper intepolation factors
real(r8) v_coord ! vertical coordinate
real(r8) AER_diff ! temp var for difference between aerosol masses
character(len=*), parameter :: subname = 'cam3_aero_data.vert_interpolate'
!-----------------------------------------------------------------------
call t_startf ('vert_interpolate')
!
! Initialize index array
!
do i=1,ncol
kupper(i) = 1
end do
!
! assign total mass to topmost level
!
aerosol(:,1,:) = aerosolc(:,1,c,:,n)
!
! At every pressure level, interpolate onto that pressure level
!
do k=2,pver
!
! Top level we need to start looking is the top level for the previous k
! for all longitude points
!
kkstart = paerlev+1
do i=1,ncol
kkstart = min0(kkstart,kupper(i))
end do
kount = 0
!
! Store level indices for interpolation
!
! for the pressure interpolation should be comparing
! pint(column,lev) with M_hybi(lev)*M_ps_cam_col(month,column,chunk)
!
lev_interp_comp = .false.
do kk=kkstart,paerlev
if(.not.lev_interp_comp) then
do i=1,ncol
v_coord = pint(i,k)
if (M_hybi(kk)*Match_ps(i) .lt. v_coord .and. v_coord .le. M_hybi(kk+1)*Match_ps(i)) then
kupper(i) = kk
kount = kount + 1
end if
end do
!
! If all indices for this level have been found, do the interpolation and
! go to the next level
!
! Interpolate in pressure.
!
if (kount.eq.ncol) then
do m=1,naer
do i=1,ncol
dpu = pint(i,k) - M_hybi(kupper(i))*Match_ps(i)
dpl = M_hybi(kupper(i)+1)*Match_ps(i) - pint(i,k)
aerosol(i,k,m) = &
(aerosolc(i,kupper(i) ,c,m,n)*dpl + &
aerosolc(i,kupper(i)+1,c,m,n)*dpu)/(dpl + dpu)
enddo !i
end do
lev_interp_comp = .true.
end if
end if
end do
!
! If we've fallen through the kk=1,levsiz-1 loop, we cannot interpolate and
! must extrapolate from the bottom or top pressure level for at least some
! of the longitude points.
!
if(.not.lev_interp_comp) then
do m=1,naer
do i=1,ncol
if (pint(i,k) .lt. M_hybi(1)*Match_ps(i)) then
aerosol(i,k,m) = aerosolc(i,1,c,m,n)
else if (pint(i,k) .gt. M_hybi(paerlev+1)*Match_ps(i)) then
aerosol(i,k,m) = 0.0_r8
else
dpu = pint(i,k) - M_hybi(kupper(i))*Match_ps(i)
dpl = M_hybi(kupper(i)+1)*Match_ps(i) - pint(i,k)
aerosol(i,k,m) = &
(aerosolc(i,kupper(i) ,c,m,n)*dpl + &
aerosolc(i,kupper(i)+1,c,m,n)*dpu)/(dpl + dpu)
end if
end do
end do
if (kount.gt.ncol) then
call endrun (subname//': Bad data: non-monotonicity suspected in dependent variable')
end if
end if
end do
! call t_startf ('vi_checks')
!
! aerosol mass beneath lowest interface (pverp) must be 0
!
aerosol(1:ncol,pverp,:) = 0._r8
!
! Set mass in layer to zero whenever it is less than
! 1.e-40 kg/m^2 in the layer
!
do m = 1, naer
do k = 1, pver
do i = 1, ncol
if (aerosol(i,k,m) < 1.e-40_r8) aerosol(i,k,m) = 0._r8
end do
end do
end do
!
! Set mass in layer to zero whenever it is less than
! 10^-15 relative to column total mass
!
error_found = .false.
do m = 1, naer
do k = 1, pver
do i = 1, ncol
AER_diff = aerosol(i,k,m) - aerosol(i,k+1,m)
if( abs(AER_diff) < 1e-15_r8*aerosol(i,1,m)) then
AER_diff = 0._r8
end if
aerosol_mass(i,k,m)= AER_diff
if (aerosol_mass(i,k,m) < 0) error_found = .true.
end do
end do
end do
if (error_found) then
do m = 1, naer
do k = 1, pver
do i = 1, ncol
if (aerosol_mass(i,k,m) < 0) then
write(iulog,*) subname//': mass < 0, m, col, lev, mass',m, i, k, aerosol_mass(i,k,m)
write(iulog,*) subname//': aerosol(k),(k+1)',aerosol(i,k,m),aerosol(i,k+1,m)
write(iulog,*) subname//': pint(k+1),(k)',pint(i,k+1),pint(i,k)
write(iulog,*)'n,c',n,c
call endrun()
end if
end do
end do
end do
end if
call t_stopf ('vert_interpolate')
return
end subroutine vert_interpolate
!================================================================================================
subroutine aerint (phys_state) 1,9
type(physics_state), intent(in) :: phys_state(begchunk:endchunk)
integer :: ntmp ! used in index swapping
integer :: start(4) ! start vector for netcdf calls
integer :: kount(4) ! count vector for netcdf calls
integer :: i,j,k ! spatial indices
integer :: m ! constituent index
integer :: cols, cole
integer :: lchnk, ncol
real(r8) :: caldayloc ! calendar day of current timestep
real(r8) :: aerosol_data(naerlon,naerlat,paerlev) ! aerosol field read in from dataset
real(r8) :: aerosol_field(naerlon,paerlev+1,naerlat) ! aerosol field to be scattered
integer latidx,lonidx
real(r8) closelat,closelon
character(len=*), parameter :: subname = 'cam3_aero_data.aerint'
!-----------------------------------------------------------------------
if (single_column) &
call shr_scam_GetCloseLatLon(aernid,scmlat,scmlon,closelat,closelon,latidx,lonidx)
!
! determine if need to read in next month data
! also determine time interpolation factors
!
caldayloc = get_curr_calday ()
!
! If model time is past current forward timeslice, then
! masterproc reads in the next timeslice for time interpolation. Messy logic is
! for interpolation between December and January (mo_nxt == 1). Just like
! ozone_data_timestep_init, sstint.
!
if (caldayloc > cdayp .and. .not. (mo_nxt == 1 .and. caldayloc >= cdaym)) then
mo_nxt = mod(mo_nxt,12) + 1
cdaym = cdayp
cdayp = Mid(mo_nxt)
!
! Check for valid date info
!
if (.not. (mo_nxt == 1 .or. caldayloc <= cdayp)) then
call endrun (subname//': Non-monotonicity suspected in input aerosol data')
end if
ntmp = nm
nm = np
np = ntmp
if(aerosol_datan%isncol) then
do lchnk=begchunk,endchunk
ncol=phys_state(lchnk)%ncol
cols=1
cole=cols+aerosol_datan%count(cols,lchnk)-1
do while(cole<=ncol)
start=(/aerosol_datan%start(cols,lchnk),mo_nxt,1,-1/)
kount=(/aerosol_datan%count(cols,lchnk),1,-1,-1/)
call handle_ncerr( nf90_get_var(aerosol_datan%ncid, aerosol_datan%psid , &
aerosol_datan%ps(cols:cole,lchnk,np), start(1:2), &
kount(1:2)),&
subname, __LINE__)
start(2)=1
start(3)=mo_nxt
kount(2)=paerlev
kount(3)=1
do m=1,naer
call handle_ncerr( nf90_get_var(aerosol_datan%ncid, aerosol_datan%dataid(m) , &
aerosol_datan%fields(cols:cole,:,lchnk,m,np), &
start(1:3), kount(1:3)),&
subname, __LINE__)
end do
if(cols==ncol) exit
cols=cols+aerosol_datan%count(cols,lchnk)
cole=cols+aerosol_datan%count(cols,lchnk)-1
end do
end do
aerosolc(:,1:paerlev,:,:,np)=aerosol_datan%fields(:,:,:,:,np)
else
do m=1,naer
if (masterproc) then
if (single_column) then
naerlon=1
naerlat=1
start(:) = (/lonidx,latidx,1,mo_nxt/)
else
start(:) = (/1,1,1,mo_nxt/)
endif
kount(:) = (/naerlon,naerlat,paerlev,1/)
call handle_ncerr( nf90_get_var (aernid, species_id(m), aerosol_data, start, kount),&
subname, __LINE__)
do j=1,naerlat
do k=1,paerlev
aerosol_field(:,k,j) = aerosol_data(:,j,k)
end do
aerosol_field(:,paerlev+1,j) = 0._r8 ! value at bottom
end do
end if
call scatter_field_to_chunk (1, paerlev+1, 1, naerlon, aerosol_field, &
aerosolc(:,:,:,m,np))
end do
!
! Retrieve PS from Match
!
if (masterproc) then
if (single_column) then
naerlon=1
naerlat=1
start(:) = (/lonidx,latidx,mo_nxt,-1/)
else
start(:) = (/1,1,mo_nxt,-1/)
endif
kount(:) = (/naerlon,naerlat,1,-1/)
call handle_ncerr( nf90_get_var (aernid, Mpsid, M_ps, start, kount),&
subname, __LINE__)
write(iulog,*) subname//': Read aerosols data for julian day', Mid(mo_nxt)
end if
call scatter_field_to_chunk (1, 1, 1, naerlon, M_ps(:,:), M_ps_cam_col(:,:,np))
end if
end if
end subroutine aerint
end module cam3_aero_data