%This is working example for the MUSIVA MetOp/IASI retrieval simulator. It has
%been written by Matthias Schneider on 19.1.2017. It requires a netcdf File
%(exampledata.nc) with the model exampe data and an apriori data File
%(apri.dat) with the H2O and delD apriori data used for the retrieval.

clear variables;

%****************************************************************
%* Reading of Example Data files (all states of whole ensemble) *
%****************************************************************
nc_infile='example_ECHAM5wiso_20140212hh00.nc';
%nc_infile='example_ECHAM5wiso_20140812hh00.nc';
lon=ncread(nc_infile,'lon'); % dim: L1x1
lat=ncread(nc_infile,'lat'); % dim: L2x1 
alt=ncread(nc_infile,'height'); % dim: L1xL2xA1x1; unit m
H2O=ncread(nc_infile,'WVMR'); % dim: L1xL2xA1x1; unit ppmv
delD=ncread(nc_infile,'deltaD'); % dim: L1xL2xA1x1; unit permil 
temp=ncread(nc_infile,'temperature'); % dim: L1xL2xA1x1; unit K
press=ncread(nc_infile,'pressure'); % dim: L1xL2xA1x1; unit hPa 
srf_skinT=ncread(nc_infile,'skin_temp'); % dim: L1xL2x1; unit K
srf_emi=ncread(nc_infile,'emissivity');  % dim: L1xL2
cld=ncread(nc_infile,'cloudflag');  % dim: L1xL2x1
%Dimensions for example file: L1=320 (longitude); L2=160 (latitude); A1=31 (altitude)
L1=size(lon,1);
L2=size(lat,1);
%****************************************************************

%***********************************
%* Reading of H2O and delD apriori *
%***********************************
fidAPR=fopen('apri.dat');
apr_i=0;
lineAPR = fgetl(fidAPR);
while 1,
   lineAPR = fgetl(fidAPR);
   if ~isstr(lineAPR), break, end
   A = sscanf(lineAPR, '%f');
   apr_i=apr_i+1;
   apralt(apr_i)=A(1);
   aprH2O(apr_i)=A(2);
   aprdelD(apr_i)=A(3);
end
fclose(fidAPR);
%***********************************

%*******************
%* Working Example *
%*******************
r=0;
for i=1:+1:L1,
    for j=1:+1:L2,
        if floor(alt(i,j,1)/10)*10<=5000, %ground height <=5000m a.s.l.
        %if cld(i,j,1)==0 && floor(alt(i,j,1)/10)*10<=5000, %no cloud and ground height <=5000m a.s.l.
            r=r+1;
            %***Fixing of IASI gridding
            A1=size(alt,3);
            nol=28; 
            alt_lev(1)=floor(alt(i,j,1)/10)*10;
            if alt_lev(1)<0, alt_lev(1)=0; end
            if alt_lev(1)>=390.0, nol=27; end 
            if alt_lev(1)>=830.0, nol=26; end 
            if alt_lev(1)>=1300.0, nol=25; end 
            if alt_lev(1)>=1820.0, nol=24; end 
            if alt_lev(1)>=2370.0, nol=23; end 
            if alt_lev(1)>=2950.0, nol=22; end 
            if alt_lev(1)>=3570.0, nol=21; end 
            if alt_lev(1)>=4220.0, nol=20; end 
            if alt_lev(1)>=4900.0, nol=19; end 
            if alt_lev(1)>=5620.0, nol=18; end 
            alt_lev(2:nol)=apralt(2+(28-nol):28);
            aprH2O_lev(1)=exp(interp1(apralt,log(aprH2O),alt_lev(1),'linear','extrap'));              
            aprH2O_lev(2:nol)=aprH2O(2+(28-nol):28);
            aprdelD_lev(1)=interp1(apralt,aprdelD,alt_lev(1),'linear','extrap');              
            aprdelD_lev(2:nol)=aprdelD(2+(28-nol):28);
                
            %***Individual states
            altij(1:A1)=alt(i,j,1:A1,1);
            pressij(1:A1)=press(i,j,1:A1,1);
            H2Oij(1:A1)=H2O(i,j,1:A1,1);
            delDij(1:A1)=delD(i,j,1:A1,1);
            tempij(1:A1)=temp(i,j,1:A1,1);
            %***Interpolation on retrieval levels
            for lev=1:+1:nol,
                press_lev(lev)=exp(interp1(altij(1:A1),log(pressij(1:A1)),alt_lev(lev),'linear','extrap'));
                H2O_lev(lev)=exp(interp1(altij(1:A1),log(H2Oij(1:A1)),alt_lev(lev),'linear','extrap'));
                delD_lev(lev)=interp1(altij(1:A1),delDij(1:A1),alt_lev(lev),'linear','extrap');
                temp_lev(lev)=interp1(altij(1:A1),tempij(1:A1),alt_lev(lev),'linear','extrap');
                %apriori proxy state
                apr_pstate(lev)=log(aprH2O_lev(lev))+log(aprdelD_lev(lev)+1)/2;
                apr_pstate(lev+nol)=log(aprdelD_lev(lev)+1);
                %proxy state
                pstate(lev)=log(H2O_lev(lev))+log(delD_lev(lev)+1)/2;
                pstate(lev+nol)=log(delD_lev(lev)+1);
            end
            %***CALL FUNCTION FOR AVK SIMULATION 
            [kern_theo_T1(1:2*nol,1:2*nol),kern_theo_T2(1:2*nol,1:2*nol)]=AVKsimulator(25,srf_emi(i,j),srf_skinT(i,j,1),press_lev(1:nol),temp_lev(1:nol),log(H2O_lev(1:nol)),alt_lev(1:nol),nol); %call function
            %*************

            %***kernel multiplication
            sm_pstate(1:2*nol)=kern_theo_T2(1:2*nol,1:2*nol)*(pstate(1:2*nol)-apr_pstate(1:2*nol))'+apr_pstate(1:2*nol)'; %'smoothed' proxy state
            for lev=1:+1:nol,
                sm_H2O_lev(lev)=exp(sm_pstate(lev)-sm_pstate(lev+nol)/2);
                sm_delD_lev(lev)=(exp(sm_pstate(lev+nol))-1);
            end 
            
            %***quality check
            IDM(1:2*nol,1:2*nol)=0;
            SadD(1:nol,1:nol)=0;
            corrlength=5000;
            corrlengthbl=500; %weak coupling with boundary layer
            for levi=1:+1:nol,
                IDM(levi,levi)=1;
                adD(levi)=0.1;
            end
            for levi=1:+1:nol,
                for levj=1:+1:nol,
                    SadD(levi,levj)=adD(levi)*adD(levj)*exp(-(((alt_lev(levi)-alt_lev(levj))^2)/(2*corrlength^2)));
                end
            end
            for levi=1:+1:3,  %weak coupling between free troposphere and boundary layer
                for levj=4:+1:nol,
                    SadD(levi,levj)=adD(levi)*adD(levj)*exp(-(((alt_lev(levi)-alt_lev(levj))^2)/(2*corrlengthbl^2)));
                    SadD(levj,levi)=adD(levi)*adD(levj)*exp(-(((alt_lev(levi)-alt_lev(levj))^2)/(2*corrlengthbl^2)));
                end
            end
            Serr(1:nol,1:nol)=(IDM(1:nol,1:nol)-kern_theo_T2(nol+1:2*nol,nol+1:2*nol))*SadD(1:nol,1:nol)*(IDM(1:nol,1:nol)-kern_theo_T2(nol+1:2*nol,nol+1:2*nol))';
            DOFS_T1=trace(kern_theo_T1(1:nol,1:nol)); %DOFS product Type 1
            DOFS_T2=trace(kern_theo_T2(1:nol,1:nol)); %DOFS product Type 2
            for lev=1:+1:nol,
                sensErr(lev)=sqrt(Serr(lev,lev));
            end
            [sensErrmin,index]=min(sensErr); %best detectable level (level with minimum of sensErr)
            
            %***output  
            %The example_output.dat is an ascii-File that contains
            %information from the exampledata.nc file (lon, lat, ground
            %altitude, skin temperature) as well as results from the
            %retrieval simulator: DOFS, sensitivity estimations, original 
            %model data and model data that has passed through the retrieval 
            %simulator for the altitude with the estimated best sensitivity
            %and for the altitude that is closest to 5km (altitude for which 
            %on average the best sensitivity is achieved).
            if r==1,
               fod=fopen('example_output.dat','wt');
               fprintf(fod,'lon lat GNDalt SKINtemp cld ');
               fprintf(fod,'DOFS_T1 DOFS_T2 ');
               fprintf(fod,'altbest Sensbest ');
               fprintf(fod,'MH2Obest MdelDbest ');
               fprintf(fod,'smMH2Obest smMdelDbest ');
               fprintf(fod,'alt5km Sens5km ');
               fprintf(fod,'MH2O5km MdelD5km ');
               fprintf(fod,'smMH2O5km smMdelD5km ');
               fprintf(fod,'\n');
               fclose(fod);
            end
            fod=fopen('example_output.dat','at');
            fprintf(fod,'%.3f %.3f %.3f %.2f %i ', lon(i,1), lat(j,1), alt_lev(1), srf_skinT(i,j,1), cld(i,j,1));
            fprintf(fod,'%.3f %.3f ', DOFS_T1, DOFS_T2); 
            fprintf(fod,'%.3f %.3f ', alt_lev(index), sensErrmin);  %output for level with best sensitivity
            fprintf(fod,'%.3f %.3f ', H2O_lev(index), delD_lev(index));
            fprintf(fod,'%.3f %.3f ', sm_H2O_lev(index), sm_delD_lev(index));
            fprintf(fod,'%.3f %.3f ', alt_lev((10-(28-nol))), sensErr(10-(28-nol)));  %output for level (10-(28-nol)), which is at 4900m a.s.l.; it is level 19 starting with uppermost level
            fprintf(fod,'%.3f %.3f ', H2O_lev(10-(28-nol)), delD_lev(10-(28-nol)));
            fprintf(fod,'%.3f %.3f ', sm_H2O_lev(10-(28-nol)), sm_delD_lev(10-(28-nol)));
            fprintf(fod,'\n');
            fclose(fod);
        end
    end
end
%*******************

