source: trunk/procs/macros/make_linfit.pro @ 164

;+
;
; compute point-wise slope of linear regression of two 3D (x,y,t) fields
;
; @version
; $Id$
;
;-
FUNCTION make_linfit,  file_name,  ncdf_db, BOXZOOM = boxzoom,  TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data

@common
@com_eg

   IF debug_w THEN print, '   ENTER make_linfit... '

   var_name1 =  (strsplit(macro_base_fld, ',', ESCAPE = ' ', /EXTRACT))[0]
   var_name2 =  (strsplit(macro_base_fld, ',', ESCAPE = ' ', /EXTRACT))[1]

;  Build file_name2 if different

   IF strpos (cmd1_back.grid, '#') NE -1 THEN BEGIN 
      varpos = strpos(file_name,  var_name1)
      file_name2= strmid(file_name, 0, varpos)+var_name2+'.nc'
      IF debug_w THEN print, '   file_name2 = ', file_name2
   ENDIF ELSE file_name2 = file_name

;  Read the variables in the correspondant netcdf file
   var1 =  nc_read(file_name, var_name1,  ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2)
   var2 =  nc_read(file_name2, var_name2,  ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1,  TIME_2 =  time_2)

   varname =  var1.name+'/'+var2.name+' corr.'
   varunit =  '[-1/1]'

   nxa = (size(var1.data))[1]
   nya = (size(var1.data))[2]

   IF debug_w THEN print, '   nxa, nya', nxa, nya

   pt_linfit = fltarr(nxa, nya)
   pt_err = fltarr(nxa, nya)
   pt_corr = fltarr(nxa, nya)
   pt_linfit[*, *] = 0
   pt_err[*, *] = 1.1
   pt_corr[*, *] = 0.

; Sampling of data and computation of new numbers of values
   @mth_decode

   FOR imth = 0, nmth-1 DO BEGIN 

      IF debug_w THEN print, '   month idx/value: ', imth, strd(imth)

      data1 =  (var1.data)[*, *, reform(idxm(imth,*), njpt)]        
      data2 =  (var2.data)[*, *, reform(idxm(imth,*), njpt)]        
      nval =  njpt
   
; Compute linear regression


      FOR idx = 0, nxa-1 DO FOR idy = 0, nya-1 DO BEGIN

         IF data1(idx, idy, 0) LT valmask/10. THEN BEGIN 

            x1i = reform(data1(idx, idy, *), nval)
            x2i = reform(data2(idx, idy, *), nval)
            x1 = x1i-mean(x1i)
            x2 = x2i-mean(x2i)

            no_compute = 0

            IF linfit_map NE '' THEN BEGIN
            ; compute regression for value above (p - default) or below (m) linfit_sep
               CASE linfit_map OF
                  'm': idt = where (x2 LE linfit_sep)
                  ELSE: idt = where (x2 GE linfit_sep)
               ENDCASE 
               IF (size(idt))[1] GE 3  THEN BEGIN 
                  IF idt[0] NE -1 THEN BEGIN 
                     x1 = x1(idt)
                     x2 = x2(idt)
                  ENDIF 
               ENDIF ELSE BEGIN 
                  no_compute = 1
               ENDELSE 
            ENDIF 
            
;           BEGIN 
;             print, '  idx, idy, linfit_map ', idx, idy, linfit_map
;             print, '  size(idt)',size(idt)
;             stop
;          ENDIF 
         
            IF no_compute EQ 0 THEN BEGIN 

               coeffl = linfit(x2, x1, CHISQ = linerrl, PROB = proberrl, SIGMA = sigmaerrl)
               correl = c_timecorrelate(x2,x1)
               pt_linfit(idx, idy) = pt_linfit(idx, idy)+coeffl(1)
               pt_err(idx, idy) = min(pt_err(idx, idy), proberrl)
               pt_corr(idx, idy) = pt_corr(idx, idy) + correl
               
               IF debug_w THEN BEGIN  
                  IF idx EQ 30 AND idy EQ 15 THEN BEGIN 
                     print, '  idx, idy, linfit_map ', idx, idy, linfit_map
                     print, '  size(idt)',size(idt)
                     print, '  pt_linfit, correl = ',coeffl(1), correl
                  ENDIF 
               ENDIF 
            ENDIF ELSE BEGIN 
               pt_linfit(idx, idy) = valmask
               pt_err(idx, idy) = 0.
               pt_corr(idx, idy) = 0.
            ENDELSE 

         ENDIF ELSE BEGIN 
            pt_linfit(idx, idy) = valmask
            pt_err(idx, idy) = 0.
            pt_corr(idx, idy) = 0.
         ENDELSE 

      ENDFOR 
   ENDFOR 

  ; make mean of period

   idm = where(pt_linfit GE valmask/10.)
   pt_linfit = pt_linfit/float(nmth)
   pt_corr =  pt_corr/float(nmth)
   IF idm[0] NE -1 THEN pt_linfit(idm) = valmask

  ; only take points where correlation larger than 0.2

   idc = where(abs(pt_corr) LT corr_thres)
   IF idc[0] NE -1 THEN pt_linfit(idc) = !values.f_nan

   corr_txt = '[r>'+strtrim(string(corr_thres, format = '(f5.2)'), 2)+']'

   varname =  varname+' '+ntxt

   print,  '    Linfit min/max pt_corr ', min(pt_corr), max(pt_corr)
   
;  Define the outputs of the function   
   field = {name: varname, data: pt_linfit, legend: '', units: '', origin: '', dim: 0,  direc:'xy'}   
   field.origin = var1.origin
   field.dim = var1.dim

   CASE linfit_map OF
      'p': BEGIN
         print,  '    Linear fit computed for anomalies ABOVE ',linfit_sep 
         lintxt = ' > '+ strtrim(string(linfit_sep), 2)
      END 
      'm': BEGIN 
         print,  '    Linear fit computed for anomalies BELOW ',linfit_sep 
         lintxt = ' < '+ strtrim(string(linfit_sep), 2)
      END 
      ELSE: lintxt = ''
   ENDCASE 

   field.legend = lintxt+' for '+ntxt+' in ['+cmdm.date1+'-'+cmdm.spec+'] '+corr_txt+' -'

; additional computations (pac_5 nino_3 nino_4 averages) *** requires whole domain ****

   IF nxt EQ jpi AND nyt EQ jpj THEN BEGIN 
      
      n3_linfit = moyenne(pt_linfit,'xy', boite = [210,  270,  -5, 5], NaN = valmask)
      n4_linfit = moyenne(pt_linfit,'xy', boite = [160,  210,  -5, 5], NaN = valmask)
      zl_linfit = moyenne(pt_linfit,'xy', boite = [130,  280,  -5, 5], NaN = valmask)
      
      print, '    Nino 3 average of slope of linear fit = ', n3_linfit
      print, '    Nino 4 average of slope of linear fit = ', n4_linfit
      print, '    Zonal 5N/5S average of slope of linear fit = ', zl_linfit

   ENDIF 

 ;  jpt=nval
 ;  meants=grossemoyenne(data2, 't')
 ;  meantsp=reform(meants,nxt*nyt)  
 ;  meanlf=reform(pt_linfit,nxt*nyt)
 ;  IF min (meantsp) GE 100 THEN meantsp=meantsp-273.15
 ;  pltsc, meanlf,meantsp,-70,30,20,34,"SST "+mth[strd-1], window=2  
 ;  pltsc, meanlf,meantsp,-70,30,20,34,"SST "+mth[strd-1], window=2, /noerase, /ov1d, col1d=2

 ;  stop

   IF debug_w THEN print, '   ... EXIT make_linfit'
   
   return, field

END