A set of tool chain to analyze the CMEPS mediator history files and and converts 1d FV3 related fields in it to 2d representations and writes into a new file as tiles

plot_cubed.ncl

; Note for user:
; The orog.tile*.nc files can be generated by simply using following command
; for each tile of surface output.
; ncks -v orog sfcf000.tile1.nc orog.tile1.nc
begin
  ;--- arguments ---
  ntiles = 6
  did = 0
  diff = True
  mask_land = False

  ;--- datasets: label, directory, prefix, variable name, time (-1 for mediator history) ---
  ;--- the max number of dataset in the list needs to be two ---
  dsets = (/ (/ "MED", "./", "ufs.cpld.cpl.hi", "Med_aoflux_atm_Faox_sen", "2016-10-04", "0" /), \
             (/ "ATM", "./", "ufs.cpld.cpl.hi", "atmImp_Faxa_sen"        , "2016-10-04", "0" /) /)
  ndsets = dimsizes(dsets(:,0))
  
  ;--- list files for each dataset ---
  fl = new((/ ndsets, ntiles /), "string")
  do j = 0, ndsets-1
    time = tointeger(dsets(j,5))
    if (isStrSubset(dsets(j,2), "cpl.hi")) then
      time_str = sprinti("%05d", time)
    else
      time_str = sprinti("%03d", time)
    end if
    fl(j,:) = systemfunc("ls "+dsets(j,1)+dsets(j,2)+"*"+dsets(j,4)+"*"+time_str+"*.nc")
  end do
  print(fl)

  ;--- loop over tiles ---
  do j = 0, ndsets-1
    do k = 0, ntiles-1
      ;--- open file ---
      nc = addfile(fl(j,k), "r")

      ;--- read variables ---
      var = nc->$dsets(j,3)$(0,:,:)
      if (isfilevar(nc, "orog")) then 
        orog = nc->orog(0,:,:)
      else
        tmp = addfile("orog/orog.tile"+(k+1)+".nc", "r")
        orog = tmp->orog(0,:,:)
        delete(tmp)
      end if
      if (isfilevaratt(nc, dsets(j,3), "long_name")) then
        name = var@long_name
      else
        name = dsets(j,3)
      end if

      ;--- query dimension of the variable ---
      dim_names = getvardims(var)

      ;--- read coordinate data ---
      if (isfilevar(nc, dim_names(1))) then
        lon2d = (/ nc->$dim_names(1)$ /)
        lat2d = (/ nc->$dim_names(0)$ /)
      else
        lon2d = (/ nc->atmImp_lon(0,:,:) /)
        lat2d = (/ nc->atmImp_lat(0,:,:) /)
      end if

      ;--- create data array ---
      if (j .eq. 0 .and. k .eq. 0) then
        dims = dimsizes(lon2d)
        nlat = dims(0)
        nlon = dims(1)
        data = new((/ ndsets, ntiles, nlat, nlon /), "double", var@_FillValue)
        lsmk = new((/ ndsets, ntiles, nlat, nlon /), "double", 999)
      end if

      ;--- fill common variable ---
      data(j,k,:,:) = (/ var /)
      lsmk(j,k,:,:) = (/ where(orog .gt. 0.0, 1.0, 0.0) /)

      ;--- special care for following tiles ---
      if any((k .eq. (/ 3, 4 /))) then
        data(j,k,:,:) = (/ transpose(data(j,k,:,:)) /)
        data(j,k,:,:) = data(j,k,::-1,:)
        lsmk(j,k,:,:) = (/ transpose(lsmk(j,k,:,:)) /)
        lsmk(j,k,:,:) = lsmk(j,k,::-1,:)
      end if
      if any((k .eq. (/ 2 /))) then
        data(j,k,:,:) = (/ transpose(data(j,k,:,:)) /)
        data(j,k,:,:) = data(j,k,:,::-1)
        lsmk(j,k,:,:) = (/ transpose(lsmk(j,k,:,:)) /)
        lsmk(j,k,:,:) = lsmk(j,k,:,::-1)
      end if

      ;--- mask out large values ---
      data(j,k,:,:) = where(data(j,k,:,:) .gt. 1.0e20, 0.0, data(j,k,:,:))

      ;--- mask out land ---
      if (mask_land) then
        data(j,k,:,:) = mask(data(j,k,:,:), lsmk(j,k,:,:) .gt. 0.0, False)
      end if

      ;--- change sign if it is required ---
      if (any(dsets(j,3) .eq. (/ "Med_aoflux_atm_Faox_lat", "Med_aoflux_atm_Faox_sen" /))) then
        data(j,k,:,:) = -1.0*data(j,k,:,:)
      end if

      print("tile "+k+" "+dsets(j,3)+" min = "+min(data(j,k,:,:))+" max = "+max(data(j,k,:,:)))

      ;--- delete temporary variables ---
      delete([/ var, dim_names /])
    end do
  end do

  ;--- calculate diff ---
  if (diff) then
    data(1,:,:,:) = data(1,:,:,:)-data(0,:,:,:)
    do k = 0, ntiles-1
      print("diff tile "+k+" "+dsets(1,3)+" min = "+min(data(1,k,:,:))+" max = "+max(data(1,k,:,:))+" avg = "+avg(data(1,k,:,:)))
    end do
  end if

  ;--- create plot ---
  wks_type = "pdf"
  wks_type@wkOrientation = "landscape"
  if (diff) then
    wks = gsn_open_wks(wks_type, "plot_cubed_"+dsets(0,3)+"_diff")
  else
    wks = gsn_open_wks(wks_type, "plot_cubed_"+dsets(did,3))
  end if

  res1 = True
  res1@gsnDraw = False
  res1@gsnFrame = False
  res1@gsnAddCyclic = False
  res1@cnFillOn = True
  res1@cnLinesOn = False
  res1@cnLineLabelsOn = False
  res1@gsnLeftString = ""
  res1@gsnRightString = ""
  size = 0.24
  res1@vpWidthF = size
  res1@vpHeightF = size
  res1@sfXArray := lon2d
  res1@sfYArray := lat2d
  res1@cnInfoLabelOn = False
  res1@lbLabelBarOn = False
  res1@cnFillMode = "RasterFill"
  res1@cnLevelSelectionMode = "ExplicitLevels"

  res2 = res1
  res2@cnFillOn = False
  res2@cnLinesOn = True
  res2@cnLineLabelsOn = False
  res2@cnLevels = (/ 0.999999 /)
  res2@cnInfoLabelOn = False
  res2@cnLineThicknessF = 0.5

  if (diff) then
    res1@cnFillPalette = "NCV_blu_red"
    if (any(dsets(did,3) .eq. (/ "tmpsfc", "atmImp_Sa_t2m"/))) then
      res1@cnLevels = fspan(-5, 5, 21)
    end if
    if (dsets(did,3) .eq. "pressfc") then
      res1@cnLevels = fspan(-20, 20, 41)
    end if
    if (any(dsets(did,3) .eq. (/ "Med_aoflux_atm_Faox_lat", "atmImp_Faxa_lat" /))) then
      res1@cnLevels = fspan(-100, 100, 21)
    end if
    if (any(dsets(did,3) .eq. (/ "Med_aoflux_atm_Faox_sen", "atmImp_Faxa_sen" /))) then
      res1@cnLevels = fspan(-100, 100, 21)
    end if
  else
    res1@cnFillPalette = "MPL_viridis"
    if (any(dsets(did,3) .eq. (/ "tmpsfc", "atmImp_Sa_t2m"/))) then
      data = data-273.15 ; K -> degC
      res1@cnLevels = fspan(-30, 30, 61)
    end if
    if (dsets(did,3) .eq. "pressfc") then
      data = data*0.01 ; Pa -> hPa or mb
      res1@cnLevels = fspan(1000.0, 1030.0, 31)
    end if
    if (any(dsets(did,3) .eq. (/ "Med_aoflux_atm_Faox_lat", "atmImp_Faxa_lat" /))) then
      res1@cnLevels = fspan(0, 400, 41)
    end if
    if (any(dsets(did,3) .eq. (/ "Med_aoflux_atm_Faox_sen", "atmImp_Faxa_sen" /))) then
      res1@cnLevels = fspan(-50, 150, 21)
    end if
  end if

  txres = True
  txres@gsnFrame = False
  txres@txJust = "CenterCenter"

  ;--- tile order --- 
  ;   2
  ; 4 0 1 3
  ;   5
  xpos = tofloat((/ 1, 2, 1, 3, 0, 1 /))*size+0.05
  ypos = 0.9-tofloat((/ 1, 1, 0, 1, 1, 2 /))*size

  do k = 0, ntiles-1
    ;--- set position ---
    res1@vpXF = xpos(k)
    res1@vpYF = ypos(k) 
    res2@vpXF = xpos(k)
    res2@vpYF = ypos(k) 
    print(k+" "+res1@vpXF+" "+res1@vpYF)

    if (k .eq. 1) then
      res1@lbLabelBarOn = True
      res1@pmLabelBarWidthF = size*1.8
      res1@pmLabelBarHeightF = 0.03
      res1@pmLabelBarParallelPosF = 0.97
    else
      res1@lbLabelBarOn = False
    end if

    res1@gsnStringFont = 22
    res1@gsnStringFontHeightF = 0.013
    res1@gsnLeftString = ""
    if (k .eq. 4) then
      if (time .gt. 1000) then
        time = time/60
      end if
      if (diff) then
        res1@gsnLeftString = dsets(1,0)+"-"+dsets(0,0)+" @ "+dsets(0,4)+" "+time+"h"
      else
        res1@gsnLeftString = dsets(did,0)+" @ "+dsets(did,4)+" "+time+"h"
      end if
    end if

    plot1 = gsn_csm_contour(wks, data(did,k,:,:), res1)
    draw(plot1)
    if (.not. mask_land) then
      plot2 = gsn_csm_contour(wks, lsmk(did,k,:,:), res2)
      draw(plot2)
    end if

    txres@txPosXF = xpos(k)+0.01
    txres@txPosYF = ypos(k)-0.01
    txres@txAngleF = 0
    txres@txFont = 22
    txres@txFontHeightF = 0.015
    text = gsn_create_text(wks, sprinti("%d", k), txres)
    draw(text)
  end do

  frame(wks)
end

plot_cubed_dif.sh

#!/bin/bash
module load ncl

# field mapping
declare -A variables
variables=( \
            ["atmImp_Faxa_lat"]="Med_aoflux_atm_Faox_lat" \
            ["atmImp_Faxa_sen"]="Med_aoflux_atm_Faox_sen" \
            ["atmImp_Faxa_taux"]="Med_aoflux_atm_Faox_taux" \
            ["atmImp_Faxa_tauy"]="Med_aoflux_atm_Faox_tauy" \
          )

# loop over fields
for key in "${!variables[@]}"
do
  echo "$key - ${variables[$key]}"
  ncl plot_cubed.ncl "vname1=\"$key\"" "vname2=\"${variables[$key]}\""
done

plot_cubed_ref.sh

(base) turuncu@cheyenne1/glade/work/turuncu/NOAHMP/ufs-weather-model (feature/lnd_noahmp) $ cat ../../XGRID/proc/plot_cubed_ref.sh 
#!/bin/bash
module load ncl
lst=`ncdump -c ref/ufs.cpld.cpl.hi.2016-10-04-00000.tile1.nc | grep "double atmImp_" | awk '{print $2}' | awk -F\( '{print $1}'`
for i in $lst
do
  echo $i
  ncl plot_cubed.ncl "vname1=\"$i\"" "vname2=\"$i\""
done

plot_cubed_ts.ncl

load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
begin
  ;--- arguments ---
  ntiles = 6
  dir = "./upd_aoflux_new_flux/"
  vname1 = "atmImp_Faxa_lat"
  vname2 = "Med_aoflux_atm_Faox_lat"

  ;--- list files ---
  fl = systemfunc("ls "+dir+"*.nc")
  fl2d = new((/ ntiles, dimsizes(fl)/ntiles /), "string") 
  do i = 0, ntiles-1
    fl2d(i,:) = fl(i:dimsizes(fl)-1:ntiles)
  end do
  ntime = dimsizes(fl2d(0,:))

  ;--- create data array ---
  data = new((/ ntiles, 2, 3, ntime /), "double")

  ;--- loop over tiles and prepare data ---
  do i = 0, ntiles-1
    ;--- open file ---
    print("[debug] -- processing tile :: "+(i+1))
    nc = addfiles(fl2d(i,:), "r")
    ListSetType(nc, "join")

    ;--- read fields ---
    dum1 = nc[:]->$vname1$(:,0,:,:)
    dum2 = nc[:]->$vname2$(:,0,:,:)

    ;--- get time ---
    time = nc[:]->time
    stime = cd_string(time, "%d %c - %H:%M")

    ;--- create mask ---
    lmask = (/ where(dum1(0,:,:) .gt. 1.0e20, False, True) /)

    ;--- mask data and calculate spatial average ---
    dum1 = mask(dum1, conform(dum1, lmask, (/ 1,2 /)), True)    
    dum2 = mask(dum2, conform(dum2, lmask, (/ 1,2 /)), True)    
    
    ;--- calculate min, max and avg ---
    data(i,0,0,:) = dim_min_n_Wrap(dum1, (/1,2/))
    data(i,0,1,:) = dim_max_n_Wrap(dum1, (/1,2/))
    data(i,0,2,:) = dim_avg_n_Wrap(dum1, (/1,2/))
    data(i,1,0,:) = dim_min_n_Wrap(dum2, (/1,2/))
    data(i,1,1,:) = dim_max_n_Wrap(dum2, (/1,2/))
    data(i,1,2,:) = dim_avg_n_Wrap(dum2, (/1,2/))

    ;--- delete temporary variables ---
    delete([/ nc, dum1, dum2, lmask /])
  end do

  ;--- calculate RMS for each tile ---
  error = dim_rmsd_n_Wrap(data(:,0,2,:), data(:,1,2,:), (/ 1 /))

  ;--- plot data ---
  wks  = gsn_open_wks("newpdf", "plot_cubed_ts")
  gsn_define_colormap(wks, "matlab_jet")

  res = True
  res@gsnDraw = False
  res@gsnFrame = False
  res@vpHeightF = 0.24
  res@vpWidthF = 0.75

  res@tmXBMode = "Explicit"
  res@tmXBValues = ispan(0, ntime-1, 12)
  res@tmXBLabels = stime(::12)

  res@tmXTOn = False
  res@tmXTLabelsOn = False
  res@tmYROn = False
  res@tmYRLabelsOn = False
  res@tmYLMinorOn = False
  res@tmXBLabelAngleF = 90.0
  res@tmXBLabelFontHeightF = 0.01
  res@tmYLLabelFontHeightF = 0.01
  colors = (/ "black", "red", "blue", "green", "gold3", "darkorchid" /)
  res@xyLineColors = colors
  zeros = new((/ntiles/), "integer")
  zeros = 0
  res@xyDashPatterns = zeros
  res@xyLineThicknessF = 2.0
  ;res@gsnYRefLine = 0.0 
  ;res@gsnYRefLineDashPatterns = 16
  res@tiXAxisFontHeightF = 0.012
  res@tiYAxisFontHeightF = 0.012
  res@tiYAxisString = "Heat Flux [W/m2]"
  res@tiXAxisString = "Time"
  res@trXMaxF = max(ispan(0, ntime-1, 12)) 
  res@trYMaxF = max(data(:,0,2,:))*1.05

  res@vpXF = 0.1
  res@vpYF = 0.9

  plot1 = gsn_csm_xy(wks, ispan(0, ntime-1, 1), data(:,0,2,:), res)
  res@xyDashPatterns = zeros+1
  plot2 = gsn_csm_xy(wks, ispan(0, ntime-1, 1), data(:,1,2,:), res)
  overlay(plot1, plot2)
  draw(plot1)
  
  print(data(0,0,2,:)+" "+data(0,1,2,:)+" "+fl2d(0,:))

  ;--- add legend ---
  res_text = True
  res_text@txFontColor = "black"
  res_text@txFontHeightF = 0.01
  res_text@txJust = "CenterLeft"

  res_line = True
  res_line@gsLineThicknessF  = 2.0

  xp = 0.02
  yp = 0.95
  do i = 0, ntiles-1
    xp = xp+0.11

    res_line@gsLineColor = colors(i) 
    gsn_text_ndc(wks, "Tile "+sprinti("%d", i+1)+" ("+sprintf("%5.2f",error(i))+")", xp-0.006, yp-0.035, res_text)
    res_line@gsLineDashPattern = 0
    gsn_polyline_ndc(wks, (/ xp-0.006, xp+0.06 /), (/ yp-0.02, yp-0.02 /), res_line)
    if (i .eq. 0) then
      gsn_text_ndc(wks, "FV3", xp+0.065, yp-0.02, res_text)
      gsn_text_ndc(wks, "MED", xp+0.065, yp, res_text)
      res_line@gsLineDashPattern = 1
      gsn_polyline_ndc(wks, (/ xp-0.006, xp+0.06 /), (/ yp, yp /), res_line)
    end if
  end do

  frame(wks)
end

proc_med_his.ncl

; This function converts 1d array to its tiled representation
; adapted from code originally written by 
; Denise Worthen, IMSG at NOAA/NWS/NCEP/EMC
undef("totile")
function totile(var1d, ntile, res)
begin
  ;--- create new array ---
  var3d = new((/ntile,res,res/), typeof(var1d), default_fillvalue(typeof(var1d)))

  ;--- loop over tiles and assgin data ---
  do i = 0, ntile-1
    istr = (i*res*res)
    iend = istr+res*res-1
    var3d(i,:,:) = onedtond(var1d(istr:iend), (/res,res/))
  end do

  return var3d
end

; This procedure creates mediator history file for FV3
; there will be a different file for each tile
undef("create_tiled_output")
procedure create_tiled_output(thefile, ntiles, atm_only)
begin
  ;--- create file names ---
  ofiles = new((/ ntiles /), "string")
  do i = 0, ntiles-1
    ofiles(i) = str_sub_str(systemfunc("basename "+thefile), ".nc", ".tile"+(i+1)+".nc")
    system("rm -f "+ofiles(i))
  end do

  ;--- open file/s ---
  nc = addfile(thefile, "r")

  ;--- loop over tile files and create them ---
  do i = 0, ntiles-1
    ;--- open output file ---
    print(""+ofiles(i))
    fout = addfile(ofiles(i), "c")

    ;--- enable define mode ---
    setfileoption(fout, "DefineMode", True)

    ;--- query dimensions ---
    dimNames = getfiledims(nc)
    dimSizes = getfiledimsizes(nc)
    ndims = dimsizes(dimNames)

    ;--- keep only FV3 related dimensions ---
    if (atm_only) then
      indx = str_match_ind(dimNames, "atm")
      dimNames_atm = dimNames(indx)
      dimSizes_atm = dimSizes(indx)
      delete([/ dimNames, dimSizes /])
      dimNames = array_append_record(dimNames_atm, "time", 0)
      dimSizes = array_append_record(dimSizes_atm, 1, 0)
      ndims = dimsizes(dimNames)
      delete([/ dimNames_atm, dimSizes_atm, indx /])
    end if
  
    ;--- fix dimensions for one-dimensional variables ---
    do j = 0, ndims-1
      if (isStrSubset(dimNames(j), "atm")) then
        if (isStrSubset(dimNames(j), "_nx")) then
          res = tointeger(sqrt(dimSizes(j)/ntiles))
        end if
        dimSizes(j) = res
      end if
    end do

    dimUnlim = new((/ ndims /), "logical")
    dimUnlim = False

    ;--- define dimensions ---
    filedimdef(fout, dimNames, dimSizes, dimUnlim)

    ;--- query list of variables ---
    varNames = getfilevarnames(nc)

    ;--- keep only FV3 related variables ---
    if (atm_only) then
      indx = str_match_ind(varNames, "atm")
      varNames_atm = varNames(indx)
      delete(varNames)
      varNames = array_append_record(varNames_atm, "time", 0)
      delete([/ varNames_atm, indx /])
    end if    

    do j = 0, dimsizes(varNames)-1
      ;--- query variable dimensions and type ---
      varDims = getfilevardims(nc, varNames(j))
      varType = getfilevartypes(nc, varNames(j))

      ;--- create variable ---
      filevardef(fout, varNames(j), varType, varDims)

      ;--- copy attributes and fill variables --
      filevarattdef(fout, varNames(j), nc->$varNames(j)$)
      if (isStrSubset(varNames(j), "atm")) then
        var = totile(nc->$varNames(j)$(0,0,:), ntiles, res)
        fout->$varNames(j)$ = (/ var(i,:,:) /)
      else
        var = nc->$varNames(j)$
        fout->$varNames(j)$ = (/ var /)
      end if

      ;--- delete temporary variables ---
      delete([/ varDims, varType, var /])
    end do

    ;--- delete temporary variables ---
    delete([/ fout, dimNames, dimSizes, varNames /])
  end do
end

begin
  ;--- arguments ---
  dir = "../cpld_control_upd_aoflux/"
  ntiles = 6

  ;--- list mediator history files ---
  ;mhf = systemfunc("ls "+dir+"ufs.cpld.cpl.hi.2016-10-03-03600*.nc | grep -v fixed | grep -v tile")
  mhf = systemfunc("ls "+dir+"ufs.cpld.cpl.hi.2016-10-04-00000*.nc | grep -v fixed | grep -v tile")

  ;--- process files and convert 1d FV3 data to 2d ---
  do i = 0, dimsizes(mhf)-1
    create_tiled_output(mhf(i), ntiles, True)
  end do
end