#!/bin/bash
#	
#	SCRIPTS USED:
#	ROUTINES USED:
#	GMT PROGRAMS USED:
#
#   Example script to prepare the data and force grids for the itration flexure program
#   files in Mercator coordinates will be called _xy.grd
#   in mercator coordinates
#
#   first set the boundaries of the area
# 

	PARFILE=$1
    source $PARFILE
    
 	RLBOUNDS="-R${WLON}/${ELON}/${SLAT}/${NLAT}"   
    
#	SWITCH GMT VERSION MANUALLY
	gmtswitch /Users/esg006/gmt5    
    gmtset PROJ_ELLIPSOID = Sphere

#
#   get the gravity and topography grids using img2grd
#   and the age grid using grdcut
#
	rm grav.img.23.1 topo_18.1.img
	ln -s /radar3/flex/sup_data/global_grav_1_min/grav.img.23.1 .
	ln -s /radar3/flex/sup_data/global_topo_1min/topo_18.1.img .
#	ln -s /radar3/flex/sup_data/srtm15_plus/topo15.grd .

	gmt img2grd grav.img.23.1 -T1 -M -I1m -D -S.1 $RLBOUNDS -Ggrav_merc.grd -V
	gmt img2grd topo_18.1.img -T2 -M -I1m -D $RLBOUNDS -Gtopo_merc.grd -V 
	gmt img2grd topo_18.1.img -T1 -M -I1m -D $RLBOUNDS -Gtopo_all.grd -V 
	
	gmt grdcut $RLBOUNDS topo15.grd -Gsrtm_topo.grd
	
	
	rm age.cpt age.grd
	ln -s /radar3/flex/age/age.cpt
#	ln -s /radar3/flex/age/age.3.6.nc age.grd
	cp /radar3/flex/age/age.3.6.nc age.grd
#
# 	generate mask for age grid based on trench axis planform
#

	if [ -n "$TRSYS" ]
	then
		echo $TRSYS
		prep_age_mask.com $PARFILE $TRSYS $MASKDIR
		
		gmt mapproject age_mask_ll.xy $RLBOUNDS -Jm1c --PROJ_ELLIPSOID=Sphere > age_mask_merc.xy
		RMBOUNDS=`gmt grdinfo -I- grav_merc.grd`
	
 	 	gmt grdmask $RMBOUNDS -I1m/1m age_mask_merc.xy -N1/1/NaN -r -Gage_mask_merc.grd
	fi
	
#
# 	OR generate mask based on boundaries of plate
#

	if [ -n "$PLMASK" ]
	then
		ln -s /radar3/flex/age/age_mask/${PLMASK} ${PLMASK}
		gmt gmtspatial ${PLMASK} $RLBOUNDS -T > age_mask_ll.xy
		
		gmt mapproject age_mask_ll.xy $RLBOUNDS -Jm1c --PROJ_ELLIPSOID=Sphere > age_mask_merc.xy
		RMBOUNDS=`gmt grdinfo -I- grav_merc.grd`
 	 	gmt grdmask $RMBOUNDS -I1m/1m age_mask_merc.xy -NNaN/1/1 -r -Gage_mask_merc.grd
	fi


 	gmt grdmath age.grd 100 DIV = age_scaled.grd 	
	gmt grdcut age_scaled.grd $RLBOUNDS -Gage_cut.grd
	gmt grdsample age_cut.grd -I1m/1m -nl -Gage.grd
	gmt grdproject age.grd -Jm1c -D1m/1m $RLBOUNDS -r -nl -V --PROJ_ELLIPSOID=Sphere -Gage_proj.grd
	# filter
	gmt grdedit age_proj.grd `grdinfo age_mask_merc.grd -I-`
	gmt grdfilter age_proj.grd -D5 -Fg60 -Gage_f1.grd -V

	# reset ages
	gmt grdmath age_proj.grd age_f1.grd AND = age_f1_reset.grd 
	# filter again
	gmt grdfilter age_f1_reset.grd -D5 -Fg30 -Gage_f2.grd -V
	# reset again
	gmt grdmath age_proj.grd age_f2.grd AND = age_f2_reset.grd
	gmt grdmath age_f2_reset.grd age_mask_merc.grd MUL = age_masked.grd
	gmt grdclip age_masked.grd -Sb1.0/1.0 -Gage_clipped.grd
	gmt grdmath age_clipped.grd 1.0 AND = age_merc.grd

#
#	ESTIMATE OBSERVED CURVATURE FROM TOPO15
#


#	gmt grdmath topo_all.grd CURV 9.9022E-11 MUL = curv_init.grd
#	gmt grdmath curv_init.grd age_mask_merc.grd MUL = curv_masked.grd
#	gmt grdfilter curv_masked.grd -Fg100 -D3 -V -Gcurv_filt.grd
#	gmt grdmath curv_filt.grd age_mask_merc.grd MUL = curv_masked2.grd

#	gmt grdmath topo_all.grd age_mask_merc.grd MUL = topo_init.grd
#	gmt grdfilter topo_init.grd -Fg100 -D3 -V -Gtopo_filt.grd
#	gmt grdmath topo_filt.grd age_mask_merc.grd MUL = topo_masked.grd
#	gmt grdmath topo_masked.grd CURV 9.9022E-11 MUL = curv_init.grd
#	gmt grdfilter curv_init.grd -Fg30 -D3 -V -Gcurv_filt.grd
#	gmt grdmath curv_filt.grd age_mask_merc.grd MUL = curv_masked.grd

#	gmt dimfilter srtm_topo.grd -Fg30 -D3 -Nm8 -V -Gtopo_dim.grd
#	gmt grdfilter srtm_topo.grd -Fg30 -D3 -V -Gtopo_filt.grd
#	gmt grdmath topo_dim.grd CURV 9.9022E-11 MUL = curv_init.grd
# 	gmt grdfilter curv_init.grd -I1m/1m -Fg100 -D3 -V -Gcurv_filt.grd
# 	gmt grdproject curv_filt.grd -Jm1c -D1m/1m $RLBOUNDS -r -nl -V --PROJ_ELLIPSOID=Sphere -Gcurv_proj.grd
#	gmt grdmath curv_proj.grd age_mask_merc.grd MUL = curv_masked.grd
#	gmt grdmath curv_masked.grd -1.0E-9 AND = curv0_xy.grd
#	gmt grdedit curv0_xy.grd `grdinfo -I- fz_xy.grd`
   
	prep_trench_segs.com ${TSEGSRC} ${NSEGSINC} ${REGNAME}
	gmt mapproject ${TSEGSRC} $RLBOUNDS -Jm1c --PROJ_ELLIPSOID=Sphere > tsegs_merc.xy
	gmt mapproject ${TSGISRC} $RLBOUNDS -Jm1c --PROJ_ELLIPSOID=Sphere > tsegs_int_merc.xy
# 	    
  	awk '{print $1,$2,$3,$4," 1.00e14 1.00e17"}' $ALLSEGS > load.dat
	
#   make a load for the trench model and force the grid into the gravity coordinate system

 	rm fz_xy.grd
 	vload $WLON $SLAT load.dat fz_xy.grd
    cp fz_xy.grd fz_merc.grd
	gmt grdedit fz_merc.grd `grdinfo -I- grav_merc.grd`
# 
#   keep the fz.grd file in the Mercator space and the 
#   fz_xy.grd in the xy space
# 
	
	gmt grdmath fz_xy.grd 0.0 MUL -1.0e-9 ADD = curv0_xy.grd 


#	grdmath fz_xy.grd 0.0 MUL -3.5e-8 ADD = curv0_xy.grd 
#	grdmath fz_xy.grd 0.0 MUL -1.0e-7 ADD = curv0_xy.grd 
# 
	cp age_merc.grd age_xy.grd
	gmt grdedit age_xy.grd `grdinfo -I- fz_xy.grd`
 
# 	age2rig $AGEGRD $RIGINIT 
 
#	REPROJECT TRENCH DISTANCE GRIDS

	echo $RLBOUNDS
	rm tdist_full.grd
	ln -s /radar3/flex/trench_distance_grids/${TRENAME}.grd tdist_full.grd
	gmt grdcut tdist_full.grd $RLBOUNDS -Gtdist_cut.grd
	gmt grdsample tdist_cut.grd -I1m/1m -nl -Gtdist.grd
	gmt grdproject tdist.grd -Jm1c -D1m/1m $RLBOUNDS -nl -V --PROJ_ELLIPSOID=Sphere -Gtdist_proj.grd
	gmt grdmath tdist_proj.grd -1000.0 AND -1.0 MUL = tdist_adj.grd
	gmt grdclip tdist_adj.grd -Sb1.0/1.0 -Gtdist_wgt.grd
	
#	GENERATE MASKS

	rm *ext.xy *int.xy 

	awk '{print $1, $2}' < b_box > ${REGNAME}_ext.xy
	awk '{print $1, $2}' < c_box > ${REGNAME}_int.xy
	
	prep_mask.com $PARFILE