Primary Data Description: Li & von Storch, 2019
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1. Model Data : a simulation with the Max-Planck-Institut Ocean Model (MPIOM) on the 0.1-degree grid TP6ML40
- the code of the model is mpiom-1.6.3 (r4550)
- the run script for spinup is hel17089-LTO_TP6ML40_mpiom-trunk_OMIP.job
- the run script for production run is hel17153-DZA_TP6ML40_mpiom-trunk_NCEP.job
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2. Analyses of the model data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for harmonic analysis (The inputs for Fig.2,3,4,6-12 are amplitudes and phases obtained from
the harmonic analysis using the script given below)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.1: Scatter diagram (Matlab)
clear all;
clc;
rhoref = 1025;
Tm2 = 12.42*3600;
wm2 = 2*pi/Tm2;
omega = 7.29*10^(-5);
tm = 0:1800:12.42*3600;
grv = 9.8;
plon = [190 199];
plat = [22 26];
regstr = 'Hawaii';
% --1 -- %
% -- import data -- %
disp('import uprm data')
dir = '/work/mh0256/m300212/ITgeneration_Oct2018/proc1_BartBarc_UV/data/proc1_BartBarc_AmpPhs/';
uamp = ncread([dir 'Ubarcl_2012Jan_stormtide-ncep_m2_amp.nc'],'m2amp');
uphs = ncread([dir 'Ubarcl_2012Jan_stormtide-ncep_m2_phase.nc'],'m2phase');
ulon = ncread([dir 'Ubarcl_2012Jan_stormtide-ncep_m2_amp.nc'],'lon');
ulat = ncread([dir 'Ubarcl_2012Jan_stormtide-ncep_m2_amp.nc'],'lat');
disp('import vprm data')
vamp = ncread([dir 'Vbarcl_2012Jan_stormtide-ncep_m2_amp.nc'],'m2amp');
vphs = ncread([dir 'Vbarcl_2012Jan_stormtide-ncep_m2_phase.nc'],'m2phase');
vlon = ncread([dir 'Vbarcl_2012Jan_stormtide-ncep_m2_amp.nc'],'lon');
vlat = ncread([dir 'Vbarcl_2012Jan_stormtide-ncep_m2_amp.nc'],'lat');
clear dir;
disp('import grid distance')
dir = '/work/mh0256/m300212/origin_data/tp6ml40_info/';
fname = [dir 'TP6ML40_fx.nc'];
tem_dlxu = ncread(fname,'dlxu');
tem_dlyv = ncread(fname,'dlyv');
dlxu = tem_dlxu(2:end-1,:);
dlyv = tem_dlyv(2:end-1,:);
disp('import wate depth')
tem_ddpo = ncread(fname,'ddpo');
tem_ddpo(:,:,1) = tem_ddpo(:,:,1)+8;
ddpo = tem_ddpo(2:end-1,:,:);
sddpo = nansum(ddpo,3);
sddpo(sddpo==0) = nan;
clear dir fname tem_*;
disp('import density perturb')
dir = '/work/mh0256/m300212/STORMTIDE_run/Tides_NCEP_SAL2/Data_analysis/proc2_harmonics_velocity_density/rho_total/';
fname = [dir 'rho_harmonics_removeHalo.nc'];
ramp = ncread(fname,'m2amp');
rphs = ncread(fname,'m2phase');
rlon = ncread(fname,'lon');
rlat = ncread(fname,'lat');
rdep = ncread(fname,'lev');
clear dir fname;
% -- 2-- %
% -- make calculation -- %
disp('cal of pprm')
nlon = size(rlon,1);
nlat = size(rlon,2);
ndep = length(rdep);
itm = 1;
new_po = zeros(nlon,nlat,ndep).*nan;
for ilon = 1 : nlon
disp(['pprm cal: ilon=' num2str(ilon)])
for ilat = 1 : nlat
if rlon(ilon,ilat)<plon(1) || rlon(ilon,ilat)>plon(2) || rlat(ilon,ilat)<plat(1) || rlat(ilon,ilat)>plat(2)
continue;
end
rho1d = squeeze(ramp(ilon,ilat,:)) .* cos(wm2*tm(itm)-degtorad(squeeze(rphs(ilon,ilat,:))));
dzw = squeeze(ddpo(ilon,ilat,:));
tempos = find(~isnan(rho1d));
po = zeros(size(rdep)) .* nan;
po(1) = grv *rho1d(1)*dzw(1)/2;
if length(tempos)>1
for idep = 2 : tempos(end)
po(idep) = po(idep-1) + grv*rho1d(idep-1)*dzw(idep-1)/2+grv*rho1d(idep)*dzw(idep)/2;
end
end
po_vMean(ilon,ilat) = nansum(po.*squeeze(ddpo(ilon,ilat,:)))/sddpo(ilon,ilat);
new_po(ilon,ilat,:) = po - po_vMean(ilon,ilat);
clear rho1d dzw tempos po po_vMean;
end
end
disp('u: linear theory cal')
ueq_LHS = zeros(nlon,nlat,ndep).*nan;
ueq_RHS = zeros(nlon,nlat,ndep).*nan;
for ilon = 1 : nlon-1
disp(['ueq cal, ilon-' num2str(ilon)])
for ilat = 1 : nlat
if isnan(new_po(ilon,ilat,1))
continue;
end
for idep = 1 : ndep
f = 2*omega*sin(degtorad(vlat(ilon,ilat)));
part_u = -wm2*uamp(ilon,ilat,idep)*sin(wm2*tm(itm)-degtorad(uphs(ilon,ilat,idep)));
part_v = -f * vamp(ilon,ilat,idep)*cos(wm2*tm(itm)-degtorad(vphs(ilon,ilat,idep)));
ueq_LHS(ilon,ilat,idep) = part_u + part_v;
ueq_RHS(ilon,ilat,idep) = -(new_po(ilon+1,ilat,idep)-new_po(ilon,ilat,idep))/(rhoref*dlxu(ilon,ilat));
clear f part_u part_v;
end
end
end
disp('v: linear equaiton')
veq_LHS = zeros(nlon,nlat,ndep).*nan;
veq_RHS = zeros(nlon,nlat,ndep).*nan;
for ilon = 1 : nlon
disp(['veq cal, ilon-' num2str(ilon)])
for ilat = 1 : nlat-1
if isnan(new_po(ilon,ilat,1))
continue;
end
for idep = 1 : ndep
f = 2*omega*sin(degtorad(ulat(ilon,ilat)));
part_u = f*uamp(ilon,ilat,idep)*cos(wm2*tm(itm)-degtorad(uphs(ilon,ilat,idep)));
part_v = -wm2*vamp(ilon,ilat,idep)*sin(wm2*tm(itm)-degtorad(vphs(ilon,ilat,idep)));
veq_LHS(ilon,ilat,idep) = part_u+part_v;
veq_RHS(ilon,ilat,idep) = -(new_po(ilon,ilat,idep)-new_po(ilon,ilat+1,idep))/(rhoref*dlyv(ilon,ilat));
clear f part_u part_v;
end
end
end
% -- 3 -- %
% -- make plots -- %
temH = reshape(sddpo,[],1);
depnum = [9 24];
rdep(depnum)
for idep = depnum
disp(['plot of linear IW theory: idep=' num2str(idep)])
temx = reshape(squeeze(ueq_LHS(:,:,idep)),[],1);
temy = reshape(squeeze(ueq_RHS(:,:,idep)),[],1);
u_xdata = temx(~isnan(temx));
u_ydata = temy(~isnan(temx));
Hu = temH(~isnan(temx));
clear temx temy;
temx = reshape(squeeze(veq_LHS(:,:,idep)),[],1);
temy = reshape(squeeze(veq_RHS(:,:,idep)),[],1);
v_xdata = temx(~isnan(temx));
v_ydata = temy(~isnan(temx));
Hv = temH(~isnan(temx));
clear temx temy;
figure(1);clf;
plot(u_xdata,u_ydata,'.b');
xlabel('LHS','fontsize',15);
ylabel('RHS','fontsize',15);
set(gca,'fontsize',15);
hold on; grid on;
xmin = floor(min(u_xdata.*10^5)); xmax = ceil(max(u_xdata.*10^5));
ymin = floor(min(u_ydata.*10^5)); ymax = ceil(max(u_ydata.*10^5));
ax_max = max([xmax ymax])
ax_min = min([xmin ymin])
clear xmin xmax ymin ymax;
lmin = ax_min;
lmax = ax_max;
plot([lmin:lmax].*10^(-5),[lmin:lmax].*10^(-5),'--r','linewidth',1.5);
if idep == 9
axis([-2 4 -2 4].*10^(-5));
else
axis([-1 1.5 -1 1.5].* 10^(-5));
end
saveas(gcf,['figures/proc3_linear_theory_paper/ueq_balance_linear_IW_theory_' num2str(round(rdep(idep))) 'm_' regstr '.jpg']);
clear ax_max ax_min lmin lmax; % ymin ymax;
figure(1);clf;
plot(v_xdata,v_ydata,'.b');
xlabel('LHS','fontsize',15);
ylabel('RHS','fontsize',15);
set(gca,'fontsize',15);
hold on; grid on;
xmin = floor(min(v_xdata*10^5)); xmax = ceil(max(v_xdata*10^5));
ymin = floor(min(v_ydata*10^5)); ymax = ceil(max(v_ydata*10^5));
ax_max = max([xmax ymax])
ax_min = min([xmin ymin])
clear xmin xmax ymin ymax;
lmin = ax_min;
lmax = ax_max;
plot([lmin:lmax].*10^(-5),[lmin:lmax].*10^(-5),'--r','linewidth',1.5);
if idep == 9
axis([-2.5 4 -2.5 4].*10^(-5));
else
axis([-1 1.5 -1 1.5].*10^(-5));
end
saveas(gcf,['figures/proc3_linear_theory_paper/veq_balance_linear_IW_theory_' num2str(round(rdep(idep))) 'm_' regstr '.jpg']);
clear ax_max ax_min lmin lmax;
clear u_*data v_*data;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.2: Amplitude of M2 barotropic tide
........................................................
Script for Fig.2a
run script:
........................................................
#!/bin/ksh
# get the speed of (U, V)
cd /scratch/m/m221050
id=/work/mh0256/m300212/ITgeneration_Oct2018/proc1_BartBarc_UV/data/proc1_BartBarc_AmpPhs
u=${id}/Ubart_m2_2012Jan_stormtide-ncep.nc
v=${id}/Vbart_m2_2012Jan_stormtide-ncep.nc
g=/pool/data/MPIOM/TP6M/TP6Ms.nc
# 1. get amplitude and phase of barotropic velocity
cdo -f ext -b f32 selvar,m2amp $u au
cdo -f ext -b f32 mulc,0.0174532 -selvar,m2phase $u pu # 0.0174532=pi/180
cdo -f ext -b f32 selvar,m2amp $v av
cdo -f ext -b f32 mulc,0.0174532 -selvar,m2phase $v pv # 0.0174532=pi/180
# 2. interpolate these amplitudes and phases of barotropic velocities (au, pu, av, pv) to scalar points
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=1
nt=1
ifile=au
ofile=ampu
/
EOF
/home/mpim/m221050/zhuhua/paper/u2s.x
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=1
nt=1
ifile=pu
ofile=phiu
/
EOF
/home/mpim/m221050/zhuhua/paper/u2s.x
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=1
nt=1
ifile=av
ofile=ampv
/
EOF
/home/mpim/m221050/zhuhua/paper/v2s.x
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=1
nt=1
ifile=pv
ofile=phiv
/
EOF
/home/mpim/m221050/zhuhua/paper/v2s.x
cdo -f nc setgrid,${g} -setgrid,r3600x2394 ampu ampu.nc
cdo -f nc setgrid,${g} -setgrid,r3600x2394 ampv ampv.nc
cdo -f nc setgrid,${g} -setgrid,r3600x2394 phiu phiu.nc
cdo -f nc setgrid,${g} -setgrid,r3600x2394 phiv phiv.nc
# 3. get amplitude of barotropic tide
cdo sqrt -add -sqr ampu.nc -sqr ampv.nc ampuv.nc
exit
...............................................
Fortran progran: u2s.f90
(for interpolating u to scale points)
...............................................
PROGRAM u2s
! 13.5.19
! interpolating u to scalar point
!
!
IMPLICIT NONE
! INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12,307)
CHARACTER(180) :: ifile,ofile
INTEGER :: i,j,k,t,nc
INTEGER :: ie,je,ke,nt,id(4)
REAL, ALLOCATABLE :: u(:,:),us(:,:)
real,parameter :: dft=-9.e+33
!
! read parameters, input and output, and allocate the variables
!
NAMELIST /CTL/ ie, je, ke, nt, ifile, ofile
OPEN(10,FILE='INPUT',STATUS='UNKNOWN', &
& ACCESS='SEQUENTIAL',FORM='FORMATTED')
READ(10,CTL)
close(10)
write(0,*) ie, je, ke, nt, ifile,ofile
ALLOCATE (u(ie,je),us(ie,je))
OPEN (unit=11,file=ifile,form='unformatted',status='old')
OPEN (unit=21,file=ofile,form='unformatted',status='unknown')
!
!
do t=1,nt
do k=1,ke
read(11) id
read(11) u
print *, ' data read:', id
us=dft
do j=1,je
do i=2,ie
if(u(i,j).ne.dft .and. u(i-1,j).ne.dft) then
us(i,j)=(u(i,j)+u(i-1,j))/2
endif
enddo
us(1,j)=us(2,j)
enddo
write(21) id
write(21) us
enddo
enddo
end PROGRAM u2s
...............................................
Fortran program: v2s.f90
(for interpolating v to scale points)
...............................................
PROGRAM v2s
! 13.5.19
! interpolating v to scalar point
!
!
IMPLICIT NONE
! INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12,307)
CHARACTER(180) :: ifile,ofile
INTEGER :: i,j,k,t,nc
INTEGER :: ie,je,ke,nt,id(4)
REAL, ALLOCATABLE :: v(:,:),vs(:,:)
real,parameter :: dft=-9.e+33
!
! read parameters, input and output, and allocate the variables
!
NAMELIST /CTL/ ie, je, ke, nt, ifile, ofile
OPEN(10,FILE='INPUT',STATUS='UNKNOWN', &
& ACCESS='SEQUENTIAL',FORM='FORMATTED')
READ(10,CTL)
close(10)
write(0,*) ie, je, ke, nt, ifile,ofile
ALLOCATE (v(ie,je),vs(ie,je))
OPEN (unit=11,file=ifile,form='unformatted',status='old')
OPEN (unit=21,file=ofile,form='unformatted',status='unknown')
!
!
do t=1,nt
do k=1,ke
read(11) id
read(11) v(:,:)
vs=dft
do i=1,ie
do j=2,je
if(v(i,j).ne.dft .and. v(i,j-1).ne.dft) then
vs(i,j)=(v(i,j)+v(i,j-1))/2
endif
enddo
vs(i,1)=vs(i,2)
enddo
write(21) id
write(21) vs
enddo
enddo
end PROGRAM v2s
...............................................
Script for Fig.2b (Matlab)
...............................................
clear all; clc;
fname = 'ampuv_setctomiss-0.nc';
UH = ncread(fname,'var1');
clear fname;
fname = 'depto.nc';
temH = ncread(fname,'depto');
depto = temH(2:end-1,:);
clear fname temH;
Hrng = 0 :100 :6000;
for iH = 1 :length(Hrng)-1
H1 = Hrng(iH);
H2 = Hrng(iH+1);
temH = depto;
temH(temH<H1) = nan;
temH(temH>=H2) = nan;
temUH = UH;
temUH(isnan(temH)) = nan;
avg_UH(iH)= nanmean(nanmean(temUH));
avg_H(iH) = nanmean(nanmean(temH));
clear H1 H2 temH temUH;
end
clf;
plot(avg_UH,-avg_H,'-k','linewidth',1.5);
grid on;
set(gca,'fontsize',15);
xlabel('amplitude of U_H (m/s)','fontsize',15);
ylabel('depth (m)','fontsize',15);
pbaspect([1 3 1]);
saveas(gcf,'averaged-ampuv_over_depth.pdf');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.3:
..........................................................................................
script: for selecting a small region from a global field
- the same script is used to get for other variables in the Emperor Trench region
- for Fig.3a: apply the script below on the bottom topography available from the model description files
- for Fig.3b: apply the script below on the global field obtained from the script for Fig.2
- for Fig.3c: apply the script below on the global field obtained from the script for Fig.5
..........................................................................................
..........................................
Select regional data
..........................................
#!/bin/ksh
# select data in Hawaii Ridge rgion
#
cd /scratch/m/m221050
od=/work/mh0256/m221050/zhuhua/paper
# set longitudes and latitudes of the selected region
lo1=190
lo2=199
la1=22
la2=26
cdo selvar,depto /pool/data/MPIOM/TP6M/TP6ML40_fx.nc depto.nc
f=depto.nc
x=${od}/${f}.nc
of=${od}/${f}_box_${lo1}_${lo2}_${la1}_${la2}.nc
cdo sellonlatbox,${lo1},${lo2},${la1},${la2} ${x} ${of}
exit
..........................................
Script for Fig.3a (NCL)
..........................................
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin
;;read u,v,t from the data at 500hPa
Hfile = addfile("data/depto_box_190_199_22_26.nc","r")
depto = Hfile->depto(0,0,:,:)
depto@lon2d = Hfile->lon
depto@lat2d = Hfile->lat
printVarSummary(depto)
;;create plots;;
wks = gsn_open_wks("pdf","figures/depto_Hawaii")
res = True
res@gsnDraw = False
res@gsnFrame = False
res@gsnLeftString = ""
res@gsnRightString = ""
;;set map;;
mpres = res
mpres@mpDataSetName = "Earth..4"
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpOutlineOn = True
mpres@mpGeophysicalLineThicknessF = 2
mpres@mpFillDrawOrder = "PostDraw"
;;set area;;
mpres@mpMinLatF = 22
mpres@mpMaxLatF = 26
mpres@mpMinLonF = 190
mpres@mpMaxLonF = 199
mpres@gsnMajorLatSpacing = 1
mpres@gsnMajorLonSpacing = 1
;;set contour;;
cnres = res
cnres@cnFillDrawOrder = "PreDraw"
cnres@cnFillOn = True
cnres@cnLinesOn = False
cnres@lbLabelBarOn = True
cnres@lbOrientation = "Vertical"
cnres@lbLabelFontHeightF = 0.008
cnres@lbLabelStride = 2
cnres@pmLabelBarOrthogonalPosF = 0.02
cnres@pmLabelBarWidthF = 0.06
cnres@cnLevelSelectionMode = "ManualLevels"
cnres@cnMinLevelValF = 0
cnres@cnMaxLevelValF = 6400
cnres@cnLevelSpacingF = 200
cnres@cnLabelBarEndStyle = "ExcludeOuterBoxes"
cmap = read_colormap_file("temp_19lev")
cnres@cnFillPalette = cmap(::-1,:)
cnres@gsnLeftString = ""
;; set contour lines for depto
clres = res
clres@cnFillOn = False
clres@cnLinesOn = True
clres@cnLineThicknesses = 2.5
clres@cnLabelMasking = True
clres@cnLineLabelBackgroundColor = "transparent"
clres@cnLineLabelFontHeightF = 0.008
clres@cnInfoLabelOn = False
clres@cnLineColor = "gray10"
clres@cnLevelSelectionMode = "ManualLevels"
clres@cnMinLevelValF = 200
clres@cnMaxLevelValF = 6000
clres@cnLevelSpacingF = 500
;;plot;;
map = gsn_csm_map(wks,mpres)
contour = gsn_csm_contour(wks,depto,cnres)
cnlines = gsn_csm_contour(wks,depto,clres)
overlay(map,contour)
overlay(map,cnlines)
draw(map)
frame(wks)
end
..........................................
Script for Fig.3b (NCL)
..........................................
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin
;;read u,v,t from the data at 500hPa
in = addfile("data/ampuv_box_190_199_22_26.nc","r")
amp = in->var1
amp@lon2d=in->lon
amp@lat2d=in->lat
printVarSummary(amp)
Hfile = addfile("data/depto_box_190_199_22_26.nc","r")
depto = Hfile->depto(0,0,:,:)
depto@lon2d = Hfile->lon
depto@lat2d = Hfile->lat
printVarSummary(depto)
;;create plots;;
wks = gsn_open_wks("pdf","figures/amp-UV_Hawaii_new_v2")
res = True
res@gsnDraw = False
res@gsnFrame = False
res@gsnLeftString = ""
res@gsnRightString = ""
;;set map;;
mpres = res
mpres@mpDataSetName = "Earth..4"
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpOutlineOn = True
mpres@mpGeophysicalLineThicknessF = 2
mpres@mpFillDrawOrder = "PostDraw"
;;set area;;
mpres@mpMinLatF = 22
mpres@mpMaxLatF = 26
mpres@mpMinLonF = 190
mpres@mpMaxLonF = 199
mpres@gsnMajorLatSpacing = 1
mpres@gsnMajorLonSpacing = 1
;;set contour;;
cnres = res
cnres@cnFillDrawOrder = "PreDraw"
cnres@cnFillOn = True
cnres@cnLinesOn = False
cnres@lbLabelBarOn = True
cnres@lbOrientation = "Vertical"
cnres@lbLabelFontHeightF = 0.008
cnres@lbLabelStride = 2
cnres@pmLabelBarOrthogonalPosF = 0.02
cnres@pmLabelBarWidthF = 0.06
cnres@cnLevelSelectionMode = "ManualLevels"
cnres@cnMinLevelValF = 0.02
cnres@cnMaxLevelValF = 0.4;35
cnres@cnLevelSpacingF = 0.02
cmap = read_colormap_file("WhiteBlueGreenYellowRed.rgb")
cnres@cnFillPalette = cmap
cnres@gsnLeftString = ""
;; set contour lines for depto
clres = res
clres@cnFillOn = False
clres@cnLinesOn = True
clres@cnLineThicknesses = 2.5 ;1.5
clres@cnLabelMasking = True
clres@cnLineLabelBackgroundColor = "transparent"
clres@cnLineLabelFontHeightF = 0.008
clres@cnInfoLabelOn = False
clres@cnLineColor = "gray10" ;"slategray" ;"skyblue4" ;"steelblue4" ;"slategray" ;"snow4" ; "black" ; "white"
clres@cnLevelSelectionMode = "ManualLevels"
clres@cnMinLevelValF = 200
clres@cnMaxLevelValF = 6000
clres@cnLevelSpacingF = 500
;;plot;;
map = gsn_csm_map(wks,mpres)
contour = gsn_csm_contour(wks,amp,cnres)
cnlines = gsn_csm_contour(wks,depto,clres)
overlay(map,contour)
overlay(map,cnlines)
draw(map)
frame(wks)
end
..........................................
Script for Fig.3c (NCL)
..........................................
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin
;;read u,v,t from the data at 500hPa
in = addfile("data/abs_nabla_d_box_190_199_22_26.nc","r")
slp = in->var193(0,:,:)
slp@lon2d=in->lon
slp@lat2d=in->lat
printVarSummary(slp)
Hfile = addfile("data/depto_box_190_199_22_26.nc","r")
depto = Hfile->depto(0,0,:,:)
depto@lon2d = Hfile->lon
depto@lat2d = Hfile->lat
printVarSummary(depto)
;;create plots;;
wks = gsn_open_wks("pdf","figures/abs-nabla-d_Hawaii_new")
res = True
res@gsnDraw = False
res@gsnFrame = False
res@gsnLeftString = ""
res@gsnRightString = ""
;;set map;;
mpres = res
mpres@mpDataSetName = "Earth..4"
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpOutlineOn = True
mpres@mpGeophysicalLineThicknessF = 2
mpres@mpFillDrawOrder = "PostDraw"
;;set area;;
mpres@mpMinLatF = 22
mpres@mpMaxLatF = 26
mpres@mpMinLonF = 190
mpres@mpMaxLonF = 199
mpres@gsnMajorLatSpacing = 1
mpres@gsnMajorLonSpacing = 1
;;set contour;;
cnres = res
cnres@cnFillDrawOrder = "PreDraw"
cnres@cnFillOn = True
cnres@cnLinesOn = False
cnres@lbLabelBarOn = True
cnres@lbOrientation = "Vertical"
cnres@lbLabelFontHeightF = 0.008
cnres@lbLabelStride = 3
cnres@pmLabelBarOrthogonalPosF = 0.02
cnres@pmLabelBarWidthF = 0.06
cnres@cnLevelSelectionMode = "ExplicitLevels"
cnres@cnLevels = (/0.001,0.005,0.007,0.009,0.011,0.013,0.015,0.017,0.019,0.021,0.023,0.025,0.027,0.029,0.031,0.033,0.035,0.037,0.039,0.041,0.043,0.045,0.047,0.049,0.051,0.053,0.055,0.057,0.059,0.061,0.063,0.065,0.067,0.069,0.071,0.073,0.075,0.077,0.079,0.081,0.083,0.085,0.087,0.089,0.091,0.093,0.095,0.099,0.1,0.5/)
cnres@cnLabelBarEndStyle = "ExcludeOuterBoxes"
cmap = read_colormap_file("MPL_RdBu")
cnres@cnFillPalette = cmap(::-1,:)
cnres@gsnLeftString = ""
;; set contour lines for depto
clres = res
clres@cnFillOn = False
clres@cnLinesOn = True
clres@cnLineThicknesses = 2.5 ;1.5
clres@cnLabelMasking = True
clres@cnLineLabelBackgroundColor = "transparent"
clres@cnLineLabelFontHeightF = 0.008
clres@cnInfoLabelOn = False
clres@cnLineColor = "gray10" ;"slategray" ;"skyblue4" ;"steelblue4" ;"slategray" ;"snow4" ; "black" ; "white"
clres@cnLevelSelectionMode = "ManualLevels"
clres@cnMinLevelValF = 200
clres@cnMaxLevelValF = 6000
clres@cnLevelSpacingF = 500
;;plot;;
map = gsn_csm_map(wks,mpres)
contour = gsn_csm_contour(wks,slp,cnres)
cnlines = gsn_csm_contour(wks,depto,clres)
overlay(map,contour)
overlay(map,cnlines)
draw(map)
frame(wks)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.4:
.................................................................
- for Fig.4a: use the script for Fig.3a, but for the Emperor Trench region with
lo1=174
lo2=180
la1=33
la2=45
- for Fig.4b: use the script for Fig.3b, but for the Emperor region
- Fig.4c: apply the following script for the Emperor Trench region on the global field obtained from the script of Fig.10
................................................................
..........................................
Script for Fig.4c (NCL)
..........................................
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin
;;read u,v,t from the data at 500hPa
in = addfile("data/Px+Py_minus_mulc-half_box_174_180_33_45.nc","r")
pgen = in->var71(0,0,:,:)
pgen@lon2d=in->lon
pgen@lat2d=in->lat
printVarSummary(pgen)
Hfile = addfile("data/depto_box_174_180_33_45.nc","r")
depto = Hfile->depto(0,0,:,:)
depto@lon2d = Hfile->lon
depto@lat2d = Hfile->lat
printVarSummary(depto)
;;create plots;;
wks = gsn_open_wks("pdf","figures/M2-ITgen_Emperor_Trough_paper_new")
res = True
res@gsnDraw = False
res@gsnFrame = False
res@gsnLeftString = ""
res@gsnRightString = ""
;;set map;;
mpres = res
mpres@mpDataSetName = "Earth..4"
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpOutlineOn = True
mpres@mpGeophysicalLineThicknessF = 2
mpres@mpFillDrawOrder = "PostDraw"
;;set area;;
mpres@mpMinLatF = 33
mpres@mpMaxLatF = 45
mpres@mpMinLonF = 174
mpres@mpMaxLonF = 180
mpres@gsnMajorLatSpacing = 5
mpres@gsnMajorLonSpacing = 5
;;set contour;;
cnres = res
cnres@cnFillDrawOrder = "PreDraw"
cnres@cnFillOn = True
cnres@cnLinesOn = False
cnres@lbLabelBarOn = True
cnres@lbOrientation = "Vertical"
cnres@lbLabelFontHeightF = 0.02;0.008
cnres@lbLabelStride = 2
cnres@pmLabelBarOrthogonalPosF = 0.05 ;0.02
cnres@pmLabelBarWidthF = 0.1 ;0.06
cnres@cnLevelSelectionMode = "ManualLevels"
cnres@cnMinLevelValF = -0.016
cnres@cnMaxLevelValF = 0.016
cnres@cnLevelSpacingF = 0.002
cmap = read_colormap_file("MPL_RdBu")
cnres@cnFillPalette = cmap(::-1,:)
cnres@gsnLeftString = ""
;; set contour lines for depto
clres = res
clres@cnFillOn = False
clres@cnLinesOn = True
clres@cnLineThicknesses = 2.5 ;1.5
clres@cnLabelMasking = True
clres@cnLineLabelBackgroundColor = "transparent"
clres@cnLineLabelFontHeightF = 0.012
clres@cnInfoLabelOn = False
clres@cnLineColor = "gray30" ;"slategray" ;"skyblue4" ;"steelblue4" ;"slategray" ;"snow4" ; "black" ; "white"
clres@cnLevelSelectionMode = "ManualLevels"
clres@cnMinLevelValF = 200
clres@cnMaxLevelValF = 6000
clres@cnLevelSpacingF = 500
;;plot;;
map = gsn_csm_map(wks,mpres)
contour = gsn_csm_contour(wks,pgen,cnres)
cnlines = gsn_csm_contour(wks,depto,clres)
overlay(map,contour)
overlay(map,cnlines)
draw(map)
frame(wks)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.5:
..................................................................................
script: using dudx.f90 and dudy.f90 to calculate the gradient of bottom topography
..................................................................................
#!/bin/ksh
# calculate the nabla of the bottom topography
#
cd /scratch/m/m221050
id=/work/mh0256/m300212/ITgeneration_Oct2018/proc1_BartBarc_UV/data/proc1_BartBarc_AmpPhs
od=/work/mh0256/m221050/zhuhua/paper
odx=${od}/dddx.nc
ody=${od}/dddy.nc
of=${od}/abs_nabla_d.nc
g=/pool/data/MPIOM/TP6M/TP6Ms.nc
cdo selvar,deuto /pool/data/MPIOM/TP6M/TP6ML40_fx.nc deuto.nc
cdo selvar,deute /pool/data/MPIOM/TP6M/TP6ML40_fx.nc deute.nc
cdo selvar,dlxu /pool/data/MPIOM/TP6M/TP6ML40_fx.nc dlxu.nc
cdo selvar,dlyv /pool/data/MPIOM/TP6M/TP6ML40_fx.nc dlyv.nc
# calculate the 2-dimensional nablace of depth
cdo -f ext copy deuto.nc deuto
cdo -f ext copy deute.nc deute
cdo -f ext copy dlxu.nc dlxu
cdo -f ext copy dlyv.nc dlyv
cat >INPUT<<EOF
&CTL
ie=3602
je=2394
nt=1
ifile=dlxu,deuto
ofile=dddx
/
EOF
/home/mpim/m221050/zhuhua/paper/dudx.x
cat >INPUT<<EOF
&CTL
ie=3602
je=2394
nt=1
ifile=dlyv,deute
ofile=dddy
/
EOF
/home/mpim/m221050/zhuhua/paper/dudy.x
cdo -f nc setgrid,${g} -selindexbox,2,3601,1,2394 -setgrid,r3602x2394 dddx dddx.nc
cdo -f nc setgrid,${g} -selindexbox,2,3601,1,2394 -setgrid,r3602x2394 dddy dddy.nc
cp dddx.nc ${odx}
cp dddy.nc ${ody}
cdo sqrt -add -sqr ${odx} -sqr ${ody} ${of}
.........................
fortran programm dudx.f90
.........................
PROGRAM u_gradient
! 13.5.19
! calculating dudx. u can be any variable (e.g. a scalar variable or a component of vector u).
!
!
!
!
IMPLICIT NONE
! INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12,307)
CHARACTER(180) :: ifile(2),ofile
INTEGER :: i,j,nc,t
INTEGER :: ie,je,nt,id(4)
REAL, ALLOCATABLE :: u(:,:)
REAL, ALLOCATABLE :: dx(:,:),dudx(:,:)
real,parameter :: dft=-9.e+33
!
! read parameters, input and output, and allocate the variables
!
NAMELIST /CTL/ ie, je, nt, ifile, ofile
OPEN(10,FILE='INPUT',STATUS='UNKNOWN', &
& ACCESS='SEQUENTIAL',FORM='FORMATTED')
READ(10,CTL)
close(10)
write(0,*) ie, je, nt, ifile, ofile
ALLOCATE (u(ie,je),dudx(ie,je),dx(ie,je))
OPEN (unit=11,file=ifile(1),form='unformatted',status='old')
OPEN (unit=12,file=ifile(2),form='unformatted',status='old')
OPEN (unit=21,file=ofile,form='unformatted',status='unknown')
!
! read grid distance
!
read(11) id
read(11) dx
print *, ' read grid x-distance:', id
do t=1,nt
read(12) id
read(12) u
dudx=dft
do j=1,je
do i=2,ie
if(u(i,j).ne.dft .and. u(i-1,j).ne.dft) then
if(dx(i,j).ne.dft) then
dudx(i,j)=(u(i,j)-u(i-1,j))/dx(i,j)
! print '(2i5,4e15.6)', i,j,u(i,j),u(i-1,j), u(i,j)-u(i-1,j),dx(i,j)
end if
endif
enddo
dudx(1,j)=dudx(2,j)
enddo
write(21) id
write(21) dudx
enddo
end PROGRAM u_gradient
..........................
fortran programm dudy.f90
..........................
PROGRAM u_gradient
! 13.5.19
! calculating dudy. u can be any variables (e.g. a scalar variable or a component of vector u).
!
! FROM div_u.f90: caution: grid distances from arcgri are positive, i.e. they result from
! northern latitude - southern latitude
! Since the model output start from the most northern latitude, the meridional gradients
! must be calculated from difference defined by
! u(i,j-1,k)-u(i,j,k)
!
!
IMPLICIT NONE
! INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12,307)
CHARACTER(180) :: ifile(2),ofile
INTEGER :: i,j,t,nc
INTEGER :: ie,je,nt,id(4)
REAL, ALLOCATABLE :: u(:,:)
REAL, ALLOCATABLE :: dy(:,:),dudy(:,:)
real,parameter :: dft=-9.e+33
!
! read parameters, input and output, and allocate the variables
!
NAMELIST /CTL/ ie, je, nt, ifile, ofile
OPEN(10,FILE='INPUT',STATUS='UNKNOWN', &
& ACCESS='SEQUENTIAL',FORM='FORMATTED')
READ(10,CTL)
close(10)
write(0,*) ie, je, nt, ifile,ofile
ALLOCATE (u(ie,je),dudy(ie,je),dy(ie,je))
OPEN (unit=11,file=ifile(1),form='unformatted',status='old')
OPEN (unit=12,file=ifile(2),form='unformatted',status='old')
OPEN (unit=21,file=ofile,form='unformatted',status='unknown')
!
! read grid distance
!
read(11) id
read(11) dy
print *, ' read grid x-distance:', id
do t=1,nt
read(12) id
read(12) u
dudy=dft
do i=1,ie
do j=2,je
if(u(i,j).ne.dft .and. u(i,j-1).ne.dft) then
if(dy(i,j).ne.dft) then
dudy(i,j)=(u(i,j-1)-u(i,j))/dy(i,j)
end if
endif
enddo
dudy(i,1)=dudy(i,2)
enddo
write(21) id
write(21) dudy
enddo
end PROGRAM u_gradient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.6:
.................................................................................................
script for calculating pressure perturbation p' from density perturbation following Eq.(7)
.................................................................................................
#!/bin/ksh
# calculating amplitude and phase of p-prime using:
# 1. harmonic analysis result of rho (i.e. amplitude and phase)
# 2. surface elevation, which is the upper bound of the vertical integral of the hydrostatic relation
#
cd /scratch/m/m221050
drho=/work/mh0256/m300212/STORMTIDE_run/Tides_NCEP_SAL2/Data_analysis/proc2_harmonics_velocity_density/rho_total
frho=rho_harmonics_removeHalo.nc
dssh=/work/mh0256/m300212/STORMTIDE_run/Tides_NCEP_SAL2/experiments/hel17153-DZA_TP6ML40_mpiom-trunk_2012Jan/outdata
fssh=hel17153-DZA_TP6ML40_mpiom-trunk_mpiom_data_zos_1hr_2012-01-01_2012-01-31.nc
od=/work/mh0256/m221050/zhuhua/paper
of1=${od}/amp_p.nc
of2=${od}/phi_p.nc
of3=${od}/amp_p_eta.nc
of4=${od}/phi_p_eta.nc
g=/pool/data/MPIOM/TP6M/TP6Ms.nc
cdo -f ext selvar,m2amp ${drho}/${frho} ampr0
cdo -f ext mulc,0.0174532 -selvar,m2phase ${drho}/${frho} phir0 # 0.0174532=pi/180
cdo setctomiss,0 ampr0 ampr
cdo setctomiss,0 phir0 phir
cdo info ampr
cdo info phir
## get zo at one time step
cdo -f ext selindexbox,2,3601,1,2394 -seltimestep,1 ${dssh}/${fssh} zo1
cdo setctomiss,0 zo1 zo
# using fortran programm ftn90/pprim.f90
# ifort -o pprim.x ../ftn90/pprim.f90
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=40
ifile=ampr,phir,zo
ofile=ampp,phip,ampe,phie
/
&DEPTH
tiestw = 0, 20, 30, 40, 50, 60, 70, 83, 98, 118,
143, 173, 208, 248, 293, 343, 398, 458, 528, 608,
698, 798, 908, 1028, 1158, 1298, 1448, 1618, 1798, 1988,
2188, 2408, 2658, 2928, 3228, 3578, 3978, 4428, 4928, 5428, 6028,
/
EOF
/home/mpim/m221050/zhuhua/paper/pprim.x
cdo -f nc setgrid,${g} -setgrid,r3600x2394 ampp ampp.nc
cp ampp.nc ${of1}
cdo -f nc setgrid,${g} -setgrid,r3600x2394 phip phip.nc
cp phip.nc ${of2}
cdo -f nc setgrid,${g} -setgrid,r3600x2394 ampe ampe.nc
cp ampe.nc ${of3}
cdo -f nc setgrid,${g} -setgrid,r3600x2394 phie phie.nc
cp phie.nc ${of4}
exit
..................................
fortran program for calculating p'
..................................
program p
! 7.5.2019
!
! calculating the pressure perturbation p' in form of the wave ansatz
! p(x,y,z,t)=ampp*cos(2*pi*omega*t+phip)
! by vertically integrating
! rho(x,y,z,t)=ampr*cos(2*pi*omega*t+phir)
! fowlloing the hydrostatic relation from z=-d to z=\eta
!
! ampp^2 = a^2+b^2,
! where a = g int_z^0 ampr *cos(phir) dz
! b = g int_z^0 ampr * sin(phir) dz
!
! phip=arctan (b/a)
!
! Two amplitudes /phases are calculated: the first one (ampp, phip) is obtained by integrating rho' from 0 to z,
! the second (ampe, phie) from 0 to eta.
!
! Note: the second is obtained to one given instance of eta
!
! input : ampr on fort.11
! phir on fort.12
! output : ampp, ampe on fort.21
! phip, phie on fort.22
!
! The program is based on mo_convection.f90
!
IMPLICIT NONE
CHARACTER(128) :: ifile(3),ofile(4)
integer :: j,je,i,ie,k,ke,kep,id1(4),id2(4)
real, allocatable :: ampr(:,:,:),phir(:,:,:),ampp(:,:,:),phip(:,:,:),ampe(:,:),phie(:,:)
real, allocatable :: a(:,:,:),b(:,:,:),eta(:,:)
real, allocatable :: dz(:), dzw(:), tiestu(:), tiestw(:)
real :: w
real,parameter :: dft=-9.e+33,g=9.81,pi=3.1416
NAMELIST /CTL/ ie, je, ke, ifile, ofile
NAMELIST /DEPTH/ tiestw
OPEN(10,FILE='INPUT',STATUS='UNKNOWN', &
& ACCESS='SEQUENTIAL',FORM='FORMATTED')
READ(10,CTL)
write(0,*) ie, je, ke, ifile, ofile
kep=ke+1
ALLOCATE(tiestw(kep))
READ(10,DEPTH)
write(0,*) tiestw
close(10)
ALLOCATE (DZ(KEP),TIESTU(KEP),dzw(kep))
ALLOCATE(ampr(ie,je,ke),phir(ie,je,ke),ampp(ie,je,ke),phip(ie,je,ke))
allocate (a(ie,je,ke),b(ie,je,ke),eta(ie,je),ampe(ie,je),phie(ie,je))
OPEN (unit=11,file=ifile(1),form='unformatted', status='old')
OPEN (unit=12,file=ifile(2),form='unformatted', status='old')
OPEN (unit=13,file=ifile(3),form='unformatted', status='old')
OPEN (unit=21,file=ofile(1),form='unformatted', status='unknown')
OPEN (unit=22,file=ofile(2),form='unformatted', status='unknown')
OPEN (unit=23,file=ofile(3),form='unformatted', status='unknown')
OPEN (unit=24,file=ofile(4),form='unformatted', status='unknown')
! calculating dz using dzw
do k=1,ke
dzw(k)=tiestw(k+1)-tiestw(k)
enddo
TIESTU(1) = 0.5 * DZW(1)
DO K=2,KE
TIESTU(K) = 0.5 * ( TIESTW(K+1) + TIESTW(K) )
END DO
TIESTU(KEP)=9000.
DZ(1) = TIESTU(1)
print*, 1, TIESTU(1),tiestw(1),dz(1),dzw(1)
DO K=2,KEP
DZ(K) = TIESTU(K) - TIESTU(K-1)
print*, k, TIESTU(k),tiestw(k),dz(k),dzw(k)
END DO
! calculating pou and pow at each time step
ampp=dft
phip=dft
do k=1,ke
read(11) id1
read(11) ampr(:,:,k)
read(12) id2
read(12) phir(:,:,k)
enddo
read(13) id2
read(13) eta
print *, ' data considered:', id1, id2
! vertical integrating the hydrostatic relation
a=dft
b=dft
ampp=dft
phip=dft
ampe=dft
phie=dft
DO j=1,je
DO i=1,ie
! p' obtained by integrating rho' from z=0 to z=eta
if(eta(i,j).ne.dft .and. phir(i,j,1).ne.dft .and. ampr(i,j,1).ne.dft) then
ampe(i,j)=g*eta(i,j)*ampr(i,j,1)
phie(i,j)=phir(i,j,1)
endif
! calculating the vertical integral from z=z to z=0
if(ampr(i,j,1).ne.dft .and. phir(i,j,1).ne.dft) then
a(i,j,1)=g*tiestu(1)*ampr(i,j,1)*cos(phir(i,j,1))
b(i,j,1)=g*tiestu(1)*ampr(i,j,1)*sin(phir(i,j,1))
ampp(i,j,1)=sqrt(a(i,j,1)*a(i,j,1)+b(i,j,1)*b(i,j,1))
endif
do k=2,ke
if (ampr(i,j,k).ne.dft .and. a(i,j,k-1).ne.dft .and. phir(i,j,k).ne.dft .and. b(i,j,k-1).ne.dft) then
a(i,j,k) = a(i,j,k-1) + 0.5*g*(ampr(i,j,k)*cos(phir(i,j,k))+ampr(i,j,k-1)*cos(phir(i,j,k-1)))*dz(k)
b(i,j,k) = b(i,j,k-1) + 0.5*g*(ampr(i,j,k)*sin(phir(i,j,k))+ampr(i,j,k-1)*sin(phir(i,j,k-1)))*dz(k)
ampp(i,j,k)=sqrt(a(i,j,k)*a(i,j,k)+b(i,j,k)*b(i,j,k))
endif
enddo ! do k
do k=1,ke
if(a(i,j,k).ne.dft .and. b(i,j,k).ne.dft) then
w=atan2(b(i,j,k),a(i,j,k))
if(b(i,j,k).gt.0 .and. a(i,j,k).gt.0) then
phip(i,j,k)=w
!print '(a,3i5,4f10.2)','1.quadrat',i,j,k,a(i,j,k),b(i,j,k),w,phip(i,j,k)
endif
if(b(i,j,k).lt.0 .and. a(i,j,k).gt.0) then
phip(i,j,k)=w+2*pi
!print '(a,3i5,4f10.2)','4.quadrat',i,j,k,a(i,j,k),b(i,j,k),w,phip(i,j,k)
endif
if(b(i,j,k).gt.0 .and. a(i,j,k).lt.0) then
phip(i,j,k)=w
!print '(a,3i5,4f10.2)', '2.quadrat',i,j,k,a(i,j,k),b(i,j,k),w,phip(i,j,k)
endif
if(b(i,j,k).lt.0 .and. a(i,j,k).lt.0) then
phip(i,j,k)=w+2*pi
!print '(a,3i5,4f10.2)', '3.quadrat',i,j,k,a(i,j,k),b(i,j,k),w,phip(i,j,k)
endif
if(a(i,j,k).eq.0) then
if(b(i,j,k).gt.0) phip(i,j,k)=pi*0.5
if(b(i,j,k).lt.0) phip(i,j,k)=-pi*0.5
endif
if(b(i,j,k).eq.0) then
if(a(i,j,k).gt.0) phip(i,j,k)=0
if(a(i,j,k).lt.0) phip(i,j,k)=-pi
endif
endif
enddo
enddo ! do i
enddo ! do j
! write out
do k=1,ke
write(21) id1(1),61,INT(tiestu(k)),id1(4)
write(21) ampp(:,:,k)
write(22) id1(2),71,INT(tiestu(k)),id2(4)
write(22) phip(:,:,k)
enddo
write(23) id1(1),62,0,id1(4)
write(23) ampe
write(24) id1(2),72,0,id2(4)
write(24) phie
end program p
..................................................................................................
script for calcultating depth-averaged p'
..........................................
#!/bin/ksh
# calculating depth average / integral
#
#
cd /scratch/m/m221050
d=/work/mh0256/m221050/zhuhua/paper
if1=${d}/amp_p.nc
if2=${d}/phi_p.nc
of1=${d}/amp_depthavg-p.nc
of2=${d}/phi_depthavg-p.nc
g=/pool/data/MPIOM/TP6M/TP6Ms.nc
cdo -f ext copy ${if1} ampp
cdo -f ext copy ${if2} phip
cdo info ampp
cdo info phip
# using fortran programm ftn90/depthavg.f90.f90
# ifort -o depthavg.x ../ftn90/depthavg.f90
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=40
lavg=.true.
ifile=ampp,phip
ofile=avgamp,avgphi
/
&DEPTH
tiestw = 0, 20, 30, 40, 50, 60, 70, 83, 98, 118,
143, 173, 208, 248, 293, 343, 398, 458, 528, 608,
698, 798, 908, 1028, 1158, 1298, 1448, 1618, 1798, 1988,
2188, 2408, 2658, 2928, 3228, 3578, 3978, 4428, 4928, 5428, 6028,
/
EOF
/home/mpim/m221050/zhuhua/paper/depthavg.x
cdo -f nc setgrid,${g} -setgrid,r3600x2394 avgamp avgamp.nc
cp avgamp.nc ${of1}
cdo -f nc setgrid,${g} -setgrid,r3600x2394 avgphi avgphi.nc
cp avgphi.nc ${of2}
..............................................
fortran program for calculating depth average
...............................................
program p
! 28.2.2019
!
! Calculating depth integral /average of the complex quantity
! p(x,y,z,t)=amp*exp(i*omega*t+i*phi)
! The result is ampp*exp(i*omega*t+phip) with
! ampp^2 = a^2+b^2,
! where a = int_d^0 amp * cos(phi) dz
! b = int_d^0 amp * sin(phi) dz
!
! phip=arctan (b/a)
!
! input : amp on fort.11
! phi on fort.12
! output : ampp on fort.21
! phip on fort.22
!
!
!
IMPLICIT NONE
CHARACTER(128) :: ifile(2),ofile(2)
integer :: j,je,i,ie,k,ke,kep,id1(4),id2(4)
real, allocatable :: amp(:,:,:),phi(:,:,:),ampp(:,:),phip(:,:)
real, allocatable :: dz(:), dzw(:), tiestu(:), tiestw(:)
real :: w,a,b,d
logical :: lavg
real,parameter :: dft=-9.e+33,g=9.81,pi=3.1416
NAMELIST /CTL/ ie, je, ke, lavg, ifile, ofile
NAMELIST /DEPTH/ tiestw
OPEN(10,FILE='INPUT',STATUS='UNKNOWN', &
& ACCESS='SEQUENTIAL',FORM='FORMATTED')
READ(10,CTL)
write(0,*) ie, je, ke, ifile, ofile
kep=ke+1
ALLOCATE(tiestw(kep))
READ(10,DEPTH)
write(0,*) tiestw
close(10)
ALLOCATE (DZ(KEP),TIESTU(KEP),dzw(kep))
ALLOCATE(amp(ie,je,ke),phi(ie,je,ke),ampp(ie,je),phip(ie,je))
OPEN (unit=11,file=ifile(1),form='unformatted', status='old')
OPEN (unit=12,file=ifile(2),form='unformatted', status='old')
OPEN (unit=21,file=ofile(1),form='unformatted', status='unknown')
OPEN (unit=22,file=ofile(2),form='unformatted', status='unknown')
! calculating dz using dzw
do k=1,kep
dzw(k)=tiestw(k+1)-tiestw(k)
enddo
TIESTU(1) = 0.5 * DZW(1)
DO K=2,KE
TIESTU(K) = 0.5 * ( TIESTW(K+1) + TIESTW(K) )
END DO
TIESTU(KEP)=9000.
DZ(1) = TIESTU(1)
print*, 1, TIESTU(1),tiestw(1),dz(1),dzw(1)
DO K=2,KEP
DZ(K) = TIESTU(K) - TIESTU(K-1)
print*, k, TIESTU(k),tiestw(k),dz(k),dzw(k)
END DO
print '(10f4.0)', dzw
! read in data
do k=1,ke
read(11) id1
read(11) amp(:,:,k)
read(12) id2
read(12) phi(:,:,k)
enddo
print *, ' data considered:', id1, id2
! vertical integration
ampp=dft
phip=dft
DO j=1,je
DO i=1,ie
a=0
b=0
d=0
do k=1,ke
if (amp(i,j,k).ne.dft .and. phi(i,j,k).ne.dft) then
a = a + amp(i,j,k)*cos(phi(i,j,k))*dzw(k)
b = b + amp(i,j,k)*sin(phi(i,j,k))*dzw(k)
d = d + dzw(k)
endif
enddo ! do k
if(d.ne.0) then
if(lavg) then
ampp(i,j)=sqrt(a*a+b*b)/d
else
ampp(i,j)=sqrt(a*a+b*b)
endif
endif
if(d.ne.0) then
w=atan2(b,a)
if(b.gt.0 .and. a.gt.0) then
phip(i,j)=w
endif
if(b.lt.0 .and. a.gt.0) then
phip(i,j)=w+2*pi
endif
if(b.gt.0 .and. a.lt.0) then
phip(i,j)=w
endif
if(b.lt.0 .and. a.lt.0) then
phip(i,j)=w+2*pi
endif
if(a.eq.0) then
if(b.gt.0) phip(i,j)=pi*0.5
if(b.lt.0) phip(i,j)=-pi*0.5
endif
if(b.eq.0) then
if(a.gt.0) phip(i,j)=0
if(a.lt.0) phip(i,j)=-pi
endif
endif
enddo ! do i
enddo ! do j
! write out
write(21) id1(1),61,INT(tiestu(ke)),id1(4)
write(21) ampp
write(22) id1(2),62,INT(tiestu(ke)),id2(4)
write(22) phip
end program p
.....................................................................................................
script for calculating the internal tide pressure p''
..................................................................................
#!/bin/ksh
# calculating ppprim using dif.f90
#
#
cd /scratch/m/m221050
d=/work/mh0256/m221050/zhuhua/paper
if1=${d}/amp_p.nc
if2=${d}/phi_p.nc
if3=${d}/amp_depthavg-p.nc
if4=${d}/phi_depthavg-p.nc
of1=${d}/amp_ppprim.nc
of2=${d}/phi_ppprim.nc
of3=${d}/ppprim-reconstructed_11-22.nc
g=/pool/data/MPIOM/TP6M/TP6Ms.nc
cdo -f ext copy ${if1} ampx
cdo -f ext copy ${if2} phix
cdo -f ext copy ${if3} ampy
cdo -f ext copy ${if4} phiy
cdo info ampx
cdo info ampy
# using fortran programm ftn90/dif.f90
# ifort -o dif.x ../ftn90/dif.f90
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=40
nt=1
ifile=ampx,phix,ampy,phiy
ofile=ampd,phid
/
EOF
/home/mpim/m221050/zhuhua/paper/dif.x
cdo -f nc setgrid,${g} -setgrid,r3600x2394 ampd ampd.nc
cp ampd.nc ${of1}
cdo -f nc setgrid,${g} -setgrid,r3600x2394 phid phid.nc
cp phid.nc ${of2}
exit
.......................
fortran program dif.f90
.......................
PROGRAM diff
! 10.5.2019
! calculating amplitude and phase of the difference of a 3d complex variable x
! and a 2d complex variable y
!
!
!
IMPLICIT NONE
CHARACTER(180) :: ifile(4),ofile(2)
INTEGER :: i,j,k,t,nc
INTEGER :: ie,je,ke,nt,id(4)
integer, allocatable :: idepth(:)
REAL, ALLOCATABLE :: ampx(:,:,:),phix(:,:,:),ampy(:,:),phiy(:,:),ampd(:,:,:),phid(:,:,:)
real,parameter :: dft=-9.e+33,pi=3.1416
real :: a,b,w
!
! read parameters, input and output, and allocate the variables
!
NAMELIST /CTL/ ie, je, ke, nt, ifile, ofile
OPEN(10,FILE='INPUT',STATUS='UNKNOWN', &
& ACCESS='SEQUENTIAL',FORM='FORMATTED')
READ(10,CTL)
close(10)
write(0,*) ie, je, ke, nt, ifile,ofile
ALLOCATE (ampx(ie,je,ke),phix(ie,je,ke),ampy(ie,je),phiy(ie,je),idepth(ke),ampd(ie,je,ke),phid(ie,je,ke))
OPEN (unit=11,file=ifile(1),form='unformatted',status='old')
OPEN (unit=12,file=ifile(2),form='unformatted',status='old')
OPEN (unit=13,file=ifile(3),form='unformatted',status='old')
OPEN (unit=14,file=ifile(4),form='unformatted',status='old')
OPEN (unit=21,file=ofile(1),form='unformatted',status='unknown')
OPEN (unit=22,file=ofile(2),form='unformatted',status='unknown')
do t=1,nt
do k=1,ke
read(11) id
idepth(k)=id(3)
read(11) ampx(:,:,k)
read(12) id
read(12) phix(:,:,k)
enddo
read(13) id
read(13) ampy
read(14)
read(14) phiy
ampd=dft
phid=dft
do i=1,ie
do j=1,je
do k=1,ke
if(ampx(i,j,k).ne.dft .and. phix(i,j,k).ne.dft .and. ampy(i,j).ne.dft .and. phiy(i,j).ne.dft) then
a=ampx(i,j,k)*cos(phix(i,j,k))-ampy(i,j)*cos(phiy(i,j))
b=ampx(i,j,k)*sin(phix(i,j,k))-ampy(i,j)*sin(phiy(i,j))
ampd(i,j,k)=sqrt(a*a+b*b)
w=atan2(b,a)
if(b.gt.0 .and. a.gt.0) then
phid(i,j,k)=w
endif
if(b.lt.0 .and. a.gt.0) then
phid(i,j,k)=w+2*pi
endif
if(b.gt.0 .and. a.lt.0) then
phid(i,j,k)=w
endif
if(b.lt.0 .and. a.lt.0) then
phid(i,j,k)=w+2*pi
endif
if(a.eq.0) then
if(b.gt.0) phid(i,j,k)=pi*0.5
if(b.lt.0) phid(i,j,k)=-pi*0.5
endif
if(b.eq.0) then
if(a.gt.0) phid(i,j,k)=0
if(a.lt.0) phid(i,j,k)=-pi
endif
endif
enddo ! k-loop
enddo ! j-loop
enddo ! i-loop
do k=1,ke
write(21) id(1),61,idepth(k),id(4)
write(21) ampd(:,:,k)
write(22) id(1),71,idepth(k),id(4)
write(22) phid(:,:,k)
enddo
enddo ! t-loop
end PROGRAM diff
..........................................................................................
script for calculating the internal tide bottom pressure
.......................................................................
#!/bin/ksh
# calculating bottom pressure
#
#
cd /scratch/m/m221050
d=/work/mh0256/m221050/zhuhua/paper
if1=${d}/amp_ppprim.nc
if2=${d}/phi_ppprim.nc
of1=${d}/amp_ppprim_bottom.nc
of2=${d}/phi_ppprim_bottom.nc
of3=${d}/ppprim_bottom-reconstructed_11-22.nc
g=/pool/data/MPIOM/TP6M/TP6Ms.nc
cdo -f ext copy ${if1} amp
cdo -f ext copy ${if2} phi
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=40
nt=1
ifile=amp
ofile=amp_bottom
/
EOF
/home/mpim/m221050/zhuhua/paper/bottom_p.x
cdo -f nc setgrid,${g} amp_bottom amp_bottom.nc
cp amp_bottom.nc ${of1}
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=40
nt=1
ifile=phi
ofile=phi_bottom
/
EOF
/home/mpim/m221050/zhuhua/paper/bottom_p.x
cdo -f nc setgrid,${g} phi_bottom phi_bottom.nc
cp phi_bottom.nc ${of2}
## reconstruction amplitude*cos(p0) with p0=w*t+phase, w=2*pi/12.42=0.50589 for M2
## at t=11,...,22
##
w=0.50589
time0=:00:00
for t in $(seq 11 22)
do
c=$(echo "$t * $w" | bc)
echo $w $t $c
cdo subc,$c phi_bottom.nc p0_$t
time=${t}${time0}
cdo settime,${time} -mul amp_bottom.nc -cos p0_$t p_$t
done
cdo -f nc mergetime p_?? ${of3}
#cdo copy p_11 ${of3}
exit
.................................
fortran program bottom_p.f90
................................
PROGRAM bottom_tho
! 10.5.2019
! calculating p in the bottom layer
!
! level k is identified as the bottom when p(k)=dft
!
!
IMPLICIT NONE
CHARACTER(180) :: ifile,ofile
INTEGER :: i,j,k,t,nc
INTEGER :: ie,je,ke,nt,id(4)
REAL, ALLOCATABLE :: p(:,:,:),pb(:,:)
real,parameter :: dft=-9.e+33
!
! read parameters, input and output, and allocate the variables
!
NAMELIST /CTL/ ie, je, ke, nt, ifile, ofile
OPEN(10,FILE='INPUT',STATUS='UNKNOWN', &
& ACCESS='SEQUENTIAL',FORM='FORMATTED')
READ(10,CTL)
close(10)
write(0,*) ie, je, ke, nt, ifile,ofile
ALLOCATE (p(ie,je,ke),pb(ie,je))
OPEN (unit=11,file=ifile,form='unformatted',status='old')
OPEN (unit=21,file=ofile,form='unformatted',status='unknown')
do t=1,nt
do k=1,ke
read(11) id
read(11) p(:,:,k)
enddo
pb=dft
do i=1,ie
do j=1,je
do k=1,ke
if(p(i,j,k).ne.dft) then
pb(i,j)=p(i,j,k)
endif
enddo ! k-loop
enddo ! j-loop
enddo ! i-loop
write(21) id(1),id(2),1,id(4)
write(21) pb
enddo
end PROGRAM bottom_tho
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.7:
..............................................................................................................
applying the following script to the global barotropic velocity and the internal tide bottom pressure
reconstructured using the respective amplitude and the phase
...............................................................................................................
#!/bin/ksh
# reconstruct the real component from amplitude and phase of a complex variable
#
cd /scratch/m/m221050
# reconstruct pporim_bottom
od=/work/mh0256/m221050/zhuhua/paper
ap=${od}/amp_ppprim_bottom.nc
pp=${od}/phi_ppprim_bottom.nc
ofp=${od}/ppprim_bottom_reconstructed_1-744.nc
## reconstruct a time series of a internal-tide variable from ampp (amplitude) and phip (phase in degrees)
## of an internal tide using
## p-prim=ampp*cos(p0), p0=w*t+phip, w=2*pi/12.42=0.50589 for M2 internal tide
##
w=0.50589
time0=:00:00
for t in $(seq 1 744)
do
c=$(echo "$t * $w" | bc)
echo $w $t $c
cdo subc,$c ${pp} pu0_$t
time=${t}${time0}
cdo settime,${time} -mul ${ap} -cos pu0_$t p_$t
done
cdo mergetime p_?? ${ofp}
# reconstruct barotropic U,V
id=/work/mh0256/m300212/ITgeneration_Oct2018/proc1_BartBarc_UV/data/proc1_BartBarc_AmpPhs
od=/work/mh0256/m221050/zhuhua/paper
u=${id}/Ubart_m2_2012Jan_stormtide-ncep.nc
v=${id}/Vbart_m2_2012Jan_stormtide-ncep.nc
ofu=${od}/Ubart_reconstructed_11-22_u2s.nc
ofv=${od}/Vbart_reconstructed_11-22_u2s.nc
g=/pool/data/MPIOM/TP6M/TP6Ms.nc
cdo -f ext -b f32 selvar,m2amp $u au
cdo -f ext -b f32 mulc,0.0174532 -selvar,m2phase $u phu # 0.0174532=pi/180
cdo -f ext -b f32 selvar,m2amp $v av
cdo -f ext -b f32 mulc,0.0174532 -selvar,m2phase $v phv # 0.0174532=pi/180
# 2. interpolate these amplitudes and phases of barotropic velocities (au, pu, av, pv) to scalar points
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=1
nt=1
ifile=au
ofile=ampu
/
EOF
/home/mpim/m221050/zhuhua/paper/u2s.x
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=1
nt=1
ifile=phu
ofile=phiu
/
EOF
/home/mpim/m221050/zhuhua/paper/u2s.x
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=1
nt=1
ifile=av
ofile=ampv
/
EOF
/home/mpim/m221050/zhuhua/paper/v2s.x
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=1
nt=1
ifile=phv
ofile=phiv
/
EOF
/home/mpim/m221050/zhuhua/paper/v2s.x
## reconstruct p-prim at t=11,...,22 from ampp (amplitude) and phip (phase in degrees) of an internal tide using
## p-prim=ampp*cos(p0), p0=w*t+phip, w=2*pi/12.42=0.50589 for M2 internal tide
##
cdo -f nc setgrid,${g} -setgrid,r3600x2394 ampu ampu.nc
cdo -f nc setgrid,${g} -setgrid,r3600x2394 ampv ampv.nc
cdo -f nc setgrid,${g} -setgrid,r3600x2394 phiu phiu.nc
cdo -f nc setgrid,${g} -setgrid,r3600x2394 phiv phiv.nc
w=0.50589
time0=:00:00
for t in $(seq 11 22)
do
c=$(echo "$t * $w" | bc)
echo $w $t $c
cdo subc,$c phiu.nc pu0_$t
time=${t}${time0}
cdo settime,${time} -mul ampu.nc -cos pu0_$t u_$t
cdo subc,$c phiv.nc pv0_$t
time=${t}${time0}
cdo settime,${time} -mul ampv.nc -cos pv0_$t v_$t
done
cdo mergetime u_?? ${ofu}
cdo mergetime v_?? ${ofv}
exit
exit
...................................................
Script for plot (NCL)
...................................................
begin
;;read data
itm = 9
in = addfile("data/ppprim_bottom-reconstructed_11-22_box_190_199_22_26.nc","r")
pprm = in->var61(itm-1,0,:,:)
pprm@lon2d=in->lon
pprm@lat2d=in->lat
printVarSummary(pprm)
Ufile = addfile("data/Ubart_reconstructed_11-22_u2s_box_190_199_22_26.nc","r")
Ubt = Ufile->var1(itm-1,::2,::2)
Ubt@lon2d = Ufile->lon(::2,::2)
Ubt@lat2d = Ufile->lat(::2,::2)
printVarSummary(Ubt)
Vfile = addfile("data/Vbart_reconstructed_11-22_u2s_box_190_199_22_26.nc","r")
Vbt = Vfile->var1(itm-1,::2,::2)
Vbt@lon2d = Vfile->lon(::2,::2)
Vbt@lat2d = Vfile->lat(::2,::2)
printVarSummary(Vbt)
Hfile = addfile("data/depto_box_190_199_22_26.nc","r")
depto = Hfile->depto(0,0,:,:)
depto@lon2d = Hfile->lon
depto@lat2d = Hfile->lat
printVarSummary(depto)
;;create plots
wks = gsn_open_wks("pdf","figures/UV-pprm-topo_Hawaii_tm-"+itm)
res = True
res@gsnDraw = False
res@gsnFrame = False
res@gsnLeftString = ""
res@gsnRightString = ""
;;set map
mpres = res
mpres@mpDataSetName = "Earth..4"
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpOutlineOn = True
mpres@mpGeophysicalLineThicknessF = 2
mpres@mpFillDrawOrder = "PostDraw"
;;set area
mpres@mpMinLatF = 22
mpres@mpMaxLatF = 26
mpres@mpMinLonF = 190
mpres@mpMaxLonF = 199
mpres@gsnMajorLatSpacing = 1
mpres@gsnMajorLonSpacing = 1
;;set contour
cnres = res
cnres@cnFillDrawOrder = "PreDraw"
cnres@cnFillOn = True
cnres@cnLinesOn = False
cnres@lbLabelBarOn = True
cnres@lbOrientation = "Vertical"
cnres@lbLabelFontHeightF = 0.008
cnres@lbLabelStride = 2
cnres@pmLabelBarOrthogonalPosF = 0.02
cnres@pmLabelBarWidthF = 0.06
cnres@cnLevelSelectionMode = "ManualLevels"
cnres@cnMinLevelValF = -250
cnres@cnMaxLevelValF = 250
cnres@cnLevelSpacingF = 25
cmap = read_colormap_file("cmp_b2r")
cnres@cnFillPalette = cmap
cnres@gsnLeftString = ""
;; set contour lines for depto
clres = res
clres@cnFillOn = False
clres@cnLinesOn = True
clres@cnLineThicknesses = 2.5
clres@cnLabelMasking = True
clres@cnLineLabelBackgroundColor = "transparent"
clres@cnLineLabelFontHeightF = 0.008
clres@cnInfoLabelOn = False
clres@cnLineColor = "gray10"
clres@cnLevelSelectionMode = "ManualLevels"
clres@cnMinLevelValF = 200
clres@cnMaxLevelValF = 6000
clres@cnLevelSpacingF = 500
;;set vector
res_vc = res
res_vc@vcGlyphStyle = "LineArrow"
res_vc@vcMinDistanceF = 0.01
res_vc@vcRefLengthF = 0.03
res_vc@vcMinFracLengthF = 0
res_vc@vcMinLevelValF = 0.05
res_vc@vcLineArrowThicknessF = 1.5
res_vc@vcFillArrowsOn = True
res_vc@vcMonoLineArrowColor = True
res_vc@vcMonoFillArrowFillColor = True
res_vc@vcMonoLineArrowColor = True
res_vc@vcLineArrowColor = "black"
res_vc@vcFillArrowFillColor = "black"
res_vc@vcFillArrowEdgeColor = "black"
res_vc@vcRefAnnoOn = True
res_vc@vcRefMagnitudeF = 0.2
res_vc@vcRefAnnoString1 = "0.2"
res_vc@vcRefAnnoSide = "Top"
res_vc@vcRefAnnoString2On = False
res_vc@vcRefAnnoPerimOn = False
res_vc@vcVectorDrawOrder = "PostDraw"
res_vc@gsnRightString = ""
;;plot
map = gsn_csm_map(wks,mpres)
contour = gsn_csm_contour(wks,pprm,cnres)
cnlines = gsn_csm_contour(wks,depto,clres)
vector = gsn_csm_vector(wks,Ubt,Vbt,res_vc)
overlay(map,contour)
overlay(map,cnlines)
overlay(map,vector)
draw(map)
frame(wks)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.8:
.................................................
Script for Fig.8a (NCL)
- Applying the following script to the phase of internal tide bottom pressure p''_b obtained from the script of Fig.6
.................................................
begin
;;read u,v,t from the data at 500hPa
in = addfile("data/phi_ppprim_bottom_mulc-57_box_190_199_22_26.nc","r")
pphs = in->var71(0,0,:,:)
pphs@lon2d=in->lon
pphs@lat2d=in->lat
printVarSummary(pphs)
Hfile = addfile("data/depto_box_190_199_22_26.nc","r")
depto = Hfile->depto(0,0,:,:)
depto@lon2d = Hfile->lon
depto@lat2d = Hfile->lat
printVarSummary(depto)
;;create plots
wks = gsn_open_wks("pdf","figures/phase_pprm_bottom_Hawaii")
res = True
res@gsnDraw = False
res@gsnFrame = False
res@gsnLeftString = ""
res@gsnRightString = ""
;;set map
mpres = res
mpres@mpDataSetName = "Earth..4"
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpOutlineOn = True
mpres@mpGeophysicalLineThicknessF = 2
mpres@mpFillDrawOrder = "PostDraw"
;;set area
mpres@mpMinLatF = 22
mpres@mpMaxLatF = 26
mpres@mpMinLonF = 190
mpres@mpMaxLonF = 199
mpres@gsnMajorLatSpacing = 1
mpres@gsnMajorLonSpacing = 1
;;set contour
cnres = res
cnres@cnFillDrawOrder = "PreDraw"
cnres@cnFillOn = True
cnres@cnLinesOn = False
cnres@lbLabelBarOn = True
cnres@lbOrientation = "Vertical"
cnres@lbLabelFontHeightF = 0.008
cnres@lbLabelStride = 2
cnres@pmLabelBarOrthogonalPosF = 0.02
cnres@pmLabelBarWidthF = 0.06
cnres@cnLevelSelectionMode = "ManualLevels"
cnres@cnMinLevelValF = 0
cnres@cnMaxLevelValF = 360
cnres@cnLevelSpacingF = 20
cnres@cnLabelBarEndStyle = "ExcludeOuterBoxes"
cmap = read_colormap_file("NCV_bright")
cnres@cnFillPalette = cmap
cnres@gsnLeftString = ""
;; set contour lines for depto
clres = res
clres@cnFillOn = False
clres@cnLinesOn = True
clres@cnLineThicknesses = 2.5
clres@cnLabelMasking = True
clres@cnLineLabelBackgroundColor = "transparent"
clres@cnLineLabelFontHeightF = 0.008
clres@cnInfoLabelOn = False
clres@cnLineColor = "gray10"
clres@cnLevelSelectionMode = "ManualLevels"
clres@cnMinLevelValF = 200
clres@cnMaxLevelValF = 6000
clres@cnLevelSpacingF = 500
;;plot
map = gsn_csm_map(wks,mpres)
contour = gsn_csm_contour(wks,pphs,cnres)
cnlines = gsn_csm_contour(wks,depto,clres)
overlay(map,contour)
overlay(map,cnlines)
draw(map)
frame(wks)
end
.................................................
Script for Fig.8b and Fig.8c (Matlab)
.................................................
clear all; clc;
slon = 192;
slat = 22;
elon = 193.7;
elat = 26;
disp('import pbot')
fname = 'data/phi_ppprim_bottom.nc';
lon2d = ncread(fname,'lon');
lat2d = ncread(fname,'lat');
temphs = ncread(fname,'var71');
phs = radtodeg(temphs);
clear fname temphs;
disp('import depto')
fname = 'data/depto_box_190_199_22_26.nc';
hlon = ncread(fname,'lon');
hlon(hlon<0) = hlon(hlon<0)+360;
hlat = ncread(fname,'lat');
depto = ncread(fname,'depto');
H1d = reshape(depto,[],1);
len = length(find(abs(H1d)>10^9))
clear fname;
disp('get points from sectin')
plon = slon+0.1 : 0.1 : elon;-0.1
plat = slat + (elat-slat)/(elon-slon) .* (plon-slon);
disp('interp')
psec = griddata(double(lon2d),double(lat2d),double(phs),double(plon),double(plat));
hsec = griddata(double(hlon),double(hlat),double(depto),double(plon),double(plat));
disp('output data')
data = [reshape(plon,[],1) reshape(plat,[],1) reshape(psec,[],1) reshape(hsec,[],1)]
dlmwrite('interpoated_phase_depth_section.dat',data,' ');
clf;
disp('make plot')
subplot(211);
plot(plat,psec,'-k','linewidth',1.5);
grid on;
ylabel('phase (degree)','fontsize',15);
pbaspect([4 1 1]);
subplot(212);
plot(plat,-hsec,'-k','linewidth',1.5);
ylabel('depth (m)','fontsize',15);
xlabel('latitude (degree)','fontsize',15);
axis([22 26 -6000 0]);
grid on;
pbaspect([4 1 1]);
saveas(gcf,'phase_section.pdf');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.9:
..................................................................
script: calculating cosine of the phase difference
...............
#!/bin/ksh
#
# calculate the difference in phase
#
cd /scratch/m/m221050
od=/work/mh0256/m221050/zhuhua/paper
phip=phi_ppprim_bottom
phiu=phiu
phiv=phiv
cdo sub -mulc,57.296 ${od}/${phip}.nc -mulc,57.296 ${od}/${phiu}.nc ${od}/sub_${phip}_${phiu}.nc
cdo sub -mulc,57.296 ${od}/${phip}.nc -mulc,57.296 ${od}/${phiv}.nc ${od}/sub_${phip}_${phiv}.nc
exit
................................................
Script for plot (NCL)
......................................
begin
;;read u,v,t from the data at 500hPa
in = addfile("data/cosphix_box_190_199_22_26.nc","r")
cosphi = in->var71(0,0,:,:)
cosphi@lon2d=in->lon
cosphi@lat2d=in->lat
printVarSummary(cosphi)
Hfile = addfile("data/depto_box_190_199_22_26.nc","r")
depto = Hfile->depto(0,0,:,:)
depto@lon2d = Hfile->lon
depto@lat2d = Hfile->lat
printVarSummary(depto)
;;create plots
wks = gsn_open_wks("pdf","figures/cos-phidiff_pphi-uphi_Hawaii")
res = True
res@gsnDraw = False
res@gsnFrame = False
res@gsnLeftString = ""
res@gsnRightString = ""
;;set map
mpres = res
mpres@mpDataSetName = "Earth..4"
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpOutlineOn = True
mpres@mpGeophysicalLineThicknessF = 2
mpres@mpFillDrawOrder = "PostDraw"
;;set area
mpres@mpMinLatF = 22
mpres@mpMaxLatF = 26
mpres@mpMinLonF = 190
mpres@mpMaxLonF = 199
mpres@gsnMajorLatSpacing = 1
mpres@gsnMajorLonSpacing = 1
;;set contour
cnres = res
cnres@cnFillDrawOrder = "PreDraw"
cnres@cnFillOn = True
cnres@cnLinesOn = False
cnres@lbLabelBarOn = True
cnres@lbOrientation = "Vertical"
cnres@lbLabelFontHeightF = 0.008
cnres@lbLabelStride = 2
cnres@pmLabelBarOrthogonalPosF = 0.02
cnres@pmLabelBarWidthF = 0.06
cnres@cnLevelSelectionMode = "ManualLevels"
cnres@cnMinLevelValF = -1
cnres@cnMaxLevelValF = 1
cnres@cnLevelSpacingF = 0.1
cnres@cnLabelBarEndStyle = "ExcludeOuterBoxes"
cmap = read_colormap_file("MPL_coolwarm")
cnres@cnFillPalette = cmap
cnres@gsnLeftString = ""
;; set contour lines for depto
clres = res
clres@cnFillOn = False
clres@cnLinesOn = True
clres@cnLineThicknesses = 2.5
clres@cnLabelMasking = True
clres@cnLineLabelBackgroundColor = "transparent"
clres@cnLineLabelFontHeightF = 0.008
clres@cnInfoLabelOn = False
clres@cnLineColor = "gray10"
clres@cnLevelSelectionMode = "ManualLevels"
clres@cnMinLevelValF = 200
clres@cnMaxLevelValF = 6000
clres@cnLevelSpacingF = 500
;;plot
map = gsn_csm_map(wks,mpres)
contour = gsn_csm_contour(wks,cosphi,cnres)
cnlines = gsn_csm_contour(wks,depto,clres)
overlay(map,contour)
overlay(map,cnlines)
draw(map)
frame(wks)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.10:
..........................................................................................
script: calculating the work done by the internal-tide form drag at the seafloor
.........
#!/bin/ksh
#
# calculate the conversion rate P from p'' at the bottom and barotropic velocities
#
cd /scratch/m/m221050
id=/work/mh0256/m300212/ITgeneration_Oct2018/proc1_BartBarc_UV/data/proc1_BartBarc_AmpPhs
od=/work/mh0256/m221050/zhuhua/paper
ampp=${od}/amp_ppprim_bottom.nc
phip=${od}/phi_ppprim_bottom.nc
u=${id}/Ubart_m2_2012Jan_stormtide-ncep.nc
v=${id}/Vbart_m2_2012Jan_stormtide-ncep.nc
of3a=${od}/Px+Py.nc
of3b=${od}/Px+Py_minus.nc
of3c=${od}/Px+Py_minus_mulc-0.5.nc
odx=${od}/dddx.nc
ody=${od}/dddy.nc
g=/pool/data/MPIOM/TP6M/TP6Ms.nc
# 1. get amplitude and phase of barotropic velocity
cdo -f ext -b f32 selvar,m2amp $u au
cdo -f ext -b f32 mulc,0.0174532 -selvar,m2phase $u phu # 0.0174532=pi/180
cdo -f ext -b f32 selvar,m2amp $v av
cdo -f ext -b f32 mulc,0.0174532 -selvar,m2phase $v phv # 0.0174532=pi/180
# 2. interpolate these amplitudes and phases of barotropic velocities (au, pu, av, pv) to scalar points
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=1
nt=1
ifile=au
ofile=ampu
/
EOF
/home/mpim/m221050/zhuhua/paper/u2s.x
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=1
nt=1
ifile=phu
ofile=phiu
/
EOF
/home/mpim/m221050/zhuhua/paper/u2s.x
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=1
nt=1
ifile=av
ofile=ampv
/
EOF
/home/mpim/m221050/zhuhua/paper/v2s.x
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ke=1
nt=1
ifile=phv
ofile=phiv
/
EOF
/home/mpim/m221050/zhuhua/paper/v2s.x
cdo -f nc setgrid,${g} -setgrid,r3600x2394 ampu ampu.nc
cdo -f nc setgrid,${g} -setgrid,r3600x2394 ampv ampv.nc
cdo -f nc setgrid,${g} -setgrid,r3600x2394 phiu phiu.nc
cdo -f nc setgrid,${g} -setgrid,r3600x2394 phiv phiv.nc
cp ampu.nc ${od}/.
cp ampv.nc ${od}/.
cp phiu.nc ${od}/.
cp phiv.nc ${od}/.
# 3. calculating the conversion rate
cdo -f ext cos -sub ${phip} phiu.nc cosphix
cdo -f ext cos -sub ${phip} phiv.nc cosphiy
cdo -f ext mul ${ampp} ampu.nc pu
cdo -f ext mul ${ampp} ampv.nc pv
cdo -f ext copy ${odx} dddx
cdo -f ext copy ${ody} dddy
# get joined mask
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ifile=pu,dddx,cosphix
ofile=pun,dddxn,cosphixn
/
EOF
/home/mpim/m221050/zhuhua/paper/mask3.x
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
ifile=pv,dddy,cosphiy
ofile=pvn,dddyn,cosphiyn
/
EOF
/home/mpim/m221050/zhuhua/paper/mask3.x
cdo mul pun cosphixn PU
cdo mul pvn cosphiyn PV
cdo mul PU dddxn Px
cdo mul PV dddyn Py
cdo add Px Py PxPy
cdo -f nc setgrid,${g} -setgrid,r3600x2394 PxPy PxPy.nc
cp PxPy.nc ${of3a}
cdo mulc,-1 ${of3a} ${of3b}
cdo mulc,0.5 ${of3b} ${of3c}
exit
.......................................
fortran programs:
- u2s.f90 and v2s.f90 are given above
- mask3.f90 is given below
.......................................
PROGRAM mask
! 19.5.2019
! remask the 3 input variables using joined mask
!
!
IMPLICIT NONE
CHARACTER(180) :: ifile(3),ofile(3)
INTEGER :: i,j,nc
INTEGER :: ie,je,id(4)
REAL, ALLOCATABLE :: x(:,:),y(:,:),z(:,:),xn(:,:),yn(:,:),zn(:,:)
real,parameter :: dft=-9.e+33,pi=3.1416
!
! read parameters, input and output, and allocate the variables
!
NAMELIST /CTL/ ie, je, ifile, ofile
OPEN(10,FILE='INPUT',STATUS='UNKNOWN', &
& ACCESS='SEQUENTIAL',FORM='FORMATTED')
READ(10,CTL)
close(10)
write(0,*) ie, je, ifile,ofile
ALLOCATE (x(ie,je),y(ie,je),z(ie,je),xn(ie,je),yn(ie,je),zn(ie,je))
OPEN (unit=11,file=ifile(1),form='unformatted',status='old')
OPEN (unit=12,file=ifile(2),form='unformatted',status='old')
OPEN (unit=13,file=ifile(3),form='unformatted',status='old')
OPEN (unit=21,file=ofile(1),form='unformatted',status='unknown')
OPEN (unit=22,file=ofile(2),form='unformatted',status='unknown')
OPEN (unit=23,file=ofile(3),form='unformatted',status='unknown')
read(11) id
read(11) x
print *,' read in x:', id
read(12) id
read(12) y
print *,' read in y:', id
read(13) id
read(13) z
print *,' read in z:', id
xn=dft
yn=dft
zn=dft
do i=1,ie
do j=1,je
if(x(i,j).ne.dft .and. y(i,j).ne.dft .and. z(i,j).ne.dft) then
xn(i,j)=x(i,j)
yn(i,j)=y(i,j)
zn(i,j)=z(i,j)
endif
enddo ! j-loop
enddo ! i-loop
write(21) id
write(21) xn
write(22) id
write(22) yn
write(23) id
write(23) zn
end PROGRAM mask
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.11:
.......................................
Script for plot (NCL)
- Applying the following script to the global field obtained from the script for Fig.10
.......................................
begin
;;read u,v,t from the data at 500hPa
in = addfile("data/Px+Py_minus_box_190_199_22_26.nc","r")
tem = in->var71(0,0,:,:)
pgen = tem/2
pgen@lon2d=in->lon
pgen@lat2d=in->lat
printVarSummary(pgen)
Hfile = addfile("data/depto_box_190_199_22_26.nc","r")
depto = Hfile->depto(0,0,:,:)
depto@lon2d = Hfile->lon
depto@lat2d = Hfile->lat
printVarSummary(depto)
;;create plots
wks = gsn_open_wks("pdf","figures/M2-ITgen_Hawaii_new")
res = True
res@gsnDraw = False
res@gsnFrame = False
res@gsnLeftString = ""
res@gsnRightString = ""
;;set map
mpres = res
mpres@mpDataSetName = "Earth..4"
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpOutlineOn = True
mpres@mpGeophysicalLineThicknessF = 2
mpres@mpFillDrawOrder = "PostDraw"
;;set area
mpres@mpMinLatF = 22
mpres@mpMaxLatF = 26
mpres@mpMinLonF = 190
mpres@mpMaxLonF = 199
mpres@gsnMajorLatSpacing = 1
mpres@gsnMajorLonSpacing = 1
;;set contour
cnres = res
cnres@cnFillDrawOrder = "PreDraw"
cnres@cnFillOn = True
cnres@cnLinesOn = False
cnres@lbLabelBarOn = True
cnres@lbOrientation = "Vertical"
cnres@lbLabelFontHeightF = 0.008
cnres@lbLabelStride = 2
cnres@pmLabelBarOrthogonalPosF = 0.02
cnres@pmLabelBarWidthF = 0.06
cnres@cnLevelSelectionMode = "ManualLevels"
cnres@cnMinLevelValF = -0.16
cnres@cnMaxLevelValF = 0.16
cnres@cnLevelSpacingF = 0.01
cnres@cnLabelBarEndStyle = "ExcludeOuterBoxes"
cmap = read_colormap_file("MPL_RdBu")
cnres@cnFillPalette = cmap(::-1,:)
cnres@gsnLeftString = ""
;; set contour lines for depto
clres = res
clres@cnFillOn = False
clres@cnLinesOn = True
clres@cnLineThicknesses = 2.5
clres@cnLabelMasking = True
clres@cnLineLabelBackgroundColor = "transparent"
clres@cnLineLabelFontHeightF = 0.008
clres@cnInfoLabelOn = False
clres@cnLineColor = "gray10"
clres@cnLevelSelectionMode = "ManualLevels"
clres@cnMinLevelValF = 200
clres@cnMaxLevelValF = 6000
clres@cnLevelSpacingF = 500
;;plot
map = gsn_csm_map(wks,mpres)
contour = gsn_csm_contour(wks,pgen,cnres)
cnlines = gsn_csm_contour(wks,depto,clres)
overlay(map,contour)
overlay(map,cnlines)
draw(map)
frame(wks)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.12:
.......................................................
script 1:
sort the bottom generation through the generation depth
........................................................
#!/bin/ksh
# calculating area-integral for selected basin
#
# IndoPac: 7+8
# Atlantic: 4+5
# SO: 6
cd /scratch/m/m221050
od=/work/mh0256/m221050/zhuhua/paper
if1=${od}/P-reconstructed_minus.nc
d0=0.1
if1=${od}/Px+Py_minus_mulc-0.5.nc
#if1=${od}/Px+Py_using-depthavg-p_minus_mulc-0.5.nc
ibek1=4
ibek2=5
of1=${od}/Px+Py_minus_mulc-0.5_bek${ibek1}${ibek2}_z_pro-unit-area.ext
of2=${od}/Px+Py_minus_mulc-0.5_bek${ibek1}${ibek2}_d.ext
#of1=${od}/Px+Py_minus_deeper-than-1000.nc
g=/pool/data/MPIOM/TP6M/TP6Ms.nc
#cdo selvar,dlxp /pool/data/MPIOM/TP6M/TP6ML40_fx.nc dlxp.nc
#cdo selvar,dlyp /pool/data/MPIOM/TP6M/TP6ML40_fx.nc dlyp.nc
#cdo -f ext copy -selindexbox,2,3601,1,2394 -setgrid,r3602x2394 dlxp.nc dlxp
#cdo -f ext copy -selindexbox,2,3601,1,2394 -setgrid,r3602x2394 dlyp.nc dlyp
#cdo -f ext copy -selindexbox,2,3601,1,2394 -setgrid,r3602x2394 /work/mh0256/m221050/zhuhua/paper/depto.nc depto
#cdo -f ext copy -selindexbox,2,3601,1,2394 -setgrid,r3602x2394 /pool/data/MPIOM/TP6M/RIBEK.nc ribek
cdo -f ext copy ${if1} x
cat >INPUT<<EOF
&CTL
ie=3600
je=2394
nt=1
ibek=${ibek1},${ibek2}
lavg=true
ifile=dlxp,dlyp,depto,ribek,x
ofile=z,d
/
EOF
/home/mpim/m221050/zhuhua/paper/area-int_basin.x
cp z ${of1}
cp d ${of2}
exit
...............................................................
fortran program 1:
- for the global basin
......................
PROGRAM area_avg
! 13-5-19
! for a given variable x, calculating the average / integral
!
! input (fort.11): dlxp, dlyp
! input (fort.12): x (quantity to be averaged / integrated)
!
! output: z: generation P*dx*dy
! d: depth of x
! z and d will be used to calculate the generation as a function of depth
IMPLICIT NONE
CHARACTER(180) :: ifile(4),ofile(2)
INTEGER :: i,j,id(4),t,n
INTEGER :: ie,je,nt
REAL, ALLOCATABLE :: x(:,:),z(:),d(:)
real, allocatable :: dlxp(:,:),dlyp(:,:),depto(:,:)
real,parameter :: dft=-9.e+33
real :: area,xs,xp,xn,p0
logical :: lavg
!
! read parameters, input and output, and allocate the variables
!
NAMELIST /CTL/ ie, je, nt, lavg, ifile, ofile
OPEN(10,FILE='INPUT',STATUS='UNKNOWN', &
& ACCESS='SEQUENTIAL',FORM='FORMATTED')
READ(10,CTL)
write(0,*) ie, je, nt, lavg, ifile, ofile
close(10)
ALLOCATE (x(ie,je),z(ie*je),d(ie*je))
ALLOCATE (dlxp(ie,je),dlyp(ie,je),depto(ie,je))
OPEN (unit=11,file=ifile(1),form='unformatted',status='old')
OPEN (unit=12,file=ifile(2),form='unformatted',status='old')
OPEN (unit=13,file=ifile(3),form='unformatted',status='old')
OPEN (unit=14,file=ifile(4),form='unformatted',status='old')
OPEN (unit=21,file=ofile(1),form='unformatted',status='unknown')
OPEN (unit=22,file=ofile(2),form='unformatted',status='unknown')
!
! read grid distance
!
read(11) id
read(11) dlxp
print *, ' read grid x-distance:', id
read(12) id
read(12) dlyp
print *, ' read grid y-distance:', id
read(13) id
read(13) depto
!
! time loop
!
do t=1,nt
read(14) id
read(14) x
!
! calculating the area
!
if(t.eq.1) then
area=0.
do i=1,ie
do j=1,je
if(x(i,j).ne.dft .and. dlxp(i,j).ne.dft .and. dlyp(i,j).ne.dft) then
area=area+dlxp(i,j)*dlyp(i,j)
endif
enddo
enddo
print *, ' area [m*m]:',area
endif
xs=0
xp=0
xn=0
n=0
do i=1,ie
do j=1,je
if(x(i,j).ne.dft .and. dlxp(i,j).ne.dft .and. dlyp(i,j).ne.dft) then
n=n+1
z(n)=x(i,j)*dlxp(i,j)*dlyp(i,j)
d(n)=depto(i,j)
xs=xs+z(n)
if(x(i,j).gt.0) then
xp=xp+1
else
xn=xn+1
endif
! print '(2i6,e12.3)',i,j,x(i,j)
endif
enddo
enddo
if(lavg) then
xs=xs/area
do i=1,n
z(i)=z(i)/area
enddo
endif
print *, xp,xn
enddo ! time loop
print *, ' area, area avg /int:', area, xs
write(21) id(1),id(2),id(3),n
write(21) (z(i),i=1,n)
write(22) id(1),id(2),id(3),n
write(22) (d(i),i=1,n)
END PROGRAM area_avg
...........................................................
fortran program 2:
- for the the Atlantic and Indo-Pacific basin
......................
PROGRAM area_avg
! 19.11.19
! same as area-int.f90, but for different basins using /pool/data/MPIOM/TP6M/RIBEK.nc
! Indo-Pacific: 7+8
! Atlantic: 4+5
! Southern Ocean: 6
! Note: two basins (determined by ibek(2)) can be consiered together.
! If only one basin is considered, ibek(2)=0
!
! 13-5-19
! for a given variable x, calculating the average / integral
!
! input (fort.11): dlxp, dlyp
! input (fort.12): x (quantity to be averaged / integrated)
!
! output: z: generation P*dx*dy
! d: depth of x
! z and d will be used to calculate the generation as a function of depth
IMPLICIT NONE
CHARACTER(180) :: ifile(5),ofile(2)
INTEGER :: i,j,id(4),t,n,ibek(2)
INTEGER :: ie,je,nt
REAL, ALLOCATABLE :: x(:,:),z(:),d(:)
real, allocatable :: dlxp(:,:),dlyp(:,:),depto(:,:),bek(:,:)
real,parameter :: dft=-9.e+33
real :: area,xs,p0
logical :: lavg
!
! read parameters, input and output, and allocate the variables
!
NAMELIST /CTL/ ie, je, nt, ibek, lavg, ifile, ofile
OPEN(10,FILE='INPUT',STATUS='UNKNOWN', &
& ACCESS='SEQUENTIAL',FORM='FORMATTED')
READ(10,CTL)
write(0,*) ie, je, nt, ibek, lavg, ifile, ofile
close(10)
ALLOCATE (x(ie,je),z(ie*je),d(ie*je))
ALLOCATE (dlxp(ie,je),dlyp(ie,je),depto(ie,je),bek(ie,je))
OPEN (unit=11,file=ifile(1),form='unformatted',status='old')
OPEN (unit=12,file=ifile(2),form='unformatted',status='old')
OPEN (unit=13,file=ifile(3),form='unformatted',status='old')
OPEN (unit=14,file=ifile(4),form='unformatted',status='old')
OPEN (unit=15,file=ifile(5),form='unformatted',status='old')
OPEN (unit=21,file=ofile(1),form='unformatted',status='unknown')
OPEN (unit=22,file=ofile(2),form='unformatted',status='unknown')
!
! read grid distance
!
read(11) id
read(11) dlxp
print *, ' read grid x-distance:', id
read(12) id
read(12) dlyp
print *, ' read grid y-distance:', id
read(13) id
read(13) depto
read(14) id
read(14) bek
!
! time loop
!
do t=1,nt
read(15) id
read(15) x
!
! calculating the area
!
if(t.eq.1) then
area=0.
if(ibek(2).eq.0) then
do i=1,ie
do j=1,je
if(x(i,j).ne.dft .and. dlxp(i,j).ne.dft .and. dlyp(i,j).ne.dft) then
if(bek(i,j).eq.ibek(1)) then
area=area+dlxp(i,j)*dlyp(i,j)
endif
endif
enddo
enddo
else
do i=1,ie
do j=1,je
if(x(i,j).ne.dft .and. dlxp(i,j).ne.dft .and. dlyp(i,j).ne.dft) then
if(bek(i,j).eq.ibek(1) .or. bek(i,j).eq.ibek(2)) then
area=area+dlxp(i,j)*dlyp(i,j)
endif
endif
enddo
enddo
endif
print *, ' area [m*m]:',area
endif
xs=0
n=0
if(ibek(2).eq.0) then
do i=1,ie
do j=1,je
if(x(i,j).ne.dft .and. dlxp(i,j).ne.dft .and. dlyp(i,j).ne.dft) then
if(bek(i,j).eq.ibek(1)) then
n=n+1
z(n)=x(i,j)*dlxp(i,j)*dlyp(i,j)
d(n)=depto(i,j)
xs=xs+z(n)
endif
endif
enddo
enddo
else
do i=1,ie
do j=1,je
if(x(i,j).ne.dft .and. dlxp(i,j).ne.dft .and. dlyp(i,j).ne.dft) then
if(bek(i,j).eq.ibek(1) .or. bek(i,j).eq.ibek(2)) then
n=n+1
z(n)=x(i,j)*dlxp(i,j)*dlyp(i,j)
d(n)=depto(i,j)
xs=xs+z(n)
endif
endif
enddo
enddo
endif
if(lavg) then
xs=xs/area
do i=1,n
z(i)=z(i)/area
enddo
endif
enddo ! time loop
print *, ' area, area avg /int:', area, xs
write(21) id(1),id(2),id(3),n
write(21) (z(i),i=1,n)
write(22) id(1),id(2),id(3),n
write(22) (d(i),i=1,n)
END PROGRAM area_avg
........................................................................
script 2:
integrating according to the generation depth
..........................................
#!/bin/ksh
# calculating the generation as a function of depth
#
cd /scratch/m/m221050
od=/work/mh0256/m221050/zhuhua/paper
nk=60
dd=100
if1=${od}/Px+Py_minus_mulc-0.5_bek78_z_pro-unit-area.ext
if2=${od}/Px+Py_minus_mulc-0.5_bek78_d.ext
of1=${od}/Px+Py_minus_mulc-0.5_bek78_fct-of-d_dd${dd}_pro-unit-area.ext
of2=${od}/Px+Py_minus_mulc-0.5_bek78_depth_dd${dd}_pro-unit-area.ext
# global dim: 5501287
# pacific: 2394871
# Atlantic: 833722
# SO: 1930622
rm -f z
rm -f d
cp ${if1} z
cp ${if2} d
cat >INPUT<<EOF
&CTL
idim=2394871
nk=${nk}
dd=${dd}
ifile=z,d
ofile=x,dp
/
EOF
/home/mpim/m221050/zhuhua/paper/fct-of-d.x
cp x ${of1}
cp dp ${of2}
exit
......................................................................................
fortran program needed: function-of-d.f90
.......................
PROGRAM function_of_d
! 4-6-19
! Given the internal-tide generation z at respective depth z, obtained as output z,d from area_int.f90,
! calculate the total generation as a function of depth, which equals the sum of all z at a given depth range dp(i)
!
! The total depth range from 22 m to 6020 m is devided into 120 intervals dp(i),i=1,...,120 with dp(i+1)-dp(i)=dd
!
! input (fort.11): z
! input (fort.12): d
!
! output: z sumed up at depth dp(i)
! dp: depths considered
!
IMPLICIT NONE
CHARACTER(180) :: ifile(2),ofile(2)
INTEGER :: i,nk,k,id(4)
INTEGER :: idim
REAL, ALLOCATABLE :: z(:),d(:),x(:),dp(:)
real,parameter :: dft=-9.e+33
real :: dd
!
! read parameters, input and output, and allocate the variables
!
NAMELIST /CTL/ idim, nk, dd, ifile, ofile
OPEN(10,FILE='INPUT',STATUS='UNKNOWN', &
& ACCESS='SEQUENTIAL',FORM='FORMATTED')
READ(10,CTL)
write(0,*) idim, nk, dd, ifile, ofile
close(10)
ALLOCATE (z(idim),d(idim),x(nk),dp(nk+1))
OPEN (unit=11,file=ifile(1),form='unformatted',status='old')
OPEN (unit=12,file=ifile(2),form='unformatted',status='old')
OPEN (unit=21,file=ofile(1),form='unformatted',status='unknown')
OPEN (unit=22,file=ofile(2),form='unformatted',status='unknown')
read(11) id
read(11) z
read(12) id
read(12) d
dp(1)=0
do k=1,nk
dp(k+1)=dp(k)+dd
enddo
print *, '########## dp #################'
print '(10f8.0)', dp
x=0
do i=1,idim
do k=1,nk
if(d(i).ge.dp(k) .and. d(i).lt.dp(k+1)) then
x(k)=x(k)+z(i)
goto 1
endif
enddo
1 continue
enddo
print *, '########## x #################'
print '(10e12.2)', x
write(21) id(1),id(2),id(2),nk
write(21) x
!write(21) x*1.e-9
write(22) id(1),id(2),id(2),nk
write(22) (dp(i),i=1,nk)
END PROGRAM function_of_d
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.A1 (NCL):
................................................
obtain the MOC data from monthly-mean data
Script 1:
- go_plot.sh
................................................
#!/bin/bash -ex
YEAR_START=2003
YEAR_END=2012
INTERVAL=10
export EXPID=hel17153-DZA_TP6ML40_mpiom-trunk
#export DATADIR=/work/mh0256/m300212/STORMTIDE_run/Tides_NCEP_SAL2/experiments/hel17153-DZA_TP6ML40_mpiom-trunk_12yr/outdata/
export DATADIR=/work/mh0256/m300212/STORMTIDE_run/Tides_NCEP_SAL2/Data_analysis/proc5_MOC/data/proc1_moc_yearMean
export CDO='cdo -s -P 4'
export CHUNK=1
export GRID=TP6M
export LEV=L40
export POOL=/pool/data/MPIOM/input/r0008
OQSDIR=$(pwd) #add here path to the OQs scripts
# Create plots
YEAR=$YEAR_START
while [ $YEAR -le $(( YEAR_END - INTERVAL + 1)) ]
do
Y1=$YEAR Y2=$(( YEAR + INTERVAL - 1 )) ${OQSDIR}/plot_mpiom.sh
YEAR=$(( YEAR + INTERVAL ))
done
exit
.................
Script 2:
- plot_mpiom.sh
.................
#! /usr/bin/env bash
set -e
set -x # DEBUG
export DATADIR=$DATADIR
export WORKDIR=$WORKDIR
export EXPID=$EXPID
export CDO=${CDO:-cdo -s -P 4}
export Y1=$Y1
export Y2=$Y2
export CHUNK=${CHUNK:-1}
export GRID=$GRID
export LEV=$LEV
export POOL=$POOL
export BINDIR=$(dirname $0)
run_bg () {(
trap 'status=$?; [ $status != 0 ] && echo $1 $status >> status' EXIT
time "$@"
) &}
rm -f status
#run_bg $BINDIR/OQs_mpiom_surf.sh
#run_bg $BINDIR/OQs_mpiom_subsurf.sh
run_bg $BINDIR/OQs_mpiom_moc.sh
wait
if [ -e status ]
then
echo "Sorry: errors during plotting: $(<status)" >&2
exit 1
fi
.................
script 3:
- OQs_mpiom_moc.sh
.................
#! /usr/bin/env bash
set -e
set -x # DEBUG
PROGRAM=$(basename $0)
TMP_DIR=${PROGRAM}_$Y1-$Y2
[ -e $TMP_DIR ] && rm -rf $TMP_DIR
mkdir -p $TMP_DIR
cd $TMP_DIR
CATFILE=moc.nc
CODE='100,101,102'
for YEAR in $(seq $Y1 $CHUNK $Y2)
do
INFILE=$(printf $DATADIR/${EXPID}_mpiom_data_moc_mm_%04.0f-01-01_%04.0f-12-31.nc $YEAR $((YEAR + CHUNK - 1)))
if [ -e $INFILE ]
then
$CDO cat -timmean -selcode,$CODE $INFILE $CATFILE
else
echo Warning ! $INFILE is missing, skipping to next file ...
fi
done
$CDO sinfo $CATFILE # DEBUG
inbase=$(printf ${EXPID}_mpiom_moc_%04.0f-01-01_%04.0f-12-31 $Y1 $Y2)
infile=$inbase.nc
outfile=$inbase.pdf
$CDO timmean $CATFILE $infile
$CDO -f nc -setunit,"Sv" -mulc,1.025e-9 -selcode,100 $infile glo_moc.$infile
ncrename -d depth_2,depth -v depth_2,depth glo_moc.$infile
$BINDIR/plotsec_ncl --rstring=glo --code=100 --min=-24 --max=24 --inc=2 --pal=BlueWhiteOrangeRed --title="" glo_moc.$infile
$CDO -f nc -setunit,"Sv" -mulc,1.025e-9 -selcode,101 $infile atl_moc.$infile
ncrename -d depth_2,depth -v depth_2,depth atl_moc.$infile
$BINDIR/plotsec_ncl --rstring=atl --code=101 --min=-24 --max=24 --inc=2 --pal=BlueWhiteOrangeRed --title="" atl_moc.$infile
$CDO -f nc -setunit,"Sv" -mulc,1.025e-9 -selcode,102 $infile indopacific_moc.$infile
ncrename -d depth_2,depth -v depth_2,depth indopacific_moc.$infile
$BINDIR/plotsec_ncl --rstring=indopac --code=102 --min=-24 --max=24 --inc=2 --pal=BlueWhiteOrangeRed --title="" indopacific_moc.$infile
cd -
mv $TMP_DIR/*.$outfile .
%rm -rf $TMP_DIR
................................................
Script for Fig. A1(a)
................................................
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
; =======;
tem=addfile("data/proc1_mass_sfbarotr/time-mean_mass_streamfunc_barotr_2003-2012.nc","r")
sf=tem->msftbarot
sf=sf/(1025*10^6)
sf@lon2d=tem->lon
sf@lat2d=tem->lat
sf@_FillValue = 9.97e+36
sf=where(sf.lt.-10^9,sf@_FillValue,sf)
printVarSummary(sf)
; =======
wks = gsn_open_wks("pdf","figures/proc2_streamfunc/barotr_streamfunc_volume_Sv_10yrMean_2003-2012_ncl")
gsn_define_colormap(wks,"cmp_b2r")
; =======
resa = True
resa@cnFillOn = True
resa@cnMissingValFillPattern = 0
resa@cnMissingValFillColor = "grey" ;;Color of miss values
resa@vpWidthF = 0.5
resa@vpHeightF = 0.3
resa@gsnSpreadColors = True
resa@lbLabelBarOn = True
resa@lbOrientation = "Vertical"
resa@lbLabelFontHeightF = 0.012
resa@lbLabelStride = 2
resa@lbLabelAutoStride = True
resa@pmLabelBarWidthF = 0.06
resa@pmLabelBarOrthogonalPosF = 0.02
resa@cnLevelSelectionMode = "ManualLevels"
resa@cnMinLevelValF = -300.
resa@cnMaxLevelValF = 300.
resa@cnLevelSpacingF = 20.
resa@gsnAddCyclic = False
resa@cnSmoothingOn = True
resa@mpMaxLatF = 90
resa@mpMaxLonF = 360
resa@mpMinLatF = -90
resa@mpMinLonF = 0
resa@mpCenterLonF = 180
resa@cnLinesOn = True
resa@cnLineLabelsOn = True
resa@cnLabelMasking = True
resa@cnLineLabelBackgroundColor = "transparent"
resa@cnLineLabelPlacementMode = "constant"
resa@cnLineLabelInterval = 2
resa@cnLineColor = "black"
resa@cnLineLabelFontHeightF = 0.006
resa@cnInfoLabelOn = False
resa@gsnContourZeroLineThicknessF = 2.
resa@cnLineThicknessF = 1.
resa@gsnCenterString = ""
resa@gsnLeftString = ""
resa@gsnRightString = ""
base = gsn_csm_contour_map_ce(wks,sf(0,0,:,:),resa)
end
................................................
Script for Fig. A1(b)
................................................
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
; =======
tem=addfile("data/proc2_moc_final/atl_moc_hel17153-DZA_TP6ML40_mpiom-trunk_mpiom_moc_2003-01-01_2012-12-31.nc","r")
AMOC=tem->atlantic_moc
AMOC@dep=tem->depth
AMOC@lat1d=tem->lat
AMOC@_FillValue = 9.97e+36
AMOC = where(AMOC.lt.-1000000000,AMOC@_FillValue,AMOC)
tem2=addfile("data/proc2_moc_final/glo_moc_hel17153-DZA_TP6ML40_mpiom-trunk_mpiom_moc_2003-01-01_2012-12-31.nc/","r")
GMOC=tem2->global_moc
GMOC@dep=tem2->depth
GMOC@lat1d=tem2->lat
GMOC = where(GMOC.lt.-1000000000,GMOC@_FillValue,GMOC)
printVarSummary(AMOC)
;=======
wks = gsn_open_wks("pdf","figures/AMOC_NCL")
gsn_define_colormap(wks,"BlueDarkRed18")
gray = NhlNewColor(wks,0.83,0.83,0.83)
resa = True
resa@gsnMaximize = True
resa@cnFillOn = True
resa@cnLinesOn = True
resa@cnLineLabelsOn = True
resa@cnInfoLabelOn = False
resa@cnLabelMasking = True
resa@cnLineLabelBackgroundColor = "transparent"
resa@cnLineLabelPlacementMode = "constant"
resa@lbLabelBarOn = True ;False
resa@lbOrientation = "Vertical"
resa@lbLabelStride = 3
resa@gsnSpreadColors = True
resa@cnMissingValFillPattern = 0
resa@cnMissingValFillColor = "grey"
resa@cnFillMode = "RasterFill"
resa@trGridType = "TriangularMesh"
resa@cnLevelSelectionMode = "ManualLevels"
resa@cnMinLevelValF = -24.
resa@cnMaxLevelValF = 24.
resa@cnLevelSpacingF = 2.
resa@gsnCenterString = ""
resa@gsnLeftString = ""
resa@gsnRightString = ""
resa@sfXArray = AMOC@lat1d
resa@sfYArray = AMOC@dep
resa@gsnYAxisIrregular2Linear = True ; converts irreg depth to linear
resa@tiYAxisString = "depth (m)"
resa@trYReverse = True ; reverse y-axis
resa@gsnContourZeroLineThicknessF = 2.
resa@cnLineThicknessF = 1.
plot = gsn_csm_contour(wks,AMOC(0,:,:,0),resa)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.A2 (NCL):
................................................
script to obtain the 10-year mean & zonal mean of
T/S from STORMTIDE2 simulation
................................................
#!/bin/bash
dd=(31 28 31 30 31 30 31 31 30 31 30 31)
dd2=(31 29 31 30 31 30 31 31 30 31 30 31)
dir=/work/mh0256/m300212/STORMTIDE_run/Tides_NCEP_SAL2/experiments/hel17153-DZA_TP6ML40_mpiom-trunk_12yr/outdata/data_mm_3d_2d/
for yr in {10..12}; do
for mm in {1..12..3}; do
mm1=$(printf "%02d" $mm)
mm2=$(printf "%02d" $((mm+2)))
echo "*****************" ${yr} " -- " ${mm1} "---" ${mm2}
fname=${dir}hel17153-DZA_TP6ML40_mpiom-trunk_20${yr}_mpiom_data_3d_mm_20${yr}-${mm1}-01_20${yr}-${mm2}-${dd[$((mm+1))]}.nc
echo $fname
outname=data/proc1_tp6m_TS/hel17153-DZA_STORMTIDE_NCEP_T-S_mm_20${yr}-${mm1}-20${yr}-${mm2}.nc
echo $outname
cdo selcode,2,5 $fname $outname
done
done
# combine all years
echo "merge files"
cdo mergetime data/proc1_tp6m_TS/hel17153-DZA_STORMTIDE_NCEP*.nc data/proc1_tp6m_TS/TS_STORMTIDE_NCEP_mergetime_2003-2012.nc
# time mean
echo "time mean"
cdo timmean data/proc1_tp6m_TS/TS_STORMTIDE_NCEP_mergetime_2003-2012.nc data/proc1_tp6m_TS/TS_STORMTIDE_NCEP_2003-2012_timmean.nc
# remap to lonlat grid
echo "remap"
cdo remapbil,griddes data/proc1_tp6m_TS/TS_STORMTIDE_NCEP_2003-2012_timmean.nc data/proc1_tp6m_TS/TS_1x1_STORMTIDE_NCEP_2003-2012_timmean.nc
# zonal mean
echo "zonal mean"
targetLevels=$(cdo -s showlevel /work/mh0256/m300212/origin_data/Observations/PHC/PHC__3.0__TempO__1x1__annual.nc | tr ' ' ',')
cdo zonmean -intlevel$targetLevels data/proc1_tp6m_TS/TS_1x1_STORMTIDE_NCEP_2003-2012_timmean.nc data/proc1_tp6m_TS/TS_STORMTIDE_NCEP_tmmean_zonmean_2003-2012_33levels.nc
echo "finished!!!"
................................................
script for A2 plot (NCL)
................................................
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "/sw/rhel6-x64/ncl-6.3.0-gccsys/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin
;=======
; import PHC3.0 data
temT=addfile("/work/mh0256/m300212/origin_data/Observations/PHC/PHC__3.0__TempO__1x1__annual_zonmean.nc","r")
temS=addfile("/work/mh0256/m300212/origin_data/Observations/PHC/PHC__3.0__SO__1x1__annual_zonmean.nc","r")
PHC_T = temT->TempO
PHC_T@_FillValue = 9.27e+36
PHC_T = where(PHC_T.eq.-999,PHC_T@_FillValue,PHC_T)
PHC_T@lat = temT->lat
PHC_T@depth = temT->lev
PHC_S = temS->SO
PHC_S@_FillValue = 9.27e+36
PHC_S = where(PHC_S.eq.-999,PHC_S@_FillValue,PHC_S)
PHC_S@lat = temS->lat
PHC_S@depth = temS->lev
printVarSummary(PHC_T)
printVarSummary(PHC_S)
data=addfile("data/proc1_tp6m_TS/TS_STORMTIDE_NCEP_tmmean_zonmean_2003-2012_33levels.nc","r")
Mod_T=data->thetao
Mod_T@lat=data->lat
Mod_T@depth=data->depth
Mod_T@_FillValue = 9.27e+36
Mod_T = where(Mod_T.lt.-9e10,Mod_T@_FillValue,Mod_T)
Mod_S=data->so
Mod_S@lat=data->lat
Mod_S@depth=data->depth
Mod_S@_FillValue = 9.27e+36
Mod_S = where(Mod_S.lt.-9e10,Mod_S@_FillValue,Mod_S)
printVarSummary(Mod_T)
printVarSummary(Mod_S)
diffT = Mod_T(0,:,:,0) - PHC_T(:,:,0)
diffT@lat=data->lat
diffT@depth=data->depth
printVarSummary(diffT)
diffS = Mod_S(0,:,:,0) - PHC_S(:,:,0)
diffS@lat=data->lat
diffS@depth=data->depth
printVarSummary(diffS)
wks = gsn_open_wks("pdf","figures/T-S-bias_STORMTIDE-VS-PHC_NCL")
gsn_define_colormap(wks,"BlWhRe")
plot = new(2,graphic)
resa = True
resa@gsnDraw = False
resa@gsnFrame = False
resa@cnFillOn = True
resa@cnLinesOn = True
resa@cnLineLabelsOn = False ;True
resa@lbLabelBarOn = True
resa@lbOrientation = "Vertical"
resa@lbLabelStride = 2
resa@gsnSpreadColors = True
resa@gsnSpreadColorEnd = -2 ; don't include gray in the end
resa@mpFillDrawOrder = "PostDraw" ; draw map fill last
resa@sfXArray = data->lat
resa@sfYArray = data->depth
resa@vpWidthF = 0.5 ; change aspect ratio of plot
resa@vpHeightF = 0.3
resa@cnMissingValFillPattern = 0
resa@cnMissingValFillColor = "grey" ;;Color of miss values
resa@gsnYAxisIrregular2Linear = True
resa@tiYAxisString = "depth (m)"
resa@trYReverse = True ; reverses y-axis
resa@tmXBMode = "Explicit"
resa@tmXBValues = (/-60,-30,0,30,60/)
resa@tmXBLabels =(/"60S","30S","0","30N","60N"/)
resa@tmXBMinorValues = ispan(-90,90,10)
resa@gsnCenterString = ""
resa@gsnLeftString = ""
resa@gsnRightString = ""
resDT = resa
resDT@tiXAxisOn = False
resDT@tmXBLabelsOn = False
resDT@cnLevelSelectionMode = "ManualLevels"
resDT@cnMinLevelValF = -4.25
resDT@cnMaxLevelValF = 4.25
resDT@cnLevelSpacingF = 0.5
plot(0) = gsn_csm_contour(wks,diffT,resDT)
resDS = resa
resDS@cnLevelSelectionMode = "ManualLevels"
resDS@cnMinLevelValF = -0.5
resDS@cnMaxLevelValF = 0.5
resDS@cnLevelSpacingF = 0.04
plot(1) = gsn_csm_contour(wks,diffS,resDS)
;; create panel
resP = True
resP@gsnFrame = False
resP@gsnPanelLabelBar = False
gsn_panel(wks,plot,(/2,1/),resP)
frame(wks)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Fig.B1 (NCL):
................................................
Script for Fig. B1(a)
................................................
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
; =====
harm = (/"m2","s2","n2","k2","k1","o1","p1","q1"/)
rng = (/80,35,20,10,35,25,11,5/)
intv = (/5,2,2,1,2,2,1,0.5/)
in = addfile("data/proc1_harmonics_LSF/harmonics_2012yr_zos_SAL2.nc","r")
lon2d = in->lon
lat2d = in->lat
do i = 1,1 ; 8
print(harm(i-1))
ampname = harm(i-1)+"amp"
phaname = harm(i-1)+"phase"
amp = in->$ampname$
pha = in->$phaname$
amp@lon2d = lon2d
amp@lat2d = lat2d
pha@lon2d = lon2d
pha@lat2d = lat2d
amp = amp*100 ;;; unit in cm ;;;
amp@_FillValue = -9999.0
pha@_FillValue = -9999.0
amp = where(amp.gt.0,amp,amp@_FillValue)
pha = where(pha.gt.0,pha,pha@_FillValue)
wks = gsn_open_wks("pdf","figures/proc3_harmonics/"+harm(i-1)+"AmpPhase_LSF_zos_2012yr_sal2_stormtide-ncep")
gsn_define_colormap(wks,"WhBlGrYeRe") ; define colormap type
gray = NhlNewColor(wks,0.83,0.83,0.83)
; ======
resp = True ; resource for phase plot
resp@cnLevelSelectionMode = "ManualLevels" ; max, min and spacing for plot
resp@cnMinLevelValF = 0.
resp@cnMaxLevelValF = 360.
resp@cnLevelSpacingF = 45.
resp@cnInfoLabelOn = False ; display information label
resp@cnLinesOn = True ; display the lines and their labels
resp@cnLineLabelsOn = False ; True
resp@cnLineLabelInterval = 1 ; interval for displaying labels
resp@cnLineColor = "black" ; color for lines and labels
resp@cnLineLabelFontColor = "black"
resp@cnLineLabelFontHeightF = 0.006 ; font size for labels
resp@cnLabelMasking = True ; make a white box for font of labels
resp@gsnDraw = False ; cancle displaying draw
resp@gsnFrame = False ; cancle displaying frame
resp@gsnLeftString = ""
resp@gsnRightString = "" ;"STORMTIDE/NCEP, 2012, SAL=0.09"
plot = gsn_csm_contour(wks,pha(0,0,:,:),resp) ; plot for phase
;=====
resa = True ; resources for amplitude
resa@cnFillOn = True ; fill in the contour
resa@cnLinesOn = False ; cancle displaying lines
resa@gsnSpreadColors = True ;
resa@lbLabelAutoStride = "True" ; automatically set colormap stride
resa@cnLevelSelectionMode = "ManualLevels" ; max, min and spacing for plot
resa@cnMinLevelValF = 0.
resa@cnMaxLevelValF = rng(i-1)
resa@cnLevelSpacingF = intv(i-1)
resa@gsnLeftString = ""
resp@gsnRightString = ""
resa@gsnAddCyclic = False
resa@cnSmoothingOn = True
resa@mpMaxLatF = 90
resa@mpMaxLonF = 360
resa@mpMinLatF = -90
resa@mpMinLonF = 0
resa@mpCenterLonF = 180
resa@gsnSpreadColorEnd = -2 ; don't include gray in the end
resa@mpFillDrawOrder = "PostDraw" ; draw map fill last
resa@gsnDraw = False ; cancle displaying draw
resa@gsnFrame = False ; cancle displaying frame
base = gsn_csm_contour_map_ce(wks,amp(0,0,:,:),resa) ; plot for amplitude
overlay(base,plot) ; overlay two pictures
draw(base) ; display draw
frame(wks) ; display frame
end do
end
................................................
Script for Fig. B1(b)
................................................
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
; ======
harm = (/"m2","s2","n2","k2","k1","o1","p1","q1"/)
rng = (/80,35,20,10,35,25,11,5/)
intv = (/5,2,2,1,2,2,1,0.5/)
do i = 1,1 ; 8
print(harm(i-1))
in = addfile("/work/mh0256/m300212/origin_data/Observations/HAMTIDE/"+harm(i-1)+".hamtide11a.nc","r")
amp = in->AMPL ;;; unit in cm ;;;
amp@lon1d = in->LON
amp@lat1d = in->LAT
amp@_FillValue = -9999.0
amp = where(amp.gt.0,amp,amp@_FillValue)
pha = in->PHAS
pha@lon1d = in->LON
pha@lat1d = in->LAT
pha@_FillValue = -9999.0
pha = where(pha.gt.0,pha,pha@_FillValue)
wks = gsn_open_wks("pdf","figures/proc3_harmonics/AmpPhs_HAMTIDE_zos_"+harm(i-1))
gsn_define_colormap(wks,"WhBlGrYeRe") ; define colormap type
gray = NhlNewColor(wks,0.83,0.83,0.83)
; ======
resp = True ; resource for phase plot
resp@cnLevelSelectionMode = "ManualLevels" ; max, min and spacing for plot
resp@cnMinLevelValF = 0.
resp@cnMaxLevelValF = 360.
resp@cnLevelSpacingF = 45.
resp@cnInfoLabelOn = False ; display information label
resp@cnLinesOn = True ; display the lines and their labels
resp@cnLineLabelsOn = True
resp@cnLineLabelInterval = 1 ; interval for displaying labels
resp@cnLineColor = "black" ; color for lines and labels
resp@cnLineLabelFontColor = "black"
resp@cnLineLabelFontHeightF = 0.006 ; font size for labels
resp@cnLabelMasking = True ; make a white box for font of labels
resp@gsnDraw = False ; cancle displaying draw
resp@gsnFrame = False ; cancle displaying frame
resp@gsnLeftString = ""
resp@gsnRightString = "" ;"STORMTIDE/NCEP, 2012, SAL=0.09"
plot = gsn_csm_contour(wks,pha,resp) ; plot for phase
resa = True ; resources for amplitude
resa@cnFillOn = True ; fill in the contour
resa@cnLinesOn = False ; cancle displaying lines
resa@gsnSpreadColors = True ;
resa@lbLabelAutoStride = "True" ; automatically set colormap stride
resa@cnLevelSelectionMode = "ManualLevels" ; max, min and spacing for plot
resa@cnMinLevelValF = 0.
resa@cnMaxLevelValF = rng(i-1)
resa@cnLevelSpacingF = intv(i-1)
resa@gsnLeftString = ""
resa@gsnRightString = ""
resa@gsnAddCyclic = False
resa@cnSmoothingOn = True
resa@mpMaxLatF = 90
resa@mpMaxLonF = 360
resa@mpMinLatF = -90
resa@mpMinLonF = 0
resa@mpCenterLonF = 180
resa@gsnSpreadColorEnd = -2 ; don't include gray in the end
resa@mpFillDrawOrder = "PostDraw" ; draw map fill last
resa@gsnDraw = False ; cancle displaying draw
resa@gsnFrame = False ; cancle displaying frame
base = gsn_csm_contour_map_ce(wks,amp,resa) ; plot for amplitude
overlay(base,plot) ; overlay two pictures
draw(base) ; display draw
frame(wks) ; display frame
end do
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Script for Tab.1 (Matlab):
................................................
% main program
................................................
clear all; clc;
harm = {'m2','s2','n2','k2','k1','o1','p1','q1'};
dir102 = '/work/mh0256/m300212/origin_data/Observations/st102/';
temdir = '/work/mh0256/m300212/STORMTIDE_run/Tides_NCEP_SAL2/Data_analysis/proc1_harmonics/NCEP_SAL_2012yr/';
fname_NCEP = [temdir 'data/proc1_harmonics_LSF/harmonics_2012yr_zos_SAL2.nc'];
clear temdir;
disp('import NCEP data')
for iharm = 1 : numel(harm)
temamp = ncread(fname_NCEP,[harm{iharm} 'amp']);
temamp(temamp==0) = nan;
amp_ncep{iharm} = temamp;
clear temamp;
temphs = ncread(fname_NCEP,[harm{iharm} 'phase']);
temphs(temphs==0) = nan;
phs_ncep{iharm} = temphs;
clear temphs;
if iharm == 1
Nlon = ncread(fname_NCEP,'lon');
Nlon(Nlon<0) = Nlon(Nlon<0) + 360;
Nlat = ncread(fname_NCEP,'lat');
end
end
disp('import 102 gauges')
for iharm = 1 : numel(harm)
fname_102 = [dir102 'data/Pegel_' harm{iharm} '_102.asci'];
st102 = load(fname_102);
if iharm == 1
lat102 = st102(:,1); % [-61.5 60.5]
lon102 = st102(:,2);
lon102(lon102<0) = lon102(lon102<0) + 360; % [0 360]
end
temA = sqrt(st102(:,3).^2 + st102(:,4).^2)./100; % [m]
temA(temA==0) = nan;
amp102(iharm,:) = temA;
temphs = atan2(st102(:,4),st102(:,3))*180/pi;
temphs(temphs<0) = temphs(temphs<0) + 360;
temphs(isnan(temA)) = nan;
phs102(iharm,:) = temphs;
clear temA temphs;
clear fname_102 st102;
end
dir = '/work/mh0256/m300212/origin_data/Observations/st102/';
fname = [dir 'stationij.mat'];
data = load(fname);
row102 = data.jpos_st;
col102 = data.ipos_st;
clear dir fname data;
disp('cal of RMS')
for iharm = 1 : numel(harm)
[sigN02(iharm),rmsN02(iharm),rat1_N02(iharm),relD1_N02(iharm),rms2N02(iharm),rat2_N02(iharm),relD2_N02(iharm),rms2A_N02(iharm),rat3_N02(iharm),relD3_N02(iharm)] = calcRMS(lon102,lat102,amp102(iharm,:),phs102(iharm,:),Nlon,Nlat,amp_ncep{iharm},phs_ncep{iharm},row102,col102);
end
rssN_g02 = 100 * (1 - (sqrt(sum(rms2N02.^2))^2) / (sqrt(sum(sigN02.^2))^2) )
disp('output file')
clear ofname;
ofname = 'Signal_RMSdiff_Ratio_cmp_tidal-gauges_VS_model_sal2_2012yr.dat';
fid = fopen(ofname,'w');
fprintf(fid,'%s\n','STORMTIDE_NCEP with sal=0.09 (2012yr)')
fprintf(fid,'%s\n',' ')
fprintf(fid,'%s\n',[' ' harm{1} ' ' harm{2} ' ' harm{3} ' ' harm{4} ' ' harm{5} ' ' harm{6} ' ' harm{7} ' ' harm{8}]);
fprintf(fid,'%s ','signal-102[cm] : ');
fprintf(fid,'%.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n',sigN02.*100); % output of [cm]
fprintf(fid,'%s\n',' ');
fprintf(fid,'%s\n','NCEP run VS 102 gauges');
fprintf(fid,'%s','RMS:total[cm] : ');
fprintf(fid,'%.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n',rms2N02.*100); % [cm]
fprintf(fid,'%s','RSSerror (full amp/phase): ');
fprintf(fid,'%.2f\n',rssN_g02);
................................................
% sub-function for calcRMS.m
................................................
function [signal,RMSD1,ratio1,relD1,RMSD2,ratio2,relD2,RMSD2_amp,ratio3,relD3,out_rownum,out_colnum] = calcRMS(lon_st,lat_st,amp_st,phs_st,lon2d,lat2d,amp_md,phs_md,rownum,colnum)
% DO it for one tidal constituent at each time
%
% make sure lon_st&lon2d, lat_st&lat2d are in the same range with the same
% units
% -- INPUT -- %
% lon_st, lat_st: lon/lat of the station locations
% amp_st, phs_st: tidal constants (amp/phase) of the respective tidal
% constituents
%
% lon2d, lat2d: 2D lon/lat data of the model
% amp_md, phs_md: tidal constants of different tidal constituents by
% harmonic analysis of the simulation
%
% rownum, colnum: [optional] location of obs stations in 3600x2394 model grides.
if exist('rownum','var')
disp('rownum exists')
else
disp('get rownum of stations')
rownum = zeros(length(lon_st),1);
colnum = zeros(length(lon_st),1);
for ipnt = 1 : length(lon_st)
dist = sqrt( (lon2d-lon_st(ipnt)).^2 + (lat2d-lat_st(ipnt)).^2 );
[val1,pos1] = min(dist);
[~,pos2] = min(val1);
rownum(ipnt) = pos1(pos2);
colnum(ipnt) = pos2;
end
out_rownum = rownum;
out_colnum = colnum;
end
Am = zeros(length(lon_st),1);
As = zeros(length(lon_st),1);
Pm = zeros(length(lon_st),1);
Ps = zeros(length(lon_st),1);
for ipnt = 1 : length(lon_st)
Am(ipnt) = amp_md(rownum(ipnt),colnum(ipnt));
Pm(ipnt) = phs_md(rownum(ipnt),colnum(ipnt));
As(ipnt) = amp_st(ipnt);
Ps(ipnt) = phs_st(ipnt);
end
% ======
% estimate of signal from tidal gauge observations (Stammer et al., 2014)
temsig = As.^2;
signal = sqrt( 0.5*nansum(temsig)/length(find(~isnan(temsig))) );
% ======
% estimate of RMS differences (malte's method)
tem_rms = 0.5*(Am.^2+As.^2) - Am.*As .*cos(degtorad(Ps-Pm));
RMSD1 = sqrt( nansum(tem_rms)/length(find(~isnan(tem_rms))) );
ratio1 = 100 * (1-RMSD1^2/signal^2);
relD1 = 100*RMSD1/signal;
% ======
% estimate of RMS differences (Stammer et al., 2014)
RMSD_cycle = zeros(length(lon_st),1);
for ipnt = 1 : length(lon_st)
wt = 0 : 0.01 : 2*pi;
RMSD_cycle(ipnt) = 0;
for iwt = 1 : length(wt)
RMSD_cycle(ipnt) = RMSD_cycle(ipnt) + ( As(ipnt)*cos(wt(iwt)-degtorad(Ps(ipnt))) - Am(ipnt) * cos(wt(iwt)-degtorad(Pm(ipnt))) )^2;
end
RMSD_cycle(ipnt) = RMSD_cycle(ipnt)/length(wt);
end
RMSD2 = sqrt( nansum(RMSD_cycle)/length(find(~isnan(RMSD_cycle))) );
ratio2 = 100 * (1-RMSD2^2/signal^2);
relD2 = 100*RMSD2/signal;
% ======
% estimate RMS using Amp only
temdata = (As-Am).^2;
RMSD2_amp = sqrt(nansum(temdata)/length(find(~isnan(temdata))));
ratio3 = 100 * (1-RMSD2_amp^2/signal^2);
relD3 = 100 * RMSD2_amp/signal;