clear all
close all

% nolake.m

% Trevor Alcott - Oct 2011

% sample matlab script to modify WRF input files (met_em)
% does not cover soil temp/moisture mods, but will get you started

filepath='/uufs/chpc.utah.edu/common/home/steen_data2/trevor/research_wrf/WPS/';

for domain=3:-1:3
    domstr=['0' num2str(domain)];
    
    namelist=dir([filepath 'met_em.d' domstr '*']);
    
    for f=1:length(namelist)
        
        filename=namelist(f).name;
        file=[filepath filename];
        
        ncid = netcdf.open(file,'NC_WRITE');
        [numdims, numvars, numglobalatts, unlimdimID] = netcdf.inq(ncid);
        for k=0:(numvars-1);
            [varname, xtype, dimids, numatts]=netcdf.inqVar(ncid,k);
            %
            % optional: displays all variables and their IDs
            if f==1 && domain==3
              disp([num2str(k) ' ' varname]);
            end
        end

        lu = uint8(netcdf.getVar(ncid,51));
        landusef = uint8(netcdf.getVar(ncid,52));
        landsea = uint8(netcdf.getVar(ncid,11));
        landmask = uint8(netcdf.getVar(ncid,53));
        lat = double(netcdf.getVar(ncid,74));
        lon = double(netcdf.getVar(ncid,73));
        elev = double(netcdf.getVar(ncid,50));
        
        gsl = find(lat>40.6 & lat<41.7 & lon>-113 & lon<-111.9 & landmask==0);

        % 3-d fields
        temp(:,:,1)=landusef(:,:,16);
        temp(gsl)=0;
        landusef(:,:,16)=temp;
        temp(:,:,1)=landusef(:,:,8);
        temp(gsl)=8;
        landusef(:,:,8)=temp;

        % 2-d fields
        landsea=landsea*0+1;
        lu(gsl)=8;
        landmask(gsl)=1;

        ncwrite(file,'LANDMASK',landmask)
        ncwrite(file,'LANDSEA',landsea)
        ncwrite(file,'LU_INDEX',lu)
        netcdf.close(ncid)
    end
end


