#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 endrun50           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)