ccsm4_nitdep.110610.ncl

;=========================================================;
; for nitrogen deposition CN forcing:
; read in modern template clm nit dep forcing file
; read in paleo landmask
; read in landmask for default (modern) grid
; all assumed to be on the same grid, if not, add regridding code or do offline
; make modifications; use zonal averages lnd/ocean specific
; write out new files
; using modern number of timeslices to compare well to modern; time is arbitrary unless doing transient
; use ONLY 1849 (first) timeslice for values
; NOTE: mods are paleo period dependent. only use this code as a GUIDE
; originally written for PT (pangean style continent, symmetric about equator, assuming little dep at poles)
;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

;====================user defined parameters

filename_land= "/your/path/here/yourtopo_1.9x2.5_topodepth_0-360.yymmdd.nc" ; edit for you filename

filename_defland= "/your/inputdata/path/here/for/default_1850/LANDFRAC_cami_0000-09-01_1.9x2.5_L26_c040809.nc"

f_template_nd= "/your/inputpath/here/for/default_1850/fndep_clm_hist_simyr1849-2006_1.9x2.5_c100428.nc"

new_cdf         = "fndep_clm_hist_PALEOCASE_1.9x2.5_cyymmdd.nc"  ;edit paleocase and yymmdd


;=============== read in atm template netcdf files
; note: same grid for all
 
fnd = addfile(f_template_nd,"r")
lat  = fnd->lat
lon  = fnd->lon
time = fnd->time
YEAR = fnd->YEAR
NDEP_year = fnd->NDEP_year
NDEP_NOy_year = fnd->NDEP_NOy_year
NDEP_NHx_year = fnd->NDEP_NHx_year

nlat = dimsizes(lat)
nlon = dimsizes(lon)
ntimes=  dimsizes(time)
print("nlon = " + nlon)
print("nlat = " + nlat)

;=============== read in land/topo on paleo grid
;                read in  modern landmask on same grid as paleo

;paleo
pland = addfile(filename_land,"r")
topo = pland->topo
plat = pland->lat
plon = pland->lon
nplat = dimsizes(plat)
nplon = dimsizes(plon)

;modern
mland = addfile(filename_defland,"r")
mlf = mland->LANDFRAC(0,:,:)
mlat = mland->lat
mlon = mland->lon

if(dimsizes(plat).ne.dimsizes(mlat))then
 print("paleo and modern should be on same grid")
 exit
end if


;================= make mods

;-----computing 1850 (default) avg values over land and ocean per latitude

;conform landmask to default array
lfC = conform(NDEP_year,mlf,(/1,2/))

;mask over ocean, land pts only
lndep = mask(NDEP_year,lfC.ne.0,True)
lnoy = mask(NDEP_NOy_year,lfC.ne.0,True)
lnhx = mask(NDEP_NHx_year,lfC.ne.0,True)
print("ocean pt should be missing: " +lndep(0,48,78))

;mask over land, ocean pts only
ondep = mask(NDEP_year,lfC.eq.0,True)
onoy = mask(NDEP_NOy_year,lfC.eq.0,True)
onhx = mask(NDEP_NHx_year,lfC.eq.0,True)
print("ocean pt should have value: " +ondep(0,48,78))

;total modern zonal avg vars
ndepav     = new((/ntimes,nplat/),float)
noyav     = new((/ntimes,nplat/),float)
nhxav     = new((/ntimes,nplat/),float)

;land modern zonal avg
lndepav = new((/ntimes,nplat/),float)
lnoyav = new((/ntimes,nplat/),float)
lnhxav = new((/ntimes,nplat/),float)

;ocean modern zonal avg
ondepav = new((/ntimes,nplat/),float)
onoyav = new((/ntimes,nplat/),float)
onhxav = new((/ntimes,nplat/),float)

do nlt = 0,nplat-1
;total
 ndepav(:,nlt) = dim_avg(NDEP_year(:,nlt,:))
 noyav(:,nlt) = dim_avg(NDEP_NOy_year(:,nlt,:))
 nhxav(:,nlt) = dim_avg(NDEP_NHx_year(:,nlt,:))
;land
 lndepav(:,nlt) = dim_avg(lndep(:,nlt,:))
 lnoyav(:,nlt) = dim_avg(lnoy(:,nlt,:))
 lnhxav(:,nlt) = dim_avg(lnhx(:,nlt,:))
;ocean
 ondepav(:,nlt) = dim_avg(ondep(:,nlt,:))
 onoyav(:,nlt) = dim_avg(onoy(:,nlt,:))
 onhxav(:,nlt) = dim_avg(onhx(:,nlt,:))
end do

;global wtg for filling in msg data
rad    = 4.0*atan(1.0)/180.0
clat   = cos(plat*rad)
ndep_g = wgt_areaave(NDEP_year, clat, 1.0, 0)
noy_g = wgt_areaave(NDEP_NOy_year, clat, 1.0, 0)
nhx_g = wgt_areaave(NDEP_NHx_year, clat, 1.0, 0)

print("modern zonal averages computed")


;full paleo grid
pndep     = new((/ntimes,nplat,nplon/),float)
pndep_noy = new((/ntimes,nplat,nplon/),float)
pndep_nhx = new((/ntimes,nplat,nplon/),float)

;zero out
pndep = 0.
pndep_noy = 0.
pndep_nhx = 0.

;--- assign paleo values over land/ocean per latitude
;modify over land/ocean
do nln = 0,nplon-1
 do nlt = 0,nplat-1
 if(topo(nlt,nln).ge.0)then ; fill in land
  pndep(:,nlt,nln) = (/lndepav(:,nlt)/)
  pndep_noy(:,nlt,nln) = (/lnoyav(:,nlt)/)
  pndep_nhx(:,nlt,nln) = (/lnhxav(:,nlt)/)
 else
  pndep(:,nlt,nln) = (/ondepav(:,nlt)/)
  pndep_noy(:,nlt,nln) = (/onoyav(:,nlt)/)
  pndep_nhx(:,nlt,nln) = (/onhxav(:,nlt)/)
 end if
 end do
end do
print("assigned to paleo over respective land/ocean pts")

; also test to see if both land and ocean data area missing, if so, fill in
; with zonal average
ndepmiss = ismissing(pndep)
do nlt = 0,nplat-1
do nln = 0,nplon-1
 if(ndepmiss(0,nlt,nln).eq.True)then
  pndep(:,nlt,nln) = (/ndepav(:,nlt)/)
  pndep_noy(:,nlt,nln) = (/noyav(:,nlt)/)
  pndep_nhx(:,nlt,nln) = (/nhxav(:,nlt)/)
 end if
end do
end do
print("check/correction completed for missing where land/ocean inconsistenties")

;=== specific to permian/triassic (pangean),
;  correct over 39S (roaring 40s modern not comparable to 40S pangea)
;  apply some hemispheric symmetric

topoC = conform(pndep,topo,(/1,2/))
topoC!1 = "lat"
topoC&lat = plat

; hemisphere symmetry fix
lndepav!1 = "lat"
lndepav&lat = plat
lnoyav!1 = "lat"
lnoyav&lat = plat
lnhxav!1 = "lat"
lnhxav&lat = plat

do nlt = 0,nplat-1
fliplat = lat(nlt) * -1.
do nln = 0,nplon -1
if(topoC(0,nlt,nln).ge.0.and.lat(nlt).lt.-30.and.lat(nlt).gt.-80.)then
 pndep(:,nlt,nln) = (/lndepav(:,{fliplat})/)
 pndep_noy(:,nlt,nln) = (/lnoyav(:,{fliplat})/)
 pndep_nhx(:,nlt,nln) = (/lnhxav(:,{fliplat})/)
else
if(topoC(0,nlt,nln).gt.0.and.lat(nlt).lt.-30.and.lat(nlt).gt.-80.)then
 pndep(:,nlt,nln) = (/ondepav(:,{fliplat})/)
 pndep_noy(:,nlt,nln) = (/onoyav(:,{fliplat})/)
 pndep_nhx(:,nlt,nln) = (/onhxav(:,{fliplat})/)
end if
end if
end do
end do


; == copy 1849 (first time slice to all timeslices)

do nt = 0,ntimes-1
do nlt = 0,nplat-1
do nln = 0,nplon-1
  pndep(:,nlt,nln) = pndep(0,nlt,nln)
  pndep_noy(:,nlt,nln) = pndep_noy(0,nlt,nln)
  pndep_nhx(:,nlt,nln) = pndep_nhx(0,nlt,nln)
end do
end do
end do


print("Making Pangean (P/T) specific changes")


;=== apply smoother to fire carbons
p = .5
q = .25   ;  -.25 = light local smoothing... +.25 =heavy
cyclic = True
pndep = smth9_Wrap(pndep,p,q,cyclic)
pndep_noy = smth9_Wrap(pndep_noy,p,q,cyclic)
pndep_nhx = smth9_Wrap(pndep_nhx,p,q,cyclic)

pndep!0 = "time"
pndep!1 = "lat"
pndep!2 = "lon"
pndep&time = time
pndep&lat = lat
pndep&lon = lon
pndep@long_name = NDEP_year@long_name
pndep@unitsNDEP_year@units

pndep_noy!0 = "time"
pndep_noy!1 = "lat"
pndep_noy!2 = "lon"
pndep_noy&time = time
pndep_noy&lat = lat
pndep_noy&lon = lon
pndep_noy@long_name = NDEP_NOy_year@long_name
pndep_noy@unitsNDEP_NOy_year@units

pndep_nhx!0 = "time"
pndep_nhx!1 = "lat"
pndep_nhx!2 = "lon"
pndep_nhx&time = time
pndep_nhx&lat = lat
pndep_nhx&lon = lon
pndep_nhx@long_name = NDEP_NHx_year@long_name
pndep_nhx@unitsNDEP_NHx_year@units

;========create new file surface file

;----> write new .nc file

system("rm "+new_cdf)
cdf_nd = addfile(new_cdf, "c")
print("Creating new ndep file...")

; define dimensions
dimNames = (/"lat","lon","time"/)
dimSizes = (/ nlat, nlon,ntimes/)
dimUnlim = (/ False, False,True/)
filedimdef(cdf_nd,dimNames,dimSizes,dimUnlim)

;define variable names
filevardef(cdf_nd,"lat","double",(/"lat"/))
filevardef(cdf_nd,"lon","double",(/"lon"/))
filevardef(cdf_nd,"time","double",(/"time"/))
filevardef(cdf_nd,"YEAR","integer",(/"time"/))
filevardef(cdf_nd,"NDEP_year","float",(/"time","lat","lon"/))
filevardef(cdf_nd,"NDEP_NOy_year","float",(/"time","lat","lon"/))
filevardef(cdf_nd,"NDEP_NHx_year","float",(/"time","lat","lon"/))


; write variables
cdf_nd->lat    = lat
cdf_nd->lon    = lon
cdf_nd->time   = time
cdf_nd->YEAR   = YEAR
cdf_nd->NDEP_year=  pndep
cdf_nd->NDEP_NOy_year=  pndep_noy
cdf_nd->NDEP_NHx_year=  pndep_nhx

; create file attributes
cdf_nd@template = "template = " + f_template_nd
cdf_nd@topo = "topo = " + filename_land
cdf_nd@info = " given 1850 zonal avg (delineated by landmask), hem symmetry" ;edit for you info

print(new_cdf + " has been written")


end