chemistry.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#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 does all the work.  Sends needed info to slaves
!
    if (masterproc) then
       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
          write(6,*)'CHEM_INIT_LOSS: longitude dependence not implemented'
          call endrun
       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
          write(6,*)'(CHEM_INIT_LOSS): error when cycling data: 1st month must be 1'
          call endrun
       end if
       if (mod(date_tr(ptrtim),10000)/100 /= 12) then
          write(6,*)'(CHEM_INIT_LOSS): error when cycling data: last month must be 12'
          call endrun
       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
             write (6,*) '(CHEM_INIT_LOSS): Error in vertical interpolation'
             call endrun
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
!
! 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,timsiz-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)