;
;  mkscripgrid.ncl
;
;  Create SCRIP grid and mapping file for a land-only point or region.
;  Requires NCL 6.1.0 or later for the ESMF regridding functions
;
;  Erik Kluzek
;  Dec/07/2011
;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
begin
  ; ===========================================================================================================
  ; Set a few constants needed later
  cdate  = systemfunc( "date +%y%m%d" );
  ldate  = systemfunc( "date" );
  ;
  ; IMPORTANT NOTE: EDIT THE FOLLOWING TO CUSTOMIZE or use ENV VARIABLE SETTINGS
  ; Edit the following as needed to interpolate to a new resolution.
  ;
  ; Input resolution and position
  ;
  name = getenv("PTNAME");                    ; Get name of this point

  latS = stringtodouble( getenv("S_LAT") );   ; Get south   latitude from env variable
  latN = stringtodouble( getenv("N_LAT") );   ; Get north   latitude from env variable
  lonE = stringtodouble( getenv("E_LON") );   ; Get east   longitude from env variable
  lonW = stringtodouble( getenv("W_LON") );   ; Get west   longitude from env variable

  nx   = stringtointeger( getenv("NX" )  );   ; Get number of grids along longitude lines
  ny   = stringtointeger( getenv("NY" )  );   ; Get number of grids along latitude lines

  imask = stringtointeger( getenv("IMASK") ); ; Get imask to use     from env variable

  print_str  = getenv("PRINT");               ; Do Extra printing for debugging

  outfilename = getenv("GRIDFILE");           ; Get filename from env variable

  gitdescribe = getenv("GITDES");             ; Git describe from the source clone

  if ( ismissing(nx) )then
     nx = 1;
  end if
  if ( ismissing(ny) )then
     ny = 1;
  end if
  if ( ismissing(imask) )then
     imask = 1;
  end if
  if ( ismissing(name) )then
     name = nx+"x"+ny+"pt_US-UMB";
  end if
  if ( ismissing(latS) )then
     latS = 45.5098;
  end if
  if ( ismissing(latN) )then
     latN = 45.6098;
  end if
  if ( ismissing(lonW) )then
     lonW = 275.2362;
  end if
  if ( ismissing(lonE) )then
     lonE = 275.3362;
  end if
  if ( ismissing(print_str) )then
     printn = False;
  else
     if ( print_str .eq. "TRUE" )then
        printn = True;
     else
        printn = False;
     end if
  end if

  if ( ismissing(outfilename) )then
     if ( imask .eq. 1 )then
        outfilename = "SCRIPgrid_"+name+"_nomask_c"+cdate+".nc";
     else
        if ( imask .eq. 0 )then
           outfilename = "SCRIPgrid_"+name+"_noocean_c"+cdate+".nc";
        else
           outfilename = "SCRIPgrid_"+name+"_mask_c"+cdate+".nc";
        end if
     end if
  end if

  if ( ismissing(gitdescribe) )then
     gitdescribe = systemfunc( "git describe" )
  end if

  system( "/bin/rm -f "+outfilename );
  if ( printn )then
     print( "output file: "+outfilename );
  end if

function fspan1up( fbegin [*]:double, fend [*]:double, number:integer )
;
; An "fspan" that can handle size of 1 and up.
; Do fspan for arrays of two or more, or average of end points for array of one.
;
local farray;
begin
   if ( number .eq. 1) then
      farray = (/ (fbegin + fend) / 2.0d00 /);
   else
      farray = fspan( fbegin, fend, number );
   end if
   return( farray );
end

  ;
  ; Compute derived quantities
  ;

  delX       = (lonE - lonW) / int2dble(nx);
  delY       = (latN - latS) / int2dble(ny);

  lonCenters = fspan1up( (lonW + delX/2.d0), (lonE - delX/2.d0), nx)
  latCenters = fspan1up( (latS + delY/2.d0), (latN - delY/2.d0), ny)
  lon = new( (/ny, nx/), "double" );
  lat = new( (/ny, nx/), "double" );
  if ( (nx .eq. 1) .or. (ny .eq. 1) )then
     if ( printn )then
        print( "Calculate corners" )
     end if
     lonCorners = new( (/ny, nx, 4/), "double" );
     latCorners = new( (/ny, nx, 4/), "double" );
  else
     if ( printn )then
        print( "Have NCL calculate corners" )
     end if
  end if
  do i = 0, nx-1
    lat(:,i)          =  latCenters;
    if ( (nx .eq. 1) .or. (ny .eq. 1) )then
       latCorners(:,i,0) =  latCenters - delY/2.d0;
       latCorners(:,i,1) =  latCenters - delY/2.d0;
       latCorners(:,i,2) =  latCenters + delY/2.d0;
       latCorners(:,i,3) =  latCenters + delY/2.d0;
    end if
  end do
  do j = 0, ny-1
    lon(j,:)          =  lonCenters;
    if ( (nx .eq. 1) .or. (ny .eq. 1) )then
       lonCorners(j,:,0) =  lonCenters - delX/2.d0;
       lonCorners(j,:,1) =  lonCenters + delX/2.d0;
       lonCorners(j,:,2) =  lonCenters + delX/2.d0;
       lonCorners(j,:,3) =  lonCenters - delX/2.d0;
    end if
  end do

  ; for some reason, "No_FillValue" isn't working in the case where imask=1
  Mask2D = new( (/ny,nx/), "integer", "No_FillValue" )
  Mask2D(:,:) = imask
  gridSize = delX+"x"+delY

  ;
  ; Create SCRIP grid file
  ;
  
  Opt = True
  Opt@Mask2D   = Mask2D
  if ( (nx .eq. 1) .or. (ny .eq. 1) )then
    Opt@GridCornerLat = latCorners
    Opt@GridCornerLon = lonCorners
  end if
  Opt@Title = "SCRIP grid file for "+name
  if (printn) then
    Opt@Debug = True
  end if
  curvilinear_to_SCRIP(outfilename, lat, lon, Opt)

  ;
  ; Add global attributes to file
  ;

  nc = addfile( outfilename, "w" );
  nc@history = ldate+": create using mkscripgrid.ncl";
  nc@comment = "Ocean is assumed to be non-existant in this region";
  nc@Version = gitdescribe;
  if ( printn )then
    print( "================================================================================================" );
    print( "Successfully created SCRIP grid file: "+outfilename);
  end if

  ; ===========================================================================================================

end