module soil_drv 4,3
!@sum soil_drv contains variables and routines for the ground
!@+   hydrology driver
!@auth I. Alienov/F. Abramopolous
      use model_com, only : im,jm
      use DOMAIN_DECOMP, only : GRID, GET
      use veg_drv, only : cosday,sinday
      implicit none
      private
      save

      public daily_earth, ground_e, init_gh, earth, conserv_wtg
     $     ,conserv_htg

      !real*8 cosday,sinday
      !real*8 cosdaym1, sindaym1               !nyk TEMPORARY for jday-1
      real*8 adlmass          ! accumulator for dleafmass in daily_earth

      real*8 spgsn !@var spgsn specific gravity of snow
!@dbparam snow_cover_coef coefficient for topography variance in
!@+       snow cover parameterisation for albedo
      real*8 :: snow_cover_coef = .15d0

      contains


      subroutine earth (ns,moddsf,moddd) 1,31
!@sum EARTH calculates surface fluxes of sensible heat,
!@+   evaporation, thermal radiation, and momentum drag over earth.
!@auth I. Alienov/F. Abramopolous
c****
      use constant, only : grav,rgas,lhe,lhs
     *     ,sha,tf,rhow,deltx
      use model_com, only : t,p,q,dtsrc,nisurf,dsig,qcheck,jdate
     *     ,jday,jhour,nday,itime,jeq,fearth,modrd,itearth
     *     ,u,v
      use DOMAIN_DECOMP, only : HALO_UPDATE, CHECKSUM, NORTH
      use geom, only : imaxj,dxyp,bydxyp
      use dynamics, only : pmid,pk,pek,pedn,pdsig,am,byam
      use somtq_com, only : mz
      use radncb, only : trhr,fsf, cosz1

      use surf_albedo, only: albvnh   ! added 5/23/03 from RADIATION.f
      !albvnh(9,6,2)=albvnh(sand+8veg,6bands,2hemi) - only need 1st band
      use sle001
     &    , only : advnc,evap_limits,
     &    pr,htpr,prs,htprs,w,ht,snowd,tp,fice,
     &    fv,fb,atrg,ashg,alhg,
     &    abetad,abetav,abetat,
     &    abetap,abetab,abeta,
     &    acna,acnc,agpp,
     &    aevap,aevapw,aevapd,aevapb,
     &    aruns,arunu,aeruns,aerunu,
     &    aepc,aepb,aepp,af0dt,af1dt,zw,tbcs,
     &    qm1,qs,
     &    pres,rho,ts,vsm,ch,srht,trht, !cd,snht,
     &    nsn,dzsn,wsn,hsn,fr_snow
     &     ,ghy_debug


      use veg_drv, only: veg_save_cell,veg_set_cell
      use vegetation, only: update_veg_locals

      use dagcom , only : aij,tsfrez,tdiurn,aj,areg,adiurn,jreg,hdiurn,
     *     ij_rune, ij_arunu, ij_pevap, ij_shdt, ij_beta, ij_trnfp0,
     *     ij_srtr, ij_neth, ij_ws, ij_ts, ij_us, ij_vs, ij_taus,
     *     ij_tauus, ij_tauvs, ij_qs, ij_tg1, ij_evap, j_trhdt, j_shdt,
     *     j_evhdt,j_evap,j_erun,j_run,j_tsrf,j_type,j_tg1,j_tg2,ij_g05
     *     ,ij_g06,ij_g11,ij_g12,ij_g13,ij_g14,ij_g15,ij_g16,ij_g17
     *     ,ij_gpp, ij_dleaf,ij_pblht
     *     ,ij_g18,ij_g19,ij_g20,ij_g21,ij_g22,ij_g23,ij_g24,ij_g25
     *     ,ij_g26,ij_g27,ijdd,idd_ts,idd_tg1,idd_qs,idd_qg,idd_swg
     *     ,idd_lwg,idd_sh,idd_lh,idd_hz0,idd_ug,idd_vg,idd_wg,idd_us
     *     ,idd_vs,idd_ws,idd_cia,idd_cm,idd_ch,idd_cq,idd_eds,idd_dbl
     *     ,idd_ev,tf_day1,tf_last,ndiupt

      use fluxes, only : dth1,dq1,uflux1,vflux1,e0,e1,evapor,prec,eprec
     *     ,runoe,erunoe,gtemp,precss
      use ghycom, only : ngm,nlsn,
     &     gdeep,wbare,wvege,htbare,htvege,snowbv,
     &     nsn_ij,dzsn_ij,wsn_ij,hsn_ij,fr_snow_ij,
     *     canopy_temp_ij,snowe,tearth,wearth,aiearth,
     &     evap_max_ij, fr_sat_ij, qg_ij, fr_snow_rad_ij,top_dev_ij

!#ifdef TRACERS_WATER
!     *     ,trvege,trbare,trsnowbv
!#endif
      use vegetation, only :
     &    veg_srht=>srht,veg_pres=>pres,veg_ch=>ch,veg_vsm=>vsm !ia

      USE SOCPBL, only : dtsurf         ! zgs,     ! global
     &     ,zs1,tgv,tkv,qg_sat,qg_aver,hemi,pole     ! rest local
     &     ,us,vs,ws,wsm,wsh,tsv,qsrf,psi,dbl    ! ,edvisc=>kms
     &     ,khs,ug,vg,wg,wint   ! ,kq=>kqs ,ppbl
      use pblcom, only : ipbl,cmgs,chgs,cqgs,tsavg,qsavg
      use pbl_drv, only : pbl, evap_max,fr_sat,uocean,vocean,psurf,trhr0


      use snow_drvm, only : snow_cover_same_as_rad


      implicit none

      integer, intent(in) :: ns,moddsf,moddd
      integer i,j,kr,jr,itype,ih,ihm,ibv
      real*8 shdt,qsats,evap,evhdt,tg2av,ace2av,trhdt,rcdmws,rcdhws
     *     ,cdq,cdm,cdh,elhx,tg,srheat,tg1,ptype,trheat,wtr2av    !,dhgs
     *     ,rhosrf,ma1,tfs,th1,thv1,p1k,psk,ps,pij,psoil,pearth
     *     ,warmer,timez,spring,zs1co,q1

!@var rhosrf0 estimated surface air density
      real*8 rhosrf0


      real*8 qsat
      real*8 srhdt
      real*8 aregij(9,im,jm)
!@var qg rel. humidity at the ground, defined: total_evap = Cq V (qg-qs)
!@var qg_nsat rel. humidity at non-saturated fraction of soil
      real*8 qg, qg_nsat
c****
c**** fearth    soil covered land fraction (1)
c****
c**** snowe     earth snow amount (kg/m**2)
c**** tearth    earth temperature of first layer (c)
c**** wearth    earth water of first layer (kg/m**2)
c**** aiearth   earth ice of first layer (kg/m**2)
c****
c**** wbare  1-6 water of bare soil layer 1-6 (m)
c**** wvege   0  water of vegetation canopy (m)
c****        1-6 water of vegetated soil layer 1-6 (m)
c**** htbare  0  bare soil layer 0 is unused
c****        1-6 heat content of bare soil layer 1-6 (j m-2)
c**** htvege  0  heat content of vegetation canopy (j m-2)
c****        1-6 heat content of vegetated soil layer 1-6 (j m-2)
c**** snowbv  1  snow depth over bare soil (m)
c****         2  snow depth over vegetated soil (m)
c****

C****   define local grid
      integer J_0, J_1

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT=J_0, J_STOP=J_1)

      dtsurf=dtsrc/nisurf
      zs1co=.5*dsig(1)*rgas/grav

      spring=-1.
      if((jday.ge.32).and.(jday.le.212)) spring=1.
      ih=1+jhour
      ihm = ih+(jdate-1)*24
c****
c**** outside loop over time steps, executed nisurf times every hour
c****
      timez=jday+(mod(itime,nday)+(ns-1.)/nisurf)/nday ! -1 ??
      if(jday.le.31) timez=timez+365.

ccc set i,j - independent stuff for tracers


      aregij = 0.
c****
c**** outside loop over j and i, executed once for each grid point
c****
C**** halo update u and v for distributed parallelization
       call checksum   (grid, U, 263, "GHY_DRV.f")
       call halo_update(grid, U, from=NORTH)
       call checksum   (grid, V, 265, "GHY_DRV.f")
       call halo_update(grid, V, from=NORTH)

!$OMP  PARALLEL DO PRIVATE
!$OMP*  (ACE2AV, ELHX,EVAP,EVHDT, CDM,CDH,CDQ,
!$OMP*   I,ITYPE,ibv, J, KR, MA1,PIJ,PSK,PEARTH,PSOIL,PS,P1K,PTYPE, QG,
!$OMP*   QG_NSAT,QSATS, RHOSRF,RHOSRF0,RCDMWS,RCDHWS, SRHDT,SRHEAT,SHDT,
!$OMP*   TRHEAT, TH1,TFS,THV1,TG1,TG,TRHDT,TG2AV, WARMER,WTR2AV,q1

!$OMP*   )
!$OMP*   SCHEDULE(DYNAMIC,2)

      loop_j: do j=J_0,J_1
      hemi=1.
      if(j.le.jm/2) hemi=-1.
c**** conditions at the south/north pole
      pole= ( j.eq.1 .or. j.eq.jm )

ccc   if(j.lt.jeq) warmer=-spring
ccc   if(j.ge.jeq) warmer=spring
      if(j.lt.jeq)  then
         warmer=-spring
       else
         warmer=spring
      end if
      loop_i: do i=1,imaxj(j)
c****
c**** determine surface conditions
c****
      pearth=fearth(i,j)
      psoil=pearth
      pij=p(i,j)
      ps=pedn(1,i,j)
      psk=pek(1,i,j)
      p1k=pk(1,i,j)
      th1=t(i,j,1)
      q1=q(i,j,1)
      thv1=th1*(1.+q1*deltx)
      tkv=thv1*psk
      tfs=tf*psoil
      ma1=am(1,i,j)
      qm1=q1*ma1
c     rhosrf=100.*ps/(rgas*tsv)
c     rhosrf=100.*ps/(rgas*tkv)
ccc   jr=jreg(i,j)
c****
c**** earth
c****
      if (pearth.le.0.) then
        ipbl(i,j,4)=0
        cycle loop_i
      endif


      itype=4
      ptype=pearth
      pr=prec(i,j)/(dtsrc*rhow)
C**** This variable was originally assoicated with super-saturated
C**** large-scale precip, but has not worked for many moons.
C**** If you want to reinstate it, uncomment this calculation.
c      prs=precss(i,j)/(dtsrc*rhow)
      prs=0.
      htpr=eprec(i,j)/dtsrc
!!! insert htprs here
      w(1:ngm,1) =  wbare(1:ngm,i,j)
      w(0:ngm,2) =  wvege(0:ngm,i,j)
      ht(0:ngm,1) = htbare(0:ngm,i,j)
      ht(0:ngm,2) = htvege(0:ngm,i,j)
      snowd(1:2)  = snowbv(1:2,i,j)
ccc extracting snow variables
      nsn(1:2)          = nsn_ij    (1:2, i, j)
      !isn(1:2)          = isn_ij    (1:2, i, j)
      dzsn(1:nlsn, 1:2) = dzsn_ij   (1:nlsn, 1:2, i, j)
      wsn(1:nlsn, 1:2)  = wsn_ij    (1:nlsn, 1:2, i, j)
      hsn(1:nlsn, 1:2)  = hsn_ij    (1:nlsn, 1:2, i, j)
      fr_snow(1:2)      = fr_snow_ij(1:2, i, j)
ccc tracers variables

      tg1 = tearth(i,j)
      srheat=fsf(itype,i,j)*cosz1(i,j)
! need this
      srhdt=srheat*dtsurf
c****
c**** boundary layer interaction
c****
      zs1=zs1co*tkv*pij/pmid(1,i,j)
c**** loop over ground time steps
      tg=tg1+tf
      elhx=lhe
      if(tg1.lt.0.)  elhx=lhs
      qg_sat=qsat(tg,elhx,ps)  !  replacing with qs from prev step
      qg = qg_ij(i,j)
      ! if ( qg > 999.d0 ) qg = qg_sat
      qg_aver = qg
      tgv=tg*(1.+qg*deltx)
      psurf=ps
      trhr0 = TRHR(0,I,J)
      rhosrf0=100.*ps/(rgas*tgv) ! estimated surface density
C**** Obviously there are no ocean currents for earth points, but
C**** variables set for consistency with surfce
      uocean=0 ; vocean=0

c***********************************************************************
c****
ccc actually PBL needs evap (kg/m^2*s) / rho_air
      evap_max = evap_max_ij(i,j) * 1000.d0 / rhosrf0
      fr_sat = fr_sat_ij(i,j)


      call pbl(i,j,itype,ptype)
c****
      cdm = cmgs(i,j,itype)
      cdh = chgs(i,j,itype)
      cdq = cqgs(i,j,itype)
c***********************************************************************
c**** calculate qs
      qs=qsrf
      ts=tsv/(1.+qs*deltx)
c**** calculate rhosrf*cdm*ws
      rhosrf=100.*ps/(rgas*tsv)
      rcdmws=cdm*wsm*rhosrf
      rcdhws=cdh*wsh*rhosrf
c**** calculate fluxes of sensible heat, latent heat, thermal
c****   radiation, and conduction heat (watts/m**2)
c      snht=-sha*rcdhws*(ts-tg)  ! -not used
      trheat=trhr(0,i,j)
c***********************************************************************
c****
c  define extra variables to be passed in surfc:
      pres  =ps
      veg_pres = ps
      rho   =rhosrf
      vsm   =ws
      veg_vsm = ws
      ch    =cdh
      veg_ch = cdh
      srht  =srheat
      veg_srht = srheat
      trht  =trheat
  !    zs    =zgs  !!! will not need after qsbal is replaced
  !    eddy  =khs   !!! will not need after qsbal is replaced
c     tspass=ts
c***********************************************************************
c****
c**** calculate ground fluxes
c     call qsbal

      call ghinij (i,j)
      call veg_set_cell(i,j)
      !call init_localveg
      call advnc
      call evap_limits( .false., evap_max_ij(i,j), fr_sat_ij(i,j) )

      call veg_save_cell(i,j)

      tg1=tbcs
      !qg_ij(i,j) = qs  !!! - this seemed to work ok
      !! trying more precise value for qg :  qsat(tg1+tf,elhx,ps)
      qg_sat = qsat(tg1+tf,elhx,ps) ! saturated soil
      qg_nsat = qs              ! non-sat soil, no evap
      if ( rcdhws > 1.d-30 ) then   ! correction to non-sat, due to evap
        qg_nsat = qg_nsat + evap_max_ij(i,j)/(0.001*rcdhws)
      endif
      qg_nsat = min( qg_nsat, qg_sat )
      qg_ij(i,j) = fr_sat_ij(i,j) * qg_sat
     &     + (1.d0 -fr_sat_ij(i,j)) * qg_nsat

      wbare(1:ngm,i,j) = w(1:ngm,1)
      wvege(0:ngm,i,j) = w(0:ngm,2)
      htbare(0:ngm,i,j) = ht(0:ngm,1)
      htvege(0:ngm,i,j) = ht(0:ngm,2)
      snowbv(1:2,i,j)   = snowd(1:2)

!!! test test test
!      aevap = 0.d0

cddd#ifdef TRACERS_WATER_OLD
cdddC**** reset tracer variables
cdddc      wsoil_tot=wsoil_tot-(aevapw+aevapd+aevapb+aruns+arunu)/rhow
cddd      wsoil_tot=wsoil_tot-(aevap+aruns+arunu)/rhow
cddd      do nx=1,ntx
cddd        n=ntix(nx)
cddd        if (tr_wd_TYPE(n).eq.nWATER) THEN
cdddc**** fix outputs to mean ratio (TO BE REPLACED WITHIN SOIL TRACERS)
cddd        trruns(nx)=aruns * trsoil_rat(nx) ! kg/m^2
cddd        trrunu(nx)=arunu * trsoil_rat(nx)
cdddc#ifdef TRACERS_SPECIAL_O18
cdddc        tevapw(nx)=aevapw*trsoil_rat(nx)*fracvl(tp(1,1),trname(n))
cdddc#else
cdddc        tevapw(nx)=aevapw * trsoil_rat(nx)
cdddc#endif
cdddc        tevapd(nx)=aevapd * trsoil_rat(nx)
cdddc        tevapb(nx)=aevapb * trsoil_rat(nx)
cddd        tevapw(nx)=aevap * trsoil_rat(nx)
cdddc**** update ratio
cc       trsoil_tot(nx)=trsoil_tot(nx)-(tevapw(nx)+tevapd(nx)+tevapb(nx)
cdddc     *       +trruns(nx)+trrunu(nx))
c       trsoil_tot(nx)=trsoil_tot(nx)-(tevapw(nx)+trruns(nx)+trrunu(nx))
cddd        trsoil_rat(nx)=trsoil_tot(nx)/(rhow*wsoil_tot)
cddd        trbare(n,1:ngm,i,j) = trsoil_rat(nx)*w(1:ngm,1)*rhow
cddd        trvege(n,:,i,j) = trsoil_rat(nx)*w(:,2)*rhow
cddd        trsnowbv(n,:,i,j)=trsoil_rat(nx)*snowd(:)*rhow
cdddc       trbare(n,:,i,j) = trw(nx,:,1)
cdddc       trvege(n,:,i,j) = trw(nx,:,2)
cdddc       trsnowbv(n,:,i,j) = trsnowd(nx,:)
cddd        gtracer(n,itype,i,j)=trsoil_rat(nx)
cddd        end if
cddd      end do
cddd#endif
ccc copy snow variables back to storage
      nsn_ij    (1:2, i, j)         = nsn(1:2)
      !isn_ij    (1:2, i, j)         = isn(1:2)
      dzsn_ij   (1:nlsn, 1:2, i, j) = dzsn(1:nlsn,1:2)
      wsn_ij    (1:nlsn, 1:2, i, j) = wsn(1:nlsn,1:2)
      hsn_ij    (1:nlsn, 1:2, i, j) = hsn(1:nlsn,1:2)
      fr_snow_ij(1:2, i, j)         = fr_snow(1:2)
ccc Save canopy temperature.
      canopy_temp_ij(i,j) = tp(0,2)  !nyk
ccc tracers


c**** set snow fraction for albedo computation (used by RAD_DRV.f)
      if ( snow_cover_same_as_rad == 0 ) then
        do ibv=1,2
          call snow_cover(fr_snow_rad_ij(ibv,i,j),
     &         snowbv(ibv,i,j), top_dev_ij(i,j) )
          fr_snow_rad_ij(ibv,i,j) = min (
     &         fr_snow_rad_ij(ibv,i,j), fr_snow_ij(ibv, i, j) )
        enddo
      else
        do ibv=1,2
          fr_snow_rad_ij(ibv,i,j) = fr_snow_ij(ibv, i, j)
        enddo
      endif

      aij(i,j,ij_g18)=aij(i,j,ij_g18)+aevapb
      aij(i,j,ij_g19)=aij(i,j,ij_g19)+aevapd
      aij(i,j,ij_g20)=aij(i,j,ij_g20)+aevapw
      aij(i,j,ij_g05)=aij(i,j,ij_g05)+abetab/nisurf
      aij(i,j,ij_g06)=aij(i,j,ij_g06)+abetap/nisurf
      aij(i,j,ij_g11)=aij(i,j,ij_g11)+abeta/nisurf
      aij(i,j,ij_g12)=aij(i,j,ij_g12)+acna/nisurf
      aij(i,j,ij_g13)=aij(i,j,ij_g13)+acnc/nisurf
      aij(i,j,ij_gpp)=aij(i,j,ij_gpp)+agpp
      aij(i,j,ij_g26)=aij(i,j,ij_g26)+abetav/nisurf
      aij(i,j,ij_g27)=aij(i,j,ij_g27)+abetat/nisurf
      aij(i,j,ij_g14)=aij(i,j,ij_g14)+aepp
      if (moddsf.eq.0) then
        aij(i,j,ij_g15)=aij(i,j,ij_g15)+tp(1,1)
        aij(i,j,ij_g16)=aij(i,j,ij_g16)+tp(2,1)
        aij(i,j,ij_g17)=aij(i,j,ij_g17)+tp(3,1)
        aij(i,j,ij_g21)=aij(i,j,ij_g21)+tp(0,2)
        aij(i,j,ij_g22)=aij(i,j,ij_g22)+tp(1,2)
        aij(i,j,ij_g23)=aij(i,j,ij_g23)+tp(2,2)
        aij(i,j,ij_g24)=aij(i,j,ij_g24)+tp(3,2)
        aij(i,j,ij_g25)=aij(i,j,ij_g25)+fb*zw(1)+fv*zw(2)
      end if
      trhdt=trheat*dtsurf-atrg
c           for radiation find composite values over earth
c           for diagnostic purposes also compute gdeep 1 2 3
      snowe(i,j)=1000.*(snowd(1)*fb+snowd(2)*fv)
      tearth(i,j)=tg1
      wearth(i,j)=1000.*( fb*w(1,1)*(1.-fice(1,1)) +
     &     fv*(      aiearth(i,j)=1000.*( fb*w(1,1)*fice(1,1) +
     &     fv*(      call retp2 (tg2av,wtr2av,ace2av)
      gdeep(i,j,1)=tg2av
      gdeep(i,j,2)=wtr2av
      gdeep(i,j,3)=ace2av
      gtemp(1,4,i,j)=tearth(i,j)
c**** calculate fluxes using implicit time step for non-ocean points
      uflux1(i,j)=uflux1(i,j)+ptype*rcdmws*(us-uocean)
      vflux1(i,j)=vflux1(i,j)+ptype*rcdmws*(vs-vocean)
c**** accumulate surface fluxes and prognostic and diagnostic quantities
      !evap=aevapw+aevapd+aevapb
      evap = aevap
      evapor(i,j,4)=evapor(i,j,4)+evap
      evhdt=-alhg
C**** hack to correct energy flux
ccc      evhdt=-evap*lhe ! hopefully not needed any more
      shdt=-ashg
      dth1(i,j)=dth1(i,j)-shdt*ptype/(sha*ma1*p1k)
      dq1(i,j) =dq1(i,j)+evap*ptype/ma1
      qsavg(i,j)=qsavg(i,j)+qs*ptype
      qsats=qsat(ts,elhx,ps)
c**** save runoff for addition to lake mass/energy resevoirs
      runoe (i,j)=runoe (i,j)+ aruns+ arunu
      erunoe(i,j)=erunoe(i,j)+aeruns+aerunu
c****
      e0(i,j,4)=e0(i,j,4)+af0dt
      e1(i,j,4)=e1(i,j,4)+af1dt




c****
c**** accumulate diagnostics
c****

ccc not sure about the code below. hopefully that''s what is meant above

!      print '(a,10(e12.4))','trevapor_trunoe',
!     &     trevapor(1,itype,i,j), trunoe(1,i,j)
!     &     ,aevap,atr_evap(1),aruns,arunu,atr_rnff(1),pr,trpr


      aij(i,j,ij_rune)=aij(i,j,ij_rune)+aruns
      aij(i,j,ij_arunu)=aij(i,j,ij_arunu)+arunu
      aij(i,j,ij_pevap)=aij(i,j,ij_pevap)+(aepc+aepb)

      if ( warmer >= 0 ) then
        if(ts.lt.tf) tsfrez(i,j,tf_day1)=timez
        tsfrez(i,j,tf_last)=timez
      else
        if ( tsfrez(i,j,tf_last)+.03 >= timez .and. ts >= tf )
     $       tsfrez(i,j,tf_last)=timez
      endif

      if(tg1.lt.tdiurn(i,j,1)) tdiurn(i,j,1)=tg1
      if(tg1.gt.tdiurn(i,j,2)) tdiurn(i,j,2)=tg1
      if(ts.lt.tdiurn(i,j,3)) tdiurn(i,j,3)=ts
      if(ts.gt.tdiurn(i,j,4)) tdiurn(i,j,4)=ts

c**** quantities accumulated for regions in diagj
ccc   areg(jr,j_trhdt)=areg(jr,j_trhdt)+trhdt*ptype*dxyp(j)
ccc   areg(jr,j_shdt )=areg(jr,j_shdt )+shdt*ptype*dxyp(j)
ccc   areg(jr,j_evhdt)=areg(jr,j_evhdt)+evhdt*ptype*dxyp(j)
ccc   areg(jr,j_evap )=areg(jr,j_evap )+evap*ptype*dxyp(j)
ccc   areg(jr,j_erun)=areg(jr,j_erun)+(aeruns+aerunu)*pearth*dxyp(j)
ccc   areg(jr,j_run )=areg(jr,j_run )+(aruns+arunu)*pearth*dxyp(j)
ccc   if ( moddsf == 0 )
ccc  $     areg(jr,j_tsrf )=areg(jr,j_tsrf )+(ts-tf)*ptype*dxyp(j)
      AREGIJ(1,I,J)  = trhdt*ptype*dxyp(j)
      AREGIJ(2,I,J)  = shdt*ptype*dxyp(j)
      AREGIJ(3,I,J)  = evhdt*ptype*dxyp(j)
      AREGIJ(4,I,J)  = evap*ptype*dxyp(j)
      AREGIJ(5,I,J)  = (aeruns+aerunu)*pearth*dxyp(j)
      AREGIJ(6,I,J)  = (aruns+arunu)*pearth*dxyp(j)
      if ( moddsf == 0 ) THEN
        AREGIJ(7,I,J)  = (ts-tf)*ptype*dxyp(j)
        AREGIJ(8,I,J)  = tg1    *ptype*dxyp(j)
        AREGIJ(9,I,J)  = tg2av  *ptype*dxyp(j)
      end if
c**** quantities accumulated for latitude-longitude maps in diagij
      aij(i,j,ij_shdt)=aij(i,j,ij_shdt)+shdt*ptype
      aij(i,j,ij_beta)=aij(i,j,ij_beta)+abetad/nisurf
      if(modrd.eq.0)aij(i,j,ij_trnfp0)=aij(i,j,ij_trnfp0)+trhdt*ptype
     *     /dtsrc
      aij(i,j,ij_srtr)=aij(i,j,ij_srtr)+(srhdt+trhdt)*ptype
      aij(i,j,ij_neth)=aij(i,j,ij_neth)+(srhdt+trhdt+shdt+evhdt)*ptype
      aij(i,j,ij_evap)=aij(i,j,ij_evap)+evap*ptype
      if ( moddsf == 0 ) then
        aij(i,j,ij_ws)=aij(i,j,ij_ws)+ws*ptype
        aij(i,j,ij_ts)=aij(i,j,ij_ts)+(ts-tf)*ptype
        aij(i,j,ij_us)=aij(i,j,ij_us)+us*ptype
        aij(i,j,ij_vs)=aij(i,j,ij_vs)+vs*ptype
        aij(i,j,ij_taus)=aij(i,j,ij_taus)+rcdmws*wsm*ptype
        aij(i,j,ij_tauus)=aij(i,j,ij_tauus)+rcdmws*(us-uocean)*ptype
        aij(i,j,ij_tauvs)=aij(i,j,ij_tauvs)+rcdmws*(vs-vocean)*ptype
        aij(i,j,ij_qs)=aij(i,j,ij_qs)+qs*ptype
        aij(i,j,ij_tg1)=aij(i,j,ij_tg1)+tg1*ptype
        aij(i,j,ij_pblht)=aij(i,j,ij_pblht)+dbl*ptype
chyd       aij(i,j,ij_arunu)=aij(i,j,ij_arunu)
chyd      *  +   (40.6*psoil+.72*(2.*(tss-tfs)-(qsatss-qss)*lhe/sha))
c**** quantities accumulated hourly for diagDD
      endif
      if ( moddd == 0 ) then
        do kr=1,ndiupt
          if(i.eq.ijdd(1,kr).and.j.eq.ijdd(2,kr)) then
            adiurn(ih,idd_ts,kr)=adiurn(ih,idd_ts,kr)+ts*ptype
            adiurn(ih,idd_tg1,kr)=adiurn(ih,idd_tg1,kr)+(tg1+tf)*ptype
            adiurn(ih,idd_qs,kr)=adiurn(ih,idd_qs,kr)+qs*ptype
            adiurn(ih,idd_qg,kr)=adiurn(ih,idd_qg,kr)+qg*ptype
            adiurn(ih,idd_swg,kr)=adiurn(ih,idd_swg,kr)+srhdt*ptype
            adiurn(ih,idd_lwg,kr)=adiurn(ih,idd_lwg,kr)+trhdt*ptype
            adiurn(ih,idd_sh,kr)=adiurn(ih,idd_sh,kr)+shdt*ptype
            adiurn(ih,idd_lh,kr)=adiurn(ih,idd_lh,kr)+evhdt*ptype
            adiurn(ih,idd_hz0,kr)=adiurn(ih,idd_hz0,kr)
     *           +(srhdt+trhdt+shdt+evhdt)*ptype
            adiurn(ih,idd_ug,kr)=adiurn(ih,idd_ug,kr)+ug*ptype
            adiurn(ih,idd_vg,kr)=adiurn(ih,idd_vg,kr)+vg*ptype
            adiurn(ih,idd_wg,kr)=adiurn(ih,idd_wg,kr)+wg*ptype
            adiurn(ih,idd_us,kr)=adiurn(ih,idd_us,kr)+us*ptype
            adiurn(ih,idd_vs,kr)=adiurn(ih,idd_vs,kr)+vs*ptype
            adiurn(ih,idd_ws,kr)=adiurn(ih,idd_ws,kr)+ws*ptype
            adiurn(ih,idd_cia,kr)=adiurn(ih,idd_cia,kr)+psi*ptype
            adiurn(ih,idd_cm,kr)=adiurn(ih,idd_cm,kr)+cdm*ptype
            adiurn(ih,idd_ch,kr)=adiurn(ih,idd_ch,kr)+cdh*ptype
            adiurn(ih,idd_cq,kr)=adiurn(ih,idd_cq,kr)+cdq*ptype
            adiurn(ih,idd_eds,kr)=adiurn(ih,idd_eds,kr)+khs*ptype
            adiurn(ih,idd_dbl,kr)=adiurn(ih,idd_dbl,kr)+dbl*ptype
            adiurn(ih,idd_ev,kr)=adiurn(ih,idd_ev,kr)+evap*ptype
            hdiurn(ihm,idd_ts,kr)=hdiurn(ihm,idd_ts,kr)+ts*ptype
            hdiurn(ihm,idd_tg1,kr)=hdiurn(ihm,idd_tg1,kr)+(tg1+tf)*ptype
            hdiurn(ihm,idd_qs,kr)=hdiurn(ihm,idd_qs,kr)+qs*ptype
            hdiurn(ihm,idd_qg,kr)=hdiurn(ihm,idd_qg,kr)+qg*ptype
            hdiurn(ihm,idd_swg,kr)=hdiurn(ihm,idd_swg,kr)+srhdt*ptype
            hdiurn(ihm,idd_lwg,kr)=hdiurn(ihm,idd_lwg,kr)+trhdt*ptype
            hdiurn(ihm,idd_sh,kr)=hdiurn(ihm,idd_sh,kr)+shdt*ptype
            hdiurn(ihm,idd_lh,kr)=hdiurn(ihm,idd_lh,kr)+evhdt*ptype
            hdiurn(ihm,idd_hz0,kr)=hdiurn(ihm,idd_hz0,kr)
     *           +(srhdt+trhdt+shdt+evhdt)*ptype
            hdiurn(ihm,idd_ug,kr)=hdiurn(ihm,idd_ug,kr)+ug*ptype
            hdiurn(ihm,idd_vg,kr)=hdiurn(ihm,idd_vg,kr)+vg*ptype
            hdiurn(ihm,idd_wg,kr)=hdiurn(ihm,idd_wg,kr)+wg*ptype
            hdiurn(ihm,idd_us,kr)=hdiurn(ihm,idd_us,kr)+us*ptype
            hdiurn(ihm,idd_vs,kr)=hdiurn(ihm,idd_vs,kr)+vs*ptype
            hdiurn(ihm,idd_ws,kr)=hdiurn(ihm,idd_ws,kr)+ws*ptype
            hdiurn(ihm,idd_cia,kr)=hdiurn(ihm,idd_cia,kr)+psi*ptype
            hdiurn(ihm,idd_cm,kr)=hdiurn(ihm,idd_cm,kr)+cdm*ptype
            hdiurn(ihm,idd_ch,kr)=hdiurn(ihm,idd_ch,kr)+cdh*ptype
            hdiurn(ihm,idd_cq,kr)=hdiurn(ihm,idd_cq,kr)+cdq*ptype
            hdiurn(ihm,idd_eds,kr)=hdiurn(ihm,idd_eds,kr)+khs*ptype
            hdiurn(ihm,idd_dbl,kr)=hdiurn(ihm,idd_dbl,kr)+dbl*ptype
            hdiurn(ihm,idd_ev,kr)=hdiurn(ihm,idd_ev,kr)+evap*ptype
          end if
        end do
      endif
c**** quantities accumulated for surface type tables in diagj
      aj(j,j_evap ,itearth)=aj(j,j_evap ,itearth)+ evap*pearth
      aj(j,j_trhdt,itearth)=aj(j,j_trhdt,itearth)+trhdt*pearth
      aj(j,j_evhdt,itearth)=aj(j,j_evhdt,itearth)+evhdt*pearth
      aj(j,j_shdt ,itearth)=aj(j,j_shdt ,itearth)+ shdt*pearth
      aj(j,j_erun ,itearth)=aj(j,j_erun ,itearth)+(aeruns+aerunu)*pearth
      aj(j,j_run  ,itearth)=aj(j,j_run  ,itearth)+(aruns+arunu)*pearth
      if(moddsf.eq.0) then
        aj(j,j_tsrf,itearth)=aj(j,j_tsrf,itearth)+(ts-tf)*pearth
        aj(j,j_tg1 ,itearth)=aj(j,j_tg1 ,itearth)+    tg1*pearth
        aj(j,j_tg2 ,itearth)=aj(j,j_tg2 ,itearth)+  tg2av*pearth
        aj(j,j_type,itearth)=aj(j,j_type,itearth)+        pearth
      end if

      end do loop_i
      end do loop_j
!$OMP  END PARALLEL DO

      DO 825 J=J_0,J_1
      DO 825 I=1,IMAXJ(J)
         IF(FEARTH(I,J).LE.0.0)  GO TO 825
         JR=JREG(I,J)
         areg(jr,j_trhdt)=areg(jr,j_trhdt)+AREGIJ(1,I,J)
         areg(jr,j_shdt )=areg(jr,j_shdt )+AREGIJ(2,I,J)
         areg(jr,j_evhdt)=areg(jr,j_evhdt)+AREGIJ(3,I,J)
         areg(jr,j_evap )=areg(jr,j_evap )+AREGIJ(4,I,J)
         areg(jr,j_erun )=areg(jr,j_erun )+AREGIJ(5,I,J)
         areg(jr,j_run  )=areg(jr,j_run  )+AREGIJ(6,I,J)
         if( moddsf == 0 ) then
           areg(jr,j_tsrf)=areg(jr,j_tsrf)+AREGIJ(7,I,J)
           areg(jr,j_tg1 )=areg(jr,j_tg1 )+AREGIJ(8,I,J)
           areg(jr,j_tg2 )=areg(jr,j_tg2 )+AREGIJ(9,I,J)
         end if
  825 CONTINUE
C
      return
      end subroutine earth


      subroutine snow_cover( fract_snow, snow_water, top_dev ) 2,1
!@sum computes snow cover from snow water eq. and topography
!@var fract_snow snow cover fraction (0-1)
!@var snow_water snow water equivalent (m)
!@var top_dev standard deviation of the surface elevation
      use constant, only : teeny
      real*8, intent(out) :: fract_snow
      real*8, intent(in) :: snow_water, top_dev

      ! using formula from the paper by A. Roesch et al
      ! (Climate Dynamics (2001), 17: 933-946)
      fract_snow =
ccc     $     .95d0 * tanh( 100.d0 * snow_water ) *
ccc                               currently using only topography part
     $     sqrt ( 1000.d0 * snow_water /
     $     (1000.d0 * snow_water + teeny + snow_cover_coef * top_dev) )

      end subroutine snow_cover



      subroutine init_gh(dtsurf,redogh,inisnow,istart) 2,33
c**** modifications needed for split of bare soils into 2 types
      use filemanager
      use param
      use constant, only : twopi,rhow,edpery,sha,lhe,tf
      use model_com, only : fearth,itime,nday,jeq,jyear
      use dagcom, only : npts,icon_wtg,icon_htg,conpt0
      use sle001

      use fluxes, only : gtemp
      use ghycom
      use dynamics, only : pedn
      use snow_drvm, only : snow_cover_coef2=>snow_cover_coef
     &     ,snow_cover_same_as_rad
      use veg_drv, only : init_vegetation
      use veg_com, only : vdata

      implicit none

      real*8, intent(in) :: dtsurf
      integer, intent(in) :: istart
      logical, intent(in) :: redogh, inisnow
      integer iu_soil,iu_top_index
      integer jday
      real*8 snowdp,wtr1,wtr2,ace1,ace2,tg1,tg2
      logical :: qcon(npts)
      integer i, j, ibv
      real*8 pearth
      logical ghy_data_missing
      character conpt(npts)*10

c****
c**** contents of sdata(i,j,k):
c****       1 -   ngm   dz(ngm)
c****   ngm+1 - 6*ngm   q(is,ngm)
c**** 6*ngm+1 - 11*ngm   qk(is,ngm)
c**** 11*ngm+1           sl
      real*8, external :: qsat
!@dbparam ghy_default_data if == 1 reset all GHY data to defaults
!@+ (do not read it from files)
      integer :: ghy_default_data = 0

C**** define local grid
      integer J_0, J_1

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT=J_0, J_STOP=J_1)

c**** set conservation diagnostics for ground water mass and energy
      conpt=conpt0
      conpt(4)="EARTH"
      qcon=(/ .false., .false., .false., .true., .false., .false.
     $     , .false., .false., .true., .false., .false./)
      call set_con(qcon,conpt,"GRND WTR","(kg/m^2)        ",
     *     "(10^-9 kg/s/m^2)",1d0,1d9,icon_wtg)
      qcon=(/ .false., .false., .false., .true., .false., .false.
     $     , .false., .false., .true., .false., .false./)
      call set_con(qcon,conpt,"GRND ENG","(10**6 J/m^2)   ",
     *     "(10^-3 W/m^2)   ",1d-6,1d3,icon_htg)

c**** read rundeck parameters
      call sync_param( "snow_cover_coef", snow_cover_coef )
! hack. snow_cover_coef should be moved to snow_drvm
      snow_cover_coef2 = snow_cover_coef
      call sync_param( "snow_cover_same_as_rad", snow_cover_same_as_rad)
      call sync_param( "snoage_def", snoage_def )
      call sync_param( "ghy_default_data", ghy_default_data )

c**** read land surface parameters or use defaults
      if ( ghy_default_data == 0 ) then ! read from files

        !!!if (istart.le.0) return ! avoid reading unneeded files
c**** read soils parameters
        call openunit("SOIL",iu_SOIL,.true.,.true.)
        call dread (iu_SOIL,dz_ij,im*jm*(11*ngm+1),dz_ij)
        call closeunit (iu_SOIL)
c**** read topmodel parameters
        call openunit("TOP_INDEX",iu_TOP_INDEX,.true.,.true.)
        call readt(iu_TOP_INDEX,0,top_index_ij,im*jm,top_index_ij,1)
        call readt(iu_TOP_INDEX,0,top_dev_ij,  im*jm,top_dev_ij,  1)
        call closeunit (iu_TOP_INDEX)
      else  ! reset to default data
        if ( istart>0 .and. istart<10 ) then ! reset all
          call reset_gh_to_defaults( .true. )
        else   ! do not reset ghy prognostic variables
          call reset_gh_to_defaults( .false. )
        endif
        !!!if (istart.le.0) return
      endif


      if (istart.le.0) return
c****
c**** initialize constants
c****
c**** time step for ground hydrology
      dt=dtsurf
c spgsn is the specific gravity of snow
      spgsn=.1d0

c**** check whether ground hydrology data exist at this point.
      ghy_data_missing = .false.
      do j=J_0,J_1
        do i=1,im
          if (fearth(i,j).gt.0) then
            if ( top_index_ij(i,j).eq.-1. ) then
              print *,"No top_index data: i,j=",i,j,top_index_ij(i,j)
              ghy_data_missing = .true.
            end if
            if ( sum(dz_ij(i,j,1:ngm)).eq.0 ) then
              print *, "No soil data: i,j=",i,j,dz_ij(i,j,1:ngm)
              ghy_data_missing = .true.
            endif
            if (wbare(1,i,j) < 1.d-10 .and. wvege(1,i,j) < 1.d-10) then
              print*,"No gh data in restart file: i,j=",i,j,
     &             wbare(:,i,j),wvege(:,i,j)
              ghy_data_missing = .true.
            endif
          end if
        enddo
      enddo
      if ( ghy_data_missing ) then
        write(6,*) 'Ground Hydrology data is missing at some pts'
        write(6,*) 'If you have a non-standard land mask, please'
        write(6,*) 'consider using extended GH data and rfs file.'
        call stop_model(
     &       'Ground Hydrology data is missing at some cells',255)
      endif

ccc read and initialize vegetation here
      call init_vegetation(redogh,istart)

      ! no need to continue computations for postprocessing
      if (istart.le.0) return

      call hl0

c****
c      print *,' '
c      print *,'soils parameters'
c      sdstnc=100.
c      print *,'interstream distance (m) sdstnc:',sdstnc
c      c1=90.
c      print *,'canopy conductance related parameter c1:',c1
c      prfr=.1d0
c      print *,'fraction (by area) of precipitation prfr:',prfr
c      print *,' '
c****
c code transplanted from subroutine input
c**** recompute ground hydrology data if necessary (new soils data)
      if (redogh) then
        jday=1+mod(itime/nday,365)
        cosday=cos(twopi/edpery*jday)
        sinday=sin(twopi/edpery*jday)

        do j=J_0,J_1
        do i=1,im
          pearth=fearth(i,j)
          if(pearth.le.0.) then

            wbare(:,i,j)=0.
            wvege(:,i,j)=0.
            htbare(:,i,j)=0.
            htvege(:,i,j)=0.
            snowbv(:,i,j)=0.

          else
ccc   ??? remove next 5 lines? -check the old version
            w(1:ngm,1) =   wbare(1:ngm,i,j)
            w(0:ngm,2) =   wvege(0:ngm,i,j)
            ht(0:ngm,1) = htbare(0:ngm,i,j)
            ht(0:ngm,2) = htvege(0:ngm,i,j)
            snowd(1:2) =  snowbv(1:2,i,j)

c**** compute soil heat capacity and ground water saturation gws
            call ghinij (i,j)
c**** fill in soils common blocks
            snowdp=snowe(i,j)/rhow
            wtr1=wearth(i,j)
            ace1=aiearth(i,j)
            tg1 =tearth(i,j)
            wtr2=wtr1
            ace2=ace1
            tg2 =tg1
c           wtr2=gdata(i,j,9)   ! this cannot be right
c           ace2=gdata(i,j,10)
c           tg2 =gdata(i,j,8)
            call ghinht (snowdp, tg1,tg2, wtr1,wtr2, ace1,ace2)

c**** copy soils prognostic quantities to model variables
            wbare(1:ngm,i,j) = w(1:ngm,1)
            wvege(0:ngm,i,j) = w(0:ngm,2)
            htbare(0:ngm,i,j) = ht(0:ngm,1)
            htvege(0:ngm,i,j) = ht(0:ngm,2)
            snowbv(1:2,i,j)   = snowd(1:2)
          end if
        end do
        end do
        write (*,*) 'ground hydrology data was made from ground data'
      end if
c**** set gtemp array
      do j=J_0,J_1
        do i=1,im
          if (fearth(i,j).gt.0) then
            gtemp(1,4,i,j)=tearth(i,j)
          end if
        end do
      end do

C GISS-ESMF EXCEPTIONAL CASE
C-BMP Global sum on evap_max_ij

ccc if not initialized yet, set evap_max_ij, fr_sat_ij, qg_ij
ccc to something more appropriate
      if ( sum(evap_max_ij(:,:)) > im*jm-1.d0 ) then ! old default
        do j=J_0,J_1
          do i=1,im
            if ( fearth(i,j) .le. 0.d0 ) cycle
            qg_ij(i,j) = qsat(tearth(i,j)+tf,lhe,pedn(1,i,j))
          enddo
        enddo
        fr_sat_ij(:,:) = 0.d0
        evap_max_ij(:,:) = 0.d0
      endif

ccc   init snow here
ccc hope this is the right place to split first layer into soil
ccc and snow  and to set snow arrays
ccc!!! this should be done only when restarting from an old
ccc!!! restart file (without snow model data)

      if (inisnow) then
        do j=J_0,J_1
        do i=1,im
          pearth=fearth(i,j)
          if(pearth.le.0.) then
            nsn_ij(:,i,j)     = 0
            !isn_ij(:,i,j)     = 0
            dzsn_ij(:,:,i,j)  = 0.
            wsn_ij(:,:,i,j)   = 0.
            hsn_ij(:,:,i,j)   = 0.
            fr_snow_ij(:,i,j) = 0.
          else
            jday=1+mod(itime/nday,365)
            cosday=cos(twopi/edpery*jday)
            sinday=sin(twopi/edpery*jday)

            w(1:ngm,1) =   wbare(1:ngm,i,j)
            w(0:ngm,2) =   wvege(0:ngm,i,j)
            ht(0:ngm,1) = htbare(0:ngm,i,j)
            ht(0:ngm,2) = htvege(0:ngm,i,j)
            snowd(1:2) =  snowbv(1:2,i,j)

            call ghinij (i,j)
            call set_snow

            nsn_ij    (1:2, i, j)         = nsn(1:2)
            !isn_ij    (1:2, i, j)         = isn(1:2)
            dzsn_ij   (1:nlsn, 1:2, i, j) = dzsn(1:nlsn,1:2)
            wsn_ij    (1:nlsn, 1:2, i, j) = wsn(1:nlsn,1:2)
            hsn_ij    (1:nlsn, 1:2, i, j) = hsn(1:nlsn,1:2)
            fr_snow_ij(1:2, i, j)         = fr_snow(1:2)

c****     copy soils prognostic quantities to model variables
             wbare(1:ngm,i,j) = w(1:ngm,1)
             wvege(0:ngm,i,j) = w(0:ngm,2)
            htbare(0:ngm,i,j) = ht(0:ngm,1)
            htvege(0:ngm,i,j) = ht(0:ngm,2)
            snowbv(1:2,i,j)   = snowd(1:2)

          end if
        end do
        end do
      end if

c**** set snow fraction for albedo computation (used by RAD_DRV.f)
      fr_snow_rad_ij(:,:,:) = 0.d0
      do j=J_0,J_1
        do i=1,im
          if ( fearth(i,j) > 0.d0 ) then
            do ibv=1,2
              call snow_cover(fr_snow_rad_ij(ibv,i,j),
     &             snowbv(ibv,i,j), top_dev_ij(i,j) )
              fr_snow_rad_ij(ibv,i,j) = min (
     &             fr_snow_rad_ij(ibv,i,j), fr_snow_ij(ibv, i, j) )
            enddo
          endif
        enddo
      enddo



      return
      end subroutine init_gh



      subroutine reset_gh_to_defaults( reset_prognostic ) 2,4
      !use model_com, only: vdata
      use ghycom
      use veg_drv, only : reset_veg_to_defaults
      logical, intent(in) :: reset_prognostic
      integer i,j

C**** define local grid
      integer J_0, J_1

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT=J_0, J_STOP=J_1)

ccc ugly, should fix later
      call reset_veg_to_defaults( reset_prognostic )

      do j=J_0,J_1
      do i=1,im

      dz_ij(i,j,1:ngm)= (/  0.99999964d-01,  0.17254400d+00,
     &     0.29771447d+00,  0.51368874d+00,  0.88633960d+00,
     &     0.15293264d+01 /)
      q_ij(i,j,1:imt,1:ngm)=
     &     reshape(     &     0.13549370d+00,  0.00000000d+00,  0.00000000d+00,
     &     0.32995611d+00,  0.52192056d+00,  0.14812243d+00,
     &     0.00000000d+00,  0.00000000d+00,  0.32145596d+00,
     &     0.48299056d+00,  0.19555295d+00,  0.00000000d+00,
     &     0.00000000d+00,  0.47638881d+00,  0.40400982d+00,
     &     0.11959970d+00,  0.00000000d+00,  0.00000000d+00,
     &     0.99985123d-01,  0.95771909d-01,  0.41175738d-01,
     &     0.00000000d+00,  0.76306665d+00,  0.00000000d+00,
     &     0.00000000d+00,  0.00000000d+00,  0.00000000d+00,
     &     0.10000000d+01 /), (/imt,ngm/) )
      qk_ij(i,j,1:imt,1:ngm)=
     &     reshape(     &     0.12878728d+00,  0.00000000d+00,  0.00000000d+00,
     &     0.32943058d+00,  0.52857041d+00,  0.14199871d+00,
     &     0.00000000d+00,  0.00000000d+00,  0.30698991d+00,
     &     0.52528000d+00,  0.16772974d+00,  0.00000000d+00,
     &     0.00000000d+00,  0.39890009d+00,  0.43742162d+00,
     &     0.16367787d+00,  0.00000000d+00,  0.00000000d+00,
     &     0.46536058d+00,  0.39922065d+00,  0.13541836d+00,
     &     0.00000000d+00,  0.00000000d+00,  0.00000000d+00,
     &     0.00000000d+00,  0.00000000d+00,  0.00000000d+00,
     &     0.10000000d+01 /), (/imt,ngm/) )
      sl_ij(i,j)= 0.22695422d+00
      top_index_ij(i,j)= 0.10832934d+02
      top_dev_ij(i,j)= 0.21665636d+03

      if ( .not. reset_prognostic ) cycle

      snowe(i,j)= 0.65458111d-01
      tearth(i,j)= -0.12476520d+00
      wearth(i,j)=  0.29203081d+02
      aiearth(i,j)=  0.93720329d-01
      wbare(:,i,j) = (/  0.17837750d-01,  0.40924843d-01,
     &     0.77932012d-01,  0.11919649d+00,  0.57237469d-01,
     &     0.10000000d-11 /)
      wvege(:,i,j) = (/  0.10000000d-11,  0.29362259d-01,
     &     0.50065177d-01,  0.82533140d-01,  0.10383620d+00,
     &     0.31552459d-01,  0.10000000d-11 /)
      htbare(:,i,j)= (/  0.00000000d+00, -0.15487181d+07,
     &     -0.50720067d+07,  0.18917623d+07,  0.77174974d+07,
     &     0.21716040d+08,  0.44723067d+08 /)
      htvege(:,i,j)= (/ -0.13991376d+05, -0.53165599d+05,
     &     0.65443775d+06,  0.29276050d+07,  0.81455096d+07,
     &     0.21575081d+08,  0.45952255d+08 /)
      snowbv(:,i,j)= (/  0.00000000d+00,  0.65458111d-04 /)

      enddo
      enddo

      end subroutine reset_gh_to_defaults



      subroutine ghinij (i0,j0) 4,6
c**** input:
c**** avh(i,j) - array of vegetation heights
c**** spgsn - specific gravity of snow
c**** output:
c**** vh - vegetation height
c**** snowm - snow masking depth
c**** wfcap - water field capacity of top soil layer, m
c****
      use snow_model, only : i_earth,j_earth
      use sle001, only : dz,qk,ng,zb,zc,q,sl,xklh0 !spgsn,
     *     ,fb,fv,prs,ijdebug,n
     *     ,thets,thetm,ws,thm,nth,shc,shw,htprs,pr !shcap,shtpr,
     *     ,htpr
     *     ,top_index,top_stdev
      use ghycom, only : ngm,imt,dz_ij,sl_ij,q_ij,qk_ij
     *     ,top_index_ij,top_dev_ij
      use veg_com, only: afb
!      use veg_drv, only : veg_set_cell


      implicit none
      integer i0,j0
!      real*8 wfcap
      integer k,ibv,i
      real*8 shtpr
!----------------------------------------------------------------------!
      real*8, parameter :: shcap(imt) = (/2d6,2d6,2d6,2.5d6,2.4d6/)


      ijdebug=i0*1000+j0
      i_earth = i0
      j_earth = j0

ccc setting vegetation
 !     call veg_set_cell(i0,j0)

ccc passing topmodel parameters
      top_index = top_index_ij(i0, j0)
      top_stdev = top_dev_ij(i0, j0)
c**** set up layers
      dz(1:ngm)=dz_ij(i0,j0,1:ngm)
      q(1:imt,1:ngm)=q_ij(i0,j0,1:imt,1:ngm)
      qk(1:imt,1:ngm)=qk_ij(i0,j0,1:imt,1:ngm)
      sl=sl_ij(i0,j0)

      n=0
      do k=1,ngm
        if(dz(k).le.0.) exit
        n=k                
      end do
          
      if(n.le.0) then
         write (99,*) 'ghinij:  n <= 0:  i,j,n=',i0,j0,n,(dz(k),k=1,43)
         call stop_model('stopped in GHY_DRV.f',255)
      end if

c**** calculate the boundaries, based on the thicknesses.
      zb(1)=0.
      do k=1,n
        zb(k+1)=zb(k)-dz(k)
      end do
c**** calculate the layer centers, based on the boundaries.
      do k=1,n
        zc(k)=.5*(zb(k)+zb(k+1))
      end do
c**** fb,fv: bare, vegetated fraction (1=fb+fv)
      fb=afb(i0,j0)
      fv=1.-fb

c****
      do ibv=1,2
        do k=1,n
          thets(k,ibv)=0.
          thetm(k,ibv)=0.
          do i=1,imt-1
            thets(k,ibv)=thets(k,ibv)+q(i,k)*thm(0,i)
            thetm(k,ibv)=thetm(k,ibv)+q(i,k)*thm(nth,i)
          end do
          ws(k,ibv)=thets(k,ibv)*dz(k)
        end do
      end do
!veg      ws(0,2)=.0001d0*alai  ! from call veg_set_cell above
  !    wfcap=fb*ws(1,1)+fv*(ws(0,2)+ws(1,2))
c****
      call xklh0
c****
      do ibv=1,2
        do k=1,n
          shc(k,ibv)=0.
          do i=1,imt
            shc(k,ibv)=shc(k,ibv)+q(i,k)*shcap(i)
          end do
          shc(k,ibv)=(1.-thets(k,ibv))*shc(k,ibv)*dz(k)
        end do
      end do
c****
c shc(0,2) is the heat capacity of the canopy
!veg      aa=ala(1,i0,j0)
!veg      shc(0,2)=(.010d0+.002d0*aa+.001d0*aa**2)*shw
c****
c htpr is the heat of precipitation.
c shtpr is the specific heat of precipitation.
      shtpr=0.
      if(pr.gt.0.)shtpr=htpr/pr
c htprs is the heat of large scale precipitation
      htprs=shtpr*prs
c****
      return
      end subroutine ghinij


      subroutine ghinht (snowdp,tg1,tg2,wtr1,wtr2,ace1,ace2) 1,1
c**** initializes new ground (w,ht,snw) from old (t,w,ice,snw)
c**** evaluates the heat in the soil layers based on the
c**** temperatures.
c**** input:
c**** w - water in soil layers, m
c**** tp - temperature of layers, c
c**** fice - fraction of ice of layers
c**** fsn - heat of fusion of water
c**** shc - specific heat capacity of soil
c**** shi - specific heat capacity of ice
c**** shw - specific heat capcity of water
c**** snowd - snow depth, equivalent water m
c**** output:
c**** ht - heat in soil layers
c**** add calculation of wfc2
c**** based on combination of layers 2-n, as in retp2
      use sle001
      implicit none

      real*8 snowdp,tg1,tg2,wtr1,wtr2,ace1,ace2
      real*8 wfc1, wfc2, wet1, wet2, wmin, fbv
      integer k, ibv, ll

      wfc1=fb*ws(1,1)+fv*(ws(0,2)+ws(1,2))
      wfc2=0.
      fbv=fb
      do 30 ibv=1,2
      do 20 k=2,n
      wfc2=wfc2+fbv*ws(k,ibv)
   20 continue
      fbv=fv
   30 continue
      wfc1=1000.*wfc1
      wfc2=1000.*wfc2
      fice(0,2)=1.
      fice(1,1)=(ace1+snowdp*1000.)/(wtr1+ace1+snowdp*1000.+1.d-20)
      fice(1,2)=fice(1,1)
      tp(0,2)=tg1
c**** w = snow(if top layer) + wmin + (wmax-wmin)*(wtr+ice)/wfc
      w(0,2)=0.
      do ibv=1,2
        w(1,ibv)=snowdp
        wmin=thetm(1,ibv)*dz(1)
        wet1=(wtr1+ace1)/(wfc1+1.d-20)
        if(wet1.gt.1.) wet1=1.
        w(1,ibv)=w(1,ibv)+wmin+(ws(1,ibv)-wmin)*wet1
        snowd(ibv)=snowdp
        tp(1,ibv)=tg1
        do k=2,n
          fice(k,ibv)=ace2/(wtr2+ace2+1.d-20)
          wmin=thetm(k,ibv)*dz(k)
          wet2=(wtr2+ace2)/(wfc2+1.d-20)
          if(wet2.gt.1.) wet2=1.
          w(k,ibv)=wmin+(ws(k,ibv)-wmin)*wet2
          tp(k,ibv)=tg2
        end do
      end do
c****

      entry ghexht
c****
c**** compute ht (heat w/m+2)
      do ibv=1,2
        ll=2-ibv
        do k=ll,n
          if(tp(k,ibv)) 2,4,6
 2        ht(k,ibv)=tp(k,ibv)*(shc(k,ibv)+w(k,ibv)*shi)-w(k,ibv)*fsn
          cycle
 4        ht(k,ibv)=-fice(k,ibv)*w(k,ibv)*fsn
          cycle
 6        ht(k,ibv)=tp(k,ibv)*(shc(k,ibv)+w(k,ibv)*shw)
        end do
      end do
      if(ijdebug.eq.0)then
       write(99,*)'ghinht id check',ijdebug
       write(99,*)'tg1,tg2',tg1,tg2
       write(99,*)'tp',tp
       write(99,*)'ht',ht
       write(99,*)'w',w
       write(99,*)'wtr1,wtr2',wtr1,wtr2
       write(99,*)'ace1,ace2',ace1,ace2
       write(99,*)'wfc1,wfc2',wfc1,wfc2
       write(99,*)'shc',shc
       write(99,*)'fice',fice
      endif
      return
      end subroutine ghinht


      subroutine retp2 (tg2av,wtr2av,ace2av) 1,1
c**** evaluates the mean temperature in the soil layers 2-ngm
c**** as well as the water and ice content.
c**** input:
c**** w - water in soil layers, m
c**** ht - heat in soil layers
c**** fsn - heat of fusion of water
c**** shc - specific heat capacity of soil
c**** shi - specific heat capacity of ice
c**** shw - specific heat capcity of water
c**** output:
c**** tg2av - temperature of layers 2 to ngm, c
c**** ice2av - ice amount in layers 2 to ngm, kg/m+2
c**** wtr2av - water in layers 2 to ngm, kg/m+2
      use sle001
      implicit none
      real*8 tg2av,wtr2av,ace2av, wc,htc,shcc,tpc,ficec,ftp
      integer k, ibv
      tg2av=0.
      wtr2av=0.
      ace2av=0.
      do 3500 ibv=1,2
      wc=0.
      htc=0.
      shcc=0.
      do k=2,n
        wc=wc+w(k,ibv)
        htc=htc+ht(k,ibv)
        shcc=shcc+shc(k,ibv)
      end do
      tpc=0.
      ficec=0.
      if(wc.ne.0.)  ficec=-htc/(fsn*wc)
      if(fsn*wc+htc.ge.0.)go to 3430
      tpc=(htc+wc*fsn)/(shcc+wc*shi)
      ficec=1.
      go to 3440
 3430 if(htc.le.0.) go to 3440
      tpc=htc/(shcc+wc*shw)
      ficec=0.
 3440 continue
      ftp=fb
      if(ibv.eq.2) ftp=fv
      tg2av=tg2av+tpc*ftp
      wtr2av=wtr2av+wc*ftp*1000.*(1.-ficec)
      ace2av=ace2av+wc*ftp*1000.*ficec
 3500 continue
      return
      end subroutine retp2


      subroutine checke(subr),9
!@sum  checke checks whether arrays are reasonable over earth
!@auth original development team
!@ver  1.0
      use model_com, only : fearth,itime,wfcs
      use geom, only : imaxj
      use ghycom, only : tearth,wearth,aiearth,snowe,wbare,wvege,htbare
     *     ,htvege,snowbv,ngm
      implicit none

      real*8 x,tgl,wtrl,acel
      integer i,j
!@var subr identifies where check was called from
      character*6, intent(in) :: subr

C**** define local grid
      integer I_0, I_1
      integer J_0, J_1

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, I_STRT=I_0, I_STOP=I_1, J_STRT=J_0, J_STOP=J_1)

c**** check for nan/inf in earth data
      call check3(wbare(1:ngm,I_0:I_1,J_0:J_1) ,ngm  ,
     *                 (I_1-I_0+1),(J_1-J_0+1),subr,'wb')
      call check3(wvege(0:ngm,I_0:I_1,J_0:J_1) ,ngm+1,
     *                 (I_1-I_0+1),(J_1-J_0+1),subr,'wv')
      call check3(htbare(0:ngm,I_0:I_1,J_0:J_1),ngm+1,
     *                 (I_1-I_0+1),(J_1-J_0+1),subr,'hb')
      call check3(htvege(0:ngm,I_0:I_1,J_0:J_1),ngm+1,
     *                 (I_1-I_0+1),(J_1-J_0+1),subr,'hv')
      call check3(snowbv(1:ngm,I_0:I_1,J_0:J_1),2    ,
     *                 (I_1-I_0+1),(J_1-J_0+1),subr,'sn')

c**** check for reasonable temperatures over earth
      x=1.001
      do j=J_0,J_1
        do i=1,imaxj(j)
          if (fearth(i,j).gt.0.) then
            tgl=tearth(i,j)
            wtrl=wearth(i,j)
            acel=aiearth(i,j)
            if ((tgl+60.)*(60.-tgl).le.0.) write (6,901) subr,i,j,itime
     *           ,fearth(i,j),'tg1 ',snowe(i,j),tgl,wtrl,acel
            if (wtrl.lt.0..or.acel.lt.0..or.(wtrl+acel).gt.x*wfcs(i
     *           ,j)) write(6,901) subr,i,j,itime,fearth(i,j),'wtr '
     *           ,snowe(i,j),tgl,wtrl,acel,wfcs(i,j)
          end if
        end do
      end do

      return
 901  format ('0gdata off, subr,i,j,i-time,pearth,',a7,2i4,i10,f5.2,1x
     *     ,a4/' snw,x,tg1,wtr1,ice1, wfc1 ',6f12.4)

      end subroutine checke


      subroutine daily_earth(end_of_day) 2,17
!@sum  daily_earth performs daily tasks for earth related functions
!@auth original development team
!@ver  1.0
!@calls RDLAI
      use constant, only : rhow,twopi,edpery,tf
      use model_com, only : nday,nisurf,jday,jyear,fearth,wfcs
      use veg_com, only : vdata                 !nyk
      use geom, only : imaxj
      use dagcom, only : aij,tdiurn,ij_strngts,ij_dtgdts,ij_tmaxe
     *     ,ij_tdsl,ij_tmnmx,ij_tdcomp, ij_dleaf
      use ghycom, only : snoage, snoage_def
      use veg_com, only : almass,aalbveg       !nyk
      use vegetation, only: crops_yr,cond_scheme !nyk
      use surf_albedo, only: albvnh  !nyk
      use sle001, only : fb,fv,ws
      use veg_drv, only : veg_set_cell

      implicit none
      real*8 tsavg,wfc1
      real*8 aleafmass, aalbveg0, fvp, sfv  !nyk veg ! , aleafmasslast
      integer i,j,itype
      integer northsouth,iv  !nyk
      logical, intent(in) :: end_of_day

C**** define local grid
      integer J_0, J_1

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT=J_0, J_STOP=J_1)

C**** Update vegetation file if necessary  (i.e. if crops_yr=0)
      if(crops_yr.eq.0) call updveg(jyear,.true.)
      if(cond_scheme.eq.2) call updsur (0,jday)
c****
c**** find leaf-area index & water field capacity for ground layer 1
c****
      cosday=cos(twopi/edpery*jday)
      sinday=sin(twopi/edpery*jday)
      do j=J_0,J_1
        if(j.le.jm/2) then      !nyk added northsouth
          northsouth=1.d0       !southern hemisphere
        else
          northsouth=2.d0       !northern hemisphere
        end if
        do i=1,im
          wfcs(i,j)=24.
          if (fearth(i,j).gt.0.) then
            !-----------------------------------------------------------
            !nyk Update vegetation albedos.
            !WARNING:  ALBVNH is not saved in a restart file.
            !          Must run without restarting.
            !albvnh(9,6,2)=albvnh(1+8veg,6bands,2hemi), band 1 is VIS.
            !if RUNDECK selects new conductance scheme
            if (cond_scheme.eq.2) then
              aalbveg0 = 0.d0
              sfv=0.d0
              do iv=1,8
                fvp=vdata(i,j,iv+1)
                sfv=sfv+fvp
                aalbveg0 = aalbveg0 + fvp*(ALBVNH(iv+1,1,northsouth))
                !write (99,*) 'fvp',fvp
                !write (99,*) 'ALBVNH',ALBVNH(iv+1,1,northsouth)
              end do
              aalbveg(i,j) = 0.08D0
              if(sfv.gt.0.) aalbveg(i,j) = aalbveg0/sfv !nyk
             !write (99,*) 'daily aalbveg', aalbveg(i,j)
            end if
            !-----------------------------------------------------------

            call ghinij(i,j)
            call veg_set_cell(i,j,.true.)
            wfc1=fb*ws(1,1)+fv*(ws(0,2)+ws(1,2))
            wfcs(i,j)=rhow*wfc1 ! canopy part changes

            !-----------------------------------------------------------
            !nyk - TEMPORARY calculate change in leaf mass per day
            !get aleafmass(i,j) at jday
            aleafmass=
     $           almass(1,i,j)+cosday*almass(2,i,j)+sinday*almass(3,i,j)

            !Calculate dlmass(i,j) increment from last jday
            !cosdaym1=cos(twopi/edpery*(jday-1))
            !sindaym1=sin(twopi/edpery*(jday-1))
            !aleafmasslast=almass(1,i,j)+cosdaym1*almass(2,i,j)+
!     $      !     sindaym1*almass(3,i,j)
            !accumulate dlmass
            !adlmass = aleafmass - aleafmasslast
            adlmass = aleafmass
            !aij(i,j,ij_dleaf)=aij(i,j,ij_dleaf)+adlmass
            aij(i,j,ij_dleaf)=adlmass  !accumulate just instant. value
            !PRINT '(F4.4)',adlmass                            !DEBUG
            !call stop_model('Just did adlmass',255)           !DEBUG
          end if
        end do
      end do

      if (end_of_day) then
        do j=J_0,J_1
        do i=1,imaxj(j)
c****
c**** increase snow age depending on snoage_def
c****
          if (snoage_def.eq.0) then ! update indep. of ts
            do itype=1,3
              snoage(itype,i,j)=1.+.98d0*snoage(itype,i,j)
            end do
          elseif (snoage_def.eq.1) then ! update if max T>0
            if (tdiurn(i,j,7).gt.0) snoage(1,i,j)=1.+.98d0
     *           *snoage(1,i,j) ! ocean ice (not currently used)
            if (tdiurn(i,j,8).gt.0) snoage(2,i,j)=1.+.98d0
     *           *snoage(2,i,j) ! land ice
            if (tdiurn(i,j,2).gt.0) snoage(3,i,j)=1.+.98d0
     *           *snoage(3,i,j) ! land
          else
            write(6,*) "This snoage_def is not defined: ",snoage_def
            write(6,*) "Please use: 0 (update indep of T)"
            write(6,*) "            1 (update if T>0)"
            call stop_model('stopped in GHY_DRV.f',255)
          end if
          tsavg=tdiurn(i,j,5)/(nday*nisurf)
          if(32.+1.8*tsavg.lt.65.)
     *         aij(i,j,ij_strngts)=aij(i,j,ij_strngts)+(33.-1.8*tsavg)
          aij(i,j,ij_dtgdts)=aij(i,j,ij_dtgdts)+18.*((tdiurn(i,j,2)-
     *         tdiurn(i,j,1))/(tdiurn(i,j,4)-tdiurn(i,j,3)+1.d-20)-1.)
          aij(i,j,ij_tdsl)=aij(i,j,ij_tdsl)+
     *         (tdiurn(i,j,4)-tdiurn(i,j,3))
          aij(i,j,ij_tdcomp)=aij(i,j,ij_tdcomp)+
     *         (tdiurn(i,j,6)-tdiurn(i,j,9))
          aij(i,j,ij_tmaxe)=aij(i,j,ij_tmaxe)+(tdiurn(i,j,4)-tf)
          if (tdiurn(i,j,6).lt.aij(i,j,ij_tmnmx))
     *         aij(i,j,ij_tmnmx)=tdiurn(i,j,6)
        end do
        end do
      end if


      return
      end subroutine daily_earth


      subroutine ground_e 1,7
!@sum  ground_e driver for applying surface fluxes to land fraction
!@auth original development team
!@ver  1.0
      use model_com, only : fearth,itearth
      use geom, only : imaxj,dxyp
      use ghycom, only : snowe, tearth,wearth,aiearth,wbare,wvege,snowbv
     *     ,fr_snow_ij,fr_snow_rad_ij, gdeep
      use veg_com, only : afb
      use dagcom, only : aj,areg,aij,jreg,ij_evap,ij_f0e,ij_evape
     *     ,ij_gwtr,ij_tg1,j_wtr1,j_ace1,j_wtr2,j_ace2
     *     ,j_snow,j_evap,j_type,ij_g01,ij_g07,ij_g28
     *     ,ij_g29,j_rsnow,ij_rsnw,ij_rsit,ij_snow
      use fluxes, only : e0,e1,evapor,eprec
      implicit none

      real*8 snow,tg1,tg2,f0dt,f1dt,evap,wtr1,wtr2,ace1,ace2
     *     ,pearth,enrgp,scove
      integer i,j,jr,k

C**** define local grid
      integer J_0, J_1

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT=J_0, J_STOP=J_1)

      do j=J_0,J_1
      do i=1,imaxj(j)
      pearth=fearth(i,j)
      jr=jreg(i,j)
      if (pearth.gt.0) then

        snow=snowe(i,j)
        tg1 = tearth(i,j)
        wtr1= wearth(i,j)
        ace1=aiearth(i,j)
        tg2=gdeep(i,j,1)
        wtr2=gdeep(i,j,2)
        ace2=gdeep(i,j,3)
        f0dt=e0(i,j,4)
        f1dt=e1(i,j,4)
        evap=evapor(i,j,4)
        enrgp=eprec(i,j)      ! including latent heat

c**** accumulate diagnostics
c**** the following is the actual snow cover of the snow model
c        scove = pearth *
c     *       ( afb(i,j)*fr_snow_ij(1,i,j)
c     *       + (1.-afb(i,j))*fr_snow_ij(2,i,j) )
c**** the following computes the snow cover as it is used in RAD_DRV.f
        scove = pearth *
     *       ( afb(i,j)*fr_snow_rad_ij(1,i,j)
     *       + (1.-afb(i,j))*fr_snow_rad_ij(2,i,j) )

        !if (snowe(i,j).gt.0.) scove=pearth
        aj(j,j_rsnow,itearth)=aj(j,j_rsnow,itearth)+scove
        areg(jr,j_rsnow)=areg(jr,j_rsnow)+scove*dxyp(j)
        aij(i,j,ij_rsnw)=aij(i,j,ij_rsnw)+scove
        aij(i,j,ij_snow)=aij(i,j,ij_snow)+snow*pearth
        aij(i,j,ij_rsit)=aij(i,j,ij_rsit)+scove

        aj(j,j_wtr1,itearth)=aj(j,j_wtr1,itearth)+wtr1*pearth
        aj(j,j_ace1,itearth)=aj(j,j_ace1,itearth)+ace1*pearth
        aj(j,j_wtr2,itearth)=aj(j,j_wtr2,itearth)+wtr2*pearth
        aj(j,j_ace2,itearth)=aj(j,j_ace2,itearth)+ace2*pearth
        aj(j,j_snow,itearth)=aj(j,j_snow,itearth)+snow*pearth
        areg(jr,j_snow)=areg(jr,j_snow)+snow*pearth*dxyp(j)
        areg(jr,j_wtr1)=areg(jr,j_wtr1)+wtr1*pearth*dxyp(j)
        areg(jr,j_ace1)=areg(jr,j_ace1)+ace1*pearth*dxyp(j)
        areg(jr,j_wtr2)=areg(jr,j_wtr2)+wtr2*pearth*dxyp(j)
        areg(jr,j_ace2)=areg(jr,j_ace2)+ace2*pearth*dxyp(j)

        aij(i,j,ij_f0e)  =aij(i,j,ij_f0e)  +f0dt+enrgp
        aij(i,j,ij_gwtr) =aij(i,j,ij_gwtr)+(wtr1+ace1+wtr2+ace2)
        aij(i,j,ij_evape)=aij(i,j,ij_evape)+evap
        do k=1,4
          aij(i,j,ij_g01+k-1)=aij(i,j,ij_g01+k-1)+wbare(k,i,j)
          aij(i,j,ij_g07+k-1)=aij(i,j,ij_g07+k-1)+wvege(k-1,i,j)
        end do
        aij(i,j,ij_g28)=aij(i,j,ij_g28)+snowbv(1,i,j)
        aij(i,j,ij_g29)=aij(i,j,ij_g29)+snowbv(2,i,j)
      end if
c****
      end do
      end do
      end subroutine ground_e


      subroutine conserv_wtg(waterg),6
!@sum  conserv_wtg calculates zonal ground water incl snow
!@auth Gavin Schmidt
!@ver  1.0
      use constant, only : rhow
      use model_com, only : fim,fearth
      use geom, only : imaxj
      use ghycom, only : ngm,wbare,wvege,snowbv
      use veg_com, only : afb
      implicit none
!@var waterg zonal ground water (kg/m^2)
      real*8, dimension(grid%j_strt:grid%j_stop) :: waterg
      integer i,j,n
      real*8 wij,fb

C**** define local grid
      integer :: J_0, J_1
      logical :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT=J_0, J_STOP=J_1,
     &               HAVE_SOUTH_POLE = HAVE_SOUTH_POLE,
     &               HAVE_NORTH_POLE = HAVE_NORTH_POLE)

      do j=J_0,J_1
        waterg(j)=0
        do i=1,imaxj(j)
          if (fearth(i,j).gt.0) then
            fb=afb(i,j)
            wij=fb*snowbv(1,i,j)+(1.-fb)*(wvege(0,i,j)+snowbv(2,i,j))
            do n=1,ngm
              wij=wij+fb*wbare(n,i,j)+(1.-fb)*wvege(n,i,j)
            end do
            waterg(j)=waterg(j)+fearth(i,j)*wij*rhow
          end if
        end do
      end do
      if (HAVE_SOUTH_POLE) waterg(1) =fim*waterg(1)
      if (HAVE_NORTH_POLE) waterg(jm)=fim*waterg(jm)
c****
      end subroutine conserv_wtg


      subroutine conserv_htg(heatg),5
!@sum  conserv_htg calculates zonal ground energy incl. snow energy
!@auth Gavin Schmidt
!@ver  1.0
      use model_com, only : fim,fearth
      use geom, only : imaxj, dxyp
      use ghycom, only : ngm,htbare,htvege,fr_snow_ij,nsn_ij,hsn_ij
      use veg_com, only : afb
      implicit none
!@var heatg zonal ground heat (J/m^2)
      real*8, dimension(grid%j_strt:grid%j_stop) :: heatg
      integer i,j
      real*8 hij,fb,fv

C**** define local grid
      integer J_0, J_1
      logical :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT=J_0, J_STOP=J_1,
     &               HAVE_SOUTH_POLE = HAVE_SOUTH_POLE,
     &               HAVE_NORTH_POLE = HAVE_NORTH_POLE)

      do j=J_0,J_1
        heatg(j)=0
        do i=1,imaxj(j)
          if (fearth(i,j).le.0) cycle
          fb=afb(i,j)
          fv=(1.d0-fb)
          hij=fb*sum( htbare(1:ngm,i,j) )
     &       +  fv*sum( htvege(0:ngm,i,j) )
     &       +  fb*fr_snow_ij(1,i,j)*sum( hsn_ij(1:nsn_ij(1,i,j),1,i,j))
     &       +  fv*fr_snow_ij(2,i,j)*sum( hsn_ij(1:nsn_ij(2,i,j),2,i,j))
          heatg(j)=heatg(j)+fearth(i,j)*hij
        end do
      end do
      if (HAVE_SOUTH_POLE) heatg(1) =fim*heatg(1)
      if (HAVE_NORTH_POLE) heatg(jm)=fim*heatg(jm)
c****
ccc debugging ...
ccc      print *,'conserv_htg energy ',
ccc     &     sum(heatg(1:jm)*dxyp(1:jm))/(sum(dxyp(1:jm))*im)
      end subroutine conserv_htg


      end module soil_drv


      subroutine check_ghy_conservation( flag ),10
ccc debugging program: cam be put at the beginning and at the end
ccc of the 'surface' to check water conservation
      use constant, only : rhow
      use geom, only : imaxj
      use model_com, only : im,jm,fearth
      use DOMAIN_DECOMP, only : GRID, GET
      use fluxes, only : prec,evapor,runoe
      use ghycom, only : ngm,wbare,wvege,htbare,htvege,snowbv,dz_ij
      use veg_com, only : afb
      implicit none
      integer flag
      real*8 total_water(im,jm), error_water
      real*8, save :: old_total_water(im,jm)
!      real*8 total_energy(im,jm), error_energy
!      real*8, save :: old_total_energy(im,jm)
      integer i,j,n
      real*8 fb,fv
ccc enrgy check not implemented yet ...

C**** define local grid
      integer J_0, J_1

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT=J_0, J_STOP=J_1)

      do j=J_0,J_1
        do i=1,imaxj(j)
          if ( fearth(i,j) <= 0.d0 ) cycle

ccc just checking ...
          do n = 1,ngm
            if ( dz_ij(i,j,n) .le. 0.d0 )
     &           call stop_model('incompatible dz',255)
          enddo

          fb = afb(i,j)
          fv = 1.d0 - fb
          total_water(i,j) = fb*sum( wbare(1:ngm,i,j) )
     &         + fv*sum( wvege(0:ngm,i,j) )
     &         + fb*snowbv(1,i,j) + fv*snowbv(2,i,j)
        end do
      end do

      ! call stop_model('just testing...',255)

      if ( flag == 0 ) then
        old_total_water(:,:) = total_water(:,:)
        return
      endif

      do j=J_0,J_1
        do i=1,imaxj(j)

          !print *,'fearth = ', i, j, fearth(i,j)

          if ( fearth(i,j) <= 0.d0 ) cycle
          fb = afb(i,j)
          fv = 1.d0 - fb
          error_water = ( total_water(i,j) - old_total_water(i,j) )*rhow
     &         - prec(i,j) + evapor(i,j,4) + runoe(i,j)

          !print *, 'err H2O: ', i, j, error_water

  !        if ( abs( error_water ) > 1.d-9 ) print *, 'error'
          if ( abs( error_water ) > 1.d-9 ) call stop_model(  ! was -15
     &         'check_ghy_conservation: water conservation problem',255)

        end do
      end do

      end subroutine check_ghy_conservation