roms preprocessing tools (d_mercator2roms.m)

Message

corvianawatie · Tue Apr 23, 2019 3:23 am

%
%  D_MERCATOR2ROMS:  Driver script to create a ROMS initial conditions
%
%  This a user modifiable script that can be used to prepare ROMS
%  initial conditions NetCDF file from Mercator dataset. It sets-up
%  all the necessary parameters and variables. USERS can use this
%  as a prototype for their application.
%
% svn $Id: d_mercator2roms.m 938 2019-01-28 06:35:10Z arango $
%=========================================================================%
%  Copyright (c) 2002-2019 The ROMS/TOMS Group                            %
%    Licensed under a MIT/X style license                                 %
%    See License_ROMS.txt                           Hernan G. Arango      %
%=========================================================================%
clear all;
clc
%% Set file names.
OPR_Dir = 'E:/matlab/initial/data_source/';
RTR_Dir = 'E:/matlab/initial/data_source/';
Tfile   = 'woa18_ann_temp.nc';
Sfile   = 'woa18_ann_salt.nc';
GRDname = 'E:/matlab/initial/data_source/roms_grd.nc';
INIname = 'E:/matlab/initial/data_source/roms_ini.nc';
CREATE = true;                   % logical switch to create NetCDF
report = false;                  % report vertical grid information
%% Get number of grid points.
[Lr,Mr]=size(nc_read(GRDname,'h'));
Lu = Lr-1;   Lv = Lr;
Mu = Mr;     Mv = Mr-1;
%%
%--------------------------------------------------------------------------
%  Create initial conditions NetCDF file.
%--------------------------------------------------------------------------
% Set full path of Mercator files assigned as initial conditions.
fileT=fullfile(RTR_Dir,Tfile); lenT=1;
fileS=fullfile(RTR_Dir,Sfile); lenS=1;
%%
%--------------------------------------------------------------------------
%  Set application parameters in structure array, S.
%--------------------------------------------------------------------------
S.ncname      = INIname;     % output NetCDF file
S.spherical   = 1;           % spherical grid
S.Lm          = Lr-2;        % number of interior RHO-points, X-direction
S.Mm          = Mr-2;        % number of interior RHO-points, Y-direction
S.N           = 50;          % number of vertical levels at RHO-points
S.NT          = 2;           % total number of tracers
S.Vtransform  = 2;           % vertical transfomation equation
S.Vstretching = 4;           % vertical stretching function
S.theta_s     = 7.0;         % S-coordinate surface control parameter
S.theta_b     = 0.1;         % S-coordinate bottom control parameter
S.Tcline      = 200.0;       % S-coordinate surface/bottom stretching width
S.hc          = S.Tcline;    % S-coordinate stretching width
%%
%--------------------------------------------------------------------------
%  Set grid variables.
%--------------------------------------------------------------------------
S.h           = nc_read(GRDname, 'h');            % bathymetry
S.lon_rho     = nc_read(GRDname, 'lon_rho');      % RHO-longitude
S.lat_rho     = nc_read(GRDname, 'lat_rho');      % RHO-latitude
S.lon_u       = nc_read(GRDname, 'lon_u');        % U-longitude
S.lat_u       = nc_read(GRDname, 'lat_u');        % U-latitude
S.lon_v       = nc_read(GRDname, 'lon_v');        % V-longitude
S.lat_v       = nc_read(GRDname, 'lat_v');        % V-latitude
S.mask_rho    = nc_read(GRDname, 'mask_rho');     % RHO-mask
S.mask_u      = nc_read(GRDname, 'mask_u');       % U-mask
S.mask_v      = nc_read(GRDname, 'mask_v');       % V-mask
S.angle       = nc_read(GRDname, 'angle');        % curvilinear angle
%%  Set vertical grid variables.
kgrid=0;                                          % RHO-points
[S.s_rho, S.Cs_r]=stretching(S.Vstretching, ...
                             S.theta_s, S.theta_b, S.hc, S.N,         ...
                             kgrid, report);
kgrid=1;                                          % W-points			 
[S.s_w,   S.Cs_w]=stretching(S.Vstretching, ...
                             S.theta_s, S.theta_b, S.hc, S.N,         ...
                             kgrid, report);
%%
%--------------------------------------------------------------------------
%  Interpolate initial conditions from Mercator data to application grid.
%--------------------------------------------------------------------------
disp(' ')
disp('Interpolating from Mercator to ROMS grid ...');
disp(' ')
%%  Uncompress input Mercator files.
% s=unix(['gunzip ',fileT]);
% s=unix(['gunzip ',fileU]);
% s=unix(['gunzip ',fileV]);
%%  Read Mercator data has a time coordinate counter (seconds) that
% starts on 11-Oct-2006.
% time   = nc_read(fileT(1:lenT),'time_counter');
% mydate = datestr(datenum('11-Oct-2006')+time/86400-0.5,0);
% disp([ '    Processing: ', mydate]);
% disp(' ')
%%  Get Mercator grid.
Tlon    =nc_read(fileT,'lon');
Tlat    =nc_read(fileT,'lat');
Tdepth  =nc_read(fileT,'depth');
%%  In the western Pacific, the level 50 (z=5727.9 m) of the Mercator data
%  is all zeros. Our grid needs depth of Zr=-5920 m.  Therefore, the depths
%  are modified in level 49 (z=5274.7 m) to bound the vertical interpolation.
% Tdepth(49)=6000; Udepth(49)=6000; Vdepth(49)=6000;
% Tdepth(50)=6100; Udepth(50)=6100; Vdepth(50)=6100;
%  Read in initial conditions fields.
Temp=nc_read(fileT,'t_mn');
Salt=nc_read(fileS,'s_mn');
Uvel=zeros(size(Temp));
Vvel=zeros(size(Temp));
Zeta=zeros(size(Temp));
Zeta=Zeta(:,:,1);
% Uvel=nc_read(fileU(1:lenU),'vozocrtx');
% Vvel=nc_read(fileV(1:lenV),'vomecrty');
%%  Determine Mercator Land/Sea mask.  Since Mercator is a Z-level model, the mask is 3D.
Rmask3d=ones(size(Temp));
ind=find(Temp >= 40.0);
Rmask3d(ind)=0;
clear ind
%%  Compress input Mercator files.
% s=unix(['gzip ',fileT(1:lenT)]);
% s=unix(['gzip ',fileU(1:lenU)]);
% s=unix(['gzip ',fileV(1:lenV)]);
%%  Set initial conditions time (seconds). 
%  The time coordinate for this
%  ROMS application is "seconds since 2007-01-01 00:00:00". The 0.5
%  coefficient here is to account Mecator daily average.
%  MyTime=time/86400-(datenum('1-Jan-2007')-datenum('11-Oct-2006'))-0.5;
%  S.ocean_time = MyTime*86400;
   S.ocean_time = 86400;               % set to Jan 1, because of forcing
%%  Interpolate free-surface initial conditions.
zeta= mercator2roms('zeta',S,Zeta,Tlon,Tlat,Rmask3d(:,:,1));
%  Compute ROMS model depths.  Ignore free-sruface contribution
%  so interpolation is bounded below mean sea level.
igrid=1;
[S.z_r]=set_depth(S.Vtransform, S.Vstretching,                        ...
                  S.theta_s, S.theta_b, S.hc, S.N,                    ...
		  igrid, S.h, zeta);     
igrid=3;
[S.z_u]=set_depth(S.Vtransform, S.Vstretching,                        ...
                  S.theta_s, S.theta_b, S.hc, S.N,                    ...
		  igrid, S.h,zeta);
igrid=4;
[S.z_v]=set_depth(S.Vtransform, S.Vstretching,                        ...
                  S.theta_s, S.theta_b, S.hc, S.N,                    ...
		  igrid, S.h,zeta);
%  Compute ROMS vertical level thicknesses (m).
N=S.N;
igrid=5;
[S.z_w]=set_depth(S.Vtransform, S.Vstretching,                        ...
                  S.theta_s, S.theta_b, S.hc, S.N,                    ...
		  igrid, S.h, zeta);

S.Hz=S.z_w(:,:,2:N+1)-S.z_w(:,:,1:N);
	      
%  Interpolate temperature and salinity.
temp = mercator2roms('temp',S,Temp,Tlon,Tlat,Rmask3d,Tdepth);
salt = mercator2roms('salt',S,Salt,Tlon,Tlat,Rmask3d,Tdepth);
Urho = zeros(size(temp));
Vrho = zeros(size(temp));
% Urho=mercator2roms('u'   ,S,Uvel,Tlon,Tlat,Rmask3d,Udepth);
% Vrho=mercator2roms('v'   ,S,Vvel,Tlon,Tlat,Rmask3d,Vdepth);
% Urho=mercator2roms('u'   ,S,Uvel,Tlon,Tlat,Rmask3d,Udepth);
% Vrho=mercator2roms('v'   ,S,Vvel,Tlon,Tlat,Rmask3d,Vdepth);

%  Process velocity: rotate and/or average to staggered C-grid locations.
[u,v]=roms_vectors(Urho,Vrho,S.angle,S.mask_u,S.mask_v);

%  Compute barotropic velocities by vertically integrating (u,v).
[ubar,vbar]=uv_barotropic(u,v,S.Hz);

%--------------------------------------------------------------------------
%  Create initial condition Netcdf file.
%--------------------------------------------------------------------------
if (CREATE)
  [~]=c_initial(S);

%%  Set attributes for "ocean_time".
%   avalue='seconds since 2007-01-01 00:00:00';
%   [status]=nc_attadd(INIname,'units',avalue,'ocean_time');
%   
%   avalue='gregorian';
%   [status]=nc_attadd(INIname,'calendar',avalue,'ocean_time');

%%  Set global attributes.
%   avalue='Philippine Archipelago Straits, ~5.5 km resolution, Grid b';
%   [status]=nc_attadd(INIname,'title',avalue);
% 
%   avalue='Mercator system PSY3V2 daily average, 0.25 degree resolution';
%   [status]=nc_attadd(INIname,'data_source',avalue);
% 
%   [status]=nc_attadd(INIname,'grd_file',GRDname);
end

%--------------------------------------------------------------------------
%  Write out initial conditions.
%--------------------------------------------------------------------------

if (CREATE)
  disp(' ')
  disp([ 'Writing initial conditions ...']);
  disp(' ')

  [~]=nc_write(INIname, 'spherical',   S.spherical);
  [~]=nc_write(INIname, 'Vtransform',  S.Vtransform);
  [~]=nc_write(INIname, 'Vstretching', S.Vstretching);
  [~]=nc_write(INIname, 'theta_s',     S.theta_s);
  [~]=nc_write(INIname, 'theta_b',     S.theta_b);
  [~]=nc_write(INIname, 'Tcline',      S.Tcline);
  [~]=nc_write(INIname, 'hc',          S.hc);
  [~]=nc_write(INIname, 's_rho',       S.s_rho);
  [~]=nc_write(INIname, 's_w',         S.s_w);
  [~]=nc_write(INIname, 'Cs_r',        S.Cs_r);
  [~]=nc_write(INIname, 'Cs_w',        S.Cs_w);

  [~]=nc_write(INIname, 'h',           S.h);
  [~]=nc_write(INIname, 'lon_rho',     S.lon_rho);
  [~]=nc_write(INIname, 'lat_rho',     S.lat_rho);
  [~]=nc_write(INIname, 'lon_u',       S.lon_u);
  [~]=nc_write(INIname, 'lat_u',       S.lat_u);
  [~]=nc_write(INIname, 'lon_v',       S.lon_v);
  [~]=nc_write(INIname, 'lat_v',       S.lat_v);
  
  IniRec = 1;

  [~]=nc_write(INIname, 'ocean_time', S.ocean_time, IniRec);

  [~]=nc_write(INIname, 'zeta', Zeta, IniRec);
  [~]=nc_write(INIname, 'ubar', ubar, IniRec);
  [~]=nc_write(INIname, 'vbar', vbar, IniRec);
  [~]=nc_write(INIname, 'u',    u,    IniRec);
  [~]=nc_write(INIname, 'v',    v,    IniRec);
  [~]=nc_write(INIname, 'temp', temp, IniRec);
  [~]=nc_write(INIname, 'salt', salt, IniRec);
end

%--------------------------------------------------------------------------
%  Set masking indices to facilitate plotting.  They can be used to
%  replace ROMS Land/Sea mask with NaNs.
%--------------------------------------------------------------------------

inr2d=find(S.mask_rho == 0);
inr3d=find(repmat(S.mask_rho,[1,1,N]) == 0);

Results:

Interpolated temp number of NaNs replaced with nearest value = 37860258
Min= NaN Max= NaN
Interpolated salt number of NaNs replaced with nearest value = 37860258
Min= NaN Max= NaN

Writing initial conditions ...

Wrote h Min= 6.00000e+01 Max= 5.00000e+03
Wrote lon_rho Min= 1.07531e+02 Max= 1.39895e+02
Wrote lat_rho Min=-1.50192e+01 Max= 1.34821e+01
Wrote lon_u Min= 1.07553e+02 Max= 1.39873e+02
Wrote lat_u Min=-1.50192e+01 Max= 1.34821e+01
Wrote lon_v Min= 1.07531e+02 Max= 1.39895e+02
Wrote lat_v Min=-1.49971e+01 Max= 1.34597e+01
Wrote zeta into record: 0001, Min= 0.00000e+00 Max= 0.00000e+00
Wrote ubar into record: 0001, Min= NaN Max= NaN
Wrote vbar into record: 0001, Min= NaN Max= NaN
Wrote u into record: 0001, Min= 0.00000e+00 Max= 0.00000e+00
Wrote v into record: 0001, Min= 0.00000e+00 Max= 0.00000e+00
Wrote temp into record: 0001, Min= 0.00000e+00 Max= 0.00000e+00
Wrote salt into record: 0001, Min= 0.00000e+00 Max= 0.00000e+00

I have some problems with these codes, the results of interpolation are 0.0 for temperature and salinity. anyone can help me to solve this problem?

wilkin · #2 Unread post by **wilkin** » Tue Apr 23, 2019 10:54 am

Looks to me like you're using a World Ocean Atlas climatology as input.

Code: Select all

Tfile   = 'woa18_ann_temp.nc';
Sfile   = 'woa18_ann_salt.nc';

This script is designed to process Mercator-Océan data from Copernicus.eu.

corvianawatie · Wed Apr 24, 2019 12:28 am

wilkin wrote:Looks to me like you're using a World Ocean Atlas climatology as input.
Code: Select all
Tfile   = 'woa18_ann_temp.nc';
Sfile   = 'woa18_ann_salt.nc';
This script is designed to process Mercator-Océan data from Copernicus.eu.

Thank you for your replies. Isn't it has similar data format?

By the way, I interpolated the WOA data (1 x 1 degree) to my grid (5 x 5 km), which might be the problem. It couldn't interpolate to a very fine resolution. Am I right?

wilkin · #4 Unread post by **wilkin** » Wed Apr 24, 2019 12:56 am

Thank you for your replies. Isn't it has similar data format?

Similar, yes, but the script is expecting certain names and conventions for the coordinates and dimensions. Code doesn't do similar.

Code: Select all

By the way, I interpolated the WOA data (1 x 1 degree) to my grid (5 x 5 km), which might be the problem. It couldn't interpolate to a very fine resolution. Am I right?

1 degree to 5 km isn't necessarily a problem for initialization. ROMS will take the smooth initial condition and evolve the higher resolution dynamics. But pay attention to the land/sea mask near the coast because you likely have points where WOA must be extrapolated into coastal bays. You can't have NaNs or missing values in ROMS initial conditions. If you interpolate over land you might get unexpected results.

Also, if you have ROMS points that are deeper that WOA valid data you might need a thoughtful extrapolation to depth.

corvianawatie · Fri Apr 26, 2019 12:53 am

wilkin wrote:
Thank you for your replies. Isn't it has similar data format?
Similar, yes, but the script is expecting certain names and conventions for the coordinates and dimensions. Code doesn't do similar.
Code: Select all
By the way, I interpolated the WOA data (1 x 1 degree) to my grid (5 x 5 km), which might be the problem. It couldn't interpolate to a very fine resolution. Am I right?
1 degree to 5 km isn't necessarily a problem for initialization. ROMS will take the smooth initial condition and evolve the higher resolution dynamics. But pay attention to the land/sea mask near the coast because you likely have points where WOA must be extrapolated into coastal bays. You can't have NaNs or missing values in ROMS initial conditions. If you interpolate over land you might get unexpected results.

Also, if you have ROMS points that are deeper that WOA valid data you might need a thoughtful extrapolation to depth.

Thank you for your suggestion. I am still trying to modify the code, to make it suitable for my datasets. And I will take more attention to the land/sea masking near the coast and bottom water..

Ocean Modeling Discussion

roms preprocessing tools (d_mercator2roms.m)

roms preprocessing tools (d_mercator2roms.m)

Re: roms preprocessing tools (d_mercator2roms.m)

Re: roms preprocessing tools (d_mercator2roms.m)

Re: roms preprocessing tools (d_mercator2roms.m)

Re: roms preprocessing tools (d_mercator2roms.m)