[ncl-talk] How to convert NCL script to a NetCDF file?

[Setareh Rahimi]

Dear all NCL users,
Using R software, I want to plot aerosol optical depth (AOD) from the
WRF-chem model. Once import the model output in the R, the software can not
recognize the lat/lon. So I need to convert my NCL script to a NetCDF file
in an easily readable way for the R. How can I convert my NCL script (which
includes converted time, lat/lon, and calculate AOD) to a new  NetCDF file,
please? The script and the model output have been attached.
 wef_chem.nc.gz
<https://drive.google.com/file/d/1K7lEFyBtr3B88PHF6QWhMRB9yc_ORn3e/view?usp=drive_web>

-- 
S.Rahimi
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20221030/dc3341f4/attachment.htm>
-------------- next part --------------
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
;***************************************************************************************************************************



;;f = addfile ("wrfout_d01_2011-07-01_00:00:12", "r")
  f = addfile ("setareh_chem.nc", "r")
      print(f)                              ; same as "nmcdump -h"

  Times = f->Times                          ; [Time | 41] x [DateStrLen | 19]  (character)
  Time_0 = wrf_times_c( Times, 0 )          ; "hours since" initial time on file   (double)
                                            ;  Time_0 is recognized by cd_calendar
  Time_3 = wrf_times_c( Times, 3 )          ; yyyymmddhh  (type "integer") 
  Time_3 = Time_3/100                       ; remove 'ss' part of Time_3
  Time_3 at units = "yyyymmdd"      
 
  print("********************************************************************************************************************")
  
  wks = gsn_open_wks("pdf" ,"WRF_aod_Daily-Average")                   ; ps,pdf,x11,ncgm,eps  
  gsn_define_colormap(wks,"BlAqGrYeOrReVi200")                         ; select color map
    
  res = True    ; plot mods desired
  res at cnConstFEnableFill =  True
  res at gsnMaximize = True
  res at gsnSpreadColors = True
  res at cnFillOn = True
  res at cnLinesOn = False
  res at cnLineLabelsOn = False
  WRF_map_c(f, res, 0)
  res at tfDoNDCOverlay = True
  res at pmTickMarkDisplayMode = "Always"
 
  
 ;***************************************************************************************************************************
  
  
  ymd = cd_calendar(Time_0, -2)
  print(ymd)    
  YMD = 20110701 
  mt  = ind(YMD.eq.ymd) 
;;print(mt) 
 
  A2           = f->TAUAER2      ; (Time, south_north, west_east)
  A2&Time      = Time_0          ; associate the 'Time' coordinate with the variable using  standard NCL & syntax

  a2           = A2(mt,:,:,:)    ; (Time, south_north, west_east)
  printVarSummary(a2)            ; [Time | 8] x [bottom_top | 35] x [south_north | 199] x [west_east | 199]
  delete(A2)

  B3           = f->TAUAER3      ; (Time, south_north, west_east)
  B3&Time      = Time_0          ; associate the 'Time' coordinate with the variable using  standard NCL & syntax

  b3           = B3(mt,:,:,:)    ; (Time, south_north, west_east)
  printVarSummary(b3)            ; [Time | 8] x [bottom_top | 35] x [south_north | 199] x [west_east | 199]
  delete(B3)

  time_3       = Time_3(mt)      ; (Time | 8)
  print(time_3)

  print("********************************************************************************************************************")
  

  opt_cdv = True                 ; option for calculate_daily_values  (cdv)
  opt_cdv at nval_crit = 4          ; default is 1

  a2Day = calculate_daily_values (a2, "avg", 0, opt_cdv)
  printVarSummary(a2Day)         ; [Time | 1] x [bottom_top | 35] x [south_north | 199] x [west_east | 199]        
  printMinMax (a2Day,1)
  
  b3Day = calculate_daily_values (b3, "avg", 0, opt_cdv)
  printVarSummary(b3Day)         ; [Time | 1] x [bottom_top | 35] x [south_north | 199] x [west_east | 199]
  printMinMax (b3Day,1)
  
  print("")
  print("***********************************************************************************************************")
  print("")
  
  angstrom_exponent = -(log(a2Day)-log(b3Day))/log(400./600.)
  printVarSummary(angstrom_exponent)          ; [1] x [35] x [199] x [199]
  printMinMax((angstrom_exponent)          ; [1] x [35] x [199] x [199]
  
  AOD550_3D = a2Day * ((400.0/550.0)^angstrom_exponent)
 ;AOD550_3D at long_name   =  "AOD550_3D"
  AOD550_3D at description =  "AOD550_3D"
  copy_VarCoords(a2Day,AOD550_3D )
  printVarSummary(AOD550_3D)      ; [Time | 1] x [bottom_top | 35] x [south_north | 199] x [west_east | 199]
  
  AOD550_2D = dim_sum_n_Wrap(AOD550_3D,0)
 ;AOD550_2D at long_name   =  "AOD550_2D"
  AOD550_2D at description =  "AOD550_2D"
  printVarSummary(AOD550_2D)      ; [bottom_top | 35] x [south_north | 199] x [west_east | 199]
  
  if (any(isnan_ieee(AOD550_2D))) then
     value = 9.96921e+36  
     replace_ieeenan (AOD550_2D,value,0)
     AOD550_2D at _FillValue = value
  end if
  
  printVarSummary(AOD550_2D)
  printMinMax (AOD550_2D,1)
  print("")
  print("***********************************************************************************************************")
  print("")

  res at tiMainString = "WRF-CHEM (aod_550) " + time_3(0)
;;res at gsnLeftString = a at description                            ; ?????????????????
  res at gsnLeftString = AOD550_2D at description
  print(res)

;;plot = gsn_csm_contour_map(wks,AOD550_2D( ???,:,:),res)
;;;====================================
;;;  Possible changes

;;;
;;;   AOD550_2D = dim_sum_n_Wrap(AOD550_3D,1)   ; sum all the levels
;;;   ;AOD550_2D at long_name   =  "AOD550_2D"
;;;   AOD550_2D at description =  "AOD550_2D"
;;;   printVarSummary(AOD550_2D)      ; [Time | 1] x [south_north | 199] x [west_east | 199]
 
;;;  if (any(isnan_ieee(AOD550_2D))) then
;;;     value = 9.96921e+36
;;;     replace_ieeenan (AOD550_2D,value,0)
;;;     AOD550_2D at _FillValue = value
;;;  end if
 
;;;  printVarSummary(AOD550_2D)
;;;  printMinMax (AOD550_2D,1)
;;;  print("")
;;; ; ;print("***********************************************************************************************************")
;;;  print("")

;;;  res at tiMainString = "WRF-CHEM (aod_550) " + time_3(0)
;;res at gsnLeftString = a at description                            ; ?????????????????
;;;  res at gsnLeftString = AOD550_2D at description
print(res)

;;plot = gsn_csm_contour_map(wks,AOD550_2D(0,:,:),res)