;Annualaverage_MOPITTCO_nyg_modified.ncl
;
; This code loads MOPITT CO measurements
; and applies averaging kernels to model results
; an annual average in 2013
;
; To use type at command line:
;        ncl Annualaverage_MOPITTCO_rrb_modified.ncl
;
;                                            rrb 20201210
;********************************************************
undef("MOP_to_MDL_order")
function MOP_to_MDL_order(fMOP:list, varMOP[1]:string)  ; function MOP_to_MDL_order(fMOP:list, varMOP[5]:string)
begin

  x = fMOP[:]->$varMOP$
printVarSummary(x)
 psurf = fMOP[:]->SurfacePressureDay_MOP03   ; psurf = fMOP[:]->$varMOP$
printVarSummary(psurf)
 cocol_a_priori = fMOP[:]->APrioriCOTotalColumnDay_MOP03
printVarSummary(cocol_a_priori)
  coprof_a_priori = fMOP[:]->APrioriCOMixingRatioProfileDay_MOP03
printVarSummary(coprof_a_priori)
  avkernel = fMOP[:]->TotalColumnAveragingKernelDay_MOP03
printVarSummary(avkernel)
; MOPITT are ordered (12,lon,lat)      : lon==>XDim_MOP03 , lat==>YDim_MOP03
; match model dimension ordering (...,lat,lon)
; Use := syntax to replace the variable

  x := x(ncl_join|:,YDim_MOP03|:,XDim_MOP03|:)
  printVarSummary(x)   ; [ncl_join | 12] x [YDim_MOP03 | 180] x [XDim_MOP03 | 360]
  print("----------")
  psurf := psurf(ncl_join|:,YDim_MOP03|:,XDim_MOP03|:)
  printVarSummary(psurf)  ; [ncl_join | 12] x [YDim_MOP03 | 180] x [XDim_MOP03 | 360]
  cocol_a_priori := cocol_a_priori(ncl_join|:,YDim_MOP03|:,XDim_MOP03|:)
  printVarSummary(cocol_a_priori)
;  coprof_a_priori := coprof_a_priori(ncl_join|:,YDim_MOP03|:,XDim_MOP03|:,Prs_MOP03|:)
  coprof_a_priori := coprof_a_priori(ncl_join|:,YDim_MOP03|:,Prs_MOP03|:,XDim_MOP03|:)
  printVarSummary(coprof_a_priori)
;  avkernel := avkernel(ncl_join|:,YDim_MOP03|:,XDim_MOP03|:,Prs1_MOP03|:)
  avkernel := avkernel(ncl_join|:,YDim_MOP03|:,Prs1_MOP03|:,XDim_MOP03|:)
  printVarSummary(avkernel)
; model longitudes are 0->358.75
; make MOPITT longitudes be 0.5->359.5 ; convenience / consistency

  x = lonFlip(x)
  printVarSummary(x)   ; [ncl_join | 12] x [YDim_MOP03 | 180] x [XDim_MOP03 | 360]
  print("----------")

  psurf = lonFlip(psurf)
  printVarSummary(psurf)

  cocol_a_priori = lonFlip(cocol_a_priori)
;  printVarSummary(cocol_a_priori)

  coprof_a_priori = lonFlip(coprof_a_priori)
  printVarSummary(coprof_a_priori)

  avkernel = lonFlip(avkernel)
  printVarSummary(avkernel)

; rename ncl_join ==> time ..............; for other variables it could not display the summary and no error.
  x!0 = "time"
  printVarSummary(x)   ; [time | 12] x [YDim_MOP03 | 180] x [XDim_MOP03 | 360]
  print("----------")
   return(x)

  psurf!0 = "time"
  printVarSummary(psurf)
;   return(psurf)

  cocol_a_priori!0 = "time"
  printVarSummary(cocol_a_priori)
;   return(cocol_a_priori)

  coprof_a_priori!0 = "time"
  printVarSummary(coprof_a_priori)
;   return(coprof_a_priori)

  avkernel!0 = "time"
  printVarSummary(avkernel)
;   return(avkernel)

end

;***********************************************************************;
;                       SPECIAL
; The following is only for Najib's usage of coordinates subscripting
;***********************************************************************;
undef("MOPITT_extra_XDIM")
function MOPITT_extra_XDIM(data[*][*][*]:numeric)
local dims, newdata, nt, ny, mx, mx1, lon0, lon1, delta, lon_Last, lonExtra, \
      XDIM, XDIM1
begin
    dims    = dimsizes(data)
    nt      = dims(0)
    ny      = dims(1)
    mx      = dims(2)
    mx1     = mx+1

    dimNames= getvardimnames(data)
   ;print(dimNames)    ; ncl_join, YDim_MOP03, XDim_MOP03

    newdata = new((/nt,ny, mx1/),typeof(data))
    newdata(:,:,0:mx-1) = data               ; pass everything
    newdata(:,:,mx)     = (/ data(:,:,0) /)  ; values only
    newdata!0           = dimNames(0)        ; ncl_join
    newdata!1           = dimNames(1)        ; YDim_MOP03
    newdata!2           = dimNames(2)        ; XDim_MOP03
    printVarSummary(newdata)
    print("+++++++++++++++")

 ;; newdata&$dimNames(0)$ = data&$dimNames(1)$  ; there are no 'time'/'ncl_join' coordinates
    newdata&$dimNames(0)$ = ispan(1,nt,1)       ; time coordinate filler (not necessary)
    newdata&$dimNames(1)$ = data&$dimNames(1)$  ; newdata&YDim_MOP03  ... latitude coordinate variable

; add extra 'bogus' longitude point so that coordinate subscripting can be used

    lon0         = data&XDim_MOP03(0)
    lon1         = data&XDim_MOP03(1)
    delta        = lon1-lon0
    lonLast      = data&XDim_MOP03(mx-1)
    lonExtra     = lonLast+delta
    XDIM         = data&XDim_MOP03

    XDIM1        = new(mx1,typeof(XDIM))
    XDIM1(0:mx-1)= (/ XDIM /)
    XDIM1(  mx  )= lonExtra
    copy_VarAtts(XDIM,XDIM1)
    print(XDIM1)
    print("lon0="+lon0+"   lon1="+lon1+"  delta="+delta+"   lonLast="+lonLast+"  lonExtra="+lonExtra)

  ;;newdata&$dimNames(1)$ = XDIM1
    newdata&XDim_MOP03    = XDIM1
;    psurf&XDim_MOP03 = XDIM1
;       printVarSummary(psurf)
    printVarSummary(newdata)

    exit

    end

 begin

;-------------------------------------------------------------------------
; Setup
;-------------------------------------------------------------------------
; mopitt input

  wkdir  = "/glade/scratch/najiby/MOPITT_temp/"
  diri   = "/glade/scratch/najiby/MOPITT_temp/"
; model input
  dir_mo = "/glade/work/najiby/f.e20.FCSD_DICE_WAFRICA/f.e20.FCSD_DICE_WAFRICA/atm/hist/"

;-------------------------------------------------------------------------
;-- read MOPITT data and define
;--------------------------------------------
; find all MOP data:
  all_mop= systemfunc("ls -1 "+diri+"MOP03JM*")
  print(all_mop)
  nmop = dimsizes(all_mop)
  date_all = new(nmop,"string")

;--------------------------------------------
; open all MOPITT files
;--------------------------------------------
fileIn = all_mop

; Read files and collect variables

  fMOP = addfiles(fileIn,"r")
  ListSetType (fMOP, "join")           ; concatenate (=default)
  printVarSummary(fMOP)

  vnames = getfilevarnames(fMOP[0])
;; print(vnames)
  lat_mop = fMOP[0]->YDim_MOP03        ;  [-89.5..89.5]
  lat_mop@long_name = "Latitude"       ; latitude  ==> Latitude
  lon_mop = fMOP[0]->XDim_MOP03        ;  [-179.5..179.5]
  lon_mop@long_name = "Longitude"      ; longitude ==> Longitude

 x  = MOP_to_MDL_order(fMOP, "RetrievedCOTotalColumnDay_MOP03")
    printVarSummary(x)
    print("----------")
  x := MOPITT_extra_XDIM(x)
  printVarSummary(x)
  print("----------")

  psurf = MOP_to_MDL_order(fMOP, "SurfacePressureDay_MOP03")
    printVarSummary(psurf)
  psurf := MOPITT_extra_XDIM(psurf)
   printVarSummary(psurf)

  cocol_a_priori = MOP_to_MDL_order(fMOP, "APrioriCOTotalColumnDay_MOP03")
    printVarSummary(cocol_a_priori)
  cocol_a_priori := MOPITT_extra_XDIM(cocol_a_priori)
  printVarSummary(cocol_a_priori)

  coprof_a_priori = MOP_to_MDL_order(fMOP, "APrioriCOMixingRatioProfileDay_MOP03")
    printVarSummary(coprof_a_priori)
  coprof_a_priori := MOPITT_extra_XDIM(coprof_a_priori)
  printVarSummary(coprof_a_priori)

  avkernel = MOP_to_MDL_order(fMOP, "TotalColumnAveragingKernelDay_MOP03")
    printVarSummary(avkernel)
  avkernel := MOPITT_extra_XDIM(avkernel)
  printVarSummary(avkernel)
  print("----------")

exit
;-------------------------------------------------------------------------
; Load and collect equivalent Model data
;-------------------------------------------------------------------------
  dir_mo  = "/glade/work/najiby/f.e20.FCSD_DICE_WAFRICA/f.e20.FCSD_DICE_WAFRICA/atm/hist/"
  file_mo = systemfunc("ls -1 "+dir_mo+"f.e20.FCSD_DICE_WAFRICA.cam.h0.2013*.nc")

  print(file_mo)
  filein = file_mo
  fMDL = addfiles(filein, "r")     ; MODEL
  ListSetType (fMDL, "cat")        ; concatenate (=default)

  vnames1 = getfilevarnames(fMDL[0])
; print(vnames1)

  time     = fMDL[:]->time
  lev1     = fMDL[0]->lev
  lat1     = fMDL[0]->lat
  lon1     = fMDL[0]->lon
  ntim1    = dimsizes(time)
  nlat1    = dimsizes(lat1)
  nlon1    = dimsizes(lon1)
  nlev1    = dimsizes(lev1)
  co_model = fMDL[:]->CO
  printVarSummary(co_model)
  print("----------")

  mod_yyyymmdd = cd_calendar(co_model&time, 2)
  model_lev = co_model&lev
  nlev1 = dimsizes(model_lev)

; load level information for interpolation

  ps   = fMDL[:]->PS
  hyam = fMDL[0]->hyam
  hybm = fMDL[0]->hybm
  p0   = fMDL[0]->P0
  mod_press = pres_hybrid_ccm(ps,p0,hyam,hybm)
  copy_VarCoords(co_model, mod_press)
  printVarSummary(mod_press)
  print("----------")

; initialize array to collect model data

co_model_select = new((/dimsizes(time), nlev1/),float)
altp_model_select = new((/dimsizes(time), nlev1/),float)
co_model_smooth = new(dimsizes(newdata),float)
ln_co_mo_ak = new(nlev1,float)

; natural log of a priori
ln_cocol_a_priori = log(cocol_a_priori)

; MOPITT Prior
MOPITT_ap_annavg = dim_avg_n(cocol_a_priori,0)

; Loop through Measurement data and find equivalent model data and Collect base model data

  time_collect = fMOP[:]@StartTime_MOP03
      ;units value presumes use of TAI93 (International Atomic Time) format
      time_collect@units = "seconds since 1993-1-1 00:00:00"
  date_collect  := cd_calendar(time_collect, 0)

;-- Collect equivalent profile
  do n=0,nmop-1
  time_ind = ind(mod_yyyymmdd.eq.time_collect(n))
   co_model_select(n,:) = co_model(time_ind:,{lat_mop(n)},{lon_mop(n)})
   altp_model_select(n,:) = mod_press(time_ind,:,{lat_mop(n)},{lon_mop(n)})
  end do

; Regrid to MOPITT grid
   co_mo_regrid = int2p_n_Wrap(altp_model_select(n,:)/100,co_model_select(n,:),psurf(n,:),0,0)
   ln_co_mo_regrid = log(co_mo_regrid)

; Read AK
;  avkernel = fMOP[:]->TotalColumnAveragingKernelDay_MOP03
;  coprof_a_priori = fMOP[:]->APrioriCOMixingRatioProfileDay_MOP03

;-- Apply averaging kernel
;--------------------------------------------
; Need to use the usual multiplication because matrix multiplication can't handle missing values
   do m=0,nlev1-1
     ln_co_mo_ak(m) = ln_cocol_a_priori(n,m) + sum((avkernel(n,m,:)) * (ln_co_mo_regrid - ln_cocol_a_priori(n,:)))
   end do
   co_mo_ak_dummy = exp(ln_co_mo_ak)
   co_mo_ak = where(ismissing(co_mo_regrid), co_mo_ak_dummy@_FillValue, co_mo_ak_dummy)

; Collect smoothed model data
do n=0,nmop-1
 co_model_smooth(n,:) = co_mo_ak
end do ; end go through all MOPITT data profiles and create smoothed model value

; Perform Annual Averages
; MOPITT
  MOPITT_co_annavg = dim_avg_n(newdata,0)
  MOPITT_co_annavg!0 = "lon"
  MOPITT_co_annavg&lon = lon_mop
  MOPITT_co_annavg!1 = "lat"
  MOPITT_co_annavg&lat = lat_mop
  MOPITT_co_annsd = dim_stddev_n(newdata,0)
  MOPITT_alt_annavg = dim_avg_n(psurf,0)

printVarSummary(MOPITT_co_annavg)
printVarSummary(MOPITT_co_annsd)
; MODEL
 model_co_model_annavg = dim_avg_n(co_model,0)
 model_co_model_annsd = dim_stddev_n(co_model,0)
 model_altp_annavg = dim_avg_n(lev1,0)/100

printVarSummary(model_co_model_annavg)
printVarSummary(model_co_model_annsd)
; convert pressure to meters

convert_altp_annavg = stdatmus_p2tdz(model_altp_annavg)
model_alt_annavg = convert_altp_annavg(2,:)

; Smoothed Model
smooth_model_co_annavg = dim_avg_n(co_model_smooth,0)
smooth_model_co_annsd = dim_stddev_n(co_model_smooth,0)


end