c  @(#)rdflux.f	1.3


      subroutine rdflux

      include'mladv.inc'

      parameter(c=1005.0,rho=1.240,ch=0.83e-3,conv=2.3901e-5)
      parameter(ce=1.5e-3,vl0=2.45e6,qsat0=0.013)
      parameter(bbrc=0.985,sigma=5.67e-8,aback=0.39,bback=0.05)
      parameter(cc0=0.5,rv=461.5)

      common/sbldef/qsw0,qsurf0,wflux0,u10m,angu10m,hdpth0
      common /unit/iin,iout,idff,idif,iorr,iowr,iplt,ifndd
      common/parm/kz,m,mm1,mm2,mp1,msq,tim,tmax,tstart,hy,rea,dt,
     +   xbasin,yy(mx0),icalc,kzp1,kzm1,iftop,ifbot,iftvv,ifrel,theta,
     +   ifsbl,ifeva

      common/inform/ifdiff,ifpert,titl(20,2),runid(6),restid(9),
     +              idstring,ident
c
      character*70 idstring
      character*80 ident,fname
      character*4 restid,runid

      data icall/0/

1     continue

      if (iftvf .eq. 1) then

          read(ifunit,*)tflx,qsw0,qsurf0,u10m,angu10m

         else if (iftvf .eq. 2) then

          read(ifunit,*)tflx,u10m,angu10m,airtem,relhum,qsw0

         else

          read(ifunit,*)tflx,u10m,angu10m,airtem,relhum,clfac,qsw0

        endif

       if (iftvf .gt. 1) then

          sstave = 0.0
          do 2 ip=1,msq
             sstave = sstave + temp(ip,1,mlcur) + tqgofk(ip,1)
  2        continue

          sstave = 273. + sstave / float(msq)

          if (iftvf .eq. 2) then
             qsat = qsat0
             vl = vl0
             cc = cc0
             qa = (relhum/100.) * qsat
            else
             qsat = (.622 * .6112 / 100.) * exp (17.67 * 
     .              (airtem - 273.16) / ( airtem - 29.66) )
             qa = (relhum/100.) * qsat
             qsat = (.622 * .6112 / 100.) * exp (17.67 * 
     .              (sstave - 273.16) / (sstave - 29.66) )
             vl = (2.501 - 0.00237 * (airtem - 273.)) * 1.0e6
             cc = clfac
           endif

          rhov = qa * rho
          ea = rhov * rv * airtem / 100.0

          hsens = c * rho * ch * u10m * (airtem - sstave) * conv
          hlat = - rho * ce * u10m * (qsat - qa) * vl * conv
          hback =  - bbrc * sigma * sstave**4
     .                       * (aback - bback * ea **.5)
     .                       * (1.0 - 0.6 * cc**2) * conv
c         hback = 0.96 * sigma * airtem**4
c    .              * (0.254 - 0.00495 * ea) * (1.0 - .63 * cc)
c    .              + 4.0 * 0.96 * sigma * airtem**3
c    .              * (sstave - airtem)
c         hback = -hback * conv
          qsurf0 = hsens + hlat + hback

        endif
c
      qsw0 = qsw0 * qswfac + qswadd
      qsurf0 = qsurf0 * qsurfac + qsuradd
      u10m = u10m * u10mfac + u10madd
      angu10m = angu10m * angfac + angadd

        if (icall.eq.0) then
           icall = 1
           iflen=lnblk(ident,80)
           fname=ident(:iflen)//'.sfc'
           open(unit=iusfc,file=fname,form='formatted',status='new')
         endif

      write(iusfc,4)tflx,u10m,angu10m,qsw0,qsurf0,hsens,hlat,hback
  4   format(10e12.4)

      write(iout,*)'rdflux',tflx,mltime
      if (abs(tflx-mltime)*86400..gt.(mlstep/2.)) then
         write(iout,*)'Surface forcing - model time mismatch'
         write(iout,*)tflx,mltime,mlstep
         goto 1
        endif

      return
      end
