c    @(#)readsbl.f	1.9
      subroutine readsbl
c
c -- read SBL control parameters and simple atmospheric fluxes
c      
      character*80 flxfil
      include'mladv.inc'
c ---------------------------------------------------------------------
c -- qg common
      common/user/pram(16)
      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 /threed/hz(k0),scisg,sigz(k01),eigval(k0),ainv(k0,k0),
     .               amat(k0,k0)
      common /extpar/ rlat0, rlng0, t0, v0, dhor, ht, time0,
     +                r0,f0star
      common /cnn/n,nm1,nm2,np1
      common /unit/iin,iout,idff,idif,iorr,iowr,iplt,ifndd

c
c---------------------------------------------------------------------
      common/sbldef/qsw0,qsurf0,wflux0,u10m,angu10m,hdpth0

c---------------------------------------------------------------------

c
c -- defaults
c
c
c -- atmospheric fluxes
c
      iftvf = 0
      qsw0 = 0.0
      qsurf0 = 0.0
      wflux0 = 0.0
      u10m = pram(4)
      angu10m = pram(5)
      qswfac  = 1.0
      qsurfac = 1.0
      u10mfac = 1.0
      angfac  = 1.0
      momflx = 1
      ifwbg = 0
c
c -- number of output fields
c
      mout = 0
c
c -- number of tracers
c
      ntrace = 0
c
c -- tracer physics turned off
c
      itphys = 0
c
c -- tracer flux computations off
c
      dtrflx = 0.
c
c -- tracer budgets turned off
c
      ibudf = 0
      ibudi = 0
c
c -- positive definite tracer flags and sinking rate
c
      do 69 i=1,maxtrace
         ipdt(i) = ipdtn
         wsink(i) = 0.0
  69   continue
c
c -- no diffusion
c
      dc = 0.0
c
c -- vertical advections turned on
c
      ifvadv = 1
c
c -- horizontal advections turned on
c
      ifadvqg = 1
      ifadvml = 1
c
c -- boundary conditions not persisted
c
      ifpers = 0
c
c -- mixing layer depth not persisted
c
      ifperh = 0
c
c -- mixed layer filtering
c
      nord = 2
      nfrq = 1
      ncnt = 1
c
c -- number of time steps to skip in computing w
c
      iwskip = 0
c
c --  sbl diagnostics
c
      ifsbld = 0
c
c -- artificial top density perturbation
c
      artop = 0.0
c
c -- use only top qg level
c
      ifmqg = 0
c
c -- do not read background profile
c
      ifrbg = 0
c
c -- minimum and maximum excursion of mixed layer
c -- if zero use the defaults set at the end of this subroutine
c
      hmin0  = 0
      hmax0  = 0
c
c --- ag is the thermal expansion coefficent *g
c --- used in heat.F and mlbflux.F
c
      ag=9.8*.000246
c
c --- temperature threshhold used in mix.F
c
      tthresh=.1
c
c --- set default decay is the extinction coefficent for sw radiation.
c --- used in heat.F
c
      decay=0.08333333

c
c
c --- set up constants needed for the garwood entrainment model.
c
      p1 = 1.0
      p2 = 1.0
      m3 = 4.5
      azr = 1.0

      f0=f0star
      m1=1.0
      m2=p2
      m4=p1

      m5=4.6
      p3=m5/m1
c
c
      idone = 0
      ilevelsin = 0
 1    continue
      read(iin,*)icard
      if(icard.lt.0) then
c
c -- exit flag
c
         idone = 1
      elseif(icard.lt.100 .or. icard .gt. 199) then
c
c -- invalid card number stop
c
         write(iout,"(' Error in SBL input definition')")
         write(iout,"(' Card ',i5,' is not supported ')")
         write(iout,"(' fatal error in RDSBL')")
         stop
      elseif(icard.eq.100) then
c
c -- define vertical grid
c         
         read(iin,*)l,amp,dzml(1)
         if(l.lt.2 .or. l.gt.lmax) then
            write(iout,"(' Number of vertical levels out of range ')")
            write(iout,"(' requested:',i5,' max ',i5)")l,lmax
            write(iout,"(' error in RDSBL: No. of vertical levels')")
            stop
         endif
         zml(1)=-0.5*dzml(1)
c

         do 10 k=2,l+1
            if (amp.lt.0.0) then
               read(iin,*)dzml(k)
              else
               dzml(k)=amp*dzml(k-1)
             endif
            zml(k)=zml(k-1)-.5*(dzml(k-1)+dzml(k))
 10      continue
         ilevelsin = 1
c
      elseif(icard.eq.101) then
c
c   How many tracers? DEFAULT = 0
c
         read(iin,*)ntrace
         if (ntrace.gt.maxtrace) then
            print*,'Tracer arrays not dimensioned large enough'
            stop
           endif

      elseif(icard.eq.102) then
c
c
c --- now read in the field output data.
c
         read(iin,*)mout
         if (mout.ge.0) then
            nout=mout
            if (nout.gt.0) then
               do 25 ii=1,nout
                  read(iin,*)(dout(jj,ii),jj=1,3)
 25            continue
            endif
         endif
      elseif(icard .eq.103) then
c
c -- mixed layer initial depth, minimum and max depth

         read(iin,*)hdpth0

      elseif(icard.eq.104) then
c
c -- atmospheric input: heat fluxes
c
         read(iin,*)qsw0,qsurf0
      elseif(icard.eq.105) then
c
c -- atmospheric input: wind
c
         read(iin,*)u10m,angu10m
      elseif(icard.eq.106) then
c
c -- threshold T and g*alpha
c
         read(iin,*)tthresh,ag
      elseif(icard.eq.107) then
c
c -- reset min and max excursion of mixed layer depth
c
         read(iin,*)hmin0,hmax0
      elseif(icard.eq.108) then
c
c -- reset Garwood's constants
c
         read(iin,*)m3,azr
         m5=((m2-.3333333*m1)*(m3/3.0/m2)**.3333333)*3.0/azr

      elseif(icard.eq.109) then
c
c -- read sbl restart file name
c
         read(iin,*)fmlrrs

      else if (icard.eq.110) then
c
c -- Time varying surface fluxes
c
         read(iin,*)iftvf
         if (iftvf.gt.0) then
            read(iin,*)flxfil
            open(unit=ifunit,file=flxfil,status='old',form='formatted')
          endif
      elseif(icard.eq.111) then

c
c -- vertical advection switch
c
         read(iin,*) ifvadv
      else if (icard.eq.112) then
c
c -- mixed layer filtering
c
         read(iin,*) nord,nfrq,ncnt
      elseif(icard.eq.113) then
c
c -- sbl boundary conditions
c
         read(iin,*) ifpers
      elseif(icard.eq.114) then
c
c -- iwskip
c
         read(iin,*) iwskip
      elseif(icard.eq.115) then
c
c -- ifperh
c
         read(iin,*) ifperh
      elseif(icard.eq.116) then
c
c -- positive definite tracer flags
c
         read(iin,*) (ipdt(i),i=1,ntrace)
      elseif(icard.eq.117) then
c
c -- itphys
c
         read(iin,*) itphys
      elseif(icard.eq.118) then
c
c -- diffusion coefficient
c
         read(iin,*) dc
      elseif(icard.eq.119) then
c
c -- tracer sinking rates
c
         read(iin,*) (wsink(i),i=1,ntrace)
      elseif(icard.eq.120) then
c
c -- artificial top density perturbation
c
         read(iin,*) artop
      elseif(icard.eq.121) then
c
c -- surface flux multipliers
c
         read(iin,*) qswfac,qsurfac,u10mfac,angfac
      elseif(icard.eq.122) then
c
c -- horizontal advection switches
c
         read(iin,*) ifadvqg,ifadvml
      elseif(icard.eq.123) then
c
c -- sbl diagnostics
c
         read(iin,*) ifsbld,idl
      elseif(icard.eq.124) then
c
c -- multiple qg level option
c
         read(iin,*) ifmqg

         zqg(1) = ht * 0.5 * hz(1)
         do 900 k=2,kz
            zqg(k) = zqg(k-1) + ht * 0.5 * (hz(k) + hz(k-1))
 900      continue
      elseif(icard.eq.125) then
c
c -- read in background profile
c
         read(iin,*) ifrbg
         if (ifrbg .eq. 1) then
            read(iin,*)bgfile
          endif

      elseif(icard.eq.126) then
c
c -- momentum flux parameterization
c
         read(iin,*) momflx

      elseif(icard.eq.127) then
c
c -- background vertical velocity
c
         read(iin,*) ifwbg
         if (ifwbg .ge. 1) then
            read(iin,*)wbgbar,wbgamp,wbgpha
            if (ifwbg .eq. 2) then
               read(iin,*)wbgref
             endif
           endif
      elseif(icard.eq.128) then
c
c -- additive surface fluxes
c
         read(iin,*)qswadd,qsuradd,u10madd,angadd
      elseif(icard.eq.129) then
c
c -- tracer flux computation interval
c
         read(iin,*)dtrflx
      elseif(icard.eq.130) then
c        
c -- tracer budget computations
c
         read(iin,*)ibudf,ibudi
         if ((ibudf .ne. 0) .or. (ibudi .ne. 0)) then
            read(iin,*)dtbudf,dtbudi
           endif
         if (ibudi .ne. 0) then
            read(iin,*)istrip
           endif

      endif
c
c
c -- read another card until done
c
      if (idone .eq. 0) goto 1
c
      if (ilevelsin .eq. 0) then
         write(iout,231)
         write(iout,*)' SBL levels undefined'
         stop
      endif
c
c --- the default  minimum depth to which the mixed layer can 
c --- shallow is  one half of the first layer or hmin=.5*dz(1)
c --- and the maximum is the depth of the last layer.
c
      if(hmin0 .eq.0) then
         hmin=.5*dzml(1)
      else
         hmin = max(0.5*dzml(1),hmin0)
      endif
      if(hmax0 .eq.0) then
         maxh=-zml(l)
      else
         maxh = min(-zml(l),hmax0)
      endif
      if(hdpth0 .lt. hmin) then
         hdpth0 = hmin
      elseif(hdpth0.gt.maxh) then
         hdpth0 = maxh
      endif
c
c -- echo some parameters
c
      write(iout,101)ag,decay,qgts,p1,p2,p3,azr,m3,(-zml(k),k=1,l)
101   format(///,'  mixed layer model constants  ',///,
     .          '  ag:              ',e13.5,//,
     .          '  decay:           ',e13.5,' (1/m)',//,
     .          '  time step        ',f7.0,' (s)',//,
     .          '  p1:              ',f7.4,//,
     .          '  p2:              ',f7.4,//,
     .          '  p3, azr          ',f7.4,', ',f7.4,//,
     .          '  m3               ',f7.4,///,
     .          '  level depths: (m)',//,10('  ',f7.3),//,
     .          10('  ',f7.3))
 231  format(' SBL card 100: SBL levels must be entered')
 
c
c -- read in tracer information
c
      if (ntrace .gt. 0) call readtrc

      end
