; WA 11/20/09 for KNMI testbed from JH build_forcing_file.ncl
; WA 1/11/10 major revision to separate advective and large-scale forcing
; WA 12/26/12 for ARM case in EDMF comparison, full implementation with
;     specified advection and radiation.  From convert_kpt_to_wrf.ncl.
; WA 4/18/16 for RICO case in WFIP testing.
; WA 10/5/16 for LASSO cases in WFIP testing.
;
; Builds a forcing file for the SCM from LASSO case specs.
; Input files for different cases are handled by linking the source directory
; to local directory lasso_config.  
; This is done in the run script (run_lasso.sh).
;
; NOTES: 

begin

; fixed parameters
P0 = 100000.0
Rd = 287.0
p00 = 1.0e5
t00 = 290.0
l00 = 50.0
g = 9.81

Lv = 2.5e6
Rv = 461.6
Rd = 287.0
e0 = 0.6112
T0 = 273.15
eps = Rd/Rv

nli = 751 ; number of levels in the input driver file
nlo = 76  ; number of levels to put in the forcing file 

fName = "lasso_force.nc"
print("Filling forcing file "+fName)

; open the forcing file (product of this script)
forceFile = addfile(fName,"rw")

; first modify some global attributes
dx = 13000      ; Set to (arbitrary) value used in namelist
forceFile@DX = dx ; CHECK THIS!!!!!!!!!!!!!
forceFile@DY = dx

; open LASSO-provided large-scale forcing file here
inF = addfile("lasso_config/input_ls_forcing.nc","r")

; Times copied from input file
forceFile->Times = inF->Times
forceFile@START_DATE = (/ forceFile->Times(0,:) /)
forceFile@SIMULATION_START_DATE = (/ forceFile->Times(0,:) /)

; Heights
; print(dimsizes(forceFile->Z_FORCE))
; print(dimsizes(inF->Z_LS))
print(dimsizes(inF->Z_LS(:,0:nli-1:10)))
forceFile->Z_FORCE = (/ inF->Z_LS(:,0:nli-1:10) /) + 320.0   ; Make it ASL
forceFile->Z_FORCE_TEND = (/ inF->Z_LS_TEND(:,0:nli-1:10) /)
 
print("Done with z_force")

;----------------------------------------------------
; U_G, V_G, and associated tendencies
;----------------------------------------------------
; WA Using provided U_LS and its TEND for U_G, OK?
; Has the effect of instantaneous forcing as desired????
; WA 10/19/16 try using initial values only for U_G, V_G and 
; using largescale vars for forcing instead.

; forceFile->U_G = (/ inF->U_LS(:,0:nli-1:10) /)
; forceFile->V_G = (/ inF->V_LS(:,0:nli-1:10) /)
; forceFile->U_G_TEND = (/ inF->U_LS_TEND(:,0:nli-1:10) /)
; forceFile->V_G_TEND = (/ inF->V_LS_TEND(:,0:nli-1:10) /)
forceFile->U_G = conform_dims( dimsizes(forceFile->U_G),(/ inF->U_LS(1,0:nli-1:10) /),1 )
forceFile->V_G = conform_dims( dimsizes(forceFile->V_G),(/ inF->V_LS(1,0:nli-1:10) /),1 )
forceFile->U_G_TEND = (/ 0.0 /)
forceFile->V_G_TEND = (/ 0.0 /)
forceFile->U_LARGESCALE = (/ inF->U_LS(:,0:nli-1:10) /)
forceFile->V_LARGESCALE = (/ inF->V_LS(:,0:nli-1:10) /)
forceFile->U_LARGESCALE_TEND = (/ inF->U_LS_TEND(:,0:nli-1:10) /)
forceFile->V_LARGESCALE_TEND = (/ inF->V_LS_TEND(:,0:nli-1:10) /)

print("Done with u_g,v_g")

;----------------------------------------------------
; W_SUBS and associated tendency
;----------------------------------------------------

forceFile->W_SUBS = (/ inF->W_LS(:,0:nli-1:10) /) ; WA CHECK SIGN!!!!!
forceFile->W_SUBS_TEND = (/ inF->W_LS_TEND(:,0:nli-1:10) /)

;----------------------------------------------------
; All advection/relaxation terms and tendencies
; Using "largescale" variables which are simpler.
; Only one scale is in use in LASSO.
; Input TH_ADV is already a tendency (K/s), same for Q
; Ignoring the second-order tendency provided for now

; forceFile->TH_LARGESCALE = (/ inF->TH_ADV(:,0:nli-1:10) /)
; forceFile->QV_LARGESCALE = (/ inF->TH_ADV(:,0:nli-1:10) /)
forceFile->TH_LARGESCALE = (/ 0.0 /)
forceFile->QV_LARGESCALE = (/ 0.0 /)
; U,V already taken care of by putting them in geostrophic terms above

; forceFile->TH_LARGESCALE_TEND = (/ inF->TH_ADV_TEND(:,0:nli-1:10) /)
; forceFile->QV_LARGESCALE_TEND = (/ inF->QV_ADV_TEND(:,0:nli-1:10) /)
forceFile->TH_LARGESCALE_TEND = (/ inF->TH_ADV(:,0:nli-1:10) /)
forceFile->QV_LARGESCALE_TEND = (/ inF->QV_ADV(:,0:nli-1:10) /)

nForceTimes = 31
tau_x = new((/nForceTimes+1,nlo/),"float")
tau_x = 60  ; WA should this be matched to WRF timestep?
; tau_x = 3600
tau_x_tend = new((/nForceTimes,nlo/),"float")
tau_x_tend  = 0    ; Does not change
; print(dimsizes(tau_x))
; print(dimsizes(forceFile->TAU_LARGESCALE))
forceFile->TAU_LARGESCALE = (/ tau_x(0:nForceTimes-1,0:nlo-1)/)
forceFile->TAU_LARGESCALE_TEND = (/ tau_x_tend(0:nForceTimes-1,0:nlo-1)/)
delete(tau_x)
delete(tau_x_tend)

print("Done with advection")

;-------------------------------------------------------------
; Surface flux forcing
;-------------------------------------------------------------

; open LASSO-provided surface forcing file here
sfcF = addfile("lasso_config/input_sfc_forcing.nc","r")
forceFile->HFX_FORCE = (/ sfcF->PRE_SH_FLX /)
forceFile->LH_FORCE = (/ sfcF->PRE_LH_FLX /)
forceFile->TSK_FORCE = (/ sfcF->PRE_TSK /)
forceFile->HFX_FORCE_TEND = (/ sfcF->PRE_SH_FLX_TEND /)
forceFile->LH_FORCE_TEND = (/ sfcF->PRE_LH_FLX_TEND /)
forceFile->TSK_FORCE_TEND = (/ sfcF->PRE_TSK_TEND /)


end