c    @(#)initial.f	1.4
      subroutine initial
c
c--------------------------------------------------------------------
c
      include'mladv.inc'
c
c
c--------------------------------------------------------------------
c
c -- qg common
      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),scsig,sigz(k01),eigval(k0),ainv(k0,k0),
     .               amat(k0,k0)
      common/fields/vort(mn0,k0),stream(mn0,k0),prvjac(mn0,k0),
     +     topden(mn0),botden(mn0),relief(mn0),prvtop(mn0),
     +     prvbot(mn0),topfct(k0),botfct(k0),pvort(mn0)
      common /extpar/ rlat0, rlng0, t0, v0, dhor, ht, time0,
     +                r0,f0star
      common /vert/ rho0, tsurf, xc(k01), salsur, eta(k01)
      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

      if( ifsbl.eq.1) then
c
c --- simple initialization
c
         write(iout,500)
c
c --- init mixed layer depth.
c
         do 10 j=1,n
            do 9 i=1,m
               ip=i+(j-1)*m
               hdpth(ip,mlcur)=hdpth0
               if (ifperh .eq. 1) hdpth(ip,mlpast) = hdpth0
 9          continue
 10      continue

c
c --- now bouyancy and temperature.
c
         do 20 k=1,l
             do 15 ip=1,msq
               buoy(ip,k,mlcur)=0.0
               temp(ip,k,mlcur)=0.0
 15         continue
 20      continue
         write(iout,211)hdpth0
c
      elseif(ifsbl.eq.2) then
         write(iout,501)
c
c -- read initial b T and h
c
         call rdbth(ier)
         if(ier.gt.0) then
            write(iout,101)
            print*, ' error while reading sbl.th.dat file '
            stop
         endif
c
c -- convert from  total b and T to SBL q and T anomalies
c         
         do 40 k=1,l
             q = -zml(k)
             do 35 ip=1,msq
               buoy(ip,k,mlcur)= buoy(ip,k,mlcur) - bqgofk(ip,q)
               temp(ip,k,mlcur)= temp(ip,k,mlcur) - tqgofk(ip,q)
 35         continue
 40      continue

      else if (ifsbl.eq.3) then
c
c --- simple initialization - option 2
c
         write(iout,500)
c
c --- init mixed layer depth.
c
         do 12 j=1,n
            do 11 i=1,m
               ip=i+(j-1)*m
               hdpth(ip,mlcur)=hdpth0
               if (ifperh .eq. 1) hdpth(ip,mlpast) = hdpth0
 11         continue
 12      continue

c
c --- now bouyancy and temperature.
c
         dbdz=(f0star*dhor/ht)**2./(scsig*sigz(1))
         bsurf2 = buoysur
         buoysur = buoysur + (hdpth0 * dbdz)
         tsurf2 = tsurf
         tsurf = tsurf + (hdpth0 * xc(1))

         call findkml(hdpth0,zml,lmax,kml)
         
         do 21 k=1,l
             do 16 ip=1,msq
               if (k .le. kml) then
                  buoy(ip,k,mlcur) = bsurf2 - bqgofk(ip,k)
                  temp(ip,k,mlcur) = tsurf2 - tqgofk(ip,k)
                 else
                  buoy(ip,k,mlcur)=0.0
                  temp(ip,k,mlcur)=0.0
                endif
 16         continue
 21      continue

         write(iout,212)hdpth0

      else if (ifsbl.ge.4) then
c
c --- simple initialization - option 3
c
         write(iout,500)
c
c --- init mixed layer depth.
c
         do 14 j=1,n
            do 13 i=1,m
               ip=i+(j-1)*m
               hdpth(ip,mlcur)=hdpth0
               if (ifperh .eq. 1) hdpth(ip,mlpast) = hdpth0
 13         continue
 14      continue

c
c --- now bouyancy and temperature.
c
         bscale=f0star*v0*dhor/(ht*scsig*sigz(1))
         dbdz=(f0star*dhor/ht)**2./(scsig*sigz(1))
         tscale=xc(1)*v0*ht/(f0star*dhor)

         call findkml(hdpth0,zml,lmax,kml)
         
         do 22 ip=1,msq
            delta = topden(ip) + artop
            bml = buoysur + bscale * delta - hdpth0 * dbdz
            tml = tsurf   + tscale * delta - hdpth0 * xc(1)
            do 17 k=1,l
               if (k .le. kml) then
                  buoy(ip,k,mlcur) = bml - bqgofk(ip,k)
                  temp(ip,k,mlcur) = tml - tqgofk(ip,k)
                 else
                  buoy(ip,k,mlcur)=0.0
                  temp(ip,k,mlcur)=0.0
                endif
 17         continue
 22      continue

         if (ifrbg .eq. 1) call readbg

         write(iout,213)hdpth0
c
      else
         write(iout,502)ifsbl
         print*, 'wrong sbl key'
         stop
      endif
c
c -- tracer initialization
c
      if (ntrace.gt.0) call trcinit

      write(iout,222)wflux0,qsw0,qsurf0

      write(iout,201)u10m,angu10m

 211  format
     .   (' the initial conditions for the mixed layer model are',
     . /,'  mixed layer depth of ',f4.1,' meters and no buoyancy or',
     .     /,'  temperature deviations other than those needed to  ',
     .     /,'  make the mixed layer well mixed.',///)
 212  format
     .   (' the initial conditions for the mixed layer model are',
     . /,'mixed layer depth of ',f4.1,' meters with ML buoyancy',
     . /,'and temperature matching that set by (tsurf,salsur) ',
     . /,'and no discontinuity at the base of the ML.',///)
 213  format
     .   (' the initial conditions for the mixed layer model are',
     . /,'mixed layer depth of ',f4.1,' meters with ML buoyancy',
     . /,'and temperature matching that set by the qg values at',
     . /,'the base of the ML.',///)
 201  format(
     .       '  wind speed of                 ',f6.2,' m/s',/,
     .       '  wind direction                ',f6.2,
     .     ' deg cclk from x-positive'//)
 222  format('  water flux           ',f8.4,' m/s ',/
     .      , '  short wave radiation ',f8.4,
     .          ' cal/cm^2/sec pos. down',/
     .      , '  heat surface flux    ',f8.4,
     .        ' cal/cm^2/sec pos. down ',/)

 101  format(' error while reading sbl.th.dat file ')
 500  format(/,'   SURFACE BOUNDARY LAYER INPUT SECTION',/,
     *'   Initial  b and T anomalies set to zero, h is a constant')
 501  format(/,'   SURFACE BOUNDARY LAYER INPUT SECTION',/
     *        ,'   Initial   T and h read: form b & T anomalies  ')
 502  format(/,' Fatal error: wrong ifsbl selection: ',i5)
      end
