;  Script to plot and overlay contours on a map plot (using gsn_csm_map) 


begin

;  - - - - -  Reading the file  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

dir = "/users/home/models/csm/inputdata/ocn/nemo/tn0.25v1/forcing/"
key = "chlorophyll.nc"

f = addfile(dir+key,"r")

lat = f->nav_lat                 ; [nj | 1020] x [ni | 1440]
lon = f->nav_lon
nav_lat1 = lat(:,1060)
nav_lon1 = lon(500,:)

nav_lona = lon(500,1150:1441)
nav_lonb = lon(500,:149)
dims = dimsizes(lon)
;print(max(lat(:,1060)))
;print(lat(:100,1060))

chla = f->chla(0,:,:)
chla!0 = "Lat"
chla!1 = "Long"
chla&Lat = nav_lat1
chla&Long = nav_lon1

chla = chla
buff = chla
buff!0 = "Lat"
buff!1 = "Long"
buff&Lat = nav_lat1
buff&Long = nav_lon1

; for later projection purposes (possibly)
lat0=50.
lat1=90.
lon0=40.
lon1=80.
lev0=10.
lev1=26.
deltalev=2.

printVarSummary(chla)
printVarSummary(buff)


; - - - - - - - - -  LOOP for masking a particular area - - - - - - - - - - - - - - - -- - - - - 

 
  do j=0, dimsizes(nav_lat1)-1
   do i=0, dimsizes(nav_lon1)-1

    chla(j,i)  = where ( lat(j,i).ge.66 .and. lat(j,i) .le. 80 .and. lon(j,i) .ge. 25 .and. lon(j,i) .le. 100, 1., 0.)
    buff(j,i) = where ( lat(j,i).ge.64 .and. lat(j,i) .le. 82 .and. lon(j,i) .ge. 24 .and. lon(j,i) .le. 101, 0.8, 0.)
  ; buff2(j,i) = where ( lat(j,i).ge.65 .and. lat(j,i) .le. 81 .and. lon(j,i) .ge. 24 .and. lon(j,i) .le. 101, 0.8, 0.)

    end do
  end do


 ;  PLOT ..... ------------------------------------------

 pltDir = "./"
 pltName = "overl_mask"
 pltType = "png"
 pltPath = pltDir+pltName

 wks = gsn_open_wks("png","map_overlay")

 gsn_define_colormap(wks, "BlueDarkRed18")

 mpres = True		; Map resources
 mpres@mpProjection        = "Satellite"  ; choose map projection
 mpres@mpCenterLonF        =  45.         ; choose center lon
 mpres@mpCenterLatF        =  70.         ; choose center lat
 mpres@mpLimitMode         = "LatLon"     ; required
 mpres@mpMinLatF           =  50.         ; min lat
 mpres@mpMaxLatF           =  120.        ; max lat
 mpres@mpMinLonF           =  10.         ; min lon
 mpres@mpMaxLonF           =  90.         ; max lon
 mpres@sfXArray            =  lon
 mpres@sfYArray            =  lat
 mpres@mpGridSpacingF        =  5.
 mpres@mpGridAndLimbOn       =  True         ; turn on lat/lon lines
 mpres@pmTickMarkDisplayMode = "Always"     ; turn on automatic tickmarks

 res                       =  True         ; plot mods desired
 res		           =  mpres
 res@gsnDraw               =  False
 res@gsnFrame              =  False
 res@cnFillOn              =  True         ; turn on color
 res@cnLinesOn             =  False        ; turn off contour lines
 res@cnFillDrawOrder       =  "PreDraw"    ; fill first
 res@cnLevelSelectionMode  =  "ExplicitLevels"   ; set explicit contour levels
 res@cnLevels    	   =  (/.1,.2,.3,.4,.5,.6,.7,.8,.9/)   ; set levels
 res@sfXArray              =  lon
 res@sfYArray              =  lat
 res@lbBoxLinesOn          =  False        ; turn off box between colors
 res@gsnSpreadColors       =  True         ; use full colormap
 res@gsnAddCyclic          = False 

 res2		=  True
 res2 		=  res

 map_plot  =  gsn_csm_map(wks, mpres)
 chla_plot =  gsn_csm_contour(wks, chla, res2)
 buff_plot = gsn_csm_contour(wks, buff, res)

 overlay(map_plot, buff_plot)
; overlay(map_plot, chla_plot) 

 draw(map_plot)
 frame(wks)

end