source: trunk/CONFIG/ORCA2_LIM/IDL_scripts/fig_ht_fromQ.pro @ 2151

pro fig_ht_fromQ,  Q1, Q2, sEXP1, sEXP2, start_end, POSTSCRIPT = postscript, MASK_FILENAME = mask_filename, _extra = ex
  
  compile_opt idl2, strictarrsubs

@common

  tmask_save = tmask  
                                ; suppress 1/2 north ligne dans TMASK
  tmask[jpi/2:jpi-1,jpj-1,0] = 0b ; 0 byte

  Qave1 = moyenne(Q1,'xy') 
  Qave2 = moyenne(Q2,'xy')

  Qnet1 = Q1 - Qave1
  Qnet2 = Q2 - Qave2

  index=where(gphit eq max(gphit))
  nx=index[0] mod jpi
  gphi_save=gphit
  gphit[0,*]=gphit[nx,*]
                                ; 
  msk = read_ncdf( 'atlmsk_nomed', filename = mask_filename, _extra = ex)
  
                                ; from 2D array to 3D array 
  msk = msk.arr[*]#replicate(1., nzt)
  
                                ; x average (zonal mean)
                                ; *1.E-15 to have PetaWatt
  Q1x = moyenne(Qnet1*e2t,'x',/integration)*1.E-15
  Q2x = moyenne(Qnet2*e2t,'x',/integration)*1.E-15
  Q1x_atl = moyenne(Qnet1*e2t*msk,'x',/integration)*1.E-15
  Q2x_atl = moyenne(Qnet2*e2t*msk,'x',/integration)*1.E-15
                                ; northward heat flux transport from antartic  
  htr1 = total(Q1x, /cumulative)
  htr2 = total(Q2x, /cumulative)
  htr1_atl = total(reverse(Q1x_atl), /cumulative)
  htr1_atl = -reverse(htr1_atl)
  htr2_atl = total(reverse(Q2x_atl), /cumulative)
  htr2_atl = -reverse(htr2_atl)

  bad_gphit = where(gphit[0,*] lt -30.)
  htr1_atl[bad_gphit] = !Values.F_NaN   ; we take from values north till 30° South
  htr2_atl[bad_gphit] = !Values.F_NaN   ; we take from values north till 30° South

  if KEYWORD_SET(POSTSCRIPT) then begin
     openps, '004_MHT_'+sEXP1+'-'+sEXP2+'_y'+start_end+ '.ps', /landscape
  endif
  
; update data informations
  varname = 'MHT'
  varunit = 'Pw'
  vargrid = 'T' 

  ibce1 = STRTRIM(Qave1, 1)
  ibce2 = STRTRIM(Qave2, 1)

  title = sEXP1+' (Red) MHT (PetaW) - (Blue) Atlantic MHT (PetaW) / year '+start_end
  subtitle = sEXP1+' (Qnet='+ibce1+' W/m2) -  (Qnet='+ibce2+' W/m2)'

  plt1d, htr1,'y',min = -2., max = 2.5, TITLE=title, SUBTITLE = subtitle $ 
        , SMALL = [2,2,1], YTITLE = 'PW',  XGRIDSTYLE = 2, TICKLEN = 1, YGRIDSTYLE = 2,  _extra = ex   
  plt1d, htr1_atl,'y', COLOR = 50, /ov1d,  _extra = ex

 if sEXP1 ne sEXP2 then begin
  
  title = sEXP1+' - '+sEXP2+' MHT (TeraW) differences / year '+start_end
  plt1d, (htr1-htr2)*1.e3,'y',min = 0., max = 500., TITLE = title $
         , SMALL = [2,2,3], TICKLEN = 1, XGRIDSTYLE = 2 $
         , YTITLE = 'PW', YGRIDSTYLE = 2,   /NOERASE, _extra = ex
  plt1d, (htr1_atl - htr2_atl),'y', COLOR = 50, /ov1d,  _extra = ex
 
 endif 

  title = 'Qnet'
  text = [   title, '<img width="80%" src=/Users/sflod/TOOLS/my_IDL/PS/v32colib_v32colib/ps/004_MHT_'+sEXP1+'-'+sEXP2+'_y'+start_end+ '.png  />  ' ]
  putfile, '/tmp/prova.txt', text

  if KEYWORD_SET(POSTSCRIPT) then begin
     closeps
  endif
  
  tmask=tmask_save
  gphit=gphi_save
  
  return
end