Ocean Modeling Discussion

ROMS/TOMS

Search for:
It is currently Mon Jul 22, 2019 7:07 pm




Post new topic Reply to topic  [ 5 posts ] 

All times are UTC

Author Message
PostPosted: Tue Apr 23, 2019 3:23 am 
Offline

Joined: Thu May 14, 2015 4:50 pm
Posts: 11
Location: Indonesian Institute of Science (LIPI)
Code:
%
%  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?


Top
 Profile  
Reply with quote  
PostPosted: Tue Apr 23, 2019 10:54 am 
Offline
User avatar

Joined: Mon Apr 28, 2003 5:44 pm
Posts: 491
Location: Rutgers University
Looks to me like you're using a World Ocean Atlas climatology as input.

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


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

_________________
John Wilkin: DMCS Rutgers University
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu


Top
 Profile  
Reply with quote  
PostPosted: Wed Apr 24, 2019 12:28 am 
Offline

Joined: Thu May 14, 2015 4:50 pm
Posts: 11
Location: Indonesian Institute of Science (LIPI)
wilkin wrote:
Looks to me like you're using a World Ocean Atlas climatology as input.

Code:
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?


Top
 Profile  
Reply with quote  
PostPosted: Wed Apr 24, 2019 12:56 am 
Offline
User avatar

Joined: Mon Apr 28, 2003 5:44 pm
Posts: 491
Location: Rutgers University
Quote:
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:
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.

_________________
John Wilkin: DMCS Rutgers University
71 Dudley Rd, New Brunswick, NJ 08901-8521, USA. ph: 609-630-0559 jwilkin@rutgers.edu


Top
 Profile  
Reply with quote  
PostPosted: Fri Apr 26, 2019 12:53 am 
Offline

Joined: Thu May 14, 2015 4:50 pm
Posts: 11
Location: Indonesian Institute of Science (LIPI)
wilkin wrote:
Quote:
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:
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..


Top
 Profile  
Reply with quote  
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 5 posts ] 

All times are UTC


Who is online

Users browsing this forum: No registered users and 4 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
cron
Powered by phpBB® Forum Software © phpBB Group