! This script computes the total wood harvested biomass per land-related unit area.
! It is computed as the total wood biomass (primf_bioh+primn_bioh+secmf_bioh+secyf_bioh+secnf_bioh) divided by the related land cell area ((primf+primn+secdf+secdn)*area).

! Data from LUH2v2 - http://luh.umd.edu/data.shtml
use "/home/basc/dnarayan/LandUsemap_ORCHIDEE/LUH2v2h/staticData_quarterdeg.nc"


let file1={SPAWN: "find /home/surface7/vuichard/LandUseMap_CMIP6/UofMD-landState-$1/land/yr/multiple-states/gn/ -name multiple-states_input4MIPs_landState_ScenarioMIP_UofMD-$1_gn_2100-2300.nc"}

list file1

use "`file1`"

let file2={SPAWN: "find /home/surface7/vuichard/LandUseMap_CMIP6/UofMD-landState-$1/land/yr/multiple-transitions/gn/ -name multiple-transitions_input4MIPs_landState_ScenarioMIP_UofMD-$1_gn_2100-2300.nc"}

list file2

use "`file2`"

! Bioh
! Total harvest wood biomass (in kg)
! Conversion from kg to g
let bioh=(primf_bioh+primn_bioh+secmf_bioh+secyf_bioh+secnf_bioh)*1000

! Land
! Total land fraction over which wood harvest acts
let land=primf+primn+secdf+secdn


! Total Grid-cell Area 
! Conversion from km2 to m2
let area=carea*1e6

! Define harvest wood biomass per m2 of "forest" 
let woodharvest=if(land[d=2]) then bioh[d=3]/(land[d=2]*area[d=1]) else 0.




repeat/l=1:200 save/file="/home/surface7/vuichard/LandUseMap_CMIP6/UofMD-landState-$1/outputs/$1_woodharvest_`l+2099`.nc" woodharvest[l=`l`]