SUBROUTINE surfce(ape,zint,cklq)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  IMPLICIT NONE
  INCLUDE "parmeta.f90"
  INTEGER, PARAMETER :: lp1=lm+1
  REAL, INTENT(in)   , DIMENSION(  im  ,   jm  , lm ) :: ape
  REAL, INTENT(in)   , DIMENSION(0:im+1, 0:jm+1, lp1) :: zint  
  REAL, INTENT(inout), DIMENSION(  im  ,   jm       ) :: cklq
  !
  REAL, PARAMETER :: epswet=.001e0
  REAL, PARAMETER :: pq0=379.90516e0
  REAL, PARAMETER :: seafc=.98e0
  REAL, PARAMETER :: tresh=.95e0
  REAL, PARAMETER :: a2=17.2693882e0
  REAL, PARAMETER :: a3=273.16e0
  REAL, PARAMETER :: a4=35.86e0
  REAL, PARAMETER :: t0=273.16e0
  REAL, PARAMETER :: t1=274.16e0
  REAL, PARAMETER :: capa=0.28589641
  REAL, PARAMETER :: cp=1004.6e0
  REAL, PARAMETER :: stbol=5.67e-8
  REAL, PARAMETER :: r=287.04e0
  REAL, PARAMETER :: row=1.e3
  REAL, PARAMETER :: elwv=2.50e6
  REAL, PARAMETER :: eliv=2.834e6
  REAL, PARAMETER :: eliw=.334e6
  REAL, PARAMETER :: a23m4=a2*(a3-a4)
  REAL, PARAMETER :: pq0sea=pq0*seafc
  REAL, PARAMETER :: pq0c=pq0*tresh
  REAL, PARAMETER :: rlivwv=eliv/elwv
  REAL, PARAMETER :: rowliw=row*eliw
  REAL, PARAMETER :: rowliv=row*eliv
  INCLUDE "parm.tbl.f90"
  INTEGER, PARAMETER :: la=13
  INTEGER, PARAMETER :: imjm=im*jm-jm/2
  INTEGER, PARAMETER :: jam=6+2*(jm-10)
  INTEGER, PARAMETER :: lscrch=4*lm+1+la +1
  INTEGER, PARAMETER :: l2=la+lm+1
  INTEGER, PARAMETER :: l4=la+3*lm+1
  REAL, DIMENSION(im,jm) :: zlm
  REAL, DIMENSION(im,jm) :: ps
  REAL, DIMENSION(im,jm) :: apes
  REAL, DIMENSION(im,jm) :: etalm
  REAL, DIMENSION(im,jm) :: plm
  REAL, DIMENSION(im,jm) :: apelm
  REAL, DIMENSION(im,jm) :: rdsin
  REAL, DIMENSION(im,jm) :: tlm
  REAL, DIMENSION(im,jm) :: thlm
  REAL, DIMENSION(im,jm) :: qlm
  REAL, DIMENSION(im,jm) :: qlms
  REAL, DIMENSION(im,jm) :: dqsdt
  REAL, DIMENSION(im,jm) :: ffs
  REAL, DIMENSION(im,jm) :: qfc1
  !
  !  additional soil model constants and variables
  !
  INTEGER, PARAMETER :: nsoil=2
  INTEGER, PARAMETER :: nroot=2
  REAL, DIMENSION(nsoil) ::  smck
  REAL, DIMENSION(nsoil) ::  stck
  !  include common blocks.
  !     common blocks soil, acmpre, acmsfc were added for gcip.
  !     common block options was added for the post.
  INCLUDE "ctlblk.f90"
  INCLUDE "loops.f90"
  INCLUDE "masks.f90"
  INCLUDE "phys.f90"
  INCLUDE "vrbls.f90"
  INCLUDE "pvrbls.f90"
  INCLUDE "soil.f90"
  INCLUDE "acmpre.f90"
  INCLUDE "acmsfc.f90"
  INCLUDE "acmrds.f90"
  INCLUDE "acmrdl.f90"
  !  include "options"

  REAL :: fdtliw, fdtliv, fdtw, ffsk, factrs, factrl, tlmh
  REAL :: satflg, dtk, z, fk, soldn, sfcprs, prcp, q2k, q2sat
  REAL :: dqsdtk, tbot, chk, chkff, vgfrck, q1k, sfctmp, sfcth2
  REAL :: t1k, cmck, snodpk, plflx, elflx, hflx, gflx
  REAL :: rnof1k, rnof2k, smeltk, soilqw, soilqm, scheck, tsfc, tsfc2
  INTEGER :: i, j, ns, llmh, ice, ivgtpk, isltpk
  !
  !                      start surfce here
  !
  !  set constants calculated here for clarity.
  !
  fdtliw=dtq2/rowliw
  fdtliv=dtq2/rowliv
  fdtw=dtq2/2.5e9
  !
  !  set osu model constants and time independent variables
  !  initialize osu model historical variables
  !

  IF(ntsd.LT.nphs)THEN
     DO j=1,jm
        DO i=1,im
           ps(i,j)=pd(i,j)+pt
           apes(i,j)=(1.e5/ps(i,j))**capa
           pctsno(i,j)=-999.0

           ! set zero-value for some output diagnostic arrays
!chou
           IF(sm(i,j).LT.0.5)THEN
              IF(sice(i,j).EQ.1.) THEN
                 !        sea-ice case
                 smstav(i,j)=1.0
                 smstot(i,j)=1.0
                 ssroff(i,j)=0.0
                 bgroff(i,j)=0.0
                 cmc(i,j)=0.0
                 DO ns=1,nsoil
                    smc(i,j,ns)=1.0
                 ENDDO
              ENDIF
           ELSE
              !       water case
              smstav(i,j)=1.0
              smstot(i,j)=1.0
              ssroff(i,j)=0.0
              bgroff(i,j)=0.0
              soiltb(i,j)=273.16
              grnflx(i,j)=0.
              subshx(i,j)=0.0
              acsnow(i,j)=0.0
              acsnom(i,j)=0.0
              snopcx(i,j)=0.0
              cmc(i,j)=0.0
              sno(i,j)=0.0
              DO ns=1,nsoil
                 smc(i,j,ns)=1.0
                 stc(i,j,ns)=273.16
              ENDDO
           ENDIF

        END DO
     END DO
  ENDIF

  !
  !  set lowest model layer variables.
  !
  DO  j=3,jm-2
     DO  i=1,im
        llmh=lmh(i,j)
        etalm(i,j)=aeta(llmh)
        apelm(i,j)=ape(i,j,llmh)
        tlm(i,j)=t(i,j,llmh)
        qlm(i,j)=q(i,j,llmh)
        zlm(i,j)=(zint(i,j,llmh)-zint(i,j,llmh+1))*0.50
     END DO
  END DO

  DO j=3,jm-2
     DO i=1,im
        ps(i,j)=pd(i,j)+pt
        apes(i,j)=(1.e5/ps(i,j))**capa
        plm(i,j)=etalm(i,j)*pd(i,j)*res(i,j)+pt
        qlms(i,j)=((1.-sm(i,j))*pq0+sm(i,j)*pq0sea) &
             /plm(i,j)*EXP(a2*(tlm(i,j)-a3)/(tlm(i,j)-a4))
        dqsdt(i,j)=qlms(i,j)*a23m4/(tlm(i,j)-a4)**2
        ffsk=akhs(i,j)*plm(i,j)*hbm2(i,j)/((qlm(i,j)*.608+1.)*tlm(i,j)*r)
        qfc1(i,j)=apes(i,j)*ffsk*elwv
        ffs(i,j)=ffsk*cp
     END DO
  END DO

  DO j=3,jm-2
     DO i=1,im
        !
        !  compute radin and rdsin for this timestep
        !  czen is in phys common and is current from call to rdtemp
        !
        IF(czmean(i,j).GT.0.)THEN
           factrs=czen(i,j)/czmean(i,j)
        ELSE
           factrs=0.
        ENDIF

        IF(sigt4(i,j).GT.0.)THEN
           tlmh=tlm(i,j)
           factrl=stbol*tlmh*tlmh*tlmh*tlmh/sigt4(i,j)
        ELSE
           factrl=0.
        ENDIF

        radin(i,j)=((rswin(i,j)-rswout(i,j))*factrs+ &
             rlwin(i,j)*factrl)*hbm2(i,j)
        rdsin(i,j)= rswin(i,j)*hbm2(i,j)*factrs
        !
        !  diagnostic radiation accumulation
        !
        aswin (i,j)=aswin (i,j)+rswin (i,j)*hbm2(i,j)*factrs
        aswout(i,j)=aswout(i,j)-rswout(i,j)*hbm2(i,j)*factrs
        aswtoa(i,j)=aswtoa(i,j)+rswtoa(i,j)*hbm2(i,j)*factrs
        alwin (i,j)=alwin (i,j)+rlwin (i,j)*hbm2(i,j)*factrl
        alwout(i,j)=alwout(i,j)-radot (i,j)*hbm2(i,j)
        alwtoa(i,j)=alwtoa(i,j)+rlwtoa(i,j)*hbm2(i,j)
        !
        !  check for saturation at the lowest model level
        !
        IF((qlm(i,j).GE.qlms(i,j)*tresh).AND.(qlm(i,j).LT.qz0(i,j)))THEN
           cklq(i,j)=0.
        ELSE
           cklq(i,j)=hbm2(i,j)
        ENDIF
     END DO
  END DO
  !
  !  ths, thlm, cheating wet for profs
  !
  DO  j=3,jm-2
     DO  i=1,im
        thlm(i,j)=tlm(i,j)*apelm(i,j)
        qfc1(i,j)=qfc1(i,j)*cklq(i,j)
     END DO
  END DO
  !
  DO j=3,jm-2
     DO i=1,im
        IF(hbm2(i,j).GE.0.5)THEN
!chou
           IF(sm(i,j).GT.0.5)THEN
              ths(i,j)=sst(i,j)*apes(i,j)
              qs(i,j)=hbm2(i,j)*pq0sea/ps(i,j) &
                   *EXP(a2*(ths(i,j)-a3*apes(i,j))/(ths(i,j)-a4*apes(i,j)))
           ENDIF
           !
           !  loading and unloading nmc/osu land soil variables
           !
!chou
           IF(sm(i,j).LT.0.5)THEN
              ice=INT(sice(i,j)+0.3)
              satflg=cklq(i,j)
              dtk=dtq2
              z=zlm(i,j)
              fk=radin(i,j)
              soldn=rdsin(i,j)
              sfcprs=plm(i,j)
              prcp=prec(i,j)*row/dtq2
              q2k=qlm(i,j)
              q2sat=qlms(i,j)
              dqsdtk=dqsdt(i,j)
              tbot=tg(i,j)
              chk=akhs(i,j)
              chkff=ffs(i,j)
              ivgtpk=ivgtyp(i,j)
              isltpk=isltyp(i,j)
              !  meb  prevent routines in sflx from going out of bounds
              IF (ivgtpk.EQ.0) ivgtpk=13
              IF (isltpk.EQ.0) isltpk=9
              !  meb  prevent routines in sflx from going out of bounds
              vgfrck=vegfrc(i,j)
              q1k=qs(i,j)
              sfctmp=thlm(i,j)/apelm(i,j)
              sfcth2=thlm(i,j)/apes(i,j)
              t1k=ths(i,j)/apes(i,j)
              cmck=cmc(i,j)
              snodpk=sno(i,j)
              !
              DO ns=1,nsoil
                 smck(ns)=smc(i,j,ns)
                 stck(ns)=stc(i,j,ns)
              END DO
              CALL sflx (ice   ,satflg,dtk   ,z, &
                   fk    ,soldn ,sfcprs,prcp  ,sfctmp,sfcth2, &
                   q2k   ,q2sat ,dqsdtk,tbot  ,chk,   chkff, &
                   ivgtpk,isltpk,vgfrck, &
                   plflx ,elflx ,hflx  ,gflx  ,rnof1k,rnof2k, &
                   q1k   ,smeltk,t1k   ,cmck  ,smck  ,stck  ,snodpk, &
                   soilqw,soilqm )
              scheck=z*chk
              IF(scheck.LE.1.3e-3)THEN
                 plflx=0.
                 elflx=0.
              ENDIF
              !
              !  gcip diagnostics & modification of qfc1 over snow
              !
              ssroff(i,j)=ssroff(i,j)+rnof1k*dtq2
              bgroff(i,j)=bgroff(i,j)+rnof2k*dtq2
              smstav(i,j)=soilqw
              soiltb(i,j)=tbot
              sfcexc(i,j)=chk
              grnflx(i,j)=gflx
              IF(sno (i,j).GT.0..OR.sice(i,j).GT.0.)THEN
                 qfc1(i,j)=qfc1(i,j)*rlivwv
              ENDIF
              IF(sno(i,j).GT.0.)THEN
                 acsnom(i,j)=acsnom(i,j)+smeltk
                 snopcx(i,j)=snopcx(i,j)-smeltk/fdtliw
              ENDIF
              potevp(i,j)=potevp(i,j)+plflx*fdtw
              potflx(i,j)=potflx(i,j)-plflx
              subshx(i,j)=subshx(i,j)+gflx
              !
              !  eta model lower boundary conditions
              !
              ths(i,j)=t1k*apes(i,j)
              IF(qfc1(i,j).GT.0.)qs(i,j)=qlm(i,j)+elflx*apes(i,j)/qfc1(i,j)
              !
              !  historical variables
              !
              sno(i,j)=snodpk
              cmc(i,j)=cmck
              !
              smstot(i,j)=soilqm
              DO  ns=1,nsoil
                 smc(i,j,ns)=smck(ns)
                 stc(i,j,ns)=stck(ns)
              END DO
           ENDIF
           !
        END IF
     END DO
  END DO

  !  variables twbs and qwbs computed here for gcip.
  !  accumulate surface heat fluxes here.
  !  for gcip accumulate actual and potential evaporation.
  !  for gcip accumulate total snow melt and
  !  the associated net heat flux.

  DO j=3,jm-2
     DO i=1,im
        twbs(i,j)=(thlm(i,j)-ths(i,j)*(1.-sm(i,j))-thz0(i,j)*sm(i,j)) &
             *ffs (i,j)/apes(i,j)
        qwbs(i,j)=(qlm (i,j)-qs (i,j)*(1.-sm(i,j))-qz0 (i,j)*sm(i,j)) &
             *qfc1(i,j)/apes(i,j)
        sfcshx(i,j)=sfcshx(i,j)+twbs(i,j)
        sfclhx(i,j)=sfclhx(i,j)+qwbs(i,j)
        sfcevp(i,j)=sfcevp(i,j)-qwbs(i,j)*fdtw
        potevp(i,j)=potevp(i,j)-qwbs(i,j)*sm(i,j)*fdtw
        potflx(i,j)=potflx(i,j)+qwbs(i,j)*sm(i,j)

        ! if cold enough, it snows (in osu model)...
        ! for gcip accumulate total snowfall.

!chou
        IF(thlm(i,j)/apelm(i,j).LE.t0.AND.sice(i,j)+sm(i,j).LT.0.5)THEN
           acsnow(i,j)=acsnow(i,j)+prec(i,j)

           !  ... otherwise it rains.

        ELSE
           accliq(i,j)=accliq(i,j)+prec(i,j)
        ENDIF

        prec(i,j)=0.
     END DO
  END DO

  !  longwave outgoing radiation

  DO j=3,jm-2
     DO i=1,im
        tsfc=ths(i,j)/apes(i,j)
        tsfc2=tsfc*tsfc
        radot(i,j)=hbm2(i,j)*epsr(i,j)*stbol*tsfc2*tsfc2
     END DO
  END DO

  !   increment time step counters for use in computing time ave values

  aphtim = aphtim + 1.
  ardsw  = ardsw  + 1.
  ardlw  = ardlw  + 1.
  asrfc  = asrfc  + 1.
END SUBROUTINE surfce
