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 :: t_grnd,t_soilc,t_soilb,t_soila REAL*8, DIMENSION(:), ALLOCATABLE :: t_grnd_eronly,t_soilc_eronly,t_soilb_eronly,t_soila_eronly REAL*8, DIMENSION(:,:), ALLOCATABLE :: t_grnd_sb,t_soilc_sb,t_soilb_sb,t_soila_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=8783 !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 (t_grnd(ndt),t_grnd_eronly(ndt),t_grnd_sb(nsub,ndt)) ALLOCATE (t_soilc(ndt),t_soilc_eronly(ndt),t_soilc_sb(nsub,ndt)) ALLOCATE (t_soilb(ndt),t_soilb_eronly(ndt),t_soilb_sb(nsub,ndt)) ALLOCATE (t_soila(ndt),t_soila_eronly(ndt),t_soila_sb(nsub,ndt)) cell_sub=0.0 !Domain and east river watershed only OPEN(53,FILE='2.5_wy2020.t_grnd_timeseries.txt') OPEN(54,FILE='2.5_wy2020.t_soilc_timeseries.txt') OPEN(55,FILE='2.5_wy2020.t_soilb_timeseries.txt') OPEN(56,FILE='2.5_wy2020.t_soila_timeseries.txt') OPEN(57,FILE='2.5_wy2020.t_grnd_timeseries_ERonly_all.txt') OPEN(58,FILE='2.5_wy2020.t_soilc_timeseries_ERonly_all.txt') OPEN(59,FILE='2.5_wy2020.t_soilb_timeseries_ERonly_all.txt') OPEN(60,FILE='2.5_wy2020.t_soila_timeseries_ERonly_all.txt') !Subbasins only OPEN(61,FILE='2.5_wy2020.t_grnd_subbasins.txt') OPEN(62,FILE='2.5_wy2020.t_soilc_subbasins.txt') OPEN(63,FILE='2.5_wy2020.t_soilb_subbasins.txt') OPEN(64,FILE='2.5_wy2020.t_soila_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='2.5_wy2020.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 !10-28-22 these have 17 layers !evap_tot,5] !evap_grnd,6] !evap_soi,7] !evap_veg,8] !tran_veg,9] !infl,10] !swe,11] !t_grnd,12 !tsoil_i=1,13 !tsoil_i=2,14 !tsoil_i=3,15 !tsoil_i=4,16 !tsoil_i=5,17 !Domain average t_grnd(i)=t_grnd(i)+clm(k,j,12) t_soilc(i)=t_soilc(i)+clm(k,j,15)+clm(k,j,16)+clm(k,j,17) t_soilb(i)=t_soilb(i)+clm(k,j,16)+clm(k,j,17) t_soila(i)=t_soila(i)+clm(k,j,17) !Watershed specific IF (sub_ind(k,j)>=1) THEN t_grnd_eronly(i)=t_grnd_eronly(i)+clm(k,j,12) t_soilc_eronly(i)=t_soilc_eronly(i)+clm(k,j,15)+clm(k,j,16)+clm(k,j,17) t_soilb_eronly(i)=t_soilb_eronly(i)+clm(k,j,16)+clm(k,j,17) t_soila_eronly(i)=t_soila_eronly(i)+clm(k,j,17) END IF !Subbasin average here DO kk=1,nsub IF (sub_ind(k,j)==kk) THEN t_grnd_sb(kk,i)=t_grnd_sb(kk,i)+clm(k,j,12) t_soilc_sb(kk,i)=t_soilc_sb(kk,i)+clm(k,j,15)+clm(k,j,16)+clm(k,j,17) t_soilb_sb(kk,i)=t_soilb_sb(kk,i)+clm(k,j,16)+clm(k,j,17) t_soila_sb(kk,i)=t_soila_sb(kk,i)+clm(k,j,17) END IF END DO END DO END DO WRITE(53,14) i*1.0,t_grnd(i) WRITE(54,14) i*1.0,t_soilc(i) WRITE(55,14) i*1.0,t_soilb(i) WRITE(56,14) i*1.0,t_soila(i) WRITE(57,14) i*1.0,t_grnd_eronly(i) WRITE(58,14) i*1.0,t_soilc_eronly(i) WRITE(59,14) i*1.0,t_soilb_eronly(i) WRITE(60,14) i*1.0,t_soila_eronly(i) WRITE(61,14) i*1.0,t_grnd_sb(1:nsub,i) WRITE(62,14) i*1.0,t_soilc_sb(1:nsub,i) WRITE(63,14) i*1.0,t_soilb_sb(1:nsub,i) WRITE(64,14) i*1.0,t_soila_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