| 1 |
gforget |
1.1 |
|
| 2 |
|
|
%%0) setup: |
| 3 |
|
|
|
| 4 |
|
|
%add matlab path and load ECCO v4 grid: |
| 5 |
|
|
p = genpath([pwd filesep 'gcmfaces/']); addpath(p); |
| 6 |
|
|
p = genpath([pwd filesep 'MITprof/']); addpath(p); |
| 7 |
|
|
p = genpath([pwd filesep 'm_map/']); addpath(p); |
| 8 |
|
|
grid_load; gcmfaces_global; |
| 9 |
|
|
m=mygrid.mskC(:,:,1);%standard land mask |
| 10 |
|
|
|
| 11 |
|
|
%input directories: |
| 12 |
|
|
dirNctiles='release1/nctiles_climatology/'; |
| 13 |
|
|
dirBinary='release1/binary_inputs/'; |
| 14 |
|
|
|
| 15 |
|
|
%constant parameters: |
| 16 |
|
|
parms.rhoconst =1029; %sea water density [in kg/m^3] |
| 17 |
|
|
parms.rhoi = 910; %sea ice density [in kg/m^3] |
| 18 |
|
|
parms.rhosn = 330; %snow density [in kg/m^3] |
| 19 |
|
|
|
| 20 |
|
|
return; |
| 21 |
|
|
|
| 22 |
|
|
|
| 23 |
|
|
%%1) sea surface temperature and related variables: |
| 24 |
|
|
|
| 25 |
|
|
%THETA: potential temperature [in degree C] |
| 26 |
|
|
%SALT: salinity [in psu] |
| 27 |
|
|
%SIarea: fractional ice-covered area [0 to 1] |
| 28 |
|
|
%MXLDEPTH: Mixed-Layer Depth [in m] |
| 29 |
|
|
% |
| 30 |
|
|
%note: MXLDEPTH is based on the Kara, Rochford, and Hurlburt 2000 formula |
| 31 |
|
|
|
| 32 |
|
|
if isdir(dirNctiles); |
| 33 |
|
|
sst_ecco_clim=read_nctiles([dirNctiles 'THETA/THETA'],'THETA',[],1); |
| 34 |
|
|
sss_ecco_clim=read_nctiles([dirNctiles 'SALT/SALT'],'SALT',[],1); |
| 35 |
|
|
mld_ecco_clim=read_nctiles([dirNctiles 'MXLDEPTH/MXLDEPTH'],'MXLDEPTH'); |
| 36 |
|
|
ice_ecco_clim=read_nctiles([dirNctiles 'SIarea/SIarea'],'SIarea'); |
| 37 |
|
|
end; |
| 38 |
|
|
|
| 39 |
|
|
%sst_reynolds: Reynolds monthly mean sea surface temperature maps |
| 40 |
|
|
%iicecover_nsidc: NSIDC monthly mean fractional ice-covered area [0 to 1] |
| 41 |
|
|
% |
| 42 |
|
|
%note: references are provided in Forget et al 2015 (GMD) |
| 43 |
|
|
|
| 44 |
|
|
if isdir(dirBinary); |
| 45 |
|
|
sst_reynolds_1992=read_bin([dirBinary 'sst_reynolds/reynolds_oiv2_r1_1992']); |
| 46 |
|
|
sst_reynolds_1992(sst_reynolds_1992==0)=NaN;%0's denote missing values in these files |
| 47 |
|
|
ice_nsidc_1992=read_bin([dirBinary 'icecover_nsidc/nsidc79_monthly_1992']); |
| 48 |
|
|
ice_nsidc_1992(ice_nsidc_1992==-9999)=NaN;%-9999's denote missing values in these files |
| 49 |
|
|
end; |
| 50 |
|
|
|
| 51 |
|
|
|
| 52 |
|
|
%%2) heat and freshwater fluxes into the ocean (and ice pack) |
| 53 |
|
|
|
| 54 |
|
|
%oceFWflx: net surface Fresh-Water flux into the ocean [in kg/m^2/s] |
| 55 |
|
|
%oceQnet: net surface heat flux into the ocean [in W/m^2] |
| 56 |
|
|
%SIatmFW: net freshwater flux from atmosphere & land [in kg/m^2/s] |
| 57 |
|
|
%SIatmQnt: net heat flux to atmosphere & land [in W/m^2] |
| 58 |
|
|
% |
| 59 |
|
|
%note: the key difference between oceQnet and SIatmQnt (or oceFWflx and SIatmFW) is |
| 60 |
|
|
% that the former is computed at the bottom of the ice pack when there is sea-ice |
| 61 |
|
|
%note: the MITgcm sign conventions are opposite for SIatmQnt and oceQnet so below |
| 62 |
|
|
% I flip the sign of SIatmQnt upon reading it from file. |
| 63 |
|
|
|
| 64 |
|
|
if isdir(dirNctiles); |
| 65 |
|
|
qnet_to_oce=read_nctiles([dirNctiles 'oceQnet/oceQnet'],'oceQnet'); |
| 66 |
|
|
qnet_to_surf=-read_nctiles([dirNctiles 'SIatmQnt/SIatmQnt'],'SIatmQnt'); |
| 67 |
|
|
qnet_to_ice=qnet_to_surf-qnet_to_oce; |
| 68 |
|
|
|
| 69 |
|
|
fwf_to_oce=read_nctiles([dirNctiles 'oceFWflx/oceFWflx'],'oceFWflx'); |
| 70 |
|
|
fwf_to_surf=read_nctiles([dirNctiles 'SIatmFW/SIatmFW'],'SIatmFW'); |
| 71 |
|
|
fwf_to_ice=fwf_to_surf-fwf_to_oce; |
| 72 |
|
|
end; |
| 73 |
|
|
|
| 74 |
|
|
|
| 75 |
|
|
%%3) sea surface height and bottom pressure: |
| 76 |
|
|
|
| 77 |
|
|
%ETAN: Surface Height Anomaly [in m] |
| 78 |
|
|
%sIceLoad: sea-ice loading [in km/m^2] |
| 79 |
|
|
%PHIBOT: Bottom Pressure Anomaly [in m^2/s^2] |
| 80 |
|
|
% |
| 81 |
|
|
%note: ETAN corresponds to the ocean-ice interface when there is seaice. The |
| 82 |
|
|
% height contribution of seaice+snow (sIceLoad/parms.rhoconst) therefore |
| 83 |
|
|
% needs to be accounted for to allow comparison with altimetry. |
| 84 |
|
|
%note: bp is converted to equivalent sea level meters by dividing by g. |
| 85 |
|
|
|
| 86 |
|
|
if isdir(dirNctiles); |
| 87 |
|
|
ETAN=read_nctiles([dirNctiles 'ETAN/ETAN'],'ETAN'); |
| 88 |
|
|
sIceLoad=read_nctiles([dirNctiles 'sIceLoad/sIceLoad'],'sIceLoad'); |
| 89 |
|
|
PHIBOT=read_nctiles([dirNctiles 'PHIBOT/PHIBOT'],'PHIBOT'); |
| 90 |
|
|
|
| 91 |
|
|
ssh=ETAN+sIceLoad/parms.rhoconst; |
| 92 |
|
|
bp=PHIBOT/9.81; |
| 93 |
|
|
end; |
| 94 |
|
|
|
| 95 |
|
|
%%4) near surface velocity and related vector variables: |
| 96 |
|
|
|
| 97 |
|
|
%UVELMASS: first component of horizontal velocity [in m/s] |
| 98 |
|
|
%VVELMASS: second component of horizontal velocity [in m/s] |
| 99 |
|
|
% |
| 100 |
|
|
%note: the loaded vector components are along grid line directions. To |
| 101 |
|
|
% allow for comparison with e.g. current meters they need to be |
| 102 |
|
|
% converted to zonal and meridional components using calc_UEVNfromUXVY. |
| 103 |
|
|
|
| 104 |
|
|
UVELMASS=read_nctiles([dirNctiles 'UVELMASS/UVELMASS'],'UVELMASS',[],1); |
| 105 |
|
|
VVELMASS=read_nctiles([dirNctiles 'VVELMASS/VVELMASS'],'VVELMASS',[],1); |
| 106 |
|
|
oceTAUX=read_nctiles([dirNctiles 'oceTAUX/oceTAUX'],'oceTAUX'); |
| 107 |
|
|
oceTAUY=read_nctiles([dirNctiles 'oceTAUY/oceTAUY'],'oceTAUY'); |
| 108 |
|
|
|
| 109 |
|
|
[Ue,Vn]=calc_UEVNfromUXVY(UVELMASS,VVELMASS); |
| 110 |
|
|
[Te,Tn]=calc_UEVNfromUXVY(oceTAUX,oceTAUY); |
| 111 |
|
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph