pro std_plots, std_EXP1, std_EXP2, ystart, _extra=ex
 
 ; scripts for nemo v3_3


PRINT, 'usage: postscript '
PRINT, '========================================================================'
PRINT, 'std_plots, ''std_EXP1'', ''std_EXP2'', ystart, /postscript'
PRINT, '========================================================================'
PRINT, 'std_plots, ''core2000'', ''core2000'', 1991, /postscript'
PRINT, '========================================================================'


  compile_opt idl2, strictarrsubs

  PRINT, ''
  PRINT, ' 	############################################'
  PRINT, ''
  PRINT, '                    LAUNCH of std_plots'
  PRINT, ''
  PRINT, ' 	############################################'
  PRINT, ''


@common
@initorca2
@std_com

  year_start = ystart * 10000L + 101L                                                 ; year of start (corresponding to filename)
  year_end = (ystart + 9L) * 10000L + 1231L                                           ; year of end (corresponding to filename)
  start_end = string(year_start,format= '(I08)')+'_'+string(year_end,format ='(I08)') ; variable string for filename
;  std_iodir_data    = '/Users/sflod/idl_PLOTS/DATA_STORE/RUN_INTERAN/brodie/'        ; path of data in NetCdf format
  std_iodir_data    = '/Users/sflod/idl_PLOTS/DATA_STORE/RUN_CLIMATO/brodie/'        ; path of data in NetCdf format
  std_iodir_climato = '/Users/sflod/idl_PLOTS/CLIMATOLOGIES/'                         ; path of climatological data
  std_iodir_mask    = '/Users/sflod/idl_PLOTS/MASK/'                                  ; path of mask files (ex: subbasins)
  
;;;;; ==================================================== WARNING ================================================================= ;;;;;;;;;;
  d1 = (ystart + 9L) * 10000L + 101L ;  date to select last record of file with 10 years of record
  d2 = year_end  ; date to select last record of file with 10 years of record
  id1 = (ystart + 9L) * 10000L + 1201L ; date to select starting date for (monthly) ice output 
  id2 = d2 ; date to select last record for (monthly) ice output
  id_mars_1 = (ystart + 9L) * 10000L + 301L ; date to select mars ice output 
  id_mars_2 = (ystart + 9L) * 10000L + 331L  ; date to select mars ice output
  id_sept_1 = (ystart + 9L) * 10000L + 901L ; date to select september ice output 
  id_sept_2 = (ystart + 9L) * 10000L + 930L ; date to select september ice output
;
; Levitus 98
  std_file_Levitus_T =  'data_1y_potential_temperature_nomask.nc'
  std_file_Levitus_S =  'data_1y_salinity_nomask.nc'

; Output run experience1
  std_file1_T     = std_EXP1+'_'+start_end+'_1Y_grid_T.nc' ; grid T input file
  std_file1_U     = std_EXP1+'_'+start_end+'_1Y_grid_U.nc' ; grid U input file
  std_file1_V     = std_EXP1+'_'+start_end+'_1Y_grid_V.nc' ; grid V input file
;  std_file1_W    = std_EXP1+'_'+start_end+'_1Y_grid_W.nc' ; grid W input file
  std_file1_I     = std_EXP1+'_'+start_end+'_1M_icemod.nc' ; ice    input file
  
; Output run experience2
  std_file2_T     = std_EXP2+'_'+start_end+'_1Y_grid_T.nc' ; grid T input file
  std_file2_U     = std_EXP2+'_'+start_end+'_1Y_grid_U.nc' ; grid U input file
  std_file2_V     = std_EXP2+'_'+start_end+'_1Y_grid_V.nc' ; grid V input file
;  std_file2_W    = std_EXP2+'_'+start_end+'_1Y_grid_W.nc' ; grid W
  std_file2_I     = std_EXP2+'_'+start_end+'_1M_icemod.nc' ; ice    input file


  PRINT, ''
  PRINT, '	std_EXP1       : ', std_EXP1
  PRINT, '	std_EXP2       : ', std_EXP2
  PRINT, '	std_iodir_data : ', std_iodir_data 
  PRINT, '	std_file1T : ', std_file1_T
  PRINT, '	std_file1U : ', std_file1_U
  PRINT, '	std_file1V : ', std_file1_V
;  PRINT, '	std_file1W : ', std_file1_W
  PRINT, '	std_file2I : ', std_file1_I
  PRINT, '	std_file2T : ', std_file2_T
  PRINT, '	std_file2U : ', std_file2_U
  PRINT, '	std_file2V : ', std_file2_V
;  PRINT, '	std_file2W : ', std_file2_W
  PRINT, '	std_file2I : ', std_file2_I
  PRINT, ''

;#########################################################################
; reading variables

;;; 3D ;;;
; temperature
  T1   = read_ncdf('thetao', d1, d2, filename = std_iodir_data+std_file1_T                          ) 
  T2   = read_ncdf('thetao', d1, d2, filename = std_iodir_data+std_file2_T                          ) 
  TLev = read_ncdf('votemper', filename = std_iodir_climato+std_file_Levitus_T                ) 
  TRey = read_ncdf('sst'     , filename = std_iodir_climato+'NewREY_ORCA2_1991_2000_1y.nc'    ) 

help,T1,T2,TLev,TRey

; salinity
  S1   = read_ncdf('so', d1, d2, filename = std_iodir_data+std_file1_T             ) 
  S2   = read_ncdf('so', d1, d2, filename = std_iodir_data+std_file2_T             ) 
  SLev = read_ncdf('vosaline', filename = std_iodir_climato+std_file_Levitus_S      ) 

; Net Downward heat flux
  Q1   = read_ncdf('qt', d1, d2, filename = std_iodir_data+std_file1_T               ) 
  Q2   = read_ncdf('qt', d1, d2, filename = std_iodir_data+std_file2_T               )
  geo  = read_ncdf('heatflow', filename=  std_iodir_climato+'geothermal_heating.nc'   )
  geo  = geo.arr*1.e-3      ; convert into W/m2

;climatology 
  QNET  = read_ncdf('qnet', filename = std_iodir_climato+'OAFlux_1my_01_12_1984_2004_orca2_qnet.nc'   )

; erp (evaporation damping)
  ERP1   = read_ncdf('wfcorr', d1, d2, filename = std_iodir_data+std_file1_T           )
  ERP2   = read_ncdf('wfcorr', d1, d2, filename = std_iodir_data+std_file2_T           )
  ERP1   = ERP1.arr * 86400
  ERP2   = ERP2.arr * 86400

; emp (evaporation minus precipitation)
  EMP1   = read_ncdf('wfo', d1, d2, filename = std_iodir_data+std_file1_T           )
  EMP2   = read_ncdf('wfo', d1, d2, filename = std_iodir_data+std_file2_T           )
  EMP1   = EMP1.arr * 86400
  EMP2   = EMP2.arr * 86400
 
;mixed layer depth
  MXL1 = read_ncdf('mldr10_1', d1, d2, filename = std_iodir_data+std_file1_T                 )  ; 10 m
  MXL2 = read_ncdf('mldr10_1', d1, d2, filename = std_iodir_data+std_file2_T                 )  ; 10 m
;climatology 
  MLD  = read_ncdf('mld', filename = std_iodir_climato+'mld_DR003_c1m_ORCA2_1y.nc'   )

; velocities
  U1   = read_ncdf('uocetr_eff', d1, d2, filename = std_iodir_data+std_file1_U             )
  U2   = read_ncdf('uocetr_eff', d1, d2, filename = std_iodir_data+std_file2_U             )
  e23 = e3u_3d(/e2)
  e23 = e23[*]#replicate(1., jpt)
  U1.arr = U1.arr /e23
  U2.arr = U2.arr /temporary(e23)
  U1.unit = 'm/s'
  U2.unit = 'm/s'
;
  V1   = read_ncdf('vocetr_eff', d1, d2, filename = std_iodir_data+std_file1_V             )
  V2   = read_ncdf('vocetr_eff', d1, d2, filename = std_iodir_data+std_file2_V             )
  e13 = e3u_3d(/e1)
  e13 = e13[*]#replicate(1., jpt)
  V1.arr = V1.arr /e13
  V2.arr = V2.arr /temporary(e13)
  V1.unit = 'm/s'
  V2.unit = 'm/s'

; Take into account the Gent & McWilliams "Eddy Induced velocity"
;  Veiv1 = read_ncdf('vomeeivv', filename = std_iodir_data+std_file2_V )
;  Veiv2 = read_ncdf('vomeeivv', filename = std_iodir_data+std_file2_V )

; ice
  Ithi_mars_1   = read_ncdf('sit', id_mars_1, id_mars_2, filename = std_iodir_data+std_file1_I ) 
  Ithi_mars_2   = read_ncdf('sit', id_mars_1, id_mars_2, filename = std_iodir_data+std_file2_I )

  Ithi_sept_1   = read_ncdf('sit', id_sept_1, id_sept_2, filename = std_iodir_data+std_file1_I ) 
  Ithi_sept_2   = read_ncdf('sit', id_sept_1, id_sept_2, filename = std_iodir_data+std_file2_I )

  Ifra_mars_1   = read_ncdf('sic', id_mars_1, id_mars_2, filename = std_iodir_data+std_file1_I )
  Ifra_mars_2   = read_ncdf('sic', id_mars_1, id_mars_2, filename = std_iodir_data+std_file2_I )

  Ifra_sept_1   = read_ncdf('sic', id_sept_1, id_sept_2, filename = std_iodir_data+std_file1_I )
  Ifra_sept_2   = read_ncdf('sic', id_sept_1, id_sept_2, filename = std_iodir_data+std_file2_I )

 
;#########################################################################

;#########################################################################
;######################  STANDARD PLOTS   ################################
;#########################################################################

; fixed color tabled
  lct,64

  PRINT, '  ' 
  PRINT, '				Experience treated:', std_EXP1
  PRINT, '  ' 
  PRINT, '  ' 
  PRINT, '  ' 

       ;;;;;;;; ErP & Emp salinity damping term
       ;;;;;;;; =============================================
  PRINT, '============================================= '
  PRINT, ' PLOTS 2D '
  PRINT, '============================================= '
  PRINT, '                                ErP salinity damping term'
  PRINT, '1'
  xxx_Erp, ERP1, ERP2, std_EXP1, std_EXP2, start_end, _extra = ex
  PRINT, '2'
  xxx_Emp, EMP1, EMP2, std_EXP1, std_EXP2, start_end, _extra = ex
 
       ;;;;; Net heat flux
        ;;;;; =============
  PRINT, ' '
  PRINT, '                               Net heat flux'
  PRINT, '3 '
  xxx_Qnet, Q1, Q2, QNET, std_EXP1, std_EXP2, start_end, _extra = ex

       ;;; Meridionnal Heat Transport
        ;;; ===========================
  PRINT, ' '
  PRINT, '                               Meridionnal Heat Transport'
  PRINT, '4 '
  fig_ht_fromQ, Q1.arr+geo, Q2.arr+geo, std_EXP1, std_EXP2, start_end, MASK_FILENAME = std_iodir_mask+'subbasins_orca21_nored.nc', _extra = ex

        ;;; Global Barotropic Function
        ;;; ===========================
  PRINT, ' '
  PRINT, ' 				Global Barotropic Function ' 
  PRINT, '5 '
  plt_bsf, U1, U2, std_EXP1, std_EXP2, start_end, _extra=ex

  ;; mean Temperature diff with New Reynolds
  ;; =======================================
  PRINT, ' '
  PRINT, ' 				mean Temperature diff with New Reynolds'
  PRINT, '6 '
  xxx_SST, T1[*,*,0], T2[*,*,0], TRey, std_EXP1, std_EXP2, start_end, _extra=ex


        ;;;; mean Salinity diff with Levitus
        ;;;; ===============================
  PRINT, ' '
  PRINT, '                               mean Salinity diff with Levitus'
  PRINT, '7 '
  xxx_SSS, S1[*,*,0], S2[*,*,0], SLev[*,*,0], std_EXP1, std_EXP2, start_end, _extra=ex

      ;;;;;;; Arctic mean Salinity diff with Levitus at z=105 meters
        ;;;;;;; ======================================================
  PRINT, ' '
  PRINT, '                               Arctic mean Salinity diff with Levitus at z=105 meters'
  PRINT, '8 exp1-levitus'
  xxx_EXP1_S_ARC, S1, SLev, std_EXP1, std_EXP2, start_end, _extra = ex       
  PRINT, '9 exp2-levitus'
  xxx_EXP2_S_ARC, S1, S2, SLev, std_EXP1, std_EXP2, start_end, _extra = ex
  PRINT, '10 exp1-levitus at 105 meters'
  xxx_EXP1_S_ARC_z105, S1, SLev, std_EXP1, std_EXP2, start_end, _extra = ex       
  PRINT, '11 exp2-levitus at 105 meters'
  xxx_EXP2_S_ARC_z105, S1, S2, SLev, std_EXP1, std_EXP2, start_end, _extra = ex


        ;;;;;;;; mean Temperature diff with Levitus at z=105 meters
        ;;;;;;;; ==================================================
  PRINT, ' '
  PRINT, '                               mean Temperature diff with Levitus at z=105 meters'
  PRINT, '12'
  xxx_T_z105, T1, T2, TRey, std_EXP1, std_EXP2, start_end, _extra=ex
                              

        ;;;;;;;; mean Salinity diff with Levitus at z=105 meters
        ;;;;;;;; ==================================================
  PRINT, ' '
  PRINT, '                               mean Salinity diff with Levitus at z=105 meters'
  PRINT, '13'
  xxx_S_z105, S1, S2, SLev, std_EXP1, std_EXP2, start_end, _extra=ex
                            

  PRINT, '============================================= '
  PRINT, ' PLOTS 3D '
  PRINT, '============================================= '
 
  PRINT, '                         Mixed layer depth'
  PRINT, '14 '
  xxx_mxl10, MXL1, MXL2, MLD, std_EXP1, std_EXP2, start_end,  _extra = ex
  PRINT, '15'
  xxx_ZonMxl, MXL1, MXL2, MLD, std_EXP1, std_EXP2, start_end, _extra = ex
 

        ;;;;;;;;; Zonal mean Temperature diff with Levitus
        ;;;;;;;;; ========================================
  PRINT, ' '
  PRINT, '                               Zonal mean Temperature diff with Levitus'
  PRINT, '16 '
  xxx_zonal_T, T1, T2, TLev, std_EXP1, std_EXP2, start_end, MASK_FILENAME = std_iodir_mask+'subbasins_orca21_nored.nc', _extra = ex

  PRINT, '16 '                      
  xxx_zonal_T, T1, T2, TLev, std_EXP1, std_EXP2, start_end, SUBBASIN = 'Atl', MASK_FILENAME = std_iodir_mask+'subbasins_orca21_nored.nc', _extra = ex
  
  PRINT, '16 '
  xxx_zonal_T, T1, T2, TLev, std_EXP1, std_EXP2, start_end, SUBBASIN = 'Ind', MASK_FILENAME = std_iodir_mask+'subbasins_orca21_nored.nc', _extra = ex
  
  PRINT, '16 '
  xxx_zonal_T, T1, T2, TLev, std_EXP1, std_EXP2, start_end, SUBBASIN = 'Pac', MASK_FILENAME = std_iodir_mask+'subbasins_orca21_nored.nc', _extra = ex
 

        ;;;;;;;;;; Zonal mean Salinity diff with Levitus
        ;;;;;;;;;; ========================================
  PRINT, ' '
  PRINT, '                               Zonal mean Salinity diff with Levitus'
  PRINT, '17 '
  xxx_zonal_S, S1, S2, SLev, std_EXP1, std_EXP2, start_end, MASK_FILENAME = std_iodir_mask+'subbasins_orca21_nored.nc', _extra = ex
                               
  PRINT, '17 '
  xxx_zonal_S, S1, S2, SLev, std_EXP1, std_EXP2, start_end, SUBBASIN = 'Atl', MASK_FILENAME = std_iodir_mask+'subbasins_orca21_nored.nc', _extra = ex
                              
  PRINT, '17 '
  xxx_zonal_S, S1, S2, SLev, std_EXP1, std_EXP2, start_end, SUBBASIN = 'Ind', MASK_FILENAME = std_iodir_mask+'subbasins_orca21_nored.nc', _extra = ex 

  PRINT, '17 '
  xxx_zonal_S, S1, S2, SLev, std_EXP1, std_EXP2, start_end, SUBBASIN = 'Pac', MASK_FILENAME = std_iodir_mask+'subbasins_orca21_nored.nc', _extra = ex 
  
        ;;;;;;; Seasonal & Annual mean Arctic/Antarctic Ice Thickness & Ice Concentration
        ;;;;;; ==========================================================================
  PRINT, ' '
  PRINT, '                               Arctic/Antarctic Ice Thickness'
  PRINT, ' 18'
  xxx_ARC_Icethick_MARS, Ithi_mars_1, Ithi_mars_2, std_EXP1, std_EXP2, start_end, _extra = ex
 
  PRINT, ' 19'
  xxx_ARC_Icethick_SEPT, Ithi_sept_1, Ithi_sept_2, std_EXP1, std_EXP2, start_end, _extra = ex
 
  PRINT, ' 20'
  xxx_ANT_Icethick_MARS, Ithi_mars_1, Ithi_mars_2, std_EXP1, std_EXP2, start_end, _extra = ex

  PRINT, ' 21'
  xxx_ANT_Icethick_SEPT, Ithi_sept_1, Ithi_sept_2, std_EXP1, std_EXP2, start_end, _extra = ex

  PRINT, '                               Arctic/Antarctic Ice Concentration (area of sea-ice per grid cell area'
  PRINT, ' 22'
  xxx_Iceleadfrac_MARS, Ifra_mars_1, Ifra_mars_2, std_EXP1, std_EXP2, start_end, _extra = ex

  PRINT, ' 23'
  xxx_Iceleadfrac_SEPT, Ifra_sept_1, Ifra_sept_2, std_EXP1, std_EXP2, start_end, _extra = ex

                       
        ;;; Global and Atlantic Meridional Function
        ;;; =======================================
  PRINT, ' '
  PRINT, ' 				Global and Atlantic Meridional Function'
  PRINT, ' 24 '
  plt_msf, V1, V2, std_EXP1, std_EXP2, start_end, SUBBASIN = 'Glo', MASK_FILENAME = std_iodir_mask+'subbasins_orca21_nored.nc', _extra = ex
  PRINT, ' 24 '
  plt_msf, V1, V2, std_EXP1, std_EXP2, start_end, SUBBASIN = 'Atl', MASK_FILENAME = std_iodir_mask+'subbasins_orca21_nored.nc', _extra = ex
  PRINT, ' 24 '
  plt_msf, V1, V2, std_EXP1, std_EXP2, start_end, SUBBASIN = 'Ind', MASK_FILENAME = std_iodir_mask+'subbasins_orca21_nored.nc', _extra = ex
  PRINT, ' 24 '
  plt_msf, V1, V2, std_EXP1, std_EXP2, start_end, SUBBASIN = 'Pac', MASK_FILENAME = std_iodir_mask+'subbasins_orca21_nored.nc', _extra = ex


  PRINT, '============================================= '
  PRINT, ' Local PLOTS '
  PRINT, '============================================= '

        ;;;;;;;;; Equatorial Temperature 
        ;;;;;;;;; ======================
  PRINT, ' '
  PRINT, '                               Equatorial Temperature'
  PRINT, ' '
  PRINT, ' 25 '
  xxx_Eq_T, T1, T2, Tlev, std_EXP1, std_EXP2, start_end, _extra = ex

        ;;;;;;;; Equatorial Salinity 
        ;;;;;;;; ===================
  PRINT, ' '
  PRINT, '                               Equatorial Salinity'
  PRINT, ' 26 '
  xxx_Eq_S, S1, S2, SLev, std_EXP1, std_EXP2, start_end, _extra = ex

        ;;; Equatorial zonal velocity
        ;;; =========================
  PRINT, ' '
  PRINT, '                               Equatorial zonal velocity'
  PRINT, ' 27 '
  xxx_Eq_U, U1, U2, std_EXP1, std_EXP2, start_end, _extra = ex
  
       ;;;;;;;;;;;; Mediterranean salt tongue 
       ;;;;;;;;;;;; =============================================== 
  
  PRINT, ' '
  PRINT, '                               Mediterranean salt tongue at depth=700 and lat 40°N'
  PRINT, ' 28 '
  xxx_Med_S_tongue, S1, S2, SLev, std_EXP1, std_EXP2, start_end, DEPTH = 700,  LAT = 40, _extra = ex
  ;
  PRINT, '                               Mediterranean salt tongue at depth=1000 and lat 38°N'
  PRINT, ' 28 '
  xxx_Med_S_tongue, S1, S2, SLev, std_EXP1, std_EXP2, start_end, DEPTH = 1000, LAT = 38, _extra = ex

        ;;;;;;;;; Mediterranean 
        ;;;;;;;;; ===============================================
  PRINT, ' '
  PRINT, '                               Mediterranean water at 38°N'
  PRINT, ' 29 '
  xxx_Med_S_depth, S1, S2, SLev, std_EXP1, std_EXP2, start_end, DEPTH = 1000, LAT = 38, _extra = ex
  PRINT, '                               Mediterranean water at 38°N'
  PRINT, ' 29 '
  xxx_Med_S_depth, S1, S2, SLev, std_EXP1, std_EXP2, start_end, DEPTH = 700,  LAT = 40, _extra = ex

;
     
        ;;;;;;;;;;; Vertical Global mean T & S
        ;;;;;;;;;;; ===========================
  PRINT, ' '
  PRINT, '                                Vertical Global mean T & S'
  PRINT, ' 30 '
  xxx_TS_ZGlobm, T1, T2, TLev, S1, S2, SLev, std_EXP1, std_EXP2, start_end, _extra = ex

 
  return
END