      program lugrid
      use M3UTILIO

c needed use M3UTILIO with intel17 ioapi libs.

c add lufrac info from LUFRAC_CRO file with INGRID data to combined OUTGRID.
c this program assumes:
c 1. only one variable on input LUFRAC file called LUFRAC with units=percent and var_desc="fractional land use"
c 2. each timestep in LUFRAC is the same data as each other timestep's data, so only one timestep is read for inputs.
c 3. only one level input for the ingrid file

      implicit none

c atmos intel17.0.3 netcdf4.4.1 ioapi3.2:

c17 ifort -fPIC -check -traceback -extend_source -zero -O3 -mp1 -o lugrid.exe.atm lugrid.f -L/home/local-rhel7/apps/netcdf-4.4.1/intel-17.0/lib -lnetcdf -lnetcdff -L/home/local-rhel7/apps/ioapi-3.2/intel-17.0/lib -lioapi -module /home/local-rhel7/apps/ioapi-3.2/intel-17.0/Linux2_x86_64ifort

      real,    allocatable :: lufrac(:,:)  ! x,y lufrac of LUFRAC_CRO 1st timestep, layer-n going to lufrac_n in GRIDCRO2D.
      real,    allocatable :: ingrid(:,:,:)  ! col,row,var
      character*16, allocatable :: luvar(:),ingridvar(:),outgridvar(:)
      character*16, allocatable :: luuni(:),ingriduni(:),outgriduni(:)
      character*80, allocatable :: ludes(:),ingriddes(:),outgriddes(:)

      integer,      allocatable :: lutyp(:),ingridtyp(:),outgridtyp(:)

      integer i,j,k,n,nfrac,nvars3
      integer jtime,jdate,tstep,mxts
      character*16 vname
      character*2 chnfrac

      INTEGER LOGUNIT,ichk

c logical netcdf file names

c LUFRAC_CRO
      integer  FTYPE1 
      integer  SDATE1
      integer  STIME1
      integer  TSTEP1
      integer  NTHIK1
      integer  NCOLS1
      integer  NROWS1
      integer  NLAYS1
      integer  NVARS1
      integer  GDTYP1
      real*8   P_ALP1
      real*8   P_BET1
      real*8   P_GAM1
      real*8   XCENT1
      real*8   YCENT1
      real*8   XORIG1
      real*8   YORIG1
      real*8   XCELL1
      real*8   YCELL1
      integer  VGTYP1
      real     VGTOP1
      integer  mxrec1

c INGRID
      integer  FTYPE2
      integer  SDATE2
      integer  STIME2
      integer  TSTEP2
      integer  NTHIK2
      integer  NCOLS2
      integer  NROWS2
      integer  NLAYS2
      integer  NVARS2
      integer  GDTYP2
      real*8   P_ALP2
      real*8   P_BET2
      real*8   P_GAM2
      real*8   XCENT2
      real*8   YCENT2
      real*8   XORIG2
      real*8   YORIG2
      real*8   XCELL2
      real*8   YCELL2
      integer  VGTYP2
      real     VGTOP2
      integer  mxrec2

      LOGUNIT = INIT3( )

c open stack params input file and get its dimensions and variable attributes
      IF ( .not. OPEN3( 'LUFRAC', FSread3, upnam3d ) ) THEN
        write(6,*) 'failure at open3 of lufrac'
        goto 999
      END IF

      if ( .not. desc3( 'LUFRAC' ) ) then
        write(6,*) 'failure to read stacks header'
        goto 999
      end if

c save lufrac input file description
      FTYPE1 = FTYPE3D
      SDATE1 = SDATE3D
      STIME1 = STIME3D
      TSTEP1 = TSTEP3D
      NTHIK1 = NTHIK3D
      NCOLS1 = NCOLS3D
      NROWS1 = NROWS3D
      NLAYS1 = NLAYS3D
      NVARS1 = NVARS3D
      GDTYP1 = GDTYP3D
      P_ALP1 = P_ALP3D
      P_BET1 = P_BET3D
      P_GAM1 = P_GAM3D
      XCENT1 = XCENT3D
      YCENT1 = YCENT3D
      XORIG1 = XORIG3D
      YORIG1 = YORIG3D
      XCELL1 = XCELL3D
      YCELL1 = YCELL3D
      VGTYP1 = VGTYP3D
      VGTOP1 = VGTOP3D
      mxrec1 = mxrec3d

      if ( nvars1 .ne. 1 ) then
         print*,'lufrac input must have just one variable: abort'
         goto 999
      end if

      nfrac = NLAYS3D
      allocate ( luvar(nvars1))
      allocate ( luuni(nvars1))
      allocate ( ludes(nvars1))
      allocate ( lutyp(nvars1))
      allocate ( lufrac(ncols1,nrows1))

      luvar(1:nvars1) = vname3d(1:nvars1)
      luuni(1:nvars1) = units3d(1:nvars1)
      ludes(1:nvars1) = vdesc3d(1:nvars1)
      lutyp(1:nvars1) = vtype3d(1:nvars1)

c open ingrid values input file and get its dimensions and variable attributes

      IF ( .not. OPEN3( 'INGRID', FSread3, upnam3d ) ) THEN
        write(6,*) 'failure at open3 of ingrid'
        goto 999
      END IF

      if ( .not. desc3( 'INGRID' ) ) then
        write(6,*) 'failure to read inline ingrid'
        goto 999
      end if

c save ingrid input file description
      FTYPE2 = FTYPE3D
      SDATE2 = SDATE3D
      STIME2 = STIME3D
      TSTEP2 = TSTEP3D
      NTHIK2 = NTHIK3D
      NCOLS2 = NCOLS3D
      NROWS2 = NROWS3D
      NLAYS2 = NLAYS3D
      NVARS2 = NVARS3D
      GDTYP2 = GDTYP3D
      P_ALP2 = P_ALP3D
      P_BET2 = P_BET3D
      P_GAM2 = P_GAM3D
      XCENT2 = XCENT3D
      YCENT2 = YCENT3D
      XORIG2 = XORIG3D
      YORIG2 = YORIG3D
      XCELL2 = XCELL3D
      YCELL2 = YCELL3D
      VGTYP2 = VGTYP3D
      VGTOP2 = VGTOP3D
      mxrec2 = mxrec3d

      allocate ( ingrid(ncols2,nrows2,nvars2))
      allocate ( ingridvar(nvars2))
      allocate ( ingriduni(nvars2))
      allocate ( ingriddes(nvars2))
      allocate ( ingridtyp(nvars2))

      ingridvar(1:nvars2) = vname3d(1:nvars2)
      ingriduni(1:nvars2) = units3d(1:nvars2)
      ingriddes(1:nvars2) = vdesc3d(1:nvars2)
      ingridtyp(1:nvars2) = vtype3d(1:nvars2)

c compare lufrac and gridcro file headers

      if ( FTYPE1 .ne. FTYPE2 ) stop 'FTYPE3D inconsistent'
c     if ( SDATE1 .ne. SDATE2 ) print*,'SDATE3D inconsistent'
c     if ( STIME1 .ne. STIME2 ) print*,'STIME3D inconsistent'
c     if ( TSTEP1 .ne. TSTEP2 ) print*,'TSTEP3D inconsistent'
c     if ( NTHIK1 .ne. NTHIK2 ) print*,'NTHIK3D inconsistent'
      if ( NCOLS1 .ne. NCOLS2 ) stop 'NCOLS3D inconsistent'
      if ( NROWS1 .ne. NROWS2 ) stop 'NROWS3D inconsistent'
c     if ( NVARS1 .ne. NVARS2 ) print*,'NVARS3D inconsistent'
      if ( GDTYP1 .ne. GDTYP2 ) stop 'GDTYP3D inconsistent'
      if ( P_ALP1 .ne. P_ALP2 ) stop 'P_ALP3D inconsistent'
      if ( P_BET1 .ne. P_BET2 ) stop 'P_BET3D inconsistent'
      if ( P_GAM1 .ne. P_GAM2 ) stop 'P_GAM3D inconsistent'
      if ( XCENT1 .ne. XCENT2 ) stop 'XCENT3D inconsistent'
      if ( YCENT1 .ne. YCENT2 ) stop 'YCENT3D inconsistent'
      if ( XORIG1 .ne. XORIG2 ) stop 'XORIG3D inconsistent'
      if ( YORIG1 .ne. YORIG2 ) stop 'YORIG3D inconsistent'
      if ( XCELL1 .ne. XCELL2 ) stop 'XCELL3D inconsistent'
      if ( YCELL1 .ne. YCELL2 ) stop 'YCELL3D inconsistent'
c     if ( VGTYP1 .ne. VGTYP2 ) print*,'VGTYP3D inconsistent'
c     if ( VGTOP1 .ne. VGTOP2 ) print*,'VGTOP3D inconsistent'
                       
c create output arrays

      nvars3 = nvars2 + nfrac  ! ingrid + lufrac
      allocate ( outgridvar(nvars3))
      allocate ( outgriduni(nvars3))
      allocate ( outgriddes(nvars3))
      allocate ( outgridtyp(nvars3))

      outgridvar(1:nvars2) = ingridvar(1:nvars2)
      outgriduni(1:nvars2) = ingriduni(1:nvars2)
      outgriddes(1:nvars2) = ingriddes(1:nvars2)
      outgridtyp(1:nvars2) = ingridtyp(1:nvars2)

      outgriduni(nvars2+1:nvars3) = luuni(1)  ! assumes just one lufrac file var input
      outgriddes(nvars2+1:nvars3) = ludes(1)
      outgridtyp(nvars2+1:nvars3) = lutyp(1)

      do n = 1,nfrac
         k = nvars2 + n
         write(chnfrac,'(i2)')n ! should range from 0-99
         if ( n .lt. 10 ) chnfrac = '0'//chnfrac(2:2)
         outgridvar(k) = 'LUFRAC_'//chnfrac
      enddo

c set up the output file by resetting some parameters of the ingrid file to outgrid specs.

      if ( .not. desc3( 'INGRID' ) ) then
        write(6,*) 'failure to read inline ingrid'
        goto 999
      end if

      nvars3d = nvars3
      vname3d(1:nvars3) = outgridvar
      units3d(1:nvars3) = outgriduni
      vdesc3d(1:nvars3) = outgriddes
      vtype3d(1:nvars3) = outgridtyp

      print*,'pre open out'

      IF ( .not. OPEN3( 'OUTGRID', FSnew3, upnam3d ) ) THEN
        write(6,*) 'failure at open3 of outgrid'
        goto 999
      END IF

      print*,'aft open out'

c first copy gridin to gridout

      jdate = sdate2
      jtime = stime2

      print*,'pre io ingrid'

      do k = 1,nvars2
         vname = ingridvar(k)
         if ( .not.  read3('INGRID', vname,1,jdate,jtime,ingrid))goto 999
         if ( .not. write3('OUTGRID',vname,jdate,jtime,ingrid))goto 999
      enddo

      print*,'aft io ingrid'

c second transfrer the landfrac from lufrac to outgrid

      print*,' pre io lufrac '

      do n = 1,nfrac
         k = nvars2 + n
         if ( .not. read3('LUFRAC', luvar,n,sdate1,stime1,lufrac))goto 999
         vname = outgridvar(k)
         if ( .not. write3('OUTGRID',vname,sdate2,stime2,lufrac))goto 999
      enddo

      print*,' aft io lufrac '

 999  continue
      IF ( .NOT. SHUT3( ) ) THEN
        write(20,*)'failure at shut3'
      END IF

      print*,' end '

 998  continue
      stop
      end
