function oct_interp_OGCM(OGCM_dir,OGCM_prefix,year,month,Roa,interp_method,... lonU,latU,lonV,latV,lonT,latT,Z,tin,... nc_clm,ncid_bry,lon,lat,angle,h,tout,obc,vtransform) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Read the local OGCM files and perform the interpolations % Adapted for reducing computationnal time for bry file % by S. Illig, IRD-LEGOS % % % Further Information: % http://www.croco-ocean.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 % % Copyright (c) 2005-2006 by Pierrick Penven % e-mail:Pierrick.Penven@ird.fr % % Updated 6-Sep-2006 by Pierrick Penven : Nothing special for the bry file % Updated 5-Nov-2006 by Pierrick Penven : A bit of cleaning... % Update 13-Sep-2009 by Gildas Cambon : Begin treatments case for the bry % Update March-2011 by Serena Illig : optimisation bry % Update Nov-2011 by Pierrick Penven : cleaning % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % conserv=1; % same barotropic velocities as the OGCM % disp([' Horizontal interpolation: ',... OGCM_prefix,'Y',num2str(year),'M',num2str(month),'.cdf']) % % % CROCO grid angle % cosa=cos(angle); sina=sin(angle); % % Open the OGCM file % ncid = netcdf.open([OGCM_dir,OGCM_prefix,'Y',num2str(year),'M',num2str(month),'.cdf'], 'NC_NOWRITE'); % % Interpole data on the OGCM Z grid and CROCO horizontal grid % % Get zeta because it is needed to compute vertical levels of CROCO grid zeta=oct_ext_data_OGCM(ncid,lonT,latT,'ssh',tin,lon,lat,1,Roa,interp_method); % if ~isempty(nc_clm) % % Read and extrapole the 2D variables % u2d=oct_ext_data_OGCM(ncid,lonU,latU,'ubar',tin,lon,lat,1,Roa,interp_method); v2d=oct_ext_data_OGCM(ncid,lonV,latV,'vbar',tin,lon,lat,1,Roa,interp_method); ubar=oct_rho2u_2d(u2d.*cosa+v2d.*sina); vbar=oct_rho2v_2d(v2d.*cosa-u2d.*sina); % % Read and extrapole the 3D variables % NZ=length(Z); [M,L]=size(lon); dz=gradient(Z); temp=zeros(NZ,M,L); salt=zeros(NZ,M,L); u=zeros(NZ,M,L-1); v=zeros(NZ,M-1,L); for k=1:NZ if rem(k,10)==0 disp([' Level ',num2str(k),' of ',num2str(NZ)]) end u2d=oct_ext_data_OGCM(ncid,lonU,latU,'u',tin,lon,lat,k,Roa,interp_method); v2d=oct_ext_data_OGCM(ncid,lonV,latV,'v',tin,lon,lat,k,Roa,interp_method); u(k,:,:)=oct_rho2u_2d(u2d.*cosa+v2d.*sina); v(k,:,:)=oct_rho2v_2d(v2d.*cosa-u2d.*sina); temp(k,:,:)=oct_ext_data_OGCM(ncid,lonT,latT,'temp',tin,lon,lat,k,Roa,interp_method); salt(k,:,:)=oct_ext_data_OGCM(ncid,lonT,latT,'salt',tin,lon,lat,k,Roa,interp_method); end end % %Initialisation in case of bry files % if ~isempty(ncid_bry) NZ=length(Z); [M,L]=size(lon); % Read and extrapole the 2D variables if obc(1)==1 zeta_south=squeeze(zeta(1,:)); ubar1d=oct_ext_data_OGCM(ncid,lonU,latU,'ubar',tin,... squeeze(lon(1:2,:)),squeeze(lat(1:2,:)),1,Roa,interp_method); vbar1d=oct_ext_data_OGCM(ncid,lonV,latV,'vbar',tin,... squeeze(lon(1:2,:)),squeeze(lat(1:2,:)),1,Roa,interp_method); ubar_south=squeeze(oct_rho2u_2d(ubar1d(1,:).*cosa(1,:)+vbar1d(1,:).*sina(1,:))); vbar_south=squeeze(oct_rho2v_2d(vbar1d.*cosa(1:2,:)-ubar1d.*sina(1:2,:))); end if obc(2)==1 zeta_east=squeeze(zeta(:,end)); ubar1d=oct_ext_data_OGCM(ncid,lonU,latU,'ubar',tin,... squeeze(lon(:,end-1:end)),squeeze(lat(:,end-1:end)),1,Roa,interp_method); vbar1d=oct_ext_data_OGCM(ncid,lonV,latV,'vbar',tin,... squeeze(lon(:,end-1:end)),squeeze(lat(:,end-1:end)),1,Roa,interp_method); ubar_east=squeeze(oct_rho2u_2d(ubar1d.*cosa(:,end-1:end)+vbar1d.*sina(:,end-1:end))); vbar_east=squeeze(oct_rho2v_2d(vbar1d(:,end).*cosa(:,end)-ubar1d(:,end).*sina(:,end))); end if obc(3)==1 zeta_north=squeeze(zeta(end,:)); ubar1d=oct_ext_data_OGCM(ncid,lonU,latU,'ubar',tin,... squeeze(lon(end-1:end,:)),squeeze(lat(end-1:end,:)),1,Roa,interp_method); vbar1d=oct_ext_data_OGCM(ncid,lonV,latV,'vbar',tin,... squeeze(lon(end-1:end,:)),squeeze(lat(end-1:end,:)),1,Roa,interp_method); ubar_north=squeeze(oct_rho2u_2d(ubar1d(end,:).*cosa(end,:)+vbar1d(end,:).*sina(end,:))); vbar_north=squeeze(oct_rho2v_2d(vbar1d.*cosa(end-1:end,:)-ubar1d.*sina(end-1:end,:))); end if obc(4)==1 zeta_west=squeeze(zeta(:,1)); ubar1d=oct_ext_data_OGCM(ncid,lonU,latU,'ubar',tin,... squeeze(lon(:,1:2)),squeeze(lat(:,1:2)),1,Roa,interp_method); vbar1d=oct_ext_data_OGCM(ncid,lonV,latV,'vbar',tin,... squeeze(lon(:,1:2)),squeeze(lat(:,1:2)),1,Roa,interp_method); ubar_west=squeeze(oct_rho2u_2d(ubar1d.*cosa(:,1:2)+vbar1d.*sina(:,1:2))); vbar_west=squeeze(oct_rho2v_2d(vbar1d(:,1).*cosa(:,1)-ubar1d(:,1).*sina(:,1))); end % % Read and extrapole the 3D variables % if obc(1)==1 temp_south=zeros(NZ,L); salt_south=zeros(NZ,L); u_south=zeros(NZ,L-1); v_south=zeros(NZ,L); end if obc(2)==1 temp_east=zeros(NZ,M); salt_east=zeros(NZ,M); u_east=zeros(NZ,M); v_east=zeros(NZ,M-1); end if obc(3)==1 temp_north=zeros(NZ,L); salt_north=zeros(NZ,L); u_north=zeros(NZ,L-1); v_north=zeros(NZ,L); end if obc(4)==1 temp_west=zeros(NZ,M); salt_west=zeros(NZ,M); u_west=zeros(NZ,M); v_west=zeros(NZ,M-1); end % for k=1:NZ if rem(k,10)==0 disp([' Level bry ',num2str(k),' of ',num2str(NZ)]) end if obc(1)==1 % Southern boundary t1d=squeeze(oct_ext_data_OGCM(ncid,lonT,latT,'temp',tin,... squeeze(lon(1:2,:)),squeeze(lat(1:2,:)),... k,Roa,interp_method)); s1d=squeeze(oct_ext_data_OGCM(ncid,lonT,latT,'salt',tin,... squeeze(lon(1:2,:)),squeeze(lat(1:2,:)),... k,Roa,interp_method)); temp_south(k,:)=squeeze(t1d(1,:)); salt_south(k,:)=squeeze(s1d(1,:)); u1d=oct_ext_data_OGCM(ncid,lonU,latU,'u',tin,... squeeze(lon(1:2,:)),squeeze(lat(1:2,:)),... k,Roa,interp_method); v1d=oct_ext_data_OGCM(ncid,lonV,latV,'v',tin,... squeeze(lon(1:2,:)),squeeze(lat(1:2,:)),... k,Roa,interp_method); u_south(k,:)=squeeze(oct_rho2u_2d(u1d(1,:).*cosa(1,:)+... v1d(1,:).*sina(1,:))); v_south(k,:)=squeeze(oct_rho2v_2d(v1d.*cosa(1:2,:)-... u1d.*sina(1:2,:))); end if obc(2)==1 % Eastern boundary t1d=squeeze(oct_ext_data_OGCM(ncid,lonT,latT,'temp',tin,... squeeze(lon(:,end-1:end)),squeeze(lat(:,end-1:end)),... k,Roa,interp_method)); s1d=squeeze(oct_ext_data_OGCM(ncid,lonT,latT,'salt',tin,... squeeze(lon(:,end-1:end)),squeeze(lat(:,end-1:end)),... k,Roa,interp_method)); temp_east(k,:)=squeeze(t1d(:,end)'); salt_east(k,:)=squeeze(s1d(:,end)'); u1d=oct_ext_data_OGCM(ncid,lonU,latU,'u',tin,... squeeze(lon(:,end-1:end)),squeeze(lat(:,end-1:end)),... k,Roa,interp_method); v1d=oct_ext_data_OGCM(ncid,lonV,latV,'v',tin,... squeeze(lon(:,end-1:end)),squeeze(lat(:,end-1:end)),... k,Roa,interp_method); cff=squeeze(oct_rho2u_2d(u1d.*cosa(:,end-1:end)+... v1d.*sina(:,end-1:end))); u_east(k,:)=cff'; cff=squeeze(oct_rho2v_2d(v1d(:,end).*cosa(:,end)-... u1d(:,end).*sina(:,end))); v_east(k,:)=cff'; end if obc(3)==1 % Northern boundary t1d=squeeze(oct_ext_data_OGCM(ncid,lonT,latT,'temp',tin,... squeeze(lon(end-1:end,:)),squeeze(lat(end-1:end,:)),... k,Roa,interp_method)); s1d=squeeze(oct_ext_data_OGCM(ncid,lonT,latT,'salt',tin,... squeeze(lon(end-1:end,:)),squeeze(lat(end-1:end,:)),... k,Roa,interp_method)); temp_north(k,:)=squeeze(t1d(end,:)); salt_north(k,:)=squeeze(s1d(end,:)); u1d=oct_ext_data_OGCM(ncid,lonU,latU,'u',tin,... squeeze(lon(end-1:end,:)),squeeze(lat(end-1:end,:)),... k,Roa,interp_method); v1d=oct_ext_data_OGCM(ncid,lonV,latV,'v',tin,... squeeze(lon(end-1:end,:)),squeeze(lat(end-1:end,:)),... k,Roa,interp_method); u_north(k,:)=squeeze(oct_rho2u_2d(u1d(end,:).*cosa(end,:)+... v1d(end,:).*sina(end,:))); v_north(k,:)=squeeze(oct_rho2v_2d(v1d.*cosa(end-1:end,:)-... u1d.*sina(end-1:end,:))); end if obc(4)==1 % Western boundary t1d=squeeze(oct_ext_data_OGCM(ncid,lonT,latT,'temp',tin,... squeeze(lon(:,1:2)),squeeze(lat(:,1:2)),... k,Roa,interp_method)); s1d=squeeze(oct_ext_data_OGCM(ncid,lonT,latT,'salt',tin,... squeeze(lon(:,1:2)),squeeze(lat(:,1:2)),... k,Roa,interp_method)); temp_west(k,:)=squeeze(t1d(:,1)'); salt_west(k,:)=squeeze(s1d(:,1)'); u1d=oct_ext_data_OGCM(ncid,lonU,latU,'u',tin,... squeeze(lon(:,1:2)),squeeze(lat(:,1:2)),... k,Roa,interp_method); v1d=oct_ext_data_OGCM(ncid,lonV,latV,'v',tin,... squeeze(lon(:,1:2)),squeeze(lat(:,1:2)),... k,Roa,interp_method); cff=squeeze(oct_rho2u_2d(u1d.*cosa(:,1:2)+... v1d.*sina(:,1:2))); u_west(k,:)=cff'; cff=squeeze(oct_rho2v_2d(v1d(:,1).*cosa(:,1)-... u1d(:,1).*sina(:,1))); v_west(k,:)=cff'; end end end % % Close the OGCM file % netcdf.close(ncid); % % % Get the CROCO vertical grid % disp(' Vertical interpolations') if ~isempty(nc_clm) theta_s=nc_clm{'theta_s'}(:); theta_b=nc_clm{'theta_b'}(:); hc=nc_clm{'hc'}(:); N=length(nc_clm('s_rho')); end if ~isempty(ncid_bry) theta_s=netcdf.getVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'theta_s')); theta_b=netcdf.getVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'theta_b')); hc=netcdf.getVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'hc')); N=netcdf.inqDimLen(ncid_bry, netcdf.inqDimID(ncid_bry, 's_rho')); end % % Add an extra bottom layer (-100000m) and an extra surface layer (+100m) % to prevent vertical extrapolations % Z=[100;Z;-100000]; % % CROCO vertical grid % zr=oct_zlevs(h,zeta,theta_s,theta_b,hc,N,'r',vtransform); zu=oct_rho2u_3d(zr); zv=oct_rho2v_3d(zr); zw=oct_zlevs(h,zeta,theta_s,theta_b,hc,N,'w',vtransform); dzr=zw(2:end,:,:)-zw(1:end-1,:,:); dzu=oct_rho2u_3d(dzr); dzv=oct_rho2v_3d(dzr); % % % Vertical interpolation in case of clim file % % if ~isempty(nc_clm) % % Add a level on top and bottom with no-gradient % u=cat(1,u(1,:,:),u); u=cat(1,u,u(end,:,:)); v=cat(1,v(1,:,:),v); v=cat(1,v,v(end,:,:)); temp=cat(1,temp(1,:,:),temp); temp=cat(1,temp,temp(end,:,:)); salt=cat(1,salt,salt(end,:,:)); salt=cat(1,salt(1,:,:),salt); % % Perform the vertical interpolations % u=oct_ztosigma(flipdim(u,1),zu,flipud(Z)); v=oct_ztosigma(flipdim(v,1),zv,flipud(Z)); temp=oct_ztosigma(flipdim(temp,1),zr,flipud(Z)); salt=oct_ztosigma(flipdim(salt,1),zr,flipud(Z)); % % Correct the horizontal transport % i.e. remove the interpolated tranport and add % the OGCM transport % if conserv==1 u=u-oct_tridim(squeeze(sum(u.*dzu)./sum(dzu)),N); v=v-oct_tridim(squeeze(sum(v.*dzv)./sum(dzv)),N); u=u+oct_tridim(ubar,N); v=v+oct_tridim(vbar,N); end % % Barotropic velocities % ubar=squeeze(sum(u.*dzu)./sum(dzu)); vbar=squeeze(sum(v.*dzv)./sum(dzv)); % end %~isempty(nc_clm) % % % Vertical interpolation in case of bry files % % if ~isempty(ncid_bry) if obc(1)==1 [u_south,v_south,ubar_south,vbar_south,... temp_south,salt_south]=oct_vinterp_OGCM_bry(zr(:,1,:),zu(:,1,:),zv(:,1,:),... dzr(:,1,:),dzu(:,1,:),dzv(:,1,:),... u_south,v_south,ubar_south,vbar_south,... temp_south,salt_south,... N,Z,conserv); end if obc(2)==1 [u_east,v_east,ubar_east,vbar_east,... temp_east,salt_east]=oct_vinterp_OGCM_bry(zr(:,:,end),zu(:,:,end),zv(:,:,end),... dzr(:,:,end),dzu(:,:,end),dzv(:,:,end),... u_east,v_east,ubar_east,vbar_east,... temp_east,salt_east,... N,Z,conserv); end if obc(3)==1 [u_north,v_north,ubar_north,vbar_north,... temp_north,salt_north]=oct_vinterp_OGCM_bry(zr(:,end,:),zu(:,end,:),zv(:,end,:),... dzr(:,end,:),dzu(:,end,:),dzv(:,end,:),... u_north,v_north,ubar_north,vbar_north,... temp_north,salt_north,... N,Z,conserv); end if obc(4)==1 [u_west,v_west,ubar_west,vbar_west,... temp_west,salt_west]=oct_vinterp_OGCM_bry(zr(:,:,1),zu(:,:,1),zv(:,:,1),... dzr(:,:,1),dzu(:,:,1),dzv(:,:,1),... u_west,v_west,ubar_west,vbar_west,... temp_west,salt_west,... N,Z,conserv); end end %~isempty(ncid_bry) %-------------------------------------------------------------- % % fill the files % % Climatology file % if ~isempty(nc_clm) nc_clm{'zeta'}(tout,:,:)=zeta; nc_clm{'SSH'}(tout,:,:)=zeta; nc_clm{'temp'}(tout,:,:,:)=temp; nc_clm{'salt'}(tout,:,:,:)=salt; nc_clm{'u'}(tout,:,:,:)=u; nc_clm{'v'}(tout,:,:,:)=v; nc_clm{'ubar'}(tout,:,:,:)=ubar; nc_clm{'vbar'}(tout,:,:,:)=vbar; end % % Boundary file % if ~isempty(ncid_bry) if obc(1)==1 netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'zeta_south'), tout,:-1, 1, zeta_south); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'temp_south'), tout,:,:-1, 1, temp_south); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'salt_south'), tout,:,:-1, 1, salt_south); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'u_south'), tout,:,:-1, 1, u_south); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'v_south'), tout,:,:-1, 1, v_south); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'ubar_south'), tout,:,:-1, 1, ubar_south); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'vbar_south'), tout,:,:-1, 1, vbar_south); % [conv] 0-based end if obc(2)==1 netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'zeta_east'), tout,:-1, 1, zeta_east); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'temp_east'), tout,:,:-1, 1, temp_east); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'salt_east'), tout,:,:-1, 1, salt_east); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'u_east'), tout,:,:-1, 1, u_east); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'v_east'), tout,:,:-1, 1, v_east); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'ubar_east'), tout,:,:-1, 1, ubar_east); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'vbar_east'), tout,:,:-1, 1, vbar_east); % [conv] 0-based end if obc(3)==1 netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'zeta_north'), tout,:-1, 1, zeta_north); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'temp_north'), tout,:,:-1, 1, temp_north); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'salt_north'), tout,:,:-1, 1, salt_north); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'u_north'), tout,:,:-1, 1, u_north); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'v_north'), tout,:,:-1, 1, v_north); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'ubar_north'), tout,:,:-1, 1, ubar_north); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'vbar_north'), tout,:,:-1, 1, vbar_north); % [conv] 0-based end if obc(4)==1 netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'zeta_west'), tout,:-1, 1, zeta_west); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'temp_west'), tout,:,:-1, 1, temp_west); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'salt_west'), tout,:,:-1, 1, salt_west); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'u_west'), tout,:,:-1, 1, u_west); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'v_west'), tout,:,:-1, 1, v_west); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'ubar_west'), tout,:,:-1, 1, ubar_west); % [conv] 0-based netcdf.putVar(ncid_bry, netcdf.inqVarID(ncid_bry, 'vbar_west'), tout,:,:-1, 1, vbar_west); % [conv] 0-based end end end