clear all 
close all
format shortg

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Chris Foster
% 11/12/2015
% University of Utah
%
% Purpose: Edit geo_em file land use to remove portions of the GSL, the
% seasonality of which is not captured by static land use.
%
% How: This is done by reading all the static variables that are expected 
% to be affected by the presence of water into matlab, editing them based
% qualitatively on satellite imagery, and overwriting the original geo_em  
% with the new (edited) variables. The edits are done seperately to each 
% portion of the lake that requires modification on a pixel by pixel basis.
%
% Note on editing: This is done in loops that loop through and change each
% pixel in a box that encompasses an area of the inner domain that is
% being changed. The pixels' properties are being changed to that of barren
% land, a pixel chosen near the lake and assumed to be represntative on
% non-water conditions around the lake.
%
% Important: This script needs to be modified and rerun for each domain you
% are modifying the lake in.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% ncread the appropriate variables into matlab arrays.

run_path = '../kingspeak/wrf-nlcd2011/';
domain = 'd01'; % Change this to modify each of your domains. (d01, d02, d03)


% Variables that won't be modified.
lat = ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'XLAT_M');
lon = ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'XLONG_M');
hgt = ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'HGT_M');

% Variables to modify.
lnd =   ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'LU_INDEX');
flg40 = ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'LANDMASK');
frc =   ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'LANDUSEF');
flg =   ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'LANDMASK');
alb =   ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'ALBEDO12M');
slt =   ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'SOILTEMP');
stp =   ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'SOILCTOP');
sct =   ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'SCT_DOM');
sbt =   ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'SOILCBOT');
scb =   ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'SCB_DOM');
grn =   ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'GREENFRAC');
snalb = ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'SNOALB');
lkdep = ncread(strcat(run_path,'WPS/geo_em_new.',domain,'.nc'),'LAKE_DEPTH');

lon  = double(lon(:,:,1));
lat  = double(lat(:,:,1));

% For plotting purposes.

latlim = [min(lat(:)) max(lat(:))];
lonlim = [min(lon(:)) max(lon(:))];

%% Set boundaries for GLS modification.

% These min/max values are the grid points within your model domains where
% you want the modifications to take place. 

% Domain 3, 1.33 km

% x_min = 159+25; 
% x_max = 239+15;
% y_min = 305+35;
% y_max = 425+6;

% Domain 2, 4.0 km

% x_min = 105;
% x_max = 130;
% y_min = 154;
% y_max = 185;

% Domain 1, 12.0 km

x_min = 79;
x_max = 88;
y_min = 105;
y_max = 117;

%% Initialize new, smaller domain for purposes of changing land use vars.

flg_orig = flg40(x_min:x_max,y_min:y_max);

new_lnd = lnd(x_min:x_max,y_min:y_max);
new_frc = frc(x_min:x_max,y_min:y_max,:);
new_flg = flg40(x_min:x_max,y_min:y_max);
new_slt = slt(x_min:x_max,y_min:y_max);
new_stp = stp(x_min:x_max,y_min:y_max,:);
new_sct = sct(x_min:x_max,y_min:y_max);
new_sbt = sbt(x_min:x_max,y_min:y_max,:);
new_scb = scb(x_min:x_max,y_min:y_max);
new_grn = grn(x_min:x_max,y_min:y_max,:);
new_alb = alb(x_min:x_max,y_min:y_max,:);
new_snalb = snalb(x_min:x_max,y_min:y_max);
new_lkdep = lkdep(x_min:x_max,y_min:y_max);

% Check to make sure your 'lake domain' is the correct size and fully
% encompasses the GSL.

figure(1)
pcolor(new_lnd')
title('Dominant Land Use Before Modification')

%% Each loop is a different edit, described as a comment beforehand.
%  Domain 3

% % Lose north west seperate water body.

for ns = 1:length(new_flg(1,:))     % ns (north to south)
    for ew = 1:13                   % ew (east to west)
        if flg_orig(ew,ns) == 0
            new_flg(ew,ns) = new_flg(27,30);
            new_lnd(ew,ns) = new_lnd(27,30);   
            new_frc(ew,ns,:) = new_frc(27,30,:);
            new_alb(ew,ns,:) = new_alb(27,30,:);
            new_slt(ew,ns) = new_slt(27,30);
            new_stp(ew,ns,:) = new_stp(27,30,:);
            new_sct(ew,ns) = new_sct(27,30);
            new_sbt(ew,ns,:) = new_sbt(27,30,:);
            new_scb(ew,ns) = new_scb(27,30);
            new_grn(ew,ns,:) = new_grn(27,30,:);
            new_snalb(ew,ns) = new_snalb(27,30);
            new_lkdep(ew,ns) = new_lkdep(27,30);
        end
    end
end

% % Lose Farmington Bay.
% % Make barren land classification (pixel we know is barren land).

for ns = 17:42     % ns (north to south)
    for ew = 50:71                   % ew (east to west)
        if flg_orig(ew,ns) == 0 
            new_flg(ew,ns) = new_flg(27,30);
            new_lnd(ew,ns) = new_lnd(27,30);   
            new_frc(ew,ns,:) = new_frc(27,30,:);
            new_alb(ew,ns,:) = new_alb(27,30,:);
            new_slt(ew,ns) = new_slt(27,30);
            new_stp(ew,ns,:) = new_stp(27,30,:);
            new_sct(ew,ns) = new_sct(27,30);
            new_sbt(ew,ns,:) = new_sbt(27,30,:);
            new_scb(ew,ns) = new_scb(27,30);
            new_grn(ew,ns,:) = new_grn(27,30,:);
            new_snalb(ew,ns) = new_snalb(27,30);
            new_lkdep(ew,ns) = new_lkdep(27,30);
        end
    end
end

% % Lose lower southwest portion of seperate water body.

for ns = 1:20     
    for ew = 1:33                 
        if flg_orig(ew,ns) == 0
            new_flg(ew,ns) = new_flg(27,30);
            new_lnd(ew,ns) = new_lnd(27,30);   
            new_frc(ew,ns,:) = new_frc(27,30,:);
            new_alb(ew,ns,:) = new_alb(27,30,:);
            new_slt(ew,ns) = new_slt(27,30);
            new_stp(ew,ns,:) = new_stp(27,30,:);
            new_sct(ew,ns) = new_sct(27,30);
            new_sbt(ew,ns,:) = new_sbt(27,30,:);
            new_scb(ew,ns) = new_scb(27,30);
            new_grn(ew,ns,:) = new_grn(27,30,:);
            new_snalb(ew,ns) = new_snalb(27,30);
            new_lkdep(ew,ns) = new_lkdep(27,30);
        end
    end
end

% % Fix Willard Bay.

for sn = 43:80
    for ew = 45:70
        if flg_orig(ew,sn) == 0 
            new_flg(ew,sn) = new_flg(27,30);
            new_lnd(ew,sn) = new_lnd(27,30);   
            new_frc(ew,sn,:) = new_frc(27,30,:);
            new_alb(ew,sn,:) = new_alb(27,30,:);
            new_slt(ew,sn) = new_slt(27,30);
            new_stp(ew,sn,:) = new_stp(27,30,:);
            new_sct(ew,sn) = new_sct(27,30);
            new_sbt(ew,sn,:) = new_sbt(27,30,:);
            new_scb(ew,sn) = new_scb(27,30);
            new_grn(ew,sn,:) = new_grn(27,30,:);
            new_snalb(ew,sn) = new_snalb(27,30);
            new_lkdep(ew,sn) = new_lkdep(27,30);
        end
    end
end

%% Domain 2

% % Lose Farmington Bay.
% % Make barren land classification (pixel we know is barren land).

for ns = 5:30     % ns (north to south)
    for ew = 18:25                   % ew (east to west)
        if flg_orig(ew,ns) == 0 
            new_flg(ew,ns) = new_flg(10,10);
            new_lnd(ew,ns) = new_lnd(10,10);   
            new_frc(ew,ns,:) = new_frc(10,10,:);
            new_alb(ew,ns,:) = new_alb(10,10,:);
            new_slt(ew,ns) = new_slt(10,10);
            new_stp(ew,ns,:) = new_stp(10,10,:);
            new_sct(ew,ns) = new_sct(10,10);
            new_sbt(ew,ns,:) = new_sbt(10,10,:);
            new_scb(ew,ns) = new_scb(10,10);
            new_grn(ew,ns,:) = new_grn(10,10,:);
            new_snalb(ew,ns) = new_snalb(10,10);
            new_lkdep(ew,ns) = new_lkdep(10,10);
        end
    end
end

% % Lose Willard Bay.
% % Make barren land classification (pixel we know is barren land).

for ns = 15:30     % ns (north to south)
    for ew = 16:25                   % ew (east to west)
        if flg_orig(ew,ns) == 0 
            new_flg(ew,ns) = new_flg(10,10);
            new_lnd(ew,ns) = new_lnd(10,10);   
            new_frc(ew,ns,:) = new_frc(10,10,:);
            new_alb(ew,ns,:) = new_alb(10,10,:);
            new_slt(ew,ns) = new_slt(10,10);
            new_stp(ew,ns,:) = new_stp(10,10,:);
            new_sct(ew,ns) = new_sct(10,10);
            new_sbt(ew,ns,:) = new_sbt(10,10,:);
            new_scb(ew,ns) = new_scb(10,10);
            new_grn(ew,ns,:) = new_grn(10,10,:);
            new_snalb(ew,ns) = new_snalb(10,10);
            new_lkdep(ew,ns) = new_lkdep(10,10);
        end
    end
end

% % Lose Northwest Seperate Body of Water
% % Make barren land classification (pixel we know is barren land).

for ns = 15:30     % ns (north to south)
    for ew = 1:4                   % ew (east to west)
        if flg_orig(ew,ns) == 0 
            new_flg(ew,ns) = new_flg(10,10);
            new_lnd(ew,ns) = new_lnd(10,10);   
            new_frc(ew,ns,:) = new_frc(10,10,:);
            new_alb(ew,ns,:) = new_alb(10,10,:);
            new_slt(ew,ns) = new_slt(10,10);
            new_stp(ew,ns,:) = new_stp(10,10,:);
            new_sct(ew,ns) = new_sct(10,10);
            new_sbt(ew,ns,:) = new_sbt(10,10,:);
            new_scb(ew,ns) = new_scb(10,10);
            new_grn(ew,ns,:) = new_grn(10,10,:);
            new_snalb(ew,ns) = new_snalb(10,10);
            new_lkdep(ew,ns) = new_lkdep(10,10);
        end
    end
end

for ns = 4:6     % ns (north to south)
    for ew = 9:11                   % ew (east to west)
        if flg_orig(ew,ns) == 0 
            new_flg(ew,ns) = new_flg(10,10);
            new_lnd(ew,ns) = new_lnd(10,10);   
            new_frc(ew,ns,:) = new_frc(10,10,:);
            new_alb(ew,ns,:) = new_alb(10,10,:);
            new_slt(ew,ns) = new_slt(10,10);
            new_stp(ew,ns,:) = new_stp(10,10,:);
            new_sct(ew,ns) = new_sct(10,10);
            new_sbt(ew,ns,:) = new_sbt(10,10,:);
            new_scb(ew,ns) = new_scb(10,10);
            new_grn(ew,ns,:) = new_grn(10,10,:);
            new_snalb(ew,ns) = new_snalb(10,10);
            new_lkdep(ew,ns) = new_lkdep(10,10);
        end
    end
end

%% Domain 1

% % Lose Farmington Bay.
% % Make barren land classification (pixel we know is barren land).

for ns = 5:9     % ns (north to south)
    for ew = 7:10                   % ew (east to west)
        if flg_orig(ew,ns) == 0 
            new_flg(ew,ns) = new_flg(10,10);
            new_lnd(ew,ns) = new_lnd(10,10);   
            new_frc(ew,ns,:) = new_frc(10,10,:);
            new_alb(ew,ns,:) = new_alb(10,10,:);
            new_slt(ew,ns) = new_slt(10,10);
            new_stp(ew,ns,:) = new_stp(10,10,:);
            new_sct(ew,ns) = new_sct(10,10);
            new_sbt(ew,ns,:) = new_sbt(10,10,:);
            new_scb(ew,ns) = new_scb(10,10);
            new_grn(ew,ns,:) = new_grn(10,10,:);
            new_snalb(ew,ns) = new_snalb(10,10);
            new_lkdep(ew,ns) = new_lkdep(10,10);
        end
    end
end

for ns = 8:9     % ns (north to south)
    for ew = 1:2                   % ew (east to west)
        if flg_orig(ew,ns) == 0 
            new_flg(ew,ns) = new_flg(10,10);
            new_lnd(ew,ns) = new_lnd(10,10);   
            new_frc(ew,ns,:) = new_frc(10,10,:);
            new_alb(ew,ns,:) = new_alb(10,10,:);
            new_slt(ew,ns) = new_slt(10,10);
            new_stp(ew,ns,:) = new_stp(10,10,:);
            new_sct(ew,ns) = new_sct(10,10);
            new_sbt(ew,ns,:) = new_sbt(10,10,:);
            new_scb(ew,ns) = new_scb(10,10);
            new_grn(ew,ns,:) = new_grn(10,10,:);
            new_snalb(ew,ns) = new_snalb(10,10);
            new_lkdep(ew,ns) = new_lkdep(10,10);
        end
    end
end

%% Plot some stuff to check how things went.

figure(2)
subplot(1,2,1)
title('Landmask Before Modification')
ax=usamap([lat(x_min,y_min),lat(x_max,y_max)],[lon(x_min,y_max),lon(x_max,y_min)])
surfm(lat(x_min:x_max,y_min:y_max),lon(x_min:x_max,y_min:y_max),flg(x_min:x_max,y_min:y_max))
subplot(1,2,2)
title('Landmask After Modification')
ax=usamap([lat(x_min,y_min),lat(x_max,y_max)],[lon(x_min,y_max),lon(x_max,y_min)])
surfm(lat(x_min:x_max,y_min:y_max),lon(x_min:x_max,y_min:y_max),new_flg)

figure(3)
subplot(1,2,1)
title('NLCD 2011 Before Modification')
ax=usamap([lat(x_min,y_min),lat(x_max,y_max)],[lon(x_min,y_max),lon(x_max,y_min)])
surfm(lat(x_min:x_max,y_min:y_max),lon(x_min:x_max,y_min:y_max),lnd(x_min:x_max,y_min:y_max))
subplot(1,2,2)
title('NLCD 2011 After Modification')
ax=usamap([lat(x_min,y_min),lat(x_max,y_max)],[lon(x_min,y_max),lon(x_max,y_min)])
surfm(lat(x_min:x_max,y_min:y_max),lon(x_min:x_max,y_min:y_max),new_lnd)

figure(4)
pcolor(new_lnd')
title('Dominant Land Use After Modification')

%% Implement changes into entire domain.

fin_flg = flg;
fin_flg(x_min:x_max,y_min:y_max) = new_flg; % Implement changes.
fin_lnd = lnd;
fin_lnd(x_min:x_max,y_min:y_max) = new_lnd; 
fin_frc = frc;
fin_frc(x_min:x_max,y_min:y_max,:) = new_frc;
fin_alb = alb;
fin_alb(x_min:x_max,y_min:y_max,:) = new_alb;
fin_slt = slt;
fin_slt(x_min:x_max,y_min:y_max) = new_slt;
fin_stp = stp;
fin_stp(x_min:x_max,y_min:y_max,:) = new_stp;
fin_sct =sct;
fin_sct(x_min:x_max,y_min:y_max) = new_sct;
fin_sbt = sbt;
fin_sbt(x_min:x_max,y_min:y_max,:) = new_sbt;
fin_scb = scb;
fin_scb(x_min:x_max,y_min:y_max) = new_scb;
fin_grn = grn;
fin_grn(x_min:x_max,y_min:y_max,:) = new_grn;
fin_snalb = snalb;
fin_snalb(x_min:x_max,y_min:y_max,:) = new_snalb;
fin_lkdep = lkdep;
fin_lkdep(x_min:x_max,y_min:y_max) = new_lkdep;

% Modify existing met_em files with new land use.

ncwrite(strcat(run_path,'WPS/geo_em.',domain,'.nc'),'LU_INDEX',fin_lnd);
ncwrite(strcat(run_path,'WPS/geo_em.',domain,'.nc'),'LANDUSEF',fin_frc);
ncwrite(strcat(run_path,'WPS/geo_em.',domain,'.nc'),'LANDMASK',fin_flg);
ncwrite(strcat(run_path,'WPS/geo_em.',domain,'.nc'),'ALBEDO12M',fin_alb);
ncwrite(strcat(run_path,'WPS/geo_em.',domain,'.nc'),'SOILTEMP',fin_slt);
ncwrite(strcat(run_path,'WPS/geo_em.',domain,'.nc'),'SOILCTOP',fin_stp);
ncwrite(strcat(run_path,'WPS/geo_em.',domain,'.nc'),'SCT_DOM',fin_sct);
ncwrite(strcat(run_path,'WPS/geo_em.',domain,'.nc'),'SOILCBOT',fin_sbt);
ncwrite(strcat(run_path,'WPS/geo_em.',domain,'.nc'),'SCB_DOM',fin_scb);
ncwrite(strcat(run_path,'WPS/geo_em.',domain,'.nc'),'GREENFRAC',fin_grn);
ncwrite(strcat(run_path,'WPS/geo_em.',domain,'.nc'),'SNOALB',fin_snalb);
ncwrite(strcat(run_path,'WPS/geo_em.',domain,'.nc'),'LAKE_DEPTH',fin_lkdep);