! The datum is NAVD88.

PROGRAM slr_program
implicit none

REAL, PARAMETER :: m = 0.0014, tauc = 0.4, g = 9.82, fw = 0.25, pho=1000, sigma=100, gamma= 0.15 !gamma=0.15 !sigma = 100  ! 0.53 (0.4, 0.3) 0.52 (0.8, 0.5) pho = 1000 sigma = 1000 sigma=30 fw = 0.1 sigma = 100, tau = 0.4 tau = 0.03 sigma 100 fw of 0.1 0.15 !tauc=0.03 !m=0.0014, sigma=100, 0.0014 50 100 tauc=0.4 03
!REAL, PARAMETER :: q = 0.00009, k2 = 0.00000405 !0.0009
!REAL, PARAMETER :: q = 1360.0, k2 = 3.3, k3= 0.00001 ! 0.00003 0.00006  0.000405 !0.0000405  0.0009 0.0000105 q = 0.00009, k2 = 0.00001
!REAL, PARAMETER :: q = 113.33, k2=3.3
REAL, PARAMETER :: q = 0.00009, k2=0.00003
INTEGER :: nrow, ncol, row, col, cellsize, nodata_land, tot_cell, year, land_type, k !, land_temp
REAL :: x, y, nodata_elev, nodata_biom, nodata_tgvi1, nodata_tgvi2
CHARACTER (LEN=14) :: temp1, temp2, temp3, temp4, temp5, temp6
INTEGER, DIMENSION (:,:), allocatable :: land, land_temp, land_init, land_c
REAL, DIMENSION(:,:),allocatable :: elv, biom, tgvi1, tgvi2, elv_temp, dep, elv_prev, U, Us, elv_init
REAL :: slr = 0.0041, taub, taubs ! 0.005973, 0.0041 0.0025 0.005 (0.42) 0.0112 (0.45) 0.015+(-3) 0.0112+0.3 0.0110 0.0123 0.0133 0.0123 0.005 0.015 0.025 0.0025 slr=0.01 not working better, slr=0.02 0.002 erosion, 0.01, 0.008, 0.009, 0.005, 0.015 biom
!slr=0.009 0.0025

CHARACTER(LEN = 45) :: landfile

REAL :: prob_temp1, prob_temp2, prob_now, prob_prev

integer :: nu=0, ndep=0, nero=0, nbio = 0

real :: usum=0.0, depsum=0.0, erosum=0.0, uavg, depavg, eroavg, biosum=0.0, bioavg

! 0.005 SLR, tau 0.40: 0.420809835 (0.536426842, 0.784468234)
! 0.006 SLR, tau 0.40: 0.439275503 (0.573443174, 0.766031444)
! 0.008 SLR, tau 0.40: 0.459106684 (0.65634620, 0.699488580)
! 0.009 SLR, tau 0.40: 0.464079142 (0.718149424, 0.646215320)
! 0.010 SLR, tau 0.40: 0.464563280 (0.754749596, 0.615519762)
! 0.0112 SLR, tau 0.40: 0.459742695 (0.8397, 0.5475)
! 0.005 SLR, tau 0.03: 0.420826674 (0.536500454, 0.784392059)
! 0.008 SLR, tau 0.03: 0.459114969 (0.656676769, 0.699149132)

INTEGER :: winrow, wincol, water, workrow, workcol, saltmarsh
!REAL, DIMENSION (3240, 2719) :: erosion
!REAL, DIMENSION (5311, 3970) :: erosion
REAL, DIMENSION (2767, 3655) :: erosion, ero_s
INTEGER, DIMENSION (23) :: howmany1, howmany2
REAL :: min_elv, l, TSS
INTEGER :: i, scenario


!open (10, FILE="F:/Wu/SLR_GB/data/elv_sub3.asc", status = "OLD")
!open (20, FILE="F:/Wu/SLR_GB/data/wet88_2.asc", status = "OLD")

!open (10, FILE="F:/Wu/pub/GB_SLR2016/Data/GIS/elev_sub.asc", status = "OLD")
!open (10, FILE="F:/Wu/pub/GB_SLR2016/Data/GIS/elv_sub_afKbd.asc", status = "OLD")
!open (20, FILE="F:/Wu/SLR_GB/biomass/wet88_2_sub.asc", status = "OLD")
!open (20, FILE="F:/Wu/pub/GB_SLR2016/Data/GIS/wet_sub.asc", status = "OLD") !1988
!open(20, file="F:/Wu/pub/GB_SLR2016/Data/GIS/nwi07_sub.asc",status="old") !2007

open (10, FILE="E:/Wu/2018/GB_SLR/elv_sub_afKbd.asc", status = "OLD")
open (20, FILE="E:/Wu/2018/GB_SLR/wet_sub.asc", status = "OLD") !1988





!open (25, FILE="F:/Wu/SLR_GB/biomass/biom08_3.asc", status = "OLD")
!open (26, FILE="F:/Wu/SLR_GB/data/tgvi1.asc", status = "OLD")
!open (27, FILE="F:/Wu/SLR_GB/data/tgvi2.asc", status = "OLD")

!open (30, FILE="F:/Wu/SLR_GB/result_2016/new_land10_biom.asc", status = "UNKNOWN")
!open (30, FILE="F:/Wu/pub/GB_SLR2016/small/new_land2010_slr4.asc", status = "UNKNOWN")


!open (30, FILE="F:/Wu/pub/GB_SLR2016/small/new_land2100.asc", status = "UNKNOWN")
!open (30, FILE="F:/Wu/pub/GB_SLR2016/result/new_land2050.asc", status = "UNKNOWN")

!open (30, FILE="F:/Wu/SLR_GB/result/new_land10_biom.asc", status = "UNKNOWN") !0.0025 SLR
!open (30, FILE="F:/Wu/SLR_GB/result/new_land10_biom.asc", status = "UNKNOWN") ! 0.025 SLR
!open (30, FILE="F:/Wu/SLR_GB/result/new_land10_newbiom2.asc", status = "UNKNOWN")
!open (30, FILE="F:/Wu/SLR_GB/result/new_land10_newbiom3.asc", status = "UNKNOWN") 
!open (30, FILE="F:/Wu/SLR_GB/result/new_land10_newbiom3_try2.asc", status = "UNKNOWN")  ! 0.0123 and -0.3 kappa = 0.25
!open (30, FILE="F:/Wu/SLR_GB/result/new_land10_newbiom3_try4.asc", status = "UNKNOWN")  ! 20 years 2008
!open (30, FILE="F:/Wu/SLR_GB/result/new_land28.asc", status = "UNKNOWN")   ! 40 years 2028
!open (30, FILE="F:/Wu/SLR_GB/result/new_land48.asc", status = "UNKNOWN")    ! 60 years 2048
!open (30, FILE="F:/Wu/SLR_GB/result/new_land68.asc", status = "UNKNOWN")    ! 80 years 2068
!open (30, FILE="F:/Wu/SLR_GB/result/new_land88.asc", status = "UNKNOWN")    ! 100 years 2088
!open (30, FILE="F:/Wu/SLR_GB/result/new_land2108.asc", status = "UNKNOWN")    ! 120 years 2108
 

!open (40, FILE="F:/Wu/SLR_GB/result/dep_biom.asc", status = "UNKNOWN") 
!open (40, FILE="F:/Wu/pub/GB_SLR2016/small/dep_2010_slr4.asc", status = "UNKNOWN")
!open (50, FILE="F:/Wu/pub/GB_SLR2016/small/ero_2010_slr4.asc", status = "UNKNOWN")

!open (30, FILE="F:/Wu/pub/GB_SLR2016/small/new_land2100_slr7_newinit.asc", status = "UNKNOWN")
!open (40, FILE="F:/Wu/pub/GB_SLR2016/small/dep_2100_slr7_newinit.asc", status = "UNKNOWN")
!open (50, FILE="F:/Wu/pub/GB_SLR2016/small/ero_2100_slr7_newinit.asc", status = "UNKNOWN") 
!open (60, FILE="F:/Wu/pub/GB_SLR2016/small/U_2100_slr7_newinit.asc", status = "UNKNOWN") 
!open (70, FILE="F:/Wu/pub/GB_SLR2016/small/average_2100_newinit.csv", status = "UNKNOWN")
!open (80, FILE="F:/Wu/pub/GB_SLR2016/small/saltmarsh_2100_newinit.csv", status = "UNKNOWN")


open (30, FILE="E:/Wu/2018/GB_SLR/threshold2/new_land2100_mult.asc", status = "UNKNOWN")
open (40, FILE="E:/Wu/2018/GB_SLR/threshold2/dep_2100_mult.asc", status = "UNKNOWN")
open (50, FILE="E:/Wu/2018/GB_SLR/threshold2/ero_2100_mult.asc", status = "UNKNOWN") 
open (60, FILE="E:/Wu/2018/GB_SLR/threshold2/U_2100_mult.asc", status = "UNKNOWN") 
open (70, FILE="E:/Wu/2018/GB_SLR/threshold2/average_2100_mult.csv", status = "UNKNOWN")
open (80, FILE="E:/Wu/2018/GB_SLR/threshold2/saltmarsh_2100_mult.csv", status = "UNKNOWN")
open (90, FILE="E:/Wu/2018/GB_SLR/threshold2/slr13_2100_fixed_mult.csv", status = "UNKNOWN")

!test2 for 2100 starting from 1988, and test4 for 2100 starting from 2007 

!open (40, FILE="F:/Wu/pub/GB_SLR2016/small/dep_2100.asc", status = "UNKNOWN")
!open (50, FILE="F:/Wu/pub/GB_SLR2016/small/ero_2100.asc", status = "UNKNOWN") 

landfile(1:38) = "E:/Wu/2018/GB_SLR/threshold2/land2100_"
landfile (42:45) = ".asc"

read (10, *) temp1, ncol
write (*, *) temp1, ncol
read (10, *) temp2, nrow
write (*, *) temp2, nrow
read (10, *) temp3, x
read (10, *) temp4, y
read (10, *) temp5, cellsize
write (*,*) cellsize
read (10, *) temp6, nodata_elev

read (20, *) temp1, ncol
write (*, *) temp1, ncol
read (20, *) temp2, nrow
read (20, *) temp3, x
read (20, *) temp4, y
read (20, *) temp5, cellsize
read (20, *) temp6, nodata_land

!read (25, *) temp1, ncol
!write (*, *) temp1, ncol
!read (25, *) temp2, nrow
!write (*, *) temp2, nrow
!read (25, *) temp3, x
!read (25, *) temp4, y
!read (25, *) temp5, cellsize
!read (25, *) temp6, nodata_elev

!read (26, *) temp1, ncol
!write (*, *) temp1, ncol
!read (26, *) temp2, nrow
!write (*, *) temp2, nrow
!read (26, *) temp3, x
!read (26, *) temp4, y
!read (26, *) temp5, cellsize
!read (26, *) temp6, nodata_tgvi1

!read (27, *) temp1, ncol
!write (*, *) temp1, ncol
!read (27, *) temp2, nrow
!write (*, *) temp2, nrow
!read (27, *) temp3, x
!read (27, *) temp4, y
!read (27, *) temp5, cellsize
!read (27, *) temp6, nodata_tgvi2

allocate(land(nrow,ncol))
allocate(elv(nrow,ncol))
allocate(biom(nrow,ncol))
!allocate(tgvi1(nrow,ncol))
!allocate(tgvi2(nrow,ncol))
allocate(land_temp(nrow,ncol))
allocate(elv_temp(nrow,ncol))
allocate(dep(nrow,ncol))
allocate(elv_prev(nrow,ncol))
allocate(U(nrow,ncol))
allocate(Us(nrow,ncol))
allocate(elv_init(nrow,ncol))
allocate(land_init(nrow,ncol))
allocate(land_c(nrow,ncol))

write (30, *) temp1, ncol
write (30, *) temp2, nrow
write(30, *) temp3, x
write (30, *) temp4, y
write (30, *) temp5, cellsize
write (30, *) temp6, nodata_land

write (40, *) temp1, ncol
write (40, *) temp2, nrow
write(40, *) temp3, x
write (40, *) temp4, y
write (40, *) temp5, cellsize
write (40, *) temp6, nodata_elev

write (50, *) temp1, ncol
write (50, *) temp2, nrow
write(50, *) temp3, x
write (50, *) temp4, y
write (50, *) temp5, cellsize
write (50, *) temp6, nodata_elev

write (60, *) temp1, ncol
write (60, *) temp2, nrow
write(60, *) temp3, x
write (60, *) temp4, y
write (60, *) temp5, cellsize
write (60, *) temp6, nodata_elev

do row=1, nrow
        read(10, *) (elv(row,col),col=1, ncol)
        read(20, *) (land(row,col),col=1,ncol)
 !       read(25, *) (biom(row,col), col=1,ncol)
 !       read(26, *) (tgvi1(row,col), col=1,ncol)
 !       read(27, *) (tgvi2(row,col), col=1,ncol)
end do

elv_init = elv
land_init = land

do land_type = 1, 23
    howmany1(land_type) = 0
end do

do row = 1, nrow
    do col = 1, ncol
        if (land(row,col) /= nodata_land ) then
           k = land(row,col)
           howmany1(k) = howmany1(k) + 1
        end if
    end do
end do


!do row = 1, nrow
!    do col = 1, ncol
!        erosion(row,col) = nodata_elev
!        dep(row,col) = nodata_elev
!        biom(row,col) = 0.0
  !      if (tgvi1(row,col) /= nodata_tgvi1) then
  !         biom(row,col) = 759.7688+ 1895.2964*tgvi1(row,col)
  !         if (biom(row,col) < 0) then
  !             biom (row,col) = 0
  !         end if
  !      end if
!    end do
!end do
do row = 1, nrow
    do col = 1, ncol
        erosion(row,col) = nodata_elev
        dep(row,col) = nodata_elev
        biom(row,col) = nodata_elev
        U(row,col) = nodata_elev
        Us(row,col) = nodata_elev
    end do
end do

write(70,*) "year, u, cells, dep, cells, ero, cells, biom, cells"

do scenario = 40, 200, 5
    
  if (scenario < 100) then
      write (landfile(39:41), 10) 0, scenario
  else
      write (landfile(39:41), 11) scenario
  end if
  
  open (100+scenario, file=landfile, status="UNKNOWN")
   
   slr = scenario/10000.0   
   elv = elv_init
   land = land_init
   
   do row = 1, nrow
    do col = 1, ncol
        erosion(row,col) = nodata_elev
        dep(row,col) = nodata_elev
        biom(row,col) = nodata_elev
        U(row,col) = nodata_elev
    end do
end do

do year = 1, 112! 112 !62 !43 93
    
!slr = slr !*1.02
    
write (*,*) "Year", year, slr
!do year = 1, 39
!do year = 1, 59   
!do year = 1, 79   
!do year = 1, 99   
!do year = 1, 119   
  do row = 1, nrow
    do col = 1, ncol            
        if (elv(row,col) /= nodata_elev) then
            elv(row,col) = elv(row,col) - slr
            
            if (land(row,col) == 20) then  !land(row,col) /= nodata_land .and. 
            
        if (elv(row,col) < 0.323) then !  0.323 0.273 = MHW (0.753 m) - NAVD (0.480 m) 
    !https://tidesandcurrents.noaa.gov/datums.html?units=1&epoch=0&id=8740166&name=Grand+Bay+NERR%2C+Mississippi+Sound&state=MS
          
            biom(row,col) = 864.2825 - 1022.8792 * elv(row,col) * elv(row,col)
            
            if (biom(row,col) <= 0 .and. biom(row,col) > nodata_elev) then
                biom (row,col) = 0
            end if
            
            if (biom(row,col) >= 4000.0) then
                biom (row,col) = 4000.0
            end if
            
                TSS = 12.22 + 0.002307*row !SST_coord.csv, SST.r
!dep(row,col) = (TSS*(q+k2*biom(row,col))*(0.273 - elv(row,col))*(0.273 - elv(row,col))/0.417 + 4.6*biom(row,col)*0.1)/(716.7*1000) !0.273  4.6*biom(row,col)*0.1/716.7 !k3*0.4*biom(row,col)      ! 0.053 0.17 0.273 0.352
            if (elv(row,col) < 0.273) then
                dep(row,col) = (TSS*q+k2*biom(row,col))*(0.273-elv(row,col))+4.6*biom(row,col)*0.15/(716.7*1000) !refractory 0.1 0.273
            ! 716.7 is the bulk density, see Bulk Density_pant.xls
            else
                dep(row,col)=4.6*biom(row,col)*0.15/(716.7*1000)
            end if
    
 !   write(*,*) dep(row,col)
    
 !               dep(row,col) =  0.26*biom(row,col)*0.05/716.7 !0.273  4.6*biom(row,col)*0.1/716.7 !k3*0.4*biom(row,col)      ! 0.053 0.17 0.273 0.352
        
             if (dep(row,col) <0 .and. dep(row,col) > nodata_elev) then
                 dep(row,col) = 0
             end if
             
!             dep(row,col) = dep(row,col) + 0.26*biom(row,col)*0.05/716.7
             
        !     if (dep(row,col) > 0.10) then !0.05
        !         dep(row,col) = 0.10
        !     end if
             
            if (elv(row,col) < 0.273) then 
               U(row,col) = gamma/2*sqrt(g*((0.273-elv(row, col))))  ! 0.273 U=0.071 0.302 0.62 0.27 0.31 0.92 0.31 abs(elv()) 0.302
            else
               U(row,col) = 0.0
            end if
            
            if (year == 16) then
               if (elv(row,col) < 0.273) then
                   Us(row,col) = 10*gamma/2*sqrt(g*(0.273-elv(row, col)))
               else
                   Us(row,col) = 0.0
               end if
            end if
   
           if (U(row,col) > nodata_elev .and. U(row,col) <= 0.0) U(row,col) = 0.0
  
               taub = pho*fw*U(row,col)*U(row,col)/8.0
  !             write (*,*) taub, U(row,col)
               if(year == 16) taubs = pho*fw*Us(row,col)*Us(row,col)/8.0
  
                water = 0
                min_elv = 9999.0
                do winrow = (row-2), (row+2)   ! 5 2 3
                    do wincol = (col-2), (col+2)
                        if (winrow <= 0) then
                            workrow = 1
                        else if (winrow >= nrow) then
                            workrow = nrow
                        else
                            workrow = winrow
                        end if
                        if (wincol <= 0) then
                            workcol = 1
                        else if (wincol >= ncol) then
                            workcol = ncol
                        else
                            workcol = wincol
                        end if
                        if(land(row,col) /= 17 .and. land(workrow,workcol) == 17) then
                            water = water + 1
                            
                            if (elv(workrow,workcol) < min_elv) then
                                min_elv = elv(workrow,workcol)
                            end if
                        end if
                    end do
                end do
          
          if (year /= 16) then
               if (taub < tauc ) then !.or. water == 0
                   erosion(row,col) = 0
               else 
                   erosion(row,col) = (m*(taub-tauc)/tauc)*sigma/716.7*365 !300.0*365  391.5*365
               end if
          else
               if (taub < tauc ) then !.or. water == 0
                   erosion(row,col) = 0
               else 
                   erosion(row,col) = (m*(taub-tauc)/tauc)*sigma/716.7*363 !300.0*365  391.5*365
               end if
               
               if (taubs < tauc) then
                    ero_s(row,col) = 0
               else 
                    ero_s(row,col) = (m*(taubs-tauc)/tauc)*sigma/716.7*2
               end if
               erosion(row,col) = erosion(row,col)+ero_s(row,col)                
          end if
               
               
               if (erosion(row,col) < 0 .and. erosion(row,col) > nodata_elev) erosion(row,col) = 0
    
   !         if (erosion(row,col) > 0.1) then
   !             erosion(row,col) = 0.1
   !         end if
            
          if (erosion(row,col) >= 0 .and. dep(row,col) >= 0) then 
            elv_temp (row,col) = elv(row,col) -erosion (row,col) + dep(row,col) !- slr !- erosion(row,col)/300 * 3600 * 24 * 365
           else if (erosion(row,col) >= 0) then
            elv_temp(row,col) = elv(row,col) - erosion(row,col) !- slr
           else if (dep(row,col) >= 0) then
            elv_temp(row,col) = elv(row,col) + dep(row,col) !- slr
          end if
        else
           elv_temp(row,col) = elv(row,col) !- slr
        END IF
            end if 
            end if
        
   !      if (elv(row,col) /= nodata_elev) elv_temp(row, col) = elv_temp(row,col) - slr
         
        END DO
  END DO
  elv_prev = elv
  elv = elv_temp
  
  land_temp = land
  
  DO ROW = 1, NROW
     DO COL = 1, NCOL
   !      if (land(row,col) /= nodata_land .and. elv(row,col) /= nodata_elev) then
 !           if (elv(row,col) < 0.065 .and. (land(row,col) == 20 )) then ! -0.17 (2.5 percentile of land 20), 0.277 0.065 0.277 -0.183 -0.144 0.053 0.065 .or. land(row,col)==10  0.302 0.16, 0.28, 0.3.or. land(row,col) == 6
            if ((land(row,col) == 20 )) then
                
  !      if ((land(row,col) == 20 .or. land(row,col)==10 )) then    
  !             land(row,col)= 17
  
                  water = 0
                min_elv = 9999.0
                do winrow = (row-2), (row+2)   !2 5 3
                    do wincol = (col-2), (col+2)
                        if (winrow <= 0) then
                            workrow = 1
                        else if (winrow >= nrow) then
                            workrow = nrow
                        else
                            workrow = winrow
                        end if
                        if (wincol <= 0) then
                            workcol = 1
                        else if (wincol >= ncol) then
                            workcol = ncol
                        else
                            workcol = wincol
                        end if
                        if(land(workrow,workcol) == 17) then
                            water = water + 1
                            
                            if (elv(workrow,workcol) < min_elv) then
                                min_elv = elv(workrow,workcol)
                            end if
                        end if
                    end do
                end do
                    
 !                 prob_temp1 = exp(1.9464508-0.4466820*elv(row,col)-0.0401783*elv(row,col)*elv(row,col))
 !                 prob_now = prob_temp1/(1+prob_temp1)
                  
 !                prob_temp2 = exp(1.9464508-0.4466820*elv_prev(row,col)-0.0401783*elv_prev(row,col)*elv_prev(row,col))
 !                 prob_prev = prob_temp2/(1+prob_temp2)
                  
                  prob_temp1 = exp(1.158594+4.051589*elv(row,col))
                  prob_now = prob_temp1/(1+prob_temp1)
                  
                 prob_temp2 = exp(1.158594+4.051589*elv_prev(row,col))
                  prob_prev = prob_temp2/(1+prob_temp2)
 
 !                 if (prob_now < prob_prev .and. prob_now <0.5 .and. water > 0) then  !0.7612432 .and. elv(row,col) <= min_elv
                  if (elv(row,col) < -0.144) then   ! .and. water > 0 -0.17 -0.144 -0.384 .and. water > 0, -0.457: -0.73 below MHW, -0.384: -0.24 m below MLW, -0.183 (MLLW) 0.7612432 .and. elv(row,col) <= min_elv
 !                MLW (4.336 m) - NAVD (4.480 m) = -0.144 m !0.384
                    land_temp (row,col) = 17
                     biom(row,col) = 0.0  ! new addition 08/22/2016
                     U(row,col) = 0.0
                     Us(row,col) = 0.0
                     erosion(row,col) = 0.0
                     dep(row,col) = 0.0
   !                  land(row,col) = land_temp
                  end if
         end if
    end do
  end do
land = land_temp

nu=0
ndep=0
nero=0
nbio = 0

usum=0.0
depsum=0.0
erosum=0.0
biosum = 0.0

do row = 1, nrow
    do col = 1, ncol
      if (land(row,col) == 20) then
        if (U(row,col) >  nodata_elev .and. elv(row,col) < 0.323) then
            nu = nu+1
            usum = usum+U(row,col)
  !      else 
  !          U(row,col) = nodata_elev
        end if
        if (dep(row,col) > nodata_elev .and. elv(row,col) < 0.323 ) then
            ndep=ndep+1
            depsum=depsum+dep(row,col)
  !      else
  !          dep(row,col) = nodata_elev
        end if
        if(erosion(row,col) > nodata_elev .and. elv(row,col) < 0.323) then
            nero = nero+1
            erosum=erosum+erosion(row,col)
  !      else
  !          erosion(row,col) = nodata_elev
        end if
        if(biom(row,col) >  nodata_elev .and. elv(row,col) < 0.323) then
            nbio = nbio+1
            biosum=biosum+biom(row,col)
  !      else
  !          biom(row,col) = nodata_elev
        end if
      end if
    end do
end do

uavg = usum/nu
depavg = depsum/ndep
eroavg = erosum/nero
bioavg = biosum/nbio

write(70,*) year, ",", uavg, ",",nu, ",",depavg, ",",ndep, ",",eroavg, ",", nero,",", bioavg, ",", nbio

do land_type = 1, 23
    howmany2(land_type) = 0
end do

do row = 1, nrow
    do col = 1, ncol
        if (land(row,col) /= nodata_land ) then
           k = land(row,col)
           howmany2(k) = howmany2(k) + 1
        end if
    end do
end do

write (80, *) year, ",", howmany2(20)

!if (year == 1) then 
!    do row = 1, nrow
!        write (40, *) (dep(row,col), col=1, ncol)
!    end do
!end if
        

end do ! end year

write (90, *) slr, ",", howmany2(20)

!do row = 1, nrow
!    do col = 1, ncol
!        if (land(row,col) == 20) then
!            land_c(row,col) = 1
!        else
!            land_c(row,col) = nodata_land
!        end if
!    end do
!end do

do row=1, nrow
!    write (100+scenario, *) (land_c(row,col), col=1, ncol)
    write (100+scenario, *) (land(row,col), col=1, ncol)
end do

close(100+scenario)

end do ! scenario


do row=1, nrow
    write (30, *) (land(row,col), col=1, ncol)
!    write (40, *) (dep(row,col), col=1, ncol)
 !   write (50, *) (erosion(row,col), col=1, ncol)
 !   write (60, *) (U(row,col), col=1, ncol)
end do
            

do land_type = 1, 23
    howmany2(land_type) = 0
end do

do row = 1, nrow
    do col = 1, ncol
        if (land(row,col) /= nodata_land ) then
           k = land(row,col)
           howmany2(k) = howmany2(k) + 1
        end if
    end do
end do

do land_type = 1, 23
    write(*,*) land_type, howmany1(land_type), howmany2(land_type)
end do

deallocate(land)
deallocate(elv)
deallocate(biom)
deallocate(land_temp)
deallocate(elv_temp)
deallocate(dep)
deallocate(elv_prev)
deallocate(U)
deallocate(land_c)

close(10)
close(20)
close(25)
close(30)
close(40)
close(50)
close(60)
close(70)
close(80)
close(90)

10 format(I1, I2)
11 format(I3)

STOP
end program slr_program