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
  REAL :: tph0di,tph0do
  REAL :: pi,r2d,rsumh,rsumv,rtboco,sphi0,tboco,tlm0di,tlm0do
  REAL :: h1,h2,h3,h4,htot,v1,v2,v3,v4,vtot
  INTEGER :: i1,i2,i3,i4,j1,j2,j3,j4,ih,iv,jh,jv,imi,jmi,lmi,lmip
  INTEGER :: in,jn,k,kbi,kntbc,l,n,nwgts,iutim0,iutim1
 

  !*** filled E-grid

  INTEGER, PARAMETER :: imt=(im*2)-1, jmt=jm        
  INTEGER, PARAMETER :: imjm=(im*jm-jm/2)
  INTEGER :: mskfil(imt,jmt)
  INTEGER :: jt,it,i
  REAL:: tfil(imt,jmt)
  REAL:: qfil(imt,jmt)
  REAL:: ufil(imt,jmt)
  REAL:: vfil(imt,jmt)
  REAL:: q2fil(imt,jmt)
  REAL:: cwmfil(imt,jmt)

  REAL::  t1d(imjm)
  REAL::  q1d(imjm)
  REAL::  q21d(imjm)
  REAL::  cwm1d(imjm)
  REAL::  u1d(imjm)
  REAL::  v1d(imjm)


  !*** 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 :: aux(kbi,lm)
  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
     IF(mod(jouthb(k),2).ne.0) THEN
     iouthb(k)=(iouthb(k)+1)/2
     ionear(k)=iouthb(k)
     ELSE
     iouthb(k)=iouthb(k)/2
     ionear(k)=iouthb(k)
     ENDIF
     IF(mod(joutvb(k),2).ne.0) THEN
     ioutvb(k)=ioutvb(k)/2
     ELSE
     ioutvb(k)=(ioutvb(k)+1)/2
     ENDIF
  ENDDO
      
  DO k=1,kbi
     in=ionear(k)
     jn=jonear(k)
     pdbo(k)=pd(in,jn)
  ENDDO

  DO l=1,lm
       ! Filling E-grid

        call hij2k(t(0,0,l),t1d,l)
        call fillv(t1d,tfil,0,imt,jmt)

        call hij2k(q(0,0,l),q1d,l)
        call fillv(q1d,qfil,0,imt,jmt)

        call hij2k(q2(0,0,l),q21d,l)
        call fillv(q21d,q2fil,0,imt,jmt)

        call hij2k(cwm(0,0,l),cwm1d,l)
        call fillv(cwm1d,cwmfil,0,imt,jmt)

        call vij2k(u(0,0,l),u1d,l)
        call fillh(u1d,ufil,0,imt,jmt)

        call vij2k(v(0,0,l),v1d,l)
        call fillh(v1d,vfil,0,imt,jmt)

       ! create new topo mask
         DO jt=1,jmt  
         DO it=1,imt,2
            i=(it+1)/2
            IF (mod(jt,2).ne.0) then
               mskfil(it,jt)=htm(i,jt,l)
               mskfil(it+1,jt)=vtm(i,jt,l)
            ELSE
               mskfil(it,jt)=vtm(i,jt,l)
               mskfil(it+1,jt)=htm(i,jt,l)
            ENDIF
         ENDDO
         ENDDO

        if ((l.eq.30).and.(kntbc.eq.3)) then
            write(90)tfil
            write(90)q2fil
            write(90)cwmfil
            write(90)ufil
            write(90)vfil
            write(90)mskfil
        endif
        
     DO k=1,kbi

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

        in=iouthb(k)
        jn=jouthb(k)
 
        IF (mod(jn,2).ne.0) THEN
            ih=in*2-1
        ELSE
            ih=in*2
        ENDIF
        jh=jouthb(k)

        i1=ih
        j1=jh
        i2=ih+1
        j2=jh
        i3=ih        !chou   filled E-grid
        j3=jh+1      !chou   filled E-grid
        i4=ih+1      !chou   filled E-grid
        j4=jh+1

        h1=mskfil(i1,j1)
        h2=mskfil(i2,j2)
        h3=mskfil(i3,j3)
        h4=mskfil(i4,j4)

        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)*tfil(i1,j1)*h1 &
             +hwgts(k,2)*tfil(i2,j2)*h2 &
             +hwgts(k,3)*tfil(i3,j3)*h3 &
             +hwgts(k,4)*tfil(i4,j4)*h4)*rsumh

        qbo(k,l)=(hwgts(k,1)*qfil(i1,j1)*h1 &
             +hwgts(k,2)*qfil(i2,j2)*h2 &
             +hwgts(k,3)*qfil(i3,j3)*h3 &
             +hwgts(k,4)*qfil(i4,j4)*h4)*rsumh

        q2bo(k,l)=(hwgts(k,1)*q2fil(i1,j1)*h1 &
             +hwgts(k,2)*q2fil(i2,j2)*h2 &
             +hwgts(k,3)*q2fil(i3,j3)*h3 &
             +hwgts(k,4)*q2fil(i4,j4)*h4)*rsumh

        cwmbo(k,l)=(hwgts(k,1)*cwmfil(i1,j1)*h1 &
             +hwgts(k,2)*cwmfil(i2,j2)*h2 &
             +hwgts(k,3)*cwmfil(i3,j3)*h3 &
             +hwgts(k,4)*cwmfil(i4,j4)*h4)*rsumh

!         aux(k,l)=tbo(k,l)-t(in,jn,l)
!        if((aux(k,l).le.-1.0).or.(aux(k,l).ge.1.0))then
!           write(93,*)tbo(k,l),t(in,jn,l),tbo(k,l)-t(in,jn,l),k,l
!        endif

!         write(80,*)tbo(k,l),t(in,jn,l),tbo(k,l)-t(in,jn,l),k,l
        
!         tbo(k,l)=t(in,jn,l)

!         write(81,*)qbo(k,l),q(in,jn,l),qbo(k,l)-q(in,jn,l),k,l
!         qbo(k,l)=q(in,jn,l)

!         write(82,*)q2bo(k,l),q2(in,jn,l),q2bo(k,l)-q2(in,jn,l),k,l
!         q2bo(k,l)=q2(in,jn,l)

!         write(83,*)cwmbo(k,l),cwm(in,jn,l),cwmbo(k,l)-cwm(in,jn,l),k,l
!         cwmbo(k,l)=cwm(in,jn,l)

        !***  the v points

        in=ioutvb(k)
        jn=joutvb(k)

        IF (mod(jn,2).ne.0) THEN
            iv=in*2
        ELSE
            iv=in*2-1
        ENDIF
        jv=joutvb(k)

        i1=iv
        j1=jv
        i2=iv+1
        j2=jv
        i3=iv        !chou   filled E-grid
        j3=jv+1      !chou   filled E-grid
        i4=iv+1      !chou   filled E-grid
        j4=jv+1

        v1=mskfil(i1,j1)
        v2=mskfil(i2,j2)
        v3=mskfil(i3,j3)
        v4=mskfil(i4,j4)
        
        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)*ufil(i1,j1)*v1 &
             +vwgts(k,2)*ufil(i2,j2)*v2 &
             +vwgts(k,3)*ufil(i3,j3)*v3 &
             +vwgts(k,4)*ufil(i4,j4)*v4)*rsumv
        vbo(k,l)=(vwgts(k,1)*vfil(i1,j1)*v1 &
             +vwgts(k,2)*vfil(i2,j2)*v2 &
             +vwgts(k,3)*vfil(i3,j3)*v3 &
             +vwgts(k,4)*vfil(i4,j4)*v4)*rsumv


         write(84,*)ubo(k,l),u(in,jn,l),ubo(k,l)-u(in,jn,l),k,l,kntbc
!         ubo(k,l)=u(in,jn,l)

         write(85,*)vbo(k,l),v(in,jn,l),vbo(k,l)-v(in,jn,l),k,l,kntbc
!         vbo(k,l)=v(in,jn,l)

     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)then
    tboco= tbcex(2)-tbcex(1) 
    WRITE(iutim1)run,idat,ihrst,tboco
  RETURN
  ENDIF

  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)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)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
!chou
  do  l=1,lm
  READ(iun1,'(16x,F18.16)')detao(l)
  enddo
  do  l=1,lmi
  READ(iun2,'(16x,F18.16)')detai(l)
  enddo
!chou

  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
!     write(88,*) ' >> l:',l,' dflo:',dflo(l)
  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)
!     write(88,*) 'kbi,l:',k,l,'htmo,ionear,jonear:',htmo(k,l),ionear(k),jonear(k)
     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
       write(94,*)"lmho reso ",lmho(k),reso(k)
     END DO
  END DO

  DO k=1,kbi
     IF(lmho(k).EQ.0)THEN
        lmho(k)=lm
        reso(k)=1.
     ENDIF
      write(94,*)"if lmh k ",lmho(k),reso(k),k
  ENDDO

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

  DO k=1,kbi
!ROZANTE
     pdslo(k)=reso(k)*pdbo(k)
!ROZANTE
     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
!ROZANTE
     pdij(k)=pdslo(k)
!ROZANTE

     !***  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)
       write(94,*)"pdbi reso ",pdbi(k,1),refi(k),k
     q2bi(k,lmhi,1)=q2sfc(k)
  ENDDO
END SUBROUTINE ptetae
!
      SUBROUTINE FILLV(ARR1,ARR2,IFLAG,IDIM,JDIM) 
!$$$  SUBPROGRAM DOCUMENTATION BLOCK                                    
!                .      .    .                                          
! SUBPROGRAM:    FILLV       FILL V POINTS ON E-GRID H FIELD            
!   PRGRMMR: TREADON         ORG: W/NMC2     DATE: 92-12-23             
!                                                                       
! ABSTRACT:                                                             
!     THIS ROUTINE FILLS VELOCITY (V) POINTS ON A MASS                  
!     POINT E-GRID.  PASSED VARIABLE IFLAG TELLS THE ROUTINE            
!     WHAT TO DO WITH NEGATIVE VALUES OBTAINED DURING THE               
!     FILLING PROCESS                                                   
!        IFLAG=0, TAKE WHATEVER THE INTERPOLATION GIVES.                
!        IFLAG=1, TRUNCATE AT ZERO.                                     
!        IFLAG=2, TRUNCATE AT A VERY SMALL POSITIVE VALUE (EPS).        
!     THE INTERPOLATION CAN USE A BILINEAR SCHEME OR CUBIC              
!     INTERPOLATING POLYMONIALS.  PARAMETERS WGTCUB AND WGTLIN          
!     CONTROL THIS OPTION.  SETTING WGTCUB TO 1.0 ACTIVATES             
!     INTERPOLATION USING CUBIC INTERPOLATING POLYNOMIALS.              
!     SETTING WGTCUB TO 0.0 ACTIVATES BILINEAR INTERPOLATION.           
!                                                                       
!   .                                                                   
!                                                                       
! PROGRAM HISTORY LOG:                                                  
!   ??-??-??  ???                                                       
!   93-12-23  RUSS TREADON - ADDED COMMENTS AND VARIABLE                
!                            GRID DIMENSIONS TO                         
!                            SUBROUTINE CALL.                           
!                                                                       
! USAGE:    CALL FILLV(ARR1,ARR2,IFLAG,IDIM,JDIM)                       
!   INPUT ARGUMENT LIST:                                                
!     ARR2     - VELOCITY POINT DATA ON STAGGERED E-GRID.               
!     IFLAG    -                                                        
!     IDIM     - FIRST DIMENSION OF FILLED E-GRID.                      
!     JDIM     - SECOND DIMENSION OF FILLED E-GRID.                     
!                                                                       
!   OUTPUT ARGUMENT LIST:                                               
!     ARR1     - VELOCITY POINT DATA ON FILLED E-GRID.                  
!                                                                       
!   OUTPUT FILES:                                                       
!     NONE                                                              
!                                                                       
!   SUBPROGRAMS CALLED:                                                 
!     UTILITIES:                                                        
!       EFILL - FILL MISSING VALUES ON STAGGERED                        
!               E-GRID WITH FIELD MEAN.                                 
!       COEFI - COMPUTES ROW-WISE INTERPOLATION VALUES.                 
!       COEFJ - COMPUTES COLUMN-WISE INTERPOLATION VALUES.              
!     LIBRARY:                                                          
!       COMMON                                                          
!                                                                       
!   ATTRIBUTES:                                                         
!     LANGUAGE: FORTRAN                                                 
!     MACHINE : CRAY Y-MP                                               
!$$$                                                                    
!                                                                       
!                                                                       
!WAS      PARAMETER (IDIM=131,JDIM=131)                                 
!WAS      PARAMETER (KCYB=JDIM*(IDIM+1)/2-JDIM/2)                       
      PARAMETER (WGTCUB=0.0,WGTLIN=1.-WGTCUB,EPS=1.E-5) 
      PARAMETER (P5C=.125,P15C=3.375,P25C=15.625                        &
     &,          P35C=42.875,P45C=91.125)                               
      PARAMETER (P5S=.25,P15S=2.25,P25S=6.25                            &
     &     ,          P35S=12.25,P45S=20.25)                            
      DIMENSION ARR1(JDIM*(IDIM+1)/2-JDIM/2), ARR2(IDIM,JDIM) 
!WAS      DIMENSION ARR1(KCYB),ARR2(IDIM,JDIM)                          
!                                                                       
!********************************************************************** 
!     START FILLV HERE.                                                 
!                                                                       
!     REPLACE MISSING VALUES ON EGRID WITH FIELD MEAN                   
      IMJM = JDIM*(IDIM+1)/2-JDIM/2 
!hou      CALL EFILL(ARR1,IMJM)                                         
!                                                                       
      IDIMM=IDIM-1 
      IDIMM2=IDIM-2 
      IDIMM3=IDIM-3 
      IDIMM8=IDIM-8 
      IDIMM9=IDIM-9 
      JDIMM=JDIM-1 
      JDIMM2=JDIM-2 
      JDIMM3=JDIM-3 
      JDIMM8=JDIM-8 
      JDIMM9=JDIM-9 
      IADDP=IDIM/2+1 
      IADDM=IDIM/2 
      DO 3 J=1,JDIM 
      DO 3 I=1,IDIM 
      ARR2(I,J)=0. 
    3 CONTINUE 
!***                                                                    
!***  REASSIGN VALUES OF MODEL GRID TO IDIM X JDIM GRID ALONG BOTTOM    
!***  AND TOP ROWS.                                                     
!***                                                                    
      DO 5 I=1,IDIM,2 
      IND1B=(I+1)/2 
      IND1T=IND1B+IDIM*(JDIM/2) 
      ARR2(I,1)=ARR1(IND1B) 
      ARR2(I,JDIM)=ARR1(IND1T) 
    5 END DO 
!***                                                                    
!***  INTERPOLATE VALUES TO THOSE CORRESPONDING TO V POINTS NOW ON THE  
!***  IDIM X JDIM GRID FOR THE BOTTOM AND TOP ROWS.                     
!***                                                                    
      DO 10 I=2,IDIMM,2 
      IND1B=I/2 
      IND1T=IND1B+IDIM*(JDIM/2) 
      ARR2(I,1)=(ARR1(IND1B)+ARR1(IND1B+1))*0.5 
      ARR2(I,JDIM)=(ARR1(IND1T)+ARR1(IND1T+1))*0.5 
   10 END DO 
!***                                                                    
!***  REASSIGN VALUES FROM MODEL GRID TO IDIM X JDIM GRID ON ALL ROWS   
!***  OTHER THAN BOTTOM AND TOP THAT ARE BOUNDED BY H POINTS.           
!***                                                                    
   15 CONTINUE 
      DO 50 J=3,JDIMM2,2 
      IND1=((J-1)/2)*IDIM 
      DO 20 I=1,IDIM,2 
      IND1=IND1+1 
      ARR2(I,J)=ARR1(IND1) 
   20 END DO 
!***                                                                    
!***  INTERPOLATE TO CORRESPONDING V POINTS ON THESE ROWS.              
!***                                                                    
      IND1=((J-1)/2)*IDIM 
      DO 40 I=2,IDIMM,2 
      IND1=IND1+1 
      ARR2(I,J)=0.25*(ARR1(IND1)+ARR1(IND1+1)+ARR1(IND1+IADDP)+         &
     & ARR1(IND1-IADDM))                                                
   40 END DO 
   50 END DO 
!***                                                                    
!***  REASSIGN VALUES FROM MODEL GRID TO IDIM X JDIM GRID ON ALL ROWS   
!***  THAT ARE BOUNDED BY V POINTS.                                     
!***                                                                    
      DO 100 J=2,JDIMM,2 
      IND1=((J-2)/2)*IDIM+IADDP 
      DO 70 I=2,IDIMM,2 
      IND1=IND1+1 
      ARR2(I,J)=ARR1(IND1) 
   70 END DO 
!***                                                                    
!***  INTERPOLATE TO CORRESPONDING V POINTS ON THESE ROWS EXCEPT ON     
!***  THE SIDES.                                                        
!***                                                                    
      IND1=((J-2)/2)*IDIM+IADDP 
      DO 90 I=3,IDIMM2,2 
      IND1=IND1+1 
      ARR2(I,J)=0.25*(ARR1(IND1)+ARR1(IND1+1)+ARR1(IND1+IADDP)+         &
     & ARR1(IND1-IADDM))                                                
   90 END DO 
  100 END DO 
!***                                                                    
!***  INTERPOLATE TO CORRESPONDING V POINTS ON THE SIDES.               
!***                                                                    
      IND1L=IDIM/2+2-IDIM 
      DO 125 J=2,JDIMM,2 
      IND1L=IND1L+IDIM 
      IND1R=IND1L+IDIM/2-1 
      ARR2(1,J)=(ARR1(IND1L-IADDP)+ARR1(IND1L+IADDM))*0.5 
      ARR2(IDIM,J)=(ARR1(IND1R-IADDM)+ARR1(IND1R+IADDP))*0.5 
  125 END DO 
                                                                        
!hou removed cubic interpolation from here                              
                                                                        
  825 IF(IFLAG.EQ.1)THEN 
        DO 900 J=1,JDIM,2 
        DO 900 I=2,IDIMM,2 
        IF(ARR2(I,J).LT.0.)ARR2(I,J)=0. 
  900   CONTINUE 
        DO 925 J=2,JDIMM,2 
        DO 925 I=1,IDIM,2 
        IF(ARR2(I,J).LT.0.)ARR2(I,J)=0. 
  925   CONTINUE 
      ELSEIF(IFLAG.EQ.2)THEN 
        DO 950 J=1,JDIM,2 
        DO 950 I=2,IDIMM,2 
        IF(ARR2(I,J).LT.0.)ARR2(I,J)=EPS 
  950   CONTINUE 
        DO 975 J=2,JDIMM,2 
        DO 975 I=1,IDIM,2 
        IF(ARR2(I,J).LT.0.)ARR2(I,J)=EPS 
  975   CONTINUE 
      ENDIF 
!                                                                       
!     END OF ROUTINE.                                                   
!                                                                       
      RETURN 
      END      
!                                    
      SUBROUTINE FILLH(ARR1,ARR2,IFLAG,IDIM,JDIM) 
!$$$  SUBPROGRAM DOCUMENTATION BLOCK                                    
!                .      .    .                                          
! SUBPROGRAM:    FILLH       FILL H POINTS ON E-GRID V FIELD            
!   PRGRMMR: TREADON         ORG: W/NMC2     DATE: 92-12-23             
!                                                                       
! ABSTRACT:                                                             
!     THIS ROUTINE FILLS MASS (H) POINTS ON A VELOCITY (V)              
!     POINT E-GRID.                                                     
!   .                                                                   
!                                                                       
! PROGRAM HISTORY LOG:                                                  
!   ??-??-??  ???                                                       
!   93-12-23  RUSS TREADON - ADDED COMMENTS AND VARIABLE                
!                            GRID DIMENSIONS TO                         
!                            SUBROUTINE CALL.                           
!                                                                       
! USAGE:    CALL FILLH(ARR1,ARR2,IDIM,JDIM)                             
!   INPUT ARGUMENT LIST:                                                
!     ARR2     - VELOCITY POINT DATA ON STAGGERED E-GRID.               
!     IDIM     - FIRST DIMENSION OF FILLED E-GRID.                      
!     JDIM     - SECOND DIMENSION OF FILLED E-GRID.                     
!                                                                       
!   OUTPUT ARGUMENT LIST:                                               
!     ARR1     - VELOCITY POINT DATA ON FILLED E-GRID.                  
!                                                                       
!   OUTPUT FILES:                                                       
!     NONE                                                              
!                                                                       
!   SUBPROGRAMS CALLED:                                                 
!     UTILITIES:                                                        
!       EFILL - FILL MISSING VALUES ON STAGGERED                        
!               E-GRID WITH FIELD MEAN.                                 
!     LIBRARY:                                                          
!       COMMON                                                          
!                                                                       
!   ATTRIBUTES:                                                         
!     LANGUAGE: FORTRAN                                                 
!     MACHINE : CRAY Y-MP                                               
!$$$                                                                    
!                                                                       
!                                                                       
!WAS      PARAMETER (IDIM=131,JDIM=131)                                 
!WAS      PARAMETER (KCYB=JDIM*(IDIM+1)/2-JDIM/2)                       
!WAS      DIMENSION ARR1(KCYB),ARR2(IDIM,JDIM)                          
      DIMENSION ARR1(JDIM*(IDIM+1)/2-JDIM/2), ARR2(IDIM,JDIM) 
!                                                                       
!********************************************************************   
!     START FILLH HERE.                                                 
!                                                                       
!     FILL MISSING VALUES ON EGRID WITH FIELD MEAN                      
      IMJM = JDIM*(IDIM+1)/2-JDIM/2 
!hou      CALL EFILL(ARR1,IMJM)                                         
!                                                                       
      IDIMM=IDIM-1 
      IDIMM2=IDIM-2 
      IDIMM3=IDIM-3 
      JDIMM=JDIM-1 
      JDIMM2=JDIM-2 
      JDIMM3=JDIM-3 
      IADDP=IDIM/2+1 
      IADDM=IDIM/2 
      DO 3 J=1,JDIM 
      DO 3 I=1,IDIM 
      ARR2(I,J)=0. 
    3 CONTINUE 
!***                                                                    
!***  REASSIGN VALUES OF MODEL GRID TO V POINTS ON IDIM X JDIM GRID     
!***  ALONG BOTTOM AND TOP ROWS.                                        
!***                                                                    
      DO 5 I=2,IDIMM,2 
      IND1B=I/2 
      IND1T=IND1B+IDIM*(JDIM/2) 
      ARR2(I,1)=ARR1(IND1B) 
      ARR2(I,JDIM)=ARR1(IND1T) 
    5 END DO 
!***                                                                    
!***  INTERPOLATE VALUES TO THOSE CORRESPONDING TO H POINTS ON THE      
!***  IDIM X JDIM GRID FOR THE BOTTOM AND TOP ROWS.                     
!***                                                                    
      DO 10 I=3,IDIMM2,2 
      IND1B=(I-1)/2 
      IND1T=IND1B+IDIM*(JDIM/2) 
      ARR2(I,1)=(ARR1(IND1B)+ARR1(IND1B+1))*0.5 
      ARR2(I,JDIM)=(ARR1(IND1T)+ARR1(IND1T+1))*0.5 
   10 END DO 
!***                                                                    
!***  REASSIGN VALUES FROM MODEL GRID TO V POINTS ON THE                
!***  IDIM X JDIM GRID ON ALL ROWS OTHER THAN BOTTOM AND                
!***  TOP THAT ARE BOUNDED BY H POINTS.                                 
!***                                                                    
   15 CONTINUE 
      DO 50 J=3,JDIMM2,2 
      IND1=((J-1)/2)*IDIM 
      DO 20 I=2,IDIMM,2 
      IND1=IND1+1 
      ARR2(I,J)=ARR1(IND1) 
   20 END DO 
!***                                                                    
!***  INTERPOLATE TO CORRESPONDING H POINTS ON THESE ROWS               
!***  EXCEPT ON THE SIDES.                                              
!***                                                                    
      IND1=((J-1)/2)*IDIM 
      DO 40 I=3,IDIMM2,2 
      IND1=IND1+1 
      ARR2(I,J)=0.25*(ARR1(IND1)+ARR1(IND1+1)+ARR1(IND1+IADDP)+         &
     & ARR1(IND1-IADDM))                                                
   40 END DO 
   50 END DO 
!***                                                                    
!***  REASSIGN VALUES FROM MODEL GRID TO V POINTS ON IDIM X JDIM GRID   
!***  ON ALL ROWS THAT ARE BOUNDED BY V POINTS.                         
!***                                                                    
      DO 100 J=2,JDIMM,2 
      IND1=((J-2)/2)*IDIM+IADDM 
      DO 70 I=1,IDIM,2 
      IND1=IND1+1 
      ARR2(I,J)=ARR1(IND1) 
   70 END DO 
!***                                                                    
!***  INTERPOLATE TO CORRESPONDING H POINTS ON THESE ROWS.              
!***                                                                    
      IND1=((J-2)/2)*IDIM+IADDM 
      DO 90 I=2,IDIMM,2 
      IND1=IND1+1 
      ARR2(I,J)=0.25*(ARR1(IND1)+ARR1(IND1+1)+ARR1(IND1+IADDP)+         &
     & ARR1(IND1-IADDM))                                                
   90 END DO 
  100 END DO 
!***                                                                    
!***  INTERPOLATE TO CORRESPONDING H POINTS ON THE SIDES.               
!***                                                                    
      IND1L=1 
      DO 125 J=3,JDIMM2,2 
      IND1L=IND1L+IDIM 
      IND1R=IND1L+IDIM/2-1 
      ARR2(1,J)=(ARR1(IND1L-IADDP)+ARR1(IND1L+IADDM))*0.5 
      ARR2(IDIM,J)=(ARR1(IND1R-IADDM)+ARR1(IND1R+IADDP))*0.5 
  125 END DO 
!***                                                                    
!***  EXTRAPOLATE TO THE FOUR CORNER H POINTS.                          
!***                                                                    
      ARR2(1,1)=0.5*(1.5*(ARR2(2,1)+ARR2(1,2))                          &
     &              -0.5*(ARR2(4,1)+ARR2(1,4)))                         
      ARR2(IDIM,1)=0.5*(1.5*(ARR2(IDIMM,1)+ARR2(IDIM,2))                &
     &              -0.5*(ARR2(IDIMM3,1)+ARR2(IDIM,4)))                 
      ARR2(1,JDIM)=0.5*(1.5*(ARR2(2,JDIM)+ARR2(1,JDIMM))                &
     &              -0.5*(ARR2(4,JDIM)+ARR2(1,JDIMM3)))                 
      ARR2(IDIM,JDIM)=0.5*(1.5*(ARR2(IDIMM,JDIM)+ARR2(IDIM,JDIMM))      &
     &              -0.5*(ARR2(IDIMM3,JDIM)+ARR2(IDIM,JDIMM3)))         
      RETURN 
      END                                           


      SUBROUTINE HIJ2K(AIJ,AK,l)
!----------------------------------------------------
      INCLUDE "parmeta.f90"
      PARAMETER (IMJM=IM*JM-JM/2)
      REAL AK(IMJM)
      INTEGER :: l
      REAL, DIMENSION(0:im+1, 0:jm+1 ) :: AIJ
!----------------------------------------------------
      IMT=2*IM-1
      DO K=1,IMJM
        JX=(K-1)/IMT+1
        IX=K-(JX-1)*IMT
        IF(IX.LE.IM)THEN
          I=IX
          J=2*JX-1
        ELSE
          I=IX-IM
          J=2*JX
        ENDIF
        AK(K)=AIJ(I,J)
      ENDDO
!
      RETURN
      END

      SUBROUTINE VIJ2K(AIJ,AK,l)
!----------------------------------------------------
      INCLUDE "parmeta.f90"
      PARAMETER (IMJM=IM*JM-JM/2)
      REAL AK(IMJM)
      INTEGER :: l
      REAL, DIMENSION(0:im+1, 0:jm+1 ) :: AIJ
!----------------------------------------------------
      IMT=2*IM-1
      DO K=1,IMJM-1
        JX=(K-1)/IMT+1
        IX=K-(JX-1)*IMT
        IF(IX.LE.IM-1)THEN
          I=IX
          J=2*JX-1
        ELSE
          I=IX-IM+1
          J=2*JX
        ENDIF
        AK(K)=AIJ(I,J)
      ENDDO
!
      RETURN
      END
