#include <misc.h>
#include <params.h>
module volcanicmass 3,8
!-----------------------------------------------------------------------
!
! Purpose: read, store, interpolate, and retrieve mass fields describing
! volcanic emissions
!
! Author: A. Conley, but mostly copied from ozone code base.
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid, only: plat, plon, plond, masterproc, plev, plevp
use commap, only: latdeg
use ppgrid, only: begchunk, endchunk, pcols, pver, pverp
use phys_grid, only: scatter_field_to_chunk, get_ncols_p
use filenames, only: bndtvvolc
use abortutils, only: endrun
use timeinterp, only: validfactors
#if ( defined SPMD )
use mpishorthand
#endif
implicit none
private
save
public read_volcanic_mass
public get_volcanic_mass
public volcanic_initialize
real(r8) cdayprev ! dataset calendar day previous month
real(r8) cdaynext ! dataset calendar day next month
integer nprev ! Array indices for previous month volcanic data
integer nnext ! Array indices for next month volcanic data
integer ntmp ! temp storage for swaping nprev, nnext
integer massid ! netcdf id for volcanic variable
integer levsiz ! size of level dimension on volcanic dataset
integer latsiz ! size of latitude dimension on volcanic dataset
integer timesiz ! size of time dimension on volcanic dataset
integer np1 ! current forward time index of volcanic dataset
integer ncid_volc ! netcdf file id for volcanic file
type volcanic_pointers
real(r8), dimension(:,:,:), pointer :: mass ! (pcols,levsiz,begchunk:endchunk)
end type volcanic_pointers
type (volcanic_pointers) :: volcanic(2) ! pointers to volcanic masses in 2 bounding months
type (volcanic_pointers) :: volcanic2(2) ! pointers to volcanic masses in 2 bounding months on CAM grid
real(r8), allocatable :: volcanic_mass(:,:,:)! (pcols,plev,begchunk:endchunk)
real(r8), allocatable :: pin(:) ! volcanic pressure level (levsiz)
real(r8), allocatable :: volclat(:) ! Latitudes of volcanic dataset (latsiz)
integer, allocatable :: date_volc(:) ! Date on volcanic dataset (YYYYMMDD)
! (timesiz)
integer, allocatable :: sec_volc(:) ! seconds of date (0-86399)
! (timesiz)
contains
subroutine volcanic_initialize 1,45
!-----------------------------------------------------------------------
!
! Purpose: Read appropriate portion of time-variant volcanic boundary
! dataset, containing volcanic mixing ratios as a
! function of latitude and pressure. Read two
! consecutive months between which the current date lies.
!
! Another routine
! then evaluates the two path length integrals (with and without
! pressure weighting) from zero to the interfaces between the input
! levels. It also stores the contribution to the integral from each
! layer.
!
! Method: Call appropriate netcdf wrapper routines and interpolate to model grid
!
! Author: CCM Core Group
! Modified: P. Worley, August 2003, for chunking and performance optimization
! Modified: A. Conley, Sept 2003, for volcanic aerosols
!
!-----------------------------------------------------------------------
use commap, only: latdeg, londeg
use rgrid, only: nlon
use time_manager, only: get_curr_date, get_perp_date, get_curr_calday, &
is_perpetual
use ioFileMod, only: getfil
!-----------------------------------------------------------------------
include 'netcdf.inc'
!-----------------------------------------------------------------------
!
! Local workspace
!
integer dateid ! netcdf id for date variable
integer secid ! netcdf id for seconds variable
integer latdimid ! netcdf id for latitude dimension
integer levdimid ! netcdf id for level dimension
integer bindimid ! netcdf id for bin dimension
integer latid ! netcdf id for latitude variable
integer levid ! netcdf id for level variable
integer dimids(nf_max_var_dims) ! variable shape
integer dim_cnt(4) ! array of counts for each dimension
integer dim_strt(4) ! array of starting indices
integer ilon, ilev, ilat, time ! longitude, level, latitude, time indices
integer ichunk ! index into chunk
integer :: yr, mon, day ! components of a date
integer :: ncdate ! current date in integer format [yyyymmdd]
integer :: ncsec ! current time of day [seconds]
real(r8) :: calday ! current calendar day
real(r8) caldayloc ! calendar day (includes yr)
real(r8), allocatable :: volc2D(:,:) ! temporary volcanic arrays
real(r8), allocatable :: volc3D(:,:,:)
real(r8), allocatable :: volcbdyprev(:,:,:,:) ! volcanic data previous time sample
real(r8), allocatable :: volcbdynext(:,:,:,:) ! volcanic data next time sample
character(len=256) :: locfn ! netcdf local filename to open
!
! find and open file; abort if fail (getfil(,,0)).
!
nprev = 1
nnext = 2
!
! SPMD: Master does all the work. Sends needed info to slaves
!
if (masterproc) then
calday = get_curr_calday()
call get_curr_date(yr, mon, day, ncsec)
ncdate = yr*10000 + mon*100 + day
caldayloc = calday + yr*365.
! call bnddyi(ncdate, 0, caldayloc) !!!!
!
! Get and check dimension info
!
call getfil(bndtvvolc, locfn, 0)
call wrap_open(locfn, 0, ncid_volc)
CALL WRAP_INQ_DIMID( ncid_volc, 'lev', levdimid )
CALL WRAP_INQ_DIMID( ncid_volc, 'lat', latdimid )
CALL WRAP_INQ_DIMID( ncid_volc, 'bin', bindimid )
CALL WRAP_INQ_DIMID( ncid_volc, 'date', dateid )
CALL WRAP_INQ_DIMLEN( ncid_volc, levdimid, levsiz )
CALL WRAP_INQ_DIMLEN( ncid_volc, latdimid, latsiz )
CALL WRAP_INQ_DIMLEN( ncid_volc, dateid, timesiz )
! CALL WRAP_INQ_VARID( ncid_volc, 'date', dateid )
! CALL WRAP_INQ_VARID( ncid_volc, 'datesec', secid )
CALL WRAP_INQ_VARID( ncid_volc, 'MVOLC', massid )
CALL WRAP_INQ_VARID( ncid_volc, 'lat', latid )
CALL WRAP_INQ_VARID( ncid_volc, 'lev', levid )
CALL WRAP_INQ_VARDIMID (ncid_volc, massid, dimids)
if (dimids(1) /= levdimid .or. dimids(2) /= latdimid .or. dimids(3) /= bindimid .or. dimids(4) /= dateid ) then
call endrun ('VOLCANIC_INITIALIZE: Data must be ordered lev, lat, time')
endif
endif
#if (defined SPMD )
call mpibcast( latsiz, 1, mpiint, 0, mpicom )
call mpibcast( levsiz, 1, mpiint, 0, mpicom )
call mpibcast( timesiz, 1, mpiint, 0, mpicom )
#endif
!
! Dynamically allocated memory for module
!
allocate (date_volc(timesiz))
allocate (sec_volc(timesiz))
allocate (pin(levsiz))
allocate (volcanic(nprev)%mass(pcols,levsiz,begchunk:endchunk))
allocate (volcanic(nnext)%mass(pcols,levsiz,begchunk:endchunk))
allocate (volcanic_mass(pcols,levsiz,begchunk:endchunk))
!
! Locally dynamic that will be deallocated before "return"
!
allocate (volc3D(plond,levsiz,plat))
if (masterproc) then
!
! More dynamically allocated memory for module
! (for masterproc only)
!
allocate (volclat(latsiz))
!
! More locally dynamic that will be deallocated before return
!
allocate (volcbdyprev(levsiz,latsiz,1,1))
allocate (volcbdynext(levsiz,latsiz,1,1))
allocate (volc2D(levsiz,plat))
!
! Retrieve longitude, latitude and level arrays for interpolation.
!
call wrap_get_var_realx (ncid_volc, latid, volclat)
call wrap_get_var_realx (ncid_volc, levid, pin)
!
! Retrieve entire date and sec variables.
!
call wrap_get_var_int (ncid_volc,dateid,date_volc)
! call wrap_get_var_int (ncid_volc,secid,sec_volc)
sec_volc = 0
dim_strt(1) = 1
dim_strt(2) = 1
dim_strt(3) = 1
dim_cnt(1) = levsiz
dim_cnt(2) = latsiz
dim_cnt(3) = 1 ! only 1 bin
dim_cnt(4) = 1
!
! Normal interpolation between consecutive time slices.
!
do time=1,timesiz-1
np1 = time + 1
call bnddyi(date_volc(time), sec_volc(time), cdayprev)
call bnddyi(date_volc(np1 ), sec_volc(np1 ), cdaynext)
yr = date_volc(time)/10000
cdayprev = cdayprev + yr*365.
yr = date_volc(np1)/10000
cdaynext = cdaynext + yr*365.
if (caldayloc > cdayprev .and. caldayloc <= cdaynext) then
dim_strt(4) = time
call wrap_get_vara_realx (ncid_volc,massid,dim_strt,dim_cnt,volcbdyprev)
dim_strt(4) = np1
call wrap_get_vara_realx (ncid_volc,massid,dim_strt,dim_cnt,volcbdynext)
goto 10
endif
end do
write(6,*)'volcanic_initialize: Failed to find dates bracketing ncdate, ncsec=', ncdate, ncsec
call endrun
10 continue
write(6,*)'volcanic_initialize: Read volcanic data for dates ',date_volc(time), &
sec_volc(time),' and ',date_volc(np1),sec_volc(np1)
!
! Latitude interpolation for "prev" data. Expand to 3D after interpolation.
!
call lininterp (volcbdyprev(:,:,1,1),volclat ,levsiz ,latsiz ,volc2D(:,:), &
latdeg ,plat )
deallocate (volcbdyprev)
do ilat=1,plat
do ilev=1,levsiz
do ilon=1,nlon(ilat)
volc3D(ilon,ilev,ilat) = volc2D(ilev,ilat)
end do
end do
end do
call scatter_field_to_chunk(1,levsiz,1,plond,volc3D,volcanic(nprev)%mass)
#if (defined SPMD )
else
call scatter_field_to_chunk(1,levsiz,1,plond,volc3D,volcanic(nprev)%mass)
endif
if (masterproc) then
#endif
!
! Latitude interpolation for "next" data. Expand to 3D after interpolation.
!
call lininterp (volcbdynext(:,:,1,1),volclat ,levsiz ,latsiz ,volc2D(:,:), &
latdeg ,plat )
deallocate (volcbdynext)
do ilat=1,plat
do ilev=1,levsiz
do ilon=1,nlon(ilat)
volc3D(ilon,ilev,ilat) = volc2D(ilev,ilat)
end do
end do
end do
deallocate (volc2D)
call scatter_field_to_chunk(1,levsiz,1,plond,volc3D,volcanic(nnext)%mass)
#if (defined SPMD )
else
call scatter_field_to_chunk(1,levsiz,1,plond,volc3D,volcanic(nnext)%mass)
#endif
endif
!
! Deallocate dynamic memory for local workspace. NOT for pointers in common.
!
deallocate (volc3D)
#if (defined SPMD )
call mpibcast (np1, 1, mpiint, 0, mpicom)
call mpibcast (pin, levsiz, mpir8, 0, mpicom)
call mpibcast (date_volc, timesiz, mpiint, 0, mpicom )
call mpibcast (sec_volc, timesiz, mpiint, 0, mpicom )
call mpibcast (cdayprev, 1, mpir8, 0, mpicom)
call mpibcast (cdaynext, 1, mpir8, 0, mpicom)
#endif
!
! accumulate masses from surface before interpolating
! allocate memory for boundary datasets on cam vertical grid
! interpolate from file levels to cam levels
! deallocate memory for boundary dataset (previous) on file levels
!
do ichunk = begchunk, endchunk
do ilev = levsiz-1,1,-1
volcanic(nprev)%mass(:,ilev,ichunk)=volcanic(nprev)%mass(:,ilev,ichunk) &
+ volcanic(nprev)%mass(:,ilev+1,ichunk)
volcanic(nnext)%mass(:,ilev,ichunk)=volcanic(nnext)%mass(:,ilev,ichunk) &
+ volcanic(nnext)%mass(:,ilev+1,ichunk)
enddo
enddo
allocate (volcanic2(nprev)%mass(pcols,plevp,begchunk:endchunk))
allocate (volcanic2(nnext)%mass(pcols,plevp,begchunk:endchunk))
call vert_interpolate(nprev)
call vert_interpolate(nnext)
! deallocate (volcanic(nprev)%mass)
return
end subroutine volcanic_initialize
subroutine read_volcanic_mass() 1,12
use rgrid, only: nlon
use time_manager, only: get_curr_date, get_curr_calday
integer :: dim_cnt(4) ! array of counts for each dimension
integer :: dim_strt(4) ! array of starting indices
integer :: ilon,ilev,ilat,icol,ncols ! indices
integer :: lchnk ! chunk index
integer :: yr,mon,day ! components of date
integer :: ncdate ! formated date [yyyymmdd]
integer :: ncsec ! seconds past midnight of current ncdate
real(r8):: calday ! current calendar day
real(r8):: caldayloc ! calendar day (including yr)
real(r8):: deltat ! time (days) between times of dataset
real(r8), allocatable :: volc2D(:,:) ! temporary volcanic arrays
real(r8), allocatable :: volc3D(:,:,:) ! temporary volcanic arrays
real(r8), allocatable :: volcbdynext(:,:,:,:) ! volcanic data next time sample
include 'netcdf.inc'
call get_curr_date(yr, mon, day, ncsec)
ncdate = yr*10000 + mon*100 + day
calday = get_curr_calday()
caldayloc = calday + yr*365.
if (caldayloc > cdaynext) then
write(6,*) 'caldayloc',caldayloc,'cdaynext',cdaynext
np1 = np1 + 1
if (np1 > timesiz) then
call endrun ('READ_VOLCANIC_MASS: Attempt to read past end of Volcanic dataset')
end if
cdayprev = cdaynext
call bnddyi(date_volc(np1), sec_volc(np1), cdaynext)
yr = date_volc(np1)/10000
cdaynext = cdaynext + yr*365.
if(caldayloc <= cdaynext) then
ntmp = nprev
nprev = nnext
nnext = ntmp
!
! Allocate memory for dynamic local workspace
!
allocate (volc3D(plond,levsiz,plat))
if (masterproc) then
dim_strt(1) = 1
dim_strt(2) = 1
dim_strt(3) = 1
dim_strt(4) = np1
dim_cnt(1) = levsiz
dim_cnt(2) = latsiz
dim_cnt(3) = 1 ! only 1 bin
dim_cnt(4) = 1
!
! Allocate memory for more dynamic local workspace
!
allocate (volcbdynext(levsiz,latsiz,1,1))
allocate (volc2D(levsiz,plat))
call wrap_get_vara_realx(ncid_volc,massid,dim_strt,dim_cnt,volcbdynext)
write(6,*)'get_volcanic_mass: read volcanic masses for date (yyyymmdd) ', date_volc(np1), ' sec ',sec_volc(np1)
!
! Spatial interpolation.
!
call lininterp (volcbdynext(:,:,1,1),volclat ,levsiz ,latsiz ,volc2D(:,:), &
latdeg ,plat )
do ilat=1,plat
do ilev=1,levsiz
do ilon=1,nlon(ilat)
volc3D(ilon,ilev,ilat) = volc2D(ilev,ilat)
end do
end do
end do
call scatter_field_to_chunk(1,levsiz,1,plond,volc3D,volcanic(nnext)%mass)
deallocate (volcbdynext)
deallocate (volc2D)
deallocate (volc3D)
else ! masterproc
call scatter_field_to_chunk(1,levsiz,1,plond,volc3D,volcanic(nnext)%mass)
deallocate (volc3D)
endif ! masterproc
do lchnk = begchunk,endchunk
do ilev = levsiz-1,1,-1
volcanic(nnext)%mass(:,ilev,lchnk)=volcanic(nnext)%mass(:,ilev,lchnk) &
+ volcanic(nnext)%mass(:,ilev+1,lchnk)
enddo
enddo
call vert_interpolate(nnext)
else ! if(caldayloc <= cdaynext)
if (masterproc) then
write(6,*)'READ_VOLCANIC_MASS: date from dataset ',date_volc(np1),' does not exceed model date ', ncdate
endif
call endrun
endif ! if(caldayloc <= cdaynext)
endif ! (caldayloc > cdaynext)
end subroutine read_volcanic_mass
subroutine get_volcanic_mass(lchnk, aerosol) 1,6
use rgrid, only: nlon
use time_manager, only: get_curr_date, get_curr_calday
integer, intent(in) :: lchnk ! chunk id for local aerosol array
real(r8), intent(inout) :: aerosol(:,:) ! aerosol array to be filled
integer :: dim_cnt(4) ! array of counts for each dimension
integer :: dim_strt(4) ! array of starting indices
integer :: ilon,ilev,ilat,icol,ncols ! indices
integer :: yr,mon,day ! components of date
integer :: ncdate ! formated date [yyyymmdd]
integer :: ncsec ! seconds past midnight of current ncdate
real(r8):: fact1, fact2 ! time interpolation factors
real(r8):: calday ! current calendar day
real(r8):: caldayloc ! calendar day (including yr)
real(r8):: deltat ! time (days) between times of dataset
real(r8), allocatable :: volc2D(:,:) ! temporary volcanic arrays
real(r8), allocatable :: volc3D(:,:,:) ! temporary volcanic arrays
real(r8), allocatable :: volcbdynext(:,:,:,:) ! volcanic data next time sample
!
! time interpolate
!
call get_curr_date(yr, mon, day, ncsec)
calday = get_curr_calday()
caldayloc = calday + yr*365.
deltat = cdaynext - cdayprev
fact1 = (cdaynext - caldayloc)/deltat
fact2 = (caldayloc - cdayprev)/deltat
!
! Check sanity of time interpolation calculation to within 32-bit roundoff
!
if (.not. validfactors (fact1, fact2)) then
write(6,*)'get_volcanic_mass: Bad fact1 and/or fact2=',fact1,fact2
call endrun ()
end if
!
! Interpolate (in time)!
!
ncols = get_ncols_p(lchnk)
do ilev=1,plev
do icol=1,ncols
aerosol(icol,ilev) = volcanic2(nprev)%mass(icol,ilev,lchnk)*fact1 + volcanic2(nnext)%mass(icol,ilev,lchnk)*fact2
enddo
enddo
return
end subroutine get_volcanic_mass
subroutine vert_interpolate(month) 5,6
!-----------------------------------------------------------------------------
!
! volcanic2(month) is assigned vertically interpolated values of
! volcanic(month). volcanic is on vertical grid defined by file.
! volcanic2 is on vertical grid of CAM simulation.
!
!-----------------------------------------------------------------------------
use pmgrid, only: plat, plon, plond, masterproc, plev, plevp
#include <comhyb.h>
integer, intent(in):: month ! which month boundary should be interpolated?
real(r8) :: model_level(pverp)
real(r8) :: AER_diff(pverp)
integer :: ncols ! number of columns in chunk
integer :: lchnk ! chunk index
integer :: kk ! boundary data index into file data
integer :: ilat, ilon, ilev, icol ! indices
real(r8) :: dpl, dpu ! delta p across lower and upper parts of interpolant
model_level(:) = 1000._r8 *(hyai(:) + hybi(:))
do lchnk = begchunk, endchunk
ncols = get_ncols_p(lchnk)
volcanic2(month)%mass(1:ncols,1,lchnk) = volcanic(month)%mass(1:ncols,1,lchnk)
!
! At every pressure level, interpolate onto that pressure level
!
kk = 2
do ilev = 2, pver
!
! 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(masterproc) then
! write(6,*) 'model_level:',model_level
! write(6,*) 'pin:',pin
! endif
if (model_level(ilev) .lt. pin(1)) then
volcanic2(month)%mass(1:ncols,ilev,lchnk) = volcanic(month)%mass(1:ncols,1,lchnk)
else if (model_level(ilev) .ge. pin(levsiz)) then
volcanic2(month)%mass(1:ncols,ilev,lchnk) = 0._r8
else
do while (model_level(ilev) .ge. pin(kk) )
kk = kk + 1
enddo
dpu = model_level(ilev) - pin(kk-1)
dpl = pin(kk) - model_level(ilev)
volcanic2(month)%mass(1:ncols,ilev,lchnk) = &
(volcanic(month)%mass(1:ncols,kk-1,lchnk)*dpl + &
volcanic(month)%mass(1:ncols,kk ,lchnk)*dpu)/(dpl + dpu)
end if
end do
! do ilev = 1, pver
! do icol = 1, ncols
! if(volcanic2(month)%mass(icol,ilev,lchnk) .lt. 0) then
! write(6,*)'vert interp error'
! write(6,*)'lchnk,ilev,icol',lchnk,ilev,icol
! write(6,*)'volcanic', volcanic(month)%mass(icol,:,lchnk)
! write(6,*)'volcanic2',volcanic2(month)%mass(icol,:,lchnk)
! call endrun
! endif
! enddo
! enddo
!
! aerosol mass beneath lowest interface (pverp) must be 0
!
volcanic2(month)%mass(1:ncols, pverp, lchnk)=0._r8
!
! Set mass in layer to zero whenever it is less than
! 1.e-40 kg/m^2 in the layer
!
! do ilev = 1, pver
! do icol = 1, ncols
! if(volcanic2(month)%mass(icol,ilev,lchnk) < 1.e-40) volcanic2(month)%mass(icol,ilev,lchnk) = 0._r8
! enddo
! enddo
!!
!! Extract mass per layer from cumulative mass.
!!
do ilev = 1, pver
do icol = 1, ncols
volcanic2(month)%mass(icol,ilev,lchnk) = volcanic2(month)%mass(icol,ilev,lchnk)-volcanic2(month)%mass(icol,ilev+1,lchnk)
if(volcanic2(month)%mass(icol,ilev,lchnk) .lt. 0._r8) then
volcanic2(month)%mass(icol,ilev,lchnk) = 0._r8
endif
enddo
enddo
do ilev = 1, pver
do icol = 1, ncols
if(volcanic2(month)%mass(icol,ilev,lchnk) .lt. 0) then
write(6,*)'mass per layer error'
write(6,*)'lchnk,ilev,icol',lchnk,ilev,icol
write(6,*)'volcanic', volcanic(month)%mass(icol,:,lchnk)
write(6,*)'volcanic2',volcanic2(month)%mass(icol,:,lchnk)
call endrun ('VERT_INTERPOLATE')
endif
enddo
enddo
enddo ! lchnk
return
end subroutine vert_interpolate
end module volcanicmass