PRO plot_strmfz_fesom_jrpaper

   close, /all
   ;; Domain
   jpj=180L
   nk=33L
   PS=1
   ps_dir='/Net/gleam/abozec/HYCOM/JRA55-res_PAPER/PS_PAPER/FESOM/'

   ;; PATH
   io='/Net/gleam/abozec/HYCOM/JRA55-res_PAPER/FESOM/NEW/OMIP3_LR/loop1/'

   ;; read amoc files
   tmp=ncdf_lec(io+'msftyz_variant2_Omon_FESOM1.4_historical_loop1_gr_198001-201812.nc',var='msftyz_variant2')
   mean_amoc=transpose(reform(tmp[1,*,*]))
   tmp=ncdf_lec(io+'msftyz_var_variant2_Omon_FESOM1.4_historical_loop1_gr_198001-201812.nc',var='msftyz_var_variant2')
   vari_amoc=transpose(reform(tmp[1,*,*]))
   stdv_amoc=sqrt(vari_amoc)

   ;; plot
   file='msftyz_variant2_Omon_FESOM1.4_historical_loop1_gr_198001-201812.nc'
   lati = ncdf_lec(io+file, var = 'lat' )
   depth = ncdf_lec(io+file, var = 'depth_coord' )

   ;; let's plot the Atlantic streamfunction
   levels = 20
   Maxss = 22 & Minss = -8.
   step = (Maxss - Minss) / levels
   userLevels = IndGen(levels) * step + Minss
   LCT, 64, NColors=levels, Bottom=3

   ind = where(mean_amoc LE Minss )
   IF (ind NE [-1]) THEN mean_amoc(ind) = Minss


   if (PS EQ 1) then openps,ps_dir+'meanvar_amoc_yz_FESOMLR_1980-2018.ps'
   scontour, mean_amoc,lati, depth,  xstyle = 1, ystyle = 1, cell_fill=2, petit = [1, 2, 1],$
    Levels=userLevels, Color=!P.Color, Background=!P.Background, charsize = 1, /rempli, $
    C_Colors=IndGen(levels) + 3, yrange = [6500, 0], xrange = [-30., 80.], xtickinterval = 20, $
    ytitle = 'Depth', title = ' FESOM low res. Mean Atlantic Streamfunction JRA55 1980-2018', /portrait, window = 1
   scontour, mean_amoc,lati, depth, /overplot, Levels=userLevels, color = 0.

   ColorBar2, Divisions=levels, Range=[Minss, Maxss], Format='(f4.1)', $
    Position = [0.1365, 0.54, 0.92, 0.56], Color=!P.Color, $
    NColors=levels, Bottom=3   ;; let's plot the global streamfunction

   ;; let's plot the Atlantic streamfunction
   levels = 20
   Maxss = 6 & Minss = 0.
   step = (Maxss - Minss) / levels
   userLevels = IndGen(levels) * step + Minss
   LCT, 64, NColors=levels, Bottom=3

   ;ind = where(vari_amoc GE  Maxss )
   ;IF (ind NE [-1]) THEN vari_amoc(ind) = Maxss

   scontour, stdv_amoc,lati, depth,  xstyle = 1, ystyle = 1, cell_fill=2, petit = [1, 2, 2],$
    Levels=userLevels, Color=!P.Color, Background=!P.Background, charsize = 1, /rempli, $
    C_Colors=IndGen(levels) + 3, yrange = [6500, 0], xrange = [-30., 80.], xtickinterval = 20, $
    ytitle = 'Depth', title = 'FESOM low res. Standard Dev. Atlantic Streamfunction JRA55 1980-2018', /portrait, window = 1,/noerase
   scontour, stdv_amoc,lati, depth, /overplot, Levels=userLevels, color = 0.

   ColorBar2, Divisions=levels, Range=[Minss, Maxss], Format='(f3.1)', $
    Position = [0.1365, 0.08, 0.92, 0.10], Color=!P.Color, $
    NColors=levels,Bottom=3  ;; let's plot the global streamfunction
    if (PS EQ 1) then  closeps

stop

end