SUBROUTINE bcex(kntbc,tbcex,nwgts,iutim0,iutim1,tph0do,tlm0do, &
     imi,jmi,lmi,kbi,tph0di,tlm0di)
  !
  ! $Author: CPTEC $
  ! $Date: 2000/04/17 12:13:43 $
  ! $Revision: 1.0 $
  !
  IMPLICIT NONE
  INCLUDE "parmeta.f90"
  INCLUDE "ctlblk.f90"
  INCLUDE "vrbls.f90"
  INCLUDE "pvrbls.f90"
  INCLUDE "cldwtr.f90"
  INCLUDE "masks.f90"
  INCLUDE "indx.f90"

  REAL :: alpha,cphi0,d2r,erl0,erlam0,erphi0,etop,h1,h2
  REAL :: h3,h4,htot
  REAL :: pi,r2d,rsumh,rsumv,rtboco,sphi0,tboco,tlm0di,tlm0do
  REAL :: tph0di,tph0do,v1,v2,v3,v4,vtot
  INTEGER :: i1,i2,i3,i4,ih,imi,in,iutim0,iutim1,iv,j1,j2,j3,j4
  INTEGER :: jh,jmi,jn,jv,k,kbi,kntbc,l,lmi,lmip,n,nwgts

  !*** the inner grid 

  REAL :: pdbi(kbi,2)
  REAL :: tbi(kbi,lmi,2)
  REAL :: qbi(kbi,lmi,2)
  REAL :: ubi(kbi,lmi,2)
  REAL :: vbi(kbi,lmi,2)
  REAL :: q2bi(kbi,lmi,2)
  REAL :: cwmbi(kbi,lmi,2)

  !*** the outer grid 

  REAL :: pdbo(kbi)
  REAL :: tbo(kbi,lm)
  REAL :: qbo(kbi,lm)
  REAL :: ubo(kbi,lm)
  REAL :: vbo(kbi,lm)
  REAL :: q2bo(kbi,lm)
  REAL :: cwmbo(kbi,lm)
  REAL :: hlati(kbi)
  REAL :: hloni(kbi)
  REAL :: vlati(kbi)
  REAL :: vloni(kbi)
  REAL :: tvloni(kbi)
  REAL :: tvlati(kbi)
  INTEGER :: iouthb(kbi)
  INTEGER :: jouthb(kbi)
  INTEGER :: ioutvb(kbi)
  INTEGER :: joutvb(kbi)
  INTEGER :: ionear(kbi)
  INTEGER :: jonear(kbi)
  REAL :: hwgts(kbi,4)
  REAL :: vwgts(kbi,4)
  REAL :: uerth(kbi,lm)
  REAL :: verth(kbi,lm)
  REAL :: cosalp(kbi)
  REAL :: sinalp(kbi)
  REAL :: pdbi0(kbi)
  REAL :: tbi0(kbi,lmi)
  REAL :: qbi0(kbi,lmi)
  REAL :: ubi0(kbi,lmi)
  REAL :: vbi0(kbi,lmi)
  REAL :: q2bi0(kbi,lmi)
  REAL :: cwmbi0(kbi,lmi)
  REAL :: tbcex(50)

  !***  read in etop via namelist

  etop=25.
10 CONTINUE
  lmip=lmi+1

  !***  read geodetic lat/lon at boundary of inner grid

  REWIND nwgts
  READ(nwgts)hlati,hloni,vlati,vloni

  !***  read in locations of the inner boundary within the outer grid
  !***  as well as interpolation weights

  READ(nwgts)iouthb,jouthb,ioutvb,joutvb,ionear,jonear
  READ(nwgts)hwgts,vwgts

  pi=3.141592654
  d2r=pi/180.
  r2d=1./d2r

  !***  compute rotation angle for winds (transformed grid (outer) to
  !***  earth lat/lon at inner grid points)

  erlam0=(360.-tlm0do)*d2r
  erphi0=tph0do*d2r
  erl0=erlam0/d2r
  cphi0=COS(erphi0)
  sphi0=SIN(erphi0)
  CALL tll(vloni,vlati,kbi,erl0,d2r,cphi0,sphi0,tvloni,tvlati)

  DO k=1,kbi
     alpha=ASIN((sphi0*SIN(tvloni(k)*d2r)) &
          /(COS(vlati(k)*d2r)))
     cosalp(k)=COS(alpha)
     sinalp(k)=SIN(alpha)
  ENDDO

  !***  perform interpolation from outer grid to inner grid 
  !***  boundary points

  !***  take pd on the inner grid boundary to be equal to the
  !***  value of pd on the outer grid point that is nearest to
  !***  the inner grid point, i.e., no horizontal interpolation.

  DO k=1,kbi
     in=ionear(k)
     jn=jonear(k)
     pdbo(k)=pd(in,jn)
  ENDDO

  DO l=1,lm

     DO k=1,kbi

        !***  the other h point arrays are interpolated.

        ih=iouthb(k)
        jh=jouthb(k)
        i1=ih
        j1=jh
        i2=ih+1
        j2=jh
        i3=ih+ihe(jh)
        j3=jh-1
        i4=ih+ihe(jh)
        j4=jh+1

        h1=htm(i1,j1,l)
        h2=htm(i2,j2,l)
        h3=htm(i3,j3,l)
        h4=htm(i4,j4,l)
        htot=h1+h2+h3+h4
        IF(htot.GT.0.)THEN
           rsumh=1./(h1*hwgts(k,1)+h2*hwgts(k,2) &
                +h3*hwgts(k,3)+h4*hwgts(k,4))
        ELSE
           rsumh=1.
        ENDIF

        tbo(k,l)=(hwgts(k,1)*t(i1,j1,l)*h1 &
             +hwgts(k,2)*t(i2,j2,l)*h2 &
             +hwgts(k,3)*t(i3,j3,l)*h3 &
             +hwgts(k,4)*t(i4,j4,l)*h4)*rsumh
        qbo(k,l)=(hwgts(k,1)*q(i1,j1,l)*h1 &
             +hwgts(k,2)*q(i2,j2,l)*h2 &
             +hwgts(k,3)*q(i3,j3,l)*h3 &
             +hwgts(k,4)*q(i4,j4,l)*h4)*rsumh
        q2bo(k,l)=(hwgts(k,1)*q2(i1,j1,l)*h1 &
             +hwgts(k,2)*q2(i2,j2,l)*h2 &
             +hwgts(k,3)*q2(i3,j3,l)*h3 &
             +hwgts(k,4)*q2(i4,j4,l)*h4)*rsumh
        cwmbo(k,l)=(hwgts(k,1)*cwm(i1,j1,l)*h1 &
             +hwgts(k,2)*cwm(i2,j2,l)*h2 &
             +hwgts(k,3)*cwm(i3,j3,l)*h3 &
             +hwgts(k,4)*cwm(i4,j4,l)*h4)*rsumh

        !***  the v points

        iv=ioutvb(k)
        jv=joutvb(k)

        i1=iv
        j1=jv
        i2=iv+1
        j2=jv
        i3=iv+ive(jv)
        j3=jv-1
        i4=iv+ive(jv)
        j4=jv+1

        v1=vtm(i1,j1,l)
        v2=vtm(i2,j2,l)
        v3=vtm(i3,j3,l)
        v4=vtm(i4,j4,l)
        vtot=v1+v2+v3+v4
        IF(vtot.GT.0.)THEN
           rsumv=1./(v1*vwgts(k,1)+v2*vwgts(k,2) &
                +v3*vwgts(k,3)+v4*vwgts(k,4))
        ELSE
           rsumv=1.
        ENDIF

        ubo(k,l)=(vwgts(k,1)*u(i1,j1,l)*v1 &
             +vwgts(k,2)*u(i2,j2,l)*v2 &
             +vwgts(k,3)*u(i3,j3,l)*v3 &
             +vwgts(k,4)*u(i4,j4,l)*v4)*rsumv
        vbo(k,l)=(vwgts(k,1)*v(i1,j1,l)*v1 &
             +vwgts(k,2)*v(i2,j2,l)*v2 &
             +vwgts(k,3)*v(i3,j3,l)*v3 &
             +vwgts(k,4)*v(i4,j4,l)*v4)*rsumv
     END DO
  END DO

  !***  begin rotation of wind components.
  !***  first, rotate from outer grid coord to earth coordinates.

  DO l=1,lm
     DO k=1,kbi
        uerth(k,l)=ubo(k,l)*cosalp(k)+vbo(k,l)*sinalp(k)
        verth(k,l)=vbo(k,l)*cosalp(k)-ubo(k,l)*sinalp(k)
     ENDDO
  ENDDO

  !***  rotate from earth coordinates to the target grid

  CALL rotlle(uerth,verth,vlati,vloni,tph0di,tlm0di,kbi,ubo,vbo)

  !***  perform vertical interpolation of boundary point fields 

  CALL ptetae(etop,ionear,jonear, & 
       imi,jmi,lmi,lmip,kbi, &
       pdbi,tbi,qbi,ubi,vbi,q2bi,cwmbi, &
       pdbo,tbo,qbo,ubo,vbo,q2bo,cwmbo,ntsd)

  !***  read in inner grid boundary values from tboco hours ago,
  !***  form the tendencies, and write out to the boundary file
  !***  in the standard format

  REWIND iutim0
  IF(kntbc.GT.1)THEN
     READ(iutim0)pdbi0,tbi0,qbi0,ubi0,vbi0,q2bi0,cwmbi0
  ENDIF

  REWIND iutim0
  WRITE(iutim0)(pdbi(n,1),n=1,kbi), &
       ((tbi(n,l,1),n=1,kbi),l=1,lmi), &
       ((qbi(n,l,1),n=1,kbi),l=1,lmi), &
       ((ubi(n,l,1),n=1,kbi),l=1,lmi), &
       ((vbi(n,l,1),n=1,kbi),l=1,lmi), &
       ((q2bi(n,l,1),n=1,kbi),l=1,lmi), &
       ((cwmbi(n,l,1),n=1,kbi),l=1,lmi)

  IF(kntbc.EQ.1)RETURN
  tboco=tbcex(kntbc)-tbcex(kntbc-1)
  rtboco=1./(tboco*3600.)

  DO n=1,kbi
     pdbi(n,2)=(pdbi(n,1)-pdbi0(n))*rtboco
     pdbi(n,1)=pdbi0(n)
  ENDDO

  DO l=1,lmi
     DO n=1,kbi
        tbi(n,l,2)=(tbi(n,l,1)-tbi0(n,l))*rtboco
        qbi(n,l,2)=(qbi(n,l,1)-qbi0(n,l))*rtboco
        ubi(n,l,2)=(ubi(n,l,1)-ubi0(n,l))*rtboco
        vbi(n,l,2)=(vbi(n,l,1)-vbi0(n,l))*rtboco
        q2bi(n,l,2)=(q2bi(n,l,1)-q2bi0(n,l))*rtboco
        cwmbi(n,l,2)=(cwmbi(n,l,1)-cwmbi0(n,l))*rtboco
        tbi(n,l,1)=tbi0(n,l)
        qbi(n,l,1)=qbi0(n,l)
        ubi(n,l,1)=ubi0(n,l)
        vbi(n,l,1)=vbi0(n,l)
        q2bi(n,l,1)=q2bi0(n,l)
        cwmbi(n,l,1)=cwmbi0(n,l)
     ENDDO
  ENDDO

  !***  for use in identifying where we are in the nest's
  !***  bc file, write out the hour at the start of each set

  WRITE(iutim1)tbcex(kntbc-1)

  WRITE(iutim1)pdbi
  WRITE(iutim1)tbi
  WRITE(iutim1)qbi
  WRITE(iutim1)ubi
  WRITE(iutim1)vbi
  WRITE(iutim1)q2bi
  WRITE(iutim1)cwmbi

  !***  since new bc conditions are needed for the last timestep
  !***  of the nest's forecast, duplicate the last set

  IF(tbcex(kntbc+1).LT.0.)THEN
     WRITE(iutim1)tbcex(kntbc)

     WRITE(iutim1)pdbi
     WRITE(iutim1)tbi
     WRITE(iutim1)qbi
     WRITE(iutim1)ubi
     WRITE(iutim1)vbi
     WRITE(iutim1)q2bi
     WRITE(iutim1)cwmbi
  ENDIF
END SUBROUTINE bcex
SUBROUTINE tll (almd,aphd,kb,tlmod,dtr,ctph0,stph0,tlm,tph)
  IMPLICIT NONE
  REAL :: anum,aph,cc,cph,crlm,ctph0,denom,dtr
  INTEGER :: k,kb
  REAL :: relm,sph,srlm,stph0,tlmod
  REAL almd(kb),aphd(kb),tlm(kb),tph(kb)
  DO k=1,kb
     relm=(almd(k)-tlmod)*dtr
     srlm=SIN(relm)
     crlm=COS(relm)
     aph=aphd(k)*dtr
     sph=SIN(aph)
     cph=COS(aph)
     cc=cph*crlm
     anum=cph*srlm
     denom=ctph0*cc+stph0*sph
     tlm(k)=-ATAN2(anum,denom)/dtr
     tph(k)=ASIN(ctph0*sph-stph0*cc)/dtr
  ENDDO
END SUBROUTINE tll
SUBROUTINE rotlle(u,v,vlat,vlon,tph0di,tlm0di,kbi,uout,vout)
  !
  ! subprogram: rotlle        rotate winds on lat/long grid to e-grid
  !   prgmmr: t.black         org: w/nmc22    date: ??-??-??
  !
  ! abstract: rotates winds on the lat/long grid to the eta model grid
  !
  ! program history log:
  !   ??-??-??  t.black
  !
  ! usage     call rotlle(u,v,vlat,vlon)
  !   input argument list:
  !     u        - lat/long u-component
  !     v        - lat/long v-component
  !     vlat     - latitude of e-grid v points (degrees)
  !     vlon     - longitude of e-grid v points (degrees)
  !
  !   output argument list:
  !     u        - eta grid u-component
  !     v        - eta grid v-component
  !
  !
  IMPLICIT NONE

  !*** rotate the lat-lon winds to the e-grid

  !*** note : input lat/long must be converted to radians !!!

  INCLUDE "parmeta.f90"
  REAL :: cc,cph,cphi0,crlm
  INTEGER :: k,kbi,l
  REAL :: orlam0,orphi0,rctph,relm,ru,rv,sph,sphi0,srlm
  REAL :: tlat,tlm0di,tlon,tph,tph0di
  INTEGER, PARAMETER :: lmax=lm
  REAL   , PARAMETER :: pi=3.141592654
  REAL   , PARAMETER :: d2r=pi/180.
  REAL :: u(kbi,lmax)
  REAL :: v(kbi,lmax)
  REAL :: vlat(kbi)
  REAL :: vlon(kbi)
  REAL :: uout(kbi,lmax)
  REAL :: vout(kbi,lmax)
  REAL :: cray(kbi)
  REAL :: dray(kbi)

  orphi0=tph0di*d2r
  orlam0=(360.+tlm0di)*d2r
  sphi0=SIN(orphi0)
  cphi0=COS(orphi0)

  DO k=1,kbi
     tlat=vlat(k)*d2r
     tlon=-vlon(k)*d2r
     relm=tlon-orlam0
     srlm=SIN(relm)
     crlm=COS(relm)
     sph=SIN(tlat)
     cph=COS(tlat)
     cc=cph*crlm
     tph=ASIN(cphi0*sph-sphi0*cc)
     rctph=1./COS(tph)
     cray(k)=sphi0*srlm*rctph
     dray(k)=(cphi0*cph+sphi0*sph*crlm)*rctph
  ENDDO

  DO l=1,lmax
     DO k=1,kbi
        ru=dray(k)*u(k,l)-cray(k)*v(k,l)
        rv=cray(k)*u(k,l)+dray(k)*v(k,l)
        uout(k,l)=ru
        vout(k,l)=rv
     ENDDO
  END DO
END SUBROUTINE rotlle
SUBROUTINE ptetae(etop,ionear,jonear, &
     imi,jmi,lmi,lmip,kbi, &
     pdbi,tbi,qbi,ubi,vbi,q2bi,cwmbi, &
     pdbo,tbo,qbo,ubo,vbo,q2bo,cwmbo,ntsd)

  !       subprogram documentaion block
  !
  !    subroutine:  ptetae              eta-to-eta interpolation program
  !    author: t black                  org: w/np22         date: 95-07-14
  !
  !
  !  abstract:  vertically interpolates eta data at one vertical 
  !             resolution to eta layers at another resolution.
  !             program assumes horizontal interpolation has been
  !             performed earlier.
  !
  !  reference:
  !  black, t. l., 1988: the step-mountain, eta coordinate model: a
  !          documentation. developmemt divison, national meteorological
  !          center.
  !
  !  program history log:
  !    95-07-14  t black - wrote code for eta restrt-->eta restrt ability
  !    96-10-29  t black - modified for general nesting capability
  !
  !  usage:
  !
  !    called from program intrst
  !
  !---------------------------------------------------------------------
  IMPLICIT NONE
  INCLUDE "parmeta.f90"
  INCLUDE "masks.f90"
  INCLUDE "srfdsp.f90"
  INTEGER, PARAMETER :: lmop=lm+1
  INTEGER, PARAMETER :: lmom=lm-1
  REAL   , PARAMETER :: cm1=2937.4
  REAL   , PARAMETER :: cm2=4.9283
  REAL   , PARAMETER :: cm3=23.5518
  REAL   , PARAMETER :: eps=0.622
  REAL   , PARAMETER :: rd=287.04
  REAL   , PARAMETER :: g=9.80
  REAL   , PARAMETER :: rg=1./g
  REAL   , PARAMETER :: epsq=2.e-12
  REAL   , PARAMETER :: epsq2=0.2
  REAL   , PARAMETER :: b1=11.87799326

  REAL :: alpdif,alpet,alpetk,alpt,apel,arglog,b,c,capa
  REAL :: cf,cloges,cp,cpog,cwmld,denom2,deta1,dpde,ese,esl,ess,etop
  REAL :: fi,gamma,gorg,h1,h2,h3,h3m1,h3m2
  INTEGER :: ibb,ibbp,ilb,ilbm,ilbp,ilbp2,imi,imt,imtp,in,itb
  INTEGER :: iter,jmi,jmt,jn,k,k1k,k2k,kbi,kbm,kntin,kntx
  INTEGER :: kp1,kstrt,l,l1,l2,l3,ld,lki,lko,llk,lmhi
  INTEGER :: lmi,lmip,mnknt,n,nk,nn,ntsd,num0,numk
  REAL :: p00,petak,phlb,pintoa,pintob,prf0,pt,q2dif,q2ij,q2ld,qetaij
  REAL :: qld,qsl,qsldm,qsmx,qsx,rgog,rhldm,rog,rtop,t0,udif
  REAL :: ukb,uld,vdif,vkb,vld

  !*** the inner grid 

  REAL   :: pdbi(kbi,2)
  REAL   :: tbi(kbi,lmi,2)
  REAL   :: qbi(kbi,lmi,2)
  REAL   :: ubi(kbi,lmi,2)
  REAL   :: vbi(kbi,lmi,2)
  REAL   :: q2bi(kbi,lmi,2)
  REAL   :: cwmbi(kbi,lmi,2)

  !*** the outer grid 

  REAL   :: pdbo(kbi)
  REAL   :: tbo(kbi,lm)
  REAL   :: qbo(kbi,lm)
  REAL   :: ubo(kbi,lm)
  REAL   :: vbo(kbi,lm)
  REAL   :: q2bo(kbi,lm)
  REAL   :: cwmbo(kbi,lm)

  INTEGER:: ionear(kbi)
  INTEGER:: jonear(kbi)
  REAL   :: etao(lmop)
  REAL   :: etai(lmip)
  REAL   :: detao(lm)
  REAL   :: aetao(lm)
  REAL   :: detai(lmi)
  REAL   :: aetai(lmi)
  REAL   :: pinto(kbi,lmop)
  REAL   :: alpr(kbi,lm)
  REAL   :: alp(kbi,lmop)
  REAL   :: alpsq(kbi,lmop)
  REAL   :: pdslo(kbi)
  REAL   :: zetao(kbi,lmop)
  REAL   :: zetai(kbi,lmip)
  INTEGER:: ldt1(kbi,lmi)
  REAL   :: dlt(kbi,lmi)
  REAL   :: pmido(kbi,lm)
  REAL   :: zsfci(kbi)
  REAL   :: refi(kbi)
  REAL   :: dflo(lmop)
  REAL   :: htmo(kbi,lm)
  INTEGER:: lmho(kbi)
  REAL   :: reso(kbi)
  REAL   :: psfco(kbi)
  REAL   :: tx(kbi)
  REAL   :: pdij(kbi)
  REAL   :: peta(kbi)
  INTEGER:: idx(kbi)
  INTEGER:: idxx(kbi)
  REAL   :: pdvp(kbi)
  REAL   :: alpeta(kbi)
  REAL   :: pmvp(kbi,2)
  REAL   :: pmidi(kbi)
  REAL   :: utop(kbi)
  REAL   :: vtop(kbi)
  REAL   :: q2top(kbi)
  REAL   :: alptu(kbi)
  REAL   :: phub(kbi)
  REAL   :: alpsfc(kbi)
  REAL   :: q2sfc(kbi)
  INTEGER:: k1(kbi)
  INTEGER:: k2(kbi)
  INTEGER:: ll(kbi)
  INTEGER :: kstrtx(kbi)

  INTEGER, PARAMETER :: iun1=24
  INTEGER, PARAMETER :: iun2=25

  !*** definition of constants needed for vertical interpolation

  cp=1004.6
  p00=1.e5
  gamma=0.0065
  prf0=101325.
  t0=288.
  gorg=g/(rd*gamma)
  rgog=1./gorg
  rog=rd/g
  cpog=cp/g
  capa=rd/cp
  pt=etop*100.
  alpt=alog(pt)

  REWIND iun1
  REWIND iun2
  READ(iun1)detao
  READ(iun2)detai
  DO k=1,kbi
     zsfci(k)=0.
     refi(k)=1.
  ENDDO

  !***  since we are on the nest's boundary, lmh always equals lmi

  lmhi=lmi


  !***  compute eta at the interfaces of the model layers

  !***  first on the outer grid

  etao(1)=0.

  DO l=2,lmop
     etao(l)=etao(l-1)+detao(l-1)
  ENDDO

  etao(lmop)=1.
  detao(lm)=1.-etao(lm)

  !***  then on the inner grid

  etai(1)=0.

  DO l=2,lmip
     etai(l)=etai(l-1)+detai(l-1)
  ENDDO

  etai(lmip)=1.
  detai(lmi)=1.-etai(lmi)

  !*** compute aeta in the middles of the layers

  DO l=1,lmi
     aetai(l)=0.5*(etai(l)+etai(l+1))
  ENDDO
  DO l=1,lm
     aetao(l)=0.5*(etao(l)+etao(l+1))
  ENDDO

  !***  compute the heights of the standard interfaces
  !***  for the outer grid

  DO l=1,lmop
     dflo(l)=g*t0*(1.-((pt+etao(l)*(prf0-pt))/prf0)**rgog)/gamma
  ENDDO

  !***  obtain htmo for the outer grid domain on the inner boundary.
  !***  simply use the height of the topography on the outer grid
  !***  at the point nearest to the inner grid point.

  DO l=1,lm
     DO k=1,kbi
        htmo(k,l)=htm(ionear(k),jonear(k),l)
     ENDDO
  ENDDO

  DO k=1,kbi
     lmho(k)=0
  ENDDO

  DO l=2,lm
     DO k=1,kbi
        IF(htmo(k,l).LT.0.1.AND.lmho(k).EQ.0)THEN
           lmho(k)=l-1
           reso(k)=1./etao(l)
        ENDIF
     END DO
  END DO

  DO k=1,kbi
     IF(lmho(k).EQ.0)THEN
        lmho(k)=lm
        reso(k)=1.
     ENDIF
  ENDDO

  !***  compute heights of the outer grid eta interfaces

  DO k=1,kbi
     pdslo(k)=reso(k)*pdbo(k)
     zetao(k,lmop)=0.0
     !       zetao(k,lmop)=fiso(k)*rg
  ENDDO

  DO l=lm,1,-1
     DO k=1,kbi
        dpde=detao(l)*pdslo(k)
        apel=pt+aetao(l)*pdslo(k)
        rtop=rd*tbo(k,l)*(1.+0.608*qbo(k,l))/apel
        fi=g*zetao(k,l+1)+rtop*dpde
        zetao(k,l)=((fi-dflo(l))*htmo(k,l)+dflo(l))*rg
     ENDDO
  ENDDO

  !***  find the pressure on the outer grid interfaces

  DO k=1,kbi
     pinto(k,1)=pt
  ENDDO

  DO l=2,lmop
     DO k=1,kbi
        pinto(k,l)=pdslo(k)*etao(l)+pt
     ENDDO
  ENDDO

  !***  find the pressure at the outer grid midlayer locations

  DO l=1,lm
     DO k=1,kbi
        pmido(k,l)=pinto(k,l)+0.5*detao(l)*pdslo(k)
     ENDDO
  ENDDO

  !***  compute the logs of pressure for upcoming interpolations

  DO ld=1,lm
     DO k=1,kbi
        alpr(k,ld)=alog(pinto(k,ld+1)/pinto(k,ld))
     ENDDO
  ENDDO

  DO ld=1,lmop
     DO k=1,kbi
        alp(k,ld)=alog(pinto(k,ld))
        alpsq(k,ld)=alp(k,ld)**2
     ENDDO
  ENDDO

  !***  vertical interpolation: quadratic interpolation of
  !***  heights (equivalent to temperature being linear in ln p)
  !***  and linear interpolation of winds


  !***  assign to pd on the nest boundary the value of pdsl at the
  !***  nearest h point in the parent eta run

  DO k=1,kbi
     pdij(k)=pdslo(k)

     !***  the following line is only true if pt is the same
     !***  for both the input and the output

     zetai(k,1)=zetao(k,1)
  ENDDO

  !***  computation of heights at upper boundaries of eta layers 
  !***  and temperature and spec humidity at middle of eta layers

  !***  ldt1 holds the index of the outer grid eta interface that is
  !***  directly above the inner grid eta interface in question.

  DO k=1,kbi
     peta(k)=pt
     ldt1(k,1)=1
  ENDDO

  mnknt=1
  DO l=2,lmi
     deta1=detai(l-1)

     numk=kbi
     kntx=mnknt
     kstrt=kntx+1

     DO k=1,kbi
        peta(k)=peta(k)+deta1*pdij(k)
        idxx(k)=k
        kstrtx(k)=1
     ENDDO

50   CONTINUE
     nk=numk
     DO nn=1,nk
        idx(nn)=idxx(nn)
     ENDDO
     kntin=kntx
     numk=0
     DO n=1,nk
        k=idx(n)
        petak=peta(k)
        pintoa=pinto(k,kntin)
        pintob=pinto(k,kntin+1)
        IF(petak.GT.pintoa.AND.petak.LE.pintob.OR. &
             petak.LT.pinto(k,2))THEN
           ldt1(k,l)=MIN(kntin,lm-1)
           IF(petak+detai(l)*pdij(k).LE.pintob)THEN
              kstrtx(k)=MIN(kstrt,kntin)
           ELSE
              kstrtx(k)=MIN(kstrt,kntin+1)
           ENDIF
        ELSE
           numk=numk+1
           idxx(numk)=k
           kntx=kntin+1
        ENDIF
     ENDDO

     IF(numk.GT.0)go to 50
     IF(numk.EQ.0)THEN
        mnknt=mnknt+1
        DO k=1,kbi
           mnknt=MIN(mnknt,kstrtx(k))
        ENDDO
     ENDIF
  END DO

  DO l=1,lmi
     DO k=1,kbi
        l1=ldt1(k,l)
        dlt(k,l)=alpr(k,l1+1)*(-alpr(k,l1))*alog(pinto(k,l1+2)/pinto(k,l1))
     ENDDO
  END DO

  !***  now compute the heights of the inner grid interfaces

  DO l=2,lmi
     DO k=1,kbi

        l1=ldt1(k,l)
        l2=l1+1
        l3=l1+2
        h1=zetao(k,l1)
        h2=zetao(k,l2)
        h3=zetao(k,l3)
        h3m2=h3-h2
        h3m1=h3-h1

        !***  the coefficients for quadratic interpolation of heights.

        b=(h3m2*(alpsq(k,l3)-alpsq(k,l1))-h3m1*(alpsq(k,l3)-alpsq(k,l2))) &
             /dlt(k,l)
        c=(h3m1*(alp(k,l3)-alp(k,l2))-h3m2*(alp(k,l3)-alp(k,l1))) &
             /dlt(k,l)

        arglog=pt+pdij(k)*etai(l)
        alpet=alog(arglog)

        !***  the heights of the inner grid eta interfaces.

        zetai(k,l)=h2+b*(alpet-alp(k,l2))+c*(alpet**2-alpsq(k,l2))

     END DO
  END DO

  DO k=1,kbi
     pdbi(k,1)=pdij(k)
  ENDDO

  !***  computation of winds within eta layers

  imt=2*imi-1
  jmt=jmi/2+1

  !***  arrangement of locations in the v boundary file

  !                       2nd
  !           ibbp--------->----------itb
  !           ilb                     kbi
  !            |                       |
  !       3rd /|\                     /|\ 4th
  !            |                       |
  !           imt                     ilbp
  !            1----------->----------ibb
  !                       1st
  ibb=imt/2
  ibbp=ibb+1
  itb=imt-1
  ilb=imt+jmt-2
  ilbp=ilb+1
  ilbp2=ilb+2
  ilbm=ilb-1
  imtp=imt+1
  kbm=kbi-1

  !***  find average pd over velocity boundary points

  DO k=1,ibb
     k1(k)=k
     k2(k)=k+1
  ENDDO
  DO k=ibbp,itb
     k1(k)=k+1
     k2(k)=k+2
  ENDDO
  DO k=imtp,ilbm
     k1(k)=k+1
     k2(k)=k+2
  ENDDO
  DO k=ilbp2,kbm
     k1(k)=k
     k2(k)=k+1
  ENDDO
  k1(imt)=1
  k2(imt)=imt+2
  k1(ilb)=ilb+1
  k2(ilb)=ibb+2
  k1(ilbp)=ibb+1
  k2(ilbp)=ilbp2
  k1(kbi)=kbi
  k2(kbi)=itb+2

  !***  average the pressure at velocity points.

  DO k=1,kbi
     k1k=k1(k)
     k2k=k2(k)
     pdvp(k)=0.5*(pdbi(k1k,1)+pdbi(k2k,1))
     pmvp(k,1)=0.5*(pmido(k1k,1)+pmido(k2k,1))
     alp(k,1)=alog(pmvp(k,1))
  END DO

  DO l=2,lm
     DO k=1,kbi
        k1k=k1(k)
        k2k=k2(k)
        pmvp(k,2)=0.5*(pmido(k1k,l)+pmido(k2k,l))
        alpr(k,l-1)=alog(pmvp(k,2)/pmvp(k,1))
        alp(k,l)=alog(pmvp(k,2))
        pmvp(k,1)=pmvp(k,2)
     ENDDO
  END DO

  !***  interpolate the winds

  mnknt=1
  DO l=1,lmi
     DO k=1,kbi
        alpeta(k)=alog(pt+pdvp(k)*aetai(l))
        idxx(k)=k
     ENDDO

     numk=kbi
     kntx=mnknt

230  CONTINUE
     nk=numk
     DO nn=1,nk
        idx(nn)=idxx(nn)
     ENDDO
     kntin=kntx
     kp1=MIN(lm,kntin+1)
     numk=0

     DO n=1,nk
        k=idx(n)
        alpetk=alpeta(k)
        IF(alpetk.GT.alp(k,kntin).AND.alpetk.LE.alp(k,kp1) &
             .AND.ABS(ubo(k,kp1)).GT.0..OR.alpetk.LT.alp(k,1))THEN
           uld=ubo(k,kp1)
           vld=vbo(k,kp1)
           cf=(alp(k,kp1)-alpetk)/alpr(k,kntin)
           ukb=uld+(ubo(k,kntin)-uld)*cf
           vkb=vld+(vbo(k,kntin)-vld)*cf
           utop(k)=ukb
           vtop(k)=vkb
           alptu(k)=alpetk
           ll(k)=l
           ubi(k,l,1)=ukb
           vbi(k,l,1)=vkb

           !***  if the new eta ground surface is below the old ground, then
           !***  assume the wind components decrease to zero at the surface.
           !***  the same is true for any new eta level below the lowest
           !***  old eta surface.

        ELSEIF(alpetk.GT.alp(k,lm).OR.ubo(k,kp1).EQ.0.)THEN
           udif=utop(k)
           vdif=vtop(k)
           alpdif=alptu(k)-alog(pt+pdvp(k))
           ll(k)=ll(k)+1
           llk=ll(k)
           ubi(k,llk,1)=utop(k)-udif*(alptu(k)-alpetk)/alpdif
           vbi(k,llk,1)=vtop(k)-vdif*(alptu(k)-alpetk)/alpdif
        ELSEIF(alpetk.GT.alp(k,kp1))THEN
           numk=numk+1
           idxx(numk)=k
           kntx=kntin+1
        ENDIF
     ENDDO
     IF(numk.GT.0)THEN
        IF(numk.EQ.kbi)mnknt=mnknt+1
        go to 230
     ENDIF
     !     if(numk.eq.0)mnknt=kntx

  END DO

  !***  interpolate q and cld water linearly in ln(p)

  DO l=lmi,1,-1
     DO k=1,kbi
        IF(l.EQ.lmi)phub(k)=0.
        phlb=phub(k)
        phub(k)=g*zetai(k,l)
        tbi(k,l,1)=-(phlb-phub(k))*(pt+aetai(l)*pdbi(k,1)) &
             /(rd*detai(l)*pdbi(k,1))
        IF(tbi(k,l,1).LT.150.)tbi(k,l,1)=150.
     ENDDO
  END DO

  DO k=1,kbi
     alp(k,1)=alog(pmido(k,1))
     psfco(k)=pdbo(k)+pt
  ENDDO

  DO ld=2,lm
     DO k=1,kbi
        alpr(k,ld-1)=alog(pmido(k,ld)/pmido(k,ld-1))
        alp(k,ld)=alog(pmido(k,ld))
     ENDDO
  ENDDO

  mnknt=1
  DO l=1,lmhi

     DO k=1,kbi
        pmidi(k)=aetai(l)*pdbi(k,1)+pt
        alpeta(k)=alog(pmidi(k))
        idxx(k)=k
     ENDDO

     numk=kbi
     kntx=mnknt

310  CONTINUE
     nk=numk
     DO nn=1,nk
        idx(nn)=idxx(nn)
     ENDDO
     kntin=kntx
     kp1=kntin+1
     num0=0
     numk=0

     DO n=1,nk
        k=idx(n)
        lko=lmho(k)
        alpetk=alpeta(k)
        IF(alpetk.GT.alp(k,kntin).AND.alpetk.LE.alp(k,kp1) &
             .AND.kp1.LE.lko.OR.alpetk.LT.alp(k,1) &
             .OR.alpetk.GT.alp(k,kp1).AND.pmidi(k).LE.psfco(k) &
             .AND.kp1.EQ.lko)THEN
           qld=qbo(k,kp1)
           cwmld=cwmbo(k,kp1)
           cf=(alp(k,kp1)-alpeta(k))/alpr(k,kntin)
           cwmbi(k,l,1)=cwmld+(cwmbo(k,kntin)-cwmld)*cf
           qetaij=qld+(qbo(k,kntin)-qld)*cf
           cloges=-cm1/tbi(k,l,1)-cm2*alog10(tbi(k,l,1))+cm3
           ese=10.**(cloges+2.)
           qsx=eps*ese/(pmidi(k)-ese*(1.-eps))
           qsmx=0.98*qsx
           qbi(k,l,1)=amin1(qetaij,qsmx,20.e-3)

           !***  if the new eta ground surface is below the old ground then
           !***  assume the rh (given t) is equal to that on the lowest
           !***  old surface and that the cloud water is zero.
           !***  note:  tbi inside the next block is actually the virtual
           !***         temperature yet it is being used as the true t

        ELSEIF(pmidi(k).GT.psfco(k))THEN
           cloges=-cm1/tbo(k,lko)-cm2*alog10(tbo(k,lko))+cm3
           ess=10.**(cloges+2.)
           qsldm=eps*ess/(pmido(k,lko)-ess*(1.-eps))
           rhldm=qbo(k,lko)/qsldm
           rhldm=amin1(rhldm,0.98)
           cloges=-cm1/tbi(k,l,1)-cm2*alog10(tbi(k,l,1))+cm3
           esl=10.**(cloges+2.)
           qsl=eps*esl/(aetai(l)*pdbi(k,1)+pt-esl*(1.-eps))
           qbi(k,l,1)=rhldm*qsl
           qbi(k,l,1)=amin1(qbi(k,l,1),20.e-3)
           cwmbi(k,l,1)=0.
        ELSEIF(alpetk.GT.alp(k,kp1))THEN
           numk=numk+1
           idxx(numk)=k
           kntx=kntin+1
        ENDIF
     ENDDO
     IF(numk.GT.0)THEN
        IF(numk.EQ.kbi)mnknt=mnknt+1
        go to 310
     ENDIF
     !     if(numk.eq.0)mnknt=kntx
  END DO

  DO l=1,lmi
     DO k=1,kbi
        qbi(k,l,1)=amax1(qbi(k,l,1),epsq)
        cwmbi(k,l,1)=amax1(cwmbi(k,l,1),0.)
     ENDDO
  ENDDO

  !***  convert from virtual to true temperature.

  DO l=1,lmi

     DO k=1,kbi
        tbi(k,l,1)=amin1(tbi(k,l,1),350.)
        tbi(k,l,1)=amax1(tbi(k,l,1),150.)
        tx(k)=tbi(k,l,1)
        idxx(k)=k
     ENDDO
     numk=kbi

     iter=0
360  CONTINUE
     iter=iter+1
     nk=numk
     numk=0
     DO nn=1,nk
        idx(nn)=idxx(nn)
     ENDDO

     DO n=1,nk
        k=idx(n)
        tbi(k,l,1)=tx(k)/(1.+0.608*qbi(k,l,1))
        cloges=-cm1/tbi(k,l,1)-cm2*alog10(tbi(k,l,1))+cm3
        ese=10.**(cloges+2.)
        denom2=aetai(l)*pdbi(k,1)+pt-ese*(1.-eps)
        qsx=eps*ese/denom2
        qsmx=0.98*qsx
        IF(qbi(k,l,1).GT.qsmx)THEN
           qbi(k,l,1)=qsmx
           numk=numk+1
           idxx(numk)=k
        ENDIF
     ENDDO
     IF(numk.GT.0.AND.iter.LE.2)go to 360

  END DO

  !***  interpolate q2 linearly in ln(p) (but on interfaces)

  DO k=1,kbi
     in=ionear(k)
     jn=jonear(k)
     q2sfc(k)=amax1(b1**(2./3.)*ustar(in,jn)*ustar(in,jn),epsq2)
     alpsfc(k)=alog(pdbi(k,1)*refi(k)+pt)
     alp(k,lmop)=alog(pinto(k,lmop))
  ENDDO

  DO ld=2,lm
     DO k=1,kbi
        alpr(k,ld)=alog(pinto(k,ld+1)/pinto(k,ld))
        alp(k,ld)=alog(pinto(k,ld))
     ENDDO
  END DO

  mnknt=2
  lki=lmhi-1
  DO l=1,lki

     DO k=1,kbi
        alpeta(k)=alog(etai(l+1)*pdbi(k,1)+pt)
        idxx(k)=k
     ENDDO

     numk=kbi
     kntx=mnknt

410  CONTINUE
     nk=numk
     DO nn=1,nk
        idx(nn)=idxx(nn)
     ENDDO

     kntin=kntx
     kp1=kntin+1
     numk=0


     DO n=1,nk
        k=idx(n)
        lko=lmho(k)
        alpetk=alpeta(k)
        IF(alpetk.GT.alp(k,kntin).AND.alpetk.LE.alp(k,kp1) &
             .AND.kp1.LE.lko.OR.alpetk.LT.alp(k,2))THEN
           q2ld=q2bo(k,kp1)
           cf=(alp(k,kp1)-alpetk)/alpr(k,kntin)
           q2ij=q2ld+(q2bo(k,kntin)-q2ld)*cf
           q2top(k)=amax1(q2ij,epsq2)
           alptu(k)=alpetk
           ll(k)=l
           q2bi(k,l,1)=q2top(k)

           !***  if the new q2 interface is below the top of the lowest 
           !***  outer grid layer, assume q2 decreases linearly to q2 at
           !***  the new ground.

        ELSEIF(alpetk.GT.alp(k,lko))THEN
           q2dif=q2top(k)-q2sfc(k)
           alpdif=alptu(k)-alpsfc(k)
           ll(k)=ll(k)+1
           llk=ll(k)
           q2bi(k,llk,1)=q2top(k)-q2dif*(alptu(k)-alpetk)/alpdif
        ELSEIF(alpetk.GT.alp(k,kp1))THEN
           numk=numk+1
           idxx(numk)=k
           kntx=kntin+1
        ENDIF
     ENDDO
     IF(numk.GT.0)THEN
        IF(numk.EQ.kbi)mnknt=mnknt+1
        go to 410
     ENDIF
     !     if(numk.eq.0)mnknt=kntx

  END DO

  !***  now redefine pd to have it equal to psfc-pt.

  DO k=1,kbi
     pdbi(k,1)=pdbi(k,1)*refi(k)
     q2bi(k,lmhi,1)=q2sfc(k)
  ENDDO
END SUBROUTINE ptetae
