#include <misc.h>
#include <params.h>
module chemistry 8,8
!---------------------------------------------------------------------------------
! Module to parameterized greenhouse gas chemical loss frequencies from
! Portmann and Solomon
!---------------------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid, only: plat, plev, plevp, plon, masterproc
use ppgrid, only: pcols, pver
use physconst, only: mwdry, mwch4, mwn2o, mwf11, mwf12, mwh2o, gravit
use constituents, only: ppcnst, cnst_add, cnst_name, advected
use ghg_surfvals, only: ch4vmr, n2ovmr, f11vmr, f12vmr
use abortutils, only: endrun
use timeinterp, only: getfactors
implicit none
private
save
!
! Public interfaces
!
public chem_register ! register consituents
public chem_implements_cnst ! returns true if consituent is implemented by this package
public chem_init_cnst ! initialize mixing ratios if not read from initial file
public chem_init ! initialize (history) variables
public chem_timestep_init ! time interpolate chemical loss frequencies
public chem_timestep_tend ! interface to tendency computation
! Namelist variable
logical, public :: trace_gas=.false. ! set true to activate this module
integer, parameter :: ptrlon=01 ! number of longitudes in input dataset
integer, parameter :: ptrlat=36 ! number of latitudes in input dataset
integer, parameter :: ptrlev=56 ! number of levels in input dataset
integer, parameter :: ptrtim=12 ! number of times(months) in input dataset
! Ratios of molecular weights
real(r8), parameter :: rmwn2o = mwn2o/mwdry ! ratio of molecular weight n2o to dry air
real(r8), parameter :: rmwch4 = mwch4/mwdry ! ratio of molecular weight ch4 to dry air
real(r8), parameter :: rmwf11 = mwf11/mwdry ! ratio of molecular weight cfc11 to dry air
real(r8), parameter :: rmwf12 = mwf12/mwdry ! ratio of molecular weight cfc12 to dry air
real(r8), parameter :: rh2och4= mwh2o/mwch4 ! ratio of molecular weight h2o to ch4
! Time/space interpolation of loss frequencies
real(r8) :: tch4i (plat,plev,ptrtim) ! input data ch4 loss rate interp. in lat and lev
real(r8) :: tn2oi (plat,plev,ptrtim) ! input data n2o loss rate interp. in lat and lev
real(r8) :: tcfc11i(plat,plev,ptrtim) ! input data cfc11 loss rate interp. in lat and lev
real(r8) :: tcfc12i(plat,plev,ptrtim) ! input data cfc12 loss rate interp. in lat and lev
real(r8) :: tch4m (plat,plev,2) ! input data ch4 loss rate interp. in lat and lev
real(r8) :: tn2om (plat,plev,2) ! input data n2o loss rate interp. in lat and lev
real(r8) :: tcfc11m(plat,plev,2) ! input data cfc11 loss rate interp. in lat and lev
real(r8) :: tcfc12m(plat,plev,2) ! input data cfc12 loss rate interp. in lat and lev
real(r8) :: tch4 (plat,plev) ! instantaneous ch4 loss rate
real(r8) :: tn2o (plat,plev) ! instantaneous ch4 loss rate
real(r8) :: tcfc11 (plat,plev) ! instantaneous ch4 loss rate
real(r8) :: tcfc12 (plat,plev) ! instantaneous ch4 loss rate
real(r8) :: cdaytrm ! calendar day for previous month data
real(r8) :: cdaytrp ! calendar day for next month data
integer :: np ! array index for previous month tracer data
integer :: nm ! array index for next month tracer data
integer :: np1 ! current forward time index of tracer dataset
integer :: date_tr(ptrtim) ! date on tracer dataset (YYYYMMDD)
integer :: sec_tr(ptrtim) ! seconds of date on tracer dataset (0-86399)
! dummy values for specific heats at constant pressure
real(r8), parameter:: cpch4 = 666.
real(r8), parameter:: cpn2o = 666.
real(r8), parameter:: cpf11 = 666.
real(r8), parameter:: cpf12 = 666.
integer, parameter :: ncnst=4 ! number of constituents
character(len=8), dimension(ncnst), parameter :: & ! constituent names
cnst_names = (/'N2O ', 'CH4 ', 'CFC11', 'CFC12'/)
character(len=8) :: srcnam(ncnst) ! names of source/sink tendencies
integer :: ixghg ! index of 1st constituent (N2O)
contains
!===============================================================================
subroutine chem_register 1,4
!-----------------------------------------------------------------------
!
! Purpose: register advected constituents for parameterized greenhouse gas chemistry
!
! Method:
! <Describe the algorithm(s) used in the routine.>
! <Also include any applicable external references.>
!
! Author: B.A. Boville
!
!-----------------------------------------------------------------------
!---------------------------Local workspace-----------------------------
integer :: m ! tracer index
!-----------------------------------------------------------------------
! Set names of diffused variable tendencies and declare them as history variables
call cnst_add(cnst_names(1), advected, mwn2o, cpn2o, 0._r8, ixghg, longname='Nitrous Oxide')
call cnst_add(cnst_names(2), advected, mwch4, cpch4, 0._r8, m, longname='Methane')
call cnst_add(cnst_names(3), advected, mwf11, cpf11, 0._r8, m)
call cnst_add(cnst_names(4), advected, mwf12, cpf12, 0._r8, m)
return
end subroutine chem_register
!===============================================================================
function chem_implements_cnst(name)
!-----------------------------------------------------------------------
!
! Purpose: return true if specified constituent is implemented by this package
!
! Author: B. Eaton
!
!-----------------------------------------------------------------------
implicit none
!-----------------------------Arguments---------------------------------
character(len=*), intent(in) :: name ! constituent name
logical :: chem_implements_cnst ! return value
!---------------------------Local workspace-----------------------------
integer :: m
!-----------------------------------------------------------------------
chem_implements_cnst = .false.
do m = 1, ncnst
if (name == cnst_names(m)) then
chem_implements_cnst = .true.
return
end if
end do
end function chem_implements_cnst
!===============================================================================
subroutine chem_init 1,4
!-----------------------------------------------------------------------
!
! Purpose: initialize parameterized greenhouse gas chemistry
! (declare history variables)
!
! Method:
! <Describe the algorithm(s) used in the routine.>
! <Also include any applicable external references.>
!
! Author: NCAR CMS
!
!-----------------------------------------------------------------------
use history, only: addfld, add_default, phys_decomp
!---------------------------Local workspace-----------------------------
integer :: m ! tracer index
!-----------------------------------------------------------------------
! Set names of diffused variable tendencies and declare them as history variables
do m = 1, 4
srcnam(m) = trim(cnst_name(ixghg-1+m)) // 'SRC'
call addfld (srcnam(m),'kg/kg/s ',pver, 'A',trim(cnst_name(ixghg-1+m))//' source/sink',phys_decomp)
call add_default (srcnam(m), 1, ' ')
end do
call chem_init_loss
return
end subroutine chem_init
!===============================================================================
subroutine chem_timestep_tend (state, ptend, cflx, dt, fh2o) 1,7
!-----------------------------------------------------------------------
!
! Purpose:
! Interface to parameterized greenhouse gas chemisty (source/sink).
!
! Method:
! <Describe the algorithm(s) used in the routine.>
! <Also include any applicable external references.>
!
! Author: B.A. Boville
!
!-----------------------------------------------------------------------
use history, only: outfld
use physics_types, only: physics_state, physics_ptend, physics_ptend_init
use phys_grid, only: get_lat_all_p
!-----------------------------------------------------------------------
implicit none
!-----------------------------------------------------------------------
!
! Arguments:
!
real(r8), intent(in) :: dt ! time step
type(physics_state), intent(in ) :: state ! Physics state variables
type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies
real(r8), intent(inout) :: cflx(pcols,ppcnst) ! Surface constituent flux (kg/m^2/s)
real(r8), intent(out) :: fh2o(pcols) ! H2O flux required to balance the H2O source from methane chemistry (kg/m^2/s)
!
! Local variables
!
integer :: i, k, m ! indices
integer :: ioff ! offset for ghg indices
integer :: lchnk ! chunk identifier
integer :: ncol ! number of atmospheric columns
integer :: lat(pcols) ! latitude index for S->N storage
!
!-----------------------------------------------------------------------
! Initialize output tendency structure
call physics_ptend_init(ptend)
ioff = ixghg - 1
lchnk = state%lchnk
ncol = state%ncol
! get latitude indices
call get_lat_all_p(lchnk, ncol, lat)
! compute tendencies and surface fluxes
call ghg_chem ( lchnk, ncol, lat, &
state%q(:,:,1), state%q(:,:,ixghg:ixghg+3), &
ptend%q(:,:,1), ptend%q(:,:,ixghg:ixghg+3), cflx(:,ixghg:ixghg+3), dt)
! set flags for tracer tendencies (water and 4 ghg's)
ptend%lq(1) = .TRUE.
ptend%lq(ioff+1:ioff+4) = .TRUE.
!
! record tendencies on history files
do m = 1, 4
call outfld (srcnam(m),ptend%q(:,:,ioff+m),pcols,lchnk)
end do
! Compute water vapor flux required to make the conservation checker happy
fh2o = 0
do k = 1, pver
do i = 1, ncol
fh2o(i) = fh2o(i) + ptend%q(i,k,1)*state%pdel(i,k)/gravit
end do
end do
return
end subroutine chem_timestep_tend
!===============================================================================
subroutine ghg_chem (lchnk, ncol, lat, qh2o, qghg, dqh2o, dqghg, fghg, dt) 1
!-----------------------------------------------------------------------
!
! Purpose:
! Apply the interpolated chemical loss rates from the input data to
! N2O, CH4, CFC11 and CFC12. Set the surface values to a constant.
!
! Method:
! <Describe the algorithm(s) used in the routine.>
! <Also include any applicable external references.>
!
! Author: NCAR CMS
!
!-----------------------------------------------------------------------
implicit none
!-----------------------------------------------------------------------
!
! Arguments:
!
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncol ! number of atmospheric columns
integer, intent(in) :: lat(pcols) ! latitude index for S->N storage
real(r8), intent(in) :: dt ! time step
real(r8), intent(in) :: qh2o(pcols,pver) ! mass mixing ratios of water vapor
real(r8), intent(in) :: qghg(pcols,pver,4) ! mass mixing ratios of greenhouse gases
real(r8), intent(out) :: dqh2o(pcols,pver) ! tendency of mass mixing ratios (water)
real(r8), intent(out) :: dqghg(pcols,pver,4) ! tendency of mass mixing ratios (ghg's)
real(r8), intent(out) :: fghg(pcols,4) ! Surface constituent flux (kg/m^2/s)
!
! Local variables
!
integer i,k ! loop indexes
real(r8) xch4 ! new methane mass mixing ratio
real(r8) xn2o ! new nitrous oxide mass mixing ratio
real(r8) xcfc11 ! new cfc11 mass mixing ratio
real(r8) xcfc12 ! new cfc12 mass mixing ratio
!
!-----------------------------------------------------------------------
!
! Apply chemical rate coefficient using time split implicit method. The
! turn the new value back into a tendency. NOTE that water tendency is
! twice methane tendency. Water is specific humidity which is in mass
! mixing ratio units. Note that
! o 1 => indx of n2o
! o 2 => indx of ch4
! o 3 => indx of cfc11
! o 4 => indx of cfc12
!
do k=1,pver-2
do i=1,ncol
xn2o = qghg(i,k,1) / (1. + tn2o (lat(i),k) * dt)
xch4 = qghg(i,k,2) / (1. + tch4 (lat(i),k) * dt)
xcfc11 = qghg(i,k,3) / (1. + tcfc11(lat(i),k) * dt)
xcfc12 = qghg(i,k,4) / (1. + tcfc12(lat(i),k) * dt)
dqghg(i,k,1) =(xn2o - qghg(i,k,1)) / dt
dqghg(i,k,2) =(xch4 - qghg(i,k,2)) / dt
dqghg(i,k,3) =(xcfc11 - qghg(i,k,3)) / dt
dqghg(i,k,4) =(xcfc12 - qghg(i,k,4)) / dt
dqh2o(i,k) = -2. * rh2och4 * dqghg(i,k,2)
end do
end do
!
! Set the "surface" tendencies (bottom 2 levels) to maintain specified
! tropospheric concentrations.
!
do k = pver-1, pver
do i=1,ncol
dqghg(i,k,1) =((rmwn2o*n2ovmr) - qghg(i,k,1)) / dt
dqghg(i,k,2) =((rmwch4*ch4vmr) - qghg(i,k,2)) / dt
dqghg(i,k,3) =((rmwf11*f11vmr) - qghg(i,k,3)) / dt
dqghg(i,k,4) =((rmwf12*f12vmr) - qghg(i,k,4)) / dt
dqh2o(i,k) = 0.
end do
end do
!
! For now set all tracer fluxes to 0
!
do i=1,ncol
fghg(i,1) = 0.
fghg(i,2) = 0.
fghg(i,3) = 0.
fghg(i,4) = 0.
end do
return
end subroutine ghg_chem
!===============================================================================
subroutine chem_init_loss 1,53
!-----------------------------------------------------------------------
!
! Purpose:
! Do initial read of time-variant chemical loss rate frequency dataset, containing
! loss rates as a function of latitude and pressure. Determine the two
! consecutive months between which the current date lies.
!
! Method:
! <Describe the algorithm(s) used in the routine.>
! <Also include any applicable external references.>
!
! Author: NCAR CMS
!-----------------------------------------------------------------------
use ioFileMod
use commap
use time_manager, only: get_curr_date, get_perp_date, get_curr_calday, &
is_perpetual
use filenames, only: bndtvg
#if ( defined SPMD )
use mpishorthand
#endif
implicit none
#include <comctl.h>
#include <comhyb.h>
#include <comlun.h>
include 'netcdf.inc'
!
! Local variables
!
integer dateid ! netcdf id for date variable
integer secid ! netcdf id for seconds variable
integer lonid ! netcdf id for longitude variable
integer latid ! netcdf id for latitude variable
integer levid ! netcdf id for level variable
integer timid ! netcdf id for time variable
integer tch4id ! netcdf id for ch4 loss rate
integer tn2oid ! netcdf id for n2o loss rate
integer tcfc11id ! netcdf id for cfc11 loss rate
integer tcfc12id ! netcdf id for cfc12 loss rate
integer lonsiz ! size of longitude dimension on tracer dataset
integer levsiz ! size of level dimension on tracer dataset
integer latsiz ! size of latitude dimension on tracer dataset
integer timsiz ! size of time dimension on tracer dataset
integer j,n,k,nt ! indices
integer ki,ko,ji,jo ! indices
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) lato(plat) ! cam model latitudes (degrees)
real(r8) zo(plev) ! cam model heights (m)
real(r8) lati(ptrlat) ! input data latitudes (degrees)
real(r8) pin(ptrlev) ! input data pressure values (mbars)
real(r8) zi(ptrlev) ! input data heights (m)
real(r8) xch4i (ptrlon,ptrlev,ptrlat,ptrtim) ! input ch4 loss rate coeff
real(r8) xn2oi (ptrlon,ptrlev,ptrlat,ptrtim) ! input n2o loss rate coeff
real(r8) xcfc11i(ptrlon,ptrlev,ptrlat,ptrtim) ! input cfc11 loss rate coeff
real(r8) xcfc12i(ptrlon,ptrlev,ptrlat,ptrtim) ! input cfc12 loss rate coeff
real(r8) xch4(ptrlat, ptrlev) ! input ch4 loss rate coeff indices changed
real(r8) xn2o(ptrlat, ptrlev) ! input n2o loss rate coeff indices changed
real(r8) xcfc11(ptrlat, ptrlev) ! input cfc11 loss rate coeff indices changed
real(r8) xcfc12(ptrlat, ptrlev) ! input cfc12 loss rate coeff indices changed
real(r8) xch4lv(ptrlat, plev) ! input ch4 loss rate coeff interp to cam levels
real(r8) xn2olv(ptrlat, plev) ! input n2o loss rate coeff interp to cam levels
real(r8) xcfc11lv(ptrlat, plev) ! input cfc11 loss rate coeff interp to cam levels
real(r8) xcfc12lv(ptrlat, plev) ! input cfc12 loss rate coeff interp to cam levels
character(len=256) :: locfn ! netcdf local filename to open
!
!-----------------------------------------------------------------------
!
! Initialize
!
nm = 1
np = 2
!
! SPMD: Master reads dataset and does spatial interpolation on all time samples.
! Spatial interpolents are broadcast. All subsequent time interpolation is
! done in every process.
!
if (masterproc) then
write(6,*)'CHEM_INIT_LOSS: greenhouse loss rates from: ', trim(bndtvg)
call getfil(bndtvg, locfn)
call wrap_open(locfn, 0, ncid_trc)
write(6,*)'CHEM_INIT_LOSS: NCOPN returns id ',ncid_trc,' for file ',trim(locfn)
!
!------------------------------------------------------------------------
! Read tracer data
!------------------------------------------------------------------------
!
! Get dimension info
!
call wrap_inq_dimid(ncid_trc, 'lat' , latid)
call wrap_inq_dimid(ncid_trc, 'lev' , levid)
call wrap_inq_dimid(ncid_trc, 'lon' , lonid)
call wrap_inq_dimid(ncid_trc, 'time', timid)
call wrap_inq_dimlen(ncid_trc, lonid, lonsiz)
call wrap_inq_dimlen(ncid_trc, levid, levsiz)
call wrap_inq_dimlen(ncid_trc, latid, latsiz)
call wrap_inq_dimlen(ncid_trc, timid, timsiz)
!
! Check dimension info
!
if (ptrlon/=1) then
call endrun ('CHEM_INIT_LOSS: longitude dependence not implemented')
endif
if (lonsiz /= ptrlon) then
write(6,*)'CHEM_INIT_LOSS: lonsiz=',lonsiz,' must = ptrlon=',ptrlon
call endrun
end if
if (levsiz /= ptrlev) then
write(6,*)'CHEM_INIT_LOSS: levsiz=',levsiz,' must = ptrlev=',ptrlev
call endrun
end if
if (latsiz /= ptrlat) then
write(6,*)'CHEM_INIT_LOSS: latsiz=',latsiz,' must = ptrlat=',ptrlat
call endrun
end if
if (timsiz /= ptrtim) then
write(6,*)'CHEM_INIT_LOSS: timsiz=',timsiz,' must = ptrtim=',ptrtim
call endrun
end if
!
! Determine necessary dimension and variable id's
!
call wrap_inq_varid(ncid_trc, 'lat' , latid)
call wrap_inq_varid(ncid_trc, 'lev' , levid)
call wrap_inq_varid(ncid_trc, 'date' , dateid)
call wrap_inq_varid(ncid_trc, 'datesec', secid)
call wrap_inq_varid(ncid_trc, 'TCH4' , tch4id)
call wrap_inq_varid(ncid_trc, 'TN2O' , tn2oid)
call wrap_inq_varid(ncid_trc, 'TCFC11' , tcfc11id)
call wrap_inq_varid(ncid_trc, 'TCFC12' , tcfc12id)
!
! Obtain entire date and sec variables. Assume that will always
! cycle over 12 month data.
!
call wrap_get_var_int(ncid_trc, dateid, date_tr)
call wrap_get_var_int(ncid_trc, secid , sec_tr)
if (mod(date_tr(1),10000)/100 /= 1) then
call endrun ('(CHEM_INIT_LOSS): error when cycling data: 1st month must be 1')
end if
if (mod(date_tr(ptrtim),10000)/100 /= 12) then
call endrun ('(CHEM_INIT_LOSS): error when cycling data: last month must be 12')
end if
!
! Obtain input data latitude and level arrays.
!
call wrap_get_var_realx(ncid_trc, latid, lati)
call wrap_get_var_realx(ncid_trc, levid, pin )
!
! Convert input pressure levels to height (m).
! First convert from millibars to pascals.
!
do k=1,ptrlev
pin(k) = pin(k)*100.
zi(k) = 7.0e3 * log (1.0e5 / pin(k))
end do
!
! Convert approximate cam pressure levels to height (m).
!
do k=1,plev
zo (k) = 7.0e3 * log (1.0e5 / hypm(k))
end do
!
! Convert cam model latitudes to degrees.
! Input model latitudes already in degrees.
!
do j=1,plat
lato(j) = clat(j)*45./atan(1.)
end do
!
! Obtain all time samples of tracer data.
!
call wrap_get_var_realx(ncid_trc, tch4id , xch4i )
call wrap_get_var_realx(ncid_trc, tn2oid , xn2oi )
call wrap_get_var_realx(ncid_trc, tcfc11id, xcfc11i)
call wrap_get_var_realx(ncid_trc, tcfc12id, xcfc12i)
!
! Close netcdf file
!
call wrap_close(ncid_trc)
!
!------------------------------------------------------------------------
! Interpolate tracer data to model grid
!------------------------------------------------------------------------
!
! Loop over all input times.
!
do nt = 1, ptrtim
!
! Remove longitude and time index and switch level and latitude indices
! for the loss coefficients.
!
do j=1,ptrlat
do k=1,ptrlev
xch4 (j,k) = xch4i (1,k,j,nt)
xn2o (j,k) = xn2oi (1,k,j,nt)
xcfc11(j,k) = xcfc11i(1,k,j,nt)
xcfc12(j,k) = xcfc12i(1,k,j,nt)
end do
end do
!
! Interpolate input data to model levels.
! If the CAM level is outside the range of the input data (this
! can happen only in troposphere) put zero for every latitude.
! Otherwise determine the input data levels bounding the current
! CAM level and interpolate.
!
do ko=1,plev
if (zo(ko) < zi(ptrlev)) then
do j=1,ptrlat
xch4lv (j,ko) = 0.0
xn2olv (j,ko) = 0.0
xcfc11lv(j,ko) = 0.0
xcfc12lv(j,ko) = 0.0
end do
goto 50
end if
do ki=1,ptrlev-1
if (zo(ko) < zi(ki) .and. zo(ko) >= zi(ki+1)) then
do j=1,ptrlat
xch4lv(j,ko) = xch4(j,ki) + (xch4(j,ki+1) - xch4(j,ki)) &
/ (zi(ki+1) - zi(ki)) * (zo(ko) - zi(ki))
xn2olv(j,ko) = xn2o(j,ki) + (xn2o(j,ki+1) - xn2o(j,ki)) &
/ (zi(ki+1) - zi(ki)) * (zo(ko) - zi(ki))
xcfc11lv(j,ko) = xcfc11(j,ki) + (xcfc11(j,ki+1) - xcfc11(j,ki)) &
/ (zi(ki+1) - zi(ki)) * (zo(ko) - zi(ki))
xcfc12lv(j,ko) = xcfc12(j,ki) + (xcfc12(j,ki+1) - xcfc12(j,ki)) &
/ (zi(ki+1) - zi(ki)) * (zo(ko) - zi(ki))
end do
goto 50
endif
end do
call endrun ('(CHEM_INIT_LOSS): Error in vertical interpolation')
50 continue
end do
!
! Interpolate input data to model latitudes.
! Determine the input data latitudes bounding the current CAM latitude and
! interpolate. Use last value from input data if the cam latitude is
! outside the range of the input data latitudes.
!
do jo=1,plat
if (lato(jo) <= lati(1)) then
do k = 1, plev
tch4i(jo,k,nt) = xch4lv(1,k)
tn2oi(jo,k,nt) = xn2olv(1,k)
tcfc11i(jo,k,nt) = xcfc11lv(1,k)
tcfc12i(jo,k,nt) = xcfc12lv(1,k)
end do
else if (lato(jo) >= lati(ptrlat)) then
do k = 1, plev
tch4i(jo,k,nt) = xch4lv(ptrlat,k)
tn2oi(jo,k,nt) = xn2olv(ptrlat,k)
tcfc11i(jo,k,nt) = xcfc11lv(ptrlat,k)
tcfc12i(jo,k,nt) = xcfc12lv(ptrlat,k)
end do
else
do ji=1,ptrlat-1
if ( (lato(jo) > lati(ji)) .and. (lato(jo) <= lati(ji+1))) then
do k=1,plev
tch4i(jo,k,nt) = xch4lv(ji,k) + (xch4lv(ji+1,k) - xch4lv(ji,k)) &
/ (lati(ji+1) - lati(ji)) * (lato(jo) - lati(ji))
tn2oi(jo,k,nt) = xn2olv(ji,k) + (xn2olv(ji+1,k) - xn2olv(ji,k)) &
/ (lati(ji+1) - lati(ji)) * (lato(jo) - lati(ji))
tcfc11i(jo,k,nt) = xcfc11lv(ji,k) + (xcfc11lv(ji+1,k) - xcfc11lv(ji,k)) &
/ (lati(ji+1) - lati(ji)) * (lato(jo) - lati(ji))
tcfc12i(jo,k,nt) = xcfc12lv(ji,k) + (xcfc12lv(ji+1,k) - xcfc12lv(ji,k)) &
/ (lati(ji+1) - lati(ji)) * (lato(jo) - lati(ji))
end do
goto 90
endif
end do
end if
write (6,*)'(CHEM_INIT_LOSS): Error in horizontal interpolation'
90 continue
end do
end do ! end loop over time samples
endif ! end of masterproc
#if ( defined SPMD )
call mpibcast (date_tr, ptrtim, mpiint, 0, mpicom)
call mpibcast (sec_tr , ptrtim, mpiint, 0, mpicom)
call mpibcast (tch4i , plev*plat*ptrtim, mpir8, 0, mpicom)
call mpibcast (tn2oi , plev*plat*ptrtim, mpir8, 0, mpicom)
call mpibcast (tcfc11i, plev*plat*ptrtim, mpir8, 0, mpicom)
call mpibcast (tcfc12i, plev*plat*ptrtim, mpir8, 0, mpicom)
#endif
!
! Initial time interpolation between December and January
!
calday = get_curr_calday()
if ( is_perpetual() ) then
call get_perp_date(yr, mon, day, ncsec)
else
call get_curr_date(yr, mon, day, ncsec)
end if
ncdate = yr*10000 + mon*100 + day
n = 12
np1 = 1
call bnddyi(date_tr(n ), sec_tr(n ), cdaytrm)
call bnddyi(date_tr(np1), sec_tr(np1), cdaytrp)
if (calday <= cdaytrp .or. calday > cdaytrm) then
do j=1,plat
do k=1,plev
tch4m (j,k,nm) = tch4i (j,k,n)
tn2om (j,k,nm) = tn2oi (j,k,n)
tcfc11m(j,k,nm) = tcfc11i(j,k,n)
tcfc12m(j,k,nm) = tcfc12i(j,k,n)
tch4m (j,k,np) = tch4i (j,k,np1)
tn2om (j,k,np) = tn2oi (j,k,np1)
tcfc11m(j,k,np) = tcfc11i(j,k,np1)
tcfc12m(j,k,np) = tcfc12i(j,k,np1)
end do
end do
goto 10
end if
!
! Initial normal interpolation between consecutive time slices.
!
do n=1,ptrtim-1
np1 = n + 1
call bnddyi(date_tr(n ), sec_tr(n ), cdaytrm)
call bnddyi(date_tr(np1), sec_tr(np1), cdaytrp)
if (calday > cdaytrm .and. calday <= cdaytrp) then
do j=1,plat
do k=1,plev
tch4m (j,k,nm) = tch4i (j,k,n)
tn2om (j,k,nm) = tn2oi (j,k,n)
tcfc11m(j,k,nm) = tcfc11i(j,k,n)
tcfc12m(j,k,nm) = tcfc12i(j,k,n)
tch4m (j,k,np) = tch4i (j,k,np1)
tn2om (j,k,np) = tn2oi (j,k,np1)
tcfc11m(j,k,np) = tcfc11i(j,k,np1)
tcfc12m(j,k,np) = tcfc12i(j,k,np1)
end do
end do
goto 10
end if
end do
write(6,*)'CHEM_INIT_LOSS: Failed to find dates bracketing ncdate, ncsec=', ncdate, ncsec
call endrun
!
! Data positioned correctly
!
10 continue
return
end subroutine chem_init_loss
!===============================================================================
subroutine chem_timestep_init 1,9
!-----------------------------------------------------------------------
!
! Purpose:
! Time interpolate chemical loss rates to current time, updating the
! bounding month arrays with new data if necessary
!
! Method:
! <Describe the algorithm(s) used in the routine.>
! <Also include any applicable external references.>
!
! Author: NCAR CMS
!
!-----------------------------------------------------------------------
use commap
use time_manager, only: get_curr_date, get_perp_date, get_curr_calday, &
is_perpetual
implicit none
#include <comctl.h>
!-----------------------------------------------------------------------
!
! Local workspace
!
integer ntmp ! temporary
integer j,k ! indices
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) fact1, fact2 ! time interpolation factors
!-----------------------------------------------------------------------
!
! If model time is past current forward timeslice, obtain the next
! timeslice for time interpolation. Messy logic is for
! interpolation between December and January (np1.eq.1).
!
calday = get_curr_calday()
if ( is_perpetual() ) then
call get_perp_date(yr, mon, day, ncsec)
else
call get_curr_date(yr, mon, day, ncsec)
end if
ncdate = yr*10000 + mon*100 + day
if (calday > cdaytrp .and. .not. (np1 == 1 .and. calday > cdaytrm)) then
np1 = mod(np1,12) + 1
if (np1 > ptrtim) then
call endrun ('CHEMINT: Attempt to access bad month')
end if
cdaytrm = cdaytrp
call bnddyi(date_tr(np1), sec_tr(np1), cdaytrp)
if (np1 == 1 .or. calday <= cdaytrp) then
ntmp = nm
nm = np
np = ntmp
do j=1,plat
do k=1,plev
tch4m (j,k,np) = tch4i (j,k,np1)
tn2om (j,k,np) = tn2oi (j,k,np1)
tcfc11m(j,k,np) = tcfc11i(j,k,np1)
tcfc12m(j,k,np) = tcfc12i(j,k,np1)
end do
end do
else
write(6,*)'CHEMINT: Input data for date',date_tr(np1), &
' sec ',sec_tr(np1), 'does not exceed model date', &
ncdate,' sec ',ncsec,' Stopping.'
call endrun
end if
end if
!
! Determine time interpolation factors. 1st arg says we are cycling yearly data
!
call getfactors (.true., np1, cdaytrm, cdaytrp, calday, &
fact1, fact2, 'SULFINT:')
!
! Do time interpolation
!
do j=1,plat
do k=1,plev
tch4(j,k) = tch4m (j,k,nm)*fact1 + tch4m (j,k,np)*fact2
tn2o(j,k) = tn2om (j,k,nm)*fact1 + tn2om (j,k,np)*fact2
tcfc11(j,k) = tcfc11m(j,k,nm)*fact1 + tcfc11m(j,k,np)*fact2
tcfc12(j,k) = tcfc12m(j,k,nm)*fact1 + tcfc12m(j,k,np)*fact2
end do
end do
return
end subroutine chem_timestep_init
!===============================================================================
subroutine chem_init_cnst(name, q) 1
!-----------------------------------------------------------------------
!
! Purpose:
! Set initial mass mixing ratios of CH4, N2O, CFC11 and CFC12.
!
!-----------------------------------------------------------------------
implicit none
!-----------------------------Arguments---------------------------------
character(len=*), intent(in) :: name ! constituent name
real(r8), intent(out) :: q(plon,plev,plat) ! mass mixing ratio
!-----------------------------------------------------------------------
if ( name == 'N2O' ) then
q = rmwn2o * n2ovmr
return
else if ( name == 'CH4' ) then
q = rmwch4 * ch4vmr
return
else if ( name == 'CFC11' ) then
q = rmwf11 * f11vmr
return
else if ( name == 'CFC12' ) then
q = rmwf12 * f12vmr
return
end if
end subroutine chem_init_cnst
end module chemistry