Warning, /verification/tutorial_global_oce_latlon/diags_matlab/ncep2global_ocean.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
051ee7f715 Jean*0001 function ncep2global_ocean
0002 % function ncep2global_ocean
0003
0004 % read various fluxes, adds them up to have net heat flux and net
0005 % freshwater flux (E-P only)
0006 % interpolates the fluxes onto 90x40 grid of global_ocean.90x40x15
0007 %
0008 % constants
0009 ql_evap = 2.5e6; % latent heat due to evaporation
0010 rho_fresh = 1000; % density of fresh water
0011 onemonth = 60*60*24*30;
0012 % grid parameter
0013 grd = mit_loadgrid(DIRECTORY WHERE THE GRID INFORMATION IS HELD ...
0014 (XC*, YC*, HFAC* etc.));
0015 lonc = grd.lonc;
0016 latc = grd.latc;
0017 cmask = grd.cmask; % 1=wet point, NaN = dry point
0018 rac = grd.rac;
0019
0020 bdir = DIRECTORY WHERE THE NCEP DATA FILES ARE STORED;
0021 % taux
0022 ncload(fullfile(bdir,'ncep_taux_monthly.cdf'),'taux','X','Y','T')
0023 % tauy
0024 ncload(fullfile(bdir,'ncep_tauy_monthly.cdf'),'tauy')
0025 % there appears to be a different east west convention in the data
0026 taux = - taux;
0027 tauy = - tauy;
0028
0029 % qnet
0030 ncload(fullfile(bdir,'ncep_netsolar_monthly.cdf'),'solr')
0031 ncload(fullfile(bdir,'ncep_netlw_monthly.cdf'),'lwflx') % long wave radiation
0032 ncload(fullfile(bdir,'ncep_latent_monthly.cdf'),'heat_flux');
0033 lhflx = heat_flux;
0034 ncload(fullfile(bdir,'ncep_sensible_monthly.cdf'),'heat_flux');
0035 shflx = heat_flux;
0036
0037 % emp
0038 ncload(fullfile(bdir,'ncep_precip_monthly.cdf'),'prate')
0039 ncload(fullfile(bdir,'ncep_runoff_monthly.cdf'),'runoff')
0040 evap = lhflx/ql_evap/rho_fresh; % evaporation rate estimated from latent
0041 % heat flux
0042 precip = prate/rho_fresh; % change the units
0043
0044 % add the fields
0045 % net heat flux: sum of solar radiation, long wave flux, latent heat flux,
0046 % and sensible heat flux
0047 qnet = solr + lwflx + lhflx + shflx;
0048 % net fresh water flux: difference between evaporation and precipitation
0049 emp = evap - precip;
0050 % substract the runoff (doesn't lead anywhere, unfortunately)
0051 empmr = emp - change(runoff,'==',NaN,0)/rho_fresh/onemonth;
0052 runoff(find(isnan(runoff)))=0;
0053 empmr = emp - runoff/rho_fresh/onemonth;
0054
0055 % interpolate the three fields onto 90x40 grid
0056 ir = interp2global(X,Y,lonc,latc); % this sets up the
0057 % interpolator ir!
0058 for k = 1:length(T);
0059 tauxg(:,:,k) = interp2global(X,Y,squeeze(taux(k,:,:)),lonc,latc,ir);
0060 tauyg(:,:,k) = interp2global(X,Y,squeeze(tauy(k,:,:)),lonc,latc,ir);
0061 qnetg(:,:,k) = interp2global(X,Y,squeeze(qnet(k,:,:)),lonc,latc,ir);
0062 empg(:,:,k) = interp2global(X,Y,squeeze(emp(k,:,:)),lonc,latc,ir);
0063 empmrg(:,:,k) = interp2global(X,Y,squeeze(empmr(k,:,:)),lonc,latc,ir);
0064 end
0065
0066 % apply landmask
0067 for k = 1:length(T);
0068 tauxg(:,:,k) = tauxg(:,:,k).*cmask(:,:,1);
0069 tauyg(:,:,k) = tauyg(:,:,k).*cmask(:,:,1);
0070 qnetg(:,:,k) = qnetg(:,:,k).*cmask(:,:,1);
0071 empg(:,:,k) = empg(:,:,k).*cmask(:,:,1);
0072 empmrg(:,:,k) = empmrg(:,:,k).*cmask(:,:,1);
0073 end
0074 % balance the fluxes
0075 effrac = rac.*cmask(:,:,1);
0076 mqnet = mit_mean(qnetg,effrac);
0077 qnetb = qnetg-mean(mqnet);
0078
0079 memp = mit_mean(empg,effrac);
0080 empb = empg-mean(memp);
0081 mempmr = mit_mean(empmrg,effrac);
0082 empmrb = empmrg-mean(mempmr);
0083
0084 % apply the landmasks one more time (because we have added something
0085 % to what should have been zero == land)
0086 for k = 1:length(T);
0087 qnetb(:,:,k) = qnetb(:,:,k).*smask;
0088 empb(:,:,k) = empb(:,:,k).*smask;
0089 empmrb(:,:,k) = empmrb(:,:,k).*smask;
0090 end
0091
0092 % save the fields
0093 acc = 'real*4';
0094 mit_writefield('ncep_taux.bin',change(tauxg,'==',NaN,0),acc);
0095 mit_writefield('ncep_tauy.bin',change(tauyg,'==',NaN,0),acc);
0096 mit_writefield('ncep_qnet.bin',change(qnetb,'==',NaN,0),acc);
0097 mit_writefield('ncep_emp.bin',change(empb,'==',NaN,0),acc);
0098 mit_writefield('ncep_empmr.bin',change(empmrb,'==',NaN,0),acc);
0099
0100 return
0101
0102 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0103
0104 function Y = mit_mean(X,area);
0105 %function Y = mit_mean(X,area);
0106
0107 na = ndims(area);
0108 n = ndims(X);
0109 full_area = nansum(area(:));
0110 if n == 2 & na == 2
0111 [nx ny] = size(X);
0112 Xa = X.*area;
0113 Y = nansum(Xa(:))/full_area;
0114 elseif n == 3 & na == 2
0115 [nx ny nt] = size(X);
0116 Xa = X.*repmat(area,[1 1 nt]);
0117 for k=1:nt
0118 Y(k) = nansum(nansum(Xa(:,:,k)))/full_area;
0119 end
0120 elseif n == 3 & na == 3
0121 [nx ny nz] = size(X);
0122 Xa = X.*area;
0123 Y = nansum(Xa(:))/full_area;
0124 else
0125 error('X and area have to have consistent dimensions')
0126 end
0127
0128 return
0129
0130 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0131
0132 function count = mit_writefield(fname,h,accuracy)
0133 %function count = mit_writefield(fname,h,accuracy)
0134
0135 ieee='ieee-be';
0136
0137 [fid message] = fopen(fname,'w',ieee);
0138 if fid <= 0
0139 error([message ', filename: ', [fname]])
0140 end
0141
0142 clear count
0143 dims = size(h);
0144 if length(dims) == 2
0145 count = fwrite(fid,h,accuracy);
0146 elseif length(dims) == 3
0147 for k=1:dims(3);
0148 count(k) = fwrite(fid,squeeze(h(:,:,k)),accuracy);
0149 end
0150 elseif length(dims) == 4
0151 for kt=1:dims(4)
0152 for k=1:dims(3)
0153 count(k,kt) = fwrite(fid,squeeze(h(:,:,k,kt)),accuracy);
0154 end
0155 end
0156 end
0157
0158 fclose(fid);
0159
0160 return
0161
0162 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0163
0164 function [zi] = interp2global(varargin);
0165 %function [zi] = interp2global(x,y,z,xi,yi,ir);
0166
0167 if nargin == 4
0168 [xx yy] = meshgrid(varargin{1},varargin{2});
0169 xx = reshape(xx,prod(size(xx)),1);
0170 yy = reshape(yy,prod(size(yy)),1);
0171 [xxi yyi] = meshgrid(varargin{3},varargin{4});
0172 xxv = reshape(xxi,prod(size(xxi)),1);
0173 yyv = reshape(yyi,prod(size(yyi)),1);
0174 % create map
0175 if ~exist('ir','var');
0176 r = 230e3; %m
0177 for k = 1:length(xxv);
0178 % find all points within radius r
0179 radius = lonlat2dist(xxv(k),yyv(k),xx,yy);
0180 ir{k} = find(radius<=r);
0181 end
0182 end
0183 zi = ir;
0184 else
0185 % zi = interp2_global(x,y,z,ix,iy','spline')';
0186 z = varargin{3};
0187 xi = varargin{4};
0188 yi = varargin{5};
0189 [xxi yyi] = meshgrid(varargin{4},varargin{5});
0190 ir = varargin{6};
0191 zzv = repmat(NaN,length(ir),1);
0192 % interpolate
0193 for k = 1:length(ir);
0194 % find all points within radius r
0195 zzv(k) = mean(z(ir{k}));
0196 end
0197 zi = reshape(zzv,size(xxi))';
0198 end
0199
0200 return
0201
0202 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0203
0204 function r = lonlat2dist(lon,lat,lons,lats)
0205 %function r = lonlat2dist(lon,lat,lons,lats)
0206 %
0207 % returns distance (in meters) between postion (lon,lat) and positions
0208 % (lons,lats) on the earth (sphere). length(r) = length(lons)
0209
0210 earth=6371000;
0211 deg2rad=pi/180;
0212 nx=length(lon);
0213 lt = deg2rad*lat;
0214 ln = deg2rad*lon;
0215 lts = deg2rad*lats;
0216 lns = deg2rad*lons;
0217 alpha = ...
0218 acos( ( cos(lt)*cos(lts) ).*cos(ln-lns) ...
0219 + ( sin(lt)*sin(lts) ) );
0220 r = earth*abs(alpha');
0221
0222 return
0223
