% WB_2021_.m (underscore) works with WBloop_2021_.m % 14th January 2021
% revised 30 Dec 2020 to reassign WsD as White Flag WBGT and have Yellow YsD
% also remove condition of Durbin-Watson to report the trend
% also base enthalpy
% v17B revised 22nd March 2017
% be sure to add to function % close(gcf); % close current figure
% WB17.m webbulb 17th edition 30th January 2017 from local disk copy of ftp://ftp.ncdc.noaa.gov/pub/data/gsod/ ;
% WB17local.n 3th February set unix=1 for run on euramoo
% Be sure to include function psych.m 
% [wbC, H2O] = psych(dbC, dpC, stPa)

% and function WBloop.m
% savings = WBloop(WMO)
%% indirectly control WMO
% Enter country index IDEA
versionstem=num2str(2021);
suffix='_';
% vector{...
switch suffix
    case 'a'
        framingvector=1:3;
    case 'b'
        framingvector=4:6;
    case 'c'
        framingvector=7:9;
    case 'd'
        framingvector=10:13;
    case 'e'
        framingvector=14:17;
    case 'f'
        framingvector=18:21;
    case 'g'
        framingvector=22:25;
    case 'h'
        framingvector=26:30;
    otherwise
        framingvector=1:30; % everything
end
% for IDEA=1:5 % WB2021a  
% for IDEA=6:10; % WB2021b   
% for IDEA=11:15; % WB2021c  
% for IDEA=16:20; % WB2021d  
% for IDEA=21:25; % WB2021e  
% for IDEA=26:30; % WB2021f  
for IDEA=framingvector
    beenthere=0; 
startyear=1987; % up to 34 years
%%
unix=1; % switch to 1 for unix, zero for windows
% unix=1; % for unix
%%
average_fname=strcat('RESULTSout_',versionstem,suffix,'.txt');
rates_fname=strcat('RESULTSrates_',versionstem,suffix,'.txt');
stats_fname=strcat('RESULTSstats_',versionstem,suffix,'.txt');
if unix==1
	% keypath='/home/xxepeter/WB/NOAAdata/'; % running unix
    keypath='/home/home01/fbsepet/nb/WB/'; % running unix on arc 3
else
    % keypath='C:\Users\xxepeter\Documents\MATLAB\NOAAdata\'; % running windows
    % keypath='C:\Users\eric\Documents\MATLAB\WB\'; % 'D:\WB\'; % running windows
    keypath='D:\WB\'; % 'D:\WB\'; % running windows
end
summary=strcat(keypath,average_fname); % summary output file
rates=strcat(keypath,rates_fname); % rates output file
stats=strcat(keypath,stats_fname); % statistics output file
if beenthere==0
    % txt = [date 'Appending\n'];
    fid = fopen(summary,'at');
    fprintf(fid,'WMOstn WBANo Acz HVAC meanPa WBref hBASE dpy [ avHDD  avCDD  avChD ] low min ave max best [ WsD YsD ] yrs dys [ Icy Hea Hum Pas Vtn EvC  DX ] 24Cd 67Fd  latitude longitude    elev base stop CC ST Placename\n');
    %            WMOstn WBANo Acz HVAC meanPa WBref hBASE dpy [ avHDD  avCDD  avChD ] low min ave max best [ WsD YsD ] yrs dys [ Icy Hea Hum Pas Vtn EvC  DX ] 24Cd 67Fd  latitude longitude    elev base stop CC ST Placename
    %              8260 99999   3 type 101325  19.4  55.4 [  1709     11     80 ]   0   3  15  25  31  [   8   0 ]   4 274 [  25 170 NaN  62  16   0   1 ]    9   67    0.0000    0.0000       0 2006 2009       WXPOD8270
    %             10010 99999   7 type 101052  19.4  55.5 [  6544      0      0 ]   0 NaN NaN NaN NaN  [   0   0 ]  32 364 [ 282  82 NaN   0   0   0   0 ]    0   67   70.9330   -8.6670       9 1987 2019 NO    JAN MAYEN(NOR-NAVY)
    %             10013 99999   5 type    NaN  19.4   NaN [  4038      0      0 ]   0 NaN NaN NaN NaN  [   0   0 ]   1 302 [ 109 192 NaN   1   0   0   0 ]    0   67       NaN       NaN     NaN 1987 1987 NO    ROST
    fclose(fid);

    fid = fopen(stats,'at');
    fprintf(fid,'WMOstn WBANo [ adjR2_H ob_H pVal_H DuWa_H pVDWH ][ adjR2_C obC pVal_C DuWaC pVDWC ][ adjR2_E ob_E   pVal_E DuWaE pVDWE ][ adjR2_D obD pValD DuWaD pVDWD ][ adjR2_W obW pValW DuWaW pVDWW ][ adjR2_Y obY pValY DuWaY pVDWY ] base stop Placename\n');
    %            WMOstn WBANo [ adjR2_H ob_H pVal_H DuWa_H pVDWH ][ adjR2_C obC pVal_C DuWaC pVDWC ][ adjR2_E ob_E   pVal_E DuWaE pVDWE ][ adjR2_D obD pValD DuWaD pVDWD ][ adjR2_W obW pValW DuWaW pVDWW ][ adjR2_Y obY pValY DuWaY pVDWY ] base stop Placename
    %              8260 99999 [  0.9938    4  0.002 0.9370 0.350 ][  0.3908   4 0.2290 1.648 0.978 ][  0.8920    4    0.037 2.060 0.710 ][ -0.3342   4 0.668 1.831 0.887 ][  0.8561   4 0.049 1.505 0.873 ][ -0.3342   4 0.668 1.831 0.887 ] 2006 2009 WXPOD8270
    %             10010 99999 [  0.6522   33  0.000 1.5550 0.133 ][     NaN  33    NaN   NaN 0.000 ][     NaN   32      NaN   NaN 0.000 ][     NaN  32   NaN   NaN 0.000 ][     NaN  32   NaN   NaN 0.000 ][     NaN  32   NaN   NaN 0.000 ] 1987 2019 JAN MAYEN(NOR-NAVY)
    %             10013 99999 [     NaN    1    NaN    NaN 0.000 ][     NaN   1    NaN   NaN 0.000 ][     NaN    1      NaN   NaN 0.000 ][     NaN   1   NaN   NaN 0.000 ][     NaN   1   NaN   NaN 0.000 ][     NaN   1   NaN   NaN 0.000 ] 1987 1987 ROST
    
    fid = fopen(rates,'at');
    fprintf(fid,'WMOstn WBANo [ HdD18 std_H  slope_H SEE_H HDD18 ][ CdD24 stdC  slope_C SEEC CDD24 ][ EnthDys stdE   slope_E SEEE EnthD ][ DXd std_D  slope_D SEE_D  DXD ][ Wsd std_W  slope_W SEE_W  WsD ][ Ysd std_Y  slope_Y SEE_Y  YsD ] base stop Placename\n');
    %            WMOstn WBANo [ HdD18 std_H  slope_H SEE_H HDD18 ][ CdD24 stdC  slope_C SEEC CDD24 ][ EnthDys stdE   slope_E SEEE EnthD ][ DXd std_D  slope_D SEE_D  DXD ][ Wsd std_W  slope_W SEE_W  WsD ][ Ysd std_Y  slope_Y SEE_Y  YsD ] base stop Placename
    %              8260 99999 [  1709  1142     1.14    90  2851 ][    11    9    0.010    7    20 ][   80.00   55     0.050   18   135 ][   1   1.0    0.000   1.0    2 ][   8     5     0.00     2   13 ][   0     1     0.00     1    1 ] 2006 2009 WXPOD8270
    %             10010 99999 [  6544   386   -32.50   228  6252 ][     0    0   -0.000    0     0 ][    0.00    0    -0.000    0     0 ][   0   0.0   -0.000   0.0    0 ][   0     0    -0.00     0    0 ][   0     0    -0.00     0    0 ] 1987 2019 JAN MAYEN(NOR-NAVY)
    %             10013 99999 [  4038     0     2.03     0  4038 ][     0    0    0.000    0     0 ][    0.00    0     0.000    0     0 ][   0   0.0    0.000   0.0    0 ][   0     0     0.00     0    0 ][   0     0     0.00     0    0 ] 1987 1987 ROST
    fclose(fid);  
    beenthere=0;
end
listfile='isd-history_2021.csv';
filekey = strcat(keypath,listfile);
delimiter = ',';
startRow = 2;
formatSpec3 = '%f%f%s%s%s%s%f%f%f%f%f%[^\n\r]';
fileID = fopen(filekey,'r');
dataArray0 = textscan(fileID, formatSpec3, 'Delimiter', delimiter, 'EmptyValue' ,NaN,'HeaderLines' ,startRow-1, 'ReturnOnError', false);
fclose(fileID);

USAF = dataArray0{:, 1}; % WMO
WBAN = dataArray0{:, 2}; % WBAN
STATIONNAME = dataArray0{:, 3}; % station name
LAT = dataArray0{:, 7};
LON = dataArray0{:, 8};
ELEVM=dataArray0{:, 9};
CTRY = dataArray0{:, 4};
STATE = dataArray0{:, 5};

vector{1} = [1:1000];
vector{2} = [1001:2000];
vector{3} = [2001:3000];
vector{4} = [3001:4000];
vector{5} = [4001:5000];
vector{6} = [5001:6000];
vector{7} = [6001:7000];
vector{8} = [7001:8000];
vector{9} = [8001:9000];
vector{10} = [9001:10000];
vector{11} = [10001:11000];
vector{12} = [11001:12000];
vector{13} = [12001:13000];
vector{14} = [13001:14000];
vector{15} = [14001:15000];
vector{16} = [15001:16000];
vector{17} = [16001:17000];
vector{18} = [17001:18000];
vector{19} = [18001:19000];
vector{20} = [19001:20000];
vector{21} = [20001:21000];
vector{22} = [21001:22000];
vector{23} = [22001:23000];
vector{24} = [23001:24000];
vector{25} = [24001:25000];
vector{26} = [25001:26000];
vector{27} = [26001:27000];
vector{28} = [27001:28000];
vector{29} = [28001:29000];
vector{30} = [29001:29555];
                                

%%	

Victor=vector{IDEA};
for j = Victor  % j=40 is starting at WMO=10000
try
    %%
    WMO = USAF(j); % work down the list based on index j
    WBANo = WBAN(j);
    country=CTRY{j};
    state=STATE{j};
    place=STATIONNAME{j};
    latitude=LAT(j); 
    longitude=LON(j); 
    elev=ELEVM(j);
    skipthissite=0;
        if isnan(elev)
            elev=0; % as good a guess as any
        end
    % savings = WBlooplocal(WMO)
    % savings = WBloop16local(WMO)
    [ASHRAEcz, timeframe, HDDstats, CDDstats, ChDstats, WBGTstatsW, WBGTstatsY, DXstats, PC, savings, Basis] = WBloop_2021_(WMO,WBANo,startyear); 
    % Basis [mean(stPa) Thumid BASEenthalpy]: Given mean station pressure (Pascals) and Thumid (wb °C), Basis of cooling enthalpy (kJ/kg-da)
    Nobs=int32(   min([  HDDstats(6), CDDstats(6), ChDstats(6), WBGTstatsW(6), WBGTstatsY(6), DXstats(6)  ])   ); % observations (minimum)
    % note there are 7 statistics in HDDstats, CDDstats, ChDstats, WBGTstats, and DXstats (1:9) [mean stdev slope SEE R2 n=obs design pV pVautocor]
    % and 9 columns in PC (1:9) [total icing heating humid passive ventn evapC DX allC], 
    % and 5 columns in savings (1:5) [least min mean max best]
    %%
    HVAC='????'; % be sure to have exactly four characters
    winter='__'; summer='__'; % assume nothing
    floatdays=10; % nominally 10 days per year
    if PC(2)>=floatdays % Too Icy in winter for heat pump
    	winter='Gs'; % gas or groundsource winter
    else
        if (PC(2)+PC(3))>=floatdays % airsource heat pump will work
            winter='Hp'; % heat pump works
        else
            if (PC(2)+PC(3)+PC(4))>=floatdays % humidification
                winter='Hu'; % heat pump could serve
            end
        end
    end
    if PC(8)>=floatdays % DXX
        summer='Ac'; % dx air-conditioning for cooling
    else
        if (PC(7)+PC(8))>=floatdays % evap cooler
            summer='Ev'; % evaporative cool summer
        else
            if (PC(6)+PC(7)+PC(8))>=floatdays % ventilation
                summer='Vt';
            end
        end
    end
    HVAC=strcat(winter,summer);
    %%
    fid = fopen(summary,'at'); % summary output file
    fprintf(fid,'%6.0f %5.0f %3.0f %s %6.0f %5.1f %5.1f %3.0f [%6.0f %6.0f %6.0f ] %3.0f %3.0f %3.0f %3.0f %3.0f  [ %3.0f %3.0f ] %3.0f %3.0f [ %3.0f %3.0f %3.0f %3.0f %3.0f %3.0f %3.0f ] %4.0f %4.0f %9.4f %9.4f %7.0f %4.0f %4.0f %2s %2s %s\n' ...
        , WMO,  WBANo, ASHRAEcz, HVAC, Basis, HDDstats(1), CDDstats(1), ChDstats(1), savings, WBGTstatsW(1), WBGTstatsY(1), Nobs, PC, latitude, longitude, elev, timeframe, country, state, place);  
    % WMOstn WBANo Acz HVAC meanPa WBref hBASE dpy [ avHDD  avCDD  avChD ] low min ave max best [ WsD YsD ] yrs dys [ Icy Hea Hum Pas Vtn EvC  DX ] 24Cd 67Fd latitude longitude    elev base stop CC ST Placename
    %   8260 99999   3 type 101325  19.4  55.4     [  1709     11     80 ]   0   3  15  25  31  [   8   0 ]   4 274 [  25 170 NaN  62  16   0   1 ]    9   67  0.0000    0.0000       0 2006 2009       WXPOD8270
    %  10010 99999   7 type 101052  19.4  55.5     [  6544      0      0 ]   0 NaN NaN NaN NaN  [   0   0 ]  32 364 [ 282  82 NaN   0   0   0   0 ]    0   67 70.9330   -8.6670       9 1987 2019 NO    JAN MAYEN(NOR-NAVY)
    %  10013 99999   5 type    NaN  19.4   NaN     [  4038      0      0 ]   0 NaN NaN NaN NaN  [   0   0 ]   1 302 [ 109 192 NaN   1   0   0   0 ]    0   67     NaN       NaN     NaN 1987 1987 NO    ROST
    fclose(fid);
    
    fid = fopen(rates,'at'); % rates output file
    Vselect=[1 2 3 4 7]; % calling for the 1st, 2nd, 3rd, 4th, 5th and 7th statistics
    fprintf(fid,'%6.0f %5.0f [%6.0f %5.0f %8.2f %5.0f %5.0f ][%6.0f %4.0f %8.3f %4.0f %5.0f ][%8.0f %4.0f %9.3f %4.0f %5.0f ][%4.0f %5.1f %8.3f %5.1f %4.0f ][%4.0f %5.1f %8.3f %5.1f %4.0f ][%4.0f %5.1f %8.3f %5.1f %4.0f ] %4.0f %4.0f %s\n' ...
        , WMO,  WBANo,  HDDstats(Vselect), CDDstats(Vselect), ChDstats(Vselect), DXstats(Vselect), WBGTstatsW(Vselect), WBGTstatsY(Vselect), timeframe, place);
    % WMOstn WBANo [ HdD18 std_H  slope_H SEE_H HDD18 ][ CdD24 stdC  slope_C SEEC CDD24 ][ EnthDys stdE   slope_E SEEE EnthD ][ DXd std_D  slope_D SEE_D  DXD ][ Wsd std_W  slope_W SEE_W  WsD ][ Ysd std_Y  slope_Y SEE_Y  YsD ] base stop Placename
    %   8260 99999 [  1709  1142     1.14    90  2851 ][    11    9    0.010    7    20 ][      80   55     0.050   18   135 ][   1   1.0    0.000   1.0    2 ][   8   5.0    0.000   2.0   13 ][   0   1.0    0.000   1.0    1 ] 2006 2009 WXPOD8270
    %  10010 99999 [  6544   386   -32.50   228  6252 ][     0    0   -0.000    0     0 ][       0    0    -0.000    0     0 ][   0   0.0   -0.000   0.0    0 ][   0   0.0   -0.000   0.0    0 ][   0   0.0   -0.000   0.0    0 ] 1987 2019 JAN MAYEN(NOR-NAVY)
    %  10013 99999 [  4038     0     2.03     0  4038 ][     0    0    0.000    0     0 ][       0    0     0.000    0     0 ][   0   0.0    0.000   0.0    0 ][   0   0.0    0.000   0.0    0 ][   0   0.0    0.000   0.0    0 ] 1987 1987 ROST
    fclose(fid);
    
    fid = fopen(stats,'at'); % statistics output file
    Vselect=[5 6 8 9 10]; % calling for the 3rd and +6/8/9/10th nobs (%2.0f); pV (%5.3f); DuWat (%5.3f); pDW (%5.3f)
    fprintf(fid,'%6.0f %5.0f [%8.4f %4.0f %6.4f %6.3f %5.3f ][%8.4f %3.0f %6.4f %5.3f %5.3f ][%8.4f %4.0f %8.4f %5.3f %5.3f ][%8.4f %3.0f %5.3f %5.3f %5.3f ][%8.4f %3.0f %5.3f %5.3f %5.3f ][%8.4f %3.0f %5.3f %5.3f %5.3f ] %4.0f %4.0f %s\n' ...
        , WMO,  WBANo,  HDDstats(Vselect), CDDstats(Vselect), ChDstats(Vselect), DXstats(Vselect), WBGTstatsW(Vselect), WBGTstatsY(Vselect), timeframe, place);
    % WMOstn WBANo [ adjR2_H ob_H pVal_H DuWa_H pVDWH ][ adjR2_C obC pVal_C DuWaC pVDWC ][ adjR2_E ob_E   pVal_E DuWaE pVDWE ][ adjR2_D obD pValD DuWaD pVDWD ][ adjR2_W obW pValW DuWaW pVDWW ][ adjR2_Y obY pValY DuWaY pVDWY ] base stop Placename
    %   8260 99999 [  0.9938    4 0.0020  0.937 0.350 ][  0.3908   4 0.2290 1.648 0.978 ][  0.8920    4   0.0370 2.060 0.710 ][ -0.3342   4 0.668 1.831 0.887 ][  0.8561   4 0.049 1.505 0.873 ][ -0.3342   4 0.668 1.831 0.887 ] 2006 2009 WXPOD8270
    %  10010 99999 [  0.6522   33 0.0000  1.555 0.133 ][     NaN  33    NaN   NaN 0.000 ][     NaN   32      NaN   NaN 0.000 ][     NaN  32   NaN   NaN 0.000 ][     NaN  32   NaN   NaN 0.000 ][     NaN  32   NaN   NaN 0.000 ] 1987 2019 JAN MAYEN(NOR-NAVY)
    %  10013 99999 [     NaN    1    NaN    NaN 0.000 ][     NaN   1    NaN   NaN 0.000 ][     NaN    1      NaN   NaN 0.000 ][     NaN   1   NaN   NaN 0.000 ][     NaN   1   NaN   NaN 0.000 ][     NaN   1   NaN   NaN 0.000 ] 1987 1987 ROST
    fclose(fid);
        
        sprintf('Results for this station in %s, %s.',STATE{j},CTRY{j})
        %%

catch
    sprintf('There are no results for this station in %s, %s.',STATE{j},CTRY{j})
end
end
% loop of IDEA
end

%{
% lister=[7005 7011 7018 7025 7026 7034 7037 7044 7047 7052 7059 7064 7070 7072];
lister=[10000 10010 10013 10014 10015 10016 10017 10020 10030];
for j=1:length(lister)
    WMO=lister(j);
    savings = WBloop(WMO)
    fid = fopen(average_fname,'at');
    txt = [average_fname ' : ' WMO '\nDate of creation:' savings]; 
    fprintf(fid,'%6.0f %7.2f\n',WMO, savings);
    fclose(fid);
end
%}