      Program merge3d
 
C --- MERGE3D    Version: 1.04      Level: 090328                  MAIN
C
c --- Zhong-Xiang Wu
C     Developed: 1/12/2009

C --- Purpose:
C       Create new 3D.DAT from existing 3D.DAT files, which may have
C       different grid sizes. All 3D.DAT files were created from MM5
C       raw output.

c --- Updates:
c
c --- v1.03 L090320 to v1.04 L090328
c     1.   Fix a bug of initial n3dd variables
C
C     Zhong Wu
C     3/20/2009

c --- v1.02 L090316 to v1.03 L090320
c     1.   Fix a bug to check sigma levels when extracted sigma levels are
C          less than MM5 settings
C     2.   Add checks for extracted domain ranges
C
C     Zhong Wu
C     3/20/2009

c --- v1.01 L090312 to v1.02 L090316 (DGS)
c     1.   Allow spaces in file pathnames
c
c --- v1.0 L090112 to v1.01 L090312 (DGS)
c     1.   Augment/refine QA-related output information
c     2.   Fix logic to determine domain overlap (CHKRANGE)

C --- Rules for multiple source 3D.DAT files:
C       1. Map projection must be same
C       2. Grid size ratio must be integer numbers
C       3. New 3D.DAT data created from source 3D.DAT files 
C          in following order:
C          a. Use source data with same grid size
C          b. Bi-linear interpolate coarse source data
C          c. Use source data in the order specifed by user if more than one 
C             source data with same grid size exist.

C --- Input files are in the format of 3D.DAT Version 2.1
      include 'merge3d.par'
      include 'merge3d.cm1'

      parameter(iecho=1)

C ------------------------------------------------------------------------

c --- Set Code Version and Level number of code
      cver='1.04'
      clevel='090328'
c --- Open control file with user input information
      fct='merge3d.inp'

c     call comline(fct)
      numargs=IARGC()
      if(numargs.ge.1)then
         call GETARG(1,fct)
      endif

C --- Setup model info
      call getcontrol

C --- Get 3D configuration info
      call setup

C --- Check consistency between MM5 data
      if(nmm5.gt.1) call chkmap

C --- Convert lat/lon range to X/Y range if needed
      if(iselect.eq.2) call convert

C --- Check coverage
      call chkrange

C     Adjust X/Y beg/end to nearest finest grid location
      call adjust

C --- Check coverage again after adjustment
      call chkrange

C --- Mark source data
      call marksrc

C --- Finish setup, start to process data

C --- Read header info
      istart=1
      do k=1,nmm5
         fin=flnames(k)
         in=inb+k
         open(in,file=fin,status='old',action='read')
         if(k.eq.1) then
            iput=1
         else
            iput=0
         endif
         call rdhdv3(in,k,iput)
      enddo

      call getgeo
      istart=0

C --- Main Loop over All Hours  ----------------------------
C     Read one hour data
      ihrout=0
      idate0=idateb

c --- loop over all source datasets
 3000 continue

C     Read one hour data for all MM5 sets and all I/J grids
      do k=1,nmm5
         in=inb+k
         call rdv3(in,k)
      enddo

c --- Data heared; Surface variables
C     Order: psldot,rain,isnow,radswd,radlwd,t2,qq2,wd10,ws10,sst
      do ivar=1,n2d
         if(ivar.eq.3) then   !isnow
          call pass2di(ivar,iwrk)
          call getflag(iwrk,ix2tmp,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)
          call save2di(ivar,ix2tmp)
         else
          call pass2df(ivar,xwrk)

          if(ivar.eq.2) then   ! rain fall
           if(id_water.eq.1) then  
           call getrain1(xwrk,x2tmp,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)
           else
           call getrain2(xwrk,x2tmp,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)
           endif
          else
           call getdatf(xwrk,x2tmp,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)
          endif

          call save2df(ivar,x2tmp)
         endif
      enddo

C --- 3D profile
      do k=nz1,nz2
         ilvl=k
         do ivar=1,n3dd
          call pass3df(ivar,ilvl,xwrk)

          if(ivar.le.8) then
           call getdatf(xwrk,x2tmp,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)
          else   ! water contents
           if(id_water.eq.1) then
            call getrain1(xwrk,x2tmp,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)
           else
            call getrain2(xwrk,x2tmp,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)
           endif
          endif

          call save3df(ivar,ilvl,x2tmp)
         enddo
      enddo

C --- Output one hour 3D
      do j=1,nyn
         do i=1,nxn
c           Data header record
            pp=x2d(1,i,j)
            pre=x2d(2,i,j)
            isnow=nint(x2d(3,i,j))
            radsw=x2d(4,i,j)
            radlw=x2d(5,i,j)

            t2=x2d(6,i,j)
            qq2=x2d(7,i,j)
            wd10=x2d(8,i,j)
            ws10=x2d(9,i,j)
            sst=x2d(10,i,j)

            write(io,1030)idate,i,j,pp,pre,isnow,radsw,radlw
     &            ,t2,qq2,wd10,ws10,sst
 1030       format(i10.10,2i3,f7.1,f5.2,i2,3f8.1,f8.2,3f8.1)

            do iz=nz1,nz2

               do ii=1,n3dd
                  xcmp(ii)=x3d(ii,iz,i,j)
               enddo

C              Convert U/V to WD/WS
               uu=xcmp(4)
               vv=xcmp(5)
               call uv2ws(uu,vv,wdd,wss)
               xcmp(4)=wdd
               xcmp(5)=wss

C              Re-compute RH from q
               call getrh(xcmp,n3d)

C              Output one-level data
               call wrtcmp(xcmp,n3dd)

            enddo
         enddo
      enddo

C --- Finish one hour, goto next hour
      ihrout=ihrout+1
      print *,' Output:',ihrout,idate0

      if(idate0.lt.idatee) then
         call getnext(idate0,idaten,1)
         idate0=idaten
         goto 3000 
      endif

C --- Finish all hours
      write(ilg,*)' Total hours output: ',ihrout
      print *,' Total hours output: ',ihrout
      
      close(io)

      do k=1,nmm5
         in=inb+k
         close(in)
      enddo
      
      stop
      end

c      include 'merge3d.blk'
      include 'coordlib.for'

C ---------------------------------------------------------------------
      Subroutine Getcontrol

C --- MERGE3D    Version: 1.04       Level: 090312
C
C     Purpose:
C     Read control input file to setup
C      
c --- Updates:
c --- v1.0 L090112 to v1.01 L090312 (DGS)
c     1.   Augment/refine QA-related output information
c
C     Zhong-Xiang Wu
C     1/12/2009

      include 'merge3d.par'
      include 'merge3d.cm1'

c --- Read title line from control file
      open(ict,file=fct,status='old')

      read(ict,101)title
 101  format(a)
c
c --- Read output list file name 
      read(ict,101)flg
      call getstr(flg)
c
c --- Open output list file
      open(ilg,file=flg,status='unknown')
c
c --- Write version and level number of code to list file and
c --- echo back control file inputs
      write(ilg,102)cver,clevel
 102  format(//1x,'MERGE3D -- Version: ',a16,' Level: ',a16//)

      write(ilg,103)trim(title)
 103  format(1x,'Title: ',a)
      write(ilg,104)trim(flg)
 104  format(/1x,'List file: ',a)

      read(ict,101)fio
      call getstr(fio)
      write(ilg,105)trim(fio)
 105  format(/1x,'Output file: ',a)

      open(io,file=fio,status='unknown',action='write')

c --- Read beginning and ending dates/times (UTC) for 3D.DAT
c     (Format: yyyymmddhh)
      read(ict,*)idateb
      read(ict,*)idatee

      call gethours(ilg,idateb,idatee,ntotal)
      write(ilg,106)idateb,idatee,ntotal
 106  format(' Beg/End/Hours (UTC - YYYYMMDDHH):',3i14)

c --- Read/write method of selecting domain (1= X/Y; 2=Lon/Lat)
      read(ict,*)iselect
      write(ilg,107)iselect
 107  format(/1x,'Method for selecting domain: ',i6,
     1 ' (1=X/Y; 2=Lon/Lat)')

      if(iselect.eq.1) then
         read(ict,*)xbeg,xend
         read(ict,*)ybeg,yend
c         write(ilg,108)xbeg,xend,ybeg,yend
         write(ilg,*)' Beg/End X (km): ',xbeg,xend
         write(ilg,*)' Beg/End Y (km): ',ybeg,yend
      else if(iselect.eq.2)then
         read(ict,*)flonbeg,flonend
         read(ict,*)flatbeg,flatend
c         write(ilg,109)flonbeg,flonend,flatbeg,flatend
         write(ilg,*)' Beg/End Lon (deg): ',flonbeg,flonend
         write(ilg,*)' Beg/End Lat (deg): ',flatbeg,flatend
      else
         write(ilg,110)iselect
         write(*,110)iselect
         stop
      endif
      write(ilg,*)
 108  format(' Beg/End X/Y (km): ',2f14.3)
 109  format(' Beg/End Lon/Lat (deg): ',2f14.4)
 110  format(/1x,'Invalid selection method for defining domain: ',i4/
     1   1x,'Selection method must be: 1 (X/Y) or 2 (Lon/Lat)')

C --- MM5 Source Datasets used to create new 3D.DAT
      read(ict,*)nmm5

      if(nmm5.lt.1 .or.nmm5.gt.mxset) then
         write(ilg,*)' Error: Illegal source dataset specified: ',nmm5
         print *,' Error: Illegal source dataset specified: ',nmm5
         stop 
      endif

      do k=1,nmm5
         read(ict,101)fin
         call getstr(fin)
         flnames(k)=fin
         write(ilg,111)k,trim(flnames(k))
 111     format(1x,'Input file: ',i4,2x,a)
      enddo
      write(ilg,*)

C --- End of reading ict

      if(id_water.eq.1) then
         write(ilg,1113)id_water
      else
         write(ilg,1114)id_water
      endif
 1113 format(' Equal frequency Rain/Water content interpolation:',i3)
 1114 format(' Threshold Rain/Water content interpolation:',i3)
      write(ilg,*)

      close(ict)

      return
      end

C -------------------------------------------------------------------
      Subroutine convert

C --- MERGE3D    Version: 1.04       Level: 090312
C
C --- Purpose:
C       Convert lon/lat to X/Y 
C 
c --- Updates:
c --- v1.0 L090112 to v1.01 L090312 (DGS)
c     1.   Augment/refine QA-related output information
c     2.   TMSONE=1.0  (not used here so no impact)
c
C --- Zhong-Xiang Wu
C     1/12/2009

      include 'merge3d.par'
      include 'merge3d.cm1'
      
c --- Initialize GLOBE1
      cmapi='LL'
      iutmi=999.9
      tmsone=1.0
      rlat1i=config(3,1)  
      rlat2i=config(4,1)
      rlat0i=config(1,1)
      rlon0i=config(2,1)
      feasti = 0.
      fnorti = 0.
      
      cmapo=cmaps(1)
      iutmo=999.9
      tmsone=1.0
      rlat1=config(3,1)
      rlat2=config(4,1)
      rlat0=config(1,1)
      rlon0=config(2,1)
      feasto = 0.
      fnorto = 0.

c --  Set conversion vectors
      call GLOBE1(
     &  cmapi,iutmi,tmsone,rlat1i,rlat2i,rlat0i,rlon0i,feasti,fnorti,
     &  cmapo,iutmo,tmsone,rlat1,rlat2,rlat0,rlon0,feasto,fnorto,
     &  caction,vecti,vecto)

c --- Convert station coordinates to model map projection and datm,
c     Translate coordinates 
      datumi='WGS-84  '
      datumo='NWS-84  '

      call GLOBE(ilg,caction,datumi,vecti,datumo,vecto,
     &     flonbeg,flatbeg,xbeg,ybeg,iutmo,c4hem)
      call GLOBE(ilg,caction,datumi,vecti,datumo,vecto,
     &     flonend,flatend,xend,yend,iutmo,c4hem)

c      write(ilg,122)flonbeg,flonend,flatbeg,flatend,xbeg,xend,ybeg,yend
 122  format('Lon/Lat & X/Y beg-end:',8f12.4)
      write(ilg,*)'Extraction range converted from Lon/Lat to X/Y ---'
      write(ilg,*)' Beg/End X (km): ',xbeg,xend
      write(ilg,*)' Beg/End Y (km): ',ybeg,yend

      return
      end
 
C ---------------------------------------------------------------------
      Subroutine Chkmap

C --- MERGE3D    Version: 1.04       Level: 090112
C
C     Purpose:
C     Check map projections
C      
C     Zhong-Xiang Wu
C     1/12/2009

      include 'merge3d.par'
      include 'merge3d.cm1'

      k1=1
      cmap1=cmaps(k1)
      rlatc1=config(1,k1)
      rlonc1=config(2,k1)
      rstd11=config(3,k1)
      rstd21=config(4,k1)

      do k=2,nmm5
         cmap2=cmaps(k)
         rlatc2=config(1,k)
         rlonc2=config(2,k)
         rstd12=config(3,k)
         rstd22=config(4,k)

         if(trim(cmap1).ne.trim(cmap2) .or.
     &        rlatc1.ne.rlatc2 .or. rlonc1.ne.rlonc2 .or.   
     &        rstd11.ne.rstd12 .or. rstd21.ne.rstd22) then 
            write(ilg,120)
            write(ilg,121)k1,cmap1,rlatc1,rlonc1,rstd11,rstd21
            write(ilg,121)k,cmap2,rlatc2,rlonc2,rstd12,rstd22
            print 120
            print 121,k1,cmap1,rlatc1,rlonc1,rstd11,rstd21
            print 121,k,cmap2,rlatc2,rlonc2,rstd12,rstd22
 120        format(' Error: Map Projections do not Match:')
 121        format(' DataSet:',i4,1x,a,1x,4f12.4)
            stop
         endif

      enddo

      return
      end

C -------------------------------------------------------------------
      Subroutine adjust

C --- MERGE3D    Version: 1.04      Level: 090312
C
C --- Purpose:
C     Adjust X/Y beg/end to nearest finest grid location
c
c --- Updates:
c --- v1.0 L090112 to v1.01 L090312 (DGS)
c     1.   Augment/refine QA-related output information
c

c --- Zhong-Xiang Wu
C     Developed: 1/12/2009

      include 'merge3d.par'
      include 'merge3d.cm1'

      x1dmn=config(5,nmm5)  ! coarse SW grid Xo
      y1dmn=config(6,nmm5)  ! coarse SW grid Yo
      dxyn=config(7,1)      ! fine grid size used as new 3D grid size

      fi1=(xbeg-x1dmn)/dxyn
      fi2=(xend-x1dmn)/dxyn
      fj1=(ybeg-y1dmn)/dxyn
      fj2=(yend-y1dmn)/dxyn

      i1=int(fi1)
      i2=int(fi2)
      j1=int(fj1)
      j2=int(fj2)

      if(float(i2).lt.fi2) i2=i2+1
      if(float(j2).lt.fj2) j2=j2+1

      xbego=xbeg
      xendo=xend
      ybego=ybeg
      yendo=yend
      
      xbeg=x1dmn+i1*dxyn
      xend=x1dmn+i2*dxyn
      ybeg=y1dmn+j1*dxyn
      yend=y1dmn+j2*dxyn

      nxn=nint((xend-xbeg)/dxyn)+1
      nyn=nint((yend-ybeg)/dxyn)+1
      ibeg=1
      iend=nxn
      jbeg=1
      jend=nyn

      write(ilg,*)' Aligning range to model grid-points ---'
      write(ilg,123)xbego,xendo,ybego,yendo,xbeg,xend,ybeg,yend
      write(ilg,124)ibeg,iend,jbeg,jend,nxn,nyn,dxyn
      print *,' Aligning range to model grid-points ---'
      print 123,xbego,xendo,ybego,yendo,xbeg,xend,ybeg,yend
      print 124,ibeg,iend,jbeg,jend,nxn,nyn,dxyn
 123  format('  X/Y ranges adjusted from: ',4f12.3,/,
     &       '                        to: ',4f12.3)
 124  format('  New ibeg,iend,jbeg,jend,nxn,nyn,dxyn: ',
     &        6i5,f12.3)

C --- Check with mxnx,mxny ranges
      if(ibeg.lt.1 .or. iend.gt.mxnx) then
         write(ilg,201)
         write(ilg,202)ibeg,iend,nxn
         write(ilg,203)1,mxnx,mxnx
         write(*,201)
         write(*,202)ibeg,iend,nxn
         write(*,203)1,mxnx,mxnx
 201     format(' Extracted grid number outrange in X: - Reduce size')
 202     format(' Required  Beg/End/NX: ',3i5)
 203     format(' Range of Beg/End/MXNX: ',3i5)
         stop
      endif

      if(jbeg.lt.1 .or. jend.gt.mxny) then
         write(ilg,301)
         write(ilg,302)jbeg,jend,nyn
         write(ilg,303)1,mxny,mxny
         write(*,301)
         write(*,302)jbeg,jend,nyn
         write(*,303)1,mxny,mxny
 301     format(' Extracted grid number outrange in Y: - Reduce size')
 302     format(' Required  Beg/End/NX: ',3i5)
 303     format(' Range of Beg/End/MXNY: ',3i5)
         stop
      endif

      return
      end

C -------------------------------------------------------------------
      Subroutine chkrange

C --- MERGE3D    Version: 1.04      Level: 090312
C
C --- Purpose:
C       Check coverage of source data
c
c --- Updates:
c --- v1.0 L090112 to v1.01 L090312 (DGS)
c     1.   Augment/refine QA-related output information
c     2.   Fix logic to determine domain overlap
c

c --- Zhong-Xiang Wu
C     Developed: 1/12/2009

      include 'merge3d.par'
      include 'merge3d.cm1'

      write(ilg,*)'CHKRANGE:'
      write(ilg,*)' Output Beg/End X (km): ',xbeg,xend
      write(ilg,*)' Output Beg/End Y (km): ',ybeg,yend

      do k=nmm5,1,-1
         x1dmn=config(5,k)
         y1dmn=config(6,k)
         dxy=config(7,k)
         nx1=config(11,k) 
         ny1=config(12,k) 
         nx2=config(13,k) 
         ny2=config(14,k) 

         x1=(nx1-1)*dxy+x1dmn
         x2=(nx2-1)*dxy+x1dmn
         y1=(ny1-1)*dxy+y1dmn
         y2=(ny2-1)*dxy+y1dmn

         write(ilg,*)' File number ',k
         write(ilg,*)'   x range = ',x1,x2
         write(ilg,*)'   y range = ',y1,y2

         if(k.eq.nmm5) then
c ---       Outer (coarse grid) domain requirement
            if(xbeg.lt.x1 .or. xend.gt.x2 .or.
     &         ybeg.lt.y1 .or. yend.gt.y2) then
               write(ilg,101)xbeg,xend,ybeg,yend,x1,x2,y1,y2
               write(*,101)xbeg,xend,ybeg,yend,x1,x2,y1,y2
 101           format(' Error: Output grid ',4f12.3,/,
     &                '     is beyond data ',4f12.3)
               stop
            endif
         endif

c         ineed=0
c         if(xbeg.ge.x1 .and. xbeg.le.x2 .and.
c     &      ybeg.ge.y1 .and. ybeg.le.y2) ineed=1
c         if(xend.ge.x1 .and. xend.le.x2 .and.
c     &      yend.ge.y1 .and. yend.le.y2) ineed=1

         ineedx=0
         ineedy=0
c ---    Overlap of output x-range and input x-range?
         if(xend.GE.x1 .AND. x2.GE.xbeg) ineedx=1
c ---    Overlap of output y-range and input y-range?
         if(yend.GE.y1 .AND. y2.GE.ybeg) ineedy=1
c ---    Both must overlap
         ineed=ineedx*ineedy
         if(ineed.eq.0) then
            write(ilg,102)k,trim(flnames(k))
            write(*,102)k,trim(flnames(k))
 102        format(' No data in this file are in output grid:',i5,2x,a)
            write(*,*) ' Review list file for details'
            stop
         endif

      enddo

      return
      end

C ----------------------------------------------------------------
      Subroutine marksrc

C --- MERGE3D    Version: 1.04      Level: 090112
C
C --- Purpose:
C       Mark source data used in new 3D.DAT
C 
C --- Zhong-Xiang Wu
C     1/12/2009

      include 'merge3d.par'
      include 'merge3d.cm1'

C --- Search areas needed for new data
      do 1000 jj=1,nyn               ! new data J 
      yyn=ybeg+(jj-1)*dxyn

      do 1000 ii=1,nxn            ! new dat I
         xxn=xbeg+(ii-1)*dxyn

c        print *,jj,ii,xxn,yyn

C        Search source datasets starting from finest domain
         do k=1,nmm5            ! start from finest domain
            x1dmn=config(5,k)
            y1dmn=config(6,k)
            dxy=config(7,k)
            nx1=config(11,k)
            ny1=config(12,k)
            nx2=config(13,k)
            ny2=config(14,k)

            xb=x1dmn+(nx1-1)*dxy
            yb=y1dmn+(ny1-1)*dxy
            xe=x1dmn+(nx2-1)*dxy
            ye=y1dmn+(ny2-1)*dxy

c            print *,k,xb,xe,yb,ye

            if(xxn.ge.xb .and. xxn.le.xe .and.
     &         yyn.ge.yb .and. yyn.le.ye) then

               fi=(xxn-x1dmn)/dxy+1    ! source grid I  
               fj=(yyn-y1dmn)/dxy+1    ! source grid J

               i=int(fi)  
               j=int(fj)

               if(i.lt.nx1 .or. j.lt.ny1 .or.
     &            i.gt.nx2 .or. j.gt.ny2) then
                  write(ilg,125)i,j,nx1,nx1,ny1,ny2
                  print 125,i,j,nx1,nx1,ny1,ny2
 125              format(' Error: I/J out range - ',6i5,/,
     &                      ' Reset Domains')
               endif

               wx=1.0-(fi-i)
               wy=1.0-(fj-j)
               if(wx.lt.0.999 .or. wy.lt.0.999) then
                  itp=1
               else
                  itp=0
               endif

C ---          Check possible out range for needed interpolation 
               if(itp.eq.1) then
                  i1=i+1
                  j1=j+1
                  if(i1.lt.nx1 .or. j1.lt.ny1 .or.
     &               i1.gt.nx2 .or. j1.gt.ny2) then
                     write(ilg,126)i1,j1,nx1,nx1,ny1,ny2
                     print 126,i1,j1,nx1,nx1,ny1,ny2
 126                 format(' Error: I1/J1 out range - ',6i5,/,
     &                      ' Reset Domains')
                     stop
                  endif
               endif

C              Save flags
               idflag(1,ii,jj)=i    ! source SW I 
               idflag(2,ii,jj)=j    ! source SW I 
               idflag(3,ii,jj)=k    ! source data array index
               idflag(4,ii,jj)=itp  ! Interp flag
               wgt(1,ii,jj)=wx      ! I-weight
               wgt(2,ii,jj)=1.0-wx  ! I-weight
               wgt(3,ii,jj)=wy      ! J-weight
               wgt(4,ii,jj)=1.0-wy  ! J-weight

               goto 1000

            endif
         enddo

         write(ilg,201)ii,jj,xxn,yyn
         print 201,ii,jj,xxn,yyn
 201     format(' Error: No source data for new grid:',2i5,2f12.3)

         stop

 1000 continue

      return
      end

C ---------------------------------------------------------------------
      Subroutine getgeo

C --- MERGE3D    Version: 1.04      Level: 090112
C
C     Purpose:
C     Construct header geo info
C      
C     Zhong-Xiang Wu
C     1/12/2009

      include 'merge3d.par'
      include 'merge3d.cm1'

C --- Contruct geo data for new data
      do k=kbeg,kend
         siglvl(k)=siglvlm(k,1)

         if(nmm5.gt.1) then
            do m=2,nmm5
               if(siglvlm(k,m).ne.siglvl(k)) then
                  write(ilg,201)m,k,siglvlm(k,m),siglvl(k)
 201              format(' Error: Sigma levels differ:',2i4,2f12.4)
                  stop
               endif
            enddo
         endif
      enddo

      call getdati(ihgtm,ihgt,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)
      call getdati(ihgtcrsm,ihgtcrs,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)

      call getdatf(rlatm,rlat,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)
      call getdatf(rlonm,rlon,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)
      call getdatf(rlatcrsm,rlatcrs,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)
      call getdatf(rloncrsm,rloncrs,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)
      call getflag(ilandm,iland,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)

      call outgeo

      return
      end

C ----------------------------------------------------------------------
      subroutine getnext(idate,idaten,idt)

c --- MERGE3D    Version: 1.04      Level: 090112               
! --- Zhong-Xiang Wu           
!
! --- PURPOSE:
!     Get next timestamp 

      call getdate(idate,iy,im,id,ih)
      call chgtim(iy,im,id,ih,idt)
      call timestamp(iy,im,id,ih,idaten)

      return
      end

C --------------------------------------------------------------
      subroutine chgtim(iyr,imon,iday,ihour,idt)

c --- MERGE3D    Version: 1.04      Level: 090112               
! --- Zhong-Xiang Wu           
!
! --- PURPOSE:
!     Increase or decrease hours. 
!     Used to convert Local Stand Time to GMT or to setup time stamp for loop

      parameter(nmonth=12)

      dimension ndays(nmonth)

      data ndays/31,28,31,30,31,30,31,31,30,31,30,31/

      if(mod(iyr,4).eq.0) then
         ndays(2)=29 
      else
         ndays(2)=28
      endif

      ihour=ihour+idt

      if(idt.lt.0) goto 2000

 1000 if(ihour.gt.23) then
         ihour=ihour-24
         iday=iday+1
         
         if(iday.gt.ndays(imon)) then
            iday=1
            imon=imon+1
            
            if(imon.gt.12) then
               imon=1
               iyr=iyr+1
               if(iyr/4 .eq. iyr/4.0) then
                  ndays(2)=29 
               else
                  ndays(2)=28
               endif
            endif
         endif
         goto 1000
      else
         return
      endif

 2000 continue

 3000 if(ihour.lt.0) then
         ihour=ihour+24
         iday=iday-1
         
         if(iday.le.0) then
            imon=imon-1
            if(imon.le.0) then
               iyr=iyr-1
               imon=12
               if(iyr/4 .eq. iyr/4.0) then
                  ndays(2)=29 
               else
                  ndays(2)=28
               endif
            endif
   
            iday=ndays(imon)
         endif
         goto 3000
      else
         return
      endif

      end

c --------------------------------------------------------------
	subroutine gethours(ilog,ndate1,ndate2,nhrs)

c --- MM5DECODE    Version: 2.5       Level: 061228               
C Purpose:
c       Get hours between two dates
C Zhong Wu
C (3/22/1999)

C Changes
C  1. Fix an error for total hours when beginning year is 
C     leap year
C  Zhong Wu      
C  5/23/2003

c	Get hours between to dates

	if(ndate2.lt.ndate1) then
	   print *,'Error in gethours: ndate2 < ndate1'
	   print *,ndate1,ndate2
	   stop
	endif

	call getdate(ndate1,iyr1,imon1,iday1,ihour1)
	call julday(iyr1,imon1,iday1,jday1)

	call getdate(ndate2,iyr2,imon2,iday2,ihour2)
	call julday(iyr2,imon2,iday2,jday2)

        call deltt(iyr1,jday1,ihour1,iyr2,jday2,ihour2,nhrs)

        nhrs=nhrs+1   ! deltt is one hour less total hours
 
	return
	end

c ----------------------------------------------------------------------
      subroutine rdcmp(in,xcmp,n3dd)  

c --- MERGE3D    Version: 1.04      Level: 090112               WRTCMP
C
C     Purpose:
C     Read one level 3D profile in compressed form
C     (Compressed output of MM5.dat data record for ioutform >= 93)
C
C     Zhong-Xiang Wu
C     1/12/2009

      parameter (xzero=0.0049)

      include 'merge3d.par'

      dimension xcmp(n3d)
      character*200 buff

      read(in,101)buff
 101  format(a)

      xcmp=0
      read(buff,201)(xcmp(ii),ii=1,9)
 201  format( f4.0,f6.0,f6.1,f4.0,f5.1,f6.2,f3.0,f5.2,5f6.3)
c 201  format( f4.0,f6.0,f6.1,f4.0,f5.1,f6.2,f3.0,f5.2,5f5.2)

      if(xcmp(9).lt.0) then
         nn=-nint(xcmp(9))
         do i1=9,9+nn-1
            xcmp(i1)=0
         enddo
      else
         read(buff,201)(xcmp(i1),i1=1,n3dd)
      endif

c     Convert WS/WD to U/V
      wdd=xcmp(4)
      wss=xcmp(5)

      call ws2uv(wdd,wss,uu,vv)

      xcmp(4)=uu
      xcmp(5)=vv

      return
      end

c ----------------------------------------------------------------------
      subroutine wrtcmp(xcmp,n3dd)  

c --- MERGE3D    Version: 1.04      Level: 090112               WRTCMP
C
C     Purpose:
C     Write one level 3D profile in compressed form
C     (Compressed output of MM5.dat data record for ioutform >= 93)
C
C     Zhong-Xiang Wu
C     1/12/2009

      parameter (xzero=0.0049)

      include 'merge3d.par'

      dimension xcmp(n3d)

      ipre=nint(xcmp(1))
      izh=nint(xcmp(2))
      atk=xcmp(3)
      iwd=nint(xcmp(4))
      aws=xcmp(5)
      awz=xcmp(6)
      irh=nint(xcmp(7))
      aq=xcmp(8)

      if(iwd.eq.0 .and. aws.gt.0.01) iwd=360 

C --- Compress
      idcmp=1
      ncmp=0
      do i=9,n3dd
         if(xcmp(i).ge.xzero) then
            idcmp=0
            goto 1000
         else
            ncmp=ncmp+1
         endif
      enddo

 1000 continue

      if(idcmp.eq.1) then
         write(io,95)ipre,izh,atk,iwd,aws,awz,irh,aq,float(-ncmp)
      else
         write(io,95)ipre,izh,atk,iwd,aws,awz,irh,aq,
     &       (xcmp(i),i=9,n3dd)
      endif
 95   format( i4,i6,f6.1,i4,f5.1,f6.2,i3,f5.2,5f6.3)

      return
      end

C -----------------------------------------------------------------
      Subroutine save3df(ivar,ilvl,x2tmp)

C --- MERGE3D    Version 1.0   Level 090112                 PASS2D
C
C Purpose:
C     Save one real variable in 2D array to its master array x2d
C
C Zhong Wu
C 1/12/2009

      include 'merge3d.par'
      include 'merge3d.cm1'

      do j=1,nyn
         do i=1,nxn
            x3d(ivar,ilvl,i,j)=x2tmp(i,j)
         enddo
      enddo

      return
      end

C -----------------------------------------------------------------
      Subroutine save2df(ivar,x2tmp)

C --- MERGE3D    Version 1.0   Level 090112                 PASS2D
C
C Purpose:
C     Save one real variable in 2D array to its master array x2d
C
C Zhong Wu
C 1/12/2009

      include 'merge3d.par'
      include 'merge3d.cm1'

      do j=1,nyn
         do i=1,nxn
            x2d(ivar,i,j)=x2tmp(i,j)
         enddo
      enddo

      return
      end

C -----------------------------------------------------------------
      Subroutine save2di(ivar,ix2tmp)

C --- MERGE3D    Version 1.0   Level 090112                 PASS2D
C
C Purpose:
C     Save one integer variable in 2D array to its master array x2d
C
C Zhong Wu
C 1/12/2009

      include 'merge3d.par'
      include 'merge3d.cm1'

      do j=1,nyn
         do i=1,nxn
            x2d(ivar,i,j)=float(ix2tmp(i,j))
         enddo
      enddo

      return
      end

C -----------------------------------------------------------------
      Subroutine pass3df(ivar,ilvl,xdata)

C --- MERGE3D    Version 1.0   Level 090112                 PASS2D
C
C Purpose:
C     Pass one real variable in the 3D data header master array to a
C     3D array
C
C Zhong Wu
C 1/12/2009

      include 'merge3d.par'
      include 'merge3d.cm1'
      dimension xdata(mxnx,mxny,mxset)

      do k=1,nmm5
         nxb=config(11,k)
         nyb=config(12,k)
         nxe=config(13,k)
         nye=config(14,k)

         do j=nyb,nye
            do i=nxb,nxe
               xdata(i,j,k)=x3dm(ivar,ilvl,i,j,k)
            enddo
         enddo

      enddo

      return
      end

C -----------------------------------------------------------------
      Subroutine pass2df(ivar,xdata)

C --- MERGE3D    Version 1.0   Level 090112                 PASS2D
C
C Purpose:
C     Pass one real variable in the 2D data header master array 
C     to a 3D array
C
C Zhong Wu
C 1/12/2009

      include 'merge3d.par'
      include 'merge3d.cm1'
      dimension xdata(mxnx,mxny,mxset)

      do k=1,nmm5
         nxb=config(11,k)
         nyb=config(12,k)
         nxe=config(13,k)
         nye=config(14,k)

         do j=nyb,nye
            do i=nxb,nxe
               xdata(i,j,k)=x2dm(ivar,i,j,k)
            enddo
         enddo

      enddo

      return
      end

C -----------------------------------------------------------------
      Subroutine pass2di(ivar,ixdata)

C --- MERGE3D    Version 1.0   Level 090112                 PASS2D
C
C Purpose:
C     Pass one integer variable in the 2D data header master array
C     to a 3D array
C
C Zhong Wu
C 1/12/2009

      include 'merge3d.par'
      include 'merge3d.cm1'
      dimension ixdata(mxnx,mxny,mxset)

      do k=1,nmm5
         nxb=config(11,k)
         nyb=config(12,k)
         nxe=config(13,k)
         nye=config(14,k)

         do j=nyb,nye
            do i=nxb,nxe
               ixdata(i,j,k)=nint(x2dm(ivar,i,j,k))
            enddo
         enddo

      enddo

      return
      end

c ----------------------------------------------------------------
      Subroutine rdhdv3(in,is,iput)

C --- MERGE3D    Version: 1.04      Level: 090328               
C
C Purpose:
C     Read 3D.DAT V2.0 file header
C
C Zhong-Xiang Wu
C 09/11/2007

c --- v1.03 L090320 to v1.04 L090328
c     1.   Fix a bug of initial n3dd variables
C
C     Zhong Wu
C     3/20/2009

      include 'merge3d.par'
      include 'merge3d.cm1'

c     Records #1
      read(in,1010)dataset,dataver,datamod
 1010 format(2a16,a64)
c      print *,' dataset,dataver,datamod: ',dataset,dataver,datamod 

      read(in,1011)ncomm
 1011 format(i4)

      if(iput.eq.1) then 
         write(io,1010)dataset,dataver,datamod
         write(io,1011)ncomm+2
         write(io,1112)cver,clevel
         if(id_water.eq.1) then
            write(io,1113)id_water
         else
            write(io,1114)id_water
         endif
      endif
 1112 format('Merged by MERGE3D version: ',a,'  Level: ',a)
 1113 format('Equal frequency Rain/Water content interpolation:',i3)
 1114 format('Threshold Rain/Water content interpolation:',i3)

      do i=1,ncomm
         read(in,1111)comment
         if(iput.eq.1) then 
c            write(io,1111)trim(comment)
            write(io,1111)comment
         endif
      enddo
c 1111 format(a)
 1111 format(a80)
c 1111 format(a132)

c --- header record #3 (mm5.dat only)
      read(in,1012)ioutw,ioutq,ioutc,iouti,ioutg,iosrf
 1012 format (6i3)
      if(ioutc.eq.0) then
         print *,' Error: Rain/cloud water must be included:',ioutc
         stop
      endif

      n3dd=n3d
      if(ioutg.eq.0) n3dd=n3d-1
      if(ioutg.eq.0 .and. iouti.eq.0) n3dd=n3d-3

      if(iput.eq.1) 
     &     write(io,1012)ioutw,ioutq,ioutc,iouti,ioutg,iosrf

c     Record #4 - projection parameters
      read(in,1013)cmap,rlatc,rlonc,truelat1,truelat2,
     &        x1dmn,y1dmn,dxy,nx,ny,nz
 1013 format(a4,f9.4,f10.4,2f7.2,2f10.3,f8.3,2i4,i3)

      if(iput.eq.1) 
     &     write(io,1013)cmap,rlatc,rlonc,truelat1,truelat2,
     &         xbeg,ybeg,dxyn,nxn,nyn,nz

c     Record #5 - physical options
      read(in,1014)inhyd,imphys,icupa,ibltyp,ifrad,isoil,
     &	    ifddaan,ifddaob,igrdt,ipbl,ishf,ilhf,iustr,iswdn,
     &	    ilwdn,ist1,ist2,ist3,ist4,ist5,ist6

      if (iput.eq.1)
     &	   write(io,1014)inhyd,imphys,icupa,ibltyp,ifrad,isoil,
     &           ifddaan,ifddaob,igrdt,ipbl,ishf,ilhf,iustr,iswdn,
     &           ilwdn,ist1,ist2,ist3,ist4,ist5,ist6
 1014 format (30(i3))
	
c     Record #6 - Extracted domain Stamp 1
      read(in,1015)idateout,nhours,nxsub,nysub,nzsub

      if(iput.eq.1)
     &    write(io,1015)idateb,ntotal,nxn,nyn,nzsub
 1015 format(i10,i5,3i4)

C --- Following records will change, output them later
c     Record #7: - Extracted domain
      read(in,1016)nx1,ny1,nx2,ny2,nz1,nz2,
     &         rxmin,rxmax,rymin,rymax
 1016 format(6i4,2f10.4,2f9.4)

      if(nz1.lt.1 .or. nz2.gt.mxnz) then
         print *,'Error: NZ1/NZ2 outrange: ',nz1,nz2,mxnz
         stop
      endif

      kbeg=nz1
      kend=nz2

      do k=nz1,nz2
         read(in,1022)siglvlm(k,is)
      enddo
 1022 format(f6.3)	       

      do j=ny1,ny2
         do i=nx1,nx2
            read(in,1023)ii,jj,rlatm(i,j,is),rlonm(i,j,is)
     &           ,ihgtm(i,j,is),ilandm(i,j,is),rlatcrsm(i,j,is)
     &           ,rloncrsm(i,j,is),ihgtcrsm(i,j,is)
            if(i.ne.ii .or. j.ne.jj) then
               print *,'I/J match error:',ii,jj,i,j
               stop
            endif
         enddo
      enddo
 1023 format(2i4,f9.4,f10.4,i5,i3,1x,f9.4,f10.4,i5)

      return
      end

C ------------------------------------------------------------------
      subroutine rdv3(in,idm)

C --- MERGE3D    Version: 1.04      Level: 090112               
C
C Purpose:
C     Read 3D.DAT V2.0 data records
C
C Zhong-Xiang Wu
C 09/11/2007

      include 'merge3d.par'
      include 'merge3d.cm1'

      nxb=config(11,idm)
      nyb=config(12,idm)
      nxe=config(13,idm)
      nye=config(14,idm)
      nz1=config(15,idm)
      nz2=config(16,idm)

 1000 continue
      do j=nyb,nye
         do i=nxb,nxe
            read(in,101)idate,ii,jj,
     &           (x2dm(ivar,ii,jj,idm),ivar=1,n2d)
 101        format(i10.10,2i3,f7.1,f5.2,f2.0,3f8.1,f8.2,3f8.1)

            if(i.ne.ii .or. j.ne.jj) then
               print *,'I/J match error:',ii,jj,i,j
               stop
            endif

            do iz=nz1,nz2
               call rdcmp(in,xcmp,n3dd)
               do i1=1,n3dd       
                  x3dm(i1,iz,i,j,idm)=xcmp(i1) ! U/V in x3d data now
               enddo
            enddo

         enddo
      enddo

      if(idate.lt.idate0) then
         goto 1000
      elseif(idate.gt.idate0) then
         write(ilg,102)idate,idate0
         print 102,idate,idate0
 102     format(' Required Date not exist:',2i14)
         stop
      endif   

      return
      end

c --------------------------------------------------------------------
      subroutine getdate(ndate,iyr,imon,iday,ihour)

C --- MERGE3D    Version: 1.04      Level: 090112               
C
C Purpose:
c     Get year,month,day and hour
C
C Zhong-Xiang Wu
C 09/11/2007

      include 'merge3d.par'

      iyr=ndate/isix
      imon=ndate/ifour-iyr*100
      iday=ndate/100-iyr*ifour-imon*100
      ihour=ndate-iyr*isix-imon*ifour-iday*100
      
c      print *,ndate
c      print *,iyr,imon,iday,ihour

      return
      end


c --------------------------------------------------------------------
      subroutine timestamp(iyr,imon,iday,ihour,idatestmp) 

      include 'merge3d.par'

C --- MERGE3D    Version: 1.04      Level: 090112               
C
C Purpose:
c     Get time stamp
C
C Zhong-Xiang Wu
C 09/11/2007

      idatestmp=ihour+iday*100+imon*ifour+iyr*isix

      return
      end

c----------------------------------------------------------------------
      subroutine julday(iyr,imo,iday,ijuldy)
c
c --- calmm5   Version: 1.04      level: 901130                  julday
c ---          j. scire, src
c
c --- purpose:  compute the julian day number from the gregorian
c               date (month, day)
c
c --- inputs:
c           iyr - integer      - year
c           imo - integer      - month
c          iday - integer      - day
c
c --- output:
c          ijul - integer      - julian day
c
c --- julday called by:  main
c --- julday calls:      none
c----------------------------------------------------------------------
      integer kday(12)
      data kday/0,31,59,90,120,151,181,212,243,273,304,334/

      ijuldy=kday(imo)+iday
      if(imo.le.2)return
      if(mod(iyr,4).eq.0)ijuldy=ijuldy+1

      return
      end

c----------------------------------------------------------------------
      subroutine grday(iyr,ijul,imo,iday)
c----------------------------------------------------------------------
c
c --- CALMET   Version: 5.1       Level: 901130                   GRDAY
c ---          J. Scire, SRC
c
c --- Compute the Gregorian date (month, day) from the Julian day
c
c --- INPUTS:
c           IYR - integer      - Year
c          IJUL - integer      - Julian day
c        Parameters: IO6
c
c --- OUTPUT:
c           IMO - integer      - Month
c          IDAY - integer      - Day
c
c --- GRDAY called by:  COMP
c --- GRDAY calls:      none
c----------------------------------------------------------------------
c
c --- include parameters
c
      integer kday(12,2)
      data kday/31,59,90,120,151,181,212,243,273,304,334,365,
     1          31,60,91,121,152,182,213,244,274,305,335,366/
c
      ileap=1
      if(mod(iyr,4).eq.0)ileap=2
      if(ijul.lt.1.or.ijul.gt.kday(12,ileap))go to 11
c
      do 10 i=1,12
      if(ijul.gt.kday(i,ileap))go to 10
      imo=i
      iday=ijul
      if(imo.ne.1)iday=ijul-kday(imo-1,ileap)
      return
10    continue
c
11    continue
      print 12,iyr,ijul
12    format(//2x,'ERROR in SUBR. GRDAY -- invalid Julian day '//2x,
     1 'iyr = ',i5,3x,'ijul = ',i5)
      stop
      end

C --------------------------------------------------------------------	
      subroutine ws2uv(wd,ws,uu,vv)

C --- MERGE3D    Version: 1.04      Level: 090112               
C
C Purpose:
c     Convert wind direction/speed to U/V
C
C Zhong-Xiang Wu
C 09/11/2007

      include 'merge3d.par'

C     Convert WD/WS to U/V
      wdrad = deg2rad*wd
      uu = -ws * sin(wdrad)
      vv = -ws * cos(wdrad)
      
      return
      end

C --------------------------------------------------------------------	

      subroutine uv2ws(uu,vv,wd,ws)

C --- MERGE3D    Version: 1.04      Level: 090112               
C
C Purpose:
c     Convert U/V to wind direction/speed
C
C Zhong-Xiang Wu
C 09/11/2007

C     Convert U/V to WD/WS
      ws=sqrt(uu**2+vv**2)
      if(ws.gt.1.e-9) then
         xwd=270.-57.295778*atan2(vv,uu)
         wd=amod(xwd,360.)
      else
         wd=0.
      endif

      return
      end

C ----------------------------------------------------------------------
      Subroutine getstr(cstr)

C --- MERGE3D    Version 1.0   Level 090112                 GETSTR
C
C Purpose:
C    Extract string content up to the '!' delimiter
c
c --- Updates:
c --- V1.0 Level 090112 to V1.02 Level 090316 (DGS)
c     - Allow spaces in strings by accepting all characters up to the
c       delimiter
C
C Zhong Wu
C 1/12/2009

      character*120 cstr

c      nt=index(cstr,' ')-1
      nt=index(cstr,'!')-1
      cstr=TRIM(cstr(1:nt))
      
      return
      end

C -------------------------------------------------------------------
      subroutine setup

C --- MERGE3D    Version 1.0   Level 090312                  SETUP
C
C Purpose:
C    Extract 3D configuration from headers 
c
c --- Updates:
c --- v1.0 L090112 to v1.01 L090312 (DGS)
c     1.   Augment/refine QA-related output information
c
C
C Zhong Wu
C 1/12/2009

      parameter(iecho=1)
      include 'merge3d.par'
      include 'merge3d.cm1'

      character*200 buff

      if(iecho.eq.1) then
c ---    Write variable names
         write(ilg,'(a)')'Map   NLat0    ELon0    NLat1  NLat2 '//
     &        '   X1dmn     Y1dmn     Dxy   nx  ny  nz nx1 ny1'//
     &        ' nx2 ny2 nz1 nz2   RXmin     RXmax    RYmin    '//
     &        'RYmax  Id_reg'
      endif

C --- Read header info
      istart=1
      do 1000 k=1,nmm5

      fin=flnames(k)
      in=inb
      open(in,file=fin,status='old',action='read')

C     Read Header Record #1
      read(in,201)dataset,dataver,datamod
 201  format(2a16,a64)

C     Read Header Record #2
      read(in,*)ncomm

C     Read Header Record #3  ! read/write original comment lines
      do i=1,ncomm
        read(in,202)comment
      enddo
 202  format(a)

C     Read Header Record #4
      read(in,203)ioutw,ioutq,ioutc,iouti,ioutg,iosrf
 203  format (30(i3))

C     Read Header Record #5
      read(in,204)cmap,rlatch,rlonch,truelat1h,truelat2h,
     &        x1dmn,y1dmn,dxy,nx,ny,nz
 204  format(a4,f9.4,f10.4,2f7.2,2f10.3,f8.3,2i4,i3)

C     Read Header Record #6
      read(in,203)inhyd,imphys,icupa,ibltyp,ifrad,isoil,
     &         ifddaan,ifddaob,igrdt,ipbl,ishf,ilhf,iustr,iswdn,
     &         ilwdn,ist1,ist2,ist3,ist4,ist5,ist6

C     Read Header Record #7
      read(in,205)idatebmm5,nhours,nxsub,nysub,nz
 205  format(i10.10,i5,3i4)

C     Read Header Record #8
      read(in,206)nx1,ny1,nx2,ny2,nz1,nz2,rxmin,rxmax,rymin,rymax
 206  format(6i4,2f10.4,2f9.4)

      id_reg=k
      write(buff,301)cmap,rlatch,rlonch,truelat1h,truelat2h,
     &        x1dmn,y1dmn,dxy,nx,ny,nz,nx1,ny1,nx2,ny2,nz1,nz2,
     &        rxmin,rxmax,rymin,rymax,id_reg
 301  format(a4,f9.4,f10.4,2f7.2,2f10.3,f8.3,2i4,i3,6i4
     &       ,2f10.4,2f9.4,i4)
      read(buff(1:4),'(a4)')cmaps(k)
      read(buff(5:),*)(config(i,k),i=1,21)

      if(iecho.eq.1) then

         read(buff,301)cmap,rlatch,rlonch,truelat1h,truelat2h,
     &        x1dmn,y1dmn,dxy,nx,ny,nz,nx1,ny1,nx2,ny2,nz1,nz2,
     &        rxmin,rxmax,rymin,rymax,id_reg
         write(ilg,301)cmap,rlatch,rlonch,truelat1h,truelat2h,
     &        x1dmn,y1dmn,dxy,nx,ny,nz,nx1,ny1,nx2,ny2,nz1,nz2,
     &        rxmin,rxmax,rymin,rymax,id_reg
c         print 302,cmap,rlatch,rlonch,truelat1h,truelat2h,
c     &        x1dmn,y1dmn,dxy,nx,ny,nz,nx1,ny1,nx2,ny2,nz1,nz2,
c     &        rxmin,rxmax,rymin,rymax,id_reg
c 302     format(1x,a4,f9.4,f10.4,2f7.2,2f10.3,f8.3,2i4,i3,6i4
c     &        ,2f10.4,2f9.4,i4)

      endif

      close(in)

 1000 continue
      write(ilg,*)
      
      return
      end

C ---------------------------------------------------------------------
      subroutine getdatf(xs,xd,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)

C --- MERGE3D    Version 1.0   Level 090112
C
C Purpose:
C     Get real data for new 3D.DAT from source datasets
C
C Zhong Wu
C 1/12/2009

      real xs(mxnx,mxny,mxset),xd(mxnx,mxny),wgt(4,mxnx,mxny)
      integer idflag(4,mxnx,mxny)

      do jj=1,nyn       !JJ in new data
         do ii=1,nxn    !II in new data

            i=idflag(1,ii,jj)    ! source SW I 
            j=idflag(2,ii,jj)    ! source SW I 
            k=idflag(3,ii,jj)    ! source data array index
            itp=idflag(4,ii,jj)  ! Interp flag (1,0)

            if(itp.eq.0) then
               xd(ii,jj)=xs(i,j,k)
            elseif(itp.eq.1) then
               i1=i+1
               j1=j+1

               a1=xs(i,j,k)
               a2=xs(i1,j,k)
               b1=xs(i,j1,k)
               b2=xs(i1,j1,k)

               wx1=wgt(1,ii,jj) ! I-weight
               wx2=wgt(2,ii,jj) ! I-weight
               wy1=wgt(3,ii,jj) ! J-weight
               wy2=wgt(4,ii,jj) ! J-weight

               aa=a1*wx1+a2*wx2
               bb=b1*wx1+b2*wx2
               cc=aa*wy1+bb*wy2
               xd(ii,jj)=cc
            else
               print *,' Error: itp flag must be 0/1 -',itp
               stop
            endif
         enddo
      enddo

      return
      end

C ---------------------------------------------------------------------
      subroutine getdati(ixs,ixd,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)

C --- MERGE3D    Version 1.0   Level 090112
C
C Purpose:
C     Get integer data for new 3D.DAT from source datasets
C
C Zhong Wu
C 1/12/2009

      integer ixs(mxnx,mxny,mxset),ixd(mxnx,mxny),idflag(4,mxnx,mxny)
      real wgt(4,mxnx,mxny)

      do jj=1,nyn       !JJ in new data
         do ii=1,nxn    !II in new data

            i=idflag(1,ii,jj)    ! source SW I 
            j=idflag(2,ii,jj)    ! source SW I 
            k=idflag(3,ii,jj)    ! source data array index
            itp=idflag(4,ii,jj)  ! Interp flag (1,0)

            if(itp.eq.0) then
               ixd(ii,jj)=ixs(i,j,k)
            elseif(itp.eq.1) then
               i1=i+1
               j1=j+1

               a1=ixs(i,j,k)
               a2=ixs(i1,j,k)
               b1=ixs(i,j1,k)
               b2=ixs(i1,j1,k)

               wx1=wgt(1,ii,jj) ! I-weight
               wx2=wgt(2,ii,jj) ! I-weight
               wy1=wgt(3,ii,jj) ! J-weight
               wy2=wgt(4,ii,jj) ! J-weight

               aa=a1*wx1+a2*wx2
               bb=b1*wx1+b2*wx2
               cc=aa*wy1+bb*wy2
               ixd(ii,jj)=nint(cc)
            else
               print *,' Error: itp flag must be 0/1 -',itp
               stop
            endif
         enddo
      enddo

      return
      end

C ---------------------------------------------------------------------
      subroutine getflag(ixs,ixd,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)

C --- MERGE3D    Version 1.0   Level 090112
C
C Purpose:
C     Get flag type data (landuse(1-25), snow cover (0/1) integer 
C     data for new 3D.DAT from source datasets
C
C Zhong Wu
C 1/12/2009

      integer ixs(mxnx,mxny,mxset),ixd(mxnx,mxny),idflag(4,mxnx,mxny)
      real wgt(4,mxnx,mxny)

      do jj=1,nyn       !JJ in new data
         do ii=1,nxn    !II in new data
            i=idflag(1,ii,jj)    ! source SW I 
            j=idflag(2,ii,jj)    ! source SW I 
            k=idflag(3,ii,jj)    ! source data array index
            itp=idflag(4,ii,jj)  ! Interp flag (1,0)

            if(itp.eq.0) then
               ixd(ii,jj)=ixs(i,j,k)
            elseif(itp.eq.1) then
               i1=i+1
               j1=j+1

               a1=ixs(i,j,k)
               a2=ixs(i1,j,k)
               b1=ixs(i,j1,k)
               b2=ixs(i1,j1,k)

               wx1=wgt(1,ii,jj) ! I-weight
               wx2=wgt(2,ii,jj) ! I-weight
               wy1=wgt(3,ii,jj) ! J-weight
               wy2=wgt(4,ii,jj) ! J-weight

               if(wx1.ge.wx2) then
                  wx1=1
                  wx2=0
               else
                  wx1=0
                  wx2=1
               endif

               if(wy1.ge.wy2) then
                  wy1=1
                  wy2=0
               else
                  wy1=0
                  wy2=1
               endif

               aa=a1*wx1+a2*wx2
               bb=b1*wx1+b2*wx2
               cc=aa*wy1+bb*wy2
               ixd(ii,jj)=nint(cc)
            else
               print *,' Error: itp flag must be 0/1 -',itp
               stop
            endif
         enddo
      enddo

      return
      end

C ---------------------------------------------------------------------
      subroutine getrain1(xs,xd,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)

C --- MERGE3D    Version 1.0   Level 090112
C
C Purpose:
C     Get rail fall based on frequency conservation
C     1. Determine Yes(1)/No(0) flag for rain
C     2. Interpolate rainfall if rain flag=1. Otherwise
C        Set rainfall=0
C
C Zhong Wu
C 1/12/2009

      real xs(mxnx,mxny,mxset),xd(mxnx,mxny),wgt(4,mxnx,mxny)
      integer idflag(4,mxnx,mxny)

      do jj=1,nyn       !JJ in new data
         do ii=1,nxn    !II in new data

            i=idflag(1,ii,jj)    ! source SW I 
            j=idflag(2,ii,jj)    ! source SW I 
            k=idflag(3,ii,jj)    ! source data array index
            itp=idflag(4,ii,jj)  ! Interp flag (1,0)

            if(itp.eq.0) then
               xd(ii,jj)=xs(i,j,k)
            elseif(itp.eq.1) then
               i1=i+1
               j1=j+1

               a1=xs(i,j,k)
               a2=xs(i1,j,k)
               b1=xs(i,j1,k)
               b2=xs(i1,j1,k)

               ia1=0
               ia2=0
               ib1=0
               ib2=0

               if(a1.gt.0) ia1=1
               if(a2.gt.0) ia2=1
               if(b1.gt.0) ib1=1
               if(b2.gt.0) ib2=1

               wx1=wgt(1,ii,jj) ! I-weight
               wx2=wgt(2,ii,jj) ! I-weight
               wy1=wgt(3,ii,jj) ! J-weight
               wy2=wgt(4,ii,jj) ! J-weight

               raa=ia1*wx1+ia2*wx2
               rbb=ib1*wx1+ib2*wx2
               rcc=raa*wy1+rbb*wy2
               icc=nint(rcc)

               if(icc.ge.1) then
                  aa=a1*wx1+a2*wx2
                  bb=b1*wx1+b2*wx2
                  cc=aa*wy1+bb*wy2
                  xd(ii,jj)=cc
               else
                  xd(ii,jj)=0.
               endif
            else
               print *,' Error: itp flag must be 0/1 -',itp
               stop
            endif
         enddo
      enddo

      return
      end

C ---------------------------------------------------------------------
      subroutine getrain2(xs,xd,idflag,wgt,mxnx,mxny,mxset,nxn,nyn)

C --- MERGE3D    Version 1.0   Level 090112
C
C Purpose:
C     Get rail fall based on threshold value
C     1. Interpolate rainfall
C     2. If rainfall < rainthld, set rainfall=0
C
C Zhong Wu
C 1/12/2009

      real xs(mxnx,mxny,mxset),xd(mxnx,mxny),wgt(4,mxnx,mxny)
      integer idflag(4,mxnx,mxny)

      parameter(rainthld=0.01)

      do jj=1,nyn       !JJ in new data
         do ii=1,nxn    !II in new data

            i=idflag(1,ii,jj)    ! source SW I 
            j=idflag(2,ii,jj)    ! source SW I 
            k=idflag(3,ii,jj)    ! source data array index
            itp=idflag(4,ii,jj)  ! Interp flag (1,0)

            if(itp.eq.0) then
               xd(ii,jj)=xs(i,j,k)
            elseif(itp.eq.1) then
               i1=i+1
               j1=j+1

               a1=xs(i,j,k)
               a2=xs(i1,j,k)
               b1=xs(i,j1,k)
               b2=xs(i1,j1,k)

               wx1=wgt(1,ii,jj) ! I-weight
               wx2=wgt(2,ii,jj) ! I-weight
               wy1=wgt(3,ii,jj) ! J-weight
               wy2=wgt(4,ii,jj) ! J-weight

               aa=a1*wx1+a2*wx2
               bb=b1*wx1+b2*wx2
               cc=aa*wy1+bb*wy2

               if(cc.lt.rainthld) cc=0.0

               xd(ii,jj)=cc
            else
               print *,' Error: itp flag must be 0/1 -',itp
               stop
            endif
         enddo
      enddo

      return
      end

C ------------------------------------------------------------------
      Subroutine outgeo

c --- MERGE3D    Version: 1.04      Level: 090112                     

C --- Purpose:
C     Output geo part in 3D.DAT header
C
C Zhong-Xiang Wu
C 11/28/2008

      include 'merge3d.par'
      include 'merge3d.cm1'

C     Min/max lat/lon for selected domain
      rxmin=99999
      rxmax=-99999
      rymin=99999
      rymax=-99999

      do j=jbeg,jend
         do i=ibeg,iend
            alat=rlat(i,j)
            alon=rlon(i,j)
            if(alat.lt.rymin) rymin=alat
            if(alat.gt.rymax) rymax=alat
            if(alon.lt.rxmin) rxmin=alon
            if(alon.gt.rxmax) rxmax=alon
         enddo
      enddo

C --- Record #7
c      write(io,1015)idateb,ntotal,nxn,nyn,nz
c 1015 format(i10,i5,3i4)

C --- Record #8
      write(io,1016)ibeg,jbeg,iend,jend,kbeg,kend,rxmin
     &      ,rxmax,rymin,rymax
 1016 format(6i4,2f10.4,2f9.4)

C     Sigma level
      write(io,1022)(siglvl(k),k=kbeg,kend)
 1022 format(f6.3)

C     Lat/lon
      do j=jbeg,jend
         do i=ibeg,iend
            alat=rlat(i,j)
            alon=rlon(i,j)
            alatcrs=rlatcrs(i,j)
            aloncrs=rloncrs(i,j)
            ih=ihgt(i,j)
            land=iland(i,j)
            ihcrs=ihgtcrs(i,j)

            write(io,1023)i,j,alat,alon,ih,land,alatcrs
     &           ,aloncrs,ihcrs
1023        format(2i4,f9.4,f10.4,i5,i3,1x,f9.4,f10.4,i5)
         enddo
      enddo

      return
      end

c------------------------------------------------------------------------------
      subroutine deltt(j1yr,j1jul,j1hr,j2yr,j2jul,j2hr,jleng)
c------------------------------------------------------------------------------
c
c --- CALUTILS   Version: 2.56     Level: 941215                  DELTT
c ---            J. Scire, SRC
c
c --- Compute the difference (in hours) between two dates & times
c ---    (time #2 - time #1)
c
c --- INPUTS:
c              J1YR - integer    - Year of date/time #1
c             J1JUL - integer    - Julian day of date/time #1
c              J1HR - integer    - Hour of date/time #1
c              J2YR - integer    - Year of date/time #2
c             J2JUL - integer    - Julian day of date/time #2
c              J2HR - integer    - Hour of date/time #2
c
c --- OUTPUT:
c             JLENG - integer    - Difference (#2 - #1) in hours
c
c --- DELTT called by:  host subroutines
c --- DELTT calls:      none
c------------------------------------------------------------------------------
c
      jmin=min0(j1yr,j2yr)
c
c --- find the number of hours between Jan. 1 of the "base" year and
c --- the first date/hour
      if(j1yr.eq.jmin)then
         j1=0
      else
         j1=0
         j1yrm1=j1yr-1
         do 10 i=jmin,j1yrm1
         if(mod(i,4).eq.0)then
            j1=j1+8784
         else
            j1=j1+8760
         endif
10       continue
      endif
      j1=j1+(j1jul-1)*24+j1hr
c
c --- find the number of hours between Jan. 1 of the "base" year and
c --- the second date/hour
      if(j2yr.eq.jmin)then
         j2=0
      else
         j2=0
         j2yrm1=j2yr-1
         do 20 i=jmin,j2yrm1
         if(mod(i,4).eq.0)then
            j2=j2+8784
         else
            j2=j2+8760
         endif
20       continue
      endif
      j2=j2+(j2jul-1)*24+j2hr
c
c --- compute the time difference (in hours)
      jleng=j2-j1
c
      return
      end

C ----------------------------------------------------------------------
      Subroutine getrh(xcmp,n3d)

c --- MERGE3D    Version: 1.04      Level: 090112                     

C --- Purpose:
C     Replace interpolated RH with re-computed RH from q
C
C Zhong-Xiang Wu
C 1/12/2009

      dimension xcmp(n3d)

C Oder in xcmp
C     'PRESS','GH','TK','WD','WS','WW','RH','QQ'
C    ,'WCLD','WRAN','WICE','WSNW','WGRP'

c --- p in mb, t in k, qsat and q in kg/kg, rh in %
      pp=xcmp(1)
      tk=xcmp(3)
      qq=xcmp(8)/1000.
      qq=max(qq,0.0)

      qs=0.622 *6.112/pp*exp(17.67* (tk-273.15)/(tk-29.65))
      rh=max(100.*qq/qs,0.)
      if(rh.gt.100.) rh=100.

      xcmp(7)=rh

      return
      end
