#include <misc.h>
#include <preproc.h>
      subroutine lsm_ccmintr (nstep  ,londim   ,pcnst    ,pnats    ,
     &                        pgcmxy ,thgcmxy  ,ugcmxy   ,
     &                        vgcmxy ,tgcmxy   ,qgcmxy   ,flwdsxy  ,
     &                        zgcmxy ,prcxy    ,prlxy    ,snowcxy  ,
     &                        snowlxy,solsxy   ,sollxy   ,solsdxy  ,
     &                        solldxy,tsxy     ,shxy     ,lhxy     ,
     &                        cfxy   ,tauxxy   ,tauyxy   ,trefxy   ,
     &                        asdirxy,asdifxy  ,aldirxy  ,aldifxy  ,
     &                        lwupxy ,snowxy   ,eccen    ,obliqr   ,
     &                        lambm0 ,mvelpp  )

#if (! defined COUP_CSM) 
! ------------------------ code history ---------------------------
! source file:       lsm_ccmintr
! purpose:           ccm 2-d lsm interface 
! author:            Mariana Vertenstein
! $Id: lsm_ccmintr.F,v 1.4 2000/04/13 21:17:13 erik Exp $
! -----------------------------------------------------------------

! ------------------------ notes ----------------------------------
! the land surface model works by gathering all the land points on a
! [lsmlon] by [lsmlat] grid into a vector of [lpt] land points. this is
! then expanded into a "big" vector of [kpt] subgrid points, allowing for
! up to [msub] subgrid points per land point. the [kpt] subgrid points
! are processed as [numlv] "little" vectors of [numkpt] points for [numlv] 
! calls to subroutine lsm. this subroutine:
! -----------------------------------------------------------------

      use precision
      use lsmpar
      use lsmvar
      implicit none
#include <phycon.h>
#include <lsmctl.h>

! ------------------- arguments -----------------------------------
! input variables

      integer londim    !atm number of longitudes
      integer pcnst     !atm number of adv constituents (water is first)
      integer pnats     !atm number of non-adv constituents
      integer nstep     !atm time step index

      real(r8) eccen    !Earth's orbital eccentricity
      real(r8) obliqr   !Earth's obliquity in radians
      real(r8) lambm0   !Mean longitude of perihelion at the vernal
                        !equinox (radians)
      real(r8) mvelpp   !Earth's moving vernal equinox longitude of 
                        !perihelion plus pi (radians)

      real(r8) pgcmxy(londim,beglat:endlat) ! atm bottom level pressure (pa)
      real(r8) thgcmxy(londim,beglat:endlat)! atm btm lvl pot temp (kelvin)
      real(r8) ugcmxy(londim,beglat:endlat) ! atm btm level zonal wind (m/s)
      real(r8) vgcmxy(londim,beglat:endlat) ! atm btm lev merid wind (m/s)
      real(r8) tgcmxy(londim,beglat:endlat) ! atm bottom level temp (kelvin)
      real(r8) qgcmxy(londim,beglat:endlat) ! atm btm lvl spc humidity(kg/kg)
      real(r8) zgcmxy(londim,beglat:endlat) ! atm btm lev hght above srf (m)
      real(r8) prcxy(londim,beglat:endlat)  ! conv precipitation rate (m h2o/s) 
      real(r8) prlxy(londim,beglat:endlat)  ! large-scale precipitation rate(m h2o/s)
      real(r8) snowcxy(londim,beglat:endlat)! conv snowfall rate (m h2o/s) 
      real(r8) snowlxy(londim,beglat:endlat)! large-scale snowfall rate(m h2o/s)
      real(r8) flwdsxy(londim,beglat:endlat)! downward longwave rad onto surface (w/m**2)
      real(r8) solsxy (londim,beglat:endlat)! vis direct beam solar rad onto srf (w/m**2)
      real(r8) sollxy (londim,beglat:endlat)! nir direct beam solar rad onto srf (w/m**2)
      real(r8) solsdxy(londim,beglat:endlat)! vis diffuse solar rad onto srf (w/m**2)
      real(r8) solldxy(londim,beglat:endlat)! nir diffuse solar rad onto srf(w/m**2)

! output variables

      real(r8) shxy(londim,beglat:endlat)    ! sens heat flux (w/m**2) [+to atm]
      real(r8) lhxy(londim,beglat:endlat)    ! latent heat flux (w/m**2) [+to atm]
      real(r8) cfxy(londim,pcnst+pnats,beglat:endlat)! constituent fluxes (kg/m**2/s)[+to atm]
      real(r8) tauxxy(londim,beglat:endlat)  ! zonal surf stress (kg/m/s**2)
      real(r8) tauyxy(londim,beglat:endlat)  ! merid surf stress (kg/m/s**2)
      real(r8) tsxy   (londim,beglat:endlat) ! surface radiative temperature (kelvin)
      real(r8) lwupxy (londim,beglat:endlat) ! emitted longwave radiation (w/m**2)
      real(r8) trefxy(londim,beglat:endlat)  ! 2 m height air temp (kelvin)
      real(r8) asdirxy(londim,beglat:endlat) ! albedo -- visible waveband, direct
      real(r8) asdifxy(londim,beglat:endlat) ! albedo -- visible waveband, diffuse
      real(r8) aldirxy(londim,beglat:endlat) ! albedo -- near infrared waveband, direct
      real(r8) aldifxy(londim,beglat:endlat) ! albedo -- near infrared waveband, diffuse
      real(r8) snowxy (londim,beglat:endlat) ! snow, liquid water equivalent (m)
! -----------------------------------------------------------------

! -----------------------------------------------------------------
! local variables

      integer i,j,k,m              !loop/array indices
      logical doalb                !true if surface albedo calculation time step
      real(r8) :: wt               !subgrid weight
      real(r8) :: po2  = 0.209     !partial pressure  o2 (mol/mol)
      real(r8) :: pco2 = 355.e-06  !partial pressure co2 (mol/mol)
      integer  :: kstrt            !starting subgrid index
      integer  :: kstop            !ending subgrid index

! -----------------------------------------------------------------

CC#ifdef SPMD
CC! return of no kpt points for this latitude band
CC      if (numks == 0) return
CC#endif 

! determine if albedo calculation is done

      doalb = irad.eq.1 .or.(mod(nstep,irad).eq.0 .and. nstep+1.ne.1)

! determine starting and ending subgrid indices

#if ( defined SPMD )
      kstrt = begkpt
      kstop = endkpt
#else
      kstrt = 1
      kstop = kpt
#endif

! map atmospheric fields to force lsm: [lsmlon] x [lsmlat] grid ->
! [lpt] vector of land points -> [kpt] vector of subgrid points
! note: atmospheric co2 and o2 are currently constants


      call t_startf ('lsm_single')
      do k = kstrt,kstop
        i = ixy(k)
        j = jxy(k)
        tgcm(k)    = tgcmxy(i,j)      
        hgcm(k)    = zgcmxy(i,j)     
        ugcm(k)    = ugcmxy(i,j)    
        vgcm(k)    = vgcmxy(i,j)   
        qgcm(k)    = qgcmxy(i,j)  
        pgcm(k)    = pgcmxy(i,j)      
        thgcm(k)   = thgcmxy(i,j)  
        qsnowfc(k) = snowcxy(i,j)*1000. ! convert from m/sec to mm/sec
        qsnowfl(k) = snowlxy(i,j)*1000. ! convert from m/sec to mm/sec
        qrainc(k)  = (prcxy(i,j) - snowcxy(i,j))*1000.0 ! convert from m/sec to mm/sec
        qrainl(k)  = (prlxy(i,j) - snowlxy(i,j))*1000.0 ! convert from m/sec to mm/sec
        solad(1,k) = solsxy(i,j)     
        solad(2,k) = sollxy(i,j)    
        solai(1,k) = solsdxy(i,j)      
        solai(2,k) = solldxy(i,j)     
        firgcm(k)  = flwdsxy(i,j)     
        egcm(k)    = qgcm(k)*pgcm(k) / (0.622+0.378*qgcm(k))
        rhogcm(k)  = (pgcm(k)-0.378*egcm(k)) / (rair*tgcm(k))
        co2gcm(k)  = pco2*pgcm(k)
        o2gcm(k)   = po2*pgcm(k) 
      end do
      call t_stopf ('lsm_single')

! call lsm driver  

      call lsmdrv (nstep, doalb, eccen, obliqr, lambm0, mvelpp)

! initialize 2-d fields to zero for subgrid averaging - but only 
! but only for land points because other surface models may have 
! already calculated fields for non-land points.

      call t_startf ('lsm_single')
      do j = beglat,endlat
        do i = 1,numlon(j)
          if (surf2d(i,j) > 0.) then 
            shxy(i,j)    = 0.
            lhxy(i,j)    = 0.
            cfxy(i,1,j)  = 0.
            tauxxy(i,j)  = 0.
            tauyxy(i,j)  = 0.
            tsxy(i,j)    = 0.
            trefxy(i,j)  = 0.
            lwupxy(i,j)  = 0.
            asdirxy(i,j) = 0.
            asdifxy(i,j) = 0.
            aldirxy(i,j) = 0.
            aldifxy(i,j) = 0.
            snowxy(i,j)  = 0.
          end if
        end do
        do m = 2,pcnst+pnats
          do i = 1,numlon(j)
            if (surf2d(i,j) > 0.) cfxy(i,m,j) = 0.
          end do
        end do
      end do

! map fields from subgrid vector of length [kpt] to [lsmlon]x[lsmlat] grid,
! NOTE: only use points where weight is nonzero so SPMD code will 
! not use uninitialized stack memory values for arrays like taux.
! NOTE: set cfxy = co2 fluxes for pcnst > 1, multiplying by 4.4e-08
! to convert from umol co2/m**2/s to kg co2/m**2/s

      do k=kstrt,kstop
         wt = wtk2xy(k)
         if (wt /= 0.) then
            i = ixy(k)
            j = jxy(k)
            shxy(i,j)    = shxy(i,j)   + fsh(k)*wt
            lhxy(i,j)    = lhxy(i,j)   + (fgev(k)+fcev(k)+fctr(k))*wt
            cfxy(i,1,j)  = cfxy(i,1,j) + (qsoil(k)+qvege(k)+qvegt(k))*wt
            tauxxy(i,j)  = tauxxy(i,j) + taux(k)*wt
            tauyxy(i,j)  = tauyxy(i,j) + tauy(k)*wt
            tsxy(i,j)    = tsxy(i,j)   + trad(k)*wt
            trefxy(i,j)  = trefxy(i,j) + tsa(k)*wt
            lwupxy(i,j)  = lwupxy(i,j) + fire(k)*wt
            asdirxy(i,j) = asdirxy(i,j)+ albd(1,k)*wt
            asdifxy(i,j) = asdifxy(i,j)+ albi(1,k)*wt
            aldirxy(i,j) = aldirxy(i,j)+ albd(2,k)*wt
            aldifxy(i,j) = aldifxy(i,j)+ albi(2,k)*wt
            snowxy(i,j)  = snowxy(i,j) + h2osno(k)*wt/1000. ! convert mm to m
!           cfxy(i,2,j)  = cfxy(i,2,j) + fco2(k)*wt * 4.4e-08
!           cfxy(i,3,j)  = cfxy(i,3,j) + fpsn(k)*wt * 4.4e-08
!           cfxy(i,4,j)  = cfxy(i,4,j) + (frm(k)+frg(k))*wt * 4.4e-08
!           cfxy(i,5,j)  = cfxy(i,5,j) + fmicr(k)*wt * 4.4e-08
         end if
      end do
      call t_stopf ('lsm_single')
#endif

      return
      end

 
