SUBROUTINE init
  !
  ! $Author: CPTEC $
  ! $Date: 2000/02/24 13:15:43 $
  ! $Revision: 1.0 $
  !
  !                .      .    .     
  ! subprogram:    init        initialize variable for model run
  !   prgrmmr: janjic          org: w/np2      date: ??-??-??
  !     
  ! abstract:  init reads in primary and auxiliary variables and constants
  !            and sets initial values for others
  !     
  ! program history log:
  !   87-06-??  janjic  -
  !   92-10-27  deaven  - changed reads of nhb, nfc, and nbc to
  !                       accomodate shortened record lengths
  !   95-03-27  black   - conversion from 1-d to 2-d in horizontal
  !   96-10-31  black   - added namelist bcexdata for the nests
  !
  !     
  ! usage:    call init from main program ebu
  !
  !   input argument list:
  !     none     
  !
  !   output argument list: 
  !     none
  !     
  !   input files:
  !     nfc - the initial values of sfc pressure, t, q, u, and v
  !     nhb - a large variety of array and scalar constants
  !     nbc - the boundary conditions and tendencies
  !
  !                              or
  !
  !     restrt - a restart file with all necessary quantities
  !
  !   output files:
  !     none
  !     
  !   subprograms called:
  !     utilities:
  !       none
  !     library:
  !       common   - ctlblk
  !                  loops
  !                  masks
  !                  dynam
  !                  phys2
  !                  mapot1
  !                  vrbls
  !                  pvrbls
  !                  boco
  !                  grids
  !                  acmclh
  !                  acmcld
  !                  acmpre
  !                  acmrdl
  !                  acmrds
  !                  acmsfc
  !                  cldwtr
  !                  cnvcld
  !                  srfdsp
  !                  soil
  !                  indx
  !                  bcnest
  !    
  IMPLICIT NONE
  INCLUDE "parmeta.f90"
  INCLUDE "parm.tbl.f90"
  REAL, PARAMETER :: d00=0.0
  REAL, PARAMETER :: d50=.50
  REAL, PARAMETER :: h1=1.0
  REAL, PARAMETER :: g=9.8
  REAL, PARAMETER :: cm1=2937.4
  REAL, PARAMETER :: cm2=4.9283
  REAL, PARAMETER :: cm3=23.5518
  REAL, PARAMETER :: eps=0.622
  REAL, PARAMETER :: q2ini=.50
  REAL, PARAMETER :: epsq2=0.2
  REAL, PARAMETER :: epsq=2.e-12
  REAL, PARAMETER :: epswet=0.0
  REAL, PARAMETER :: z0land=.10
  REAL, PARAMETER :: z0sea=.001
  REAL, PARAMETER :: fcm=.00001
  REAL, PARAMETER :: dtr=0.1745329e-1
  REAL, PARAMETER :: h360=360.0
  REAL, PARAMETER :: h1905=190.5
  REAL, PARAMETER :: h105=105.0
  INTEGER, PARAMETER :: lda=lm+9
  INTEGER, PARAMETER :: la=13
  INTEGER, PARAMETER :: imjm=im*jm-jm/2
  INTEGER, PARAMETER :: jmp1=jm+1
  INTEGER, PARAMETER :: jam=6+2*(jm-10)
  INTEGER, PARAMETER :: lb=2*im+jm-3
  INTEGER, PARAMETER :: lm1=lm-1
  INTEGER, PARAMETER :: lp1=lm+1
  INTEGER, PARAMETER :: l3=lda+lm+1
  INTEGER, PARAMETER :: l4=lda+2*lm+1
  INTEGER, PARAMETER :: l7=l4+3
  INTEGER, PARAMETER :: imt=2*im-1
  INTEGER, PARAMETER :: nsoil=2

  !*** declare variables

  REAL :: cloges
  REAL :: dummy
  REAL :: ese
  INTEGER :: iday
  INTEGER :: ihrstb
  INTEGER :: imnth
  INTEGER :: itboco
  INTEGER :: iyr
  INTEGER :: jbcday
  INTEGER :: jbchr
  INTEGER :: jdifhr
  INTEGER :: jmdday
  INTEGER :: jmdhr
  INTEGER :: k
  INTEGER :: khh2
  INTEGER :: khha
  INTEGER :: khl2
  INTEGER :: khla
  INTEGER :: knt
  INTEGER :: kvh2
  INTEGER :: kvha
  INTEGER :: kvl2
  INTEGER :: kvla
  INTEGER :: l1
  INTEGER :: l2
  INTEGER :: llmh
  INTEGER :: nhibu
  INTEGER :: nmap
  INTEGER :: nposit
  INTEGER :: nradlh
  INTEGER :: nradsh
  INTEGER :: nsolar
  INTEGER :: nstart
  INTEGER :: ntim1
  INTEGER :: ntim2
  REAL :: pdomg
  REAL :: pshltr
  REAL :: resomg
  REAL :: tboco
  REAL :: tcp
  REAL :: tddamp
  REAL :: tend
  REAL :: tstart
  LOGICAL runb
  CHARACTER *32 label
  INTEGER idatb(3)
  REAL pslp(im,jm)
  INTEGER i, j, l, n, nrd, ns

  !*** include common blocks.

  INCLUDE "ctlblk.f90"
  INCLUDE "loops.f90"
  INCLUDE "masks.f90"
  INCLUDE "dynam.f90"
  INCLUDE "phys2.f90"
  INCLUDE "mapot1.f90"
  INCLUDE "vrbls.f90"
  INCLUDE "contin.f90"
  INCLUDE "pvrbls.f90"
  INCLUDE "boco.f90"
  INCLUDE "acmclh.f90"
  INCLUDE "acmcld.f90"
  INCLUDE "acmpre.f90"
  INCLUDE "acmrdl.f90"
  INCLUDE "acmrds.f90"
  INCLUDE "acmsfc.f90"
  INCLUDE "cldwtr.f90"
  INCLUDE "cnvcld.f90"
  INCLUDE "srfdsp.f90"
  INCLUDE "solar.f90"
  INCLUDE "soil.f90"
  INCLUDE "indx.f90"
  INCLUDE "z0efft.f90"
  INCLUDE "bcnest.f90"

  !***  declare namelists.

  NAMELIST /fcstdata/ tstart,tend,tcp,restrt,subpost,nmap,tshde,spl, &
       nphs,ncnvc,nradsh,nradlh,tprec,theat,tclod,trdsw,trdlw,tsrfc
  NAMELIST /bcexdata/ nest,tbegp,tbcex, tph0do,tlm0do, &
       imi1,jmi1,lmi1,tph0di1,tlm0di1,imi2,jmi2,lmi2,tph0di2,tlm0di2

  !*** calculate the i-index east-west increments

  DO j=1,jm
     ihe(j)=MOD(j+1,2)
     ihw(j)=ihe(j)-1
     ive(j)=MOD(j,2)
     ivw(j)=ive(j)-1
  ENDDO

  !*** calculate the indirect i indices for radtn

  knt=0
  DO i=1,im
     knt=knt+1
     irad(knt)=i
  ENDDO
  DO i=1,im-1
     knt=knt+1
     irad(knt)=im+2+i
  ENDDO

  !*** zero the external edge of arrays in common.

  CALL zero2(pdsl)
  CALL zero3(t,lm)
  CALL zero3(q,lm)
  CALL zero3(u,lm)
  CALL zero3(v,lm)
  CALL zero2(res)
  CALL zero3(rtop,lm)
  CALL zero3(omgalf,lm)
  CALL zero3(div,lm)
  CALL zero3(etadt,lm-1)
  CALL zero3(htm,lm)
  CALL zero3(vtm,lm)
  CALL zero2(hbm2)
  CALL zero2(akms)
  CALL zero2(uz0)
  CALL zero2(vz0)
  CALL zero2(fad)

  !*** read z0 effective

  READ(42) zeffij

  !*** read "constant" data from unit connected to nhibu.

  nhibu=12
  lsl  =lsm
  READ(nhibu) nfcst,nbc,list,dt,idtad,sigma, &
       khla,khha,kvla,kvha,khl2,khh2,kvl2,kvh2
  list=6
  READ(nhibu)lmh
  READ(nhibu)lmv
  READ(nhibu)((hbm2(i,j),i=1,im),j=1,jm)
  READ(nhibu)vbm2
  READ(nhibu)vbm3
  READ(nhibu)sm
  READ(nhibu)sice
  DO l=1,lm
     READ(nhibu) ((htm(i,j,l),i=1,im),j=1,jm)
  END DO
  DO l=1,lm
     READ(nhibu) ((vtm(i,j,l),i=1,im),j=1,jm)
  END DO
  READ(nhibu) dy,cpgfv,en,ent,r,pt,tddamp, &
       f4d,f4q,ef4t,deta,rdeta,aeta,f4q2,eta,dfl, &
       em,emt
  READ(nhibu) dx
  READ(nhibu) wpdar
  READ(nhibu) cpgfu
  READ(nhibu) curv
  READ(nhibu) fcp
  READ(nhibu) fdiv
  READ(nhibu) ((fad(i,j),i=1,im),j=1,jm)
  READ(nhibu) f
  READ(nhibu) ddmpu
  READ(nhibu) ddmpv
  READ(nhibu) pt2,glat
  READ(nhibu) glon
  READ(nhibu) plq,rdpq,rdtheq,stheq,the0q
  READ(nhibu) ros,cs,ds,roi,ci,di, &
       pl,thl,rdq,rdth,rdp,rdthe, &
       deta2,aeta2,dfrlg, &
       qs0,sqs,sthe,the0
  READ(nhibu) wfk
  READ(nhibu) epsr
  READ(nhibu) tg
  READ(nhibu) gffc
  READ(nhibu) sst
  READ(nhibu) alb
  READ(nhibu) hdac
  READ(nhibu) hdacv
  READ(nhibu) ttblq
  READ(nhibu) ptbl,ttbl,r1,pt1,tsph, &
       wbd,sbd,tlm0d,tph0d,dlmd,dphd,cmld,dp30, &
       x1p,y1p,ixm,iym,deta1,aeta1,eta1
  READ(nhibu) ivgtyp
  READ(nhibu) isltyp
  READ(nhibu) islope
  READ(nhibu) vegfrc
  READ(nhibu) sldpth
  READ(nhibu) rtdpth

  !*** read annual values of the greenwich mean sidereal time
  !*** for 0000 utc december 31 and of the geometric mean
  !*** celestial longitude of the sun for 0000 tdt december 31

  nsolar=19
  l2=0
  DO n=1,6
     l1=l2+1
     l2=l1+4
     READ(nsolar,25)(csidx(l),l=l1,l2)
25   FORMAT(5f9.6)
  ENDDO
  l2=0
  DO n=1,6
     l1=l2+1
     l2=l1+4
     READ(nsolar,30)(clonx(l),l=l1,l2)
30   FORMAT(5f11.6)
  ENDDO

  nhibu = 12
  WRITE(list,*)'init:  read constants file'

  !*** read namelist fcstdata which controls timesteps, 
  !*** accumulation periods, standard output

  restrt = .FALSE.
  REWIND 11
  READ(11,fcstdata)
  WRITE(list,*)'init:  read namelist fcstdata - listed below'
  WRITE(list,*)'  tstart,tend  :  ',tstart,tend
  WRITE(list,*)'  tcp          :  ',tcp
  WRITE(list,*)'  restrt       :  ',restrt
  WRITE(list,*)'  subpost      :  ',subpost
  WRITE(list,*)'  nmap,nphs    :  ',nmap,nphs
  WRITE(list,*)'  ncnvc        :  ',ncnvc
  WRITE(list,*)'  nradsh,nradlh:  ',nradsh,nradlh
  WRITE(list,*)'  tprec,theat  :  ',tprec,theat
  WRITE(list,*)'  tclod,trdsw  :  ',tclod,trdsw
  WRITE(list,*)'  trdlw,tsrfc  :  ',trdlw,tsrfc
  WRITE(list,*)'  tshde (posted forecast hours) below:  '
  WRITE(list,75) (tshde(k),k=1,49)
  WRITE(list,*)'  spl (posted pressure levels) below: '
  WRITE(list,80) (spl(l),l=1,lsm)
75 FORMAT(14(f4.1,1x))
80 FORMAT(8(f8.1,1x))

  !*** set time stepping related constants.

  first  = .TRUE.
  nstart = INT(tstart*tsph+d50)
  ntstm  = INT(tend  *tsph+d50)+1
  ncp    = INT(tcp   *tsph+d50)
  nddamp = INT(tddamp*tsph+d50)
  nprec  = INT(tprec *tsph+d50)
  nheat  = INT(theat *tsph+d50)
  nclod  = INT(tclod *tsph+d50)
  nrdsw  = INT(trdsw *tsph+d50)
  nrdlw  = INT(trdlw *tsph+d50)
  nsrfc  = INT(tsrfc *tsph+d50)

  !*** set various physics package timestep variables.

  nrads = NINT(tsph)*nradsh
  nradl = NINT(tsph)*nradlh
  dtq2  = nphs * dt
  tdtq2 = dtq2 + dtq2
  dtd   = d50  * dtq2
  tdtd  = dtd  + dtd
  ktm   = INT(dtq2/dtd+d50)
  WRITE(list,*)' '
  WRITE(list,*)'set time stepping constants'
  WRITE(list,*)' first             :  ',first
  WRITE(list,*)' nstart,nstsm,ncp  :  ',nstart,ntstm,ncp
  WRITE(list,*)' nddamp,nprec,nheat:  ',nddamp,nprec,nheat
  WRITE(list,*)' nclod,nrdsw,nrdlw :  ',nclod,nrdsw,nrdlw
  WRITE(list,*)' nsrfc             :  ',nsrfc
  WRITE(list,*)' nrads,nradl,ktm   :  ',nrads,nradl,ktm
  WRITE(list,*)' dtq2,tdtq2        :  ',dtq2,tdtq2
  WRITE(list,*)' dtd,tdtd          :  ',dtd,tdtd
  WRITE(list,*)' '

  !*** read namelist bcexdata which controls the generation
  !*** of extraction of boundary conditions for the
  !*** nested run.

  REWIND 17
  READ(17,bcexdata)
  WRITE(list,*)'init:  read namelist bcexdata - listed below'
  WRITE(list,*)'  nest,tbegp  :  ',nest,tbegp
  WRITE(list,*)'  tbcex (nested fcst bc times) below:  '
  WRITE(list,75) (tbcex(k),k=1,50)
  kntbc=1

  !*** compute derived map output constants.

  DO l = 1,lsl
     alsl(l) = LOG(spl(l))
  END DO
  DO i=1,nmap
     ishde(i)=INT(tshde(i)*tsph+d50)+1
  END DO

  !*** read initial conditions or restart file.
  !*** first, the .not. restart file case.

  IF (.NOT.restrt) THEN
     READ(nfcst) run,idat,ihrst,ntsd
     READ(nfcst) pd
     READ(nfcst) ((res(i,j),i=1,im),j=1,jm)
     READ(nfcst) fis
     DO l=1,lm
        READ(nfcst) ((u(i,j,l),i=1,im),j=1,jm)
     END DO
     DO l=1,lm
        READ(nfcst) ((v(i,j,l),i=1,im),j=1,jm)
     END DO
     DO l=1,lm
        READ(nfcst) ((t(i,j,l),i=1,im),j=1,jm)
     END DO
     DO l=1,lm
        READ(nfcst) ((q(i,j,l),i=1,im),j=1,jm)
     END DO
     READ(nfcst) wet
     READ(nfcst) sno
     READ(nfcst) smc
     READ(nfcst) cmc
     READ(nfcst) stc
     iyr   = idat(3)-1900
     imnth = idat(1)
     iday  = idat(2)
     WRITE(list,*)'init:  read initial condition file'
  ELSE

     !*** second, the restart file case.

     WRITE(list,*)'init:  read restart file'
     READ(nfcst) run,idat,ihrst,ntsd,label
     READ(nfcst) pdomg,resomg
     DO l = 1,lm
        READ(nfcst) ((omgalf(i,j,l),i=1,im),j=1,jm)
     END DO
     WRITE(list,*)'  read  ',label
     READ(nfcst) run,idat,ihrst,ntsd,label, &
          first,iout,nshde
     first=.TRUE.
     READ(nfcst) pd,((res(i,j),i=1,im),j=1,jm),fis
     READ(nfcst) pdb,tb,qb,ub,vb
     DO l = 1,lm
        READ(nfcst) ((t(i,j,l),i=1,im),j=1,jm)
        READ(nfcst) ((q(i,j,l),i=1,im),j=1,jm)
        READ(nfcst) ((u(i,j,l),i=1,im),j=1,jm)
        READ(nfcst) ((v(i,j,l),i=1,im),j=1,jm)
        READ(nfcst) ((q2(i,j,l),i=1,im),j=1,jm)
        !*** dummy read of the total radiative temperature tendencies
        !*** which are not used explicitly in the integration
        READ(nfcst)
        READ(nfcst) ((cwm(i,j,l),i=1,im),j=1,jm)
        READ(nfcst) ((train(i,j,l),i=1,im),j=1,jm)
        READ(nfcst) ((tcucn(i,j,l),i=1,im),j=1,jm)
     END DO
     WRITE(list,*)'  read  ',label
     READ(nfcst) run,idat,ihrst,ntsd,label, &
          rswin,rswout,tg,z0,((akms(i,j),i=1,im),j=1,jm),czen
     READ(nfcst) akhs,ths,qs,twbs,qwbs,hbot,cfracl
     READ(nfcst) thz0,qz0, & 
          ((uz0(i,j),i=1,im),j=1,jm), &
          ((vz0(i,j),i=1,im),j=1,jm), &
          ustar,htop,cfracm
     READ(nfcst) sno,wet,cldefi,rf,pslp,cuppt,cfrach
     READ(nfcst) soiltb,sfcexc,smstav,smstot,grnflx,pctsno
     READ(nfcst) rlwin,radot,czmean,sigt4
     READ(nfcst) u00,ul,lc,sr
     READ(nfcst) run,idat,ihrst,ntsd,label, &
          prec,acprec,accliq,cuprec
     READ(nfcst) acfrcv,ncfrcv,acfrst,ncfrst
     READ(nfcst) acsnow,acsnom,ssroff,bgroff
     READ(nfcst) sfcshx,sfclhx,subshx,snopcx, &
          sfcuvx,sfcevp,potevp
     READ(nfcst) aswin,aswout,aswtoa,alwin,alwout,alwtoa
     READ(nfcst) ardsw,ardlw,asrfc,avrain,avcnvc
     READ(nfcst) th10,q10,u10,v10,tshltr,qshltr,pshltr
     READ(nfcst) smc
     READ(nfcst) cmc
     READ(nfcst) stc
     WRITE(list,*)'  read  ',label
     !*** call radiation to obtain the short and longwave
     !*** temperature tendencies
     CALL radtn
  ENDIF

  !*** done reading initial conditions or a restart file.
  !*** if not running the model, print date of initial conditions
  !*** just read and stop.  otherwise, continue.

  IF (.NOT. run) THEN
     WRITE(list,165) ihrst,idat
     WRITE(list,166)
     STOP 2
165  FORMAT('0*** ',i2,' gmt ',2(i2,'/'),i4,' ***')
166  FORMAT('0f*** no initial conditions. run terminated.')
  END IF

  !*** if the timestep counter (ntsd) exceeds the "stop model" t
  !*** timestep,continue, stop execution.  otherwise, continue.

  IF (ntsd.GE.ntstm) THEN
     WRITE(list,165) ihrst,idat
     WRITE(list,195)
     STOP 3
195  FORMAT('0f*** forecast already done. run terminated.')
     !*** read boundary conditions.
  END IF
  REWIND nbc
  READ (nbc) runb,idatb,ihrstb,tboco
  nboco = INT(tboco*tsph+d50)
  CALL w3fs13(idatb(3),idatb(1),idatb(2),jbcday)
  CALL w3fs13(idat(3),idat(1),idat(2),jmdday)
  jbchr  = jbcday*24 + ihrstb
  jmdhr  = jmdday*24 + ihrst
  jdifhr = jmdhr     - jbchr
  WRITE(list,*)'init:  read boundary conditions'
  IF (jdifhr.NE.0) THEN
     itboco = NINT(tboco)
     nposit = jdifhr/itboco+1
     WRITE(list,*)'  loop nposit=',nposit,' times to read bc'
     IF (restrt) THEN
        WRITE(list,*)'  ',nposit,' dummy reads into nbc=', &
             nbc,'.  use bc in ic.'
     ELSE
        WRITE(list,*)'init:  ',nposit,' reads into nbc=', &
             nbc,' for bc'
     ENDIF
     DO nrd = 1,nposit
        IF (restrt) THEN
           READ(nbc) dummy
           READ(nbc) dummy
           READ(nbc) dummy
           READ(nbc) dummy
           READ(nbc) dummy
        ELSE
           READ(nbc) pdb
           READ(nbc) tb
           READ(nbc) qb
           READ(nbc) ub
           READ(nbc) vb
        ENDIF
     END DO
  ELSE 
     READ(nbc) pdb
     READ(nbc) tb
     READ(nbc) qb
     READ(nbc) ub
     READ(nbc) vb
     WRITE(list,*)'  read unit nbc=',nbc
  ENDIF
  IF (MOD(ntsd,nboco).EQ.0) THEN
     BACKSPACE nbc
     BACKSPACE nbc
     BACKSPACE nbc
     BACKSPACE nbc
     BACKSPACE nbc
     WRITE(list,*)'  backspace unit nbc=',nbc
  ENDIF

  !*** set arrays controlling post processing.

  WRITE(list,*)'init:  read iout,nshde,ntsd=',iout,nshde,ntsd
  DO i = 1,nmap
     iout = i
     IF (ishde(i).GE.ntsd) EXIT
  END DO
  nshde = ishde(iout)
  WRITE(list,*)'init:  set iout,nshde =',iout,nshde, &
       ' for ishde,ntsd=',ishde(iout),ntsd
  aphtim=0.0
  aratim=0.0
  acutim=0.0

  !*** initialize physics variables if starting this run from scratch.

  IF(.NOT.restrt)THEN
     DO j=1,jm
        DO i=1,im
           prec(i,j)   = d00
           acprec(i,j) = d00
           cuprec(i,j) = d00
           z0(i,j)     = sm(i,j)*z0sea+(h1-sm(i,j))* &
                (fis(i,j)*fcm+z0land)
           qs(i,j)     = d00
           akms(i,j)   = d00
           akhs(i,j)   = d00
           twbs(i,j)   = d00
           qwbs(i,j)   = d00
           cldefi(i,j) = h1
           htop(i,j)   = REAL(lmh(i,j))
           hbot(i,j)   = REAL(lmh(i,j))
           !*** this was moved from surfce to here
           !*** reduce the high gdas soil moisture 
           DO ns=1,nsoil
              IF(smc(i,j,ns).GT.0.31)THEN
                 smc(i,j,ns)=MIN(0.714*smc(i,j,ns)+0.089,0.36)
              ENDIF
           ENDDO
           !***  end of gdas soil moisture adjustment
        ENDDO
     ENDDO
     !*** initialize cloud fields
     DO l=1,lm
        DO j=1,jm
           DO i=1,im
              cwm(i,j,l)=d00
           ENDDO
        ENDDO
     ENDDO
     ntim1=0.*tsph
     ntim2=24.*tsph
     !*** initialize accumulator arrays to zero.
     ardsw=0.0
     ardlw=0.0
     asrfc=0.0
     avrain=0.0
     avcnvc=0.0
     DO j=1,jm
        DO i=1,im
           acfrcv(i,j) = d00
           ncfrcv(i,j) = 0
           acfrst(i,j) = d00
           ncfrst(i,j) = 0
           acsnow(i,j) = d00
           acsnom(i,j) = d00
           ssroff(i,j) = d00
           bgroff(i,j) = d00
           alwin(i,j)  = d00
           alwout(i,j) = d00
           alwtoa(i,j) = d00
           aswin(i,j)  = d00
           aswout(i,j) = d00
           aswtoa(i,j) = d00
           sfcshx(i,j) = d00
           sfclhx(i,j) = d00
           subshx(i,j) = d00
           snopcx(i,j) = d00
           sfcuvx(i,j) = d00
           sfcevp(i,j) = d00
           potevp(i,j) = d00
           potflx(i,j) = d00
        ENDDO
     ENDDO
     !*** initialize saturation specific humidity over the water.
     DO j=1,jm
        DO i=1,im
           cloges = -cm1/sst(i,j)-cm2*alog10(sst(i,j))+cm3
           ese    = 10.**(cloges+2.)
           qs(i,j)= sm(i,j)*eps*ese/(pd(i,j)+pt-ese*(1.-eps))
        ENDDO
     ENDDO
     !*** pad ground wetness if it is too small.
     DO j=1,jm
        DO i=1,im
           wet(i,j)=amax1(wet(i,j),epswet)
        ENDDO
     ENDDO
     !*** initialize turbulent kinetic energy (tke) to a small
     !*** value (epsq2) above ground.  set tke to zero in the
     !*** the lowest model layer.  in the lowest two atmospheric
     !*** eta layers set tke to a small value (q2ini).
     DO l=1,lm1
        DO j=1,jm
           DO i=1,im
              q2(i,j,l)=htm(i,j,l+1)*hbm2(i,j)*epsq2
           ENDDO
        ENDDO
     ENDDO
     DO j=1,jm
        DO i=1,im
           q2(i,j,lm)    = d00
           llmh          = lmh(i,j)
           q2(i,j,llmh-2)= hbm2(i,j)*q2ini
           q2(i,j,llmh-1)= hbm2(i,j)*q2ini
        ENDDO
     ENDDO
     !*** pad above ground specific humidity if it is too small.
     !*** initialize latent heating accumulation arrays.
     DO l=1,lm
        DO j=1,jm
           DO i=1,im
              IF(q(i,j,l).LT.epsq)q(i,j,l)=epsq*htm(i,j,l)
              train(i,j,l)=d00
              tcucn(i,j,l)=d00
           ENDDO
        ENDDO
     ENDDO
     !*** end of scratch start initialization block.
     WRITE(list,*)'init:  initialized arrays for clean start'
  ENDIF

  !*** restart initializing.  check to see if we need to zero
  !*** accumulation arrays.

  IF(restrt)THEN
     !***   set ntim1 to number of hours in preceding edas
     ntim1=3.*tsph
     ntim2=24.*tsph
     !*** average cloud amount array
     IF(MOD(ntsd,nclod).LT.nphs)THEN
        WRITE(list,*)'  zero avg cld amt array'
        DO j=1,jm
           DO i=1,im
              acfrcv(i,j)=d00
              ncfrcv(i,j)=0
              acfrst(i,j)=d00
              ncfrst(i,j)=0
           ENDDO
        ENDDO
     ENDIF
     !*** grid-scale and convective latent heating arrays.
     IF(MOD(ntsd,nheat).LT.ncnvc)THEN
        WRITE(list,*)'  zero accum latent heating arrays'
        avrain=d00
        avcnvc=d00
        DO l=1,lm
           DO j=1,jm
              DO i=1,im
                 train(i,j,l)=d00
                 tcucn(i,j,l)=d00
              ENDDO
           ENDDO
        ENDDO
     ENDIF
     !*** total and convective precipitation arrays.
     !*** total snow and snow melt arrays.
     !*** storm surface and base ground run off arrays.
     IF(MOD(ntsd,nprec).LT.nphs)THEN
        WRITE(list,*)'  zero accum precip arrays'
        DO j=1,jm
           DO i=1,im
              acprec(i,j)=d00
              cuprec(i,j)=d00
              acsnow(i,j)=d00
              acsnom(i,j)=d00
              ssroff(i,j)=d00
              bgroff(i,j)=d00
           ENDDO
        ENDDO
     ENDIF
     !*** long wave radiation arrays.
     IF(MOD(ntsd,nrdlw).LT.nphs)THEN
        WRITE(list,*)'  zero accum lw radtn arrays'
        ardlw=d00
        DO j=1,jm
           DO i=1,im
              alwin(i,j) =d00
              alwout(i,j)=d00
              alwtoa(i,j)=d00
           ENDDO
        ENDDO
     ENDIF
     !*** short wave radiation arrays.
     IF(MOD(ntsd,nrdsw).LT.nphs)THEN
        WRITE(list,*)'  zero accum sw radtn arrays'
        ardsw=d00
        DO j=1,jm
           DO i=1,im
              aswin(i,j) =d00
              aswout(i,j)=d00
              aswtoa(i,j)=d00
           ENDDO
        ENDDO
     ENDIF
     !*** surface sensible and latent heat flux arrays.
     IF(MOD(ntsd,nsrfc).LT.nphs)THEN
        WRITE(list,*)'  zero accum sfc flux arrays'
        asrfc=d00
        DO j=1,jm
           DO i=1,im
              sfcshx(i,j)=d00
              sfclhx(i,j)=d00
              subshx(i,j)=d00
              snopcx(i,j)=d00
              sfcuvx(i,j)=d00
              sfcevp(i,j)=d00
              potevp(i,j)=d00
              potflx(i,j)=d00
           ENDDO
        ENDDO
     ENDIF
     !*** endif for restart file accumulation zero block.
     WRITE(list,*)'init:  initialized arrays for restart start'
  ENDIF

  !*** initialize cloud constants

  DO j=1,jm
     DO i=1,im
        u00(i,j)=(h1-sm(i,j))*0.75+sm(i,j)*0.80
     END DO
  END DO
  DO l=1,2*lm
     IF(l.GE.lm-10.AND.l.LE.lm)THEN
        ul(l)=0.1*float(l-lm+10)
     ELSE
        ul(l)=d00
     ENDIF
  END DO

  !*** set index arrays for upstream advection

  knt=0
  DO j=3,5
     knt=knt+1
     ihla(knt)=2
     ihha(knt)=im-1-MOD(j+1,2)
     ivla(knt)=2 
     ivha(knt)=im-1-MOD(j,2)
     jra(knt)=j
  ENDDO
  DO j=jm-4,jm-2
     knt=knt+1
     ihla(knt)=2
     ihha(knt)=im-1-MOD(j+1,2)
     ivla(knt)=2 
     ivha(knt)=im-1-MOD(j,2)
     jra(knt)=j
  ENDDO
  DO j=6,jm-5
     knt=knt+1
     ihla(knt)=2
     ihha(knt)=2+MOD(j,2)
     ivla(knt)=2 
     ivha(knt)=2+MOD(j+1,2)
     jra(knt)=j
  ENDDO
  DO j=6,jm-5
     knt=knt+1
     ihla(knt)=im-2
     ihha(knt)=im-2+MOD(j,2)
     ivla(knt)=im-2
     ivha(knt)=im-2+MOD(j+1,2)
     jra(knt)=j
  ENDDO

  !*** end of subroutine init.

  WRITE(list,*)'init:  exit init and start model integration'
  WRITE(list,*)' '
END SUBROUTINE init
