begin

 err = NhlGetErrorObjectId()
 setvalues err
    "errLevel" : "Fatal"          ; only report Fatal errors
 end setvalues

 type         = "png"
 iLead        = ILEAD
 Szn4Month    = SZNLENLGIC    ; False ==  3 month season forecast data ;
 STATE        = LOGICSTATE 
 iMonth       = OMONS
 currentYr    = CURRENTYR  
 SEASON       = SEASONNAME
 latS         = LATS
 latN         = LATN
 lonW         = LONW
 lonE         = LONE
 PRDCTR       = SLCTDPRD ;"Nino34" ;"WVG"  
 SigPlev      = SIGLOGIC 
 COUNTRY      = MYCOUNTRY ;one of the countries below as displayed.
 resl         = RESL
 ZoomIn       = ZOOMIN
                        ; The order is the same as maskNameNCL
 StateNames =     (/ "FJI",    "FSM",    "PLW",   "PNG", \
                     "SLB",    "VUT",   "WSM", "TUV", \
                     "COK",    "KIR" /)
 maskNameNCL := (/"Fiji",   "Federated States of Micronesia",  "Palau", "Papua New Guinea",\
                  "Solomon Islands" , "Vanuatu", "American Samoa","Tuvalu", \
                  "Cook Islands","Kiribati" /)


 dircolor   =  "../colors/"
 dirSHP     =  "../shpfiles/"
 regInd    := (/0,1,2,3,4,5,6,7,8,9 /)
  maskcPa    = maskNameNCL(regInd)
           ;  FJI   FSM   PLW   PNG    SLB   VUT   WSM   TUV   COK  KIR/-17.8
 sLatMin = (/-19.5, 1.0,  6.5,  -11., -11.0, -20., -14.5,  -10., -22., -3.5 /)
 sLatMax = (/-16.0, 9.0,  8.0,  -1.,  -5.5,  -13., -13.0,  -5.,  -8.,  2.5/)
 sLonMin = (/177.0, 137., 134., 141., 156.0, 166., -173., 176., -165., -174.5/)
 sLonMax = (/180.5, 163., 135., 157., 163.0, 170., -171., 179., -157., 177.0 /)
 allMask = (/True, False, True, True, True, True, False, False, False, False /)

;;;;;;;;;;;;;
 dirModelRoot     = "./"
 gcmnames         = (/"CanSIP","NASA","GFDLS","CCSM4", "CFSv2"/)
 smodStart        = 0
 smodLast         = 4

 RASTER           = True 
 plotcombined     = True
 plotProbSeparate = False 

; maskPlot         = False
; maskLev          = 5
 shaderaindry     = False
;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;
;;; Set lat-lon of domain

  ;;; directory where colors and shpfiles dirs  are located

 useNCLmask       = True
 if (STATE) then
     StateInd = ind(StateNames .eq. COUNTRY)
     StName = StateNames(StateInd)
     print(" ")
     print("... working for "+ str_capital(maskNameNCL(StateInd)) )  
      f =addfile(dirSHP+"gadm41_"+StName+"_shp/"+"gadm41_"+StName+"_0.shp","r")
      state_lon             = f->x
      state_lat             = f->y
      delete(f)
      state_lon         = where(state_lon .lt. 0., state_lon+360.0,state_lon)
      useNCLmask = allMask(StateInd)

 else
     print(" ")
     print("... working for "+ str_capital(maskcPa) )  
 end if

 landocean = (/"ocean", "land"/)
 if (lonW .lt. 0.) then
    lonW = lonW + 360.
 end if

 if (lonE .lt. 0.) then
    lonE = lonE + 360.
 end if

;;;; Need to have the lower r/wesolution at 0.5 for real-time verification 
     mlatLres            = toint((latN - latS)/resl + 1)
     mlonLres            = toint((lonE - lonW)/resl + 1)
     LresLat        = fspan(latS, latN, mlatLres)
     LresLon        = fspan(lonW, lonE, mlonLres)
     LresLat@units  = "degrees_north"
     LresLon@units  = "degrees_east"

   month      = (/"Jan","Feb","Mar","Apr","May","Jun", \
               "Jul","Aug","Sep","Oct","Nov","Dec" /)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

   FTYPES   =  "_PRCP_Logit"+PRDCTR+"_3monthSeasonal.nc"

  if ( Szn4Month) then
     FTYPES   = "_PRCP_Logit"+PRDCTR+"_4monthSeasonal.nc"
  end if
 FileSeasonal = FTYPES

 figmethods = "_Prob_Logit"+PRDCTR

;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;
; useNCLmask       = True    
 pthr  = 3.  ; probability difference (abs) b/n above and below 
              ; categories below which both categories are set to 
              ; missing.  
 fthr  = 3.  ; same as pthr

 VARN   = "Prec"
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;;;;; Obsfiles for plotting mean, making dry areas

 minProbLev   = 33   ; cut off probability

 countmincat  = 10    ; GHA must have at least 30 point to contour probability


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;; Create directories for figures
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   if (STATE) then
      dir0 = "./epsfiles/" + COUNTRY +"/"
   else
      dir0 = "./epsfiles/"
   end if
   system("if ! test -d " + dir0 +" ; then mkdir -p " + dir0 + " ; fi")
   plotdir = dir0
   delete(dir0)

;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;
; Set some basic resources for all

  res                                 = True
  res@gsnFrame                        = False
  res@gsnDraw                         = False
  res@gsnLeftString                   = ""
  res@gsnRightString                  = ""
  res@tfDoNDCOverlay                  = False
  res@gsnMaximize                     = False
  res@gsnAddCyclic                    = False
  res@mpDataBaseVersion               = "MediumRes"
  res@mpDataSetName                   = "Earth..4"
  res@mpGeophysicalLineThicknessF     = 4
  res@mpNationalLineThicknessF        = 4
  res@mpOutlineBoundarySets           = "Allboundaries"
  res@pmTickMarkDisplayMode           = "Always"
  res@mpLimitMode                     = "LatLon"
  res@mpFillOn                        = False
  res@mpOutlineOn                     = False
  res@cnInfoLabelOn                   = False
  res@cnConstFLabelOn                 = False
  res@cnFillOn                        = True
  res@cnLinesOn                       = False
  res@cnConstFEnableFill              = False
  res@cnNoDataLabelOn                 = False
  res@mpCenterLonF                    = 180

  if (STATE) then
     if (ZoomIn) then

        min_lat               = sLatMin(StateInd)
        max_lat               = sLatMax(StateInd)
        min_lon               = sLonMin(StateInd)
        max_lon               = sLonMax(StateInd)

        if (min_lon .lt. 0.) then
            min_lon = min_lon + 360.
        end if

        if (max_lon .lt. 0.) then
            max_lon = max_lon + 360.
        end if

        res@mpMinLatF            = min_lat     ; range to zoom in on
        res@mpMaxLatF            = max_lat
        res@mpMinLonF            = min_lon
        res@mpMaxLonF            = max_lon

     else

        min_lat               = min(state_lat)
        max_lat               = max(state_lat)
        min_lon               = min(state_lon)
        max_lon               = max(state_lon)
        res@mpMinLatF            = min_lat     ; range to zoom in on
        res@mpMaxLatF            = max_lat
        res@mpMinLonF            = min_lon
        res@mpMaxLonF            = max_lon
     end if

     f1shp = dirSHP + "gadm41_"+StName+"_shp/" + "gadm41_"+StName+"_0.shp"

    if (useNCLmask) then
         res@mpFillOn                = True
         res@mpOutlineOn             = False
         res@mpFillBoundarySets      = "NoBoundaries"
         res@mpOutlineBoundarySets   = "National"
         res@mpFillAreaSpecifiers    = landocean
         res@mpSpecifiedFillColors   = (/"white","white" /)
         res@mpAreaMaskingOn         = 1
         res@mpMaskAreaSpecifiers    = (/maskNameNCL(StateInd)/)
         res@cnFillDrawOrder         = "Predraw"
         res@cnLineDrawOrder         = "Predraw"
    end if

  else
     res@mpMinLatF            = latS      ; range to zoom in on
     res@mpMaxLatF            = latN
     res@mpMinLonF            = lonW
     res@mpMaxLonF            = lonE

     f1shp                  = dirSHP + "ppacific.shp"
     f1shp                  = dirSHP + "allPacific.shp"
     useNCLmask  = False
     if (useNCLmask) then
         res@mpFillOn                = True
         res@mpOutlineOn             = False
         res@mpFillBoundarySets      = "NoBoundaries"
         res@mpOutlineBoundarySets   = "National"
         res@mpFillAreaSpecifiers    = landocean
         res@mpSpecifiedFillColors   = (/"white","white" /)
         res@mpAreaMaskingOn         = 1
         res@mpMaskAreaSpecifiers    = maskcPa 
         res@cnFillDrawOrder         = "Predraw"
         res@cnLineDrawOrder         = "Predraw"
     end if
  end if

  resproball             = res
  resprob                = res

  res@cnLevelSelectionMode = "ExplicitLevels"

  res@vpWidthF           = 0.65
  res@vpHeightF          = 0.8
  res@vpXF               = 0.13
  res@vpYF               = 0.925

  resprob               = res
  resproball            = res

  enspall                        = resproball
  enspall@cnLevelSelectionMode     = "ExplicitLevels"
  enspall@cnLevels                 = (/20,30,40,50,60,70,80,90,100/)
  enspall@cnLineLabelsOn           = False
  enspall@lbLabelBarOn             = True
  enspall@lbOrientation            = "Vertical"

;;;;;;;;;;;;;;;;;;;;;;
;;; Read color files 
;;;;;;;;;;;;;;;;;;;;;;

  colorabove   = asciiread(dircolor+"colAbovePrec.txt",(/11,3/),"float")
  colorbelow   = asciiread(dircolor+"colBelowPrec.txt",(/11,3/),"float")
  colornormal  = asciiread(dircolor+"colNormalPrec.txt",(/11,3/),"float")
  coloranomaly = asciiread(dircolor+"colAnomalyPrec.txt",(/23,3/),"float")

  problevels   = (/40,45, 50,60,70,80,90/)
  problabels   = problevels + ""
  pcolind      = (/0,2,4,5,6,8,9,10/)

  ;colornormal  := (/"grey80","grey81","grey82","grey83","grey84","grey85", \
  ;                 "grey86","grey87","grey88","grey89","grey90" /)

  lnres                  = True
  lnres@gsLineColor      = "black"
  lnres@gsLineThicknessF = 2.0

  ensp  = resprob
  ensp@cnLevelSelectionMode     = "ExplicitLevels"
  ensp@cnLevels                 = problevels
  ensp@cnLineLabelsOn           = False
  ensp@lbLabelBarOn             = False
  ensp@cnFillMode                = "RasterFill"
;  ensp@lbOrientation            = "Vertical"

   bensp                           = True
   bensp@gsnDraw                   = False
   bensp@gsnFrame                  = False
   bensp@cnInfoLabelOn             = False
   bensp@cnFillOn                  = True
   bensp@cnFillMode                = "RasterFill"
   bensp@cnLinesOn                 = False
   bensp@cnLineLabelsOn            = False
   bensp@cnNoDataLabelOn           = False
   bensp@cnConstFEnableFill        = True
   bensp@cnConstFLabelOn           = False
   bensp@lbLabelBarOn              = False
   bensp@mpCentrLonF               = 180.

   bensp@cnLevelSelectionMode      = "ExplicitLevels"
   bensp@cnLevels                  = problevels

   if (useNCLmask) then
       bensp@cnFillDrawOrder         = "Predraw"
       bensp@cnLineDrawOrder         = "Predraw"
   end if

   if (RASTER) then
      ensp@cnFillMode                = "RasterFill"
      bensp@cnFillMode               = "RasterFill"
   end if

;;;;;;;
;;; for labelbar
  plbres                      = True

; Set initial width and height.
  plbres@vpWidthF             = 0.03

; Allow more control over labelbars.
  plbres@lbAutoManage         = False

; No margins around labelbar.
  plbres@lbBottomMarginF      = 0.0
  plbres@lbLeftMarginF        = 0.0
  plbres@lbRightMarginF       = 0.0
  plbres@lbTopMarginF         = 0.0

; Turn various features on and off.
  plbres@lbLabelsOn           = True
  plbres@lbPerimOn            = False
  plbres@lbTitleOn            = True
  plbres@lbMonoFillPattern    = True
  plbres@lbOrientation        = "Vertical"
  plbres@lbLabelFontHeightF   = 0.013

  plbres@lbJustification      = "TopLeft"

  plbres@lbLabelAlignment     = "InteriorEdges"

  plbres@lbLabelJust          = "CenterLeft"
  plbres@lbLabelOffsetF       = 0.5


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;; Read Data and plot      
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  modDatDir  = dirModelRoot + "analdat/"
  filesReg = modDatDir
  allfiles = filesReg + gcmnames + FTYPES
  fdata = addfiles(allfiles,"r")
  ListSetType (fdata, "join")
  ;print(allfiles)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;ENSMEAN PROBABILITY
;;;;;
     if (plotcombined) then
           tbelow  := fdata[:]->below
           tabove  := fdata[:]->above
           tnormal := fdata[:]->normal
           bsig    := fdata[:]->blwsig
           asig    := fdata[:]->abvsig
           nsig    := fdata[:]->nrmsig
           tbelow  := dim_avg_n_Wrap(tbelow,0)
           tabove  := dim_avg_n_Wrap(tabove,0)
           tnormal := dim_avg_n_Wrap(tnormal,0)
           tnormal  = (/ 100. - (tbelow + tabove) /)
           bsig    := dim_avg_n_Wrap(bsig,0) 
           asig    := dim_avg_n_Wrap(asig,0) 
           nsig    := dim_avg_n_Wrap(nsig,0)
 
           if (SigPlev) then
             tbelow = (/ tbelow * bsig/)
             tabove = (/ tabove * asig/)
             tnormal= (/ tnormal* nsig/)
           end if

        ;printVarSummary(tbelow) 
        tbelow@units = ""
        tabove@units = ""
        tnormal@units = ""

;;;;;;;;
        t0 := tbelow
        t0 = t0@_FillValue
        t1 := t0
        t2 := t0

        t0Raw := t0
        t1Raw := t1
        t2Raw := t2
;;;;;
;;;;;ENSEMBLE PROCESSED PROBABILITY DATA FOR SEPARATE PLOTING
;;;;;

        dim := dimsizes(t0)
        ydim = dim(0) - 1
        xdim = dim(1) - 1
        do j = 0,ydim
           do i = 0, xdim
              vv   = new(3, float , -9999 )

              vv(0) = tbelow(j,i)
              vv(1) = tnormal(j,i)
              vv(2) = tabove(j,i)

              t0Raw(j,i) = (/ vv(0) /)
              t1Raw(j,i) = (/ vv(1) /)
              t2Raw(j,i) = (/ vv(2) /)

              ix = maxind(vv)
              if (ismissing(ix)) then
                 t0(j,i) = (/t0@_FillValue/)
                 t1(j,i) = (/t1@_FillValue/)
                 t2(j,i) = (/t2@_FillValue/)
              else
                 if ( ix .eq. 0 .and. vv(0) .ge. minProbLev) then
                    t0(j,i) = (/vv(0)/)
                 end if
                 if (ix .eq. 1 .and. vv(1) .ge. minProbLev) then
                    t1(j,i) = (/vv(1)/)
                 end if
                 if (ix .eq. 2 .and. vv(2) .ge. minProbLev)  then
                    t2(j,i) = (/vv(2)/)
                 end if
              end if

              v02 := abs(vv(0)   - vv(2))
              v01 := abs(vv(0)   - vv(1))
              v12 := abs(vv(1)   - vv(2))

              if (.not. ismissing(ix)) then
                 if (ix .eq. 0 .or. ix .eq. 2) then
                    if ( .not. ismissing(v02) .and. v02 .le. pthr) then
                       t0(j,i) = (/t0@_FillValue/) ;white if no signf. diff
                       t2(j,i) = (/t2@_FillValue/) ;white if no signf. diff
                       if (.not. ismissing(vv(1)) .and. vv(1) .gt. minProbLev) then
                          t1(j,i) = (/ vv(1) /)
                       else
                          t1(j,i) = (/t1@_FillValue/)
                       end if
                       ix  = 102  ; Normal favored for A or B when A-B < pthr
                    end if
                 end if

                 if (ix .eq. 0 .or. ix .eq. 1) then   ; favor normal
                    if ( .not. ismissing(v01) .and. v01 .le. pthr) then
                       t0(j,i) = (/t0@_FillValue/)  ;white if no signf. diff
                       t2(j,i) = (/t2@_FillValue /) ;white if no signf. diff
                       if (.not. ismissing(vv(1)) .and. vv(1) .gt. minProbLev) then
                          t1(j,i) = (/ vv(1) /)
                       else
                          t1(j,i) = (/t1@_FillValue/)
                       end if
                       ix  = 101 ; Normal favored for B or N when B-N < pthr
                    end if
                 end if

                 if (ix .eq. 1 .or. ix .eq. 2) then   ; favor normal
                    if ( .not. ismissing(v12) .and. v12 .le. pthr) then
                       t0(j,i) = (/t0@_FillValue/)  ;white if no signf. diff
                       t2(j,i) = (/t2@_FillValue /) ;white if no signf. diff
                       if (.not. ismissing(vv(1)) .and. vv(1) .gt. minProbLev) then
                          t1(j,i) = (/ vv(1) /)
                       else
                          t1(j,i) = (/t1@_FillValue/)
                       end if
                       ix  = 112  ;Normal favored for A or N when A-N < pthr
                    end if
                 end if
              else
                 t0(j,i) = (/t0@_FillValue/)  ; white  no sig.
                 t2(j,i) = (/t2@_FillValue /) ; white if no signf.
                 t1(j,i) = (/ vv(1) /)
                 ix  = 199  ;Normal favored for indeterminate
              end if

              delete([/ix,vv/])
;              vv = maskPercent(j,i)
;              if (.not. ismissing(vv) .and. vv .lt. maskLev .and. maskPlot) then
;                 t0(j,i) = t0@_FillValue
;                 t1(j,i) = t0@_FillValue
;                 t2(j,i) = t0@_FillValue
;              end if
;              delete(vv)
           end do ; i xdiim
        end do ; j ydim
;           t0 = (/t0 * mskrgn   /)
;           t1 = (/t1 * mskrgn   /)
;           t2 = (/t2 * mskrgn   /)

;        if (STATE) then
;           fmask = addfile(dirSHP+COUNTRY+"MaskP1deg.nc","r")
;          xmask    := fmask->mask
;          t0 = (/t0 * xmask   /)
;          t1 = (/t1 * xmask   /)
;          t2 = (/t2 * xmask   /)
;        end if


        t0file := t0
        t1file := t1
        t2file := t2

;;;;;
;;;;;PLOT ENSEMBBLE BELOW PROBABILITY
;;;;;
        if (SigPlev) then
            xname  = plotdir +  "Ens_"+VARN + "_Signf"  + \
                     "_"+iLead+"monLead" + figmethods
        else
            xname  = plotdir +  "Ens_"+VARN + "_NoSignf" +  \
                     "_"+iLead+"monLead" +  figmethods
        end if

        wks = gsn_open_wks(type,xname)

        tmp = t0
        ensp@cnFillPalette             = colorbelow(pcolind,:)
        cntrb = gsn_csm_contour_map(wks,tmp,ensp)
        if (type .eq. "png" .or. type .eq. "eps" .or. STATE) then
          poly1 = gsn_add_shapefile_polylines(wks,cntrb,f1shp,lnres)
        end if
        delete([/tmp,ensp@cnFillPalette/])

;;;;;
;;;;;PLOT ENSEMBLE NEAR NORMAL PROBABILITY
;;;;;
        tmp = t1
        bensp@cnFillPalette            = colornormal(pcolind,:)
        cntrn = gsn_csm_contour(wks,tmp,bensp)
        delete([/tmp,bensp@cnFillPalette/])

;;;;;
;;;;;PLOT ENSEMBLE ABOVE PROBABILITY
;;;;;
        tmp = t2
        bensp@cnFillPalette             = colorabove(pcolind,:)
        cntra = gsn_csm_contour(wks,tmp,bensp)
        delete([/tmp,bensp@cnFillPalette/])

;;;;;
;;;;;ENSEMBLE SHADED AREA BY PERCENT
;;;;;
;        if (shaderaindry) then
;            x := maskPercent
;            tmp := x
;            tmp = (/ where(x .ge. maskLev, x@_FillValue, 0.) /)
;            cntrx = gsn_csm_contour(wks,tmp,xensp)
;           ; delete([/tmp,x/])
;        end if

;;;;;
;;;;;ENSEMBLE LABELBAR
;;;;;
        txres               = True
        txres@txFontHeightF = 0.013
        txres@txAngleF      = 90.

        getvalues cntrb                      ; Get plot size for use in
           "vpHeightF" : vph                 ; creating labelbar.
           "vpWidthF"  : vpw
           "vpXF"      : vpx
           "vpYF"      : vpy
        end getvalues
        plbres@vpHeightF            = vph/3. +0.04

        xpos = vpw+vpx+0.01
        ypos = (1.-vpy)+ vph/3 + 0.100 -0.02
        ytxt = ypos - 0.08
        plbres@lbBoxMinorExtentF    = 1.0
        plbres@lbFillColors         = colorbelow(pcolind,:) ;plbresB
        lbid1 = gsn_create_labelbar_ndc(wks,dimsizes(problevels)+1,problabels, \
                                   xpos,ypos,plbres)
        gsn_text_ndc(wks,"     Below (%) " ,xpos+0.11 ,ytxt,txres)

        ypos = ypos + vph/3 + 0.015
        ytxt = ypos - 0.07
        delete(plbres@lbFillColors)
        plbres@lbFillColors         = colornormal(pcolind,:)
        lbid2 = gsn_create_labelbar_ndc(wks,dimsizes(problevels)+1,problabels, \
                                   xpos,ypos,plbres)
        gsn_text_ndc(wks,"    Normal (%) " ,xpos+0.11,ytxt,txres)
        delete(plbres@lbFillColors)


        ypos = ypos + vph/3 + 0.015
        ytxt = ypos - 0.07
        plbres@lbFillColors         = colorabove(pcolind,:)
        lbid3 = gsn_create_labelbar_ndc(wks,dimsizes(problevels)+1,problabels, \
                                   xpos,ypos,plbres)
        gsn_text_ndc(wks,"   Above (%) " ,xpos+0.11,ytxt,txres)

;;;;;
;;;;;ENSEMBLE OVERLAY
;;;;;
        overlay(cntrb,cntrn)
        overlay(cntrb,cntra)
        if (shaderaindry) then
           overlay(cntrb,cntrx)
           draw(cntrb)
        else
           draw(cntrb)
        end if

        draw(lbid1)
        draw(lbid2)
        draw(lbid3)

        frame(wks)
        delete(wks)
        delete([/txres,lbid1,lbid2,lbid3/])

;;;;;
;;;;;ENSEMBLE PLOT SEPARATE PROBABILITY
;;;;;
        if (plotProbSeparate) then ; Ensemble  separate

           if (SigPlev) then
               figName  = "Ens_"+VARN + "_Signf"  + \
                     "_"+iLead+"monLead" + figmethods
           else
               figName  = "Ens_"+VARN + "_NoSignf" +  \
                     "_"+iLead+"monLead" + figmethods
           end if

           wks = gsn_open_wks(type,plotdir+figName+"Above")

           tmp := tabove 
           ;tmp = (/ tmp * mskrgn /)
           enssep = enspall
           enssep@cnLevels    := (/20, 25,30,35,40, 45, 50,55,60/)
           enssep@cnLevels    := (/40,41,42,43,44, 45, 50,55,60,70,80/)
           enssep@cnFillPalette  := colorabove

           cntra = gsn_csm_contour_map(wks,tmp,enssep)
           poly1 = gsn_add_shapefile_polylines(wks,cntra,f1shp,lnres)
           draw(cntra)
           frame(wks)
           delete(wks)

           wks = gsn_open_wks(type,plotdir+figName+"Normal")
           tmp := tnormal
           ;tmp = (/ tmp * mskrgn /)
           enssep = enspall
           enssep@cnLevels    := (/20, 25,30,35,40, 45, 50,55,60/)
           enssep@cnLevels    := (/40,41,42,43,44, 45, 50,55,60,70,80/)
           enssep@cnFillPalette  = colornormal
           cntra = gsn_csm_contour_map(wks,tmp,enssep)
           poly1 = gsn_add_shapefile_polylines(wks,cntra,f1shp,lnres)
           draw(cntra)
           frame(wks)
           delete(wks)

           wks = gsn_open_wks(type,plotdir+ figName+"Below")
           tmp := tbelow 
           ;tmp = (/ tmp * mskrgn /)
           enssep = enspall
           enssep@cnLevels    := (/40,41,42,43,44, 45, 50,55,60,70,80/)
           enssep@cnFillPalette  := colorbelow

           cntra = gsn_csm_contour_map(wks,tmp,enssep)
           poly1 = gsn_add_shapefile_polylines(wks,cntra,f1shp,lnres)
           draw(cntra)
           frame(wks)
           delete(wks)

           delete(tmp)
        end if ; ensemble separate plots

        delete([/ tbelow, tabove,tnormal/])

     end if  ; plotcombined

end