Material for Noah-MP Workshop & Tutorial 2023

README.txt

# install docker
https://www.docker.com

# create container and prepare for development
docker run -it ubuntu:latest

apt-get -qq update
apt-get -qq install tar unzip file gringo
apt-get -qq install build-essential binutils-dev gfortran
apt-get -qq install python3-dev python3-pip python3-boto3 python3-yaml
apt-get -qq install wget awscli ca-certificates gh libxml2-dev
pip3 install --target ~/.local botocore

# install dependenciies with spack
mkdir /test
cd /test
git clone https://github.com/spack/spack.git
cd spack
git checkout b2c3973
cd -
. spack/share/spack/setup-env.sh
spack compiler find
spack external find

# create ~/.spack/config.yaml
config:
  url_fetch_method: curl
  connect_timeout: 60

# create spack.yaml (under /test directory) with following content
spack:
  concretizer:
    targets:
      granularity: generic
      host_compatible: false
    unify: when_possible
  specs:
  - esmf@8.5.0b10+external-parallelio %gcc@11.3.0 target=x86_64
  packages:
    all:
      target: ['x86_64']
      compiler: [gcc@11.3.0]
  view: /test/.spack-ci/view
  config:
    source_cache: /test/.spack-ci/source_cache
    misc_cache: /test/.spack-ci/misc_cache
    test_cache: /test/.spack-ci/test_cache
    install_tree:
      root: /test/.spack-ci/opt
    install_missing_compilers: true

# install dependencies (these commands must be run in same folder with spack.yaml)
spack --color always -e . concretize -f
spack --color always -e . install -j3 --deprecated --no-checksum

# build and install data component
export PATH=/test/.spack-ci/view/bin:$PATH
export ESMFMKFILE=/test/.spack-ci/view/lib/esmf.mk
export FC=gfortran
cd /test
git clone -b hotfix/std_build https://github.com/uturuncoglu/CDEPS.git cdeps
cd cdeps
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=../install \
  -DPIO_C_LIBRARY=/test/.spack-ci/view/lib \
  -DPIO_C_INCLUDE_DIR=/test/.spack-ci/view/include \
  -DPIO_Fortran_LIBRARY=/test/.spack-ci/view/lib \
  -DPIO_Fortran_INCLUDE_DIR=/test/.spack-ci/view/include \
  -DCMAKE_Fortran_FLAGS="-ffree-line-length-none -fallow-argument-mismatch -fallow-invalid-boz" \
  -DDISABLE_FoX=ON ../
make
make install

# build and install noahmp
export NetCDF_ROOT=/test/.spack-ci/view
export FC=gfortran
cd /test
git clone -b feature/nuopc_cap https://github.com/esmf-org/noahmp.git
cd noahmp
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=../install ../
make
make install

# create /test/esmxBuild.yaml with following content
components:
  datm:
    cmake_config: /test/cdeps/install/lib/cmake/datm-esmx.cmake
    fort_module: cdeps_datm_comp
  noahmp:
    cmake_config: /test/noahmp/install/lib/cmake/noahmp-esmx.cmake
    fort_module: lnd_comp_nuopc

# create executable using NUOPC provided generic driver - ESMX
export ESMF_ESMXDIR=/test/.spack-ci/view/include/ESMX
cd /test
cmake -H$ESMF_ESMXDIR -Bbuild
cd build
make

# create run directory and download input
export PYTHONPATH=~/.local:$PYTHONPATH
mkdir /test/run
cd /test/run
python3 /test/noahmp/.github/workflows/scripts/get_input.py --ifile /test/noahmp/.github/workflows/tests/test_datm_lnd.yaml 
python3 /test/noahmp/.github/workflows/scripts/gen_config.py --ifile /test/noahmp/.github/workflows/tests/test_datm_lnd.yaml

# run application (you might need to add extra options to mpirun like --mca btl_tcp_if_include eth0)
mpirun --oversubscribe -np 6 /test/build/esmx

# you could copy files outside of the Docker container using docker cp command!

plot_cubed.ncl

load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"

procedure add_box(wks,plot,ll[2],width,cname,oval)
local xbox, ybox, gsres, dumstr
begin
  xbox = (/ll(0),ll(0)+width,ll(0)+width,ll(0),ll(0)/)
  ybox = (/ll(1),ll(1),ll(1)+width,ll(1)+width,ll(1)/)

  gsres                = True
  gsres@gsFillColor    = cname
  gsres@gsFillOpacityF = oval
  dumstr               = unique_string("gon")
;
; Adding it as an attribute is a sneaky way to
; make sure it "lives" outside this procedure.
;
  plot@$dumstr$ = gsn_add_polygon(wks,plot,xbox,ybox,gsres)
end

begin
  ;--- parameters ---
  ntiles = 6

  ;--- arguments (suffix: (1) FV3 and (2) Noah-MP) ---
  ;vname1 = "lhtfl"
  ;vname2 = "evap"
  ;vname1 = "tmpsfc"
  ;vname2 = "tskin"
  vname1 = "shtfl"
  vname2 = "hflx" 
  
  tid1 = 24 ;1
  tid2 = 23 ;0
  diff = True 
  dataset = 1 ; 0 for FV3 and 1 for Noah-MP (only for diff = False)
  percentage_diff = False

  ;--- loop over tiles ---
  do i = 0, ntiles-1
    ;--- open file ---
    lsA = systemfunc ("ls sfcf*.tile"+sprinti("%d", (i+1))+".nc")
    ncA = addfiles(lsA, "r")
    ListSetType(ncA, "cat")
    lsB = systemfunc ("ls ufs.cpld.lnd.out.*.tile"+sprinti("%d", (i+1))+".nc")
    ncB = addfiles(lsB, "r")
    ListSetType(ncB, "cat")
    ncC = addfile("INPUT/oro_data.tile"+sprinti("%d", (i+1))+".nc", "r")

    ;--- read fields ---
    dum1 = ncA[:]->$vname1$(tid1,:,:)
    dum2 = ncB[:]->$vname2$(tid2,:,:)
    land_frac = ncC->land_frac

    ;--- read coordinate data ---
    lon2d = (/ ncA[0]->grid_xt /)
    lat2d = (/ ncA[0]->grid_yt /)

    ;--- string for label bar ---
    if (diff) then
      lstr = "Diff: "+vname2+"-"+vname1+" ["+dum2@units+"]"
    else
      if (dataset .eq. 0) then
        lstr = vname1+" ["+dum1@units+"]"
      else
        lstr = vname2+" ["+dum2@units+"]"
      end if
    end if

    time = ncA[:]->time(tid1)
    time@calendar = "standard"
    time_str = ut_string(time, "%D-%N-%Y_%H:%M")

    ;--- create data array ---
    if (i .eq. 0) then
      dims = dimsizes(lon2d)
      nlat = dims(0)
      nlon = dims(1)
      data = new((/ 2, ntiles, nlat, nlon /), "double", dum1@_FillValue)
      lmsk = new((/ ntiles, nlat, nlon /), "integer")
      lon  = new((/ ntiles, nlat, nlon /), typeof(lon2d))
      lat  = new((/ ntiles, nlat, nlon /), typeof(lat2d))
    end if

    ;--- store coordinate arrays ---
    lon(i,:,:) = (/ lon2d /)
    lat(i,:,:) = (/ lat2d /)
    if (i .eq. 0 .or. i .eq. 5) then
      lon(i,:,:) = where(lon(i,:,:) .ge. 180. .or. lon(i,:,:) .eq.360., lon(i,:,:)-360, lon(i,:,:) )
    end if

    ;--- fill common variable ---
    data(0,i,:,:) = (/ dum1 /)
    data(1,i,:,:) = (/ dum2 /)
    lmsk(i,:,:) = (/ where(land_frac .gt. 0, 1, 0) /)
    lmsk(i,:,:) = (/ mask(lmsk(i,:,:), lmsk(i,:,:) .gt. 0, True) /)

    ;--- special care for following tiles ---
    if any((i .eq. (/ 3, 4 /))) then
      data(0,i,:,:) = (/ transpose(data(0,i,:,:)) /)
      data(0,i,:,:) = data(0,i,::-1,:)
      data(1,i,:,:) = (/ transpose(data(1,i,:,:)) /)
      data(1,i,:,:) = data(1,i,::-1,:)
      lmsk(i,:,:) = (/ transpose(lmsk(i,:,:)) /)
      lmsk(i,:,:) = lmsk(i,::-1,:)
      lon(i,:,:) = (/ transpose(lon(i,:,:)) /)
      lon(i,:,:) = lon(i,::-1,:)
      lat(i,:,:) = (/ transpose(lat(i,:,:)) /)
      lat(i,:,:) = lat(i,::-1,:)
    end if
    if any((i .eq. (/ 2 /))) then
      data(0,i,:,:) = (/ transpose(data(0,i,:,:)) /)
      data(0,i,:,:) = data(0,i,:,::-1)
      data(1,i,:,:) = (/ transpose(data(1,i,:,:)) /)
      data(1,i,:,:) = data(1,i,:,::-1)
      lmsk(i,:,:) = (/ transpose(lmsk(i,:,:)) /)
      lmsk(i,:,:) = lmsk(i,:,::-1)
      lon(i,:,:) = (/ transpose(lon(i,:,:)) /)
      lon(i,:,:) = lon(i,:,::-1)
      lat(i,:,:) = (/ transpose(lat(i,:,:)) /)
      lat(i,:,:) = lat(i,:,::-1)
    end if

    ;--- maskout data ---
    data(0,i,:,:) = mask(data(0,i,:,:), lmsk(i,:,:) .eq. 1, True)
    data(0,i,:,:) = mask(data(0,i,:,:), .not. ismissing(lmsk(i,:,:)), True)

    ;--- delete temporary variables ---
    delete([/ dum1, dum2, lat2d, lon2d, land_frac /])
  end do

  ;--- calculate diff ---
  if (diff) then
    did = 1
    if (percentage_diff) then
      data(0,:,:,:) = mask(data(0,:,:,:), data(0,:,:,:) .eq. 0.0, False) 
      data(1,:,:,:) = (data(1,:,:,:)-data(0,:,:,:))/data(0,:,:,:)
    else
      data(1,:,:,:) = data(1,:,:,:)-data(0,:,:,:)
    end if
  end if

  ;--- create plot ---
  wks_type = "pdf"
  wks_type@wkOrientation = "landscape"
  if (diff) then
    wks = gsn_open_wks(wks_type, "plot_cubed_diff_"+vname1+"_"+vname2+"_"+time_str)
  else
    wks = gsn_open_wks(wks_type, "plot_cubed_"+time_str)
  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.23
  res1@vpWidthF = size
  res1@vpHeightF = size
  res1@cnInfoLabelOn = False
  res1@lbLabelBarOn = False
  res1@cnFillMode = "RasterFill"
  res1@cnLevelSelectionMode = "ExplicitLevels"
  res1@cnMissingValFillColor = "white"
  res1@tmXBOn = False
  res1@tmXTOn = False
  res1@tmYLOn = False
  res1@tmYROn = False

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

  if (diff) then
    res1@cnFillPalette = "NCV_blu_red"

    if (percentage_diff) then
      res1@cnLevels = fspan(-0.25, 0.25, 51)
    else
      if (any(vname2 .eq. (/ "xcoszin", "xlatin" /))) then
        res1@cnLevels = fspan(-0.1, 0.1, 21)
      end if
      if (any(vname2 .eq. (/ "zorl", "sigmaf" /))) then
        res1@cnLevels = fspan(-0.01, 0.01, 21)
      end if
      if (any(vname2 .eq. (/ "rainn_mp", "rainc_mp", "snow_mp", "ice_mp", "graupel_mp", "zvfun" /))) then
        res1@cnLevels = fspan(-0.001, 0.001, 21)
      end if
      if (any(vname2 .eq. (/ "snet", "dlwflx", "dswsfc", "ps", "weasd" /))) then
        res1@cnLevels = fspan(-10, 10, 21)
      end if
      if (any(vname2 .eq. (/ "fm1", "fh1", "fm101", "fh21" /))) then
        res1@cnLevels = fspan(-5, 5, 21)
      end if
      if (any(vname2 .eq. (/ "fh1", "ustar1", "stress1", "fm101", "sfalb", "prslki", "prslk1", "prsik1", "wind" /))) then
        res1@cnLevels = fspan(-1, 1, 21)
      end if
      if (any(vname2 .eq. (/ "cm", "ch", "q1", "zf", "sfalb", "tsurf", "tskin", "t1", "tg3" /))) then
        res1@cnLevels = fspan(-1, 1, 21)
      end if
      if (any(vname2 .eq. (/ "u1", "v1", "rb1" /))) then
        res1@cnLevels = fspan(-5, 5, 21)
      end if
      if (any(vname2 .eq. (/ "garea" /))) then
        res1@cnLevels = fspan(-500, 500, 21)
      end if
      if (any(vname2 .eq. (/ "lhtfl", "hflx", "evap", "snwdph" /))) then
        res1@cnLevels = fspan(-50, 50, 21)
      end if
    end if
  else
    res1@cnFillPalette = "NCV_banded"
    if (any(vname1 .eq. (/ "soilt1", "tg3"/))) then
      data = data-273.15 ; K -> degC
      res1@cnLevels = fspan(-30, 30, 61)
    end if
    if (any(vname1 .eq. (/ "evap"/)) .or. any(vname2 .eq. (/ "evap"/))) then
      res1@cnLevels = fspan(-50, 200, 25)
    end if
    if (any(vname1 .eq. (/ "hflx"/)) .or. any(vname2 .eq. (/ "hflx"/))) then
      res1@cnLevels = fspan(-100, 300, 41)
    end if
  end if

  ;--- text resource ---
  txres = True
  txres@gsnFrame = False
  txres@txJust = "CenterCenter"

  ;--- box resource ---
  gnres = True
  gnres@gsFillColor = -1
  gnres@gsEdgesOn = True
  gnres@gsEdgeColor = "grey"

  ;--- 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)

    ;--- add labelbar ---
    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

    ;--- add plot label ---
    res1@gsnStringFont = 22
    res1@gsnStringFontHeightF = 0.013
    res1@gsnLeftString = ""
    if (k .eq. 4) then
      if (diff) then
        res1@gsnLeftString = vname2+"-"+vname1
      else
        if (dataset .eq. 0) then
          res1@gsnLeftString = vname1
        else
          res1@gsnLeftString = vname2
        end if
      end if
    end if

    ;--- plot ---
    if (diff) then
      plot1 = gsn_csm_contour(wks, data(1,k,:,:), res1)
    else
      plot1 = gsn_csm_contour(wks, data(dataset,k,:,:), res1)
    end if
    draw(plot1)

    ;--- add box around each plot ---
    x_box = (/ 0.0, size, size, 0.0, 0.0 /)+xpos(k)
    y_box = (/ 0.0, 0.0, size, size, 0.0 /)+ypos(k)-size
    gsn_polygon_ndc(wks, x_box, y_box, gnres)

    ;--- add tile number ---
    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+1), txres)
    draw(text)

    ;--- add labelbar string ---
    if (k .eq. 3) then
      txres@txPosXF = xpos(k)+0.01
      txres@txPosYF = ypos(k)-size-0.04
      txres@txFontHeightF = 0.011
      text = gsn_create_text(wks, lstr, txres)
      draw(text)
    end if
  end do

 ;--- add date string ---
 txres@txPosXF = xpos(3)+size
 txres@txPosYF = ypos(3)+0.01
 txres@txJust = "centerright"
 text = gsn_create_text(wks, time_str, txres)
 draw(text)

 ;--- create plot ---
 frame(wks)
end

plot_cubed_ts.ncl

load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
begin
  ;--- arguments ---
  ntiles = 6
  vname1 = "lhtfl"  ; latent heat
  vname2 = "evap"
  ;vname1 = "shtfl"  ; sensible heat total
  ;vname2 = "hflx"
  ;vname1 = "tmpsfc" ; temperature radiative
  ;vname2 = "tskin"

  ;--- coordinates of predefined land grid point for each tile ---
  ;--- different set of grid points can be defined in here --- 
  xc = (/ 79, 46, 78, 78, 10, 54 /)
  yc = (/ 47, 83, 12, 12, 42, 47 /)

  ;--- loop over tiles and prepare data ---
  do i = 0, ntiles-1
    ;--- open file ---
    lsA = systemfunc ("ls sfcf*.tile"+sprinti("%d", (i+1))+".nc")
    ncA = addfiles(lsA, "r")
    ListSetType(ncA, "cat")
    lsB = systemfunc ("ls ufs.cpld.lnd.out.*.tile"+sprinti("%d", (i+1))+".nc")
    ncB = addfiles(lsB, "r")
    ListSetType(ncB, "cat")

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

    ;--- read fields ---
    dum1 = ncA[:]->$vname1$(1:ntime,yc(i),xc(i))
    dum2 = ncB[:]->$vname2$(:ntime-1,yc(i),xc(i))

    ;--- create data array ---
    if (i .eq. 0) then
      data = new((/ 2, ntiles, ntime /), "double")
    end if

    data(0,i,:) = (/ dum1 /)
    data(1,i,:) = (/ dum2 /)

    ;--- delete temporary variables ---
    delete([/ ncA, ncB, dum1, dum2 /])
  end do

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

  ;--- plot data ---
  wks  = gsn_open_wks("newpdf", "plot_cubed_ts_"+vname1+"_"+vname2)
  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 = "Data"
  res@tiXAxisString = "Time"
  res@trXMaxF = max(ispan(0, ntime-1, 12)) 
  minF = min(data)
  if (minF .lt. 0.0) then
    res@trYMinF = min(data)*1.1
  else
    res@trYMinF = min(data)*0.9
  end if
  maxF = max(data)
  if (maxF .lt. 0.0) then
    res@trYMaxF = max(data)*0.9
  else
    res@trYMaxF = max(data)*1.1
  end if

  res@vpXF = 0.1
  res@vpYF = 0.9

  plot1 = gsn_csm_xy(wks, ispan(0, ntime-1, 1), data(0,:,:), res)
  res@xyDashPatterns = zeros+1
  plot2 = gsn_csm_xy(wks, ispan(0, ntime-1, 1), data(1,:,:), res)
  overlay(plot1, plot2)
  draw(plot1)
  
  ;--- 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.126

    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, "Noah-MP", 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