Warning, /utils/matlab/cs_grid/latloncap/plot_llc.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
6a75a3b8ef Chri*0001 % Make some color scale and vector plots for a gmt mesh
0002 % Functions used -
0003 % rdnctiles - to read in mnc generated tile files. these are read into a
0004 % structure containing each tile separately.
0005 %
0006 % croptiled - this function takes a rdnctiles structure and crops the
0007 % data held according to a set of crop specifications
0008 % given as an argument. it returns an updated rdnctiles style
0009 % structure containing each tile separately.
0010 %
0011 % calcFacetAngles - calculate cos and sin terms needed to rotate vectors
0012 % onto underlying coordinate system.
0013 %
0014 % zmean - routine for calculating vertical average.
0015 %
0016 % vel2latlon - rotate velocity vectors according to angle factors
0017 % previously calculated.
0018 %
0019 % llc_pcol - plots a scalar map.
0020 %
0021 % llc_vec - plots a vector map.
0022
0023 % Set what to get and where to get it from
0024 % gfpat - where to get grid files
0025 % dfpat - where to get output files
0026 % flist - fields to get
0027 % tlev - timestep number
0028 % klo, khi - restrict indices in vertical to this range
0029 gfpat='../llc_grid_72tiles/*nc'; dfpat='../llc_48month/*nc';
0030 flist={'dyG', 'dxG', 'XG','YG','XC', 'YC', 'Z','Zl','Ttave','uVeltave','vVeltave'};
0031 tlev=172800;
0032 klo=1; khi=11;
0033
0034 % Load tiles
0035 % ( The example below includes a special mod to get rid of tiles in the llc 72 tile mesh. )
0036 % ( The Antarctic grid tiles (69,70,71,72) contain junk, so get rid of them. )
0037 t=rdnctiles({dfpat,gfpat},flist,[tlev],'bytile',[]);
0038 t=t(1:68);
0039
0040 % Useful commands!
0041 % This command gives the max over all separate tiles.
0042 % max(cellfun(@(x)max(x(:)),(arrayfun(@(x)(x.uVeltave),(arrayfun(@(x)(x.var),t)),'UniformOutput',0))))
0043 % This command extracts a great big long vector from all the separate tiles.
0044 % a=cell2mat((cellfun(@(x)(x(:)),(arrayfun(@(x)(x.uVeltave),(arrayfun(@(x)(x.var),t)),'UniformOutput',0)),'UniformOutput',0)));
0045 % max(a(:));
0046
0047 % Crop fields to desired set of points
0048 % (cspec should be incorporated into rdnctiles to do cropping on the fly).
0049 % croptiled() is a function that is used to crop tiles datastructures returned by rdnctiles.
0050 cspec.crops(1).ranks=3;
0051 cspec.crops(1).indexlo=klo; cspec.crops(1).indexhi=khi;
0052 cspec.crops(1).fields={'uVeltave', 'vVeltave'};
0053 cspec.crops(2).ranks=3;
0054 cspec.crops(2).indexlo=1; cspec.crops(2).indexhi=1;
0055 cspec.crops(2).fields={'Ttave'};
0056 t=croptiled(t,cspec);
0057
0058 %
0059 % Calculate rotation angle to get latlon flow
0060 for tn=1:length(t)
0061 [t(tn).var.ca,t(tn).var.sa,t(tn).var.ua,t(tn).var.va] = ...
0062 calcFacetAngles( t(tn).var.YG, t(tn).var.dxG, t(tn).var.dyG );
0063 end
0064
0065 % Calculate vertical means
0066 for tn=1:length(t)
0067 t(tn).var.ubz=zmean(t(tn).var.uVeltave,t(tn).var.Zl);
0068 t(tn).var.vbz=zmean(t(tn).var.vVeltave,t(tn).var.Zl);
0069 end
0070
0071 % Calculate lat-lon projection of vertical mean flow
0072 for tn=1:length(t)
0073 ca=t(tn).var.ca; sa=t(tn).var.sa;
0074 uvec=t(tn).var.ubz; vvec=t(tn).var.vbz;
0075 [u_ll, v_ll] = vel2latlon(ca, sa, uvec, vvec);
0076 t(tn).var.ubz_ll=u_ll; t(tn).var.vbz_ll=v_ll;
0077 end
0078
0079 % Plot a scalar
0080 figure(1);clf
0081 hc=llc_pcol(t,'ubz_ll',1,'sphere');
0082
0083 % Plot some vector fields
0084 % == Arctic view
0085 figure(2);clf
0086 fspec.proj='cyl';
0087 fspec.lathi= 90; fspec.latlo= 45;
0088 fspec.lonlo=-180; fspec.lonhi=+180;
0089 fspec.az=+90; fspec.el=+90;
0090 fspec.alx = 3.8e3; fspec.aly =-3.8e3;
0091 fspec.axmag=0.0; fspec.aymag=0.5;
0092 fspec.scal=3.;
0093 hv=llc_vec(t,'ubz_ll','vbz_ll',fspec);
0094
0095 % == North Atlantic sector
0096 figure(3);clf
0097 fspec.proj='cyl';
0098 fspec.lathi= 70; fspec.latlo= 20;
0099 fspec.lonlo=-110; fspec.lonhi= -20;
0100 fspec.az= 0; fspec.el=+90;
0101 fspec.alx = 4.2e3; fspec.aly = -5.0e3;
0102 fspec.axmag=0.5; fspec.aymag=0.;
0103 fspec.scal=3.;
0104 hv=llc_vec(t,'ubz_ll','vbz_ll',fspec);
0105
0106
0107 % == South Atlantic sector
0108 figure(4);clf
0109 fspec.proj='cyl';
0110 fspec.lathi= -25; fspec.latlo= -70;
0111 fspec.lonlo= -80; fspec.lonhi= 40;
0112 fspec.az=+90; fspec.el=-90;
0113 fspec.alx =1.0e3; fspec.aly = 2.0e3;
0114 fspec.axmag=0.0; fspec.aymag=0.5;
0115 fspec.scal=3.;
0116 hv=llc_vec(t,'ubz_ll','vbz_ll',fspec);
0117
0118 % == Equatorial Pacific (define in two pieces)
0119 figure(5);clf
0120 fspec.proj='cart';
0121 fspec.lathi= [25 25]; fspec.latlo= [ -25 -25];
0122 fspec.lonlo= [-180 100]; fspec.lonhi=[ -65 180];
0123 fspec.az=+180; fspec.el=0;
0124 fspec.alx =130; fspec.aly = -20;
0125 fspec.axmag=0.5; fspec.aymag=0.;
0126 fspec.scal=3.;
0127 fspec.skip=1;
0128 hv=llc_vec(t,'ubz_ll','vbz_ll',fspec);
0129
0130 % == Antarctic view
0131 figure(6);clf
0132 fspec.proj='cyl';
0133 fspec.lathi= -22; fspec.latlo= -85;
0134 fspec.lonlo=-180; fspec.lonhi=180;
0135 fspec.az=+90; fspec.el=-90;
0136 fspec.alx =5.0e3; fspec.aly = -5.0e3;
0137 fspec.axmag=0.0; fspec.aymag=0.5;
0138 fspec.scal=3.;
0139 hv=llc_vec(t,'ubz_ll','vbz_ll',fspec);
