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(2) || 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<INPUT<INPUT<INPUT<=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<INPUT<INPUT<INPUT<INPUT<INPUT<INPUT<INPUT<INPUT<INPUT<INPUT<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<INPUT<INPUT<INPUT<INPUT<INPUT<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<INPUT<> 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: $(&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;