Program Indicator_File
implicit none
INTEGER i,j,k,z,iaa,nsub,nx,ny,nz,nz1,ndt,ii,ind1,ind2,ind3,ind4,ind,nwrite,kk,i_final
REAL*8 dx, dy,x0,y0,z0
REAL*8, DIMENSION(:,:,:), ALLOCATABLE :: clm,slopex,slopey,slopet,manning
REAL*8, DIMENSION(:), ALLOCATABLE :: et,evap_grnd,evap_soi,evap_veg,tran_veg,infl,swe
REAL*8, DIMENSION(:), ALLOCATABLE :: et_eronly,evap_grnd_eronly,evap_soi_eronly,evap_veg_eronly,tran_veg_eronly,infl_eronly,swe_eronly
REAL*8, DIMENSION(:,:), ALLOCATABLE :: et_sb,evap_grnd_sb,evap_soi_sb,evap_veg_sb,tran_veg_sb,infl_sb,swe_sb
INTEGER,DIMENSION(:,:), ALLOCATABLE :: sub_ind
INTEGER, DIMENSION(:), ALLOCATABLE :: cell_sub
INTEGER,DIMENSION(:), ALLOCATABLE :: out_k,out_j
Character (LEN=200) Nom_Fich,fichier,Nom_Fich1,fichier1,Nom_Fich2,fichier2,Nom_Fich3,fichier3,Nom_Fich4,fichier4,fid
  character*200 fname01
    nx=150
    ny=170
    !nz is number of CLM layers here
    nz=17
    ndt=8759
    !ndt=17496
    nwrite=0 
    dx=100.0
    dy=100.0
    nsub=8
    nz1=1

    ALLOCATE (cell_sub(nsub)) 
    ALLOCATE (clm(nx,ny,nz),slopex(nx,ny,nz1),slopey(nx,ny,nz1),slopet(nx,ny,nz1),manning(nx,ny,nz1))
    ALLOCATE (out_k(nsub),out_j(nsub),sub_ind(nx,ny))
    ALLOCATE (et(ndt),et_eronly(ndt),et_sb(nsub,ndt))
    ALLOCATE (evap_grnd(ndt),evap_grnd_eronly(ndt),evap_grnd_sb(nsub,ndt))
    ALLOCATE (evap_soi(ndt),evap_soi_eronly(ndt),evap_soi_sb(nsub,ndt))
    ALLOCATE (evap_veg(ndt),evap_veg_eronly(ndt),evap_veg_sb(nsub,ndt))
    ALLOCATE (tran_veg(ndt),tran_veg_eronly(ndt),tran_veg_sb(nsub,ndt))
    ALLOCATE (infl(ndt),infl_eronly(ndt),infl_sb(nsub,ndt))
    ALLOCATE (swe(ndt),swe_eronly(ndt),swe_sb(nsub,ndt))


    cell_sub=0.0


    !Domain and east river watershed only
      OPEN(29,FILE='4.0_wy2015.et_timeseries_all.txt')
      OPEN(30,FILE='4.0_wy2015.et_timeseries_ERonly_all.txt')

      OPEN(31,FILE='4.0_wy2015.evap_grnd_timeseries_all.txt')
      OPEN(32,FILE='4.0_wy2015.evap_grnd_timeseries_ERonly_all.txt')
      
      OPEN(33,FILE='4.0_wy2015.evap_soi_timeseries_all.txt')
      OPEN(34,FILE='4.0_wy2015.evap_soi_timeseries_ERonly_all.txt')

      OPEN(35,FILE='4.0_wy2015.evap_veg_timeseries_all.txt')
      OPEN(36,FILE='4.0_wy2015.evap_veg_timeseries_ERonly_all.txt')

      OPEN(37,FILE='4.0_wy2015.tran_veg_timeseries_all.txt')
      OPEN(38,FILE='4.0_wy2015.tran_veg_timeseries_ERonly_all.txt')

      OPEN(39,FILE='4.0_wy2015.infl_timeseries_all.txt')
      OPEN(40,FILE='4.0_wy2015.infl_timeseries_ERonly_all.txt')

      OPEN(41,FILE='4.0_wy2015.swe_timeseries.txt')
      OPEN(42,FILE='4.0_wy2015.swe_timeseries_ERonly_all.txt')



!Subbasins only
      OPEN(43,FILE='4.0_wy2015.et_subbasins.txt')
      OPEN(44,FILE='4.0_wy2015.evap_grnd_subbasins.txt')
      OPEN(45,FILE='4.0_wy2015.evap_soi_subbasins.txt')
      OPEN(46,FILE='4.0_wy2015.evap_veg_subbasins.txt')
      OPEN(47,FILE='4.0_wy2015.tran_veg_subbasins.txt')
      OPEN(48,FILE='4.0_wy2015.infl_subbasins.txt')
      OPEN(49,FILE='4.0_wy2015.swe_subbasins.txt')




     OPEN(20,FILE='subbasin_indicator.txt')
      READ(20,*)
      DO j=1,ny
                READ(20,*) sub_ind(1:nx,j)
      END DO
      CLOSE (20)

        DO kk=1,nsub
                DO j=1,ny
                        DO k=1,nx
                                IF (sub_ind(k,j)==kk) THEN
                                        cell_sub(kk)=cell_sub(kk)+1
                                END IF
                        END DO
                 END DO
        write(*,*) 'Number of cells in subbasin ',kk,' is ',cell_sub(kk)
        END DO




      !Main time loop here
      DO i=1,ndt,1
        
        IF ( mod(i,500) == 0 ) THEN
	        write(*,*) i
        END IF

        write(fid, '(I5.5)') i
        fname01='4.0_wy2015.out.clm_output.' // trim(adjustl(fid)) // '.C.pfb'
        write(*,*) fname01

	call pfb_read(clm,fname01,nx,ny,nz)	

            DO j=1,ny
                DO k=1,nx
                !evap_tot,5]
                !evap_grnd,6]
                !evap_soi,7]
                !evap_veg,8]
                !tran_veg,9]
                !infl,10]
                !swe,11]



                !Domain average
                        et(i)=et(i)+clm(k,j,5)
                        evap_grnd(i)=evap_grnd(i)+clm(k,j,6)
                        evap_soi(i)=evap_soi(i)+clm(k,j,7)
                        evap_veg(i)=evap_veg(i)+clm(k,j,8)
                        tran_veg(i)=tran_veg(i)+clm(k,j,9)
                        infl(i)=infl(i)+clm(k,j,10)
                        swe(i)=swe(i)+clm(k,j,11)

                        !Watershed specific
                        IF (sub_ind(k,j)>=1) THEN
                                        et_eronly(i)=et_eronly(i)+clm(k,j,5)
                                        evap_grnd_eronly(i)=evap_grnd_eronly(i)+clm(k,j,6)
                                        evap_soi_eronly(i)=evap_soi_eronly(i)+clm(k,j,7)
                                        evap_veg_eronly(i)=evap_veg_eronly(i)+clm(k,j,8)
                                        tran_veg_eronly(i)=tran_veg_eronly(i)+clm(k,j,9)
                                        infl_eronly(i)=infl_eronly(i)+clm(k,j,10)
                                        swe_eronly(i)=swe_eronly(i)+clm(k,j,11)
                        END IF

                        !Subbasin average here
                        DO kk=1,nsub
                                IF (sub_ind(k,j)==kk) THEN
                                        et_sb(kk,i)=et_sb(kk,i)+clm(k,j,5)
                                        evap_grnd_sb(kk,i)=evap_grnd_sb(kk,i)+clm(k,j,6)
                                        evap_soi_sb(kk,i)=evap_soi_sb(kk,i)+clm(k,j,7)
                                        evap_veg_sb(kk,i)=evap_veg_sb(kk,i)+clm(k,j,8)
                                        tran_veg_sb(kk,i)=tran_veg_sb(kk,i)+clm(k,j,9)
                                        infl_sb(kk,i)=infl_sb(kk,i)+clm(k,j,10)
                                        swe_sb(kk,i)=swe_sb(kk,i)+clm(k,j,11)
                                END IF
                        END DO

                END DO
            END DO


                        WRITE(29,14) i*1.0,et(i)
                        WRITE(30,14) i*1.0,et_eronly(i)

                        WRITE(31,14) i*1.0,evap_grnd(i)
                        WRITE(32,14) i*1.0,evap_grnd_eronly(i)

                        WRITE(33,14) i*1.0,evap_soi(i)
                        WRITE(34,14) i*1.0,evap_soi_eronly(i)

                        WRITE(35,14) i*1.0,evap_veg(i)
                        WRITE(36,14) i*1.0,evap_veg_eronly(i)

                        WRITE(37,14) i*1.0,tran_veg(i)
                        WRITE(38,14) i*1.0,tran_veg_eronly(i)

                        WRITE(39,14) i*1.0,infl(i)
                        WRITE(40,14) i*1.0,infl_eronly(i)

                        WRITE(41,14) i*1.0,swe(i)
                        WRITE(42,14) i*1.0,swe_eronly(i)

                        WRITE(43,14) i*1.0,et_sb(1:nsub,i)
                        WRITE(44,14) i*1.0,evap_grnd_sb(1:nsub,i)
                        WRITE(45,14) i*1.0,evap_soi_sb(1:nsub,i)
                        WRITE(46,14) i*1.0,evap_veg_sb(1:nsub,i)
                        WRITE(47,14) i*1.0,tran_veg_sb(1:nsub,i)
                        WRITE(48,14) i*1.0,infl_sb(1:nsub,i)
                        WRITE(49,14) i*1.0,swe_sb(1:nsub,i)

 END DO

 14 FORMAT(2000(1x,e13.7))

END
 
 subroutine pfb_read(value,fname,nx,ny,nz) 
  implicit none
  real*8 value(nx,ny,nz)
  real*8 dx, dy, dz, x1, y1, z1								
  integer*4 i,j,k, nni, nnj, nnk, ix, iy, iz,			&
            ns,  rx, ry, rz,nx,ny,nz, nnx, nny, nnz,is
  character*100 fname
  
  open(100,file=trim(adjustl(fname)),form='unformatted',   &
                 recordtype='stream',convert='BIG_ENDIAN',status='old') 
  ! Read in header info

! Start: reading of domain spatial information
  read(100) x1 !X
  read(100) y1 !Y 
  read(100) z1 !Z

  read(100) nx !NX
  read(100) ny !NY
  read(100) nz !NZ

  read(100) dx !DX
  read(100) dy !DY
  read(100) dz !DZ

  read(100) ns !num_subgrids
! End: reading of domain spatial information

! Start: loop over number of sub grids
  do is = 0, (ns-1)

! Start: reading of sub-grid spatial information
   read(100) ix
   read(100) iy
   read(100) iz

   read(100) nnx
   read(100) nny
   read(100) nnz
   read(100) rx
   read(100) ry
   read(100) rz

! End: reading of sub-grid spatial information

! Start: read in saturation data from each individual subgrid
  do  k=iz +1 , iz + nnz
   do  j=iy +1 , iy + nny
    do  i=ix +1 , ix + nnx
     read(100) value(i,j,k)
    end do
   end do
  end do
! End: read in saturation data from each individual subgrid

! End: read in saturation data from each individual subgrid

  end do
! End: loop over number of sub grids



  close(100)
  end subroutine


