%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Create a CROCO run-off forcing file % % % Further Information: % http://www.crocoagrif.org % % This file is part of CROCOTOOLS % % CROCOTOOLS is free software; you can redistribute it and/or modify % it under the terms of the GNU General Public License as published % by the Free Software Foundation; either version 2 of the License, % or (at your option) any later version. % % CROCOTOOLS is distributed in the hope that it will be useful, but % WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program; if not, write to the Free Software % Foundation, Inc., 59 Temple Place, Suite 330, Boston, % MA 02111-1307 USA % % % July 2013: G. Cambon (IRD/LEGOS) & M. Herrmann (IRD/LEGOS) % July 2023 : G. Cambon (IRD/LOPS) % clear all close all %%%%%%%%%%%%%%%%%%% USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%%% %% crocotools_param %========================================================================================= % Choose the grid level into which you ant to set up the runoffs gridlevel=0; if ( gridlevel == 0 ) % #-> Parent / zoom #O grdname = [CROCO_files_dir,'croco_grd.nc']; rivname = [CROCO_files_dir,'croco_runoff.nc']; clmname = [CROCO_files_dir,'croco_clm.nc']; % <- climato file for runoff else % # -> Child / zoom #XX grdname = [CROCO_files_dir,'croco_grd.nc.',num2str(gridlevel)]; rivname = [CROCO_files_dir,'croco_runoff.nc.',num2str(gridlevel)]; clmname = [CROCO_files_dir,'croco_clm.nc.',num2str(gridlevel)]; % <- climato file for runoff end % Choose the monthly runoff forcing time and cycle in days clim_run=1; if (clim_run == 1) qbar_time=[15:30:345]; qbar_cycle=360; else qbar_time=[15.2188:30.4375:350.0313]; qbar_cycle=365.25; end % - times and cycles for runoff conditions: % - clim_run = 1 % climato forcing experiments with climato calendar % qbar_time=[15:30:365]; % qbar_cycle=360; % % - clim_run = 0 % interanual forcing experiments with real calendar % qbar_time=[15.2188:30.4375:350.0313]; % qbar_cycle=365.25; %========================================================================================= % Choose if you process variable tracer concentration(temp, salt, NO3, ...) psource_ncfile_ts=1; if psource_ncfile_ts psource_ncfile_ts_auto=1 ; psource_ncfile_ts_manual=0; end % - psource_ncfile_ts = 0 => Constant analytical runoff tracers concentration no processing % It reads analytical values in croco.in % or use default value defined in % analytical.F % % - pource_ncfile_ts = 1 => Variable runoff tracers % concentration processing is activated. % It needs the climatology % file created with oct_make_clim.m % % In this case, either choose: % - psource_ts_auto : auto definition using % the nearest point in the climatlogy file % % - psource_ts_manual : manually definition the % variable tracer concentration %========================================================================================= % Fancy plots plotting=1; plotting_zoom=0; % %========================================================================================= % Add biogeochemical variables if (makenpzd | makepisces | makebioebus | makequota) makebio = 1; else makebio = 0; end %% disp(' ') disp(['Create runoff forcing from Dai and Trenberth''s global monthly climatological run-off dataset']) disp(' ') title_name='runoff forcing file (Dai and Trenberth, 2002 dataset)'; %========================================================================================= define_dir=0 ; %%->flag to define directly the orientation / direction of the runoff % if define_dir==1 %% Define orientation/direction of the flow. First column is the u- (0) or v- (1) %% orientation. Second column is the direction of the flow in the choosen orientation % $$$ %%EXAMPLE GIGATl FAMILY for river # 1 to 9 without #7 and #8 % $$$ dir(1,:) = [0 , 1]; % # Amazon % $$$ dir(2,:) = [0 , -1]; % # Congo % $$$ dir(3,:) = [0 , 1]; % # Orinoco % $$$ dir(4,:) = [1 , -1]; % # Mississipi % $$$ dir(5,:) = [0 , 1]; % # Parana % $$$ dir(6,:) = [1 , 1]; % # Tocantins % $$$ dir(7,:) = [0 , 1]; % # Tapajos same dir/sense as Amazon % $$$ dir(8,:) = [0 , 1]; % # Saint Lawrence % $$$ dir(9,:) = [0 , 1]; % # Xingu same dir/sense as Amazon % $$$ 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 %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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % [latriv,lonriv,my_flow,myrivername,rivernumber]=oct_runoff_glob_extract(grdname,global_clim_rivername); %% => rivernumber is the total number of river in the domain if rivernumber == 0 %%at least a river disp(['Create a "fictive" runoff forcing file : river with no discharge']) oct_create_runoff(rivname,grdname,title_name,... qbar_time,qbar_cycle, ... 'fictiveriver_@nest',1,18,1,psource_ncfile_ts,makebio,makepisces,makequota) my_flow=0; nw_id = netcdf.open(rivname, 'NC_WRITE'); disp(['Write in runoff file']) netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'Qbar'), my_flow'); if psource_ncfile_ts==1 netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'temp_src'), my_temp_src'); netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'salt_src'), my_salt_src'); if makebio netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'NO3_src'), my_no3_src'); if makepisces netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'PO4_src'), my_po4_src'); netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'Si_src'), my_sil_src'); netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'DIC_src'), my_dic_src'); netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'DOC_src'), my_doc_src'); netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'TALK_src'), my_talk_src'); if makequota netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'DON_src'), my_don_src'); netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'DOP_src'), my_dop_src'); end end end end netcdf.close(nw_id); disp([' ']) disp(['Line to enter in the croco.in file in the psource_ncfile section :']) disp(['-----------------------------------------------------------------']) disp(['psource_ncfile: Nsrc Isrc Jsrc Dsrc qbardir Lsrc Tsrc runoff file name']) disp([' CROCO_FILES/croco_runoff.nc(.#nestlevel)']) disp([' ',num2str(0)']) %% else latriv=latriv'; lonriv=lonriv'; %% rivername=strvcat(myrivername); rivernumber=size(rivername,1); rivname_StrLen=size(rivername,2); %% %% Determine the positions of the river, from its lon/lat position %% extract j and i to put in croco.in / croco.in.1 for that use of croco_grd.nc %% %% Read the grid %% [lat,lon,mask]=read_latlonmask(grdname,'r'); [latu,lonu,masku]=read_latlonmask(grdname,'u'); [latv,lonv,maskv]=read_latlonmask(grdname,'v'); lonmin=min(lon(:)); lonmax=max(lon(:)); latmin=min(lat(:)); latmax=max(lat(:)); %% %% Perform a first pass on each river %% disp(' ') disp(['First guess:']) disp(['============']) % for k= 1:rivernumber disp([' ']) disp(['- Process river #',num2str(k),': ',char(myrivername(k,:))]) indomain(k)=oct_check_domain_runoff(lon,lat,lonriv(k),latriv(k)); [j,i]=oct_runoff_grid_pos(lon,lat,lonriv(k),latriv(k)); J(k)=j-1; 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)))]) 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... %% rivertoprocess=find(indomain==1); %%================= %%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 'lon',[lonmin-dl/2 lonmax-dl/2],... 'lat',[latmin-dl/2 latmax-dl/2]); m_pcolor(lon-dl/2,lat-dl/2,mask) shading flat m_grid('box','fancy','xtick',5,'ytick',5,'tickdir','out'); set(findobj('tag','m_grid_color'),'facecolor','white'); hold on for k=1:number_rivertoprocess lon_src=lon(J(k),I(k)); lat_src=lat(J(k),I(k)); [px,py]=m_ll2xy(lon_src,lat_src); h1=plot(px,py,'ro'); set(h1,'Clipping','off') legend(h1,'first guess position'); h2=m_text(lon(J(k),I(k)),lat(J(k),I(k))+0.1,myrivername(k,:)); set(h2,'fontweight','demi','fontsize',13); end %pause %% if plotting_zoom %% 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; figure(nfig) lon_src=lon(J(k),I(k)); lat_src=lat(J(k),I(k)); 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 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); % %% -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 %% %% Choose which river you really want to process... %% indomain_last=indomain-indomain; for k0=1:number_rivertoprocess k=rivertoprocess(k0); indomain_last(k)=input(['Do you want to use river (Yes[1], No[0]) ? ', rivername(k0,:)]); if indomain_last(k)==0 && plotting_zoom % Close the (zoom) figure if the river is not used close(nfig0+k0); end end rivertoprocess=find(indomain_last==1); if isempty(rivertoprocess) error(['No river selected !']) end number_rivertoprocess=length(rivertoprocess); rivername=strvcat(myrivername(rivertoprocess,:)); rivname_StrLen=size(rivername,2); my_flow=my_flow(:,find(indomain_last==1)); %% %% Define the orientation/direction of the flow. First column is the u- (0) or v- (1) %% orientation. Second column is the direction of the flow in the choosen orientation %% if define_dir==0 for k0=1:number_rivertoprocess k=rivertoprocess(k0); disp(['=====']) disp(['River ',rivername(k0,:)]) disp(['Choose the orientation of the flow']) dir11=NaN; while ~(dir11==0 | dir11==1) dir11 = input('0=zonal or 1=meridional. '); end disp(['Choose the direction of the flow. ']) dir12=NaN; while ~(dir12==1 | dir12==-1) dir12= input('1 is positive [S-N or W-E], -1 negative [N-S or E-W]. '); end dir(k,:)=[dir11 dir12]; disp(['k=',num2str(k)]) disp(['k0=',num2str(k0)]) disp(['dir(k,:)=',num2str(dir(k,:))]) end end %% %% Create the runoff forcing file %% disp(' ') disp(' Create the runoff file...') oct_create_runoff(rivname,grdname,title_name,... qbar_time,qbar_cycle, ... rivername,number_rivertoprocess,... rivname_StrLen,dir,psource_ncfile_ts,... makebio,makepisces,makequota) %% %% Adjust the rivers positions relative to the mask %% disp(['Find the real positions of the rivers in the grid: ']) disp(['===================================================']) for k0=1:number_rivertoprocess k=rivertoprocess(k0); disp(['Process final position for river ',myrivername(k,:)]) disp(['Choose the orientation']) jj=J(k); ii=I(k); dir2=dir(k,:); % % Compute the u and v cells indexes in CROCO convention [jj2for_out,ii2for_out, jj2for,ii2for,jj2,ii2]=oct_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,:)),' 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([' 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)); if makebio==1 my_no3_src0=zeros(rivernumber,length(woa_time)); if makepisces==1 my_po4_src0=zeros(rivernumber,length(woa_time)); my_sil_src0=zeros(rivernumber,length(woa_time)); my_dic_src0=zeros(rivernumber,length(woa_time)); my_doc_src0=zeros(rivernumber,length(woa_time)); my_talk_src0=zeros(rivernumber,length(woa_time)); if makequota==1 my_don_src0=zeros(rivernumber,length(woa_time)); my_dop_src0=zeros(rivernumber,length(woa_time)); end end end if psource_ncfile_ts_auto %%============================================================== ncidclim = netcdf.open(clmname, 'NC_NOWRITE'); N=netcdf.inqDimLen(ncidclim, netcdf.inqDimID(ncidclim, 's_rho')); for k0=1:number_rivertoprocess k=rivertoprocess(k0); %% %% For temperature, use the target surface rho_point in the clim file %% to reduce any heat flux induced by the rivers %% if dir(k,1)==0 & dir(k,2)==1 %Zonal flux W->E i_rho=I2(k); j_rho=J2(k); elseif dir(k,1)==0 & dir(k,2)==-1 %Zonal flux E->W i_rho=I2(k)-1; j_rho=J2(k); elseif dir(k,1)==1 & dir(k,2)==1 %Meridional flux S->N i_rho=I2(k); j_rho=J2(k); elseif dir(k,1)==1 & dir(k,2)==-1 %Meridional flux N->S i_rho=I2(k); j_rho=J2(k)-1; end T=squeeze(netcdf.getVar(ncidclim, netcdf.inqVarID(ncidclim, 'temp'))); my_temp_src0(k,:)=T'; %% %% For salinity ... ? %% S=squeeze(netcdf.getVar(ncidclim, netcdf.inqVarID(ncidclim, 'salt')))-10; % hum... S(S<2)=2.0001; % to prevent negative salinities in the % equation of state disp([' Use psource_ncfile_ts_auto using S = sclim -10 ']) disp([' Check line 464 in oct_make_runoff.m to change ' ... 'this arbitrary runoff salinity']) %%S=2; my_salt_src0(k,:)=S'; if makebio==1 NO3=squeeze(netcdf.getVar(ncidclim, netcdf.inqVarID(ncidclim, 'NO3'))); my_no3_src0(k,:)=NO3'; if makepisces==1 PO4=squeeze(netcdf.getVar(ncidclim, netcdf.inqVarID(ncidclim, 'PO4'))); my_po4_src0(k,:)=PO4'; Si=squeeze(netcdf.getVar(ncidclim, netcdf.inqVarID(ncidclim, 'Si'))); my_sil_src0(k,:)=Si'; DIC=squeeze(netcdf.getVar(ncidclim, netcdf.inqVarID(ncidclim, 'DIC'))); my_dic_src0(k,:)=DIC'; DOC=squeeze(netcdf.getVar(ncidclim, netcdf.inqVarID(ncidclim, 'DOC'))); my_doc_src0(k,:)=DOC'; TALK=squeeze(netcdf.getVar(ncidclim, netcdf.inqVarID(ncidclim, 'TALK'))); my_talk_src0(k,:)=TALK'; if makequota==1 DON=squeeze(netcdf.getVar(ncidclim, netcdf.inqVarID(ncidclim, 'DON'))); my_don_src0(k,:)=DON'; DOP=squeeze(netcdf.getVar(ncidclim, netcdf.inqVarID(ncidclim, 'DOP'))); my_dop_src0(k,:)=DOP'; end end end end netcdf.close(ncidclim); %%============================================================== elseif psource_ncfile_ts_manual %% Alternativaly : Define all mannually the tracer %% t, s, and eventually biogeochemical tracer concentration %% my_temp/salt/_src0 : line are rivers / rows are month %% size(my_temp_src0)=[number_rivertoprocess 12] %% my_temp/salt/_src : line are rivers effectivelly processed / rows are month %% Example for GIGATL FAMILY my_temp_src0(:,:)=[28.1677 27.7407 27.7290 27.9662 28.2562 ... 28.3312 28.3832 28.4892 28.8444 28.9035 ... 28.8618 28.6118 ; 27.4397 28.2505 28.7115 28.3359 26.6957 ... 24.3065 22.5866 22.1818 23.2263 24.9318 ... 26.2492 26.8046 ; 26.9818 27.0154 27.2094 27.5274 27.7536 ... 27.8428 28.1083 28.5623 28.7645 28.6072 ... 28.1333 27.3927 ; 18.4467 18.1003 19.3736 21.6955 24.9263 27.6829 ... 29.2091 29.5412 28.3085 25.5634 22.6769 ... 20.2277 ; 22.6310 23.2607 22.4562 20.2833 17.2070 13.8076 ... 11.6466 11.2372 12.5740 15.1383 18.1340 ... 20.6393 ; 27.9732 27.7875 27.9070 28.2647 28.5391 28.5458 ... 28.3633 28.1903 28.2218 28.2947 28.3242 ... 28.1359 ; 28.1677 27.7407 27.7290 27.9662 28.2562 28.3312 ... 28.3832 28.4892 28.8444 28.9035 28.8618 ... 28.6118 ; -0.3210 -0.5013 -0.1819 1.5011 5.6784 10.6697 ... 14.2692 14.6227 11.4167 7.0775 3.4114 ... 1.1195 ; 28.1677 27.7407 27.7290 27.9662 28.2562 28.3312 ... 28.3832 28.4892 28.8444 28.9035 28.8618 ... 28.6118 ; 26.9245 26.7369 26.8849 27.2396 27.8931 28.3006 ... 28.1773 28.4447 28.7895 28.8083 28.3689 ... 27.6033 ; 22.6310 23.2607 22.4562 20.2833 17.2070 13.8076 ... 11.6466 11.2372 12.5740 15.1383 18.1340 ... 20.6393 ; 9 9 9 9 9 9 9 9 9 9 9 9 ; %% <= not used 28.6263 28.9582 29.3015 29.4045 29.1542 28.3245 ... 27.3215 26.8242 26.8997 27.4027 28.0562 ... 28.4404 ]; % $$$ my_salt_src0(:,:)=[ 31.8610 30.2483 28.0273 25.8570 ... % $$$ 25.5184 26.1775 27.5553 28.7191 29.7456 ... % $$$ 30.7106 31.5538 32.1978 ; % $$$ % $$$ 30.1895 30.2550 30.4903 30.8359 32.1160 33.4259 ... % $$$ 33.8053 33.5998 32.8721 32.1840 31.4867 ... % $$$ 30.8132 ; % $$$ % $$$ 33.1099 33.8338 34.3348 34.4912 34.0427 33.1361 ... % $$$ 31.6499 30.3496 30.2760 31.2658 32.1016 ... % $$$ 32.5288 ; % $$$ % $$$ 31.4802 30.8317 29.5911 28.0779 27.4586 27.6335 ... % $$$ 28.4153 29.0401 29.9319 31.0356 31.6389 ... % $$$ 31.6811 ; % $$$ % $$$ 28.3315 28.3112 28.3833 28.5506 28.9386 29.3127 ... % $$$ 29.1937 28.8813 28.5683 28.3198 28.2123 ... % $$$ 28.2293 ; % $$$ % $$$ 30.6917 28.2332 25.6524 23.7380 24.2345 26.3038 ... % $$$ 28.2613 29.4522 30.6267 31.8787 32.6207 ... % $$$ 32.4137 ; % $$$ % $$$ 31.8610 30.2483 28.0273 25.8570 25.5184 26.1775 ... % $$$ 27.5553 28.7191 29.7456 30.7106 31.5538 ... % $$$ 32.1978 ; % $$$ % $$$ 30.5800 31.1429 31.3964 30.9931 29.8125 28.4991 ... % $$$ 27.7825 27.6972 27.8768 28.8042 29.5877 ... % $$$ 30.1445 ; % $$$ % $$$ 31.8610 30.2483 28.0273 25.8570 25.5184 26.1775 ... % $$$ 27.5553 28.7191 29.7456 30.7106 31.5538 ... % $$$ 32.1978 ; % $$$ % $$$ 35.5176 35.7714 35.9853 35.8987 35.6900 35.5734 ... % $$$ 35.7705 35.4552 35.1000 34.6942 34.7754 ... % $$$ 35.1856 ; % $$$ % $$$ 28.3315 28.3112 28.3833 28.5506 28.9386 29.3127 ... % $$$ 29.1937 28.8813 28.5683 28.3198 28.2123 ... % $$$ 28.2293 ; % $$$ % $$$ 33.1350 33.2516 33.2083 33.0608 32.9854 33.0303 ... % $$$ 33.3633 33.4677 33.3594 33.2315 33.0206 ... % $$$ 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 %% my_temp_src= my_temp_src0(find(indomain_last==1),:); my_salt_src= my_salt_src0(find(indomain_last==1),:); if makebio my_no3_src = my_no3_src0(find(indomain_last==1),:); if makepisces my_po4_src = my_po4_src0(find(indomain_last==1),:); my_sil_src = my_sil_src0(find(indomain_last==1),:); my_dic_src = my_dic_src0(find(indomain_last==1),:); my_doc_src = my_doc_src0(find(indomain_last==1),:); my_talk_src = my_talk_src0(find(indomain_last==1),:); if makequota my_don_src = my_don_src0(find(indomain_last==1),:); my_dop_src = my_dop_src0(find(indomain_last==1),:); end end end end %% %%============================================================== %% Continue the figure if plotting==1 figure(1) hold on m_proj('mercator',... % dl/2 shift to have cells centered on lon/lat 'lon',[lonmin-dl/2 lonmax-dl/2],... '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),I2for_out(k)); % Adjusted position lat_src=latu(J2(k),I2for_out(k)); % + nice plot else %meridional flow (v-grid) 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 h3=plot(px,py,'k>'); elseif dir(k,1)==0 & dir(k,2)==-1 h3=plot(px,py,'k<'); elseif dir(k,1)==1 & dir(k,2)==1 h3=plot(px,py,'k^'); elseif dir(k,1)==1 & dir(k,2)==-1 h3=plot(px,py,'kv'); end set(h3,'Clipping','off') end legend([h1,h3],{'Approximative first guess river location','final adjusted river location'}); title({'\bf Location of river in the croco grid';'(from Dai and Trenberth dataset)'}); if plotting_zoom nfig0=5; for k0=1:number_rivertoprocess k=rivertoprocess(k0); nfig=nfig0+k; 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_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(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 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',... % 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-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 [px0,py0]=m_ll2xy(lon_src0,lat_src0); % First guess position 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 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 %% Fill the river discharge and eventually %% t/s concentration, no3 concentration nw_id = netcdf.open(rivname, 'NC_WRITE'); disp(['Write in runoff file']) %%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; % debugflow=0; if debugflow cff=100; disp(['================== WARNING ============================']) disp(['debugflow is on See line 633 in oct_make_runoff.m ']) disp(['Use a cff coef to increase the discharge= ',num2str(cff)]) disp(['=======================================================']) end my_flow=cff.*my_flow; netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'Qbar'), my_flow'); disp(['... discharges']) if psource_ncfile_ts==1 %% take care : no transpostion needed !! netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'temp_src'), my_temp_src); disp(['... temperature concentration']) netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'salt_src'), my_salt_src); disp(['... salt concentration']) if makebio netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'NO3_src'), my_no3_src'); disp(['... NO3 concentration']) if makepisces netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'PO4_src'), my_po4_src'); disp(['... PO4 concentration']) netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'Si_src'), my_sil_src'); disp(['... Si concentration']) netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'DIC_src'), my_dic_src'); disp(['... DIC concentration']) netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'DOC_src'), my_doc_src'); disp(['... DOC concentration']) netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'TALK_src'), my_talk_src'); disp(['... TALK concentration']) if makequota netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'DON_src'), my_don_src'); disp(['... DON concentration']) netcdf.putVar(nw_id, netcdf.inqVarID(nw_id, 'DOP_src'), my_dop_src'); disp(['... DOP concentration']) end end end end if psource_ncfile_ts == 1 disp([' ==>']) disp([' PSOURCE_NCFILE_TS = 1']) disp([' Variable river discharge in m3/s + variable tracer ' ... '(temp/salt) concentration ']) if psource_ncfile_ts_auto disp([' PSOURCE_NCFILE_TS_AUTO = 1']) disp([' auto definition of the variable runoff tracer ']) disp([' concentration using the nearest point']) disp([' in the climatlogy file']) else disp([' PSOURCE_NCFILE_TS_MANUAL = 1']) disp([' manual definition of the variable runoff tracer']) disp([' concentration (see example in oct_make_runoff.m directly)']) end disp([' ...Note : ']) disp([' ... The Tsrc value reported in croco.in are the annual-mean tracer value']) disp([' ... It''s just for information !']) disp([' ... The Tsrc used are read in the runoff netCDF ' ... 'file created']) else disp([' ==> ']) disp([' PSOURCE_NCFILE_TS = 0']) disp([' Variable river discharge in m3/s + constant tracer (temp/salt) concentration ' ... 'fluxes']) disp([' ... The Tsrc value reported in croco.in are the ' ... 'constant value imposed manually']) end netcdf.close(nw_id); %% %% Line to enter in the croco.in file in the psource section %% disp([' ']) disp(['Line to enter in the croco.in file in the psource_ncfile section :']) disp(['-----------------------------------------------------------------']) disp(['psource_ncfile: Nsrc Isrc Jsrc Dsrc qbardir Lsrc Tsrc runoff file name']) disp([' CROCO_FILES/croco_runoff.nc(.#nestlevel)']) disp([' ',num2str(number_rivertoprocess)]) for k0=1:number_rivertoprocess k=rivertoprocess(k0); if psource_ncfile_ts==1 T=mean(my_temp_src0(k,:)); S=mean(my_salt_src0(k,:)); if makebio N=mean(my_no3_src0(k,:)); if makepisces P=mean(my_po4_src0(k,:)); S=mean(my_sil_src0(k,:)); A=mean(my_talk_src0(k,:)); Di=mean(my_dic_src0(k,:)); Dc=mean(my_doc_src0(k,:)); if makequota Dn=mean(my_don_src0(k,:)); Dp=mean(my_dop_src0(k,:)); end end end else T=20; S=15; if makebio N=10; if makepisces P=1; S=10; A=2300; Di=2100; Dc=28; if makequota Dn=10; Dp=1; end end end end %Print the listing if makebio==0 disp([' ',num2str(I2for_out(k)),' ',num2str(J2for_out(k)),... ' ',num2str(dir(k,1)),' ',num2str(dir(k,2)),' T T ',... num2str(T),' ',num2str(S)]) else if makepisces==0 disp([' ',num2str(I2for_out(k)),' ',num2str(J2for_out(k)),... ' ',num2str(dir(k,1)),' ',num2str(dir(k,2)),' T T T ',... num2str(T),' ',num2str(S),' ',num2str(N)]) else if makequota==0 disp([' ',num2str(I2for_out(k)),' ',num2str(J2for_out(k)),... ' ',num2str(dir(k,1)),' ',num2str(dir(k,2)),' 4*T 2*F T F T 2*F T 12*F T 2*F']) else disp([' ',num2str(I2for_out(k)),' ',num2str(J2for_out(k)),... ' ',num2str(dir(k,1)),' ',num2str(dir(k,2)),' 4*T 2*F T F T 2*F T 12*F T 2*F 2*T 13*F ']) end end end end disp(['-----------------------------------------------------------------']) %% %% Plot the seasonal cycle %% figure(30) if psource_ncfile_ts==1 subplot(3,2,1) end hold on plot([1:12],my_flow) legend(rivername) box on, grid on title(['\bf Monthly clim of the domain run off']) xlabel(['\bf Month']);ylabel(['\bf Discharge in m3/s']) set(gca,'Xtick',[0.5:11.5],'XtickLabel',['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D']); if psource_ncfile_ts==1 subplot(3,2,3) plot([1:12],my_temp_src) box on, grid on xlabel(['\bf Month']);ylabel(['\bf Temp [C]']) set(gca,'Xtick',[0.5:11.5],'XtickLabel',['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D']); % subplot(3,2,5) plot([1:12],my_salt_src) box on, grid on xlabel(['\bf Month']);ylabel(['\bf Salt [psu]']) set(gca,'Xtick',[0.5:11.5],'XtickLabel',['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D']); % if makebio subplot(3,2,2) plot([1:12],my_no3_src) box on, grid on xlabel(['\bf Month']);ylabel(['\bf NO3 [mmol/m3]']) set(gca,'Xtick',[0.5:11.5],'XtickLabel',['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D']); % if makepisces subplot(3,2,4) plot([1:12],my_doc_src) box on, grid on xlabel(['\bf Month']);ylabel(['\bf DOC [mmol/m3]']) set(gca,'Xtick',[0.5:11.5],'XtickLabel',['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D']); % subplot(3,2,6) plot([1:12],my_sil_src) box on, grid on xlabel(['\bf Month']);ylabel(['\bf Si [mmol/m3]']) set(gca,'Xtick',[0.5:11.5],'XtickLabel',['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D']); end end end end %% end at least one river !