% Iceberg Model -----------------------------------------------------------

%

% NOTE: To run this script you'll need the 5 .mat files that are

% found in this folder:

https://www.dropbox.com/sh/9qc7l3yacyqgn3c/AAC69yLJUH39VsjbEwAdRRMOa?dl=0

%

% This model accompanies the paper "An Analytical Model of Iceberg Drift"

% Wagner, Dell, Eisenman, JPO (2017), henceforth WDE17

%

% Other articles making use of this model:

% Gone with the wind: How model biases skew iceberg meltwater distributions

% Wagner & Eisenman, GRL (2017).

%

% "On the representation of capsizing in iceberg models"

% Wagner, Stern, Dell, Eisenman, Ocean Model (2017)

%

% "Wave inhibition by sea ice enables trans-Atlantic ice rafting of debris

% during Heinrich Events"

% Wagner, Dell, Eisenman, Keeling, Severinghaus, EPSL, (submitted)

%

% The model computes the drift and decay of icebergs, by solving their

% motion analytically for given velocity and SST conditions. Those

% conditions are loaded first (see "input fields").

%

% Then model parameters and analytical expressions are set.

%

% Then run parameters (time, domain, space domain etc) are specified.

%

% Then the release locations are loaded (this is a .mat-file, constructed

% in "construct_seeding.m".

%

% Finally the different bergsizes are loaded and specified.

%

% Then we compute the Lagrangian trajectories, looping over every iceberg

% size class. These loops consist of 2 components for every timestep.

%

% 1) drifting - computes iceberg velocity and translation

% 2) melting  - computes change in iceberg dimensions and volume

%

% Trajectories for every iceberg size class are saved as individual files

% in the output trajectory.

%

% Till Wagner & Ian Eisenman, May 2017,

% tjwagner@ucsd.edu, eisenman@ucsd.edu

 

% -------------------------------------------------------------------------

% pick type of input ------------------------------------------------------

model = 'ECCO2';    %you can either run this with ECCO2 input fields or

% model = 'CCSM4';  %CCSM4 input fields

% -------------------------------------------------------------------------

% load input fields -------------------------------------------------------

% -------------------------------------------------------------------------

% ECCO2 data can be downloaded from

% ftp://ecco2.jpl.nasa.gov/data1/cube/cube92/lat_lon/quart_90S_90N/

% with JRA-25 surface wind reanalysis data

% (https://rda.ucar.edu/datasets/ds625.0/)

% -------------------------------------------------------------------------

% CCSM4 data can be downloaded from

% -------------------------------------------------------------------------

% the demo data set here is for the year 1992, and for the Atlantic only

% here we load surface ocean currents, winds, SST, and landmask from 1 file

load(strcat(model,'_1992_Atlantic.mat'))

fprintf('model data loaded \n')

% -------------------------------------------------------------------------

% Iceberg Parameters and analytical expressions of

% gamma, Lambda, alpha, beta (eqns 7,8,9)

% -------------------------------------------------------------------------

R = 6378*1e3;       % earth radius in m

rhow = 1027;        % density of water (kg/m^3)

rhoa = 1.2;         % density of air   (kg/m^3)

rhoi = 850;         % density of shelf ice (kg/m^3), Silva et al (2006)

drho = rhow-rhoi;

Cw = 0.9;           % bulk coefficient water (Bigg et al 1997)

Ca = 1.3;           % bulk coefficient air   (Bigg et al 1997)

Om = 7.2921*1e-5;   % rotation rate of earth (rad/s)

ff = @(lat) 2*Om*sin(abs(lat)*pi/180);  % latitude in degrees

ga  = sqrt(rhoa*drho/rhow/rhoi*Ca/Cw);  % gamma = sqrt(ca/cw), Eq. 7

S  = @(l,w) l.*w./(l+w);                % harmonic mean length

La = @(u,lat,S) Cw*ga/ff(lat).*u/S/pi;  % Lambda, Eq. 9

% alpha and beta (eq 8 in WDE17)-------------------------------------------

a = @(La) 1./(2*La.^3).*(sqrt(1+4*La.^4)-1);

b = @(La) 1./(sqrt(2)*La.^3).*sqrt((1+La.^4).*sqrt(1+4*La.^4)-3*La.^4-1);

% Melt parameters (WDE17 Appendix)-----------------------------------------

Ti = -4;

a1 = 8.7e-6; a2 = 5.8e-7;

b1 = 8.8e-8; b2 = 1.5e-8;

c  = 6.7e-6;

% -------------------------------------------------------------------------

% specify the space domain ------------------------------------------------

LAT = double(input.latw); LON = double(input.lonw);

if strcmp(model,'CCSM4')

    LAT2 = LAT(:); LON2 = LON(:);

    vec = reshape(1:length(LON2),151,192);

end

minLAT = min(LAT(:)); maxLAT = max(LAT(:));

minLON = min(LON(:)); maxLON = max(LON(:));

msk = input.landmask;

% -------------------------------------------------------------------------

% set run parameters ------------------------------------------------------

if strcmp(model,'ECCO2') % Input fields time step in days

    DT = 3;

else

    DT = 1;

end

trajnum = 25;            % total number of iceberg trajectories to compute

final_t = 366/DT;        % time span of input field

startrange = round(linspace(1,final_t/2,trajnum)); %evenly space seeding

dt = 24*3600;            % model timestep in seconds (dt = 1 day)

dtR = dt/R*180/pi;       % need this ratio for distances

% -------------------------------------------------------------------------

t = 1:final_t;                  %input time

nt= length(t)*DT;               %number of model timesteps

tt = linspace(1,length(t),nt);  %model time

% -------------------------------------------------------------------------

% Load Seeding fields -----------------------------------------------------

load(strcat(model,'_Laurentide_Seed'))

seed_X = repmat(Seed_X(:),[100,1]); %cycle through each location 100x

seed_Y = repmat(Seed_Y(:),[100,1]); %i.e. this can run 3600 icebergs

% -------------------------------------------------------------------------

% these are the circulation fields-----------------------------------------

uwF = input.uw(:,:,t); vwF = input.vw(:,:,t);   %water vels input

uaF = input.ua(:,:,t); vaF = input.va(:,:,t);   %air vels input

sst = input.sst(:,:,t);                         %sst vels input

% -------------------------------------------------------------------------

% loop over individual initial iceberg size classes -----------------------

% -------------------------------------------------------------------------

bvec = 1:10;   %vector of which size classes to compute - has to be [1,10]

bergdims=[100          67          67  % [L W H]

    200         133         133

    300         200         200

    400         267         267

    500         333         300

    600         400         300

    750         500         300

    900         600         300

    1200         800         300

    1500        1000         300];

% -------------------------------------------------------------------------

for bb = bvec

    % ---------------------------------------------------------------------

    bergsize = bb;   % current berg size class

    fprintf('run bergsize B%d \n',bergsize)

    % ---------------------------------------------------------------------

    XIL = nan(trajnum,nt); YIL = XIL; VOL = XIL;        %set output  arrays

    % ---------------------------------------------------------------------

    L = bergdims(bergsize,1); %initialize iceberg dimensions

    W = bergdims(bergsize,2);

    H = bergdims(bergsize,3);

    % ---------------------------------------------------------------------

    % run drift and melt---------------------------------------------------

    mm=0; ss=0; ob=0; %set counters for melted, survived, and escaped bergs

    for j = 1:trajnum

        if mod(j,10)==0; fprintf('%d trajectories computed \n',j); end

        

        ts = startrange(j);            %pick seeding time (of input field)

        tts= ts*DT;                    %trajectory start time (of model)

        lt = nt-tts;                   %trajectory run length

        

        xil = nan(1,lt); yil = xil;  v = xil;   %initialize output vectors

        yig = seed_Y(j); xig = seed_X(j);   % cycle through start locations

        if strcmp(model,'ECCO2')

            xil(1) = LON(xig); yil(1) = LAT(yig);   %initial lon and lat

        else

            xil(1) = LON(yig,xig); yil(1) = LAT(yig,xig);

        end

        l = L*ones(1,lt); w = l*W/L; h = l*H/L; %initial berg dimensions

        v(1) = L*W*H;                           %initial volume and dvol

 

        % integrate as long as the iceberg is in the domain and not

        % melted and over the time period specified above

        i = 0; outofbound = 0; melted = 0;

        while outofbound == 0 && melted == 0 && i<lt-1

            i = i+1;

            % -------------------------------------------------------------

            % Computes iceberg drift component, WDE17, Section 3

            % -------------------------------------------------------------

            % find nearest neighbor of surface condition field

            if strcmp(model,'ECCO2') % this only works on a regular grid

                YI = dsearchn(LAT,yil(i));

                XI = dsearchn(LON,xil(i));

            else  %CCSM4 is on an irregular grid

                indic = dsearchn([LAT2,LON2],[yil(i),xil(i)]);

                [XI,YI]=find(indic==vec);

            end

            % now interpolate fields linearly between timesteps------------

            timestep = tt(tts+i);

            t1  = floor(timestep); t2 = t1+1;

            dt1 = timestep-t1; dt2 = t2-timestep;

            ua = uaF(XI,YI,t1)*dt1+uaF(XI,YI,t2)*dt2;

            va = vaF(XI,YI,t1)*dt1+vaF(XI,YI,t2)*dt2;

            uw = uwF(XI,YI,t1)*dt1+uwF(XI,YI,t2)*dt2;

            vw = vwF(XI,YI,t1)*dt1+vwF(XI,YI,t2)*dt2;

            SST= sst(XI,YI,t1)*dt1+sst(XI,YI,t2)*dt2;

            % compute wind speed and Lambda at location (for given iceberg 

            % size)--------------------------------------------------------

            Ua = sqrt(ua^2+va^2);

            LA = La(Ua,yil(i),S(l(i),w(i)));

            % now compute iceberg velocity Eq (6)--------------------------

            ui = uw + ga*( a(LA)*va + b(LA)*ua);

            vi = vw + ga*(-a(LA)*ua + b(LA)*va);

            % iceberg translation (convert from m to deg lat/lon)----------

            dlon = ui*dtR;

            dlat = vi*dtR;

            yil(i+1) = yil(i) + dlat;

            xil(i+1) = xil(i) + dlon/cos((yil(i+1)+yil(i))/2*pi/180);

            % check you haven't gone out of bounds-------------------------

            if xil(i+1)>maxLON || xil(i+1)<minLON || ...

                    yil(i+1)>maxLAT || yil(i+1)<minLAT

                outofbound = 1;

                ob = ob+1;

                fprintf('iceberg %d left domain at timestep %d \n',j,i);

            else % now check iceberg not on land---------------------------

                if strcmp(model,'ECCO2')

                    yi2(1) = find(LAT<=yil(i+1),1,'last');

                    yi2(2) = find(LAT>yil(i+1),1,'first');

                    xi2(1) = find(LON<=xil(i+1),1,'last');

                    xi2(2) = find(LON>xil(i+1),1,'first');

                    % when new position within one grid box of land:

                    if any(find(msk(xi2,yi2)==0))

                        yil(i+1) = yil(i);  %iceberg reset to time step i

                        xil(i+1) = xil(i);

                    end

                else % similar for CCSM4

                    indic = dsearchn([LAT2,LON2],[yil(i+1),xil(i+1)]);

                    [XI,YI]=find(indic==vec);

                    if input.landmask(XI,YI)==1

                        yil(i+1) = yil(i);

                        xil(i+1) = xil(i);

                    end

                end

            end

            % -------------------------------------------------------------

            % Here we compute the iceberg melting component, WDE15 Appendix

            % -------------------------------------------------------------

            % Compute Melt Terms Eq A1 in WDE17----------------------------

            Me = a1*Ua^0.5 + a2*Ua;

            Mv = b1*SST + b2*SST^2;

            Mb = c*sqrt((ui - uw).^2+(vi-vw).^2)^0.8*(SST-Ti)*l(i)^(-0.2);

            % Apply Melt Rates --------------------------------------------

            dldt = - Mv - Me; dhdt = - Mb;

            l(i+1) = l(i)+dldt*dt;

            w(i+1) = w(i)+dldt*dt;

            h(i+1) = h(i)+dhdt*dt;

            % -------------------------------------------------------------

            % Make sure the berg is not negative size

            if l(i+1)<0 || w(i+1)<0 || h(i+1)<0

                l(i+1)=0; w(i+1)=0; h(i+1)=0;

                melted = 1;

                mm = mm+1;

            end

            % Rollovers (Eq A2 in WDE17)

            if w(i+1)/h(i+1) < sqrt(6*rhoi/rhow*(1-rhoi/rhow))

                hn = w(i+1); w(i+1) = h(i+1); h(i+1) = hn;

            end

            % Make sure length>width:

            if w(i+1)>l(i+1)

                wn = l(i+1); l(i+1)=w(i+1); w(i+1) = wn;

            end

            % Compute new volume

            v(i+1) = l(i+1)*w(i+1)*h(i+1);

            % Check whether iceberg survived-------------------------------

            if i == lt-1 && v(i+1) > 0

                ss = ss+1;

            end

        end

        ind = 1:i+1;

        XIL(j,ind)=xil(ind); YIL(j,ind)=yil(ind);   %store trajectory

        VOL(j,ind)=v(ind);                          %store volume

    end

    % ---------------------------------------------------------------------

    fprintf('%d icebergs died, %d lived, %d left the domain \n',mm,ss,ob)

    % ---------------------------------------------------------------------

    save(sprintf('%s_Size%d',model,bb),'XIL','YIL','VOL'); %save output

end

% -------------------------------------------------------------------------

% Plot Iceberg Trajectories------------------------------------------------

%

% This loads in output files from "iceberg_model.m" for different iceberg

% sizes and plots them (colorcoded). 

%

% Requires the matlab mapping toolbox. I made the atl.mat file to quickly

% plot a low-resolution landmask of the north atlantic. 

% -------------------------------------------------------------------------

 

figure(1); clf

% make North Atlantic Map plot (requires mapping toolbox)------------------

load atl

latlim = [33 70];

lonlim = [-80 0];

worldmap(latlim,lonlim)

gri = gridm('off');

geo = geoshow(atl(1:10), 'FaceColor', [.9 .9 .9]);

setm(gca,'FFaceColor',[1 1 1])

hold on

hm = mlabel('on');

set(hm(:),'fontsize',13);

set(hm(:),'verticalalignment','baseline');

pm = plabel('on');

set(pm(:),'fontsize',13);

set(pm(:),'verticalalignment','baseline');

 

%uncomment the following if you've precomputed output files and -----------

%want to call just the plotting cell --------------------------------------

% model = 'ECCO2'; 

 

N = 25;   %number of trajectories to be plotted

bvec = [9 7 5 3 1];  %specify different iceberg sizes you're interested in

p = cell(1,length(bvec)); pl = p; leg = p;

ii = 0;

for i  = bvec

    load(sprintf('%s_Size%d',model,i))

    bb = i; bs = bb;

    XL = XIL;

    YL = YIL;

    ind = 1:N;

    ii = ii+1;

    for tind = 1:length(ind)

        vend = find(VOL(ind(tind),:)<.1*VOL(ind(tind),1),1,'first'); 

        %plot until iceberg is 90% decayed

        if isempty(vend)

            vend = size(VOL,2);

        end

        if tind == 1

            p{ii} = plotm(YL(ind(1),1:vend)',XL(ind(1),1:vend)'...

                ,'-','col',1-[bb/15 1-bb/15 1-bb/15]);

        else

            plotm(YL(ind(tind),1:vend)',XL(ind(tind),1:vend)'...

                ,'-','col',1-[bb/15 1-bb/15 1-bb/15]);

        end

    end

end

 

for i = 1:ii; pp = p{i}; pl{i}=pp(1); leg{i} = sprintf('%d',bvec(i)); end

ll = legend([pl{:}],leg(:),'location','northeastoutside',...

    'interpreter','latex');

legendtitle(ll,'Size Class','fontweight','normal');

title(sprintf('%s Forcing',model))