;;  Adapted from idl script for PVWave  12 Sept 2005
;;  It reads MFT model output files and plots contours of the
;;  tracer concentration in a ps file.
;;
;; First created: For 4 views     J G Li   26 Nov 2008
;; Last modified: Color grids     J G Li   26 Aug 2009
;; Adapted for 25km SMC grids     J G Li    3 Feb 2010
;; Adapted for 6-25km SMC grids   J G Li   18 Feb 2010
;; Add Arctic boundary rings.     J G Li   10 Jun 2011
;; S4R sterographic projection.   J G Li   30 Jun 2011
;; Modified for SMC625 grid.      J G Li   15 Dec 2011
;; Add 33 buoy locations.         J G Li   11 Apr 2012
;

 nc=0L
 ng=0L
 n1=0L
 n2=0L
 n4=0L
 n8=0L
 n9=0L

 openr, 9, 'DatGMC/SMC625Cels.dat'
 readf, 9, ng, n1, n2, n4 
 cel=intarr(5,ng)
 readf, 9, cel
 close, 9

 Print, ' GloMCels read done'
 Print, ng, n1, n2, n4 
 Print, cel(*,0)
 Print, cel(*,ng-1)

;; Maximum j row number in Global part
 jmxglb = max( cel(1,*) )
 Print, ' Maximum j row =', jmxglb

;;  Arctic part from SMC625Arcl.dat
 na=0L

 openr, 8, 'DatGMC/SMC625BArc.dat'
 readf, 8, na, n8, n9
 ael=intarr(5,na)
 readf, 8, ael
 close, 8

 Print, ' SMC625BArc read done'
 Print, na, n8, n9
 Print, ael(*,0)
 Print, ael(*,na-1)

;;  Merge Arctic part but excluding boundary cells together
 nb=n8+n9  
 nc=ng + na - nb 
 arl=ael(*,nb:na-1)
 cel=transpose([transpose(cel), transpose(arl)])
 Print, nc, n1, n2, n4, n8, n9 


 d2rad=!Pi/180.0
 radius=10.0

;; Number of radius of projection distance
 np=4.0

;; Projected Equator radius * np to be (np-1)*radius
 PRadus=(np-1.0)*radius

 dxlon=0.087890635D
 dylat=0.05859375D
;dxlon=0.3515625D
;dylat=0.2343750D

;; Outline circle at equator
 ciran=findgen(1081)*d2rad/3.0
 xcirc=radius*cos(ciran)
 ycirc=radius*sin(ciran)

 colorlevls=[(indgen(11)+1)*23]
 heights=[-0.1, 0.1, 0.5, 0.9, 1.1, 2.0, 3.0, 4.0, 4.9, 5.1, 6.0]

 nx=4097
 ny=3073
;;  Note full grid from south to north poles requires 3072*dlat

 xlon=(Dindgen(nx))*dxlon
 ylat=(Dindgen(ny)-1536.0D)*dylat
 jeqt=1536
;;  Equator is now at ylat(1536)=0.0 

;; Color levels to be used
  mclr=128
  colrs=[indgen(mclr)]
  depth=[10000,1000,100,10,1]
  cdeps=TRANSPOSE( cel(4,*) )
  clrdep=FIX( (4.0-ALOG10(FLOAT(cdeps)))*31.5 )
  marks=FIX( (4.0-ALOG10(FLOAT(depth)))*31.5 )

;; Key sympol as filled vertical tangular
 xsymb=[0, 0, 0.8, 0.8, 0]
 ysymb=[0, 5.5, 5.5, 0, 0]
 usersym, xsymb, ysymb, /fill

 xkeys=[findgen(mclr)/float(mclr-1)-1.52]*radius
 ykeys= findgen(mclr)*0.0+radius*0.88

;; Read in 33 buoys id lat lon data from buoy33.dat
   ndbuoy = 'DatGMC/buoy33.dat'
    nnumb = 33
   buoids = LONARR(nnumb)
   buolat = FLTARR(nnumb)
   buolon = FLTARR(nnumb)
   status = DC_READ_FREE( ndbuoy, buoids, buolat, buolon, /Column )

;; New Pole as projection direction (central point)
 file4=[ '0',  '1',  '2',  '3' ]
 plon4=[  0.0,  90.0,  180.0,  270.0 ]
 plat4=[ 23.5,  40.0,   40.0,   40.0 ]

;; Open two files to store converted quadrilaterals coordinates
;; They contain one loop of the 4 trapezoid corners (5 points)
;; x and y coornates (5 each) for each cell, respectively.
;; Equatorial lat lon are also stored to identify two hemispheres
 openw, 10, 'S6GSTLaL.dat'
 openw, 11, 'S6GSTSpX.dat'
 openw, 12, 'S6GSTSpY.dat'
 printf,10, nc, 2,  Format="(2I8)"
 printf,11, nc, 5,  Format="(2I8)"
 printf,12, nc, 5,  Format="(2I8)"

 sxcel=fltarr(5,nc+1)
 sycel=fltarr(5,nc+1)
 seqll=fltarr(2,nc)

;set_plot, 'PS'
;;  Repeat for 2x2 plots 
;; FOR kk=0, 3, 2 DO Begin
   kk=0
 
 plon= plon4(kk)
 plat= plat4(kk)

;PR,/ps,/landscape,/color,file="g25grdcxy"+file4(kk)+".ps", $
;  xsize=20.0, ysize=39.0, xoffset=2.5, yoffset=1.5

 SET_PLOT, 'PS'
;DEVICE, /Encaps, /Portrait, Bits_Per_Pixel=8, /Color, File="g6kstrgrd.eps",  $
 DEVICE,          /Portrait, Bits_Per_Pixel=8, /Color, File="smc625grd.ps",   $
         YSize=21.0, XSize=40.0, xoffset=1.0, yoffset=3.0

;;*************** Set colours                *******************
  RESTORE_COLOURS,"/home/h05/frjl/PVWave/palettes/linuxspectrum.clr"

;; Set back and fore ground color to be white and black, respectively
    tvlct, 255, 255, 255,  !P.Background
; PRINT, ' !P.Background=',!P.Background
    tvlct,  0,  0,  0, !P.Color
; PRINT, ' !P.Color=', !P.Color

;***************  *******************  *******************
  !p.charthick = 2  &  !p.charsize = 1.0

   !P.Multi = [0, 2, 0, 0, 0]

;  Northern hemisphere

   plot, xcirc, ycirc, xrange=[-10,10], yrange=[-10,10],  $
       xstyle=4, ystyle=4, thick=2, $
       xmargin=[2,1], ymargin=[1,1]
;      title=' GMC Grid North NC='+string(nc, Format='(I6)')+ $
;            ' N4='+string(n4, Format='(I4)')
;            ' N2='+string(n2, Format='(I5)')+ $
;            ' N1='+string(n1, Format='(I6)')

;; Exclude last cell, the North Polar Cell
 for i=0L,nc-2L do BEGIN

   xc=[cel(0,i),cel(0,i)+cel(2,i),cel(0,i)+cel(2,i),cel(0,i),cel(0,i)]
   yc=[cel(1,i),cel(1,i),cel(3,i)+cel(1,i),cel(3,i)+cel(1,i),cel(1,i)]
   slat=ylat(yc+jeqt)
   slon=xlon(xc)

;; Convert slat slon to elat elon with given new pole

   LL2EqDeg, slat, slon, elat, elon, PoLat=Plat, PoLon=Plon

   if(elat(0) ge 0.0) then BEGIN
     pradmp=PRadus*cos(elat*d2rad)/(np - 1.0 + sin(elat*d2rad))
     syc = pradmp*cos(elon*d2rad)
     sxc =-pradmp*sin(elon*d2rad)
   endif else BEGIN
     pradmp=PRadus*cos(elat*d2rad)/(np - 1.0 - sin(elat*d2rad))
     syc = pradmp*cos(elon*d2rad)
     sxc = pradmp*sin(elon*d2rad)
   endelse

;  syc= radius*cos(elat*d2rad)*cos(elon*d2rad)
;  sxc0=-radius*cos(elat*d2rad)*sin(elon*d2rad)
;  if(elat(0) ge 0.0) then sxc= sxc0
;  if(elat(0) lt 0.0) then sxc=-sxc0

   seq=[elat(0), elon(0)]

   sxcel(*,i)=sxc
   sycel(*,i)=syc
   seqll(*,i)=seq

;; Store converted x and y values
   printf,10, seq,  Format="(2F9.3)"
   printf,11, sxc,  Format="(5F9.4)"
   printf,12, syc,  Format="(5F9.4)"

   if(elat(0) ge 0) then Begin
;  oplot, sxc, syc, Color=0
   oplot, sxc, syc, Color=colrs(clrdep(i))
   xyouts, (sxc(0)+sxc(2))*0.5, (syc(0)+syc(2))*0.5, '.', $
           charsize=0.005*cel(3,i)*Abs(sin(elat(0)*d2rad)), alignment=0.5, Color=colrs(clrdep(i)) 
;; Mark the arctical map-east reference region by first ring of its boundary cells at jmxglb - 3*4
     if(cel(1,i) eq jmxglb - 12) then polyfill, sxc, syc, Color=236
;; Mark the end of local-east reference region by last row of global part at jmxglb
     if(cel(1,i) eq jmxglb     ) then polyfill, sxc, syc, Color=136
   endif
 endfor

;;  Polar cell as a octagon, note polar cell size set as size-128 for flux calculation
;;  As it mergers 8 size-128 cells, each side of the octagon should be size-128
;;  10 apexes are calculated to conform with other cells.  To plot it, drop the last apex.
   i=nc-1
   xc=cel(0,i)+indgen(9)*cel(2,i)
   yc=xc*0+cel(1,i)
   slat=ylat(yc+jeqt)
   slon=xlon(xc)

;; Convert slat slon to elat elon with given new pole

   LL2EqDeg, slat, slon, elat, elon, PoLat=Plat, PoLon=Plon

   if(elat(0) ge 0.0) then BEGIN
     pradmp=PRadus*cos(elat*d2rad)/(np - 1.0 + sin(elat*d2rad))
     syc = pradmp*cos(elon*d2rad)
     sxc =-pradmp*sin(elon*d2rad)
   endif else BEGIN
     pradmp=PRadus*cos(elat*d2rad)/(np - 1.0 - sin(elat*d2rad))
     syc = pradmp*cos(elon*d2rad)
     sxc = pradmp*sin(elon*d2rad)
   endelse

;  syc= radius*cos(elat*d2rad)*cos(elon*d2rad)
;  sxc0=-radius*cos(elat*d2rad)*sin(elon*d2rad)
;  if(elat(0) ge 0.0) then sxc= sxc0
;  if(elat(0) lt 0.0) then sxc=-sxc0

   seq=[elat(0), elon(0)]

   seqll(*,i)=seq
   sxcel(*,i)=sxc(0:4)
   sycel(*,i)=syc(0:4)
   sxcel(*,i+1)=sxc(4:8)
   sycel(*,i+1)=syc(4:8)

;; Store converted x and y values
   printf,10, seq,  Format="(2F9.3)"
   printf,11, sxc(0:4),  Format="(5F9.4)"
   printf,12, syc(0:4),  Format="(5F9.4)"
   printf,11, sxc(4:8),  Format="(5F9.4)"
   printf,12, syc(4:8),  Format="(5F9.4)"

   if(elat(0) ge 0) then Begin
   oplot, sxc, syc, Color=colrs(clrdep(i))
   xyouts, (sxc(0)+sxc(4))*0.5, (syc(0)+syc(4))*0.5, '.', $
           charsize=1.0*Abs(sin(elat(0)*d2rad)), alignment=0.5, Color=254
;          charsize=0.6*Abs(sin(elat(0)*d2rad)), alignment=0.5, Color=colrs(clrdep(i))
   endif

;;  Convert buoy locations and overlay on grid map
   openw, 15, 'S6GSBuoy.dat'
   printf,15, nnumb, 7,  Format="(2I8)"
;   nnumb = 33
;  buoids = LONARR(nnumb)
   buosxc = FLTARR(nnumb)
   buosyc = FLTARR(nnumb)
   slat = buolat
   slon = buolon

   LL2EqDeg, slat, slon, elat, elon, PoLat=Plat, PoLon=Plon

   buelat = elat
   buelon = elon

  for i=0, nnumb-1 do BEGIN
   if(elat(i) ge 0.0) then BEGIN
     pradmp=PRadus*cos(elat(i)*d2rad)/(np - 1.0 + sin(elat(i)*d2rad))
     syc = pradmp*cos(elon(i)*d2rad)
     sxc =-pradmp*sin(elon(i)*d2rad)
;; Mark buoy position on map
     txtsiz=Abs(sin(elat(i)*d2rad))
     xyouts, sxc, syc, '.', charsize=4.0*txtsiz, alignment=0.5, Color=234
     xyouts, sxc, syc, 'Y', charsize=2.0*txtsiz, alignment=0.5, Color=234
   endif else BEGIN
     pradmp=PRadus*cos(elat(i)*d2rad)/(np - 1.0 - sin(elat(i)*d2rad))
     syc = pradmp*cos(elon(i)*d2rad)
     sxc = pradmp*sin(elon(i)*d2rad)
   endelse

     buosxc(i)=sxc
     buosyc(i)=syc

     printf,15, buoids(i), slat(i), slon(i), elat(i), elon(i), sxc, syc,  Format="(I8,4F9.3,2F9.4)"

  endfor
  
  close, 15


;  Southern hemisphere

   plot, xcirc, ycirc, xrange=[-10,10], yrange=[-10,10],  $
       xstyle=4, ystyle=4, thick=2, $
       xmargin=[1,2], ymargin=[1,1]

;; Exclude last cell, the North Polar Cell
 for i=0L,nc-2L do BEGIN

   seq=seqll(*,i)

   if(seq(0) lt 0) then Begin
   sxc=sxcel(*,i)
   syc=sycel(*,i)

;  oplot, sxc, syc, Color=0
   oplot, sxc, syc, Color=colrs(clrdep(i))
   xyouts, (sxc(0)+sxc(2))*0.5, (syc(0)+syc(2))*0.5, '.', $
           charsize=0.01*cel(3,i)*Abs(sin(elat(0)*d2rad)), alignment=0.5, Color=colrs(clrdep(i))
   endif
 endfor

;; Mark buoy position on map if visible
 for i=0, nnumb-1 do BEGIN
   if(buelat(i) lt 0.0) then BEGIN
     txtsiz=Abs(sin(buelat(i)*d2rad))
     xyouts, buosxc(i), buosyc(i), '.', charsize=4.0*txtsiz, alignment=0.5, Color=234
     xyouts, buosxc(i), buosyc(i), 'Y', charsize=2.0*txtsiz, alignment=0.5, Color=234
   endif
 endfor

;;  Put cell information inside plot
   xyouts, -10.2, -10.0, '!5Global SMC 6-25 km Grid',  $
           Alignment=0.5, Color=0,   Charthick=2,CharSize=1.5
   xyouts, -10.2, - 9.0, 'NC='+string(nc, Format='(I7)'),$
           Alignment=0.5, Color=254, Charthick=2,CharSize=1.3
   xyouts, -10.2, - 8.4, 'N1='+string(n1, Format='(I7)'),$
           Alignment=0.5, Color=254, Charthick=2,CharSize=1.3
   xyouts, -10.2, - 7.8, 'N2='+string(n2, Format='(I7)'),$
           Alignment=0.5, Color=254, Charthick=2,CharSize=1.3
   xyouts, -10.2, - 7.2, 'N4='+string(n4, Format='(I7)'),$
           Alignment=0.5, Color=254, Charthick=2,CharSize=1.3
;  xyouts, -10.2, - 6.6, 'N8='+string(n8, Format='(I5)'),$
;          Alignment=0.5, Color=254, Charthick=2,CharSize=1.3
;  xyouts, -10.2, - 6.0, 'N16='+string(n9-169, Format='(I4)'),$
;          Alignment=0.5, Color=254, Charthick=2,CharSize=1.3
;  xyouts, -10.2, - 5.4, 'N32='+string(128, Format='(I4)'),$
;          Alignment=0.5, Color=254, Charthick=2,CharSize=1.3

;  xyouts, -10.2,   5.8, 'N64='+string(32, Format='(I3)'),$
;          Alignment=0.5, Color=254, Charthick=2,CharSize=1.3
;  xyouts, -10.2,   6.4, 'N128='+string(8, Format='(I2)'),$
;          Alignment=0.5, Color=254, Charthick=2,CharSize=1.3
;  xyouts, -10.2,   7.0, 'NPol='+string(1, Format='(I2)'),$
;          Alignment=0.5, Color=254, Charthick=2,CharSize=1.3

;; Add a color key
 plots, xkeys, ykeys, psym=8, color=colrs, symsize=1.0
 plots, xkeys(marks), ykeys(marks), psym=8, color=colrs(marks), symsize=1.2
 for j=0,4 do xyouts, xkeys(marks(j)), ykeys(marks(j))+0.8, $
         string(fix(depth(j)), Format='(I5)'), charsize=1.0, $
         charthick=2.0, color=0, alignment=0.8 
 xyouts, -10.2, 8.2, 'Depth (m)', Alignment=0.5, Color=0, Charthick=2,CharSize=1.2

;PREND,'B2_CL13',/noprint,/view,/keep

;; End of kk loop
;; Endfor

 DEVICE, /Close

 end