kr=0.25
ke=0.8
ks=1.5
kc=0.1

for(i=0;i<size[argument0];i+=1)
{
    for(j=0;j<size[argument0];j+=1)
    {
        water[i,j]=random(1)
        sediment[i,j]=0
        altitude[i,j]=ds_grid_get(argument0,i,j)+water[i,j]
    }
}

repeat(10)
{
    for(i=0;i<size[argument0];i+=1)
    {
        for(j=0;j<size[argument0];j+=1)
        {
            ds_grid_set(argument0,i,j,ds_grid_get(argument0,i,j)-ks*water[i,j])

            if i!=0 && j!=0 && i!=size[argument0]-1 && j!=size[argument0]-1
            {
                water[i,j]=water[i,j]+kr
                altitude[i,j]=ds_grid_get(argument0,i,j)+water[i,j]
                sediment[i,j]=sediment[i,j]+ks*water[i,j]

                dtotal=0
                if altitude[i,j]-altitude[i+1,j]>0 then dtotal+=altitude[i,j]-altitude[i+1,j]
                if altitude[i,j]-altitude[i-1,j]>0 then dtotal+=altitude[i,j]-altitude[i-1,j]
                if altitude[i,j]-altitude[i,j+1]>0 then dtotal+=altitude[i,j]-altitude[i,j+1]
                if altitude[i,j]-altitude[i,j-1]>0 then dtotal+=altitude[i,j]-altitude[i,j-1]

                delta_a=altitude[i,j]-mean(altitude[i+1,j],altitude[i-1,j],altitude[i,j+1],altitude[i,j-1])
                if delta_a<0 then delta_a=0

                if dtotal>0
                {
                    for(k=0;k<4;k+=1;)
                    {
                        delta_w=0
                        delta_s=0
                        if k=0 then di=altitude[i,j]-altitude[i+1,j]
                        if k=1 then di=altitude[i,j]-altitude[i-1,j]
                        if k=2 then di=altitude[i,j]-altitude[i,j+1]
                        if k=3 then di=altitude[i,j]-altitude[i,j-1]

                        if k=0 && di>0 && water[i,j]>0
                        {
                            delta_w=min(water[i,j],delta_a)*di/dtotal
                            delta_s=sediment[i,j]*delta_w/water[i,j]
                            water[i+1,j]+=delta_w
                            sediment[i+1,j]+=delta_s
                        }
                        if k=1 && di>0 && water[i,j]>0
                        {
                            delta_w=min(water[i,j],delta_a)*di/dtotal
                            delta_s=sediment[i,j]*delta_w/water[i,j]
                            water[i-1,j]+=delta_w
                            sediment[i-1,j]+=delta_s
                        }
                        if k=2 && di>0 && water[i,j]>0
                        {
                            delta_w=min(water[i,j],delta_a)*di/dtotal
                            delta_s=sediment[i,j]*delta_w/water[i,j]
                            water[i,j+1]+=delta_w
                            sediment[i,j+1]+=delta_s

                        }
                        if k=3 && di>0 && water[i,j]>0
                        {
                            delta_w=min(water[i,j],delta_a)*di/dtotal
                            delta_s=sediment[i,j]*delta_w/water[i,j]
                            water[i,j-1]+=delta_w
                            sediment[i,j-1]+=delta_s
                        }

                        if delta_w!=0 && delta_s!=0 && water[i,j]!=0
                        {
                            water[i,j]-=delta_w
                            sediment[i,j]-=delta_s
                        }
                    }
                }

            }
            water[i,j]*=ke
            m_max=kc*water[i,j]
            delta_m=max(0,sediment[i,j]-m_max)

            sediment[i,j]=sediment[i,j]-delta_m
            ds_grid_set(argument0,i,j,ds_grid_get(argument0,i,j)+delta_m)
        }
    }