% generate forcing data from 3DCVA output
%
% advections need to be re-calculated from domain mean fields.
% set x_1=x_2 and y_1=y_2 for one grid forcing

clear all;
close all;
addpath '../';  % need resize_lev_time.m

%% settings 
% set the input and output filenames
ipt_file = [ '../3DCVAout_ruc.nc'];
opt_file = [ '../forcing_ruc_10by9_onefile.nc'];

% set the domain of the forcing data
% you have two input options: lat/lon or x/y
% lat/lon grids are centered at *.25 and *.75
lat_1=34.5;   lat_2=39;
lon_1=-100;  lon_2=-95;
% x_1=5;x_2=5;
% y_1=5;y_2=5;

% function generate_forcing_3DCVA(ipt_file,opt_file,x_1,x_2,y_1,y_2)
disp(opt_file);

% set time parameters
time_interval = 1200;    % output time intervals, in second (default 20min)
base_date = 20110422;   % yyyymmdd

% set constants
cpd=1005.7;
Lv0=2.5e6;
missing_value = -9999;

%% read in data
ncid = netcdf.open(ipt_file, 'nc_nowrite');
varname='time';
varid = netcdf.inqVarID(ncid,varname);
time = netcdf.getVar(ncid,varid,'double');
timeunit = netcdf.getAtt(ncid,varid,'units');
varname='lat';
varid = netcdf.inqVarID(ncid,varname);
lat = netcdf.getVar(ncid,varid,'double');
varname='lon';
varid = netcdf.inqVarID(ncid,varname);
lon = netcdf.getVar(ncid,varid,'double');
varname='lev';
varid = netcdf.inqVarID(ncid,varname);
lev = netcdf.getVar(ncid,varid,'double');
varname='lat_stag';
varid = netcdf.inqVarID(ncid,varname);
lat_stag = netcdf.getVar(ncid,varid,'double');
varname='lon_stag';
varid = netcdf.inqVarID(ncid,varname);
lon_stag = netcdf.getVar(ncid,varid,'double');


initial_time=['seconds ' timeunit(9:end)];
disp(['time starts from: ' initial_time]);
ttt=strsplit(timeunit,':');
ttt=ttt{1};
initial_hour=str2double(ttt(end-1:end));
tsec=time+initial_hour*3600;     % change into second
tsec_new = tsec(1):time_interval:tsec(length(tsec));
bdate = base_date+floor(tsec/86400);

% if input is longitude and latitude, change to grid index
if ~exist('x_1','var')
    x_1=find(lon>=lon_1,1);
    x_2=find(lon<=lon_2,1,'last');
    y_1=find(lat>=lat_1,1);
    y_2=find(lat<=lat_2,1,'last');
end
x_domain=x_1:x_2;
y_domain=y_1:y_2;

lat_new=mean(lat(y_domain));lon_new=mean(lon(x_domain));

lev = lev*100; % hPa to Pa
lev = flip(lev,1);
lev_new = 100*[100:25:1000];     % reshape vertical levels

timelen=length(tsec);
timelen_new=length(tsec_new);
levlen=length(lev);
levlen_new=length(lev_new);
latlen=length(lat);
lonlen=length(lon);
latlen_new=length(y_domain);
lonlen_new=length(x_domain);
clear lat lon;

%% read in variables
varname_input={'avg_T','avg_r','avg_u','avg_v','avg_z','omega_c','Q1','m_Q2', ...
    'avg_E','h_advr','h_advE','h_advT','v_advT','v_advr','v_advE','m_divuv',...
    'm_Dps_dt','rsd_water','m_sdivuvr','evapor','rsd_heat','m_sdivuvh',...
    'radiation','Lprecip',...
    'zs','lh','sh','p','T','rh','us','vs' ...
    ,'srflwdn','srflwup','srfswdn','srfswup'...
    ,'liq','vap','lwt','swt','lwc','swc','ins'...
    ,'clddz','cldtz','cldhgh','cldmid','cldlow','cldtot' };

% corresponding to the output varname 
% T,q,u,v,omega,Q1,Q2,
% s,divT,divq,divs,vertdivT,vertdivq,vertdivs,WD
% Ptend,dcolH2Odt,colH20Adv,evapsrf,dcolDSEdt,dcolDSEAdt,
% colRH,colLH
% phis,lhflx,shflx,Ps,Tsair,RHsair,usrf,vsrf,
% srflwdn,srflwup,srfswdn,srfswup,
% LWP,PWV, TOA fluxes from satellite,
% cloud information from satellite,


Var_input=cell(1,length(varname_input));
for ii=1:length(varname_input)
    varid = netcdf.inqVarID(ncid,varname_input{ii});
    var=netcdf.getVar(ncid,varid,'double');
    var(var<-9000)=NaN;
    switch ndims(var)
        case 4
            var=flip(var,3);
            var_new = resize_lev_time(var,lev,tsec,lev_new,tsec_new);
            Var_input{ii}=mean(mean(var_new(x_domain,y_domain,:,:),2,'omitnan'),1,'omitnan');
        case 3
           	if length(tsec)==length(tsec_new)
                var_new=var;
            else
                var_new = zeros(lonlen,latlen,timelen_new);
                for yy=1:latlen
                    for xx=1:lonlen
                        var_tmp=interp1(tsec,reshape(var(xx,yy,:),1,[]),tsec_new,'linear');
                        var_new(xx,yy,:)=reshape(var_tmp,1,1,[]);
                    end
                end
            end
            Var_input{ii}=mean(mean(var_new(x_domain,y_domain,:),2,'omitnan'),1,'omitnan');
    end
end
% change some unit of the variables
idx = (strcmp(varname_input,'avg_r'));
Var_input{idx}=Var_input{idx}*cpd/Lv0;  % K to kg/kg
idx = (strcmp(varname_input,'h_advr'));
Var_input{idx}=Var_input{idx}*cpd/Lv0;  % K/s to kg/kg/s
idx = (strcmp(varname_input,'v_advr'));
Var_input{idx}=Var_input{idx}*cpd/Lv0;  % K/s to kg/kg/s
idx = (strcmp(varname_input,'m_divuv'));
Var_input{idx}=Var_input{idx}*-1;  % 
idx = (strcmp(varname_input,'m_Q2'));
Var_input{idx}=Var_input{idx}*-1;  % 
idx = (strcmp(varname_input,'m_Dps_dt'));
Var_input{idx}=Var_input{idx}*-1;
idx = (strcmp(varname_input,'rsd_water'));
Var_input{idx}=Var_input{idx}*cpd/Lv0/1000;        % K*kg/s/m^2 to m/s
idx = (strcmp(varname_input,'m_sdivuvr'));
Var_input{idx}=Var_input{idx}*(-1)*cpd/Lv0/1000;   % K*kg/s/m^2 to m/s
idx = (strcmp(varname_input,'evapor'));
Var_input{idx}=Var_input{idx}*cpd/Lv0/1000;        % K*kg/s/m^2 to m/s
idx = (strcmp(varname_input,'rsd_heat'));
Var_input{idx}=Var_input{idx}*cpd;            % K/s/m^2 to W/m^2
idx = (strcmp(varname_input,'m_sdivuvh'));
Var_input{idx}=Var_input{idx}*(-1)*cpd;       % K/s/m^2 to W/m^2
idx = (strcmp(varname_input,'radiation'));
Var_input{idx}=Var_input{idx}*cpd;       % K/s/m^2 to W/m^2
idx = (strcmp(varname_input,'Lprecip'));
Var_input{idx}=Var_input{idx}*cpd;       % K/s/m^2 to W/m^2
idx = (strcmp(varname_input,'lh'));
Var_input{idx} = Var_input{idx}*-1;   
idx = (strcmp(varname_input,'sh'));
Var_input{idx} = Var_input{idx}*-1;   
idx = (strcmp(varname_input,'zs'));
Var_input{idx}=Var_input{idx}*9.8;    % z(gpm) to phi (m^2/s^2)
idx = (strcmp(varname_input,'liq'));
Var_input{idx}=Var_input{idx}*0.01; % cm to m
idx = (strcmp(varname_input,'T'));
Var_input{idx}=Var_input{idx}+273.15; % C to K


% change surface RH into qs
idx = (strcmp(varname_input,'T'));
Tmp=Var_input{idx};    % temperature, used to calculate q
idx = (strcmp(varname_input,'rh'));
rh_sfc = Var_input{idx};   % RH
idx = (strcmp(varname_input,'p'));
Var_input{idx}=Var_input{idx}*100; % mb to Pa
ps = Var_input{idx};  
es=100*1.003*exp(53.67957-(6743.769./Tmp)-4.8451*log(Tmp));
ev=es.*rh_sfc*0.01;
qv=0*rh_sfc;
for tt=1:timelen_new
    qs_tmp = 287.04/461.5*ev(:,:,tt)./(ps(:,:,tt)-ev(:,:,tt));
    if squeeze(qs_tmp)<0
        disp(qs_tmp);
        error('s');
    end
    qv(:,:,tt)=qs_tmp;
end 
idx = (strcmp(varname_input,'rh'));
Var_input{idx}=qv;

%% re-calculate advections and Q1Q2
% read in omega
varid = netcdf.inqVarID(ncid,'omega_c');
w=netcdf.getVar(ncid,varid,'double');
w(w<-9000)=NaN;
w=mean(mean(w(x_domain,y_domain,:,:),2,'omitnan'),1,'omitnan');
w = flip(w,3); 
varid = netcdf.inqVarID(ncid,'KB');
kb = netcdf.getVar(ncid,varid,'double');

% read in state
varname_state={'u','v','q','s','z','u_xstag','q_xstag','s_xstag',...
    'v_ystag','q_ystag','s_ystag'};

varid = netcdf.inqVarID(ncid,'q');
q = netcdf.getVar(ncid,varid,'double');
varid = netcdf.inqVarID(ncid,'s');
s = netcdf.getVar(ncid,varid,'double');
varid = netcdf.inqVarID(ncid,'z');
z = netcdf.getVar(ncid,varid,'double');
varid = netcdf.inqVarID(ncid,'u');
u = netcdf.getVar(ncid,varid,'double');
varid = netcdf.inqVarID(ncid,'v');
v = netcdf.getVar(ncid,varid,'double');
varid = netcdf.inqVarID(ncid,'u_xstag');
u_x = netcdf.getVar(ncid,varid,'double');
varid = netcdf.inqVarID(ncid,'q_xstag');
q_x = netcdf.getVar(ncid,varid,'double');
varid = netcdf.inqVarID(ncid,'s_xstag');
s_x = netcdf.getVar(ncid,varid,'double');
varid = netcdf.inqVarID(ncid,'v_ystag');
v_y = netcdf.getVar(ncid,varid,'double');
varid = netcdf.inqVarID(ncid,'q_ystag');
q_y = netcdf.getVar(ncid,varid,'double');
varid = netcdf.inqVarID(ncid,'s_ystag');
s_y = netcdf.getVar(ncid,varid,'double');
q = mean(mean(q(x_domain,y_domain,:,:),2,'omitnan'),1,'omitnan');
s = mean(mean(s(x_domain,y_domain,:,:),2,'omitnan'),1,'omitnan');
z = mean(mean(z(x_domain,y_domain,:,:),2,'omitnan'),1,'omitnan');
u = mean(mean(u(x_domain,y_domain,:,:),2,'omitnan'),1,'omitnan');
v = mean(mean(v(x_domain,y_domain,:,:),2,'omitnan'),1,'omitnan');
T = s-z*9.8/cpd;
u_x1 = mean(u_x(x_domain(1),y_domain,:,:),2,'omitnan');
u_x2 = mean(u_x(x_domain(end)+1,y_domain,:,:),2,'omitnan');
q_x1 = mean(q_x(x_domain(1),y_domain,:,:),2,'omitnan');
q_x2 = mean(q_x(x_domain(end)+1,y_domain,:,:),2,'omitnan');
s_x1 = mean(s_x(x_domain(1),y_domain,:,:),2,'omitnan');
s_x2 = mean(s_x(x_domain(end)+1,y_domain,:,:),2,'omitnan');
v_y1 = mean(v_y(x_domain,y_domain(1),:,:),1,'omitnan');
v_y2 = mean(v_y(x_domain,y_domain(end)+1,:,:),1,'omitnan');
q_y1 = mean(q_y(x_domain,y_domain(1),:,:),1,'omitnan');
q_y2 = mean(q_y(x_domain,y_domain(end)+1,:,:),1,'omitnan');
s_y1 = mean(s_y(x_domain,y_domain(1),:,:),1,'omitnan');
s_y2 = mean(s_y(x_domain,y_domain(end)+1,:,:),1,'omitnan');

netcdf.close(ncid);
    
% get dqdt and dsdt
q=flip(q,3);
s=flip(s,3);
z=flip(z,3);
T=flip(T,3);
u=flip(u,3);
v=flip(v,3);
u_x1=flip(u_x1,3);  u_x2=flip(u_x2,3);
q_x1=flip(q_x1,3);  q_x2=flip(q_x2,3);
s_x1=flip(s_x1,3);  s_x2=flip(s_x2,3);
v_y1=flip(v_y1,3);  v_y2=flip(v_y2,3);
q_y1=flip(q_y1,3);  q_y2=flip(q_y2,3);
s_y1=flip(s_y1,3);  s_y2=flip(s_y2,3);
dqdt = zeros(1,1,levlen,timelen);
dsdt = zeros(1,1,levlen,timelen);
dtu = (time(2)-time(1))*60;
for ll=2:timelen
    dqdt(1,1,:,ll)=(q(1,1,:,ll)-q(1,1,:,ll-1))/dtu;
    dsdt(1,1,:,ll)=(s(1,1,:,ll)-s(1,1,:,ll-1))/dtu;
end

% recalculate horizontal and vertical advections
dx = (lon_stag(x_2+1)-lon_stag(x_1))*pi/180*6371000 * cosd(lat_new);
dy = (lat_stag(y_2+1)-lat_stag(y_1))*pi/180*6371000;
dp0 = zeros(1,1,levlen,timelen);
lev1=flip(lev); % from downside up
kb_mean = zeros(1,timelen);
for ll=1:timelen
 	kk=round(mean(mean(kb(:,:,ll))));
    while ps(1,1,ll)<lev1(kk)+10    % ensure ps-lev1(kk) is not too small
        kk=kk+1;
    end
    kb_mean(ll)=kk;
    dp0(1,1,kk,ll)=(ps(1,1,ll)-lev1(kk));
	for kk=kb_mean(ll)+1:(levlen-1)
    	dp0(1,1,kk,ll)=(lev1(kk)-lev1(kk+1));
	end
end
dp0=flip(dp0,3);    % from 2 to levlen-kb_mean(ll)+1
dp=dp0; 

WD = (u_x2-u_x1)./dx + (v_y2-v_y1)./dy;
advq = -u.*(q_x2-q_x1)./dx - v.*(q_y2-q_y1)./dy;
advs = -u.*(s_x2-s_x1)./dx - v.*(s_y2-s_y1)./dy;
advT = advs;

vertadvq = zeros(1,1,levlen,timelen);
vertadvs = zeros(1,1,levlen,timelen);
vertadvT = zeros(1,1,levlen,timelen);

% forward difference
% for ll=2:timelen
%     for kk=2:(levlen-kb_mean(ll)+1)
%         vertadvq(1,1,kk,ll) = -w(1,1,kk,ll)* ...
%             (q(1,1,kk,ll)-q(1,1,kk-1,ll))./dp(1,1,kk,ll);
%         vertadvs(1,1,kk,ll) = -w(1,1,kk,ll)* ...
%             (s(1,1,kk,ll)-s(1,1,kk-1,ll))./dp(1,1,kk,ll);
%         vertadvT(1,1,kk,ll) = -w(1,1,kk,ll)* ...
%             (T(1,1,kk,ll)-T(1,1,kk-1,ll))./dp(1,1,kk,ll);
%     end
% end
% central difference
for ll=2:timelen
    for kk=2:(levlen-kb_mean(ll)+1)
        vertadvq(1,1,kk,ll) = -0.5*w(1,1,kk,ll)* ...
            (q(1,1,kk+1,ll)-q(1,1,kk-1,ll))./dp(1,1,kk,ll);
        vertadvs(1,1,kk,ll) = -0.5*w(1,1,kk,ll)* ...
            (s(1,1,kk+1,ll)-s(1,1,kk-1,ll))./dp(1,1,kk,ll);
        vertadvT(1,1,kk,ll) = -0.5*w(1,1,kk,ll)* ...
            (T(1,1,kk+1,ll)-T(1,1,kk-1,ll))./dp(1,1,kk,ll);
    end
end
Q1 = dsdt - advs - vertadvs;
Q2 = ( - dqdt + advq + vertadvq)*Lv0/cpd;

T = resize_lev_time(T,lev,tsec,lev_new,tsec_new);
Q1 = resize_lev_time(Q1,lev,tsec,lev_new,tsec_new);
Q2 = resize_lev_time(Q2,lev,tsec,lev_new,tsec_new);
advq = resize_lev_time(advq,lev,tsec,lev_new,tsec_new);
advs = resize_lev_time(advs,lev,tsec,lev_new,tsec_new);
advT = resize_lev_time(advT,lev,tsec,lev_new,tsec_new);
vertadvq = resize_lev_time(vertadvq,lev,tsec,lev_new,tsec_new);
vertadvs = resize_lev_time(vertadvs,lev,tsec,lev_new,tsec_new);
vertadvT = resize_lev_time(vertadvT,lev,tsec,lev_new,tsec_new);
WD = resize_lev_time(WD,lev,tsec,lev_new,tsec_new);

idx = (strcmp(varname_input,'avg_T'));
Var_input{idx}=T; 
idx = (strcmp(varname_input,'Q1'));
Var_input{idx}=Q1;
idx = (strcmp(varname_input,'m_Q2'));
Var_input{idx}=Q2;
idx = (strcmp(varname_input,'h_advr'));
Var_input{idx}=advq;
idx = (strcmp(varname_input,'h_advE'));
Var_input{idx}=advs;
idx = (strcmp(varname_input,'h_advT'));
Var_input{idx}=advT;
idx = (strcmp(varname_input,'v_advT'));
Var_input{idx}=vertadvT;
idx = (strcmp(varname_input,'v_advr'));
Var_input{idx}=vertadvq;
idx = (strcmp(varname_input,'v_advE'));
Var_input{idx}=vertadvs;
idx = (strcmp(varname_input,'m_divuv'));
Var_input{idx}=WD;


%% ---------- write data ---------------
ncid = netcdf.create(opt_file,'NC_WRITE');
% dimension
londimid = netcdf.defDim(ncid,'lon',1);
latdimid = netcdf.defDim(ncid,'lat',1);
levdimid = netcdf.defDim(ncid,'lev',levlen_new);
tdimid  = netcdf.defDim(ncid,'time',timelen_new);
% variable
lonid = netcdf.defVar(ncid,'lon','double',londimid);
latid = netcdf.defVar(ncid,'lat','double',latdimid);
levid = netcdf.defVar(ncid,'lev','double',levdimid);
time2id = netcdf.defVar(ncid,'time','double',tdimid);
timeid = netcdf.defVar(ncid,'tsec','int',tdimid);
bdateid = netcdf.defVar(ncid,'bdate','int',[]);
phiid = netcdf.defVar(ncid,'phis','double',[]);
lhid = netcdf.defVar(ncid,'lhflx','double',[londimid,latdimid,tdimid]);
shid = netcdf.defVar(ncid,'shflx','double',[londimid,latdimid,tdimid]);
precid = netcdf.defVar(ncid,'Prec','double',[londimid,latdimid,tdimid]);
psid = netcdf.defVar(ncid,'Ps','double',[londimid,latdimid,tdimid]);
Tsid = netcdf.defVar(ncid,'Tsair','double',[londimid,latdimid,tdimid]);
qsid = netcdf.defVar(ncid,'qs','double',[londimid,latdimid,tdimid]);
usid = netcdf.defVar(ncid,'usrf','double',[londimid,latdimid,tdimid]);
vsid = netcdf.defVar(ncid,'vsrf','double',[londimid,latdimid,tdimid]);
slwdnid = netcdf.defVar(ncid,'srflwdn','double',[londimid,latdimid,tdimid]);
slwupid = netcdf.defVar(ncid,'srflwup','double',[londimid,latdimid,tdimid]);
sswdnid = netcdf.defVar(ncid,'srfswdn','double',[londimid,latdimid,tdimid]);
sswupid = netcdf.defVar(ncid,'srfswup','double',[londimid,latdimid,tdimid]);
lcldid = netcdf.defVar(ncid,'lowcld','double',[londimid,latdimid,tdimid]);
mcldid = netcdf.defVar(ncid,'midcld','double',[londimid,latdimid,tdimid]);
hcldid = netcdf.defVar(ncid,'hghcld','double',[londimid,latdimid,tdimid]);
tcldid = netcdf.defVar(ncid,'totcld','double',[londimid,latdimid,tdimid]);
cldthkid = netcdf.defVar(ncid,'cldthk','double',[londimid,latdimid,tdimid]);
cldhtid = netcdf.defVar(ncid,'cldht','double',[londimid,latdimid,tdimid]);
liqid = netcdf.defVar(ncid,'LWP','double',[londimid,latdimid,tdimid]);
PWVid = netcdf.defVar(ncid,'PWV','double',[londimid,latdimid,tdimid]);
lwtid = netcdf.defVar(ncid,'lwt','double',[londimid,latdimid,tdimid]);
swtid = netcdf.defVar(ncid,'swt','double',[londimid,latdimid,tdimid]);
lwcid = netcdf.defVar(ncid,'lwc','double',[londimid,latdimid,tdimid]);
swcid = netcdf.defVar(ncid,'swc','double',[londimid,latdimid,tdimid]);
insid = netcdf.defVar(ncid,'ins','double',[londimid,latdimid,tdimid]);

ptendid = netcdf.defVar(ncid,'Ptend','double',[londimid,latdimid,tdimid]);
dsqdtid = netcdf.defVar(ncid,'dcolH2Odt','double',[londimid,latdimid,tdimid]);
sdivqid = netcdf.defVar(ncid,'colH20Adv','double',[londimid,latdimid,tdimid]);
evapid = netcdf.defVar(ncid,'evapsrf','double',[londimid,latdimid,tdimid]);
dssdtid = netcdf.defVar(ncid,'dcolDSEdt','double',[londimid,latdimid,tdimid]);
sdivsid = netcdf.defVar(ncid,'dcolDSEAdt','double',[londimid,latdimid,tdimid]);
radid = netcdf.defVar(ncid,'colRH','double',[londimid,latdimid,tdimid]);
LPid = netcdf.defVar(ncid,'colLH','double',[londimid,latdimid,tdimid]);

Tid = netcdf.defVar(ncid,'T','double',[londimid,latdimid,levdimid,tdimid]);
qid = netcdf.defVar(ncid,'q','double',[londimid,latdimid,levdimid,tdimid]);
uid = netcdf.defVar(ncid,'u','double',[londimid,latdimid,levdimid,tdimid]);
vid = netcdf.defVar(ncid,'v','double',[londimid,latdimid,levdimid,tdimid]);
zid = netcdf.defVar(ncid,'geoz','double',[londimid,latdimid,levdimid,tdimid]);
wid = netcdf.defVar(ncid,'omega','double',[londimid,latdimid,levdimid,tdimid]);
q1id = netcdf.defVar(ncid,'Q1','double',[londimid,latdimid,levdimid,tdimid]);
q2id = netcdf.defVar(ncid,'Q2','double',[londimid,latdimid,levdimid,tdimid]);
sid = netcdf.defVar(ncid,'s','double',[londimid,latdimid,levdimid,tdimid]);
hadvTid = netcdf.defVar(ncid,'divT','double',[londimid,latdimid,levdimid,tdimid]);
hadvqid = netcdf.defVar(ncid,'divq','double',[londimid,latdimid,levdimid,tdimid]);
hadvsid = netcdf.defVar(ncid,'divs','double',[londimid,latdimid,levdimid,tdimid]);
vadvTid = netcdf.defVar(ncid,'vertdivT','double',[londimid,latdimid,levdimid,tdimid]);
vadvqid = netcdf.defVar(ncid,'vertdivq','double',[londimid,latdimid,levdimid,tdimid]);
vadvsid = netcdf.defVar(ncid,'vertdivs','double',[londimid,latdimid,levdimid,tdimid]);
wdid = netcdf.defVar(ncid,'WD','double',[londimid,latdimid,levdimid,tdimid]);
netcdf.endDef(ncid);

netcdf.putVar(ncid,lonid,lon_new);
netcdf.putVar(ncid,latid,lat_new);
netcdf.putVar(ncid,levid,lev_new);
netcdf.putVar(ncid,time2id,tsec_new);
netcdf.putVar(ncid,timeid,tsec_new);
netcdf.putVar(ncid,bdateid,bdate(1));
idx=(strcmp(varname_input,'zs'));
netcdf.putVar(ncid,phiid,mean(squeeze(Var_input{idx})));
idx=(strcmp(varname_input,'lh'));
netcdf.putVar(ncid,lhid,Var_input{idx});
idx=(strcmp(varname_input,'sh'));
netcdf.putVar(ncid,shid,Var_input{idx});
idx=(strcmp(varname_input,'p'));
netcdf.putVar(ncid,psid,Var_input{idx});
idx=(strcmp(varname_input,'T'));
netcdf.putVar(ncid,Tsid,Var_input{idx});
idx=(strcmp(varname_input,'rh'));
netcdf.putVar(ncid,qsid,Var_input{idx});
idx=(strcmp(varname_input,'us'));
netcdf.putVar(ncid,usid,Var_input{idx});
idx=(strcmp(varname_input,'vs'));
netcdf.putVar(ncid,vsid,Var_input{idx});
idx=(strcmp(varname_input,'srflwdn'));
netcdf.putVar(ncid,slwdnid,Var_input{idx});
idx=(strcmp(varname_input,'srflwup'));
netcdf.putVar(ncid,slwupid,Var_input{idx});
idx=(strcmp(varname_input,'srfswdn'));
netcdf.putVar(ncid,sswdnid,Var_input{idx});
idx=(strcmp(varname_input,'srfswup'));
netcdf.putVar(ncid,sswupid,Var_input{idx});
idx=(strcmp(varname_input,'cldlow'));
netcdf.putVar(ncid,lcldid,Var_input{idx});    
idx=(strcmp(varname_input,'cldmid'));
netcdf.putVar(ncid,mcldid,Var_input{idx});  
idx=(strcmp(varname_input,'cldhgh'));
netcdf.putVar(ncid,hcldid,Var_input{idx});  
idx=(strcmp(varname_input,'cldtot'));
netcdf.putVar(ncid,tcldid,Var_input{idx});  
idx=(strcmp(varname_input,'clddz'));
netcdf.putVar(ncid,cldthkid,Var_input{idx});  
idx=(strcmp(varname_input,'cldtz'));
netcdf.putVar(ncid,cldhtid,Var_input{idx});   
idx=(strcmp(varname_input,'liq'));
netcdf.putVar(ncid,liqid,Var_input{idx});
idx=(strcmp(varname_input,'vap'));
netcdf.putVar(ncid,PWVid,Var_input{idx});
idx=(strcmp(varname_input,'lwt'));
netcdf.putVar(ncid,lwtid,Var_input{idx});
idx=(strcmp(varname_input,'swt'));
netcdf.putVar(ncid,swtid,Var_input{idx});
idx=(strcmp(varname_input,'lwc'));
netcdf.putVar(ncid,lwcid,Var_input{idx});  
idx=(strcmp(varname_input,'swc'));
netcdf.putVar(ncid,swcid,Var_input{idx});  
idx=(strcmp(varname_input,'ins'));
netcdf.putVar(ncid,insid,Var_input{idx});

idx=(strcmp(varname_input,'m_Dps_dt'));
netcdf.putVar(ncid,ptendid,Var_input{idx});
idx=(strcmp(varname_input,'rsd_water'));
netcdf.putVar(ncid,dsqdtid,Var_input{idx});
idx=(strcmp(varname_input,'m_sdivuvr'));
netcdf.putVar(ncid,sdivqid,Var_input{idx});
idx=(strcmp(varname_input,'evapor'));
netcdf.putVar(ncid,evapid,Var_input{idx});
idx=(strcmp(varname_input,'rsd_heat'));
netcdf.putVar(ncid,dssdtid,Var_input{idx});
idx=(strcmp(varname_input,'m_sdivuvh'));
netcdf.putVar(ncid,sdivsid,Var_input{idx});
idx=(strcmp(varname_input,'radiation'));
netcdf.putVar(ncid,radid,Var_input{idx});
idx=(strcmp(varname_input,'Lprecip'));
netcdf.putVar(ncid,LPid,Var_input{idx});
netcdf.putVar(ncid,precid,Var_input{idx}/Lv0/1000);   % precip: m/s 

idx=(strcmp(varname_input,'avg_T'));
netcdf.putVar(ncid,Tid,Var_input{idx});
idx=(strcmp(varname_input,'avg_r'));
netcdf.putVar(ncid,qid,Var_input{idx});
idx=(strcmp(varname_input,'avg_u'));
netcdf.putVar(ncid,uid,Var_input{idx});
idx=(strcmp(varname_input,'avg_v'));
netcdf.putVar(ncid,vid,Var_input{idx});
idx=(strcmp(varname_input,'avg_z'));
netcdf.putVar(ncid,zid,Var_input{idx});
idx=(strcmp(varname_input,'omega_c'));
netcdf.putVar(ncid,wid,Var_input{idx});
idx=(strcmp(varname_input,'Q1'));
netcdf.putVar(ncid,q1id,Var_input{idx});
idx=(strcmp(varname_input,'m_Q2'));
netcdf.putVar(ncid,q2id,Var_input{idx});
idx=(strcmp(varname_input,'avg_E'));
netcdf.putVar(ncid,sid,Var_input{idx});
idx=(strcmp(varname_input,'h_advT'));
netcdf.putVar(ncid,hadvTid,Var_input{idx});
idx=(strcmp(varname_input,'h_advr'));
netcdf.putVar(ncid,hadvqid,Var_input{idx});
idx=(strcmp(varname_input,'h_advE'));
netcdf.putVar(ncid,hadvsid,Var_input{idx});
idx=(strcmp(varname_input,'v_advT'));
netcdf.putVar(ncid,vadvTid,Var_input{idx});
idx=(strcmp(varname_input,'v_advr'));
netcdf.putVar(ncid,vadvqid,Var_input{idx});
idx=(strcmp(varname_input,'v_advE'));
netcdf.putVar(ncid,vadvsid,Var_input{idx});
idx=(strcmp(varname_input,'m_divuv'));
netcdf.putVar(ncid,wdid,Var_input{idx});

% parameters
netcdf.reDef(ncid);
netcdf.putAtt(ncid,lonid,'units','degrees_east');
netcdf.putAtt(ncid,latid,'units','degrees_north');
netcdf.putAtt(ncid,levid,'units','Pa');
netcdf.putAtt(ncid,time2id,'units',initial_time);
netcdf.putAtt(ncid,timeid,'long_name','Time in seconds after 0Z on base date');
netcdf.putAtt(ncid,timeid,'units','s');
netcdf.putAtt(ncid,bdateid,'long_name','base date as 8 digit integer');
netcdf.putAtt(ncid,bdateid,'units','yyyymmdd');
netcdf.putAtt(ncid,phiid,'long_name','surface geopotential');
netcdf.putAtt(ncid,phiid,'units','m^2/s^2');
netcdf.putAtt(ncid,lhid,'long_name','surface latent heat flux');
netcdf.putAtt(ncid,lhid,'units','W/m2');
netcdf.putAtt(ncid,lhid,'missing_value',missing_value);
netcdf.putAtt(ncid,shid,'long_name','surface sensible heat flux');
netcdf.putAtt(ncid,shid,'units','W/m2');
netcdf.putAtt(ncid,shid,'missing_value',missing_value);
netcdf.putAtt(ncid,precid,'long_name','precipitation');
netcdf.putAtt(ncid,precid,'units','m/s');
netcdf.putAtt(ncid,precid,'missing_value',missing_value);
netcdf.putAtt(ncid,psid,'long_name','area mean surface pressure');
netcdf.putAtt(ncid,psid,'units','Pa');
netcdf.putAtt(ncid,psid,'missing_value',missing_value);
netcdf.putAtt(ncid,Tsid,'long_name','surface temperature');
netcdf.putAtt(ncid,Tsid,'units','K');
netcdf.putAtt(ncid,Tsid,'missing_value',missing_value);
netcdf.putAtt(ncid,qsid,'long_name','surface specific humidity');
netcdf.putAtt(ncid,qsid,'units','kg/kg');
netcdf.putAtt(ncid,qsid,'missing_value',missing_value);
netcdf.putAtt(ncid,usid,'long_name','surface u wind');
netcdf.putAtt(ncid,usid,'units','m/s');
netcdf.putAtt(ncid,usid,'missing_value',missing_value);
netcdf.putAtt(ncid,vsid,'long_name','surface v wind');
netcdf.putAtt(ncid,vsid,'units','m/s');
netcdf.putAtt(ncid,vsid,'missing_value',missing_value);
netcdf.putAtt(ncid,slwdnid,'long_name','surface longwave downward');
netcdf.putAtt(ncid,slwdnid,'units','W/m2');
netcdf.putAtt(ncid,slwdnid,'missing_value',missing_value);
netcdf.putAtt(ncid,slwupid,'long_name','surface longwave upward');
netcdf.putAtt(ncid,slwupid,'units','W/m2');
netcdf.putAtt(ncid,slwupid,'missing_value',missing_value);
netcdf.putAtt(ncid,sswdnid,'long_name','surface shortwave downward');
netcdf.putAtt(ncid,sswdnid,'units','W/m2');
netcdf.putAtt(ncid,sswdnid,'missing_value',missing_value);
netcdf.putAtt(ncid,sswupid,'long_name','surface shortwave upward');
netcdf.putAtt(ncid,sswupid,'units','W/m2');
netcdf.putAtt(ncid,sswupid,'missing_value',missing_value);
netcdf.putAtt(ncid,lcldid,'long_name','low cloud fraction');
netcdf.putAtt(ncid,lcldid,'units','%');
netcdf.putAtt(ncid,lcldid,'missing_value',missing_value);
netcdf.putAtt(ncid,mcldid,'long_name','mid cloud fraction');
netcdf.putAtt(ncid,mcldid,'units','%');
netcdf.putAtt(ncid,mcldid,'missing_value',missing_value);
netcdf.putAtt(ncid,hcldid,'long_name','high cloud fraction');
netcdf.putAtt(ncid,hcldid,'units','%');
netcdf.putAtt(ncid,hcldid,'missing_value',missing_value);
netcdf.putAtt(ncid,tcldid,'long_name','total cloud fraction');
netcdf.putAtt(ncid,tcldid,'units','%');
netcdf.putAtt(ncid,tcldid,'missing_value',missing_value);
netcdf.putAtt(ncid,cldthkid,'long_name','cloud thickness');
netcdf.putAtt(ncid,cldthkid,'units','km');
netcdf.putAtt(ncid,cldthkid,'missing_value',missing_value);
netcdf.putAtt(ncid,cldhtid,'long_name','cloud top height');
netcdf.putAtt(ncid,cldhtid,'units','km');
netcdf.putAtt(ncid,cldhtid,'missing_value',missing_value);
netcdf.putAtt(ncid,liqid,'long_name','cloud liquid water path');
netcdf.putAtt(ncid,liqid,'units','m');
netcdf.putAtt(ncid,liqid,'missing_value',missing_value);
netcdf.putAtt(ncid,PWVid,'long_name','total precipitable water');
netcdf.putAtt(ncid,PWVid,'units','cm');
netcdf.putAtt(ncid,PWVid,'missing_value',missing_value);
netcdf.putAtt(ncid,lwtid,'long_name','total sky longwave flux at TOA');
netcdf.putAtt(ncid,lwtid,'units','W/m^2');
netcdf.putAtt(ncid,lwtid,'missing_value',missing_value);
netcdf.putAtt(ncid,swtid,'long_name','total sky shortwave flux at TOA');
netcdf.putAtt(ncid,swtid,'units','W/m^2');
netcdf.putAtt(ncid,swtid,'missing_value',missing_value);
netcdf.putAtt(ncid,lwcid,'long_name','clear sky longwave flux at TOA');
netcdf.putAtt(ncid,lwcid,'units','W/m^2');
netcdf.putAtt(ncid,lwcid,'missing_value',missing_value);
netcdf.putAtt(ncid,swcid,'long_name','clear sky shortwave flux at TOA');
netcdf.putAtt(ncid,swcid,'units','W/m^2');
netcdf.putAtt(ncid,swcid,'missing_value',missing_value);
netcdf.putAtt(ncid,insid,'long_name','total insolation flux at TOA');
netcdf.putAtt(ncid,insid,'units','W/m^2');
netcdf.putAtt(ncid,insid,'missing_value',missing_value);

netcdf.putAtt(ncid,ptendid,'long_name','surface pressure tendency');
netcdf.putAtt(ncid,ptendid,'units','Pa/s');
netcdf.putAtt(ncid,ptendid,'missing_value',missing_value);
netcdf.putAtt(ncid,dsqdtid,'long_name','column total tendency of q');
netcdf.putAtt(ncid,dsqdtid,'units','m/s');
netcdf.putAtt(ncid,dsqdtid,'missing_value',missing_value);
netcdf.putAtt(ncid,sdivqid,'long_name','column total advection of q');
netcdf.putAtt(ncid,sdivqid,'units','m/s');
netcdf.putAtt(ncid,sdivqid,'missing_value',missing_value);
netcdf.putAtt(ncid,evapid,'long_name','surface evaporation');
netcdf.putAtt(ncid,evapid,'units','m/s');
netcdf.putAtt(ncid,evapid,'missing_value',missing_value);
netcdf.putAtt(ncid,dssdtid,'long_name','column total tendency of s');
netcdf.putAtt(ncid,dssdtid,'units','W/m2');
netcdf.putAtt(ncid,dssdtid,'missing_value',missing_value);
netcdf.putAtt(ncid,sdivsid,'long_name','column total advection of s');
netcdf.putAtt(ncid,sdivsid,'units','W/m2');
netcdf.putAtt(ncid,sdivsid,'missing_value',missing_value);
netcdf.putAtt(ncid,radid,'long_name','column total radiative heating');
netcdf.putAtt(ncid,radid,'units','W/m2');
netcdf.putAtt(ncid,radid,'missing_value',missing_value);
netcdf.putAtt(ncid,LPid,'long_name','column total latend heating');
netcdf.putAtt(ncid,LPid,'units','W/m2');
netcdf.putAtt(ncid,LPid,'missing_value',missing_value);

netcdf.putAtt(ncid,Tid,'long_name','temperature');
netcdf.putAtt(ncid,Tid,'units','K');
netcdf.putAtt(ncid,Tid,'missing_value',missing_value);
netcdf.putAtt(ncid,qid,'long_name','water vapor mixing ratio');
netcdf.putAtt(ncid,qid,'units','kg/kg');
netcdf.putAtt(ncid,qid,'missing_value',missing_value);
netcdf.putAtt(ncid,uid,'long_name','u wind');
netcdf.putAtt(ncid,uid,'units','m/s');
netcdf.putAtt(ncid,uid,'missing_value',missing_value);
netcdf.putAtt(ncid,vid,'long_name','v wind');
netcdf.putAtt(ncid,vid,'units','m/s');
netcdf.putAtt(ncid,vid,'missing_value',missing_value);
netcdf.putAtt(ncid,wid,'long_name','vertical velocity');
netcdf.putAtt(ncid,wid,'units','Pa/s');
netcdf.putAtt(ncid,wid,'missing_value',missing_value);
netcdf.putAtt(ncid,zid,'long_name','geopotential height');
netcdf.putAtt(ncid,zid,'units','gpm');
netcdf.putAtt(ncid,zid,'missing_value',missing_value);
netcdf.putAtt(ncid,q1id,'long_name','apparent heating source');
netcdf.putAtt(ncid,q1id,'units','K/s');
netcdf.putAtt(ncid,q1id,'missing_value',missing_value);
netcdf.putAtt(ncid,q2id,'long_name','apparent moisture sinks');
netcdf.putAtt(ncid,q2id,'units','K/s');
netcdf.putAtt(ncid,q2id,'missing_value',missing_value);
netcdf.putAtt(ncid,sid,'long_name','dry static energy');
netcdf.putAtt(ncid,sid,'units','K');
netcdf.putAtt(ncid,sid,'missing_value',missing_value);
netcdf.putAtt(ncid,hadvTid,'long_name','horizontal temperature advections');
netcdf.putAtt(ncid,hadvTid,'units','K/s');
netcdf.putAtt(ncid,hadvTid,'missing_value',missing_value);
netcdf.putAtt(ncid,hadvqid,'long_name','horizontal water vapor advections');
netcdf.putAtt(ncid,hadvqid,'units','kg/kg/s');
netcdf.putAtt(ncid,hadvqid,'missing_value',missing_value);
netcdf.putAtt(ncid,hadvsid,'long_name','horizontal s advections');
netcdf.putAtt(ncid,hadvsid,'units','K/s');
netcdf.putAtt(ncid,hadvsid,'missing_value',missing_value);
netcdf.putAtt(ncid,vadvTid,'long_name','vertical temperature advections');
netcdf.putAtt(ncid,vadvTid,'units','K/s');
netcdf.putAtt(ncid,vadvTid,'missing_value',missing_value);
netcdf.putAtt(ncid,vadvqid,'long_name','vertical water vapor advections');
netcdf.putAtt(ncid,vadvqid,'units','kg/kg/s');
netcdf.putAtt(ncid,vadvqid,'missing_value',missing_value);
netcdf.putAtt(ncid,vadvsid,'long_name','vertical s advections');
netcdf.putAtt(ncid,vadvsid,'units','K/s');
netcdf.putAtt(ncid,vadvsid,'missing_value',missing_value);
netcdf.putAtt(ncid,wdid,'long_name','horizontal wind divergence');
netcdf.putAtt(ncid,wdid,'units','s-1');
netcdf.putAtt(ncid,wdid,'missing_value',missing_value);

% global attributes
varid = netcdf.getConstant('GLOBAL');
netcdf.putAtt(ncid,varid,'title','Forcing data generated from 3DCVA');
netcdf.putAtt(ncid,varid,'author','Shuaiqi Tang: tang32@llnl.gov');
netcdf.putAtt(ncid,varid,'program name','generate_forcing_3DCVA.m');
netcdf.putAtt(ncid,varid,'input data',ipt_file);
netcdf.putAtt(ncid,varid,'references','Tang and Zhang (2015). https://doi.org/10.1002/2015JD023621');
netcdf.close(ncid);
