; Computes LAI and CHLR arrays per land type from MODIS data and Olson landmap.
; Generates twelve monthly LAI files per year, and twelve monthly CHLR
; files per year that is available.
; Each file contains 73 variables, one for each land type. Values are 0 
; except where that land type is present, in which case the original LAI
; or CHLR data is used. 
; (L. Lundgren, 11/3/16)

; NCL routines
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"  

begin

; Read Olson variables
inFile_olson = "/n/holylfs/EXTERNAL_REPOS/GEOS-CHEM/gcgrid/gcdata/GEOS_NATIVE/Olson_Land_Map_201203/Olson_2001_Land_Map.025x025.generic.nc"
fIn_olson = addfile( inFile_olson, "r" )
vNames_olson = getfilevarnames(fIn_olson)
olson = fIn_olson->OLSON
dims_olson  = dimsizes(olson)
nTime_olson = dims_olson(0)
nLat_olson  = dims_olson(1)
nLon_olson  = dims_olson(2)

;=========================================================================
; LAI or CHLR, loop over files (one per year
;=========================================================================
do datavar = 0, 0    ; 0 == LAI, 1 == CHLR
do yr = 2000, 2020

;  ; if LAI, only go from 2005 to 2016
;  if ( datavar .eq. 0 .and. yr .gt. 2016 )
;    continue
;  end if

  ; Source data ( 0 = LAI, 1 = CHLR )
  if ( datavar .eq. 0 )
    ; for Odyssey:
    inFile_modis   = "../MODIS.LAI.BNUv6.generic.025x025."+yr+".nc"
    print(inFile_modis)
  else
    ; for Odyssey:
    inFile_modis   = "/n/holylfs/EXTERNAL_REPOS/GEOS-CHEM/gcgrid/gcdata/GEOS_NATIVE/MODIS_LAI_201204/For_Olson_2001/MODIS.CHLRv.V5.generic.025x025."+yr+".nc"
    print(inFile_modis)
  end if    
  
  ; Open and inspect input files
  fIn_modis = addfile( inFile_modis, "r" )
  
  ; Read MODIS data
  modis      = fIn_modis->MODIS
  lon_modis  = fIn_modis->lon
  lat_modis  = fIn_modis->lat
  time_modis = fIn_modis->time
  
  ; Get dimensions of MODIS data
  dims_modis  = dimsizes(modis)
  nTime = dims_modis(0)
  nLat  = dims_modis(1)
  nLon  = dims_modis(2)
  
  ; Debugging (lai)
  ;vNames_modis = getfilevarnames(fIn_modis)
  ;print(vNames_modis)
  ;print(time_modis)
  ;print(dims_modis)
  ;print(dims_olson)
  ;print(nTime_modis)
  ;print(nLat_modis)
  ;print(nLon_modis)
  ;print(modis(1,719,1439))
  ;print(olson(0,719,1439))

  ;=========================================================================
  ; Loop over months and set month-specific information
  ; NOTE: There will be one output file per month in order to fit the data
  ;=========================================================================
  do iTime = 0, nTime-1

    nTimeM1 = nTime-1
    print("Processing slice " +iTime+ " of " +nTimeM1)

    ; Number of days
    iDay = 1 + ( time_modis(iTime) / 1440 )

    if ( iDay .lt. 10 ) then
       day = "00"+iDay 
    else if ( iDay .lt. 100 ) then
       day = "0"+iDay
    else
       day = ""+iDay
    end if
    end if 
    print(day)

    if ( datavar .eq. 0 )
      outFile = "../v2021-06/"+yr+"/MODIS_XLAI_BNUv6.025x025."+yr+"_"+day+".nc"
    else
      outFile = "./XCHLR_for_GCHP/"+yr+"/MODIS_XCHLR.025x025."+yr+mo+".nc"
    end if
    print(outFile)
  
    ; Skip this file if it already exists
    ; NOTE: disable this in future. It is used since need to create
    ; the files across multiple slurm runs, but don't want to customize
    ; based on what's already been created
    if ( isfilepresent( outFile ) )
      continue
    end if
  
    ;=======================================================================
    ; Define variables and attributes
    ;=======================================================================
    
    ;; Copy lat and lon	             
    lat = lat_modis
    lon = lon_modis
  
    ;%%% Create time dimension and index array %%%	             
    time_new                    = new( (/ 1 /), integer )
    time_new(0)                 = time_modis(iTime)
    time_new!0                  = "time"
    time_new@long_name          = "time"
    time_new@units              = "minutes since "+yr+"-01-01 00:00:00"
    time_new@_ChunkSizes        = 1
    time_new&time               = time_new
  
    ;%%% Create template variable for modis data
    modisxx_template            = new( (/ 1, nLat, nLon /), float )
    modisxx_template!0          = "time"
    modisxx_template!1          = "lat"
    modisxx_template!2          = "lon"
    modisxx_template&time       = time_new
    modisxx_template&lat        = lat
    modisxx_template&lon        = lon
    if ( datavar .eq. 0 )
       modisxx_template@long_name     = "MODIS BNU leaf area index"
       modisxx_template@units         = "cm2 cm-2"
    else
       modisxx_template@long_name     = "MODIS-Aqua Chlorophyll-a"
       modisxx_template@units         = "mg m-3"
    end if
    modisxx_template@_ChunkSizes      = (/ 1, nLat, nLon /)
  
    ;%%% Initialize new land-type-specific MODIS variables
    modis00 = modisxx_template
    modis01 = modisxx_template
    modis02 = modisxx_template
    modis03 = modisxx_template
    modis04 = modisxx_template
    modis05 = modisxx_template
    modis06 = modisxx_template
    modis07 = modisxx_template
    modis08 = modisxx_template
    modis09 = modisxx_template
    modis10 = modisxx_template
    modis11 = modisxx_template
    modis12 = modisxx_template
    modis13 = modisxx_template
    modis14 = modisxx_template
    modis15 = modisxx_template
    modis16 = modisxx_template
    modis17 = modisxx_template
    modis18 = modisxx_template
    modis19 = modisxx_template
    modis20 = modisxx_template
    modis21 = modisxx_template
    modis22 = modisxx_template
    modis23 = modisxx_template
    modis24 = modisxx_template
    modis25 = modisxx_template
    modis26 = modisxx_template
    modis27 = modisxx_template
    modis28 = modisxx_template
    modis29 = modisxx_template
    modis30 = modisxx_template
    modis31 = modisxx_template
    modis32 = modisxx_template
    modis33 = modisxx_template
    modis34 = modisxx_template
    modis35 = modisxx_template
    modis36 = modisxx_template
    modis37 = modisxx_template
    modis38 = modisxx_template
    modis39 = modisxx_template
    modis40 = modisxx_template
    modis41 = modisxx_template
    modis42 = modisxx_template
    modis43 = modisxx_template
    modis44 = modisxx_template
    modis45 = modisxx_template
    modis46 = modisxx_template
    modis47 = modisxx_template
    modis48 = modisxx_template
    modis49 = modisxx_template
    modis50 = modisxx_template
    modis51 = modisxx_template
    modis52 = modisxx_template
    modis53 = modisxx_template
    modis54 = modisxx_template
    modis55 = modisxx_template
    modis56 = modisxx_template
    modis57 = modisxx_template
    modis58 = modisxx_template
    modis59 = modisxx_template
    modis60 = modisxx_template
    modis61 = modisxx_template
    modis62 = modisxx_template
    modis63 = modisxx_template
    modis64 = modisxx_template
    modis65 = modisxx_template
    modis66 = modisxx_template
    modis67 = modisxx_template
    modis68 = modisxx_template
    modis69 = modisxx_template
    modis70 = modisxx_template
    modis71 = modisxx_template
    modis72 = modisxx_template

    ; Set values in new arrays to original MODIS data for this month
    ; NOTE: for some reason the regridded 025 x 025 data have
    ; the variables in double and not float.  Convert back here.
    modis00(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis01(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis02(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis03(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis04(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis05(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis06(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis07(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis08(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis09(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis10(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis11(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis12(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis13(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis14(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis15(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis16(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis17(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis18(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis19(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis20(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis21(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis22(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis23(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis24(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis25(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis26(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis27(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis28(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis29(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis30(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis31(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis32(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis33(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis34(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis35(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis36(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis37(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis38(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis39(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis40(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis41(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis42(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis43(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis44(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis45(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis46(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis47(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis48(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis49(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis50(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis51(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis52(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis53(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis54(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis55(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis56(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis57(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis58(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis59(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis60(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis61(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis62(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis63(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis64(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis65(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis66(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis67(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis68(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis69(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis70(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis71(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
    modis72(0,:,:) = doubletofloat(modis(iTime, 0:nLat-1, 0:nLon-1))
  
    ; Loop over latitudes and longitudes to apply land mask per type
    do ilat  = 0, nLat-1
      do ilon  = 0, nLon-1
    
        if ( olson(0,ilat,ilon) .ne. 0 )
          modis00(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 1 )
          modis01(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 2 )
          modis02(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 3 )
          modis03(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 4 )
          modis04(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 5 )
          modis05(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 6 )
          modis06(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 7 )
          modis07(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 8 )
          modis08(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 9 )
          modis09(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 10 )
          modis10(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 11 )
          modis11(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 12 )
          modis12(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 13 )
          modis13(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 14 )
          modis14(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 15 )
          modis15(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 16 )
          modis16(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 17 )
          modis17(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 18 )
          modis18(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 19 )
          modis19(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 20 )
          modis20(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 21 )
          modis21(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 22 )
          modis22(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 23 )
          modis23(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 24 )
          modis24(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 25 )
          modis25(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 26 )
          modis26(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 27 )
          modis27(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 28 )
          modis28(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 29 )
          modis29(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 30 )
          modis30(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 31 )
          modis31(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 32 )
          modis32(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 33 )
          modis33(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 34 )
          modis34(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 35 )
          modis35(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 36 )
          modis36(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 37 )
          modis37(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 38 )
          modis38(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 39 )
          modis39(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 40 )
          modis40(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 41 )
          modis41(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 42 )
          modis42(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 43 )
          modis43(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 44 )
          modis44(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 45 )
          modis45(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 46 )
          modis46(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 47 )
          modis47(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 48 )
          modis48(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 49 )
          modis49(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 50 )
          modis50(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 51 )
          modis51(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 52 )
          modis52(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 53 )
          modis53(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 54 )
          modis54(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 55 )
          modis55(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 56 )
          modis56(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 57 )
          modis57(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 58 )
          modis58(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 59 )
          modis59(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 60 )
          modis60(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 61 )
          modis61(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 62 )
          modis62(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 63 )
          modis63(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 64 )
          modis64(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 65 )
          modis65(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 66 )
          modis66(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 67 )
          modis67(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 68 )
          modis68(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 69 )
          modis69(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 70 )
          modis70(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 71 )
          modis71(0,ilat,ilon) = 0.0
        end if
        if ( olson(0,ilat,ilon) .ne. 72 )
          modis72(0,ilat,ilon) = 0.0
        end if
    
      end do
    end do
  
    ;=========================================================================
    ; Save to output file for this month and year
    ;========================================================================= 
    print("Saving out "+outFile)
  
    ; Open output file (remove prior version if it exists)
    system( "rm -f " + outFile )
    fOut           = addfile( outFile, "c" )
    if ( datavar .eq. 0 )
       fOut->XLAI00  = modis00
       fOut->XLAI01  = modis01
       fOut->XLAI02  = modis02
       fOut->XLAI03  = modis03
       fOut->XLAI04  = modis04
       fOut->XLAI05  = modis05
       fOut->XLAI06  = modis06
       fOut->XLAI07  = modis07
       fOut->XLAI08  = modis08
       fOut->XLAI09  = modis09
       fOut->XLAI10  = modis10
       fOut->XLAI11  = modis11
       fOut->XLAI12  = modis12
       fOut->XLAI13  = modis13
       fOut->XLAI14  = modis14
       fOut->XLAI15  = modis15
       fOut->XLAI16  = modis16
       fOut->XLAI17  = modis17
       fOut->XLAI18  = modis18
       fOut->XLAI19  = modis19
       fOut->XLAI20  = modis20
       fOut->XLAI21  = modis21
       fOut->XLAI22  = modis22
       fOut->XLAI23  = modis23
       fOut->XLAI24  = modis24
       fOut->XLAI25  = modis25
       fOut->XLAI26  = modis26
       fOut->XLAI27  = modis27
       fOut->XLAI28  = modis28
       fOut->XLAI29  = modis29
       fOut->XLAI30  = modis30
       fOut->XLAI31  = modis31
       fOut->XLAI32  = modis32
       fOut->XLAI33  = modis33
       fOut->XLAI34  = modis34
       fOut->XLAI35  = modis35
       fOut->XLAI36  = modis36
       fOut->XLAI37  = modis37
       fOut->XLAI38  = modis38
       fOut->XLAI39  = modis39
       fOut->XLAI40  = modis40
       fOut->XLAI41  = modis41
       fOut->XLAI42  = modis42
       fOut->XLAI43  = modis43
       fOut->XLAI44  = modis44
       fOut->XLAI45  = modis45
       fOut->XLAI46  = modis46
       fOut->XLAI47  = modis47
       fOut->XLAI48  = modis48
       fOut->XLAI49  = modis49
       fOut->XLAI50  = modis50
       fOut->XLAI51  = modis51
       fOut->XLAI52  = modis52
       fOut->XLAI53  = modis53
       fOut->XLAI54  = modis54
       fOut->XLAI55  = modis55
       fOut->XLAI56  = modis56
       fOut->XLAI57  = modis57
       fOut->XLAI58  = modis58
       fOut->XLAI59  = modis59
       fOut->XLAI60  = modis60
       fOut->XLAI61  = modis61
       fOut->XLAI62  = modis62
       fOut->XLAI63  = modis63
       fOut->XLAI64  = modis64
       fOut->XLAI65  = modis65
       fOut->XLAI66  = modis66
       fOut->XLAI67  = modis67
       fOut->XLAI68  = modis68
       fOut->XLAI69  = modis69
       fOut->XLAI70  = modis70
       fOut->XLAI71  = modis71
       fOut->XLAI72  = modis72
    else
       fOut->XCHLR00  = modis00
       fOut->XCHLR01  = modis01
       fOut->XCHLR02  = modis02
       fOut->XCHLR03  = modis03
       fOut->XCHLR04  = modis04
       fOut->XCHLR05  = modis05
       fOut->XCHLR06  = modis06
       fOut->XCHLR07  = modis07
       fOut->XCHLR08  = modis08
       fOut->XCHLR09  = modis09
       fOut->XCHLR10  = modis10
       fOut->XCHLR11  = modis11
       fOut->XCHLR12  = modis12
       fOut->XCHLR13  = modis13
       fOut->XCHLR14  = modis14
       fOut->XCHLR15  = modis15
       fOut->XCHLR16  = modis16
       fOut->XCHLR17  = modis17
       fOut->XCHLR18  = modis18
       fOut->XCHLR19  = modis19
       fOut->XCHLR20  = modis20
       fOut->XCHLR21  = modis21
       fOut->XCHLR22  = modis22
       fOut->XCHLR23  = modis23
       fOut->XCHLR24  = modis24
       fOut->XCHLR25  = modis25
       fOut->XCHLR26  = modis26
       fOut->XCHLR27  = modis27
       fOut->XCHLR28  = modis28
       fOut->XCHLR29  = modis29
       fOut->XCHLR30  = modis30
       fOut->XCHLR31  = modis31
       fOut->XCHLR32  = modis32
       fOut->XCHLR33  = modis33
       fOut->XCHLR34  = modis34
       fOut->XCHLR35  = modis35
       fOut->XCHLR36  = modis36
       fOut->XCHLR37  = modis37
       fOut->XCHLR38  = modis38
       fOut->XCHLR39  = modis39
       fOut->XCHLR40  = modis40
       fOut->XCHLR41  = modis41
       fOut->XCHLR42  = modis42
       fOut->XCHLR43  = modis43
       fOut->XCHLR44  = modis44
       fOut->XCHLR45  = modis45
       fOut->XCHLR46  = modis46
       fOut->XCHLR47  = modis47
       fOut->XCHLR48  = modis48
       fOut->XCHLR49  = modis49
       fOut->XCHLR50  = modis50
       fOut->XCHLR51  = modis51
       fOut->XCHLR52  = modis52
       fOut->XCHLR53  = modis53
       fOut->XCHLR54  = modis54
       fOut->XCHLR55  = modis55
       fOut->XCHLR56  = modis56
       fOut->XCHLR57  = modis57
       fOut->XCHLR58  = modis58
       fOut->XCHLR59  = modis59
       fOut->XCHLR60  = modis60
       fOut->XCHLR61  = modis61
       fOut->XCHLR62  = modis62
       fOut->XCHLR63  = modis63
       fOut->XCHLR64  = modis64
       fOut->XCHLR65  = modis65
       fOut->XCHLR66  = modis66
       fOut->XCHLR67  = modis67
       fOut->XCHLR68  = modis68
       fOut->XCHLR69  = modis69
       fOut->XCHLR70  = modis70
       fOut->XCHLR71  = modis71
       fOut->XCHLR72  = modis72
    end if
    fOut->time     = time_new
    fOut->lat      = lat
    fOut->lon      = lon
    
    ; Add global file attributes
    ; Redefine the title string
    if ( datavar .eq. 0 )
      fOut@Title              = "MODIS BNU leaf area indices per land type" 
    else
      fOut@Title              = "MODIS-Aqua Chlorophyll-a" 
    end if
    fOut@Conventions          = "COARDS" 
    fOut@History              = "Last modified on: " + systemfunc( "date" )
    fOut@ProductionDateTime   = "Last modified on: " + systemfunc( "date" )
    fOut@ModificationDateTime = "Last modified on: " + systemfunc( "date" )
    fOut@format               = "netCDF-3" ;
    ;fOut@model                = "GEOS5" ;
    ;fOut@nlayers              = "72" ;
    ;fOut@start_date           = yr+mo+"15";
    ;fOut@start_time           = "00:00:00" ;
    ;fOut@end_date             = yr+mo+"15";
    ;fOut@end_time             = "23:59:59" ;
    ;fOut@delta_lon            = "0.25" ;
    ;fOut@delta_lat            = "0.25" ;
    ;fOut@delta_time           = "000000" ;
  
    ; clean up
    delete(lat)
    delete(lon)
    delete(time_new)
    delete(modisxx_template)
    delete(modis00)
    delete(modis01)
    delete(modis02)
    delete(modis03)
    delete(modis04)
    delete(modis05)
    delete(modis06)
    delete(modis07)
    delete(modis08)
    delete(modis09)
    delete(modis10)
    delete(modis11)
    delete(modis12)
    delete(modis13)
    delete(modis14)
    delete(modis15)
    delete(modis16)
    delete(modis17)
    delete(modis18)
    delete(modis19)
    delete(modis20)
    delete(modis21)
    delete(modis22)
    delete(modis23)
    delete(modis24)
    delete(modis25)
    delete(modis26)
    delete(modis27)
    delete(modis28)
    delete(modis29)
    delete(modis30)
    delete(modis31)
    delete(modis32)
    delete(modis33)
    delete(modis34)
    delete(modis35)
    delete(modis36)
    delete(modis37)
    delete(modis38)
    delete(modis39)
    delete(modis40)
    delete(modis41)
    delete(modis42)
    delete(modis43)
    delete(modis44)
    delete(modis45)
    delete(modis46)
    delete(modis47)
    delete(modis48)
    delete(modis49)
    delete(modis50)
    delete(modis51)
    delete(modis52)
    delete(modis53)
    delete(modis54)
    delete(modis55)
    delete(modis56)
    delete(modis57)
    delete(modis58)
    delete(modis59)
    delete(modis60)
    delete(modis61)
    delete(modis62)
    delete(modis63)
    delete(modis64)
    delete(modis65)
    delete(modis66)
    delete(modis67)
    delete(modis68)
    delete(modis69)
    delete(modis70)
    delete(modis71)
    delete(modis72)
    
  end do ; end of month loop
  
  ; clean up
  delete(lat_modis)
  delete(lon_modis)
  delete(time_modis)
  delete(inFile_modis)
  delete(modis)

end do   ; end of year loop
end do   ; end of type loop (LAI or CHLR)
end