chemistry.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

                   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

    endif                     ! end of masterproc

    return
  end subroutine chem_init_loss

!===============================================================================
  subroutine chem_time_interp
!----------------------------------------------------------------------- 
! 
! Purpose: 
! Time interpolate chemical loss rates to current time, reading
! in new monthly 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
#if ( defined SPMD )
    use mpishorthand
#endif

    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
    real(r8) deltat            ! time (days) between interpolating ozone data

!-----------------------------------------------------------------------
!
! SPMD: Master does all the work.  Sends needed info to slaves
!
    if ( masterproc) then
!
! 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
             write(6,*)'CHEMINT: Attempt to access bad month'
             call endrun
          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 factors for time interpolation.
!
       if (np1 == 1) then     ! Dec-Jan interpolation
          deltat = cdaytrp + 365. - cdaytrm
          if (calday > cdaytrp) then ! We're in December
             fact1 = (cdaytrp + 365. - calday)/deltat
             fact2 = (calday - cdaytrm)/deltat
          else                ! We're in January
             fact1 = (cdaytrp - calday)/deltat
             fact2 = (calday + 365. - cdaytrm)/deltat
          end if
       else                   ! Non Dec-Jan interpolation
          deltat = cdaytrp - cdaytrm
          fact1 = (cdaytrp - calday)/deltat
          fact2 = (calday - cdaytrm)/deltat
       end if
!
! Check sanity of time interpolation factors to within 32-bit roundoff.
!
       if (abs(fact1+fact2-1.) > 1.e-6 .or. &
            fact1 > 1.000001 .or. fact1 < -1.e-6 .or. &
            fact2 > 1.000001 .or. fact2 < -1.e-6) then
          write(6,*)'CHEMINT: Bad fact1 and/or fact2=',fact1,fact2
          call endrun
       end if
!
! 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
    end if			! end of if-masterproc
#if ( defined SPMD )
    call mpibcast (tch4  , plev*plat,mpir8,0,mpicom)
    call mpibcast (tn2o  , plev*plat,mpir8,0,mpicom)
    call mpibcast (tcfc11, plev*plat,mpir8,0,mpicom)
    call mpibcast (tcfc12, plev*plat,mpir8,0,mpicom)
#endif

    return
  end subroutine chem_time_interp

!===============================================================================
  subroutine chem_init_mix(lat, ps, q, nlon)
!----------------------------------------------------------------------- 
! 
! Purpose: 
! Specify initial mass mixing ratios of CH4, N2O, CFC11 and CFC12.
! Distributions assume constant mixing ratio in the troposphere
! and a decrease of mixing ratio in the stratosphere. Tropopause
! defined by ptrop. The scale height of the particular trace gas
! depends on latitude. This assumption produces a more realistic
! stratospheric distribution of the various trace gases.
! 
! Method: 
! <Describe the algorithm(s) used in the routine.> 
! <Also include any applicable external references.> 
! 
! Author: J.T. Kiehl
! 
!-----------------------------------------------------------------------

    use commap

    implicit none

#include <comhyb.h>
!-----------------------------Arguments---------------------------------
!
! Input
!
    integer, intent(in) :: lat                  ! current latitude index
    integer, intent(in) :: nlon
    real(r8), intent(in) :: ps(plond)               ! surface pressure
!
! Input/Output
!
    real(r8), intent(inout) :: q(plond,plev,ppcnst)           !  mass mixing ratio

!
!--------------------------Local Variables------------------------------
!
    integer i                ! longitude loop index
    integer k                ! level index
!
    real(r8) coslat              ! cosine of latitude
    real(r8) dlat                ! latitude in degrees
    real(r8) pmid(plond,plev)    ! model pressures
    real(r8) ptrop               ! pressure level of tropopause
    real(r8) pratio              ! pressure divided by ptrop
!
    real(r8) xn2o                ! pressure scale height for n2o
    real(r8) xch4                ! pressure scale height for ch4
    real(r8) xcfc11              ! pressure scale height for cfc11
    real(r8) xcfc12              ! pressure scale height for cfc12
!
    real(r8) ch40                ! tropospheric mass mixing ratio for ch4
    real(r8) n2o0                ! tropospheric mass mixing ratio for n2o
    real(r8) cfc110              ! tropospheric mass mixing ratio for cfc11
    real(r8) cfc120              ! tropospheric mass mixing ratio for cfc12
!
!-----------------------------------------------------------------------
!   For 1990, I get the following:
!   
!   CFC-11 mixing ratio = 270.1 pptv
!   CFC-12 mixing ratio = 465.3 pptv
!   
!   CFC-11 radiative forcing = 0.0594 W/m2
!   CFC-12 radiative forcing = 0.1300 W/m2
!   Other CFC radiative forcing = 0.0595 W/m2
!   
!   Hence, the magic number is "TWO".
!   
!   These values compare reasonably well to IPCC.   They have slightly higher
!   CFC-11 and 12 mixing ratios in the l990 report (280 and 484), but I think
!   mine are better.   They had 0.062 and 0.14 for the corresponding forcings.
!   They also had higher CCl4 and CFC-113 mixing ratios (in one case by quite a
!   bit) and had radiative forcing by other CFCs at 0.085.  I think the above
!   number can be defended as better.
!   
!   Just in case you want it, my estimates for l990 CH4 and N2O would be:
!   
!   CH4 - 1.715 ppmv
!   N2O - 310.0 ppbv
!   
!   but you probably want to use Tom's numbers instead for those.
!   
!   Thanks,
!   Susan
!-----------------------------------------------------------------------
!
! 1990 ch4 vmr from Wigley   : 1.722e-6
! 1990 n2o vmr from Wigley   : 308.4e-9
! 1990 cfc11 vmr from Solomon: 270.1e-12 * 2.
! factor of 2 on cfc11 to account for all other halocarbons
! 1900 cfc12 vmr from Solomon: 465.3e-12
!
    ch40   = rmwch4 * ch4vmr
    n2o0   = rmwn2o * n2ovmr
    cfc110 = rmwf11 * f11vmr
    cfc120 = rmwf12 * f12vmr
!
! Set stratospheric scale height factor for gases
!
    dlat = abs(latdeg(lat))
    if(dlat.le.45.0) then
       xn2o = 0.3478 + 0.00116 * dlat
       xch4 = 0.2353
       xcfc11 = 0.7273 + 0.00606 * dlat
       xcfc12 = 0.4000 + 0.00222 * dlat
    else
       xn2o = 0.4000 + 0.013333  * (dlat - 45)
       xch4 = 0.2353 + 0.0225489 * (dlat - 45)
       xcfc11 = 1.00 + 0.013333  * (dlat - 45)
       xcfc12 = 0.50 + 0.024444  * (dlat - 45)
    end if
!
! Set pressure of tropopause and model layer pressures
!
    coslat = cos(clat(lat))
    ptrop = 250.0e2 - 150.0e2*coslat**2.0

    do k=1,plev
       do i=1,nlon
          pmid(i,k) = hyam(k)*ps0 + hybm(k)*ps(i)
       end do
    end do
!
! Determine initial mixing ratios
!
    do k = 1,plev
       do i = 1,nlon
          if (pmid(i,k) >= ptrop) then
             q(i,k,ixghg  ) = n2o0
             q(i,k,ixghg+1) = ch40
             q(i,k,ixghg+2) = cfc110
             q(i,k,ixghg+3) = cfc120
          else
             pratio = pmid(i,k)/ptrop
             q(i,k,ixghg  ) = n2o0   * (pratio)**xn2o
             q(i,k,ixghg+1) = ch40   * (pratio)**xch4
             q(i,k,ixghg+2) = cfc110 * (pratio)**xcfc11
             q(i,k,ixghg+3) = cfc120 * (pratio)**xcfc12
          end if
       end do
    end do
!
! Adjust water using methane mass mxiing ratio
!
    do k=1,plev
       do i=1,nlon
          q(i,k,1) = q(i,k,1) + 2.*(ch40-q(i,k,ixghg+1))
       end do
    end do
!
    return
  end subroutine chem_init_mix

end module chemistry