function[] = femao2lake()

disp('Starting the script of conversion from FEMAO output to the LAKE input')

Kelvin0 = 273.15;
mmm2msec = 1e-3/(30*24*60*60); % mm/month to m/s
relhum = 0.9;
expname = 'essea1948-2011';

numpoint = 2; %number of point
M = 40; % number of layers in LAKE
%depth = 20; % the depth (m) in the point
width = 1e+4; 

% point 1
phi(1) = 71;
lambda(1) = 165;
depth(1) = 20;
% point 2
phi(2) = 79; 
lambda(2) = 165;
depth(2) = 2100;
%point 3
phi(3) = 85; 
lambda(3) = 165;
depth(3) = 2100; % check the value

dir = 'd:\Models\LAKE\2013\March\lake\data\essea\';
%dir = '/home/victor/Files/main/models/LAKE/2013/March/data/essea/';
% Atmospheric variables
ta = importdata(strcat(dir,'ta_methane',num2str(numpoint),'.dat'));
pa = importdata(strcat(dir,'pa_methane',num2str(numpoint),'.dat'));
cloud = importdata(strcat(dir,'cloud_methane',num2str(numpoint),'.dat'));
wind = importdata(strcat(dir,'wind_methane',num2str(numpoint),'.dat'));
precip = importdata(strcat(dir,'precip_methane',num2str(numpoint),'.dat'));

% Ocean variables
zin = [ 0 2.5 5 7.5 10 15 20 27.5 37.5 50 65 82.5 100 125 ...
    150 175 200 225 250 275 300 325 350 375 410 450 500 600 700 850 ...
    1050 1250 1500 1750 2100 2500 3000 3500 4000 4500]';
ncols = size(zin,1) + 1;
format = '%f ';
for i = 2:ncols
    format = [format '%f '];
end

% Water temperature
fid = fopen(strcat(dir,'T3Dmodel',num2str(numpoint),'.dat'),'r');
temp = fscanf(fid,'%f');
nraws = numel(temp)/ncols;
temp = reshape(temp,ncols,nraws);
temp = temp';

% Salinity
fid = fopen(strcat(dir,'S3Dmodel',num2str(numpoint),'.dat'),'r');
sal = fscanf(fid,'%f');
sal = reshape(sal,ncols,nraws);
sal = sal'*1e-3; %converting to promille

% u-speed
fid = fopen(strcat(dir,'u3Dmodel',num2str(numpoint),'.dat'),'r');
u = fscanf(fid,'%f');
u = reshape(u,ncols,nraws);
u = u'*1e-2; %converting to m/s

% v-speed
fid = fopen(strcat(dir,'v3Dmodel',num2str(numpoint),'.dat'),'r');
v = fscanf(fid,'%f');
v = reshape(v,ncols,nraws);
v = v'*1e-2;  %converting to m/s

dzLAKE = depth(numpoint)/M;
zLAKE = [0:dzLAKE:depth(numpoint)];

% Interpolation at LAKE grid
work = zeros(M+1,7); % 7 -number of columns in "inflow file"
work1 = zeros(M+1,3); % 3 - number of columns in "outflow file"
for i = 1:nraws
    tempint = interp1(zin,temp(i,2:ncols),zLAKE);
    salint  = interp1(zin,sal(i,2:ncols),zLAKE);
    uint = interp1(zin,u(i,2:ncols),zLAKE);
    vint = interp1(zin,v(i,2:ncols),zLAKE);    
    work2(1:M+1) = width;    
    work = [[1:1:M+1]' work2' [(uint(:).^2 + vint(:).^2).^(0.5)] tempint' salint' uint' vint']; 
    work1 = [[1:1:M+1]' work2' [(uint(:).^2 + vint(:).^2).^(0.5)]];
    if i == 1
        dlmwrite(strcat(dir,expname,'infl_',num2str(numpoint),'.dat'), work, 'delimiter', ' ')
        dlmwrite(strcat(dir,expname,'outfl_',num2str(numpoint),'.dat'), work1, 'delimiter', ' ')
    end
    dlmwrite(strcat(dir,expname,'infl_',num2str(numpoint),'.dat'), work, '-append', 'delimiter', ' ')
    dlmwrite(strcat(dir,expname,'outfl_',num2str(numpoint),'.dat'), work1, '-append', 'delimiter', ' ')
end
clear work work1 work2

nyear_begin = 1948;
nyear_end = 2011;

ncols = 10;
N_Uspeed    = 9;
N_Vspeed    = 10;
N_Temp       = 4;
N_Hum        = 6;
N_Pres        = 5;
N_SWdown  = 3;
N_Precip     =  7;
N_Cloud       =  8;

[m,~]=size(ta);
work = zeros(m,ncols);

n = 0;
for i = nyear_begin:nyear_end
    nm = 1; % month number
    for j = 1:365
        n = n + 1;
        work(n,N_Uspeed) = wind(n,3);
        work(n,N_Vspeed) = wind(n,4);
        work(n,N_Temp) = ta(n,3) + Kelvin0;
        work(n,N_Pres) = pa(n,3);
        work(n,N_Hum) = 0.622*relhum*satvap(ta(n,3))/pa(n,3);
        if (mod(j,30) == 0 && j ~= 360)
            nm = nm + 1;          
        end
        work(n,N_Precip) = precip(nm,2)*mmm2msec;
        work(n,N_Cloud) = cloud(nm,2);
        work(n,N_SWdown) = solar(j,phi(numpoint),lambda(numpoint),relhum*satvap(ta(n,3)),cloud(nm,2));
    end    
end

dlmwrite(strcat(dir,expname,'_',num2str(numpoint),'.dat'), work, ' ')

disp('The script finished successfully')

end

function[SW] = solar(m1day,fi,Xout,evapor,cloud)

Solar = 1367; %solar constant
%mhour = 12; % For daily timestep

% For daily timestep
SW = 0;
for mhour = 0:23

% Solar time at the location
stime=(mhour - 0.5); % mhour - GMT time
stime=stime +Xout/15.;

day = m1day + stime/24.;
if (stime > 24.) 
    stime = stime - 24.;
end

% Declination angle
delta= 23.44*cos((172.-day)*pi/180.);
sd=sin(delta*pi/180.);
cd=cos(delta*pi/180.);

CZ= max([0. sin(fi*pi/180.)*sd+cos(fi*pi/180.)*cd*cos((12.-stime)*pi/12.)]);
Q0=Solar*CZ^2/((CZ+2.7)*evapor*1.e-5 +1.085*CZ+0.10);

% Clouds
SW = SW + Q0*(1.-0.6*cloud^3);

end
SW = SW/24; % averaging over the diurnal cycle

end

function[ea]=satvap(TC)

% The -12C limit is by E.V. Volodin. AMIP. See QMAX function.
if (TC > 0.) % Warm atmosphere.
    % In AOMIP - only warm atmosphere.
    % ea= 0.9* 611.*(10.**(7.5*TC/(273.16-35.86+TC)))
    
    ea = 10.^(2.+(0.7859+0.03477*TC)/(1.+0.00412*TC)) ;
 
else % Atmosphere with snow and ice. 
    ea = 10.^(0.00422*TC+2.+(0.7859+0.03477*TC)/(1.+0.00412*TC)) ;
    % ea= 0.9* 611.*10.**(9.5*TC/(273.16- 7.66+TC))
end % TA > 0C

end