To edit HTML, one can use the ASHE by
typing
ashe
from any machine. You can also get the documentation.
Also, see the beginner's
guide to HTML.
Hydrosys/Neutral contains the neutral surfaces stuff (see this section).
Up: General Menu
The structure of the data is relatively staightforward and described in the routines getsdat.m (for the station data) and getpdat.m (for the pair data).
Example: loads the data for section A36N, in the Atlantic.
path(path,'/data4/ganacho/web/Hydrosys')
datadir='/data4/ganacho/web/Hydrosys/Hdata/';
secid='a36n';
getsdat
g_pltopo(Nstat, Slon,Pres, Maxd(:,Itemp), Botd)
hold on;title(['Potential Temperature, Station data, ' secid ])
path(path,'/data4/ganacho/web/Matlab')
cint=[0 1 2 3 4 5 10 15 20];
[c,h]=contour(Slon,-Presctd,temp,cint)'
clabel(c,h)
xlabel('longitude');set(gca,'ylim',[-6000,0]);
Example: to plot the potential temperature, section A36N
path(path,'/data4/ganacho/web/Hydrosys')
datadir='/data4/ganacho/web/Hydrosys/Hdata/';
secid='a36n';
getpdat
g_pltopo(Nstat, Slon,Pres, Maxd(:,Itemp), Botd)
hold on;title(['Potential Temperature, Pair data, ' secid ])
path(path,'/data4/ganacho/web/Matlab')
cint=[0 1 2 3 4 5 10 15 20];
[c,h]=contour(Plon,-Presctd,temp,cint);
clabel(c,h) xlabel('longitude');set(gca,'ylim',[-6000,0]);
To read data, one wants to read the header file first: datadir='/data4/ganacho/web/Hydrosys/Hdata/';
secid='a36n';
hdrname = [secid '_pair.hdr.mat']; eval(['load ' datadir hdrname
])
For example, to read the temperature only, type : temp=rhydro([datadir
Pairfiles(Itemp,:)], Precision(Itemp,:),...
MPres(Itemp), Npair, Maxdp(:,Itemp));
To read all the properties (indice iprop), one by one: for iprop = 1:Nprop pname = Propnm(iprop,:); %property name punit = Propunits(iprop,:); %property unit pfile=[datadir Pairfiles(iprop,:)];%property file name prop=rhydro(pfile,Precision(iprop,:),MPres(iprop),Npair,Maxdp(:,iprop)); eval([Propnm(iprop,:) '=prop;']); end So that each property is loaded in each variable (temp, sali, ...).
To read the velocity: gvel = rhydro([datadir Velfile], Velprec, ... MPres(iprop), Npair, Maxdp(:,iprop));
To save the properties: pprop=temp; %save the temperature ovw=0; %do not overwrite existing file whydro(pprop,[OPdir Pairfiles(Itemp,:)],Precision(Itemp,:),Maxdp,ovw)
The examples are taken from getpdat.m and g_pairset.m
Up: General Menu
Mergech output files |
||||||||
|---|---|---|---|---|---|---|---|---|
| Section | root name | nstat | first | last | nvar | ndep | CTD | |
| A36NA24Nflst | at109 | 213 | 1 102192 | 101191202 | 5 | 37 | yes | |
| A48N | hudson | 78 | * | * | 5 | 37 | yes | |
| I32S | i5 | 106 | * | * | 5 | 37 | yes | |
| I18S | at93 | 68 | * | * | 5 | 37 | no | |
| MZ_NMZ_S | mozam | 348 | * | * | 5 | 37 | no | |
The data from a particular cruise or set of cruises are in the "root
name" file, in the directory
/data39/alison/phd/data/<rootname>.<ext>,
where ext='.hdr' for the header, and '.std' for the data. The data file
size must be nstat*nvar(+1)*ndep*4. In the case of ctd data, then it
includes
the dynamic height, so there is one more variable (+1). The record
length
is nvar(+1)*ndep*4, as there is one record for each station. To convert
a file to the new format, copy the example file conv_at109.m,
that does the conversion for the North Atlantic, and modify it to the
new
set of parameters. For converting the station data only, set the
parameter
p_conv_stat
to 1, and p_conv_pair to zero. Verify afterwards how the
fields
look like. The latter routine calls automatically conv_sdata.m,
that actually does the conversion.
Most of the geovel output files sit in the /data39/alison/phd/data/
directory. The names of the pair data files are organized as follows.
The root name takes the extension .hdr, .dat for station data and _pair.hdr, _pair.dat for pair data, .hdr standing for "header information" and .dat for the data fields. Information about the missing parameters can be found in the mkequats parameter file (see hydroinv ) as the mkequats program reads the pair data files. nstat stands for the total number of station in the station header file; npair stands for the total number of pairs in header and data files; first and last pairs are given for each section, when one file contains several sections. nvar is the number of properties in each file (temperature, ...). The velocity is always read after the nvar properties.
ndep contains the total number of standard depths plus the depth at
the bottom. The "37" standard depths are in std.37
the "38" standard depths are in std.38
In the routine read_pairdata, the values at the bottom are removed
from the data because they contain no actual iinformation and they
intriduce
complications in the plotting routines.
The record lenght depends on the cruise. It is usually
(nvar+1)*ndep*4.
See the comments on the first lines of geovel.f
for a complete description, although there are unclear things in the
comments.
To check, verify the size of the pair data file: it should be
npairtotal*(nvar+1)*ndep*4*npair.
Note: it is also possible to access directly the data from the old system (without translation into the hydrosys format) with the routines read_pairdata.m and read_stathdr.m. They read all the information in the pair header and data files, and in the station header file.
Up: General Menu
....
i.e., when I spot an outlier, (in the last step of the program) I zoom
to get the station number, and then I write down the depth range.
After w2o_spotoutliers is run, I set gis=[
sili( idepth_outlier , istat_outlier ) = NaN
for each outlier. These lines are copied in the woce2obs.log filein
the run directory for documentation.
Finally, the data are saved: w2o_saveobsfile (and plotted)
Up: General Menu
Parameter file example: ctd2std_input_p21w.m Copy and adapt this file.
Parameters:
To run ctd2std, first type the name of the adapted input file at
the matlab prompt; for example: ctd2std_input_p21w. The
program
will then load the header file of the specified section, in this case
p21w,
load the standard depths, and read the ctd data from that section. Then
ctd_get.m averages
the
data around each standard depth point, interp/extrapolating in the
intervals
if necessary. Next, ctd_step3.m
fills in missing data and stations, first by automatically copying near
the surface and near deep bottoms if necessary. Any remaining missing
data
or stations are plotted for the user, who can then choose from a number
of different filling methods including copy top/bottom and vertical and
horizontal interp/extrapolation for each continuous hole. The user can
also remove data and refill it. If p_ooxyg is activated, the
oxygen
data for each station will be plotted along with that station's oxygen
bottle data. These stations can either be accepted as they are or
altered
graphically with a mouse. Any missing or altered data is recorded in
the
cell variable ctd2std_rec, which is then saved in the header
file.
This variable contains the indices of the previously missing or altered
data for each station and property (for example, ctd2std_rec{57,3}
would contain indices of filled or altered oxygen data for station 57).
These can then be used as input for cpres to find the
corresponding
pressure ranges or prop(station number,property index) to
find
the filled/altered data values. After ctd_step3, the program saves the
treated data, plots contours of the section for each property, and
prints
them if desired.
Up: General Menu
Parameter file example: obs2std_input_p21w.m Copy and adapt this file.
Parameters:
To run obs2std, first type the name of the adapted input file at
the matlab prompt; for example: obs2std_input_p21w. The
program
will then load the header file of the specified section, in this case
p21w,
load the standard depths, and read the observed data from that section.
obs2_step1.m
vertically
interpolates automatically within certain adjustable limits and removes
closely spaced data. obs2_step2.m
vertically extrapolates automatically at deep enough bottoms and at the
surface. obs2_step3.m
fills holes and missing stations interactively, using user input
choices
such as vertical inter/extapolation, horizontal inter/extrapolation
along
isopycnals, and copy to the top/bottom. Any user-altered data are
stored
by their corresponding indices in stdd (the standard pressure vector)
in
the cell variable obs2std_rec, which is later saved in the
header
file. This variable contains the indices of the previously missing or
altered
data for each station and property (for example, obs2std_rec{57,1}
would contain indices of filled or altered oxygen data (or whichever
was
the first observation property treated for that section) for station
57).
These can then be used as input for stdd to find the
corresponding
pressure ranges or bprop(station number,property index) to
find
the filled/altered data values. obs2_step4.m
looks for data anomalies and plots graphs with red diamonds around
suspicious
points. The user can either OK these plots by double clicking or remove
data around these areas, which will be filled later on with obs2_step3.
Finally, the standardized data is saved in a memory-efficient format
and
contours of the section for each property can be plotted and printed if
desired.
Up: General Menu
Up: General Menu
Parameter file example: gv_input_a36n.m.
Copy and adapt this file. The parameter secid must correspond
to the file names.
Parameters:
Up: General Menu
Up: General Menu
Up: General Menu
The mkequats.m routine first call the subroutine mk_set_layprops.m that computes layer interface position, integrated properties/transports through each section and layer, and averages over the whole box.
To compute horizontal areas of layer interfaces:
Section 1
*<------- length 1 --------->
*oooooooooooooooooooooooooooo*
*
*
******ooooooooooooooooooooooooo*
<-zc-><-------length 2-------->
Section 2
ydist=distance between the two sections
zc=zonal coast
The area is computed as in a trapeze:
A=ydist*(sum(lay. interface. lenght)+zc)/2
So that the decrease of the area with depth (because of topography)
is taken into account by lay. interface. lenght. Of course we do not
have
a trapeze everywhere. This is why the parameter ydist is not provided
but
computed from the harea parameter and ydist*(sum(lay. interface.
lenght)+zc)/2
so that it gives the right area at the surface. Zonal coast is used to
increase the area from that of a trapezoidal calculation if needed. The
area does not need to be very accurate in the inverse calculation.
To compute a horizontal area enclosed by a section/continent
polygon:
areas.f, binary SUN areas
(from /data39/alison/phd/progs/bin/), or the Matlab equivalent areas.m
which, given a set of (latitude, longitude) locations returns the area
of the polygon they define (as long as all interior angles measure less
than 180 degrees). Input file example: natl.inp
An alternative to areas.f and its equivalents is neutarea.m, which finds the area of various neutral (gamma) surfaces within a longitude/latitude box. The program reads gridded Levitus data with a program called readlev.m. The inputs are the northern and southern Latitude boundaries, the eastern and western Longitude boundaries, and Input directory for the levitus data file, the Output directory for the ascii file containing the computed areas, and the values of gamma at which to find the areas (an input gamma of 0 will produce the surface area of all ocean within the box that is both recognized by the Levitus data set and defined in the program gamman.m, which is also called by neutarea.m).
To make inversions involving more than one box, these matrix equations will have to be combined to form a global matrix (script weldquats.m).
The model will allow an adjustment of a priori size dEkstd
to the Ekman transport. The adjustment is solved for throught the
inversion.
The solution vector contains the Ekman adjustment, 1 unknown for each
section,
just after the ref. velocities of each section. isb=section
number.
Allows to select only some pairs to enter the conservation equation.
Convenient if the section goes through and island for example, or to
define
sub-boxes. The output matrix has the same structure, but pairs that are
not selected have zero columns. All other parameters (gsecs.npair,
gsecs.binit, gsecs.bstd, flxeq.fpair, flxeq.lpair, etc.) must be
precised
with the original pair indices, not with indices within the
subselection. isb=section
number.
Only tracers are mixed, not mass (e.g. Pedlosky, page 6-7 ?). One
a priori value of Kz is given for each interface.
ieq is the equation indice ex: boxi.conseq.anom=[0,1,1,1,1,0,1,1,1];.
Mass is nerver an anomaly. Anomaly equations are of the form
div(vC)-avgC
* div(v), that is the mass conservation equation times the average
tracer
property is substracted from the tracer conservation equation. To
weight
those equation, one can use mass weight/std. deviation of the property
within the layer. In this case, use the Norm, e.g., boxi.conseq.norm={'one','std','std','std','std','one',...
'std','std','std'};
Norm std is used only for conservation equations. If std is asked,
norm 1 (no norm) is used for all flux equations. If rms is asked, rms
is
used for conservation and flux equations.
Correction to initial freshwater flux (0 by default). If there is
outcropping, the average of the percentage of total Ekman transport in
each layer (paramameter gsecs.perEk(ilayer,isection) ) is
taken
to determine the percentage of freshwater that goes into each layer.
This
may be not the most accurate representation... The freshwater
correction
(as for initial freshwater, positive is precipitation) comes right
after
dianeutral transfers in the solution vector.
Example for Gauss-Markov:
caOptions for SVD:
load mendocino_1_equn %load A matrix and others
p_scale_row_norm=0; %no row normalization
p_scale_col_norm=0; %no column normalization
p_plots=1; %plot the results
p_resobs=1; %plot the resolution and observation matrices
invmethod='Gauss'; %Gauss-Markov estimator
boxinvert
p_scale_row_norm=1; %row normalizationUp: General Menu
p_scale_col_norm=1; %column normalization
invmethod='svd';
gKs=[10:20]; %cut-off ranks for which to plot the solution
hydrotrans.m computes transports and uncertainties. Plot them (by calling hydrotrans_plot.m). Procedure: set parameters/type hydrotrans.Main parameters:
p_plt_cumtrans=0; %computes cumulative transports and uncertainties from beginning of the sectionUp: General Menu
p_calculate_Ekman=0; %(if adjustement to Ekman transport was not used)
%re-calculate the net flux without the initial Ekman transport
%it can be interpreted as a new Ekman flux if net mass flux is 0 or Bering
gprop2plot=[1,2,6]; %results to plot (MASS,HEAT,SILI) (default=all)
p_plots=1; %do the plots (default=1)
p_plotstream=1; %plot streamfunctions
p_bw=1; %plot black and white (default=0)
p_posterplot=0; % 0=plot relative and absolute transports and residuals
% 2=plot absolute transports only
% 1=plot same as 2 but on a very small page
p_resonly=0; %plot only the residuals
p_res_units=1; %convert nutrient/O2 residuals to mol/m2/yr
To compute global transports and residuals:
1) run_global for each box. In the loop, include save_transec, get_sumresiduals and hydrotrans (p_plt_streamf=0;p_pltstream=0).
2) save global transports and residuals :
eval(['save ' ModelidG '_Transec.mat Transec Amatnets_glb'])
eval(['save ' ModelidG '_lay2sum.mat lay2sum'])
3) make global calculation based on bhat_global, P, Transec, Amatnets_glb and lay2sum. See glbdiagnostic.m
Main parameters:
Glbfigs/glbcumtransmap.m: cumulative transports on map
Up: General Menu