diff --git a/Rivers/locate_runoff.m b/Rivers/locate_runoff.m
index cb5e978a43d3e34f9c5ab5caa177f82249955dfa..547004132174a70c538562e724d27a8d62541cc0 100644
--- a/Rivers/locate_runoff.m
+++ b/Rivers/locate_runoff.m
@@ -1,8 +1,46 @@
-function [j2,i2]=locate_runoff(dir,j,i,mask,masku,maskv)
-%
+function [j2for_out,i2for_out, j2for,i2for,j2,i2]=locate_runoff(dir,j,i,mask,masku,maskv)
+%=======================================================================================================
+% input :
% j,i : first guess river position
% dir : direction and sense of the flow
+%
+% output :
+% j2,i2 : index position (eta, xi) at rho-point for runoff, in matlab convention
+% in MATLAB, rho, u and v indexing start at 1 !
+%
+% MATLAB
+% |-------- v(i,j) -------|
+% | |
+% | |
+% u(i-1,j) r(i,j) u(i,j)
+% | |
+% | |
+% |-------- v(i,j-1) -----|
+%
+% j2for, i2for : index position (eta, xi) at rho-point for runoff, in CROCO convention
+% in CROCO fortran, rho start at 0, u and v indexing start at 1 !,
+%
+% zero index is a ghost cell, with no staggerd u-point on right and v point below
+%
+%
+% CROCO fortran
+% |-------- v(i,j+1) -----|
+% | |
+% | |
+% u(i,j) r(i,j) u(i+1,j)
+% | |
+% | |
+% |-------- v(i,j) -------|
+%
+% j2for_out, i2for_out : index for runoff positions, at u and v cells, in CROCO fortran index convention.
+% i- and j- are of the u or v cells that are flowing into the target wet cell
+%
+% It is computed depending :
+% - the runoff direction (zonal=0, meridian=1)
+% - the sense (east-west/south-north =1 ; west-east/north-south =-1
+% - the landmask around ...
%
+%=======================================================================================================
if mask(j,i)==1
disp('River positionned in sea')
insea=1;
@@ -12,76 +50,119 @@ else
end
if dir(1)==0
if insea
- if dir(2)==-1 %est - west => TESTED
- while masku(j,i)==1
- i=i+1;
- %disp(['MASKU:',num2str(masku(j,i))])
- end
- %disp(['--'])
- %disp(['MASKU:',num2str(masku(j,i))])
- elseif dir(2)==1 % west - est => TESTED
- while masku(j,i)==1
+ if dir(2)==1 % >> : west - east => TESTED
+ while mask(j,i) == 1
i=i-1;
- %disp(['MASKU:',num2str(masku(j,i))])
+ disp(['MASK:',num2str(mask(j,i))])
+ end
+ disp(['--'])
+ disp(['MASK:',num2str(mask(j,i))])
+ disp(['MASKU:',num2str(masku(j,i))])
+ elseif dir(2)==-1 % << : east - west => TESTED
+ while mask(j,i) == 1
+ i=i+1;
+ disp(['MASK:',num2str(mask(j,i))])
end
- %disp(['--'])
- %disp(['MASKU:',num2str(masku(j,i))])
+ %in this case add +1 in i to get back into last wet cell
+ j=j+1;
+ disp(['--'])
+ disp(['MASK:',num2str(mask(j,i))])
+ disp(['MASKU:',num2str(masku(j,i))])
end
else %inland
- if dir(2)==1 %west-est => TESTED
- while masku(j,i)~=1
+ if dir(2) == 1 % >> : west-est => TESTED
+ while mask(j,i) ~= 1
i=i+1;
- %disp(['MASKU:',num2str(masku(j,i))])
+ disp(['MASK:',num2str(mask(j,i))])
end
disp(['--'])
- i=i-1
- %disp(['MASK:',num2str(masku(j,i))])
- elseif dir(2)==-1 %est-west => TESTED
- while masku(j,i)~=1
+ disp(['MASK:',num2str(mask(j,i))])
+ disp(['MASKU:',num2str(masku(j,i))])
+ elseif dir(2) == -1 % << : east-west => TESTED
+ while mask(j,i) ~= 1
i=i-1;
- %disp(['MASKU:',num2str(masku(j,i))])
+ disp(['MASK:',num2str(mask(j,i))])
end
- %disp(['--'])
- i=i+1;
- %disp(['MASKU:',num2str(masku(j,i))])
+ disp(['--'])
+ disp(['MASK:',num2str(mask(j,i))])
+ disp(['MASKU:',num2str(masku(j,i))])
end
end
else %dir(k,1)=1
if insea
- if dir(2)==-1 % nord - sud => TESTED
- while maskv(j,i)==1
- j=j+1;
- %disp(['MASKV:',num2str(maskv(j,i))])
- end
- %disp(['--'])
- %disp(['MASKV:',num2str(maskv(j,i))])
- elseif dir(2)==1 % sud - nord => TESTED
- while maskv(j,i)==1
+ if dir(2) == 1 % ^ : sud - nord => TESTED
+ while mask(j,i) == 1
j=j-1;
- %disp(['MASKV:',num2str(maskv(j,i))])
+ disp(['MASK:',num2str(mask(j,i))])
end
- %disp(['--'])
- %disp(['MASKV:',num2str(maskv(j,i))])
+ %in this case add +1 in j to get back into last wet cell
+ j=j+1;
+ %
+ disp(['-- ^^ start insea'])
+ disp(['MASK:',num2str(mask(j,i))])
+ disp(['MASKV:',num2str(maskv(j,i))])
+
+ elseif dir(2) == -1 % v : nord - sud => TESTED
+ while mask(j,i) == 1
+ j=j+1;
+ disp(['MASK:',num2str(mask(j,i))])
+ end
+ disp(['--'])
+ disp(['MASK:',num2str(mask(j,i))])
+ disp(['MASKV:',num2str(maskv(j,i))])
+
end
else %inland
- if dir(2)==1 %sud-nord => TESTED
- while maskv(j,i)~=1
+ if dir(2) == 1 % ^: sud-nord => TESTED
+ while mask(j,i) ~= 1
j=j+1;
- %disp(['MASKV:',num2str(maskv(j,i))])
+ disp(['MASK:',num2str(mask(j,i))])
end
- %disp(['--'])
- j=j-1;
- %disp(['MASKV:',num2str(maskv(j,i))])
- elseif dir(2)==-1 %nord-sud => TESTED
- while maskv(j,i)~=1
+ disp(['--^^ start inland'])
+ disp(['MASK:',num2str(mask(j,i))])
+ disp(['MASKV:',num2str(maskv(j,i))])
+
+ elseif dir(2) == -1 % v : nord-sud => TESTED
+ while mask(j,i) ~= 1
j=j-1;
- %disp(['MASKV:',num2str(maskv(j,i))])
+ disp(['MASK:',num2str(mask(j,i))])
end
- %disp(['--'])
- j=j+1;
- %disp(['MASKV:',num2str(maskv(j,i))])
+ disp(['--'])
+ disp(['MASK:',num2str(mask(j,i))])
+ disp(['MASKV:',num2str(maskv(j,i))])
end
end
end
+
+% from matlab to roms index convention
j2=j; i2=i;
+j2for=j-1; i2for=i-1;
+
+% regarding the flow direction, sense and landmask surrounding
+% -- toward west >>
+if ( dir(1) == 0 & dir(2) == 1 ) ;
+ j2for_out = j2for ;
+ i2for_out = i2for ;
+end
+% -- toward est <<
+if ( dir(1) == 0 & dir(2) == -1 );
+ j2for_out = j2for ;
+ i2for_out = i2for + 1 ;
+end
+% -- toward north ^^
+if ( dir(1) == 1 & dir(2) == 1 );
+ j2for_out = j2for ;
+ i2for_out = i2for ;
+end
+% -- toward south vv
+if ( dir(1) == 1 & dir(2) == -1 );
+ j2for_out = j2for + 1 ;
+ i2for_out = i2for ;
+end
+
+% j2mat_out is equal to j2for_out
+% i2mat_out is equal to j2for_out % I think ...
+% j2mat_out = j2for_out
+% i2mat_out = i2for_out
+
return
diff --git a/Rivers/make_runoff.m b/Rivers/make_runoff.m
index 8d42e3c332ff61de40d56ef791c10fadebfa2794..3a532fc5dd1904e55f2f40395cbbb8120b3b0e39 100644
--- a/Rivers/make_runoff.m
+++ b/Rivers/make_runoff.m
@@ -25,6 +25,7 @@
%
%
% July 2013: G. Cambon (IRD/LEGOS) & M. Herrmann (IRD/LEGOS)
+% July 2023 : G. Cambon (IRD/LOPS)
%
clear all
close all
@@ -134,7 +135,16 @@ if define_dir==1
% $$$ dir(10,:)= [1 , 1]; % # Magdalena
% $$$ dir(11,:)= [1, -1]; % # Urugay same dir/sense as Parana
% $$$ dir(12,:)= [0 , 0]; % # Danube (not used)
-% $$$ dir(13,:)= [0 , -1]; % # Niger
+% $$$ dir(13,:)= [0 , -1]; % # Niger
+
+%ARVOR FAMILY for river # 1 to 6 without Danube
+ dir(1,:) = [0 , 1]; % # Orinoco
+ dir(2,:) = [0 , 1]; % # Mississipi
+ dir(3,:) = [0 , 1]; % # Saint Lawrence
+ dir(4,:) = [1 , 1]; % # Magdalena
+ dir(5,:) = [0 , 0]; % # Danube (not used)
+ dir(6,:) = [1 , -1]; % # Niger
+
end
%=========================================================================================
%%%%%%%%%%%%%%%%%%% END USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -204,11 +214,22 @@ else
I(k)=i-1;
I(I<1)=1;
J(J<1)=1;
-
+
disp([' Position is approximetly J=',num2str(J(k)),' and I=',num2str(I(k))])
- disp(['lon src in grid (rho point) ~',num2str(lon(J(k),I(k)))])
- disp(['lat src in grid (rho point) ~',num2str(lat(J(k),I(k)))])
+ disp([' lon src in grid (rho point) ~',num2str(lon(J(k),I(k)))])
+ disp([' lat src in grid (rho point) ~',num2str(lat(J(k),I(k)))])
+ end
+ %==========================================
+ mod_manual =0
+ if mod_manual
+ % Modify position runoff manually
+ % River number to modify number to modify
+ kmod = 3;
+ J(kmod) = J(kmod);
+ I(kmod) = I(kmod) - 10 ;
+ disp(['River #',num2str(kmod),': ',char(myrivername(kmod,:)),' modified manually'])
end
+ %==========================================
%%
%% Check the river you really have to process.
%% Remove rivers out of the domain, if any...
@@ -218,15 +239,15 @@ else
%%debug
%%rivertoprocess=1
%=================
-
number_rivertoprocess=length(rivertoprocess);
rivername=strvcat(myrivername(rivertoprocess,:));
rivernumber=number_rivertoprocess;
rivname_StrLen=size(rivername,2);
%% Make a figure
-
+
if plotting==1
+ dl = mean(lat(2:end,1) - lat(1:end-1,1));
fig1=figure(1);
set(fig1,'Position',[659 34 1334 1307],'units','pixels');
m_proj('mercator',... % dl/2 shift to have cells centered on lon/lat
@@ -250,17 +271,17 @@ else
pause
%%
if plotting_zoom
- %% to be adapeted regarding the configuration
-% $$$ extendpointY=[500 500 500 500 ...
-% $$$ 500 600 900 500 ...
-% $$$ 500 500 500 50 ...
-% $$$ 500 ];
-% $$$ extendpointX=[500 500 500 500 ...
-% $$$ 500 600 900 500 ...
-% $$$ 500 500 500 50 ...
-% $$$ 500 ];
- extendpointY=50*ones(number_rivertoprocess,1);
- extendpointX=50*ones(number_rivertoprocess,1);
+ %% to be adapted regarding the configuration
+ % extendpointY=[500 500 500 500 ...
+ % 500 600 900 500 ...
+ % 500 500 500 50 ...
+ % 500 ];
+ % extendpointX=[500 500 500 500 ...
+ % 500 600 900 500 ...
+ % 500 500 500 50 ...
+ % 500 ];
+ extendpointY=20*ones(number_rivertoprocess,1);
+ extendpointX=20*ones(number_rivertoprocess,1);
nfig0=5;
for k=1:number_rivertoprocess
nfig=nfig0+k;
@@ -282,37 +303,25 @@ else
mylat_point=lat(Jzoom_start : Jzoom_end,Izoom_start : Izoom_end);
mymask_point=mask(Jzoom_start : Jzoom_end,Izoom_start : Izoom_end);
-% $$$ %%-Replace by specific case for GIGATL1
-% $$$ if (k ~= 12 )
-% $$$ mylon_point=lon(J(k)+1 - extendpointYkk : J(k)+1 + extendpointYkk,...
-% $$$ I(k)+1 - extendpointXkk : I(k)+1 + extendpointXkk);
-% $$$ mylat_point=lat(J(k)+1 - extendpointYkk : J(k)+1 + extendpointYkk,...
-% $$$ I(k)+1 - extendpointXkk : I(k)+1 + extendpointXkk);
-% $$$ mymask_point=mask(J(k)+1 - extendpointYkk : J(k)+1 + extendpointYkk,...
-% $$$ I(k)+1 - extendpointXkk : I(k)+1 + extendpointXkk);
-% $$$ else
-% $$$ mylon_point=lon(J(k)+1 - extendpointYkk : J(k)+1,...
-% $$$ I(k)+1 - extendpointXkk : I(k)+1 ) ;
-% $$$
-% $$$ mylat_point=lat(J(k)+1 - extendpointYkk : J(k)+1,...
-% $$$ I(k)+1 - extendpointXkk : I(k)+1) ;
-% $$$
-% $$$ mymask_point=mask(J(k)+1 - extendpointYkk : J(k)+1,...
-% $$$ I(k)+1 - extendpointXkk : I(k)+1);
-% $$$ end
-% $$$ %%- end Specific plotting
-
-
- m_proj('mercator',... % dl/2 shift to have cells centered on lon/lat
- 'lon',[min(min(mylon_point-dl/2)) max(max(mylon_point-dl/2)) ],...
- 'lat',[min(min(mylat_point-dl/2)) max(max(mylat_point-dl/2)) ] );
- m_pcolor(mylon_point-dl/2,mylat_point-dl/2,mymask_point) ; shading flat
- m_grid('box','fancy','xtick',5,'ytick',5,'tickdir','out');
- set(findobj('tag','m_grid_color'),'facecolor','white');
- hold on
- [px,py]=m_ll2xy(lon_src,lat_src);
- h11=plot(px,py,'ro');
- title(['\bf', myrivername(k,:)]);
+ % %% -Replace by specific case for GIGATL1
+ % if (k ~= 12 )
+ % mylon_point=lon(J(k)+1 - extendpointYkk : J(k)+1 + extendpointYkk,...
+ % I(k)+1 - extendpointXkk : I(k)+1 + extendpointXkk);
+ % mylat_point=lat(J(k)+1 - extendpointYkk : J(k)+1 + extendpointYkk,...
+ % I(k)+1 - extendpointXkk : I(k)+1 + extendpointXkk);
+ % mymask_point=mask(J(k)+1 - extendpointYkk : J(k)+1 + extendpointYkk,...
+ % I(k)+1 - extendpointXkk : I(k)+1 + extendpointXkk);
+ % else
+ % mylon_point=lon(J(k)+1 - extendpointYkk : J(k)+1,...
+ % I(k)+1 - extendpointXkk : I(k)+1 ) ;
+
+ % mylat_point=lat(J(k)+1 - extendpointYkk : J(k)+1,...
+ % I(k)+1 - extendpointXkk : I(k)+1) ;
+
+ % mymask_point=mask(J(k)+1 - extendpointYkk : J(k)+1,...
+ % I(k)+1 - extendpointXkk : I(k)+1);
+ % end
+ % %%- end Specific plotting
end
end
end
@@ -384,101 +393,33 @@ else
jj=J(k);
ii=I(k);
dir2=dir(k,:);
-
-% $$$ %% Example manual modification for GIGATL family runs
-% $$$ %% GIGATL6 [BASE]: Amazon Tocantins Orenoque Mississipi
-% $$$ %% Tapajos Xingu Magdalena
-% $$$
-% $$$ if k0 == 1 % Amazon
-% $$$ jj=jj+5
-% $$$
-% $$$ end
-% $$$ if k0 == 6 % Tocantins
-% $$$ ii=ii+1
-% $$$ end
-% $$$ if k0 == 3 %Orenoque
-% $$$ jj=jj-12
-% $$$ end
-% $$$ if k0 == 4 %Mississipi
-% $$$ ii=ii-6
-% $$$ end
-% $$$ if k0 == 7 %Tapajos
-% $$$ jj=jj+14
-% $$$ end
-% $$$ if k0 == 9 %Xingu
-% $$$ jj=jj+24
-% $$$ end
-% $$$ if k0 == 10 %Magdalena
-% $$$ ii=ii+7
-% $$$ end
-% $$$
-% $$$ %%====================================================================
-% $$$ %%
-% $$$ %% GIGATL3 (Multiply base by 2 for the jj and ii offset]
-% $$$ %% => Amazon Tocantins Orenoque Mississipi Tapajos Xingu Magdalena
-% $$$ if k0 == 1 % Amazon
-% $$$ jj=jj+10 - 2
-% $$$ end
-% $$$ if k0 == 6 % Tocantins
-% $$$ ii=ii+2
-% $$$ end
-% $$$ if k0 == 3 %Orenoque
-% $$$ jj=jj-24 -2
-% $$$ end
-% $$$ if k0 == 4 %Mississipi
-% $$$ ii=ii-12
-% $$$ end
-% $$$ if k0 == 7 %Tapajos
-% $$$ jj=jj+28
-% $$$ end
-% $$$ if k0 == 9 %Xingu
-% $$$ jj=jj+48
-% $$$ end
-% $$$ if k0 == 10 %Magdalena
-% $$$ ii=ii+14
-% $$$ end
-% $$$
-% $$$ %%====================================================================
-% $$$ %%
-% $$$ %% GIGATL1 (Multiply base by 6 for the jj and ii offset]
-% $$$ %% => Amazon Tocantins Orenoque Mississipi Tapajos Xingu Magdalena
-% $$$ if k0 == 1 % Amazon
-% $$$ jj=jj+24
-% $$$ end
-% $$$ if k0 == 6 % Tocantins
-% $$$ ii=ii+6
-% $$$ end
-% $$$ if k0 == 3 %Orenoque
-% $$$ jj=jj-78 - 21
-% $$$ end
-% $$$ if k0 == 4 %Mississipi
-% $$$ ii=ii-36
-% $$$ end
-% $$$ if k0 == 7 %Tapajos
-% $$$ jj=jj+84
-% $$$ end
-% $$$ if k0 == 9 %Xingu
-% $$$ jj=jj+144
-% $$$ end
-% $$$ if k0 == 10 %Magdalena
-% $$$ ii=ii+42
-% $$$ end
-% $$$
-% $$$ %% End Example manual modifications for GIGATL family runs
-% $$$ %%$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
-
- [jj2,ii2]=locate_runoff(dir2,jj,ii,mask,masku,maskv);
+ %
+ % Compute the u and v cells indexes in CROCO convention
+ [jj2for_out,ii2for_out, jj2for,ii2for,jj2,ii2]=locate_runoff(dir2,jj,ii,mask,masku,maskv);
+ %
+ J2for_out(k)=jj2for_out;
+ I2for_out(k)=ii2for_out;
+ %
+ J2for(k)=jj2for;
+ I2for(k)=ii2for;
+ %
J2(k)=jj2;
I2(k)=ii2;
- disp([char(myrivername(k,:)),' is J=',num2str(J2(k)),' and I=',num2str(I2(k))])
- disp([' '])
+ %
+ disp([char(myrivername(k,:)),' runoff wet cell (rho-point) indices are :'])
+ disp(['matlab J=',num2str(J2(k)),' and I=',num2str(I2(k))])
+ disp(['CROCO fortran conv J=',num2str(J2for(k)),' and I=',num2str(I2for(k))])
+ disp([''])
+ disp(['For u or v cells indices are (CROCO fortran convention)'])
+ disp(['Jsrc=',num2str(J2for_out(k)),' and Isrc=',num2str(I2for_out(k))])
+ disp(['===='])
end
%%
%% Adjust the rivers temperature and salinity
%%
if psource_ncfile_ts==1
disp([' '])
- disp([' Adjust the rivers temperature and salinity '])
+ disp(['Adjust the rivers temperature and salinity '])
%%disp([' Use the closest surface point in the climatology file '])
my_temp_src0=zeros(rivernumber ,length(woa_time));
my_salt_src0=zeros(rivernumber ,length(woa_time));
@@ -513,7 +454,7 @@ else
%%
S=squeeze(ncclim{'salt'}(:,N,j_rho,i_rho))-10; % hum...
S(S<2)=2.0001; % to prevent negative salinities in the
- % equation of state
+ % equation of state
disp([' Use psource_ncfile_ts_auto using S = sclim -10 '])
disp([' Check line 464 in make_runoff.m to change ' ...
'this arbitrary runoff salinity'])
@@ -636,13 +577,12 @@ else
% $$$ 33.0312 ];
%%
my_salt_src0(:,:)=ones(rivernumber,12)*5;
- %%
if makebio
% $$$ my_no3_src0=[0 0 0 0 0 0 0 0 0 0 0 0];
end
end
%%
- %%Effectivelly processed rivers
+ %% Effectivelly processed rivers
%%
my_temp_src= my_temp_src0(find(indomain_last==1),:);
my_salt_src= my_salt_src0(find(indomain_last==1),:);
@@ -662,12 +602,13 @@ else
'lat',[latmin-dl/2 latmax-dl/2]);
for k0=1:number_rivertoprocess
k=rivertoprocess(k0);
+
if dir(k,1)==0 %zonal flow (u-grid)
- lon_src=lonu(J2(k),I2(k));
- lat_src=latu(J2(k),I2(k));
+ lon_src=lonu(J2(k),I2for_out(k)); % Adjusted position
+ lat_src=latu(J2(k),I2for_out(k)); % + nice plot
else %meridional flow (v-grid)
- lon_src=lonv(J2(k),I2(k));
- lat_src=latv(J2(k),I2(k));
+ lon_src=lonv(J2for_out(k),I2(k)); % Adjusted position
+ lat_src=latv(J2for_out(k),I2(k)); % + nice plot
end
[px,py]=m_ll2xy(lon_src,lat_src);
if dir(k,1)==0 & dir(k,2)==1
@@ -693,37 +634,43 @@ else
close(nfig);
figure(nfig);
+ lon_src0=lon(J(k),I(k)); % First guess position
+ lat_src0=lat(J(k),I(k));
if dir(k,1)==0 %zonal flow (u-grid)
- lon_src0=lonu(J(k),I(k)); % First guess position
- lat_src0=latu(J(k),I(k));
- lon_src=lonu(J2(k),I2(k)); % Adjusted position
- lat_src=latu(J2(k),I2(k));
+ lon_src=lonu(J2(k),I2for_out(k)); % Adjusted position
+ lat_src=latu(J2(k),I2for_out(k)); % + nice plot
else %meridional flow (v-grid)
- lon_src0=lonv(J(k),I(k)); % First guess position
- lat_src0=latv(J(k),I(k));
- lon_src=lonv(J2(k),I2(k)); % Adjusted position
- lat_src=latv(J2(k),I2(k));
+ lon_src=lonv(J2for_out(k),I2(k)); % Adjusted position
+ lat_src=latv(J2for_out(k),I2(k)); % + nice plot
end
-
+ %
+ lon_src_r=lon(J2(k),I2(k));
+ lat_src_r=lat(J2(k),I2(k));
+ dlzoom = lat(J2(k),I2(k)) - lat(J2(k)-1,I2(k));
+ %
+ disp(['Final position of the target wet cell (rho-point) is J=',num2str(J2for_out(k)),' and I=',num2str(I2for_out(k))])
+ disp([' lon src in grid ',num2str(lon_src_r)])
+ disp([' lat src in grid ',num2str(lat_src_r)])
+ disp([' '])
+ %
extendpointYkk=extendpointY(k);
extendpointXkk=extendpointX(k);
-
+ %
Izoom_start=max(1,I(k) - extendpointXkk);
Izoom_end=min(size(mask,2),I(k) + extendpointXkk);
Jzoom_start=max(1,J(k) - extendpointYkk);
Jzoom_end=min(size(mask,1),J(k) + extendpointYkk);
- %% Generic case
+ %% Generic case for zoom
mylon_point=lon(Jzoom_start : Jzoom_end,Izoom_start : Izoom_end);
mylat_point=lat(Jzoom_start : Jzoom_end,Izoom_start : Izoom_end);
mymask_point=mask(Jzoom_start : Jzoom_end,Izoom_start : Izoom_end);
-
%%
- m_proj('mercator',... % dl/2 shift to have cells centered on lon/lat
- 'lon',[min(min(mylon_point-dl/2)) max(max(mylon_point-dl/2)) ],...
- 'lat',[min(min(mylat_point-dl/2)) max(max(mylat_point-dl/2)) ] );
+ m_proj('mercator',... % dlzoom/2 shift to have cells centered on lon/lat
+ 'lon',[min(min(mylon_point-dlzoom/2)) max(max(mylon_point-dlzoom/2)) ],...
+ 'lat',[min(min(mylat_point-dlzoom/2)) max(max(mylat_point-dlzoom/2)) ] );
- m_pcolor(mylon_point-dl/2,mylat_point-dl/2,mymask_point) ; shading flat
+ m_pcolor(mylon_point-dlzoom/2,mylat_point-dlzoom/2,mymask_point) ; %shading flat
m_grid('box','fancy','xtick',5,'ytick',5,'tickdir','out');
set(findobj('tag','m_grid_color'),'facecolor','white');
hold on
@@ -731,17 +678,31 @@ else
h11=plot(px0,py0,'ro');
[px,py]=m_ll2xy(lon_src,lat_src); % Adjusted position
if dir(k,1)==0 & dir(k,2)==1
- h33=plot(px,py,'k>');
- elseif dir(k,1)==0 & dir(k,2)==-1
- h33=plot(px,py,'k<');
- elseif dir(k,1)==1 & dir(k,2)==1
- h33=plot(px,py,'k^');
- elseif dir(k,1)==1 & dir(k,2)==-1
- h33=plot(px,py,'kv');
- end
+ h33=plot(px,py,'k>');
+ elseif dir(k,1)==0 & dir(k,2)==-1
+ h33=plot(px,py,'k<');
+ elseif dir(k,1)==1 & dir(k,2)==1
+ h33=plot(px,py,'k^');
+ elseif dir(k,1)==1 & dir(k,2)==-1
+ h33=plot(px,py,'kv');
+ end
legend([h11,h33],{'Approximative first guess river location','final adjusted river location'});
-
title(['\bf', myrivername(k,:)]);
+ %
+ % pcolor(mylon_point-dlzoom/2,mylat_point-dlzoom/2,mymask_point) ; %shading flat
+ % hold on
+ % h11=plot(lon_src0,lat_src0,'ro');
+ % if dir(k,1)==0 & dir(k,2)==1
+ % h33=plot(lon_src,lat_src,'k>');
+ % elseif dir(k,1)==0 & dir(k,2)==-1
+ % h33=plot(lon_src,lat_src,'k<');
+ % elseif dir(k,1)==1 & dir(k,2)==1
+ % h33=plot(lon_src,lat_src,'k^');
+ % elseif dir(k,1)==1 & dir(k,2)==-1
+ % h33=plot(lon_src,lat_src,'kv');
+ % end
+ % legend([h11,h33],{'Approximative first guess river location','final adjusted river location'});
+ % title(['\bf', myrivername(k,:)]);
end
end
end
@@ -754,17 +715,27 @@ else
%%nw{'runoff_position'}(:,:)=[I2(find(indomain_last==1)) J2(find(indomain_last==1))];
%%nw{'runoff_direction'}(:,:)=dir(find(indomain_last==1),:);
%%disp(['... river positions'])
- %% Write qbar, temp,salt and bgc variables conc.
- cff=1;
- %%
- %%debug cff=150;
- %%disp(['================== WARNING ============================'])
- %%disp(['Use a cff coef to increase the discharge= ',num2str(cff)])
- %%disp(['======================================================='])
- %% end debug
+ %%Write qbar, temp,salt and bgc variables conc.
+ cff = 1;
+ %
+ debugflow=0;
+ if debugflow
+ cff=10;
+ disp(['================== WARNING ============================'])
+ disp(['debugflow is on See line 633 in make_runoff.m '])
+ disp(['Use a cff coef to increase the discharge= ',num2str(cff)])
+ disp(['======================================================='])
+ end
my_flow=cff.*my_flow;
nw{'Qbar'}(:) = my_flow';
+ if debugflow
+ disp(['debugflow is on'])
+ disp(['SIZE MY_FLOW=',num2str(size(my_flow))]);
+ nw{'Qbar'}(4,:) = my_flow(:,1);
+ nw{'Qbar'}(5,:) = my_flow(:,1);
+ end
+
disp(['... discharges'])
if psource_ncfile_ts==1
%% take care : no transpostion needed !!
@@ -776,8 +747,6 @@ else
nw{'NO3_src'}(:) = my_no3_src';
disp(['... NO3 concentration'])
end
-
-
end
if psource_ncfile_ts == 1
disp([' ==>'])
@@ -829,8 +798,8 @@ else
T=20;
S=15;
end
- %% listing...
- disp([' ',num2str(I2(k)),' ',num2str(J2(k)),...
+ %Print the listing
+ disp([' ',num2str(I2for_out(k)),' ',num2str(J2for_out(k)),...
' ',num2str(dir(k,1)),' ',num2str(dir(k,2)),' T T ',...
num2str(T),' ',num2str(S)])
end
diff --git a/Rivers/runoff_glob_extract.m b/Rivers/runoff_glob_extract.m
index 1f05a15a60ff25524c76b0fedbdb050414288ce6..ca71602c81453104c081ed08b53931506c08d3e3 100644
--- a/Rivers/runoff_glob_extract.m
+++ b/Rivers/runoff_glob_extract.m
@@ -42,12 +42,9 @@ close(ncriv)
maxmaxflow=(max(max(FLOW_clm)));
maxflow=(max(FLOW_clm)'); % attention transpose
-meanflowmax=mean(FLOW_clm(:,1)); % pour amazon
-meanflow=nanmean(FLOW_clm(:,:),1)'; % pour les autres
-
-%size(meanflow)
-%plot(meanflow)
-%size(maxflow)
+meanflowmax=mean(FLOW_clm(:,1)); % for the biggest
+meanflow=nanmean(FLOW_clm(:,:),1)'; % for others
+%
%==========================================================================================================
%% => by default : all rivers in the domain
@@ -55,12 +52,18 @@ rivdetectype='DEFAULT';
my_riv=find(latriv_mou>=minlat & latriv_mou<=maxlat & lonriv_mou >= minlon & lonriv_mou <= maxlon);
%==
-
% => megatl : rivers in the boxes + max value >= 20%*max val of the flow amazonia peak)
%rivdetectype='MEGATL';
%thold=0.1;
%thold=0.03;
%my_riv=find(latriv_mou>=minlat & latriv_mou<=maxlat & lonriv_mou >= minlon & lonriv_mou <= maxlon & meanflow >= thold.*meanflowmax);
+
+% % => arvor : rivers in the boxes + max value >= 20%*max val of the flow orinoco peak)
+% rivdetectype='ARVOR';
+% %thold=0.1;
+% thold=0.03;
+% my_riv=find(latriv_mou>=minlat & latriv_mou<=maxlat & lonriv_mou >= minlon & lonriv_mou <= maxlon & meanflow >= thold.*meanflowmax);
+
%========================================
%
disp(['There are ',num2str(length(my_riv)),' rivers in the domain : '])