[MITgcm-support] runoff in coupled runs
David Ferreira
dfer at mit.edu
Thu Jun 26 07:21:13 EDT 2014
Hi Maura,
Here is an example of matlab code that produces a run-off map file. This
is for an idealized configuration (2 NH ridges and a circular Antarctic
at the South pole), but it gives the method.
It also produces a bathy file, a land fraction file, an albedo file and
a vegetation fraction file. All are prescribed fields in the coupled
model. Other fields (e.g. ground water content) you leave aside as they
are prognostics (Except if you want to initialize them with some
specific value otherwise the model does it with some default value).
cheers,
david
On 6/25/14 9:40 PM, Maura Brunetti wrote:
> Hello
>
> does anyone know how to construct the runoff file (as for example runOff_cs32_3644.bin in MITgcm/verification/cpl_aim+ocn/input_cpl)?
> I need to run a coupled atm/ocean simulation for the jurassic period, thus I need to give this input field…
> Any help is greatly appreciated!
>
>
> Maura
>
>
>
> --
> Maura Brunetti
> Institute for Environmental Sciences (ISE)
> University of Geneva
> Site de Battelle / D, 7 route de Drize
> CH-1227 Carouge / GE -- Switzerland
>
>
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support
-------------- next part --------------
kg=1; kwr=0;
if kg < 2
grdA=rdmnc('/net/ds-08/scratch-0/dfer/Double_Dr_Ref_bio/GridAtm/grid.*');
arc=grdA.rA;
[ncx nc]=size(arc); npts=ncx*nc;
xcs=grdA.XC;
ycs=grdA.YC;
xcg=grdA.XG(1:ncx,1:nc);
ycg=grdA.YG(1:ncx,1:nc);
ar=reshape(arc,npts,1);
nh=nc/2; np=nh+1;
end
if kg < 3
lndFrc = zeros(24,6,24);
lndFrc(np:nc,3,nh)=1;
lndFrc(nh,3,1:nh)=1;
lndFrc(1:21,4,nh)=1;
lndFrc(nh,2,4:nc)=1;
lndFrc=reshape(lndFrc,24*6,24);
%lndFrc(grdA.YC<=-71) = 2;
lndFrc(grdA.YC<=-75) = 1;
figure('position',[66 109 519 710],'Paperposition',[.25 .5 8 10])
subplot(2,2,1)
plotcube(xcg,ycg,lndFrc)
set(gca,'PlotBoxAspectRatio',[1 1 1])
set(gca,'visible','off')
set(gca,'CameraPosition',[-10 10 9])
set(gca,'CameraViewAngle',6)
caxis([-.2 1.7])
subplot(2,2,2)
plotcube(xcg,ycg,lndFrc)
set(gca,'PlotBoxAspectRatio',[1 1 1])
set(gca,'visible','off')
set(gca,'CameraPosition',[-10 10 -9])
set(gca,'CameraViewAngle',6)
caxis([-.2 1.7])
subplot(2,1,2)
merccube_mod(xcg,ycg,lndFrc)
caxis([-.2 1.7])
set(gca,'yaxislocation','right','fontsize',18,'tickdir','out')
title('Miocene')
%print -deps Miocene.eps
Depth = lndFrc;
Depth( lndFrc == 0 ) = -4000;
Depth( lndFrc == 1 ) = 0;
%- run-Off catchment area:
mapL=zeros(npts,1); mapO=zeros(npts,1);
clear ijA ijO roff
n=0;
%- treat face 2: lndFrc(nh,2,4:nc)=1;
f=2;
i=nh;
for j=5:nc
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j-1)*ncx+(f-1)*nc+i+1;
roff(n)=ar(ijA(n))/2;
end
for j=5:nc
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j-1)*ncx+(f-1)*nc+i-1;
roff(n)=ar(ijA(n))/2;
end
j = 4;
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j-1)*ncx+(f-1)*nc+i+1;
roff(n)=ar(ijA(n))/3;
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j-1)*ncx+(f-1)*nc+i-1;
roff(n)=ar(ijA(n))/3;
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j-1-1)*ncx+(f-1)*nc+i;
roff(n)=ar(ijA(n))/3;
%--
%- treat face 3:
% lndFrc(nh,3,1:nh)=1;
f=3;
i=nh;
for j=1:nh-1
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j-1)*ncx+(f-1)*nc+i+1;
roff(n)=ar(ijA(n))/2;
end
for j=1:nh-1
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j-1)*ncx+(f-1)*nc+i-1;
roff(n)=ar(ijA(n))/2;
end
% lndFrc(np:nc,3,nh)=1;
f=3;
j=nh;
for i=np:nc
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j-1-1)*ncx+(f-1)*nc+i;
roff(n)=ar(ijA(n))/2;
end
for i=np:nc
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j+1-1)*ncx+(f-1)*nc+i;
roff(n)=ar(ijA(n))/2;
end
% Joining point of the 2 ridges at NP
f=3;
i=nh;
j=nh;
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j+1-1)*ncx+(f-1)*nc+i;
roff(n)=ar(ijA(n))/3;
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j-1)*ncx+(f-1)*nc+i-1;
roff(n)=ar(ijA(n))/3;
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j+1-1)*ncx+(f-1)*nc+i-1;
roff(n)=ar(ijA(n))/3;
%- treat face 4: lndFrc(1:21,4,nh)=1;
f=4;
j=nh;
for i=1:20
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j-1-1)*ncx+(f-1)*nc+i;
roff(n)=ar(ijA(n))/2;
end
for i=1:20
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j+1-1)*ncx+(f-1)*nc+i;
roff(n)=ar(ijA(n))/2;
end
i = 21;
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j-1-1)*ncx+(f-1)*nc+i;
roff(n)=ar(ijA(n))/3;
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=(j+1-1)*ncx+(f-1)*nc+i;
roff(n)=ar(ijA(n))/3;
n=n+1;
ijA(n)=(j-1)*ncx+(f-1)*nc+i;
ijO(n)=ijA(n)+1;
roff(n)=ar(ijA(n))/3;
%--
%- treat face 6:
% lndFrc(np,6,1:nh)=1;
f=6;
Yctmp = reshape(grdA.YC,24,6,24); YC6 = squeeze(Yctmp(:,6,:));
Xctmp = reshape(grdA.XC,24,6,24); XC6 = squeeze(Xctmp(:,6,:));
FaceL6 = zeros(24,24);
FaceL6(YC6 <= -75)=1;
FaceO6 = zeros(24,24);
FaceO6(YC6 > -75 & YC6 <= -71)=1;
IL = find(FaceL6==1);
IO = find(FaceO6==1);
figure
pcol((0:23),(0:23),FaceL6-FaceO6); hold on; grid on
set(gca,'xtick',[0:24],'ytick',[0:24])
for l=1:length(IL)
n=n+1;
clear dist
jl = floor(IL(l)/24)+1;
il = IL(l)-(jl-1)*24;
plot(il-0.5,jl-0.5,'w+','linewidth',2)
ijA(n)=(jl-1)*ncx+(f-1)*nc+il;
for o=1:length(IO)
%dist(o) = m_lldist([XC6(IL(l)) XC6(IO(o))],[YC6(IL(l)) YC6(IO(o))]);
jo = floor(IO(o)/24)+1;
io = IO(o)-(jo-1)*24;
dist(o) = sqrt((il-io)^2+(jl-jo)^2);
end
[a o]=min(dist);
jo = floor(IO(o)/24)+1;
io = IO(o)-(jo-1)*24;
plot(io-0.5,jo-0.5,'k+','linewidth',2)
ijO(n)=(jo-1)*ncx+(f-1)*nc+io;
roff(n)=ar(ijA(n));
pause(.5)
end
nRO=n;
for n=1:nRO,
mapL(ijA(n))=mapL(ijA(n))+roff(n);
mapO(ijO(n))=mapO(ijO(n))+roff(n);
end
mapL=reshape(mapL,ncx,nc);
mapO=reshape(mapO,ncx,nc);
end
%----------
% verification
figure
merccube_mod(xcg,ycg,mapL./arc-lndFrc)
colorbar
figure;clf
cbV=2; shift=-1; ccB=[0 0]; AxBx=[-180 180 -90 90];
ipt=79; jpt=11;
var=mapL+mapO;
%var=mapO;
var=var./arc;
[fac]=grph_CS(var,xcs,ycs,xcg,ycg,ccB(1),ccB(2),shift,cbV,AxBx);
if ipt*jpt ~= 0, hold on; plot(fac*xcs(ipt,jpt),fac*ycs(ipt,jpt),'rx'); hold off; end
%title(['month = ',int2str(nm)]);
%nResp= input('Continue: 1 ; Stop: 0 ; answer ? ');
%end ; end
colorbar
figure
plotcube(xcg,ycg,var)
set(gca,'PlotBoxAspectRatio',[1 1 1])
set(gca,'visible','off')
set(gca,'CameraPosition',[-13. -10. -3.])
set(gca,'CameraViewAngle',6)
colorbar
figure
merccube_mod(xcg,ycg,var)
colorbar
%----------
nameC = 'Antarct_DDr';
if kwr > 0
file_name=['landFrc_cs24.' nameC '.bin'];
fid=fopen(file_name,'w','b');
fwrite(fid,lndFrc,'real*8'); fclose(fid);
fprintf(['write Land Fraction on file: ',file_name,'\n']);
%--
tmp=zeros(3,nRO); tmp(1,:)=ijA; tmp(2,:)=ijO; tmp(3,:)=roff;
file_name=['runOff_cs24_',int2str(nRO),'.bin'];
fid=fopen(file_name,'w','b');
fwrite(fid,tmp,'real*8'); fclose(fid);
fprintf(['write RunOff map on file: ',file_name,'\n']);
%- Other input file:
albedo=25; vegetFrac=75;
var=albedo*lndFrc;
var(grdA.YC <= -75) = 50;
file_name=['albedo_',int2str(albedo),'pc.' nameC '.bin'];
fid=fopen(file_name,'w','b');
fwrite(fid,var,'real*8'); fclose(fid);
fprintf(['write Albedo on file: ',file_name,'\n']);
var=vegetFrac*lndFrc;
file_name=['vegFrc_',int2str(vegetFrac),'pc.' nameC '.bin'];
fid=fopen(file_name,'w','b');
fwrite(fid,var,'real*8'); fclose(fid);
fprintf(['write Vegetation Frac. on file: ',file_name,'\n']);
file_name=[nameC '.c24.Bathy.4p0km.bin'];
fid=fopen(file_name,'w','b');
fwrite(fid,Depth,'real*8'); fclose(fid);
fprintf(['write Bathy. on file: ',file_name,'\n']);
end
%---
return
More information about the MITgcm-support
mailing list