source: trunk/src/flux_evaluation_tpr_timeseries.pro @ 106

;------------------------------------------------------------
pro flux_evaluation_tpr_timeseries, $
;                    var,    $ ;;  flux variable (swr, lwr, lhf, shf) to calculate the statistics
                    date1,  $ ;;  start date (in julian date. eg. 20000101)
                    date2     ;;  end date (in julian date eg. 20091231)

@common
;------------------------------------------------------------
reinitplt, /z,/invert
key_portrait = 0
openps, FILENAME = 'idl.ps'
;------------------------------------------------------------
;; part to change

min_obs=10.  ;; this will allow to calculate statistics at locations with more than 180. valid observation

;; choose the appropriate min and max values for the following.

bias_mi=-5      &  bias_ma=5     & bias_int=0.5
std_mi=0.5      &  std_ma=1.5    & std_int=0.05
rmsd_mi=0       &  rmsd_ma=15    & rmsd_int=1.5
cor_mi=0.5      &  cor_ma=1.     & cor_int=0.025

;------------------------------------------------------------
;; Before running this program, you have to compile the following subroutines
;; .r extract_tpr_location
;; .r read_tpr_netflux

;; TPR locations.  This needs to be updated with time since more locations are added to the array.

sitelist=['5n165e','8s67e','12s55e', '8s55e', '8s80.5e', '1.5s80.5e', '0n80.5e', '1.5n80.5e', '1.5s90e', $
           '0n90e', '1.5n90e', '4n90e','8n90e','12n90e', '15n90e', '5s95e', $
           '8s165e', '8s180w',  '8s155w', '8s125w', '8s110w', '8s95w',  '5s156e', '5s165e', '5s180w', '5s170w', $
          '5s155w', '5s140w', '5s125w', '5s110w', '5s95w', '2s156e', '2s165e', '2s180w', '2s170w', '2s155w', '2s140w', $
          '2s125w', '2s110w', '2s95w', '0n147e', '0n156e', '0n165e', '0n180w', '0n170w', '0n155w', '0n140w', '0n125w', $
          '0n110w', '0n95w', '2n147e', '2n156e', '2n165e', '2n180w', '2n170w', '2n155w', '2n140w', '2n125w', '2n110w', $
          '2n95w', '5n147e', '5n156e', '5n170w', '5n155w', '5n140w', '5n125w', '5n110w', '5n95w', $
          '8n156e', '8n165e', '8n180w', '8n170w', '9n140w', '8n125w', '8n110w', '8n95w', $
          '0n0e', '0n10w', '0n23w', '0n35w', '10s10w', '12n23w', '12n38w', '14s32w', '15n38w', '19s34w', '20n38w', $
          '21n23w', '4n23w', '4n38w', '6s10w', '8n38w', '8s30w']

;------------------------------------------------------------------------------------------------------------------------
;;   This program will create the following text files with statistics of respective variables
;------------------------------------------------------------------------------------------------------------------------
close,/all

fi='/Users/pkb/work/MY_SAXO/TropFlux_update/flux_stat.txt'
openw,1,fi

printf,1, 'x     y      cor    bias     std ratio     rmsd'
;------------------------------------------------------------------------------------------------------------------------
;; First, this program reads the full TropFlux data and later extract it at specific TPR locations

file="/Volumes/PAYASAM/TropFlux/TropFlux/shf_tropflux_1d_1989_2010.nc"
initncdf, file
var=-1*read_ncdf("shf", date1, date2, file=file,/nostr)
help, var

;------------------------------------------------------------------------------------------------------------------------
nn=n_elements(sitelist)
date1=date1
date2=date2
nsmooth=1. ;; this will return daily TPR values
mooring=0 & product=0

for n=0, nn-1 do begin

;; reading data from mooring

    site=sitelist(n) & csite=site
    print, csite
    x=x_site_location(site)
    y=y_site_location(site)
    if (y ge 0. and y le 30.) then y=y+360.
    dx=0.5 & dy=0.5 & box=[y-dy, y+dy, x-dx, x+dx]
     
    read_tpr_netflux, csite,date1,date2,nsmooth, $
            sw,lw,sh,lh

;; select the appropriate variables for evaluation (trp = sw or lw or sh or lh)

   tpr=sh

   ind=where(finite(tpr)) & no_valid=n_elements(ind)

   if (no_valid ge min_obs) then begin


           extract_tpr_location,var,box, $
                var_tpr
           var_tpr=reform(var_tpr)

           stats_5d, tpr,var_tpr, $  ;; tpr=TPR observation and var_tpr=gridded product extracted at TPR location
              cor, bias, std, rmsd

           printf, 1, x, y, cor, bias, std, rmsd, format='(f6.2, 3x, f6.2, 3x, f4.2, 3x, f7.2, 3x, f4.2, 3x, f5.2)'
           cstat=string(cor, bias, std, rmsd, format='(f4.2,3x,f7.2,3x,f4.2,3x,f5.2)')
           print, cstat
 
;;         PLOTTING THE TIME-SERIES
           array=[tpr, var_tpr] & mi=min(array,/nan) & ma=max(array,/nan) & int=(ma-mi)/3.
           pltt, ts_smooth(tpr,5,/nan), "t",/rempl, small=[1,3,1], lct=65, $
                 title='Five day stats are shown below.  TPR (black) and Product (red) at'+csite+' ', charsize=1., $
                 subtitle=cstat
           ind=where(finite(tpr,/nan)) & var_tpr(ind)=!Values.f_nan
           pltt, ts_smooth(var_tpr,5,/nan), "t",/ov1d, color=250
           erase
           mooring=[mooring,tpr] & product=[product,var_tpr]
    endif
endfor

close,/all
;----------------------------------------------------------
closeps

fig='flux_evaluation_tpr_timeseries.ps'
spawn, 'mv '+psdir+'idl.ps '+updatedir+fig
spawn, 'gv '+updatedir+fig
return
end
;--------------------------------------------------------------------------
function x_site_location, site
    n1=strpos(site, 's')
if (n1 gt -1) then begin
    ns=-1.
    x=strmid(site, 0, n1)
    x=float(x)*ns
endif else begin
    n1=strpos(site, 'n')
    x=strmid(site, 0, n1)
    ny=1.
    x=float(x)*ny
endelse
return, float(x)
end
;--------------------------------------------------------------------------
function y_site_location, site
    n1=strpos(site, 'e')
if (n1 gt -1) then begin
    n=strpos(site, 's')
    if (n gt -1) then begin
        y=strmid(site, n+1, n1-n-1)
    endif else begin
        n=strpos(site, 'n')
        y=strmid(site, n+1, n1-n-1)
    endelse

endif else begin
    n1=strpos(site, 'w')
    n=strpos(site, 's')
    if (n gt -1) then begin
        y=strmid(site, n+1, n1-n-1)
        y=180+(180-float(y))
    endif else begin
        n=strpos(site, 'n')
        y=strmid(site, n+1, n1-n-1)
        y=180+(180-float(y))
   endelse
endelse
return,float(y)
end

;--------------------------------------------------------------------------