| 1 | molod | 1.1 | % | 
| 2 |  |  | % Read in CS 510 input fields | 
| 3 |  |  | xc1 = zeros(511,511,6); | 
| 4 |  |  | yc1 = zeros(511,511,6); | 
| 5 |  |  | dxf1 = zeros(511,511,6); | 
| 6 |  |  | dyf1 = zeros(511,511,6); | 
| 7 |  |  | ra1 = zeros(511,511,6); | 
| 8 |  |  | xg1 = zeros(511,511,6); | 
| 9 |  |  | yg1 = zeros(511,511,6); | 
| 10 |  |  | dxv1 = zeros(511,511,6); | 
| 11 |  |  | dyu1 = zeros(511,511,6); | 
| 12 |  |  | raz1 = zeros(511,511,6); | 
| 13 |  |  | dxc1 = zeros(511,511,6); | 
| 14 |  |  | dyc1 = zeros(511,511,6); | 
| 15 |  |  | raw1 = zeros(511,511,6); | 
| 16 |  |  | ras1 = zeros(511,511,6); | 
| 17 |  |  | dxg1 = zeros(511,511,6); | 
| 18 |  |  | dyg1 = zeros(511,511,6); | 
| 19 |  |  | for iface=1:6 | 
| 20 |  |  | fid = fopen(['tile00' num2str(iface) '.mitgrid'],'r','b'); | 
| 21 |  |  | xc1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 22 |  |  | yc1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 23 |  |  | dxf1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 24 |  |  | dyf1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 25 |  |  | ra1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 26 |  |  | xg1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 27 |  |  | yg1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 28 |  |  | dxv1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 29 |  |  | dyu1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 30 |  |  | raz1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 31 |  |  | dxc1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 32 |  |  | dyc1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 33 |  |  | raw1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 34 |  |  | ras1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 35 |  |  | dxg1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 36 |  |  | dyg1(:,:,iface) = fread(fid,[511 511],'real*8'); | 
| 37 |  |  | end | 
| 38 |  |  | % | 
| 39 |  |  | % Now output (approx 1 deg) grid | 
| 40 |  |  | ratio = 5; | 
| 41 |  |  | xc2 = zeros(103,103,6); | 
| 42 |  |  | yc2 = zeros(103,103,6); | 
| 43 |  |  | dxf2 = zeros(103,103,6); | 
| 44 |  |  | dyf2 = zeros(103,103,6); | 
| 45 |  |  | ra2 = zeros(103,103,6); | 
| 46 |  |  | xg2 = zeros(103,103,6); | 
| 47 |  |  | yg2 = zeros(103,103,6); | 
| 48 |  |  | dxv2 = zeros(103,103,6); | 
| 49 |  |  | dyu2 = zeros(103,103,6); | 
| 50 |  |  | raz2 = zeros(103,103,6); | 
| 51 |  |  | dxc2 = zeros(103,103,6); | 
| 52 |  |  | dyc2 = zeros(103,103,6); | 
| 53 |  |  | raw2 = zeros(103,103,6); | 
| 54 |  |  | ras2 = zeros(103,103,6); | 
| 55 |  |  | dxg2 = zeros(103,103,6); | 
| 56 |  |  | anglecos2 = zeros(103,103,6); | 
| 57 |  |  | anglesin2 = zeros(103,103,6); | 
| 58 |  |  | dyg2 = zeros(103,103,6); | 
| 59 |  |  | % Create left, x-center, bottom and y-center indeces into 510 grid | 
| 60 |  |  | ileft=[1:5:510]; | 
| 61 |  |  | icent=ileft+2; | 
| 62 |  |  | jbott=[1:5:510]; | 
| 63 |  |  | jcent=jbott+2; | 
| 64 |  |  | % First do the interior of each face | 
| 65 |  |  | for iface=1:6 | 
| 66 |  |  | % lats and lons | 
| 67 |  |  | xc2(1:102,1:102,iface)=xc1(icent,jcent,iface); | 
| 68 |  |  | yc2(1:102,1:102,iface)=yc1(icent,jcent,iface); | 
| 69 |  |  | xg2(1:102,1:102,iface)=xg1(ileft,jbott,iface); | 
| 70 |  |  | yg2(1:102,1:102,iface)=yg1(ileft,jbott,iface); | 
| 71 |  |  | for ipnt=0:ratio-1 | 
| 72 |  |  | dxf2(1:102,1:102,iface)=dxf2(1:102,1:102,iface) + dxf1(ileft+ipnt,jcent,iface); | 
| 73 |  |  | dxg2(1:102,1:102,iface)=dxg2(1:102,1:102,iface) + dxg1(ileft+ipnt,jbott,iface); | 
| 74 |  |  | dxv2(2:102,1:102,iface)=dxv2(2:102,1:102,iface) + dxv1(icent(1:101)+ipnt+1,jbott(1:102),iface); | 
| 75 |  |  | dxc2(2:102,1:102,iface)=dxc2(2:102,1:102,iface) + dxc1(icent(1:101)+ipnt+1,jcent(1:102),iface); | 
| 76 |  |  | end | 
| 77 |  |  | % dx's and dy's | 
| 78 |  |  | for jpnt=0:ratio-1 | 
| 79 |  |  | dyf2(1:102,1:102,iface)=dyf2(1:102,1:102,iface) + dyf1(icent,jbott+jpnt,iface); | 
| 80 |  |  | dyg2(1:102,1:102,iface)=dyg2(1:102,1:102,iface) + dyg1(ileft,jbott+jpnt,iface); | 
| 81 |  |  | dyu2(1:102,2:102,iface)=dyu2(1:102,2:102,iface) + dyu1(ileft(1:102),jcent(1:101)+jpnt+1,iface); | 
| 82 |  |  | dyc2(1:102,2:102,iface)=dyc2(1:102,2:102,iface) + dyc1(icent(1:102),jcent(1:101)+jpnt+1,iface); | 
| 83 |  |  | end | 
| 84 |  |  | % Construct values that we don't know yet - use grid symmetries for edge dx's and dy's | 
| 85 |  |  | dxc2(1,1:102,iface)=dxc1(1,jcent(1:102),iface) + 2.*dxc1(2,jcent(1:102),iface) + 2.*dxc1(3,jcent(1:102),iface); | 
| 86 |  |  | dxv2(1,1:102,iface)=dxv1(1,jbott(1:102),iface) + 2.*dxv1(2,jbott(1:102),iface) + 2.*dxv1(3,jbott(1:102),iface); | 
| 87 |  |  | dyc2(1:102,1,iface)=dyc1(icent(1:102),1,iface) + 2.*dyc1(icent(1:102),2,iface) + 2.*dyc1(icent(1:102),3,iface); | 
| 88 |  |  | dyu2(1:102,1,iface)=dyu1(ileft(1:102),1,iface) + 2.*dyu1(ileft(1:102),2,iface) + 2.*dyu1(ileft(1:102),3,iface); | 
| 89 |  |  | % Areas | 
| 90 |  |  | for ipnt=0:ratio-1 | 
| 91 |  |  | for jpnt=0:ratio-1 | 
| 92 |  |  | ra2(1:102,1:102,iface)=ra2(1:102,1:102,iface) + ra1(ileft+ipnt,jbott+jpnt,iface); | 
| 93 |  |  | ras2(1:102,2:102,iface)=ras2(1:102,2:102,iface) + ras1(ileft(1:102)+ipnt,jcent(1:101)+jpnt+1,iface); | 
| 94 |  |  | raw2(2:102,1:102,iface)=raw2(2:102,1:102,iface) + raw1(icent(1:101)+ipnt+1,jbott(1:102)+jpnt,iface); | 
| 95 |  |  | raz2(2:102,2:102,iface)=raz2(2:102,2:102,iface) + raz1(icent(1:101)+ipnt+1,jcent(1:101)+jpnt+1,iface); | 
| 96 |  |  | end | 
| 97 |  |  | end | 
| 98 |  |  | % Construct values that we don't know yet - use grid symmetries for edge areas | 
| 99 |  |  | for ipnt=0:ratio-1 | 
| 100 |  |  | ras2(1:102,1,iface) = ras2(1:102,1,iface) + ... | 
| 101 |  |  | ras1(ileft(1:102)+ipnt,1,iface) + 2.*ras1(ileft(1:102)+ipnt,2,iface) + 2.*ras1(ileft(1:102)+ipnt,3,iface); | 
| 102 |  |  | end | 
| 103 |  |  | for jpnt=0:ratio-1 | 
| 104 |  |  | raw2(1,1:102,iface) = raw2(1,1:102,iface) + ... | 
| 105 |  |  | raw1(1,jbott(1:102)+jpnt,iface) + 2.*raw1(2,jbott(1:102)+jpnt,iface) + 2.*raw1(3,jbott(1:102)+jpnt,iface); | 
| 106 |  |  | end | 
| 107 |  |  | % Construct values that we don't know yet - vorticity point area, along face edges than at face corner | 
| 108 |  |  | for jpnt=0:ratio-1 | 
| 109 |  |  | raz2(1,2:102,iface) = raz2(1,2:102,iface) + ... | 
| 110 |  |  | raz1(1,jcent(1:101)+jpnt+1,iface) + 2.*raz1(2,jcent(1:101)+jpnt+1,iface) + 2.*raz1(3,jcent(1:101)+jpnt+1,iface); | 
| 111 |  |  | end | 
| 112 |  |  | for ipnt=0:ratio-1 | 
| 113 |  |  | raz2(2:102,1,iface) = raz2(2:102,1,iface) + ... | 
| 114 |  |  | raz1(icent(1:101)+ipnt+1,1,iface) + 2.*raz1(icent(1:101)+1+ipnt,2,iface) + 2.*raz1(icent(1:101)+1+ipnt,3,iface); | 
| 115 |  |  | end | 
| 116 |  |  | raz2(1,1,iface) = sum(raz1(1:3,1,iface),1) + 2.*sum(raz1(1,2:3,iface),2) + ... | 
| 117 |  |  | 3.*sum(raz1(2:3,3,iface),1) + sum(raz1(2:3,2,iface),1); | 
| 118 |  |  | end | 
| 119 |  |  | % Now the exchanges to fill up the extra column and row | 
| 120 |  |  | % Exchange for lon and lat's at center and corners - no directions (xc is always xc) | 
| 121 |  |  | xc2(1:102,103,[1 3 5])=xc2(1,102:-1:1,[3 5 1]); | 
| 122 |  |  | xc2(103,1:102,[1 3 5])=xc2(1,1:102,[2 4 6]); | 
| 123 |  |  | xc2(1:102,103,[2 4 6])=xc2(1:102,1,[3 5 1]); | 
| 124 |  |  | xc2(103,1:102,[2 4 6])=xc2(102:-1:1,1,[4 6 2]); | 
| 125 |  |  | xc2(103,103,[1 2 3 4 5 6])=xc2(1,1,[3 4 5 6 1 2]); | 
| 126 |  |  | yc2(1:102,103,[1 3 5])=yc2(1,102:-1:1,[3 5 1]); | 
| 127 |  |  | yc2(103,1:102,[1 3 5])=yc2(1,1:102,[2 4 6]); | 
| 128 |  |  | yc2(1:102,103,[2 4 6])=yc2(1:102,1,[3 5 1]); | 
| 129 |  |  | yc2(103,1:102,[2 4 6])=yc2(102:-1:1,1,[4 6 2]); | 
| 130 |  |  | yc2(103,103,[1 2 3 4 5 6])=yc2(1,1,[3 4 5 6 1 2]); | 
| 131 |  |  | xg2(1:102,103,[1 3 5])=xg2(1,102:-1:1,[3 5 1]); | 
| 132 |  |  | xg2(103,1:102,[1 3 5])=xg2(1,1:102,[2 4 6]); | 
| 133 |  |  | xg2(1:102,103,[2 4 6])=xg2(1:102,1,[3 5 1]); | 
| 134 |  |  | xg2(103,1:102,[2 4 6])=xg2(102:-1:1,1,[4 6 2]); | 
| 135 |  |  | xg2(103,103,[1 2 3 4 5 6])=xg2(1,1,[3 4 5 6 1 2]); | 
| 136 |  |  | yg2(1:102,103,[1 3 5])=yg2(1,102:-1:1,[3 5 1]); | 
| 137 |  |  | yg2(103,1:102,[1 3 5])=yg2(1,1:102,[2 4 6]); | 
| 138 |  |  | yg2(1:102,103,[2 4 6])=yg2(1:102,1,[3 5 1]); | 
| 139 |  |  | yg2(103,1:102,[2 4 6])=yg2(102:-1:1,1,[4 6 2]); | 
| 140 |  |  | yg2(103,103,[1 2 3 4 5 6])=yg2(1,1,[3 4 5 6 1 2]); | 
| 141 |  |  | % Exchange for areas at center - no directions | 
| 142 |  |  | ra2(1:102,103,[1 3 5])=ra2(1,102:-1:1,[3 5 1]); | 
| 143 |  |  | ra2(103,1:102,[1 3 5])=ra2(1,1:102,[2 4 6]); | 
| 144 |  |  | ra2(1:102,103,[2 4 6])=ra2(1:102,1,[3 5 1]); | 
| 145 |  |  | ra2(103,1:102,[2 4 6])=ra2(102:-1:1,1,[4 6 2]); | 
| 146 |  |  | ra2(103,103,[1 2 3 4 5 6])=ra2(1,1,[3 4 5 6 1 2]); | 
| 147 |  |  | % Exchange for areas at south and west edges - no directions | 
| 148 |  |  | ras2(1:102,103,[1 3 5])=raw2(1,102:-1:1,[3 5 1]); | 
| 149 |  |  | ras2(103,1:102,[1 3 5])=ras2(1,1:102,[2 4 6]); | 
| 150 |  |  | ras2(1:102,103,[2 4 6])=ras2(1:102,1,[3 5 1]); | 
| 151 |  |  | ras2(103,1:102,[2 4 6])=raw2(102:-1:1,1,[4 6 2]); | 
| 152 |  |  | ras2(103,103,[1 2 3 4 5 6])=ras2(1,1,[3 4 5 6 1 2]); | 
| 153 |  |  | raw2(1:102,103,[1 3 5])=ras2(1,102:-1:1,[3 5 1]); | 
| 154 |  |  | raw2(103,1:102,[1 3 5])=raw2(1,1:102,[2 4 6]); | 
| 155 |  |  | raw2(1:102,103,[2 4 6])=raw2(1:102,1,[3 5 1]); | 
| 156 |  |  | raw2(103,1:102,[2 4 6])=ras2(102:-1:1,1,[4 6 2]); | 
| 157 |  |  | raw2(103,103,[1 2 3 4 5 6])=raw2(1,1,[3 4 5 6 1 2]); | 
| 158 |  |  | % Exchange for areas at vorticity points - no directions | 
| 159 |  |  | % Special case: define upper edge assuming symmetry in y | 
| 160 |  |  | %   Done because the edge values are ambiguously defined | 
| 161 |  |  | %   for quantities defined on a corner | 
| 162 |  |  | %%%%raz2(1:102,103,[1 3 5])=raz2(1,102:-1:1,[3 5 1]); | 
| 163 |  |  | %%%%raz2(1:102,103,[2 4 6])=raz2(1:102,1,[3 5 1]); | 
| 164 |  |  | %%%%raz2(103,1:102,[1 3 5])=raz2(1,1:102,[2 4 6]); | 
| 165 |  |  | %%%%raz2(103,1:102,[2 4 6])=raz2(102:-1:1,1,[4 6 2]); | 
| 166 |  |  | raz2(1:102,103,:)=raz2(1:102,1,:); | 
| 167 |  |  | raz2(103,1:102,:)=raz2(1,1:102,:); | 
| 168 |  |  | raz2(103,103,[1 2 3 4 5 6])=raz2(1,1,[3 4 5 6 1 2]); | 
| 169 |  |  | % | 
| 170 |  |  | % Exchange for dx's and dy's at center - direction but no sign (dx is sometimes from dy) | 
| 171 |  |  | dxf2(1:102,103,[1 3 5])=dyf2(1,102:-1:1,[3 5 1]); | 
| 172 |  |  | dxf2(103,1:102,[1 3 5])=dxf2(1,1:102,[2 4 6]); | 
| 173 |  |  | dxf2(1:102,103,[2 4 6])=dxf2(1:102,1,[3 5 1]); | 
| 174 |  |  | dxf2(103,1:102,[2 4 6])=dyf2(102:-1:1,1,[4 6 2]); | 
| 175 |  |  | dxf2(103,103,[1 2 3 4 5 6])=dxf2(1,1,[3 4 5 6 1 2]); | 
| 176 |  |  | dyf2(1:102,103,[1 3 5])=dxf2(1,102:-1:1,[3 5 1]); | 
| 177 |  |  | dyf2(103,1:102,[1 3 5])=dyf2(1,1:102,[2 4 6]); | 
| 178 |  |  | dyf2(1:102,103,[2 4 6])=dyf2(1:102,1,[3 5 1]); | 
| 179 |  |  | dyf2(103,1:102,[2 4 6])=dxf2(102:-1:1,1,[4 6 2]); | 
| 180 |  |  | dyf2(103,103,[1 2 3 4 5 6])=dyf2(1,1,[3 4 5 6 1 2]); | 
| 181 |  |  | dxg2(1:102,103,[1 3 5])=dyg2(1,102:-1:1,[3 5 1]); | 
| 182 |  |  | dxg2(103,1:102,[1 3 5])=dxg2(1,1:102,[2 4 6]); | 
| 183 |  |  | dxg2(1:102,103,[2 4 6])=dxg2(1:102,1,[3 5 1]); | 
| 184 |  |  | dxg2(103,1:102,[2 4 6])=dyg2(102:-1:1,1,[4 6 2]); | 
| 185 |  |  | dxg2(103,103,[1 2 3 4 5 6])=dxg2(1,1,[3 4 5 6 1 2]); | 
| 186 |  |  | dyg2(1:102,103,[1 3 5])=dxg2(1,102:-1:1,[3 5 1]); | 
| 187 |  |  | dyg2(103,1:102,[1 3 5])=dyg2(1,1:102,[2 4 6]); | 
| 188 |  |  | dyg2(1:102,103,[2 4 6])=dyg2(1:102,1,[3 5 1]); | 
| 189 |  |  | dyg2(103,1:102,[2 4 6])=dxg2(102:-1:1,1,[4 6 2]); | 
| 190 |  |  | dyg2(103,103,[1 2 3 4 5 6])=dyg2(1,1,[3 4 5 6 1 2]); | 
| 191 |  |  | % Exchange for dx's and dy's at edges - direction but no sign (dx is sometimes from dy) | 
| 192 |  |  | dxc2(1:102,103,[1 3 5])=dyc2(1,102:-1:1,[3 5 1]); | 
| 193 |  |  | dxc2(103,1:102,[1 3 5])=dxc2(1,1:102,[2 4 6]); | 
| 194 |  |  | dxc2(1:102,103,[2 4 6])=dxc2(1:102,1,[3 5 1]); | 
| 195 |  |  | dxc2(103,1:102,[2 4 6])=dyc2(102:-1:1,1,[4 6 2]); | 
| 196 |  |  | dxc2(103,103,[1 2 3 4 5 6])=dxc2(1,1,[3 4 5 6 1 2]); | 
| 197 |  |  | dyc2(1:102,103,[1 3 5])=dxc2(1,102:-1:1,[3 5 1]); | 
| 198 |  |  | dyc2(103,1:102,[1 3 5])=dyc2(1,1:102,[2 4 6]); | 
| 199 |  |  | dyc2(1:102,103,[2 4 6])=dyc2(1:102,1,[3 5 1]); | 
| 200 |  |  | dyc2(103,1:102,[2 4 6])=dxc2(102:-1:1,1,[4 6 2]); | 
| 201 |  |  | dyc2(103,103,[1 2 3 4 5 6])=dyc2(1,1,[3 4 5 6 1 2]); | 
| 202 |  |  | % Special case: define upper edge assuming symmetry in y | 
| 203 |  |  | %   Done because the edge values are ambiguously defined | 
| 204 |  |  | %   for quantities defined on a corner | 
| 205 |  |  | %%%%dxv2(1:102,103,[1 3 5])=dyu2(1,102:-1:1,[3 5 1]); | 
| 206 |  |  | %%%%dxv2(103,1:102,[1 3 5])=dxv2(1,1:102,[2 4 6]); | 
| 207 |  |  | %%%%dxv2(1:102,103,[2 4 6])=dxv2(1:102,1,[3 5 1]); | 
| 208 |  |  | %%%%dxv2(103,1:102,[2 4 6])=dyu2(102:-1:1,1,[4 6 2]); | 
| 209 |  |  | %%%%dyu2(1:102,103,[1 3 5])=dxv2(1,102:-1:1,[3 5 1]); | 
| 210 |  |  | %%%%dyu2(103,1:102,[1 3 5])=dyu2(1,1:102,[2 4 6]); | 
| 211 |  |  | %%%%dyu2(1:102,103,[2 4 6])=dyu2(1:102,1,[3 5 1]); | 
| 212 |  |  | %%%%dyu2(103,1:102,[2 4 6])=dxv2(102:-1:1,1,[4 6 2]); | 
| 213 |  |  | dxv2(1:102,103,:)=dxv2(1:102,1,:); | 
| 214 |  |  | dxv2(103,1:102,:)=dxv2(1,1:102,:); | 
| 215 |  |  | dxv2(103,103,[1 2 3 4 5 6])=dxv2(1,1,[3 4 5 6 1 2]); | 
| 216 |  |  | dyu2(1:102,103,:)=dyu2(1:102,1,:); | 
| 217 |  |  | dyu2(103,1:102,:)=dyu2(1,1:102,:); | 
| 218 |  |  | dyu2(103,103,[1 2 3 4 5 6])=dyu2(1,1,[3 4 5 6 1 2]); | 
| 219 |  |  | % | 
| 220 |  |  | % All done with grid values - now compute angles | 
| 221 |  |  | ygarg1=reshape(permute(yg2(1:102,1:102,:),[1 3 2]),[612 102]); | 
| 222 |  |  | racarg2=reshape(permute(ra2(1:102,1:102,:),[1 3 2]),[612 102]); | 
| 223 |  |  | dxgarg3=reshape(permute(dxg2(1:102,1:102,:),[1 3 2]),[612 102]); | 
| 224 |  |  | dygarg4=reshape(permute(dyg2(1:102,1:102,:),[1 3 2]),[612 102]); | 
| 225 |  |  | [AngleCS,AngleSN] = cubeCalcAngle(ygarg1,racarg2,dxgarg3,dygarg4); | 
| 226 |  |  | anglecos2(1:102,1:102,:)=permute(reshape(AngleCS,[102 6 102]),[1 3 2]); | 
| 227 |  |  | anglesin2(1:102,1:102,:)=permute(reshape(AngleSN,[102 6 102]),[1 3 2]); | 
| 228 |  |  | % | 
| 229 |  |  | % And now write it all out | 
| 230 |  |  | for iface=1:6 | 
| 231 |  |  | fid = fopen(['onedegcube.face00' num2str(iface) '.bin'],'w','b'); | 
| 232 |  |  | fwrite(fid,xc2(:,:,iface),'real*8'); | 
| 233 |  |  | fwrite(fid,yc2(:,:,iface),'real*8'); | 
| 234 |  |  | fwrite(fid,dxf2(:,:,iface),'real*8'); | 
| 235 |  |  | fwrite(fid,dyf2(:,:,iface),'real*8'); | 
| 236 |  |  | fwrite(fid,ra2(:,:,iface),'real*8'); | 
| 237 |  |  | fwrite(fid,xg2(:,:,iface),'real*8'); | 
| 238 |  |  | fwrite(fid,yg2(:,:,iface),'real*8'); | 
| 239 |  |  | fwrite(fid,dxv2(:,:,iface),'real*8'); | 
| 240 |  |  | fwrite(fid,dyu2(:,:,iface),'real*8'); | 
| 241 |  |  | fwrite(fid,raz2(:,:,iface),'real*8'); | 
| 242 |  |  | fwrite(fid,dxc2(:,:,iface),'real*8'); | 
| 243 |  |  | fwrite(fid,dyc2(:,:,iface),'real*8'); | 
| 244 |  |  | fwrite(fid,raw2(:,:,iface),'real*8'); | 
| 245 |  |  | fwrite(fid,ras2(:,:,iface),'real*8'); | 
| 246 |  |  | fwrite(fid,dxg2(:,:,iface),'real*8'); | 
| 247 |  |  | fwrite(fid,dyg2(:,:,iface),'real*8'); | 
| 248 |  |  | fwrite(fid,anglecos2(:,:,iface),'real*8'); | 
| 249 |  |  | fwrite(fid,anglesin2(:,:,iface),'real*8'); | 
| 250 |  |  | fclose(fid); | 
| 251 |  |  | end |