begin
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; This ncl reads in pre downloaded MERRA files and generates both the 2-dimensional and 3-dimensional IVT files   ;;;
;;; that were used in ARTMIP Tier 1. 										    ;;;
;;;														    ;;;
;;; Author: Brian Kawzenuk (Center for Western Weather and Water Extremes; University of California, San Diego)     ;;;
;;;														    ;;;
;;; For this script to run several filenames (MERRAfile (line 34), outfile_2D (line 124), and outfile_3D (line 136) ;;;
;;; must be updated to reflect the local computing environment. Three date variables must also be passed to the     ;;;
;;; script (YYYY representing year, MM represnting month, and DD representing day). The accompanied shell script    ;;;
;;; loops through these dates and passes the variables to this script.						    ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Add leading 0 to month and day variables if necessary ;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
yyyy = tostring(YYYY)
if(MM.lt.10) then
 mm = "0"+tostring(MM)
else
 mm = tostring(MM)
end if

if(DD.lt.10) then
 dd = "0"+tostring(DD)
else
 dd = tostring(DD)
end if

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Call in variables from MERRA data files ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
MERRAfile="/path/to/MERRA/data/file/MERRA2_400.inst3_3d_asm_Np."+yyyy+mm+dd+".SUB.nc4"		; MERRA file name
if(isfilepresent(MERRAfile))
 
in = addfile(Mfile,"r")

p = doubletofloat(in->lev*100)		; Pressure levels
QV = in->QV(:,:,:,:)			; Sepcific humidity
U = in->U(:,:,:,:)			; U-component of wind
V = in->V(:,:,:,:)			; V-component of wind
PS = in->PS(:,:,:)			; Surface pressure



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Generate arrays to be used in the IVT and IWV calculations ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    uflux     = U(:,0:dimsizes(U(0,:,0,0))-2,:,:)
    uflux     = uflux*0
    vflux     = uflux
    QFL	      = uflux
    IVT       = U(:,0,:,:)
    IVT       = IVT*0    
    uIVT      = IVT
    vIVT      = IVT
    IWV       = IVT   

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Loop through each pressure level to calculate IVT and IWV ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
m=-1
do while m.lt.(dimsizes(p)-2)
   m=m+1
   layer_diff     = p(m)-p(m+1)					; Depth of layer in pressure (dp)
   ql             = (QV(:,m,:,:)+QV(:,m+1,:,:))/2		; Average specific humidity in layer (q)
   ul             = (U(:,m,:,:)+U(:,m+1,:,:))/2			; Average u-component of wind in layer (u)
   vl             = (V(:,m,:,:)+V(:,m+1,:,:))/2			; Average v-component of wind in layer (v)
   qfl            = (ql/9.8)*layer_diff				; Divide average q by gravity and multiply by dp
   uflux(:,m,:,:) = doubletofloat(qfl*ul)			; Mulitple value in previous line by average u
   vflux(:,m,:,:) = doubletofloat(qfl*vl)			; Same as previous line but for v
   QFL(:,m,:,:)   = (/qfl/)
end do

uIVT            = dim_sum_n(uflux,1)				; Sum array along the pressure levels to calculate u-component of IVT
vIVT            = dim_sum_n(vflux,1)				; Sum array along the pressure levels to calculate v-component of IVT
IVT             = sqrt((uIVT*uIVT)+(vIVT*vIVT))			; Calculate the magnitude of IVT using the u and v components of IVT
IWV             = dim_sum_n(QFL,1)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Add metadata to IWV and IVT variables ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
IWV@long_name = "Integrated Water Vapor"
IWV@standard_name = "Integrated Water Vapor"
IWV@short_name = "IWV"
IWV@units = "mm"
IWV&lat@units = "degrees_north"
IWV&lon@units = "degrees_east"

uIVT@long_name = "U-component of IVT"
uIVT@standard_name = "U-component of IVT"
uIVT@short_name = "uIVT"
uIVT@units = "kg/m/s"
uIVT&lat@units = "degrees_north"
uIVT&lon@units = "degrees_east"

vIVT@long_name = "V-component of IVT"
vIVT@standard_name = "V-component of IVT"
vIVT@short_name = "vIVT"
vIVT@units = "kg/m/s"
vIVT&lat@units = "degrees_north"
vIVT&lon@units = "degrees_east"

IVT@long_name = "Integrated Water Vapor Transport"
IVT@standard_name = "Integrated Water Vapor Transport"
IVT@short_name = "IVT"
IVT@units = "kg/m/s"
IVT&lat@units = "degrees_north"
IVT&lon@units = "degrees_east"

lev2 = U&lev(0:21)
uflux!1 = "Lev"
vflux!1 = "Lev"
uflux&Lev = lev2
vflux&Lev = lev2

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Loop through each timestep to write netcdf files ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
hours = (/"00","03","06","09","12","15","18","21"/)
do time=0, 7
hh = hours(time)
outfile_2D = "/full/path/to/2D/file/ARTMIP_MERRA_2DIVT_"+yyyy+mm+dd+"_"+hh+".nc"		; 2-dimensional output filename
if (isfilepresent(outfile_2D))
 system("/bin/rm "+outfile_2D)
end if
fout1 = addfile(outfile_2D,"c")
fout1 ->PS = PS(time,:,:)
fout1 ->IVT = IVT(time,:,:)
fout1 ->uIVT = uIVT(time,:,:)
fout1 ->vIVT = vIVT(time,:,:)
fout1 ->IWV = IWV(time,:,:)


outfile_3D = "/full/path/to/3D/file/ARTMIP_MERRA_3DIVT_"+yyyy+mm+dd+"_"+hh+".nc"		; 3-dimensional output filename
if (isfilepresent(outfile_3D))
 system("/bin/rm "+outfile_3D)
end if
fout2 = addfile(outfile_3D,"c")
fout2 ->QV = QV(time,:,:,:)
fout2 ->u = U(time,:,:,:)
fout2 ->v = V(time,:,:,:)
fout2 ->uIVT_levels = uflux(time,:,:,:)
fout2 ->vIVT_levels = vflux(time,:,:,:)
end do

end if
end