[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