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_mixend module chemistry