Generate river's forcing

[vleguenn]

Following the suggestions from Sarah W. (and Anna K.) it was decided to implement rivers using the JRA55 datasets (available here on livljobs : _/projectsa/NEMO/slwa/globalrivers/JRA55-do-1-5-0). Each yearly file provide daily values of river flow in kg/m2/s.

[vleguenn]

1. Exploration of the JRA55 datasets for the river flow

The extraction of the Caribbean area from the Global river datasets (JRA55 in kg m-2 s-1) looks like this for the 1st January 2010:

Note that in the figure above, all the river flow values on the west coast of Mexico (in the Pacific side) are not of interest for our case. We will discard those values when generating river files. To obtain the flow in m3/s : area (m2) flow (kg m-2 s-1) 0.001

If we only show river flow with values above the threshold of 0.0015 kg m-2 s-1, we can highlight river mouths of some of the main rivers :

The river flow from JRA55 data spread alongside the entire coast and therefore does not show only one single value at the location of the river mouth. We can assess which are the main rivers represented inside JRA55 datasets. In total, I obtained a total of 13 rivers: Orinoco, Atrato, Magdalena, Mississipi, Mobile, Essequibo, Coppename, Corantijn, San Juan, Papaloapan, Grijalva, Panuco and Escondido. The river mouth of those 13 rivers is represented in the map below (alongside the bathymetry) :

[vleguenn]

2. Processing of the data

The raw data from JRA55 need to be processed so that it is in the correct format for NEMO inputs.

The script comes from Anna K., it is located in livljobs (_/projectsa/accord/SEAsiaR36/RIVERS) and it is called _rivers_JRA55spread.m. This Matlab script reads the river files, interpolates to the new grid by taking into account where is the coastline, then spread the rivers, and finally create a new netcdf file. The function _remaprivers.m just interpolates in the file.

Note that the script allows to spread the rivers and the threshold value used is arbitrary. It may need a lot of tries and errors to find a "good" value. If you do NOT want to spread the rivers, you have to comment all the lines inside %% spread large rivers to adjacent sea points.

The script _rivers_JRA55spread.m is : (Note here that here the threshold has the value of 2000)

clear;clc;close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Anna Katavouta NOC 07/2020
% create river forcing for NEMO from 
% MRI-JRA55 input MIPs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% choose year
year=2009;

%% read Seasia_R12 domain
cd C:\Users\User\Desktop\path
file='domain_cfg_Caribbean.nc';
lat_R12=ncread(file,'nav_lat');
lon_R12=ncread(file,'nav_lon');
e1t_R12=ncread(file,'e1t');
e2t_R12=ncread(file,'e2t');
Land_R12=single(ncread(file,'top_level'));Land_R12(Land_R12==0)=nan;

%% read MRI rivers
cd C:\Users\User\Desktop\path
file=['friver_input4MIPs_atmosphericState_OMIP_MRI-JRA55-do-1-5-0_gr_',num2str(year),'0101-',num2str(year),'1231.nc'];

ilat=380; nlat=108; ilon=1040; nlon=182;

x1=ncread(file,'lon',[ilon],[nlon]);
lat=ncread(file,'lat',[ilat],[nlat])';
lon=x1;lon(x1>180)=lon(x1>180)-360;
friver=ncread(file,'friver',[ilon ilat 1],[nlon nlat Inf]); % in kg m-2 s-1
% friver=squeeze(mean(friver,3));
 friver(friver==0)=NaN;

 % cell areas
cd C:\Users\User\Desktop\path
file='areacello_input4MIPs_atmosphericState_OMIP_MRI-JRA55-do-1-5-0_gr.nc';
a1=ncread(file,'areacello',[ilon ilat],[nlon nlat]);

friver=friver.*a1; % kg/s
temp=friver*0.001; % m3/s

mask=friver;mask(~isnan(mask))=1;
%% find coast in our domain - Jason way
Re=6367456.0*pi/180;
dr=mean([e1t_R12(:) ; e2t_R12(:)]);
dx=mean(e1t_R12(:))./(Re.*cos(mean(lon_R12(:)).*180/pi));
dy=mean(e2t_R12(:))./Re;

%get mask coast_R12
mask_coast= isnan(Land_R12);
mask_sea= Land_R12==1;

mask_coast = [mask_coast(1,:); mask_coast; mask_coast(end,:)] ;
mask_coast = [mask_coast(:,1), mask_coast, mask_coast(:,end)] ;

mask_coast = mask_coast(1:end-2,2:end-1) + mask_coast(3:end,2:end-1)      ...
          + mask_coast(2:end-1,1:end-2) + mask_coast(2:end-1,3:end) ;
mask_R12_C = (mask_coast>0 & +mask_sea==1) ;

% No rivers on boundaries
mask_R12_C(1,:)=0; mask_R12_C(end,:)=0;mask_R12_C(:,1)=0;mask_R12_C(:,end)=0;
mask_R12_C(2,:)=0; mask_R12_C(end-1,:)=0;mask_R12_C(:,2)=0;mask_R12_C(:,end-1)=0;

mask_R12=double(mask_R12_C);
mask_R12(isnan(Land_R12))=nan;

%% loop over time as different mask and sometimes rivers only in some dates
[LLAT LLON]=meshgrid(lat,lon);
%River=zeros(size(lon_R36,1),size(lon_R36,2),size(temp,3));

cd C:\Users\chepe\Desktop\from_sarah\Anna_script
for tt=1:size(temp,3)
    River_remap=remap_river(mask(:,:,tt),LLON,LLAT,mask_R12,lon_R12,lat_R12,temp(:,:,tt));
    River(:,:,tt)=River_remap;
    tt
end

%% spread large rivers to adjacent sea points
%spead it equally
River_uns=River;
mask_spread=mask_R12; mask_spread(~isnan(mask_spread))=1;
for tt=1:size(River,3)
    [Indi,Indj]=find(squeeze(River(:,:,tt))>2000);
    Rad_s=2;%choose how many points to spread intitially, with each iteration this 
    %radisu will grow until all large rivers have been spread

    while Indi>0
    for ii=1:length(Indi)
         Dummy2=mask_spread(Indi(ii)-Rad_s:Indi(ii)+Rad_s,Indj(ii)-Rad_s:Indj(ii)+Rad_s);
         np=nansum(nansum(Dummy2,2),1);
         [Di, Dj]=find(~isnan(Dummy2));
         Ri_SUM=nansum(nansum(River(Indi(ii)-Rad_s:Indi(ii)+Rad_s,Indj(ii)-Rad_s:Indj(ii)+Rad_s,tt),2),1);
         for jj=1:length(Di)
             River(Indi(ii)-(Rad_s+1)+Di(jj),Indj(ii)-(Rad_s+1)+Dj(jj),tt)=Ri_SUM/np;
         end
         %River(Indi(ii),Indj(ii),tt)=Ri_mid/np;
    end
    % repat until all the large rivers have been spread
    Rad_s=Rad_s+1;
    display([num2str(Rad_s),' with time ',num2str(tt)])
    [Indi,Indj]=find(squeeze(River(:,:,tt))>2000);
    % cap it to 45 Km radius
%     if Rad_s>=15
%         Indi=0;
%     end
    end

    % check that you included all river output
    if (nansum(squeeze(nansum(River_uns(:,:,tt),1)),2)-nansum(squeeze(nansum(River(:,:,tt),1)),2))>0.000001
        display (['Error total runoff in t=',num2str(tt)])
    end
end

%% change to Kg m^(-2) s^(-1) that nemo wants
rho0=1000;
for tt=1:size(River,3)
    rorunoff(:,:,tt)=River(:,:,tt).*rho0./(e1t_R12.*e2t_R12);
end

rorunoff(isnan(rorunoff))=0;
%% create depth information
DEPTH_R=10.*ones(size(River,1),size(River,2));%,size(River,3));

%% write data

cd C:\Users\User\Desktop\path

output=(['rivers_T2000_y',num2str(year),'.nc']);
nx=size(lon_R12,1);
ny=size(lon_R12,2);
tt=1:size(rorunoff,3);

nccreate(output,'time_counter','Dimensions',{'time_counter',Inf},'Format','netcdf4_classic');
nccreate(output,'nav_lat','Dimensions',{'x',nx,'y',ny},'Format','netcdf4_classic');
nccreate(output,'nav_lon','Dimensions',{'x',nx,'y',ny},'Format','netcdf4_classic');
nccreate(output,'rorunoff','Dimensions',{'x',nx,'y',ny,'time_counter',Inf},'Format','netcdf4_classic');
nccreate(output,'DEPTHrunoff','Dimensions',{'x',nx,'y',ny},'Format','netcdf4_classic');

ncwrite(output,'time_counter',tt);
ncwrite(output,'nav_lon',lon_R12);
ncwrite(output,'nav_lat',lat_R12);
ncwrite(output,'rorunoff',rorunoff);
ncwriteatt(output, 'rorunoff', 'long_name','river runoff' );
ncwriteatt(output, 'rorunoff','units','Kg m^-2 s^-1');
ncwrite(output,'DEPTHrunoff',DEPTH_R);

The script for the function _remaprivers.m is:

function River=remap_river(mask,LLON,LLAT,mask_R12,lon_R12,lat_R12,river_in)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Remap rivers
% Anna Katavouta, NOC, Liverpool 07/2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    River=zeros(size(lon_R12,1),size(lon_R12,2));
    %% move everything to our coast
    mask_log=(~isnan(mask));
    lon_log=LLON(mask_log);
    lat_log=LLAT(mask_log);

    %keep only what is in your domain;
    lat_log(lat_log>nanmax(lat_R12(:)))=nan;lat_log(lat_log<nanmin(lat_R12(:)))=nan;
    lon_log(isnan(lat_log))=nan;
    lon_log(lon_log<nanmin(lon_R12(:)))=nan;lon_log(lon_log>nanmax(lon_R12(:)))=nan;
    lat_log(isnan(lon_log))=nan;
    lon_save=lon_log;
    lat_log(isnan(lat_log))=[];
    lon_log(isnan(lon_log))=[];

    %find nearest neighbour on the R12 coast grid
    lon_coast=lon_R12;lat_coast=lat_R12;
    lon_coast(find(mask_R12~=1))=nan;lat_coast(find(mask_R12~=1))=nan;
    lon_coast(isnan(lon_coast))=[];lat_coast(isnan(lat_coast))=[];
    for ii=1:length(lon_log)
       dxf = [lon_log(ii), lat_log(ii)];
       G = [lon_coast(:),lat_coast(:)];
       d =  sum( (G-dxf).^2, 2);
       [minDist, idxMinDist] = min(d);
%        if minDist>0.4 %change it as you want for far from coast
%           CC(ii,1)=nan;CC(ii,2)=nan;
%           ind_ci(ii)=nan;
%           ind_cj(ii)=nan;
%        else
          solution = G(idxMinDist,:);
          CC(ii,:)=solution;
          [Ic Jc]=find(lon_R12==solution(1));
          ind_ci(ii)=Ic(1);
          [Ic Jc]=find(lat_R12==solution(2));
          ind_cj(ii)=Jc(1);
       %end
    end    

    %% recreate a 2x2 matrix on domain grid 
    %ensure that you add catchmens that go to the same point   
    R_log1=river_in; 
    R_log=R_log1(mask_log);
    R_log(isnan(lon_save))=[];
    for ii=1:length(ind_ci)
       if ~isnan(ind_ci(ii)) 
           River(ind_ci(ii),ind_cj(ii))=River(ind_ci(ii),ind_cj(ii))+(R_log(ii));
        end      
    end

    % check that you included all river output
    R_log(isnan(ind_ci))=nan;R_log(isnan(ind_cj))=nan;
    if (nansum(R_log)-nansum(squeeze(nansum(River,1)),2))>0.000001
       display (['Error total runoff in t=',num2str(tt)])
    end

end

[vleguenn]

3. Visualisation of the netcdf file generated

After using Anna K.'s script which prepare the data as input for NEMO grid, we obtain the following maps for the 1st of January 2013 using a) WITH SPREADING (treshold of 2000)

a) WITH SPREADING (threshold of 2000) Threshold of 2000 is used to spread flow to the adjacent cells.

Note on spreading : The lower the value of the threshold, the more it spreads. It spread around a radius (like equally along the coast and towards the ocean). it is better to always spread as a radius (so both coast and outcoast) rather than laterally only along the coast. The spreading is a more like trial and error, so doing some first test you will see from results if you need to do more spreading.

[vleguenn]

4. Location (path) of the outputs

• Latest run available with rivers : (2009 is the spin-up), 2010 to 2016.

• Year 2017 is currently running on Archer 2. So far, no errors linked with the implementation of the rivers.

The outputs have been transferred from Archer2 (_/work/n01/n01/valegu/CARIBBEAN_NEMO_4_0_4/nemo/cfgs/Caribbean/EXPFullForcings/) into livljobs (_/projectsa/CME/3DCaribbeanNEMO/OUTPUT/RivSPREADMSK/). The folder permission has been open to Marta (marpay@noc.ac.uk) and Gaby (gmaya@noc.ac.uk).

[vleguenn]

5. List of variables and frequencies available

Table 1 : Summary of the outputs. See Table 2 for the meaning of the variables.

Hourly means (extract only surface values)
Grid_T: tos, sos Grid_U: uos Grid_V: vos

25h means Grid_T: thetao, so, tos,tossq, sos, zos, zossq, wfo, rsntds, tohfls, taum, mldkz5, mldr10_1 Grid_U: uo, uos, tauuo Grid_V: vo, vos, tauvo Grid_W: wo, difvho

Daily means
Grid_T: thetao, so, tos,tossq, sos, zos, zossq, wfo, rsntds, tohfls, taum, mldkz5, mldr10_1 Grid_U: uo, uos, tauuo Grid_V: vo, vos, tauvo Grid_W: wo, difvho, difvmo

5 days means
Grid_T: thetao, so, tos,tossq, sos, zos, zossq, wfo, rsntds, tohfls, taum, mldkz5, mldr10_1 Grid_U: uo, uos, tauuo Grid_V: vo, vos, tauvo Grid_W: wo, difvho

Monthly means Grid_T: thetao, so, tos,tossq, sos, zos, zossq, wfo, rsntds, tohfls, taum, mldkz5, mldr10_1 Grid_U: uo, uos, tauuo Grid_V: vo, vos, tauvo Grid_W: wo, difvho

Yearly means
Grid_T: thetao, so, tos,tossq, sos, zos, zossq, wfo, rsntds, tohfls, taum, mldkz5, mldr10_1 Grid_U: uo, uos, tauuo Grid_V: vo, vos, tauvo Grid_W: wo, difvho

Table 2: Details of the variables used in Table 1.

_GridT thetao 3D temperature (°C) tos 2D Surface temperature (°C) tossq 2D Square of sea surface temperature (°C²) so 3D Salinity
sos 2D Surface salinity
zos 2D Surface height (m) zossq 2D Square of sea surface height (m2) wfo 2D Net upward water flux (kg/m2/s) rsntds 2D Shortwave radiation (w/m2) Tohfls 2D Net downward heat flux (w/m2) taum 2D Wind stress module (N/m2) Mldkz5 2D Turbocline depth (m) Mldr10_1 2D Mixed layer depth (m)

_GridU uo 3D Ocean current along I axis (m/s) uos 2D Ocean surface current along I axis (m/s) tauuo 2D Wind stress along I axis (N/m2)

_GridV vo 3D Ocean current along j axis (m/s) vos 2D Ocean surface current along j axis (m/s) tauvo 2D Wind stress along j axis (N/m2)

_GridW wo 3D Ocean vertical velocity (m/s) difvho 3D Vertical eddy diffusivity (m2/s) difvmo ? Ocean vertical momentum diffusivity (m2/s)

[vleguenn]

6. Initial results of a run with rivers

[vleguenn]

7. Comparison run with/without rivers

[vleguenn]

8. Next steps

• Release the Caribbean (+rivers) version on GitHub.
• Generate a Zenodo DOI